CreateRunOptions.R 15.5 KB
Newer Older
1
2
3
CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndPeriod_Run,
                             IniStates = NULL, IniResLevels = NULL, 
                             Outputs_Cal = NULL, Outputs_Sim = "all",
4
                             RunSnowModule, MeanAnSolidPrecip = NULL, IsHyst = FALSE, 
5
                             warnings = TRUE, verbose = TRUE) {
6
7
  
  if (!missing(RunSnowModule)) {
8
    warning("argument 'RunSnowModule' is deprecated; please adapt 'FUN_MOD' instead.", call. = FALSE)
9
  }
10
11
12
  if (!is.logical(IsHyst) | length(IsHyst) != 1L) {
    stop("'IsHyst' must be a 'logical' of length 1")
  }
13
  ObjectClass <- NULL
14
  
15
16
  FUN_MOD <- match.fun(FUN_MOD)
  
Delaigue Olivier's avatar
Delaigue Olivier committed
17
  ##check_FUN_MOD
18
  BOOL <- FALSE
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
  if (identical(FUN_MOD, RunModel_GR4H)) {
    ObjectClass <- c(ObjectClass, "GR", "hourly")
    BOOL <- TRUE
  }
  if (identical(FUN_MOD, RunModel_GR4J) | identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_GR6J)) {
    ObjectClass <- c(ObjectClass, "GR", "daily")
    BOOL <- TRUE
  }
  if (identical(FUN_MOD, RunModel_GR2M)) {
    ObjectClass <- c(ObjectClass, "GR", "monthly")
    BOOL <- TRUE
  }
  if (identical(FUN_MOD, RunModel_GR1A)) {
    ObjectClass <- c(ObjectClass, "GR", "yearly")
    BOOL <- TRUE
  }
  if (identical(FUN_MOD, RunModel_CemaNeige)) {
    ObjectClass <- c(ObjectClass, "CemaNeige", "daily")
    BOOL <- TRUE
  }
  if (identical(FUN_MOD, RunModel_CemaNeigeGR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
    ObjectClass <- c(ObjectClass, "GR", "CemaNeige", "daily")
    BOOL <- TRUE
  }
43
44
45
  if (IsHyst) {
    ObjectClass <- c(ObjectClass, "hysteresis")
  }
46
  if (!BOOL) {
47
    stop("incorrect 'FUN_MOD' for use in 'CreateRunOptions'")
48
49
  }
  
Delaigue Olivier's avatar
Delaigue Olivier committed
50
  ##check_InputsModel
51
  if (!inherits(InputsModel, "InputsModel")) {
52
    stop("'InputsModel' must be of class 'InputsModel'")
53
54
  }
  if ("GR" %in% ObjectClass & !inherits(InputsModel, "GR")) {
55
    stop("'InputsModel' must be of class 'GR'")
56
57
58
  }
  if ("CemaNeige" %in% ObjectClass &
      !inherits(InputsModel, "CemaNeige")) {
59
    stop("'InputsModel' must be of class 'CemaNeige'")
60
61
62
  }
  if ("hourly" %in% ObjectClass &
      !inherits(InputsModel, "hourly")) {
63
    stop("'InputsModel' must be of class 'hourly'")
64
65
  }
  if ("daily" %in% ObjectClass & !inherits(InputsModel, "daily")) {
66
    stop("'InputsModel' must be of class 'daily'")
67
68
69
  }
  if ("monthly" %in% ObjectClass &
      !inherits(InputsModel, "monthly")) {
70
    stop("'InputsModel' must be of class 'monthly'")
71
72
73
  }
  if ("yearly" %in% ObjectClass &
      !inherits(InputsModel, "yearly")) {
74
    stop("'InputsModel' must be of class 'yearly'")
75
76
77
  }
  
  
Delaigue Olivier's avatar
Delaigue Olivier committed
78
  ##check_IndPeriod_Run
79
  if (!is.vector(IndPeriod_Run)) {
80
    stop("'IndPeriod_Run' must be a vector of numeric values")
81
82
  }
  if (!is.numeric(IndPeriod_Run)) {
83
    stop("'IndPeriod_Run' must be a vector of numeric values")
84
  }
85
  if (!identical(as.integer(IndPeriod_Run), as.integer(seq(from = IndPeriod_Run[1], to = tail(IndPeriod_Run, 1), by = 1)))) {
86
    stop("'IndPeriod_Run' must be a continuous sequence of integers")
87
88
  }
  if (storage.mode(IndPeriod_Run) != "integer") {
89
    stop("'IndPeriod_Run' should be of type integer")
90
91
92
  }
  
  
Delaigue Olivier's avatar
Delaigue Olivier committed
93
  ##check_IndPeriod_WarmUp
94
95
  WTxt <- NULL
  if (is.null(IndPeriod_WarmUp)) {
96
    WTxt <- paste0(WTxt,"\t Model warm up period not defined -> default configuration used")
97
    ##If_the_run_period_starts_at_the_very_beginning_of_the_time_series
98
    if (IndPeriod_Run[1L] == 1L) {
99
      IndPeriod_WarmUp <- as.integer(0)
100
      WTxt <- paste0(WTxt,"\t    No data were found for model warm up!")
Delaigue Olivier's avatar
Delaigue Olivier committed
101
      ##We_look_for_the_longest_period_preceeding_the_run_period_with_a_maximum_of_one_year
102
103
104
105
106
107
108
    } else {
      TmpDateR0 <- InputsModel$DatesR[IndPeriod_Run[1]]
      TmpDateR  <- TmpDateR0 - 365 * 24 * 60 * 60
      ### minimal date to start the warmup
      if (format(TmpDateR, format = "%d") != format(TmpDateR0, format = "%d")) {
        ### leap year
        TmpDateR <- TmpDateR - 1 * 24 * 60 * 60
Delaigue Olivier's avatar
Delaigue Olivier committed
109
      }
110
111
112
      IndPeriod_WarmUp <- which(InputsModel$DatesR == max(InputsModel$DatesR[1], TmpDateR)):(IndPeriod_Run[1] - 1)
      if ("hourly"  %in% ObjectClass) {
        TimeStep <- as.integer(60 * 60)
113
      }
114
115
      if ("daily"   %in% ObjectClass) {
        TimeStep <- as.integer(24 * 60 * 60)
116
      }
117
118
      if ("monthly" %in% ObjectClass) {
        TimeStep <- as.integer(30.44 * 24 * 60 * 60)
119
      }
120
121
      if ("yearly"  %in% ObjectClass) {
        TimeStep <- as.integer(365.25 * 24 * 60 * 60)
122
      }
123
124
125
126
127
      if (length(IndPeriod_WarmUp) * TimeStep / (365 * 24 * 60 * 60) >= 1) {
        WTxt <- paste0(WTxt, "\t    The year preceding the run period is used \n")
      } else {
        WTxt <- paste0(WTxt, "\t    Less than a year (without missing values) was found for model warm up: \n")
        WTxt <- paste0(WTxt, "\t    (", length(IndPeriod_WarmUp), " time-steps are used for initialisation) \n")
128
      }
129
130
131
132
    }
  }
  if (!is.null(IndPeriod_WarmUp)) {
    if (!is.vector(IndPeriod_WarmUp)) {
133
      stop("'IndPeriod_WarmUp' must be a vector of numeric values")
134
135
    }
    if (!is.numeric(IndPeriod_WarmUp)) {
136
      stop("'IndPeriod_WarmUp' must be a vector of numeric values")
137
138
    }
    if (storage.mode(IndPeriod_WarmUp) != "integer") {
139
      stop("'IndPeriod_WarmUp' should be of type integer")
140
    }
141
142
    if (identical(IndPeriod_WarmUp, as.integer(0)) & verbose) {
      message(paste0(WTxt, "\t No warm up period is used \n"))
143
144
145
146
147
    }
    if ((IndPeriod_Run[1] - 1) != tail(IndPeriod_WarmUp, 1) & !identical(IndPeriod_WarmUp, as.integer(0))) {
      WTxt <- paste0(WTxt, "\t Model warm up period is not directly before the model run period \n")
    }
  }
148
  if (!is.null(WTxt) & warnings) {
149
150
151
152
153
154
155
156
    warning(WTxt)
  }
  
  
  ## check IniResLevels    
  if ("GR" %in% ObjectClass & ("monthly" %in% ObjectClass | "daily" %in% ObjectClass | "hourly" %in% ObjectClass)) {
    if (!is.null(IniResLevels)) {
      if (!is.vector(IniResLevels) | !is.numeric(IniResLevels) | any(is.na(IniResLevels))) {
157
        stop("'IniResLevels' must be a vector of numeric values")
Delaigue Olivier's avatar
Delaigue Olivier committed
158
      }
159
160
161
162
163
      if ((identical(FUN_MOD, RunModel_GR4H) |
           identical(FUN_MOD, RunModel_GR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR4J) |
           identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) |
           identical(FUN_MOD, RunModel_GR2M)) &
          length(IniResLevels) != 2) {
164
        stop("The length of 'IniResLevels' must be 2 for the chosen 'FUN_MOD'")
165
      }
166
167
      if ((identical(FUN_MOD,RunModel_GR6J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)) &
          length(IniResLevels) != 3) {
168
        stop("The length of 'IniResLevels' must be 3 for the chosen 'FUN_MOD'")
169
      }
170
171
172
173
174
    } else if (is.null(IniStates)) {
      if (identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
        IniResLevels <- as.double(c(0.3, 0.5, 0))
      } else {
        IniResLevels <- as.double(c(0.3, 0.5, NA))
175
      }
Delaigue Olivier's avatar
Delaigue Olivier committed
176
    }
177
178
  } else {
    if (!is.null(IniResLevels)) {
179
      stop("'IniResLevels' can only be used with monthly or daily or hourly GR models")
180
181
182
    }
  }
  ## check IniStates
