From 313d00cc74a950d6eba3ae0ddd5c475f1b9baa45 Mon Sep 17 00:00:00 2001 From: Cinzia Mazzetti Date: Thu, 18 Apr 2024 14:33:17 +0000 Subject: [PATCH] mct using average discharge and closing water balance on catchment --- src/lisflood/Lisflood_dynamic.py | 7 ++- src/lisflood/Lisflood_initial.py | 4 +- src/lisflood/hydrological_modules/routing.py | 55 ++++++++++--------- .../hydrological_modules/waterbalance.py | 1 + 4 files changed, 34 insertions(+), 33 deletions(-) diff --git a/src/lisflood/Lisflood_dynamic.py b/src/lisflood/Lisflood_dynamic.py index 29548348..d547ce81 100644 --- a/src/lisflood/Lisflood_dynamic.py +++ b/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 diff --git a/src/lisflood/Lisflood_initial.py b/src/lisflood/Lisflood_initial.py index 76d2bfa5..f238251e 100644 --- a/src/lisflood/Lisflood_initial.py +++ b/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() diff --git a/src/lisflood/hydrological_modules/routing.py b/src/lisflood/hydrological_modules/routing.py index 9b94da8b..8e03898e 100644 --- a/src/lisflood/hydrological_modules/routing.py +++ b/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)) diff --git a/src/lisflood/hydrological_modules/waterbalance.py b/src/lisflood/hydrological_modules/waterbalance.py index aed1dedd..65d39b71 100755 --- a/src/lisflood/hydrological_modules/waterbalance.py +++ b/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]