.Rbuildignore

@@ -2,3 +2,14 @@
 ^\.Rproj\.user$
 ^\.Rprofile$
 ^packrat/
+^tests/tmp/
+^\.gitlab-ci.yml$
+^\.regressionignore$
+^\.vignettechunkignore$
+^\.gitlab-ci\.yml$
+^\.vscode$
+^Rplots\.pdf$
+^ci$
+^data-raw$
+^revdep$
+^\.devcontainer$

.devcontainer/devcontainer.json

+{
+	"image": "rocker/geospatial:devel",
+	"customizations": {
+		"vscode": {
+			"extensions": [
+				"eamodio.gitlens",
+				"REditorSupport.r"
+			]
+		}
+	},
+	// Use 'postCreateCommand' to run commands after the container is created.
+	"postCreateCommand": "R -q -e 'install.packages(\"languageserver\");remotes::install_deps(dep = TRUE)'",
+	"postStartCommand": "R -q -e 'devtools::install()'"
+}

.gitignore

-.Rproj.user
+# Specific files for airGR
+packrat/lib*/
+
+# Compiled files
+/src/*.o
+/src/*.so
+/src/*.dll
+/src-*
+
+# Test temporary files
+/tests/tmp/
+*.pdf
+!man/figures/*.pdf
+
+# revdep
+/revdep/
+
+######################################################################################################
+### Generic .gitignore for R (source: https://github.com/github/gitignore/blob/master/R.gitignore) ###
+######################################################################################################
+
+# History files
 .Rhistory
+.Rapp.history
+
+# Session Data files
 .RData
-airGR.Rproj
-packrat/lib*/
+
+# User-specific files
+.Ruserdata
+
+# Example code in package build process
+*-Ex.R
+
+# Output files from R CMD build
+/*.tar.gz
+
+# Output files from R CMD check
+/*.Rcheck/
+
+# RStudio files
+.Rproj.user
+
+# produced vignettes
+vignettes/*.html
+vignettes/*.pdf
+
+# OAuth2 token, see https://github.com/hadley/httr/releases/tag/v0.3
+.httr-oauth
+
+# knitr and R markdown default cache directories
+*_cache/
+/cache/
+
+# Temporary files created by R markdown
+*.utf8.md
+*.knit.md
+
+# R Environment Variables
+.Renviron
+.Rprofile
+
+# pkgdown site
+docs/
+
+# vscode IDE
+.vscode/*
+*.code-workspace
+.history/

.gitlab-ci.yml

+stages:
+  - check
+  - scheduled_tests
+  - revdepcheck
+
+variables:
+  R_LIBS_USER: "$CI_PROJECT_DIR/ci/lib"
+
+default:
+  tags: [docker]
+
+cache:
+  key: "airGR-${R_VERSION}"
+  paths:
+    - $R_LIBS_USER
+
+before_script:
+  - mkdir -p $R_LIBS_USER
+  - echo "R_LIBS='$R_LIBS_USER'" > .Renviron
+  - echo 'options(repos = c(CRAN = "https://packagemanager.posit.co/cran/__linux__/focal/latest"))' > .Rprofile
+  - apt-get update && apt-get install -y libstdc++6
+  - R -q -e 'remotes::install_deps(dependencies = TRUE)'
+
+.R-devel:
+  image: rocker/geospatial:devel
+  variables:
+    R_VERSION: devel
+
+.R-patched:
+  only:
+    refs:
+      - tags
+      - schedules
+  image: rocker/geospatial:latest
+  variables:
+    R_VERSION: patched
+
+.R-oldrel:
+  only:
+    refs:
+      - tags
+      - schedules
+  image: rocker/geospatial:4.1
+  variables:
+    R_VERSION: oldrel
+
+.scheduled_tests:
+  only:
+    refs:
+      - tags
+      - schedules
+      - master
+  stage: scheduled_tests
+  script:
+    - R -q -e 'devtools::install()'
+    - Rscript tests/scheduled_tests/scheduled.R
+    - Rscript tests/scheduled_tests/regression.R stable
+    - R CMD INSTALL .
+    - Rscript tests/scheduled_tests/regression.R dev
+    - Rscript tests/scheduled_tests/regression.R compare
+
+.check:
+  stage: check
+  script:
+  - R -q -e 'tinytex::tlmgr("option repository https://ftp.tu-chemnitz.de/pub/tug/historic/systems/texlive/2021/tlnet-final")'
+  - tlmgr update --self && tlmgr install ec epstopdf-pkg amsmath
+  - R -q -e 'remotes::update_packages("rcmdcheck")'
+  - R -q -e 'rcmdcheck::rcmdcheck(args = c("--as-cran"), error_on = "warning")'
+
+.test_all:
+  stage: check
+  script:
+  - R -q -e 'testthat::test_local()'
+
+benchmark_devel:
+  extends: .R-devel
+  stage: check
+  allow_failure: true
+  script:
+  - R -q -e 'remotes::update_packages("microbenchmark", repos = "http://cran.r-project.org")'
+  - R -q -e 'install.packages("airGR", repos = "http://cran.r-project.org")'
+  - Rscript tests/scheduled_tests/benchmarkRunModel.R
+  - R CMD INSTALL .
+  - Rscript tests/scheduled_tests/benchmarkRunModel.R
+  artifacts:
+    paths:
+      - tests/tmp/benchmark.tsv
+      - tests/tmp/mean_execution_time.tsv
+
+scheduled_tests_patched:
+  extends:
+    - .scheduled_tests
+    - .R-patched
+
+scheduled_tests_devel:
+  extends:
+    - .scheduled_tests
+    - .R-devel
+
+scheduled_tests_oldrel:
+  extends:
+    - .scheduled_tests
+    - .R-oldrel
+
+check_patched:
+  extends:
+    - .check
+    - .R-patched
+
+check_devel:
+  extends:
+    - .check
+    - .R-devel
+
+test_all_patched:
+  extends:
+    - .test_all
+    - .R-patched
+
+test_all_devel:
+  extends:
+    - .test_all
+    - .R-devel
+
+test_all_oldrel:
+  extends:
+    - .test_all
+    - .R-oldrel
+
+revdepcheck_devel:
+  stage: revdepcheck
+  only:
+    refs:
+      - tags
+      - schedule
+      - master
+  extends: .R-devel
+  script:
+    - R -q -e 'remotes::install_github("https://github.com/r-lib/revdepcheck")'
+    - R -q -e 'revdepcheck::revdep_check(timeout = as.difftime(20, units = "mins"))'
+    - R -q -e 'stopifnot(all("+" == sapply(revdepcheck::revdep_summary(), "[[", "status")))'
+    - R -q -e 'if (any(sapply(revdepcheck::revdep_summary(), function(x) {any(x$cmp$change == 1)}))) stop()'
+  artifacts:
+    when: on_failure
+    paths:
+      - revdep/README.md
+      - revdep/problems.md
+      - revdep/failures.md
+      - revdep/cran.md
+      - revdep/checks/*.log

.regressionignore

+# .test-regression.ignore contains the list of topic/variables produces by
+# documentation examples that should be ignore in the regression test
+# The format of this file is: 5 lines of comments followed by one line by
+# ignored variable : [Topic]<SPACE>[Variable] or *<SPACE>[Variable] for every variable whatever the topic
+# Example for ignoring OutputsModel variable produced by example("RunModel_GR2M"): RunModel_GR2M OutputsModel

.vignettechunkignore

+# This file is used by the script tests/testthat/test-vignettes which test all
+# chunks including those with `eval=FALSE`
+# It serves to ignore chunks that should not be tested anyway
+# Format: `vignette file name`[space]`id of the chunk`
+V02.1_param_optim.Rmd hydroPSO1
+V02.1_param_optim.Rmd hydroPSO2
+V02.1_param_optim.Rmd resGLOB

DESCRIPTION

 Package: airGR
 Type: Package
 Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.3.0.17
-Date: 2019-05-21
+Version: 1.7.6.9000
+Date: 2023-10-25
 Authors@R: c(
   person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
-  person("Olivier", "Delaigue", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@irstea.fr"),
-  person("Guillaume", "Thirel", role = c("aut"), comment = c(ORCID = "0000-0002-1444-1830")),
+  person("Olivier", "Delaigue", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@inrae.fr"),
+  person("Guillaume", "Thirel", role = c("aut", "ths"), comment = c(ORCID = "0000-0002-1444-1830")),
+  person("David", "Dorchies", role = c("aut"), comment = c(ORCID = "0000-0002-6595-7984")),
   person("Charles", "Perrin", role = c("aut", "ths"), comment = c(ORCID = "0000-0001-8552-1881")),
   person("Claude", "Michel", role = c("aut", "ths")),
   person("Vazken", "Andréassian", role = c("ctb", "ths"), comment = c(ORCID = "0000-0001-7124-9303")),
-  person("François", "Bourgin", role = c("ctb"), comment = c(ORCID = "0000-0002-2820-7260", vignette = "'Parameter estimation' vignettes")),
+  person("François", "Bourgin", role = c("ctb"), comment = c(ORCID = "0000-0002-2820-7260")),
   person("Pierre", "Brigode", role = c("ctb"), comment = c(ORCID = "0000-0001-8257-0741")),
   person("Nicolas", "Le Moine", role = c("ctb")),
-  person("Thibaut", "Mathevet", role = c("ctb")),
+  person("Thibaut", "Mathevet", role = c("ctb"), comment = c(ORCID = "0000-0002-4142-4454")),
   person("Safouane", "Mouelhi", role = c("ctb")),
-  person("Ludovic", "Oudin", role = c("ctb")),
+  person("Ludovic", "Oudin", role = c("ctb"), comment = c(ORCID = "0000-0002-3712-0933")),
   person("Raji", "Pushpalatha", role = c("ctb")),
   person("Audrey", "Valéry", role = c("ctb"))
   )
-Depends: R (>= 3.0.1)
-Suggests: knitr, rmarkdown, coda, DEoptim, dplyr, FME, ggmcmc, hydroPSO, Rmalschains
-Description: Hydrological modelling tools developed at Irstea-Antony (HYCAR Research Unit, France). The package includes several conceptual rainfall-runoff models (GR4H, GR4J, GR5J, GR6J, GR2M, GR1A), a snow accumulation and melt model (CemaNeige) and the associated functions for their calibration and evaluation. Use help(airGR) for package description and references.
+Depends: R (>= 3.1.0)
+Imports:
+  graphics,
+  grDevices,
+  stats,
+  utils
+Suggests:
+  knitr, markdown, rmarkdown,
+  caRamel, coda, DEoptim, FME, ggmcmc, Rmalschains,
+  GGally, ggplot2,
+  testthat
+Description: Hydrological modelling tools developed at INRAE-Antony (HYCAR Research Unit, France). The package includes several conceptual rainfall-runoff models (GR4H, GR5H, GR4J, GR5J, GR6J, GR2M, GR1A) that can be applied either on a lumped or semi-distributed way. A snow accumulation and melt model (CemaNeige) and the associated functions for the calibration and evaluation of models are also included. Use help(airGR) for package description and references.
 License: GPL-2
 URL: https://hydrogr.github.io/airGR/
+BugReports: https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues
 NeedsCompilation: yes
 Encoding: UTF-8
 VignetteBuilder: knitr
+RoxygenNote: 7.1.1

NAMESPACE

@@ -8,7 +8,13 @@ useDynLib(airGR, .registration = TRUE)
 #####################################
 ##            S3 methods           ##
 #####################################
-S3method("plot", "OutputsModel")
+S3method('[', InputsModel)
+#S3method('[', OutputsModel) ### to add in version 2.0
+S3method(plot, OutputsModel)
+S3method(SeriesAggreg, data.frame)
+S3method(SeriesAggreg, list)
+S3method(SeriesAggreg, InputsModel)
+S3method(SeriesAggreg, OutputsModel)
 
 
 
@@ -18,8 +24,10 @@ S3method("plot", "OutputsModel")
 export(Calibration)
 export(Calibration_Michel)
 export(CreateCalibOptions)
+export(CreateErrorCrit_GAPX)
 export(CreateIniStates)
 export(CreateInputsCrit)
+export(CreateInputsCrit_Lavenne)
 export(CreateInputsModel)
 export(CreateRunOptions)
 export(DataAltiExtrapolation_Valery)
@@ -28,19 +36,24 @@ export(ErrorCrit_KGE)
 export(ErrorCrit_KGE2)
 export(ErrorCrit_NSE)
 export(ErrorCrit_RMSE)
-export(PEdaily_Oudin)
+export(Imax)
+export(PE_Oudin)
+export(plot.OutputsModel) ### to remove from version 2.0
 export(RunModel)
 export(RunModel_CemaNeige)
 export(RunModel_CemaNeigeGR4H)
+export(RunModel_CemaNeigeGR5H)
 export(RunModel_CemaNeigeGR4J)
 export(RunModel_CemaNeigeGR5J)
 export(RunModel_CemaNeigeGR6J)
 export(RunModel_GR1A)
 export(RunModel_GR2M)
 export(RunModel_GR4H)
+export(RunModel_GR5H)
 export(RunModel_GR4J)
 export(RunModel_GR5J)
 export(RunModel_GR6J)
+export(RunModel_Lag)
 export(SeriesAggreg)
 export(TransfoParam)
 export(TransfoParam_CemaNeige)
@@ -48,13 +61,13 @@ export(TransfoParam_CemaNeigeHyst)
 export(TransfoParam_GR1A)
 export(TransfoParam_GR2M)
 export(TransfoParam_GR4H)
+export(TransfoParam_GR5H)
 export(TransfoParam_GR4J)
 export(TransfoParam_GR5J)
 export(TransfoParam_GR6J)
-export(plot.OutputsModel)
-exportPattern(".FortranOutputs")
-exportPattern(".ErrorCrit")
-
+export(TransfoParam_Lag)
+#export(.ErrorCrit)
+#export(.FeatModels)
 
 
 #####################################

R/Calibration.R

-Calibration <- function(InputsModel, 
-                        RunOptions, 
-                        InputsCrit, 
+Calibration <- function(InputsModel,
+                        RunOptions,
+                        InputsCrit,
                         CalibOptions,
-                        FUN_MOD, 
+                        FUN_MOD,
                         FUN_CRIT,                      # deprecated
-                        FUN_CALIB = Calibration_Michel, 
+                        FUN_CALIB = Calibration_Michel,
                         FUN_TRANSFO = NULL,
-                        verbose = TRUE) {
-  
+                        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))
+                   verbose = verbose, ...))
 }
 

R/Calibration_Michel.R

-Calibration_Michel <- function(InputsModel, 
-                               RunOptions, 
-                               InputsCrit, 
+Calibration_Michel <- function(InputsModel,
+                               RunOptions,
+                               InputsCrit,
                                CalibOptions,
-                               FUN_MOD, 
+                               FUN_MOD,
                                FUN_CRIT,           # deprecated
-                               FUN_TRANSFO = NULL, 
-                               verbose = TRUE) {
-  
-  
+                               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'")
   }
@@ -46,100 +53,16 @@ Calibration_Michel <- function(InputsModel,
   }
   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")
-  }  
-  
-  
-  ##_check_FUN_TRANSFO
-  if (is.null(FUN_TRANSFO)) {
-    if (identical(FUN_MOD, RunModel_GR4H         )) {
-      FUN_TRANSFO <- TransfoParam_GR4H
-    }
-    if (identical(FUN_MOD, RunModel_GR4J         )) {
-      FUN_TRANSFO <- TransfoParam_GR4J
-    }
-    if (identical(FUN_MOD, RunModel_GR5J         )) {
-      FUN_TRANSFO <- TransfoParam_GR5J
-    }
-    if (identical(FUN_MOD, RunModel_GR6J         )) {
-      FUN_TRANSFO <- TransfoParam_GR6J
-    }
-    if (identical(FUN_MOD, RunModel_GR2M         )) {
-      FUN_TRANSFO <- TransfoParam_GR2M
-    }
-    if (identical(FUN_MOD, RunModel_GR1A         )) {
-      FUN_TRANSFO <- TransfoParam_GR1A
-    }
-    if (identical(FUN_MOD, RunModel_CemaNeige    )) {
-      if (inherits(CalibOptions, "hysteresis")) {
-        FUN_TRANSFO <- TransfoParam_CemaNeigeHyst
-      } else {
-        FUN_TRANSFO <- TransfoParam_CemaNeige
-      }
-    }
-    if (identical(FUN_MOD, RunModel_CemaNeigeGR4H) |
-        identical(FUN_MOD, RunModel_CemaNeigeGR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
-      if (identical(FUN_MOD, RunModel_CemaNeigeGR4H)) {
-        FUN1 <- TransfoParam_GR4H
-      }
-      if (identical(FUN_MOD, RunModel_CemaNeigeGR4J)) {
-        FUN1 <- TransfoParam_GR4J
-      }
-      if (identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
-        FUN1 <- TransfoParam_GR5J
-      }
-      if (identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
-        FUN1 <- TransfoParam_GR6J
-      }
-      if (inherits(CalibOptions, "hysteresis")) {
-        FUN2 <- TransfoParam_CemaNeigeHyst
-      } else {
-        FUN2 <- TransfoParam_CemaNeige
-      }
-      if (inherits(CalibOptions, "hysteresis")) {
-        FUN_TRANSFO <- function(ParamIn, Direction) {
-          Bool <- is.matrix(ParamIn)
-          if (!Bool) {
-            ParamIn <- rbind(ParamIn)
-          }
-          ParamOut <- NA * ParamIn
-          NParam   <- ncol(ParamIn)
-          ParamOut[,          1:(NParam-4)] <- FUN1(ParamIn[,          1:(NParam-4)], Direction)
-          ParamOut[, (NParam-3):NParam    ] <- FUN2(ParamIn[, (NParam-3):NParam    ], Direction)
-          if (!Bool) {
-            ParamOut <- ParamOut[1, ]
-          }
-          return(ParamOut)
-        }
-      } else {
-        FUN_TRANSFO <- function(ParamIn, Direction) {
-          Bool <- is.matrix(ParamIn)
-          if (!Bool) {
-            ParamIn <- rbind(ParamIn)
-          }
-          ParamOut <- NA * ParamIn
-          NParam   <- ncol(ParamIn)
-          ParamOut[,          1:(NParam-2)] <- FUN1(ParamIn[,          1:(NParam-2)], Direction)
-          ParamOut[, (NParam-1):NParam    ] <- FUN2(ParamIn[, (NParam-1):NParam    ], Direction)
-          if (!Bool) {
-            ParamOut <- ParamOut[1, ]
-          }
-          return(ParamOut)
-        }
-      }
-      
-    }
-    if (is.null(FUN_TRANSFO)) {
-      stop("'FUN_TRANSFO' was not found (in 'Calibration' function)")
-    }
   }
-  
-  ##_variables_initialisation 
+
+
+  ##_variables_initialisation
   ParamFinalR <- NULL
   ParamFinalT <- NULL
   CritFinal   <- NULL
@@ -168,20 +91,19 @@ Calibration_Michel <- function(InputsModel,
   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
+  ## use unique() to avoid duplicated values when a parameter is set
   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)
-  }    
-  
-  
+    expand.grid(lapply(seq_len(ncol(DistribParam)), function(x) unique(DistribParam[, x])))
+  }
+
+
   ##Creation_of_new_candidates_______________________________________________
   OptimParam <- is.na(CalibOptions$FixedParam)
   if (PrefilteringType == 1) {
@@ -190,7 +112,7 @@ Calibration_Michel <- function(InputsModel,
   if (PrefilteringType == 2) {
     DistribParamR <- CalibOptions$StartParamDistrib
     DistribParamR[, !OptimParam] <- NA
-    CandidatesParamR <- ProposeCandidatesGrid(DistribParamR)$NewCandidates
+    CandidatesParamR <- ProposeCandidatesGrid(DistribParamR)
   }
   ##Remplacement_of_non_optimised_values_____________________________________
   CandidatesParamR <- apply(CandidatesParamR, 1, function(x) {
@@ -202,7 +124,7 @@ Calibration_Michel <- function(InputsModel,
   } else {
     CandidatesParamR <- cbind(CandidatesParamR)
   }
-  
+
   ##Loop_to_test_the_various_candidates______________________________________
   iNewOptim <- 0
   Ncandidates <- nrow(CandidatesParamR)
@@ -221,12 +143,12 @@ Calibration_Michel <- function(InputsModel,
         if (iNew == round(k / 10 * Ncandidates)) {
           message(" ", 10 * k, "%", appendLF = FALSE)
         }
-      } 
+      }
     }
     ##Model_run
     Param <- CandidatesParamR[iNew, ]
-    OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
-    
+    OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, ...)
+
     ##Calibration_criterion_computation
     OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
     if (!is.na(OutputsCrit$CritValue)) {
@@ -245,8 +167,8 @@ Calibration_Michel <- function(InputsModel,
   if (verbose & Ncandidates > 1) {
     message(" 100%)\n", appendLF = FALSE)
   }
-  
-  
+
+
   ##End_of_first_step_Parameter_Screening____________________________________
   ParamStartR <- CandidatesParamR[iNewOptim, ]
   if (!is.matrix(ParamStartR)) {
@@ -269,13 +191,13 @@ Calibration_Michel <- function(InputsModel,
   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
@@ -299,10 +221,10 @@ Calibration_Michel <- function(InputsModel,
           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]
+            PotentialCandidateT[1, I] <- RangesT[1, I]
           }
           if (PotentialCandidateT[1, I] > RangesT[2, I]) {
-            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) {
@@ -326,11 +248,11 @@ Calibration_Michel <- function(InputsModel,
     Output$NewCandidatesT <- matrix(VECT, ncol = NParam, byrow = TRUE)
     return(Output)
   }
-  
-  
+
+
   ##Initialisation_of_variables
   if (verbose) {
-    message("Steepest-descent local search in progress") 
+    message("Steepest-descent local search in progress")
   }
   Pace <- 0.64
   PaceDiag <- rep(0, NParam)
@@ -342,18 +264,18 @@ Calibration_Michel <- function(InputsModel,
   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")
@@ -367,16 +289,16 @@ Calibration_Michel <- function(InputsModel,
     } else {
       CandidatesParamR <- cbind(CandidatesParamR)
     }
-    
-    
+
+
     ##Loop_to_test_the_various_candidates_____________________________________
     iNewOptim <- 0
     for (iNew in 1:nrow(CandidatesParamR)) {
       ##Model_run
       Param <- CandidatesParamR[iNew, ]
-      OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
+      OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, ...)
       ##Calibration_criterion_computation
-      OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)      
+      OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
       if (!is.na(OutputsCrit$CritValue)) {
         if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
           CritOptim <- OutputsCrit$CritValue * OutputsCrit$Multiplier
@@ -385,8 +307,8 @@ Calibration_Michel <- function(InputsModel,
       }
     }
     NRuns <- NRuns + nrow(CandidatesParamR)
-    
-    
+
+
     ##When_a_progress_has_been_achieved_______________________________________
     if (iNewOptim != 0) {
       ##We_store_the_optimal_set
@@ -401,7 +323,7 @@ Calibration_Michel <- function(InputsModel,
       ##We_update_PaceDiag
       VectPace <- NewParamOptimT-OldParamOptimT
       for (iC in 1:NParam) {
-        if (OptimParam[iC]) { 
+        if (OptimParam[iC]) {
           PaceDiag[iC] <- CLG * PaceDiag[iC] + (1-CLG) * VectPace[iC]
         }
       }
@@ -410,8 +332,8 @@ Calibration_Michel <- function(InputsModel,
       Pace <- Pace / 2
       Compt <- 0
     }
-    
-    
+
+
     ##Test_of_an_additional_candidate_using_diagonal_progress_________________
     if (ITER > 4 * NParam) {
       NRuns <- NRuns + 1
@@ -435,7 +357,7 @@ Calibration_Michel <- function(InputsModel,
       CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
       ##Model_run
       Param <- CandidatesParamR[iNew, ]
-      OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
+      OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, ...)
       ##Calibration_criterion_computation
       OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
       if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
@@ -447,34 +369,34 @@ Calibration_Michel <- function(InputsModel,
         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) { 
+  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) { 
+  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))
@@ -483,7 +405,7 @@ Calibration_Michel <- function(InputsModel,
       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)]
+      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, ")")
@@ -496,13 +418,13 @@ Calibration_Michel <- function(InputsModel,
   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,
@@ -511,8 +433,7 @@ Calibration_Michel <- function(InputsModel,
                        CritName = CritName, CritBestValue = CritBestValue)
   class(OutputsCalib) <- c("OutputsCalib", "HBAN")
   return(OutputsCalib)
-  
-  
-  
-}
 
+
+
+}

R/CreateErrorCrit_GAPX.R

+CreateErrorCrit_GAPX <- function(FUN_TRANSFO) {
+
+  FUN_CRIT <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
+    ## Arguments check
+    if (!inherits(OutputsModel, "OutputsModel")) {
+      stop("'OutputsModel' must be of class 'OutputsModel'")
+    }
+
+    OutputsModel$RunOptions$ParamT <- FUN_TRANSFO(OutputsModel$RunOptions$Param, "RT")
+
+    EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "GAPX", OutputsModel = OutputsModel, warnings = warnings)
+
+    CritValue <- NA
+
+    if (EC$CritCompute) {
+
+      ParamApr <- EC$VarObs[!EC$TS_ignore]
+      ParamOpt <- EC$VarSim[!EC$TS_ignore]
+
+      ## ErrorCrit
+      Crit <- 1 - sum(((ParamApr - ParamOpt) / 20)^2)^0.5
+
+      if (is.numeric(Crit) & is.finite(Crit)) {
+        CritValue <- Crit
+      }
+
+      ## Verbose
+      if (verbose) {
+        message(sprintf("Crit. %s = %.4f", EC$CritName, CritValue))
+      }
+    }
+
+    ## Output
+    OutputsCrit <- list(CritValue       = CritValue,
+                        CritName        = EC$CritName,
+                        CritBestValue   = EC$CritBestValue,
+                        Multiplier      = EC$Multiplier,
+                        Ind_notcomputed = EC$Ind_TS_ignore)
+
+    class(OutputsCrit) <- c("GAPX", "ErrorCrit")
+    return(OutputsCrit)
+  }
+
+  class(FUN_CRIT) <- c("FUN_CRIT", class(FUN_CRIT))
+
+  return(FUN_CRIT)
+}

R/CreateIniStates.R

-CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE,
-                            ProdStore = 350, RoutStore = 90, ExpStore = NULL,
+CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = FALSE,
+                            ProdStore = 350, RoutStore = 90, ExpStore = NULL, IntStore = NULL,
                             UH1 = NULL, UH2 = NULL,
                             GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
                             GthrCemaNeigeLayers = NULL, GlocmaxCemaNeigeLayers = NULL,
+                            SD = NULL,
                             verbose = TRUE) {
-  
-  
+
+
   ObjectClass <- NULL
-  
+
   UH1n <- 20L
   UH2n <- UH1n * 2L
- 
-  nameFUN_MOD <- as.character(substitute(FUN_MOD))
+
   FUN_MOD <- match.fun(FUN_MOD)
-  
-  ## check FUN_MOD
-  BOOL <- FALSE
-  if (identical(FUN_MOD, RunModel_GR4H)) {
-    ObjectClass <- c(ObjectClass, "GR", "hourly")
-    BOOL <- TRUE
-  }
-  if (identical(FUN_MOD, RunModel_GR4J) |
-      identical(FUN_MOD, RunModel_GR5J) |
-      identical(FUN_MOD, RunModel_GR6J)) {
-    ObjectClass <- c(ObjectClass, "GR", "daily")
-    BOOL <- TRUE
-  }
-  if (identical(FUN_MOD, RunModel_GR2M)) {
-    ObjectClass <- c(ObjectClass, "GR", "monthly")
-    BOOL <- TRUE
-  }
-  if (identical(FUN_MOD, RunModel_GR1A)) {
-    stop("'RunModel_GR1A' does not require 'IniStates' object")
-  }
-  if (identical(FUN_MOD, RunModel_CemaNeige)) {
-    ObjectClass <- c(ObjectClass, "CemaNeige", "daily")
-    BOOL <- TRUE
-  }
-  if (identical(FUN_MOD, RunModel_CemaNeigeGR4H)) {
-    ObjectClass <- c(ObjectClass, "GR", "CemaNeige", "hourly")
-    BOOL <- TRUE
-  }
-  if (identical(FUN_MOD, RunModel_CemaNeigeGR4J) |
-      identical(FUN_MOD, RunModel_CemaNeigeGR5J) |
-      identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
-    ObjectClass <- c(ObjectClass, "GR", "CemaNeige", "daily")
-    BOOL <- TRUE
-  }
-  if (!BOOL) {
-    stop("incorrect 'FUN_MOD' for use in 'CreateIniStates'")
-  }
+
+  FeatFUN_MOD <- .GetFeatModel(FUN_MOD = FUN_MOD, DatesR = InputsModel$DatesR)
+  ObjectClass <- FeatFUN_MOD$Class
+
   if (!"CemaNeige" %in% ObjectClass & IsHyst) {
     stop("'IsHyst' cannot be TRUE if CemaNeige is not used in 'FUN_MOD'")
   }
-  
+  if (!(identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) & IsIntStore) {
+    stop("'IsIntStore' cannot be TRUE if GR5H is not used in 'FUN_MOD'")
+  }
+
   ## check InputsModel
   if (!inherits(InputsModel, "InputsModel")) {
     stop("'InputsModel' must be of class 'InputsModel'")
@@ -65,118 +35,136 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE,
       !inherits(InputsModel, "CemaNeige")) {
     stop("'InputsModel' must be of class 'CemaNeige'")
   }
-  
-  
+
+
   ## check states
   if (any(eTGCemaNeigeLayers > 0)) {
     stop("Positive values are not allowed for 'eTGCemaNeigeLayers'")
-  }  
-  
+  }
+
   if (identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
     if (is.null(ExpStore)) {
       stop("'RunModel_*GR6J' need an 'ExpStore' value")
     }
   } else if (!is.null(ExpStore)) {
     if (verbose) {
-      warning(sprintf("'%s' does not require 'ExpStore'. Value set to NA", nameFUN_MOD))
+      warning(sprintf("'%s' does not require 'ExpStore'. Value set to NA", FeatFUN_MOD$NameFunMod))
     }
     ExpStore <- Inf
   }
-  
+
   if (identical(FUN_MOD, RunModel_GR2M)) {
     if (!is.null(UH1)) {
       if (verbose) {
-        warning(sprintf("'%s' does not require 'UH1'. Values set to NA", nameFUN_MOD))
+        warning(sprintf("'%s' does not require 'UH1'. Values set to NA", FeatFUN_MOD$NameFunMod))
       }
       UH1 <- rep(Inf, UH1n)
     }
     if (!is.null(UH2)) {
       if (verbose) {
-        warning(sprintf("'%s' does not require 'UH2'. Values set to NA", nameFUN_MOD))
+        warning(sprintf("'%s' does not require 'UH2'. Values set to NA", FeatFUN_MOD$NameFunMod))
       }
       UH2 <- rep(Inf, UH2n)
     }
   }
-  
+
   if ((identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J)) & !is.null(UH1)) {
     if (verbose) {
-      warning(sprintf("'%s' does not require 'UH1'. Values set to NA", nameFUN_MOD))
+      warning(sprintf("'%s' does not require 'UH1'. Values set to NA", FeatFUN_MOD$NameFunMod))
     }
     UH1 <- rep(Inf, UH1n)
   }
- 
+  if (!(identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) & !is.null(IntStore)) {
+    if (verbose) {
+      warning(sprintf("'%s' does not require 'IntStore'. Values set to NA", FeatFUN_MOD$NameFunMod))
+    }
+    IntStore <- Inf
+  }
+
   if ("CemaNeige" %in% ObjectClass & ! "GR" %in% ObjectClass) {
     if (!is.null(ProdStore)) {
       if (verbose) {
-        warning(sprintf("'%s' does not require 'ProdStore'. Values set to NA", nameFUN_MOD))
+        warning(sprintf("'%s' does not require 'ProdStore'. Values set to NA", FeatFUN_MOD$NameFunMod))
       }
     }
     ProdStore <- Inf
     if (!is.null(RoutStore)) {
       if (verbose) {
-        warning(sprintf("'%s' does not require 'RoutStore'. Values set to NA", nameFUN_MOD))
+        warning(sprintf("'%s' does not require 'RoutStore'. Values set to NA", FeatFUN_MOD$NameFunMod))
       }
     }
     RoutStore <- Inf
     if (!is.null(ExpStore)) {
       if (verbose) {
-        warning(sprintf("'%s' does not require 'ExpStore'. Values set to NA", nameFUN_MOD))
+        warning(sprintf("'%s' does not require 'ExpStore'. Values set to NA", FeatFUN_MOD$NameFunMod))
       }
     }
     ExpStore <- Inf
+    if (!is.null(IntStore)) {
+      if (verbose) {
+        warning(sprintf("'%s' does not require 'IntStore'. Values set to NA", FeatFUN_MOD$NameFunMod))
+      }
+    }
+    IntStore <- Inf
     if (!is.null(UH1)) {
       if (verbose) {
-        warning(sprintf("'%s' does not require 'UH1'. Values set to NA", nameFUN_MOD))
+        warning(sprintf("'%s' does not require 'UH1'. Values set to NA", FeatFUN_MOD$NameFunMod))
       }
     }
     UH1 <- rep(Inf, UH1n)
     if (!is.null(UH2)) {
       if (verbose) {
-        warning(sprintf("'%s' does not require 'UH2'. Values set to NA", nameFUN_MOD))
+        warning(sprintf("'%s' does not require 'UH2'. Values set to NA", FeatFUN_MOD$NameFunMod))
       }
     }
     UH2 <- rep(Inf, UH2n)
   }
-  if("CemaNeige" %in% ObjectClass & !IsHyst &
+  if (IsIntStore & is.null(IntStore)) {
+      stop(sprintf("'%s' need values for 'IntStore'", FeatFUN_MOD$NameFunMod))
+  }
+  if ("CemaNeige" %in% ObjectClass & !IsHyst &
      (is.null(GCemaNeigeLayers) | is.null(eTGCemaNeigeLayers))) {
-      stop(sprintf("'%s' need values for 'GCemaNeigeLayers' and 'GCemaNeigeLayers'", nameFUN_MOD))
+      stop(sprintf("'%s' need values for 'GCemaNeigeLayers' and 'GCemaNeigeLayers'", FeatFUN_MOD$NameFunMod))
   }
-  if("CemaNeige" %in% ObjectClass & IsHyst &
+  if ("CemaNeige" %in% ObjectClass & IsHyst &
      (is.null(GCemaNeigeLayers) | is.null(eTGCemaNeigeLayers) |
       is.null(GthrCemaNeigeLayers) | is.null(GlocmaxCemaNeigeLayers))) {
-    stop(sprintf("'%s' need values for 'GCemaNeigeLayers', 'GCemaNeigeLayers', 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'", nameFUN_MOD))
+    stop(sprintf("'%s' need values for 'GCemaNeigeLayers', 'GCemaNeigeLayers', 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'", FeatFUN_MOD$NameFunMod))
   }
-  if("CemaNeige" %in% ObjectClass & !IsHyst &
+  if ("CemaNeige" %in% ObjectClass & !IsHyst &
      (!is.null(GthrCemaNeigeLayers) | !is.null(GlocmaxCemaNeigeLayers))) {
     if (verbose) {
-      warning(sprintf("'%s' does not require 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'. Values set to NA", nameFUN_MOD))
+      warning(sprintf("'%s' does not require 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'. Values set to NA", FeatFUN_MOD$NameFunMod))
     }
     GthrCemaNeigeLayers    <- Inf
-    GlocmaxCemaNeigeLayers <- Inf 
+    GlocmaxCemaNeigeLayers <- Inf
   }
-  if(!"CemaNeige" %in% ObjectClass &
+  if (!"CemaNeige" %in% ObjectClass &
      (!is.null(GCemaNeigeLayers) | !is.null(eTGCemaNeigeLayers) | !is.null(GthrCemaNeigeLayers) | !is.null(GlocmaxCemaNeigeLayers))) {
     if (verbose) {
-      warning(sprintf("'%s' does not require 'GCemaNeigeLayers' 'GCemaNeigeLayers', 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'. Values set to NA", nameFUN_MOD))
+      warning(sprintf("'%s' does not require 'GCemaNeigeLayers' 'GCemaNeigeLayers', 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'. Values set to NA", FeatFUN_MOD$NameFunMod))
     }
     GCemaNeigeLayers       <- Inf
     eTGCemaNeigeLayers     <- Inf
     GthrCemaNeigeLayers    <- Inf
-    GlocmaxCemaNeigeLayers <- Inf    
+    GlocmaxCemaNeigeLayers <- Inf
   }
-  
-  
+
+
   ## set states
-  if("CemaNeige" %in% ObjectClass) {
+  if ("CemaNeige" %in% ObjectClass) {
     NLayers <- length(InputsModel$LayerPrecip)
   } else {
     NLayers <- 1
   }
-  
-  
+
+
   ## manage NULL values
   if (is.null(ExpStore)) {
-    ExpStore <- Inf 
+    ExpStore <- Inf
+  }
+  if (is.null(IntStore)) {
+    IntStore <- Inf
   }
   if (is.null(UH1)) {
     if ("hourly"  %in% ObjectClass) {
@@ -203,19 +191,24 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE,
   if (is.null(GthrCemaNeigeLayers)) {
     GthrCemaNeigeLayers <- rep(Inf, NLayers)
   }
+  if (any(is.infinite(GthrCemaNeigeLayers))) {
+    GthrCemaNeigeLayers <- rep(Inf, NLayers)
+  }
   if (is.null(GlocmaxCemaNeigeLayers)) {
     GlocmaxCemaNeigeLayers <- rep(Inf, NLayers)
   }
-  
-  
+  if (any(is.infinite(GlocmaxCemaNeigeLayers))) {
+    GlocmaxCemaNeigeLayers <- rep(Inf, NLayers)
+  }
+
   # check negative values
-  if (any(ProdStore < 0) | any(RoutStore < 0) |
+  if (any(ProdStore < 0) | any(RoutStore < 0) | any(IntStore < 0) |
       any(UH1 < 0) | any(UH2 < 0) |
       any(GCemaNeigeLayers < 0)) {
-    stop("Negative values are not allowed for any of 'ProdStore', 'RoutStore', 'UH1', 'UH2', 'GCemaNeigeLayers'")
+    stop("Negative values are not allowed for any of 'ProdStore', 'RoutStore', 'IntStore', 'UH1', 'UH2', 'GCemaNeigeLayers'")
   }
-  
-  
+
+
   ## check length
   if (!is.numeric(ProdStore) || length(ProdStore) != 1L) {
     stop("'ProdStore' must be numeric of length one")
@@ -226,6 +219,9 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE,
   if (!is.numeric(ExpStore) || length(ExpStore) != 1L) {
     stop("'ExpStore' must be numeric of length one")
   }
+  if (!is.numeric(IntStore) || length(IntStore) != 1L) {
+    stop("'IntStore' must be numeric of length one")
+  }
   if ( "hourly" %in% ObjectClass & (!is.numeric(UH1) || length(UH1) != UH1n * 24)) {
     stop(sprintf("'UH1' must be numeric of length 480 (%i * 24)", UH1n))
   }
@@ -252,23 +248,46 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE,
       stop(sprintf("'eTGCemaNeigeLayers' must be numeric of length %i", NLayers))
     }
   }
-  
+
+  # SD model state handling
+  if (!is.null(SD)) {
+    if (!inherits(InputsModel, "SD")) {
+      stop("'SD' argument provided and 'InputsModel' is not of class 'SD'")
+    }
+    if (!is.list(SD)) {
+      stop("'SD' argument must be a list")
+    }
+    lapply(SD, function(x) {
+      if (!is.numeric(x)) stop("Each item of 'SD' list argument must be numeric")
+    })
+    if (length(SD) != length(InputsModel$LengthHydro)) {
+      stop("Number of items of 'SD' list argument must be the same as the number of upstream connections",
+           sprintf(" (%i required, found %i)", length(InputsModel$LengthHydro), length(SD)))
+    }
+  }
 
   ## format output
-  IniStates <- list(Store = list(Prod = ProdStore, Rout = RoutStore, Exp = ExpStore),
+  IniStates <- list(Store = list(Prod = ProdStore, Rout = RoutStore, Exp = ExpStore, Int = IntStore),
                     UH = list(UH1 = UH1, UH2 = UH2),
                     CemaNeigeLayers = list(G = GCemaNeigeLayers, eTG = eTGCemaNeigeLayers,
                                            Gthr = GthrCemaNeigeLayers, Glocmax = GlocmaxCemaNeigeLayers))
   IniStatesNA <- unlist(IniStates)
   IniStatesNA[is.infinite(IniStatesNA)] <- NA
   IniStatesNA <- relist(IniStatesNA, skeleton = IniStates)
-  
+
+  if (!is.null(SD)) {
+    IniStatesNA$SD <- SD
+  }
+
   class(IniStatesNA) <- c("IniStates", ObjectClass)
-  if(IsHyst) {
+  if (IsHyst) {
     class(IniStatesNA) <- c(class(IniStatesNA), "hysteresis")
   }
-  
+  if (IsIntStore) {
+    class(IniStatesNA) <- c(class(IniStatesNA), "interception")
+  }
+
   return(IniStatesNA)
-  
-  
+
+
 }

R/CreateInputsCrit.R

 CreateInputsCrit <- function(FUN_CRIT,
                              InputsModel,
                              RunOptions,
-                             Qobs,              # deprecated
                              Obs,
                              VarObs = "Q",
                              BoolCrit = NULL,
                              transfo = "",
                              Weights = NULL,
-                             Ind_zeroes = NULL, # deprecated
                              epsilon = NULL,
-                             warnings = TRUE,
-                             verbose = TRUE) {
-  
-  
+                             warnings = TRUE) {
+
+
   ObjectClass <- NULL
-  
-  
+
+
   ## ---------- check arguments
-  
-  if (!missing(Qobs)) {
-    if (missing(Obs)) {
-      if (warnings) {
-        warning("argument 'Qobs' is deprecated. Please use 'Obs' and 'VarObs' instead")
-      }
-      Obs <- Qobs
-      # VarObs <- "Qobs"
-    } else {
-      warning("argument 'Qobs' is deprecated. The values set in 'Obs' will be used instead")
-    }
-  }
-  if (!missing(Ind_zeroes) & warnings) {
-    warning("deprecated 'Ind_zeroes' argument")
-  }
-  if (!missing(verbose)) {
-    warning("deprecated 'verbose' argument. Use 'warnings', instead")
-  }
-  
-  
+
   ## check 'InputsModel'
   if (!inherits(InputsModel, "InputsModel")) {
     stop("'InputsModel' must be of class 'InputsModel'")
   }
-  
-  
+
+
   ## length of index of period to be used for the model run
   LLL <- length(InputsModel$DatesR[RunOptions$IndPeriod_Run])
-  
-  
+
+
   ## check 'Obs' and definition of idLayer
-  vecObs <- unlist(Obs)
-  if (length(vecObs) %% LLL != 0 | !is.numeric(vecObs)) {
-    stop(sprintf("'Obs' must be a (list of) vector(s) of numeric values of length %i", LLL), call. = FALSE)
+  if (!is.numeric(unlist(Obs))) {
+    stop("'Obs' must be a (list of) vector(s) of numeric values")
+  }
+  Obs2 <- Obs
+  if ("ParamT" %in% VarObs) {
+    if (is.list(Obs2)) {
+      Obs2[[which(VarObs == "ParamT")]] <- NULL
+    } else {
+      Obs2 <- NULL
+    }
+  }
+  if (!is.null(Obs2)) {
+    vecObs <- unlist(Obs2)
+    if (length(vecObs) %% LLL != 0) {
+      stop(sprintf("'Obs' must be a (list of) vector(s) of numeric values of length %i", LLL), call. = FALSE)
+    }
   }
   if (!is.list(Obs)) {
     idLayer <- list(1L)
@@ -65,8 +56,8 @@ CreateInputsCrit <- function(FUN_CRIT,
     })
     Obs <- lapply(Obs, function(x) rowMeans(as.data.frame(x)))
   }
-  
-  
+
+
   ## create list of arguments
   listArgs <- list(FUN_CRIT   = FUN_CRIT,
                    Obs        = Obs,
@@ -76,8 +67,8 @@ CreateInputsCrit <- function(FUN_CRIT,
                    transfo    = as.character(transfo),
                    Weights    = Weights,
                    epsilon    = epsilon)
-  
-  
+
+
   ## check lists lengths
   for (iArgs in names(listArgs)) {
     if (iArgs %in% c("Weights", "BoolCrit", "epsilon")) {
@@ -92,11 +83,11 @@ CreateInputsCrit <- function(FUN_CRIT,
       listArgs[[iArgs]] <- list(listArgs[[iArgs]])
     }
   }
-  
+
   ## check 'FUN_CRIT'
   listArgs$FUN_CRIT <- lapply(listArgs$FUN_CRIT, FUN = match.fun)
-  
-  
+
+
   ## check 'VarObs'
   if (missing(VarObs)) {
     listArgs$VarObs <- as.list(rep("Q", times = length(listArgs$Obs)))
@@ -104,8 +95,8 @@ CreateInputsCrit <- function(FUN_CRIT,
     #   warning("'VarObs' automatically set to \"Q\"")
     # }
   }
-  
-  
+
+
   ## check 'VarObs' + 'RunOptions'
   if ("Q" %in% VarObs & !inherits(RunOptions, "GR")) {
     stop("'VarObs' cannot contain Q if a GR rainfall-runoff model is not used")
@@ -119,67 +110,76 @@ CreateInputsCrit <- function(FUN_CRIT,
   if ("SWE" %in% VarObs & inherits(RunOptions, "CemaNeige") & !"SnowPack" %in% RunOptions$Outputs_Sim) {
     stop("'SnowPack' is missing in 'Outputs_Sim' of 'RunOptions', which is necessary to output SWE with CemaNeige")
   }
-  
-  
+
+
   ## check 'transfo'
   if (missing(transfo)) {
     listArgs$transfo <- as.list(rep("", times = length(listArgs$Obs)))
     # if (warnings) {
     #   warning("'transfo' automatically set to \"\"")
     # }
-  }  
-  
+  }
+
   ## check length of each args
   if (length(unique(sapply(listArgs, FUN = length))) != 1) {
     stopListArgs <- paste(sapply(names(listArgs), shQuote), collapse = ", ")
     stop(sprintf("arguments %s must have the same length", stopListArgs))
   }
-  
-  
+
+
   ## check 'RunOptions'
   if (!inherits(RunOptions , "RunOptions")) {
     stop("'RunOptions' must be of class 'RunOptions'")
   }
-  
-  
+
+
   ## check 'Weights'
   if (length(listArgs$Weights) > 1 & sum(unlist(listArgs$Weights)) == 0 & !any(sapply(listArgs$Weights, is.null))) {
     stop("sum of 'Weights' cannot be equal to zero")
   }
-  
-  
+
+
   ## ---------- reformat
-  
+
   ## reformat list of arguments
   listArgs2 <- lapply(seq_along(listArgs$FUN_CRIT), function(i) lapply(listArgs, "[[", i))
-  
+
   ## preparation of warning messages
-  inVarObs  <- c("Q", "SCA", "SWE")
+  inVarObs  <- c("Q", "SCA", "SWE", "ParamT")
   msgVarObs <- "'VarObs' must be a (list of) character vector(s) and one of %s"
   msgVarObs <- sprintf(msgVarObs, paste(sapply(inVarObs, shQuote), collapse = ", "))
   inTransfo  <- c("", "sqrt", "log", "inv", "sort", "boxcox") # pow is not checked by inTransfo, but appears in the warning message and checkef after (see ## check 'transfo')
   msgTransfo <- "'transfo' must be a (list of) character vector(s) and one of %s, or numeric value for power transformation"
   msgTransfo <- sprintf(msgTransfo, paste(sapply(inTransfo, shQuote), collapse = ", "))
-  
-  
+
+
   ## ---------- loop on the list of inputs
-  
+
   InputsCrit <- lapply(listArgs2, function(iListArgs2) {
-    
+
+    ## define FUN_CRIT as a character string
+    iListArgs2$FUN_CRIT <- match.fun(iListArgs2$FUN_CRIT)
+
     ## check 'FUN_CRIT'
-    if (!(identical(iListArgs2$FUN_CRIT, ErrorCrit_NSE ) | identical(iListArgs2$FUN_CRIT, ErrorCrit_KGE ) |
-          identical(iListArgs2$FUN_CRIT, ErrorCrit_KGE2) | identical(iListArgs2$FUN_CRIT, ErrorCrit_RMSE))) {
+    if (!all(class(iListArgs2$FUN_CRIT) == c("FUN_CRIT", "function"))) {
       stop("incorrect 'FUN_CRIT' for use in 'CreateInputsCrit'", call. = FALSE)
     }
     if (identical(iListArgs2$FUN_CRIT, ErrorCrit_RMSE) & length(listArgs$Weights) > 1 & all(!is.null(unlist(listArgs$Weights)))) {
       stop("calculating a composite criterion with the RMSE is not allowed since RMSE is not a dimensionless metric", call. = FALSE)
     }
-    
+
     ## check 'Obs'
-    if (!is.vector(iListArgs2$Obs) | length(iListArgs2$Obs) != LLL | !is.numeric(iListArgs2$Obs)) {
-      stop(sprintf("'Obs' must be a (list of) vector(s) of numeric values of length %i", LLL), call. = FALSE)
+    if (iListArgs2$VarObs == "ParamT") {
+      # Parameter for regularisation
+      L2 <- RunOptions$FeatFUN_MOD$NbParam
+    } else {
+      # Observation time series
+      L2 <- LLL
     }
-    
+    if (!is.vector(iListArgs2$Obs) | length(iListArgs2$Obs) != L2 | !is.numeric(iListArgs2$Obs)) {
+      stop(sprintf("'Obs' must be a (list of) vector(s) of numeric values of length %i", L2), call. = FALSE)
+    }
+
     ## check 'BoolCrit'
     if (is.null(iListArgs2$BoolCrit)) {
       iListArgs2$BoolCrit <- rep(TRUE, length(iListArgs2$Obs))
@@ -187,15 +187,15 @@ CreateInputsCrit <- function(FUN_CRIT,
     if (!is.logical(iListArgs2$BoolCrit)) {
       stop("'BoolCrit' must be a (list of) vector(s) of boolean", call. = FALSE)
     }
-    if (length(iListArgs2$BoolCrit) != LLL) {
-      stop("'BoolCrit' and 'InputsModel' series must have the same length", call. = FALSE)
+    if (length(iListArgs2$BoolCrit) != L2) {
+      stop("'BoolCrit' and the period defined in 'RunOptions' must have the same length", call. = FALSE)
     }
-    
+
     ## check 'VarObs'
     if (!is.vector(iListArgs2$VarObs) | length(iListArgs2$VarObs) != 1 | !is.character(iListArgs2$VarObs) | !all(iListArgs2$VarObs %in% inVarObs)) {
       stop(msgVarObs, call. = FALSE)
     }
-    
+
     ## check 'VarObs' + 'Obs'
     if (any(iListArgs2$VarObs %in% "SCA")) {
       idSCA <- which(iListArgs2$VarObs == "SCA")
@@ -207,7 +207,7 @@ CreateInputsCrit <- function(FUN_CRIT,
       if (min(vecSCA, na.rm = TRUE) < 0 | max(vecSCA, na.rm = TRUE) > 1) {
         stop("'Obs' outside [0,1] for \"SCA\"", call. = FALSE)
       }
-    } 
+    }
     inPosVarObs <- c("Q", "SWE")
     if (any(iListArgs2$VarObs %in% inPosVarObs)) {
       idQSS <- which(iListArgs2$VarObs %in% inPosVarObs)
@@ -223,8 +223,8 @@ CreateInputsCrit <- function(FUN_CRIT,
         stop(sprintf("'Obs' outside [0,Inf[ for \"%s\"", iListArgs2$VarObs), call. = FALSE)
       }
     }
-    
-    
+
+
     ## check 'transfo'
     if (is.null(iListArgs2$transfo) | !is.vector(iListArgs2$transfo) | length(iListArgs2$transfo) != 1 | !is.character(iListArgs2$transfo)) {
       stop(msgTransfo, call. = FALSE)
@@ -239,14 +239,14 @@ CreateInputsCrit <- function(FUN_CRIT,
       }
       iListArgs2$transfo <- paste0("^", iListArgs2$transfo)
     }
-    
+
     ## check 'Weights'
     if (!is.null(iListArgs2$Weights)) {
       if (!is.vector(iListArgs2$Weights) | length(iListArgs2$Weights) != 1 | !is.numeric(iListArgs2$Weights) | any(iListArgs2$Weights < 0)) {
         stop("'Weights' must be a single (list of) positive or equal to zero value(s)", call. = FALSE)
       }
     }
-    
+
     ## check 'epsilon'
     if (!is.null(iListArgs2$epsilon)) {
       if (!is.vector(iListArgs2$epsilon) | length(iListArgs2$epsilon) != 1 | !is.numeric(iListArgs2$epsilon) | any(iListArgs2$epsilon <= 0)) {
@@ -255,7 +255,7 @@ CreateInputsCrit <- function(FUN_CRIT,
     } else if (iListArgs2$transfo %in% c("log", "inv") & any(iListArgs2$Obs %in% 0) & warnings) {
       warning("zeroes detected in Obs: the corresponding time-steps will be excluded by the 'ErrorCrit*' functions as the epsilon argument was set to NULL", call. = FALSE)
     }
-    
+
     ## check 'transfo' + 'FUN_CRIT'
     if (iListArgs2$transfo == "log" & warnings) {
       warn_log_kge <- "we do not advise using the %s with a log transformation on Obs (see the details section in the 'CreateInputsCrit' help)"
@@ -266,7 +266,7 @@ CreateInputsCrit <- function(FUN_CRIT,
         warning(sprintf(warn_log_kge, "KGE'"), call. = FALSE)
       }
     }
-    
+
 
     ## Create InputsCrit
     iInputsCrit <- list(FUN_CRIT   = iListArgs2$FUN_CRIT,
@@ -279,21 +279,14 @@ CreateInputsCrit <- function(FUN_CRIT,
                         Weights    = iListArgs2$Weights)
     class(iInputsCrit) <- c("Single", "InputsCrit", ObjectClass)
     return(iInputsCrit)
-    
+
   })
   names(InputsCrit) <- paste0("IC", seq_along(InputsCrit))
-  
-  ## define FUN_CRIT as a characater string
-  listErrorCrit <- c("ErrorCrit_KGE", "ErrorCrit_KGE2", "ErrorCrit_NSE",  "ErrorCrit_RMSE")
-  InputsCrit <- lapply(InputsCrit, function(i) {
-    i$FUN_CRIT <- listErrorCrit[sapply(listErrorCrit, function(j) identical(i$FUN_CRIT, get(j)))]
-    i
-    })
-  
+
   listVarObs <- sapply(InputsCrit, FUN = "[[", "VarObs")
   inCnVarObs <- c("SCA", "SWE")
   if (!"ZLayers" %in% names(InputsModel)) {
-    if(any(listVarObs %in% inCnVarObs)) {
+    if (any(listVarObs %in% inCnVarObs)) {
       stop(sprintf("'VarObs' can not be equal to %i if CemaNeige is not used",
                    paste(sapply(inCnVarObs, shQuote), collapse = " or ")))
     }
@@ -311,10 +304,10 @@ CreateInputsCrit <- function(FUN_CRIT,
       }
     }
   }
-  
+
   ## define idLayer as an index of the layer to use
   for (iInCnVarObs in unique(listVarObs)) {
-    if (iInCnVarObs == "Q") {
+    if (!iInCnVarObs %in% inCnVarObs) {
       for (i in which(listVarObs == iInCnVarObs)) {
         InputsCrit[[i]]$idLayer <- NA
       }
@@ -330,8 +323,8 @@ CreateInputsCrit <- function(FUN_CRIT,
       }
     }
   }
-  
-  
+
+
   ## if only one criterion --> not a list of InputsCrit but directly an InputsCrit
   if (length(InputsCrit) < 2) {
     InputsCrit <- InputsCrit[[1L]]
@@ -348,12 +341,12 @@ CreateInputsCrit <- function(FUN_CRIT,
     combInputsCrit <- combn(x = length(InputsCrit), m = 2)
     apply(combInputsCrit, MARGIN = 2, function(i) {
       equalInputsCrit <- identical(InputsCrit[[i[1]]], InputsCrit[[i[2]]])
-      if(equalInputsCrit) {
+      if (equalInputsCrit) {
         warning(sprintf("elements %i and %i of the criteria list are identical. This might not be necessary", i[1], i[2]), call. = FALSE)
       }
     })
   }
-  
+
   return(InputsCrit)
-  
+
 }

R/CreateInputsCrit_Lavenne.R

+CreateInputsCrit_Lavenne <- function(FUN_CRIT = ErrorCrit_KGE,
+                                      InputsModel,
+                                      RunOptions,
+                                      Obs,
+                                      VarObs = "Q",
+                                      AprParamR,
+                                      AprCrit = 1,
+                                      k = 0.15,
+                                      BoolCrit = NULL,
+                                      transfo = "sqrt",
+                                      epsilon = NULL) {
+
+  # Check parameters
+  if (!is.numeric(AprCrit) || length(AprCrit) != 1 || AprCrit > 1) {
+    stop("'AprCrit' must be a numeric of length 1 with a maximum value of 1")
+  }
+  if (!is.numeric(k) || length(k) != 1 || k < 0 || k > 1) {
+    stop("'k' must be a numeric of length 1 with a value between 0 and 1")
+  }
+  if (!is.null(BoolCrit) && !is.logical(BoolCrit)) {
+    stop("'BoolCrit must be logical")
+  }
+  if (!is.character(transfo)) {
+    stop("'transfo' must be character")
+  }
+  if (!is.null(epsilon) && !is.numeric(epsilon)) {
+    stop("'epsilon' must be numeric")
+  }
+  if (!is.numeric(AprParamR) || length(AprParamR) != RunOptions$FeatFUN_MOD$NbParam) {
+    stop("'AprParamR' must be a numeric vector of length ",
+         RunOptions$FeatFUN_MOD$NbParam)
+  }
+
+
+  FUN_TRANSFO <- .FunTransfo(RunOptions$FeatFUN_MOD)
+
+  AprParamT <- FUN_TRANSFO(AprParamR, "RT")
+
+  ErrorCrit_GAPX <- CreateErrorCrit_GAPX(FUN_TRANSFO)
+
+
+
+  CreateInputsCrit(FUN_CRIT = list(FUN_CRIT, ErrorCrit_GAPX),
+                   InputsModel = InputsModel,
+                   RunOptions = RunOptions,
+                   Obs = list(Obs, AprParamT),
+                   VarObs = c("Q", "ParamT"),
+                   Weights = c(1 - k, k * max(0, AprCrit)),
+                   BoolCrit = list(BoolCrit, NULL),
+                   transfo = list(transfo, ""),
+                   epsilon = list(epsilon, NULL))
+}

R/CreateRunOptions.R

-CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndPeriod_Run,
-                             IniStates = NULL, IniResLevels = NULL, 
+CreateRunOptions <- function(FUN_MOD, InputsModel,
+                             IndPeriod_WarmUp = NULL, IndPeriod_Run,
+                             IniStates = NULL, IniResLevels = NULL, Imax = NULL,
                              Outputs_Cal = NULL, Outputs_Sim = "all",
-                             RunSnowModule, MeanAnSolidPrecip = NULL, IsHyst = FALSE, 
+                             MeanAnSolidPrecip = NULL, IsHyst = FALSE,
                              warnings = TRUE, verbose = TRUE) {
-  
-  if (!missing(RunSnowModule)) {
-    warning("deprecated 'RunSnowModule' argument: please adapt 'FUN_MOD' instead.", call. = FALSE)
-  }
-  if (!is.logical(IsHyst) | length(IsHyst) != 1L) {
-    stop("'IsHyst' must be a 'logical' of length 1")
+
+  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
   }
-  ObjectClass <- NULL
-  
+
+  ## check FUN_MOD
   FUN_MOD <- match.fun(FUN_MOD)
-  
-  ##check_FUN_MOD
-  BOOL <- FALSE
-  if (identical(FUN_MOD, RunModel_GR4H)) {
-    ObjectClass <- c(ObjectClass, "GR", "hourly")
-    BOOL <- TRUE
-  }
-  if (identical(FUN_MOD, RunModel_GR4J) | identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_GR6J)) {
-    ObjectClass <- c(ObjectClass, "GR", "daily")
-    BOOL <- TRUE
-  }
-  if (identical(FUN_MOD, RunModel_GR2M)) {
-    ObjectClass <- c(ObjectClass, "GR", "monthly")
-    BOOL <- TRUE
-  }
-  if (identical(FUN_MOD, RunModel_GR1A)) {
-    ObjectClass <- c(ObjectClass, "GR", "yearly")
-    BOOL <- TRUE
-  }
-  if (identical(FUN_MOD, RunModel_CemaNeige)) {
-    ObjectClass <- c(ObjectClass, "CemaNeige", "daily")
-    BOOL <- TRUE
+  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 (identical(FUN_MOD, RunModel_CemaNeigeGR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
-    ObjectClass <- c(ObjectClass, "GR", "CemaNeige", "daily")
-    BOOL <- TRUE
+  if ("CemaNeige" %in% FeatFUN_MOD$Class) {
+    FeatFUN_MOD$IsHyst <- IsHyst
+    if (IsHyst) {
+      ObjectClass <- c(ObjectClass, "hysteresis")
+      FeatFUN_MOD$NbParam <- FeatFUN_MOD$NbParam + 2
+    }
   }
-  if (identical(FUN_MOD, RunModel_CemaNeigeGR4H)) {
-    ObjectClass <- c(ObjectClass, "GR", "CemaNeige", "hourly")
-    BOOL <- TRUE
+
+  ## SD model
+  FeatFUN_MOD$IsSD <- inherits(InputsModel, "SD")
+  if (FeatFUN_MOD$IsSD) {
+    FeatFUN_MOD$NbParam <- FeatFUN_MOD$NbParam + 1
   }
-  if (IsHyst) {
-    ObjectClass <- c(ObjectClass, "hysteresis")
+
+  if (!"CemaNeige" %in% ObjectClass & "hysteresis" %in% ObjectClass) {
+    stop("'IsHyst' cannot be TRUE for the chosen 'FUN_MOD'")
   }
-  if (!BOOL) {
-    stop("incorrect 'FUN_MOD' for use in 'CreateRunOptions'")
+  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'")
@@ -77,8 +79,8 @@ CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndP
       !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")
@@ -92,15 +94,15 @@ CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndP
   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 <- as.integer(0)
+      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 {
@@ -112,19 +114,7 @@ CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndP
         TmpDateR <- TmpDateR - 1 * 24 * 60 * 60
       }
       IndPeriod_WarmUp <- which(InputsModel$DatesR == max(InputsModel$DatesR[1], TmpDateR)):(IndPeriod_Run[1] - 1)
-      if ("hourly"  %in% ObjectClass) {
-        TimeStep <- as.integer(60 * 60)
-      }
-      if ("daily"   %in% ObjectClass) {
-        TimeStep <- as.integer(24 * 60 * 60)
-      }
-      if ("monthly" %in% ObjectClass) {
-        TimeStep <- as.integer(30.44 * 24 * 60 * 60)
-      }
-      if ("yearly"  %in% ObjectClass) {
-        TimeStep <- as.integer(365.25 * 24 * 60 * 60)
-      }
-      if (length(IndPeriod_WarmUp) * TimeStep / (365 * 24 * 60 * 60) >= 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:")
@@ -142,41 +132,79 @@ CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndP
     if (storage.mode(IndPeriod_WarmUp) != "integer") {
       stop("'IndPeriod_WarmUp' should be of type integer")
     }
-    if (identical(IndPeriod_WarmUp, as.integer(0)) & verbose) {
+    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, as.integer(0))) {
+    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    
+
+
+  ## 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))) {
-        stop("'IniResLevels' must be a vector of numeric values")
+      # 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_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 ((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_GR6J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)) &
-          length(IniResLevels) != 3) {
-        stop("the length of 'IniResLevels' must be 3 for the chosen 'FUN_MOD'")
+      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))
-      } else {
-        IniResLevels <- as.double(c(0.3, 0.5, NA))
+        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)) {
@@ -198,7 +226,7 @@ CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndP
   NState <- NULL
   if ("GR" %in% ObjectClass | "CemaNeige" %in% ObjectClass) {
     if ("hourly"  %in% ObjectClass) {
-      NState <- 7 + 3 * 24 * 20+ 4 * NLayers
+      NState <- 7 + 3 * 24 * 20 + 4 * NLayers
     }
     if ("daily"   %in% ObjectClass) {
       NState <- 7 + 3 * 20 + 4 * NLayers
@@ -211,7 +239,7 @@ CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndP
     }
   }
   if (!is.null(IniStates)) {
-    
+
     if (!inherits(IniStates, "IniStates")) {
       stop("'IniStates' must be an object of class 'IniStates'")
     }
@@ -221,25 +249,31 @@ CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndP
     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)) & !all(is.na(IniStates$UH$UH1))) { ## GR5J
+    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")) | 
+    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 (!"CemaNeige" %in% ObjectClass & "hysteresis" %in% ObjectClass) {
-      stop("'IsHyst' cannot be TRUE for the chosen 'FUN_MOD'")
-    }
     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")
@@ -256,7 +290,7 @@ CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndP
     if (!"CemaNeige" %in% ObjectClass & any(is.na(IniStates$CemaNeigeLayers$Glocmax))) {
       IniStates$CemaNeigeLayers$Glocmax <- NULL
     }
-    IniStates$Store$Rest <- rep(NA, 4)
+    IniStates$Store$Rest <- rep(NA, 3)
     IniStates <- unlist(IniStates)
     IniStates[is.na(IniStates)] <- 0
     if ("monthly" %in% ObjectClass) {
@@ -265,34 +299,16 @@ CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndP
   } else {
     IniStates <- as.double(rep(0.0, NState))
   }
-  
-  
+
+
   ##check_Outputs_Cal_and_Sim
-  
+
   ##Outputs_all
-  Outputs_all <- NULL
-  if (identical(FUN_MOD,RunModel_GR4H) | identical(FUN_MOD,RunModel_CemaNeigeGR4H)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR4H")$GR)
-  }
-  if (identical(FUN_MOD,RunModel_GR4J) | identical(FUN_MOD,RunModel_CemaNeigeGR4J)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR4J")$GR)
+  Outputs_all <- c("DatesR", unlist(FortranOutputs), "WarmUpQsim", "StateEnd", "Param")
+  if (FeatFUN_MOD$IsSD) {
+    Outputs_all <- c(Outputs_all, "QsimDown", "Qsim_m3")
   }
-  if (identical(FUN_MOD,RunModel_GR5J) | identical(FUN_MOD,RunModel_CemaNeigeGR5J)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR5J")$GR)
-  }
-  if (identical(FUN_MOD,RunModel_GR6J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR6J")$GR)
-  }
-  if (identical(FUN_MOD,RunModel_GR2M)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR2M")$GR)
-  }
-  if (identical(FUN_MOD,RunModel_GR1A)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR1A")$GR)
-  }
-  if ("CemaNeige" %in% ObjectClass) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = NULL, isCN = TRUE)$CN)
-  }
-  
+
   ##check_Outputs_Sim
   if (!is.vector(Outputs_Sim)) {
     stop("'Outputs_Sim' must be a vector of characters")
@@ -304,28 +320,26 @@ CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndP
     stop("'Outputs_Sim' must not contain NA")
   }
   if ("all" %in% Outputs_Sim) {
-    Outputs_Sim <- c("DatesR", Outputs_all, "StateEnd")
+    Outputs_Sim <- Outputs_all
   }
-  Test <- which(!Outputs_Sim %in% c("DatesR", Outputs_all, "StateEnd"))
+  Test <- which(!Outputs_Sim %in% Outputs_all)
   if (length(Test) != 0) {
     stop(paste0( "'Outputs_Sim' is incorrectly defined: ",
                  paste(Outputs_Sim[Test], collapse = ", "), " not found"))
   }
   Outputs_Sim <- Outputs_Sim[!duplicated(Outputs_Sim)]
-  
-  
+
+
   ##check_Outputs_Cal
   if (is.null(Outputs_Cal)) {
     if ("GR" %in% ObjectClass) {
-      Outputs_Cal <- c("Qsim")
-    }
-    if ("CemaNeige" %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")
     }
-    if ("GR" %in% ObjectClass &
-        "CemaNeige" %in% ObjectClass) {
-      Outputs_Cal <- c("PliqAndMelt", "Qsim")
-    }
   } else {
     if (!is.vector(Outputs_Cal)) {
       stop("'Outputs_Cal' must be a vector of characters")
@@ -338,17 +352,16 @@ CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndP
     }
   }
   if ("all" %in% Outputs_Cal) {
-    Outputs_Cal <- c("DatesR", Outputs_all, "StateEnd")
-    
+    Outputs_Cal <- Outputs_all
   }
-  Test <- which(!Outputs_Cal %in% c("DatesR", Outputs_all, "StateEnd"))
+  Test <- which(!Outputs_Cal %in% Outputs_all)
   if (length(Test) != 0) {
     stop(paste0("'Outputs_Cal' is incorrectly defined: ",
                 paste(Outputs_Cal[Test], collapse = ", "), " not found"))
   }
   Outputs_Cal <- unique(Outputs_Cal)
-  
-  
+
+
   ##check_MeanAnSolidPrecip
   if ("CemaNeige" %in% ObjectClass & is.null(MeanAnSolidPrecip)) {
     NLayers <- length(InputsModel$LayerPrecip)
@@ -397,8 +410,8 @@ CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndP
       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) {
@@ -420,8 +433,8 @@ CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndP
       Outputs_Sim <- c(Outputs_Sim, "PliqAndMelt")
     }
   }
-  
-  
+
+
   ##check_Qsim
   if ("GR" %in% ObjectClass) {
     if (!"Qsim" %in% Outputs_Cal & !"all" %in% Outputs_Cal) {
@@ -443,22 +456,27 @@ CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndP
       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)
-  
+                     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)
-  
-  
+
+
 }