RunModel_CemaNeigeGR5J.R

RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){

  
  IsHyst <- inherits(RunOptions, "hysteresis")
  NParam <- ifelse(test = IsHyst, yes = 9L, no = 7L)
  NParamCN <- NParam - 5L
  NStates <- 4L
  FortranOutputs <- .FortranOutputs(GR = "GR5J", isCN = TRUE)

    ##Arguments_check
      if(inherits(InputsModel,"InputsModel")==FALSE){ stop("'InputsModel' must be of class 'InputsModel'") }  
      if(inherits(InputsModel,"daily"      )==FALSE){ stop("'InputsModel' must be of class 'daily'      ") }  
      if(inherits(InputsModel,"GR"         )==FALSE){ stop("'InputsModel' must be of class 'GR'         ") }  
      if(inherits(InputsModel,"CemaNeige"  )==FALSE){ stop("'InputsModel' must be of class 'CemaNeige'  ") }  
      if(inherits(RunOptions,"RunOptions"  )==FALSE){ stop("'RunOptions' must be of class 'RunOptions'  ") }  
      if(inherits(RunOptions,"GR"          )==FALSE){ stop("'RunOptions' must be of class 'GR'          ") }  
      if(inherits(RunOptions,"CemaNeige"   )==FALSE){ stop("'RunOptions' must be of class 'CemaNeige'   ") }  
      if(!is.vector(Param) | !is.numeric(Param)){ stop("'Param' must be a numeric vector") }
      if(sum(!is.na(Param))!=NParam){ stop(paste("'Param' must be a vector of length ",NParam," and contain no NA",sep="")) }
      Param <- as.double(Param);
      
      Param_X1X3_threshold <- 1e-2
      Param_X4_threshold   <- 0.5
      if (Param[1L] < Param_X1X3_threshold) {
        warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
        Param[1L] <- Param_X1X3_threshold
      }
      if (Param[3L] < Param_X1X3_threshold) {
        warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
        Param[3L] <- Param_X1X3_threshold
      }
      if (Param[4L] < Param_X4_threshold) {
        warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
        Param[4L] <- Param_X4_threshold
      }      

    ##Input_data_preparation
      if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; }
      IndPeriod1     <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run);
      LInputSeries   <- as.integer(length(IndPeriod1))
      IndPeriod2     <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
      ParamCemaNeige <- Param[(length(Param)-1-2*as.integer(IsHyst)):length(Param)];
      NParamMod      <- as.integer(length(Param)-(2+2*as.integer(IsHyst)));
      ParamMod       <- Param[1:NParamMod];
      NLayers        <- length(InputsModel$LayerPrecip);
      NStatesMod     <- as.integer(length(RunOptions$IniStates)-NStates*NLayers);
      ExportDatesR   <- "DatesR"   %in% RunOptions$Outputs_Sim;
      ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;


    ##SNOW_MODULE________________________________________________________________________________##
    if(inherits(RunOptions,"CemaNeige")==TRUE){
      if("all" %in% RunOptions$Outputs_Sim){ IndOutputsCemaNeige <- as.integer(1:length(FortranOutputs$CN)); 
      } else { IndOutputsCemaNeige <- which(FortranOutputs$CN %in% RunOptions$Outputs_Sim);  }
      CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- "CemaNeigeLayers";

      
    ##Call_DLL_CemaNeige_________________________
      for(iLayer in 1:NLayers){
        if (!IsHyst) {
          StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers)]
        } else {
          StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
        }
        RESULTS <- .Fortran("frun_cemaneige",PACKAGE="airGR",
                        ##inputs
                            LInputs=LInputSeries,                                                         ### length of input and output series
                            InputsPrecip=InputsModel$LayerPrecip[[iLayer]][IndPeriod1],                   ### input series of total precipitation [mm/d]
                            InputsFracSolidPrecip=InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
                            InputsTemp=InputsModel$LayerTemp[[iLayer]][IndPeriod1],                       ### input series of air mean temperature [degC]