test

INDEX

-BasinInfo               Data sample: characteristics of a fictional
-                        catchment (L0123001, L0123002 or L0123003)
-BasinObs                Data sample: time series of observations of a
-                        fictional catchment (L0123001, L0123002 or
-                        L0123003)
-Calibration             Calibration algorithm which minimises an error
-                        criterion on the model outputs using the
-                        provided functions
-Calibration_HBAN        Calibration algorithm which minimises the error
-                        criterion using the Irstea-HBAN procedure
-Calibration_optim       Calibration algorithm which minimises the error
-                        criterion using the stats::optim function
-CreateCalibOptions      Creation of the CalibOptions object required to
-                        the Calibration functions
-CreateInputsCrit        Creation of the InputsCrit object required to
-                        the ErrorCrit functions
-CreateInputsModel       Creation of the InputsModel object required to
-                        the RunModel functions
-CreateRunOptions        Creation of the RunOptions object required to
-                        the RunModel functions
-DataAltiExtrapolation_HBAN
-                        Altitudinal extrapolation of precipitation and
-                        temperature series
-ErrorCrit               Error criterion using the provided function
-ErrorCrit_KGE           Error criterion based on the KGE formula
-ErrorCrit_KGE2          Error criterion based on the KGE' formula
-ErrorCrit_NSE           Error criterion based on the NSE formula
-ErrorCrit_RMSE          Error criterion based on the RMSE
-PEdaily_Oudin           Computation of daily series of potential
-                        evapotranspiration with Oudin's formula
-RunModel                Run with the provided hydrological model
-                        function
-RunModel_CemaNeige      Run with the CemaNeige snow module
-RunModel_CemaNeigeGR4J
-                        Run with the CemaNeigeGR4J hydrological model
-RunModel_CemaNeigeGR5J
-                        Run with the CemaNeigeGR5J hydrological model
-RunModel_CemaNeigeGR6J
-                        Run with the CemaNeigeGR6J hydrological model
-RunModel_GR4J           Run with the GR4J hydrological model
-RunModel_GR5J           Run with the GR5J hydrological model
-RunModel_GR6J           Run with the GR6J hydrological model
-TransfoParam            Transformation of the parameters using the
-                        provided function
-TransfoParam_CemaNeige
-                        Transformation of the parameters from the
-                        CemaNeige module
-TransfoParam_GR4J       Transformation of the parameters from the GR4J
-                        model
-TransfoParam_GR5J       Transformation of the parameters from the GR5J
-                        model
-TransfoParam_GR6J       Transformation of the parameters from the GR6J
-                        model
-airGR                   Modelling tools used at Irstea-HBAN (France),
-                        including GR4J, GR5J, GR6J and CemaNeige
-plot_OutputsModel       Default preview of model outputs

MD5

