RunModel_CemaNeige.R 6.74 KB
Newer Older
Delaigue Olivier's avatar
Delaigue Olivier committed
1
2
RunModel_CemaNeige <- function(InputsModel,RunOptions,Param){

3
    NParam <- 2
4
    FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp");
Delaigue Olivier's avatar
Delaigue Olivier committed
5
6
7
8
9
10
11

    ##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,"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,"CemaNeige"   )==FALSE){ stop("RunOptions must be of class 'CemaNeige'    \n"); return(NULL); }  
12
      if(!is.vector(Param) | !is.numeric(Param)){ stop("Param must be a numeric vector \n"); return(NULL); }
Delaigue Olivier's avatar
Delaigue Olivier committed
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
      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);

    ##Input_data_preparation
      if(identical(RunOptions$IndPeriod_WarmUp,0)){ RunOptions$IndPeriod_WarmUp <- NULL; }
      IndPeriod1    <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run);
      IndPeriod2    <- (length(RunOptions$IndPeriod_WarmUp)+1):length(IndPeriod1);
      ExportDatesR    <- "DatesR"   %in% RunOptions$Outputs_Sim;
      ExportStateEnd  <- "StateEnd" %in% RunOptions$Outputs_Sim;



    ##SNOW_MODULE________________________________________________________________________________##
      ParamCemaNeige  <- Param;
      NLayers         <- length(InputsModel$LayerPrecip);
      if(sum(is.na(ParamCemaNeige))!=0){ stop("Param contains missing values \n"); return(NULL); }
      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){
35
        StateStartCemaNeige <- RunOptions$IniStates[(7+20+40) + c(iLayer, iLayer+NLayers)]
36
        RESULTS <- .Fortran("frun_CemaNeige",PACKAGE="airGR",
Delaigue Olivier's avatar
Delaigue Olivier committed
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
                        ##inputs
                            LInputs=as.integer(length(IndPeriod1)),                                        ### 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=length(IndPeriod1),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;
55
        
Delaigue Olivier's avatar
Delaigue Olivier committed
56
57
58
59
60
61
62
63

        ##Data_storage
        CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
        names(CemaNeigeLayers[[iLayer]]) <- FortranOutputsCemaNeige[IndOutputsCemaNeige];
        if(ExportStateEnd){ CemaNeigeStateEnd <- c(CemaNeigeStateEnd,RESULTS$StateEnd); }
        rm(RESULTS); 
      } ###ENDFOR_iLayer
      names(CemaNeigeLayers) <- paste("Layer",formatC(1:NLayers,width=2,flag="0"),sep="");
64
65
66
67
68
69
70
71
72
73
      
      if (ExportStateEnd) { 
        CemaNeigeStateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeige, InputsModel = InputsModel,
                                             ProdStore = NULL, RoutStore = NULL, ExpStore = NULL,
                                             UH1 = NULL, UH2 = NULL,
                                             GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(2*NLayers)[seq_len(2*NLayers) %% 2 == 1]],
                                             eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(2*NLayers)[seq_len(2*NLayers) %% 2 == 0]],
                                             verbose = FALSE)
      }
      
Delaigue Olivier's avatar
Delaigue Olivier committed
74
75
76
77
78
79
80
81
82
83
      ##Output_data_preparation
      if(ExportDatesR==FALSE & ExportStateEnd==FALSE){
        OutputsModel <- list(CemaNeigeLayers);
        names(OutputsModel) <- NameCemaNeigeLayers;   }
      if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
        OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
                           list(CemaNeigeLayers) );
        names(OutputsModel) <- c("DatesR",NameCemaNeigeLayers);   }
      if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
        OutputsModel <- c( list(CemaNeigeLayers),
84
                           list(CemaNeigeStateEnd));
Delaigue Olivier's avatar
Delaigue Olivier committed
85
86
87
88
        names(OutputsModel) <- c(NameCemaNeigeLayers,"StateEnd");   }
      if(ExportDatesR==TRUE & ExportStateEnd==TRUE){
        OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
                           list(CemaNeigeLayers),
89
                           list(CemaNeigeStateEnd));
Delaigue Olivier's avatar
Delaigue Olivier committed
90
91
92
93
94
95
96
97
        names(OutputsModel) <- c("DatesR",NameCemaNeigeLayers,"StateEnd");   }

    ##End
      class(OutputsModel) <- c("OutputsModel","daily","CemaNeige");
      return(OutputsModel);

}