RunModel_CemaNeigeGR5J.R

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

    NParam <- 7;
    FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp");
    FortranOutputsMod       <- c("PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim");

    ##Arguments_check
      if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
      if(inherits(InputsModel,"daily"      )==FALSE){ stop("InputsModel must be of class 'daily'       \n"); return(NULL); }  
      if(inherits(InputsModel,"GR"         )==FALSE){ stop("InputsModel must be of class 'GR'          \n"); return(NULL); }  
      if(inherits(InputsModel,"CemaNeige"  )==FALSE){ stop("InputsModel must be of class 'CemaNeige'   \n"); return(NULL); }  
      if(inherits(RunOptions,"RunOptions"  )==FALSE){ stop("RunOptions must be of class 'RunOptions'   \n"); return(NULL); }  
      if(inherits(RunOptions,"GR"          )==FALSE){ stop("RunOptions must be of class 'GR'           \n"); return(NULL); }  
      if(inherits(RunOptions,"CemaNeige"   )==FALSE){ stop("RunOptions must be of class 'CemaNeige'    \n"); return(NULL); }  
      if(!is.vector(Param)){ stop("Param must be a vector \n"); return(NULL); }
      if(sum(!is.na(Param))!=NParam){ stop(paste("Param must be a vector of length ",NParam," and contain no NA \n",sep="")); return(NULL); }
      Param <- as.double(Param);
      
      Param_X1X3_threshold <- 1e-2
      if (Param[1L] < Param_X1X3_threshold) {
        warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n New X1 value: %.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 New X3 value: %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
        Param[3L] <- Param_X1X3_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):length(Param)];
      NParamMod      <- as.integer(length(Param)-2);
      ParamMod       <- Param[1:NParamMod];
      NLayers        <- length(InputsModel$LayerPrecip);
      NStatesMod     <- as.integer(length(RunOptions$IniStates)-2*NLayers);
      ExportDatesR   <- "DatesR"   %in% RunOptions$Outputs_Sim;
      ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;


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

    ##Call_DLL_CemaNeige_________________________
      for(iLayer in 1:NLayers){
        StateStartCemaNeige <- RunOptions$IniStates[ (NStatesMod+2*(iLayer-1)+1):(NStatesMod+2*(iLayer-1)+2) ];
        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]
                            MeanAnSolidPrecip=RunOptions$MeanAnSolidPrecip[iLayer],                        ### value of annual mean solid precip [mm/y]
                            NParam=as.integer(2),                                                          ### number of model parameter = 2
                            Param=ParamCemaNeige,                                                          ### parameter set
                            NStates=as.integer(2),                                                         ### number of state variables used for model initialising = 2
                            StateStart=StateStartCemaNeige,                                                ### state variables used when the model run starts
                            NOutputs=as.integer(length(IndOutputsCemaNeige)),                              ### number of output series
                            IndOutputs=IndOutputsCemaNeige,                                                ### indices of output series
                        ##outputs                                                               
                            Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsCemaNeige)),  ### output series [mm]
                            StateEnd=rep(as.double(-999.999),as.integer(2))                                          ### state variables at the end of the model run (reservoir levels [mm] and HU)
                         )
        RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
        RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;