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

3
  
4
  IsHyst <- inherits(RunOptions, "hysteresis")
5
  NParam <- ifelse(test = IsHyst, yes = 10L, no = 8L)
6
  NStates <- 4L
7
  FortranOutputs <- .FortranOutputs(GR = "GR6J", isCN = TRUE)
Delaigue Olivier's avatar
Delaigue Olivier committed
8
9

    ##Arguments_check
10
11
12
13
14
15
16
17
18
      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="")) }
Delaigue Olivier's avatar
Delaigue Olivier committed
19
      Param <- as.double(Param);
20
      
21
      Param_X1X3X6_threshold <- 1e-2
22
      Param_X4_threshold     <- 0.5
23
      if (Param[1L] < Param_X1X3X6_threshold) {
24
        warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
25
        Param[1L] <- Param_X1X3X6_threshold
26
      }
27
      if (Param[3L] < Param_X1X3X6_threshold) {
28
        warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
29
        Param[3L] <- Param_X1X3X6_threshold
30
      }
31
32
33
34
      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
      }      
35
      if (Param[6L] < Param_X1X3X6_threshold) {
36
        warning(sprintf("Param[6] (X6: coefficient for emptying exponential store [mm]) < %.2f\n X6 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
37
38
        Param[6L] <- Param_X1X3X6_threshold
      }         
Delaigue Olivier's avatar
Delaigue Olivier committed
39
40
41
42
43
44

    ##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;
45
46
      ParamCemaNeige <- Param[(length(Param)-1-2*as.integer(IsHyst)):length(Param)];
      NParamMod      <- as.integer(length(Param)-(2+2*as.integer(IsHyst)));
Delaigue Olivier's avatar
Delaigue Olivier committed
47
48
      ParamMod       <- Param[1:NParamMod];
      NLayers        <- length(InputsModel$LayerPrecip);
49
      NStatesMod     <- as.integer(length(RunOptions$IniStates)-NStates*NLayers);
Delaigue Olivier's avatar
Delaigue Olivier committed
50
51
52
53
54
      ExportDatesR   <- "DatesR"   %in% RunOptions$Outputs_Sim;
      ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;


    ##SNOW_MODULE________________________________________________________________________________##
55
    if(inherits(RunOptions,"CemaNeige")==TRUE){
56
57
      if("all" %in% RunOptions$Outputs_Sim){ IndOutputsCemaNeige <- as.integer(1:length(FortranOutputs$CN)); 
      } else { IndOutputsCemaNeige <- which(FortranOutputs$CN %in% RunOptions$Outputs_Sim);  }
Delaigue Olivier's avatar
Delaigue Olivier committed
58
59
      CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- "CemaNeigeLayers";

60
      
Delaigue Olivier's avatar
Delaigue Olivier committed
61
62
    ##Call_DLL_CemaNeige_________________________
      for(iLayer in 1:NLayers){
63
64
65
66
67
        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)]
        }
68
        RESULTS <- .Fortran("frun_CemaNeige",PACKAGE="airGR",
Delaigue Olivier's avatar
Delaigue Olivier committed
69
70
71
72
73
74
                        ##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]
75
76
77
                            NParam=as.integer(NParam),                                                     ### number of model parameter = 2
                            Param=as.double(ParamCemaNeige),                                               ### parameter set
                            NStates=as.integer(NStates),                                                   ### number of state variables used for model initialising = 2
Delaigue Olivier's avatar
Delaigue Olivier committed
78
                            StateStart=StateStartCemaNeige,                                                ### state variables used when the model run starts
79
                            IsHyst = as.integer(IsHyst),                       ### use of hysteresis                        
Delaigue Olivier's avatar
Delaigue Olivier committed
80
81
82
83
                            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]
84
                            StateEnd=rep(as.double(-999.999),as.integer(NStates))                                    ### state variables at the end of the model run (reservoir levels [mm] and HU)
Delaigue Olivier's avatar
Delaigue Olivier committed
85
86
87
88
89
90
                         )
        RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
        RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;

        ##Data_storage
        CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
91
        names(CemaNeigeLayers[[iLayer]]) <- FortranOutputs$CN[IndOutputsCemaNeige];
Delaigue Olivier's avatar
Delaigue Olivier committed
92
93
94
95
96
97
        IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt");
        if(iLayer==1){ CatchMeltAndPliq <- RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
        if(iLayer >1){ CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
        if(ExportStateEnd){ CemaNeigeStateEnd <- c(CemaNeigeStateEnd,RESULTS$StateEnd); }
        rm(RESULTS); 
      } ###ENDFOR_iLayer
98
      names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
Delaigue Olivier's avatar
Delaigue Olivier committed
99
    } ###ENDIF_RunSnowModule
