mct using average discharge and closing water balance on catchment

src/lisflood/Lisflood_dynamic.py

@@ -218,7 +218,6 @@ def splitlanduse(array1, array2=None, array3=None):
 
         self.ChanQAvg = self.sumDisDay/self.NoRoutSteps
         # Average channel outflow over the model computation step
-        # NOT the same as ChanQ even when if using only one routing sub-step
 
 
         # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -230,12 +229,14 @@ def splitlanduse(array1, array2=None, array3=None):
             # Set water level dynamic wave to dummy value (needed
 
         if option['InitLisflood'] or option['repAverageDis']:
-            self.CumQ += self.ChanQ
+            # self.CumQ += self.ChanQ
+            self.CumQ += self.ChanQAvg
             #cmcheck - we should use ChanQAvg here not ChanQ
             self.avgdis = self.CumQ/self.TimeSinceStart
             # to calculate average discharge over the entire simulation
 
-        self.DischargeM3Out += np.where(self.AtLastPointC ,self.ChanQ * self.DtSec,0)
+        #self.DischargeM3Out += np.where(self.AtLastPointC ,self.ChanQ * self.DtSec,0)
+        self.DischargeM3Out += np.where(self.AtLastPointC, self.ChanQAvg * self.DtSec, 0)
         # Cumulative outflow out of map
         # cmcheck - we should use ChanQAvg here not ChanQ
 

src/lisflood/Lisflood_initial.py

@@ -234,9 +234,7 @@ def __init__(self):
 
         # initialise averaged discharge
         maskinfo = MaskInfo.instance()
-
-        #cmcheck - now it's a state variable and it's initialised elsewhere
-        # self.ChanQAvg = maskinfo.in_zero()
+        self.ChanQAvg = maskinfo.in_zero()
 
         # initialise outputs once everything is set
         self.output_module.initial()

src/lisflood/hydrological_modules/routing.py

@@ -380,20 +380,19 @@ def initial(self):
             # Here it becomes Outflow (x+dx) Q at the beginning of the computation step (t) for full section (instant)
             # Used to calculated Inflow (x) from upstream pixels at the beginning of the computation step (t)
             self.var.ChanQ = np.where(PrevDischarge == -9999, self.var.ChanQKin, PrevDischarge) #np
-            # Initialise channel discharge: cold start: equal to ChanQKin [m3/s]
+            # Initialise instantaneous channel discharge: cold start: equal to ChanQKin [m3/s]
 
             # # cmcheck - deve diventare una variabile di stato
             # PrevAvgDischarge = loadmap('DischargeMaps')
-            # self.var.ChanQAvg = np.where(PrevAvgDischarge == -9999, maskinfo.in_zero(), PrevAvgDischarge) #np
+            # self.var.ChanQAvgDt = np.where(PrevAvgDischarge == -9999, maskinfo.in_zero(), PrevAvgDischarge) #np
             # cmcheck - per adesso
             self.var.ChanQAvgDt =  maskinfo.in_zero()
-
-            # Average channel outflow (x+dx) from channels during previous model computation step (dt)
+            # Initialise average channel outflow (x+dx) from channels during previous model computation step (dt)
 
             # cmcheck
             ChanQAvgDtPcr = decompress(self.var.ChanQAvgDt)  # pcr
             self.var.ChanQAvgInDt = compressArray(upstream(self.var.LddChan, ChanQAvgDtPcr))
-            # Average channel inflow (x) from channels for next model computation step (dt)
+            # Initialise average channel inflow (x) from channels for next model computation step (dt)
             # Using LddChan here because we need all pixels upstream, kin and mct
 
 
@@ -706,21 +705,23 @@ def dynamic(self, NoRoutingExecuted):
                 # Outflow (x+dx) at time t beginning of calculation step (instant)
                 # This is used to calculate inflow from upstream cells
 