-dc740898c129f840d12d222a2ead106d *DESCRIPTION
-eacf8601006227c3c02bcfaf97cbdc0e *INDEX
-f06c85e4eb89b267cfb165274067c4e8 *Meta/Rd.rds
-32a1c5de93e3b6254dbd86b07ba073ba *Meta/data.rds
-e5a11fd9f38a3bf1d0dadd7840a739f7 *Meta/hsearch.rds
-a1d82a0c2244a09e38104c219298794a *Meta/links.rds
-3b9ab8f86ffaa46a406cb5352028ef27 *Meta/nsInfo.rds
-7ecce7e1cc77345bfed1ac9c9b04e61d *Meta/package.rds
-52ca795872157b1a3e2b6f6bfbc480e0 *NAMESPACE
-ebf0fc819595d631b8bf280c4b049940 *R/airGR
-fca9fb51c6dd9ac775dbba00d895d0a7 *R/airGR.rdb
-9362eca21ff55ff6e7866a8653c5371e *R/airGR.rdx
-63a6f712183a364edfac4df460e83c4b *data/L0123001.rda
-aa6993b50e8a59fb4b37c8eb92df6760 *data/L0123002.rda
-cd27feeac8c19aad0c511b33010ea9ed *data/L0123003.rda
-c067204cb6e5b9db41af38aa96b519b0 *help/AnIndex
-a55d619bb922901e93a88bc0cd44fa4a *help/airGR.rdb
-efdee1c7837ed1b695a529e17fa09725 *help/airGR.rdx
-228f1599718d0523f25967469dcfbb61 *help/aliases.rds
-d312909d7d245fe60010b0b763b00f13 *help/paths.rds
-592202faa3a70495a0bf299af2596ac4 *html/00Index.html
-444535b9cb76ddff1bab1e1865a3fb14 *html/R.css
-0b6e15e196f42c1f59b5b95fc0fad2b9 *libs/i386/airGR.dll

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_HBAN}}, \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_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  quiet               (optional) [boolean] boolean indicating if the function is run in quiet mode or not, default=FALSE
+#_FunctionOutputs_________________________________________________________________________________
+#' @return  [list] see \code{\link{Calibration_HBAN}} or \code{\link{Calibration_optim}} 
+#**************************************************************************************************
+Calibration <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_CALIB=Calibration_HBAN,FUN_TRANSFO=NULL,quiet=FALSE){
+    return( FUN_CALIB(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO,quiet=quiet) )
+}
+

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_HBAN}},
+#'          \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_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_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{RunModel}}, \code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}}, \code{\link{CreateInputsCrit}}
+#' @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  RunOptions          (optional) [object of class \emph{RunOptions}] see \code{\link{CreateRunOptions}} for details
+#' @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 (e.g. Calibration_optim)
+#'                             \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)
+#'                             \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,RunOptions=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_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_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_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_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(identical(FUN_MOD,RunModel_GR4J) | identical(FUN_MOD,RunModel_GR5J) | identical(FUN_MOD,RunModel_GR6J) | identical(FUN_MOD,RunModel_CemaNeige)){
+        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);  }
+
+  ##check_RunOptions
+    if(!is.null(RunOptions)){ 
+      if(inherits(RunOptions,"RunOptions")==FALSE){ stop("RunOptions must be of class 'RunOptions' if not null= \n"); return(NULL); }  
+    }
+
+  ##NParam
+    if("GR4J"          %in% ObjectClass){ NParam <- 4; }
+    if("GR5J"          %in% ObjectClass){ NParam <- 5; }
+    if("GR6J"          %in% ObjectClass){ NParam <- 6; }
+    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(!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("GR4J"%in% ObjectClass){ 
+        ParamT <- matrix( c( +3.60, -2.00, +3.40, -9.10,
+                             +3.90, -0.90, +4.10, -8.70,
+                             +4.50, -0.10, +5.00, -8.10),ncol=NParam,byrow=TRUE); }
+      if("GR5J"%in% ObjectClass){ 
+        ParamT <- matrix( c( +3.60, -1.70, +3.30, -9.10, -0.70,
+                             +3.90, -0.60, +4.10, -8.70, +0.30,
+                             +4.50, -0.10, +5.00, -8.10, +0.50),ncol=NParam,byrow=TRUE); }
+      if("GR6J"%in% ObjectClass){ 
+        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=NParam,byrow=TRUE); }
+      if("CemaNeige"%in% ObjectClass){ 
+        ParamT <- matrix( c( -6.26, +0.55,
+                             -2.13, +0.92,
+                             +4.86, +1.40),ncol=NParam,byrow=TRUE); }
+      if("CemaNeigeGR4J"%in% ObjectClass){ 
+        ParamT <- matrix( c( +3.60, -2.00, +3.40, -9.10, -6.26, +0.55,
+                             +3.90, -0.90, +4.10, -8.70, -2.13, +0.92,
+                             +4.50, -0.10, +5.00, -8.10, +4.86, +1.40),ncol=NParam,byrow=TRUE); }
+      if("CemaNeigeGR5J"%in% ObjectClass){ 
+        ParamT <- matrix( c( +3.60, -1.70, +3.30, -9.10, -0.70, -6.26, +0.55,
+                             +3.90, -0.60, +4.10, -8.70, +0.30, -2.13, +0.92,
+                             +4.50, -0.10, +5.00, -8.10, +0.50, +4.86, +1.40),ncol=NParam,byrow=TRUE); }
+      if("CemaNeigeGR6J"%in% ObjectClass){ 
+        ParamT <- matrix( c( +3.60, -1.00, +3.30, -9.10, -0.90, +3.00, -6.26, +0.55,
+                             +3.90, -0.50, +4.10, -8.70, +0.10, +4.00, -2.13, +0.92,
+                             +4.50, +0.50, +5.00, -8.10, +1.10, +5.00, +4.86, +1.40),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);  }
+    }
+
+
+  ##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);
+
+
+}
+
+
+

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_RMSE) | identical(FUN_CRIT,ErrorCrit_NSE) | identical(FUN_CRIT,ErrorCrit_KGE) | identical(FUN_CRIT,ErrorCrit_KGE2)){
+      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);
+
+
+}
+

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_HBAN}}
+#' @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 daily total precipitation (catchment average) [mm], required to create the GR model and CemaNeige module inputs
+#' @param  PotEvap     [numeric] time series of daily potential evapotranspiration (catchment average) [mm], required to create the GR model inputs
+#' @param  TempMean    [numeric] time series of daily mean air temperature [degC], required to create the CemaNeige module inputs
+#' @param  TempMin     (optional) [numeric] time series of daily min air temperature [degC], possibly used to create the CemaNeige module inputs
+#' @param  TempMax     (optional) [numeric] time series of daily 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 daily total precipitation (catchment average) [mm] \cr
+#'          \emph{$PotEvap             }  \tab   [numeric] time series of daily potential evapotranspiration (catchment average) [mm], \cr\tab defined if FUN_MOD includes GR4J, GR5J or GR6J \cr \cr
+#'          \emph{$LayerPrecip         }  \tab   [list] list of time series of daily precipitation (layer average) [mm], \cr\tab defined if FUN_MOD includes CemaNeige \cr \cr
+#'          \emph{$LayerTempMean       }  \tab   [list] list of time series of daily 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 daily 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_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_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); } 
+
+  ##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]]!=TimeStep){ 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];
+      }
+    }
+
+
+  ##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]; } }
+      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); }
+    }
+
+
+  ##DataAltiExtrapolation_HBAN
+    if("CemaNeige" %in% ObjectClass){
+      RESULT <- DataAltiExtrapolation_HBAN(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="")); } }
+    }
+
+
+  ##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));    }
+
+    class(InputsModel) <- c("InputsModel",ObjectClass);
+    return(InputsModel);
+
+
+}
+

