Solving issue with numbe in discharge_avg in solve1pixel

src/lisflood/hydrological_modules/kinematic_wave_parallel_tools.py

@@ -97,6 +97,7 @@ def solve1Pixel(pix, discharge_avg, discharge, lateral_inflow, constant,\
     for ups_ix in range(num_upstream_pixels[pix]):
         upstream_inflow += discharge[upstream_lookup[pix,ups_ix]]
         upstream_inflow_avg += discharge_avg[upstream_lookup[pix, ups_ix]]
+
     const_plus_ups_infl = upstream_inflow + constant[pix] # upstream_inflow + alpha*dx/dt*Qold**beta + dx*specific_lateral_inflow
     # If old discharge, upstream inflow and lateral inflow are below accuracy: set discharge to 0 and exit
     if const_plus_ups_infl <= NEWTON_TOL:
@@ -122,15 +123,15 @@ def solve1Pixel(pix, discharge_avg, discharge, lateral_inflow, constant,\
     if discharge[pix] == NEWTON_TOL:
         discharge[pix] = 0
 
-    # cmcheck
-    # avoid negative discharge
-    if discharge[pix] < 0:
-        discharge[pix] = 0
     # volume of water in channel at end of computation step
     channel_volume_end = a_dx * discharge[pix]**beta
+
     # mass water balance to calc average outflow
-    discharge_avg[pix] = upstream_inflow_avg + lateral_inflow + (channel_volume_start - channel_volume_end) * inv_time_delta
+    discharge_avg[pix] = upstream_inflow_avg + lateral_inflow[pix] + (channel_volume_start[pix] - channel_volume_end[pix]) * inv_time_delta
 
+    # avoid negative average discharge
+    if discharge_avg[pix] < 0:
+        discharge_avg[pix] = 0
 
     # to simulate inf or nan: discharge[pix] = 1.0/0.0
     # with gil:

src/lisflood/hydrological_modules/surface_routing.py

@@ -94,6 +94,12 @@ def initial(self):
         self.var.OFQOther = ((self.var.OFM3Other * self.var.InvPixelLength * self.var.InvOFAlpha.values[self.var.dim_runoff[1].index('Other')])**(self.var.InvBeta)).astype(float)
         self.var.OFQForest = ((self.var.OFM3Forest * self.var.InvPixelLength * self.var.InvOFAlpha.values[self.var.dim_runoff[1].index('Forest')])**(self.var.InvBeta)).astype(float)
 
+        # cmcheck
+        # Initial average overland discharge [m3 s-1]
+        self.var.OFQDirect_avg = maskinfo.in_zero()
+        self.var.OFQOther_avg = maskinfo.in_zero()
+        self.var.OFQForest_avg = maskinfo.in_zero()
+
     def initialSecond(self):
         """ 2nd initialisation part of the surface routing module:
             to be called after all needed parameters are set (PixelLength, LddToChan, Alpha...)
@@ -148,13 +154,10 @@ def dynamic(self):
         SideflowOther = np.sum(self.var.SurfaceRunSoil.values[ilusevalues],self.var.SurfaceRunSoil.dims.index("landuse")) * self.var.MMtoM3 * self.var.InvPixelLength * self.var.InvDtSec
         SideflowForest = self.var.SurfaceRunSoil.values[self.var.epic_settings.soil_uses.index('Forest')] * self.var.MMtoM3 * self.var.InvPixelLength * self.var.InvDtSec
         # All surface runoff that is generated during current time step added as side flow [m3/s/m pixel-length]
-        # cmcheck
-        OFQDirect_avg = 0.
-        OFQOther_avg = 0.
-        OFQForest_avg = 0.
-        self.direct_surface_router.kinematicWaveRouting(OFQDirect_avg,self.var.OFQDirect, SideflowDirect)
-        self.other_surface_router.kinematicWaveRouting(OFQOther_avg, self.var.OFQOther, SideflowOther)
-        self.forest_surface_router.kinematicWaveRouting(OFQForest_avg, self.var.OFQForest, SideflowForest)
+
+        self.direct_surface_router.kinematicWaveRouting(self.var.OFQDirect_avg,self.var.OFQDirect, SideflowDirect)
+        self.other_surface_router.kinematicWaveRouting(self.var.OFQOther_avg, self.var.OFQOther, SideflowOther)
+        self.forest_surface_router.kinematicWaveRouting(self.var.OFQForest_avg, self.var.OFQForest, SideflowForest)
 
 # to PCRASTER