adjusting checks for small q values

Author: Cinzia Mazzetti

src/lisflood/hydrological_modules/routing.py

@@ -1381,7 +1381,7 @@ def MCTRouting_single(self, V00, q10, q01, q00, ql,q0mm, Cm0, Dm0, dt, xpix, s0,
 
         # check for negative and zero discharge values
         # zero outflow is not allowed
-        if q11 <= 0:
+        if q11 < eps: # cmcheck <=0
             q11 = eps
 
         # calc reference discharge at time t
@@ -1393,13 +1393,13 @@ def MCTRouting_single(self, V00, q10, q01, q00, ql,q0mm, Cm0, Dm0, dt, xpix, s0,
 
             # reference I discharge at x=0
             qmx0 = (q00 + q01) / 2.
-            if qmx0 == 0:
+            if qmx0 < eps: # cmcheck ==0
                 qmx0 = eps
             hmx0 = self.hoq(qmx0, s0, Balv, ANalv, Nalv)
 
             # reference O discharge at x=1
             qmx1 = (q10 + q11) / 2.
-            if qmx1 == 0:
+            if qmx1 < eps: # cmcheck ==0
                 qmx1 = eps
             hmx1 = self.hoq(qmx1, s0,Balv,ANalv,Nalv)
 
@@ -1413,7 +1413,7 @@ def MCTRouting_single(self, V00, q10, q01, q00, ql,q0mm, Cm0, Dm0, dt, xpix, s0,
             # Q(t+dt)=(I(t+dt)+O'(t+dt))/2
             qm1 = (q01 + q11) / 2.
             #cm
-            if qm1 == 0:
+            if qm1 < eps: # cmcheck ==0
                 qm1 = eps
             #cm
             hm1 = self.hoq(qm1,s0,Balv,ANalv,Nalv)
@@ -1442,7 +1442,7 @@ def MCTRouting_single(self, V00, q10, q01, q00, ql,q0mm, Cm0, Dm0, dt, xpix, s0,
             # Mass balance equation that takes into consideration the lateral flow
             q11 = c1 * q01 + c2 * q00 + c3 * q10 + c4 * ql
 
-            if q11 <= 0.:
+            if q11 < eps: # cmcheck <=0
                 q11 = eps
 
             #### end of for loop
@@ -1475,7 +1475,7 @@ def MCTRouting_single(self, V00, q10, q01, q00, ql,q0mm, Cm0, Dm0, dt, xpix, s0,
         q1mm = q0mm + ql + (V00 - V11) / dt
 
         # cmcheck
-        # q1m cannot be too small or it will cause instability
+        # q1m cannot be smaller than eps or it will cause instability
         if q1mm < eps:
             q1mm = eps
             V11 = V00 + (q0mm + ql - q1mm) * dt
@@ -1486,6 +1486,7 @@ def MCTRouting_single(self, V00, q10, q01, q00, ql,q0mm, Cm0, Dm0, dt, xpix, s0,
         return q11, V11, q1mm, Cm1, Dm1
 
 
+
     def hoq(self,q,s0,Balv,ANalv,Nalv):
         """Water depth from discharge.
         Given a generic cross-section (rectangular, triangular or trapezoidal) and a steady-state discharge q=Q*, it computes