diff --git a/.Rbuildignore b/.Rbuildignore
index acd9e68a9a2e50eff8ef783a77b63bab66e63424..94af330f02015c9a87d80c74d9a1fe07a91fb191 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -9,3 +9,5 @@
 ^\.vscode$
 ^Rplots\.pdf$
 ^ci$
+^data-raw$
+^revdep$
diff --git a/.gitignore b/.gitignore
index 5eb9ae2487a29398fd86791b85cf36ef46b8763e..b9e05842088f86b33e293868b0b28f461db25096 100644
--- a/.gitignore
+++ b/.gitignore
@@ -12,6 +12,8 @@ packrat/lib*/
 *.pdf
 !man/figures/*.pdf
 
+# revdep
+/revdep/
 
 ######################################################################################################
 ### Generic .gitignore for R (source: https://github.com/github/gitignore/blob/master/R.gitignore) ###
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index c68113ae34b581318d04a12b9f9196d997e8b9fd..938d200bee886d9916de746355d2fb2788b32956 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,6 +1,6 @@
 stages:
   - check
-  - regression
+  - scheduled_tests
   - revdepcheck
 
 default:
@@ -10,13 +10,14 @@ default:
     - PATH=~/R/sources/R-${R_VERSION}/bin:$PATH
     - R -e 'remotes::install_deps(dep = TRUE)'
 
-.regression:
-  stage: regression
+.scheduled_tests:
+  stage: scheduled_tests
   script:
-    - Rscript tests/testthat/regression_tests.R stable
+    - Rscript tests/scheduled_tests/scheduled.R
+    - Rscript tests/scheduled_tests/regression.R stable
     - R CMD INSTALL .
-    - Rscript tests/testthat/regression_tests.R dev
-    - Rscript tests/testthat/regression_tests.R compare
+    - Rscript tests/scheduled_tests/regression.R dev
+    - Rscript tests/scheduled_tests/regression.R compare
 
 .check:
   stage: check
@@ -33,26 +34,31 @@ default:
     NOT_CRAN: "false"
   extends: .check
 
-regression_patched:
+scheduled_tests_patched:
+  only:
+    refs:
+      - dev
+      - master
+      - schedules
   variables:
     R_VERSION: "patched"
-  extends: .regression
+  extends: .scheduled_tests
 
-regression_devel:
+scheduled_tests_devel:
   only:
     refs:
       - schedules
   variables:
     R_VERSION: "devel"
-  extends: .regression
+  extends: .scheduled_tests
 
-regression_oldrel:
+scheduled_tests_oldrel:
   only:
     refs:
       - schedules
   variables:
     R_VERSION: "oldrel"
-  extends: .regression
+  extends: .scheduled_tests
 
 check_not_cran_patched:
   variables:
diff --git a/.regressionignore b/.regressionignore
index 5ceae7cdb93d9480568fbd4545a4d8c7968177a1..f7eeb15b1d582fbb8eb0c291683b03a4467b6029 100644
--- a/.regressionignore
+++ b/.regressionignore
@@ -4,3 +4,62 @@
 # ignored variable : [Topic]<SPACE>[Variable].
 # Example for ignoring OutputsModel variable produced by example("RunModel_GR2M"): RunModel_GR2M OutputsModel
 
+Calibration_Michel RunOptions
+Calibration RunOptions
+CreateCalibOptions RunOptions
+CreateIniStates RunOptions
+CreateInputsCrit RunOptions
+CreateInputsModel RunOptions
+CreateRunOptions RunOptions
+ErrorCrit_KGE RunOptions
+ErrorCrit_KGE2 RunOptions
+ErrorCrit_NSE RunOptions
+ErrorCrit_RMSE RunOptions
+ErrorCrit RunOptions
+Imax RunOptions
+Param_Sets_GR4J RunOptions_Cal
+Param_Sets_GR4J RunOptions_Val
+RunModel_CemaNeige RunOptions
+RunModel_CemaNeigeGR4J RunOptions
+RunModel_CemaNeigeGR5J RunOptions
+RunModel_CemaNeigeGR6J RunOptions
+RunModel_GR1A RunOptions
+RunModel_GR2M RunOptions
+RunModel_GR4H RunOptions
+RunModel_GR4J RunOptions
+RunModel_GR5H RunOptions
+RunModel_GR5J RunOptions
+RunModel_GR6J RunOptions
+RunModel_Lag RunOptions
+RunModel RunOptions
+SeriesAggreg RunOptions
+Calibration OutputsModel
+Calibration_Michel OutputsModel
+CreateCalibOptions OutputsModel
+CreateIniStates OutputsModel
+CreateInputsCrit OutputsModel
+CreateInputsModel OutputsModel
+CreateRunOptions OutputsModel
+ErrorCrit OutputsModel
+ErrorCrit_KGE OutputsModel
+ErrorCrit_KGE2 OutputsModel
+ErrorCrit_NSE OutputsModel
+ErrorCrit_RMSE OutputsModel
+Imax OutputsModel
+RunModel OutputsModel
+RunModel_CemaNeige OutputsModel
+RunModel_CemaNeigeGR4J OutputsModel
+RunModel_CemaNeigeGR5J OutputsModel
+RunModel_CemaNeigeGR6J OutputsModel
+RunModel_GR1A OutputsModel
+RunModel_GR2M OutputsModel
+RunModel_GR4H OutputsModel
+RunModel_GR4J OutputsModel
+RunModel_GR5H OutputsModel
+RunModel_GR5J OutputsModel
+RunModel_GR6J OutputsModel
+RunModel_Lag OutputsModel
+SeriesAggreg OutputsModel
+Param_Sets_GR4J OutputsModel_Val
+RunModel_Lag OutputsModelDown
+SeriesAggreg SimulatedMonthlyRegime
diff --git a/R/Calibration_Michel.R b/R/Calibration_Michel.R
index 0dd1eec50da9aca74f4c9353c0b8e6c878d9207d..1ad475be2811e1dd7e47042362adabe31044c103 100644
--- a/R/Calibration_Michel.R
+++ b/R/Calibration_Michel.R
@@ -17,7 +17,7 @@ Calibration_Michel <- function(InputsModel,
   # Handling 'FUN_TRANSFO' from direct argument or provided by 'CaliOptions'
   if (!is.null(FUN_TRANSFO)) {
     FUN_TRANSFO <- match.fun(FUN_TRANSFO)
-  } else if(!is.null(CalibOptions$FUN_TRANSFO)) {
+  } else if (!is.null(CalibOptions$FUN_TRANSFO)) {
     FUN_TRANSFO <- CalibOptions$FUN_TRANSFO
   } else {
     stop("'FUN_TRANSFO' is not provided neither as 'FUN_TRANSFO' argument or in 'CaliOptions' argument")
diff --git a/R/CreateCalibOptions.R b/R/CreateCalibOptions.R
index 1d13f2ae29cf394512b220986941176b575b2bf5..f6d5ffb60c6ebb770d9fb0a3b0439d9df2885ad4 100644
--- a/R/CreateCalibOptions.R
+++ b/R/CreateCalibOptions.R
@@ -12,7 +12,7 @@ CreateCalibOptions <- function(FUN_MOD,
 
   FUN_MOD     <- match.fun(FUN_MOD)
   FUN_CALIB   <- match.fun(FUN_CALIB)
-  if(!is.null(FUN_TRANSFO)) {
+  if (!is.null(FUN_TRANSFO)) {
     FUN_TRANSFO <- match.fun(FUN_TRANSFO)
   }
   if (!is.logical(IsHyst) | length(IsHyst) != 1L) {
diff --git a/R/CreateIniStates.R b/R/CreateIniStates.R
index 6b311d414f2e3f5fd486056b5ace76e0b605c724..9fe48b137b3a9abbda09836d93f5907f1ec1fd39 100644
--- a/R/CreateIniStates.R
+++ b/R/CreateIniStates.R
@@ -153,19 +153,19 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
     }
     UH2 <- rep(Inf, UH2n)
   }
-  if(IsIntStore & is.null(IntStore)) {
+  if (IsIntStore & is.null(IntStore)) {
       stop(sprintf("'%s' need values for 'IntStore'", nameFUN_MOD))
   }
-  if("CemaNeige" %in% ObjectClass & !IsHyst &
+  if ("CemaNeige" %in% ObjectClass & !IsHyst &
      (is.null(GCemaNeigeLayers) | is.null(eTGCemaNeigeLayers))) {
       stop(sprintf("'%s' need values for 'GCemaNeigeLayers' and 'GCemaNeigeLayers'", nameFUN_MOD))
   }
-  if("CemaNeige" %in% ObjectClass & IsHyst &
+  if ("CemaNeige" %in% ObjectClass & IsHyst &
      (is.null(GCemaNeigeLayers) | is.null(eTGCemaNeigeLayers) |
       is.null(GthrCemaNeigeLayers) | is.null(GlocmaxCemaNeigeLayers))) {
     stop(sprintf("'%s' need values for 'GCemaNeigeLayers', 'GCemaNeigeLayers', 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'", nameFUN_MOD))
   }
-  if("CemaNeige" %in% ObjectClass & !IsHyst &
+  if ("CemaNeige" %in% ObjectClass & !IsHyst &
      (!is.null(GthrCemaNeigeLayers) | !is.null(GlocmaxCemaNeigeLayers))) {
     if (verbose) {
       warning(sprintf("'%s' does not require 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'. Values set to NA", nameFUN_MOD))
@@ -173,7 +173,7 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
     GthrCemaNeigeLayers    <- Inf
     GlocmaxCemaNeigeLayers <- Inf
   }
-  if(!"CemaNeige" %in% ObjectClass &
+  if (!"CemaNeige" %in% ObjectClass &
      (!is.null(GCemaNeigeLayers) | !is.null(eTGCemaNeigeLayers) | !is.null(GthrCemaNeigeLayers) | !is.null(GlocmaxCemaNeigeLayers))) {
     if (verbose) {
       warning(sprintf("'%s' does not require 'GCemaNeigeLayers' 'GCemaNeigeLayers', 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'. Values set to NA", nameFUN_MOD))
@@ -186,7 +186,7 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
 
 
   ## set states
-  if("CemaNeige" %in% ObjectClass) {
+  if ("CemaNeige" %in% ObjectClass) {
     NLayers <- length(InputsModel$LayerPrecip)
   } else {
     NLayers <- 1
@@ -284,17 +284,17 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
   }
 
   # SD model state handling
-  if(!is.null(SD)) {
-    if(!inherits(InputsModel, "SD")) {
+  if (!is.null(SD)) {
+    if (!inherits(InputsModel, "SD")) {
       stop("'SD' argument provided and 'InputsModel' is not of class 'SD'")
     }
-    if(!is.list(SD)) {
+    if (!is.list(SD)) {
       stop("'SD' argument must be a list")
     }
     lapply(SD, function(x) {
-      if(!is.numeric(x)) stop("Each item of 'SD' list argument must be numeric")
+      if (!is.numeric(x)) stop("Each item of 'SD' list argument must be numeric")
     })
-    if(length(SD) != length(InputsModel$LengthHydro)) {
+    if (length(SD) != length(InputsModel$LengthHydro)) {
       stop("Number of items of 'SD' list argument must be the same as the number of upstream connections",
            sprintf(" (%i required, found %i)", length(InputsModel$LengthHydro), length(SD)))
     }
@@ -309,15 +309,15 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
   IniStatesNA[is.infinite(IniStatesNA)] <- NA
   IniStatesNA <- relist(IniStatesNA, skeleton = IniStates)
 
-  if(!is.null(SD)) {
+  if (!is.null(SD)) {
     IniStatesNA$SD <- SD
   }
 
   class(IniStatesNA) <- c("IniStates", ObjectClass)
-  if(IsHyst) {
+  if (IsHyst) {
     class(IniStatesNA) <- c(class(IniStatesNA), "hysteresis")
   }
-  if(IsIntStore) {
+  if (IsIntStore) {
     class(IniStatesNA) <- c(class(IniStatesNA), "interception")
   }
 
diff --git a/R/CreateInputsCrit.R b/R/CreateInputsCrit.R
index 56e1732b03f3781842908774c2c21a2d7ab25ad2..4107b328b4f37834caacd541adbc1adf2829f5b4 100644
--- a/R/CreateInputsCrit.R
+++ b/R/CreateInputsCrit.R
@@ -293,7 +293,7 @@ CreateInputsCrit <- function(FUN_CRIT,
   listVarObs <- sapply(InputsCrit, FUN = "[[", "VarObs")
   inCnVarObs <- c("SCA", "SWE")
   if (!"ZLayers" %in% names(InputsModel)) {
-    if(any(listVarObs %in% inCnVarObs)) {
+    if (any(listVarObs %in% inCnVarObs)) {
       stop(sprintf("'VarObs' can not be equal to %i if CemaNeige is not used",
                    paste(sapply(inCnVarObs, shQuote), collapse = " or ")))
     }
@@ -348,7 +348,7 @@ CreateInputsCrit <- function(FUN_CRIT,
     combInputsCrit <- combn(x = length(InputsCrit), m = 2)
     apply(combInputsCrit, MARGIN = 2, function(i) {
       equalInputsCrit <- identical(InputsCrit[[i[1]]], InputsCrit[[i[2]]])
-      if(equalInputsCrit) {
+      if (equalInputsCrit) {
         warning(sprintf("elements %i and %i of the criteria list are identical. This might not be necessary", i[1], i[2]), call. = FALSE)
       }
     })
diff --git a/R/CreateInputsModel.R b/R/CreateInputsModel.R
index 958d93ee5d39a125f49f94e81a23abcb4da64606..cd6d26e9cfa469abf661adb7bcfbd1237bd85cbd 100644
--- a/R/CreateInputsModel.R
+++ b/R/CreateInputsModel.R
@@ -153,10 +153,10 @@ CreateInputsModel <- function(FUN_MOD,
     if (nrow(Qupstream) != LLL) {
       stop("'Qupstream' must have same number of rows as 'DatesR' length")
     }
-    if(any(is.na(Qupstream))) {
+    if (any(is.na(Qupstream))) {
       warning("'Qupstream' contains NA values: model outputs will contain NAs")
     }
-    if(any(LengthHydro > 1000)) {
+    if (any(LengthHydro > 1000)) {
       warning("The unit of 'LengthHydro' has changed from m to km in airGR >= 1.6.12: values superior to 1000 km seem unrealistic")
     }
     QupstrUnit <- tolower(QupstrUnit)
diff --git a/R/CreateRunOptions.R b/R/CreateRunOptions.R
index 4edbdd7184f0b489ebde8c6628bae9eb5a6ab998..8821935b9b971faf543f3a5f47735376a5a1edb2 100644
--- a/R/CreateRunOptions.R
+++ b/R/CreateRunOptions.R
@@ -24,12 +24,17 @@ CreateRunOptions <- function(FUN_MOD, InputsModel,
   ObjectClass <- FeatFUN_MOD$Class
   TimeStepMean <- FeatFUN_MOD$TimeStepMean
 
+  ## Model output variable list
+  FortranOutputs <- .FortranOutputs(GR = FeatFUN_MOD$CodeModHydro,
+                                    isCN = "CemaNeige" %in% FeatFUN_MOD$Class)
+
   ## manage class
   if (IsIntStore) {
     ObjectClass <- c(ObjectClass, "interception")
   }
   if (IsHyst) {
     ObjectClass <- c(ObjectClass, "hysteresis")
+    FeatFUN_MOD$NbParam <- FeatFUN_MOD$NbParam + 2
   }
 
   if (!"CemaNeige" %in% ObjectClass & "hysteresis" %in% ObjectClass) {
@@ -290,31 +295,7 @@ CreateRunOptions <- function(FUN_MOD, InputsModel,
   ##check_Outputs_Cal_and_Sim
 
   ##Outputs_all
-  Outputs_all <- NULL
-  if (identical(FUN_MOD, RunModel_GR4H) | identical(FUN_MOD, RunModel_CemaNeigeGR4H)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR4H")$GR)
-  }
-  if (identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR5H")$GR)
-  }
-  if (identical(FUN_MOD, RunModel_GR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR4J)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR4J")$GR)
-  }
-  if (identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR5J")$GR)
-  }
-  if (identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR6J")$GR)
-  }
-  if (identical(FUN_MOD, RunModel_GR2M)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR2M")$GR)
-  }
-  if (identical(FUN_MOD, RunModel_GR1A)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR1A")$GR)
-  }
-  if ("CemaNeige" %in% ObjectClass) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = NULL, isCN = TRUE)$CN)
-  }
+  Outputs_all <- c("DatesR", unlist(FortranOutputs), "WarmUpQsim", "StateEnd")
 
   ##check_Outputs_Sim
   if (!is.vector(Outputs_Sim)) {
@@ -327,9 +308,9 @@ CreateRunOptions <- function(FUN_MOD, InputsModel,
     stop("'Outputs_Sim' must not contain NA")
   }
   if ("all" %in% Outputs_Sim) {
-    Outputs_Sim <- c("DatesR", Outputs_all, "StateEnd")
+    Outputs_Sim <- Outputs_all
   }
-  Test <- which(!Outputs_Sim %in% c("DatesR", Outputs_all, "StateEnd"))
+  Test <- which(!Outputs_Sim %in% Outputs_all)
   if (length(Test) != 0) {
     stop(paste0( "'Outputs_Sim' is incorrectly defined: ",
                  paste(Outputs_Sim[Test], collapse = ", "), " not found"))
@@ -361,10 +342,9 @@ CreateRunOptions <- function(FUN_MOD, InputsModel,
     }
   }
   if ("all" %in% Outputs_Cal) {
-    Outputs_Cal <- c("DatesR", Outputs_all, "StateEnd")
-
+    Outputs_Cal <- Outputs_all
   }
-  Test <- which(!Outputs_Cal %in% c("DatesR", Outputs_all, "StateEnd"))
+  Test <- which(!Outputs_Cal %in% Outputs_all)
   if (length(Test) != 0) {
     stop(paste0("'Outputs_Cal' is incorrectly defined: ",
                 paste(Outputs_Cal[Test], collapse = ", "), " not found"))
@@ -473,7 +453,9 @@ CreateRunOptions <- function(FUN_MOD, InputsModel,
                      IniStates = IniStates,
                      IniResLevels = IniResLevels,
                      Outputs_Cal = Outputs_Cal,
-                     Outputs_Sim = Outputs_Sim)
+                     Outputs_Sim = Outputs_Sim,
+                     FortranOutputs = FortranOutputs,
+                     FeatFUN_MOD = FeatFUN_MOD)
 
   if ("CemaNeige" %in% ObjectClass) {
     RunOptions <- c(RunOptions, list(MeanAnSolidPrecip = MeanAnSolidPrecip))
diff --git a/R/RunModel_CemaNeige.R b/R/RunModel_CemaNeige.R
index ce2e4a6915b4a0b6f722294bc68b410cccfafeeb..a2aad9cec56aec3133671ea04b56a1b71c03d1d9 100644
--- a/R/RunModel_CemaNeige.R
+++ b/R/RunModel_CemaNeige.R
@@ -148,7 +148,7 @@ RunModel_CemaNeige <- function(InputsModel, RunOptions, Param) {
   
   ## End
   class(OutputsModel) <- c("OutputsModel", time_step, "CemaNeige")
-  if(IsHyst) {
+  if (IsHyst) {
     class(OutputsModel) <- c(class(OutputsModel), "hysteresis")
   }
   return(OutputsModel)
diff --git a/R/RunModel_CemaNeigeGR4H.R b/R/RunModel_CemaNeigeGR4H.R
index bda65ddbbbfffc5dde71ac8a17953c4494241657..51318c1faa81aea3bfad90c75b2eb1771b65a6b3 100644
--- a/R/RunModel_CemaNeigeGR4H.R
+++ b/R/RunModel_CemaNeigeGR4H.R
@@ -3,40 +3,12 @@ RunModel_CemaNeigeGR4H <- function(InputsModel, RunOptions, Param) {
 
   ## Initialization of variables
   IsHyst <- inherits(RunOptions, "hysteresis")
-  NParam <- ifelse(test = IsHyst, yes = 8L, no = 6L)
-  NParamCN <- NParam - 4L
+  NParamCN <- RunOptions$FeatFUN_MOD$NbParam - 4L
   NStates <- 4L
-  FortranOutputs <- .FortranOutputs(GR = "GR4H", isCN = TRUE)
 
 
-  ## Arguments check
-  if (!inherits(InputsModel, "InputsModel")) {
-    stop("'InputsModel' must be of class 'InputsModel'")
-  }
-  if (!inherits(InputsModel, "hourly")) {
-    stop("'InputsModel' must be of class 'hourly'")
-  }
-  if (!inherits(InputsModel, "GR")) {
-    stop("'InputsModel' must be of class 'GR'")
-  }
-  if (!inherits(InputsModel, "CemaNeige")) {
-    stop("'InputsModel' must be of class 'CemaNeige'")
-  }
-  if (!inherits(RunOptions, "RunOptions")) {
-    stop("'RunOptions' must be of class 'RunOptions'")
-  }
-  if (!inherits(RunOptions, "GR")) {
-    stop("'RunOptions' must be of class 'GR'")
-  }
-  if (!inherits(RunOptions, "CemaNeige")) {
-    stop("'RunOptions' must be of class 'CemaNeige'")
-  }
-  if (!is.vector(Param) | !is.numeric(Param)) {
-    stop("'Param' must be a numeric vector")
-  }
-  if (sum(!is.na(Param)) != NParam) {
-    stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
-  }
+  .ArgumentsCheckGR(InputsModel, RunOptions, Param)
+
   Param <- as.double(Param)
 
 
@@ -76,9 +48,9 @@ RunModel_CemaNeigeGR4H <- function(InputsModel, RunOptions, Param) {
   ## CemaNeige________________________________________________________________________________
   if (inherits(RunOptions, "CemaNeige")) {
     if ("all" %in% RunOptions$Outputs_Sim) {
-      IndOutputsCemaNeige <- as.integer(1:length(FortranOutputs$CN))
+      IndOutputsCemaNeige <- as.integer(1:length(RunOptions$FortranOutputs$CN))
     } else {
-      IndOutputsCemaNeige <- which(FortranOutputs$CN %in% RunOptions$Outputs_Sim)
+      IndOutputsCemaNeige <- which(RunOptions$FortranOutputs$CN %in% RunOptions$Outputs_Sim)
     }
     CemaNeigeLayers <- list()
     CemaNeigeStateEnd <- NULL
@@ -116,7 +88,7 @@ RunModel_CemaNeigeGR4H <- function(InputsModel, RunOptions, Param) {
 
       ## Data storage
       CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
-      names(CemaNeigeLayers[[iLayer]]) <- FortranOutputs$CN[IndOutputsCemaNeige]
+      names(CemaNeigeLayers[[iLayer]]) <- RunOptions$FortranOutputs$CN[IndOutputsCemaNeige]
       IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt")
       if (iLayer == 1) {
         CatchMeltAndPliq <- RESULTS$Outputs[, IndPliqAndMelt] / NLayers
@@ -142,9 +114,9 @@ RunModel_CemaNeigeGR4H <- function(InputsModel, RunOptions, Param) {
 
   ## GR model
   if ("all" %in% RunOptions$Outputs_Sim) {
-    IndOutputsMod <- as.integer(1:length(FortranOutputs$GR))
+    IndOutputsMod <- as.integer(1:length(RunOptions$FortranOutputs$GR))
   } else {
-    IndOutputsMod <- which(FortranOutputs$GR %in% RunOptions$Outputs_Sim)
+    IndOutputsMod <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
   }
 
   ## Use of IniResLevels
@@ -186,45 +158,14 @@ RunModel_CemaNeigeGR4H <- function(InputsModel, RunOptions, Param) {
   }
 
   if (inherits(RunOptions, "CemaNeige") & "Precip" %in% RunOptions$Outputs_Sim) {
-    RESULTS$Outputs[, which(FortranOutputs$GR[IndOutputsMod] == "Precip")] <- InputsModel$Precip[IndPeriod1]
-  }
-
-  ## Output data preparation
-  ## OutputsModel only
-  if (!ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(CemaNeigeLayers))
-    names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers)
-  }
-  ## DatesR and OutputsModel only
-  if (ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(CemaNeigeLayers))
-    names(OutputsModel) <- c("DatesR", FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers)
-  }
-  ## OutputsModel and StateEnd only
-  if (!ExportDatesR & ExportStateEnd) {
-    OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(CemaNeigeLayers),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers, "StateEnd")
-  }
-  ## DatesR and OutputsModel and StateEnd
-  if (ExportDatesR & ExportStateEnd) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(CemaNeigeLayers),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c("DatesR", FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers, "StateEnd")
-  }
-
-  ## End
-  rm(RESULTS)
-  class(OutputsModel) <- c("OutputsModel", "hourly", "GR", "CemaNeige")
-  if (IsHyst) {
-    class(OutputsModel) <- c(class(OutputsModel), "hysteresis")
+    RESULTS$Outputs[, which(RunOptions$FortranOutputs$GR[IndOutputsMod] == "Precip")] <-
+      InputsModel$Precip[IndPeriod1]
   }
-  return(OutputsModel)
 
+  ## OutputsModel generation
+  .GetOutputsModelGR(InputsModel,
+                     RunOptions,
+                     RESULTS,
+                     LInputSeries,
+                     CemaNeigeLayers)
 }
diff --git a/R/RunModel_CemaNeigeGR4J.R b/R/RunModel_CemaNeigeGR4J.R
index 6ff79a433e234a1aa50aa0cb692e8e9c845e729d..34f15d5ea3c3a546e30a5de3b857c4484e315194 100644
--- a/R/RunModel_CemaNeigeGR4J.R
+++ b/R/RunModel_CemaNeigeGR4J.R
@@ -1,44 +1,16 @@
 RunModel_CemaNeigeGR4J <- function(InputsModel, RunOptions, Param) {
-  
-  
+
+
   ## Initialization of variables
   IsHyst <- inherits(RunOptions, "hysteresis")
-  NParam <- ifelse(test = IsHyst, yes = 8L, no = 6L)
-  NParamCN <- NParam - 4L
+  NParamCN <- RunOptions$FeatFUN_MOD$NbParam - 4L
   NStates <- 4L
-  FortranOutputs <- .FortranOutputs(GR = "GR4J", isCN = TRUE)
-  
-  
-  ## Arguments check
-  if (!inherits(InputsModel, "InputsModel")) {
-    stop("'InputsModel' must be of class 'InputsModel'")
-  }
-  if (!inherits(InputsModel, "daily")) {
-    stop("'InputsModel' must be of class 'daily'")
-  }
-  if (!inherits(InputsModel, "GR")) {
-    stop("'InputsModel' must be of class 'GR'")
-  }
-  if (!inherits(InputsModel, "CemaNeige")) {
-    stop("'InputsModel' must be of class 'CemaNeige'")
-  }
-  if (!inherits(RunOptions, "RunOptions")) {
-    stop("'RunOptions' must be of class 'RunOptions'")
-  }
-  if (!inherits(RunOptions, "GR")) {
-    stop("'RunOptions' must be of class 'GR'")
-  }
-  if (!inherits(RunOptions, "CemaNeige")) {
-    stop("'RunOptions' must be of class 'CemaNeige'")
-  }
-  if (!is.vector(Param) | !is.numeric(Param)) {
-    stop("'Param' must be a numeric vector")
-  }
-  if (sum(!is.na(Param)) != NParam) {
-    stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
-  }
+
+
+  .ArgumentsCheckGR(InputsModel, RunOptions, Param)
+
   Param <- as.double(Param)
-  
+
 
   Param_X1X3_threshold <- 1e-2
   Param_X4_threshold   <- 0.5
@@ -53,8 +25,8 @@ RunModel_CemaNeigeGR4J <- function(InputsModel, RunOptions, Param) {
   if (Param[4L] < Param_X4_threshold) {
     warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
     Param[4L] <- Param_X4_threshold
-  }      
-  
+  }
+
   ## Input data preparation
   if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
     RunOptions$IndPeriod_WarmUp <- NULL
@@ -71,28 +43,28 @@ RunModel_CemaNeigeGR4J <- function(InputsModel, RunOptions, Param) {
   ## Output data preparation
   ExportDatesR   <- "DatesR"   %in% RunOptions$Outputs_Sim
   ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
-  
-  
+
+
   ## CemaNeige________________________________________________________________________________
   if (inherits(RunOptions, "CemaNeige")) {
     if ("all" %in% RunOptions$Outputs_Sim) {
-      IndOutputsCemaNeige <- as.integer(1:length(FortranOutputs$CN)) 
+      IndOutputsCemaNeige <- as.integer(1:length(RunOptions$FortranOutputs$CN))
     } else {
-      IndOutputsCemaNeige <- which(FortranOutputs$CN %in% RunOptions$Outputs_Sim)
+      IndOutputsCemaNeige <- which(RunOptions$FortranOutputs$CN %in% RunOptions$Outputs_Sim)
     }
     CemaNeigeLayers <- list()
     CemaNeigeStateEnd <- NULL
     NameCemaNeigeLayers <- "CemaNeigeLayers"
-    
-    
+
+
     ## Call CemaNeige Fortran_________________________
-    for(iLayer in 1:NLayers) {
+    for (iLayer in 1:NLayers) {
       if (!IsHyst) {
         StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers)]
       } else {
         StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
       }
-      RESULTS <- .Fortran("frun_cemaneige", PACKAGE = "airGR", 
+      RESULTS <- .Fortran("frun_cemaneige", PACKAGE = "airGR",
                           ## inputs
                           LInputs = LInputSeries,                                                         ### length of input and output series
                           InputsPrecip = InputsModel$LayerPrecip[[iLayer]][IndPeriod1],                   ### input series of total precipitation [mm/d]
@@ -106,16 +78,16 @@ RunModel_CemaNeigeGR4J <- function(InputsModel, RunOptions, Param) {
                           IsHyst = as.integer(IsHyst),                                                    ### use of hysteresis
                           NOutputs = as.integer(length(IndOutputsCemaNeige)),                             ### number of output series
                           IndOutputs = IndOutputsCemaNeige,                                               ### indices of output series
-                          ## outputs                                                               
+                          ## outputs
                           Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsCemaNeige)), ### output series [mm, mm/d or degC]
                           StateEnd = rep(as.double(-99e9), as.integer(NStates))                                        ### state variables at the end of the model run
       )
       RESULTS$Outputs[RESULTS$Outputs   <= -99e8] <- NA
       RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
-      
+
       ## Data storage
       CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
-      names(CemaNeigeLayers[[iLayer]]) <- FortranOutputs$CN[IndOutputsCemaNeige]
+      names(CemaNeigeLayers[[iLayer]]) <- RunOptions$FortranOutputs$CN[IndOutputsCemaNeige]
       IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt")
       if (iLayer == 1) {
         CatchMeltAndPliq <- RESULTS$Outputs[, IndPliqAndMelt] / NLayers
@@ -136,24 +108,24 @@ RunModel_CemaNeigeGR4J <- function(InputsModel, RunOptions, Param) {
     NameCemaNeigeLayers <- NULL
     CatchMeltAndPliq <- InputsModel$Precip[IndPeriod1]
   }
-  
-  
-  
+
+
+
   ## GR model______________________________________________________________________________________
   if ("all" %in% RunOptions$Outputs_Sim) {
-    IndOutputsMod <- as.integer(1:length(FortranOutputs$GR))
+    IndOutputsMod <- as.integer(1:length(RunOptions$FortranOutputs$GR))
   } else {
-    IndOutputsMod <- which(FortranOutputs$GR %in% RunOptions$Outputs_Sim)
+    IndOutputsMod <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
   }
-  
+
   ## Use of IniResLevels
   if (!is.null(RunOptions$IniResLevels)) {
     RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * ParamMod[1] ### production store level (mm)
     RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * ParamMod[3] ### routing store level (mm)
   }
-  
+
   ## Call GR model Fortan
-  RESULTS <- .Fortran("frun_gr4j", PACKAGE = "airGR", 
+  RESULTS <- .Fortran("frun_gr4j", PACKAGE = "airGR",
                       ## inputs
                       LInputs = LInputSeries,                          ### length of input and output series
                       InputsPrecip = CatchMeltAndPliq,                 ### input series of total precipitation [mm/d]
@@ -164,7 +136,7 @@ RunModel_CemaNeigeGR4J <- function(InputsModel, RunOptions, Param) {
                       StateStart = RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts
                       NOutputs = as.integer(length(IndOutputsMod)),    ### number of output series
                       IndOutputs = IndOutputsMod,                      ### indices of output series
-                      ## outputs                                        
+                      ## outputs
                       Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsMod)), ### output series [mm or mm/d]
                       StateEnd = rep(as.double(-99e9), NStatesMod)                                           ### state variables at the end of the model run
   )
@@ -173,57 +145,26 @@ RunModel_CemaNeigeGR4J <- function(InputsModel, RunOptions, Param) {
   if (ExportStateEnd) {
     RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
     idNStates <- seq_len(NStates*NLayers) %% NStates
-    RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR4J, InputsModel = InputsModel, IsHyst = IsHyst, 
-                                        ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL, 
+    RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR4J, InputsModel = InputsModel, IsHyst = IsHyst,
+                                        ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
                                         UH1 = RESULTS$StateEnd[(1:20) + 7],
-                                        UH2 = RESULTS$StateEnd[(1:40) + (7+20)], 
-                                        GCemaNeigeLayers       = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]], 
-                                        eTGCemaNeigeLayers     = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]], 
-                                        GthrCemaNeigeLayers    = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]], 
-                                        GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]], 
+                                        UH2 = RESULTS$StateEnd[(1:40) + (7+20)],
+                                        GCemaNeigeLayers       = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
+                                        eTGCemaNeigeLayers     = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
+                                        GthrCemaNeigeLayers    = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
+                                        GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
                                         verbose = FALSE)
   }
-  
+
   if (inherits(RunOptions, "CemaNeige") & "Precip" %in% RunOptions$Outputs_Sim) {
-    RESULTS$Outputs[, which(FortranOutputs$GR[IndOutputsMod] == "Precip")] <- InputsModel$Precip[IndPeriod1]
-  }
-  
-  ## Output data preparation
-  ## OutputsModel only
-  if (!ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]), 
-                      list(CemaNeigeLayers))
-    names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers)
-  }
-  ## DatesR and OutputsModel only
-  if (ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), 
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]), 
-                      list(CemaNeigeLayers))
-    names(OutputsModel) <- c("DatesR", FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers)
+    RESULTS$Outputs[, which(RunOptions$FortranOutputs$GR[IndOutputsMod] == "Precip")] <-
+      InputsModel$Precip[IndPeriod1]
   }
-  ## OutputsModel and StateEnd only
-  if (!ExportDatesR & ExportStateEnd) {
-    OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]), 
-                      list(CemaNeigeLayers), 
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers, "StateEnd")
-  }
-  ## DatesR and OutputsModel and StateEnd
-  if (ExportDatesR & ExportStateEnd) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), 
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]), 
-                      list(CemaNeigeLayers), 
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c("DatesR", FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers, "StateEnd")
-  }
-  
-  ## End
-  rm(RESULTS)
-  class(OutputsModel) <- c("OutputsModel", "daily", "GR", "CemaNeige")
-  if (IsHyst) {
-    class(OutputsModel) <- c(class(OutputsModel), "hysteresis")
-  }
-  return(OutputsModel)
-  
+
+  ## OutputsModel generation
+  .GetOutputsModelGR(InputsModel,
+                     RunOptions,
+                     RESULTS,
+                     LInputSeries,
+                     CemaNeigeLayers)
 }
diff --git a/R/RunModel_CemaNeigeGR5H.R b/R/RunModel_CemaNeigeGR5H.R
index 0a4995fb2b570e96e2ee593bc566c9812f737363..a78946b618549ee741ee2c06daf58000ea147505 100644
--- a/R/RunModel_CemaNeigeGR5H.R
+++ b/R/RunModel_CemaNeigeGR5H.R
@@ -1,51 +1,22 @@
 RunModel_CemaNeigeGR5H <- function(InputsModel, RunOptions, Param) {
-  
-  
+
+
   ## Initialization of variables
   IsHyst <- inherits(RunOptions, "hysteresis")
-  NParam <- ifelse(test = IsHyst, yes = 9L, no = 7L)
-  NParamCN <- NParam - 5L
+  NParamCN <- RunOptions$FeatFUN_MOD$NbParam - 5L
   NStates <- 4L
-  FortranOutputs <- .FortranOutputs(GR = "GR5H", isCN = TRUE)
   IsIntStore <- inherits(RunOptions, "interception")
   if (IsIntStore) {
     Imax <- RunOptions$Imax
   } else {
     Imax <- -99
   }
-  
-  
-  ## Arguments check
-  if (!inherits(InputsModel, "InputsModel")) {
-    stop("'InputsModel' must be of class 'InputsModel'")
-  }
-  if (!inherits(InputsModel, "hourly")) {
-    stop("'InputsModel' must be of class 'hourly'")
-  }
-  if (!inherits(InputsModel, "GR")) {
-    stop("'InputsModel' must be of class 'GR'")
-  }
-  if (!inherits(InputsModel, "CemaNeige")) {
-    stop("'InputsModel' must be of class 'CemaNeige'")
-  }
-  if (!inherits(RunOptions, "RunOptions")) {
-    stop("'RunOptions' must be of class 'RunOptions'")
-  }
-  if (!inherits(RunOptions, "GR")) {
-    stop("'RunOptions' must be of class 'GR'")
-  }
-  if (!inherits(RunOptions, "CemaNeige")) {
-    stop("'RunOptions' must be of class 'CemaNeige'")
-  }
-  if (!is.vector(Param) | !is.numeric(Param)) {
-    stop("'Param' must be a numeric vector")
-  }
-  if (sum(!is.na(Param)) != NParam) {
-    stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
-  }
+
+  .ArgumentsCheckGR(InputsModel, RunOptions, Param)
+
   Param <- as.double(Param)
-  
-  
+
+
   Param_X1X3_threshold <- 1e-2
   Param_X4_threshold   <- 0.5
   if (Param[1L] < Param_X1X3_threshold) {
@@ -59,8 +30,8 @@ RunModel_CemaNeigeGR5H <- function(InputsModel, RunOptions, Param) {
   if (Param[4L] < Param_X4_threshold) {
     warning(sprintf("Param[4] (X4: unit hydrograph time constant [h]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
     Param[4L] <- Param_X4_threshold
-  }     
-  
+  }
+
   ## Input data preparation
   if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
     RunOptions$IndPeriod_WarmUp <- NULL
@@ -73,24 +44,24 @@ RunModel_CemaNeigeGR5H <- function(InputsModel, RunOptions, Param) {
   ParamMod       <- Param[1:NParamMod]
   NLayers        <- length(InputsModel$LayerPrecip)
   NStatesMod     <- as.integer(length(RunOptions$IniStates) - NStates * NLayers)
-  
+
   ## Output data preparation
   ExportDatesR   <- "DatesR"   %in% RunOptions$Outputs_Sim
   ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
-  
-  
+
+
   ## CemaNeige________________________________________________________________________________
   if (inherits(RunOptions, "CemaNeige")) {
     if ("all" %in% RunOptions$Outputs_Sim) {
-      IndOutputsCemaNeige <- as.integer(1:length(FortranOutputs$CN))
+      IndOutputsCemaNeige <- as.integer(1:length(RunOptions$FortranOutputs$CN))
     } else {
-      IndOutputsCemaNeige <- which(FortranOutputs$CN %in% RunOptions$Outputs_Sim)
+      IndOutputsCemaNeige <- which(RunOptions$FortranOutputs$CN %in% RunOptions$Outputs_Sim)
     }
     CemaNeigeLayers <- list()
     CemaNeigeStateEnd <- NULL
     NameCemaNeigeLayers <- "CemaNeigeLayers"
-    
-    
+
+
     ## Call CemaNeige Fortran_________________________
     for (iLayer in 1:NLayers) {
 
@@ -113,16 +84,16 @@ RunModel_CemaNeigeGR5H <- function(InputsModel, RunOptions, Param) {
                           IsHyst = as.integer(IsHyst),                                                    ### use of hysteresis
                           NOutputs = as.integer(length(IndOutputsCemaNeige)),                             ### number of output series
                           IndOutputs = IndOutputsCemaNeige,                                               ### indices of output series
-                          ## outputs                                                               
+                          ## outputs
                           Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsCemaNeige)), ### output series [mm, mm/h or degC]
                           StateEnd = rep(as.double(-99e9), as.integer(NStates))                                        ### state variables at the end of the model run
       )
       RESULTS$Outputs[RESULTS$Outputs   <= -99e8] <- NA
       RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
-      
+
       ## Data storage
       CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
-      names(CemaNeigeLayers[[iLayer]]) <- FortranOutputs$CN[IndOutputsCemaNeige]
+      names(CemaNeigeLayers[[iLayer]]) <- RunOptions$FortranOutputs$CN[IndOutputsCemaNeige]
       IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt")
       if (iLayer == 1) {
         CatchMeltAndPliq <- RESULTS$Outputs[, IndPliqAndMelt] / NLayers
@@ -148,9 +119,9 @@ RunModel_CemaNeigeGR5H <- function(InputsModel, RunOptions, Param) {
 
   ## GR model
   if ("all" %in% RunOptions$Outputs_Sim) {
-    IndOutputsMod <- as.integer(1:length(FortranOutputs$GR))
+    IndOutputsMod <- as.integer(1:length(RunOptions$FortranOutputs$GR))
   } else {
-    IndOutputsMod <- which(FortranOutputs$GR %in% RunOptions$Outputs_Sim)
+    IndOutputsMod <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
   }
 
   ## Use of IniResLevels
@@ -191,54 +162,20 @@ RunModel_CemaNeigeGR5H <- function(InputsModel, RunOptions, Param) {
                                         UH2 = RESULTS$StateEnd[(1:(40*24)) + (7+20*24)],
                                         GCemaNeigeLayers       = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
                                         eTGCemaNeigeLayers     = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
-                                        GthrCemaNeigeLayers    = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]], 
+                                        GthrCemaNeigeLayers    = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
                                         GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
                                         verbose = FALSE)
   }
-  
+
   if (inherits(RunOptions, "CemaNeige") & "Precip" %in% RunOptions$Outputs_Sim) {
-    RESULTS$Outputs[, which(FortranOutputs$GR[IndOutputsMod] == "Precip")] <- InputsModel$Precip[IndPeriod1]
+    RESULTS$Outputs[, which(RunOptions$FortranOutputs$GR[IndOutputsMod] == "Precip")] <-
+      InputsModel$Precip[IndPeriod1]
   }
-  
-  ## Output data preparation
-  ## OutputsModel only
-  if (!ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(CemaNeigeLayers))
-    names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers)
-  }
-  ## DatesR and OutputsModel only
-  if (ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(CemaNeigeLayers))
-    names(OutputsModel) <- c("DatesR", FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers)
-  }
-  ## OutputsModel and StateEnd only
-  if (!ExportDatesR & ExportStateEnd) {
-    OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(CemaNeigeLayers),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers, "StateEnd")
-  }
-  ## DatesR and OutputsModel and StateEnd
-  if (ExportDatesR & ExportStateEnd) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(CemaNeigeLayers),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c("DatesR", FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers, "StateEnd")
-  }
-  
-  ## End
-  rm(RESULTS)
-  class(OutputsModel) <- c("OutputsModel", "hourly", "GR", "CemaNeige")
-  if (IsHyst) {
-    class(OutputsModel) <- c(class(OutputsModel), "hysteresis")
-  }
-  if (IsIntStore) {
-    class(OutputsModel) <- c(class(OutputsModel), "interception")
-  }
-  return(OutputsModel)
-  
+
+  ## OutputsModel generation
+  .GetOutputsModelGR(InputsModel,
+                     RunOptions,
+                     RESULTS,
+                     LInputSeries,
+                     CemaNeigeLayers)
 }
diff --git a/R/RunModel_CemaNeigeGR5J.R b/R/RunModel_CemaNeigeGR5J.R
index a8dabc3c97764c99f6e9c69f14dc951ba828e3fe..cdc4faaeae4263dd5a54f3559f67ec7a35b4b899 100644
--- a/R/RunModel_CemaNeigeGR5J.R
+++ b/R/RunModel_CemaNeigeGR5J.R
@@ -1,45 +1,17 @@
 RunModel_CemaNeigeGR5J <- function(InputsModel, RunOptions, Param) {
-  
-  
+
+
   ## Initialization of variables
   IsHyst <- inherits(RunOptions, "hysteresis")
-  NParam <- ifelse(test = IsHyst, yes = 9L, no = 7L)
-  NParamCN <- NParam - 5L
+  NParamCN <- RunOptions$FeatFUN_MOD$NbParam - 5L
   NStates <- 4L
-  FortranOutputs <- .FortranOutputs(GR = "GR5J", isCN = TRUE)
-  
-  
-  ## Arguments check
-  if (!inherits(InputsModel, "InputsModel")) {
-    stop("'InputsModel' must be of class 'InputsModel'")
-  }
-  if (!inherits(InputsModel, "daily")) {
-    stop("'InputsModel' must be of class 'daily'")
-  }
-  if (!inherits(InputsModel, "GR")) {
-    stop("'InputsModel' must be of class 'GR'")
-  }
-  if (!inherits(InputsModel, "CemaNeige")) {
-    stop("'InputsModel' must be of class 'CemaNeige'")
-  }
-  if (!inherits(RunOptions, "RunOptions")) {
-    stop("'RunOptions' must be of class 'RunOptions'")
-  }
-  if (!inherits(RunOptions, "GR")) {
-    stop("'RunOptions' must be of class 'GR'")
-  }
-  if (!inherits(RunOptions, "CemaNeige")) {
-    stop("'RunOptions' must be of class 'CemaNeige'")
-  }
-  if (!is.vector(Param) | !is.numeric(Param)) {
-    stop("'Param' must be a numeric vector")
-  }
-  if (sum(!is.na(Param)) != NParam) {
-    stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
-  }
+
+
+  .ArgumentsCheckGR(InputsModel, RunOptions, Param)
+
   Param <- as.double(Param)
-  
-  
+
+
   Param_X1X3_threshold <- 1e-2
   Param_X4_threshold   <- 0.5
   if (Param[1L] < Param_X1X3_threshold) {
@@ -53,8 +25,8 @@ RunModel_CemaNeigeGR5J <- function(InputsModel, RunOptions, Param) {
   if (Param[4L] < Param_X4_threshold) {
     warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
     Param[4L] <- Param_X4_threshold
-  }      
-  
+  }
+
   ## Input data preparation
   if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
     RunOptions$IndPeriod_WarmUp <- NULL
@@ -69,22 +41,22 @@ RunModel_CemaNeigeGR5J <- function(InputsModel, RunOptions, Param) {
   NStatesMod     <- as.integer(length(RunOptions$IniStates) - NStates * NLayers)
   ExportDatesR   <- "DatesR"   %in% RunOptions$Outputs_Sim
   ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
-  
-  
+
+
   ## CemaNeige________________________________________________________________________________
   if (inherits(RunOptions, "CemaNeige")) {
     if ("all" %in% RunOptions$Outputs_Sim) {
-      IndOutputsCemaNeige <- as.integer(1:length(FortranOutputs$CN))
+      IndOutputsCemaNeige <- as.integer(1:length(RunOptions$FortranOutputs$CN))
     } else {
-      IndOutputsCemaNeige <- which(FortranOutputs$CN %in% RunOptions$Outputs_Sim)
+      IndOutputsCemaNeige <- which(RunOptions$FortranOutputs$CN %in% RunOptions$Outputs_Sim)
     }
     CemaNeigeLayers <- list()
     CemaNeigeStateEnd <- NULL
     NameCemaNeigeLayers <- "CemaNeigeLayers"
-    
-    
+
+
     ## Call CemaNeige Fortran_________________________
-    for(iLayer in 1:NLayers) {
+    for (iLayer in 1:NLayers) {
       if (!IsHyst) {
         StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers)]
       } else {
@@ -104,20 +76,20 @@ RunModel_CemaNeigeGR5J <- function(InputsModel, RunOptions, Param) {
                           IsHyst = as.integer(IsHyst),                                                    ### use of hysteresis
                           NOutputs = as.integer(length(IndOutputsCemaNeige)),                             ### number of output series
                           IndOutputs = IndOutputsCemaNeige,                                               ### indices of output series
-                          ## outputs                                                               
+                          ## outputs
                           Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsCemaNeige)), ### output series [mm, mm/d or degC]
                           StateEnd = rep(as.double(-99e9), as.integer(NStates))                                        ### state variables at the end of the model run
       )
       RESULTS$Outputs[RESULTS$Outputs   <= -99e8] <- NA
       RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
-      
+
       ## Data storage
       CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
-      names(CemaNeigeLayers[[iLayer]]) <- FortranOutputs$CN[IndOutputsCemaNeige]
+      names(CemaNeigeLayers[[iLayer]]) <- RunOptions$FortranOutputs$CN[IndOutputsCemaNeige]
       IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt")
       if (iLayer == 1) {
         CatchMeltAndPliq <- RESULTS$Outputs[, IndPliqAndMelt] / NLayers
-        }
+      }
       if (iLayer > 1) {
         CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[, IndPliqAndMelt] / NLayers
       }
@@ -133,23 +105,23 @@ RunModel_CemaNeigeGR5J <- function(InputsModel, RunOptions, Param) {
     CemaNeigeStateEnd <- NULL
     NameCemaNeigeLayers <- NULL
     CatchMeltAndPliq <- InputsModel$Precip[IndPeriod1]
-    }
-  
-  
-  
+  }
+
+
+
   ## GR model______________________________________________________________________________________
   if ("all" %in% RunOptions$Outputs_Sim) {
-    IndOutputsMod <- as.integer(1:length(FortranOutputs$GR))
+    IndOutputsMod <- as.integer(1:length(RunOptions$FortranOutputs$GR))
   } else {
-    IndOutputsMod <- which(FortranOutputs$GR %in% RunOptions$Outputs_Sim)
+    IndOutputsMod <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
   }
-  
+
   ## Use of IniResLevels
   if (!is.null(RunOptions$IniResLevels)) {
     RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * ParamMod[1] ### production store level (mm)
     RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * ParamMod[3] ### routing store level (mm)
   }
-  
+
   ## Call GR model Fortan
   RESULTS <- .Fortran("frun_gr5j", PACKAGE = "airGR",
                       ## inputs
@@ -162,7 +134,7 @@ RunModel_CemaNeigeGR5J <- function(InputsModel, RunOptions, Param) {
                       StateStart = RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts
                       NOutputs = as.integer(length(IndOutputsMod)),    ### number of output series
                       IndOutputs = IndOutputsMod,                      ### indices of output series
-                      ## outputs                                        
+                      ## outputs
                       Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsMod)), ### output series [mm or mm/d]
                       StateEnd = rep(as.double(-99e9), NStatesMod)                                           ### state variables at the end of the model run
   )
@@ -177,51 +149,20 @@ RunModel_CemaNeigeGR5J <- function(InputsModel, RunOptions, Param) {
                                         UH2 = RESULTS$StateEnd[(1:40) + (7+20)],
                                         GCemaNeigeLayers       = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
                                         eTGCemaNeigeLayers     = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
-                                        GthrCemaNeigeLayers    = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]], 
+                                        GthrCemaNeigeLayers    = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
                                         GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
                                         verbose = FALSE)
   }
-  
+
   if (inherits(RunOptions, "CemaNeige") & "Precip" %in% RunOptions$Outputs_Sim) {
-    RESULTS$Outputs[, which(FortranOutputs$GR[IndOutputsMod] == "Precip")] <- InputsModel$Precip[IndPeriod1]
-  }
-  
-  ## Output data preparation
-  ## OutputsModel only
-  if (!ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(CemaNeigeLayers))
-    names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers)
-  }
-  ## DatesR and OutputsModel only
-  if (ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(CemaNeigeLayers))
-    names(OutputsModel) <- c("DatesR", FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers)
-  }
-  ## OutputsModel and StateEnd only
-  if (!ExportDatesR & ExportStateEnd) {
-    OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(CemaNeigeLayers),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers, "StateEnd")
-  }
-  ## DatesR and OutputsModel and StateEnd
-  if (ExportDatesR & ExportStateEnd) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(CemaNeigeLayers),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c("DatesR", FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers, "StateEnd")
-  }
-  
-  ## End
-  rm(RESULTS)
-  class(OutputsModel) <- c("OutputsModel", "daily", "GR", "CemaNeige")
-  if (IsHyst) {
-    class(OutputsModel) <- c(class(OutputsModel), "hysteresis")
-  }
-  return(OutputsModel)
-  
+    RESULTS$Outputs[, which(RunOptions$FortranOutputs$GR[IndOutputsMod] == "Precip")] <-
+      InputsModel$Precip[IndPeriod1]
+  }
+
+  ## OutputsModel generation
+  .GetOutputsModelGR(InputsModel,
+                     RunOptions,
+                     RESULTS,
+                     LInputSeries,
+                     CemaNeigeLayers)
 }
diff --git a/R/RunModel_CemaNeigeGR6J.R b/R/RunModel_CemaNeigeGR6J.R
index 48859cb09a837e014909452794db9a5bbad05721..6199e4178d8e02f7ffd8930eee623105e8f3b51b 100644
--- a/R/RunModel_CemaNeigeGR6J.R
+++ b/R/RunModel_CemaNeigeGR6J.R
@@ -3,40 +3,12 @@ RunModel_CemaNeigeGR6J <- function(InputsModel, RunOptions, Param) {
 
   ## Initialization of variables
   IsHyst <- inherits(RunOptions, "hysteresis")
-  NParam <- ifelse(test = IsHyst, yes = 10L, no = 8L)
-  NParamCN <- NParam - 6L
+  NParamCN <- RunOptions$FeatFUN_MOD$NbParam - 6L
   NStates <- 4L
-  FortranOutputs <- .FortranOutputs(GR = "GR6J", isCN = TRUE)
 
 
-  ## Arguments check
-  if (!inherits(InputsModel, "InputsModel")) {
-    stop("'InputsModel' must be of class 'InputsModel'")
-  }
-  if (!inherits(InputsModel, "daily")) {
-    stop("'InputsModel' must be of class 'daily'")
-  }
-  if (!inherits(InputsModel, "GR")) {
-    stop("'InputsModel' must be of class 'GR'")
-  }
-  if (!inherits(InputsModel, "CemaNeige")) {
-    stop("'InputsModel' must be of class 'CemaNeige'")
-  }
-  if (!inherits(RunOptions, "RunOptions")) {
-    stop("'RunOptions' must be of class 'RunOptions'")
-  }
-  if (!inherits(RunOptions, "GR")) {
-    stop("'RunOptions' must be of class 'GR'")
-  }
-  if (!inherits(RunOptions, "CemaNeige")) {
-    stop("'RunOptions' must be of class 'CemaNeige'")
-  }
-  if (!is.vector(Param) | !is.numeric(Param)) {
-    stop("'Param' must be a numeric vector")
-  }
-  if (sum(!is.na(Param)) != NParam) {
-    stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
-  }
+  .ArgumentsCheckGR(InputsModel, RunOptions, Param)
+
   Param <- as.double(Param)
 
 
@@ -78,9 +50,9 @@ RunModel_CemaNeigeGR6J <- function(InputsModel, RunOptions, Param) {
   ## CemaNeige________________________________________________________________________________
   if (inherits(RunOptions, "CemaNeige")) {
     if ("all" %in% RunOptions$Outputs_Sim) {
-      IndOutputsCemaNeige <- as.integer(1:length(FortranOutputs$CN))
+      IndOutputsCemaNeige <- as.integer(1:length(RunOptions$FortranOutputs$CN))
     } else {
-      IndOutputsCemaNeige <- which(FortranOutputs$CN %in% RunOptions$Outputs_Sim)
+      IndOutputsCemaNeige <- which(RunOptions$FortranOutputs$CN %in% RunOptions$Outputs_Sim)
     }
     CemaNeigeLayers <- list()
     CemaNeigeStateEnd <- NULL
@@ -88,7 +60,7 @@ RunModel_CemaNeigeGR6J <- function(InputsModel, RunOptions, Param) {
 
 
     ## Call CemaNeige Fortran_________________________
-    for(iLayer in 1:NLayers) {
+    for (iLayer in 1:NLayers) {
       if (!IsHyst) {
         StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers)]
       } else {
@@ -117,7 +89,7 @@ RunModel_CemaNeigeGR6J <- function(InputsModel, RunOptions, Param) {
 
       ## Data storage
       CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
-      names(CemaNeigeLayers[[iLayer]]) <- FortranOutputs$CN[IndOutputsCemaNeige]
+      names(CemaNeigeLayers[[iLayer]]) <- RunOptions$FortranOutputs$CN[IndOutputsCemaNeige]
       IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt")
       if (iLayer == 1) {
         CatchMeltAndPliq <- RESULTS$Outputs[, IndPliqAndMelt] / NLayers
@@ -143,9 +115,9 @@ RunModel_CemaNeigeGR6J <- function(InputsModel, RunOptions, Param) {
 
   ## GR model______________________________________________________________________________________
   if ("all" %in% RunOptions$Outputs_Sim) {
-    IndOutputsMod <- as.integer(1:length(FortranOutputs$GR))
+    IndOutputsMod <- as.integer(1:length(RunOptions$FortranOutputs$GR))
   } else {
-    IndOutputsMod <- which(FortranOutputs$GR %in% RunOptions$Outputs_Sim)
+    IndOutputsMod <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
   }
 
   ## Use of IniResLevels
@@ -188,46 +160,15 @@ RunModel_CemaNeigeGR6J <- function(InputsModel, RunOptions, Param) {
   }
 
   if (inherits(RunOptions, "CemaNeige") & "Precip" %in% RunOptions$Outputs_Sim) {
-    RESULTS$Outputs[, which(FortranOutputs$GR[IndOutputsMod] == "Precip")] <- InputsModel$Precip[IndPeriod1]
-  }
-
-  ## Output data preparation
-  ## OutputsModel only
-  if (!ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(CemaNeigeLayers))
-    names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers)
-  }
-  ## DatesR and OutputsModel only
-  if (ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(CemaNeigeLayers))
-    names(OutputsModel) <- c("DatesR", FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers)
-  }
-  ## OutputsModel and StateEnd only
-  if (!ExportDatesR & ExportStateEnd) {
-    OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(CemaNeigeLayers),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers, "StateEnd")
-  }
-  ## DatesR and OutputsModel and StateEnd
-  if (ExportDatesR & ExportStateEnd) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(CemaNeigeLayers),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c("DatesR", FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers, "StateEnd")
-  }
-
-  ## End
-  rm(RESULTS)
-  class(OutputsModel) <- c("OutputsModel", "daily", "GR", "CemaNeige")
-  if (IsHyst) {
-    class(OutputsModel) <- c(class(OutputsModel), "hysteresis")
+    RESULTS$Outputs[, which(RunOptions$FortranOutputs$GR[IndOutputsMod] == "Precip")] <-
+      InputsModel$Precip[IndPeriod1]
   }
-  return(OutputsModel)
 
+  ## OutputsModel generation
+  .GetOutputsModelGR(InputsModel,
+                     RunOptions,
+                     RESULTS,
+                     LInputSeries,
+                     CemaNeigeLayers)
 }
 
diff --git a/R/RunModel_GR1A.R b/R/RunModel_GR1A.R
index 620965f0c58c6a3ce8edcf9bd2801d3fc04433e9..87ff6aceb6d53f9b745f2542c2ed85b29f223b2f 100644
--- a/R/RunModel_GR1A.R
+++ b/R/RunModel_GR1A.R
@@ -1,33 +1,7 @@
 RunModel_GR1A <- function(InputsModel, RunOptions, Param) {
 
+  .ArgumentsCheckGR(InputsModel, RunOptions, Param)
 
-  ## Initialization of variables
-  NParam <- 1
-  FortranOutputs <- .FortranOutputs(GR = "GR1A")$GR
-
-
-  ## Arguments check
-  if (!inherits(InputsModel, "InputsModel")) {
-    stop("'InputsModel' must be of class 'InputsModel'")
-  }
-  if (!inherits(InputsModel, "yearly")) {
-    stop("'InputsModel' must be of class 'yearly'")
-  }
-  if (!inherits(InputsModel, "GR")) {
-    stop("'InputsModel' must be of class 'GR'")
-  }
-  if (!inherits(RunOptions, "RunOptions")) {
-    stop("'RunOptions' must be of class 'RunOptions'")
-  }
-  if (!inherits(RunOptions, "GR")) {
-    stop("'RunOptions' must be of class 'GR'")
-  }
-  if (!is.vector(Param) | !is.numeric(Param)) {
-    stop("'Param' must be a numeric vector")
-  }
-  if (sum(!is.na(Param)) != NParam) {
-    stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
-  }
   Param <- as.double(Param)
 
 
@@ -38,9 +12,9 @@ RunModel_GR1A <- function(InputsModel, RunOptions, Param) {
   IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
   LInputSeries <- as.integer(length(IndPeriod1))
   if ("all" %in% RunOptions$Outputs_Sim) {
-    IndOutputs <- as.integer(1:length(FortranOutputs))
+    IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
   } else {
-    IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim)
+    IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
   }
 
 
@@ -69,34 +43,9 @@ RunModel_GR1A <- function(InputsModel, RunOptions, Param) {
   RESULTS$Outputs[RESULTS$Outputs   <= -99e8] <- NA
   RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
 
-  ## Output data preparation
-  ## OutputsModel only
-  if (!ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
-    names(OutputsModel) <- FortranOutputs[IndOutputs]
-  }
-  ## DatesR and OutputsModel only
-  if (ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]))
-    names(OutputsModel) <- c("DatesR", FortranOutputs[IndOutputs])
-  }
-  ## OutputsModel and StateEnd only
-  if (!ExportDatesR & ExportStateEnd) {
-    OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c(FortranOutputs[IndOutputs], "StateEnd")
-  }
-  ## DatesR and OutputsModel and StateEnd
-  if ((ExportDatesR & ExportStateEnd) | "all" %in% RunOptions$Outputs_Sim) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c("DatesR", FortranOutputs[IndOutputs], "StateEnd")
-  }
-
-  ## End
-  class(OutputsModel) <- c("OutputsModel", "yearly", "GR")
-  return(OutputsModel)
-
+  ## OutputsModel generation
+  .GetOutputsModelGR(InputsModel,
+                     RunOptions,
+                     RESULTS,
+                     LInputSeries)
 }
diff --git a/R/RunModel_GR2M.R b/R/RunModel_GR2M.R
index e0261695d36e4b4c22b61f83c0e456df47f0a865..dd2f696d5cc30f0dd5e98f8f7a4198f4964aca6c 100644
--- a/R/RunModel_GR2M.R
+++ b/R/RunModel_GR2M.R
@@ -1,33 +1,7 @@
 RunModel_GR2M <- function(InputsModel, RunOptions, Param) {
 
+  .ArgumentsCheckGR(InputsModel, RunOptions, Param)
 
-  ## Initialization of variables
-  NParam <- 2
-  FortranOutputs <- .FortranOutputs(GR = "GR2M")$GR
-
-
-  ## Arguments check
-  if (!inherits(InputsModel, "InputsModel")) {
-    stop("'InputsModel' must be of class 'InputsModel'")
-  }
-  if (!inherits(InputsModel, "monthly")) {
-    stop("'InputsModel' must be of class 'monthly'")
-  }
-  if (!inherits(InputsModel, "GR")) {
-    stop("'InputsModel' must be of class 'GR'")
-  }
-  if (!inherits(RunOptions, "RunOptions")) {
-    stop("'RunOptions' must be of class 'RunOptions'")
-  }
-  if (!inherits(RunOptions, "GR")) {
-    stop("'RunOptions' must be of class 'GR'")
-  }
-  if (!is.vector(Param) | !is.numeric(Param)) {
-    stop("'Param' must be a numeric vector")
-  }
-  if (sum(!is.na(Param)) != NParam) {
-    stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
-  }
   Param <- as.double(Param)
 
   Param_X1X2_threshold <- 1e-2
@@ -47,9 +21,9 @@ RunModel_GR2M <- function(InputsModel, RunOptions, Param) {
   IndPeriod1   <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
   LInputSeries <- as.integer(length(IndPeriod1))
   if ("all" %in% RunOptions$Outputs_Sim) {
-    IndOutputs <- as.integer(1:length(FortranOutputs))
+    IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
   } else {
-    IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim)
+    IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
   }
 
   ## Output data preparation
@@ -91,34 +65,9 @@ RunModel_GR2M <- function(InputsModel, RunOptions, Param) {
   }
 
   ## Output data preparation
-  ## OutputsModel only
-  if (!ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
-    names(OutputsModel) <- FortranOutputs[IndOutputs]
-  }
-  ## DatesR and OutputsModel only
-  if (ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]))
-    names(OutputsModel) <- c("DatesR", FortranOutputs[IndOutputs])
-  }
-  ## OutputsModel and StateEnd only
-  if (!ExportDatesR & ExportStateEnd) {
-    OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c(FortranOutputs[IndOutputs], "StateEnd")
-  }
-  ## DatesR and OutputsModel and StateEnd
-  if ((ExportDatesR & ExportStateEnd) | "all" %in% RunOptions$Outputs_Sim) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                       lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                       list(RESULTS$StateEnd))
-    names(OutputsModel) <- c("DatesR", FortranOutputs[IndOutputs], "StateEnd")
-  }
-
-  ## End
-  rm(RESULTS)
-  class(OutputsModel) <- c("OutputsModel", "monthly", "GR")
-  return(OutputsModel)
-
+  ## OutputsModel generation
+  .GetOutputsModelGR(InputsModel,
+                     RunOptions,
+                     RESULTS,
+                     LInputSeries)
 }
diff --git a/R/RunModel_GR4H.R b/R/RunModel_GR4H.R
index 69beb4b05f3616dd0f8fb556130edb69d54ebb85..f70f3a1d0f165c51ef451e44d56ca50f83a685b6 100644
--- a/R/RunModel_GR4H.R
+++ b/R/RunModel_GR4H.R
@@ -1,33 +1,7 @@
 RunModel_GR4H <- function(InputsModel, RunOptions, Param) {
 
+  .ArgumentsCheckGR(InputsModel, RunOptions, Param)
 
-  ## Initialization of variables
-  NParam <- 4
-  FortranOutputs <- .FortranOutputs(GR = "GR4H")$GR
-
-
-  ## Arguments check
-  if (!inherits(InputsModel, "InputsModel")) {
-    stop("'InputsModel' must be of class 'InputsModel'")
-  }
-  if (!inherits(InputsModel, "hourly")) {
-    stop("'InputsModel' must be of class 'hourly'")
-  }
-  if (!inherits(InputsModel, "GR")) {
-    stop("'InputsModel' must be of class 'GR'")
-  }
-  if (!inherits(RunOptions, "RunOptions")) {
-    stop("'RunOptions' must be of class 'RunOptions'")
-  }
-  if (!inherits(RunOptions, "GR")) {
-    stop("'RunOptions' must be of class 'GR'")
-  }
-  if (!is.vector(Param) | !is.numeric(Param)) {
-    stop("'Param' must be a numeric vector")
-  }
-  if (sum(!is.na(Param)) != NParam) {
-    stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
-  }
   Param <- as.double(Param)
 
   Param_X1X3_threshold <- 1e-2
@@ -52,9 +26,9 @@ RunModel_GR4H <- function(InputsModel, RunOptions, Param) {
   IndPeriod1   <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
   LInputSeries <- as.integer(length(IndPeriod1))
   if ("all" %in% RunOptions$Outputs_Sim) {
-    IndOutputs <- as.integer(1:length(FortranOutputs))
+    IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
   } else {
-    IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim)
+    IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
   }
 
   ## Output data preparation
@@ -96,35 +70,9 @@ RunModel_GR4H <- function(InputsModel, RunOptions, Param) {
                                         verbose = FALSE)
   }
 
-  ## Output data preparation
-  ## OutputsModel only
-  if (!ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
-    names(OutputsModel) <- FortranOutputs[IndOutputs]
-  }
-  ## DatesR and OutputsModel only
-  if (ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]))
-    names(OutputsModel) <- c("DatesR", FortranOutputs[IndOutputs])
-  }
-  ## OutputsModel and StateEnd only
-  if (!ExportDatesR & ExportStateEnd) {
-    OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c(FortranOutputs[IndOutputs], "StateEnd")
-  }
-  ## DatesR and OutputsModel and StateEnd
-  if ((ExportDatesR & ExportStateEnd) | "all" %in% RunOptions$Outputs_Sim) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c("DatesR", FortranOutputs[IndOutputs], "StateEnd")
-  }
-
-  ## End
-  rm(RESULTS)
-  class(OutputsModel) <- c("OutputsModel", "hourly", "GR")
-  return(OutputsModel)
-
+  ## OutputsModel generation
+  .GetOutputsModelGR(InputsModel,
+                     RunOptions,
+                     RESULTS,
+                     LInputSeries)
 }
diff --git a/R/RunModel_GR4J.R b/R/RunModel_GR4J.R
index 0af262be0d397dcdd5a81533ac8e38bdf4746c22..fab8517e6383fbcfcf0244bd285d547ead1c3ea2 100644
--- a/R/RunModel_GR4J.R
+++ b/R/RunModel_GR4J.R
@@ -1,33 +1,7 @@
 RunModel_GR4J <- function(InputsModel, RunOptions, Param) {
 
+  .ArgumentsCheckGR(InputsModel, RunOptions, Param)
 
-  ## Initialization of variables
-  NParam <- 4
-  FortranOutputs <- .FortranOutputs(GR = "GR4J")$GR
-
-
-  ## Arguments check
-  if (!inherits(InputsModel, "InputsModel")) {
-    stop("'InputsModel' must be of class 'InputsModel'")
-  }
-  if (!inherits(InputsModel, "daily")) {
-    stop("'InputsModel' must be of class 'daily'")
-  }
-  if (!inherits(InputsModel, "GR")) {
-    stop("'InputsModel' must be of class 'GR'")
-  }
-  if (!inherits(RunOptions, "RunOptions")) {
-    stop("'RunOptions' must be of class 'RunOptions'")
-  }
-  if (!inherits(RunOptions, "GR")) {
-    stop("'RunOptions' must be of class 'GR'")
-  }
-  if (!is.vector(Param) | !is.numeric(Param)) {
-    stop("'Param' must be a numeric vector")
-  }
-  if (sum(!is.na(Param)) != NParam) {
-    stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
-  }
   Param <- as.double(Param)
 
   Param_X1X3_threshold <- 1e-2
@@ -52,16 +26,11 @@ RunModel_GR4J <- function(InputsModel, RunOptions, Param) {
   IndPeriod1   <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
   LInputSeries <- as.integer(length(IndPeriod1))
   if ("all" %in% RunOptions$Outputs_Sim) {
-    IndOutputs <- as.integer(1:length(FortranOutputs))
+    IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
   } else {
-    IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim)
+    IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
   }
 
-  ## Output data preparation
-  IndPeriod2     <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
-  ExportDatesR   <- "DatesR"   %in% RunOptions$Outputs_Sim
-  ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
-
   ## Use of IniResLevels
   if (!is.null(RunOptions$IniResLevels)) {
     RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
@@ -86,7 +55,7 @@ RunModel_GR4J <- function(InputsModel, RunOptions, Param) {
   )
   RESULTS$Outputs[RESULTS$Outputs   <= -99e8] <- NA
   RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
-  if (ExportStateEnd) {
+  if ("StateEnd" %in% RunOptions$Outputs_Sim) {
     RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
     RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
                                         ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
@@ -96,35 +65,9 @@ RunModel_GR4J <- function(InputsModel, RunOptions, Param) {
                                         verbose = FALSE)
   }
 
-  ## Output data preparation
-  ## OutputsModel only
-  if (!ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
-    names(OutputsModel) <- FortranOutputs[IndOutputs]
-  }
-  ## DatesR and OutputsModel only
-  if (ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]))
-    names(OutputsModel) <- c("DatesR", FortranOutputs[IndOutputs])
-  }
-  ## OutputsModel and StateEnd only
-  if (!ExportDatesR & ExportStateEnd) {
-    OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c(FortranOutputs[IndOutputs], "StateEnd")
-  }
-  ## DatesR and OutputsModel and StateEnd
-  if ((ExportDatesR & ExportStateEnd) | "all" %in% RunOptions$Outputs_Sim) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c("DatesR", FortranOutputs[IndOutputs], "StateEnd")
-  }
-
-  ## End
-  rm(RESULTS)
-  class(OutputsModel) <- c("OutputsModel", "daily", "GR")
-  return(OutputsModel)
-
+  ## OutputsModel generation
+  .GetOutputsModelGR(InputsModel,
+                     RunOptions,
+                     RESULTS,
+                     LInputSeries)
 }
diff --git a/R/RunModel_GR5H.R b/R/RunModel_GR5H.R
index 58ea098a765a1ee7195ad32b8fc7587172d22d8a..1d0139741fbcd7cda8bd567ee388ba7ef9805580 100644
--- a/R/RunModel_GR5H.R
+++ b/R/RunModel_GR5H.R
@@ -2,8 +2,6 @@ RunModel_GR5H <- function(InputsModel, RunOptions, Param) {
 
 
   ## Initialization of variables
-  NParam <- 5
-  FortranOutputs <- .FortranOutputs(GR = "GR5H")$GR
   IsIntStore <- inherits(RunOptions, "interception")
   if (IsIntStore) {
     Imax <- RunOptions$Imax
@@ -11,29 +9,8 @@ RunModel_GR5H <- function(InputsModel, RunOptions, Param) {
     Imax <- -99
   }
 
+  .ArgumentsCheckGR(InputsModel, RunOptions, Param)
 
-  ## Arguments check
-  if (!inherits(InputsModel, "InputsModel")) {
-    stop("'InputsModel' must be of class 'InputsModel'")
-  }
-  if (!inherits(InputsModel, "hourly")) {
-    stop("'InputsModel' must be of class 'hourly'")
-  }
-  if (!inherits(InputsModel, "GR")) {
-    stop("'InputsModel' must be of class 'GR'")
-  }
-  if (!inherits(RunOptions, "RunOptions")) {
-    stop("'RunOptions' must be of class 'RunOptions'")
-  }
-  if (!inherits(RunOptions, "GR")) {
-    stop("'RunOptions' must be of class 'GR'")
-  }
-  if (!is.vector(Param) | !is.numeric(Param)) {
-    stop("'Param' must be a numeric vector")
-  }
-  if (sum(!is.na(Param)) != NParam) {
-    stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
-  }
   Param <- as.double(Param)
 
   Param_X1X3_threshold <- 1e-2
@@ -58,9 +35,9 @@ RunModel_GR5H <- function(InputsModel, RunOptions, Param) {
   IndPeriod1   <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
   LInputSeries <- as.integer(length(IndPeriod1))
   if ("all" %in% RunOptions$Outputs_Sim) {
-    IndOutputs <- as.integer(1:length(FortranOutputs))
+    IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
   } else {
-    IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim)
+    IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
   }
 
   ## Output data preparation
@@ -107,38 +84,9 @@ RunModel_GR5H <- function(InputsModel, RunOptions, Param) {
                                         verbose = FALSE)
   }
 
-  ## Output data preparation
-  ## OutputsModel only
-  if (!ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
-    names(OutputsModel) <- FortranOutputs[IndOutputs]
-  }
-  ## DatesR and OutputsModel only
-  if (ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]))
-    names(OutputsModel) <- c("DatesR", FortranOutputs[IndOutputs])
-  }
-  ## OutputsModel and StateEnd only
-  if (!ExportDatesR & ExportStateEnd) {
-    OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c(FortranOutputs[IndOutputs], "StateEnd")
-  }
-  ## DatesR and OutputsModel and StateEnd
-  if ((ExportDatesR & ExportStateEnd) | "all" %in% RunOptions$Outputs_Sim) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c("DatesR", FortranOutputs[IndOutputs], "StateEnd")
-  }
-
-  ## End
-  rm(RESULTS)
-  class(OutputsModel) <- c("OutputsModel", "hourly", "GR")
-  if (IsIntStore) {
-    class(OutputsModel) <- c(class(OutputsModel), "interception")
-  }
-  return(OutputsModel)
-
+  ## OutputsModel generation
+  .GetOutputsModelGR(InputsModel,
+                     RunOptions,
+                     RESULTS,
+                     LInputSeries)
 }
diff --git a/R/RunModel_GR5J.R b/R/RunModel_GR5J.R
index fa5f83ac424e1df581cedcf2e02dcdb5145470df..73e17c6cb32116164131153a9e6decf7f49ce39b 100644
--- a/R/RunModel_GR5J.R
+++ b/R/RunModel_GR5J.R
@@ -1,35 +1,9 @@
 RunModel_GR5J <- function(InputsModel, RunOptions, Param) {
-  
-  
-  ## Initialization of variables
-  NParam <- 5
-  FortranOutputs <- .FortranOutputs(GR = "GR5J")$GR
-  
-  
-  ## Arguments check
-  if (!inherits(InputsModel, "InputsModel")) {
-    stop("'InputsModel' must be of class 'InputsModel'")
-  }  
-  if (!inherits(InputsModel, "daily")) {
-    stop("'InputsModel' must be of class 'daily'")
-  }  
-  if (!inherits(InputsModel, "GR")) {
-    stop("'InputsModel' must be of class 'GR'")
-  }  
-  if (!inherits(RunOptions, "RunOptions")) {
-    stop("'RunOptions' must be of class 'RunOptions'")
-  }  
-  if (!inherits(RunOptions, "GR")) {
-    stop("'RunOptions' must be of class 'GR'")
-  }  
-  if (!is.vector(Param) | !is.numeric(Param)) {
-    stop("'Param' must be a numeric vector")
-  }
-  if (sum(!is.na(Param)) != NParam) {
-    stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
-  }
+
+  .ArgumentsCheckGR(InputsModel, RunOptions, Param)
+
   Param <- as.double(Param)
-  
+
   Param_X1X3_threshold <- 1e-2
   Param_X4_threshold   <- 0.5
   if (Param[1L] < Param_X1X3_threshold) {
@@ -43,8 +17,8 @@ RunModel_GR5J <- function(InputsModel, RunOptions, Param) {
   if (Param[4L] < Param_X4_threshold) {
     warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
     Param[4L] <- Param_X4_threshold
-  }      
-  
+  }
+
   ## Input data preparation
   if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
     RunOptions$IndPeriod_WarmUp <- NULL
@@ -52,24 +26,24 @@ RunModel_GR5J <- function(InputsModel, RunOptions, Param) {
   IndPeriod1   <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
   LInputSeries <- as.integer(length(IndPeriod1))
   if ("all" %in% RunOptions$Outputs_Sim) {
-    IndOutputs <- as.integer(1:length(FortranOutputs)) 
+    IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
   } else {
-    IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim)
+    IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
   }
-  
+
   ## Output data preparation
   IndPeriod2     <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
   ExportDatesR   <- "DatesR"   %in% RunOptions$Outputs_Sim
   ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
-  
+
   ## Use of IniResLevels
   if (!is.null(RunOptions$IniResLevels)) {
     RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
     RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
   }
-  
+
   ## Call GR model Fortan
-  RESULTS <- .Fortran("frun_gr5j", PACKAGE = "airGR", 
+  RESULTS <- .Fortran("frun_gr5j", PACKAGE = "airGR",
                       ## inputs
                       LInputs = LInputSeries,                             ### length of input and output series
                       InputsPrecip = InputsModel$Precip[IndPeriod1],      ### input series of total precipitation [mm/d]
@@ -88,43 +62,17 @@ RunModel_GR5J <- function(InputsModel, RunOptions, Param) {
   RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
   if (ExportStateEnd) {
     RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
-    RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR5J, InputsModel = InputsModel, 
-                                        ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL, 
+    RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR5J, InputsModel = InputsModel,
+                                        ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
                                         UH1 = NULL,
-                                        UH2 = RESULTS$StateEnd[(1:40) + (7+20)], 
-                                        GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL, 
+                                        UH2 = RESULTS$StateEnd[(1:40) + (7+20)],
+                                        GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
                                         verbose = FALSE)
   }
-  
-  ## Output data preparation
-  ## OutputsModel only
-  if (!ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
-    names(OutputsModel) <- FortranOutputs[IndOutputs]
-  }
-  ## DatesR and OutputsModel only
-  if (ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), 
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]))
-    names(OutputsModel) <- c("DatesR", FortranOutputs[IndOutputs])
-  }
-  ## OutputsModel and StateEnd only
-  if (!ExportDatesR & ExportStateEnd) {
-    OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]), 
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c(FortranOutputs[IndOutputs], "StateEnd")
-  }
-  ## DatesR and OutputsModel and StateEnd
-  if ((ExportDatesR & ExportStateEnd) | "all" %in% RunOptions$Outputs_Sim) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), 
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]), 
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c("DatesR", FortranOutputs[IndOutputs], "StateEnd")
-  }
-  
-  ## End
-  rm(RESULTS) 
-  class(OutputsModel) <- c("OutputsModel", "daily", "GR")
-  return(OutputsModel)
-  
+
+  ## OutputsModel generation
+  .GetOutputsModelGR(InputsModel,
+                     RunOptions,
+                     RESULTS,
+                     LInputSeries)
 }
diff --git a/R/RunModel_GR6J.R b/R/RunModel_GR6J.R
index 572511f71610a75c262d34317ecfe6de8c4f7295..659a0b5cbece1c88aa883fc3a5d4247fccd0f98e 100644
--- a/R/RunModel_GR6J.R
+++ b/R/RunModel_GR6J.R
@@ -1,35 +1,9 @@
 RunModel_GR6J <- function(InputsModel, RunOptions, Param) {
-  
-  
-  ## Initialization of variables
-  NParam <- 6
-  FortranOutputs <- .FortranOutputs(GR = "GR6J")$GR
-  
-  
-  ## Arguments check
-  if (!inherits(InputsModel, "InputsModel")) {
-    stop("'InputsModel' must be of class 'InputsModel'")
-  }  
-  if (!inherits(InputsModel, "daily")) {
-    stop("'InputsModel' must be of class 'daily'")
-  }  
-  if (!inherits(InputsModel, "GR")) {
-    stop("'InputsModel' must be of class 'GR'")
-  }  
-  if (!inherits(RunOptions, "RunOptions")) {
-    stop("'RunOptions' must be of class 'RunOptions'")
-  }  
-  if (!inherits(RunOptions, "GR")) {
-    stop("'RunOptions' must be of class 'GR'")
-  }  
-  if (!is.vector(Param) | !is.numeric(Param)) {
-    stop("'Param' must be a numeric vector")
-  }
-  if (sum(!is.na(Param)) != NParam) {
-    stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
-  }
+
+  .ArgumentsCheckGR(InputsModel, RunOptions, Param)
+
   Param <- as.double(Param)
-  
+
   Param_X1X3X6_threshold <- 1e-2
   Param_X4_threshold     <- 0.5
   if (Param[1L] < Param_X1X3X6_threshold) {
@@ -47,8 +21,8 @@ RunModel_GR6J <- function(InputsModel, RunOptions, Param) {
   if (Param[6L] < Param_X1X3X6_threshold) {
     warning(sprintf("Param[6] (X6: coefficient for emptying exponential store [mm]) < %.2f\n X6 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
     Param[6L] <- Param_X1X3X6_threshold
-  }      
-  
+  }
+
   ## Input data preparation
   if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
     RunOptions$IndPeriod_WarmUp <- NULL
@@ -56,25 +30,25 @@ RunModel_GR6J <- function(InputsModel, RunOptions, Param) {
   IndPeriod1   <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
   LInputSeries <- as.integer(length(IndPeriod1))
   if ("all" %in% RunOptions$Outputs_Sim) {
-    IndOutputs <- as.integer(1:length(FortranOutputs)) 
+    IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
   } else {
-    IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim)
+    IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
   }
-  
+
   ## Output data preparation
   IndPeriod2     <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
   ExportDatesR   <- "DatesR"   %in% RunOptions$Outputs_Sim
   ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
-  
+
   ## Use of IniResLevels
   if (!is.null(RunOptions$IniResLevels)) {
     RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
     RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
     RunOptions$IniStates[3] <- RunOptions$IniResLevels[3]          ### exponential store level (mm)
   }
-  
+
   ## Call GR model Fortan
-  RESULTS <- .Fortran("frun_gr6j", PACKAGE = "airGR", 
+  RESULTS <- .Fortran("frun_gr6j", PACKAGE = "airGR",
                       ## inputs
                       LInputs = LInputSeries,                             ### length of input and output series
                       InputsPrecip = InputsModel$Precip[IndPeriod1],      ### input series of total precipitation [mm/d]
@@ -93,43 +67,17 @@ RunModel_GR6J <- function(InputsModel, RunOptions, Param) {
   RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
   if (ExportStateEnd) {
     RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
-    RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR6J, InputsModel = InputsModel, 
-                                        ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = RESULTS$StateEnd[3L], 
+    RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR6J, InputsModel = InputsModel,
+                                        ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = RESULTS$StateEnd[3L],
                                         UH1 = RESULTS$StateEnd[(1:20) + 7],
-                                        UH2 = RESULTS$StateEnd[(1:40) + (7+20)], 
-                                        GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL, 
+                                        UH2 = RESULTS$StateEnd[(1:40) + (7+20)],
+                                        GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
                                         verbose = FALSE)
   }
-  
-  ## Output data preparation
-  ## OutputsModel only
-  if (!ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
-    names(OutputsModel) <- FortranOutputs[IndOutputs]
-  }
-  ## DatesR and OutputsModel only
-  if (ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), 
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]))
-    names(OutputsModel) <- c("DatesR", FortranOutputs[IndOutputs])
-  }
-  ## OutputsModel and StateEnd only
-  if (!ExportDatesR & ExportStateEnd) {
-    OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]), 
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c(FortranOutputs[IndOutputs], "StateEnd")
-  }
-  ## DatesR and OutputsModel and StateEnd
-  if ((ExportDatesR & ExportStateEnd) | "all" %in% RunOptions$Outputs_Sim) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), 
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]), 
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c("DatesR", FortranOutputs[IndOutputs], "StateEnd")
-  }
-  
-  ## End
-  rm(RESULTS) 
-  class(OutputsModel) <- c("OutputsModel", "daily", "GR")
-  return(OutputsModel)
-  
+
+  ## OutputsModel generation
+  .GetOutputsModelGR(InputsModel,
+                     RunOptions,
+                     RESULTS,
+                     LInputSeries)
 }
diff --git a/R/RunModel_Lag.R b/R/RunModel_Lag.R
index 473e7758993b58d2c2317e29279b76424e2b25f8..7889967059b1646cc2e8c1d368897f8f3a28fdad 100644
--- a/R/RunModel_Lag.R
+++ b/R/RunModel_Lag.R
@@ -69,9 +69,21 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
       IniSD[iStart:(iStart + PT[x])]
     })
   } else {
-    IniStates <- lapply(seq_len(NbUpBasins), function(x) {
-      rep(0, floor(PT[x] + 1))
-    })
+    IniStates <- lapply(
+      seq_len(NbUpBasins),
+      function(iUpBasins) {
+        iWarmUp <- seq.int(
+          from = max(1, RunOptions$IndPeriod_WarmUp[length(RunOptions$IndPeriod_WarmUp)] - floor(PT[iUpBasins])),
+          to   = max(1, RunOptions$IndPeriod_WarmUp[length(RunOptions$IndPeriod_WarmUp)])
+        )
+        ini <- InputsModel$Qupstream[iWarmUp, iUpBasins]
+        if (length(ini) != floor(PT[iUpBasins] + 1)) {
+          # If warm-up period is not enough long complete beginning with first value
+          ini <- c(rep(ini[1], floor(PT[iUpBasins] + 1) - length(ini)), ini)
+        }
+        return(as.vector(ini))
+      }
+    )
   }
   # message("Initstates: ", paste(IniStates, collapse = ", "))
 
@@ -88,7 +100,7 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
   # message("Qsim: ", paste(OutputsModel$Qsim, collapse = ", "))
 
   # Warning for negative flows or NAs only in extended outputs
-  if(length(RunOptions$Outputs_Sim) > 2) {
+  if (length(RunOptions$Outputs_Sim) > 2) {
     if (any(OutputsModel$Qsim[!is.na(OutputsModel$Qsim)] < 0)) {
       warning(length(which(OutputsModel$Qsim < 0)), " time steps with negative flow, set to zero.")
       OutputsModel$Qsim[OutputsModel$Qsim < 0] <- 0
diff --git a/R/SeriesAggreg.OutputsModel.R b/R/SeriesAggreg.OutputsModel.R
index 3bf0d8aefd24f930c4d6c0ba36fd9786d11e640f..d065182bab683414545511a75f78aa09e6f7c58d 100644
--- a/R/SeriesAggreg.OutputsModel.R
+++ b/R/SeriesAggreg.OutputsModel.R
@@ -2,6 +2,6 @@ SeriesAggreg.OutputsModel <- function(x, Format, ...) {
   SeriesAggreg.list(x,
                     Format,
                     ConvertFun = NA,
-                    except = "StateEnd",
+                    except = c("WarmUpQsim", "StateEnd"),
                     ...)
 }
diff --git a/R/SeriesAggreg.data.frame.R b/R/SeriesAggreg.data.frame.R
index f963a0540325f30772829a6ad1c67a0aa9fbcd9d..5a467039d1d6d2ee4da3dd99be7b6a8e5d236757 100644
--- a/R/SeriesAggreg.data.frame.R
+++ b/R/SeriesAggreg.data.frame.R
@@ -60,10 +60,10 @@ SeriesAggreg.data.frame <- function(x,
   lapply(ConvertFun, function(y) {
     if (!grepl("^q\\d+$", y, ignore.case = TRUE)) {
       TestOutput <- listConvertFun[[y]](1:10)
-      if(!is.numeric(TestOutput)) {
+      if (!is.numeric(TestOutput)) {
         stop(sprintf("Returned value of '%s' function should be numeric", y))
       }
-      if(length(TestOutput) != 1) {
+      if (length(TestOutput) != 1) {
         stop(sprintf("Returned value of '%s' function should be of length 1", y))
       }
     }
diff --git a/R/Utils.R b/R/Utils.R
index 82c66da0d6410fa7e5747985c73fc44f860836a7..d5e14535fc3839bd4564168b99a56f40d24a39b2 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -3,7 +3,7 @@
 ## function to check
 ## =================================================================================
 
-# .onLoad <- function(libname, pkgname){
+# .onLoad <- function(libname, pkgname) {
 #   if (requireNamespace("airGRteaching", quietly = TRUE)) {
 #     if (packageVersion("airGRteaching") %in% package_version(c("0.2.0.9", "0.2.2.2", "0.2.3.2"))) {
 #       packageStartupMessage("In order to be compatible with the present version of 'airGR', please update your version of the 'airGRteaching' package.")
@@ -71,6 +71,7 @@
         stop("the time step of the model inputs must be ", res$TimeUnit)
       }
     }
+    res$CodeModHydro <- gsub("CemaNeige", "", res$CodeMod)
     return(res)
   }
 }
diff --git a/R/UtilsRunModel.R b/R/UtilsRunModel.R
new file mode 100644
index 0000000000000000000000000000000000000000..096aff84a47a57e68253a14b75b3d3a9b56115e3
--- /dev/null
+++ b/R/UtilsRunModel.R
@@ -0,0 +1,100 @@
+#' Create `OutputsModel` for GR non-Cemaneige models
+#'
+#' @param InputsModel output of [CreateInputsModel]
+#' @param RunOptions output of [CreateRunOptions]
+#' @param RESULTS outputs of [.Fortran]
+#' @param LInputSeries number of time steps of warm-up + run periods
+#' @param CemaNeigeLayers outputs of Cemaneige pre-process
+#'
+#' @return OutputsModel object
+#' @noRd
+#'
+.GetOutputsModelGR <- function(InputsModel,
+                               RunOptions,
+                               RESULTS,
+                               LInputSeries,
+                               CemaNeigeLayers = NULL) {
+
+  IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
+  FortranOutputs <- RunOptions$FortranOutputs$GR
+
+  IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim)
+
+  OutputsModel <- list()
+
+  if ("DatesR" %in% RunOptions$Outputs_Sim) {
+    OutputsModel$DatesR <- InputsModel$DatesR[RunOptions$IndPeriod_Run]
+  }
+
+  seqOutputs <- seq_len(RESULTS$NOutputs)
+  names(seqOutputs) <- FortranOutputs[IndOutputs]
+
+  OutputsModel <- c(OutputsModel,
+                    lapply(seqOutputs, function(i) RESULTS$Outputs[IndPeriod2, i]))
+
+  if (!is.null(CemaNeigeLayers)) {
+    OutputsModel$CemaNeigeLayers <- CemaNeigeLayers
+  }
+
+  if ("WarmUpQsim" %in% RunOptions$Outputs_Sim) {
+    OutputsModel$WarmUpQsim <- RESULTS$Outputs[seq_len(length(RunOptions$IndPeriod_WarmUp)),
+                                               which(FortranOutputs == "Qsim")]
+    class(OutputsModel$WarmUpQsim) <- c("WarmUpOutputsModelItem", class(OutputsModel$WarmUpQsim))
+  }
+
+  if ("StateEnd" %in% RunOptions$Outputs_Sim) {
+    OutputsModel$StateEnd <- RESULTS$StateEnd
+  }
+
+  class(OutputsModel) <- c("OutputsModel", class(RunOptions)[-1])
+
+  return(OutputsModel)
+}
+
+
+#' Check arguments of `RunModel_*GR*` functions
+#'
+#' @param InputsModel see [CreateInputsModel]
+#' @param RunOptions  see [CreateRunOptions]
+#' @param Param [numeric] [vector] model calibration parameters
+#'
+#' @return [NULL]
+#' @noRd
+#'
+.ArgumentsCheckGR <- function(InputsModel, RunOptions, Param) {
+  if (!inherits(InputsModel, "InputsModel")) {
+    stop("'InputsModel' must be of class 'InputsModel'")
+  }
+  if (!inherits(InputsModel, RunOptions$FeatFUN_MOD$TimeUnit)) {
+    stop("'InputsModel' must be of class '", RunOptions$FeatFUN_MOD$TimeUnit, "'")
+  }
+  if (!inherits(InputsModel, "GR")) {
+    stop("'InputsModel' must be of class 'GR'")
+  }
+  if (class(RunOptions)[1] != "RunOptions") {
+    if (!inherits(RunOptions, "RunOptions")) {
+      stop("'RunOptions' must be of class 'RunOptions'")
+    } else {
+      stop("'RunOptions' class of 'RunOptions' must be in first position")
+    }
+  }
+  if (!inherits(RunOptions, "GR")) {
+    stop("'RunOptions' must be of class 'GR'")
+  }
+
+  if ("CemaNeige" %in% RunOptions$FeatFUN_MOD$Class) {
+    if (!inherits(InputsModel, "CemaNeige")) {
+      stop("'InputsModel' must be of class 'CemaNeige'")
+    }
+    if (!inherits(RunOptions, "CemaNeige")) {
+      stop("'RunOptions' must be of class 'CemaNeige'")
+    }
+  }
+
+  if (!is.vector(Param) | !is.numeric(Param)) {
+    stop("'Param' must be a numeric vector")
+  }
+  if (sum(!is.na(Param)) != RunOptions$FeatFUN_MOD$NbParam) {
+    stop(paste("'Param' must be a vector of length", RunOptions$FeatFUN_MOD$NbParam, "and contain no NA"))
+  }
+}
diff --git a/R/UtilsSeriesAggreg.R b/R/UtilsSeriesAggreg.R
index 9b47efac0bdc70bd1cabbf3cf6db3e6ed1c0603e..d8dff978a383ded8746fb438ca4276720c5901f1 100644
--- a/R/UtilsSeriesAggreg.R
+++ b/R/UtilsSeriesAggreg.R
@@ -53,7 +53,7 @@
     iRes <- AggregConvertFunTable$ConvertFun[AggregConvertFunTable$x == iX]
     iRes <- ifelse(test = any(is.na(iRes)), yes = NA, no = iRes) # R < 4.0 compatibility
   })
-  if(Format %in% c("%d", "%m")) {
+  if (Format %in% c("%d", "%m")) {
     res <- rep("mean", length(res))
   }
   return(res)
diff --git a/data-raw/vignettes.R b/data-raw/vignettes.R
new file mode 100644
index 0000000000000000000000000000000000000000..618fb37da8b32c31d036b6dd385f324cfb51766e
--- /dev/null
+++ b/data-raw/vignettes.R
@@ -0,0 +1,25 @@
+## code to prepare datasets used in vignettes
+library(airGR)
+
+source("tests/testthat/helper_vignettes.R")
+
+# V02.1_param_optim
+RunVignetteChunks("V02.1_param_optim")
+save(resGLOB, resPORT,
+     file = "inst/vignettesData/vignetteParamOptim.rda",
+     version = 2)
+save(algo, Ind_Run, InputsCrit_inv, InputsModel, MOptimGR4J, optMO, RunOptions,
+     file = "inst/vignettesData/vignetteParamOptimCaramel.rda",
+     version = 2)
+
+# V02.2_param_mcmc
+RunVignetteChunks("V02.2_param_mcmc")
+save(gelRub, multDRAM,
+     file = "inst/vignettesData/vignetteParamMCMC.rda",
+     version = 2)
+
+# V04_cemaneige_hysteresis
+RunVignetteChunks("V04_cemaneige_hysteresis")
+save(OutputsCrit_Cal, OutputsCrit_Cal_NoHyst, OutputsCrit_Val, OutputsCrit_Val_NoHyst,
+     file = "inst/vignettesData/vignetteParamMCMC.rda",
+     version = 2)
diff --git a/inst/vignettesData/vignetteParamOptimCaramel.rda b/inst/vignettesData/vignetteParamOptimCaramel.rda
index 7d2e9fcc41220cab22fa5b1e57e7c666ef74803e..adf9461d0a5e0e4dd70f5b55034920aaa28e84c2 100644
Binary files a/inst/vignettesData/vignetteParamOptimCaramel.rda and b/inst/vignettesData/vignetteParamOptimCaramel.rda differ
diff --git a/man/CreateCalibOptions.Rd b/man/CreateCalibOptions.Rd
index c98ad748b70f67f09ada0082c21502b5358963ac..22e16bba08e7c1053f09a19bd4efeb52186f8570 100644
--- a/man/CreateCalibOptions.Rd
+++ b/man/CreateCalibOptions.Rd
@@ -5,7 +5,7 @@
 \alias{CreateCalibOptions}
 
 
-\title{Creation of the CalibOptions object required but the Calibration functions}
+\title{Creation of the CalibOptions object required but the Calibration* functions}
 
 
 \description{
diff --git a/man/SeriesAggreg.Rd b/man/SeriesAggreg.Rd
index 48d5da1c4acf52f7cbd067fa0fb7bff37e28477f..cbd43f901794561ba690f863003230a216401980 100644
--- a/man/SeriesAggreg.Rd
+++ b/man/SeriesAggreg.Rd
@@ -30,6 +30,8 @@ Warning: on the aggregated outputs, the dates correspond to the beginning of the
   E.g. use "Q90" for calculating 90th percentile in the aggregation.
   The formula used is: \code{quantile(x, probs = perc / 100, type = 8, na.rm = TRUE)}.
 
+  As there are multiple ways to take into account missing values in aggregation functions, \code{NA}s are not supported by \code{SeriesAggreg} and it provides \code{NA} values when \code{NA}s are present in the \code{x} input.
+
 }
 
 
diff --git a/tests/testthat/helper_regression.R b/tests/scheduled_tests/regression.R
similarity index 72%
rename from tests/testthat/helper_regression.R
rename to tests/scheduled_tests/regression.R
index bc05c5521a19946aa2cef6170f7f28f830dfefe2..88c7b0f6511f99953bb07a56d497735388b53d37 100644
--- a/tests/testthat/helper_regression.R
+++ b/tests/scheduled_tests/regression.R
@@ -1,3 +1,5 @@
+# Helper functions for regression
+
 StoreStableExampleResults <- function(
   package = "airGR",
   path = file.path("tests/tmp", Sys.getenv("R_VERSION"), "stable"),
@@ -77,7 +79,31 @@ StoreTopicResults <- function(topic, package, path, run.dontrun = TRUE, run.dont
 CompareStableDev <- function() {
   res <- testthat::test_file("tests/testthat/regression.R")
   dRes <- as.data.frame(res)
-  if(any(dRes[, "failed"] > 0) | any(dRes[, "error"])) {
+  if (any(dRes[, "failed"] > 0) | any(dRes[, "error"])) {
     quit(status = 1)
   }
 }
+
+###############
+# MAIN SCRIPT #
+###############
+
+# Execute Regression test by comparing RD files stored in folders /tests/tmp/ref and /tests/tmp/test
+Args <- commandArgs(trailingOnly = TRUE)
+
+lActions <- list(
+  stable = StoreStableExampleResults,
+  dev = StoreDevExampleResults,
+  compare = CompareStableDev
+)
+
+if (length(Args) == 1 && Args %in% names(lActions)) {
+  lActions[[Args]]()
+} else {
+  stop("This script should be run with one argument in the command line:\n",
+       "`Rscript tests/regression_tests.R [stable|dev|compare]`.\n",
+       "Available arguments are:\n",
+       "- stable: install stable version from CRAN, run and store examples\n",
+       "- dev: install dev version from current directory, run and store examples\n",
+       "- compare: stored results of both versions")
+}
diff --git a/tests/scheduled_tests/scheduled.R b/tests/scheduled_tests/scheduled.R
new file mode 100644
index 0000000000000000000000000000000000000000..ece2a2576562b87566add4ad697463c88f86eae7
--- /dev/null
+++ b/tests/scheduled_tests/scheduled.R
@@ -0,0 +1,46 @@
+#' Script for running scheduled test
+#'
+#' All files with the pattern /testthat/tests/scheduled-*.R are tested
+#' as testthat does for files /testthat/tests/test-*.R.
+#'
+#' This script should be started with `source` command from the root of the package.
+#' @example
+#' source("tests/scheduled.R")
+
+####################
+# Helper functions #
+####################
+
+#' Wrapper for [quit] which is only applied outside of RStudio
+#'
+#' @param status See `status` parameter of [quit]. Default `quit = 1`.
+#' @param ... Other parameters sent to [quit]
+#'
+#' @return NULL
+#' @export
+quit2 <- function(status = 1, ...) {
+  if (all(!grepl("rstudio", Sys.getenv(), ignore.case = TRUE))) {
+    quit(status, ...)
+  }
+}
+
+###############
+# MAIN SCRIPT #
+###############
+
+library(testthat)
+library(airGR)
+
+scheduled_tests <- list.files(
+  path = "tests/testthat",
+  pattern = "^scheduled-.*\\.R$",
+  full.names = TRUE
+)
+
+lRes <- lapply(scheduled_tests, test_file)
+for (res in lRes) {
+  dRes <- as.data.frame(res)
+  if (any(dRes[, "failed"] > 0) | any(dRes[, "error"])) {
+    quit2()
+  }
+}
diff --git a/tests/testthat/helper_vignettes.R b/tests/testthat/helper_vignettes.R
index 931f99a55ee90fa4417040a36b58f66220aaac71..98951fef380002bc1243e705ceaf2ac5c5d124a1 100644
--- a/tests/testthat/helper_vignettes.R
+++ b/tests/testthat/helper_vignettes.R
@@ -70,9 +70,12 @@ RunRmdChunks <- function(fileRmd,
 RunVignetteChunks <- function(vignette,
                               tmpFolder = "../tmp",
                               force.eval = TRUE) {
-  if(file.exists(sprintf("../../vignettes/%s.Rmd", vignette))) {
+  if (file.exists(sprintf("../../vignettes/%s.Rmd", vignette))) {
     # testthat context in development environnement
     RunRmdChunks(sprintf("../../vignettes/%s.Rmd", vignette), tmpFolder, force.eval)
+  } else if (file.exists(sprintf("vignettes/%s.Rmd", vignette))) {
+    # context in direct run in development environnement
+    RunRmdChunks(sprintf("vignettes/%s.Rmd", vignette), tmpFolder, force.eval)
   } else {
     # R CMD check context in package environnement
     RunRmdChunks(system.file(sprintf("doc/%s.Rmd", vignette), package = "airGR"), tmpFolder, force.eval)
diff --git a/tests/testthat/regression_tests.R b/tests/testthat/regression_tests.R
deleted file mode 100644
index 4456e09edc85aaa6b06bd36f3b518d41c192e786..0000000000000000000000000000000000000000
--- a/tests/testthat/regression_tests.R
+++ /dev/null
@@ -1,21 +0,0 @@
-# Execute Regression test by comparing RD files stored in folders /tests/tmp/ref and /tests/tmp/test
-Args <- commandArgs(trailingOnly = TRUE)
-
-source("tests/testthat/helper_regression.R")
-
-lActions <- list(
-  stable = StoreStableExampleResults,
-  dev = StoreDevExampleResults,
-  compare = CompareStableDev
-)
-
-if(Args %in% names(lActions)) {
-  lActions[[Args]]()
-} else {
-  stop("This script should be run with one argument in the command line:\n",
-       "`Rscript tests/regression_tests.R [stable|dev|compare]`.\n",
-       "Available arguments are:\n",
-       "- stable: install stable version from CRAN, run and store examples\n",
-       "- dev: install dev version from current directory, run and store examples\n",
-       "- compare: stored results of both versions")
-}
diff --git a/tests/testthat/scheduled-Calibration.R b/tests/testthat/scheduled-Calibration.R
new file mode 100644
index 0000000000000000000000000000000000000000..6b6b98b6a0426913d2b54c49283d2339c9b44f42
--- /dev/null
+++ b/tests/testthat/scheduled-Calibration.R
@@ -0,0 +1,141 @@
+context("Calibration")
+
+sModels <- c(
+  "name          IsHyst data     aggreg ParamFinalR",
+  "GR1A          FALSE  L0123001 %Y     0.91125",
+  "GR2M          FALSE  L0123001 %Y%m   259.8228;0.9975",
+  "GR4J          FALSE  L0123001 NA     223.6315877;0.5781516;97.5143942;2.2177177",
+  "GR5J          FALSE  L0123001 NA     220.3863609;0.8944531;93.5640705;1.7628720;0.4846427",
+  "GR6J          FALSE  L0123001 NA     192.8761657;0.6933087;49.1783293;2.2145422;0.5088240;6.8146261",
+  "CemaNeigeGR4J FALSE  L0123001 NA     2.043839e+02;5.781516e-01;1.025141e+02;2.217718e+00;1.501502e-03;1.432036e+01",
+  "CemaNeigeGR5J FALSE  L0123001 NA     1.983434e+02;8.747758e-01;9.849443e+01;1.768769e+00;4.829830e-01;1.501502e-03;1.432036e+01",
+  "CemaNeigeGR6J FALSE  L0123001 NA     184.9341841;0.5551637;59.7398917;2.2177177;0.4760000;6.0496475;0.0000000;14.4642868",
+  "CemaNeigeGR4J TRUE   L0123001 NA     208.5127103;0.5781516;102.5140641;2.2274775;0.0000000;6.7644613;8.0000000;1.0000000",
+  "CemaNeigeGR5J TRUE   L0123001 NA     202.350228;0.901525;98.494430;1.788288;0.483984;0.000000;7.401500;6.100000;1.000000",
+  "CemaNeigeGR6J TRUE   L0123001 NA     188.67010241;0.56662930;60.34028760;2.22747748;0.47600000;5.98945247;0.03203203;7.93816892;10.80000000;1.00000000",
+  "GR4H            FALSE  L0123003 NA     711.676649;-1.158469;150.556095;4.686093",
+  "GR5H            FALSE  L0123003 NA     804.0021672;-0.1898488;137.7524699;3.0436628;0.1951163",
+  "CemaNeigeGR4H   FALSE  L0123003 NA     1.595685e+03;-8.183484e-01;2.320697e+02;5.000000e-01;5.005005e-04;9.342369e+01",
+  "CemaNeigeGR5H   FALSE  L0123003 NA     33.34921883;-4.98925432;332.00673122;1.58534106;0.20792716;0.02214393;4.28498513",
+  "CemaNeigeGR4H   TRUE   L0123003 NA     1.766316e+03;-6.920667e-01;2.192034e+02;3.451688e+00;5.005005e-04;4.869585e+01;1.111447e+01;5.064090e-01",
+  "CemaNeigeGR5H   TRUE   L0123003 NA     66.6863310;-1.4558128;138.3795123;2.6499450;0.2325000;0.0000000;0.3017014;48.4000000;0.9914915"
+)
+dfModels <- read.table(text = paste(sModels, collapse = "\n"), header = TRUE)
+
+ModelCalibration <- function(model) {
+  sModel <- paste0("RunModel_", model$name)
+  sIM_FUN_MOD <- sModel
+
+  if (model$data == "L0123003") {
+    # hourly time step database
+    dates <- c("2004-01-01 00:00", "2004-12-31 23:00", "2005-01-01 00:00", "2008-12-31 23:00")
+    date_format = "%Y-%m-%d %H:%M"
+    TempMean <- fakeHourlyTemp()
+  } else {
+    # yearly, monthly, daily time step databases
+    dates <- c("1985-01-01", "1985-12-31", "1986-01-01", "2012-12-31")
+    date_format <- "%Y-%m-%d"
+    if (!is.na(model$aggreg)) {
+      # Aggregation on monthly and yearly databases
+      sIM_FUN_MOD <- "RunModel_GR4J" # CreateInputsModel with daily data
+      date_format <- model$aggreg
+    }
+  }
+
+  ## loading catchment data
+  data(list = model$data)
+  if (model$data != "L0123003") TempMean <- BasinObs$T
+
+  # preparation of the InputsModel object
+  InputsModel <- CreateInputsModel(FUN_MOD = sIM_FUN_MOD,
+                                   DatesR = BasinObs$DatesR,
+                                   Precip = BasinObs$P,
+                                   PotEvap = BasinObs$E,
+                                   TempMean = TempMean,
+                                   ZInputs = median(BasinInfo$HypsoData),
+                                   HypsoData = BasinInfo$HypsoData,
+                                   NLayers = 5)
+
+  if (!is.na(model$aggreg)) {
+    # conversion of InputsModel to target time step
+    InputsModel <- SeriesAggreg(InputsModel, Format = model$aggreg)
+
+    dfQobs <- SeriesAggreg(data.frame(DatesR = BasinObs$DatesR, Qmm = BasinObs$Qmm),
+                           Format = model$aggreg, ConvertFun = "sum")
+    Obs <- dfQobs$Qmm
+  } else {
+    Obs <- BasinObs$Qmm
+  }
+
+  # calibration period selection
+  dates <- sapply(dates, function(x) format(as.Date(x), format = date_format))
+  Ind_WarmUp <- seq(
+    which(format(InputsModel$DatesR, format = date_format)==dates[1]),
+    which(format(InputsModel$DatesR, format = date_format)==dates[2])
+  )
+  Ind_Run <- seq(
+    which(format(InputsModel$DatesR, format = date_format)==dates[3]),
+    which(format(InputsModel$DatesR, format = date_format)==dates[4])
+  )
+
+  # preparation of the RunOptions object
+  suppressWarnings(
+    RunOptions <- CreateRunOptions(
+      FUN_MOD = sModel,
+      InputsModel = InputsModel,
+      IndPeriod_Run = Ind_Run,
+      IndPeriod_WarmUp = Ind_WarmUp,
+      IsHyst = as.logical(model$IsHyst)
+    )
+  )
+
+  # calibration criterion: preparation of the InputsCrit object
+  InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
+                                 RunOptions = RunOptions, Obs = Obs[Ind_Run])
+  # preparation of CalibOptions object
+  CalibOptions <- CreateCalibOptions(sModel, IsHyst = as.logical(model$IsHyst))
+
+  # calibration
+  suppressWarnings(OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions,
+                                               InputsCrit = InputsCrit, CalibOptions = CalibOptions,
+                                               FUN_MOD = sModel))
+  OutputsCalib$ParamFinalR
+}
+
+
+TestModelCalibration <- function(model) {
+  model <- as.list(model)
+
+  test_that(paste(model$name, ifelse(as.logical(model$IsHyst), "Hysteresis", ""), "works"), {
+
+    ParamFinalR <- ModelCalibration(model)
+
+    expect_equal(ParamFinalR,
+                 as.numeric(strsplit(model$ParamFinalR, ";")[[1]]),
+                 tolerance = 1E-6)
+  })
+}
+
+
+#' Create Fake hourly temperature from daily temperatures in L0123001
+#'
+#' @param start_date [character] start date in format "%Y-%m-%d"
+#' @param end_date [character] end date in format "%Y-%m-%d"
+#' @return [numeric] hourly temperature time series between `start_date` and `end_date`
+fakeHourlyTemp <- function(start_date = "2004-01-01", end_date = "2008-12-31") {
+  dates <- as.POSIXct(c(start_date, end_date), tz = "UTC")
+  data(L0123002)
+  indJ <- seq.int(which(BasinObs$DatesR == as.POSIXct(dates[1])),
+                  which(BasinObs$DatesR == as.POSIXct(dates[2])))
+  TJ <- BasinObs$T[indJ]
+
+  TH <- approx((seq.int(length(TJ)) - 1) * 24,TJ,
+               seq.int(length(TJ) * 24 ) - 1,
+               rule = 2)$y
+  varT_1J <- -sin(0:23/24 * 2 * pi) # Temp min at 6 and max at 18
+  varT <- rep(varT_1J, length(TJ))
+  TH <- TH + varT * 5 # For a mean daily amplitude of 10°
+  TH
+}
+
+apply(dfModels, 1, TestModelCalibration)
diff --git a/tests/testthat/test-RunModel_Lag.R b/tests/testthat/test-RunModel_Lag.R
index c7cc53ec8c72d8b6eba1648c5977ac6e535b3ee1..81537b6bda92f61c04007ef532d60207141d7f07 100644
--- a/tests/testthat/test-RunModel_Lag.R
+++ b/tests/testthat/test-RunModel_Lag.R
@@ -165,7 +165,7 @@ Qm3GR4Only <- OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2L] * 1e3
 test_that("1 input with lag of 1 time step delay out gives an output delayed of one time step", {
   OutputsSD <- RunModel(InputsModel, RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
   Qm3SdSim <- OutputsSD$Qsim_m3
-  Qm3UpstLagObs <- c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]) * InputsModel$BasinAreas[1] * 1e3
+  Qm3UpstLagObs <- Qupstream[Ind_Run - 1] * InputsModel$BasinAreas[1] * 1e3
   expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
 })
 
@@ -174,14 +174,14 @@ test_that("1 input with lag of 0.5 time step delay out gives an output delayed o
                         Param = c(InputsModel$LengthHydro * 1e3 / (12 * 3600), Param),
                         FUN_MOD = RunModel_GR4J)
   Qm3SdSim <- OutputsSD$Qsim_m3
-  Qm3UpstLagObs <- (Qupstream[Ind_Run] + c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]])) / 2 * InputsModel$BasinAreas[1] * 1e3
+  Qm3UpstLagObs <- (Qupstream[Ind_Run] + Qupstream[Ind_Run - 1]) / 2 * InputsModel$BasinAreas[1] * 1e3
   expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
 })
 
 test_that("Qupstream in different units should return the same result", {
   OutputsSD_mm <- RunModel(InputsModel, RunOptions,
-                        Param = ParamSD,
-                        FUN_MOD = RunModel_GR4J)
+                           Param = ParamSD,
+                           FUN_MOD = RunModel_GR4J)
   InputsModel_m3 <- CreateInputsModel(
     FUN_MOD = RunModel_GR4J,
     DatesR = BasinObs$DatesR,
@@ -208,8 +208,8 @@ test_that("Qupstream in different units should return the same result", {
     QupstrUnit = "m3/s"
   )
   OutputsSD_m3s <- RunModel(InputsModel_m3s, RunOptions,
-                           Param = ParamSD,
-                           FUN_MOD = RunModel_GR4J)
+                            Param = ParamSD,
+                            FUN_MOD = RunModel_GR4J)
   expect_equal(OutputsSD_mm$Qsim, OutputsSD_m3s$Qsim)
 
   InputsModel_ls <- CreateInputsModel(
@@ -233,7 +233,7 @@ InputsCrit <- CreateInputsCrit(
   InputsModel = InputsModel,
   RunOptions = RunOptions,
   VarObs = "Q",
-  Obs = (c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]) * BasinAreas[1L] +
+  Obs = (Qupstream[Ind_Run - 1] * BasinAreas[1L] +
            BasinObs$Qmm[Ind_Run] * BasinAreas[2L]) / sum(BasinAreas)
 )
 
@@ -283,7 +283,7 @@ test_that("1 no area input with lag of 1 time step delay out gives an output del
   expect_false(any(is.na(OutputsSD$Qsim)))
 
   Qm3SdSim <- OutputsSD$Qsim_m3
-  Qm3UpstLagObs <- c(0, InputsModel$Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]])
+  Qm3UpstLagObs <- InputsModel$Qupstream[Ind_Run - 1]
 
   expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
 })
@@ -338,6 +338,21 @@ test_that("Error on Wrong length of iniState$SD", {
                                   InputsModel = IM, IndPeriod_Run = Ind_Run2,
                                   IndPeriod_WarmUp = 0L,
                                   IniStates = OutputsModel1$StateEnd)
-  expect_error(RunModel(InputsModel = IM, RunOptions = RunOptions2, Param = PSDini, FUN_MOD = RunModel_GR4J)
+  expect_error(
+    RunModel(InputsModel = IM, RunOptions = RunOptions2, Param = PSDini, FUN_MOD = RunModel_GR4J)
   )
 })
+
+test_that("First Qupstream time steps must be repeated if warm-up period is too short", {
+  IM2 <- IM[2558:3557]
+  IM2$BasinAreas[3] <- 0
+  IM2$Qupstream <- matrix(rep(1:1000, 2), ncol = 2)
+  RunOptions2 <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
+                                  InputsModel = IM2, IndPeriod_Run = seq_len(1000),
+                                  IndPeriod_WarmUp = 0L)
+  OM2 <- RunModel(InputsModel = IM2,
+                  RunOptions = RunOptions2,
+                  Param = PSDini,
+                  FUN_MOD = RunModel_GR4J)
+  expect_equal(OM2$Qsim_m3[1:3], rep(2,3))
+})
diff --git a/vignettes/V02.1_param_optim.Rmd b/vignettes/V02.1_param_optim.Rmd
index bcd3f3e8d21ef017391aaac5f94c8a61650678bc..7ded694184a5c890656f25e7c90b3ca044d9e145 100644
--- a/vignettes/V02.1_param_optim.Rmd
+++ b/vignettes/V02.1_param_optim.Rmd
@@ -269,7 +269,7 @@ ggplot() +
   theme_bw()
 ```
 
-The pameter sets can be viewed in the parameter space, illustrating different populations.
+The parameter sets can be viewed in the parameter space, illustrating different populations.
 
 ```{r fig.height=6, fig.width=6, message=FALSE, warning=FALSE}
 param_optMO <- apply(optMO$parameters, MARGIN = 1, FUN = function(x) {
diff --git a/vignettes/V05_sd_model.Rmd b/vignettes/V05_sd_model.Rmd
index 632819c13cdd8c32d751f23f4933467a030de2d8..800db1ee9ee77a6487f07122ba15040da61ca3aa 100644
--- a/vignettes/V05_sd_model.Rmd
+++ b/vignettes/V05_sd_model.Rmd
@@ -1,11 +1,11 @@
 ---
-title: "Simulating a reservoir with semi-distributed GR4J model"
+title: "Simulated vs observed upstream flows in calibration of semi-distributed GR4J model"
 author: "David Dorchies"
 bibliography: V00_airgr_ref.bib
 output: rmarkdown::html_vignette
 vignette: >
   %\VignetteEngine{knitr::rmarkdown}
-  %\VignetteIndexEntry{Simulating a reservoir with semi-distributed GR4J model}
+  %\VignetteIndexEntry{Simulated vs observed upstream flows in calibration of semi-distributed GR4J model}
   %\VignetteEncoding{UTF-8}
 ---
 
@@ -124,10 +124,13 @@ OutputsCalibDown1 <- Calibration_Michel(InputsModel = InputsModelDown1,
                                         FUN_MOD = RunModel_GR4J)
 ```
 
-To run the complete model, we should substitute the observed upstream flow by the simulated one:
+To run the complete model, we should substitute the observed upstream flow by the simulated one for the entire period of simulation (warm-up + run):
 
 ```{r}
 InputsModelDown2 <- InputsModelDown1
+# Simulated flow during warm-up period 
+InputsModelDown2$Qupstream[Ind_Run[seq_len(365)] - 365] <- OutputsModelUp$WarmUpQsim
+# Simulated flow during run period
 InputsModelDown2$Qupstream[Ind_Run] <- OutputsModelUp$Qsim
 ```