Forked from HYCAR-Hydro / airGR
Source project has a limited visibility.
RunModel_LLR.R 6.36 KiB
RunModel_LLR <- function(InputsModel, RunOptions, Param, QcontribDown) {
  NParam <- 2
  ## argument check
  if (!inherits(InputsModel, "InputsModel")) {
    stop("'InputsModel' must be of class 'InputsModel'")
  if (!inherits(InputsModel, "SD")) {
    stop("'InputsModel' must be of class 'SD'")
  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 (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")
    if (length(QcontribDown$Qsim) != length(RunOptions$IndPeriod_Run)) {
      stop("Time series Qsim in 'QcontribDown' should have the same length as 'RunOptions$IndPeriod_Run'")
    if (!identical(RunOptions$IndPeriod_WarmUp, 0L) && !identical(RunOptions$Outputs_Sim, RunOptions$Outputs_Cal)) {
      # This test is not necessary during calibration but usefull in other cases because
      # WarmUpQsim is then used for downstream sub-basins because of the delay in Qupstream
      if (is.null(QcontribDown$RunOptions$WarmUpQsim) ||
          length(QcontribDown$RunOptions$WarmUpQsim) != length(RunOptions$IndPeriod_WarmUp)) {
        stop("Time series WarmUpQsim in 'QcontribDown' should have the same length as 'RunOptions$IndPeriod_WarmUp'")
  } else if (is.vector(QcontribDown) && is.numeric(QcontribDown)) {
    if (length(QcontribDown) != length(RunOptions$IndPeriod_Run)) {
      stop("'QcontribDown' should have the same length as 'RunOptions$IndPeriod_Run'")
  } else {
    stop("'QcontribDown' must be a numeric vector or a 'OutputsModel' object")
  # data set up
  NbUpBasins <- length(InputsModel$LengthHydro)
  if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
    RunOptions$IndPeriod_WarmUp <- NULL
  IndPeriod1   <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
  LInputSeries <- as.integer(length(IndPeriod1))
  IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
  if (inherits(QcontribDown, "OutputsModel")) {
    OutputsModel <- QcontribDown
    if (is.null(OutputsModel$RunOptions$WarmUpQsim)) {
      OutputsModel$RunOptions$WarmUpQsim <- rep(NA, length(RunOptions$IndPeriod_WarmUp))
    QsimDown <- c(OutputsModel$RunOptions$WarmUpQsim, OutputsModel$Qsim)
  } else if (is.vector(QcontribDown) && is.numeric(QcontribDown)) {
    OutputsModel <- list()
    class(OutputsModel) <- c("OutputsModel", class(RunOptions)[-1])
    QsimDown <- c(rep(NA, length(RunOptions$IndPeriod_WarmUp)),
                  QcontribDown)
  ## Parameters set up
  TParam <- Param[1L]
  KParam <- Param[2L]
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
## propagation time from upstream meshes to outlet PT <- floor(InputsModel$LengthHydro * 1e3 / TParam / RunOptions$FeatFUN_MOD$TimeStep) PK <- ifelse(InputsModel$LengthHydro == 0, 0, sqrt(InputsModel$LengthHydro/max(InputsModel$LengthHydro)) * KParam) C0 <- ifelse(PK ==0, 0, exp(-1/PK)) C1 <- ((PK/1) * (1-C0))-C0 C2 <- 1 - (PK/1) * (1-C0) ## Lag model computation Qsim_m3 <- matrix(0 , nrow = length(IndPeriod1), ncol = NbUpBasins) QsimDown_input <- matrix(QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3, ncol = 1) for (upstream_basin in seq_len(NbUpBasins)) { Qupstream_m3 <- InputsModel$Qupstream[IndPeriod1, upstream_basin] Qroute <- Qllr <- Qupstream_m3 if (is.na(Qroute[1])){ Qroute[1] <- 0 } Qroute_t1 <- Qroute[1] for (q_upstream in seq_along(Qupstream_m3)[2:max(seq_along(Qupstream_m3))]){ if (!is.na(Qroute[q_upstream - 1])){ Qroute_t1 <- Qroute[q_upstream - 1] } Qroute[q_upstream] <- C0[upstream_basin] * Qroute_t1 + C1[upstream_basin] * Qupstream_m3[q_upstream - 1] + C2[upstream_basin] * Qupstream_m3[q_upstream] Qllr[q_upstream + PT[upstream_basin]] <- Qroute[q_upstream] } Qsim_m3[, upstream_basin] <- Qllr[(1 +length(Qllr) - length(IndPeriod1)) : (length(IndPeriod1) + (length(Qllr) - length(IndPeriod1)))] } Qsim_m3 <- cbind(QsimDown_input, Qsim_m3) Qsim_m3 <- rowSums(Qsim_m3) ## OutputsModel if ("Qsim_m3" %in% RunOptions$Outputs_Sim) { OutputsModel$Qsim_m3 <- Qsim_m3[IndPeriod2] } if ("Qsim" %in% RunOptions$Outputs_Sim) { # Convert back Qsim to mm OutputsModel$Qsim <- Qsim_m3[IndPeriod2] / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3 # message("Qsim: ", paste(OutputsModel$Qsim, collapse = ", ")) } if ("QsimDown" %in% RunOptions$Outputs_Sim) { # Convert back Qsim to mm OutputsModel$QsimDown <- QsimDown[IndPeriod2] } # Warning for negative flows or NAs only in extended outputs if (length(RunOptions$Outputs_Sim) > 2) { 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") } } if ("StateEnd" %in% RunOptions$Outputs_Sim) { SD <- lapply(seq(NbUpBasins), function(x) { lastTS <- RunOptions$IndPeriod_Run[length(RunOptions$IndPeriod_Run)] InputsModel$Qupstream[(lastTS - floor(PT[x])):lastTS, x] }) if (is.null(OutputsModel$StateEnd)) { OutputsModel$StateEnd <- CreateIniStates(RunModel_Lag, InputsModel, SD = SD) } else {
141142143144145146147148149150151152153154155156157
OutputsModel$StateEnd$SD <- SD } # message("StateEnd: ", paste(OutputsModel$StateEnd$SD, collapse = ", ")) } if ("WarmUpQsim" %in% RunOptions$Outputs_Sim) { OutputsModel$RunOptions$WarmUpQsim <- Qsim_m3[seq_len(length(RunOptions$IndPeriod_WarmUp))] / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3 } if ("Param" %in% RunOptions$Outputs_Sim) { OutputsModel$RunOptions$Param <- c(Param, OutputsModel$RunOptions$Param) } class(OutputsModel) <- c(class(OutputsModel), "SD") return(OutputsModel) }