.Rbuildignore

 ^.*\.Rproj$
 ^\.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
+
+# 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.0.0
-Date: 2016-04-14
+Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
+Version: 1.7.6.9000
+Date: 2023-10-25
 Authors@R: c(
-  person("Laurent", "Coron", role = c("aut", "cre", "trl")),
-  person("Charles", "Perrin", role = c("aut", "ths")),
+  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@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")),
-  person("Pierre", "Brigode", role = c("ctb")),
-  person("Olivier", "Delaigue", role = c("ctb")),
+  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")),
+  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("Guillaume", "Thirel", role = c("ctb")),
   person("Audrey", "Valéry", role = c("ctb"))
   )
-Author: Laurent Coron, Charles Perrin, with contributions
- from Vazken Andréassian, Pierre Brigode, Olivier Delaigue,
- Nicolas Le Moine, Thibaut Mathevet, Safouane Mouelhi,
- Ludovic Oudin, Raji Pushpalatha, Guillaume Thirel, Audrey Valéry.
- Based on earlier work by Claude Michel.
-Maintainer: Laurent Coron, Olivier Delaigue <airGR@irstea.fr>
-Depends: R (>= 3.0.1)
-Description: This package brings into R the hydrological modelling tools developed
-  at Irstea-Antony (HBAN Research Unit, France). The package includes several conceptual
-  rainfall-runoff models (GR4H, GR4J, GR5J, GR6J, GR2M, GR1A), a snowmelt module (Cemaneige)
-  and the associated functions for their calibration and evaluation. Use help(airGR) for package description.
+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: http://webgr.irstea.fr/modeles/?lang=en
+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

 #####################################
 ##              Load DLL           ##
 #####################################
-useDynLib(airGR)
+useDynLib(airGR, .registration = TRUE)
+
+
+
+#####################################
+##            S3 methods           ##
+#####################################
+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)
 
 
 
@@ -10,9 +23,11 @@ useDynLib(airGR)
 #####################################
 export(Calibration)
 export(Calibration_Michel)
-export(Calibration_optim)
 export(CreateCalibOptions)
+export(CreateErrorCrit_GAPX)
+export(CreateIniStates)
 export(CreateInputsCrit)
+export(CreateInputsCrit_Lavenne)
 export(CreateInputsModel)
 export(CreateRunOptions)
 export(DataAltiExtrapolation_Valery)
@@ -21,26 +36,44 @@ 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)
+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)
+export(TransfoParam_Lag)
+#export(.ErrorCrit)
+#export(.FeatModels)
+
 
+#####################################
+##               Import            ##
+#####################################
+import(stats)
+import(graphics)
+import(grDevices)
+import(utils)

NEWS.md