-                # ChanM3KinStart = self.var.ChanM3Kin.copy()
+                ChanM3KinStart = self.var.ChanM3Kin.copy()
                 # Channel storage at time t beginning of calculation step (instant)
 
                 ChanQKinOutEnd,ChanM3KinEnd = self.KINRouting(ChanQKinOutStart,SideflowChan)
                 # Outflow (x+dx) at time t+dt end of calculation step (instant)
                 # Channel storage at time t+dt end of calculation step (instant)
+                # This is in fact the same as the average outflow -> from some feature in kinematic routing function
 
                 # cmcheck
-                self.var.ChanQAvgDt = (self.var.ChanQAvgInDt * self.var.DtRouting + SideflowChanM3 - ChanM3KinEnd + self.var.ChanM3) * self.var.InvDtRouting
+                # Calculate average outflow using water balance for channel grid cell over sub-routing step
+                self.var.ChanQAvgDt = (self.var.ChanQAvgInDt * self.var.DtRouting + SideflowChanM3 - ChanM3KinEnd + ChanM3KinStart) * self.var.InvDtRouting
                 # if np.any(self.var.ChanQAvgDt < 0):
                 #     print("At least one value is less than 0.")
                 self.var.ChanQAvgDt[self.var.ChanQAvgDt < 0] = 0
                 # Average outflow (x+dx) QAvg during routing step (average)
-                # Qout_avg = Qin_avg -(Vend - Vini)
-
+                # Qout_avg = Qin_avg -(Vend - Vstart)/DtRouting
+                # Qin_avg = Qin_avg_channels + Qin_avg_sideflow
 
                 # updating variables for next routing step
                 self.var.ChanQKin = ChanQKinOutEnd.copy()
@@ -735,15 +736,16 @@ def dynamic(self, NoRoutingExecuted):
                 # Channel storage V at the end of computation step t+dt for full section (instant)
                 # same as ChanM3KinEnd for Kinematic routing only
 
-                # Update average channel inflow (x) Qavg for next step
-                ChanQAvgInDt = decompress(self.var.ChanQAvgDt)  # pcr
-                self.var.ChanQAvgInDt = compressArray(upstream(self.var.LddChan, ChanQAvgInDt))
+                # Update average channel inflow (x) Qavg for next sub-routing step
+                ChanQAvgDtPcr = decompress(self.var.ChanQAvgDt)  # pcr
+                self.var.ChanQAvgInDt = compressArray(upstream(self.var.LddChan, ChanQAvgDtPcr))
 
                 # self.var.sumDisDay += self.var.ChanQ
                 # # sum of total river outflow on model sub-step
 
                 self.var.sumDisDay += self.var.ChanQAvgDt
-                # sum of average river outflow on model routing sub-step
+                # sum of average river outflow on routing sub-step
+                # used for calculating average discharge on model time step
 
 
             # SPLIT ROUTING - no InitLisfllod
@@ -775,29 +777,27 @@ def dynamic(self, NoRoutingExecuted):
 
                 if not (option['SplitRouting']):
                     # Kinematic routing
+                    # This is calculated for every grid cell, including MCT grid cells
 
                     ChanQKinOutStart = self.var.ChanQ.copy()
                     # Outflow (x+dx) Q at time t beginning of calculation step (instant)
+                    # This is used to calculate inflow from upstream cells
 
                     ChanM3KinStart = self.var.ChanM3.copy()
-
-                    # ChanM3KinStart = self.var.ChanM3.copy()
-                    # Channel storage at time t (instant)
+                    # Channel storage at time t beginning of calculation step (instant)
 
                     ChanQKinOutEnd,ChanM3KinEnd = self.KINRouting(ChanQKinOutStart,SideflowChan)
                     # Outflow at time t+dt end of calculation step (instant)
                     # Channel storage at time t+dt end of calculation step (instant)
-
+                    # This is in fact the same as the average outflow -> from some feature in kinematic routing function
 
                     # cmcheck
-
-                    # cal average outflow discharge for channels during routing sub step dt
-
+                    # Calculate average outflow using water balance for channel grid cell
                     ChanQKinAvgDt = (self.var.ChanQAvgInDt * self.var.DtRouting + SideflowChanM3 - ChanM3KinEnd + ChanM3KinStart) * self.var.InvDtRouting
                     # if np.any(self.var.ChanQAvgDt < 0):
                     #     print("At least one value is less than 0.")
                     ChanQKinAvgDt[ChanQKinAvgDt < 0] = 0
