Failed to fetch fork details. Try again later.
-
Dorchies David authored
Refs HYCAR-Hydro/airgr#152, HYCAR-Hydro/airgr#153
057d235c
Forked from
HYCAR-Hydro / airGR
Source project has a limited visibility.
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)
}