From 9eab0a15e2f2b8c3666c9df7187eba8e8db9ebd6 Mon Sep 17 00:00:00 2001
From: Dorchies David <david.dorchies@inrae.fr>
Date: Tue, 2 Aug 2022 15:46:02 +0200
Subject: [PATCH] refactor(RunModel_LRQ): Use non simplified formula of with
 tanh from Bentura & Michel (1997)

Refs HYCAR-Hydro/airgr#152, HYCAR-Hydro/airgr#153
---
 R/RunModel_LRQ.R | 14 +++++++++-----
 1 file changed, 9 insertions(+), 5 deletions(-)

diff --git a/R/RunModel_LRQ.R b/R/RunModel_LRQ.R
index 46bcd68..8bf801b 100644
--- a/R/RunModel_LRQ.R
+++ b/R/RunModel_LRQ.R
@@ -65,10 +65,10 @@ RunModel_LRQ <- function(InputsModel, RunOptions, Param, QcontribDown) {
 
 
   ## Parameters set up
-  TParam <- Param[1L]
-  KParam <- Param[2L] # KParam en mm/pdt
+  CParam <- Param[1L] # CParam in m/s
+  KParam <- Param[2L] # KParam in mm/TS
 
-  PT <-  InputsModel$LengthHydro * 1e3 / TParam / RunOptions$FeatFUN_MOD$TimeStep
+  PT <-  InputsModel$LengthHydro * 1e3 / CParam / RunOptions$FeatFUN_MOD$TimeStep #PT in h
   PK <- sqrt(KParam * InputsModel$BasinAreas[-length(InputsModel$BasinAreas)] * 1e-3) * PT #m3/2 pdt1/2
 
   ## Lag set up
@@ -119,14 +119,18 @@ RunModel_LRQ <- function(InputsModel, RunOptions, Param, QcontribDown) {
   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 in seq_along(Qroute)[1:(length(Qroute) - 1)]){
-      Qroute[q + 1] <- ((Qroute[q]^0.5 + (((Qupstream_m3[q] + Qupstream_m3[q + 1])/2)/PK[upstream_basin])) /
-                          (1 + ((Qroute[q]^0.5) / PK[upstream_basin])))^2
+      Qroute[q + 1] <- ((Qroute[q]^0.5 + ((Qupstream_m3[q] + Qupstream_m3[q + 1])/2)^0.5 *
+                          tanh((((Qupstream_m3[q] + Qupstream_m3[q + 1])/2)^0.5)/PK[upstream_basin]) ) /
+                         (1 + (Qroute[q]/((Qupstream_m3[q] + Qupstream_m3[q + 1])/2))^0.5 *
+                          tanh((((Qupstream_m3[q] + Qupstream_m3[q + 1])/2)^0.5)/PK[upstream_basin])) )^2
     }
 
     Qllr <- Lag(PT[upstream_basin], Qroute, Qupstream_m3)[IndPeriod1]
-- 
GitLab