183
  if (is.null(IniStates) & is.null(IniResLevels) & warnings) {
184
    warning("\t Model states initialisation not defined -> default configuration used")
185
  }
186
  if (!is.null(IniStates) & !is.null(IniResLevels) & warnings) {
187
    warning("\t 'IniStates' and 'IniResLevels' are both defined -> Store levels are taken from 'IniResLevels'")
188
189
190
191
192
193
194
195
196
197
198
199
  }
  if ("CemaNeige" %in% ObjectClass) {
    NLayers <- length(InputsModel$LayerPrecip)
  } else {
    NLayers <- 0
  }
  NState <- NULL
  if ("GR" %in% ObjectClass | "CemaNeige" %in% ObjectClass) {
    if ("hourly"  %in% ObjectClass) {
      NState <- 7 + 3 * 24 * 20
    }
    if ("daily"   %in% ObjectClass) {
200
      NState <- 7 + 3 * 20 + 4 * NLayers
201
202
203
204
205
206
207
208
209
210
211
    }
    if ("monthly" %in% ObjectClass) {
      NState <- 2
    }
    if ("yearly"  %in% ObjectClass) {
      NState <- 1
    }
  }
  if (!is.null(IniStates)) {
    
    if (!inherits(IniStates, "IniStates")) {
212
      stop("'IniStates' must be an object of class 'IniStates'\n")
213
214
    }
    if (sum(ObjectClass %in% class(IniStates)) < 2) {
215
      stop(paste0("Non convenient 'IniStates' for this 'FUN_MOD'\n"))
216
217
    }
    if (identical(FUN_MOD, RunModel_GR1A) & !is.null(IniStates)) { ## GR1A
218
      stop(paste0("'IniStates' is not available for this 'FUN_MOD'\n"))
219
220
    }
    if ((identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J)) & !all(is.na(IniStates$UH$UH1))) { ## GR5J
221
      stop(paste0("Non convenient IniStates for this 'FUN_MOD.' In 'IniStates', UH1 has to be a vector of NA for GR5J"))
222
223
    }
    if ((identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) & is.na(IniStates$Store$Exp)) { ## GR6J
224
      stop(paste0("Non convenient IniStates for this 'FUN_MOD.' GR6J needs an exponential store value in 'IniStates'"))
225
226
    }
    if (!(identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) & !is.na(IniStates$Store$Exp)) { ## except GR6J
227
      stop(paste0("Non convenient IniStates for this 'FUN_MOD.' No exponential store value needed in 'IniStates'"))
228
229
    }
    # if (length(na.omit(unlist(IniStates))) != NState) {
