diff --git a/.Rbuildignore b/.Rbuildignore
index 4d699692ccac38490cb0fe2a7ef40b98cf33116a..7c1c5baa8f03668a04b8d352020ed07372de4599 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -5,3 +5,6 @@
 ^tests/tmp/
 ^\.gitlab-ci.yml$
 ^\.regressionignore$
+^\.gitlab-ci\.yml$
+^\.vscode$
+^Rplots\.pdf$
diff --git a/.regressionignore b/.regressionignore
index fb307437cc691d1933469ab72fd011cf73874771..1dd20c9de2d828cd149545d7fbc7da4326b40af2 100644
--- a/.regressionignore
+++ b/.regressionignore
@@ -3,9 +3,28 @@
 # The format of this file is: 5 lines of comments followed by one line by
 # ignored variable : [Topic]<SPACE>[Variable].
 # Example for ignoring OutputsModel variable produced by example("RunModel_GR2M"): RunModel_GR2M OutputsModel
+RunModel_GR1A BasinObs
+RunModel_GR1A ConvertFun
+RunModel_GR1A NewTabSeries
+RunModel_GR1A NewTimeFormat
+RunModel_GR1A OutputsModel
+RunModel_GR1A TabSeries
+RunModel_GR1A TimeFormat
+RunModel_GR1A YearFirstMonth
+RunModel_GR2M BasinObs
+RunModel_GR2M ConvertFun
+RunModel_GR2M NewTabSeries
+RunModel_GR2M NewTimeFormat
 RunModel_GR2M OutputsModel
 RunModel_GR2M RunOptions
 RunModel_GR1A OutputsModel
 Calibration_Michel CalibOptions
 Calibration CalibOptions
 CreateCalibOptions CalibOptions
+
+# New version of the SeriesAggreg function
+RunModel_GR2M TabSeries
+RunModel_GR2M TimeFormat
+SeriesAggreg BasinInfo
+SeriesAggreg BasinObs
+SeriesAggreg NewTabSeries
diff --git a/DESCRIPTION b/DESCRIPTION
index 674ba67d234ac3acf322291ce3a81f004811824c..32468507acef67c47ed52d79dd8d14ae743a74cd 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: airGR
 Type: Package
 Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.6.3.78
