v1.6.2.1 feat(SD): direct flow in m3 per time step for Qupstream with area at `NA`

Refs #34

DESCRIPTION

 Package: airGR
 Type: Package
 Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.6.1.16
-Date: 2020-06-04
+Version: 1.6.2.1
+Date: 2020-06-05
 Authors@R: c(
   person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
   person("Olivier", "Delaigue", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@inrae.fr"),

NEWS.md

 ## Release History of the airGR Package
 
+### 1.6.2.1 Release Notes (2020-06-05)
+
+#### New features
+
+- Add direct upstream flow not associated to an area in semi-distributed model
+
+____________________________________________________________________________________
 
 
 ### 1.6.1.11 Release Notes (2020-04-07)

R/RunModel.R

@@ -14,27 +14,30 @@ RunModel <- function (InputsModel, RunOptions, Param, FUN_MOD) {
       TimeStep <- 60 * 60
     }
 
-    # total area
-    AreaTot <- sum(InputsModel$BasinAreas)
-
     # propagation time from upstream meshes to outlet
     PT <- InputsModel$LengthHydro / Param[length(Param)] / TimeStep
     HUTRANS <- rbind(1 - (PT - floor(PT)), PT - floor(PT))
 
     NbUpBasins <- length(InputsModel$LengthHydro)
     LengthTs <- length(OutputsModelDown$QsimDown)
-    OutputsModelDown$Qsim <- OutputsModelDown$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] / AreaTot
+    OutputsModelDown$Qsim <- OutputsModelDown$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1E3
 
     for (upstream_basin in seq_len(NbUpBasins)) {
       Qupstream <- 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
+      }
       OutputsModelDown$Qsim <- OutputsModelDown$Qsim +
         c(rep(0, floor(PT[upstream_basin])),
           Qupstream[(1 + floor(PT[upstream_basin])):LengthTs]) *
-        HUTRANS[1, upstream_basin] * InputsModel$BasinAreas[upstream_basin] / AreaTot +
+        HUTRANS[1, upstream_basin] +
         c(rep(0, floor(PT[upstream_basin] + 1)),
           Qupstream[(2 + floor(PT[upstream_basin])):LengthTs]) *
-        HUTRANS[2, upstream_basin] * InputsModel$BasinAreas[upstream_basin] / AreaTot
+        HUTRANS[2, upstream_basin]
     }
+    # Convert back Qsim to mm
+    OutputsModelDown$Qsim <- OutputsModelDown$Qsim / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1E3
 
   } else {
 

tests/testthat/test-RunModel_LAG.R

@@ -63,30 +63,22 @@ OutputsGR4JOnly <-
 
 test_that("Upstream basin with nil area should return same Qdown as GR4J alone", {
   UpstBasinArea = InputsModel$BasinAreas[1]
-  InputsModel$BasinAreas[1] <<- 0
+  InputsModel$BasinAreas[1] <- 0
   OutputsSD <-
     RunModel(InputsModel,
              RunOptions,
              Param = c(Param, 1),
              FUN_MOD = RunModel_GR4J)
   expect_equal(OutputsGR4JOnly$Qsim, OutputsSD$Qsim)
-  InputsModel$BasinAreas[1] <<- UpstBasinArea
 })
 
 test_that(
   "Downstream basin with nil area and nul upstream length should return same Qdown as Qupstream alone",
   {
-    InputsModelZeroDown <- 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)
-    )
+    InputsModel$LengthHydro <- 0
+    InputsModel$BasinAreas <- c(BasinInfo$BasinArea, 0)
     OutputsSD <-
-      RunModel(InputsModelZeroDown,
+      RunModel(InputsModel,
                RunOptions,
                Param = c(Param, 1),
                FUN_MOD = RunModel_GR4J)
@@ -95,12 +87,12 @@ test_that(
 )
 
 ParamSD = c(Param, InputsModel$LengthHydro / (24 * 60 * 60)) # Speed corresponding to one time step delay
-OutputsSD <-
-  RunModel(InputsModel, RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
 
 test_that("1 input with lag of 1 time step delay out gives an output delayed of one time step", {
   QlsGR4Only <-
     OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2] * 1E6 / 86400
+  OutputsSD <-
+    RunModel(InputsModel, RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
   QlsSdSim <-
     OutputsSD$Qsim * sum(InputsModel$BasinAreas) * 1E6 / 86400
   QlsUpstLagObs <-
@@ -130,3 +122,23 @@ test_that("Params from calibration with simulated data should be similar to init
     )
   expect_equal(OutputsCalib$ParamFinalR, ParamSD, tolerance = 1E-3)
 })
+
+test_that("1 no area input with lag of 1 time step delay out gives an output delayed of one time step converted to mm", {
+  Qm3GR4Only <-
+    OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2] * 1E3
+  # Specify that upstream flow is not related to an area
+  InputsModel$BasinAreas = c(NA, BasinInfo$BasinArea)
+  # Convert upstream flow to Liter/day
+  InputsModel$Qupstream <- matrix(Qupstream, ncol = 1) * BasinInfo$BasinArea * 2
+  OutputsSD <-
+    RunModel(InputsModel, RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
+  
+  expect_false(any(is.na(OutputsSD$Qsim)))
+  
+  Qm3SdSim <-
+    OutputsSD$Qsim * sum(InputsModel$BasinAreas, na.rm = TRUE) * 1E3
+  Qm3UpstLagObs <-
+    c(0, InputsModel$Qupstream[Ind_Run[1:(length(Ind_Run) - 1)] + 1])
+  
+  expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
+})