100
    if(inherits(RunOptions,"CemaNeige")==FALSE){
Delaigue Olivier's avatar
Delaigue Olivier committed
101
102
103
104
105
106
      CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- NULL;
      CatchMeltAndPliq  <- InputsModel$Precip[IndPeriod1]; }



    ##MODEL______________________________________________________________________________________##
107
108
      if("all" %in% RunOptions$Outputs_Sim){ IndOutputsMod <- as.integer(1:length(FortranOutputs$GR)); 
      } else { IndOutputsMod <- which(FortranOutputs$GR %in% RunOptions$Outputs_Sim);  }
Delaigue Olivier's avatar
Delaigue Olivier committed
109
110

    ##Use_of_IniResLevels
111
112
113
114
      if(!is.null(RunOptions$IniResLevels)){
        RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * ParamMod[1] ### production store level (mm)
        RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * ParamMod[3] ### routing store level (mm)
        RunOptions$IniStates[3] <- RunOptions$IniResLevels[3]               ### exponential store level (mm)
Delaigue Olivier's avatar
Delaigue Olivier committed
115
116
117
      }

    ##Call_fortan
118
      RESULTS <- .Fortran("frun_GR6J",PACKAGE="airGR",
Delaigue Olivier's avatar
Delaigue Olivier committed
119
                 ##inputs
120
121
122
123
124
125
                     LInputs=LInputSeries,                          ### length of input and output series
                     InputsPrecip=CatchMeltAndPliq,                 ### input series of total precipitation [mm/d]
                     InputsPE=InputsModel$PotEvap[IndPeriod1],      ### input series potential evapotranspiration [mm/d]
                     NParam=NParamMod,                              ### number of model parameter
                     Param=ParamMod,                                ### parameter set
                     NStates=NStatesMod,                            ### number of state variables used for model initialising
Delaigue Olivier's avatar
Delaigue Olivier committed
126
                     StateStart=RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts
127
128
                     NOutputs=as.integer(length(IndOutputsMod)),    ### number of output series
                     IndOutputs=IndOutputsMod,                      ### indices of output series
Delaigue Olivier's avatar
Delaigue Olivier committed
129
                 ##outputs                                        
130
131
                     Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsMod)), ### output series [mm]
                     StateEnd=rep(as.double(-999.999),NStatesMod)                                      ### state variables at the end of the model run
Delaigue Olivier's avatar
Delaigue Olivier committed
132
133
134
                     )
      RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
      RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
135
      if (ExportStateEnd) { 
136
        idNStates <- seq_len(NStates*NLayers) %% NStates
137
        RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR6J, InputsModel = InputsModel, IsHyst = IsHyst,
138
139
                                            ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = RESULTS$StateEnd[3L],
                                            UH1 = RESULTS$StateEnd[(1:20)+7], UH2 = RESULTS$StateEnd[(1:40)+(7+20)],
140
141
142
143
                                            GCemaNeigeLayers       = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
                                            eTGCemaNeigeLayers     = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
                                            GthrCemaNeigeLayers    = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]], 
                                            GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
144
145
                                            verbose = FALSE)
      }
146
      
147
      if(inherits(RunOptions,"CemaNeige")==TRUE & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputs$GR[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; }
Delaigue Olivier's avatar
Delaigue Olivier committed
148
149
150
151
152
153

    ##Output_data_preparation
      ##OutputsModel_only
      if(ExportDatesR==FALSE & ExportStateEnd==FALSE){
        OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                           list(CemaNeigeLayers) );
154
        names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers); }
Delaigue Olivier's avatar
Delaigue Olivier committed
155
156
157
158
159
      ##DatesR_and_OutputsModel_only
      if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
        OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
                           lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                           list(CemaNeigeLayers) );
160
        names(OutputsModel) <- c("DatesR",FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers);      }
Delaigue Olivier's avatar
Delaigue Olivier committed
161
162
163
164
      ##OutputsModel_and_SateEnd_only
      if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
        OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                           list(CemaNeigeLayers),
165
                           list(RESULTS$StateEnd) );
166
        names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers,"StateEnd");      }
Delaigue Olivier's avatar
Delaigue Olivier committed
167
168
169
170
171
      ##DatesR_and_OutputsModel_and_SateEnd
      if(ExportDatesR==TRUE & ExportStateEnd==TRUE){
        OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
                           lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                           list(CemaNeigeLayers),
172
                           list(RESULTS$StateEnd) );
173
        names(OutputsModel) <- c("DatesR",FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers,"StateEnd");      }
Delaigue Olivier's avatar
Delaigue Olivier committed
174
175
176
177

    ##End
      rm(RESULTS); 
      class(OutputsModel) <- c("OutputsModel","daily","GR","CemaNeige");
178
179
180
      if(IsHyst) {
        class(OutputsModel) <- c(class(OutputsModel), "hysteresis")
      }
Delaigue Olivier's avatar
Delaigue Olivier committed
181
182
183
184
      return(OutputsModel);

}