From ac3a4626d40de94f59945603a11680553397bad5 Mon Sep 17 00:00:00 2001
From: Dorchies David <david.dorchies@inrae.fr>
Date: Tue, 2 Mar 2021 11:31:46 +0100
Subject: [PATCH] fix: RunModel_Lag initial states unit conversion issue

Refs #19, HYCAR-Hydro/airgr#104
---
 R/RunModel_Lag.R               | 13 ++++++++++++-
 tests/testthat/test-RunModel.R |  3 ++-
 2 files changed, 14 insertions(+), 2 deletions(-)

diff --git a/R/RunModel_Lag.R b/R/RunModel_Lag.R
index e3bb075..b625ac1 100644
--- a/R/RunModel_Lag.R
+++ b/R/RunModel_Lag.R
@@ -75,6 +75,7 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
       rep(0, floor(PT[x] + 1))
     })
   }
+  #message("Initstates: ",paste(IniStates, collapse = ", "))
 
   for (upstream_basin in seq_len(NbUpBasins)) {
     Qupstream <- InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin]
@@ -85,6 +86,9 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
     Qupstream1 <- c(IniStates[[upstream_basin]][-length(IniStates[[upstream_basin]])], Qupstream[1:(LengthTs - floor(PT[upstream_basin]))])
     Qupstream2 <- IniStates[[upstream_basin]]
     if(LengthTs - floor(PT[upstream_basin]) - 1 > 0) Qupstream2 <- c(Qupstream2, Qupstream[1:(LengthTs - floor(PT[upstream_basin]) - 1)])
+    #message("Qupstream1: ", paste(Qupstream1, collapse = ", "))
+    #message("Qupstream2: ", paste(Qupstream2, collapse = ", "))
+
     OutputsModel$Qsim <- OutputsModel$Qsim +
                          Qupstream1 * HUTRANS[1, upstream_basin] +
                          Qupstream2 * HUTRANS[2, upstream_basin]
@@ -96,12 +100,19 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
   }
   # 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) {
       LengthTs <- tail(RunOptions$IndPeriod_Run,1)
-      InputsModel$Qupstream[(LengthTs - floor(PT[x])):LengthTs, x]
+      Qupstream <- InputsModel$Qupstream[(LengthTs - floor(PT[x])):LengthTs, x]
+      if (!is.na(InputsModel$BasinAreas[x])) {
+        # Upstream flow with area needs to be converted to m3 by time step
+        Qupstream <- Qupstream * InputsModel$BasinAreas[x] * 1e3
+      }
+      return(Qupstream)
     })
+    #message("StateEnd: ",paste(OutputsModel$StateEnd$SD, collapse = ", "))
   }
 
   return(OutputsModel)
diff --git a/tests/testthat/test-RunModel.R b/tests/testthat/test-RunModel.R
index b12867f..e32f931 100644
--- a/tests/testthat/test-RunModel.R
+++ b/tests/testthat/test-RunModel.R
@@ -17,8 +17,9 @@ test_that("RunModelSupervisor with no regulation should returns same results as
   Qobs <- cbind(sapply(Severn$BasinsObs, function(x) {x$discharge_spec}))
   InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
   # RunOptions
+  nTS <- 365
   IndPeriod_Run <- seq(
-    length(InputsModel[[1]]$DatesR) - 365,
+    length(InputsModel[[1]]$DatesR) - nTS + 1,
     length(InputsModel[[1]]$DatesR)
   )
   IndPeriod_WarmUp = seq(IndPeriod_Run[1]-366,IndPeriod_Run[1]-1)
-- 
GitLab