diff --git a/NAMESPACE b/NAMESPACE
index 8348ba43e310606ba7183ff78740ef1d789dbab3..7769adf8b96150e8c1790d81bc2aabe724578fcb 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -55,6 +55,7 @@ export(RunModel_GR5J)
export(RunModel_GR6J)
export(RunModel_Lag)
export(RunModel_LLR)
+export(RunModel_LBLR)
export(SeriesAggreg)
export(TransfoParam)
export(TransfoParam_CemaNeige)
@@ -68,6 +69,7 @@ export(TransfoParam_GR5J)
export(TransfoParam_GR6J)
export(TransfoParam_Lag)
export(TransfoParam_LLR)
+export(TransfoParam_LBLR)
export(.ErrorCrit)
export(.FeatModels)
diff --git a/R/CreateCalibOptions.R b/R/CreateCalibOptions.R
index 9102561e3700abb14637a8546f1a4db45b3bb8b8..8c2cb2e0539d0a78878f7bc72048e0f93e2f598d 100644
--- a/R/CreateCalibOptions.R
+++ b/R/CreateCalibOptions.R
@@ -207,10 +207,16 @@ CreateCalibOptions <- function(FUN_MOD,
+5.00), ncol = 1, byrow = TRUE)
ParamT <- cbind(ParamTSD, ParamT)
}
- if (identical(RunModel_LLR, FUN_SD)){
- ParamTSD <- matrix(c(+1.25, +1.25,
- +2.50, +2.50,
- +5.00, +5.00), ncol = 2, byrow = TRUE)
+ else if (identical(RunModel_LLR, FUN_SD)){
+ ParamTSD <- matrix(c(+1.25, 0,
+ +2.50, 250,
+ +5.00, 500), ncol = 2, byrow = TRUE)
+ ParamT <- cbind(ParamTSD, ParamT)
+ }
+ else if (identical(RunModel_LBLR, FUN_SD)){
+ ParamTSD <- matrix(c(+1.25, 0,
+ +2.50, 250,
+ +5.00, 500), ncol = 2, byrow = TRUE)
ParamT <- cbind(ParamTSD, ParamT)
}
}
diff --git a/R/RunModel_LBLR.R b/R/RunModel_LBLR.R
new file mode 100644
index 0000000000000000000000000000000000000000..9d483ee70cbcdd46d7e48011bf2d28be13693e3f
--- /dev/null
+++ b/R/RunModel_LBLR.R
@@ -0,0 +1,199 @@
+RunModel_LBLR <- 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]
+
+
+ PK <- ifelse(InputsModel$LengthHydro == 0, 0,
+ sqrt(InputsModel$LengthHydro/max(InputsModel$LengthHydro)) * KParam)
+
+ PT <- ifelse(InputsModel$LengthHydro == 0, 0,
+ InputsModel$LengthHydro/max(InputsModel$LengthHydro)) * (TParam + KParam) - PK
+ PT <- ifelse(PT < 0, 0, PT)
+ rprop <- PT - floor(PT)
+
+ ## Lag set up
+ first_lag <- function(t, Qroute, Qupstream){
+ Q <- c(Qupstream[1 : (floor(t) + 1)], Qroute[1 : (length(Qroute) - (floor(t) + 1))])
+ return(Q)
+ }
+ scnd_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) {
+ 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)))
+ }
+ 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_upstream in seq_along(Qupstream_m3)[2:max(seq_along(Qupstream_m3))]){
+ Qroute[q_upstream] <- (1 - 1/PK[upstream_basin]) * Qroute[q_upstream - 1] +
+ 1/ PK[upstream_basin] * Qupstream_m3[q_upstream]
+ }
+ Qsim_m3[, upstream_basin] <- first_lag(PT[upstream_basin], Qroute, Qupstream_m3)[IndPeriod1] *
+ rprop[upstream_basin] + scnd_lag(PT[upstream_basin], Qroute, Qupstream_m3)[IndPeriod1] * (1 - rprop[upstream_basin])
+ }
+ 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_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)
+}
diff --git a/R/TransfoParam_LBLR.R b/R/TransfoParam_LBLR.R
new file mode 100644
index 0000000000000000000000000000000000000000..e92baf75630348732ac5114ded93ff9003472c76
--- /dev/null
+++ b/R/TransfoParam_LBLR.R
@@ -0,0 +1,42 @@
+TransfoParam_LBLR <- function(ParamIn, Direction) {
+
+ ## number of model parameters
+ NParam <- 2L
+
+ ## check arguments
+ isVecParamIn <- is.vector(ParamIn)
+ if (isVecParamIn) {
+ ParamIn <- matrix(ParamIn, nrow = 1)
+ }
+ if (!inherits(ParamIn, "matrix")) {
+ stop("'ParamIn' must be of class 'matrix'")
+ }
+ if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
+ stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
+ }
+ if (ncol(ParamIn) != NParam) {
+ stop(sprintf("the Lag & Route model requires %i parameters", NParam))
+ }
+
+
+ ## transformation
+ if (Direction == "TR") {
+ ParamOut <- ParamIn
+ ParamOut[, 1] <- 20 * (ParamIn[, 1] + 10) / 20.0
+ ParamOut[, 2] <- 20 * (ParamIn[, 2] + 10) / 20.0
+ }
+ if (Direction == "RT") {
+ ParamOut <- ParamIn
+ ParamOut[, 1] <- ParamIn[, 1] * 20.0 / 20 - 10
+ ParamOut[, 2] <- ParamIn[, 2] * 20.0 / 20 - 10
+ }
+
+
+ ## check output
+ if (isVecParamIn) {
+ ParamOut <- as.vector(ParamOut)
+ }
+
+ return(ParamOut)
+
+}
diff --git a/R/Utils.R b/R/Utils.R
index 87b91ce17a5458690d7fdd290f8ab39c8c265fa0..e6cfa5ee255307c6eb258b8597c717abb8ae877c 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -80,6 +80,9 @@
}else if(identical(RunModel_LLR, FUN_SD)){
res$NParamSD <- FeatMod$NbParam[15]
res$FUN_SD <- "LLR"
+ }else if(identical(RunModel_LBLR, FUN_SD)){
+ res$NParamSD <- FeatMod$NbParam[16]
+ res$FUN_SD <- "LBLR"
}
}
return(res)
diff --git a/R/UtilsCalibOptions.R b/R/UtilsCalibOptions.R
index 38d1c0337f07dde1676b0c10b9bc187580d223cc..010271e2d2a48eb9f1e908891465078e71bb8d8b 100644
--- a/R/UtilsCalibOptions.R
+++ b/R/UtilsCalibOptions.R
@@ -31,6 +31,8 @@
FUN_ROUT <- TransfoParam_Lag
} else if (identical(RunModel_LLR, FUN_SD)){
FUN_ROUT <- TransfoParam_LLR
+ } else if (identical(RunModel_LBLR, FUN_SD)){
+ FUN_ROUT <- TransfoParam_LBLR
}
}
diff --git a/inst/modelsFeatures/FeatModelsGR.csv b/inst/modelsFeatures/FeatModelsGR.csv
index d23cc8e061070307219924d82aa8e0d5297c7d23..71d2324046ab3edb6a9d7255be14252da95be258 100644
--- a/inst/modelsFeatures/FeatModelsGR.csv
+++ b/inst/modelsFeatures/FeatModelsGR.csv
@@ -14,3 +14,4 @@ CemaNeigeGR4H;CemaNeigeGR4H;6;hourly;NA;GR;airGR
CemaNeigeGR5H;CemaNeigeGR5H;7;hourly;NA;GR;airGR
Lag;Lag;1;NA;NA;SD;airGR
LLR;LLR;2;NA;NA;SD;airGR
+LBLR;LBLR;2;NA;NA;SD;airGR