+## Release History of the airGR Package
+
+
+
+### 1.7.6 Release Notes (2023-10-25)
+
+#### Bug fixes
+
+- `CreateCalibOptions()` now uses parameter screening for `RunModel_Lag()` which are now expressed in the transformed space instead of the parameter space. ([#156](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/156))
+- `RunModel_CemaNeige*()` now takes into account the case when `dG = 0`. ([#178](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/178))
+
+
+#### Minor user-visible changes
+
+- `Calibration_Michel()` now runs faster as the `ProposeCandidatesGrid()` was improved to create the propose candidates grid. ([#157](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/157)) 
+- `TransfoParamGR5J()` now returns the correct  error message when the number of parameters is incorrect. ([#168](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/168)) 
+
+
+#### CRAN-compatibility updates
+
+- `frun_*` Fortran subroutine does not use anymore the 'DLLEXPORT' command. ([#180](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/180))
+- The 'Rmalschains' package is back on CRAN and it is again suggested (cf. the 'param_optim' vignette). ([#175](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/175))
+- The 'hydroPSO' package is no longer suggested (but the code linked to its use and is always present in the 'param_optim' vignette). ([#182](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/182))
+
+____________________________________________________________________________________
+
+
+### 1.7.4 Release Notes (2023-04-11)
+
+#### CRAN-compatibility updates
+
+- The 'Rmalschains' package is no longer suggested (but the code linked to its use and is always present in the 'param_optim' vignette). ([#172](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/172))
+- `.ErrorCrit()` and `.FeatModels` function are no more exported. ([#173](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/173))
+
+____________________________________________________________________________________
+
+
+### 1.7.0 Release Notes (2022-02-21)
+
+#### New features
+
+- Semi-distributed modelling mode can now use the regularisation calibration proposed by [Lavenne et al. (2019)](https://doi.org/10.1029/2018WR024266). Added the `CreateInputsCrit_Lavenne()` to define a composite criterion based on the formula. Added the `CreateErrorCrit_GAPX()` function to compute an error criterion based on the GAPX formula. ([#111](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/111))
+- `OutputsModel`, returned by the `RunModel_*GR*()` function, gains a `RunOptions` element which is a list and contains 2 sub-elements: `WarmUpQsim` (vector series of simulated discharge on the warm-up period) and `Param` (vector of the model parameter values). ([#123](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/123))
+- `plot.OutputsModel()` gains a `AxisTS` argument in order to manage x-axis representing calendar dates and times. It avoids to display ugly x-axis. ([#122](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/122))
+
+
+#### Deprecated and defunct
+
+- The deprecated `LatRad` argument has been removed from the `PEdaily_Oudin()` function. ([#81](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/81))
+- The deprecated `Qobs` argument has been removed from the `CreateInputsCrit()` function. ([#81](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/81))
+- The deprecated `Ind_zeroes` argument has been removed from the `CreateInputsCrit()` function. ([#81](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/81))
+- The deprecated `verbose` argument has been removed from the `CreateInputsCrit()` function. ([#81](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/81))
+- The deprecated `FUN_CRIT` argument has been removed from the `ErrorCrit()` function. ([#81](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/81))
+
+
+#### Bug fixes
+
+- `SeriesAggreg()` now correctly reorders regime time series when the monthly regime is computed from a time series that does not start in January. It also keeps original `data.frame` column names. ([#133](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/133))
+- `DataAltiExtrapolation_Valery()` now correctly extract HypsoData values for each elevation layers. The selected indices were wrong (one less than expected) ([#144](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/144))
+- `CreateIniStates()` does not return anymore an error message when `IntStore` is set and `RunModel_GR5H` is used. ([#144](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/145))
+
+
+#### Major user-visible changes
+
+- `RunModel_Lag()` now handles warm-up period simulation (set in `CreateRunOptions()`). ([#132](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/132))
+
+
+#### Minor user-visible changes
+
+- `PE_Oudin()` can use inconsistent time series. It allows to mixing time series from different stations. ([#134](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/134))
+- Added the use of `.GetFeatModel()` in `CreateCalibOptions()` and `CreateIniStates()` functions in order to simplified their codes. ([#111](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/111))
+- Added the `.FunTransfo` in order to manage the parameter transformations and to simplified the code of the `CreateCalibOptions()` function ([#111](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/111))
+- Added the `.ArgumentsCheckGR()` function in order to check the arguments of the` RunModel_*()` functions and simplified their codes. ([#129](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/129))
+- Added the `.GetOutputsModelGR()` function in order to manage the outputs of the` RunModel_*()` functions and simplified their codes. ([#129](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/129))
+- The code of the `plot.OutputsModel()` function has been slightly simplified. ([#122](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/122), [#147](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/147))
+
+
+#### Version control and issue tracking
+
+- Added tests to check that the parameter sets returned by calibration algorithm do not change for any of the models. ([#120](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/120))
+- Added tests to detect Decreased performance of calibration execution time. ([#136](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/136))
+- Fixed the reverse package dependencies checked by the CI pipelines. ([#146](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/146))
+
+____________________________________________________________________________________
+
+
+### 1.6.12 Release Notes (2021-04-27)
+
+#### New features
+
+- `CreateInputsModel()` gains a `QupstrUnit` argument in order to manage the unit of the flow in the `Qupstream` argument in case of the use of a semi-distributed version of a hydrological model. ([#110](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/110))
+- `RunModel_Lag()` gains a `QcontribDown` argument containing the time series of the runoff contribution of the downstream sub-basin in case of the use of a semi-distributed version of a hydrological model. ([#109](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/109))
+
+
+#### Bug fixes
+
+- Fixed bug in `RunModel`. The `RunModel_Lag()` can now be passed to the `FUN_MOD` argument. ([#108](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/108))
+- Fixed bug in `RunModel_Lag()`. The function no longer returns two values for a single time step run. ([#102](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/102))
+- Fixed bug in `RunModel_Lag()`. The `StateEnd` value is now correct when there are more than a single upstream basin. ([#103](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/103))
+- Fixed bug in `RunModel_Lag()`. The `StateEnd` value is now correct when the upstream flow unit is mm/time step. ([#104](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/104))
+- Fixed bug in `RunModel_CemaNeigeGR5H()`. The solid precipitation are now taken into account in the GR5H model. ([#105](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/105))
+- Fixed bug in `RunModel_CemaNeige()` and `CreateInputsModel()`. `RunModel_CemaNeige()` now runs at the hourly time step. ([#106](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/106))
+- Fixed the 'param_optim' vignette. The starting points used for the multi-start approach are now in the transformed space.([#101](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/101))
+
+
+#### Major user-visible changes
+
+- `LengthHydro` must now be set in kilometers (not anymore in meters) in the `CreateInputsModel()` function. ([#112](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/112))
+- `TransfoParam_GR5H()` now use the same transformation as `TransfoParam_GR4H()` for the X1 parameter. The previous transformation set by Ficchì seems unnecessary as it provokes irrealistically high X1 values. ([#50](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/50))
+
+
+#### Minor user-visible changes
+
+- The `RunModel*()` functions now run faster. The computation times are significantly shorter for long times series with many time steps (e.g. hourly times series), due to a better management of the missing values in and out the Fortran codes. Only simulation computation times have been improved (it is largely invisible to the user for calibration computation times). ([#113](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/113))
+- The external calibration algorithms used in 'param_optim' and 'param_mcmc' vignettes now run faster. The `RunModel_*()` functions used during the parameter estimation process now run faster because the outputs contain only the simulated flows (see the `Outputs_Sim` argument in the `CreateRunOptions()` help page).
+- Added `.FeatModels()` and `.GetFeatModel()` functions in order to repectively store and get model features (e.g. name, number of parameters, time unit). Therefore the codes of the `CreateInputsModel()` and the `CreateRunOptions()` functions have been simplified. ([#106](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/106))
+
+
+#### Version control and issue tracking
+
+- The CI pipelines now fail when the checks return a warning message (and not just an error message). ([#86](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/86))
+- The reverse dependencies packages (e.g. the 'airGRteaching' or the 'airGRdatassim') are now checked by the CI pipelines. ([#86](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/86))
+
+____________________________________________________________________________________
+
+
+### 1.6.10.4 Release Notes (2021-01-29)
+
+#### New features
+
+- Added a section 'param_optim' vignette to explain how to manage with multiobjective optimization using the 'caRamel' package. ([#61](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/61))
+
+
+#### Major user-visible changes
+
+- `Imax()` now returns an error message when `IndPeriod_Run` doesn't select 24 hours by day, instead of `numeric(0)`. ([#92](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/92))
+
+
+#### Minor user-visible changes
+
+- Fixed warning returned by GCC Fortran when compiling `frun_GR5H.f90`. ([#93](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/93))
+
+
+#### CRAN-compatibility updates
+
+- Coerce `character` dates into `POSIXlt` in `RunModel_GR1A()` example and in `SeriesAggreg()` tests in order to avoid bad subsetting on time series due to mixing UTC and local time (error returned on macOS flavors). ([#94](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/94))
+
+____________________________________________________________________________________
+
+
+### 1.6.9.27 Release Notes (2021-01-18)
+
+#### New features
+
+- Added `SeriesAggreg` S3 method with functions for `InputsModel`, `OutputsModel`, `list`, `data.frame` class objects. This new version of the `SeriesAggreg()` function also allows to compute regimes. ([#25](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/25))
+- Added `.GetAggregConvertFun()` private function in order to choose automatically the `ConvertFun` to apply on each element of objects used in `SeriesAggreg.InputsModel()` and `SeriesAggreg.OutputsModel()`.
+- Added `.AggregConvertFunTable` data.frame that allows the user to see what names of list items or data.frame column names are guessed and eventually customise this correspondence table.
+- `PE_Oudin()` now presents a `RunFortran` argument to run the code in Fortran or in R. The Fortran mode is the fastest. ([#62](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/62))
+- Added `RunModel_Lag()` which allows to perform a single run for the Lag model over the test period in order to run semi-distributed GR models. ([#34](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/34))
+- Added the 'sd_model' vignette to explain how to manage the use of semi-distributed GR models. ([#34](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/34))
+- Added `[` S3 method for `InputsModel` class object in order to extract subsets of it. ([#67](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/67))
+
+
+#### Deprecated and defunct
+
+- The `TimeFormat` argument is now deprecated in `SeriesAggreg()`. ([#41](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/41))
+- The `NewTimeFormat` argument is now deprecated in `SeriesAggreg()` and replaced by the `Format` argument. ([#41](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/41))
+- The deprecated `RunSnowModule` argument has been removed from the `CreateRunOptions()` function. ([#23](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/23))
+
+
+#### Bug fixes
+
+- Fixed bug in`SeriesAggreg()`. The function now runs when `TimeLag >= 3600`.
+([#41](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/41))
+- Fixed bug in`SeriesAggreg()`. The function now runs when the time series contain some columns entirely filled with missing values. ([#43](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/43))
+- Fixed bug in `RunModel_GR1A()`. Reversed PotEvap and Precip outputs are now reordered (in the previous versions PotEvap contained the precipitation values and Precip contained the evapotranspiration values, the Qsim values were already correct). ([#65](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/65))
+
+
+#### Major user-visible changes
+
+- Added output to `RunModel_GR2M()` function (Ps). ([#51](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/51))
+- `PE_Oudin()` can now run for several locations (i.e. several latitudes) in the Fortran mode (`RunFortran = TRUE`). In this case `Lat` must be of the same length as `Temp`. ([#62](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/62))
+- `RunModel()` now allows to run semi-distributed GR models. ([#34](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/34))
+- The `ConvertFun` argument of the `SeriesAggreg()` function can now be set to names of aggregation functions that return value of length 1 (not only `"sum"` or `"mean"`, but e.g. `"min"`, `"max"`, `"Q95"`). ([#82](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/82))
+
+
+#### Minor user-visible changes
+
+- The `.FortranOutputs()` function is no longer exported in the namespace.
+- `RunModel_GR1A()` now uses the Fortran version of the model code. This code is no longer duplicated: the R version which was used was removed. ([#65](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/65))
+- Character argument verification now use partial matching in `PE_Oudin()` and `SeriesAggreg()` functions. ([#37](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/37))
+- `RunModel_*()` funcions were cleaned up, with no effect on their outputs. ([#14](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/14))
+- `.ErrorCrit()` function now returns a warning message when a criterion computed on less than 10 time-steps (whatever the unit of the time step). ([#14](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/14))
+- Added the diagram of GR5H in the `RunModel_GR5H()` documentation. ([#49](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/49))
+- The `Exch` was renames `AExch` in the `RunModel_GR2M()` output. ([#87](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/87))
+- Added 'Es' and 'Ps' on the GR2M diagram available in the `RunModel_GR2M()` help page. ([#88](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/88))
+- The `plot.OutputsModel()` function does not check anymore the time step by comparing the calculation of the difference of the last two time steps because it is already checked by the class of the `OutputsModel` object, which is therefore assumed to be necessarily valid. ([#56](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/56))
+
+
+#### Version control and issue tracking
+
+- Implement automatic tests in the package. ([#52](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/52))
+
+
+#### CRAN-compatibility updates
+
+- 'airGR' now depends on R >= 3.1.0 because of the use of the `anyNA` function.
+- The 'hydroPSO' package is back on CRAN and it is again suggested (cf. the 'param_optim' vignette). ([#38](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/38))
+- For more safety, the following "basic" packages are now imported : 'graphics', 'grDevices', 'stats', 'utils. ([#74](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/74))
+
+____________________________________________________________________________________
+
+
+### 1.4.3.65 Release Notes (2020-02-28)
+
+#### CRAN-compatibility updates
+
+- The run period is reduced in the example of the `Imax()` function in order to run faster.
+- The 'hydroPSO' package is no longer suggested (but the code linked to its use and is always present in the 'param_optim' vignette).
+
+____________________________________________________________________________________
+
+
+### 1.4.3.60 Release Notes (2020-01-29)
+
+#### New features
+
+- A digital object identifier (DOI) now allows to identify the manual of the 'airGR' package. When you use 'airGR' in your work, please always cite both the article and the manual. The last one allows to know the version of the package that is used in order to enhance reproducible research. The references can be displayed with the `citation("airGR")` command.
+
+#### Bug fixes
+
+- Fixed bug in `Imax()`. The default value of the `TestedValues` argument was wrong due to a mistyped argument name in the `seq()` function.
+
+____________________________________________________________________________________
+
+
+### 1.4.3.52 Release Notes (2020-01-21)
+
+
+#### New features
+
+- `plot.Outputsmodel()` now allows to draw actual evapotranspiration when `which = "ActEvap"` or `which = "All"` (overlaid to potential evapotranspiration if already drawn). ([#2](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/2))
+- Added `RunModel_GR5H()` and `RunModel_CemaNeigeGR5H()` functions to run the hourly model GR5H (with or without the CemaNeige module). These models present an optional additionnal interception store. ([#13](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/13))
+- Added `Imax()` which allows to estimate the maximum capacity of the GR5H interception store. ([#13](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/13))
+
+
+#### Bug fixes
+
+- Fixed bug in `TransfoParam_GR1A()`. The number of model parameters was wrong (2 instead of 1) which caused an error during the GR1A model calibration. ([#1](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/1))
+- Fixed bug in `plot.OutputsModel()`. The function does not return any error message when `log_scale = TRUE`, `Qobs = NULL` and user want to draw flows time series.
+- Fixed bug in `RunModel_*GR*()`. The functions do not return any error message anymore due to slightly negative values returned by GR4H, GR4J, GR5J or GR6J Fortran codes (the message was returned by `CreateIniStates()` when the final states were created). The `RunModel_*GR*()` functions now return zero instead of these slightly negative values, except for the ExpStore where negatives values are allowed.
+- Fixed bug in the `.ErrorCrit()` function. The Box-Cox transformation formula is now corrected when the `ErrorCrit*()` functions are used.
+
+
+#### Major user-visible changes
+
+- Added outputs to `RunModel_GR4H()` function (Pn, Ps, AExch1, AExch2).
+
+
+#### Minor user-visible changes
+
+- Added the diagram of GR2M in the `RunModel_GR2M()` documentation.
+- Fortran codes cleaned and translated from F77 to F90. ([#18](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/18))
+
+
+#### CRAN-compatibility updates
+
+- Cleaning of the Fortran codes (comment formatting).
+
+____________________________________________________________________________________
+
+
+### 1.3.2.42 Release Notes (2019-09-20)
+
+
+#### Version control and issue tracking
+
+- Users can now track changes (`https://gitlab.irstea.fr/HYCAR-Hydro/airgr`) and issues (`https://gitlab.irstea.fr/HYCAR-Hydro/airgr/issues`).
+
+
+#### Bug fixes
+
+- Fixed bug in `RunModel_CemaNeige()`. The function now runs correctly when `IndPeriod_WarmUp = 0L` in `CreateRunOptions()` in order to completely disable the warm-up period (e.g. to perform a forecast form a given initial state).
+- Fixed bug in `CreateIniStates()`. The function now returns the right number of end states when CemaNeige is used without hysteresis.
+- Fixed bug in the `RunModel_CemaNeige*()` functions. G and Gthr end states are no more inverted in the output values.
+
+
+#### Minor user-visible changes
+
+- Spurious flows set to `NA` into the `BasinObs` time series of the `L0123001` dataset.
+
+____________________________________________________________________________________
+
+
+### 1.3.2.23 Release Notes (2019-06-20)
+
+
+#### New features
+
+- `CreateInputsCrit()` now allows power (as a numeric or as a character) and the Box-Cox transformations in the `transfo` argument.
+- Added `RunModel_CemaNeigeGR4H()` function to run the hourly model GR4H with the CemaNeige module.
+- Added `PE_Oudin()` function to compute Oudin's potential evapotranspiration for hourly or daily time steps.
+- `plot.OutputsModel()` now presents a `LayoutMat` argument (and additionnal related argument: `LayoutWidths`, `LayoutHeights`) to specify complex plot arrangements.
+
+
+#### Deprecated and defunct
+
+- The `PEdaily_Oudin()` function is deprecated and his use has been replaced by the use of `PE_Oudin()`.
+
+
+#### Bug fixes
+
+- Fixed bug in `plot.OutputsModel()`. The function now runs correctly when the `which` argument contains the `"CorQQ"` value without `"CumFreq"`.
+
+
+#### Major user-visible changes
+
+- `plot.OutputsModel()` can now draw PE or error time series if the `which` argument is set to `"all"` or `"PotEvap"` or `"Error"`.
+- `plot.OutputsModel()` now allows new values for the which argument: `"all"` corresponds to all graphs, `"synth"` corresponds to the main graphs (default value; corresponding to `"all"` in the previous versions of the package) (i.e. `c("Precip", "Temp", "SnowPack", "Flows", "Regime", "CumFreq", "CorQQ")`), `"ts"` corresponds to the time series graphs (i.e. `c("Precip", "PotEvap", "Temp", "SnowPack", "Flows")`) and "perf" corresponds to the performance graphs (i.e. `c("Error", "Regime", "CumFreq", "CorQQ")`).
+
+
+#### Minor user-visible changes
+
+- `.ErrorCrit()` private function added to check inputs into `ErrorCrit_*()` functions. The `ErrorCrit_*()` functions were simplified accordingly.
+- `CreateInputsCrit()` now returns `FUN_CRIT` as a character string.
+- An example is addeed to illustred the use of the `plot.OutputsModel()` function.
+
+____________________________________________________________________________________
+
+
+### 1.2.13.16 Release Notes (2019-04-03)
+
+
+#### New features
+
+- `CreateInputsCrit()` now presents a `VarObs` argument in order to allow to prepare an `InputsCrit` object in order to run a criterion on other variables than observed discharges with the `ErrorCrit()` function (at the moment SCA and SWE).
+- `CreateInputsCrit()` can now prepare an `InputsCrit` object in order to compute a single criterion (`Single` class), multiple criteria (`Multi` class) or a composite criterion (`Compo` class) with the `ErrorCrit()` function.
+- `CreateInputsCrit()` now presents a `Weights` argument in order to allow to prepare an `InputsCrit` object in order to compute a composite criterion (`Compo` class) with `ErrorCrit()` or `Calibration_Michel()`.
+- `CreateInputsCrit()` now returns a `idLayer` element to indicate which layer(s) to use for SCA or SWE aggregation.
+- `CreateInputsCrit()` now presents a `warnings` argument to replace the verbose action (the `verbose` argument is kept to print messages).
+- In `CreateInputsCrit()`, it is now possible to set the following arguments as atomic (as before) or as list: `FUN_CRIT`, `VarObs`, `Obs`, `BoolCrit`, `transfo`, `Weights`. If the list format is chosen, all the lists must have the same length.
+- `CreateRunOptions()`, `CreateIniStates()` and `CreateCalibOptions()` now present a `IsHyst` argument to give the possibility to use the Linear Hysteresis with CemaNeige. The objects returned present an `hysteresis` class.
+- `CreateRunOptions()` now presents a `warnings` argument to replace the verbose action (the `verbose` argument is kept to print messages).
+- Added `TransfoParam_CemaNeigeHyst()` function in order to take into account transformation of the parameters of the CemaNeige module when the Linear Hysteresis is used.
+- Added the `X0310010` dataset to run the examples using the Linear Hysteresis with CemaNeige (it contains necessary SCA data).
+- Added the 'cemaneige_hysteresis' vignette to explain how to manage the use of the Linear Hysteresis with CemaNeige.
+
+
+#### Deprecated and defunct
+
+- The `Qobs` argument is now deprecated in `CreateInputsCrit()` and has been renamed `Obs`.
+- The `FUN_CRIT` argument is now deprecated in `ErrorCrit()`. This function now gets this information from the `InputsCrit` argument.
+- The `FUN_CRIT` argument is now deprecated in `Calibration_Michel()`. This function now gets this information from the `InputsCrit` argument.
+- The `plot_OutputsModel()` had been deprecated in 'airGR' 1.0.4 (it had been replaced by the use of `plot.OutputsModel()` or `plot()`) and is defunct now.
+
+
+#### Major user-visible changes
+
+- `CreateInputsCrit()` now return a list of `InputsCrit` (each element is of the `Single` class) in the cases of multiple or a composite criteria.
+- `ErrorCrit_*()` functions now return an error message if the `InputsCrit` object is of class `Multi` or `Compo`.
+- `ErrorCrit()` function can now run on a multiple or a composite `InputsCrit`. In these cases, it returns a list of `ErrorCrit`.
+- `ErrorCrit()` and `ErrorCrit_*()` functions can now assess Q, SCA or SWE simulations.
+- `Calibration_Michel()` function can now run on a composite `InputsCrit`. It returns a composite value of error and the formula used to calculate it.
+- Model diagrams added in documentations of `RunModel_GR4J()`, `RunModel_GR5J()` and `RunModel_GR6J()` functions.
+- It is now possible to be redirected to the `plot.OutputsModel()` documentation with `?plot`.
+- It is now possible to use a character vector for all `FUN_*` arguments (in addition to function objects) in the following functions: `Calibration()`, `Calibration_Michel()`, `CreateCalibOptions()`, `CreateIniStates()`, `CreateIniStates()`, `CreateInputsCrit()`, `CreateInputsModel()`, `CreateRunOptions()`, `ErrorCrit()`, `RunModel()` and `TransfoParam()`.
+
+
+#### Minor user-visible changes
+
+- `ErrorCrit_*()` functions now return objects of class `ErrorCrit` and `NSE`, `KGE`, `KGE2` or `RMSE`.
+- `.FortranOutputs()` private function added to manage Fortran outputs.
+- Outputs of `frun_GR2M` Fortran subroutine were reordered.
+- `DataAltiExtrapolation_Valery()` now returns named elements of lists relative to elevation layer.
+- `Calibration()` function now returns an error message if `FUN_CALIB` is not a function.
+- Inputs of `PEdaily_Oudin()` are now checked.
+- `PEdaily_Oudin()` example corrected (the Julian day was one day too early).
+- `plot.OutputsModel()` does not return a warning message anymore when `Qobs = NULL`.
+- Inputs of `TransfoParam*()` functions are now checked.
+- The order of authors has been updated in the DESCRIPTION and the CITATION files.
+
+
+#### CRAN-compatibility updates
+
+- Tabulations removed, unused variables removed and variable statements fixed in Fortran files.
+
+____________________________________________________________________________________
+
+
+### 1.0.15.2 Release Notes (2018-10-10)
+
+
+#### Bug fixes
+
+- Fixed bug in `CreateRunOptions()`. The function now accounts correctly for leap years when no warm-up period is defined.
+
+
+#### Minor user-visible changes
+
+- `CreateRunOptions()` was cleaned up, with no effect on its outputs.
+
+
+#### CRAN-compatibility updates
+
+- The `vignetteParam*.rda` datasets moved to the inst directory. It contains different objects needed for 'param_optim' and 'param_mcmc' vignettes.
+
+____________________________________________________________________________________
+
+
+### 1.0.14.1 Release Notes (2018-09-28)
+
+
+#### New features
+
+- `PEdaily_Oudin()` now presents a `LatUnit` argument which allows to choose the unit of the latitude (radians and degrees).
+
+
+#### Deprecated and defunct
+
+- The `LatRad` argument is now deprecated in `PEdaily_Oudin()` and replaced by the `Lat` argument.
+- The unused `Ind_zeroes` argument of the `CreateInputsCrit()` function is now deprecated.
+- The `verbose` argument is now deprecated in `CreateInputsCrit()` and replaced by the `warnings` argument.
+
+
+#### Major user-visible changes
+
+- `Calibration_Michel()` is now faster during the grid-screening step when a parameter is set using `FixedParam` in `CreateCalibOptions()`.
+- `CreateCalibOptions()` now returns an error when all the parameters are set in the `FixedParam` argument and a warning message when all the parameters are free (NA) in the `FixedParam` argument.
+- `CreateInputsCrit()` now returns an error when `epsilon` is not positive.
+- `CreateInputsCrit()` now returns a warning message in the following case: there are zeroes values in `Qobs`, `epsilon = NULL` and `transfo = log` or `inv`.
+- `ErrorCrit_*()` functions now return a warning message in the following case: there are zeroes values in `Qobs` or `Qsim`, `epsilon = NULL` and `transfo = log` or `inv`.
+
+
+#### Minor user-visible changes
+
+- Several functions of the package were cleaned up or slightly modified, with no effect on their outputs.
+- Dubious Qls and Qmm values set to NA values between 1997-01-05 and 1997-01-21 in the `L0123001` dataset.
+- ORCID numbers are now joined to the names of the authors of the package.
+
+
+#### CRAN-compatibility updates
+
+- Function name changed in a vignettes to avoid error during the check on Debian distribution
+- As recomanded by CRAN managers, the NEWS file is now at the text format and is no more just a link to the 'airGR' Website
+- Added the `Vignette_Param.` datasets in order to reduce runtime during the re-building of vignettes. It contains different objects needed for param_optim and param_mcmc vignettes.
+
+____________________________________________________________________________________
+
+
+### 1.0.10.11 Release Notes (2018-06-29)
+
+
+#### Bug fixes
+
+- Fixed bug in `RunModel_GR2M()`. The function now returns the total precipitation (P) instead of the net rainfall (P1).
+
+
+#### Major user-visible changes
+
+- `RunModel_GR2M()` now returns more explicit precipitation outputs names.
+- `CreateInputsCrit()` now returns a warning message when the KGE (or KGE') is used with a log transformation on flows.
+- The article reference is corrected.
+
+
+#### Minor user-visible changes
+
+- The documentation and help of several functions were improved.
+
+____________________________________________________________________________________
+
+
+### 1.0.9.64 Release Notes (2017-11-10)
+
+
+#### New features
+
+- An article describing the 'airGR' package has been published. Its reference has been added and will be displayed with `citation("airGR")`.
+- Added `CreateIniStates()` function in order to help user to format the `IniStates` argument for `CreateRunOptions()`.
+- Added the `Param_Sets_GR4J` dataset. It contains generalist parameter sets for the GR4J model.
+- Three vignettes have been added. They are relative to different calibration methods (including the generalist parameters sets of the GR4J model).
+
+
+#### Deprecated and defunct
+
+- The `RunSnowModule` argument is now deprecated in `CreateRunOptions()`.
+
+
+#### Bug fixes
+
+- Fixed bug in `RunModel_GR4H()`: in `frun_GR4H` Fortran subroutine, `St(2)` is now set to 0 (instead of `St(1)`) when `St(2) < 0`.
+- Fixed bug in `plot.OutputsModel()` for the regime plot when the period is less than 1 year.
+- Fixed bug in `plot.OutputsModel()` when there is no common data to plot the cumulative frequency or the correlation QQ.
+- Fixed bug in `plot.OutputsModel()` for the y-axis labelling of flows time series when `log_scale = TRUE` and `BasinArea` is used.
+
+
+#### Major user-visible changes
+
+- `RunModel_GR4J()`, `RunModel_GR5J()` and `RunModel_GR6J()` (and `CemaNeige_GR*J()`) now return Ps, Pn and actual exchanges. See the model Fortran codes for more details about the calculation of these variables.
+- `CreateInputsModel()` now returns an error when `DatesR` contains duplicated values.
+- `RunModel_GR5J` now returns `StateEnd` in the same order as the other models.
+
+
+#### Minor user-visible changes
+
+- `plot.OutputsModel()` now returns a warning message when the length of Qobs is different from the length of Qsim.
+- The X1 parameter from GR4H, GR4J, GR2M, GR5J and GR6J, the X3 parameter from GR4H, GR4J, GR5J and GR6J and the X6 parameter from GR6J are now set to 1e-2 when they are fixed to lower values. `RunModel_*()` functions now return a warning message in this case. `RunModel_*()` functions now return a warning when X4 < 0.5 and its value is set to 0.5.
+- The commands `?L0123001`, `?L0123002` and `?L0123003` now return the documentation page related to `BasinObs`.
+- Many functions of the package were cleaned up or slightly modified, with no effect on their outputs.
+- The documentation and help of several functions were improved.
+
+
+#### CRAN-compatibility updates
+
+- "airGR.c" file registers native routines.
+
+____________________________________________________________________________________
+
+
+### 1.0.5.12 Release Notes (2017-01-23)
+
+
+#### New features
+
+- `DataAltiExtrapolation_Valery()` and `CreateInputsModel()` now present a `PrecipScale` argument which allows rescaling precipitation when it is interpolated on the elevation layers when CemaNeige is used.
+
+
+#### Bug fixes
+
+- Fixed bug in `DataAltiExtrapolation_Valery()`. The elevation gradients for air temperature returned by `CreateInputsModel()` are improved.
+
+
+#### User-visible changes
+
+- `DataAltiExtrapolation_Valery()` has been improved. `DataAltiExtrapolation_Valery()` now runs faster (and by consequence `CreateInputsModel()` too, when CemaNeige is used).
+
+____________________________________________________________________________________
+
+
+### 1.0.4 Release Notes (2017-01-18)
+
+
+#### New features
+
+- `RunModel_CemaNeige()`, `RunModel_CemaNeigeGR4J()`, `RunModel_CemaNeigeGR5J()` and `RunModel_CemaNeigeGR6J()` now return air temperature for each elevation layer.
+
+
+#### Deprecated and defunct
+
+- S3 plot method defined for `OutputsModel` objects. It means that the `plot_OutputsModel()` function is deprecated and his use has been replaced by the use of `plot.OutputsModel()` or `plot()`.
+- In `plot.OutputsModel()` the `PlotChoice` argument is deprecated and has been renamed `which`.
+
+
+#### User-visible changes
+
+- `plot.OutputsModel()` displays air temperature time series for each layer when `CemaNeige` is used (argument `which = "Temp"` or `"all"`).
+
+____________________________________________________________________________________
+
+
+### 1.0.3 Release Notes (2016-12-09)
+
+
+#### New features
+
+- `ErrorCrit_*()` functions gain a `warnings` argument to replace the verbose action and the `verbose` argument now prints the criterion value(s).
+
+
+#### Bug fixes
+
+- Fixed bug in `CreateCalibOptions()` when `StartParamList` or `StartParamDistrib` arguments are used.
+
+
+#### User-visible changes
+
+- `CreateInputsModel()` now returns an error if `NLayers <= 0` when `CemaNeige` is used.
+- `plot_OutputsModel()` now displays raw values on the y-axis when the discharge time series is represented with log scale (formerly, log values of discharges were displayed on the y-axis).
+
+____________________________________________________________________________________
+
+
+### 1.0.2 Release Notes (2016-11-03)
+
+
+#### New features
+
+- `SeriesAggreg()` gains a `TimeLag` argument that corresponds to a numeric value indicating a time lag (in seconds) for the time series aggregation (useful to aggregate hourly time series to the daily time step for instance).
+ In addition, the function now accepts input dates in both `POSIXt` formats (`POSIXct` and `POSIXlt`). The output is in `POSIXct` format.
+- `plot_OutputsModel()` gains a `log_scale` argument in order to plot the flow with a log scale.
+- A tutorial is available online on the following link: https://hydrogr.github.io/airGR/.
+ It can also be displayed with the `vignette("airGR")` command.
+
+
+#### Deprecated and defunct
+
+- `CreateCalibOptions()` loses the `OptimParam` argument that was redundant with the `FixedParam` argument. The `Calibration_Michel()` was modified to take into account this change by using directly `FixedParam`, but this is transparent to the user.
+- `CreateCalibOptions()` loses the `StartParam` argument that was not used.
+
+
+#### Bug fixes
+
+- The value `sort` for the `transfo` argument of `CreateInputsCrit()` was not taken into account. It is now fixed.
+
+
+#### Major user-visible changes
+
+- The `RunModel_GR6J()` and `RunModel_CemaNeigeGR6J()` models were modified back to versions previous to 1.0.1 to prevent from unwanted efficiency criteria deterioration related to the calibration with `Calibration_Michel()`.
+ The actual model codes were not modified but the `TransfoParam_GR6J()` and `CreateCalibOptions()` functions were modified regarding the X5 parameter.
+ It is strongly advised to use airGR 1.0.2 for the `RunModel_GR6J()` and `RunModel_CemaNeigeGR6J()` functions if you are using `Calibration_Michel()`, as they are much more efficient.
+ In case you were using your own calibration algorithm, you will not notice any difference.
+
+
+#### Minor user-visible changes
+
+- `CreateInputsModel()` and `DataAltiExtrapolation_Valery()` functions now allow both `POSIXt` formats (`POSIXct` and `POSIXlt`).
+
+____________________________________________________________________________________
+
+
+### 1.0.1 Release Notes (2016-04-21)
+
+
+#### Deprecated and defunct
+
+- The `Calibration_HBAN()` and `DataAltiExtrapolation_HBAN()` functions have respectively been renamed as `Calibration_Michel()` and `DataAltiExtrapolation_Valery()` after the names of their creators.
+- The `Calibration_optim()` function has been removed from the package.
+- The silent mode is now defined by the `verbose = TRUE` argument (formerly `quiet = FALSE`) in the following functions:
+`Calibration()`, `Calibration_Michel()`, `CreateInputsModel()`, `CreateRunOptions()`, `DataAltiExtrapolation_Valery()`, `ErrorCrit()`, `ErrorCrit_KGE()`, `ErrorCrit_KGE2()`, `ErrorCrit_NSE()`, `ErrorCrit_RMSE()`, `plot_OutputsModel()`, `SeriesAggreg()`.
+
+
+#### Major user-visible changes
+- The GR5J model has been modified: previously, two unit hydrographs were used, now only one is remaining.
+ As a consequence, simulations from the GR5J (`RunModel_GR5J()` function) and CemaNeige (`RunModel_CemaNeigeGR5J()` function) models will be different.
+- An important proportion of the transformations of the parameters have been modified (`TransfoParam_*()` functions). Since this modifies the local search, calibration results will be different .
+- The quantiles of the parameters have been recalculated with the new transformations (`CreateCalibOptions()` function). Since these quantiles constitute the starting point of the calibration algorithm, calibration results will be different.
+
+
+#### Minor user-visible changes
+
+- The Fortran model core codes have been modified:
+	- optimisation of the codes for fastening of computation;
+	- simplification of the internal variables for easier reading and understanding.
+- The list of the contributors and authors is now full.
+- The references of the package has been updated; they are returned by the following R-command `citation("airGR")`.
+
+____________________________________________________________________________________
+
+
+### 0.8.1.2 Release Notes (2015-08-21)
+
+
+#### Bug fixes
+
+- Fixed bug in `CreateInputsModel()` that was related to the handling of missing values.
+- Fixed bug in `CreateRunOptions()` that prevented the correct use of the `IniResLevels` argument (to manually set the filling rate of the production and routing stores).
+
+
+#### Minor user-visible changes
+
+- Removal of an unnecessary warning when `IndPeriod_WarmUp = 0`.
+
+
+#### CRAN-compatibility updates
+
+- Modification of namespace file to ensure proper use under linux without compilation issues.
+
+____________________________________________________________________________________
+
+
+### 0.8.0.2 Release Notes (2015-04-15)
+
+
+#### New features
+
+- Three new hydrological models: `RunModel_GR4H() function for ` GR4H (hourly), `RunModel_GR2M()` function for GR2M (monthly) and `RunModel_GR1A()` function for GR1A (yearly).
+- New function `SeriesAggreg()` to easily aggreg timesteps.
+
+
+#### Bug fixes
+
+- Fixed bug in `ErrorCrit_RMSE()` which led to incorrect calibration (the criterion was maximised instead of minimised).
+
+
+#### Major user-visible changes
+
+- Update of the functions `CreateRunOptions()`, `CreateCalibOptions()` and `plot_OutputsModel()` to handle the new models.
+- Modification of CemaNeige Fortran code to add an update of Gratio after the SnowPack update (no impact on snow simulation).
+
+
+#### Minor user-visible changes
+
+- Improvement of the `plot_OutputsModel()` function to allow a selection among available plots.
+- Minor update in `ErrorCrit_KGE()` and `ErrorCrit_KGE2()` to handle case when only one values in not NA.
+- Update of the scripts in airGR-advanced-example to match the structures of the `BasinData` objects.
+- Correction of formatting issue in airGR-advanced-example regarding the "List_HypsoData.txt" file.
+
+____________________________________________________________________________________
+
+
+### 0.7.4 Release Notes (2014-11-01)
+
+
+#### New features
+
+- New argument in many functions (`quiet = TRUE` or `FALSE`) to choose if the warnings should be suppressed or not.
+
+
+#### Deprecated and defunct
+
+- The `CalibrationAlgo_*()` functions were renamed into `Calibration_*()`.
+
+
+#### Bug fixes
+
+- Fixed bug in `CreateCalibOptions()` to handle models with only one parameter.
+- Fixed bug in `Calibration_HBAN()`. The function was not working properly with models having only one parameter.
+
+
+#### Major user-visible changes
+
+- CemaNeige users must now specify one `MeanAnSolidPrecip` for each elevation layer. The `CreateRunOptions()` function is impacted.
+- CemaNeige users can now specify the mean elevation of the input series (before it was always considered equal to the catchment median elevation).
+ The impacted functions are `CreateInputsModel()` and `DataAltiExtrapolation_HBAN()`.
+- New architecture with better format verification procedure (using classes) and simpler setting of default configuration.
+- New architecture where the model, calibration and error functions are in the arguments of the functions
+ (the exotic use of "generic function" created by the users has been removed).
+- Improved documentation and examples.
+
+
+#### Minor user-visible changes
+
+- Improvements allowing the arrival of new models.
+- Improvements of the argument verifications in `CreateInputsModel()`, `CreateRunOptions()`, `CreateInputsCrit()`, `CreateCalibOptions()`.
+- Improvements of all the `ErrorCrit()` functions to better account for the cases with constant flow values or local zeros.
+- Improvement of the `plot_OutputsModel` function (to handle 0 in Qobs and Qsim).
+- Improved documentation.
+
+____________________________________________________________________________________
+
+
+### 0.6.2 Release Notes (2014-02-12)
+
+
+#### New features
+
+- Additional functions for results plotting (the 'zoo' package is required for some of them).
+- Add multi-objective calibration using `nsga2()` (the 'mco' package is required).
+- The field Multiplier has been added in the ErrorCrit() outputs, to indicate whether the criterion is an error (to minimise) or and efficiency (to maximise).
+ This allows to provide real efficiency values in the outputs e.g. NSE[Q] instead of (-1) &times; NSE[Q].
+
+
+#### Deprecated and defunct
+
+- `EfficiencyCrit()` have been replaced by `ErrorCrit()` to avoid misunderstanding (by default, the algorithms minimise the error criterion).
+
+
+#### Bug fixes
+
+- RC11 bug correction: the automatic selection of the warm-up period was not working properly when no data was available from warm-up (i.e. when the user had set the run to start at the very first index).
+- RC10 bug correction: the `CalibrationAlgo_HBAN()` function was not working in the very rare case when the diagonal search was activated and lead to a set outside the authorised range.
+- RC9 bug correction: the `CalibrationAlgo_HBAN()` function was not working properly with models having only one parameter.
+- RC8 bug correction of the `ModelDefaultIniOptions()` function (this bug was introduced in the RC7 and caused an error when `IndPeriod_WarmUp = NULL`.
+- RC7 bug correction of the `ModelDefaultIniOptions()` function (the automatic selection of one year for warm-up was not handling properly missing data).
+- RC6 correction of the help files (the description of CemaNeige parameters were inverted).
+- RC5 differs from previous releases in the way the data are read and stored (in a list instead of individual vectors).
+ The package is similar, only the examples of Main and the files in MyScriptBlocks have changed.
+ All basin data are now stored inside a list named `BasinData`. This will greatly ease the future use of Rdata files (instead of txt files) as storage format for the time series of observation.
+
+
+#### Major user-visible changes
+
+- The definition of the generic function is now made in a much simpler way (e.g. see `DefineFunctions_Model()` or `DefineFunctions_ErrorCrit()`).
+
+
+#### Minor user-visible changes
+
+- Code improvements to reduce the computation time.
+- Clearer instructions for the adding and modification of a model.
+- Improvements of the documentation.
+
+
+____________________________________________________________________________________
+
+
+### 0.5.2 Release Notes (2014-02-05)
+
+
+#### Deprecated and defunct
+
+- The `SelectPer` arguments are replaced by `IndPeriod` to ease understanding.
+- The `PE` arguments are replaced by `PotEvap()` to ease understanding.
+- The `Fsol` arguments are replaced by `FracSolidPrecip` to ease understanding.
+
+
+#### Major user-visible changes
+
+- R <= 2.15 in not supported by default.
+- The check that `SelectPer_Run()` is continuous is now made in the `CheckArg()` functions.
+- Check of the model functioning time step.
+- Name of the calibration criterion provided in `OutputsAlgo()`.
+
+
+#### Minor user-visible changes
+
+- Missing values in Fortran are now -999.999 instead of -9.999.
+
+
+____________________________________________________________________________________
+
+
+### 0.5.1 Release Notes (2014-01-27)
+
+
+#### New features
+
+- New `EfficiencyCrit_NSE_sqrtQ()` function to compute NSE criterion on sqrt flows.
+
+
+#### Bug fixes
+
+- Incorrect arguments in the call to `RunModelAndCrit` from `CalibrationAlgo_optim_stats` and `CalibrationAlgo_nlminb_stats`.
+- `CalibrationAlgo_nlminb_stats` argument was wrongly defined in `DefineFunctions_CalibrationAlgo()` (`optim` instead of `nlminb`).
+- Format checking for `RunOptions` was incorrectly made in `CheckArg()` function.

R/BasinData.R

-#' @name BasinInfo
-#' @docType data
-#' @title Data sample: characteristics of a fictional catchment (L0123001, L0123002 or L0123003)
-#' @description
-#' R-object containing the code, station's name, area and hypsometric curve of the catchment.
-#' @encoding UTF-8
-#' @format
-#' List named 'BasinInfo' containing
-#' \itemize{
-#' \item two strings: catchment's code and station's name
-#' \item one float: catchment's area in km2
-#' \item one numeric vector: catchment's hypsometric curve (min, quantiles 01 to 99 and max) in metres
-#' }
-#' @examples 
-#'    require(airGR)
-#'    data(L0123001)
-#'    str(BasinInfo)
-
-NULL
-
-
-#' @name BasinObs
-#' @docType data
-#' @title Data sample: time series of observations of a fictional catchment (L0123001, L0123002 or L0123003)
-#' @description
-#' R-object containing the times series of precipitation, temperature, potential evapotranspiration and discharges. \cr
-#' Times series for L0123001 or L0123002 are at the daily time-step for use with daily models such as GR4J, GR5J, GR6J, CemaNeigeGR4J, CemaNeigeGR5J and CemaNeigeGR6J.
-#' Times series for L0123003 are at the hourly time-step for use with hourly models such as GR4H.
-#' @encoding UTF-8
-#' @format
-#' Data frame named 'BasinObs' containing
-#' \itemize{
-#' \item one POSIXlt vector: time series dates in the POSIXlt format
-#' \item five numeric vectors: time series of catchment average precipitation [mm], catchment average air temperature [degC], catchment average potential evapotranspiration [mm], outlet discharge [l/s], outlet discharge [mm]
-#' }
-#' @examples 
-#'    require(airGR)
-#'    data(L0123001)
-#'    str(BasinObs)
-
-NULL
-

R/Calibration.R

-#*************************************************************************************************
-#' Calibration algorithm which minimises the error criterion using the provided functions. \cr
-#*************************************************************************************************
-#' @title  Calibration algorithm which minimises an error criterion on the model outputs using the provided functions
-#' @author Laurent Coron (June 2014)
-#' @seealso \code{\link{Calibration_Michel}}, \code{\link{Calibration_optim}},
-#'          \code{\link{RunModel}}, \code{\link{ErrorCrit}}, \code{\link{TransfoParam}},
-#'          \code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}}, 
-#'          \code{\link{CreateInputsCrit}}, \code{\link{CreateCalibOptions}}.
-#' @example tests/example_Calibration.R
-#' @export
-#' @encoding UTF-8
-#_FunctionInputs__________________________________________________________________________________
-#' @param  InputsModel         [object of class \emph{InputsModel}] see \code{\link{CreateInputsModel}} for details
-#' @param  RunOptions          [object of class \emph{RunOptions}] see \code{\link{CreateRunOptions}} for details
-#' @param  InputsCrit          [object of class \emph{InputsCrit}] see \code{\link{CreateInputsCrit}} for details
-#' @param  CalibOptions        [object of class \emph{CalibOptions}] see \code{\link{CreateCalibOptions}} for details
-#' @param  FUN_MOD             [function] hydrological model function (e.g. RunModel_GR4J, RunModel_CemaNeigeGR4J)
-#' @param  FUN_CRIT            [function] error criterion function (e.g. ErrorCrit_RMSE, ErrorCrit_NSE)
-#' @param  FUN_CALIB           (optional) [function] calibration algorithm function (e.g. Calibration_Michel, Calibration_optim), default=Calibration_Michel
-#' @param  FUN_TRANSFO         (optional) [function] model parameters transformation function, if the FUN_MOD used is native in the package FUN_TRANSFO is automatically defined
-#' @param  quiet               (optional) [boolean] boolean indicating if the function is run in quiet mode or not, default=FALSE
-#_FunctionOutputs_________________________________________________________________________________
-#' @return  [list] see \code{\link{Calibration_Michel}} or \code{\link{Calibration_optim}} 
-#**************************************************************************************************
-Calibration <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_CALIB=Calibration_Michel,FUN_TRANSFO=NULL,quiet=FALSE){
-    return( FUN_CALIB(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO,quiet=quiet) )
+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, ...))
 }
 

R/Calibration_Michel.R

-#*************************************************************************************************
-#' Calibration algorithm which minimises the error criterion. \cr
-#' \cr
-#' The algorithm is based on a local search procedure. 
-#' First, a screening is performed using either a rough predefined grid or a list of parameter sets
-#' and then a simple steepest descent local search algorithm is performed.
-#' 
-#' A screening is first performed either from a rough predefined grid (considering various initial 
-#' values for each paramete) or from a list of initial parameter sets. \cr
-#' The best set identified in this screening is then used as a starting point for the steepest 
-#' descent local search algorithm. \cr
-#' For this search, the parameters are used in a transformed version, to obtain uniform           
-#' variation ranges (and thus a similar pace), while the true ranges might be quite different. \cr
-#' At each iteration, we start from a parameter set of NParam values (NParam being the number of 
-#' free parameters of the chosen hydrological model) and we determine the 2*NParam-1 new candidates 
-#' by changing one by one the different parameters (+/- pace). \cr
-#' All these candidates are tested and the best one kept to be the starting point for the next    
-#' iteration. At the end of each iteration, the pace is either increased or decreased to adapt    
-#' the progression speed. A diagonal progress can occasionally be done.                        \cr
-#' The calibration algorithm stops when the pace becomes too small.                            \cr
-#'
-#' To optimise the exploration of the parameter space, transformation functions are used to convert
-#' the model parameters. This is done using the TransfoParam functions.
-#*************************************************************************************************
-#' @title  Calibration algorithm which minimises the error criterion using the Irstea-HBAN procedure
-#' @author Laurent Coron (August 2013)
-#' @references
-#'   Michel, C. (1991),
-#'        Hydrologie appliquée aux petits bassins ruraux, Hydrology handout (in French), Cemagref, Antony, France.
-#' @example tests/example_Calibration_Michel.R
-#' @seealso \code{\link{Calibration}}, \code{\link{Calibration_optim}},
-#'          \code{\link{RunModel_GR4J}}, \code{\link{TransfoParam_GR4J}}, \code{\link{ErrorCrit_RMSE}},
-#'          \code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}}, 
-#'          \code{\link{CreateInputsCrit}}, \code{\link{CreateCalibOptions}}.
-#' @encoding UTF-8
-#' @export
-#_FunctionInputs__________________________________________________________________________________
-#' @param  InputsModel         [object of class \emph{InputsModel}] see \code{\link{CreateInputsModel}} for details
-#' @param  RunOptions          [object of class \emph{RunOptions}] see \code{\link{CreateRunOptions}} for details
-#' @param  InputsCrit          [object of class \emph{InputsCrit}] see \code{\link{CreateInputsCrit}} for details
-#' @param  CalibOptions        [object of class \emph{CalibOptions}] see \code{\link{CreateCalibOptions}} for details
-#' @param  FUN_MOD             [function] hydrological model function (e.g. RunModel_GR4J, RunModel_CemaNeigeGR4J)
-#' @param  FUN_CRIT            [function] error criterion function (e.g. ErrorCrit_RMSE, ErrorCrit_NSE)
-#' @param  FUN_TRANSFO         (optional) [function] model parameters transformation function, if the FUN_MOD used is native in the package FUN_TRANSFO is automatically defined
-#' @param  quiet               (optional) [boolean] boolean indicating if the function is run in quiet mode or not, default=FALSE
-#_FunctionOutputs_________________________________________________________________________________
-#' @return  [list] list containing the function outputs organised as follows:
-#'          \tabular{ll}{
-#'          \emph{$ParamFinalR  }  \tab   [numeric] parameter set obtained at the end of the calibration \cr
-#'          \emph{$CritFinal    }  \tab   [numeric] error criterion obtained at the end of the calibration \cr
-#'          \emph{$NIter        }  \tab   [numeric] number of iterations during the calibration \cr
-#'          \emph{$NRuns        }  \tab   [numeric] number of model runs done during the calibration \cr
-#'          \emph{$HistParamR   }  \tab   [numeric] table showing the progression steps in the search for optimal set: parameter values \cr
-#'          \emph{$HistCrit     }  \tab   [numeric] table showing the progression steps in the search for optimal set: criterion values \cr
-#'          \emph{$MatBoolCrit  }  \tab   [boolean] table giving the requested and actual time steps when the model is calibrated \cr
-#'          \emph{$CritName     }  \tab   [character] name of the calibration criterion \cr
-#'          \emph{$CritBestValue}  \tab   [numeric] theoretical best criterion value \cr
-#'          }
-#**************************************************************************************************
-Calibration_Michel <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO=NULL,quiet=FALSE){
-
-
-##_____Arguments_check_____________________________________________________________________
-    if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
-    if(inherits(RunOptions,"RunOptions")==FALSE){ stop("RunOptions must be of class 'RunOptions' \n"); return(NULL); }  
-    if(inherits(InputsCrit,"InputsCrit")==FALSE){ stop("InputsCrit must be of class 'InputsCrit' \n"); return(NULL); }  
-    if(inherits(CalibOptions,"CalibOptions")==FALSE){ stop("CalibOptions must be of class 'CalibOptions' \n"); return(NULL); }  
-    if(inherits(CalibOptions,"HBAN")==FALSE){ stop("CalibOptions must be of class 'HBAN' if Calibration_Michel is used \n"); return(NULL); }  
-
-
-   ##_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    )){ FUN_TRANSFO <- TransfoParam_CemaNeige; }
-      if(identical(FUN_MOD,RunModel_CemaNeigeGR4J) | identical(FUN_MOD,RunModel_CemaNeigeGR5J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)){
-        if(identical(FUN_MOD,RunModel_CemaNeigeGR4J)){ FUN1 <- TransfoParam_GR4J; FUN2 <- TransfoParam_CemaNeige; }
-        if(identical(FUN_MOD,RunModel_CemaNeigeGR5J)){ FUN1 <- TransfoParam_GR5J; FUN2 <- TransfoParam_CemaNeige; }
-        if(identical(FUN_MOD,RunModel_CemaNeigeGR6J)){ FUN1 <- TransfoParam_GR6J; FUN2 <- TransfoParam_CemaNeige; }
-        FUN_TRANSFO <- function(ParamIn,Direction){
-          Bool <- is.matrix(ParamIn);
-          if(Bool==FALSE){ 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==FALSE){ ParamOut <- ParamOut[1,]; }
-          return(ParamOut);
+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
+  ## use unique() to avoid duplicated values when a parameter is set
+  ProposeCandidatesGrid <- function(DistribParam) {
+    expand.grid(lapply(seq_len(ncol(DistribParam)), function(x) unique(DistribParam[, x])))
+  }
+
+
+  ##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)
+  }
+  ##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)
         }
       }
-      if(is.null(FUN_TRANSFO)){ stop("FUN_TRANSFO was not found (in Calibration function) \n"); return(NULL);  }
     }
-
-    ##_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 \n"); return(NULL);  }
-    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){
-      ##Managing_matrix_sizes
-        Nvalmax <- nrow(DistribParam);
-        NParam <- ncol(DistribParam);
-        ##we_add_columns_to_MatDistrib_until_it_has_20_columns
-        DistribParam2 <- matrix(NA,nrow=Nvalmax,ncol=20);
-        DistribParam2[1:Nvalmax,1:NParam] <- DistribParam;
-        ##we_check_the_number_of_values_to_test_for_each_param
-        NbDistrib <- rep(1,20);
-        for(iC in 1:20){ NbDistrib[iC] <- max( 1 , Nvalmax-sum(is.na(DistribParam2[,iC])) ); }
-      ##Loop_on_the_various_values_to_test ###(if 4 param and 3 values for each => 3^4 sets)
-      ##NB_we_always_do_20_loops ###which_is_here_the_max_number_of_param_that_can_be_optimised
-        VECT <- NULL;
-        for(iL01 in 1:NbDistrib[01]){ for(iL02 in 1:NbDistrib[02]){ for(iL03 in 1:NbDistrib[03]){ for(iL04 in 1:NbDistrib[04]){ for(iL05 in 1:NbDistrib[05]){ 
-        for(iL06 in 1:NbDistrib[06]){ for(iL07 in 1:NbDistrib[07]){ for(iL08 in 1:NbDistrib[08]){ for(iL09 in 1:NbDistrib[09]){ for(iL10 in 1:NbDistrib[10]){
-        for(iL11 in 1:NbDistrib[11]){ for(iL12 in 1:NbDistrib[12]){ for(iL13 in 1:NbDistrib[13]){ for(iL14 in 1:NbDistrib[14]){ for(iL15 in 1:NbDistrib[15]){ 
-        for(iL16 in 1:NbDistrib[16]){ for(iL17 in 1:NbDistrib[17]){ for(iL18 in 1:NbDistrib[18]){ for(iL19 in 1:NbDistrib[19]){ for(iL20 in 1:NbDistrib[20]){
-          VECT <- c(VECT,                         
-            DistribParam2[iL01,01],DistribParam2[iL02,02],DistribParam2[iL03,03],DistribParam2[iL04,04],DistribParam2[iL05,05],
-            DistribParam2[iL06,06],DistribParam2[iL07,07],DistribParam2[iL08,08],DistribParam2[iL09,09],DistribParam2[iL10,10],
-            DistribParam2[iL11,11],DistribParam2[iL12,12],DistribParam2[iL13,13],DistribParam2[iL14,14],DistribParam2[iL15,15],
-            DistribParam2[iL16,16],DistribParam2[iL17,17],DistribParam2[iL18,18],DistribParam2[iL19,19],DistribParam2[iL20,20]);
-        } } } } }
-        } } } } }
-        } } } } }
-        } } } } }
-        MAT <- matrix(VECT,ncol=20,byrow=TRUE)[,1:NParam];
-        if(is.matrix(MAT)==FALSE){ MAT <- cbind(MAT); }
-        Output <- NULL;
-        Output$NewCandidates <- MAT;
-        return(Output);
+    ##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
 
 
-    ##Creation_of_new_candidates_______________________________________________
-    if(PrefilteringType==1){ CandidatesParamR <- CalibOptions$StartParamList; }
-    if(PrefilteringType==2){ DistribParamR <- CalibOptions$StartParamDistrib; DistribParamR[,!CalibOptions$OptimParam] <- NA; CandidatesParamR <- ProposeCandidatesGrid(DistribParamR)$NewCandidates; }
-    ##Remplacement_of_non_optimised_values_____________________________________
-    CandidatesParamR <- apply(CandidatesParamR,1,function(x){ x[!CalibOptions$OptimParam] <- CalibOptions$FixedParam[!CalibOptions$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(!quiet & Ncandidates>1){
-      if(PrefilteringType==1){ cat(paste("\t List-Screening in progress (",sep="")); }
-      if(PrefilteringType==2){ cat(paste("\t Grid-Screening in progress (",sep="")); }
-      cat("0%");
-    }
-    for(iNew in 1:nrow(CandidatesParamR)){
-      if(!quiet & Ncandidates>1){
-        for(k in c(2,4,6,8)){ if(iNew==round(k/10*Ncandidates)){ cat(paste(" ",10*k,"%",sep="")); } } 
-      }
-      ##Model_run
-      Param <- CandidatesParamR[iNew,];
-      OutputsModel <- FUN_MOD(InputsModel,RunOptions,Param);
-      ##Calibration_criterion_computation
-      OutputsCrit <- FUN_CRIT(InputsCrit,OutputsModel);      
-      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;
-      }
+
+
+  ##_____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(!quiet & Ncandidates>1){ cat(" 100%) \n"); }
-	  
-
-    ##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(!quiet){
-      if(Ncandidates> 1){ cat(paste("\t Screening completed (",NRuns," runs): \n",sep=""));  }
-      if(Ncandidates==1){ cat(paste("\t Starting point for steepest-descent local search: \n",sep="")); }
-      cat(paste("\t     Param = ",paste(formatC(ParamStartR,format="f",width=8,digits=3),collapse=" , "),"\n",sep=""));
-      cat(paste("\t     Crit ",format(CritName,width=12,justify="left")," = ",formatC(CritStart*Multiplier,format="f",digits=4),"\n",sep=""));
+    if (ncol(NewParamOptimT)!=ncol(OldParamOptimT) | ncol(NewParamOptimT) != length(OptimParam)) {
+      stop("each input set must have the same number of values")
     }
-    ##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 \n"); return(NULL); }
-      if(ncol(NewParamOptimT)!=ncol(OldParamOptimT) | ncol(NewParamOptimT)!=length(OptimParam)){ stop("each input set must have the same number of values \n"); return(NULL); }
-      ##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]==TRUE){
-          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==TRUE){ VECT <- c(VECT,PotentialCandidateT); }
+    ##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);
     }
-      
+    Output <- NULL
+    Output$NewCandidatesT <- matrix(VECT, ncol = NParam, byrow = TRUE)
+    return(Output)
+  }
 
-    ##Initialisation_of_variables
-    if(!quiet){ 
-      cat("\t Steepest-descent local search in progress \n"); 
-    }
-    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;
 
+  ##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)){
+
+  ##START_LOOP_ITER_________________________________________________________
+  for (ITER in 1:(100 * NParam)) {
 
 
     ##Exit_loop_when_Pace_becomes_too_small___________________________________
-    if(Pace<0.01){ break; }
-  
+    if (Pace < 0.01) {
+      break
+    }
+
 
     ##Creation_of_new_candidates______________________________________________
-    CandidatesParamT <- ProposeCandidatesLoc(NewParamOptimT,OldParamOptimT,RangesT,CalibOptions$OptimParam,Pace)$NewCandidatesT;
-    CandidatesParamR <- FUN_TRANSFO(CandidatesParamT,"TR");
+    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[!CalibOptions$OptimParam] <- CalibOptions$FixedParam[!CalibOptions$OptimParam]; return(x); });
-    if(NParam>1){ CandidatesParamR <- t(CandidatesParamR); } else { CandidatesParamR <- cbind(CandidatesParamR); }
+    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)){
+    iNewOptim <- 0
+    for (iNew in 1:nrow(CandidatesParamR)) {
       ##Model_run
-      Param <- CandidatesParamR[iNew,];
-      OutputsModel <- FUN_MOD(InputsModel,RunOptions,Param);
+      Param <- CandidatesParamR[iNew, ]
+      OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, ...)
       ##Calibration_criterion_computation
-      OutputsCrit <- FUN_CRIT(InputsCrit,OutputsModel);      
-      if(!is.na(OutputsCrit$CritValue)){ if(OutputsCrit$CritValue*OutputsCrit$Multiplier < CritOptim){
-        CritOptim <- OutputsCrit$CritValue*OutputsCrit$Multiplier;
-        iNewOptim <- iNew;
-      } }
+      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);
+    NRuns <- NRuns + nrow(CandidatesParamR)
 
 
     ##When_a_progress_has_been_achieved_______________________________________
-    if(iNewOptim!=0){
+    if (iNewOptim != 0) {
       ##We_store_the_optimal_set
-      OldParamOptimT <- NewParamOptimT;
-      NewParamOptimT <- matrix(CandidatesParamT[iNewOptim,1:NParam],nrow=1);
-      Compt <- Compt+1;
+      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;
+      if (Compt > 2 * NParam) {
+        Pace <- Pace * 2
+        Compt <- 0
       }
       ##We_update_PaceDiag
-      VectPace <- NewParamOptimT-OldParamOptimT;
-      for(iC in 1:NParam){ if(CalibOptions$OptimParam[iC]==TRUE){ 
-        if(VectPace[iC]!=0){ PaceDiag[iC] <- CLG*PaceDiag[iC]+(1-CLG)*VectPace[iC]; }
-        if(VectPace[iC]==0){ PaceDiag[iC] <- CLG*PaceDiag[iC]; }
-      } }
+      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;
+      ##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(CalibOptions$OptimParam[iC]==TRUE){
-          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 <- FUN_MOD(InputsModel,RunOptions,Param);
-        ##Calibration_criterion_computation
-        OutputsCrit <- FUN_CRIT(InputsCrit,OutputsModel);      
-        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);
+    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(!quiet){ 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 & !quiet){ 
-      cat("\t No progress achieved \n"); 
-    }
-    
-    ##End_of_Steepest_Descent_Local_Search____________________________________
-    ParamFinalR <- NewParamOptimR;
-    ParamFinalT <- NewParamOptimT;
-    CritFinal   <- CritOptim;
-    NIter       <- 1+ITER;
-    if(!quiet){ 
-      cat(paste("\t Calibration completed (",NIter," iterations, ",NRuns," runs): \n",sep=""));
-      cat(paste("\t     Param = ",paste(formatC(ParamFinalR,format="f",width=8,digits=3),collapse=" , "),"\n",sep=""));
-      cat(paste("\t     Crit ",format(CritName,width=12,justify="left")," = ",formatC(CritFinal*Multiplier,format="f",digits=4),"\n",sep=""));
-    }
     ##Results_archiving_______________________________________________________
-  	HistParamR <- cbind(HistParamR[1:NIter,]); colnames(HistParamR) <- paste("Param",1:NParam,sep="");
-  	HistParamT <- cbind(HistParamT[1:NIter,]); colnames(HistParamT) <- paste("Param",1:NParam,sep="");
-	  HistCrit   <- cbind(HistCrit[1:NIter,]);   ###colnames(HistCrit) <- paste("HistCrit");
+    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");
+  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(as.double(ParamFinalR),CritFinal*Multiplier,NIter,NRuns,HistParamR,HistCrit*Multiplier,MatBoolCrit,CritName,CritBestValue);
-    names(OutputsCalib) <- c("ParamFinalR","CritFinal","NIter","NRuns","HistParamR","HistCrit","MatBoolCrit","CritName","CritBestValue");
-    class(OutputsCalib) <- c("OutputsCalib","HBAN");
-    return(OutputsCalib);
+  ##_____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)
 
 
 
 }
-
-
-
-  
-

R/Calibration_optim.R

-#*************************************************************************************************
-#' Calibration algorithm which minimises the error criterion. \cr
-#' \cr
-#' The algorithm is based on the "optim" function from the "stats" R-package 
-#' (using method="L-BFGS-B", i.e. a local optimization quasi-Newton method).
-#'
-#' To optimise the exploration of the parameter space, transformation functions are used to convert
-#' the model parameters. This is done using the TransfoParam functions.
-#*************************************************************************************************
-#' @title  Calibration algorithm which minimises the error criterion using the stats::optim function
-#' @author Laurent Coron (August 2013)
-#' @example tests/example_Calibration_optim.R
-#' @seealso \code{\link{Calibration}}, \code{\link{Calibration_Michel}},
-#'          \code{\link{RunModel_GR4J}}, \code{\link{TransfoParam_GR4J}}, \code{\link{ErrorCrit_RMSE}},
-#'          \code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}}, 
-#'          \code{\link{CreateInputsCrit}}, \code{\link{CreateCalibOptions}}.
-#' @encoding UTF-8
-#' @export
-#_FunctionInputs__________________________________________________________________________________
-#' @param  InputsModel         [object of class \emph{InputsModel}] see \code{\link{CreateInputsModel}} for details
-#' @param  RunOptions          [object of class \emph{RunOptions}] see \code{\link{CreateRunOptions}} for details
-#' @param  InputsCrit          [object of class \emph{InputsCrit}] see \code{\link{CreateInputsCrit}} for details
-#' @param  CalibOptions        [object of class \emph{CalibOptions}] see \code{\link{CreateCalibOptions}} for details
-#' @param  FUN_MOD             [function] hydrological model function (e.g. RunModel_GR4J, RunModel_CemaNeigeGR4J)
-#' @param  FUN_CRIT            [function] error criterion function (e.g. ErrorCrit_RMSE, ErrorCrit_NSE)
-#' @param  FUN_TRANSFO         (optional) [function] model parameters transformation function, if the FUN_MOD used is native in the package FUN_TRANSFO is automatically defined
-#' @param  quiet               (optional) [boolean] boolean indicating if the function is run in quiet mode or not, default=FALSE
-#_FunctionOutputs_________________________________________________________________________________
-#' @return  [list] list containing the function outputs organised as follows:
-#'          \tabular{ll}{
-#'          \emph{$ParamFinalR  }  \tab   [numeric] parameter set obtained at the end of the calibration \cr
-#'          \emph{$CritFinal    }  \tab   [numeric] error criterion obtained at the end of the calibration \cr
-#'          \emph{$Nruns        }  \tab   [numeric] number of model runs done during the calibration \cr
-#'          \emph{$CritName     }  \tab   [character] name of the calibration criterion \cr
-#'          \emph{$CritBestValue}  \tab   [numeric] theoretical best criterion value \cr
-#'          }
-#**************************************************************************************************
-Calibration_optim <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO=NULL,quiet=FALSE){
-    
-
-    ##_check_class
-    if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
-    if(inherits(RunOptions,"RunOptions")==FALSE){ stop("RunOptions must be of class 'RunOptions' \n"); return(NULL); }  
-    if(inherits(InputsCrit,"InputsCrit")==FALSE){ stop("InputsCrit must be of class 'InputsCrit' \n"); return(NULL); }  
-    if(inherits(CalibOptions,"CalibOptions")==FALSE){ stop("CalibOptions must be of class 'CalibOptions' \n"); return(NULL); }  
-    if(inherits(CalibOptions,"optim")==FALSE){ stop("CalibOptions must be of class 'optim' if Calibration_optim is used \n"); return(NULL); }  
-
-
-    ##_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    )){ FUN_TRANSFO <- TransfoParam_CemaNeige; }
-      if(identical(FUN_MOD,RunModel_CemaNeigeGR4J) | identical(FUN_MOD,RunModel_CemaNeigeGR5J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)){
-        if(identical(FUN_MOD,RunModel_CemaNeigeGR4J)){ FUN1 <- TransfoParam_GR4J; FUN2 <- TransfoParam_CemaNeige; }
-        if(identical(FUN_MOD,RunModel_CemaNeigeGR5J)){ FUN1 <- TransfoParam_GR5J; FUN2 <- TransfoParam_CemaNeige; }
-        if(identical(FUN_MOD,RunModel_CemaNeigeGR6J)){ FUN1 <- TransfoParam_GR6J; FUN2 <- TransfoParam_CemaNeige; }
-        FUN_TRANSFO <- function(ParamIn,Direction){
-          Bool <- is.matrix(ParamIn);
-          if(Bool==FALSE){ 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==FALSE){ ParamOut <- ParamOut[1,]; }
-          return(ParamOut);
-        }
-      }
-      if(is.null(FUN_TRANSFO)){ stop("FUN_TRANSFO was not found (in Calibration function) \n"); return(NULL);  }
-    }
-
-
-    ##_RunModelAndCrit
-    RunModelAndCrit <- function(par,InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO){
-        ParamT <- NA*CalibOptions$FixedParam;
-        ParamT[CalibOptions$OptimParam] <- par;
-        Param <- FUN_TRANSFO(ParamIn=ParamT,Direction="TR");
-        Param[!CalibOptions$OptimParam] <- CalibOptions$FixedParam[!CalibOptions$OptimParam];
-        OutputsModel <- FUN_MOD(InputsModel=InputsModel,RunOptions=RunOptions,Param=Param);
-        OutputsCrit <- FUN_CRIT(InputsCrit=InputsCrit,OutputsModel=OutputsModel);   
-        return(OutputsCrit$CritValue*OutputsCrit$Multiplier);
-    }
-
-
-    ##_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
-    ##_screenPrint
-    if(!quiet){ 
-      cat(paste("\t Calibration in progress (function optim from the stats package) \n",sep=""));
-    }
-    
-
-    ##_lower_and_upper_limit_values (transformed)
-    RangesR <- CalibOptions$SearchRanges;
-    RangesT <- FUN_TRANSFO(RangesR,"RT");
-    lower   <- RangesT[1,CalibOptions$OptimParam];
-    upper   <- RangesT[2,CalibOptions$OptimParam];
-
-    ##_starting_values (transformed)
-    ParamStartT <- FUN_TRANSFO(CalibOptions$StartParam,"RT");
-    par_start   <- ParamStartT[CalibOptions$OptimParam]; 
-
-
-    ##_calibration
-    RESULT <- optim(par=par_start,fn=RunModelAndCrit,gr=NULL,
-                InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO,  ## arguments for the RunModelAndCrit function (other than par)
-                method="L-BFGS-B",lower=lower,upper=upper,control=list(),hessian=FALSE)
-    
-
-    ##_outputs_preparation
-    ParamFinalT <- NA*ParamStartT;
-    ParamFinalT[CalibOptions$OptimParam] <- RESULT$par;
-    ParamFinalR <- FUN_TRANSFO(ParamFinalT,"TR");
-    ParamFinalR[!CalibOptions$OptimParam] <- CalibOptions$FixedParam[!CalibOptions$OptimParam];  
-    CritFinal <- RESULT$value;
-
-    ##_storage_of_crit_info
-    OutputsModel  <- FUN_MOD(InputsModel=InputsModel,RunOptions=RunOptions,Param=ParamFinalR);
-    OutputsCrit   <- FUN_CRIT(InputsCrit=InputsCrit,OutputsModel=OutputsModel);   
-    CritName      <- OutputsCrit$CritName;
-    CritBestValue <- OutputsCrit$CritBestValue;
-    Multiplier    <- OutputsCrit$Multiplier;
-
-    ##_screenPrint
-    if(!quiet){ 
-      if(RESULT$convergence==0){ 
-        cat(paste("\t Calibration completed: \n",sep=""));
-        cat(paste("\t     Param = ",paste(formatC(ParamFinalR,format="f",width=8,digits=3),collapse=" , "),"\n",sep=""));
-        cat(paste("\t     Crit ",format(CritName,width=12,justify="left")," = ",formatC(CritFinal*Multiplier,format="f",digits=4),"\n",sep=""));
-      } else {
-        cat(paste("\t Calibration failed: \n",sep=""));
-        cat(paste("\t     ",RESULT$message,sep=""));
-      }
-    }
-    
-    
-    ##_function_output
-    OutputsCalib <- list(as.double(ParamFinalR),CritFinal*Multiplier,as.integer(RESULT$counts[1]),CritName,CritBestValue);
-    names(OutputsCalib) <- c("ParamFinalR","CritFinal","NRuns","CritName","CritBestValue");
-    class(OutputsCalib) <- c("OutputsCalib","optim");
-    return(OutputsCalib);
-
-
-}
-
-
-
-

R/CreateCalibOptions.R

-#*************************************************************************************************
-#' Creation of the CalibOptions object required to the Calibration functions.
-#'
-#' Users wanting to use FUN_MOD, FUN_CALIB or FUN_TRANSFO functions that are not included in 
-#' the package must create their own CalibOptions object accordingly.
-#*************************************************************************************************
-#' @title  Creation of the CalibOptions object required to the Calibration functions
-#' @author Laurent Coron (June 2014)
-#' @seealso \code{\link{Calibration}}, \code{\link{RunModel}}
-#' @example tests/example_Calibration.R
-#' @encoding UTF-8
-#' @export
-#_FunctionInputs__________________________________________________________________________________
-#' @param  FUN_MOD             [function] hydrological model function (e.g. RunModel_GR4J, RunModel_CemaNeigeGR4J)
-#' @param  FUN_CALIB           (optional) [function] calibration algorithm function (e.g. Calibration_HBAN, Calibration_optim), default=Calibration_HBAN
-#' @param  FUN_TRANSFO         (optional) [function] model parameters transformation function, if the FUN_MOD used is native in the package FUN_TRANSFO is automatically defined
-#' @param  OptimParam          (optional) [boolean] vector of booleans indicating which parameters must be optimised (NParam columns, 1 line) 
-#' @param  FixedParam          (optional) [numeric] vector giving the values to allocate to non-optimised parameter values (NParam columns, 1 line)
-#' @param  SearchRanges        (optional) [numeric] matrix giving the ranges of real parameters (NParam columns, 2 lines)
-#'                             \tabular{llllll}{
-#'                                           \tab [X1]   \tab [X2]   \tab [X3]   \tab [...]   \tab [Xi] \cr
-#'                                  [1,]     \tab    0   \tab   -1   \tab    0   \tab  ...    \tab  0.0 \cr
-#'                                  [2,]     \tab 3000   \tab   +1   \tab  100   \tab  ...    \tab  3.0 \cr
-#'                             }
-#' @param  StartParam          (optional) [numeric] vector of parameter values used to start global search calibration procedure (this argument is used by Calibration_optim but not by Calibration_HBAN)
-#'                             \tabular{llllll}{
-#'                                           \tab [X1]   \tab [X2]   \tab [X3]   \tab [...]   \tab [Xi] \cr
-#'                                           \tab 1000   \tab -0.5   \tab   22   \tab  ...    \tab  1.1 \cr
-#'                             }
-#' @param  StartParamList      (optional) [numeric] matrix of parameter sets used for grid-screening calibration procedure (values in columns, sets in line) (this argument is used by Calibration_HBAN but not by Calibration_optim)
-#'                             \tabular{llllll}{
-#'                                           \tab [X1]   \tab [X2]   \tab [X3]   \tab [...]   \tab [Xi] \cr
-#'                                  [set1]   \tab  800   \tab -0.7   \tab   25   \tab  ...    \tab  1.0 \cr
-#'                                  [set2]   \tab 1000   \tab -0.5   \tab   22   \tab  ...    \tab  1.1 \cr
-#'                                  [...]    \tab  ...   \tab  ...   \tab  ...   \tab  ...    \tab  ... \cr
-#'                                  [set n]  \tab  200   \tab -0.3   \tab   17   \tab  ...    \tab  1.0 \cr
-#'                             }
-#' @param  StartParamDistrib   (optional) [numeric] matrix of parameter values used for grid-screening calibration procedure (values in columns, percentiles in line) \cr
-#'                             \tabular{llllll}{
-#'                                           \tab [X1]   \tab [X2]   \tab [X3]   \tab [...]   \tab [Xi] \cr
-#'                                  [value1] \tab  800   \tab -0.7   \tab   25   \tab  ...    \tab  1.0 \cr
-#'                                  [value2] \tab 1000   \tab   NA   \tab   50   \tab  ...    \tab  1.2 \cr
-#'                                  [value3] \tab 1200   \tab   NA   \tab   NA   \tab  ...    \tab  1.6 \cr
-#'                             }
-#_FunctionOutputs_________________________________________________________________________________
-#' @return  [list] object of class \emph{CalibOptions} containing the data required to evaluate the model outputs; it can include the following:
-#'          \tabular{ll}{
-#'          \emph{$OptimParam       }  \tab   [boolean] vector of booleans indicating which parameters must be optimised \cr
-#'          \emph{$FixedParam       }  \tab   [numeric] vector giving the values to allocate to non-optimised parameter values \cr
-#'          \emph{$SearchRanges     }  \tab   [numeric] matrix giving the ranges of real parameters \cr
-#'          \emph{$StartParam       }  \tab   [numeric] vector of parameter values used to start global search calibration procedure \cr
-#'          \emph{$StartParamList   }  \tab   [numeric] matrix of parameter sets used for grid-screening calibration procedure \cr
-#'          \emph{$StartParamDistrib}  \tab   [numeric] matrix of parameter values used for grid-screening calibration procedure \cr
-#'          }
-#**************************************************************************************************
-CreateCalibOptions <- function(FUN_MOD,FUN_CALIB=Calibration_HBAN,FUN_TRANSFO=NULL,OptimParam=NULL,FixedParam=NULL,SearchRanges=NULL,
-                               StartParam=NULL,StartParamList=NULL,StartParamDistrib=NULL){
-
-
-  ObjectClass <- NULL;
-
-  ##check_FUN_MOD
-    BOOL <- FALSE;
-    if(identical(FUN_MOD,RunModel_GR4H         )){ ObjectClass <- c(ObjectClass,"GR4H"                 ); BOOL <- TRUE; }
-    if(identical(FUN_MOD,RunModel_GR4J         )){ ObjectClass <- c(ObjectClass,"GR4J"         ); BOOL <- TRUE; }
-    if(identical(FUN_MOD,RunModel_GR5J         )){ ObjectClass <- c(ObjectClass,"GR5J"         ); BOOL <- TRUE; }
-    if(identical(FUN_MOD,RunModel_GR6J         )){ ObjectClass <- c(ObjectClass,"GR6J"         ); BOOL <- TRUE; }
-    if(identical(FUN_MOD,RunModel_GR2M         )){ ObjectClass <- c(ObjectClass,"GR2M"                 ); BOOL <- TRUE; }
-    if(identical(FUN_MOD,RunModel_GR1A         )){ ObjectClass <- c(ObjectClass,"GR1A"                 ); BOOL <- TRUE; }
-    if(identical(FUN_MOD,RunModel_CemaNeige    )){ ObjectClass <- c(ObjectClass,"CemaNeige"    ); BOOL <- TRUE; }
-    if(identical(FUN_MOD,RunModel_CemaNeigeGR4J)){ ObjectClass <- c(ObjectClass,"CemaNeigeGR4J"); BOOL <- TRUE; }
-    if(identical(FUN_MOD,RunModel_CemaNeigeGR5J)){ ObjectClass <- c(ObjectClass,"CemaNeigeGR5J"); BOOL <- TRUE; }
-    if(identical(FUN_MOD,RunModel_CemaNeigeGR6J)){ ObjectClass <- c(ObjectClass,"CemaNeigeGR6J"); BOOL <- TRUE; }
-    if(!BOOL){ stop("incorrect FUN_MOD for use in CreateCalibOptions \n"); return(NULL); } 
-
-  ##check_FUN_CALIB
-    BOOL <- FALSE;
-    if(identical(FUN_CALIB,Calibration_HBAN  )){ ObjectClass <- c(ObjectClass,"HBAN"  ); BOOL <- TRUE; }
-    if(identical(FUN_CALIB,Calibration_optim )){ ObjectClass <- c(ObjectClass,"optim" ); BOOL <- TRUE; }
-    if(!BOOL){ stop("incorrect FUN_CALIB for use in CreateCalibOptions \n"); return(NULL); } 
-
-  ##check_FUN_TRANSFO
-    if(is.null(FUN_TRANSFO)){
-      ##_set_FUN1
-      if(identical(FUN_MOD,RunModel_GR4H     )                                             ){ FUN1 <- TransfoParam_GR4H     ; }
-      if(identical(FUN_MOD,RunModel_GR4J     ) | identical(FUN_MOD,RunModel_CemaNeigeGR4J) ){ FUN1 <- TransfoParam_GR4J     ; }
-      if(identical(FUN_MOD,RunModel_GR5J     ) | identical(FUN_MOD,RunModel_CemaNeigeGR5J) ){ FUN1 <- TransfoParam_GR5J     ; }
-      if(identical(FUN_MOD,RunModel_GR6J     ) | identical(FUN_MOD,RunModel_CemaNeigeGR6J) ){ FUN1 <- TransfoParam_GR6J     ; }
-      if(identical(FUN_MOD,RunModel_GR2M     )                                             ){ FUN1 <- TransfoParam_GR2M     ; }
-      if(identical(FUN_MOD,RunModel_GR1A     )                                             ){ FUN1 <- TransfoParam_GR1A     ; }
-      if(identical(FUN_MOD,RunModel_CemaNeige)                                             ){ FUN1 <- TransfoParam_CemaNeige; }
-      if(is.null(FUN1)){ stop("FUN1 was not found \n"); return(NULL);  }
-      ##_set_FUN2
-      FUN2 <- TransfoParam_CemaNeige;
-      ##_set_FUN_TRANSFO
-      if(sum(ObjectClass %in% c("GR4H","GR4J","GR5J","GR6J","GR2M","GR1A","CemaNeige"))>0){
-        FUN_TRANSFO <- FUN1; 
-      } else {
-        FUN_TRANSFO <- function(ParamIn,Direction){
-          Bool <- is.matrix(ParamIn);
-          if(Bool==FALSE){ ParamIn <- rbind(ParamIn); }
-          ParamOut <- NA*ParamIn;
-          NParam   <- ncol(ParamIn);
-          if(NParam <= 3){
-          ParamOut[,         1:(NParam-2)] <- FUN1(cbind(ParamIn[,1:(NParam-2)]),Direction);
-          } else {
-          ParamOut[,         1:(NParam-2)] <- FUN1(ParamIn[,1:(NParam-2)     ],Direction); }
-          ParamOut[,(NParam-1):NParam    ] <- FUN2(ParamIn[,(NParam-1):NParam],Direction);
-          if(Bool==FALSE){ ParamOut <- ParamOut[1,]; }
-          return(ParamOut);
-        }
-      }
-    }
-    if(is.null(FUN_TRANSFO)){ stop("FUN_TRANSFO was not found \n"); return(NULL);  }
-
-  ##NParam
-    if("GR4H"          %in% ObjectClass){ NParam <- 4; }
-    if("GR4J"          %in% ObjectClass){ NParam <- 4; }
-    if("GR5J"          %in% ObjectClass){ NParam <- 5; }
-    if("GR6J"          %in% ObjectClass){ NParam <- 6; }
-    if("GR2M"          %in% ObjectClass){ NParam <- 2; }
-    if("GR1A"          %in% ObjectClass){ NParam <- 1; }
-    if("CemaNeige"     %in% ObjectClass){ NParam <- 2; }
-    if("CemaNeigeGR4J" %in% ObjectClass){ NParam <- 6; }
-    if("CemaNeigeGR5J" %in% ObjectClass){ NParam <- 7; }
-    if("CemaNeigeGR6J" %in% ObjectClass){ NParam <- 8; }
-
-  ##check_OptimParam
-    if(is.null(OptimParam)){
-      OptimParam <- rep(TRUE,NParam);
-    } else {
-      if(!is.vector(OptimParam)    ){ stop("OptimParam must be a vector of booleans \n"); return(NULL);  }
-      if(length(OptimParam)!=NParam){ stop("Incompatibility between OptimParam length and FUN_MOD \n"); return(NULL);  }
-      if(!is.logical(OptimParam)   ){ stop("OptimParam must be a vector of booleans \n"); return(NULL);  }
-    }
-
-  ##check_FixedParam
-    if(is.null(FixedParam)){
-      FixedParam <- rep(NA,NParam);
-    } else {
-      if(!is.vector(FixedParam)              ){ stop("FixedParam must be a vector \n"); return(NULL);  }
-      if(length(FixedParam)!=NParam          ){ stop("Incompatibility between OptimParam length and FUN_MOD \n"); return(NULL);  }
-      if(sum(!OptimParam)>0){
-      if(!is.numeric(FixedParam[!OptimParam])){ stop("if OptimParam[i]==FALSE, FixedParam[i] must be a numeric value \n"); return(NULL);  } }
-    }
-
-  ##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 \n"); return(NULL);  }
-      if(!is.numeric(SearchRanges)  ){ stop("SearchRanges must be a matrix of numeric values \n"); return(NULL);  }
-      if(sum(is.na(SearchRanges))!=0){ stop("SearchRanges must not include NA values \n"); return(NULL);  }
-      if(nrow(SearchRanges)!=2      ){ stop("SearchRanges must have 2 rows \n"); return(NULL);  }
-      if(ncol(SearchRanges)!=NParam ){ stop("Incompatibility between SearchRanges ncol and FUN_MOD \n"); return(NULL);  }
-    }
-
-  ##check_StartParamList_and_StartParamDistrib__default_values
-    if( ("HBAN"  %in% ObjectClass & is.null(StartParamList) & is.null(StartParamDistrib)) |
-        ("optim" %in% ObjectClass & is.null(StartParam)) ){
-
-      if("GR4H"%in% ObjectClass){ 
-        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=NParam,byrow=TRUE); }
-      if("GR4J"%in% ObjectClass){ 
-        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=NParam,byrow=TRUE); }
-      if("GR5J"%in% ObjectClass){ 
-        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=NParam,byrow=TRUE); }
-      if("GR6J"%in% ObjectClass){ 
-        ParamT <- matrix( c( +4.41, +0.41, +2.88, -9.10, -0.13, +0.81,
-                             +5.02, +0.61, +3.45, -8.68, +1.95, +2.27,
-                             +5.58, +0.78, +4.18, -8.12, +3.59, +3.56),ncol=NParam,byrow=TRUE); }
-      if("GR2M"%in% ObjectClass){ 
-        ParamT <- matrix( c( +5.03, -7.15,
-                             +5.22, -6.74,
-                             +5.85, -6.37),ncol=NParam,byrow=TRUE); }
-      if("GR1A"%in% ObjectClass){ 
-        ParamT <- matrix( c( -1.69,
-                             -0.38,
-                             +1.39),ncol=NParam,byrow=TRUE); }
-      if("CemaNeige"%in% ObjectClass){ 
-        ParamT <- matrix( c( -9.96, +6.63,
-                             -9.14, +6.90,
-                             +4.10, +7.21),ncol=NParam,byrow=TRUE); }
-      if("CemaNeigeGR4J"%in% ObjectClass){ 
-        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=NParam,byrow=TRUE); }
-      if("CemaNeigeGR5J"%in% ObjectClass){ 
-        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=NParam,byrow=TRUE); }
-      if("CemaNeigeGR6J"%in% ObjectClass){ 
-        ParamT <- matrix( c( +4.41, +0.41, +2.88, -9.10, -0.13, +0.81, -9.96, +6.63,
-                             +5.02, +0.61, +3.45, -8.68, +1.95, +2.27, -9.14, +6.90,
-                             +5.58, +0.78, +4.18, -8.12, +3.59, +3.56, +4.10, +7.21),ncol=NParam,byrow=TRUE); }
-
-      StartParamList    <- NULL;
-      StartParamDistrib <- TransfoParam(ParamIn=ParamT,Direction="TR",FUN_TRANSFO=FUN_TRANSFO);
-      StartParam        <- StartParamDistrib[2,];
-    }
-
-  ##check_StartParamList_and_StartParamDistrib__format
-    if("HBAN" %in% ObjectClass & !is.null(StartParamList)){
-      if(!is.matrix( StartParamList)  ){ stop("StartParamList must be a matrix \n"); return(NULL);  }
-      if(!is.numeric(StartParamList)  ){ stop("StartParamList must be a matrix of numeric values \n"); return(NULL);  }
-      if(sum(is.na(StartParamList))!=0){ stop("StartParamList must not include NA values \n"); return(NULL);  }
-      if(ncol(StartParamList)!=NParam ){ stop("Incompatibility between StartParamList ncol and FUN_MOD \n"); return(NULL);  }
-    }
-    if("HBAN" %in% ObjectClass & !is.null(StartParamDistrib)){
-      if(!is.matrix( StartParamDistrib)      ){ stop("StartParamDistrib must be a matrix \n"); return(NULL);  }
-      if(!is.numeric(StartParamDistrib[1,])  ){ stop("StartParamDistrib must be a matrix of numeric values \n"); return(NULL);  }
-      if(sum(is.na(StartParamDistrib[1,]))!=0){ stop("StartParamDistrib must not include NA values on the first line \n"); return(NULL);  }
-      if(ncol(StartParamDistrib)!=NParam     ){ stop("Incompatibility between StartParamDistrib ncol and FUN_MOD \n"); return(NULL);  }
-    }
-    if("optim" %in% ObjectClass & !is.null(StartParam)){
-      if(!is.vector( StartParam)    ){ stop("StartParam must be a vector \n"); return(NULL);  }
-      if(!is.numeric(StartParam)    ){ stop("StartParam must be a vector of numeric values \n"); return(NULL);  }
-      if(sum(is.na(StartParam))!=0  ){ stop("StartParam must not include NA values \n"); return(NULL);  }
-      if(length(StartParam)!=NParam ){ stop("Incompatibility between StartParam length and FUN_MOD \n"); return(NULL);  }
+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)
 
-  ##Create_CalibOptions
-    CalibOptions <- list(OptimParam=OptimParam,FixedParam=FixedParam,SearchRanges=SearchRanges);
-    if(!is.null(StartParam       )){ CalibOptions <- c(CalibOptions,list(StartParam=StartParam)); }
-    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);
+  } 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(-8.75,
+                           -7.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)
+
+}

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, 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
+
+  FUN_MOD <- match.fun(FUN_MOD)
+
+  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'")
+  }
+  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'")
+  }
+
+
+  ## 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", 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", FeatFUN_MOD$NameFunMod))
+      }
+      UH1 <- rep(Inf, UH1n)
+    }
+    if (!is.null(UH2)) {
+      if (verbose) {
+        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", 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", FeatFUN_MOD$NameFunMod))
+      }
+    }
+    ProdStore <- Inf
+    if (!is.null(RoutStore)) {
+      if (verbose) {
+        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", 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", FeatFUN_MOD$NameFunMod))
+      }
+    }
+    UH1 <- rep(Inf, UH1n)
+    if (!is.null(UH2)) {
+      if (verbose) {
+        warning(sprintf("'%s' does not require 'UH2'. Values set to NA", FeatFUN_MOD$NameFunMod))
+      }
+    }
+    UH2 <- rep(Inf, UH2n)
+  }
+  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'", FeatFUN_MOD$NameFunMod))
+  }
+  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'", FeatFUN_MOD$NameFunMod))
+  }
+  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", FeatFUN_MOD$NameFunMod))
+    }
+    GthrCemaNeigeLayers    <- Inf
+    GlocmaxCemaNeigeLayers <- Inf
+  }
+  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", FeatFUN_MOD$NameFunMod))
+    }
+    GCemaNeigeLayers       <- Inf
+    eTGCemaNeigeLayers     <- Inf
+    GthrCemaNeigeLayers    <- Inf
+    GlocmaxCemaNeigeLayers <- Inf
+  }
+
+
+  ## set states
+  if ("CemaNeige" %in% ObjectClass) {
+    NLayers <- length(InputsModel$LayerPrecip)
+  } else {
+    NLayers <- 1
+  }
+
+
+  ## manage NULL values
+  if (is.null(ExpStore)) {
+    ExpStore <- Inf
+  }
+  if (is.null(IntStore)) {
+    IntStore <- Inf
+  }
+  if (is.null(UH1)) {
+    if ("hourly"  %in% ObjectClass) {
+      k <- 24
+    } else {
+      k <- 1
+    }
+    UH1 <- rep(Inf, UH1n * k)
+  }
+  if (is.null(UH2)) {
+    if ("hourly"  %in% ObjectClass) {
+      k <- 24
+    } else {
+      k <- 1
+    }
+    UH2 <- rep(Inf, UH2n * k)
+  }
+  if (is.null(GCemaNeigeLayers)) {
+    GCemaNeigeLayers <- rep(Inf, NLayers)
+  }
+  if (is.null(eTGCemaNeigeLayers)) {
+    eTGCemaNeigeLayers <- rep(Inf, NLayers)
+  }
+  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) | any(IntStore < 0) |
+      any(UH1 < 0) | any(UH2 < 0) |
+      any(GCemaNeigeLayers < 0)) {
+    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")
+  }
+  if (!is.numeric(RoutStore) || length(RoutStore) != 1L) {
+    stop("'RoutStore' must be numeric of length one")
+  }
+  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))
+  }
+  if (!"hourly" %in% ObjectClass & (!is.numeric(UH1) || length(UH1) != UH1n)) {
+    stop(sprintf("'UH1' must be numeric of length %i", UH1n))
+  }
+  if ( "hourly" %in% ObjectClass & (!is.numeric(UH2) || length(UH2) != UH2n * 24)) {
+    stop(sprintf("'UH2' must be numeric of length 960 (%i * 24)", UH2n))
+  }
+  if (!"hourly" %in% ObjectClass & (!is.numeric(UH2) || length(UH2) != UH2n)) {
+    stop(sprintf("'UH2' must be numeric of length %i (2 * %i)", UH2n, UH1n))
+  }
+  if (!is.numeric(GCemaNeigeLayers) || length(GCemaNeigeLayers) != NLayers) {
+    stop(sprintf("'GCemaNeigeLayers' must be numeric of length %i", NLayers))
+  }
+  if (!is.numeric(eTGCemaNeigeLayers) || length(eTGCemaNeigeLayers) != NLayers) {
+    stop(sprintf("'eTGCemaNeigeLayers' must be numeric of length %i", NLayers))
+  }
+  if (IsHyst) {
+    if (!is.numeric(GthrCemaNeigeLayers) || length(GthrCemaNeigeLayers) != NLayers) {
+      stop(sprintf("'eTGCemaNeigeLayers' must be numeric of length %i", NLayers))
+    }
+    if (!is.numeric(GlocmaxCemaNeigeLayers) || length(GlocmaxCemaNeigeLayers) != NLayers) {
+      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, 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) {
+    class(IniStatesNA) <- c(class(IniStatesNA), "hysteresis")
+  }
+  if (IsIntStore) {
+    class(IniStatesNA) <- c(class(IniStatesNA), "interception")
+  }
+
+  return(IniStatesNA)
+
+
+}

R/CreateInputsCrit.R

-#*************************************************************************************************
-#' Creation of the InputsCrit object required to the ErrorCrit functions.
-#'
-#' Users wanting to use FUN_CRIT functions that are not included in 
-#' the package must create their own InputsCrit object accordingly.
-#*************************************************************************************************
-#' @title  Creation of the InputsCrit object required to the ErrorCrit functions
-#' @author Laurent Coron (June 2014)
-#' @seealso \code{\link{RunModel}}, \code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}}, \code{\link{CreateCalibOptions}}
-#' @example tests/example_ErrorCrit.R
-#' @encoding UTF-8
-#' @export
-#_FunctionInputs__________________________________________________________________________________
-#' @param  FUN_CRIT     [function] error criterion function (e.g. ErrorCrit_RMSE, ErrorCrit_NSE)
-#' @param  InputsModel  [object of class \emph{InputsModel}] see \code{\link{CreateInputsModel}} for details
-#' @param  RunOptions   [object of class \emph{RunOptions}] see \code{\link{CreateRunOptions}} for details
-#' @param  Qobs         [numeric] series of observed discharges [mm]
-#' @param  BoolCrit     (optional) [boolean] boolean giving the time steps to consider in the computation (all time steps are consider by default)
-#' @param  transfo      (optional) [character] name of the transformation (e.g. "", "sqrt", "log", "inv", "sort")
-#' @param  Ind_zeroes   (optional) [numeric] indices of the time-steps where zeroes are observed
-#' @param  epsilon      (optional) [numeric] epsilon to add to all Qobs and Qsim if \emph{$Ind_zeroes} is not empty
-#_FunctionOutputs_________________________________________________________________________________
-#' @return  [list] object of class \emph{InputsCrit} containing the data required to evaluate the model outputs; it can include the following:
-#'          \tabular{ll}{
-#'          \emph{$BoolCrit  }  \tab   [boolean] boolean giving the time steps to consider in the computation \cr
-#'          \emph{$Qobs      }  \tab   [numeric] series of observed discharges [mm] \cr
-#'          \emph{$transfo   }  \tab   [character] name of the transformation (e.g. "", "sqrt", "log", "inv", "sort") \cr
-#'          \emph{$Ind_zeroes}  \tab   [numeric] indices of the time-steps where zeroes are observed \cr
-#'          \emph{$epsilon   }  \tab   [numeric] epsilon to add to all Qobs and Qsim if \emph{$Ind_zeroes} is not empty \cr
-#'          }
-#**************************************************************************************************
-CreateInputsCrit <- function(FUN_CRIT,InputsModel,RunOptions,Qobs,BoolCrit=NULL,transfo="",Ind_zeroes=NULL,epsilon=NULL){
-
-  ObjectClass <- NULL;
-
-  ##check_FUN_CRIT
-    BOOL <- FALSE;
-    if(identical(FUN_CRIT,ErrorCrit_NSE) | identical(FUN_CRIT,ErrorCrit_KGE) | identical(FUN_CRIT,ErrorCrit_KGE2) | 
-       identical(FUN_CRIT,ErrorCrit_RMSE)){
-      BOOL <- TRUE; }
-    if(!BOOL){ stop("incorrect FUN_CRIT for use in CreateInputsCrit \n"); return(NULL); } 
-
-  ##check_arguments
-    if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n" ); return(NULL); } 
-    if(inherits(RunOptions ,"RunOptions" )==FALSE){ stop("RunOptions must be of class 'RunOptions' \n" ); return(NULL); } 
-    LLL <- length(InputsModel$DatesR[RunOptions$IndPeriod_Run])
-
-    if(is.null(Qobs)    ){ stop("Qobs is missing \n"); return(NULL); } 
-    if(!is.vector( Qobs)){ stop(paste("Qobs must be a vector of numeric values \n",sep="")); return(NULL); }
-    if(!is.numeric(Qobs)){ stop(paste("Qobs must be a vector of numeric values \n",sep="")); return(NULL); }
-    if(length(Qobs)!=LLL){ stop("Qobs and InputsModel series must have the same length \n"); return(NULL); } 
-
-    if(is.null(BoolCrit)){ BoolCrit <- rep(TRUE,length(Qobs)); }
-    if(!is.logical(BoolCrit)){ stop("BoolCrit must be a vector of boolean \n" ); return(NULL); } 
-    if(length(BoolCrit)!=LLL){ stop("BoolCrit and InputsModel series must have the same length \n"); return(NULL); } 
-
-    if(is.null(transfo)      ){ stop("transfo must be a chosen among the following: '', 'sqrt', 'log' or 'inv' \n"); return(NULL); } 
-    if(!is.vector(transfo   )){ stop("transfo must be a chosen among the following: '', 'sqrt', 'log' or 'inv' \n"); return(NULL); } 
-    if(length(transfo)!=1    ){ stop("transfo must be a chosen among the following: '', 'sqrt', 'log' or 'inv' \n"); return(NULL); } 
-    if(!is.character(transfo)){ stop("transfo must be a chosen among the following: '', 'sqrt', 'log' or 'inv' \n"); return(NULL); } 
-    if(transfo %in% c("","sqrt","log","inv") == FALSE){
-                                stop("transfo must be a chosen among the following: '', 'sqrt', 'log' or 'inv' \n"); return(NULL); } 
-
-    if(!is.null(Ind_zeroes)){
-      if(!is.vector( Ind_zeroes)){ stop("Ind_zeroes must be a vector of integers \n" ); return(NULL); } 
-      if(!is.integer(Ind_zeroes)){ stop("Ind_zeroes must be a vector of integers \n" ); return(NULL); } 
-    }
-    if(!is.null(epsilon)){
-      if(!is.vector( epsilon) | length(epsilon)!=1 | !is.numeric(epsilon)){ 
-        stop("epsilon must be single numeric value \n" ); return(NULL); }
-      epsilon=as.double(epsilon);
-    }
-
-  ##Create_InputsCrit
-    InputsCrit <- list(BoolCrit=BoolCrit,Qobs=Qobs,transfo=transfo,Ind_zeroes=Ind_zeroes,epsilon=epsilon);
-    class(InputsCrit) <- c("InputsCrit",ObjectClass);
-    return(InputsCrit);
+CreateInputsCrit <- function(FUN_CRIT,
+                             InputsModel,
+                             RunOptions,
+                             Obs,
+                             VarObs = "Q",
+                             BoolCrit = NULL,
+                             transfo = "",
+                             Weights = NULL,
+                             epsilon = NULL,
+                             warnings = TRUE) {
 
 
-}
+  ObjectClass <- NULL
+
+
+  ## ---------- check arguments
+
+  ## 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
+  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)
+    Obs <- list(Obs)
+  } else {
+    idLayer <- lapply(Obs, function(i) {
+      if (is.list(i)) {
+        length(i)
+      } else {
+        1L
+      }
+    })
+    Obs <- lapply(Obs, function(x) rowMeans(as.data.frame(x)))
+  }
+
+
+  ## create list of arguments
+  listArgs <- list(FUN_CRIT   = FUN_CRIT,
+                   Obs        = Obs,
+                   VarObs     = VarObs,
+                   BoolCrit   = BoolCrit,
+                   idLayer    = idLayer,
+                   transfo    = as.character(transfo),
+                   Weights    = Weights,
+                   epsilon    = epsilon)
+
+
+  ## check lists lengths
+  for (iArgs in names(listArgs)) {
+    if (iArgs %in% c("Weights", "BoolCrit", "epsilon")) {
+      if (any(is.null(listArgs[[iArgs]]))) {
+        listArgs[[iArgs]] <- lapply(seq_along(listArgs$FUN_CRIT), function(x) NULL)
+      }
+    }
+    if (iArgs %in% c("FUN_CRIT", "VarObs", "transfo", "Weights") & length(listArgs[[iArgs]]) > 1L) {
+      listArgs[[iArgs]] <- as.list(listArgs[[iArgs]])
+    }
+    if (!is.list(listArgs[[iArgs]])) {
+      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)))
+    # if (warnings) {
+    #   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")
+  }
+  if (any(c("SCA", "SWE") %in% VarObs) & !inherits(RunOptions, "CemaNeige")) {
+    stop("'VarObs' cannot contain SCA or SWE if CemaNeige is not used")
+  }
+  if ("SCA" %in% VarObs & inherits(RunOptions, "CemaNeige") & !"Gratio"   %in% RunOptions$Outputs_Sim) {
+    stop("'Gratio' is missing in 'Outputs_Sim' of 'RunOptions', which is necessary to output SCA with CemaNeige")
+  }
+  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", "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 (!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 (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))
+    }
+    if (!is.logical(iListArgs2$BoolCrit)) {
+      stop("'BoolCrit' must be a (list of) vector(s) of boolean", 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")
+      if (length(idSCA) == 1L) {
+        vecSCA <- iListArgs2$Obs
+      } else {
+        vecSCA <- unlist(iListArgs2$Obs[idSCA])
+      }
+      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)
+      if (length(idQSS) == 1L) {
+        vecQSS <- iListArgs2$Obs
+      } else {
+        vecQSS <- unlist(iListArgs2$Obs[idQSS])
+      }
+      if (all(is.na(vecQSS))) {
+        stop("'Obs' contains only missing values", call. = FALSE)
+      }
+      if (min(vecQSS, na.rm = TRUE) < 0) {
+        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)
+    }
+    isNotInTransfo <- !(iListArgs2$transfo %in% inTransfo)
+    if (any(isNotInTransfo)) {
+      powTransfo <- iListArgs2$transfo[isNotInTransfo]
+      powTransfo <- gsub("\\^|[[:alpha:]]", "", powTransfo)
+      numExpTransfo <- suppressWarnings(as.numeric(powTransfo))
+      if (any(is.na(numExpTransfo))) {
+        stop(msgTransfo, call. = FALSE)
+      }
+      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)) {
+        stop("'epsilon' must be a single (list of) positive value(s)", call. = FALSE)
+      }
+    } 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)"
+      if (identical(iListArgs2$FUN_CRIT, ErrorCrit_KGE)) {
+        warning(sprintf(warn_log_kge, "KGE"), call. = FALSE)
+      }
+      if (identical(iListArgs2$FUN_CRIT, ErrorCrit_KGE2)) {
+        warning(sprintf(warn_log_kge, "KGE'"), call. = FALSE)
+      }
+    }
+
+
+    ## Create InputsCrit
+    iInputsCrit <- list(FUN_CRIT   = iListArgs2$FUN_CRIT,
+                        Obs        = iListArgs2$Obs,
+                        VarObs     = iListArgs2$VarObs,
+                        BoolCrit   = iListArgs2$BoolCrit,
+                        idLayer    = iListArgs2$idLayer,
+                        transfo    = iListArgs2$transfo,
+                        epsilon    = iListArgs2$epsilon,
+                        Weights    = iListArgs2$Weights)
+    class(iInputsCrit) <- c("Single", "InputsCrit", ObjectClass)
+    return(iInputsCrit)
+
+  })
+  names(InputsCrit) <- paste0("IC", seq_along(InputsCrit))
+
+  listVarObs <- sapply(InputsCrit, FUN = "[[", "VarObs")
+  inCnVarObs <- c("SCA", "SWE")
+  if (!"ZLayers" %in% names(InputsModel)) {
+    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 ")))
+    }
+  } else {
+    listGroupLayer0 <- sapply(InputsCrit, FUN = "[[", "idLayer")
+    listGroupLayer <- rep(listVarObs, times = listGroupLayer0)
+    tabGroupLayer  <- as.data.frame(table(listGroupLayer))
+    colnames(tabGroupLayer) <- c("VarObs", "freq")
+    nLayers <- length(InputsModel$ZLayers)
+    for (iInCnVarObs in inCnVarObs) {
+      if (any(listVarObs %in% iInCnVarObs)) {
+        if (tabGroupLayer[tabGroupLayer$VarObs %in% iInCnVarObs, "freq"] != nLayers) {
+          stop(sprintf("'Obs' must contain %i vector(s) about %s", nLayers, iInCnVarObs))
+        }
+      }
+    }
+  }
+
+  ## define idLayer as an index of the layer to use
+  for (iInCnVarObs in unique(listVarObs)) {
+    if (!iInCnVarObs %in% inCnVarObs) {
+      for (i in which(listVarObs == iInCnVarObs)) {
+        InputsCrit[[i]]$idLayer <- NA
+      }
+    } else {
+      aa <- listGroupLayer0[listVarObs == iInCnVarObs]
+      aa <- unname(aa)
+      bb <- cumsum(c(0, aa[-length(aa)]))
+      cc <- lapply(seq_along(aa), function(x) seq_len(aa[x]) + bb[x])
+      k <- 1
+      for (i in which(listVarObs == iInCnVarObs)) {
+        InputsCrit[[i]]$idLayer <- cc[[k]]
+        k <- k + 1
+      }
+    }
+  }
+
+
+  ## if only one criterion --> not a list of InputsCrit but directly an InputsCrit
+  if (length(InputsCrit) < 2) {
+    InputsCrit <- InputsCrit[[1L]]
+    InputsCrit["Weights"] <- list(Weights = NULL)
+  } else {
+    if (any(sapply(listArgs$Weights, is.null))) {
+      for (iListArgs in InputsCrit) {
+        iListArgs$Weights <- NULL
+      }
+      class(InputsCrit) <- c("Multi", "InputsCrit", ObjectClass)
+    } else {
+      class(InputsCrit) <- c("Compo", "InputsCrit", ObjectClass)
+    }
+    combInputsCrit <- combn(x = length(InputsCrit), m = 2)
+    apply(combInputsCrit, MARGIN = 2, function(i) {
+      equalInputsCrit <- identical(InputsCrit[[i[1]]], InputsCrit[[i[2]]])
+      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/CreateInputsModel.R

-#*************************************************************************************************
-#' Creation of the InputsModel object required to the RunModel functions.
-#'
-#' Users wanting to use FUN_MOD functions that are not included in 
-#' the package must create their own InputsModel object accordingly.
-#*************************************************************************************************
-#' @title  Creation of the InputsModel object required to the RunModel functions
-#' @author Laurent Coron (June 2014)
-#' @seealso \code{\link{RunModel}}, \code{\link{CreateRunOptions}}, \code{\link{CreateInputsCrit}}, \code{\link{CreateCalibOptions}}, \code{\link{DataAltiExtrapolation_Valery}}
-#' @example tests/example_RunModel.R
-#' @encoding UTF-8
-#' @export
-#_FunctionInputs__________________________________________________________________________________
-#' @param  FUN_MOD     [function] hydrological model function (e.g. RunModel_GR4J, RunModel_CemaNeigeGR4J)
-#' @param  DatesR      [POSIXlt] vector of dates required to create the GR model and CemaNeige module inputs
-#' @param  Precip      [numeric] time series of total precipitation (catchment average) [mm], required to create the GR model and CemaNeige module inputs
-#' @param  PotEvap     [numeric] time series of potential evapotranspiration (catchment average) [mm], required to create the GR model inputs
-#' @param  TempMean    (optional) [numeric] time series of mean air temperature [degC], required to create the CemaNeige module inputs
-#' @param  TempMin     (optional) [numeric] time series of min air temperature [degC], possibly used to create the CemaNeige module inputs
-#' @param  TempMax     (optional) [numeric] time series of max air temperature [degC], possibly used to create the CemaNeige module inputs
-#' @param  ZInputs     (optional) [numeric] real giving the mean elevation of the Precip and Temp series (before extrapolation) [m]
-#' @param  HypsoData   (optional) [numeric] vector of 101 reals: min, q01 to q99 and max of catchment elevation distribution [m], required to create the GR model inputs, if not defined a single elevation is used for CemaNeige
-#' @param  NLayers     (optional) [numeric] integer giving the number of elevation layers requested [-], required to create the GR model inputs, default=5
-#' @param  quiet       (optional) [boolean] boolean indicating if the function is run in quiet mode or not, default=FALSE
-#_FunctionOutputs_________________________________________________________________________________
-#' @return  [list] object of class \emph{InputsModel} containing the data required to evaluate the model outputs; it can include the following:
-#'          \tabular{ll}{
-#'          \emph{$DatesR              }  \tab   [POSIXlt] vector of dates \cr
-#'          \emph{$Precip              }  \tab   [numeric] time series of total precipitation (catchment average) [mm] \cr
-#'          \emph{$PotEvap             }  \tab   [numeric] time series of potential evapotranspiration (catchment average) [mm], \cr\tab defined if FUN_MOD includes GR4H, GR4J, GR5J, GR6J, GR2M or GR1A \cr \cr
-#'          \emph{$LayerPrecip         }  \tab   [list] list of time series of precipitation (layer average) [mm], \cr\tab defined if FUN_MOD includes CemaNeige \cr \cr
-#'          \emph{$LayerTempMean       }  \tab   [list] list of time series of mean air temperature (layer average) [degC], \cr\tab defined if FUN_MOD includes CemaNeige \cr \cr
-#'          \emph{$LayerFracSolidPrecip}  \tab   [list] list of time series of solid precip. fract. (layer average) [-], \cr\tab defined if FUN_MOD includes CemaNeige \cr \cr
-#'          }
-#**************************************************************************************************
-CreateInputsModel <- function(FUN_MOD,DatesR,Precip,PotEvap=NULL,TempMean=NULL,TempMin=NULL,TempMax=NULL,ZInputs=NULL,HypsoData=NULL,NLayers=5,quiet=FALSE){
-
-  ObjectClass <- NULL;
-
-  ##check_FUN_MOD
-    BOOL <- FALSE;
-    if(identical(FUN_MOD,RunModel_GR4H)){
-      ObjectClass <- c(ObjectClass,"hourly","GR"); 
-      TimeStep <- as.integer(60*60);
-      BOOL <- TRUE; 
-    }
-    if(identical(FUN_MOD,RunModel_GR4J) | identical(FUN_MOD,RunModel_GR5J) | identical(FUN_MOD,RunModel_GR6J)){
-      ObjectClass <- c(ObjectClass,"daily","GR"); 
-      TimeStep <- as.integer(24*60*60);
-      BOOL <- TRUE; 
-    }
-    if(identical(FUN_MOD,RunModel_GR2M)){
-      ObjectClass <- c(ObjectClass,"GR","monthly"); 
-      TimeStep <- as.integer(c(28,29,30,31)*24*60*60);
-      BOOL <- TRUE; 
-    }
-    if(identical(FUN_MOD,RunModel_GR1A)){
-      ObjectClass <- c(ObjectClass,"GR","yearly"); 
-      TimeStep <- as.integer(c(365,366)*24*60*60);
-      BOOL <- TRUE; 
-    }
-    if(identical(FUN_MOD,RunModel_CemaNeige)){
-      ObjectClass <- c(ObjectClass,"daily","CemaNeige");
-      TimeStep <- as.integer(24*60*60);
-      BOOL <- TRUE; 
-    }
-    if(identical(FUN_MOD,RunModel_CemaNeigeGR4J) | identical(FUN_MOD,RunModel_CemaNeigeGR5J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)){
-      ObjectClass <- c(ObjectClass,"daily","GR","CemaNeige");
-      TimeStep <- as.integer(24*60*60);
-      BOOL <- TRUE; 
-    }
-    if(!BOOL){ stop("incorrect FUN_MOD for use in CreateInputsModel \n"); return(NULL); } 
+CreateInputsModel <- function(FUN_MOD,
+                              DatesR,
+                              Precip, PrecipScale = TRUE,
+                              PotEvap = NULL,
+                              TempMean = NULL, TempMin = NULL, TempMax = NULL,
+                              ZInputs = NULL, HypsoData = NULL, NLayers = 5,
+                              Qupstream = NULL, LengthHydro = NULL, BasinAreas = NULL,
+                              QupstrUnit = "mm",
+                              verbose = TRUE) {
+
+
+  ObjectClass <- NULL
+
+  ## check DatesR
+  if (is.null(DatesR)) {
+    stop("'DatesR' is missing")
+  }
+  if (!"POSIXlt" %in% class(DatesR) & !"POSIXct" %in% class(DatesR)) {
+    stop("'DatesR' must be defined as 'POSIXlt' or 'POSIXct'")
+  }
+  if (!"POSIXlt" %in% class(DatesR)) {
+    DatesR <- as.POSIXlt(DatesR)
+  }
+  if (any(duplicated(DatesR))) {
+    stop("'DatesR' must not include duplicated values")
+  }
+  LLL <- length(DatesR)
+
+  ## check FUN_MOD
+  FUN_MOD <- match.fun(FUN_MOD)
+  FeatFUN_MOD <- .GetFeatModel(FUN_MOD = FUN_MOD, DatesR = DatesR)
+  ObjectClass <- FeatFUN_MOD$Class
+  TimeStep <- FeatFUN_MOD$TimeStep
+
 
   ##check_arguments
-    if("GR" %in% ObjectClass | "CemaNeige" %in% ObjectClass){
-      if(is.null(DatesR)){ stop("DatesR is missing \n"); return(NULL); } 
-      if("POSIXlt" %in% class(DatesR) == FALSE & "POSIXct" %in% class(DatesR) == FALSE){ stop("DatesR must be defined as POSIXlt or POSIXct \n"); return(NULL); }
-      if("POSIXlt" %in% class(DatesR) == FALSE){ DatesR <- as.POSIXlt(DatesR); }
-      if(difftime(tail(DatesR,1),tail(DatesR,2),units="secs")[[1]] %in% TimeStep==FALSE){ stop(paste("the time step of the model inputs must be ",TimeStep," seconds \n",sep="")); return(NULL); }    
-      LLL <- length(DatesR);
-    }
-    if("GR" %in% ObjectClass){
-      if(is.null(Precip  )){ stop("Precip is missing \n"  ); return(NULL); } 
-      if(is.null(PotEvap )){ stop("PotEvap is missing \n" ); return(NULL); } 
-      if(!is.vector( Precip) | !is.vector( PotEvap)){ stop("Precip and PotEvap must be vectors of numeric values \n"); return(NULL); } 
-      if(!is.numeric(Precip) | !is.numeric(PotEvap)){ stop("Precip and PotEvap must be vectors of numeric values \n"); return(NULL); } 
-      if(length(Precip)!=LLL | length(PotEvap)!=LLL){ stop("Precip, PotEvap and DatesR must have the same length \n"); return(NULL); } 
-    }
-    if("CemaNeige" %in% ObjectClass){
-      if(is.null(Precip  )){ stop("Precip is missing \n"  ); return(NULL); } 
-      if(is.null(TempMean)){ stop("TempMean is missing \n"); return(NULL); } 
-      if(!is.vector( Precip) | !is.vector( TempMean)){ stop("Precip and TempMean must be vectors of numeric values \n"); return(NULL); } 
-      if(!is.numeric(Precip) | !is.numeric(TempMean)){ stop("Precip and TempMean must be vectors of numeric values \n"); return(NULL); } 
-      if(length(Precip)!=LLL | length(TempMean)!=LLL){ stop("Precip, TempMean and DatesR must have the same length \n"); return(NULL); } 
-      if(is.null(TempMin)!=is.null(TempMax)){ stop("TempMin and TempMax must be both defined if not null \n"); return(NULL); }
-      if(!is.null(TempMin) & !is.null(TempMax)){ 
-        if(!is.vector( TempMin) | !is.vector( TempMax)){ stop("TempMin and TempMax must be vectors of numeric values \n"); return(NULL); } 
-        if(!is.numeric(TempMin) | !is.numeric(TempMax)){ stop("TempMin and TempMax must be vectors of numeric values \n"); return(NULL); } 
-        if(length(TempMin)!=LLL | length(TempMax)!=LLL){ stop("TempMin, TempMax and DatesR must have the same length \n"); return(NULL); }
-      }
-      if(!is.null(HypsoData)){ 
-        if(!is.vector( HypsoData)){ stop("HypsoData must be a vector of numeric values  if not null \n"); return(NULL); } 
-        if(!is.numeric(HypsoData)){ stop("HypsoData must be a vector of numeric values  if not null \n"); return(NULL); } 
-        if(length(HypsoData)!=101){ stop("HypsoData must be of length 101 if not null \n"); return(NULL); } 
-        if(sum(is.na(HypsoData))!=0 & sum(is.na(HypsoData))!=101){ stop("HypsoData must not contain any NA if not null \n"); return(NULL); } 
-      }
-      if(!is.null(ZInputs)){
-        if(length(ZInputs)!=1                   ){ stop("\t ZInputs must be a single numeric value if not null \n"); return(NULL); } 
-        if(is.na(ZInputs) | !is.numeric(ZInputs)){ stop("\t ZInputs must be a single numeric value if not null \n"); return(NULL); } 
-      }
-      if(is.null(HypsoData)){ 
-        if(!quiet){ warning("\t HypsoData is missing => a single layer is used and no extrapolation is made \n"); }
-        HypsoData <- as.numeric(rep(NA,101)); ZInputs <- as.numeric(NA); NLayers <- as.integer(1);
-      }
-      if(is.null(ZInputs)){ 
-        if(!quiet & !identical(HypsoData,as.numeric(rep(NA,101)))){ warning("\t ZInputs is missing => HypsoData[51] is used \n"); }
-        ZInputs <- HypsoData[51];
+
+  if ("GR" %in% ObjectClass) {
+    if (is.null(Precip)) {
+      stop("Precip is missing")
+    }
+    if (is.null(PotEvap)) {
+      stop("'PotEvap' is missing")
+    }
+    if (!is.vector(Precip) | !is.vector(PotEvap)) {
+      stop("'Precip' and 'PotEvap' must be vectors of numeric values")
+    }
+    if (!is.numeric(Precip) | !is.numeric(PotEvap)) {
+      stop("'Precip' and 'PotEvap' must be vectors of numeric values")
+    }
+    if (length(Precip) != LLL | length(PotEvap) != LLL) {
+      stop("'Precip', 'PotEvap' and 'DatesR' must have the same length")
+    }
+  }
+  if ("CemaNeige" %in% ObjectClass) {
+    if (is.null(Precip)) {
+      stop("'Precip' is missing")
+    }
+    if (is.null(TempMean)) {
+      stop("'TempMean' is missing")
+    }
+    if (!is.vector(Precip) | !is.vector(TempMean)) {
+      stop("'Precip' and 'TempMean' must be vectors of numeric values")
+    }
+    if (!is.numeric(Precip) | !is.numeric(TempMean)) {
+      stop("'Precip' and 'TempMean' must be vectors of numeric values")
+    }
+    if (length(Precip) != LLL | length(TempMean) != LLL) {
+      stop("'Precip', 'TempMean' and 'DatesR' must have the same length")
+    }
+    if (is.null(TempMin) != is.null(TempMax)) {
+      stop("'TempMin' and 'TempMax' must be both defined if not null")
+    }
+    if (!is.null(TempMin) & !is.null(TempMax)) {
+      if (!is.vector(TempMin) | !is.vector(TempMax)) {
+        stop("'TempMin' and 'TempMax' must be vectors of numeric values")
       }
+      if (!is.numeric(TempMin) | !is.numeric(TempMax)) {
+        stop("'TempMin' and 'TempMax' must be vectors of numeric values")
+      }
+      if (length(TempMin) != LLL | length(TempMax) != LLL) {
+        stop("'TempMin', 'TempMax' and 'DatesR' must have the same length")
+      }
+    }
+    if (!is.null(HypsoData)) {
+      if (!is.vector(HypsoData)) {
+        stop("'HypsoData' must be a vector of numeric values  if not null")
+      }
+      if (!is.numeric(HypsoData)) {
+        stop("'HypsoData' must be a vector of numeric values  if not null")
+      }
+      if (length(HypsoData) != 101) {
+        stop("'HypsoData' must be of length 101 if not null")
+      }
+      if (sum(is.na(HypsoData)) != 0 & sum(is.na(HypsoData)) != 101) {
+        stop("'HypsoData' must not contain any NA if not null")
+      }
+    }
+    if (!is.null(ZInputs)) {
+      if (length(ZInputs) != 1) {
+        stop("'ZInputs' must be a single numeric value if not null")
+      }
+      if (is.na(ZInputs) | !is.numeric(ZInputs)) {
+        stop("'ZInputs' must be a single numeric value if not null")
+      }
+    }
+    if (is.null(HypsoData)) {
+      if (verbose) {
+        warning("'HypsoData' is missing: a single layer is used and no extrapolation is made", call. = FALSE)
+      }
+      HypsoData <- as.numeric(rep(NA, 101))
+      ZInputs   <- as.numeric(NA)
+      NLayers   <- as.integer(1)
+
+    }
+    if (is.null(ZInputs)) {
+      if (verbose & !identical(HypsoData, as.numeric(rep(NA, 101)))) {
+        warning("'ZInputs' is missing: HypsoData[51] is used")
+      }
+      ZInputs <- HypsoData[51L]
+    }
+    if (NLayers <= 0) {
+      stop("'NLayers' must be a positive integer value")
+    }
+    if (NLayers != as.integer(NLayers)) {
+      warning("Coerce 'NLayers' to be of integer type (",  NLayers, ": ", as.integer(NLayers), ")")
+      NLayers <- as.integer(NLayers)
     }
+  }
 
+  ## check semi-distributed mode
+  if (!is.null(Qupstream) & !is.null(LengthHydro) & !is.null(BasinAreas)) {
+    ObjectClass <- c(ObjectClass, "SD")
+  } else if (verbose & !all(c(is.null(Qupstream), is.null(LengthHydro), is.null(BasinAreas)))) {
+    warning("Missing argument: 'Qupstream', 'LengthHydro' and 'BasinAreas' must all be set to run in a semi-distributed mode. The lumped mode will be used")
+  }
+  if ("SD" %in% ObjectClass) {
+    if (!("daily" %in% ObjectClass) & !("hourly" %in% ObjectClass)) {
+      stop("Only daily and hourly time steps can be used in a semi-distributed mode")
+    }
+    if (!is.matrix(Qupstream) | !is.numeric(Qupstream)) {
+      stop("'Qupstream' must be a matrice of numeric values")
+    }
+    if (!is.vector(LengthHydro) | !is.vector(BasinAreas) | !is.numeric(LengthHydro) | !is.numeric(BasinAreas)) {
+      stop("'LengthHydro' and 'BasinAreas' must be vectors of numeric values")
+    }
+    if (ncol(Qupstream) != length(LengthHydro)) {
+      stop("'Qupstream' number of columns and 'LengthHydro' length must be equal")
+    }
+    if (length(LengthHydro) + 1 != length(BasinAreas)) {
+      stop("'BasinAreas' must have one more element than 'LengthHydro'")
+    }
+    if (nrow(Qupstream) != LLL) {
+      stop("'Qupstream' must have same number of rows as 'DatesR' length")
+    }
+    if (any(is.na(Qupstream))) {
+      warning("'Qupstream' contains NA values: model outputs will contain NAs")
+    }
+    if (any(LengthHydro > 1000)) {
+      warning("The unit of 'LengthHydro' has changed from m to km in airGR >= 1.6.12: values superior to 1000 km seem unrealistic")
+    }
+    QupstrUnit <- tolower(QupstrUnit)
+    QupstrUnit <- match.arg(arg = QupstrUnit, choices = c("mm", "m3", "m3/s", "l/s"))
+  }
 
   ##check_NA_values
-    BOOL_NA <- rep(FALSE,length(DatesR));
-    if("GR" %in% ObjectClass){
-      BOOL_NA_TMP <- (Precip  < 0) | is.na(Precip );  if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(!quiet){ warning("\t Values < 0 or NA values detected in Precip series  \n"); } } 
-      BOOL_NA_TMP <- (PotEvap < 0) | is.na(PotEvap);  if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(!quiet){ warning("\t Values < 0 or NA values detected in PotEvap series \n"); } } 
-    }
-    if("CemaNeige" %in% ObjectClass){
-      BOOL_NA_TMP <- (Precip  < 0    ) | is.na(Precip  );  if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(!quiet){ warning("\t Values < 0 or NA values detected in Precip series       \n"); } } 
-      BOOL_NA_TMP <- (TempMean<(-150)) | is.na(TempMean);  if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(!quiet){ warning("\t Values < -150) or NA values detected in TempMean series \n"); } }
-      if(!is.null(TempMin) & !is.null(TempMax)){
-      BOOL_NA_TMP <- (TempMin<(-150)) | is.na(TempMin);  if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(!quiet){ warning("\t Values < -150) or NA values detected in TempMin series \n"); } }
-      BOOL_NA_TMP <- (TempMax<(-150)) | is.na(TempMax);  if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(!quiet){ warning("\t Values < -150) or NA values detected in TempMax series \n"); } } }
-    }
-    if(sum(BOOL_NA)!=0){
-      WTxt <- NULL;
-      WTxt <- paste(WTxt,"\t Missing values are not allowed in InputsModel \n",sep="");
-      Select <- (max(which(BOOL_NA))+1):length(BOOL_NA);
-      if(Select[1]>Select[2]){ stop(paste("time series could not be trunced since missing values were detected at the list time-step  \n",sep="")); return(NULL); }
-      if("GR" %in% ObjectClass){
-        Precip <- Precip[Select];  PotEvap <- PotEvap[Select]; }
-      if("CemaNeige" %in% ObjectClass){
-        Precip <- Precip[Select];  TempMean <- TempMean[Select]; if(!is.null(TempMin) & !is.null(TempMax)){ TempMin <- TempMin[Select]; TempMax <- TempMax[Select]; } }
-      DatesR <- DatesR[Select];
-      WTxt <- paste(WTxt,"\t -> data were trunced to keep the most recent available time-steps \n",sep="");
-      WTxt <- paste(WTxt,"\t -> ",length(Select)," time-steps were kept \n",sep="");
-      if(!is.null(WTxt) & !quiet){ warning(WTxt); }
+  BOOL_NA <- rep(FALSE, length(DatesR))
+
+  if ("GR" %in% ObjectClass) {
+    BOOL_NA_TMP <- (Precip  < 0) | is.na(Precip)
+    if (sum(BOOL_NA_TMP) != 0) {
+      BOOL_NA <- BOOL_NA |  BOOL_NA_TMP
+      if (verbose) {
+        warning("Values < 0 or NA values detected in 'Precip' series")
+      }
+    }
+    BOOL_NA_TMP <- (PotEvap < 0) | is.na(PotEvap)
+    if (sum(BOOL_NA_TMP) != 0) {
+      BOOL_NA <- BOOL_NA | BOOL_NA_TMP
+      if (verbose) {
+        warning("Values < 0 or NA values detected in 'PotEvap' series")
+      }
+    }
+  }
+  if ("CemaNeige" %in% ObjectClass) {
+    BOOL_NA_TMP <- (Precip  < 0) | is.na(Precip)
+    if (sum(BOOL_NA_TMP) != 0) {
+      BOOL_NA <- BOOL_NA | BOOL_NA_TMP
+      if (verbose) {
+        warning("Values < 0 or NA values detected in 'Precip' series")
+      }
+    }
+    BOOL_NA_TMP <- (TempMean < (-150)) | is.na(TempMean)
+    if (sum(BOOL_NA_TMP) != 0) {
+      BOOL_NA <- BOOL_NA | BOOL_NA_TMP
+      if (verbose) {
+        warning("Values < -150 or NA values detected in 'TempMean' series")
+      }
+    }
+    if (!is.null(TempMin) & !is.null(TempMax)) {
+      BOOL_NA_TMP <- (TempMin < (-150)) | is.na(TempMin)
+      if (sum(BOOL_NA_TMP) != 0) {
+        BOOL_NA <- BOOL_NA | BOOL_NA_TMP
+        if (verbose) {
+          warning("Values < -150 or NA values detected in 'TempMin' series")
+        }
+      }
+      BOOL_NA_TMP <- (TempMax < (-150)) | is.na(TempMax)
+      if (sum(BOOL_NA_TMP) != 0) {
+        BOOL_NA <- BOOL_NA | BOOL_NA_TMP
+        if (verbose) {
+          warning("Values < -150 or NA values detected in 'TempMax' series")
+        }
+      }
+    }
+  }
+  if (sum(BOOL_NA) != 0) {
+    WTxt <- NULL
+    WTxt <- paste(WTxt, "\t Missing values are not allowed in 'InputsModel'", sep = "")
+
+    Select <- (max(which(BOOL_NA)) + 1):length(BOOL_NA)
+
+    if (Select[1L] > Select[2L]) {
+      stop("time series could not be trunced since missing values were detected at the last time-step")
     }
+    if ("GR" %in% ObjectClass) {
+      Precip  <- Precip[Select]
+      PotEvap <- PotEvap[Select]
+    }
+    if ("CemaNeige" %in% ObjectClass) {
+      Precip   <- Precip[Select]
+      TempMean <- TempMean[Select]
+      if (!is.null(TempMin) & !is.null(TempMax)) {
+        TempMin <- TempMin[Select]
+        TempMax <- TempMax[Select]
+      }
+    }
+
+    DatesR <- DatesR[Select]
+
+    WTxt <- paste0(WTxt, "\t -> data were trunced to keep the most recent available time-steps")
+    WTxt <- paste0(WTxt, "\t -> ", length(Select), " time-steps were kept")
+
+    if (!is.null(WTxt) & verbose) {
+      warning(WTxt)
+    }
+  }
 
 
   ##DataAltiExtrapolation_Valery
-    if("CemaNeige" %in% ObjectClass){
-      RESULT <- DataAltiExtrapolation_Valery(DatesR=DatesR,Precip=Precip,TempMean=TempMean,TempMin=TempMin,TempMax=TempMax,ZInputs=ZInputs,HypsoData=HypsoData,NLayers=NLayers,quiet=quiet);
-      if(!quiet){ if(NLayers==1){ cat(paste("\t Input series were successfully created on 1 elevation layer for use by CemaNeige \n",sep=""));
-                         } else { cat(paste("\t Input series were successfully created on ",NLayers," elevation layers for use by CemaNeige \n",sep="")); } }
+  if ("CemaNeige" %in% ObjectClass) {
+    RESULT <- DataAltiExtrapolation_Valery(DatesR = DatesR,
+                                           Precip = Precip, PrecipScale = PrecipScale,
+                                           TempMean = TempMean, TempMin = TempMin, TempMax = TempMax,
+                                           ZInputs = ZInputs, HypsoData = HypsoData, NLayers = NLayers,
+                                           verbose = verbose)
+    if (verbose) {
+      if (NLayers == 1) {
+        message("input series were successfully created on 1 elevation layer for use by CemaNeige")
+      } else {
+        message( "input series were successfully created on ", NLayers, " elevation layers for use by CemaNeige")
+      }
     }
+  }
 
 
   ##Create_InputsModel
-    InputsModel <- list(DatesR=DatesR);
-    if("GR" %in% ObjectClass){
-      InputsModel <- c(InputsModel,list(Precip=as.double(Precip),PotEvap=as.double(PotEvap)));    }
-    if("CemaNeige" %in% ObjectClass){
-      InputsModel <- c(InputsModel,list(LayerPrecip=RESULT$LayerPrecip,LayerTempMean=RESULT$LayerTempMean,
-                                        LayerFracSolidPrecip=RESULT$LayerFracSolidPrecip,ZLayers=RESULT$ZLayers));    }
+  InputsModel <- list(DatesR = DatesR)
+  if ("GR" %in% ObjectClass) {
+    InputsModel <- c(InputsModel, list(Precip = as.double(Precip), PotEvap = as.double(PotEvap)))
+  }
+  if ("CemaNeige" %in% ObjectClass) {
+    InputsModel <- c(InputsModel, list(LayerPrecip          = RESULT$LayerPrecip,
+                                       LayerTempMean        = RESULT$LayerTempMean,
+                                       LayerFracSolidPrecip = RESULT$LayerFracSolidPrecip,
+                                       ZLayers              = RESULT$ZLayers))
+  }
+  if ("SD" %in% ObjectClass) {
+    # Qupstream is internally stored in m3/time step
+    if (QupstrUnit == "mm") {
+      iConvBasins <- which(!is.na(BasinAreas[seq.int(length(LengthHydro))]))
+      Qupstream[, iConvBasins] <- Qupstream[, iConvBasins] * rep(BasinAreas[iConvBasins], each = LLL) * 1e3
+    } else if (QupstrUnit == "m3/s") {
+      Qupstream <- Qupstream * TimeStep
+    } else if (QupstrUnit == "l/s") {
+      Qupstream <- Qupstream * TimeStep / 1e3
+    }
+    InputsModel <- c(InputsModel, list(Qupstream = Qupstream,
+                                       LengthHydro = LengthHydro,
+                                       BasinAreas = BasinAreas))
+  }
 
-    class(InputsModel) <- c("InputsModel",ObjectClass);
-    return(InputsModel);
+  class(InputsModel) <- c("InputsModel", ObjectClass)
 
+  return(InputsModel)
 
-}
 
+
+}

R/CreateRunOptions.R

-CreateRunOptions <- function(FUN_MOD,InputsModel,IndPeriod_WarmUp=NULL,IndPeriod_Run,IniStates=NULL,IniResLevels=NULL,
-                             Outputs_Cal=NULL,Outputs_Sim="all",RunSnowModule=TRUE,MeanAnSolidPrecip=NULL,quiet=FALSE){
+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
+  }
 
-  ObjectClass <- NULL;
+  ## 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
 
-  ##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; 
-    }
-    if(identical(FUN_MOD,RunModel_CemaNeigeGR4J) | identical(FUN_MOD,RunModel_CemaNeigeGR5J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)){
-      ObjectClass <- c(ObjectClass,"GR","CemaNeige","daily");
-      BOOL <- TRUE; 
+  ## 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
     }
-    if(!BOOL){ stop("incorrect FUN_MOD for use in CreateRunOptions \n"); return(NULL); } 
+  }
 
+  ## 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' \n"); return(NULL); } 
-    if("GR" %in% ObjectClass & !inherits(InputsModel,"GR")){
-      stop("InputsModel must be of class 'GR' \n"); return(NULL); } 
-    if("CemaNeige" %in% ObjectClass & !inherits(InputsModel,"CemaNeige")){
-      stop("InputsModel must be of class 'CemaNeige' \n"); return(NULL); } 
-    if("hourly" %in% ObjectClass & !inherits(InputsModel,"hourly")){
-      stop("InputsModel must be of class 'hourly' \n"); return(NULL); } 
-    if("daily" %in% ObjectClass & !inherits(InputsModel,"daily")){
-      stop("InputsModel must be of class 'daily' \n"); return(NULL); } 
-    if("monthly" %in% ObjectClass & !inherits(InputsModel,"monthly")){
-      stop("InputsModel must be of class 'monthly' \n"); return(NULL); } 
-    if("yearly" %in% ObjectClass & !inherits(InputsModel,"yearly")){
-      stop("InputsModel must be of class 'yearly' \n"); return(NULL); } 
+  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 \n"); return(NULL); } 
-    if(!is.numeric(IndPeriod_Run)){ stop("IndPeriod_Run must be a vector of numeric values \n"); return(NULL); } 
-    if(identical(as.integer(IndPeriod_Run),as.integer(seq(from=IndPeriod_Run[1],to=tail(IndPeriod_Run,1),by=1)))==FALSE){
-      stop("IndPeriod_Run must be a continuous sequence of integers \n"); return(NULL); } 
-    if(storage.mode(IndPeriod_Run)!="integer"){ stop("IndPeriod_Run should be of type integer \n"); return(NULL); } 
+  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,"\t Model warm-up period not defined -> default configuration used \n",sep="");
-      ##If_the_run_period_starts_at_the_very_beginning_of_the_time_series
-      if(IndPeriod_Run[1]==as.integer(1)){
-        IndPeriod_WarmUp <- as.integer(0);
-        WTxt <- paste(WTxt,"\t    No data were found for model warm-up! \n",sep="");
+  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 {
-        TmpDateR <- InputsModel$DatesR[IndPeriod_Run[1]] - 365*24*60*60; ### minimal date to start the warmup        
-        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){ 
-        WTxt <- paste(WTxt,"\t    The year preceding the run period is used \n",sep="");
-        } else {
-        WTxt <- paste(WTxt,"\t    Less than a year (without missing values) was found for model warm-up: \n",sep="");
-        WTxt <- paste(WTxt,"\t    (",length(IndPeriod_WarmUp)," time-steps are used for initialisation) \n",sep=""); 
-        }
+        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_Run must be a vector of numeric values \n"); return(NULL); } 
-      if(!is.numeric(IndPeriod_WarmUp)){ stop("IndPeriod_Run must be a vector of numeric values \n"); return(NULL); } 
-      if(storage.mode(IndPeriod_WarmUp)!="integer"){ stop("IndPeriod_Run should be of type integer \n"); return(NULL); } 
-      if(identical(IndPeriod_WarmUp,as.integer(0))){
-        WTxt <- paste(WTxt,"\t No warm-up period is used! \n",sep=""); }
-      if((IndPeriod_Run[1]-1)!=tail(IndPeriod_WarmUp,1) & !identical(IndPeriod_WarmUp,as.integer(0))){ 
-        WTxt <- paste(WTxt,"\t Model warm-up period is not directly before the model run period \n",sep=""); }
-    }
-    if(!is.null(WTxt) & !quiet){ warning(WTxt); }
-
-
-  ##check_IniStates_and_IniResLevels
-    if(is.null(IniStates) & is.null(IniResLevels) & !quiet){ 
-      warning("\t Model states initialisation not defined -> default configuration used \n"); }
-    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 <- 3*24*20 + 7; }
-      if("daily"   %in% ObjectClass){ if (identical(FUN_MOD,RunModel_GR5J)){NState <- 2*20 + 2*NLayers + 7; } else {NState <- 3*20 + 2*NLayers + 7; }}
-      if("monthly" %in% ObjectClass){ NState <- 2; }
-      if("yearly"  %in% ObjectClass){ NState <- 1; }
-    }
-    if(!is.null(IniStates)){
-      if(!is.vector( IniStates)   ){ stop("IniStates must be a vector of numeric values \n"); return(NULL);  }
-      if(!is.numeric(IniStates)   ){ stop("IniStates must be a vector of numeric values \n"); return(NULL);  }
-      if(length(IniStates)!=NState){ stop(paste("The length of IniStates must be ",NState," for the chosen FUN_MOD \n",sep="")); return(NULL);  }
-    } else {
-      IniStates <- as.double(rep(0.0,NState));
+  }
+  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("GR" %in% ObjectClass & ("monthly" %in% ObjectClass | "daily" %in% ObjectClass | "hourly" %in% ObjectClass)){
-      if(!is.null(IniResLevels)){
-        if(!is.vector(IniResLevels) ){ stop("IniResLevels must be a vector of numeric values \n"); return(NULL);  }
-        if(!is.numeric(IniResLevels)){ stop("IniResLevels must be a vector of numeric values \n"); return(NULL);  }
-        if(length(IniResLevels)!=2 ) { stop("The length of IniStates must be 2 for the chosen FUN_MOD \n"); return(NULL);  }
+    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 {
-        IniResLevels <- as.double(c(0.3,0.5));
+        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
+        }
       }
-    } else {
-      if(!is.null(IniResLevels)){ stop("IniResLevels can only be used with monthly or daily or hourly GR models \n") }
+      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 <- NULL;
-      if(identical(FUN_MOD,RunModel_GR4H)){
-        Outputs_all <- c(Outputs_all,"PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim"); }
-      if(identical(FUN_MOD,RunModel_GR4J) | identical(FUN_MOD,RunModel_CemaNeigeGR4J)){
-        Outputs_all <- c(Outputs_all,"PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim"); }
-      if(identical(FUN_MOD,RunModel_GR5J) | identical(FUN_MOD,RunModel_CemaNeigeGR5J)){
-        Outputs_all <- c(Outputs_all,"PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim"); }
-      if(identical(FUN_MOD,RunModel_GR6J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)){
-        Outputs_all <- c(Outputs_all,"PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QR1","Exp","QD","Qsim"); }
-      if(identical(FUN_MOD,RunModel_GR2M)){
-        Outputs_all <- c(Outputs_all,"PotEvap","Precip","AE","Perc","P3","Exch","Prod","Rout","Qsim"); }
-      if(identical(FUN_MOD,RunModel_GR1A)){
-        Outputs_all <- c(Outputs_all,"PotEvap","Precip","Qsim"); }
-      if("CemaNeige" %in% ObjectClass){
-        Outputs_all <- c(Outputs_all,"Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt"); }
-
-    ##check_Outputs_Sim
-      if(!is.vector(   Outputs_Sim)){ stop("Outputs_Sim must be a vector of characters \n"); return(NULL);  }
-      if(!is.character(Outputs_Sim)){ stop("Outputs_Sim must be a vector of characters \n"); return(NULL);  }
-      if(sum(is.na(Outputs_Sim))!=0){ stop("Outputs_Sim must not contain NA \n"); return(NULL);  }
-      if("all" %in% Outputs_Sim){ Outputs_Sim <- c("DatesR",Outputs_all,"StateEnd");  }
-      Test <- which(Outputs_Sim %in% c("DatesR",Outputs_all,"StateEnd") == FALSE); if(length(Test)!=0){ 
-        stop(paste("Outputs_Sim is incorrectly defined: ",paste(Outputs_Sim[Test],collapse=", ")," not found \n",sep="")); return(NULL); } 
-      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("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 \n"); return(NULL);  }
-        if(!is.character(Outputs_Cal)){ stop("Outputs_Cal must be a vector of characters \n"); return(NULL);  }
-        if(sum(is.na(Outputs_Cal))!=0){ stop("Outputs_Cal must not contain NA \n"); return(NULL);  }
-      }
-      if("all" %in% Outputs_Cal){ Outputs_Cal <- c("DatesR",Outputs_all,"StateEnd");  }
-      Test <- which(Outputs_Cal %in% c("DatesR",Outputs_all,"StateEnd") == FALSE); if(length(Test)!=0){ 
-        stop(paste("Outputs_Cal is incorrectly defined: ",paste(Outputs_Cal[Test],collapse=", ")," not found \n",sep="")); return(NULL); } 
-      Outputs_Cal <- Outputs_Cal[!duplicated(Outputs_Cal)];
+  ##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_RunSnowModule
-    if("CemaNeige" %in% ObjectClass){
-      if(!is.vector( RunSnowModule)){ stop("RunSnowModule must be a single boolean \n"); return(NULL);  }
-      if(!is.logical(RunSnowModule)){ stop("RunSnowModule must be either TRUE or FALSE \n"); return(NULL);  }
-      if(length(RunSnowModule)!=1  ){ stop("RunSnowModule must be either TRUE or FALSE \n"); return(NULL);  }
+  ##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' \n"); return(NULL);  }
-      MeanAnSolidPrecip <- rep(mean(SolidPrecip)*Factor,NLayers); ### default value: same Gseuil for all layers
-      if(!quiet){ warning("\t MeanAnSolidPrecip not defined -> it was automatically set to c(",paste(round(MeanAnSolidPrecip),collapse=","),") \n"); }
-    }
-    if("CemaNeige" %in% ObjectClass & !is.null(MeanAnSolidPrecip)){
-      if(!is.vector( MeanAnSolidPrecip)    ){ stop(paste("MeanAnSolidPrecip must be a vector of numeric values \n",sep="")); return(NULL);  }
-      if(!is.numeric(MeanAnSolidPrecip)    ){ stop(paste("MeanAnSolidPrecip must be a vector of numeric values \n",sep="")); return(NULL);  }
-      if(length(MeanAnSolidPrecip)!=NLayers){ stop(paste("MeanAnSolidPrecip must be a numeric vector of length ",NLayers," \n",sep="")); return(NULL);  }
+  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(RunSnowModule & "GR" %in% ObjectClass & "CemaNeige" %in% ObjectClass){
-      if("PliqAndMelt" %in% Outputs_Cal == FALSE & "all" %in% Outputs_Cal == FALSE){
-        WTxt <- NULL;
-        WTxt <- paste(WTxt,"\t PliqAndMelt was not defined in Outputs_Cal but is needed to feed the hydrological model with the snow module outputs \n",sep="");
-        WTxt <- paste(WTxt,"\t -> it was automatically added \n",sep="");
-        if(!is.null(WTxt) & !quiet){ warning(WTxt); }
-        Outputs_Cal <- c(Outputs_Cal,"PliqAndMelt"); }
-      if("PliqAndMelt" %in% Outputs_Sim == FALSE & "all" %in% Outputs_Sim == FALSE){
-        WTxt <- NULL;
-        WTxt <- paste(WTxt,"\t PliqAndMelt was not defined in Outputs_Sim but is needed to feed the hydrological model with the snow module outputs \n",sep="");
-        WTxt <- paste(WTxt,"\t -> it was automatically added \n",sep="");
-        if(!is.null(WTxt) & !quiet){ warning(WTxt); }
-        Outputs_Sim <- c(Outputs_Sim,"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);
-    if("CemaNeige" %in% ObjectClass){
-        RunOptions <- c(RunOptions,list(RunSnowModule=RunSnowModule,MeanAnSolidPrecip=MeanAnSolidPrecip));    }
-    class(RunOptions) <- c("RunOptions",ObjectClass);
-    return(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)
 
 
 }