RunModel_Lag.R 4.66 KB
Newer Older
1
RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
2
  NParam <- 1
3

4
  ##Arguments_check
5
6
  if (!inherits(InputsModel, "InputsModel")) {
    stop("'InputsModel' must be of class 'InputsModel'")
7
  }
8
9
  if (!inherits(InputsModel, "SD")) {
    stop("'InputsModel' must be of class 'SD'")
10
  }
11
12
13
14
15
16
17
18
19
  if (!inherits(RunOptions, "RunOptions")) {
    stop("'RunOptions' must be of class 'RunOptions'")
  }
  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"))
  }
20
21
22
23
24
25
26
27
28
29
30
31
  if (inherits(QcontribDown, "OutputsModel")) {
    if (is.null(QcontribDown$Qsim)) {
      stop("'QcontribDown' should contain a key 'Qsim' containing the output of the runoff of the downstream subcatchment")
    }
    OutputsModel <- QcontribDown
    OutputsModel$QsimDown <- OutputsModel$Qsim
  } else if (is.vector(QcontribDown) && is.numeric(QcontribDown)) {
    OutputsModel <- list()
    class(OutputsModel) <- c("OutputsModel", class(OutputsModel))
    OutputsModel$QsimDown <- QcontribDown
  } else {
    stop("'QcontribDown' must be a numeric vector or a 'OutputsModel' object")
32
  }
33
34
  if (length(OutputsModel$QsimDown) != length(RunOptions$IndPeriod_Run)) {
    stop("Time series in 'QcontribDown' should have the same lenght as 'RunOptions$IndPeriod_Run'")
35
  }
36

37
38
  if (inherits(InputsModel, "hourly")) {
    TimeStep <- 60 * 60
39
40
41
42
  } else if (inherits(InputsModel, "daily")) {
    TimeStep <- 60 * 60 * 24
  } else {
    stop("'InputsModel' should be of class \"daily\" or \"hourly\"")
43
44
45
  }

  # propagation time from upstream meshes to outlet
46
  PT <- InputsModel$LengthHydro * 1e3 / Param[1L] / TimeStep
47
48
49
50
  HUTRANS <- rbind(1 - (PT - floor(PT)), PT - floor(PT))

  NbUpBasins <- length(InputsModel$LengthHydro)
  LengthTs <- length(OutputsModel$QsimDown)
51
  OutputsModel$Qsim_m3 <- OutputsModel$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
52

53
  IniSD <- RunOptions$IniStates[grep("SD", names(RunOptions$IniStates))]
54
55
56
57
58
59
60
61
62
63
  if (length(IniSD) > 0) {
    if (sum(floor(PT)) + NbUpBasins != length(IniSD)) {
      stop(
        sprintf(
          "SD initial states has a length of %i and a length of %i is required",
          length(IniSD),
          sum(floor(PT)) + NbUpBasins
        )
      )
    }
64
65
    IniStates <- lapply(seq_len(NbUpBasins), function(x) {
      iStart <- 1
66
67
      if (x > 1) {
        iStart <- iStart + sum(floor(PT[1:x - 1]) + 1)
68
69
70
71
      }
      IniSD[iStart:(iStart + PT[x])]
    })
  } else {
72
73
74
75
    IniStates <- lapply(
      seq_len(NbUpBasins),
      function(iUpBasins) {
        iWarmUp <- seq.int(
76
77
          from = max(1, RunOptions$IndPeriod_WarmUp[length(RunOptions$IndPeriod_WarmUp)] - floor(PT[iUpBasins])),
          to   = max(1, RunOptions$IndPeriod_WarmUp[length(RunOptions$IndPeriod_WarmUp)])
78
79
80
81
82
83
84
85
86
        )
        ini <- InputsModel$Qupstream[iWarmUp, iUpBasins]
        if(length(ini) != floor(PT[iUpBasins] + 1)) {
          # If warm-up period is not enough long complete beginning with first value
          ini <- c(rep(ini[1], floor(PT[iUpBasins] + 1) - length(ini)), ini)
        }
        return(ini)
      }
    )
87
  }
88
  # message("Initstates: ", paste(IniStates, collapse = ", "))
89

90
  for (upstream_basin in seq_len(NbUpBasins)) {
91
92
    Qupstream <- c(IniStates[[upstream_basin]],
                   InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin])
93
    # message("Qupstream[", upstream_basin, "]: ", paste(Qupstream, collapse = ", "))
94
    OutputsModel$Qsim_m3 <- OutputsModel$Qsim_m3 +
95
96
      Qupstream[2:(1 + LengthTs)] * HUTRANS[1, upstream_basin] +
      Qupstream[1:LengthTs] * HUTRANS[2, upstream_basin]
97
  }
98
99
100
101
  # Convert back Qsim to mm
  OutputsModel$Qsim <- OutputsModel$Qsim_m3 / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
  # message("Qsim: ", paste(OutputsModel$Qsim, collapse = ", "))

102
  # Warning for negative flows or NAs only in extended outputs
103
  if(length(RunOptions$Outputs_Sim) > 2) {
104
105
106
107
108
109
110
111
    if (any(OutputsModel$Qsim[!is.na(OutputsModel$Qsim)] < 0)) {
      warning(length(which(OutputsModel$Qsim < 0)), " time steps with negative flow, set to zero.")
      OutputsModel$Qsim[OutputsModel$Qsim < 0] <- 0
    }
    # Warning for NAs
    if (any(is.na(OutputsModel$Qsim))) {
      warning(length(which(is.na(OutputsModel$Qsim))), " time steps with NA values")
    }
112
  }
113

114
  if ("StateEnd" %in% RunOptions$Outputs_Sim) {
115
    OutputsModel$StateEnd$SD <- lapply(seq(NbUpBasins), function(x) {
116
117
      lastTS <- RunOptions$IndPeriod_Run[length(RunOptions$IndPeriod_Run)]
      InputsModel$Qupstream[(lastTS - floor(PT[x])):lastTS, x]
118
    })
119
    # message("StateEnd: ", paste(OutputsModel$StateEnd$SD, collapse = ", "))
120
121
  }

122
  return(OutputsModel)
123
}