R/CreateRunOptions.R

+#*************************************************************************************************
+#' Creation of the RunOptions object required to the RunModel functions.
+#'
+#' Users wanting to use FUN_MOD functions that are not included in 
+#' the package must create their own RunOptions object accordingly.
+#' 
+#' ##### Initialisation options #####
+#' 
+#' The model initialisation options can either be set to a default configuration or be defined by the user.
+#' 
+#' This is done via three vectors: \cr \emph{IndPeriod_WarmUp}, \emph{IniStates}, \emph{IniResLevels}. \cr
+#' A default configuration is used for initialisation if these vectors are not defined.
+#' 
+#' (1) Default initialisation options:
+#' 
+#' \itemize{
+#' \item \emph{IndPeriod_WarmUp} default setting ensures a one-year warm-up using the time-steps preceding the \emph{IndPeriod_Run}. 
+#' The actual length of this warm-up might be shorter depending on data availability (no missing value being allowed on model input series).
+#' 
+#' \item \emph{IniStates} and \emph{IniResLevels} are automatically set to initialise all the model states at 0, except for the production and routing stores which are initialised at 50\% of their capacity. This initialisation is made at the very beginning of the model call (i.e. at the beginning of \emph{IndPeriod_WarmUp} or at the beginning of IndPeriod_Run if the warm-up period is disabled).
+#' }
+#' 
+#' (2) Customisation of initialisation options:
+#' 
+#' \itemize{
+#' \item \emph{IndPeriod_WarmUp} can be used to specify the indices of the warm-up period (within the time-series prepared in InputsModel). \cr
+#' - remark 1:	for most common cases, indices corresponding to one or several years preceding \emph{IndPeriod_Run} are used (e.g. \emph{IndPeriod_WarmUp <- 1000:1365} and \emph{IndPeriod_Run <- 1366:5000)}. \cr
+#' However, it is also possible to perform a long-term initialisation if other indices than the warm-up ones are set in \emph{IndPeriod_WarmUp} (e.g. \emph{IndPeriod_WarmUp <- c( 1:5000 , 1:5000 , 1:5000 ,1000:1365 )}). \cr
+#' - remark 2:	it is also possible to completely disable the warm-up period when using \emph{IndPeriod_WarmUp <- 0}.
+#' 
+#' \item \emph{IniStates} and \emph{IniResLevels} can be used to specify the initial model states. \cr
+#' - remark 1:	if \emph{IniStates} is used, all model states must be provided (e.g. 60 floats [mm] are required for GR4J, GR5J and GR6J; 60+2*NLayers floats [mm] are required for CemaNeigeGR4J, CemaNeigeGR5J and CemaNeigeGR6J; see fortran source code for details). \cr
+#' - remark 2:	in addition to \emph{IniStates}, \emph{IniResLevels} allows to set the filling rate of the production and routing stores for the GR models. For instance for GR4J, GR5J and GR6J: \emph{IniResLevels <- c(0.3,0.5)} should be used to obtain initial fillings of 30\% and 50\% for the production and routing stores, respectively.  \emph{IniResLevels} is optional and can only be used if \emph{IniStates} is also defined (the state values corresponding to these two stores in \emph{IniStates} are not used in such case). \cr \cr
+#' }
+#*************************************************************************************************
+#' @title  Creation of the RunOptions object required to the RunModel functions
+#' @author Laurent Coron (June 2014)
+#' @seealso \code{\link{RunModel}}, \code{\link{CreateInputsModel}}, \code{\link{CreateInputsCrit}}, \code{\link{CreateCalibOptions}}
+#' @example tests/example_RunModel.R
+#' @encoding UTF-8
+#' @export
+#_FunctionInputs__________________________________________________________________________________
+#' @param  FUN_MOD             [function] hydrological model function (e.g. RunModel_GR4J, RunModel_CemaNeigeGR4J)
+#' @param  InputsModel         [object of class \emph{InputsModel}] see \code{\link{CreateInputsModel}} for details
+#' @param  IndPeriod_WarmUp    (optional) [numeric] index of period to be used for the model warm-up [-]
+#' @param  IndPeriod_Run       [numeric] index of period to be used for the model run [-]
+#' @param  IniStates           (optional) [numeric] vector of initial model states [mm]
+#' @param  IniResLevels        (optional) [numeric] vector of initial filling rates for production and routing stores (2 values between 0 and 1) [-]
+#' @param  Outputs_Cal         (optional) [character] vector giving the outputs needed for the calibration \cr (e.g. c("Qsim")), the least outputs the fastest the calibration
+#' @param  Outputs_Sim         (optional) [character] vector giving the requested outputs \cr (e.g. c("DatesR","Qsim","SnowPack")), default="all"
+#' @param  RunSnowModule       (optional) [boolean] option indicating whether CemaNeige should be activated, default=TRUE
+#' @param  MeanAnSolidPrecip   (optional) [numeric] vector giving the annual mean of average solid precipitation for each layer (computed from InputsModel if not defined) [mm/y]
+#' @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{RunOptions} containing the data required to evaluate the model outputs; it can include the following:
+#'          \tabular{ll}{
+#'          \emph{IndPeriod_WarmUp }  \tab  [numeric] index of period to be used for the model warm-up [-] \cr
+#'          \emph{IndPeriod_Run    }  \tab  [numeric] index of period to be used for the model run [-] \cr
+#'          \emph{IniStates        }  \tab  [numeric] vector of initial model states [mm] \cr
+#'          \emph{IniResLevels     }  \tab  [numeric] vector of initial filling rates for production and routing stores [-] \cr
+#'          \emph{Outputs_Cal      }  \tab  [character] character vector giving only the outputs needed for the calibration \cr
+#'          \emph{Outputs_Sim      }  \tab  [character] character vector giving the requested outputs \cr
+#'          \emph{RunSnowModule    }  \tab  [boolean] option indicating whether CemaNeige should be activated \cr
+#'          \emph{MeanAnSolidPrecip}  \tab  [numeric] vector giving the annual mean of average solid precipitation for each layer [mm/y] \cr
+#'          }
+#**************************************************************************************************'
+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){
+
+
+  ObjectClass <- NULL;
+
+
+  ##check_FUN_MOD
+    BOOL <- FALSE;
+    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_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; 
+    }
+    if(!BOOL){ stop("incorrect FUN_MOD for use in CreateRunOptions \n"); return(NULL); } 
+
+
+  ##check_InputsModel
+    if(!inherits(InputsModel,"InputsModel")){
+      stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); } 
+    if("daily" %in% ObjectClass & !inherits(InputsModel,"daily")){
+      stop("InputsModel must be of class 'daily' \n"); return(NULL); } 
+
+
+  ##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); } 
+
+
+  ##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="");
+      ##We_look_for_the_longest_period_preceeding_the_run_period_with_a_maximum_of_one_year
+      } 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("daily"   %in% ObjectClass){ TimeStep <- as.integer(       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    Only ",length(IndPeriod_WarmUp)," time-steps are used!                 \n",sep=""); 
+        }
+      }
+    }
+    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)){ 
+        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("GR" %in% ObjectClass){ 
+      if("daily"   %in% ObjectClass){ NH <-    20; }
+    } else {
+      NH <- 0;
+    }
+    if("CemaNeige" %in% ObjectClass){ NLayers <- length(InputsModel$LayerPrecip); } else { NLayers <- 0; }
+    NState <- 3*NH + 2*NLayers;
+    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(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);  }
+    } else {
+      if("GR" %in% ObjectClass){ IniResLevels <- as.double(c(0.3,0.5)); }
+    }
+
+
+  ##check_Outputs_Cal_and_Sim
+
+    ##Outputs_all
+      Outputs_all <- NULL;
+      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("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)];
+
+
+  ##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_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);  }
+    }
+
+
+  ##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"); }
+    }
+
+
+  ##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);
+
+
+}
+

