RunModel.R

RunModel <- function (InputsModel, RunOptions, Param, FUN_MOD) {

  FUN_MOD <- match.fun(FUN_MOD)

  if (inherits(InputsModel, "SD")) {
    OutputsModelDown <- FUN_MOD(InputsModel = InputsModel, RunOptions = RunOptions,
                                Param = Param[-length(Param)])
    OutputsModelDown$QsimDown <- OutputsModelDown$Qsim

    if (inherits(InputsModel, "daily")) {
      TimeStep <- 60 * 60 * 24
    }
    if (inherits(InputsModel, "hourly")) {
      TimeStep <- 60 * 60
    }

    # propagation time from upstream meshes to outlet
    PT <- InputsModel$LengthHydro / Param[length(Param)] / TimeStep
    HUTRANS <- rbind(1 - (PT - floor(PT)), PT - floor(PT))

    NbUpBasins <- length(InputsModel$LengthHydro)
    LengthTs <- length(OutputsModelDown$QsimDown)
    OutputsModelDown$Qsim <- OutputsModelDown$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1E3

    for (upstream_basin in seq_len(NbUpBasins)) {
      Qupstream <- InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin]
      if(!is.na(InputsModel$BasinAreas[upstream_basin])) {
        # Upstream flow with area needs to be converted to m3 by time step
        Qupstream <- Qupstream * InputsModel$BasinAreas[upstream_basin] * 1E3
      }
      OutputsModelDown$Qsim <- OutputsModelDown$Qsim +
        c(rep(0, floor(PT[upstream_basin])),
          Qupstream[1:(LengthTs - floor(PT[upstream_basin]))]) *
        HUTRANS[1, upstream_basin] +
        c(rep(0, floor(PT[upstream_basin] + 1)),
          Qupstream[1:(LengthTs - floor(PT[upstream_basin]) - 1)]) *
        HUTRANS[2, upstream_basin]
    }
    # Warning for negative flows
    if(any(OutputsModelDown$Qsim < 0)) {
      warning(length(which(OutputsModelDown$Qsim < 0)), " time steps with negative flow, set to zero.")
      OutputsModelDown$Qsim[OutputsModelDown$Qsim < 0] <- 0
    }
    # Convert back Qsim to mm
    OutputsModelDown$Qsim <- OutputsModelDown$Qsim / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1E3

  } else {

    OutputsModelDown <- FUN_MOD(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param)

  }
  return(OutputsModelDown)
}