RunModel_Lag.R 3.74 KB
Newer Older
1
RunModel_Lag <- function(InputsModel, RunOptions, Param) {
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
20
  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"))
  }
  if (is.null(InputsModel$OutputsModel)) {
21
22
23
    stop(
      "'InputsModel' should contain an 'OutputsModel' key containing the output of the runoff of the downstream subcatchment"
    )
24
  }
25
  if (is.null(InputsModel$OutputsModel$Qsim)) {
26
27
28
    stop(
      "'InputsModel$OutputsModel' should contain a key 'Qsim' containing the output of the runoff of the downstream subcatchment"
    )
29
  }
30
  if (sum(!is.na(InputsModel$OutputsModel$Qsim)) != length(RunOptions$IndPeriod_Run)) {
31
32
33
    stop(
      "'InputsModel$OutputsModel$Qim' should have the same lenght as 'RunOptions$IndPeriod_Run' and contain no NA"
    )
34
  }
35

36
37
  OutputsModel <- InputsModel$OutputsModel
  OutputsModel$QsimDown <- OutputsModel$Qsim
38

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

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

  NbUpBasins <- length(InputsModel$LengthHydro)
  LengthTs <- length(OutputsModel$QsimDown)
53
  OutputsModel$Qsim <- OutputsModel$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
54

55
  IniSD <- RunOptions$IniStates[grep("SD", names(RunOptions$IniStates))]
56
57
58
59
60
61
62
63
64
65
  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
        )
      )
    }
66
67
    IniStates <- lapply(seq_len(NbUpBasins), function(x) {
      iStart <- 1
68
69
      if (x > 1) {
        iStart <- iStart + sum(floor(PT[1:x - 1]) + 1)
70
71
72
73
74
      }
      IniSD[iStart:(iStart + PT[x])]
    })
  } else {
    IniStates <- lapply(seq_len(NbUpBasins), function(x) {
75
76
      rep(0, floor(PT[x] + 1))
    })
77
  }
78

79
  for (upstream_basin in seq_len(NbUpBasins)) {
80
    Qupstream <- InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin]
81
    if (!is.na(InputsModel$BasinAreas[upstream_basin])) {
82
      # Upstream flow with area needs to be converted to m3 by time step
83
      Qupstream <- Qupstream * InputsModel$BasinAreas[upstream_basin] * 1e3
84
85
    }
    OutputsModel$Qsim <- OutputsModel$Qsim +
86
      c(IniStates[[upstream_basin]][-length(IniStates[[upstream_basin]])],
87
88
        Qupstream[1:(LengthTs - floor(PT[upstream_basin]))]) *
      HUTRANS[1, upstream_basin] +
89
      c(IniStates[[upstream_basin]],
90
91
92
93
        Qupstream[1:(LengthTs - floor(PT[upstream_basin]) - 1)]) *
      HUTRANS[2, upstream_basin]
  }
  # Warning for negative flows
94
  if (any(OutputsModel$Qsim < 0)) {
95
96
97
98
    warning(length(which(OutputsModel$Qsim < 0)), " time steps with negative flow, set to zero.")
    OutputsModel$Qsim[OutputsModel$Qsim < 0] <- 0
  }
  # Convert back Qsim to mm
99
  OutputsModel$Qsim <- OutputsModel$Qsim / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
100

101
  if ("StateEnd" %in% RunOptions$Outputs_Sim) {
102
    OutputsModel$StateEnd$SD <- lapply(seq(NbUpBasins), function(x) {
103
104
      Qupstream[(LengthTs - floor(PT[x])):LengthTs]
    })
105
106
  }

107
  return(OutputsModel)
108
}