diff --git a/R/Calibration.R b/R/Calibration.R
index 319dbb4b2f9d083f64631e49dc1b65953012e434..e4fc1984fad39839755ebd90ee3f900b76c7a9e5 100644
--- a/R/Calibration.R
+++ b/R/Calibration.R
@@ -3,7 +3,6 @@ Calibration <- function(InputsModel,
InputsCrit,
CalibOptions,
FUN_MOD,
- FUN_SD,
FUN_CRIT, # deprecated
FUN_CALIB = Calibration_Michel,
FUN_TRANSFO = NULL,
@@ -23,7 +22,7 @@ Calibration <- function(InputsModel,
}
return(FUN_CALIB(InputsModel = InputsModel, RunOptions = RunOptions, InputsCrit = InputsCrit,
CalibOptions = CalibOptions,
- FUN_MOD = FUN_MOD, FUN_SD, FUN_TRANSFO = FUN_TRANSFO,
+ FUN_MOD = FUN_MOD, FUN_TRANSFO = FUN_TRANSFO,
verbose = verbose, ...))
}
diff --git a/R/Calibration_Michel.R b/R/Calibration_Michel.R
index 8a59765e57dfd2f28bbd0fe2d3b0af9e6d4e554d..c7119cba1a00fbcef54e38d7a499b46418a28f75 100644
--- a/R/Calibration_Michel.R
+++ b/R/Calibration_Michel.R
@@ -3,7 +3,6 @@ Calibration_Michel <- function(InputsModel,
InputsCrit,
CalibOptions,
FUN_MOD,
- FUN_SD,
FUN_CRIT, # deprecated
FUN_TRANSFO = NULL,
verbose = TRUE,
@@ -148,7 +147,7 @@ Calibration_Michel <- function(InputsModel,
}
##Model_run
Param <- CandidatesParamR[iNew, ]
- OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, FUN_SD = FUN_SD,...)
+ OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, FUN_SD = RunOptions$FUN_SD,...)
##Calibration_criterion_computation
OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
diff --git a/R/CreateCalibOptions.R b/R/CreateCalibOptions.R
index 8a2e9a3bbeae5bc0262a719c8db547f01c0df6de..a7b0810feaecf4b4737a0d661b0d8fe9cddbb6cb 100644
--- a/R/CreateCalibOptions.R
+++ b/R/CreateCalibOptions.R
@@ -64,7 +64,7 @@ CreateCalibOptions <- function(FUN_MOD,
NParam <- NParam + 2
}
if (IsSD) {
- NParam <- NParam + FeatFUN_MOD$ParamSD
+ NParam <- NParam + FeatFUN_MOD$NParamSD
}
## check FixedParam
@@ -204,6 +204,12 @@ 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)
+ ParamT <- cbind(ParamTSD, ParamT)
+ }
}
StartParamList <- NULL
StartParamDistrib <- TransfoParam(ParamIn = ParamT, Direction = "TR", FUN_TRANSFO = FUN_TRANSFO)
diff --git a/R/CreateRunOptions.R b/R/CreateRunOptions.R
index 76daea77263c05c010abdbc74805b47120bf260a..9a8b4e9082ce32dbccdb775761f141b2ba05a199 100644
--- a/R/CreateRunOptions.R
+++ b/R/CreateRunOptions.R
@@ -44,7 +44,7 @@ CreateRunOptions <- function(FUN_MOD, InputsModel,
## SD model
FeatFUN_MOD$IsSD <- inherits(InputsModel, "SD")
if (FeatFUN_MOD$IsSD) {
- FeatFUN_MOD$NbParam <- FeatFUN_MOD$NbParam + FeatFUN_MOD$ParamSD
+ FeatFUN_MOD$NbParam <- FeatFUN_MOD$NbParam + FeatFUN_MOD$NParamSD
}
if (!"CemaNeige" %in% ObjectClass & "hysteresis" %in% ObjectClass) {
@@ -466,7 +466,8 @@ CreateRunOptions <- function(FUN_MOD, InputsModel,
Outputs_Cal = Outputs_Cal,
Outputs_Sim = Outputs_Sim,
FortranOutputs = FortranOutputs,
- FeatFUN_MOD = FeatFUN_MOD)
+ FeatFUN_MOD = FeatFUN_MOD,
+ FUN_SD = FUN_SD)
if ("CemaNeige" %in% ObjectClass) {
RunOptions <- c(RunOptions, list(MeanAnSolidPrecip = MeanAnSolidPrecip))
diff --git a/R/RunModel.R b/R/RunModel.R
index 378f87e3f7963dc3f9809c5b161f9dfb4f882ca1..685cfbe9f2c96f731e7372a2f5101f691da4a1b9 100644
--- a/R/RunModel.R
+++ b/R/RunModel.R
@@ -4,8 +4,8 @@ RunModel <- function(InputsModel, RunOptions, Param, FUN_MOD, FUN_SD = RunModel_
if (inherits(InputsModel, "SD") && !identical(FUN_MOD, FUN_SD)) {
# Lag model take one parameter at the beginning of the vector
- iFirstParamRunOffModel <- RunOptions$FeatFUN_MOD$ParamSD + 1
- RunOptions$FeatFUN_MOD$NbParam <- RunOptions$FeatFUN_MOD$NbParam - RunOptions$FeatFUN_MOD$ParamSD
+ iFirstParamRunOffModel <- RunOptions$FeatFUN_MOD$NParamSD + 1
+ RunOptions$FeatFUN_MOD$NbParam <- RunOptions$FeatFUN_MOD$NbParam - RunOptions$FeatFUN_MOD$NParamSD
} else {
# All parameters
iFirstParamRunOffModel <- 1
@@ -14,7 +14,7 @@ RunModel <- function(InputsModel, RunOptions, Param, FUN_MOD, FUN_SD = RunModel_
Param = Param[iFirstParamRunOffModel:length(Param)], ...)
if (inherits(InputsModel, "SD") && !identical(FUN_MOD, FUN_SD)) {
- OutputsModel <- RunModel_Lag(InputsModel, RunOptions, Param[1:RunOptions$FeatFUN_MOD$ParamSD],
+ OutputsModel <- RunModel_Lag(InputsModel, RunOptions, Param[1:RunOptions$FeatFUN_MOD$NParamSD],
OutputsModel)
}
return(OutputsModel)
diff --git a/R/RunModel_GR5H.R b/R/RunModel_GR5H.R
index 6c5fb4323b477b1316d51feede865e405c23206f..29d0d39a91c414fb5398a3797e0063f5789f5b02 100644
--- a/R/RunModel_GR5H.R
+++ b/R/RunModel_GR5H.R
@@ -1,93 +1,93 @@
-RunModel_GR5H <- function(InputsModel, RunOptions, Param) {
-
-
- ## Initialization of variables
- IsIntStore <- inherits(RunOptions, "interception")
- if (IsIntStore) {
- Imax <- RunOptions$Imax
- } else {
- Imax <- -99
- }
-
- .ArgumentsCheckGR(InputsModel, RunOptions, Param)
-
- Param <- as.double(Param)
-
- Param_X1X3_threshold <- 1e-2
- Param_X4_threshold <- 0.5
- if (Param[1L] < Param_X1X3_threshold) {
- warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
- Param[1L] <- Param_X1X3_threshold
- }
- if (Param[3L] < Param_X1X3_threshold) {
- warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
- Param[3L] <- Param_X1X3_threshold
- }
- if (Param[4L] < Param_X4_threshold) {
- warning(sprintf("Param[4] (X4: unit hydrograph time constant [h]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
- Param[4L] <- Param_X4_threshold
- }
-
- ## Input data preparation
- if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
- RunOptions$IndPeriod_WarmUp <- NULL
- }
- IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
- LInputSeries <- as.integer(length(IndPeriod1))
- if ("all" %in% RunOptions$Outputs_Sim) {
- IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
- } else {
- IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
- }
-
- ## Output data preparation
- IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
- ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
- ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
-
- ## Use of IniResLevels
- if (!is.null(RunOptions$IniResLevels)) {
- RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
- RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
- if (IsIntStore) {
- RunOptions$IniStates[4] <- RunOptions$IniResLevels[4] * Imax ### interception store level (mm)
- }
- }
-
- ## Call GR model Fortan
- RESULTS <- .Fortran("frun_gr5h", PACKAGE = "airGR",
- ## inputs
- LInputs = LInputSeries, ### length of input and output series
- InputsPrecip = InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/h]
- InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/h]
- NParam = as.integer(length(Param)), ### number of model parameter
- Param = Param, ### parameter set
- NStates = as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
- StateStart = RunOptions$IniStates, ### state variables used when the model run starts
- Imax = Imax, ### maximal capacity of interception store
- NOutputs = as.integer(length(IndOutputs)), ### number of output series
- IndOutputs = IndOutputs, ### indices of output series
- ## outputs
- Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm or mm/h]
- StateEnd = rep(as.double(-99e9), length(RunOptions$IniStates)) ### state variables at the end of the model run
- )
- RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
- RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
- if (ExportStateEnd) {
- RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
- RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR5H, InputsModel = InputsModel,
- ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
- IntStore = RESULTS$StateEnd[4L],
- UH1 = NULL,
- UH2 = RESULTS$StateEnd[(1:(40*24)) + (7+20*24)],
- GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
- verbose = FALSE)
- }
-
- ## OutputsModel generation
- .GetOutputsModelGR(InputsModel,
- RunOptions,
- RESULTS,
- LInputSeries,
- Param)
-}
+RunModel_GR5H <- function(InputsModel, RunOptions, Param) {
+
+
+ ## Initialization of variables
+ IsIntStore <- inherits(RunOptions, "interception")
+ if (IsIntStore) {
+ Imax <- RunOptions$Imax
+ } else {
+ Imax <- -99
+ }
+
+ .ArgumentsCheckGR(InputsModel, RunOptions, Param)
+
+ Param <- as.double(Param)
+
+ Param_X1X3_threshold <- 1e-2
+ Param_X4_threshold <- 0.5
+ if (Param[1L] < Param_X1X3_threshold) {
+ warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
+ Param[1L] <- Param_X1X3_threshold
+ }
+ if (Param[3L] < Param_X1X3_threshold) {
+ warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
+ Param[3L] <- Param_X1X3_threshold
+ }
+ if (Param[4L] < Param_X4_threshold) {
+ warning(sprintf("Param[4] (X4: unit hydrograph time constant [h]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
+ Param[4L] <- Param_X4_threshold
+ }
+
+ ## Input data preparation
+ if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
+ RunOptions$IndPeriod_WarmUp <- NULL
+ }
+ IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
+ LInputSeries <- as.integer(length(IndPeriod1))
+ if ("all" %in% RunOptions$Outputs_Sim) {
+ IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
+ } else {
+ IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
+ }
+
+ ## Output data preparation
+ IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
+ ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
+ ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
+
+ ## Use of IniResLevels
+ if (!is.null(RunOptions$IniResLevels)) {
+ RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
+ RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
+ if (IsIntStore) {
+ RunOptions$IniStates[4] <- RunOptions$IniResLevels[4] * Imax ### interception store level (mm)
+ }
+ }
+
+ ## Call GR model Fortan
+ RESULTS <- .Fortran("frun_gr5h", PACKAGE = "airGR",
+ ## inputs
+ LInputs = LInputSeries, ### length of input and output series
+ InputsPrecip = InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/h]
+ InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/h]
+ NParam = as.integer(length(Param)), ### number of model parameter
+ Param = Param, ### parameter set
+ NStates = as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
+ StateStart = RunOptions$IniStates, ### state variables used when the model run starts
+ Imax = Imax, ### maximal capacity of interception store
+ NOutputs = as.integer(length(IndOutputs)), ### number of output series
+ IndOutputs = IndOutputs, ### indices of output series
+ ## outputs
+ Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm or mm/h]
+ StateEnd = rep(as.double(-99e9), length(RunOptions$IniStates)) ### state variables at the end of the model run
+ )
+ RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
+ RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
+ if (ExportStateEnd) {
+ RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
+ RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR5H, InputsModel = InputsModel,
+ ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
+ IntStore = RESULTS$StateEnd[4L],
+ UH1 = NULL,
+ UH2 = RESULTS$StateEnd[(1:(40*24)) + (7+20*24)],
+ GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
+ verbose = FALSE)
+ }
+
+ ## OutputsModel generation
+ .GetOutputsModelGR(InputsModel,
+ RunOptions,
+ RESULTS,
+ LInputSeries,
+ Param)
+}
diff --git a/R/RunModel_LLR.R b/R/RunModel_LLR.R
new file mode 100644
index 0000000000000000000000000000000000000000..f5dff3161c35775e4d7aff98ca6be4b561fcea42
--- /dev/null
+++ b/R/RunModel_LLR.R
@@ -0,0 +1,183 @@
+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]
+
+ C0 <- exp(-1/KParam)
+ C1 <- ((KParam/1) * (1-C0))-C0
+ C2 <- 1 - (KParam/1) * (1-C0)
+
+
+ ## propagation time from upstream meshes to outlet
+ PT <- InputsModel$LengthHydro * 1e3 / TParam / RunOptions$FeatFUN_MOD$TimeStep
+ print(PT)
+ HUTRANS <- rbind(1 - (PT - floor(PT)), PT - floor(PT))
+ ## 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))
+ }
+ )
+ }
+ # message("IniStates: ", paste(IniStates, collapse = ", "))
+
+ ## Lag model computation
+ Qsim_m3 <- QsimDown *
+ InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
+
+ for (upstream_basin in seq_len(NbUpBasins)) {
+ Qupstream <- c(IniStates[[upstream_basin]],
+ InputsModel$Qupstream[IndPeriod1, upstream_basin])
+ # message("Qupstream[", upstream_basin, "]: ", paste(Qupstream, collapse = ", "))
+ Qsim_m3 <- Qsim_m3 +
+ Qupstream[2:(1 + LInputSeries)] * HUTRANS[1, upstream_basin] +
+ Qupstream[1:LInputSeries] * HUTRANS[2, upstream_basin]
+ }
+
+ ## 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 {
+ 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.R b/R/TransfoParam.R
index dcda7b088a732c5a867cb13e09d0f7d2784861ed..6173285ce7f383f71cda8cfaca7182e134005c01 100644
--- a/R/TransfoParam.R
+++ b/R/TransfoParam.R
@@ -1,4 +1,4 @@
-TransfoParam <- function(ParamIn, Direction, FUN_TRANSFO) {
- FUN_TRANSFO <- match.fun(FUN_TRANSFO)
- return(FUN_TRANSFO(ParamIn, Direction))
-}
+TransfoParam <- function(ParamIn, Direction, FUN_TRANSFO) {
+ FUN_TRANSFO <- match.fun(FUN_TRANSFO)
+ return(FUN_TRANSFO(ParamIn, Direction))
+}
diff --git a/R/TransfoParam_LLR.R b/R/TransfoParam_LLR.R
new file mode 100644
index 0000000000000000000000000000000000000000..d791f8e4f11be93350ffa4a1f782f1a34bdc95da
--- /dev/null
+++ b/R/TransfoParam_LLR.R
@@ -0,0 +1,39 @@
+TransfoParam_LLR <- 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[, 1] <- 20 * (ParamIn[, 1] + 10) / 20.0
+ ParamOut[, 2] <- 20 * (ParamIn[, 2] + 10) / 20.0
+ }
+ if (Direction == "RT") {
+ ParamOut[, 1] <- ParamIn[, 1] * 20.0 / 20 - 10
+ ParamOut[, 2] <- 20 * (ParamIn[, 2] + 10) / 20.0
+ }
+
+
+ ## check output
+ if (isVecParamIn) {
+ ParamOut <- as.vector(ParamOut)
+ }
+
+ return(ParamOut)
+
+}
diff --git a/R/Utils.R b/R/Utils.R
index 2e84abd2ee019e06ec830ce836f31cc7cdf7ed2b..f3716df2bbdf5945ed4a8e22ab45dde021ac5c0b 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -74,8 +74,11 @@
}
}
if (!is.null(FUN_SD)){
- if (identical(match.fun(RunModel_Lag), FUN_SD)){
- res$ParamSD <- FeatMod$NbParam[14]
+ if (identical(RunModel_Lag, FUN_SD)){
+ res$NParamSD <- FeatMod$NbParam[14]
+ }
+ if(identical(RunModel_LLR, FUN_SD)){
+ res$NParamSD <- FeatMod$NbParam[15]
}
}
return(res)
diff --git a/R/UtilsCalibOptions.R b/R/UtilsCalibOptions.R
index 92c9fa2f323d03293d3f35d249e83f23793e88d1..a9ecfe084e28ad51b7738a62c5aa1e6e7dbe119e 100644
--- a/R/UtilsCalibOptions.R
+++ b/R/UtilsCalibOptions.R
@@ -30,6 +30,9 @@
if (identical(RunModel_Lag, FUN_SD)){
FUN_ROUT <- TransfoParam_Lag
}
+ if (identical(RunModel_LLR, FUN_SD)){
+ FUN_ROUT <- TransfoParam_LLR
+ }
}
@@ -45,8 +48,8 @@
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
- ParamOut[, (FeatFUN_MOD$ParamSD + 1):NParam] <- FUN_GR(ParamIn[, (FeatFUN_MOD$ParamSD + 1):NParam], Direction)
- ParamOut[, FeatFUN_MOD$ParamSD ] <- FUN_ROUT(as.matrix(ParamIn[, FeatFUN_MOD$ParamSD]), Direction)
+ ParamOut[, (FeatFUN_MOD$NParamSD + 1):NParam] <- FUN_GR(ParamIn[, (FeatFUN_MOD$NParamSD + 1):NParam], Direction)
+ ParamOut[, FeatFUN_MOD$NParamSD ] <- FUN_ROUT(as.matrix(ParamIn[, FeatFUN_MOD$NParamSD]), Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
diff --git a/inst/modelsFeatures/FeatModelsGR.csv b/inst/modelsFeatures/FeatModelsGR.csv
index 38d9e5ca7bfab6eef3e8932e1f3c9ebadd5a220a..db34200d5bf8ba1c2a177215323162c319320643 100644
--- a/inst/modelsFeatures/FeatModelsGR.csv
+++ b/inst/modelsFeatures/FeatModelsGR.csv
@@ -1,15 +1,16 @@
-CodeMod;NameMod;NbParam;TimeUnit;Id;Class;Pkg
-GR1A;GR1A;1;yearly;NA;GR;airGR
-GR2M;GR2M;2;monthly;NA;GR;airGR
-GR4J;GR4J;4;daily;NA;GR;airGR
-GR5J;GR5J;5;daily;NA;GR;airGR
-GR6J;GR6J;6;daily;NA;GR;airGR
-GR4H;GR4H;4;hourly;NA;GR;airGR
-GR5H;GR5H;5;hourly;NA;GR;airGR
-CemaNeige;CemaNeige;2;NA;NA;NA;airGR
-CemaNeigeGR4J;CemaNeigeGR4J;6;daily;NA;GR;airGR
-CemaNeigeGR5J;CemaNeigeGR5J;7;daily;NA;GR;airGR
-CemaNeigeGR6J;CemaNeigeGR6J;8;daily;NA;GR;airGR
-CemaNeigeGR4H;CemaNeigeGR4H;6;hourly;NA;GR;airGR
-CemaNeigeGR5H;CemaNeigeGR5H;7;hourly;NA;GR;airGR
-Lag;Lag;1;NA;NA;SD;airGR
+CodeMod;NameMod;NbParam;TimeUnit;Id;Class;Pkg
+GR1A;GR1A;1;yearly;NA;GR;airGR
+GR2M;GR2M;2;monthly;NA;GR;airGR
+GR4J;GR4J;4;daily;NA;GR;airGR
+GR5J;GR5J;5;daily;NA;GR;airGR
+GR6J;GR6J;6;daily;NA;GR;airGR
+GR4H;GR4H;4;hourly;NA;GR;airGR
+GR5H;GR5H;5;hourly;NA;GR;airGR
+CemaNeige;CemaNeige;2;NA;NA;NA;airGR
+CemaNeigeGR4J;CemaNeigeGR4J;6;daily;NA;GR;airGR
+CemaNeigeGR5J;CemaNeigeGR5J;7;daily;NA;GR;airGR
+CemaNeigeGR6J;CemaNeigeGR6J;8;daily;NA;GR;airGR
+CemaNeigeGR4H;CemaNeigeGR4H;6;hourly;NA;GR;airGR
+CemaNeigeGR5H;CemaNeigeGR5H;7;hourly;NA;GR;airGR
+Lag;Lag;1;NA;NA;SD;airGR
+LinearLagAndRoute;Lag;2;NA;NA;SD;airGR