R/ErrorCrit.R

+#*****************************************************************************************************************
+#' Function which computes an error criterion with the provided function.
+#*****************************************************************************************************************
+#' @title  Error criterion using the provided function
+#' @author Laurent Coron (June 2014)
+#' @seealso \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}}, \code{\link{ErrorCrit_KGE}}
+#' @example tests/example_ErrorCrit.R
+#' @useDynLib airgr
+#' @encoding UTF-8
+#' @export
+#_FunctionInputs__________________________________________________________________________________________________
+#' @param  InputsCrit      [object of class \emph{InputsCrit}] see \code{\link{CreateInputsCrit}} for details
+#' @param  OutputsModel    [object of class \emph{OutputsModel}] see \code{\link{RunModel_GR4J}} or \code{\link{RunModel_CemaNeigeGR4J}} for details
+#' @param  FUN_CRIT        [function] error criterion function (e.g. ErrorCrit_RMSE, ErrorCrit_NSE)
+#' @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, see \code{\link{ErrorCrit_RMSE}} or \code{\link{ErrorCrit_NSE}} for details
+#*****************************************************************************************************************'
+ErrorCrit <- function(InputsCrit,OutputsModel,FUN_CRIT,quiet=FALSE){
+    return( FUN_CRIT(InputsCrit,OutputsModel,quiet=quiet) )
+}
+

