diff --git a/R/RunModel_LLR.R b/R/RunModel_LLR.R index c11b12d70bbcf80894105c161ef682eb2250c9db..a37601b82603d15dd3c012fcdd358a549858f892 100644 --- a/R/RunModel_LLR.R +++ b/R/RunModel_LLR.R @@ -68,19 +68,19 @@ RunModel_LLR <- function(InputsModel, RunOptions, Param, QcontribDown) { TParam <- Param[1L] KParam <- Param[2L] - ## propagation time from upstream meshes to outlet - PT <- ifelse(InputsModel$LengthHydro == 0, 0, InputsModel$LengthHydro/max(InputsModel$LengthHydro)) * - (TParam + KParam) - ifelse(InputsModel$LengthHydro == 0, 0, - sqrt(InputsModel$LengthHydro/max(InputsModel$LengthHydro))) * KParam - PT[PT < 0] <- 0 - PK <- ifelse(InputsModel$LengthHydro == 0, - 0, - sqrt(InputsModel$LengthHydro/max(InputsModel$LengthHydro)) * KParam) + PT <- InputsModel$LengthHydro * 1e3 / Param[1L] / RunOptions$FeatFUN_MOD$TimeStep + PK <- KParam * PT C0 <- ifelse(PK ==0, 0, exp(-1/PK)) C1 <- ((PK/1) * (1-C0))-C0 C2 <- 1 - (PK/1) * (1-C0) + ## Lag set up + + Lag <- function(t, Qroute, Qupstream){ + Q <- c(Qupstream[1 : floor(t)], Qroute[1 : (length(Qroute) - floor(t))]) + return(Q) + } ## set up initial states IniSD <- RunOptions$IniStates[grep("SD", names(RunOptions$IniStates))] if (length(IniSD) > 0) { @@ -135,10 +135,9 @@ RunModel_LLR <- function(InputsModel, RunOptions, Param, QcontribDown) { 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)))] + Qllr <- Lag(PT[upstream_basin], Qroute, Qupstream_m3)[IndPeriod1] + Qsim_m3[, upstream_basin] <- Qllr } Qsim_m3 <- cbind(QsimDown_input, Qsim_m3)