An error occurred while loading the file. Please try again.
-
Grand Francois authoreda8566117
RunModel_LRQ <- function(InputsModel, RunOptions, Param, QcontribDown) {
NParam <- 2
#RunModel_LRQ
## 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
CParam <- Param[1L] # CParam in m/s
KParam <- Param[2L] # KParam in mm/TS
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
PT <- InputsModel$LengthHydro * 1e3 / CParam / RunOptions$FeatFUN_MOD$TimeStep #PT in h
PK <- sqrt(KParam * InputsModel$BasinAreas[-length(InputsModel$BasinAreas)] * 1e-3) * PT #m3/2 pdt1/2
## Lag set up
Lag <- function(t, Qroute, Qupstream){
Q <- c(Qupstream[1:(1 + floor(t))], Qroute[1:(length(Qroute) - (floor(t) + 1))])
return(Q)
}
## set up initial states
IniSD <- RunOptions$IniStates[grep("SD", names(RunOptions$IniStates))]
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
)
)
}
IniStates <- lapply(seq_len(NbUpBasins), function(x) {
iStart <- 1
if (x > 1) {
iStart <- iStart + sum(floor(PT[1:x - 1]) + 1)
}
as.vector(IniSD[iStart:(iStart + PT[x])])
})
} else {
IniStates <- lapply(
seq_len(NbUpBasins),
function(iUpBasins) {
iWarmUp <- seq.int(
from = max(1, IndPeriod1[1] - floor(PT[iUpBasins]) - 1),
to = max(1, IndPeriod1[1] - 1)
)
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(as.vector(ini))
}
)
}
## Lag & route 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 <- c(IniStates[[upstream_basin]],
InputsModel$Qupstream[IndPeriod1, upstream_basin])
Qroute <- Qupstream_m3
for (q in seq_along(Qroute)[1:(length(Qroute) - 1)]){
Qroute[q + 1] <- ((Qroute[q]^0.5 + ((Qupstream_m3[q] + Qupstream_m3[q + 1])/2)^0.5 *
tanh((((Qupstream_m3[q] + Qupstream_m3[q + 1])/2)^0.5)/PK[upstream_basin]) ) /
(1 + (Qroute[q]/((Qupstream_m3[q] + Qupstream_m3[q + 1])/2))^0.5 *
tanh((((Qupstream_m3[q] + Qupstream_m3[q + 1])/2)^0.5)/PK[upstream_basin])) )^2
}
Qllr <- Lag(PT[upstream_basin], Qroute, Qupstream_m3)[IndPeriod1]
Qsim_m3[, upstream_basin] <- Qllr
}
Qsim_m3 <- cbind(QsimDown_input, Qsim_m3)
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196
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_LLR, InputsModel, SD = SD)
} else {
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)
}