R/ErrorCrit_KGE.R

+#*****************************************************************************************************************
+#' Function which computes an error criterion based on the KGE formula proposed by Gupta et al. (2009).
+#'
+#' In addition to the criterion value, the function outputs include a multiplier (-1 or +1) which allows 
+#' the use of the function for model calibration: the product CritValue*Multiplier is the criterion to be minimised 
+#' (e.g. Multiplier=+1 for RMSE, Multiplier=-1 for NSE).
+#*****************************************************************************************************************
+#' @title  Error criterion based on the KGE formula
+#' @author Laurent Coron (June 2014)
+#' @references
+#'   Gupta, H. V., Kling, H., Yilmaz, K. K. and Martinez, G. F. (2009), 
+#'       Decomposition of the mean squared error and NSE performance criteria: Implications
+#'       for improving hydrological modelling, Journal of Hydrology, 377(1-2), 80-91, doi:10.1016/j.jhydrol.2009.08.003. \cr
+#' @seealso \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}}, \code{\link{ErrorCrit_KGE2}}
+#' @examples ## see example of the ErrorCrit function
+#' @encoding UTF-8
+#' @export
+#_FunctionInputs__________________________________________________________________________________________________
+#' @param  InputsCrit      [object of class \emph{InputsCrit}] see \code{\link{CreateInputsCrit}} for details
+#' @param  OutputsModel    [object of class \emph{OutputsModel}] see \code{\link{RunModel_GR4J}} or \code{\link{RunModel_CemaNeigeGR4J}} for details
+#' @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{$CritValue      }   \tab   [numeric] value of the criterion \cr
+#'          \emph{$CritName       }   \tab   [character] name of the criterion \cr
+#'          \emph{$SubCritValues  }   \tab   [numeric] values of the sub-criteria \cr
+#'          \emph{$SubCritNames   }   \tab   [character] names of the sub-criteria \cr
+#'          \emph{$CritBestValue  }   \tab   [numeric] theoretical best criterion value \cr
+#'          \emph{$Multiplier     }   \tab   [numeric] integer indicating whether the criterion is indeed an error (+1) or an efficiency (-1) \cr
+#'          \emph{$Ind_notcomputed}   \tab   [numeric] indices of the time-steps where InputsCrit$BoolCrit=FALSE or no data is available \cr
+#'          }
+#*****************************************************************************************************************
+ErrorCrit_KGE <- function(InputsCrit,OutputsModel,quiet=FALSE){
+
+
+##Arguments_check________________________________
+  if(inherits(InputsCrit,"InputsCrit")==FALSE){ stop("InputsCrit must be of class 'InputsCrit' \n"); return(NULL); }  
+  if(inherits(OutputsModel,"OutputsModel")==FALSE){ stop("OutputsModel must be of class 'OutputsModel' \n"); return(NULL); }  
+
+
+##Initialisation_________________________________
+  CritName <- NA;
+  if(InputsCrit$transfo==""    ){ CritName <- "KGE[Q]"      ; }
+  if(InputsCrit$transfo=="sqrt"){ CritName <- "KGE[sqrt(Q)]"; }
+  if(InputsCrit$transfo=="log" ){ CritName <- "KGE[log(Q)]" ; }
+  if(InputsCrit$transfo=="inv" ){ CritName <- "KGE[1/Q]"    ; }
+  if(InputsCrit$transfo=="sort"){ CritName <- "KGE[sort(Q)]"; }
+  CritValue       <- NA;
+  CritBestValue   <- +1;
+  Multiplier      <- -1; ### must be equal to -1 or +1 only
+
+
+##Data_preparation_______________________________
+  VarObs <- InputsCrit$Qobs  ; VarObs[!InputsCrit$BoolCrit] <- NA; 
+  VarSim <- OutputsModel$Qsim; VarSim[!InputsCrit$BoolCrit] <- NA;  
+  ##Data_transformation
+  if("Ind_zeroes" %in% names(InputsCrit) & "epsilon" %in% names(InputsCrit)){ if(length(InputsCrit$Ind_zeroes)>0){
+    VarObs <- VarObs + InputsCrit$epsilon;
+    VarSim <- VarSim + InputsCrit$epsilon;
+  } }
+  if(InputsCrit$transfo=="sqrt"){ VarObs <- sqrt(VarObs); VarSim <- sqrt(VarSim); }
+  if(InputsCrit$transfo=="log" ){ VarObs <- log(VarObs) ; VarSim <- log(VarSim) ; VarSim[VarSim      < -1E100] <- NA; }
+  if(InputsCrit$transfo=="inv" ){ VarObs <- 1/VarObs    ; VarSim <- 1/VarSim    ; VarSim[abs(VarSim) > 1E+100] <- NA; }
+  if(InputsCrit$transfo=="sort"){ VarObs <- sort(VarObs); VarSim <- sort(VarSim); }
+  ##TS_ignore
+  TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit ;
+  if(sum(!TS_ignore)==0){ OutputsCrit <- list(NA); names(OutputsCrit) <- c("CritValue"); return(OutputsCrit); }
+  if(sum(!TS_ignore)<365 & !quiet){ warning("\t criterion computed on less than 365 time-steps \n"); }  
+  ##Other_variables_preparation
+  meanVarObs <- mean(VarObs[!TS_ignore]);
+  meanVarSim <- mean(VarSim[!TS_ignore]);
+  iCrit           <- 0;
+  SubCritNames    <- NULL;
+  SubCritValues   <- NULL;
+
+
+
+##SubErrorCrit_____KGE_rPearson__________________
+  iCrit <- iCrit+1;
+  SubCritNames[iCrit]  <- paste(CritName," rPEARSON(sim vs. obs)",sep="");
+  SubCritValues[iCrit] <- NA;
+  Numer <- sum( (VarObs[!TS_ignore]-meanVarObs)*(VarSim[!TS_ignore]-meanVarSim) );
+  Deno1 <- sqrt( sum((VarObs[!TS_ignore]-meanVarObs)^2) );
+  Deno2 <- sqrt( sum((VarSim[!TS_ignore]-meanVarSim)^2) );
+  if(Numer==0){ if(Deno1==0 & Deno2==0){ Crit <- 1; } else { Crit <- 0; } 
+  } else { Crit  <- Numer/(Deno1*Deno2); }
+  if(is.numeric(Crit) & is.finite(Crit)){ SubCritValues[iCrit] <- Crit; }
+
+
+##SubErrorCrit_____KGE_alpha_____________________
+  iCrit <- iCrit+1;
+  SubCritNames[iCrit]  <-  paste(CritName," STDEVsim/STDEVobs",sep="");
+  SubCritValues[iCrit] <- NA;
+  Numer <- sd(VarSim[!TS_ignore]); 
+  Denom <- sd(VarObs[!TS_ignore]); 
+  if(Numer==0 & Denom==0){ Crit <- 1; } else { Crit <- Numer/Denom ; }
+  if(is.numeric(Crit) & is.finite(Crit)){ SubCritValues[iCrit] <- Crit; }
+
+
+##SubErrorCrit_____KGE_beta______________________  
+  iCrit <- iCrit+1;
+  SubCritNames[iCrit]  <-  paste(CritName," MEANsim/MEANobs",sep="");
+  SubCritValues[iCrit] <- NA;
+  if(meanVarSim==0 & meanVarObs==0){ Crit <- 1; } else { Crit <- meanVarSim/meanVarObs ; }
+  if(is.numeric(Crit) & is.finite(Crit)){ SubCritValues[iCrit] <- Crit; }
+
+
+##ErrorCrit______________________________________
+  if(sum(is.na(SubCritValues))==0){
+  CritValue <- ( 1 - sqrt( (SubCritValues[1]-1)^2 + (SubCritValues[2]-1)^2 + (SubCritValues[3]-1)^2 ) );
+  }
+
+
+##Output_________________________________________
+  OutputsCrit <- list(CritValue,CritName,SubCritValues,SubCritNames,CritBestValue,Multiplier,which(TS_ignore));
+  names(OutputsCrit) <- c("CritValue","CritName","SubCritValues","SubCritNames","CritBestValue","Multiplier","Ind_notcomputed");
+  return(OutputsCrit);
+
+}
+
+