-                    # Average outflow (x+dx) QAvg during routing step (average)
+                    # Average outflow (x+dx) QAvg during sub-routing step (average)
                     # Qout_avg = Qin_avg -(Vend - Vini)
 
 
@@ -822,7 +822,7 @@ def dynamic(self, NoRoutingExecuted):
                 # MCT routing
 
                 ChanQMCTOutStart = self.var.ChanQ.copy()
-                # Outflow (x+dx) at time t  q10 (instant)
+                # Outflow (x+dx) at time t (instant)  q10
 
                 # debug
                 # ChanQKinOutEnd[ChanQKinOutEnd != 0] = 0
@@ -832,7 +832,7 @@ def dynamic(self, NoRoutingExecuted):
                 # Channel storage at time t V00
 
                 ChanQMCTInStart = self.var.PrevQMCTin.copy()   
-                # Inflow (x) at time t instant (instant)
+                # Inflow (x) at time t instant (instant)  q00
                 # This is coming from upstream pixels
 
                 # calling MCT routing
@@ -842,10 +842,11 @@ def dynamic(self, NoRoutingExecuted):
                 # Channel storage at time t+dt end of calculation step (instant)
 
 
-                # # update input (x) Q at t for next step (instant)
-                # ChanQMCTStartPcr = decompress(ChanQMCTOutStart)  # pcr
-                # self.var.PrevQMCTin = compressArray(upstream(self.var.LddChan, ChanQMCTStartPcr))
+                # # update input (x) Q at t for next step (instant)  q10 -> q00
+                # ChanQMCTOutStartPcr = decompress(ChanQMCTOutStart)  # pcr
+                # self.var.PrevQMCTin = compressArray(upstream(self.var.LddChan, ChanQMCTOutStartPcr))
                 # # using LddChan here because we need to input from upstream pixels to include kinematic pixels
+                # update input (x) Q at t for next step (instant)
 
 
                 # Storing MCT Courant and Reynolds numbers for state files
@@ -855,6 +856,7 @@ def dynamic(self, NoRoutingExecuted):
                 # cmcheck
 
                 # calc average discharge outflow for MCT channels during routing sub step dt
+                # Calculate average outflow using water balance for MCT channel grid cell over sub-routing step
                 # ChanQMCTAvgDt = (self.var.ChanQAvgInDt * self.var.DtRouting + SideflowChanMCT * self.var.DtRouting - ChanM3MCTEnd + ChanM3MCTStart) * self.var.InvDtRouting
                 ChanQMCTAvgDt = np.where(self.var.IsChannelKinematic, 0., (
                             self.var.ChanQAvgInDt * self.var.DtRouting + SideflowChanMCT * self.var.DtRouting - ChanM3MCTEnd + ChanM3MCTStart) * self.var.InvDtRouting)
@@ -884,7 +886,6 @@ def dynamic(self, NoRoutingExecuted):
                 self.var.PrevQMCTin = compressArray(upstream(self.var.LddChan, ChanQMCTStartPcr))
                 # using LddChan here because we need to input from upstream pixels to include kinematic pixels
 
-
                 # Update average channel inflow (x) Qavg for next step
                 ChanQAvgInDtPcr = decompress(self.var.ChanQAvgDt)  # pcr
                 self.var.ChanQAvgInDt = compressArray(upstream(self.var.LddChan, ChanQAvgInDtPcr))

src/lisflood/hydrological_modules/waterbalance.py

@@ -213,6 +213,7 @@ def dynamic(self):
 
             ##sum1 = self.var.sumDis.copy()
             sum1 = self.var.ChanQAvg.copy()
+            # averge outflow during model time step
             sum1[self.var.AtLastPointC == 0] = 0
             WaterOut = np.take(np.bincount(self.var.Catchments,weights=sum1 * self.var.DtSec),self.var.Catchments)
             # Water that goes out of the system at the channels level [m3]