Skip to content
GitLab
Select Git revision
Find file BlameHistoryPermalink
RunModel_CemaNeigeGR4J.R 9.53 KiB
1 
RunModel_CemaNeigeGR4J <- function(InputsModel, RunOptions, Param) {
2
3
4 
  ## Initialization of variables
5 
  IsHyst <- inherits(RunOptions, "hysteresis")
6 
  NParamCN <- RunOptions$FeatFUN_MOD$NbParam - 4L
7 
  NStates <- 4L
8
9
10 
  .ArgumentsCheckGR(InputsModel, RunOptions, Param)
11
12 
  Param <- as.double(Param)
13
14
15 
  Param_X1X3_threshold <- 1e-2
16 
  Param_X4_threshold   <- 0.5
17 
  if (Param[1L] < Param_X1X3_threshold) {
18 
    warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
19 
    Param[1L] <- Param_X1X3_threshold
20 
  }
21 
  if (Param[3L] < Param_X1X3_threshold) {
22 
    warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
23 
    Param[3L] <- Param_X1X3_threshold
24 
  }
25 
  if (Param[4L] < Param_X4_threshold) {
26 
    warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
27 
    Param[4L] <- Param_X4_threshold
28 
  }
29
30 
  ## Input data preparation
31 
  if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
32 
    RunOptions$IndPeriod_WarmUp <- NULL
33 
  }
34 
  IndPeriod1     <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
35 
  LInputSeries   <- as.integer(length(IndPeriod1))
36 
  IndPeriod2     <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
37 
  ParamCemaNeige <- Param[(length(Param) - 1 - 2 * as.integer(IsHyst)):length(Param)]
38 
  NParamMod      <- as.integer(length(Param) - (2 + 2 * as.integer(IsHyst)))
39 
  ParamMod       <- Param[1:NParamMod]
40 
  NLayers        <- length(InputsModel$LayerPrecip)
41 
  NStatesMod     <- as.integer(length(RunOptions$IniStates) - NStates * NLayers)
42
43 
  ## Output data preparation
44 
  ExportDatesR   <- "DatesR"   %in% RunOptions$Outputs_Sim
45 
  ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
46
47
48 
  ## CemaNeige________________________________________________________________________________
49 
  if (inherits(RunOptions, "CemaNeige")) {
50 
    if ("all" %in% RunOptions$Outputs_Sim) {
51 
      IndOutputsCemaNeige <- as.integer(1:length(RunOptions$FortranOutputs$CN))
52 
    } else {
53 
      IndOutputsCemaNeige <- which(RunOptions$FortranOutputs$CN %in% RunOptions$Outputs_Sim)
54 
    }
55 
    CemaNeigeLayers <- list()
56 
    CemaNeigeStateEnd <- NULL
57 
    NameCemaNeigeLayers <- "CemaNeigeLayers"
58
59
60 
    ## Call CemaNeige Fortran_________________________
61 
    for (iLayer in 1:NLayers) {
62 
      if (!IsHyst) {
63 
        StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers)]
64 
      } else {
65 
        StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
66 
      }
67 
      RESULTS <- .Fortran("frun_cemaneige", PACKAGE = "airGR",
68 
                          ## inputs
69 
                          LInputs = LInputSeries,                                                         ### length of input and output series
70 
                          InputsPrecip = InputsModel$LayerPrecip[[iLayer]][IndPeriod1],                   ### input series of total precipitation [mm/d]