230
    #   stop(paste0("The length of IniStates must be ", NState, " for the chosen FUN_MOD"))
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
    # }
    if (!"CemaNeige" %in% ObjectClass & any(is.na(IniStates$CemaNeigeLayers$G  ))) {
      IniStates$CemaNeigeLayers$G   <- NULL
    }
    if (!"CemaNeige" %in% ObjectClass & any(is.na(IniStates$CemaNeigeLayers$eTG))) {
      IniStates$CemaNeigeLayers$eTG <- NULL
    }
    IniStates$Store$Rest <- rep(NA, 4)
    IniStates <- unlist(IniStates)
    IniStates[is.na(IniStates)] <- 0
    if ("monthly" %in% ObjectClass) {
      IniStates <- IniStates[seq_len(NState)]
    }
  } else {
    IniStates <- as.double(rep(0.0, NState))
  }
  
  
Delaigue Olivier's avatar
Delaigue Olivier committed
249
  ##check_Outputs_Cal_and_Sim
250
251
252
253
  
  ##Outputs_all
  Outputs_all <- NULL
  if (identical(FUN_MOD,RunModel_GR4H)) {
254
    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR4H")$GR)
255
256
  }
  if (identical(FUN_MOD,RunModel_GR4J) | identical(FUN_MOD,RunModel_CemaNeigeGR4J)) {
257
    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR4J")$GR)
258
259
  }
  if (identical(FUN_MOD,RunModel_GR5J) | identical(FUN_MOD,RunModel_CemaNeigeGR5J)) {
260
    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR5J")$GR)
261
262
  }
  if (identical(FUN_MOD,RunModel_GR6J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)) {
263
    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR6J")$GR)
264
265
  }
  if (identical(FUN_MOD,RunModel_GR2M)) {
266
    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR2M")$GR)
267
268
  }
  if (identical(FUN_MOD,RunModel_GR1A)) {
269
    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR1A")$GR)
270
271
  }
  if ("CemaNeige" %in% ObjectClass) {
272
    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = NULL, isCN = TRUE)$CN)