+Version: 1.6.8.35
 Date: 2021-01-05
 Authors@R: c(
   person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
@@ -20,7 +20,7 @@ Authors@R: c(
   person("Raji", "Pushpalatha", role = c("ctb")),
   person("Audrey", "Valéry", role = c("ctb"))
   )
-Depends: R (>= 3.0.1)
+Depends: R (>= 3.1.0)
 Imports:
   graphics,
   grDevices,
diff --git a/NAMESPACE b/NAMESPACE
index a87a61cd909a547a5ac81dce8ee9ee99564fc734..5d211da10e84c0fd65df7e835c17e466ee96d4a7 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -8,7 +8,11 @@ useDynLib(airGR, .registration = TRUE)
 #####################################
 ##            S3 methods           ##
 #####################################
-S3method("plot", "OutputsModel")
+S3method(plot, OutputsModel)
+S3method(SeriesAggreg, data.frame)
+S3method(SeriesAggreg, list)
+S3method(SeriesAggreg, InputsModel)
+S3method(SeriesAggreg, OutputsModel)
 
 
 
@@ -58,11 +62,9 @@ export(TransfoParam_GR4J)
 export(TransfoParam_GR5J)
 export(TransfoParam_GR6J)
 export(TransfoParam_Lag)
-export(plot.OutputsModel)
 export(.ErrorCrit)
 
 
-
 #####################################
 ##               Import            ##
 #####################################
diff --git a/NEWS.md b/NEWS.md
index 44e9557197e3ac98962fbc96dcaa7a21591b6258..e05beedc8ab353471e66308793fd17e105b16815 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,6 +2,19 @@
 
 
 
+### 1.6.8.34 Release Notes (2021-01-05)
+
+#### New features
+
+- Added <code>SeriesAggreg</code> S3 method with functions for `InputsModel`, `OutputsModel`, `list`, `data.frame` class objects. This new version of the <code>SeriesAggreg()</code> function also allows to compute regimes. ([#25](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/25))
+- Added <code>.GetAggregConvertFun()</code> private function in order to choose automatically the <code>ConvertFun</code> to apply on each element of objects used in <code>SeriesAggreg.InputsModel</code> and <code>SeriesAggreg.OutputsModel</code>.
+- Added <code>.AggregConvertFunTable</code> data.frame that allows the user to see what names of list items or data.frame column names are guessed and eventually customise this correspondence table.
+
+#### Bug fixes
+
+- <code>TimeLag</code> of the <code>SeriesAggreg()</code> function now runs when <code>TimeLag >= 3600</code>. ([#41](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/41))
+
+____________________________________________________________________________________
 
 
 ### 1.6.3.78 Release Notes (2021-01-05)
@@ -28,9 +41,7 @@
 #### Major user-visible changes
 
 - Added output to <code>RunModel_GR2M()</code> function (Ps). ([#51](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/51))
-
 - <code>PE_Oudin()</code> can now run for several locations (i.e. several latitudes) in the Fortran mode (<code>RunFortran = TRUE</code>). In this case <code>Lat</code> must be of the same length as <code>Temp</code>. ([#62](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/62))
-
 - <code>RunModel()</code> now allows to run semi-distributed GR models. ([#34](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/34))
 
 
@@ -779,5 +790,3 @@ ________________________________________________________________________________
 - <code>CalibrationAlgo_nlminb_stats</code> argument was wrongly defined in <code>DefineFunctions_CalibrationAlgo()</code> (<code>optim</code> instead of <code>nlminb</code>).
 
 - Format checking for <code>RunOptions</code> was incorrectly made in <code>CheckArg()</code> function.
-
-
diff --git a/R/Imax.R b/R/Imax.R
index 41768866b7190655b8910de34ec7e7b3c2a4b670..3c6058f1b1c4945ee734b2a50a185c5ae3d9f8c5 100644
--- a/R/Imax.R
+++ b/R/Imax.R
@@ -32,7 +32,7 @@ Imax <- function(InputsModel,
   TabSeries <- data.frame(DatesR = InputsModel$DatesR[IndPeriod_Run], 
                           Precip = InputsModel$Precip[IndPeriod_Run], 
                           PotEvap = InputsModel$PotEvap[IndPeriod_Run])
-  daily_data <- SeriesAggreg(TabSeries, "hourly", "daily", 
+  daily_data <- SeriesAggreg(TabSeries, Format = "%Y%m%d", 
                              ConvertFun = c("sum", "sum"))
   
   ##calculate total interception of daily GR models on the period
diff --git a/R/SeriesAggreg.InputsModel.R b/R/SeriesAggreg.InputsModel.R
new file mode 100644
index 0000000000000000000000000000000000000000..a6414fb823235d3fb19904602adcd8a177b82466
--- /dev/null
+++ b/R/SeriesAggreg.InputsModel.R
@@ -0,0 +1,6 @@
+SeriesAggreg.InputsModel <- function(x, ...) {
+  SeriesAggreg.list(x,
+                    ConvertFun = .GetAggregConvertFun(names(x)),
+                    except = c("ZLayers", "LengthHydro", "BasinAreas"),
+                    ...)
+}
diff --git a/R/SeriesAggreg.OutputsModel.R b/R/SeriesAggreg.OutputsModel.R
new file mode 100644
index 0000000000000000000000000000000000000000..264bd5b7fb7f41e226004ce244d545246eefe881
--- /dev/null
+++ b/R/SeriesAggreg.OutputsModel.R
@@ -0,0 +1,6 @@
+SeriesAggreg.OutputsModel <- function(x, ...) {
+  SeriesAggreg.list(x,
+                    ConvertFun = .GetAggregConvertFun(names(x)),
+                    except = "StateEnd",
+                    ...)
+}
diff --git a/R/SeriesAggreg.R b/R/SeriesAggreg.R
index 209afb36de7f5ce2d05b6d93a08a74aa0d25550d..737ce62eb7840cc38b59fad260c1d07bf150fc45 100644
--- a/R/SeriesAggreg.R
+++ b/R/SeriesAggreg.R
@@ -1,174 +1,3 @@
-SeriesAggreg <- function(TabSeries,
-                         TimeFormat, NewTimeFormat,
-                         ConvertFun,
-                         YearFirstMonth = 1, TimeLag = 0,
-                         verbose = TRUE) {
-  
-  
-  ##_Arguments_check
-  
-  ##check_TabSeries
-  if (!is.data.frame(TabSeries)) {
-    stop("'TabSeries' must be a data.frame containing the dates and data to be aggregated")
-  }
-  if (ncol(TabSeries) < 2) {
-    stop("'TabSeries' must contain at least two columns (including the column of dates)")
-  }
-  ##check_TimeFormat
-  if (!any(class(TabSeries[, 1]) %in% "POSIXt")) {
-    stop("'TabSeries' first column must be a vector of class 'POSIXlt' or 'POSIXct'")
-  }
-  if (any(class(TabSeries[, 1]) %in% "POSIXlt")) {
-    TabSeries[, 1] <- as.POSIXct(TabSeries[, 1])
-  }
-  for (iCol in 2:ncol(TabSeries)) {
-    if (!is.numeric(TabSeries[,iCol])) {
-      stop("'TabSeries' columns (other than the first one) be of numeric class")
-    }
-  }
-  ##check TimeFormat and NewTimeFormat
-  TimeStep <- c("hourly", "daily", "monthly", "yearly")
-  TimeFormat    <- match.arg(TimeFormat   , choices = TimeStep)
-  NewTimeFormat <- match.arg(NewTimeFormat, choices = TimeStep)
-  ##check_ConvertFun
-  if (!is.vector(ConvertFun)) {
-    stop("'ConvertFun' must be a vector of character")
-  }
-  if (!is.character(ConvertFun)) {
-    stop("'ConvertFun' must be a vector of character")
-  }
-  if (length(ConvertFun) != (ncol(TabSeries) - 1)) {
-    stop(
-      paste("'ConvertFun' must be of length", ncol(TabSeries) - 1, "(length=ncol(TabSeries)-1)")
-    )
-  }
-  if (sum(ConvertFun %in% c("sum", "mean") == FALSE) != 0) {
-    stop("'ConvertFun' elements must be one of 'sum' or 'mean'")
-  }
-  ##check_YearFirstMonth
-  if (!is.vector(YearFirstMonth)) {
-    stop("'YearFirstMonth' must be an integer between 1 and 12")
-  }
-  if (!is.numeric(YearFirstMonth)) {
-    stop("'YearFirstMonth' must be an integer between 1 and 12")
-  }
-  YearFirstMonth <- as.integer(YearFirstMonth)
-  if (length(YearFirstMonth) != 1) {
-    stop("'YearFirstMonth' must be only one integer between 1 and 12")
-  }
-  if (YearFirstMonth %in% (1:12) == FALSE) {
-    stop("'YearFirstMonth' must be only one integer between 1 and 12")
-  } 
-  ##check_DatesR_integrity
-  by <- switch(TimeFormat,
-               'hourly'  = "hours",
-               'daily'   = "days",
-               'monthly' = "months",
-               'yearly'  = "years")
-  TmpDatesR <- seq(from = TabSeries[1, 1], to = tail(TabSeries[, 1], 1), by = by)
-  if (!identical(TabSeries[, 1], TmpDatesR)) {
-    stop("some dates might not be ordered or are missing in 'TabSeries'")
-  }
-  ##check_conversion_direction
-  if ((TimeFormat == "daily"   & NewTimeFormat %in% c("hourly")                  ) |
-      (TimeFormat == "monthly" & NewTimeFormat %in% c("hourly","daily")          ) |
-      (TimeFormat == "yearly"  & NewTimeFormat %in% c("hourly","daily","monthly"))) { 
-    stop("only time aggregation can be performed")
-  } 
-  ##check_if_conversion_not_needed
-  if ((TimeFormat == "hourly"  & NewTimeFormat == "hourly" ) |
-      (TimeFormat == "daily"   & NewTimeFormat == "daily"  ) |
-      (TimeFormat == "monthly" & NewTimeFormat == "monthly") |
-      (TimeFormat == "yearly"  & NewTimeFormat == "yearly" )) { 
-    if (verbose) {
-      warning("the old and new format are identical \n\t -> no time-step conversion was performed")
-      return(TabSeries)
-    }
-  }
-  
-  
-  ##_Time_step_conversion
-  
-  ##_Handle_conventional_difference_between_hourly_series_and_others
-  TmpDatesR <- TabSeries[, 1]
-  #if (TimeFormat=="hourly") { TmpDatesR <- TmpDatesR - 60*60; }
-  TmpDatesR <- TmpDatesR + TimeLag
-  Hmax <- "00"
-  if (TimeFormat == "hourly") {
-    Hmax <- "23"
-  }
-  
-  ##_Identify_the_part_of_the_series_to_be_aggregated
-  NDaysInMonth <- list("31", c("28", "29"), "31", "30", "31", "30", "31", "31", "30", "31", "30", "31")
-  NDaysAndMonth <- sprintf("%02i%s", c(1:2, 2:12), unlist(NDaysInMonth))
-  YearLastMonth <- YearFirstMonth + 11
-  if (YearLastMonth > 12) {
-    YearLastMonth <- YearLastMonth - 12
-  }
-  YearFirstMonthTxt <- formatC(YearFirstMonth, format = "d", width = 2, flag = "0")
-  YearLastMonthTxt  <- formatC(YearLastMonth , format = "d", width = 2, flag = "0")
-  if (NewTimeFormat == "daily") {
-    Ind1 <- which(format(TmpDatesR,    "%H") == "00")
-    Ind2 <- which(format(TmpDatesR,    "%H") == Hmax)
-    if (Ind2[1] < Ind1[1]) {
-      Ind2 <- Ind2[2:length(Ind2)]
-    }
-    if (tail(Ind1, 1) > tail(Ind2, 1)) {
-      Ind1 <- Ind1[1:(length(Ind1) - 1)]
-    }
-    ### Aggr <- NULL; iii <- 0; for(kkk in 1:length(Ind1)) {
-    ### iii <- iii+1; Aggr <- c(Aggr,rep(iii,length(Ind1[kkk]:Ind2[kkk]))); }
-    Aggr <- as.numeric(format(TmpDatesR[min(Ind1):max(Ind2)], "%Y%m%d"))
-    ### more efficient
-    NewDatesR <- data.frame(seq(from = TmpDatesR[min(Ind1)], to = TmpDatesR[max(Ind2)], by = "days"))
-  }
-  if (NewTimeFormat=="monthly") {
-    Ind1 <- which(format(TmpDatesR,  "%d%H") == "0100")
-    Ind2 <- which(format(TmpDatesR,"%m%d%H") %in% paste0(NDaysAndMonth, Hmax))
-    Ind2[1:(length(Ind2) - 1)][diff(Ind2) == 1] <- NA
-    Ind2 <- Ind2[!is.na(Ind2)] ### to keep only feb 29 if both feb 28 and feb 29 exists
-    if (Ind2[1] < Ind1[1]) {
-      Ind2 <- Ind2[2:length(Ind2)]
-    }
-    if (tail(Ind1, 1) > tail(Ind2, 1)) {
-      Ind1 <- Ind1[1:(length(Ind1) - 1)]
-    }
-    ### Aggr <- NULL; iii <- 0; for(kkk in 1:length(Ind1)) { 
-    ### iii <- iii+1; Aggr <- c(Aggr,rep(iii,length(Ind1[kkk]:Ind2[kkk]))); }
-    Aggr <- as.numeric(format(TmpDatesR[min(Ind1):max(Ind2)],"%Y%m"));  ### more efficient
-    NewDatesR <- data.frame(seq(from=TmpDatesR[min(Ind1)],to=TmpDatesR[max(Ind2)],by="months"))
-  }
-  if (NewTimeFormat == "yearly") {
-    Ind1 <- which(format(TmpDatesR, "%m%d%H") %in% paste0(YearFirstMonthTxt, "0100"))
-    Ind2 <- which(format(TmpDatesR, "%m%d%H") %in% paste0(YearLastMonthTxt, NDaysInMonth[[YearLastMonth]], Hmax))
-    Ind2[1:(length(Ind2) - 1)][diff(Ind2) == 1] <- NA
-    Ind2 <- Ind2[!is.na(Ind2)]
-    ### to keep only feb 29 if both feb 28 and feb 29 exists
-    if (Ind2[1] < Ind1[1]) {
-      Ind2 <- Ind2[2:length(Ind2)]
-    }
-    if (tail(Ind1, 1) > tail(Ind2, 1)) {
-      Ind1 <- Ind1[1:(length(Ind1) - 1)]
-    }
-    Aggr <- NULL
-    iii <- 0
-    for (kkk in 1:length(Ind1)) {
-      iii <- iii + 1
-      Aggr <- c(Aggr, rep(iii, length(Ind1[kkk]:Ind2[kkk])))
-    }
-    ### Aggr <- as.numeric(format(TmpDatesR[min(Ind1):max(Ind2)],"%Y")); ### not working if YearFirstMonth != 01
-    NewDatesR <- data.frame(seq(from = TmpDatesR[min(Ind1)], to = TmpDatesR[max(Ind2)], by = "years"))
-  }
-  ##_Aggreation_and_export
-  NewTabSeries <- data.frame(NewDatesR)
-  for (iCol in 2:ncol(TabSeries)) {
-    AggregData <- aggregate(TabSeries[min(Ind1):max(Ind2), iCol],
-                            by = list(Aggr),
-                            FUN = ConvertFun[iCol - 1], na.rm = FALSE)[, 2]
-    NewTabSeries <- data.frame(NewTabSeries, AggregData)
-  }
-  names(NewTabSeries) <- names(TabSeries)
-  return(NewTabSeries)
-  
-  
-}
\ No newline at end of file
+SeriesAggreg <- function(x, Format, ...) {
+  UseMethod("SeriesAggreg")
+}
diff --git a/R/SeriesAggreg.data.frame.R b/R/SeriesAggreg.data.frame.R
new file mode 100644
index 0000000000000000000000000000000000000000..68f2d625382d175a4c406e0cc44621ceb729270a
--- /dev/null
+++ b/R/SeriesAggreg.data.frame.R
@@ -0,0 +1,164 @@
+SeriesAggreg.data.frame <- function(x,
+                                    Format,
+                                    ConvertFun,
+                                    TimeFormat = NULL,
+                                    NewTimeFormat = NULL,
+                                    YearFirstMonth = 1,
+                                    TimeLag = 0,
+                                    ...) {
+  ## Arguments checks
+  if (!is.null(TimeFormat)) {
+    warning("deprecated 'TimeFormat' argument", call. = FALSE)
+  }
+  if (missing(Format)) {
+    Format <- .GetSeriesAggregFormat(NewTimeFormat)
+  } else if (!is.null(NewTimeFormat)) {
+    warning("deprecated 'NewTimeFormat' argument: 'Format' argument is used instead",
+            call. = FALSE)
+  }
+  if (is.null(Format)) {
+    stop("argument 'Format' is missing")
+  }
+
+  ## check x
+  if (!is.data.frame(x)) {
+    stop("'x' must be a data.frame containing the dates and data to be aggregated")
+  }
+  if (ncol(x) < 2) {
+    stop("'x' must contain at least two columns (including the column of dates)")
+  }
+  ## check x date column
+  if (!inherits(x[[1L]], "POSIXt")) {
+    stop("'x' first column must be a vector of class 'POSIXlt' or 'POSIXct'")
+  }
+  if (inherits(x[[1L]], "POSIXlt")) {
+    x[[1L]] <- as.POSIXct(x[[1L]])
+  }
+  ## check x other columns (boolean converted to numeric)
+  apply(x[, -1L, drop = FALSE],
+        MARGIN = 2,
+        FUN = function(iCol) {
+          if (!is.numeric(iCol)) {
+            stop("'x' columns (other than the first one) must be of numeric class")
+          }
+        })
+  ## check Format
+  listFormat <- c("%Y%m%d", "%Y%m", "%Y", "%m", "%d")
+  Format <- gsub(pattern = "[[:punct:]]+", replacement = "%", Format)
+  Format <- match.arg(Format, choices = listFormat)
+
+  ## check ConvertFun
+  listConvertFun <- list(sum = sum, mean = mean)
+  ConvertFun <- names(listConvertFun)[match(ConvertFun, names(listConvertFun))]
+  if (anyNA(ConvertFun)) {
+    stop("'ConvertFun' should be a one of 'sum' or 'mean'")
+  }
+  if (length(ConvertFun) != (ncol(x) - 1)) {
+    stop(sprintf("'ConvertFun' must be of length %i (ncol(x)-1)", ncol(x) - 1))
+  }
+  ## check YearFirstMonth
+  msgYearFirstMonth <- "'YearFirstMonth' should be a single vector of numeric value between 1 and 12"
+  YearFirstMonth <- match(YearFirstMonth, 1:12)
+  if (anyNA(YearFirstMonth)) {
+    stop(msgYearFirstMonth)
+  }
+  if (length(YearFirstMonth) != 1) {
+    stop(msgYearFirstMonth)
+  }
+  if (YearFirstMonth != 1 & Format != "%Y") {
+    warning("'YearFirstMonth' is ignored because  Format != '%Y'")
+  }
+  ## check TimeLag
+  msgTimeLag <- "'TimeLag' should be a single vector of a positive numeric value"
+  if (!is.vector(TimeLag)) {
+    stop(msgTimeLag)
+  }
+  if (!is.numeric(TimeLag)) {
+    stop(msgTimeLag)
+  }
+  if (length(TimeLag) != 1 | !any(TimeLag >= 0)) {
+    stop(msgTimeLag)
+  }
+
+  TabSeries0 <- x
+  colnames(TabSeries0)[1L] <- "DatesR"
+  TabSeries0$DatesR <- TabSeries0$DatesR + TimeLag
+
+  TabSeries2 <- TabSeries0
+
+  if (!Format %in% c("%d", "%m")) {
+    start <- sprintf("%i-01-01 00:00:00",
+                     as.numeric(format(TabSeries2$DatesR[1L], format = "%Y")) - 1)
+    stop  <- sprintf("%i-12-31 00:00:00",
+                     as.numeric(format(TabSeries2$DatesR[nrow(TabSeries2)], format = "%Y")) + 1)
+    unitTs <- format(diff(x[1:2, 1]))
+    if (gsub("[0-9]+ ", "", unitTs) == "hours") {
+      byTs <- "hours"
+    } else {
+      if (gsub(" days$", "", unitTs) == "1") {
+        byTs <- "days"
+      } else {
+        byTs <- "months"
+      }
+    }
+    fakeTs <- data.frame(DatesR = seq(from = as.POSIXct(start, tz = "UTC"),
+                                      to   = as.POSIXct(stop , tz = "UTC"),
+                                      by   = byTs) + TimeLag)
+    TabSeries2 <- merge(fakeTs, TabSeries2, by = "DatesR", all.x = TRUE)
+  }
+  TabSeries2$DatesRini <- TabSeries2$DatesR %in% TabSeries0$DatesR
+
+
+  TabSeries2$Selec2 <- format(TabSeries2$DatesR, Format)
+
+  if (nchar(Format) > 2 | Format == "%Y") {
+    # Compute aggregation
+    TabSeries2$Selec <- !duplicated(TabSeries2$Selec2)
+    if (all(TabSeries2$Selec)) {
+      warning("the requested time 'Format' is the same as the one in 'x'. No time-step conversion was performed")
+      return(x)
+    }
+    if (Format == "%Y") {
+      yfm <- sprintf("%02.f", YearFirstMonth)
+      spF1 <- "%m"
+      spF2 <- "%Y-%m"
+      TabSeries2$Selec1 <- format(TabSeries2$DatesR, spF1)
+      TabSeries2$Selec2 <- format(TabSeries2$DatesR, spF2)
+      TabSeries2$Selec <- !duplicated(TabSeries2$Selec2) & TabSeries2$Selec1 == yfm
+    }
+    TabSeries2$Fac2 <- cumsum(TabSeries2$Selec)
+  } else {
+    # Compute regime
+    if (Format == "%d") {
+      spF2 <- "%m-%d"
+      TabSeries2$Selec2 <- format(TabSeries2$DatesR, spF2)
+    }
+    TabSeries2$Fac2 <- TabSeries2$Selec2
+    TabSeries2$Selec <- !duplicated(TabSeries2$Selec2)
+    ConvertFun <- rep("mean", ncol(x) - 1)
+  }
+
+  listTsAggreg <- lapply(names(listConvertFun), function(y) {
+    if (any(ConvertFun == y)) {
+      colTsAggreg <- c("Fac2", colnames(x)[-1L][ConvertFun == y])
+      aggregate(. ~ Fac2,
+                data = TabSeries2[, colTsAggreg],
+                FUN = listConvertFun[[y]],
+                na.action = na.pass)
+    } else {
+      NULL
+    }
+  })
+  listTsAggreg <- listTsAggreg[!sapply(listTsAggreg, is.null)]
+  tsAggreg <- do.call(cbind, listTsAggreg)
+  tsAggreg <- tsAggreg[, !duplicated(colnames(tsAggreg))]
+  tsAggreg <- merge(tsAggreg,
+                    TabSeries2[, c("Fac2", "DatesR", "DatesRini", "Selec")],
+                    by = "Fac2",
+                    all.x = TRUE,
+                    all.y = FALSE)
+  tsAggreg <- tsAggreg[tsAggreg$Selec & tsAggreg$DatesRini, ]
+  tsAggreg <- tsAggreg[, colnames(TabSeries0)]
+  return(tsAggreg)
+
+}
diff --git a/R/SeriesAggreg.list.R b/R/SeriesAggreg.list.R
new file mode 100644
index 0000000000000000000000000000000000000000..8834347dbc5ac40fa35a5c19b074699173c34f15
--- /dev/null
+++ b/R/SeriesAggreg.list.R
@@ -0,0 +1,152 @@
+SeriesAggreg.list <- function(x,
+                              Format,
+                              ConvertFun,
+                              NewTimeFormat = NULL,
+                              simplify = FALSE,
+                              except = NULL,
+                              recursive = TRUE,
+                              ...) {
+  if (missing(Format)) {
+    Format <- .GetSeriesAggregFormat(NewTimeFormat)
+  } else if (!is.null(NewTimeFormat)) {
+    warning("deprecated 'NewTimeFormat' argument: 'Format' argument is used instead",
+            call. = FALSE)
+  }
+
+  # Determination of DatesR
+  if (!is.null(x$DatesR)) {
+    if (!inherits(x$DatesR, "POSIXt")) {
+      stop("'x$DatesR' should be of class 'POSIXt'")
+    }
+    DatesR <- x$DatesR
+  } else {
+    # Auto-detection of POSIXt item in Tabseries
+    itemPOSIXt <- which(sapply(x, function(x) {
+      inherits(x, "POSIXt")
+    }, simplify = TRUE))[1]
+    if (is.na(itemPOSIXt)) {
+      stop("At least one item of argument 'x' should be of class 'POSIXt'")
+    }
+    warning("Item 'DatesR' not found in 'x' argument: the item ",
+            names(x)[itemPOSIXt],
+            " has been automatically chosen")
+    DatesR <- x[[names(x)[itemPOSIXt]]]
+  }
+
+  # Selection of numeric items for aggregation
+  numericCols <- names(which(sapply(x, inherits, "numeric")))
+  arrayCols <- names(which(sapply(x, inherits, "array")))
+  numericCols <- setdiff(numericCols, arrayCols)
+  if (!is.null(except)) {
+    if (!inherits(except, "character")) {
+      stop("Argument 'except' should be a 'character'")
+    }
+    numericCols <- setdiff(numericCols, except)
+  }
+
+  cols <- x[numericCols]
+  lengthCols <- sapply(cols, length, simplify = TRUE)
+  if (any(lengthCols != length(DatesR))) {
+    sErr <- paste0(names(lengthCols)[lengthCols != length(DatesR)],
+                   " (", lengthCols[lengthCols != length(DatesR)], ")",
+                   collapse = ", ")
+    warning("The length of the following `numeric` items in 'x' ",
+            "is different than the length of 'DatesR (",
+            length(DatesR),
+            ")': they will be ignored in the aggregation: ",
+            sErr)
+    cols <- cols[lengthCols == length(DatesR)]
+  }
+  dfOut <- NULL
+  if (length(cols)) {
+    ConvertFun2 <- .GetAggregConvertFun(names(cols))
+    if (is.null(recursive)) {
+      if (missing(ConvertFun)) {
+        stop("'ConvertFun' argument should provided if 'recursive = NULL'")
+      } else if (!is.na(ConvertFun)) {
+        ConvertFun2 <- rep(ConvertFun, length(cols))
+      }
+    }
+    dfOut <- SeriesAggreg(cbind(DatesR, as.data.frame(cols)),
+                          Format,
+                          ...,
+                          ConvertFun = ConvertFun2)
+  }
+
+  if (simplify) {
+    # Returns data.frame of numeric found in the first level of the list
+    return(dfOut)
+
+  } else {
+    res <- list()
+    # Convert aggregated data.frame into list
+    if (!is.null(dfOut)) {
+      res <- as.list(dfOut)
+      ## To be consistent with InputsModel class and because plot.OutputsModel use the POSIXlt class
+      res$DatesR <- as.POSIXlt(res$DatesR)
+    }
+
+    # Exploration of embedded lists and data.frames
+    if (is.null(recursive) || recursive) {
+      listCols <- x[sapply(x, inherits, "list")]
+      dfCols <- x[sapply(x, inherits, "data.frame")]
+      dfCols <- c(dfCols, x[sapply(x, inherits, "matrix")])
+      listCols <- listCols[setdiff(names(listCols), names(dfCols))]
+      if (length(listCols) > 0) {
+        # Check for predefined ConvertFun for all sub-elements
+        ConvertFun <- .GetAggregConvertFun(names(listCols))
+        # Run SeriesAggreg for each embedded list
+        listRes <- lapply(names(listCols), function(x) {
+          listCols[[x]]$DatesR <- DatesR
+          SeriesAggreg(listCols[[x]],
+                       Format = Format,
+                       except = except,
+                       ConvertFun = ConvertFun[x],
+                       recursive = NULL,
+                       ...)
+        })
+        names(listRes) <- names(listCols)
+        if (is.null(res$DatesR)) {
+          # Copy DatesR in top level list
+          res$DatesR <- listRes[[1]]$DatesR
+        }
+        # Remove DatesR in embedded lists
+        lapply(names(listRes), function(x) {
+          listRes[[x]]$DatesR <<- NULL
+        })
+        res <- c(res, listRes)
+      }
+      if (length(dfCols) > 0) {
+        # Processing matrix and dataframes
+        for (i in length(dfCols)) {
+          key <- names(dfCols)[i]
+          m <- dfCols[[i]]
+          if (nrow(m) != length(DatesR)) {
+            warning(
+              "The number of rows of the 'matrix' item ",
+              key, " (", nrow(m),
+              ") is different than the length of 'DatesR ('", length(DatesR),
+              "), it will be ignored in the aggregation"
+            )
+          } else {
+            ConvertFun <- rep(.GetAggregConvertFun(key), ncol(m))
+            res[[key]] <- SeriesAggreg(data.frame(DatesR, m),
+                                       Format = Format,
+                                       ConvertFun = ConvertFun)
+          }
+        }
+      }
+    }
+
+    # Store elements that are not aggregated
+    res <- c(res, x[setdiff(names(x), names(res))])
+
+    class(res) <- gsub("hourly|daily|monthly|yearly",
+                       .GetSeriesAggregClass(Format),
+                       class(x))
+
+    return(res)
+
+  }
+
+}
diff --git a/R/Utils.R b/R/Utils.R
index 253ecb3f69d8f8116ecad3bf42e506dc080564c4..d1cde5eb0dc10d9f3f686ddf03dab9c24671d850 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -11,16 +11,15 @@
 #   }
 # }
 
-
 ## =================================================================================
 ## function to manage Fortran outputs
 ## =================================================================================
 
 .FortranOutputs <- function(GR = NULL, isCN = FALSE) {
-  
+
   outGR <- NULL
   outCN <- NULL
-  
+
   if (is.null(GR)) {
     GR <- ""
   }
@@ -38,7 +37,7 @@
                "AE", "EI", "ES",
                "Perc", "PR",
                "Q9", "Q1",
-               "Rout", "Exch", 
+               "Rout", "Exch",
                "AExch1", "AExch2",
                "AExch", "QR",
                "QD",
@@ -48,7 +47,7 @@
                "AE",
                "Perc", "PR",
                "Q9", "Q1",
-               "Rout", "Exch", 
+               "Rout", "Exch",
                "AExch1", "AExch2",
                "AExch", "QR",
                "QD",
@@ -66,12 +65,12 @@
                "Qsim")
   }
   if (isCN) {
-    outCN <- c("Pliq", "Psol", 
-               "SnowPack", "ThermalState", "Gratio", 
-               "PotMelt", "Melt", "PliqAndMelt", "Temp", 
+    outCN <- c("Pliq", "Psol",
+               "SnowPack", "ThermalState", "Gratio",
+               "PotMelt", "Melt", "PliqAndMelt", "Temp",
                "Gthreshold", "Glocalmax")
   }
-  
+
   res <- list(GR = outGR, CN = outCN)
-  
+
 }
diff --git a/R/UtilsSeriesAggreg.R b/R/UtilsSeriesAggreg.R
new file mode 100644
index 0000000000000000000000000000000000000000..94d7e181b9e4e37cc9c2fffc9d478b94c464348d
--- /dev/null
+++ b/R/UtilsSeriesAggreg.R
@@ -0,0 +1,57 @@
+.GetSeriesAggregFormat <- function(NewTimeFormat) {
+  errNewTimeFormat <- FALSE
+  if (missing(NewTimeFormat)) {
+    errNewTimeFormat <- TRUE
+  } else if (is.null(NewTimeFormat)) {
+    errNewTimeFormat <- TRUE
+  }
+  if (errNewTimeFormat) {
+    stop("Argument `Format` is missing")
+  }
+  if (!is.null(NewTimeFormat)) {
+    TimeStep <- c("hourly", "daily", "monthly", "yearly")
+    NewTimeFormat <- match.arg(NewTimeFormat, choices = TimeStep)
+    Format <- switch(NewTimeFormat,
+                     hourly  = "%Y%m%d%h",
+                     daily   = "%Y%m%d",
+                     monthly = "%Y%m",
+                     yearly  = "%Y")
+    msgNewTimeFormat <- sprintf("'Format' automatically set to %s", sQuote(Format))
+    warning("deprecated 'NewTimeFormat' argument: please use 'Format' instead.",
+            msgNewTimeFormat,
+            call. = FALSE)
+    return(Format)
+  }
+  return(NULL)
+}
+
+.GetSeriesAggregClass <- function(Format) {
+  Format <- substr(Format,
+                   start = nchar(Format),
+                   stop = nchar(Format))
+  switch(Format,
+         h = "hourly",
+         d = "daily",
+         m = "monthly",
+         Y = "yearly")
+}
+
+.GetAggregConvertFun <- function(Outputs) {
+  AggregConvertFunTable <- rbind(
+    data.frame(ConvertFun = "mean",
+              Outputs = c("Prod","Rout","Exp","SnowPack","ThermalState",
+                          "Gratio","Temp","Gthreshold","Glocalmax","LayerTempMean", "T"),
+              stringsAsFactors = FALSE), # R < 4.0 compatibility: avoids mixing numeric and factor into numeric
+    data.frame(ConvertFun = "sum",
+              Outputs = c("PotEvap","Precip","Pn","Ps","AE","Perc","PR","Q9",
+                          "Q1","Exch","AExch1","AExch2","AExch","QR","QRExp",
+                          "QD","Qsim","Pliq","Psol","PotMelt","Melt","PliqAndMelt",
+                          "LayerPrecip","LayerFracSolidPrecip", "Qmm", "Qls", "E", "P", "Qupstream"),
+              stringsAsFactors = FALSE) # R < 4.0 compatibility: avoids mixing numeric and factor into numeric
+  )
+  res <- sapply(Outputs, function(iOutputs) {
+    iRes <- AggregConvertFunTable$ConvertFun[AggregConvertFunTable$Outputs == iOutputs]
+    iRes <- ifelse(test = any(is.na(iRes)), yes = NA, no = iRes) # R < 4.0 compatibility
+  })
+  return(res)
+}
diff --git a/R/plot.OutputsModel.R b/R/plot.OutputsModel.R
index 8a52f9d30e7913c60c21805a5917115edddc1a74..619e5ea9618613258f6f2d0928361468b89a29bd 100644
--- a/R/plot.OutputsModel.R
+++ b/R/plot.OutputsModel.R
@@ -159,23 +159,19 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       return(filter(x, filter = rep(1 / n, n), sides = 2, circular = TRUE))
     }    
     BOOL_TS  <- FALSE
-    TimeStep <- difftime(tail(OutputsModel$DatesR, 1), tail(OutputsModel$DatesR, 2), units = "secs")[[1]]
-    if (inherits(OutputsModel, "hourly") &
-        TimeStep %in% (60 * 60)) {
+    if (inherits(OutputsModel, "hourly")) {
       BOOL_TS <- TRUE
       NameTS <- "hour"
       plotunit <- "[mm/h]"
       formatAxis <- "%m/%Y"
     }
-    if (inherits(OutputsModel, "daily") &
-        TimeStep %in% (24 * 60 * 60)) {
+    if (inherits(OutputsModel, "daily")) {
       BOOL_TS <- TRUE
       NameTS <- "day"
       plotunit <- "[mm/d]"
       formatAxis <- "%m/%Y"
     }
-    if (inherits(OutputsModel, "monthly") &
-        TimeStep %in% (c(28, 29, 30, 31) * 24 * 60 * 60)) {
+    if (inherits(OutputsModel, "monthly")) {
       BOOL_TS <- TRUE
       NameTS <- "month"
       plotunit <- "[mm/month]"
@@ -184,16 +180,15 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
         OutputsModel$DatesR <- as.POSIXlt(format(OutputsModel$DatesR, format = "%Y-%m-01"), tz = "UTC", format = "%Y-%m-%d")
       }
     }
-    if (inherits(OutputsModel, "yearly") &
-        TimeStep %in% (c(365, 366) * 24 * 60 * 60)) {
+    if (inherits(OutputsModel, "yearly")) {
       BOOL_TS <- TRUE
       NameTS <- "year"
       plotunit <- "[mm/y]"
       formatAxis <- "%Y"
     }
-    if (!BOOL_TS) {
-      stop("the time step of the model inputs could not be found")
-    }
+    # if (!BOOL_TS) {
+    #   stop("the time step of the model inputs could not be found")
+    # }
   }
   if (length(IndPeriod_Plot) == 0) {
     IndPeriod_Plot <- 1:length(OutputsModel$DatesR)
diff --git a/airGR.Rproj b/airGR.Rproj
index b1ded457422a5fb84fd0f19adb554b7b06a2a6a5..398aa1438c2c45f6eb3efd75ef2d4052d5ae7923 100644
--- a/airGR.Rproj
+++ b/airGR.Rproj
@@ -12,6 +12,9 @@ Encoding: UTF-8
 RnwWeave: knitr
 LaTeX: pdfLaTeX
 
+AutoAppendNewline: Yes
+StripTrailingWhitespace: Yes
+
 BuildType: Package
 PackageUseDevtools: Yes
 PackageInstallArgs: --no-multiarch --with-keep.source
diff --git a/man/RunModel_GR1A.Rd b/man/RunModel_GR1A.Rd
index 1bb06147d14e4ecc26a93deaa8a157e14efe9d2d..8bb42dac18820fc1f52ab3223525f025081a66ba 100644
--- a/man/RunModel_GR1A.Rd
+++ b/man/RunModel_GR1A.Rd
@@ -24,7 +24,7 @@ RunModel_GR1A(InputsModel, RunOptions, Param)
 \item{RunOptions}{[object of class \emph{RunOptions}] see \code{\link{CreateRunOptions}} for details}
 
 \item{Param}{[numeric] vector of 1 parameter
-\tabular{ll}{                                                                      
+\tabular{ll}{
 GR1A X1 \tab model parameter [-] \cr
 }}
 }
@@ -56,23 +56,21 @@ library(airGR)
 data(L0123001)
 
 ## conversion of example data from daily to yearly time step
-TabSeries       <- data.frame(BasinObs$DatesR, BasinObs$P, BasinObs$E, BasinObs$T, BasinObs$Qmm)
-TimeFormat      <- "daily"
-NewTimeFormat   <- "yearly"
-ConvertFun      <- c("sum", "sum", "mean", "sum")
-YearFirstMonth  <- 09;
-NewTabSeries    <- SeriesAggreg(TabSeries = TabSeries, TimeFormat = TimeFormat, 
-                                NewTimeFormat = NewTimeFormat, ConvertFun = ConvertFun, 
-                                YearFirstMonth = YearFirstMonth)
-BasinObs        <- NewTabSeries
-names(BasinObs) <- c("DatesR", "P", "E", "T", "Qmm")
+TabSeries <- data.frame(DatesR = BasinObs$DatesR,
+                        P = BasinObs$P,
+                        E = BasinObs$E,
+                        Qmm = BasinObs$Qmm)
+TabSeries <- TabSeries[TabSeries$DatesR < "2012-09-01", ]
+BasinObs <- SeriesAggreg(TabSeries, Format = "\%Y",
+                         YearFirstMonth = 09,
+                         ConvertFun = c("sum", "sum", "sum"))
 
 ## preparation of the InputsModel object
-InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR1A, DatesR = BasinObs$DatesR, 
+InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR1A, DatesR = BasinObs$DatesR,
                                  Precip = BasinObs$P, PotEvap = BasinObs$E)
 
 ## run period selection
-Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%Y")=="1990"), 
+Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%Y")=="1990"),
                which(format(BasinObs$DatesR, format = "\%Y")=="1999"))
 
 ## preparation of the RunOptions object
@@ -87,7 +85,7 @@ OutputsModel <- RunModel_GR1A(InputsModel = InputsModel, RunOptions = RunOptions
 plot(OutputsModel, Qobs = BasinObs$Qmm[Ind_Run])
 
 ## efficiency criterion: Nash-Sutcliffe Efficiency
-InputsCrit  <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, 
+InputsCrit  <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
                                 RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])
 OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
 }
@@ -99,8 +97,8 @@ Laurent Coron, Claude Michel, Olivier Delaigue, Guillaume Thirel
 
 
 \references{
-Mouelhi S. (2003). 
-  Vers une chaîne cohérente de modèles pluie-débit conceptuels globaux aux pas de temps pluriannuel, annuel, mensuel et journalier. 
+Mouelhi S. (2003).
+  Vers une chaîne cohérente de modèles pluie-débit conceptuels globaux aux pas de temps pluriannuel, annuel, mensuel et journalier.
   PhD thesis (in French), ENGREF - Cemagref Antony, France.
 }
 
diff --git a/man/RunModel_GR2M.Rd b/man/RunModel_GR2M.Rd
index 208f39edca9487b3ce2c953910500bbedc1fc36f..ccfd7ee83684268d90f55295afdfde06fe4ee30b 100644
--- a/man/RunModel_GR2M.Rd
+++ b/man/RunModel_GR2M.Rd
@@ -24,7 +24,7 @@ RunModel_GR2M(InputsModel, RunOptions, Param)
 \item{RunOptions}{[object of class \emph{RunOptions}] see \code{\link{CreateRunOptions}} for details}
 
 \item{Param}{[numeric] vector of 2 parameters
-\tabular{ll}{                                                                      
+\tabular{ll}{
 GR2M X1 \tab production store capacity [mm]       \cr
 GR2M X2 \tab groundwater exchange coefficient [-] \cr
 }}
@@ -32,8 +32,8 @@ GR2M X2 \tab groundwater exchange coefficient [-] \cr
 
 
 \value{
-[list] list containing the function outputs organised as follows:                                         
-  \tabular{ll}{                                                                                         
+[list] list containing the function outputs organised as follows:
+  \tabular{ll}{
     \emph{$DatesR  } \tab [POSIXlt] series of dates                                            \cr
     \emph{$PotEvap } \tab [numeric] series of input potential evapotranspiration [mm/month]    \cr
     \emph{$Precip  } \tab [numeric] series of input total precipitation [mm/month]             \cr
@@ -47,7 +47,7 @@ GR2M X2 \tab groundwater exchange coefficient [-] \cr
     \emph{$Rout    } \tab [numeric] series of routing store level [mm]                         \cr
     \emph{$Qsim    } \tab [numeric] series of simulated discharge [mm/month]                   \cr
     \emph{$StateEnd} \tab [numeric] states at the end of the run (production store level and routing store level) [mm], \cr\tab see \code{\link{CreateIniStates}} for more details \cr
-  }                                                                                                     
+  }
   (refer to the provided references or to the package source code for further details on these model outputs)
 }
 
@@ -68,23 +68,16 @@ library(airGR)
 ## loading catchment data
 data(L0123001)
 
-## conversion of example data from daily to monthly time step
-TabSeries       <- data.frame(BasinObs$DatesR, BasinObs$P, BasinObs$E, BasinObs$T, BasinObs$Qmm)
-TimeFormat      <- "daily"
-NewTimeFormat   <- "monthly"
-ConvertFun      <- c("sum", "sum", "mean", "sum")
-NewTabSeries    <- SeriesAggreg(TabSeries = TabSeries, TimeFormat = TimeFormat, 
-                                NewTimeFormat = NewTimeFormat, ConvertFun = ConvertFun)
-BasinObs        <- NewTabSeries
-names(BasinObs) <- c("DatesR", "P", "E", "T", "Qmm")
-
-## preparation of the InputsModel object
-InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR2M, DatesR = BasinObs$DatesR, 
+## preparation of the InputsModel object with daily time step data
+InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
                                  Precip = BasinObs$P, PotEvap = BasinObs$E)
 
+## conversion of InputsModel to monthly time step
+InputsModel <- SeriesAggreg(InputsModel, Format = "\%Y\%m")
+
 ## run period selection
-Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%Y-\%m")=="1990-01"), 
-               which(format(BasinObs$DatesR, format = "\%Y-\%m")=="1999-12"))
+Ind_Run <- seq(which(format(InputsModel$DatesR, format = "\%Y-\%m")=="1990-01"),
+               which(format(InputsModel$DatesR, format = "\%Y-\%m")=="1999-12"))
 
 ## preparation of the RunOptions object
 RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR2M,
@@ -94,12 +87,18 @@ RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR2M,
 Param <- c(X1 = 265.072, X2 = 1.040)
 OutputsModel <- RunModel_GR2M(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param)
 
+## conversion of observed discharge to monthly time step
+Qobs <- SeriesAggreg(data.frame(BasinObs$DatesR, BasinObs$Qmm),
+                    Format = "\%Y\%m",
+                    ConvertFun = "sum")
+Qobs <- Qobs[Ind_Run, 2]
+
 ## results preview
-plot(OutputsModel, Qobs = BasinObs$Qmm[Ind_Run])
+plot(OutputsModel, Qobs = Qobs)
 
 ## efficiency criterion: Nash-Sutcliffe Efficiency
-InputsCrit  <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, 
-                                RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])
+InputsCrit  <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
+                                RunOptions = RunOptions, Obs = Qobs)
 OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
 }
 
@@ -110,12 +109,12 @@ Laurent Coron, Claude Michel, Safouane Mouelhi, Olivier Delaigue, Guillaume Thir
 
 
 \references{
-Mouelhi S. (2003). 
-  Vers une chaîne cohérente de modèles pluie-débit conceptuels globaux aux pas de temps pluriannuel, annuel, mensuel et journalier. 
+Mouelhi S. (2003).
+  Vers une chaîne cohérente de modèles pluie-débit conceptuels globaux aux pas de temps pluriannuel, annuel, mensuel et journalier.
   PhD thesis (in French), ENGREF - Cemagref Antony, France.
 \cr\cr
-Mouelhi, S., Michel, C., Perrin, C. and Andréassian V. (2006). 
-  Stepwise development of a two-parameter monthly water balance model. 
+Mouelhi, S., Michel, C., Perrin, C. and Andréassian V. (2006).
+  Stepwise development of a two-parameter monthly water balance model.
   Journal of Hydrology, 318(1-4), 200-214, doi: \href{https://www.doi.org/10.1016/j.jhydrol.2005.06.014}{10.1016/j.jhydrol.2005.06.014}.
 }
 
diff --git a/man/SeriesAggreg.Rd b/man/SeriesAggreg.Rd
index 4498eb6e09ef89614ef7561ff74e7e061e581363..2c095b545f4e10a5fb985e31a4e55c35fec73d69 100644
--- a/man/SeriesAggreg.Rd
+++ b/man/SeriesAggreg.Rd
@@ -3,41 +3,81 @@
 
 \name{SeriesAggreg}
 \alias{SeriesAggreg}
+\alias{SeriesAggreg.list}
+\alias{SeriesAggreg.data.frame}
+\alias{SeriesAggreg.InputsModel}
+\alias{SeriesAggreg.OutputsModel}
 
 
-\title{Conversion of time series to another time step (aggregation only)}
+\title{Conversion of time series to another time step (aggregation only) and regime computation}
 
 
 \description{
-Conversion of time series to another time step (aggregation only). \cr
-Warning : on the aggregated outputs, the dates correspond to the beginning of the time step \cr
-(e.g. for daily time-series 2005-03-01 00:00 = value for period 2005-03-01 00:00 - 2005-03-01 23:59) \cr
-(e.g. for monthly time-series 2005-03-01 00:00 = value for period 2005-03-01 00:00 - 2005-03-31 23:59) \cr
-(e.g. for yearly time-series 2005-03-01 00:00 = value for period 2005-03-01 00:00 - 2006-02-28 23:59)
+Conversion of time series to another time step (aggregation only) and regime computation. \cr
+Warning: on the aggregated outputs, the dates correspond to the beginning of the time step \cr
+(e.g. for daily time series 2005-03-01 00:00 = value for period 2005-03-01 00:00 - 2005-03-01 23:59) \cr
+(e.g. for monthly time series 2005-03-01 00:00 = value for period 2005-03-01 00:00 - 2005-03-31 23:59) \cr
+(e.g. for yearly time series 2005-03-01 00:00 = value for period 2005-03-01 00:00 - 2006-02-28 23:59)
+}
+
+
+\details{
+  \code{\link{SeriesAggreg.InputsModel}} and \code{\link{SeriesAggreg.OutputsModel}}
+  call \code{\link{SeriesAggreg.list}} which itself calls \code{\link{SeriesAggreg.data.frame}}.
+  So, all arguments passed to any \code{\link{SeriesAggreg}} method will be passed to \code{\link{SeriesAggreg.data.frame}}.
 }
 
 
 \usage{
-SeriesAggreg(TabSeries, TimeFormat, NewTimeFormat, ConvertFun,
-             YearFirstMonth = 1, TimeLag = 0, verbose = TRUE)
+\method{SeriesAggreg}{data.frame}(x,
+             Format,
+             ConvertFun,
+             TimeFormat = NULL,
+             NewTimeFormat = NULL,
+             YearFirstMonth = 1,
+             TimeLag = 0,
+             \dots)
+
+\method{SeriesAggreg}{list}(x,
+             Format,
+             ConvertFun,
+             NewTimeFormat = NULL,
+             simplify = FALSE,
+             except = NULL,
+             recursive = TRUE,
+             \dots)
+
+\method{SeriesAggreg}{InputsModel}(x, \dots)
+
+\method{SeriesAggreg}{OutputsModel}(x, \dots)
 }
 
 
 \arguments{
-\item{TabSeries}{[POSIXt+numeric] data.frame containing the vector of dates (POSIXt) and the time series values numeric)}
+\item{x}{[InputsModel], [OutputsModel], [list] or [data.frame] containing the vector of dates (\emph{POSIXt}) and the time series of numeric values}
+
+\item{Format}{[character] output time step format (i.e. yearly times series: \code{"\%Y"}, monthly time series: \code{"\%Y\%m"}, daily time series: \code{"\%Y\%m\%d"}, monthly regimes: \code{"\%m"}, daily regimes: \code{"\%d"})}
+
+\item{TimeFormat}{(deprecated) [character] input time step format (i.e. \code{"hourly"}, \code{"daily"}, \code{"monthly"} or \code{"yearly"}). Use the \code{x} argument instead}
+
+\item{NewTimeFormat}{(deprecated) [character] output time step format (i.e. \code{"hourly"}, \code{"daily"}, \code{"monthly"} or \code{"yearly"}). Use the \code{Format} argument instead}
 
-\item{TimeFormat}{[character] input time-step format (i.e. \code{"hourly"}, \code{"daily"}, \code{"monthly"} or \code{"yearly"})}
+\item{ConvertFun}{[character] names of aggregation functions (e.g. for P[mm], T[degC], Q[mm]: \code{ConvertFun = c("sum", "mean", "sum"})) (default: use the name of the column (see details) or \code{"mean"} for regime calculation)}
 
-\item{NewTimeFormat}{[character] output time-step format (i.e. \code{"hourly"}, \code{"daily"}, \code{"monthly"} or \code{"yearly"})}
+\item{YearFirstMonth}{(optional) [numeric] integer used when \code{Format = "\%Y"} to set when the starting month of the year (e.g. 01 for calendar year or 09 for hydrological year starting in September)}
 
-\item{ConvertFun}{[character] names of aggregation functions (e.g. for P[mm], T[degC], Q[mm] : \code{ConvertFun = c("sum", "mean", "sum"}))}
+\item{TimeLag}{(optional) [numeric] numeric indicating a time lag (in seconds) for the time series aggregation (especially useful to aggregate hourly time series into daily time series)}
 
-\item{YearFirstMonth}{(optional) [numeric] integer used when \code{NewTimeFormat = "yearly"} to set when the starting month of the year (e.g. 01 for calendar year or 09 for hydrological year starting in September)}
+\item{simplify}{(optional) [boolean] if set to \code{TRUE}, a \code{\link{data.frame}} is returned instead of a \code{\link{list}}. Embedded lists are then ignored. (default = \code{FALSE})}
 
-\item{TimeLag}{(optional) [numeric] numeric indicating a time lag (in seconds) for the time series aggregation (especially useful to aggregate hourly time series in daily time series)}
+\item{except}{(optional) [character] the name of the items to skip in the aggregation (default = \code{NULL})}
 
-\item{verbose}{(optional) [boolean] boolean indicating if the function is run in verbose mode or not, default = \code{FALSE}}
+\item{recursive}{(optional) [boolean] if set to \code{FALSE}, embedded lists and dataframes are not aggregated (default = \code{TRUE})}
+
+\item{\dots}{Arguments passed to \code{\link{SeriesAggreg.list}} and then to \code{\link{SeriesAggreg.data.frame}}}
 }
+
+
 \value{
 [POSIXct+numeric] data.frame containing a vector of aggregated dates (POSIXct) and time series values numeric)
 }
@@ -52,20 +92,29 @@ data(L0123002)
 ## preparation of the initial time series data frame at the daily time step
 TabSeries <- BasinObs[, c("DatesR", "P", "E", "T", "Qmm")]
 
-## conversion at the monthly time step
-NewTabSeries <- SeriesAggreg(TabSeries = TabSeries,
-                             TimeFormat = "daily", NewTimeFormat = "monthly",
+## monthly time series
+NewTabSeries <- SeriesAggreg(TabSeries,
+                             Format = "\%Y\%m",
                              ConvertFun = c("sum", "sum", "mean", "sum"))
+str(NewTabSeries)
 
-## conversion at the yearly time step
-NewTabSeries <- SeriesAggreg(TabSeries = TabSeries,
-                             TimeFormat = "daily", NewTimeFormat = "yearly",
+## monthly regimes
+NewTabSeries <- SeriesAggreg(TabSeries,
+                             Format = "\%m",
                              ConvertFun = c("sum", "sum", "mean", "sum"))
+str(NewTabSeries)
+
+## conversion of InputsModel
+example("RunModel_GR2M")
+
+## monthly regimes on OutputsModel object
+SimulatedMonthlyRegime <- SeriesAggreg(OutputsModel, Format = "\%m")
+str(SimulatedMonthlyRegime)
 
 }
 
 
 \author{
-Laurent Coron
+Olivier Delaigue, David Dorchies
 }
 
diff --git a/tests/testthat/test-SeriesAggreg.R b/tests/testthat/test-SeriesAggreg.R
new file mode 100644
index 0000000000000000000000000000000000000000..f71db81857c4531093ed3c6d123dcbea053dfe45
--- /dev/null
+++ b/tests/testthat/test-SeriesAggreg.R
@@ -0,0 +1,209 @@
+context("SeriesAggreg")
+
+## load catchment data
+data(L0123002)
+
+test_that("No warning with InputsModel Cemaneige'", {
+  ## preparation of the InputsModel object
+  InputsModel <- CreateInputsModel(
+    FUN_MOD = RunModel_CemaNeige,
+    DatesR = BasinObs$DatesR,
+    Precip = BasinObs$P,
+    TempMean = BasinObs$T,
+    ZInputs = BasinInfo$HypsoData[51],
+    HypsoData = BasinInfo$HypsoData,
+    NLayers = 5
+  )
+  # Expect no warning: https://stackoverflow.com/a/33638939/5300212
+  expect_warning(SeriesAggreg(InputsModel, "%m"),
+                 regexp = NA)
+})
+
+test_that("Warning: deprecated 'TimeFormat' argument", {
+  InputsModel <- CreateInputsModel(
+    FUN_MOD = RunModel_GR4J,
+    DatesR = BasinObs$DatesR,
+    Precip = BasinObs$P,
+    PotEvap = BasinObs$E
+  )
+  expect_warning(SeriesAggreg(InputsModel, Format = "%Y%m", TimeFormat = "daily"),
+                 regexp = "deprecated 'TimeFormat' argument")
+})
+
+test_that("Warning: deprecated 'NewTimeFormat' argument: please use 'Format' instead",
+          {
+            InputsModel <- CreateInputsModel(
+              FUN_MOD = RunModel_GR4J,
+              DatesR = BasinObs$DatesR,
+              Precip = BasinObs$P,
+              PotEvap = BasinObs$E
+            )
+            expect_warning(SeriesAggreg(InputsModel, NewTimeFormat = "monthly"),
+                           regexp = "deprecated 'NewTimeFormat' argument: please use 'Format' instead")
+          })
+
+test_that("Warning: deprecated 'NewTimeFormat' argument: 'Format' argument is used instead",
+          {
+            InputsModel <- CreateInputsModel(
+              FUN_MOD = RunModel_GR4J,
+              DatesR = BasinObs$DatesR,
+              Precip = BasinObs$P,
+              PotEvap = BasinObs$E
+            )
+            expect_warning(SeriesAggreg(InputsModel, Format = "%Y%m", NewTimeFormat = "monthly"),
+                           regexp = "deprecated 'NewTimeFormat' argument: 'Format' argument is used instead")
+          })
+
+test_that("Check SeriesAggreg output values on yearly aggregation", {
+  TabSeries <- data.frame(
+    DatesR = BasinObs$DatesR,
+    P = BasinObs$P,
+    E = BasinObs$E,
+    Qmm = BasinObs$Qmm
+  )
+  GoodValues <- apply(BasinObs[BasinObs$DatesR >= "1984-09-01" &
+                                 BasinObs$DatesR < "1985-09-01",
+                               c("P", "E", "Qmm")], 2, sum)
+  TestedValues <- unlist(SeriesAggreg(TabSeries,
+                                      Format = "%Y",
+                                      YearFirstMonth = 9,
+                                      ConvertFun = rep("sum", 3))[1, c("P", "E", "Qmm")])
+  expect_equal(GoodValues, TestedValues)
+})
+
+test_that("No DatesR should warning", {
+  TabSeries <- list(
+    Dates = BasinObs$DatesR,
+    P = BasinObs$P,
+    E = BasinObs$E,
+    Qmm = BasinObs$Qmm
+  )
+  expect_warning(SeriesAggreg(TabSeries, "%Y%m"), regexp = "has been automatically chosen")
+})
+
+test_that("Check SeriesAggreg.list 'DatesR' argument", {
+  InputsModel <- CreateInputsModel(
+    FUN_MOD = RunModel_GR4J,
+    DatesR = BasinObs$DatesR,
+    Precip = BasinObs$P,
+    PotEvap = BasinObs$E
+  )
+  DatesR <- InputsModel$DatesR
+  # No InputsModel$DatesR
+  InputsModel$DatesR <- NULL
+  expect_error(SeriesAggreg(InputsModel, "%Y%m"), regexp = "'POSIXt'")
+  # Other list item chosen
+  InputsModel$SuperDates <- DatesR
+  expect_warning(SeriesAggreg(InputsModel, "%Y%m"), regexp = "SuperDates")
+  # Wrong InputsModel$DatesR
+  InputsModel$DatesR <- BasinObs$P
+  expect_error(SeriesAggreg(InputsModel, "%Y%m"), regexp = "'POSIXt'")
+
+})
+
+test_that("Check SeriesAggreg.list with embedded lists", {
+  InputsModel <-
+    CreateInputsModel(
+      FUN_MOD = RunModel_CemaNeige,
+      DatesR = BasinObs$DatesR,
+      Precip = BasinObs$P,
+      TempMean = BasinObs$T,
+      ZInputs = BasinInfo$HypsoData[51],
+      HypsoData = BasinInfo$HypsoData,
+      NLayers = 5
+    )
+  I2 <- SeriesAggreg(InputsModel, "%Y%m")
+  expect_equal(length(I2$ZLayers), 5)
+  expect_null(I2$LayerPrecip$DatesR)
+  expect_equal(length(I2$DatesR), length(I2$LayerPrecip$L1))
+})
+
+test_that("Check SeriesAggreg.outputsModel", {
+  InputsModel <-
+    CreateInputsModel(
+      FUN_MOD = RunModel_CemaNeigeGR4J,
+      DatesR = BasinObs$DatesR,
+      Precip = BasinObs$P,
+      PotEvap = BasinObs$E,
+      TempMean = BasinObs$T,
+      ZInputs = median(BasinInfo$HypsoData),
+      HypsoData = BasinInfo$HypsoData,
+      NLayers = 5
+    )
+
+  ## run period selection
+  Ind_Run <-
+    seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1990-01-01"),
+        which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1999-12-31"))
+
+  ## preparation of the RunOptions object
+  suppressWarnings(
+    RunOptions <-
+      CreateRunOptions(
+        FUN_MOD = RunModel_CemaNeigeGR4J,
+        InputsModel = InputsModel,
+        IndPeriod_Run = Ind_Run
+      )
+  )
+
+  ## simulation
+  Param <- c(
+    X1 = 408.774,
+    X2 = 2.646,
+    X3 = 131.264,
+    X4 = 1.174,
+    CNX1 = 0.962,
+    CNX2 = 2.249
+  )
+  OutputsModel <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
+                                         RunOptions = RunOptions,
+                                         Param = Param)
+
+  O2 <- SeriesAggreg(OutputsModel, "%Y%m")
+  expect_equal(length(O2$StateEnd), 3)
+  expect_equal(length(O2$DatesR),
+               length(O2$CemaNeigeLayers$Layer01$Pliq))
+})
+
+test_that("Check data.frame handling in SeriesAggreg.list", {
+  QObsUp <- imputeTS::na_interpolation(BasinObs$Qmm)
+  InputsModelDown1 <- CreateInputsModel(
+    FUN_MOD = RunModel_GR4J,
+    DatesR = BasinObs$DatesR,
+    Precip = BasinObs$P,
+    PotEvap = BasinObs$E,
+    Qupstream = matrix(QObsUp, ncol = 1),
+    # Upstream observed flow
+    LengthHydro = 100 * 1000,
+    # Distance between upstream catchment outlet and the downstream one in m
+    BasinAreas = c(180, 180) # Upstream and downstream areas in km²
+  )
+  expect_warning(SeriesAggreg(InputsModelDown1, "%Y%m"),
+                 regexp = NA)
+  I2 <- SeriesAggreg(InputsModelDown1, "%Y%m")
+  expect_equal(length(I2$DatesR), nrow(I2$Qupstream))
+  InputsModelDown1$Qupstream <-
+    InputsModelDown1$Qupstream[-1, , drop = FALSE] # https://stackoverflow.com/a/7352287/5300212
+  expect_warning(SeriesAggreg(InputsModelDown1, "%Y%m"),
+                 regexp = "it will be ignored in the aggregation")
+})
+test_that("SeriesAggreg from and to the same time step should return initial time series", {
+  InputsModel <- CreateInputsModel(
+    FUN_MOD = RunModel_GR4J,
+    DatesR = BasinObs$DatesR,
+    Precip = BasinObs$P,
+    PotEvap = BasinObs$E
+  )
+  I2 <- SeriesAggreg(InputsModel, "%Y%m")
+  expect_warning(SeriesAggreg(I2, "%Y%m"), regexp = "No time-step conversion was performed")
+  expect_equal(I2, suppressWarnings(SeriesAggreg(I2, "%Y%m")))
+})
+test_that("SeriesAggreg.data.frame with first column not named DatesR should work",
+          {
+            expect_warning(SeriesAggreg(
+              data.frame(BasinObs$DatesR, BasinObs$Qmm),
+              Format = "%Y%m",
+              ConvertFun = "sum"
+            ),
+            regexp = NA)
+          })