RunModel_CemaNeige.R 7.11 KB
Newer Older
1
RunModel_CemaNeige <- function(InputsModel, RunOptions, Param) {
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
  
  
  ## 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) {
68
69
70
71
72
73
    
    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)]
    }
74
75
76
77
78
79
80
    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]
81
                        NParam = as.integer(NParam),                                                    ### number of model parameter
82
                        Param = as.double(ParamCemaNeige),                                              ### parameter set
83
                        NStates = as.integer(NStates),                                                  ### number of state variables used for model initialising
84
                        StateStart = StateStartCemaNeige,                                               ### state variables used when the model run starts
85
                        IsHyst = as.integer(IsHyst),                                                    ### use of hysteresis
86
87
88
89
90
91
                        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)),
92
                        StateEnd = rep(-999.999, NStates)                                               ### state variables at the end of the model run (reservoir levels [mm] and HU)
93
    )
94
95
    RESULTS$Outputs[round(RESULTS$Outputs  , 3) == -999.999] <- NA
    RESULTS$StateEnd[round(RESULTS$StateEnd, 3) == -999.999] <- NA
96
97
98
99
100
101
102
103
104
105
106
107
108
109
    
    
    
    
    ## 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
  
110
  names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
111
112
  
  if (ExportStateEnd) { 
113
    idNStates <- seq_len(NStates*NLayers) %% NStates
114
    CemaNeigeStateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeige, InputsModel = InputsModel, IsHyst = IsHyst,
115
116
                                         ProdStore = NULL, RoutStore = NULL, ExpStore = NULL,
                                         UH1 = NULL, UH2 = NULL,
117
118
119
120
                                         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]],
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
                                         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")
149
150
151
  if(IsHyst) {
    class(OutputsModel) <- c(class(OutputsModel), "hysteresis")
  }
152
153
  return(OutputsModel)
}