273
274
275
276
  }
  
  ##check_Outputs_Sim
  if (!is.vector(Outputs_Sim)) {
277
    stop("Outputs_Sim must be a vector of characters")
278
279
  }
  if (!is.character(Outputs_Sim)) {
280
    stop("Outputs_Sim must be a vector of characters")
281
282
  }
  if (sum(is.na(Outputs_Sim)) != 0) {
283
    stop("Outputs_Sim must not contain NA")
284
285
286
287
  }
  if ("all" %in% Outputs_Sim) {
    Outputs_Sim <- c("DatesR", Outputs_all, "StateEnd")
  }
288
  Test <- which(!Outputs_Sim %in% c("DatesR", Outputs_all, "StateEnd"))
289
  if (length(Test) != 0) {
290
    stop(paste0( "'Outputs_Sim' is incorrectly defined: ",
291
                 paste(Outputs_Sim[Test], collapse = ", "), " not found"))
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
  }
  Outputs_Sim <- Outputs_Sim[!duplicated(Outputs_Sim)]
  
  
  ##check_Outputs_Cal
  if (is.null(Outputs_Cal)) {
    if ("GR" %in% ObjectClass) {
      Outputs_Cal <- c("Qsim")
    }
    if ("CemaNeige" %in% ObjectClass) {
      Outputs_Cal <- c("all")
    }
    if ("GR" %in% ObjectClass &
        "CemaNeige" %in% ObjectClass) {
      Outputs_Cal <- c("PliqAndMelt", "Qsim")
    }
  } else {
    if (!is.vector(Outputs_Cal)) {
310
      stop("'Outputs_Cal' must be a vector of characters")
311
312
    }
    if (!is.character(Outputs_Cal)) {
313
      stop("'Outputs_Cal' must be a vector of characters")
314
315
    }
    if (sum(is.na(Outputs_Cal)) != 0) {
316
      stop("'Outputs_Cal' must not contain NA")
317
318
319
320
321
322
    }
  }
  if ("all" %in% Outputs_Cal) {
    Outputs_Cal <- c("DatesR", Outputs_all, "StateEnd")
    
  }
323
  Test <- which(!Outputs_Cal %in% c("DatesR", Outputs_all, "StateEnd"))
324
  if (length(Test) != 0) {
325
    stop(paste0("'Outputs_Cal' is incorrectly defined: ",
326
                paste(Outputs_Cal[Test], collapse = ", "), " not found"))
327
  }
328
  Outputs_Cal <- unique(Outputs_Cal)