R/ErrorCrit_KGE2.R

+#*****************************************************************************************************************
+#' Function which computes an error criterion based on the KGE' formula proposed by Kling et al. (2012).
+#'
+#' In addition to the criterion value, the function outputs include a multiplier (-1 or +1) which allows 
+#' the use of the function for model calibration: the product CritValue*Multiplier is the criterion to be minimised 
+#' (e.g. Multiplier=+1 for RMSE, Multiplier=-1 for NSE).
+#*****************************************************************************************************************
+#' @title  Error criterion based on the KGE' formula
+#' @author Laurent Coron (June 2014)
+#' @references
+#'   Gupta, H. V., Kling, H., Yilmaz, K. K. and Martinez, G. F. (2009), 
+#'       Decomposition of the mean squared error and NSE performance criteria: Implications
+#'       for improving hydrological modelling, Journal of Hydrology, 377(1-2), 80-91, doi:10.1016/j.jhydrol.2009.08.003. \cr
+#'   Kling, H., Fuchs, M. and Paulin, M. (2012), 
+#'       Runoff conditions in the upper Danube basin under an ensemble of climate change scenarios,
+#'       Journal of Hydrology, 424-425, 264-277, doi:10.1016/j.jhydrol.2012.01.011.
+#' @seealso \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}}, \code{\link{ErrorCrit_KGE}}
+#' @examples ## see example of the ErrorCrit function
+#' @encoding UTF-8
+#' @export
+#_FunctionInputs__________________________________________________________________________________________________
+#' @param  InputsCrit      [object of class \emph{InputsCrit}] see \code{\link{CreateInputsCrit}} for details
+#' @param  OutputsModel    [object of class \emph{OutputsModel}] see \code{\link{RunModel_GR4J}} or \code{\link{RunModel_CemaNeigeGR4J}} for details
+#' @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{$CritValue      }   \tab   [numeric] value of the criterion \cr
+#'          \emph{$CritName       }   \tab   [character] name of the criterion \cr
+#'          \emph{$SubCritValues  }   \tab   [numeric] values of the sub-criteria \cr
+#'          \emph{$SubCritNames   }   \tab   [character] names of the sub-criteria \cr
+#'          \emph{$CritBestValue  }   \tab   [numeric] theoretical best criterion value \cr
+#'          \emph{$Multiplier     }   \tab   [numeric] integer indicating whether the criterion is indeed an error (+1) or an efficiency (-1) \cr
+#'          \emph{$Ind_notcomputed}   \tab   [numeric] indices of the time-steps where InputsCrit$BoolCrit=FALSE or no data is available \cr
+#'          }
+#*****************************************************************************************************************'
+ErrorCrit_KGE2 <- function(InputsCrit,OutputsModel,quiet=FALSE){
+
+
+##Arguments_check________________________________
+  if(inherits(InputsCrit,"InputsCrit")==FALSE){ stop("InputsCrit must be of class 'InputsCrit' \n"); return(NULL); }  
+  if(inherits(OutputsModel,"OutputsModel")==FALSE){ stop("OutputsModel must be of class 'OutputsModel' \n"); return(NULL); }  
+
+
+##Initialisation_________________________________
+  CritName <- NA;
+  if(InputsCrit$transfo==""    ){ CritName <- "KGE'[Q]"      ; }
+  if(InputsCrit$transfo=="sqrt"){ CritName <- "KGE'[sqrt(Q)]"; }
+  if(InputsCrit$transfo=="log" ){ CritName <- "KGE'[log(Q)]" ; }
+  if(InputsCrit$transfo=="inv" ){ CritName <- "KGE'[1/Q]"    ; }
+  if(InputsCrit$transfo=="sort"){ CritName <- "KGE'[sort(Q)]"; }
+  CritValue       <- NA;
+  CritBestValue   <- +1;
+  Multiplier      <- -1; ### must be equal to -1 or +1 only
+
+
+##Data_preparation_______________________________
+  VarObs <- InputsCrit$Qobs  ; VarObs[!InputsCrit$BoolCrit] <- NA; 
+  VarSim <- OutputsModel$Qsim; VarSim[!InputsCrit$BoolCrit] <- NA;  
+  ##Data_transformation
+  if("Ind_zeroes" %in% names(InputsCrit) & "epsilon" %in% names(InputsCrit)){ if(length(InputsCrit$Ind_zeroes)>0){
+    VarObs <- VarObs + InputsCrit$epsilon;
+    VarSim <- VarSim + InputsCrit$epsilon;
+  } }
+  if(InputsCrit$transfo=="sqrt"){ VarObs <- sqrt(VarObs); VarSim <- sqrt(VarSim); }
+  if(InputsCrit$transfo=="log" ){ VarObs <- log(VarObs) ; VarSim <- log(VarSim) ; VarSim[VarSim      < -1E100] <- NA; }
+  if(InputsCrit$transfo=="inv" ){ VarObs <- 1/VarObs    ; VarSim <- 1/VarSim    ; VarSim[abs(VarSim) > 1E+100] <- NA; }
+  if(InputsCrit$transfo=="sort"){ VarObs <- sort(VarObs); VarSim <- sort(VarSim); }
+  ##TS_ignore
+  TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit ;
+  if(sum(!TS_ignore)==0){ OutputsCrit <- list(NA); names(OutputsCrit) <- c("CritValue"); return(OutputsCrit); }
+  if(sum(!TS_ignore)<365 & !quiet){ warning("\t criterion computed on less than 365 time-steps \n"); }  
+  ##Other_variables_preparation
+  meanVarObs <- mean(VarObs[!TS_ignore]);
+  meanVarSim <- mean(VarSim[!TS_ignore]);
+  iCrit           <- 0;
+  SubCritNames    <- NULL;
+  SubCritValues   <- NULL;
+
+
+##SubErrorCrit_____KGE_rPearson__________________
+  iCrit <- iCrit+1;
+  SubCritNames[iCrit]  <- paste(CritName," rPEARSON(sim vs. obs)",sep="");
+  SubCritValues[iCrit] <- NA;
+  Numer <- sum( (VarObs[!TS_ignore]-meanVarObs)*(VarSim[!TS_ignore]-meanVarSim) );
+  Deno1 <- sqrt( sum((VarObs[!TS_ignore]-meanVarObs)^2) );
+  Deno2 <- sqrt( sum((VarSim[!TS_ignore]-meanVarSim)^2) );
+  if(Numer==0){ if(Deno1==0 & Deno2==0){ Crit <- 1; } else { Crit <- 0; } 
+  } else { Crit  <- Numer/(Deno1*Deno2); }
+  if(is.numeric(Crit) & is.finite(Crit)){ SubCritValues[iCrit] <- Crit; }
+
+
+##SubErrorCrit_____KGE_gama______________________
+  iCrit <- iCrit+1;
+  SubCritNames[iCrit]  <-  paste(CritName," CVsim/CVobs",sep="");
+  SubCritValues[iCrit] <- NA;
+  if(meanVarSim==0){ if(sd(VarSim[!TS_ignore])==0){ CVsim <- 1; } else { CVsim <- 99999; } } else { CVsim <- sd(VarSim[!TS_ignore])/meanVarSim; }
+  if(meanVarObs==0){ if(sd(VarObs[!TS_ignore])==0){ CVobs <- 1; } else { CVobs <- 99999; } } else { CVobs <- sd(VarObs[!TS_ignore])/meanVarObs; }
+  if(CVsim==0 & CVobs==0){ Crit <- 1; } else { Crit <- CVsim/CVobs ; }
+  if(is.numeric(Crit) & is.finite(Crit)){ SubCritValues[iCrit] <- Crit; }
+
+
+##SubErrorCrit_____KGE_beta______________________
+  iCrit <- iCrit+1;
+  SubCritNames[iCrit]  <-  paste(CritName," MEANsim/MEANobs",sep="");
+  SubCritValues[iCrit] <- NA;
+  if(meanVarSim==0 & meanVarObs==0){ Crit <- 1; } else { Crit <- meanVarSim/meanVarObs ; }
+  if(is.numeric(Crit) & is.finite(Crit)){ SubCritValues[iCrit] <- Crit; }
+
+
+##ErrorCrit______________________________________
+  if(sum(is.na(SubCritValues))==0){
+  CritValue <- ( 1 - sqrt( (SubCritValues[1]-1)^2 + (SubCritValues[2]-1)^2 + (SubCritValues[3]-1)^2 ) );
+  }
+
+
+##Output_________________________________________
+  OutputsCrit <- list(CritValue,CritName,SubCritValues,SubCritNames,CritBestValue,Multiplier,which(TS_ignore));
+  names(OutputsCrit) <- c("CritValue","CritName","SubCritValues","SubCritNames","CritBestValue","Multiplier","Ind_notcomputed");
+  return(OutputsCrit);
+
+}
+
+