diff --git a/R/Calibration.R b/R/Calibration.R
index d526294f62ef1d2f7a66a1e2180f72b9e506d75f..319dbb4b2f9d083f64631e49dc1b65953012e434 100644
--- a/R/Calibration.R
+++ b/R/Calibration.R
@@ -1,29 +1,29 @@
-Calibration <- function(InputsModel,
-                        RunOptions,
-                        InputsCrit,
-                        CalibOptions,
-                        FUN_MOD,
-                        FUN_CRIT,                      # deprecated
-                        FUN_CALIB = Calibration_Michel,
-                        FUN_TRANSFO = NULL,
-                        verbose = TRUE,
-                        ...) {
-
-  FUN_MOD   <- match.fun(FUN_MOD)
-
-  if (!missing(FUN_CRIT)) {
-    FUN_CRIT <- match.fun(FUN_CRIT)
-  }
-
-  FUN_CALIB <- match.fun(FUN_CALIB)
-
-  if (!is.null(FUN_TRANSFO)) {
-    FUN_TRANSFO <- match.fun(FUN_TRANSFO)
-  }
-
-  return(FUN_CALIB(InputsModel = InputsModel, RunOptions = RunOptions, InputsCrit = InputsCrit,
-                   CalibOptions = CalibOptions,
-                   FUN_MOD = FUN_MOD, FUN_TRANSFO = FUN_TRANSFO,
-                   verbose = verbose, ...))
-}
-
+Calibration <- function(InputsModel,
+                        RunOptions,
+                        InputsCrit,
+                        CalibOptions,
+                        FUN_MOD,
+                        FUN_SD,
+                        FUN_CRIT,                      # deprecated
+                        FUN_CALIB = Calibration_Michel,
+                        FUN_TRANSFO = NULL,
+                        verbose = TRUE,
+                        ...) {
+
+
+  FUN_MOD   <- match.fun(FUN_MOD)
+  if (!missing(FUN_CRIT)) {
+    FUN_CRIT <- match.fun(FUN_CRIT)
+  }
+
+  FUN_CALIB <- match.fun(FUN_CALIB)
+
+  if (!is.null(FUN_TRANSFO)) {
+    FUN_TRANSFO <- match.fun(FUN_TRANSFO)
+  }
+  return(FUN_CALIB(InputsModel = InputsModel, RunOptions = RunOptions, InputsCrit = InputsCrit,
+                   CalibOptions = CalibOptions,
+                   FUN_MOD = FUN_MOD, FUN_SD, FUN_TRANSFO = FUN_TRANSFO,
+                   verbose = verbose, ...))
+}
+
diff --git a/R/Calibration_Michel.R b/R/Calibration_Michel.R
index 4a60a4ca1732b389d175da0315c6611fac48c86d..8a59765e57dfd2f28bbd0fe2d3b0af9e6d4e554d 100644
--- a/R/Calibration_Michel.R
+++ b/R/Calibration_Michel.R
@@ -1,440 +1,440 @@
-Calibration_Michel <- function(InputsModel,
-                               RunOptions,
-                               InputsCrit,
-                               CalibOptions,
-                               FUN_MOD,
-                               FUN_CRIT,           # deprecated
-                               FUN_TRANSFO = NULL,
-                               verbose = TRUE,
-                               ...) {
-
-
-  FUN_MOD  <- match.fun(FUN_MOD)
-
-  if (!missing(FUN_CRIT)) {
-    FUN_CRIT <- match.fun(FUN_CRIT)
-  }
-
-  # 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)) {
-    FUN_TRANSFO <- CalibOptions$FUN_TRANSFO
-  } else {
-    stop("'FUN_TRANSFO' is not provided neither as 'FUN_TRANSFO' argument or in 'CaliOptions' argument")
-  }
-
-  ##_____Arguments_check_____________________________________________________________________
-  if (!inherits(InputsModel, "InputsModel")) {
-    stop("'InputsModel' must be of class 'InputsModel'")
-  }
-  if (!inherits(RunOptions, "RunOptions")) {
-    stop("'RunOptions' must be of class 'RunOptions'")
-  }
-  if (!inherits(InputsCrit, "InputsCrit")) {
-    stop("'InputsCrit' must be of class 'InputsCrit'")
-  }
-  if (inherits(InputsCrit, "Multi")) {
-    stop("'InputsCrit' must be of class 'Single' or 'Compo'")
-  }
-  if (inherits(InputsCrit, "Single")) {
-    listVarObs <- InputsCrit$VarObs
-  }
-  if (inherits(InputsCrit, "Compo")) {
-    listVarObs <- sapply(InputsCrit, FUN = "[[", "VarObs")
-  }
-  if ("SCA" %in% listVarObs & !"Gratio" %in% RunOptions$Outputs_Cal) {
-    warning("Missing 'Gratio' is automatically added to 'Output_Cal' in 'RunOptions' as it is necessary in the objective function for comparison with SCA")
-    RunOptions$Outputs_Cal <- c(RunOptions$Outputs_Cal, "Gratio")
-  }
-  if ("SWE" %in% listVarObs & !"SnowPack" %in% RunOptions$Outputs_Cal) {
-    warning("Missing 'SnowPack' is automatically added to 'Output_Cal' in 'RunOptions' as it is necessary in the objective function for comparison with SWE")
-    RunOptions$Outputs_Cal <- c(RunOptions$Outputs_Cal, "SnowPack")
-  }
-  if (!inherits(CalibOptions, "CalibOptions")) {
-    stop("'CalibOptions' must be of class 'CalibOptions'")
-  }
-  if (!inherits(CalibOptions, "HBAN")) {
-    stop("'CalibOptions' must be of class 'HBAN' if 'Calibration_Michel' is used")
-  }
-  if (!missing(FUN_CRIT)) {
-    warning("argument 'FUN_CRIT' is deprecated. The error criterion function is now automatically get from the 'InputsCrit' object")
-  }
-
-
-  ##_variables_initialisation
-  ParamFinalR <- NULL
-  ParamFinalT <- NULL
-  CritFinal   <- NULL
-  NRuns <- 0
-  NIter <- 0
-  if ("StartParamDistrib" %in% names(CalibOptions)) {
-    PrefilteringType <- 2
-  } else {
-    PrefilteringType <- 1
-  }
-  if (PrefilteringType == 1) {
-    NParam <- ncol(CalibOptions$StartParamList)
-  }
-  if (PrefilteringType == 2) {
-    NParam <- ncol(CalibOptions$StartParamDistrib)
-  }
-  if (NParam > 20) {
-    stop("Calibration_Michel can handle a maximum of 20 parameters")
-  }
-  HistParamR    <- matrix(NA, nrow = 500 * NParam, ncol = NParam)
-  HistParamT    <- matrix(NA, nrow = 500 * NParam, ncol = NParam)
-  HistCrit      <- matrix(NA, nrow = 500 * NParam, ncol = 1)
-  CritName      <- NULL
-  CritBestValue <- NULL
-  Multiplier    <- NULL
-  CritOptim     <- +1e100
-  ##_temporary_change_of_Outputs_Sim
-  RunOptions$Outputs_Sim <- RunOptions$Outputs_Cal  ### this reduces the size of the matrix exchange with fortran and therefore speeds the calibration
-
-
-
-  ##_____Parameter_Grid_Screening____________________________________________________________
-
-
-  ##Definition_of_the_function_creating_all_possible_parameter_sets_from_different_values_for_each_parameter
-  ProposeCandidatesGrid <- function(DistribParam) {
-    NewCandidates <- expand.grid(lapply(seq_len(ncol(DistribParamR)), function(x) DistribParam[, x]))
-    NewCandidates <- unique(NewCandidates) # to avoid duplicates when a parameter is set
-    Output <- list(NewCandidates = NewCandidates)
-  }
-
-
-  ##Creation_of_new_candidates_______________________________________________
-  OptimParam <- is.na(CalibOptions$FixedParam)
-  if (PrefilteringType == 1) {
-    CandidatesParamR <- CalibOptions$StartParamList
-  }
-  if (PrefilteringType == 2) {
-    DistribParamR <- CalibOptions$StartParamDistrib
-    DistribParamR[, !OptimParam] <- NA
-    CandidatesParamR <- ProposeCandidatesGrid(DistribParamR)$NewCandidates
-  }
-  ##Remplacement_of_non_optimised_values_____________________________________
-  CandidatesParamR <- apply(CandidatesParamR, 1, function(x) {
-    x[!OptimParam] <- CalibOptions$FixedParam[!OptimParam]
-    return(x)
-  })
-  if (NParam > 1) {
-    CandidatesParamR <- t(CandidatesParamR)
-  } else {
-    CandidatesParamR <- cbind(CandidatesParamR)
-  }
-
-  ##Loop_to_test_the_various_candidates______________________________________
-  iNewOptim <- 0
-  Ncandidates <- nrow(CandidatesParamR)
-  if (verbose & Ncandidates > 1) {
-    if (PrefilteringType == 1) {
-      message("List-Screening in progress (", appendLF = FALSE)
-    }
-    if (PrefilteringType == 2) {
-      message("Grid-Screening in progress (", appendLF = FALSE)
-    }
-    message("0%", appendLF = FALSE)
-  }
-  for (iNew in 1:nrow(CandidatesParamR)) {
-    if (verbose & Ncandidates > 1) {
-      for (k in c(2, 4, 6, 8)) {
-        if (iNew == round(k / 10 * Ncandidates)) {
-          message(" ", 10 * k, "%", appendLF = FALSE)
-        }
-      }
-    }
-    ##Model_run
-    Param <- CandidatesParamR[iNew, ]
-    OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, ...)
-
-    ##Calibration_criterion_computation
-    OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
-    if (!is.na(OutputsCrit$CritValue)) {
-      if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
-        CritOptim <- OutputsCrit$CritValue * OutputsCrit$Multiplier
-        iNewOptim <- iNew
-      }
-    }
-    ##Storage_of_crit_info
-    if (is.null(CritName) | is.null(CritBestValue) | is.null(Multiplier)) {
-      CritName      <- OutputsCrit$CritName
-      CritBestValue <- OutputsCrit$CritBestValue
-      Multiplier    <- OutputsCrit$Multiplier
-    }
-  }
-  if (verbose & Ncandidates > 1) {
-    message(" 100%)\n", appendLF = FALSE)
-  }
-
-
-  ##End_of_first_step_Parameter_Screening____________________________________
-  ParamStartR <- CandidatesParamR[iNewOptim, ]
-  if (!is.matrix(ParamStartR)) {
-    ParamStartR <- matrix(ParamStartR, nrow = 1)
-  }
-  ParamStartT <- FUN_TRANSFO(ParamStartR, "RT")
-  CritStart   <- CritOptim
-  NRuns       <- NRuns + nrow(CandidatesParamR)
-  if (verbose) {
-    if (Ncandidates > 1) {
-      message(sprintf("\t Screening completed (%s runs)", NRuns))
-    }
-    if (Ncandidates == 1) {
-      message("\t Starting point for steepest-descent local search:")
-    }
-    message("\t     Param = ", paste(sprintf("%8.3f", ParamStartR), collapse = ", "))
-    message(sprintf("\t     Crit. %-12s = %.4f", CritName, CritStart * Multiplier))
-  }
-  ##Results_archiving________________________________________________________
-  HistParamR[1, ] <- ParamStartR
-  HistParamT[1, ] <- ParamStartT
-  HistCrit[1, ]   <- CritStart
-
-
-
-
-  ##_____Steepest_Descent_Local_Search_______________________________________________________
-
-
-  ##Definition_of_the_function_creating_new_parameter_sets_through_a_step_by_step_progression_procedure
-  ProposeCandidatesLoc <- function(NewParamOptimT, OldParamOptimT, RangesT, OptimParam, Pace) {
-    ##Format_checking
-    if (nrow(NewParamOptimT) != 1 | nrow(OldParamOptimT) != 1) {
-      stop("each input set must be a matrix of one single line")
-    }
-    if (ncol(NewParamOptimT)!=ncol(OldParamOptimT) | ncol(NewParamOptimT) != length(OptimParam)) {
-      stop("each input set must have the same number of values")
-    }
-    ##Proposal_of_new_parameter_sets ###(local search providing 2 * NParam-1 new sets)
-    NParam <- ncol(NewParamOptimT)
-    VECT <- NULL
-    for (I in 1:NParam) {
-      ##We_check_that_the_current_parameter_should_indeed_be_optimised
-      if (OptimParam[I]) {
-        for (J in 1:2) {
-          Sign <- 2 * J - 3   #Sign can be equal to -1 or +1
-          ##We_define_the_new_potential_candidate
-          Add <- TRUE
-          PotentialCandidateT <- NewParamOptimT
-          PotentialCandidateT[1, I] <- NewParamOptimT[I] + Sign * Pace
-          ##If_we_exit_the_range_of_possible_values_we_go_back_on_the_boundary
-          if (PotentialCandidateT[1, I] < RangesT[1, I] ) {
-            PotentialCandidateT[1, I] <- RangesT[1, I]
-          }
-          if (PotentialCandidateT[1, I] > RangesT[2, I]) {
-            PotentialCandidateT[1, I] <- RangesT[2, I]
-          }
-          ##We_check_the_set_is_not_outside_the_range_of_possible_values
-          if (NewParamOptimT[I] == RangesT[1, I] & Sign < 0) {
-            Add <- FALSE
-          }
-          if (NewParamOptimT[I] == RangesT[2, I] & Sign > 0) {
-            Add <- FALSE
-          }
-          ##We_check_that_this_set_has_not_been_tested_during_the_last_iteration
-          if (identical(PotentialCandidateT, OldParamOptimT)) {
-            Add <- FALSE
-          }
-          ##We_add_the_candidate_to_our_list
-          if (Add) {
-            VECT <- c(VECT, PotentialCandidateT)
-          }
-        }
-      }
-    }
-    Output <- NULL
-    Output$NewCandidatesT <- matrix(VECT, ncol = NParam, byrow = TRUE)
-    return(Output)
-  }
-
-
-  ##Initialisation_of_variables
-  if (verbose) {
-    message("Steepest-descent local search in progress")
-  }
-  Pace <- 0.64
-  PaceDiag <- rep(0, NParam)
-  CLG <- 0.7^(1 / NParam)
-  Compt <- 0
-  CritOptim <- CritStart
-  ##Conversion_of_real_parameter_values
-  RangesR <- CalibOptions$SearchRanges
-  RangesT <- FUN_TRANSFO(RangesR, "RT")
-  NewParamOptimT <- ParamStartT
-  OldParamOptimT <- ParamStartT
-
-
-  ##START_LOOP_ITER_________________________________________________________
-  for (ITER in 1:(100 * NParam)) {
-
-
-    ##Exit_loop_when_Pace_becomes_too_small___________________________________
-    if (Pace < 0.01) {
-      break
-    }
-
-
-    ##Creation_of_new_candidates______________________________________________
-    CandidatesParamT <- ProposeCandidatesLoc(NewParamOptimT, OldParamOptimT, RangesT, OptimParam, Pace)$NewCandidatesT
-    CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
-    ##Remplacement_of_non_optimised_values_____________________________________
-    CandidatesParamR <- apply(CandidatesParamR, 1, function(x) {
-      x[!OptimParam] <- CalibOptions$FixedParam[!OptimParam]
-      return(x)
-    })
-    if (NParam > 1) {
-      CandidatesParamR <- t(CandidatesParamR)
-    } else {
-      CandidatesParamR <- cbind(CandidatesParamR)
-    }
-
-
-    ##Loop_to_test_the_various_candidates_____________________________________
-    iNewOptim <- 0
-    for (iNew in 1:nrow(CandidatesParamR)) {
-      ##Model_run
-      Param <- CandidatesParamR[iNew, ]
-      OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, ...)
-      ##Calibration_criterion_computation
-      OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
-      if (!is.na(OutputsCrit$CritValue)) {
-        if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
-          CritOptim <- OutputsCrit$CritValue * OutputsCrit$Multiplier
-          iNewOptim <- iNew
-        }
-      }
-    }
-    NRuns <- NRuns + nrow(CandidatesParamR)
-
-
-    ##When_a_progress_has_been_achieved_______________________________________
-    if (iNewOptim != 0) {
-      ##We_store_the_optimal_set
-      OldParamOptimT <- NewParamOptimT
-      NewParamOptimT <- matrix(CandidatesParamT[iNewOptim, 1:NParam], nrow = 1)
-      Compt <- Compt + 1
-      ##When_necessary_we_increase_the_pace ### if_successive_progress_occur_in_a_row
-      if (Compt > 2 * NParam) {
-        Pace <- Pace * 2
-        Compt <- 0
-      }
-      ##We_update_PaceDiag
-      VectPace <- NewParamOptimT-OldParamOptimT
-      for (iC in 1:NParam) {
-        if (OptimParam[iC]) {
-          PaceDiag[iC] <- CLG * PaceDiag[iC] + (1-CLG) * VectPace[iC]
-        }
-      }
-    } else {
-      ##When_no_progress_has_been_achieved_we_decrease_the_pace_________________
-      Pace <- Pace / 2
-      Compt <- 0
-    }
-
-
-    ##Test_of_an_additional_candidate_using_diagonal_progress_________________
-    if (ITER > 4 * NParam) {
-      NRuns <- NRuns + 1
-      iNewOptim <- 0
-      iNew <- 1
-      CandidatesParamT <- NewParamOptimT+PaceDiag
-      if (!is.matrix(CandidatesParamT)) {
-        CandidatesParamT <- matrix(CandidatesParamT, nrow = 1)
-      }
-      ##If_we_exit_the_range_of_possible_values_we_go_back_on_the_boundary
-      for (iC in 1:NParam) {
-        if (OptimParam[iC]) {
-          if (CandidatesParamT[iNew, iC] < RangesT[1, iC]) {
-            CandidatesParamT[iNew, iC] <- RangesT[1, iC]
-          }
-          if (CandidatesParamT[iNew, iC] > RangesT[2, iC]) {
-            CandidatesParamT[iNew, iC] <- RangesT[2, iC]
-          }
-        }
-      }
-      CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
-      ##Model_run
-      Param <- CandidatesParamR[iNew, ]
-      OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, ...)
-      ##Calibration_criterion_computation
-      OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
-      if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
-        CritOptim <- OutputsCrit$CritValue * OutputsCrit$Multiplier
-        iNewOptim <- iNew
-      }
-      ##When_a_progress_has_been_achieved
-      if (iNewOptim != 0) {
-        OldParamOptimT <- NewParamOptimT
-        NewParamOptimT <- matrix(CandidatesParamT[iNewOptim, 1:NParam], nrow = 1)
-      }
-
-    }
-
-
-    ##Results_archiving_______________________________________________________
-    NewParamOptimR       <- FUN_TRANSFO(NewParamOptimT, "TR")
-    HistParamR[ITER+1, ] <- NewParamOptimR
-    HistParamT[ITER+1, ] <- NewParamOptimT
-    HistCrit[ITER+1, ]   <- CritOptim
-    ### if (verbose) { cat(paste("\t     Iter ",formatC(ITER,format="d",width=3), "    Crit ",formatC(CritOptim,format="f",digits=4), "    Pace ",formatC(Pace,format="f",digits=4), "\n",sep=""))}
-
-
-
-  } ##END_LOOP_ITER_________________________________________________________
-  ITER <- ITER - 1
-
-
-  ##Case_when_the_starting_parameter_set_remains_the_best_solution__________
-  if (CritOptim == CritStart & verbose) {
-    message("\t No progress achieved")
-  }
-
-  ##End_of_Steepest_Descent_Local_Search____________________________________
-  ParamFinalR <- NewParamOptimR
-  ParamFinalT <- NewParamOptimT
-  CritFinal   <- CritOptim
-  NIter       <- 1 + ITER
-  if (verbose) {
-    message(sprintf("\t Calibration completed (%s iterations, %s runs)", NIter, NRuns))
-    message("\t     Param = ", paste(sprintf("%8.3f", ParamFinalR), collapse = ", "))
-    message(sprintf("\t     Crit. %-12s = %.4f", CritName, CritFinal * Multiplier))
-    if (inherits(InputsCrit, "Compo")) {
-      listweights  <- OutputsCrit$CritCompo$MultiCritWeights
-      listNameCrit <- OutputsCrit$CritCompo$MultiCritNames
-      msgForm <- paste(sprintf("%.2f", listweights), listNameCrit, sep = " * ", collapse = ", ")
-      msgForm <- unlist(strsplit(msgForm, split = ","))
-      msgFormSep <- rep(c(",", ",", ",\n\t\t    "), times = ceiling(length(msgForm)/3))[1:length(msgForm)]
-      msgForm <- paste(msgForm, msgFormSep, sep = "", collapse = "")
-      msgForm <- gsub("\\,\\\n\\\t\\\t    $|\\,$", "", msgForm)
-      message("\tFormula: sum(", msgForm, ")")
-    }
-  }
-  ##Results_archiving_______________________________________________________
-  HistParamR <- cbind(HistParamR[1:NIter, ])
-  colnames(HistParamR) <- paste0("Param", 1:NParam)
-  HistParamT <- cbind(HistParamT[1:NIter, ])
-  colnames(HistParamT) <- paste0("Param", 1:NParam)
-  HistCrit   <- cbind(HistCrit[1:NIter, ])
-  ###colnames(HistCrit) <- paste("HistCrit")
-
-  BoolCrit_Actual <- InputsCrit$BoolCrit
-  BoolCrit_Actual[OutputsCrit$Ind_notcomputed] <- FALSE
-  MatBoolCrit <- cbind(InputsCrit$BoolCrit, BoolCrit_Actual)
-  colnames(MatBoolCrit) <- c("BoolCrit_Requested", "BoolCrit_Actual")
-
-
-  ##_____Output______________________________________________________________________________
-  OutputsCalib <- list(ParamFinalR = as.double(ParamFinalR), CritFinal = CritFinal * Multiplier,
-                       NIter = NIter, NRuns = NRuns,
-                       HistParamR = HistParamR, HistCrit = HistCrit * Multiplier,
-                       MatBoolCrit = MatBoolCrit,
-                       CritName = CritName, CritBestValue = CritBestValue)
-  class(OutputsCalib) <- c("OutputsCalib", "HBAN")
-  return(OutputsCalib)
-
-
-
-}
+Calibration_Michel <- function(InputsModel,
+                               RunOptions,
+                               InputsCrit,
+                               CalibOptions,
+                               FUN_MOD,
+                               FUN_SD,
+                               FUN_CRIT,           # deprecated
+                               FUN_TRANSFO = NULL,
+                               verbose = TRUE,
+                               ...) {
+
+  FUN_MOD  <- match.fun(FUN_MOD)
+
+  if (!missing(FUN_CRIT)) {
+    FUN_CRIT <- match.fun(FUN_CRIT)
+  }
+
+  # 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)) {
+    FUN_TRANSFO <- CalibOptions$FUN_TRANSFO
+  } else {
+    stop("'FUN_TRANSFO' is not provided neither as 'FUN_TRANSFO' argument or in 'CaliOptions' argument")
+  }
+
+  ##_____Arguments_check_____________________________________________________________________
+  if (!inherits(InputsModel, "InputsModel")) {
+    stop("'InputsModel' must be of class 'InputsModel'")
+  }
+  if (!inherits(RunOptions, "RunOptions")) {
+    stop("'RunOptions' must be of class 'RunOptions'")
+  }
+  if (!inherits(InputsCrit, "InputsCrit")) {
+    stop("'InputsCrit' must be of class 'InputsCrit'")
+  }
+  if (inherits(InputsCrit, "Multi")) {
+    stop("'InputsCrit' must be of class 'Single' or 'Compo'")
+  }
+  if (inherits(InputsCrit, "Single")) {
+    listVarObs <- InputsCrit$VarObs
+  }
+  if (inherits(InputsCrit, "Compo")) {
+    listVarObs <- sapply(InputsCrit, FUN = "[[", "VarObs")
+  }
+  if ("SCA" %in% listVarObs & !"Gratio" %in% RunOptions$Outputs_Cal) {
+    warning("Missing 'Gratio' is automatically added to 'Output_Cal' in 'RunOptions' as it is necessary in the objective function for comparison with SCA")
+    RunOptions$Outputs_Cal <- c(RunOptions$Outputs_Cal, "Gratio")
+  }
+  if ("SWE" %in% listVarObs & !"SnowPack" %in% RunOptions$Outputs_Cal) {
+    warning("Missing 'SnowPack' is automatically added to 'Output_Cal' in 'RunOptions' as it is necessary in the objective function for comparison with SWE")
+    RunOptions$Outputs_Cal <- c(RunOptions$Outputs_Cal, "SnowPack")
+  }
+  if (!inherits(CalibOptions, "CalibOptions")) {
+    stop("'CalibOptions' must be of class 'CalibOptions'")
+  }
+  if (!inherits(CalibOptions, "HBAN")) {
+    stop("'CalibOptions' must be of class 'HBAN' if 'Calibration_Michel' is used")
+  }
+  if (!missing(FUN_CRIT)) {
+    warning("argument 'FUN_CRIT' is deprecated. The error criterion function is now automatically get from the 'InputsCrit' object")
+  }
+
+
+  ##_variables_initialisation
+  ParamFinalR <- NULL
+  ParamFinalT <- NULL
+  CritFinal   <- NULL
+  NRuns <- 0
+  NIter <- 0
+  if ("StartParamDistrib" %in% names(CalibOptions)) {
+    PrefilteringType <- 2
+  } else {
+    PrefilteringType <- 1
+  }
+  if (PrefilteringType == 1) {
+    NParam <- ncol(CalibOptions$StartParamList)
+  }
+  if (PrefilteringType == 2) {
+    NParam <- ncol(CalibOptions$StartParamDistrib)
+  }
+  if (NParam > 20) {
+    stop("Calibration_Michel can handle a maximum of 20 parameters")
+  }
+  HistParamR    <- matrix(NA, nrow = 500 * NParam, ncol = NParam)
+  HistParamT    <- matrix(NA, nrow = 500 * NParam, ncol = NParam)
+  HistCrit      <- matrix(NA, nrow = 500 * NParam, ncol = 1)
+  CritName      <- NULL
+  CritBestValue <- NULL
+  Multiplier    <- NULL
+  CritOptim     <- +1e100
+  ##_temporary_change_of_Outputs_Sim
+  RunOptions$Outputs_Sim <- RunOptions$Outputs_Cal  ### this reduces the size of the matrix exchange with fortran and therefore speeds the calibration
+
+
+
+  ##_____Parameter_Grid_Screening____________________________________________________________
+
+
+  ##Definition_of_the_function_creating_all_possible_parameter_sets_from_different_values_for_each_parameter
+  ProposeCandidatesGrid <- function(DistribParam) {
+    NewCandidates <- expand.grid(lapply(seq_len(ncol(DistribParamR)), function(x) DistribParam[, x]))
+    NewCandidates <- unique(NewCandidates) # to avoid duplicates when a parameter is set
+    Output <- list(NewCandidates = NewCandidates)
+  }
+
+
+  ##Creation_of_new_candidates_______________________________________________
+  OptimParam <- is.na(CalibOptions$FixedParam)
+  if (PrefilteringType == 1) {
+    CandidatesParamR <- CalibOptions$StartParamList
+  }
+  if (PrefilteringType == 2) {
+    DistribParamR <- CalibOptions$StartParamDistrib
+    DistribParamR[, !OptimParam] <- NA
+    CandidatesParamR <- ProposeCandidatesGrid(DistribParamR)$NewCandidates
+  }
+  ##Remplacement_of_non_optimised_values_____________________________________
+  CandidatesParamR <- apply(CandidatesParamR, 1, function(x) {
+    x[!OptimParam] <- CalibOptions$FixedParam[!OptimParam]
+    return(x)
+  })
+  if (NParam > 1) {
+    CandidatesParamR <- t(CandidatesParamR)
+  } else {
+    CandidatesParamR <- cbind(CandidatesParamR)
+  }
+
+  ##Loop_to_test_the_various_candidates______________________________________
+  iNewOptim <- 0
+  Ncandidates <- nrow(CandidatesParamR)
+  if (verbose & Ncandidates > 1) {
+    if (PrefilteringType == 1) {
+      message("List-Screening in progress (", appendLF = FALSE)
+    }
+    if (PrefilteringType == 2) {
+      message("Grid-Screening in progress (", appendLF = FALSE)
+    }
+    message("0%", appendLF = FALSE)
+  }
+  for (iNew in 1:nrow(CandidatesParamR)) {
+    if (verbose & Ncandidates > 1) {
+      for (k in c(2, 4, 6, 8)) {
+        if (iNew == round(k / 10 * Ncandidates)) {
+          message(" ", 10 * k, "%", appendLF = FALSE)
+        }
+      }
+    }
+    ##Model_run
+    Param <- CandidatesParamR[iNew, ]
+    OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, FUN_SD = FUN_SD,...)
+
+    ##Calibration_criterion_computation
+    OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
+    if (!is.na(OutputsCrit$CritValue)) {
+      if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
+        CritOptim <- OutputsCrit$CritValue * OutputsCrit$Multiplier
+        iNewOptim <- iNew
+      }
+    }
+    ##Storage_of_crit_info
+    if (is.null(CritName) | is.null(CritBestValue) | is.null(Multiplier)) {
+      CritName      <- OutputsCrit$CritName
+      CritBestValue <- OutputsCrit$CritBestValue
+      Multiplier    <- OutputsCrit$Multiplier
+    }
+  }
+  if (verbose & Ncandidates > 1) {
+    message(" 100%)\n", appendLF = FALSE)
+  }
+
+
+  ##End_of_first_step_Parameter_Screening____________________________________
+  ParamStartR <- CandidatesParamR[iNewOptim, ]
+  if (!is.matrix(ParamStartR)) {
+    ParamStartR <- matrix(ParamStartR, nrow = 1)
+  }
+  ParamStartT <- FUN_TRANSFO(ParamStartR, "RT")
+  CritStart   <- CritOptim
+  NRuns       <- NRuns + nrow(CandidatesParamR)
+  if (verbose) {
+    if (Ncandidates > 1) {
+      message(sprintf("\t Screening completed (%s runs)", NRuns))
+    }
+    if (Ncandidates == 1) {
+      message("\t Starting point for steepest-descent local search:")
+    }
+    message("\t     Param = ", paste(sprintf("%8.3f", ParamStartR), collapse = ", "))
+    message(sprintf("\t     Crit. %-12s = %.4f", CritName, CritStart * Multiplier))
+  }
+  ##Results_archiving________________________________________________________
+  HistParamR[1, ] <- ParamStartR
+  HistParamT[1, ] <- ParamStartT
+  HistCrit[1, ]   <- CritStart
+
+
+
+
+  ##_____Steepest_Descent_Local_Search_______________________________________________________
+
+
+  ##Definition_of_the_function_creating_new_parameter_sets_through_a_step_by_step_progression_procedure
+  ProposeCandidatesLoc <- function(NewParamOptimT, OldParamOptimT, RangesT, OptimParam, Pace) {
+    ##Format_checking
+    if (nrow(NewParamOptimT) != 1 | nrow(OldParamOptimT) != 1) {
+      stop("each input set must be a matrix of one single line")
+    }
+    if (ncol(NewParamOptimT)!=ncol(OldParamOptimT) | ncol(NewParamOptimT) != length(OptimParam)) {
+      stop("each input set must have the same number of values")
+    }
+    ##Proposal_of_new_parameter_sets ###(local search providing 2 * NParam-1 new sets)
+    NParam <- ncol(NewParamOptimT)
+    VECT <- NULL
+    for (I in 1:NParam) {
+      ##We_check_that_the_current_parameter_should_indeed_be_optimised
+      if (OptimParam[I]) {
+        for (J in 1:2) {
+          Sign <- 2 * J - 3   #Sign can be equal to -1 or +1
+          ##We_define_the_new_potential_candidate
+          Add <- TRUE
+          PotentialCandidateT <- NewParamOptimT
+          PotentialCandidateT[1, I] <- NewParamOptimT[I] + Sign * Pace
+          ##If_we_exit_the_range_of_possible_values_we_go_back_on_the_boundary
+          if (PotentialCandidateT[1, I] < RangesT[1, I] ) {
+            PotentialCandidateT[1, I] <- RangesT[1, I]
+          }
+          if (PotentialCandidateT[1, I] > RangesT[2, I]) {
+            PotentialCandidateT[1, I] <- RangesT[2, I]
+          }
+          ##We_check_the_set_is_not_outside_the_range_of_possible_values
+          if (NewParamOptimT[I] == RangesT[1, I] & Sign < 0) {
+            Add <- FALSE
+          }
+          if (NewParamOptimT[I] == RangesT[2, I] & Sign > 0) {
+            Add <- FALSE
+          }
+          ##We_check_that_this_set_has_not_been_tested_during_the_last_iteration
+          if (identical(PotentialCandidateT, OldParamOptimT)) {
+            Add <- FALSE
+          }
+          ##We_add_the_candidate_to_our_list
+          if (Add) {
+            VECT <- c(VECT, PotentialCandidateT)
+          }
+        }
+      }
+    }
+    Output <- NULL
+    Output$NewCandidatesT <- matrix(VECT, ncol = NParam, byrow = TRUE)
+    return(Output)
+  }
+
+
+  ##Initialisation_of_variables
+  if (verbose) {
+    message("Steepest-descent local search in progress")
+  }
+  Pace <- 0.64
+  PaceDiag <- rep(0, NParam)
+  CLG <- 0.7^(1 / NParam)
+  Compt <- 0
+  CritOptim <- CritStart
+  ##Conversion_of_real_parameter_values
+  RangesR <- CalibOptions$SearchRanges
+  RangesT <- FUN_TRANSFO(RangesR, "RT")
+  NewParamOptimT <- ParamStartT
+  OldParamOptimT <- ParamStartT
+
+
+  ##START_LOOP_ITER_________________________________________________________
+  for (ITER in 1:(100 * NParam)) {
+
+
+    ##Exit_loop_when_Pace_becomes_too_small___________________________________
+    if (Pace < 0.01) {
+      break
+    }
+
+
+    ##Creation_of_new_candidates______________________________________________
+    CandidatesParamT <- ProposeCandidatesLoc(NewParamOptimT, OldParamOptimT, RangesT, OptimParam, Pace)$NewCandidatesT
+    CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
+    ##Remplacement_of_non_optimised_values_____________________________________
+    CandidatesParamR <- apply(CandidatesParamR, 1, function(x) {
+      x[!OptimParam] <- CalibOptions$FixedParam[!OptimParam]
+      return(x)
+    })
+    if (NParam > 1) {
+      CandidatesParamR <- t(CandidatesParamR)
+    } else {
+      CandidatesParamR <- cbind(CandidatesParamR)
+    }
+
+
+    ##Loop_to_test_the_various_candidates_____________________________________
+    iNewOptim <- 0
+    for (iNew in 1:nrow(CandidatesParamR)) {
+      ##Model_run
+      Param <- CandidatesParamR[iNew, ]
+      OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, ...)
+      ##Calibration_criterion_computation
+      OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
+      if (!is.na(OutputsCrit$CritValue)) {
+        if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
+          CritOptim <- OutputsCrit$CritValue * OutputsCrit$Multiplier
+          iNewOptim <- iNew
+        }
+      }
+    }
+    NRuns <- NRuns + nrow(CandidatesParamR)
+
+
+    ##When_a_progress_has_been_achieved_______________________________________
+    if (iNewOptim != 0) {
+      ##We_store_the_optimal_set
+      OldParamOptimT <- NewParamOptimT
+      NewParamOptimT <- matrix(CandidatesParamT[iNewOptim, 1:NParam], nrow = 1)
+      Compt <- Compt + 1
+      ##When_necessary_we_increase_the_pace ### if_successive_progress_occur_in_a_row
+      if (Compt > 2 * NParam) {
+        Pace <- Pace * 2
+        Compt <- 0
+      }
+      ##We_update_PaceDiag
+      VectPace <- NewParamOptimT-OldParamOptimT
+      for (iC in 1:NParam) {
+        if (OptimParam[iC]) {
+          PaceDiag[iC] <- CLG * PaceDiag[iC] + (1-CLG) * VectPace[iC]
+        }
+      }
+    } else {
+      ##When_no_progress_has_been_achieved_we_decrease_the_pace_________________
+      Pace <- Pace / 2
+      Compt <- 0
+    }
+
+
+    ##Test_of_an_additional_candidate_using_diagonal_progress_________________
+    if (ITER > 4 * NParam) {
+      NRuns <- NRuns + 1
+      iNewOptim <- 0
+      iNew <- 1
+      CandidatesParamT <- NewParamOptimT+PaceDiag
+      if (!is.matrix(CandidatesParamT)) {
+        CandidatesParamT <- matrix(CandidatesParamT, nrow = 1)
+      }
+      ##If_we_exit_the_range_of_possible_values_we_go_back_on_the_boundary
+      for (iC in 1:NParam) {
+        if (OptimParam[iC]) {
+          if (CandidatesParamT[iNew, iC] < RangesT[1, iC]) {
+            CandidatesParamT[iNew, iC] <- RangesT[1, iC]
+          }
+          if (CandidatesParamT[iNew, iC] > RangesT[2, iC]) {
+            CandidatesParamT[iNew, iC] <- RangesT[2, iC]
+          }
+        }
+      }
+      CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
+      ##Model_run
+      Param <- CandidatesParamR[iNew, ]
+      OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, ...)
+      ##Calibration_criterion_computation
+      OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
+      if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
+        CritOptim <- OutputsCrit$CritValue * OutputsCrit$Multiplier
+        iNewOptim <- iNew
+      }
+      ##When_a_progress_has_been_achieved
+      if (iNewOptim != 0) {
+        OldParamOptimT <- NewParamOptimT
+        NewParamOptimT <- matrix(CandidatesParamT[iNewOptim, 1:NParam], nrow = 1)
+      }
+
+    }
+
+
+    ##Results_archiving_______________________________________________________
+    NewParamOptimR       <- FUN_TRANSFO(NewParamOptimT, "TR")
+    HistParamR[ITER+1, ] <- NewParamOptimR
+    HistParamT[ITER+1, ] <- NewParamOptimT
+    HistCrit[ITER+1, ]   <- CritOptim
+    ### if (verbose) { cat(paste("\t     Iter ",formatC(ITER,format="d",width=3), "    Crit ",formatC(CritOptim,format="f",digits=4), "    Pace ",formatC(Pace,format="f",digits=4), "\n",sep=""))}
+
+
+
+  } ##END_LOOP_ITER_________________________________________________________
+  ITER <- ITER - 1
+
+
+  ##Case_when_the_starting_parameter_set_remains_the_best_solution__________
+  if (CritOptim == CritStart & verbose) {
+    message("\t No progress achieved")
+  }
+
+  ##End_of_Steepest_Descent_Local_Search____________________________________
+  ParamFinalR <- NewParamOptimR
+  ParamFinalT <- NewParamOptimT
+  CritFinal   <- CritOptim
+  NIter       <- 1 + ITER
+  if (verbose) {
+    message(sprintf("\t Calibration completed (%s iterations, %s runs)", NIter, NRuns))
+    message("\t     Param = ", paste(sprintf("%8.3f", ParamFinalR), collapse = ", "))
+    message(sprintf("\t     Crit. %-12s = %.4f", CritName, CritFinal * Multiplier))
+    if (inherits(InputsCrit, "Compo")) {
+      listweights  <- OutputsCrit$CritCompo$MultiCritWeights
+      listNameCrit <- OutputsCrit$CritCompo$MultiCritNames
+      msgForm <- paste(sprintf("%.2f", listweights), listNameCrit, sep = " * ", collapse = ", ")
+      msgForm <- unlist(strsplit(msgForm, split = ","))
+      msgFormSep <- rep(c(",", ",", ",\n\t\t    "), times = ceiling(length(msgForm)/3))[1:length(msgForm)]
+      msgForm <- paste(msgForm, msgFormSep, sep = "", collapse = "")
+      msgForm <- gsub("\\,\\\n\\\t\\\t    $|\\,$", "", msgForm)
+      message("\tFormula: sum(", msgForm, ")")
+    }
+  }
+  ##Results_archiving_______________________________________________________
+  HistParamR <- cbind(HistParamR[1:NIter, ])
+  colnames(HistParamR) <- paste0("Param", 1:NParam)
+  HistParamT <- cbind(HistParamT[1:NIter, ])
+  colnames(HistParamT) <- paste0("Param", 1:NParam)
+  HistCrit   <- cbind(HistCrit[1:NIter, ])
+  ###colnames(HistCrit) <- paste("HistCrit")
+
+  BoolCrit_Actual <- InputsCrit$BoolCrit
+  BoolCrit_Actual[OutputsCrit$Ind_notcomputed] <- FALSE
+  MatBoolCrit <- cbind(InputsCrit$BoolCrit, BoolCrit_Actual)
+  colnames(MatBoolCrit) <- c("BoolCrit_Requested", "BoolCrit_Actual")
+
+
+  ##_____Output______________________________________________________________________________
+  OutputsCalib <- list(ParamFinalR = as.double(ParamFinalR), CritFinal = CritFinal * Multiplier,
+                       NIter = NIter, NRuns = NRuns,
+                       HistParamR = HistParamR, HistCrit = HistCrit * Multiplier,
+                       MatBoolCrit = MatBoolCrit,
+                       CritName = CritName, CritBestValue = CritBestValue)
+  class(OutputsCalib) <- c("OutputsCalib", "HBAN")
+  return(OutputsCalib)
+
+
+
+}
diff --git a/R/CreateCalibOptions.R b/R/CreateCalibOptions.R
index 400b4d4d7dbcfddcb33ab68f19421a78ebb43b29..8a2e9a3bbeae5bc0262a719c8db547f01c0df6de 100644
--- a/R/CreateCalibOptions.R
+++ b/R/CreateCalibOptions.R
@@ -1,254 +1,257 @@
-CreateCalibOptions <- function(FUN_MOD,
-                               FUN_CALIB = Calibration_Michel,
-                               FUN_TRANSFO = NULL,
-                               IsHyst = FALSE,
-                               IsSD = FALSE,
-                               FixedParam = NULL,
-                               SearchRanges = NULL,
-                               StartParamList = NULL,
-                               StartParamDistrib = NULL) {
-
-  ObjectClass <- NULL
-
-  FUN_MOD     <- match.fun(FUN_MOD)
-  FUN_CALIB   <- match.fun(FUN_CALIB)
-  if (!is.null(FUN_TRANSFO)) {
-    FUN_TRANSFO <- match.fun(FUN_TRANSFO)
-  }
-  if (!is.logical(IsHyst) | length(IsHyst) != 1L) {
-    stop("'IsHyst' must be a logical of length 1")
-  }
-  if (!is.logical(IsSD) | length(IsSD) != 1L) {
-    stop("'IsSD' must be a logical of length 1")
-  }
-
-  ## check FUN_MOD
-  FeatFUN_MOD <- .GetFeatModel(FUN_MOD = FUN_MOD)
-  FeatFUN_MOD$IsHyst <- IsHyst
-  FeatFUN_MOD$IsSD <- IsSD
-  ObjectClass <- FeatFUN_MOD$Class
-
-  if (identical(FUN_MOD, RunModel_Lag) && IsSD) {
-      stop("RunModel_Lag should not be used with 'isSD=TRUE'")
-  }
-  if (IsHyst) {
-    ObjectClass <- c(ObjectClass, "hysteresis")
-  }
-  if (IsSD) {
-    ObjectClass <- c(ObjectClass, "SD")
-  }
-
-  ## check FUN_CALIB
-  BOOL <- FALSE
-
-  if (identical(FUN_CALIB, Calibration_Michel)) {
-    ObjectClass <- c(ObjectClass, "HBAN")
-    BOOL <- TRUE
-  }
-  if (!BOOL) {
-    stop("incorrect 'FUN_CALIB' for use in 'CreateCalibOptions'")
-    return(NULL)
-
-  }
-
-  ## check FUN_TRANSFO
-  if (is.null(FUN_TRANSFO)) {
-    FUN_TRANSFO <- .FunTransfo(FeatFUN_MOD)
-  }
-
-  ## NParam
-  NParam <- FeatFUN_MOD$NbParam
-
-  if (IsHyst) {
-    NParam <- NParam + 2
-  }
-  if (IsSD) {
-    NParam <- NParam + 1
-  }
-
-  ## check FixedParam
-  if (is.null(FixedParam)) {
-    FixedParam <- rep(NA, NParam)
-  } else {
-    if (!is.vector(FixedParam)) {
-      stop("FixedParam must be a vector")
-    }
-    if (length(FixedParam) != NParam) {
-      stop("Incompatibility between 'FixedParam' length and 'FUN_MOD'")
-    }
-    if (all(!is.na(FixedParam))) {
-      stop("At least one parameter must be not set (NA)")
-    }
-    if (all(is.na(FixedParam))) {
-      warning("You have not set any parameter in 'FixedParam'")
-    }
-  }
-
-  ## check SearchRanges
-  if (is.null(SearchRanges)) {
-    ParamT <- matrix(c(rep(-9.99, NParam), rep(+9.99, NParam)),
-                     ncol = NParam, byrow = TRUE)
-    SearchRanges <- TransfoParam(ParamIn = ParamT, Direction = "TR", FUN_TRANSFO = FUN_TRANSFO)
-
-  } else {
-    if (!is.matrix(SearchRanges)) {
-      stop("'SearchRanges' must be a matrix")
-    }
-    if (!is.numeric(SearchRanges)) {
-      stop("'SearchRanges' must be a matrix of numeric values")
-    }
-    if (sum(is.na(SearchRanges)) != 0) {
-      stop("'SearchRanges' must not include NA values")
-    }
-    if (nrow(SearchRanges) != 2) {
-      stop("'SearchRanges' must have 2 rows")
-    }
-    if (ncol(SearchRanges) != NParam) {
-      stop("Incompatibility between 'SearchRanges' ncol and 'FUN_MOD'")
-    }
-  }
-
-  ## check StartParamList and StartParamDistrib default values
-  if (("HBAN"  %in% ObjectClass & is.null(StartParamList) & is.null(StartParamDistrib))) {
-    if ("GR4H" == FeatFUN_MOD$CodeMod) {
-      ParamT <- matrix(c(+5.12, -1.18, +4.34, -9.69,
-                         +5.58, -0.85, +4.74, -9.47,
-                         +6.01, -0.50, +5.14, -8.87), ncol = 4, byrow = TRUE)
-    }
-    if (("GR5H" == FeatFUN_MOD$CodeMod) & ("interception" %in% ObjectClass)) {
-      ParamT <- matrix(c(+3.46, -1.25, +4.04, -9.53, -9.34,
-                         +3.74, -0.41, +4.78, -8.94, -3.33,
-                         +4.29, +0.16, +5.39, -7.39, +3.33), ncol = 5, byrow = TRUE)
-    }
-    if (("GR5H" == FeatFUN_MOD$CodeMod) & !("interception" %in% ObjectClass)) {
-      ParamT <- matrix(c(+3.28, -0.39, +4.14, -9.54, -7.49,
-                         +3.62, -0.19, +4.80, -9.00, -6.31,
-                         +4.01, -0.04, +5.43, -7.53, -5.33), ncol = 5, byrow = TRUE)
-    }
-    if ("GR4J" == FeatFUN_MOD$CodeMod) {
-      ParamT <- matrix(c(+5.13, -1.60, +3.03, -9.05,
-                         +5.51, -0.61, +3.74, -8.51,
-                         +6.07, -0.02, +4.42, -8.06), ncol = 4, byrow = TRUE)
-    }
-    if ("GR5J" == FeatFUN_MOD$CodeMod) {
-      ParamT <- matrix(c(+5.17, -1.13, +3.08, -9.37, -7.45,
-                         +5.55, -0.46, +3.75, -9.09, -4.69,
-                         +6.10, -0.11, +4.43, -8.60, -0.66), ncol = 5, byrow = TRUE)
-
-    }
-    if ("GR6J" == FeatFUN_MOD$CodeMod) {
-      ParamT <- matrix(c(+3.60, -1.00, +3.30, -9.10, -0.90, +3.00,
-                         +3.90, -0.50, +4.10, -8.70, +0.10, +4.00,
-                         +4.50, +0.50, +5.00, -8.10, +1.10, +5.00), ncol = 6, byrow = TRUE)
-    }
-    if ("GR2M" == FeatFUN_MOD$CodeMod) {
-      ParamT <- matrix(c(+5.03, -7.15,
-                         +5.22, -6.74,
-                         +5.85, -6.37), ncol = 2, byrow = TRUE)
-    }
-    if ("GR1A" == FeatFUN_MOD$CodeMod) {
-      ParamT <- matrix(c(-1.69,
-                         -0.38,
-                         +1.39), ncol = 1, byrow = TRUE)
-    }
-
-
-    if ("CemaNeige" == FeatFUN_MOD$CodeMod) {
-      ParamT <- matrix(c(-9.96, +6.63,
-                         -9.14, +6.90,
-                         +4.10, +7.21), ncol = 2, byrow = TRUE)
-    }
-    if ("CemaNeigeGR4H" == FeatFUN_MOD$CodeMod) {
-      ParamT <- matrix(c(+5.12, -1.18, +4.34, -9.69, -9.96, +6.63,
-                         +5.58, -0.85, +4.74, -9.47, -9.14, +6.90,
-                         +6.01, -0.50, +5.14, -8.87, +4.10, +7.21), ncol = 6, byrow = TRUE)
-    }
-    if (("CemaNeigeGR5H" == FeatFUN_MOD$CodeMod) & ("interception" %in% ObjectClass)) {
-      ParamT <- matrix(c(+3.46, -1.25, +4.04, -9.53, -9.34, -9.96, +6.63,
-                         +3.74, -0.41, +4.78, -8.94, -3.33, -9.14, +6.90,
-                         +4.29, +0.16, +5.39, -7.39, +3.33, +4.10, +7.21), ncol = 7, byrow = TRUE)
-    }
-    if (("CemaNeigeGR5H" == FeatFUN_MOD$CodeMod) & !("interception" %in% ObjectClass)) {
-      ParamT <- matrix(c(+3.28, -0.39, +4.14, -9.54, -7.49, -9.96, +6.63,
-                         +3.62, -0.19, +4.80, -9.00, -6.31, -9.14, +6.90,
-                         +4.01, -0.04, +5.43, -7.53, -5.33, +4.10, +7.21), ncol = 7, byrow = TRUE)
-    }
-    if ("CemaNeigeGR4J" == FeatFUN_MOD$CodeMod) {
-      ParamT <- matrix(c(+5.13, -1.60, +3.03, -9.05, -9.96, +6.63,
-                         +5.51, -0.61, +3.74, -8.51, -9.14, +6.90,
-                         +6.07, -0.02, +4.42, -8.06, +4.10, +7.21), ncol = 6, byrow = TRUE)
-    }
-    if ("CemaNeigeGR5J" == FeatFUN_MOD$CodeMod) {
-      ParamT <- matrix(c(+5.17, -1.13, +3.08, -9.37, -7.45, -9.96, +6.63,
-                         +5.55, -0.46, +3.75, -9.09, -4.69, -9.14, +6.90,
-                         +6.10, -0.11, +4.43, -8.60, -0.66, +4.10, +7.21), ncol = 7, byrow = TRUE)
-    }
-    if ("CemaNeigeGR6J" == FeatFUN_MOD$CodeMod) {
-      ParamT <- matrix(c(+3.60, -1.00, +3.30, -9.10, -0.90, +3.00, -9.96, +6.63,
-                         +3.90, -0.50, +4.10, -8.70, +0.10, +4.00, -9.14, +6.90,
-                         +4.50, +0.50, +5.00, -8.10, +1.10, +5.00, +4.10, +7.21), ncol = 8, byrow = TRUE)
-    }
-
-    if (IsHyst) {
-      ParamTHyst <- matrix(c(-7.00, -7.00,
-                             -0.00, -0.00,
-                             +7.00, +7.00), ncol = 2, byrow = TRUE)
-      ParamT <- cbind(ParamT, ParamTHyst)
-    }
-    if (IsSD) {
-      ParamTSD <- matrix(c(+1.25,
-                           +2.50,
-                           +5.00), ncol = 1, byrow = TRUE)
-      ParamT <- cbind(ParamTSD, ParamT)
-    }
-
-    StartParamList    <- NULL
-    StartParamDistrib <- TransfoParam(ParamIn = ParamT, Direction = "TR", FUN_TRANSFO = FUN_TRANSFO)
-
-  }
-
-  ## check StartParamList and StartParamDistrib format
-  if ("HBAN" %in% ObjectClass & !is.null(StartParamList)) {
-    if (!is.matrix(StartParamList)) {
-      stop("'StartParamList' must be a matrix")
-    }
-    if (!is.numeric(StartParamList)) {
-      stop("'StartParamList' must be a matrix of numeric values")
-    }
-    if (sum(is.na(StartParamList)) != 0) {
-      stop("'StartParamList' must not include NA values")
-    }
-    if (ncol(StartParamList) != NParam) {
-      stop("Incompatibility between 'StartParamList' ncol and 'FUN_MOD'")
-    }
-  }
-  if ("HBAN" %in% ObjectClass & !is.null(StartParamDistrib)) {
-    if (!is.matrix(StartParamDistrib)) {
-      stop("'StartParamDistrib' must be a matrix")
-    }
-    if (!is.numeric(StartParamDistrib[1, ])) {
-      stop("'StartParamDistrib' must be a matrix of numeric values")
-    }
-    if (sum(is.na(StartParamDistrib[1, ])) != 0) {
-      stop("'StartParamDistrib' must not include NA values on the first line")
-    }
-    if (ncol(StartParamDistrib) != NParam) {
-      stop("Incompatibility between 'StartParamDistrib' ncol and 'FUN_MOD'")
-    }
-  }
-
-
-  ##Create_CalibOptions
-  CalibOptions <- list(FixedParam = FixedParam, SearchRanges = SearchRanges, FUN_TRANSFO = FUN_TRANSFO)
-
-  if (!is.null(StartParamList)) {
-    CalibOptions <- c(CalibOptions, list(StartParamList = StartParamList))
-  }
-  if (!is.null(StartParamDistrib)) {
-    CalibOptions <- c(CalibOptions, list(StartParamDistrib = StartParamDistrib))
-  }
-  class(CalibOptions) <- c("CalibOptions", ObjectClass)
-
-  return(CalibOptions)
-
-}
+CreateCalibOptions <- function(FUN_MOD,
+                               FUN_CALIB = Calibration_Michel,
+                               FUN_TRANSFO = NULL,
+                               IsHyst = FALSE,
+                               IsSD = FALSE,
+                               FUN_SD = NULL,
+                               FixedParam = NULL,
+                               SearchRanges = NULL,
+                               StartParamList = NULL,
+                               StartParamDistrib = NULL) {
+
+  ObjectClass <- NULL
+
+  FUN_MOD     <- match.fun(FUN_MOD)
+  FUN_CALIB   <- match.fun(FUN_CALIB)
+  if (!is.null(FUN_TRANSFO)) {
+    FUN_TRANSFO <- match.fun(FUN_TRANSFO)
+  }
+  if (!is.logical(IsHyst) | length(IsHyst) != 1L) {
+    stop("'IsHyst' must be a logical of length 1")
+  }
+  if (!is.logical(IsSD) | length(IsSD) != 1L) {
+    stop("'IsSD' must be a logical of length 1")
+  }
+
+  ## check FUN_MOD
+  FeatFUN_MOD <- .GetFeatModel(FUN_MOD = FUN_MOD, FUN_SD = FUN_SD)
+  FeatFUN_MOD$IsHyst <- IsHyst
+  FeatFUN_MOD$IsSD <- IsSD
+  ObjectClass <- FeatFUN_MOD$Class
+
+  if (identical(FUN_MOD, RunModel_Lag) && IsSD) {
+      stop("RunModel_Lag should not be used with 'isSD=TRUE'")
+  }
+  if (IsHyst) {
+    ObjectClass <- c(ObjectClass, "hysteresis")
+  }
+  if (IsSD) {
+    ObjectClass <- c(ObjectClass, "SD")
+  }
+
+  ## check FUN_CALIB
+  BOOL <- FALSE
+
+  if (identical(FUN_CALIB, Calibration_Michel)) {
+    ObjectClass <- c(ObjectClass, "HBAN")
+    BOOL <- TRUE
+  }
+  if (!BOOL) {
+    stop("incorrect 'FUN_CALIB' for use in 'CreateCalibOptions'")
+    return(NULL)
+
+  }
+
+  ## check FUN_TRANSFO
+  if (is.null(FUN_TRANSFO)) {
+    FUN_TRANSFO <- .FunTransfo(FeatFUN_MOD, FUN_SD)
+  }
+
+  ## NParam
+  NParam <- FeatFUN_MOD$NbParam
+
+  if (IsHyst) {
+    NParam <- NParam + 2
+  }
+  if (IsSD) {
+    NParam <- NParam + FeatFUN_MOD$ParamSD
+  }
+
+  ## check FixedParam
+  if (is.null(FixedParam)) {
+    FixedParam <- rep(NA, NParam)
+  } else {
+    if (!is.vector(FixedParam)) {
+      stop("FixedParam must be a vector")
+    }
+    if (length(FixedParam) != NParam) {
+      stop("Incompatibility between 'FixedParam' length and 'FUN_MOD'")
+    }
+    if (all(!is.na(FixedParam))) {
+      stop("At least one parameter must be not set (NA)")
+    }
+    if (all(is.na(FixedParam))) {
+      warning("You have not set any parameter in 'FixedParam'")
+    }
+  }
+
+  ## check SearchRanges
+  if (is.null(SearchRanges)) {
+    ParamT <- matrix(c(rep(-9.99, NParam), rep(+9.99, NParam)),
+                     ncol = NParam, byrow = TRUE)
+
+    SearchRanges <- TransfoParam(ParamIn = ParamT, Direction = "TR", FUN_TRANSFO = FUN_TRANSFO)
+
+  } else {
+    if (!is.matrix(SearchRanges)) {
+      stop("'SearchRanges' must be a matrix")
+    }
+    if (!is.numeric(SearchRanges)) {
+      stop("'SearchRanges' must be a matrix of numeric values")
+    }
+    if (sum(is.na(SearchRanges)) != 0) {
+      stop("'SearchRanges' must not include NA values")
+    }
+    if (nrow(SearchRanges) != 2) {
+      stop("'SearchRanges' must have 2 rows")
+    }
+    if (ncol(SearchRanges) != NParam) {
+      stop("Incompatibility between 'SearchRanges' ncol and 'FUN_MOD'")
+    }
+  }
+
+  ## check StartParamList and StartParamDistrib default values
+  if (("HBAN"  %in% ObjectClass & is.null(StartParamList) & is.null(StartParamDistrib))) {
+    if ("GR4H" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(+5.12, -1.18, +4.34, -9.69,
+                         +5.58, -0.85, +4.74, -9.47,
+                         +6.01, -0.50, +5.14, -8.87), ncol = 4, byrow = TRUE)
+    }
+    if (("GR5H" == FeatFUN_MOD$CodeMod) & ("interception" %in% ObjectClass)) {
+      ParamT <- matrix(c(+3.46, -1.25, +4.04, -9.53, -9.34,
+                         +3.74, -0.41, +4.78, -8.94, -3.33,
+                         +4.29, +0.16, +5.39, -7.39, +3.33), ncol = 5, byrow = TRUE)
+    }
+    if (("GR5H" == FeatFUN_MOD$CodeMod) & !("interception" %in% ObjectClass)) {
+      ParamT <- matrix(c(+3.28, -0.39, +4.14, -9.54, -7.49,
+                         +3.62, -0.19, +4.80, -9.00, -6.31,
+                         +4.01, -0.04, +5.43, -7.53, -5.33), ncol = 5, byrow = TRUE)
+    }
+    if ("GR4J" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(+5.13, -1.60, +3.03, -9.05,
+                         +5.51, -0.61, +3.74, -8.51,
+                         +6.07, -0.02, +4.42, -8.06), ncol = 4, byrow = TRUE)
+    }
+    if ("GR5J" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(+5.17, -1.13, +3.08, -9.37, -7.45,
+                         +5.55, -0.46, +3.75, -9.09, -4.69,
+                         +6.10, -0.11, +4.43, -8.60, -0.66), ncol = 5, byrow = TRUE)
+
+    }
+    if ("GR6J" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(+3.60, -1.00, +3.30, -9.10, -0.90, +3.00,
+                         +3.90, -0.50, +4.10, -8.70, +0.10, +4.00,
+                         +4.50, +0.50, +5.00, -8.10, +1.10, +5.00), ncol = 6, byrow = TRUE)
+    }
+    if ("GR2M" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(+5.03, -7.15,
+                         +5.22, -6.74,
+                         +5.85, -6.37), ncol = 2, byrow = TRUE)
+    }
+    if ("GR1A" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(-1.69,
+                         -0.38,
+                         +1.39), ncol = 1, byrow = TRUE)
+    }
+
+
+    if ("CemaNeige" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(-9.96, +6.63,
+                         -9.14, +6.90,
+                         +4.10, +7.21), ncol = 2, byrow = TRUE)
+    }
+    if ("CemaNeigeGR4H" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(+5.12, -1.18, +4.34, -9.69, -9.96, +6.63,
+                         +5.58, -0.85, +4.74, -9.47, -9.14, +6.90,
+                         +6.01, -0.50, +5.14, -8.87, +4.10, +7.21), ncol = 6, byrow = TRUE)
+    }
+    if (("CemaNeigeGR5H" == FeatFUN_MOD$CodeMod) & ("interception" %in% ObjectClass)) {
+      ParamT <- matrix(c(+3.46, -1.25, +4.04, -9.53, -9.34, -9.96, +6.63,
+                         +3.74, -0.41, +4.78, -8.94, -3.33, -9.14, +6.90,
+                         +4.29, +0.16, +5.39, -7.39, +3.33, +4.10, +7.21), ncol = 7, byrow = TRUE)
+    }
+    if (("CemaNeigeGR5H" == FeatFUN_MOD$CodeMod) & !("interception" %in% ObjectClass)) {
+      ParamT <- matrix(c(+3.28, -0.39, +4.14, -9.54, -7.49, -9.96, +6.63,
+                         +3.62, -0.19, +4.80, -9.00, -6.31, -9.14, +6.90,
+                         +4.01, -0.04, +5.43, -7.53, -5.33, +4.10, +7.21), ncol = 7, byrow = TRUE)
+    }
+    if ("CemaNeigeGR4J" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(+5.13, -1.60, +3.03, -9.05, -9.96, +6.63,
+                         +5.51, -0.61, +3.74, -8.51, -9.14, +6.90,
+                         +6.07, -0.02, +4.42, -8.06, +4.10, +7.21), ncol = 6, byrow = TRUE)
+    }
+    if ("CemaNeigeGR5J" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(+5.17, -1.13, +3.08, -9.37, -7.45, -9.96, +6.63,
+                         +5.55, -0.46, +3.75, -9.09, -4.69, -9.14, +6.90,
+                         +6.10, -0.11, +4.43, -8.60, -0.66, +4.10, +7.21), ncol = 7, byrow = TRUE)
+    }
+    if ("CemaNeigeGR6J" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(+3.60, -1.00, +3.30, -9.10, -0.90, +3.00, -9.96, +6.63,
+                         +3.90, -0.50, +4.10, -8.70, +0.10, +4.00, -9.14, +6.90,
+                         +4.50, +0.50, +5.00, -8.10, +1.10, +5.00, +4.10, +7.21), ncol = 8, byrow = TRUE)
+    }
+
+    if (IsHyst) {
+      ParamTHyst <- matrix(c(-7.00, -7.00,
+                             -0.00, -0.00,
+                             +7.00, +7.00), ncol = 2, byrow = TRUE)
+      ParamT <- cbind(ParamT, ParamTHyst)
+    }
+    if (IsSD) {
+      if (identical(RunModel_Lag, FUN_SD)){
+        ParamTSD <- matrix(c(+1.25,
+                             +2.50,
+                             +5.00), ncol = 1, byrow = TRUE)
+        ParamT <- cbind(ParamTSD, ParamT)
+      }
+    }
+    StartParamList    <- NULL
+    StartParamDistrib <- TransfoParam(ParamIn = ParamT, Direction = "TR", FUN_TRANSFO = FUN_TRANSFO)
+
+  }
+
+  ## check StartParamList and StartParamDistrib format
+  if ("HBAN" %in% ObjectClass & !is.null(StartParamList)) {
+    if (!is.matrix(StartParamList)) {
+      stop("'StartParamList' must be a matrix")
+    }
+    if (!is.numeric(StartParamList)) {
+      stop("'StartParamList' must be a matrix of numeric values")
+    }
+    if (sum(is.na(StartParamList)) != 0) {
+      stop("'StartParamList' must not include NA values")
+    }
+    if (ncol(StartParamList) != NParam) {
+      stop("Incompatibility between 'StartParamList' ncol and 'FUN_MOD'")
+    }
+  }
+  if ("HBAN" %in% ObjectClass & !is.null(StartParamDistrib)) {
+    if (!is.matrix(StartParamDistrib)) {
+      stop("'StartParamDistrib' must be a matrix")
+    }
+    if (!is.numeric(StartParamDistrib[1, ])) {
+      stop("'StartParamDistrib' must be a matrix of numeric values")
+    }
+    if (sum(is.na(StartParamDistrib[1, ])) != 0) {
+      stop("'StartParamDistrib' must not include NA values on the first line")
+    }
+    if (ncol(StartParamDistrib) != NParam) {
+      stop("Incompatibility between 'StartParamDistrib' ncol and 'FUN_MOD'")
+    }
+  }
+
+
+  ##Create_CalibOptions
+  CalibOptions <- list(FixedParam = FixedParam, SearchRanges = SearchRanges, FUN_TRANSFO = FUN_TRANSFO)
+
+  if (!is.null(StartParamList)) {
+    CalibOptions <- c(CalibOptions, list(StartParamList = StartParamList))
+  }
+  if (!is.null(StartParamDistrib)) {
+    CalibOptions <- c(CalibOptions, list(StartParamDistrib = StartParamDistrib))
+  }
+  class(CalibOptions) <- c("CalibOptions", ObjectClass)
+
+  return(CalibOptions)
+
+}
diff --git a/R/CreateRunOptions.R b/R/CreateRunOptions.R
index fd1694e271d03e1af37623eaf1f31dc2a8bb8a20..76daea77263c05c010abdbc74805b47120bf260a 100644
--- a/R/CreateRunOptions.R
+++ b/R/CreateRunOptions.R
@@ -1,482 +1,483 @@
-CreateRunOptions <- function(FUN_MOD, InputsModel,
-                             IndPeriod_WarmUp = NULL, IndPeriod_Run,
-                             IniStates = NULL, IniResLevels = NULL, Imax = NULL,
-                             Outputs_Cal = NULL, Outputs_Sim = "all",
-                             MeanAnSolidPrecip = NULL, IsHyst = FALSE,
-                             warnings = TRUE, verbose = TRUE) {
-
-  if (!is.null(Imax)) {
-    if (!is.numeric(Imax) | length(Imax) != 1L) {
-      stop("'Imax' must be a non negative 'numeric' value of length 1")
-    } else {
-      if (Imax < 0) {
-        stop("'Imax' must be a non negative 'numeric' value of length 1")
-      }
-    }
-    IsIntStore <- TRUE
-  } else {
-    IsIntStore <- FALSE
-  }
-
-  ## check FUN_MOD
-  FUN_MOD <- match.fun(FUN_MOD)
-  FeatFUN_MOD <- .GetFeatModel(FUN_MOD = FUN_MOD, DatesR = InputsModel$DatesR)
-  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 ("CemaNeige" %in% FeatFUN_MOD$Class) {
-    FeatFUN_MOD$IsHyst <- IsHyst
-    if (IsHyst) {
-      ObjectClass <- c(ObjectClass, "hysteresis")
-      FeatFUN_MOD$NbParam <- FeatFUN_MOD$NbParam + 2
-    }
-  }
-
-  ## SD model
-  FeatFUN_MOD$IsSD <- inherits(InputsModel, "SD")
-  if (FeatFUN_MOD$IsSD) {
-    FeatFUN_MOD$NbParam <- FeatFUN_MOD$NbParam + 1
-  }
-
-  if (!"CemaNeige" %in% ObjectClass & "hysteresis" %in% ObjectClass) {
-    stop("'IsHyst' cannot be TRUE for the chosen 'FUN_MOD'")
-  }
-  if (!(identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) & "interception" %in% ObjectClass) {
-    stop("'IMax' cannot be set for the chosen 'FUN_MOD'")
-  }
-
-  ##check_InputsModel
-  if (!inherits(InputsModel, "InputsModel")) {
-    stop("'InputsModel' must be of class 'InputsModel'")
-  }
-  if ("GR" %in% ObjectClass & !inherits(InputsModel, "GR")) {
-    stop("'InputsModel' must be of class 'GR'")
-  }
-  if ("CemaNeige" %in% ObjectClass &
-      !inherits(InputsModel, "CemaNeige")) {
-    stop("'InputsModel' must be of class 'CemaNeige'")
-  }
-  if ("hourly" %in% ObjectClass &
-      !inherits(InputsModel, "hourly")) {
-    stop("'InputsModel' must be of class 'hourly'")
-  }
-  if ("daily" %in% ObjectClass & !inherits(InputsModel, "daily")) {
-    stop("'InputsModel' must be of class 'daily'")
-  }
-  if ("monthly" %in% ObjectClass &
-      !inherits(InputsModel, "monthly")) {
-    stop("'InputsModel' must be of class 'monthly'")
-  }
-  if ("yearly" %in% ObjectClass &
-      !inherits(InputsModel, "yearly")) {
-    stop("'InputsModel' must be of class 'yearly'")
-  }
-
-
-  ##check_IndPeriod_Run
-  if (!is.vector(IndPeriod_Run)) {
-    stop("'IndPeriod_Run' must be a vector of numeric values")
-  }
-  if (!is.numeric(IndPeriod_Run)) {
-    stop("'IndPeriod_Run' must be a vector of numeric values")
-  }
-  if (!identical(as.integer(IndPeriod_Run), as.integer(seq(from = IndPeriod_Run[1], to = tail(IndPeriod_Run, 1), by = 1)))) {
-    stop("'IndPeriod_Run' must be a continuous sequence of integers")
-  }
-  if (storage.mode(IndPeriod_Run) != "integer") {
-    stop("'IndPeriod_Run' should be of type integer")
-  }
-
-
-  ##check_IndPeriod_WarmUp
-  WTxt <- NULL
-  if (is.null(IndPeriod_WarmUp)) {
-    WTxt <- paste(WTxt, "model warm up period not defined: default configuration used", sep = "")
-    ##If_the_run_period_starts_at_the_very_beginning_of_the_time_series
-    if (IndPeriod_Run[1L] == 1L) {
-      IndPeriod_WarmUp <- 0L
-      WTxt <- paste0(WTxt,"\n no data were found for model warm up!")
-      ##We_look_for_the_longest_period_preceeding_the_run_period_with_a_maximum_of_one_year
-    } else {
-      TmpDateR0 <- InputsModel$DatesR[IndPeriod_Run[1]]
-      TmpDateR  <- TmpDateR0 - 365 * 24 * 60 * 60
-      ### minimal date to start the warmup
-      if (format(TmpDateR, format = "%d") != format(TmpDateR0, format = "%d")) {
-        ### leap year
-        TmpDateR <- TmpDateR - 1 * 24 * 60 * 60
-      }
-      IndPeriod_WarmUp <- which(InputsModel$DatesR == max(InputsModel$DatesR[1], TmpDateR)):(IndPeriod_Run[1] - 1)
-      if (length(IndPeriod_WarmUp) * TimeStepMean / (365 * 24 * 60 * 60) >= 1) {
-        WTxt <- paste0(WTxt, "\n  the year preceding the run period is used \n")
-      } else {
-        WTxt <- paste0(WTxt, "\n  less than a year (without missing values) was found for model warm up:")
-        WTxt <- paste0(WTxt, "\n  (", length(IndPeriod_WarmUp), " time-steps are used for initialisation)")
-      }
-    }
-  }
-  if (!is.null(IndPeriod_WarmUp)) {
-    if (!is.vector(IndPeriod_WarmUp)) {
-      stop("'IndPeriod_WarmUp' must be a vector of numeric values")
-    }
-    if (!is.numeric(IndPeriod_WarmUp)) {
-      stop("'IndPeriod_WarmUp' must be a vector of numeric values")
-    }
-    if (storage.mode(IndPeriod_WarmUp) != "integer") {
-      stop("'IndPeriod_WarmUp' should be of type integer")
-    }
-    if (identical(IndPeriod_WarmUp, 0L) & verbose) {
-      message(paste0(WTxt, " No warm up period is used"))
-    }
-    if ((IndPeriod_Run[1] - 1) != tail(IndPeriod_WarmUp, 1) & !identical(IndPeriod_WarmUp, 0L)) {
-      WTxt <- paste0(WTxt, " Model warm up period is not directly before the model run period")
-    }
-  }
-  if (!is.null(WTxt) & warnings) {
-    warning(WTxt)
-  }
-
-
-  ## check IniResLevels
-  if ("GR" %in% ObjectClass & ("monthly" %in% ObjectClass | "daily" %in% ObjectClass | "hourly" %in% ObjectClass)) {
-    if (!is.null(IniResLevels)) {
-      # if (!is.vector(IniResLevels) | !is.numeric(IniResLevels) | any(is.na(IniResLevels))) {
-      if (!is.vector(IniResLevels) | is.character(IniResLevels) | is.factor(IniResLevels) | length(IniResLevels) != 4) {
-        stop("'IniResLevels' must be a vector of 4 numeric values")
-      }
-      # if ((identical(FUN_MOD, RunModel_GR4H) | identical(FUN_MOD, RunModel_CemaNeigeGR4H) |
-      #      # (identical(FUN_MOD, RunModel_GR5H) & !IsIntStore) |
-      #      identical(FUN_MOD, RunModel_GR5H) |
-      #      identical(FUN_MOD, RunModel_GR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR4J) |
-      #      identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) |
-      #      identical(FUN_MOD, RunModel_GR2M)) &
-      #     length(IniResLevels) != 2) {
-      #   stop("the length of 'IniResLevels' must be 2 for the chosen 'FUN_MOD'")
-      # }
-      if (any(is.na(IniResLevels[1:2]))) {
-        stop("the first 2 values of 'IniResLevels' cannot be missing values for the chosen 'FUN_MOD'")
-      }
-      # if ((identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J) |
-      #      (identical(FUN_MOD, RunModel_GR5H) & IsIntStore)) &
-      #     length(IniResLevels) != 3) {
-      #   stop("the length of 'IniResLevels' must be 3 for the chosen 'FUN_MOD'")
-      # }
-      if ((identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J))) {
-        if (is.na(IniResLevels[3L])) {
-          stop("the third value of 'IniResLevels' cannot be a missing value for the chosen 'FUN_MOD'")
-        }
-      } else {
-        if (!is.na(IniResLevels[3L])) {
-          warning("the third value of 'IniResLevels' is set to NA value for the chosen 'FUN_MOD'. Only GR6J presents an exponential store")
-          IniResLevels[3L] <- NA
-        }
-      }
-      if (identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) {
-        if (IsIntStore & is.na(IniResLevels[4L])) {
-          stop("the fourth value of 'IniResLevels' cannot be a missing value for the chosen 'FUN_MOD' (GR5H with an interception store)")
-        }
-        if (!IsIntStore & !is.na(IniResLevels[4L])) {
-          warning("the fourth value of 'IniResLevels' is set to NA value for the chosen 'FUN_MOD'. Only GR5H used with an 'Imax' value presents an interception store")
-          IniResLevels[4L] <- NA
-        }
-      } else {
-        if (!is.na(IniResLevels[4L])) {
-          warning("the fourth value of 'IniResLevels' is set to NA value for the chosen 'FUN_MOD'. Only GR5H used with an 'Imax' value presents an interception store")
-          IniResLevels[4L] <- NA
-        }
-      }
-    } else if (is.null(IniStates)) {
-      IniResLevels <- as.double(c(0.3, 0.5, NA, NA))
-      if (identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
-        IniResLevels <- as.double(c(0.3, 0.5, 0, NA))
-      }
-      if ((identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) & IsIntStore) {
-        IniResLevels <- as.double(c(0.3, 0.5, NA, 0))
-      }
-      # if (!identical(FUN_MOD, RunModel_GR6J) & !identical(FUN_MOD, RunModel_CemaNeigeGR6J) &
-      #     !identical(FUN_MOD, RunModel_GR5H) & !identical(FUN_MOD, RunModel_CemaNeigeGR5H)) {
-      # if (is.null(IniStates)) {
-      #   IniResLevels <- as.double(c(0.3, 0.5, NA, NA))
-      # }
-    }
-  } else {
-    if (!is.null(IniResLevels)) {
-      stop("'IniResLevels' can only be used with monthly or daily or hourly GR models")
-    }
-  }
-  ## check IniStates
-  if (is.null(IniStates) & is.null(IniResLevels) & warnings) {
-    warning("model states initialisation not defined: default configuration used")
-  }
-  if (!is.null(IniStates) & !is.null(IniResLevels) & warnings) {
-    warning("'IniStates' and 'IniResLevels' are both defined: store levels are taken from 'IniResLevels'")
-  }
-  if ("CemaNeige" %in% ObjectClass) {
-    NLayers <- length(InputsModel$LayerPrecip)
-  } else {
-    NLayers <- 0
-  }
-  NState <- NULL
-  if ("GR" %in% ObjectClass | "CemaNeige" %in% ObjectClass) {
-    if ("hourly"  %in% ObjectClass) {
-      NState <- 7 + 3 * 24 * 20 + 4 * NLayers
-    }
-    if ("daily"   %in% ObjectClass) {
-      NState <- 7 + 3 * 20 + 4 * NLayers
-    }
-    if ("monthly" %in% ObjectClass) {
-      NState <- 2
-    }
-    if ("yearly"  %in% ObjectClass) {
-      NState <- 1
-    }
-  }
-  if (!is.null(IniStates)) {
-
-    if (!inherits(IniStates, "IniStates")) {
-      stop("'IniStates' must be an object of class 'IniStates'")
-    }
-    if (sum(ObjectClass %in% class(IniStates)) < 2) {
-      stop(paste0("non convenient 'IniStates' for the chosen 'FUN_MOD'"))
-    }
-    if (identical(FUN_MOD, RunModel_GR1A) & !is.null(IniStates)) { ## GR1A
-      stop(paste0("'IniStates' is not available for the chosen 'FUN_MOD'"))
-    }
-    if ((identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) |
-         identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) &
-        !all(is.na(IniStates$UH$UH1))) { ## GR5J or GR5H
-      stop(paste0("non convenient 'IniStates' for the chosen 'FUN_MOD'.' In 'IniStates', 'UH1' has to be a vector of NA for GR5J"))
-    }
-    if ((identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) & is.na(IniStates$Store$Exp)) { ## GR6J
-      stop(paste0("non convenient 'IniStates' for the chosen 'FUN_MOD'.' GR6J needs an exponential store value in 'IniStates'"))
-    }
-    if ((identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) & is.na(IniStates$Store$Int)) { ## GR5H interception
-
-      stop(paste0("non convenient 'IniStates' for the chosen 'FUN_MOD'.' GR5H (with interception store) needs an interception store value in 'IniStates'"))
-    }
-    if (!(identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) & !is.na(IniStates$Store$Exp)) { ## except GR6J
-      stop(paste0("non convenient 'IniStates' for the chosen 'FUN_MOD'.' No exponential store value needed in 'IniStates'"))
-    }
-    if (!(identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) & !is.na(IniStates$Store$Int)) { ## except GR5H interception
-      stop(paste0("non convenient 'IniStates' for the chosen 'FUN_MOD'.' No interception store value needed in 'IniStates'"))
-    }
-    # if (length(na.omit(unlist(IniStates))) != NState) {
-    #   stop(paste0("the length of IniStates must be ", NState, " for the chosen FUN_MOD"))
-    # }
-    if ((!"CemaNeige" %in% ObjectClass &  inherits(IniStates, "CemaNeige")) |
-        ( "CemaNeige" %in% ObjectClass & !inherits(IniStates, "CemaNeige"))) {
-      stop("'FUN_MOD' and 'IniStates' must be both of class 'CemaNeige'")
-    }
-    if (( "hysteresis" %in% ObjectClass & !inherits(IniStates, "hysteresis")) |
-        (!"hysteresis" %in% ObjectClass &  inherits(IniStates, "hysteresis"))) {
-      stop("'IsHyst' and 'IniStates' are not consistent on the use of the hysteresis")
-    }
-    if (!"CemaNeige" %in% ObjectClass & any(is.na(IniStates$CemaNeigeLayers$G  ))) {
-      IniStates$CemaNeigeLayers$G <- NULL
-    }
-    if (!"CemaNeige" %in% ObjectClass & any(is.na(IniStates$CemaNeigeLayers$eTG))) {
-      IniStates$CemaNeigeLayers$eTG <- NULL
-    }
-    if (!"CemaNeige" %in% ObjectClass & any(is.na(IniStates$CemaNeigeLayers$Gthr))) {
-      IniStates$CemaNeigeLayers$Gthr <- NULL
-    }
-    if (!"CemaNeige" %in% ObjectClass & any(is.na(IniStates$CemaNeigeLayers$Glocmax))) {
-      IniStates$CemaNeigeLayers$Glocmax <- NULL
-    }
-    IniStates$Store$Rest <- rep(NA, 3)
-    IniStates <- unlist(IniStates)
-    IniStates[is.na(IniStates)] <- 0
-    if ("monthly" %in% ObjectClass) {
-      IniStates <- IniStates[seq_len(NState)]
-    }
-  } else {
-    IniStates <- as.double(rep(0.0, NState))
-  }
-
-
-  ##check_Outputs_Cal_and_Sim
-
-  ##Outputs_all
-  Outputs_all <- c("DatesR", unlist(FortranOutputs), "WarmUpQsim", "StateEnd", "Param")
-  if (FeatFUN_MOD$IsSD) {
-    Outputs_all <- c(Outputs_all, "QsimDown", "Qsim_m3")
-  }
-
-  ##check_Outputs_Sim
-  if (!is.vector(Outputs_Sim)) {
-    stop("'Outputs_Sim' must be a vector of characters")
-  }
-  if (!is.character(Outputs_Sim)) {
-    stop("'Outputs_Sim' must be a vector of characters")
-  }
-  if (sum(is.na(Outputs_Sim)) != 0) {
-    stop("'Outputs_Sim' must not contain NA")
-  }
-  if ("all" %in% Outputs_Sim) {
-    Outputs_Sim <- Outputs_all
-  }
-  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"))
-  }
-  Outputs_Sim <- Outputs_Sim[!duplicated(Outputs_Sim)]
-
-
-  ##check_Outputs_Cal
-  if (is.null(Outputs_Cal)) {
-    if ("GR" %in% ObjectClass) {
-      Outputs_Cal <- c("Qsim", "Param")
-      if ("CemaNeige" %in% ObjectClass) {
-        Outputs_Cal <- c("PliqAndMelt", Outputs_Cal)
-      }
-    } else if ("CemaNeige" %in% ObjectClass) {
-      Outputs_Cal <- c("all")
-    }
-  } else {
-    if (!is.vector(Outputs_Cal)) {
-      stop("'Outputs_Cal' must be a vector of characters")
-    }
-    if (!is.character(Outputs_Cal)) {
-      stop("'Outputs_Cal' must be a vector of characters")
-    }
-    if (sum(is.na(Outputs_Cal)) != 0) {
-      stop("'Outputs_Cal' must not contain NA")
-    }
-  }
-  if ("all" %in% Outputs_Cal) {
-    Outputs_Cal <- Outputs_all
-  }
-  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"))
-  }
-  Outputs_Cal <- unique(Outputs_Cal)
-
-
-  ##check_MeanAnSolidPrecip
-  if ("CemaNeige" %in% ObjectClass & is.null(MeanAnSolidPrecip)) {
-    NLayers <- length(InputsModel$LayerPrecip)
-    SolidPrecip <- NULL
-    for (iLayer in 1:NLayers) {
-      if (iLayer == 1) {
-        SolidPrecip <-
-          InputsModel$LayerFracSolidPrecip[[1]] * InputsModel$LayerPrecip[[iLayer]] /
-          NLayers
-      } else {
-        SolidPrecip <- SolidPrecip + InputsModel$LayerFracSolidPrecip[[iLayer]] * InputsModel$LayerPrecip[[iLayer]] / NLayers
-      }
-    }
-    Factor <- NULL
-    if (inherits(InputsModel, "hourly")) {
-      Factor <- 365.25 * 24
-    }
-    if (inherits(InputsModel, "daily")) {
-      Factor <- 365.25
-    }
-    if (inherits(InputsModel, "monthly")) {
-      Factor <- 12
-    }
-    if (inherits(InputsModel, "yearly")) {
-      Factor <- 1
-    }
-    if (is.null(Factor)) {
-      stop("'InputsModel' must be of class 'hourly', 'daily', 'monthly' or 'yearly'")
-    }
-    MeanAnSolidPrecip <- rep(mean(SolidPrecip) * Factor, NLayers)
-    ### default value: same Gseuil for all layers
-    if (warnings) {
-      warning("'MeanAnSolidPrecip' not defined: it was automatically set to c(",
-              paste(round(MeanAnSolidPrecip), collapse = ","), ")  from the 'InputsModel' given to the function. ",
-              "Be careful in case your application is (short-term) forecasting.\n")
-    }
-  }
-  if ("CemaNeige" %in% ObjectClass & !is.null(MeanAnSolidPrecip)) {
-    if (!is.vector(MeanAnSolidPrecip)) {
-      stop(paste0("'MeanAnSolidPrecip' must be a vector of numeric values"))
-    }
-    if (!is.numeric(MeanAnSolidPrecip)) {
-      stop(paste0("'MeanAnSolidPrecip' must be a vector of numeric values"))
-    }
-    if (length(MeanAnSolidPrecip) != NLayers) {
-      stop(paste0("'MeanAnSolidPrecip' must be a numeric vector of length ", NLayers, ""))
-    }
-  }
-
-
-  ##check_PliqAndMelt
-  if ("GR" %in% ObjectClass & "CemaNeige" %in% ObjectClass) {
-    if (!"PliqAndMelt" %in% Outputs_Cal & !"all" %in% Outputs_Cal) {
-      WTxt <- NULL
-      WTxt <- paste0(WTxt, "'PliqAndMelt' was not defined in 'Outputs_Cal' but is needed to feed the hydrological model with the snow modele outputs \n")
-      WTxt <- paste0(WTxt, ": it was automatically added \n")
-      if (!is.null(WTxt) & warnings) {
-        warning(WTxt)
-      }
-      Outputs_Cal <- c(Outputs_Cal, "PliqAndMelt")
-    }
-    if (!"PliqAndMelt" %in% Outputs_Sim & !"all" %in% Outputs_Sim) {
-      WTxt <- NULL
-      WTxt <- paste0(WTxt, "'PliqAndMelt' was not defined in 'Outputs_Sim' but is needed to feed the hydrological model with the snow modele outputs \n")
-      WTxt <- paste0(WTxt, ": it was automatically added \n")
-      if (!is.null(WTxt) & warnings) {
-        warning(WTxt)
-      }
-      Outputs_Sim <- c(Outputs_Sim, "PliqAndMelt")
-    }
-  }
-
-
-  ##check_Qsim
-  if ("GR" %in% ObjectClass) {
-    if (!"Qsim" %in% Outputs_Cal & !"all" %in% Outputs_Cal) {
-      WTxt <- NULL
-      WTxt <- paste0(WTxt, "'Qsim' was not defined in 'Outputs_Cal' \n")
-      WTxt <- paste0(WTxt, ": it was automatically added \n")
-      if (!is.null(WTxt) & warnings) {
-        warning(WTxt)
-      }
-      Outputs_Cal <- c(Outputs_Cal, "Qsim")
-    }
-    if (!"Qsim" %in% Outputs_Sim & !"all" %in% Outputs_Sim) {
-      WTxt <- NULL
-      WTxt <- paste0(WTxt, "'Qsim' was not defined in 'Outputs_Sim' \n")
-      WTxt <- paste0(WTxt, ": it was automatically added \n")
-      if (!is.null(WTxt) & warnings) {
-        warning(WTxt)
-      }
-      Outputs_Sim <- c(Outputs_Sim, "Qsim")
-    }
-  }
-
-  ##Create_RunOptions
-  RunOptions <- list(IndPeriod_WarmUp = IndPeriod_WarmUp,
-                     IndPeriod_Run = IndPeriod_Run,
-                     IniStates = IniStates,
-                     IniResLevels = IniResLevels,
-                     Outputs_Cal = Outputs_Cal,
-                     Outputs_Sim = Outputs_Sim,
-                     FortranOutputs = FortranOutputs,
-                     FeatFUN_MOD = FeatFUN_MOD)
-
-  if ("CemaNeige" %in% ObjectClass) {
-    RunOptions <- c(RunOptions, list(MeanAnSolidPrecip = MeanAnSolidPrecip))
-  }
-  if ("interception" %in% ObjectClass) {
-    RunOptions <- c(RunOptions, list(Imax = Imax))
-  }
-  class(RunOptions) <- c("RunOptions", ObjectClass)
-
-  return(RunOptions)
-
-
-}
-
+CreateRunOptions <- function(FUN_MOD, InputsModel,
+                             IndPeriod_WarmUp = NULL, IndPeriod_Run,
+                             IniStates = NULL, IniResLevels = NULL, Imax = NULL,
+                             Outputs_Cal = NULL, Outputs_Sim = "all",
+                             MeanAnSolidPrecip = NULL, FUN_SD = NULL, IsHyst = FALSE,
+                             warnings = TRUE, verbose = TRUE) {
+
+  if (!is.null(Imax)) {
+    if (!is.numeric(Imax) | length(Imax) != 1L) {
+      stop("'Imax' must be a non negative 'numeric' value of length 1")
+    } else {
+      if (Imax < 0) {
+        stop("'Imax' must be a non negative 'numeric' value of length 1")
+      }
+    }
+    IsIntStore <- TRUE
+  } else {
+    IsIntStore <- FALSE
+  }
+
+  ## check FUN_MOD
+  FUN_MOD <- match.fun(FUN_MOD)
+  FeatFUN_MOD <- .GetFeatModel(FUN_MOD = FUN_MOD, DatesR = InputsModel$DatesR,
+                               FUN_SD = FUN_SD)
+  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 ("CemaNeige" %in% FeatFUN_MOD$Class) {
+    FeatFUN_MOD$IsHyst <- IsHyst
+    if (IsHyst) {
+      ObjectClass <- c(ObjectClass, "hysteresis")
+      FeatFUN_MOD$NbParam <- FeatFUN_MOD$NbParam + 2
+    }
+  }
+
+  ## SD model
+  FeatFUN_MOD$IsSD <- inherits(InputsModel, "SD")
+  if (FeatFUN_MOD$IsSD) {
+    FeatFUN_MOD$NbParam <- FeatFUN_MOD$NbParam + FeatFUN_MOD$ParamSD
+  }
+
+  if (!"CemaNeige" %in% ObjectClass & "hysteresis" %in% ObjectClass) {
+    stop("'IsHyst' cannot be TRUE for the chosen 'FUN_MOD'")
+  }
+  if (!(identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) & "interception" %in% ObjectClass) {
+    stop("'IMax' cannot be set for the chosen 'FUN_MOD'")
+  }
+
+  ##check_InputsModel
+  if (!inherits(InputsModel, "InputsModel")) {
+    stop("'InputsModel' must be of class 'InputsModel'")
+  }
+  if ("GR" %in% ObjectClass & !inherits(InputsModel, "GR")) {
+    stop("'InputsModel' must be of class 'GR'")
+  }
+  if ("CemaNeige" %in% ObjectClass &
+      !inherits(InputsModel, "CemaNeige")) {
+    stop("'InputsModel' must be of class 'CemaNeige'")
+  }
+  if ("hourly" %in% ObjectClass &
+      !inherits(InputsModel, "hourly")) {
+    stop("'InputsModel' must be of class 'hourly'")
+  }
+  if ("daily" %in% ObjectClass & !inherits(InputsModel, "daily")) {
+    stop("'InputsModel' must be of class 'daily'")
+  }
+  if ("monthly" %in% ObjectClass &
+      !inherits(InputsModel, "monthly")) {
+    stop("'InputsModel' must be of class 'monthly'")
+  }
+  if ("yearly" %in% ObjectClass &
+      !inherits(InputsModel, "yearly")) {
+    stop("'InputsModel' must be of class 'yearly'")
+  }
+
+
+  ##check_IndPeriod_Run
+  if (!is.vector(IndPeriod_Run)) {
+    stop("'IndPeriod_Run' must be a vector of numeric values")
+  }
+  if (!is.numeric(IndPeriod_Run)) {
+    stop("'IndPeriod_Run' must be a vector of numeric values")
+  }
+  if (!identical(as.integer(IndPeriod_Run), as.integer(seq(from = IndPeriod_Run[1], to = tail(IndPeriod_Run, 1), by = 1)))) {
+    stop("'IndPeriod_Run' must be a continuous sequence of integers")
+  }
+  if (storage.mode(IndPeriod_Run) != "integer") {
+    stop("'IndPeriod_Run' should be of type integer")
+  }
+
+
+  ##check_IndPeriod_WarmUp
+  WTxt <- NULL
+  if (is.null(IndPeriod_WarmUp)) {
+    WTxt <- paste(WTxt, "model warm up period not defined: default configuration used", sep = "")
+    ##If_the_run_period_starts_at_the_very_beginning_of_the_time_series
+    if (IndPeriod_Run[1L] == 1L) {
+      IndPeriod_WarmUp <- 0L
+      WTxt <- paste0(WTxt,"\n no data were found for model warm up!")
+      ##We_look_for_the_longest_period_preceeding_the_run_period_with_a_maximum_of_one_year
+    } else {
+      TmpDateR0 <- InputsModel$DatesR[IndPeriod_Run[1]]
+      TmpDateR  <- TmpDateR0 - 365 * 24 * 60 * 60
+      ### minimal date to start the warmup
+      if (format(TmpDateR, format = "%d") != format(TmpDateR0, format = "%d")) {
+        ### leap year
+        TmpDateR <- TmpDateR - 1 * 24 * 60 * 60
+      }
+      IndPeriod_WarmUp <- which(InputsModel$DatesR == max(InputsModel$DatesR[1], TmpDateR)):(IndPeriod_Run[1] - 1)
+      if (length(IndPeriod_WarmUp) * TimeStepMean / (365 * 24 * 60 * 60) >= 1) {
+        WTxt <- paste0(WTxt, "\n  the year preceding the run period is used \n")
+      } else {
+        WTxt <- paste0(WTxt, "\n  less than a year (without missing values) was found for model warm up:")
+        WTxt <- paste0(WTxt, "\n  (", length(IndPeriod_WarmUp), " time-steps are used for initialisation)")
+      }
+    }
+  }
+  if (!is.null(IndPeriod_WarmUp)) {
+    if (!is.vector(IndPeriod_WarmUp)) {
+      stop("'IndPeriod_WarmUp' must be a vector of numeric values")
+    }
+    if (!is.numeric(IndPeriod_WarmUp)) {
+      stop("'IndPeriod_WarmUp' must be a vector of numeric values")
+    }
+    if (storage.mode(IndPeriod_WarmUp) != "integer") {
+      stop("'IndPeriod_WarmUp' should be of type integer")
+    }
+    if (identical(IndPeriod_WarmUp, 0L) & verbose) {
+      message(paste0(WTxt, " No warm up period is used"))
+    }
+    if ((IndPeriod_Run[1] - 1) != tail(IndPeriod_WarmUp, 1) & !identical(IndPeriod_WarmUp, 0L)) {
+      WTxt <- paste0(WTxt, " Model warm up period is not directly before the model run period")
+    }
+  }
+  if (!is.null(WTxt) & warnings) {
+    warning(WTxt)
+  }
+
+
+  ## check IniResLevels
+  if ("GR" %in% ObjectClass & ("monthly" %in% ObjectClass | "daily" %in% ObjectClass | "hourly" %in% ObjectClass)) {
+    if (!is.null(IniResLevels)) {
+      # if (!is.vector(IniResLevels) | !is.numeric(IniResLevels) | any(is.na(IniResLevels))) {
+      if (!is.vector(IniResLevels) | is.character(IniResLevels) | is.factor(IniResLevels) | length(IniResLevels) != 4) {
+        stop("'IniResLevels' must be a vector of 4 numeric values")
+      }
+      # if ((identical(FUN_MOD, RunModel_GR4H) | identical(FUN_MOD, RunModel_CemaNeigeGR4H) |
+      #      # (identical(FUN_MOD, RunModel_GR5H) & !IsIntStore) |
+      #      identical(FUN_MOD, RunModel_GR5H) |
+      #      identical(FUN_MOD, RunModel_GR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR4J) |
+      #      identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) |
+      #      identical(FUN_MOD, RunModel_GR2M)) &
+      #     length(IniResLevels) != 2) {
+      #   stop("the length of 'IniResLevels' must be 2 for the chosen 'FUN_MOD'")
+      # }
+      if (any(is.na(IniResLevels[1:2]))) {
+        stop("the first 2 values of 'IniResLevels' cannot be missing values for the chosen 'FUN_MOD'")
+      }
+      # if ((identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J) |
+      #      (identical(FUN_MOD, RunModel_GR5H) & IsIntStore)) &
+      #     length(IniResLevels) != 3) {
+      #   stop("the length of 'IniResLevels' must be 3 for the chosen 'FUN_MOD'")
+      # }
+      if ((identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J))) {
+        if (is.na(IniResLevels[3L])) {
+          stop("the third value of 'IniResLevels' cannot be a missing value for the chosen 'FUN_MOD'")
+        }
+      } else {
+        if (!is.na(IniResLevels[3L])) {
+          warning("the third value of 'IniResLevels' is set to NA value for the chosen 'FUN_MOD'. Only GR6J presents an exponential store")
+          IniResLevels[3L] <- NA
+        }
+      }
+      if (identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) {
+        if (IsIntStore & is.na(IniResLevels[4L])) {
+          stop("the fourth value of 'IniResLevels' cannot be a missing value for the chosen 'FUN_MOD' (GR5H with an interception store)")
+        }
+        if (!IsIntStore & !is.na(IniResLevels[4L])) {
+          warning("the fourth value of 'IniResLevels' is set to NA value for the chosen 'FUN_MOD'. Only GR5H used with an 'Imax' value presents an interception store")
+          IniResLevels[4L] <- NA
+        }
+      } else {
+        if (!is.na(IniResLevels[4L])) {
+          warning("the fourth value of 'IniResLevels' is set to NA value for the chosen 'FUN_MOD'. Only GR5H used with an 'Imax' value presents an interception store")
+          IniResLevels[4L] <- NA
+        }
+      }
+    } else if (is.null(IniStates)) {
+      IniResLevels <- as.double(c(0.3, 0.5, NA, NA))
+      if (identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
+        IniResLevels <- as.double(c(0.3, 0.5, 0, NA))
+      }
+      if ((identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) & IsIntStore) {
+        IniResLevels <- as.double(c(0.3, 0.5, NA, 0))
+      }
+      # if (!identical(FUN_MOD, RunModel_GR6J) & !identical(FUN_MOD, RunModel_CemaNeigeGR6J) &
+      #     !identical(FUN_MOD, RunModel_GR5H) & !identical(FUN_MOD, RunModel_CemaNeigeGR5H)) {
+      # if (is.null(IniStates)) {
+      #   IniResLevels <- as.double(c(0.3, 0.5, NA, NA))
+      # }
+    }
+  } else {
+    if (!is.null(IniResLevels)) {
+      stop("'IniResLevels' can only be used with monthly or daily or hourly GR models")
+    }
+  }
+  ## check IniStates
+  if (is.null(IniStates) & is.null(IniResLevels) & warnings) {
+    warning("model states initialisation not defined: default configuration used")
+  }
+  if (!is.null(IniStates) & !is.null(IniResLevels) & warnings) {
+    warning("'IniStates' and 'IniResLevels' are both defined: store levels are taken from 'IniResLevels'")
+  }
+  if ("CemaNeige" %in% ObjectClass) {
+    NLayers <- length(InputsModel$LayerPrecip)
+  } else {
+    NLayers <- 0
+  }
+  NState <- NULL
+  if ("GR" %in% ObjectClass | "CemaNeige" %in% ObjectClass) {
+    if ("hourly"  %in% ObjectClass) {
+      NState <- 7 + 3 * 24 * 20 + 4 * NLayers
+    }
+    if ("daily"   %in% ObjectClass) {
+      NState <- 7 + 3 * 20 + 4 * NLayers
+    }
+    if ("monthly" %in% ObjectClass) {
+      NState <- 2
+    }
+    if ("yearly"  %in% ObjectClass) {
+      NState <- 1
+    }
+  }
+  if (!is.null(IniStates)) {
+
+    if (!inherits(IniStates, "IniStates")) {
+      stop("'IniStates' must be an object of class 'IniStates'")
+    }
+    if (sum(ObjectClass %in% class(IniStates)) < 2) {
+      stop(paste0("non convenient 'IniStates' for the chosen 'FUN_MOD'"))
+    }
+    if (identical(FUN_MOD, RunModel_GR1A) & !is.null(IniStates)) { ## GR1A
+      stop(paste0("'IniStates' is not available for the chosen 'FUN_MOD'"))
+    }
+    if ((identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) |
+         identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) &
+        !all(is.na(IniStates$UH$UH1))) { ## GR5J or GR5H
+      stop(paste0("non convenient 'IniStates' for the chosen 'FUN_MOD'.' In 'IniStates', 'UH1' has to be a vector of NA for GR5J"))
+    }
+    if ((identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) & is.na(IniStates$Store$Exp)) { ## GR6J
+      stop(paste0("non convenient 'IniStates' for the chosen 'FUN_MOD'.' GR6J needs an exponential store value in 'IniStates'"))
+    }
+    if ((identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) & is.na(IniStates$Store$Int)) { ## GR5H interception
+
+      stop(paste0("non convenient 'IniStates' for the chosen 'FUN_MOD'.' GR5H (with interception store) needs an interception store value in 'IniStates'"))
+    }
+    if (!(identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) & !is.na(IniStates$Store$Exp)) { ## except GR6J
+      stop(paste0("non convenient 'IniStates' for the chosen 'FUN_MOD'.' No exponential store value needed in 'IniStates'"))
+    }
+    if (!(identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) & !is.na(IniStates$Store$Int)) { ## except GR5H interception
+      stop(paste0("non convenient 'IniStates' for the chosen 'FUN_MOD'.' No interception store value needed in 'IniStates'"))
+    }
+    # if (length(na.omit(unlist(IniStates))) != NState) {
+    #   stop(paste0("the length of IniStates must be ", NState, " for the chosen FUN_MOD"))
+    # }
+    if ((!"CemaNeige" %in% ObjectClass &  inherits(IniStates, "CemaNeige")) |
+        ( "CemaNeige" %in% ObjectClass & !inherits(IniStates, "CemaNeige"))) {
+      stop("'FUN_MOD' and 'IniStates' must be both of class 'CemaNeige'")
+    }
+    if (( "hysteresis" %in% ObjectClass & !inherits(IniStates, "hysteresis")) |
+        (!"hysteresis" %in% ObjectClass &  inherits(IniStates, "hysteresis"))) {
+      stop("'IsHyst' and 'IniStates' are not consistent on the use of the hysteresis")
+    }
+    if (!"CemaNeige" %in% ObjectClass & any(is.na(IniStates$CemaNeigeLayers$G  ))) {
+      IniStates$CemaNeigeLayers$G <- NULL
+    }
+    if (!"CemaNeige" %in% ObjectClass & any(is.na(IniStates$CemaNeigeLayers$eTG))) {
+      IniStates$CemaNeigeLayers$eTG <- NULL
+    }
+    if (!"CemaNeige" %in% ObjectClass & any(is.na(IniStates$CemaNeigeLayers$Gthr))) {
+      IniStates$CemaNeigeLayers$Gthr <- NULL
+    }
+    if (!"CemaNeige" %in% ObjectClass & any(is.na(IniStates$CemaNeigeLayers$Glocmax))) {
+      IniStates$CemaNeigeLayers$Glocmax <- NULL
+    }
+    IniStates$Store$Rest <- rep(NA, 3)
+    IniStates <- unlist(IniStates)
+    IniStates[is.na(IniStates)] <- 0
+    if ("monthly" %in% ObjectClass) {
+      IniStates <- IniStates[seq_len(NState)]
+    }
+  } else {
+    IniStates <- as.double(rep(0.0, NState))
+  }
+
+
+  ##check_Outputs_Cal_and_Sim
+
+  ##Outputs_all
+  Outputs_all <- c("DatesR", unlist(FortranOutputs), "WarmUpQsim", "StateEnd", "Param")
+  if (FeatFUN_MOD$IsSD) {
+    Outputs_all <- c(Outputs_all, "QsimDown", "Qsim_m3")
+  }
+
+  ##check_Outputs_Sim
+  if (!is.vector(Outputs_Sim)) {
+    stop("'Outputs_Sim' must be a vector of characters")
+  }
+  if (!is.character(Outputs_Sim)) {
+    stop("'Outputs_Sim' must be a vector of characters")
+  }
+  if (sum(is.na(Outputs_Sim)) != 0) {
+    stop("'Outputs_Sim' must not contain NA")
+  }
+  if ("all" %in% Outputs_Sim) {
+    Outputs_Sim <- Outputs_all
+  }
+  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"))
+  }
+  Outputs_Sim <- Outputs_Sim[!duplicated(Outputs_Sim)]
+
+
+  ##check_Outputs_Cal
+  if (is.null(Outputs_Cal)) {
+    if ("GR" %in% ObjectClass) {
+      Outputs_Cal <- c("Qsim", "Param")
+      if ("CemaNeige" %in% ObjectClass) {
+        Outputs_Cal <- c("PliqAndMelt", Outputs_Cal)
+      }
+    } else if ("CemaNeige" %in% ObjectClass) {
+      Outputs_Cal <- c("all")
+    }
+  } else {
+    if (!is.vector(Outputs_Cal)) {
+      stop("'Outputs_Cal' must be a vector of characters")
+    }
+    if (!is.character(Outputs_Cal)) {
+      stop("'Outputs_Cal' must be a vector of characters")
+    }
+    if (sum(is.na(Outputs_Cal)) != 0) {
+      stop("'Outputs_Cal' must not contain NA")
+    }
+  }
+  if ("all" %in% Outputs_Cal) {
+    Outputs_Cal <- Outputs_all
+  }
+  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"))
+  }
+  Outputs_Cal <- unique(Outputs_Cal)
+
+
+  ##check_MeanAnSolidPrecip
+  if ("CemaNeige" %in% ObjectClass & is.null(MeanAnSolidPrecip)) {
+    NLayers <- length(InputsModel$LayerPrecip)
+    SolidPrecip <- NULL
+    for (iLayer in 1:NLayers) {
+      if (iLayer == 1) {
+        SolidPrecip <-
+          InputsModel$LayerFracSolidPrecip[[1]] * InputsModel$LayerPrecip[[iLayer]] /
+          NLayers
+      } else {
+        SolidPrecip <- SolidPrecip + InputsModel$LayerFracSolidPrecip[[iLayer]] * InputsModel$LayerPrecip[[iLayer]] / NLayers
+      }
+    }
+    Factor <- NULL
+    if (inherits(InputsModel, "hourly")) {
+      Factor <- 365.25 * 24
+    }
+    if (inherits(InputsModel, "daily")) {
+      Factor <- 365.25
+    }
+    if (inherits(InputsModel, "monthly")) {
+      Factor <- 12
+    }
+    if (inherits(InputsModel, "yearly")) {
+      Factor <- 1
+    }
+    if (is.null(Factor)) {
+      stop("'InputsModel' must be of class 'hourly', 'daily', 'monthly' or 'yearly'")
+    }
+    MeanAnSolidPrecip <- rep(mean(SolidPrecip) * Factor, NLayers)
+    ### default value: same Gseuil for all layers
+    if (warnings) {
+      warning("'MeanAnSolidPrecip' not defined: it was automatically set to c(",
+              paste(round(MeanAnSolidPrecip), collapse = ","), ")  from the 'InputsModel' given to the function. ",
+              "Be careful in case your application is (short-term) forecasting.\n")
+    }
+  }
+  if ("CemaNeige" %in% ObjectClass & !is.null(MeanAnSolidPrecip)) {
+    if (!is.vector(MeanAnSolidPrecip)) {
+      stop(paste0("'MeanAnSolidPrecip' must be a vector of numeric values"))
+    }
+    if (!is.numeric(MeanAnSolidPrecip)) {
+      stop(paste0("'MeanAnSolidPrecip' must be a vector of numeric values"))
+    }
+    if (length(MeanAnSolidPrecip) != NLayers) {
+      stop(paste0("'MeanAnSolidPrecip' must be a numeric vector of length ", NLayers, ""))
+    }
+  }
+
+
+  ##check_PliqAndMelt
+  if ("GR" %in% ObjectClass & "CemaNeige" %in% ObjectClass) {
+    if (!"PliqAndMelt" %in% Outputs_Cal & !"all" %in% Outputs_Cal) {
+      WTxt <- NULL
+      WTxt <- paste0(WTxt, "'PliqAndMelt' was not defined in 'Outputs_Cal' but is needed to feed the hydrological model with the snow modele outputs \n")
+      WTxt <- paste0(WTxt, ": it was automatically added \n")
+      if (!is.null(WTxt) & warnings) {
+        warning(WTxt)
+      }
+      Outputs_Cal <- c(Outputs_Cal, "PliqAndMelt")
+    }
+    if (!"PliqAndMelt" %in% Outputs_Sim & !"all" %in% Outputs_Sim) {
+      WTxt <- NULL
+      WTxt <- paste0(WTxt, "'PliqAndMelt' was not defined in 'Outputs_Sim' but is needed to feed the hydrological model with the snow modele outputs \n")
+      WTxt <- paste0(WTxt, ": it was automatically added \n")
+      if (!is.null(WTxt) & warnings) {
+        warning(WTxt)
+      }
+      Outputs_Sim <- c(Outputs_Sim, "PliqAndMelt")
+    }
+  }
+
+
+  ##check_Qsim
+  if ("GR" %in% ObjectClass) {
+    if (!"Qsim" %in% Outputs_Cal & !"all" %in% Outputs_Cal) {
+      WTxt <- NULL
+      WTxt <- paste0(WTxt, "'Qsim' was not defined in 'Outputs_Cal' \n")
+      WTxt <- paste0(WTxt, ": it was automatically added \n")
+      if (!is.null(WTxt) & warnings) {
+        warning(WTxt)
+      }
+      Outputs_Cal <- c(Outputs_Cal, "Qsim")
+    }
+    if (!"Qsim" %in% Outputs_Sim & !"all" %in% Outputs_Sim) {
+      WTxt <- NULL
+      WTxt <- paste0(WTxt, "'Qsim' was not defined in 'Outputs_Sim' \n")
+      WTxt <- paste0(WTxt, ": it was automatically added \n")
+      if (!is.null(WTxt) & warnings) {
+        warning(WTxt)
+      }
+      Outputs_Sim <- c(Outputs_Sim, "Qsim")
+    }
+  }
+
+  ##Create_RunOptions
+  RunOptions <- list(IndPeriod_WarmUp = IndPeriod_WarmUp,
+                     IndPeriod_Run = IndPeriod_Run,
+                     IniStates = IniStates,
+                     IniResLevels = IniResLevels,
+                     Outputs_Cal = Outputs_Cal,
+                     Outputs_Sim = Outputs_Sim,
+                     FortranOutputs = FortranOutputs,
+                     FeatFUN_MOD = FeatFUN_MOD)
+
+  if ("CemaNeige" %in% ObjectClass) {
+    RunOptions <- c(RunOptions, list(MeanAnSolidPrecip = MeanAnSolidPrecip))
+  }
+  if ("interception" %in% ObjectClass) {
+    RunOptions <- c(RunOptions, list(Imax = Imax))
+  }
+  class(RunOptions) <- c("RunOptions", ObjectClass)
+
+  return(RunOptions)
+
+
+}
+
diff --git a/R/RunModel.R b/R/RunModel.R
index d68bf066453a2fdad3d907ff366b6a63e675a737..378f87e3f7963dc3f9809c5b161f9dfb4f882ca1 100644
--- a/R/RunModel.R
+++ b/R/RunModel.R
@@ -1,21 +1,21 @@
-RunModel <- function(InputsModel, RunOptions, Param, FUN_MOD, ...) {
-
-  FUN_MOD <- match.fun(FUN_MOD)
-
-  if (inherits(InputsModel, "SD") && !identical(FUN_MOD, RunModel_Lag)) {
-    # Lag model take one parameter at the beginning of the vector
-    iFirstParamRunOffModel <- 2
-    RunOptions$FeatFUN_MOD$NbParam <- RunOptions$FeatFUN_MOD$NbParam - 1
-  } else {
-    # All parameters
-    iFirstParamRunOffModel <- 1
-  }
-
-  OutputsModel <- FUN_MOD(InputsModel = InputsModel, RunOptions = RunOptions,
-                          Param = Param[iFirstParamRunOffModel:length(Param)], ...)
-
-  if (inherits(InputsModel, "SD") && !identical(FUN_MOD, RunModel_Lag)) {
-    OutputsModel <- RunModel_Lag(InputsModel, RunOptions, Param[1], OutputsModel)
-  }
-  return(OutputsModel)
-}
+RunModel <- function(InputsModel, RunOptions, Param, FUN_MOD, FUN_SD = RunModel_Lag,...) {
+
+  FUN_MOD <- match.fun(FUN_MOD)
+
+  if (inherits(InputsModel, "SD") && !identical(FUN_MOD, FUN_SD)) {
+    # Lag model take one parameter at the beginning of the vector
+    iFirstParamRunOffModel <- RunOptions$FeatFUN_MOD$ParamSD + 1
+    RunOptions$FeatFUN_MOD$NbParam <- RunOptions$FeatFUN_MOD$NbParam - RunOptions$FeatFUN_MOD$ParamSD
+  } else {
+    # All parameters
+    iFirstParamRunOffModel <- 1
+  }
+  OutputsModel <- FUN_MOD(InputsModel = InputsModel, RunOptions = RunOptions,
+                          Param = Param[iFirstParamRunOffModel:length(Param)], ...)
+
+  if (inherits(InputsModel, "SD") && !identical(FUN_MOD, FUN_SD)) {
+    OutputsModel <- RunModel_Lag(InputsModel, RunOptions, Param[1:RunOptions$FeatFUN_MOD$ParamSD],
+                                 OutputsModel)
+  }
+  return(OutputsModel)
+}
diff --git a/R/TransfoParam_GR5H.R b/R/TransfoParam_GR5H.R
index cd9b604d6fec293a6298bb517deaaeb7fa17bf54..72c95e4da3947e02ace027541cc844aece56b31a 100644
--- a/R/TransfoParam_GR5H.R
+++ b/R/TransfoParam_GR5H.R
@@ -1,50 +1,49 @@
-TransfoParam_GR5H <- function(ParamIn, Direction) {
-
-  ## number of model parameters
-  NParam <- 5L
-
-
-  ## check arguments
-  isVecParamIn <- is.vector(ParamIn)
-  if (isVecParamIn) {
-    ParamIn <- matrix(ParamIn, nrow = 1)
-  }
-  if (!inherits(ParamIn, "matrix")) {
-    stop("'ParamIn' must be of class 'matrix'")
-  }
-  if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
-    stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
-  }
-  if (ncol(ParamIn) != NParam) {
-    stop(sprintf("the GR5H model requires %i parameters", NParam))
-  }
-
-
-  ## transformation
-  if (Direction == "TR") {
-    ParamOut <- ParamIn
-    ParamOut[, 1] <- exp(ParamIn[, 1])                             ### GR5H X1 (production store capacity)
-    ParamOut[, 2] <- sinh(ParamIn[, 2])                            ### GR5H X2 (groundwater exchange coefficient)
-    ParamOut[, 3] <- exp(ParamIn[, 3])                             ### GR5H X3 (routing store capacity)
-    ParamOut[, 4] <- 480 + (480 - 0.01) * (ParamIn[, 4] - 10) / 20 ### GR5H X4 (unit hydrograph time constant)
-    ParamOut[, 5] <- (ParamIn[, 5] + 10) / 20                      ### GR5H X5 (groundwater exchange coefficient 2)
-  }
-  if (Direction == "RT") {
-    ParamOut <- ParamIn
-    ParamOut[, 1] <- log(ParamIn[, 1])                             ### GR5H X1 (production store capacity)
-    ParamOut[, 2] <- asinh(ParamIn[, 2])                           ### GR5H X2 (groundwater exchange coefficient)
-    ParamOut[, 3] <- log(ParamIn[, 3])                             ### GR5H X3 (routing store capacity)
-    ParamOut[, 4] <- (ParamIn[, 4] - 480) * 20 / (480 - 0.01) + 10 ### GR5H X4 (unit hydrograph time constant)
-    ParamOut[, 5] <- ParamIn[, 5] * 20 - 10                        ### GR5H X5 (groundwater exchange coefficient 2)
-  }
-
-
-  ## check output
-  if (isVecParamIn) {
-    ParamOut <- as.vector(ParamOut)
-  }
-
-  return(ParamOut)
-
-}
-
+TransfoParam_GR5H <- function(ParamIn, Direction) {
+
+  ## number of model parameters
+  NParam <- 5L
+
+  ## check arguments
+  isVecParamIn <- is.vector(ParamIn)
+  if (isVecParamIn) {
+    ParamIn <- matrix(ParamIn, nrow = 1)
+  }
+  if (!inherits(ParamIn, "matrix")) {
+    stop("'ParamIn' must be of class 'matrix'")
+  }
+  if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
+    stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
+  }
+  if (ncol(ParamIn) != NParam) {
+    stop(sprintf("the GR5H model requires %i parameters", NParam))
+  }
+
+
+  ## transformation
+  if (Direction == "TR") {
+    ParamOut <- ParamIn
+    ParamOut[, 1] <- exp(ParamIn[, 1])                             ### GR5H X1 (production store capacity)
+    ParamOut[, 2] <- sinh(ParamIn[, 2])                            ### GR5H X2 (groundwater exchange coefficient)
+    ParamOut[, 3] <- exp(ParamIn[, 3])                             ### GR5H X3 (routing store capacity)
+    ParamOut[, 4] <- 480 + (480 - 0.01) * (ParamIn[, 4] - 10) / 20 ### GR5H X4 (unit hydrograph time constant)
+    ParamOut[, 5] <- (ParamIn[, 5] + 10) / 20                      ### GR5H X5 (groundwater exchange coefficient 2)
+  }
+  if (Direction == "RT") {
+    ParamOut <- ParamIn
+    ParamOut[, 1] <- log(ParamIn[, 1])                             ### GR5H X1 (production store capacity)
+    ParamOut[, 2] <- asinh(ParamIn[, 2])                           ### GR5H X2 (groundwater exchange coefficient)
+    ParamOut[, 3] <- log(ParamIn[, 3])                             ### GR5H X3 (routing store capacity)
+    ParamOut[, 4] <- (ParamIn[, 4] - 480) * 20 / (480 - 0.01) + 10 ### GR5H X4 (unit hydrograph time constant)
+    ParamOut[, 5] <- ParamIn[, 5] * 20 - 10                        ### GR5H X5 (groundwater exchange coefficient 2)
+  }
+
+
+  ## check output
+  if (isVecParamIn) {
+    ParamOut <- as.vector(ParamOut)
+  }
+
+  return(ParamOut)
+
+}
+
diff --git a/R/Utils.R b/R/Utils.R
index 505f47f3208666db0c54d652dc54ea37ff2dae31..2e84abd2ee019e06ec830ce836f31cc7cdf7ed2b 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -1,263 +1,268 @@
-
-## =================================================================================
-## function to check
-## =================================================================================
-
-# .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.")
-#     }
-#   }
-# }
-
-
-
-## =================================================================================
-## function to extract model features
-## =================================================================================
-
-## table of feature models
-.FeatModels <- function() {
-  path <- system.file("modelsFeatures/FeatModelsGR.csv", package = "airGR")
-  read.table(path, header = TRUE, sep = ";", stringsAsFactors = FALSE)
-}
-
-
-## function to extract model features
-.GetFeatModel <- function(FUN_MOD, DatesR = NULL) {
-  FeatMod <- .FeatModels()
-  NameFunMod <- ifelse(test = FeatMod$Pkg %in% "airGR",
-                       yes  = paste("RunModel", FeatMod$NameMod, sep = "_"),
-                       no   = FeatMod$NameMod)
-  FunMod <- lapply(NameFunMod, FUN = match.fun)
-  IdMod <- which(sapply(FunMod, FUN = function(x) identical(FUN_MOD, x)))
-  if (length(IdMod) < 1) {
-    stop("'FUN_MOD' must be one of ", paste(NameFunMod, collapse = ", "))
-  } else {
-    res <- as.list(FeatMod[IdMod, ])
-    res$NameFunMod <- NameFunMod[IdMod]
-    if (!is.null(DatesR)) {
-      DiffTimeStep <- as.numeric(difftime(DatesR[length(DatesR)],
-                                          DatesR[length(DatesR)-1],
-                                          units = "secs"))
-      if (is.na(res$TimeUnit)) {
-        if (any(DiffTimeStep %in% 3600:3601)) { # 3601: leap second
-          res$TimeUnit <- "hourly"
-        } else {
-          res$TimeUnit <- "daily"
-        }
-      }
-    }
-    res$TimeStep <- switch(res$TimeUnit,
-                           hourly  =       1,
-                           daily   =       1 * 24,
-                           monthly =   28:31 * 24,
-                           yearly  = 365:366 * 24)
-    res$TimeStepMean <- switch(res$TimeUnit,
-                               hourly  =           1,
-                               daily   =           1 * 24,
-                               monthly = 365.25 / 12 * 24,
-                               yearly  =      365.25 * 24)
-    res$TimeStep     <- res$TimeStep * 3600
-    res$TimeStepMean <- as.integer(res$TimeStepMean * 3600)
-    res$Class <- c(res$TimeUnit, res$Class)
-    res$CodeModHydro <- res$CodeMod
-    if (grepl("CemaNeige", res$NameFunMod)) {
-      res$Class <- c(res$Class, "CemaNeige")
-      res$CodeModHydro <- gsub("CemaNeige", "", res$CodeMod)
-    }
-    res$Class <- res$Class[!is.na(res$Class)]
-    if (!is.null(DatesR)) {
-      if (all(DiffTimeStep != res$TimeStep)) {
-        stop("the time step of the model inputs must be ", res$TimeUnit)
-      }
-    }
-    return(res)
-  }
-}
-
-
-
-## =================================================================================
-## function to manage Fortran outputs
-## =================================================================================
-
-.FortranOutputs <- function(GR = NULL, isCN = FALSE) {
-
-  outGR <- NULL
-  outCN <- NULL
-
-  if (is.null(GR)) {
-    GR <- ""
-  }
-  if (GR == "GR1A") {
-    outGR <- c("PotEvap", "Precip",
-               "Qsim")
-  } else if (GR == "GR2M") {
-    outGR <- c("PotEvap", "Precip", "Prod", "Pn", "Ps",
-               "AE",
-               "Perc", "PR",
-               "Rout",
-               "AExch",
-               "Qsim")
-  } else if (GR == "GR5H") {
-    outGR <- c("PotEvap", "Precip", "Interc", "Prod", "Pn", "Ps",
-               "AE", "EI", "ES",
-               "Perc", "PR",
-               "Q9", "Q1",
-               "Rout", "Exch",
-               "AExch1", "AExch2",
-               "AExch", "QR",
-               "QD",
-               "Qsim")
-  } else if (GR %in% c("GR4J", "GR5J", "GR4H")) {
-    outGR <- c("PotEvap", "Precip", "Prod", "Pn", "Ps",
-               "AE",
-               "Perc", "PR",
-               "Q9", "Q1",
-               "Rout", "Exch",
-               "AExch1", "AExch2",
-               "AExch", "QR",
-               "QD",
-               "Qsim")
-  } else if (GR == "GR6J") {
-    outGR <- c("PotEvap", "Precip", "Prod", "Pn", "Ps",
-               "AE",
-               "Perc", "PR",
-               "Q9", "Q1",
-               "Rout", "Exch",
-               "AExch1", "AExch2",
-               "AExch", "QR",
-               "QRExp", "Exp",
-               "QD",
-               "Qsim")
-  }
-  if (isCN) {
-    outCN <- c("Pliq", "Psol",
-               "SnowPack", "ThermalState", "Gratio",
-               "PotMelt", "Melt", "PliqAndMelt", "Temp",
-               "Gthreshold", "Glocalmax")
-  }
-
-  res <- list(GR = outGR, CN = outCN)
-
-}
-
-
-
-## =================================================================================
-## functions to extract parts of InputsModel or OutputsModel objects
-## =================================================================================
-
-## InputsModel
-
-.ExtractInputsModel <- function(x, i) {
-  res <- lapply(x, function(x) {
-    if (is.matrix(x)) {
-      res0 <- x[i, , drop = FALSE]
-    }
-    if (is.vector(x) | inherits(x, "POSIXt")) {
-      res0 <- x[i]
-    }
-    if (is.list(x) & !inherits(x, "POSIXt")) {
-      if (inherits(x, "OutputsModel")) {
-        res0 <- .ExtractOutputsModel(x = x, i = i)
-      } else {
-        res0 <- .ExtractInputsModel(x = x, i = i)
-      }
-    }
-    return(res0)
-  })
-  if (!is.null(x$ZLayers)) {
-    res$ZLayers <- x$ZLayers
-  }
-  if (inherits(x, "SD")) {
-    res$LengthHydro <- x$LengthHydro
-    res$BasinAreas  <- x$BasinAreas
-  }
-  class(res) <- class(x)
-  res
-}
-
-'[.InputsModel' <- function(x, i) {
-  if (!inherits(x, "InputsModel")) {
-    stop("'x' must be of class 'InputsModel'")
-  }
-  if (is.factor(i)) {
-    i <- as.character(i)
-  }
-  if (is.numeric(i)) {
-    .ExtractInputsModel(x, i)
-  } else {
-    NextMethod()
-  }
-}
-
-
-## OutputsModel
-
-.ExtractOutputsModel <- function(x, i) {
-  res <- lapply(x, function(x) {
-    if (is.matrix(x)  && length(dim(x)) == 2L) {
-      res0 <- x[i, ]
-    }
-    if (is.array(x) && length(dim(x)) == 3L) {
-      res0 <- x[i, , ]
-    }
-    if (is.vector(x) | inherits(x, "POSIXt")) {
-      res0 <- x[i]
-    }
-    if (is.list(x) & !inherits(x, "POSIXt")) {
-      res0 <- .ExtractOutputsModel(x = x, i = i)
-    }
-    return(res0)
-  })
-  if (!is.null(x$RunOptions)) {
-    res$RunOptions <- x$RunOptions
-  }
-  if (!is.null(x$StateEnd)) {
-    res$StateEnd <- x$StateEnd
-  }
-  class(res) <- class(x)
-  res
-}
-
-.IndexOutputsModel <- function(x, i) {
-  # '[.OutputsModel' <- function(x, i) {
-  if (!inherits(x, "OutputsModel")) {
-    stop("'x' must be of class 'OutputsModel'")
-  }
-  if (is.factor(i)) {
-    i <- as.character(i)
-  }
-  if (is.numeric(i)) {
-    .ExtractOutputsModel(x, i)
-  } else {
-    NextMethod()
-  }
-}
-
-
-
-## =================================================================================
-## function to try to set local time in English
-## =================================================================================
-
-.TrySetLcTimeEN <- function() {
-  locale <- list("English_United Kingdom",
-                 "en_US",
-                 "en_US.UTF-8",
-                 "en_US.utf8",
-                 "en")
-  dateTest <- as.POSIXct("2000-02-15", tz = "UTC", format = "%Y-%m-%d")
-  monthTestTarget <- "February"
-  monthTest <- function() {
-    format(dateTest, format = "%B")
-  }
-  lapply(locale, function(x) {
-    if (monthTest() != monthTestTarget) {
-      Sys.setlocale(category = "LC_TIME", locale = x)
-    }
-  })
-}
+
+## =================================================================================
+## function to check
+## =================================================================================
+
+# .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.")
+#     }
+#   }
+# }
+
+
+
+## =================================================================================
+## function to extract model features
+## =================================================================================
+
+## table of feature models
+.FeatModels <- function() {
+  path <- system.file("modelsFeatures/FeatModelsGR.csv", package = "airGR")
+  read.table(path, header = TRUE, sep = ";", stringsAsFactors = FALSE)
+}
+
+
+## function to extract model features
+.GetFeatModel <- function(FUN_MOD, DatesR = NULL, FUN_SD = NULL) {
+  FeatMod <- .FeatModels()
+  NameFunMod <- ifelse(test = FeatMod$Pkg %in% "airGR",
+                       yes  = paste("RunModel", FeatMod$NameMod, sep = "_"),
+                       no   = FeatMod$NameMod)
+  FunMod <- lapply(NameFunMod, FUN = match.fun)
+  IdMod <- which(sapply(FunMod, FUN = function(x) identical(FUN_MOD, x)))
+  if (length(IdMod) < 1) {
+    stop("'FUN_MOD' must be one of ", paste(NameFunMod, collapse = ", "))
+  } else {
+    res <- as.list(FeatMod[IdMod, ])
+    res$NameFunMod <- NameFunMod[IdMod]
+    if (!is.null(DatesR)) {
+      DiffTimeStep <- as.numeric(difftime(DatesR[length(DatesR)],
+                                          DatesR[length(DatesR)-1],
+                                          units = "secs"))
+      if (is.na(res$TimeUnit)) {
+        if (any(DiffTimeStep %in% 3600:3601)) { # 3601: leap second
+          res$TimeUnit <- "hourly"
+        } else {
+          res$TimeUnit <- "daily"
+        }
+      }
+    }
+    res$TimeStep <- switch(res$TimeUnit,
+                           hourly  =       1,
+                           daily   =       1 * 24,
+                           monthly =   28:31 * 24,
+                           yearly  = 365:366 * 24)
+    res$TimeStepMean <- switch(res$TimeUnit,
+                               hourly  =           1,
+                               daily   =           1 * 24,
+                               monthly = 365.25 / 12 * 24,
+                               yearly  =      365.25 * 24)
+    res$TimeStep     <- res$TimeStep * 3600
+    res$TimeStepMean <- as.integer(res$TimeStepMean * 3600)
+    res$Class <- c(res$TimeUnit, res$Class)
+    res$CodeModHydro <- res$CodeMod
+    if (grepl("CemaNeige", res$NameFunMod)) {
+      res$Class <- c(res$Class, "CemaNeige")
+      res$CodeModHydro <- gsub("CemaNeige", "", res$CodeMod)
+    }
+    res$Class <- res$Class[!is.na(res$Class)]
+    if (!is.null(DatesR)) {
+      if (all(DiffTimeStep != res$TimeStep)) {
+        stop("the time step of the model inputs must be ", res$TimeUnit)
+      }
+    }
+    if (!is.null(FUN_SD)){
+      if (identical(match.fun(RunModel_Lag), FUN_SD)){
+        res$ParamSD <- FeatMod$NbParam[14]
+      }
+    }
+    return(res)
+  }
+}
+
+
+
+## =================================================================================
+## function to manage Fortran outputs
+## =================================================================================
+
+.FortranOutputs <- function(GR = NULL, isCN = FALSE) {
+
+  outGR <- NULL
+  outCN <- NULL
+
+  if (is.null(GR)) {
+    GR <- ""
+  }
+  if (GR == "GR1A") {
+    outGR <- c("PotEvap", "Precip",
+               "Qsim")
+  } else if (GR == "GR2M") {
+    outGR <- c("PotEvap", "Precip", "Prod", "Pn", "Ps",
+               "AE",
+               "Perc", "PR",
+               "Rout",
+               "AExch",
+               "Qsim")
+  } else if (GR == "GR5H") {
+    outGR <- c("PotEvap", "Precip", "Interc", "Prod", "Pn", "Ps",
+               "AE", "EI", "ES",
+               "Perc", "PR",
+               "Q9", "Q1",
+               "Rout", "Exch",
+               "AExch1", "AExch2",
+               "AExch", "QR",
+               "QD",
+               "Qsim")
+  } else if (GR %in% c("GR4J", "GR5J", "GR4H")) {
+    outGR <- c("PotEvap", "Precip", "Prod", "Pn", "Ps",
+               "AE",
+               "Perc", "PR",
+               "Q9", "Q1",
+               "Rout", "Exch",
+               "AExch1", "AExch2",
+               "AExch", "QR",
+               "QD",
+               "Qsim")
+  } else if (GR == "GR6J") {
+    outGR <- c("PotEvap", "Precip", "Prod", "Pn", "Ps",
+               "AE",
+               "Perc", "PR",
+               "Q9", "Q1",
+               "Rout", "Exch",
+               "AExch1", "AExch2",
+               "AExch", "QR",
+               "QRExp", "Exp",
+               "QD",
+               "Qsim")
+  }
+  if (isCN) {
+    outCN <- c("Pliq", "Psol",
+               "SnowPack", "ThermalState", "Gratio",
+               "PotMelt", "Melt", "PliqAndMelt", "Temp",
+               "Gthreshold", "Glocalmax")
+  }
+
+  res <- list(GR = outGR, CN = outCN)
+
+}
+
+
+
+## =================================================================================
+## functions to extract parts of InputsModel or OutputsModel objects
+## =================================================================================
+
+## InputsModel
+
+.ExtractInputsModel <- function(x, i) {
+  res <- lapply(x, function(x) {
+    if (is.matrix(x)) {
+      res0 <- x[i, , drop = FALSE]
+    }
+    if (is.vector(x) | inherits(x, "POSIXt")) {
+      res0 <- x[i]
+    }
+    if (is.list(x) & !inherits(x, "POSIXt")) {
+      if (inherits(x, "OutputsModel")) {
+        res0 <- .ExtractOutputsModel(x = x, i = i)
+      } else {
+        res0 <- .ExtractInputsModel(x = x, i = i)
+      }
+    }
+    return(res0)
+  })
+  if (!is.null(x$ZLayers)) {
+    res$ZLayers <- x$ZLayers
+  }
+  if (inherits(x, "SD")) {
+    res$LengthHydro <- x$LengthHydro
+    res$BasinAreas  <- x$BasinAreas
+  }
+  class(res) <- class(x)
+  res
+}
+
+'[.InputsModel' <- function(x, i) {
+  if (!inherits(x, "InputsModel")) {
+    stop("'x' must be of class 'InputsModel'")
+  }
+  if (is.factor(i)) {
+    i <- as.character(i)
+  }
+  if (is.numeric(i)) {
+    .ExtractInputsModel(x, i)
+  } else {
+    NextMethod()
+  }
+}
+
+
+## OutputsModel
+
+.ExtractOutputsModel <- function(x, i) {
+  res <- lapply(x, function(x) {
+    if (is.matrix(x)  && length(dim(x)) == 2L) {
+      res0 <- x[i, ]
+    }
+    if (is.array(x) && length(dim(x)) == 3L) {
+      res0 <- x[i, , ]
+    }
+    if (is.vector(x) | inherits(x, "POSIXt")) {
+      res0 <- x[i]
+    }
+    if (is.list(x) & !inherits(x, "POSIXt")) {
+      res0 <- .ExtractOutputsModel(x = x, i = i)
+    }
+    return(res0)
+  })
+  if (!is.null(x$RunOptions)) {
+    res$RunOptions <- x$RunOptions
+  }
+  if (!is.null(x$StateEnd)) {
+    res$StateEnd <- x$StateEnd
+  }
+  class(res) <- class(x)
+  res
+}
+
+.IndexOutputsModel <- function(x, i) {
+  # '[.OutputsModel' <- function(x, i) {
+  if (!inherits(x, "OutputsModel")) {
+    stop("'x' must be of class 'OutputsModel'")
+  }
+  if (is.factor(i)) {
+    i <- as.character(i)
+  }
+  if (is.numeric(i)) {
+    .ExtractOutputsModel(x, i)
+  } else {
+    NextMethod()
+  }
+}
+
+
+
+## =================================================================================
+## function to try to set local time in English
+## =================================================================================
+
+.TrySetLcTimeEN <- function() {
+  locale <- list("English_United Kingdom",
+                 "en_US",
+                 "en_US.UTF-8",
+                 "en_US.utf8",
+                 "en")
+  dateTest <- as.POSIXct("2000-02-15", tz = "UTC", format = "%Y-%m-%d")
+  monthTestTarget <- "February"
+  monthTest <- function() {
+    format(dateTest, format = "%B")
+  }
+  lapply(locale, function(x) {
+    if (monthTest() != monthTestTarget) {
+      Sys.setlocale(category = "LC_TIME", locale = x)
+    }
+  })
+}
diff --git a/R/UtilsCalibOptions.R b/R/UtilsCalibOptions.R
index aa4aaf5f4eeb8ea4418fbd2f86b2d8b53e17efae..92c9fa2f323d03293d3f35d249e83f23793e88d1 100644
--- a/R/UtilsCalibOptions.R
+++ b/R/UtilsCalibOptions.R
@@ -1,132 +1,136 @@
-.FunTransfo <- function(FeatFUN_MOD) {
-
-  IsHyst <- FeatFUN_MOD$IsHyst
-  IsSD <- FeatFUN_MOD$IsSD
-
-  ## set FUN_GR
-  if (FeatFUN_MOD$NameFunMod == "Cemaneige") {
-    if (IsHyst) {
-      FUN_GR <- TransfoParam_CemaNeigeHyst
-    } else {
-      FUN_GR <- TransfoParam_CemaNeige
-    }
-  } else {
-    # fatal error if the TransfoParam function does not exist
-    FUN_GR <- match.fun(sprintf("TransfoParam_%s", FeatFUN_MOD$CodeModHydro))
-  }
-
-  ## set FUN_SNOW
-  if ("CemaNeige" %in% FeatFUN_MOD$Class) {
-    if (IsHyst) {
-      FUN_SNOW <- TransfoParam_CemaNeigeHyst
-    } else {
-      FUN_SNOW <- TransfoParam_CemaNeige
-    }
-  }
-
-  ## set FUN_LAG
-  if (IsSD) {
-    FUN_LAG <- TransfoParam_Lag
-  }
-
-    ## set FUN_TRANSFO
-  if (! "CemaNeige" %in% FeatFUN_MOD$Class) {
-    if (!IsSD) {
-      FUN_TRANSFO <- FUN_GR
-    } else {
-      FUN_TRANSFO <- function(ParamIn, Direction) {
-        Bool <- is.matrix(ParamIn)
-        if (!Bool) {
-          ParamIn <- rbind(ParamIn)
-        }
-        ParamOut <- NA * ParamIn
-        NParam   <- ncol(ParamIn)
-        ParamOut[, 2:NParam] <- FUN_GR(ParamIn[, 2:NParam], Direction)
-        ParamOut[, 1       ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
-        if (!Bool) {
-          ParamOut <- ParamOut[1, ]
-        }
-        return(ParamOut)
-      }
-    }
-  } else {
-    if (IsHyst & !IsSD) {
-      FUN_TRANSFO <- function(ParamIn, Direction) {
-        Bool <- is.matrix(ParamIn)
-        if (!Bool) {
-          ParamIn <- rbind(ParamIn)
-        }
-        ParamOut <- NA * ParamIn
-        NParam   <- ncol(ParamIn)
-        ParamOut[, 1:(NParam - 4)     ] <- FUN_GR(ParamIn[, 1:(NParam - 4)], Direction)
-        ParamOut[, (NParam - 3):NParam] <- FUN_SNOW(ParamIn[, (NParam - 3):NParam], Direction)
-        if (!Bool) {
-          ParamOut <- ParamOut[1, ]
-        }
-        return(ParamOut)
-      }
-    }
-    if (!IsHyst & !IsSD) {
-      FUN_TRANSFO <- function(ParamIn, Direction) {
-        Bool <- is.matrix(ParamIn)
-        if (!Bool) {
-          ParamIn <- rbind(ParamIn)
-        }
-        ParamOut <- NA * ParamIn
-        NParam   <- ncol(ParamIn)
-        if (NParam <= 3) {
-          ParamOut[, 1:(NParam - 2)] <- FUN_GR(cbind(ParamIn[, 1:(NParam - 2)]), Direction)
-        } else {
-          ParamOut[, 1:(NParam - 2)] <- FUN_GR(ParamIn[, 1:(NParam - 2)], Direction)
-        }
-        ParamOut[, (NParam - 1):NParam] <- FUN_SNOW(ParamIn[, (NParam - 1):NParam], Direction)
-        if (!Bool) {
-          ParamOut <- ParamOut[1, ]
-        }
-        return(ParamOut)
-      }
-    }
-    if (IsHyst & IsSD) {
-      FUN_TRANSFO <- function(ParamIn, Direction) {
-        Bool <- is.matrix(ParamIn)
-        if (!Bool) {
-          ParamIn <- rbind(ParamIn)
-        }
-        ParamOut <- NA * ParamIn
-        NParam   <- ncol(ParamIn)
-        ParamOut[, 2:(NParam - 4)     ] <- FUN_GR(ParamIn[, 2:(NParam - 4)], Direction)
-        ParamOut[, (NParam - 3):NParam] <- FUN_SNOW(ParamIn[, (NParam - 3):NParam], Direction)
-        ParamOut[, 1                  ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
-        if (!Bool) {
-          ParamOut <- ParamOut[1, ]
-        }
-        return(ParamOut)
-      }
-    }
-    if (!IsHyst & IsSD) {
-      FUN_TRANSFO <- function(ParamIn, Direction) {
-        Bool <- is.matrix(ParamIn)
-        if (!Bool) {
-          ParamIn <- rbind(ParamIn)
-        }
-        ParamOut <- NA * ParamIn
-        NParam   <- ncol(ParamIn)
-        if (NParam <= 3) {
-          ParamOut[, 2:(NParam - 2)] <- FUN_GR(cbind(ParamIn[, 2:(NParam - 2)]), Direction)
-        } else {
-          ParamOut[, 2:(NParam - 2)] <- FUN_GR(ParamIn[, 2:(NParam - 2)],  Direction)
-        }
-        ParamOut[, (NParam - 1):NParam] <- FUN_SNOW(ParamIn[, (NParam - 1):NParam], Direction)
-        ParamOut[, 1                  ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
-        if (!Bool) {
-          ParamOut <- ParamOut[1, ]
-        }
-        return(ParamOut)
-      }
-    }
-  }
-  if (is.null(FUN_TRANSFO)) {
-    stop("'FUN_TRANSFO' was not found")
-  }
-  return(FUN_TRANSFO)
-}
+.FunTransfo <- function(FeatFUN_MOD, FUN_SD) {
+
+  IsHyst <- FeatFUN_MOD$IsHyst
+  IsSD <- FeatFUN_MOD$IsSD
+
+  ## set FUN_GR
+  if (FeatFUN_MOD$NameFunMod == "Cemaneige") {
+    if (IsHyst) {
+      FUN_GR <- TransfoParam_CemaNeigeHyst
+    } else {
+      FUN_GR <- TransfoParam_CemaNeige
+    }
+  } else {
+    # fatal error if the TransfoParam function does not exist
+    FUN_GR <- match.fun(sprintf("TransfoParam_%s", FeatFUN_MOD$CodeModHydro))
+  }
+
+  ## set FUN_SNOW
+  if ("CemaNeige" %in% FeatFUN_MOD$Class) {
+    if (IsHyst) {
+      FUN_SNOW <- TransfoParam_CemaNeigeHyst
+    } else {
+      FUN_SNOW <- TransfoParam_CemaNeige
+    }
+  }
+
+  ## set FUN_ROUT
+
+  if (IsSD) {
+    if (identical(RunModel_Lag, FUN_SD)){
+      FUN_ROUT <- TransfoParam_Lag
+    }
+  }
+
+
+    ## set FUN_TRANSFO
+  if (! "CemaNeige" %in% FeatFUN_MOD$Class) {
+    if (!IsSD) {
+      FUN_TRANSFO <- FUN_GR
+    } else {
+      FUN_TRANSFO <- function(ParamIn, Direction) {
+        Bool <- is.matrix(ParamIn)
+        if (!Bool) {
+          ParamIn <- rbind(ParamIn)
+        }
+        ParamOut <- NA * ParamIn
+        NParam   <- ncol(ParamIn)
+        ParamOut[, (FeatFUN_MOD$ParamSD + 1):NParam] <- FUN_GR(ParamIn[, (FeatFUN_MOD$ParamSD + 1):NParam], Direction)
+        ParamOut[, FeatFUN_MOD$ParamSD       ] <- FUN_ROUT(as.matrix(ParamIn[, FeatFUN_MOD$ParamSD]), Direction)
+        if (!Bool) {
+          ParamOut <- ParamOut[1, ]
+        }
+        return(ParamOut)
+      }
+    }
+  } else {
+    if (IsHyst & !IsSD) {
+      FUN_TRANSFO <- function(ParamIn, Direction) {
+        Bool <- is.matrix(ParamIn)
+        if (!Bool) {
+          ParamIn <- rbind(ParamIn)
+        }
+        ParamOut <- NA * ParamIn
+        NParam   <- ncol(ParamIn)
+        ParamOut[, 1:(NParam - 4)     ] <- FUN_GR(ParamIn[, 1:(NParam - 4)], Direction)
+        ParamOut[, (NParam - 3):NParam] <- FUN_SNOW(ParamIn[, (NParam - 3):NParam], Direction)
+        if (!Bool) {
+          ParamOut <- ParamOut[1, ]
+        }
+        return(ParamOut)
+      }
+    }
+    if (!IsHyst & !IsSD) {
+      FUN_TRANSFO <- function(ParamIn, Direction) {
+        Bool <- is.matrix(ParamIn)
+        if (!Bool) {
+          ParamIn <- rbind(ParamIn)
+        }
+        ParamOut <- NA * ParamIn
+        NParam   <- ncol(ParamIn)
+        if (NParam <= 3) {
+          ParamOut[, 1:(NParam - 2)] <- FUN_GR(cbind(ParamIn[, 1:(NParam - 2)]), Direction)
+        } else {
+          ParamOut[, 1:(NParam - 2)] <- FUN_GR(ParamIn[, 1:(NParam - 2)], Direction)
+        }
+        ParamOut[, (NParam - 1):NParam] <- FUN_SNOW(ParamIn[, (NParam - 1):NParam], Direction)
+        if (!Bool) {
+          ParamOut <- ParamOut[1, ]
+        }
+        return(ParamOut)
+      }
+    }
+    if (IsHyst & IsSD) {
+      FUN_TRANSFO <- function(ParamIn, Direction) {
+        Bool <- is.matrix(ParamIn)
+        if (!Bool) {
+          ParamIn <- rbind(ParamIn)
+        }
+        ParamOut <- NA * ParamIn
+        NParam   <- ncol(ParamIn)
+        ParamOut[, 2:(NParam - 4)     ] <- FUN_GR(ParamIn[, 2:(NParam - 4)], Direction)
+        ParamOut[, (NParam - 3):NParam] <- FUN_SNOW(ParamIn[, (NParam - 3):NParam], Direction)
+        ParamOut[, 1                  ] <- FUN_ROUT(as.matrix(ParamIn[, 1]), Direction)
+        if (!Bool) {
+          ParamOut <- ParamOut[1, ]
+        }
+        return(ParamOut)
+      }
+    }
+    if (!IsHyst & IsSD) {
+      FUN_TRANSFO <- function(ParamIn, Direction) {
+        Bool <- is.matrix(ParamIn)
+        if (!Bool) {
+          ParamIn <- rbind(ParamIn)
+        }
+        ParamOut <- NA * ParamIn
+        NParam   <- ncol(ParamIn)
+        if (NParam <= 3) {
+          ParamOut[, 2:(NParam - 2)] <- FUN_GR(cbind(ParamIn[, 2:(NParam - 2)]), Direction)
+        } else {
+          ParamOut[, 2:(NParam - 2)] <- FUN_GR(ParamIn[, 2:(NParam - 2)],  Direction)
+        }
+        ParamOut[, (NParam - 1):NParam] <- FUN_SNOW(ParamIn[, (NParam - 1):NParam], Direction)
+        ParamOut[, 1                  ] <- FUN_ROUT(as.matrix(ParamIn[, 1]), Direction)
+        if (!Bool) {
+          ParamOut <- ParamOut[1, ]
+        }
+        return(ParamOut)
+      }
+    }
+  }
+  if (is.null(FUN_TRANSFO)) {
+    stop("'FUN_TRANSFO' was not found")
+  }
+  return(FUN_TRANSFO)
+}