diff --git a/R/CreateInputsModel.R b/R/CreateInputsModel.R
index 3bd0f14752105a78f3010a3d8cd744d646f6fa2e..7a24cf3c98fe3bb6e9bc4d1dfa20aca275f5f664 100644
--- a/R/CreateInputsModel.R
+++ b/R/CreateInputsModel.R
@@ -5,6 +5,7 @@ CreateInputsModel <- function(FUN_MOD,
TempMean = NULL, TempMin = NULL, TempMax = NULL,
ZInputs = NULL, HypsoData = NULL, NLayers = 5,
Qupstream = NULL, LengthHydro = NULL, BasinAreas = NULL,
+ QupstrUnit = "mm",
verbose = TRUE) {
@@ -215,6 +216,8 @@ CreateInputsModel <- function(FUN_MOD,
if(any(LengthHydro > 1000)) {
warning("The unit of 'LengthHydro' has changed from m to km in v1.7 of airGR: values superior to 1000 km seem unrealistic")
}
+ QupstrUnit <- tolower(QupstrUnit)
+ QupstrUnit <- match.arg(arg = QupstrUnit, choices = c("mm", "m3", "m3/s", "l/s"))
}
##check_NA_values
@@ -330,6 +333,16 @@ CreateInputsModel <- function(FUN_MOD,
ZLayers = RESULT$ZLayers))
}
if ("SD" %in% ObjectClass) {
+ # Qupstream is internally stored in m3/time step
+ if (QupstrUnit == "mm") {
+ iConvBasins <- which(!is.na(BasinAreas[seq.int(length(LengthHydro))]))
+ Qupstream[,iConvBasins] <-
+ Qupstream[,iConvBasins] * rep(BasinAreas[iConvBasins], each = LLL) * 1e3
+ } else if (QupstrUnit == "m3/s") {
+ Qupstream <- Qupstream * TimeStep
+ } else if (QupstrUnit == "l/s") {
+ Qupstream <- Qupstream * TimeStep / 1e3
+ }
InputsModel <- c(InputsModel, list(Qupstream = Qupstream,
LengthHydro = LengthHydro,
BasinAreas = BasinAreas))
diff --git a/R/RunModel_Lag.R b/R/RunModel_Lag.R
index b7f4398897b991ff7cb54980f12aad0f930ef98f..473e7758993b58d2c2317e29279b76424e2b25f8 100644
--- a/R/RunModel_Lag.R
+++ b/R/RunModel_Lag.R
@@ -48,7 +48,7 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
NbUpBasins <- length(InputsModel$LengthHydro)
LengthTs <- length(OutputsModel$QsimDown)
- OutputsModel$Qsim <- OutputsModel$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
+ OutputsModel$Qsim_m3 <- OutputsModel$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
IniSD <- RunOptions$IniStates[grep("SD", names(RunOptions$IniStates))]
if (length(IniSD) > 0) {
@@ -78,15 +78,15 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
for (upstream_basin in seq_len(NbUpBasins)) {
Qupstream <- c(IniStates[[upstream_basin]],
InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin])
- if (!is.na(InputsModel$BasinAreas[upstream_basin])) {
- # Upstream flow with area needs to be converted to m3 by time step
- Qupstream <- Qupstream * InputsModel$BasinAreas[upstream_basin] * 1e3
- }
# message("Qupstream[", upstream_basin, "]: ", paste(Qupstream, collapse = ", "))
- OutputsModel$Qsim <- OutputsModel$Qsim +
+ OutputsModel$Qsim_m3 <- OutputsModel$Qsim_m3 +
Qupstream[2:(1 + LengthTs)] * HUTRANS[1, upstream_basin] +
Qupstream[1:LengthTs] * HUTRANS[2, upstream_basin]
}
+ # Convert back Qsim to mm
+ OutputsModel$Qsim <- OutputsModel$Qsim_m3 / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
+ # message("Qsim: ", paste(OutputsModel$Qsim, collapse = ", "))
+
# 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)) {
@@ -98,9 +98,6 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
warning(length(which(is.na(OutputsModel$Qsim))), " time steps with NA values")
}
}
- # Convert back Qsim to mm
- OutputsModel$Qsim <- OutputsModel$Qsim / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
- # message("Qsim: ", paste(OutputsModel$Qsim, collapse = ", "))
if ("StateEnd" %in% RunOptions$Outputs_Sim) {
OutputsModel$StateEnd$SD <- lapply(seq(NbUpBasins), function(x) {
diff --git a/man/Calibration_Michel.Rd b/man/Calibration_Michel.Rd
index 60439df643b4e82a52b8e7a59051b3fdfaa66d9b..28de8c4424cc84889a83fdbd9a123d5329c7af3f 100644
--- a/man/Calibration_Michel.Rd
+++ b/man/Calibration_Michel.Rd
@@ -19,7 +19,7 @@ Then a steepest descent local search algorithm is performed, starting from the r
\usage{
Calibration_Michel(InputsModel, RunOptions, InputsCrit, CalibOptions,
- FUN_MOD, FUN_CRIT, FUN_TRANSFO = NULL, verbose = TRUE)
+ FUN_MOD, FUN_CRIT, FUN_TRANSFO = NULL, verbose = TRUE, ...)
}
@@ -39,6 +39,8 @@ Calibration_Michel(InputsModel, RunOptions, InputsCrit, CalibOptions,
\item{FUN_TRANSFO}{(optional) [function] model parameters transformation function, if the \code{FUN_MOD} used is native in the package \code{FUN_TRANSFO} is automatically defined}
\item{verbose}{(optional) [boolean] boolean indicating if the function is run in verbose mode or not, default = \code{TRUE}}
+
+\item{...}{(optional) arguments to pass to \code{\link{RunModel}}}
}
diff --git a/man/CreateInputsModel.Rd b/man/CreateInputsModel.Rd
index 365955eb5d5af4700695578bcaae7e61ebfba4ae..16d8c1b400044f792baaabb51f61d9b83592a6af 100644
--- a/man/CreateInputsModel.Rd
+++ b/man/CreateInputsModel.Rd
@@ -19,7 +19,7 @@ CreateInputsModel(FUN_MOD, DatesR, Precip, PrecipScale = TRUE, PotEvap = NULL,
TempMean = NULL, TempMin = NULL, TempMax = NULL,
ZInputs = NULL, HypsoData = NULL, NLayers = 5,
Qupstream = NULL, LengthHydro = NULL, BasinAreas = NULL,
- verbose = TRUE)
+ QupstrUnit = "mm", verbose = TRUE)
\method{[}{InputsModel}(x, i)
}
@@ -50,12 +50,14 @@ CreateInputsModel(FUN_MOD, DatesR, Precip, PrecipScale = TRUE, PotEvap = NULL,
\item{verbose}{(optional) [boolean] boolean indicating if the function is run in verbose mode or not, default = \code{TRUE}}
-\item{Qupstream}{(optional) [numerical matrix] time series of upstream flows (catchment average) [mm/time step or m3/time step, see details], required to create the SD model inputs . See details}
+\item{Qupstream}{(optional) [numerical matrix] time series of upstream flows (catchment average), its unit is defined by the \code{QupstrUnit} parameter, required to create the SD model inputs. See details}
\item{LengthHydro}{(optional) [numeric] real giving the distance between the downstream outlet and each upstream inlet of the sub-catchment [km], required to create the SD model inputs . See details}
\item{BasinAreas}{(optional) [numeric] real giving the area of each upstream sub-catchment [km2] and the area of the downstream sub-catchment in the last item, required to create the SD model inputs . See details}
+\item{QupstrUnit}{(optional) [character] unit of the flow in the argument \code{Qupstream}, available units are: "mm" for mm/time-step (default), "m3" for m3/time-step, "m3/s" and "l/s". See details}
+
\item{x}{[InputsModel] object of class InputsModel}
\item{i}{[integer] of the indices to subset a time series or [character] names of the elements to extract}
@@ -80,9 +82,11 @@ Users wanting to use \code{FUN_MOD} functions that are not included in
the package must create their own InputsModel object accordingly. \cr
Please note that if CemaNeige is used, and \code{ZInputs} is different than \code{HypsoData}, then precipitation and temperature are interpolated with the \code{DataAltiExtrapolation_Valery} function.
-Users wanting to use a semi-distributed (SD) lag model should provide valid \code{Qupstream}, \code{LengthHydro}, and \code{BasinAreas} parameters. Each upstream sub-catchment is described by an upstream flow time series (one column in \code{Qupstream} matrix), a distance between the downstream outlet and the upstream inlet (one item in \code{LengthHydro}) and an area (one item in \code{BasinAreas}).
-The order of the columns or of the items should be consistent for all these parameters. \code{BasinAreas} should contain one extra information (stored in the last item) which is the area of the downstream sub-catchment.
-Upstream flows that are not related to a sub-catchment such as a release or withdraw spot are identified by an area equal to \code{NA} and an upstream flow expressed in m3/time step instead of mm/time step.
+Users wanting to use a semi-distributed (SD) model should provide valid \code{Qupstream}, \code{LengthHydro}, and \code{BasinAreas} arguments. Each upstream sub-catchment is described by an upstream flow time series (one column in \code{Qupstream} matrix), a distance between the downstream outlet and the upstream inlet (one item in \code{LengthHydro}) and an area (one item in \code{BasinAreas}).
+The order of the columns or of the items should be consistent for all these parameters.
+\code{BasinAreas} should contain one extra information (stored in the last item) which is the area of the downstream sub-catchment.
+Upstream flows that are not related to a sub-catchment such as release or withdraw spots are identified by an area equal to \code{NA}, and if \code{unit="mm"} the upstream flow must be expressed in m3/time step instead of mm/time step which is not possible in absence of a related area.
+Please note that the use of SD model requires to use the \code{\link{RunModel}} function instead of \code{\link{RunModel_GR4J}} or the other \code{RunModel_*} functions.
}
\examples{
diff --git a/man/RunModel.Rd b/man/RunModel.Rd
index 2e76d41d065a140ff0440005ab6635eff63de620..6762ed73594e76e2af745523f5f6a039bd76aa91 100644
--- a/man/RunModel.Rd
+++ b/man/RunModel.Rd
@@ -15,7 +15,7 @@ Function which performs a single model run with the provided function over the s
\usage{
-RunModel(InputsModel, RunOptions, Param, FUN_MOD)
+RunModel(InputsModel, RunOptions, Param, FUN_MOD, ...)
%
%\method{[}{OutputsModel}(x, i)
}
@@ -29,6 +29,8 @@ RunModel(InputsModel, RunOptions, Param, FUN_MOD)
\item{Param}{[numeric] vector of model parameters (See details for SD lag model)}
\item{FUN_MOD}{[function] hydrological model function (e.g. \code{\link{RunModel_GR4J}}, \code{\link{RunModel_CemaNeigeGR4J}})}
+
+\item{...}{(optional) arguments to pass to \code{FUN_MOD}}
%
%\item{x}{[InputsModel] object of class InputsModel}
%
diff --git a/tests/testthat/test-RunModel_Lag.R b/tests/testthat/test-RunModel_Lag.R
index 8b83924c280147c3a826e940a66afaffa20de40b..26e10832da4e5f5961733b469b6bef8e94889951 100644
--- a/tests/testthat/test-RunModel_Lag.R
+++ b/tests/testthat/test-RunModel_Lag.R
@@ -39,6 +39,22 @@ RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel,
IndPeriod_Run = Ind_Run))
+test_that("'QupstrUnit' must correspond to one possible value", {
+ expect_error(
+ InputsModel <- CreateInputsModel(
+ FUN_MOD = RunModel_GR4J,
+ DatesR = BasinObs$DatesR,
+ Precip = BasinObs$P,
+ PotEvap = BasinObs$E,
+ Qupstream = matrix(Qupstream, ncol = 1),
+ LengthHydro = 1,
+ BasinAreas = BasinAreas,
+ QupstrUnit = "m3/h"
+ ),
+ regexp = "'arg' should be one of \"mm\", \"m3\", \"m3/s\", \"l/s\""
+ )
+})
+
test_that("QcontribDown parameter should be a numeric vector or an OutputModel object", {
regexp = "'QcontribDown' must be a numeric vector or a 'OutputsModel' object"
expect_error(
@@ -125,8 +141,15 @@ test_that("'Qupstream' contain NA values", {
})
test_that("Upstream basin with nil area should return same Qdown as GR4J alone", {
- UpstBasinArea <- InputsModel$BasinAreas[1L]
- InputsModel$BasinAreas[1L] <- 0
+ InputsModel <- CreateInputsModel(
+ FUN_MOD = RunModel_GR4J,
+ DatesR = BasinObs$DatesR,
+ Precip = BasinObs$P,
+ PotEvap = BasinObs$E,
+ Qupstream = matrix(Qupstream, ncol = 1),
+ LengthHydro = 1,
+ BasinAreas = c(0,BasinAreas[2])
+ )
OutputsSD <- RunModel(InputsModel,
RunOptions,
Param = c(1, Param),
@@ -135,8 +158,15 @@ test_that("Upstream basin with nil area should return same Qdown as GR4J alone",
})
test_that("Downstream basin with nil area and nul upstream length should return same Qdown as Qupstream alone", {
- InputsModel$LengthHydro <- 0
- InputsModel$BasinAreas <- c(BasinInfo$BasinArea, 0)
+ InputsModel <- CreateInputsModel(
+ FUN_MOD = RunModel_GR4J,
+ DatesR = BasinObs$DatesR,
+ Precip = BasinObs$P,
+ PotEvap = BasinObs$E,
+ Qupstream = matrix(Qupstream, ncol = 1),
+ LengthHydro = 0,
+ BasinAreas = c(BasinInfo$BasinArea, 0)
+ )
OutputsSD <- RunModel(InputsModel,
RunOptions,
Param = c(1, Param),
@@ -146,22 +176,72 @@ test_that("Downstream basin with nil area and nul upstream length should return
ParamSD <- c(InputsModel$LengthHydro * 1e3 / (24 * 60 * 60), Param) # Speed corresponding to one time step delay
-QlsGR4Only <- OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2L] * 1e6 / 86400
+Qm3GR4Only <- OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2L] * 1e3
test_that("1 input with lag of 1 time step delay out gives an output delayed of one time step", {
OutputsSD <- RunModel(InputsModel, RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
- QlsSdSim <- OutputsSD$Qsim * sum(InputsModel$BasinAreas) * 1e6 / 86400
- QlsUpstLagObs <- c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]) * InputsModel$BasinAreas[1L] * 1e6 / 86400
- expect_equal(QlsSdSim - QlsGR4Only, QlsUpstLagObs)
+ Qm3SdSim <- OutputsSD$Qsim_m3
+ Qm3UpstLagObs <- c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]) * InputsModel$BasinAreas[1] * 1e3
+ expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
})
test_that("1 input with lag of 0.5 time step delay out gives an output delayed of 0.5 time step", {
OutputsSD <- RunModel(InputsModel, RunOptions,
Param = c(InputsModel$LengthHydro * 1e3 / (12 * 3600), Param),
FUN_MOD = RunModel_GR4J)
- QlsSdSim <- OutputsSD$Qsim * sum(InputsModel$BasinAreas) * 1e6 / 86400
- QlsUpstLagObs <- (Qupstream[Ind_Run] + c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]))/2 * InputsModel$BasinAreas[1L] * 1e6 / 86400
- expect_equal(QlsSdSim - QlsGR4Only, QlsUpstLagObs)
+ Qm3SdSim <- OutputsSD$Qsim_m3
+ Qm3UpstLagObs <- (Qupstream[Ind_Run] + c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]])) / 2 * InputsModel$BasinAreas[1] * 1e3
+ expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
+})
+
+test_that("Qupstream in different units should return the same result", {
+ OutputsSD_mm <- RunModel(InputsModel, RunOptions,
+ Param = ParamSD,
+ FUN_MOD = RunModel_GR4J)
+ InputsModel_m3 <- CreateInputsModel(
+ FUN_MOD = RunModel_GR4J,
+ DatesR = BasinObs$DatesR,
+ Precip = BasinObs$P,
+ PotEvap = BasinObs$E,
+ Qupstream = matrix(Qupstream, ncol = 1) * BasinAreas[1] * 1e3,
+ LengthHydro = 1,
+ BasinAreas = BasinAreas,
+ QupstrUnit = "m3"
+ )
+ OutputsSD_m3 <- RunModel(InputsModel_m3, RunOptions,
+ Param = ParamSD,
+ FUN_MOD = RunModel_GR4J)
+ expect_equal(OutputsSD_mm$Qsim, OutputsSD_m3$Qsim)
+
+ InputsModel_m3s <- CreateInputsModel(
+ FUN_MOD = RunModel_GR4J,
+ DatesR = BasinObs$DatesR,
+ Precip = BasinObs$P,
+ PotEvap = BasinObs$E,
+ Qupstream = matrix(Qupstream, ncol = 1) * BasinAreas[1] * 1e3 / 86400,
+ LengthHydro = 1,
+ BasinAreas = BasinAreas,
+ QupstrUnit = "m3/s"
+ )
+ OutputsSD_m3s <- RunModel(InputsModel_m3s, RunOptions,
+ Param = ParamSD,
+ FUN_MOD = RunModel_GR4J)
+ expect_equal(OutputsSD_mm$Qsim, OutputsSD_m3s$Qsim)
+
+ InputsModel_ls <- CreateInputsModel(
+ FUN_MOD = RunModel_GR4J,
+ DatesR = BasinObs$DatesR,
+ Precip = BasinObs$P,
+ PotEvap = BasinObs$E,
+ Qupstream = matrix(Qupstream, ncol = 1) * BasinAreas[1] * 1e6 / 86400,
+ LengthHydro = 1,
+ BasinAreas = BasinAreas,
+ QupstrUnit = "L/s"
+ )
+ OutputsSD_ls <- RunModel(InputsModel_ls, RunOptions,
+ Param = ParamSD,
+ FUN_MOD = RunModel_GR4J)
+ expect_equal(OutputsSD_mm$Qsim, OutputsSD_ls$Qsim)
})
InputsCrit <- CreateInputsCrit(
@@ -218,7 +298,7 @@ test_that("1 no area input with lag of 1 time step delay out gives an output del
expect_false(any(is.na(OutputsSD$Qsim)))
- Qm3SdSim <- OutputsSD$Qsim * sum(InputsModel$BasinAreas, na.rm = TRUE) * 1e3
+ Qm3SdSim <- OutputsSD$Qsim_m3
Qm3UpstLagObs <- c(0, InputsModel$Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]])
expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)