RunModel_CemaNeige.R 6.94 KB
Newer Older
1
RunModel_CemaNeige <- function(InputsModel, RunOptions, Param, IsHyst = FALSE) {
2
3
4

  
  ## Initialization of variables
5
6
  IsHyst <- inherits(RunOptions, "hysteresis")
  NParam <- ifelse(test = IsHyst, yes = 4L, no = 2L)
7
  NStates <- 4L
8
  FortranOutputsCemaNeige <- .FortranOutputs(GR = NULL, isCN = TRUE)$CN
9
  
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
  
  ## Arguments_check
  if (!inherits(InputsModel, "InputsModel")) {
    stop("'InputsModel' must be of class 'InputsModel'")
  }
  if (!inherits(InputsModel, "daily")) {
    stop("'InputsModel' must be of class 'daily'")
  }
  if (!inherits(InputsModel, "CemaNeige")) {
    stop("'InputsModel' must be of class 'CemaNeige'")
  }
  if (!inherits(RunOptions, "RunOptions")) {
    stop("'RunOptions' must be of class 'RunOptions'")
  }
  if (!inherits(RunOptions, "CemaNeige")) {
    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(sprintf("'Param' must be a vector of length %i and contain no NA", NParam))
  }
  
  ## 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
  
43
  
Delaigue Olivier's avatar
Delaigue Olivier committed
44
  
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
  
  
  ## SNOW_MODULE________________________________________________________________________________
  ParamCemaNeige <- Param
  NLayers        <- length(InputsModel$LayerPrecip)
  
  if (sum(is.na(ParamCemaNeige)) != 0) {
    stop("Param contains missing values")
  }
  if ("all" %in% RunOptions$Outputs_Sim) {
    IndOutputsCemaNeige <- 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[(7 + 20 + 40) + c(iLayer, iLayer + NLayers)]
    RESULTS <- .Fortran("frun_CemaNeige", PACKAGE = "airGR",
                        ## 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]
76
                        NParam = as.integer(NParam),                                                    ### number of model parameter
77
                        Param = as.double(ParamCemaNeige),                                              ### parameter set
78
                        NStates = as.integer(NStates),                                                  ### number of state variables used for model initialising
79
                        StateStart = StateStartCemaNeige,                                               ### state variables used when the model run starts
80
                        IsHyst = as.integer(IsHyst),                                                    ### use of hysteresis
81
82
83
84
85
86
                        NOutputs = as.integer(length(IndOutputsCemaNeige)),                             ### number of output series
                        IndOutputs = IndOutputsCemaNeige,                                               ### indices of output series
                        ## outputs
                        Outputs = matrix(-999.999,                                                      ### output series [mm]
                                         nrow = length(IndPeriod1),
                                         ncol = length(IndOutputsCemaNeige)),
87
                        StateEnd = rep(-999.999, NStates)                                               ### state variables at the end of the model run (reservoir levels [mm] and HU)
88
    )
89
90
    RESULTS$Outputs[round(RESULTS$Outputs  , 3) == -999.999] <- NA
    RESULTS$StateEnd[round(RESULTS$StateEnd, 3) == -999.999] <- NA
91
92
93
94
95
96
97
98
99
100
101
102
103
104
    
    
    
    
    ## 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
  
105
  names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
106
107
  
  if (ExportStateEnd) { 
108
    idNStates <- seq_len(NStates*NLayers) %% NStates
109
110
111
    CemaNeigeStateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeige, InputsModel = InputsModel,
                                         ProdStore = NULL, RoutStore = NULL, ExpStore = NULL,
                                         UH1 = NULL, UH2 = NULL,
112
113
114
115
                                         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]],
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
                                         verbose = FALSE)
  }
  
  ## Output_data_preparation
  if (!ExportDatesR & !ExportStateEnd) {
    OutputsModel <- list(CemaNeigeLayers)
    names(OutputsModel) <- NameCemaNeigeLayers
  }
  if (ExportDatesR & !ExportStateEnd) {
    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
                      list(CemaNeigeLayers))
    names(OutputsModel) <- c("DatesR", NameCemaNeigeLayers)
  }
  if (!ExportDatesR & ExportStateEnd) {
    OutputsModel <- c(list(CemaNeigeLayers),
                      list(CemaNeigeStateEnd))
    names(OutputsModel) <- c(NameCemaNeigeLayers, "StateEnd")
  }
  if (ExportDatesR & ExportStateEnd) {
    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
                      list(CemaNeigeLayers),
                      list(CemaNeigeStateEnd))
    names(OutputsModel) <- c("DatesR", NameCemaNeigeLayers, "StateEnd")
  }
  
  
  ## End
  class(OutputsModel) <- c("OutputsModel", "daily", "CemaNeige")
144
145
146
  if(IsHyst) {
    class(OutputsModel) <- c(class(OutputsModel), "hysteresis")
  }
147
148
149
  return(OutputsModel)
}