329
330
331
332
333
334
335
336
337
338
339
  
  
  ##check_MeanAnSolidPrecip
  if ("CemaNeige" %in% ObjectClass & is.null(MeanAnSolidPrecip)) {
    NLayers <- length(InputsModel$LayerPrecip)
    SolidPrecip <- NULL
    for (iLayer in 1:NLayers) {
      if (iLayer == 1) {
        SolidPrecip <-
          InputsModel$LayerFracSolidPrecip[[1]] * InputsModel$LayerPrecip[[iLayer]] /
          NLayers
Delaigue Olivier's avatar
Delaigue Olivier committed
340
      } else {
341
        SolidPrecip <- SolidPrecip + InputsModel$LayerFracSolidPrecip[[iLayer]] * InputsModel$LayerPrecip[[iLayer]] / NLayers
Delaigue Olivier's avatar
Delaigue Olivier committed
342
343
      }
    }
344
345
346
347
348
349
350
351
352
353
354
355
356
357
    Factor <- NULL
    if (inherits(InputsModel, "hourly")) {
      Factor <- 365.25 * 24
    }
    if (inherits(InputsModel, "daily")) {
      Factor <- 365.25
    }
    if (inherits(InputsModel, "monthly")) {
      Factor <- 12
    }
    if (inherits(InputsModel, "yearly")) {
      Factor <- 1
    }
    if (is.null(Factor)) {
358
      stop("'InputsModel' must be of class 'hourly', 'daily', 'monthly' or 'yearly'")
359
360
361
    }
    MeanAnSolidPrecip <- rep(mean(SolidPrecip) * Factor, NLayers)
    ### default value: same Gseuil for all layers
362
    if (warnings) {
363
364
365
      warning("\t 'MeanAnSolidPrecip' not defined -> it was automatically set to c(",
              paste(round(MeanAnSolidPrecip), collapse = ","), ")  from the 'InputsModel' given to the function. ",
              "Be careful in case your application is (short-term) forecasting.\n")
366
367
368
369
    }
  }
  if ("CemaNeige" %in% ObjectClass & !is.null(MeanAnSolidPrecip)) {
    if (!is.vector(MeanAnSolidPrecip)) {
370
      stop(paste0("'MeanAnSolidPrecip' must be a vector of numeric values"))
371
372
    }
    if (!is.numeric(MeanAnSolidPrecip)) {
373
      stop(paste0("'MeanAnSolidPrecip' must be a vector of numeric values"))
374
375
    }
    if (length(MeanAnSolidPrecip) != NLayers) {
376
      stop(paste0("'MeanAnSolidPrecip' must be a numeric vector of length ", NLayers, ""))
377
378
379
380
    }
  }
  
  
Delaigue Olivier's avatar
Delaigue Olivier committed
381
  ##check_PliqAndMelt
382
  if ("GR" %in% ObjectClass & "CemaNeige" %in% ObjectClass) {
383
    if (!"PliqAndMelt" %in% Outputs_Cal & !"all" %in% Outputs_Cal) {
384
      WTxt <- NULL
385
      WTxt <- paste0(WTxt, "\t 'PliqAndMelt' was not defined in 'Outputs_Cal' but is needed to feed the hydrological model with the snow modele outputs \n")
386
      WTxt <- paste0(WTxt, "\t -> it was automatically added \n")
387
      if (!is.null(WTxt) & warnings) {
388
389
390
        warning(WTxt)
      }
      Outputs_Cal <- c(Outputs_Cal, "PliqAndMelt")
Delaigue Olivier's avatar
Delaigue Olivier committed
391
    }
392
    if (!"PliqAndMelt" %in% Outputs_Sim & !"all" %in% Outputs_Sim) {
393
      WTxt <- NULL
394
      WTxt <- paste0(WTxt, "\t 'PliqAndMelt' was not defined in 'Outputs_Sim' but is needed to feed the hydrological model with the snow modele outputs \n")
395
      WTxt <- paste0(WTxt, "\t -> it was automatically added \n")
396
      if (!is.null(WTxt) & warnings) {
397
398
399
400
401
402
403
        warning(WTxt)
      }
      Outputs_Sim <- c(Outputs_Sim, "PliqAndMelt")
    }
  }
  
  
Delaigue Olivier's avatar
Delaigue Olivier committed
404
  ##Create_RunOptions
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
  RunOptions <- list(IndPeriod_WarmUp = IndPeriod_WarmUp,
                     IndPeriod_Run = IndPeriod_Run,
                     IniStates = IniStates,
                     IniResLevels = IniResLevels,
                     Outputs_Cal = Outputs_Cal,
                     Outputs_Sim = Outputs_Sim)
  
  if ("CemaNeige" %in% ObjectClass) {
    RunOptions <-
      c(RunOptions, list(MeanAnSolidPrecip = MeanAnSolidPrecip))
  }
  class(RunOptions) <- c("RunOptions", ObjectClass)
  
  return(RunOptions)
  
  
Delaigue Olivier's avatar
Delaigue Olivier committed
421
422
}