Bugfix/mct lake

VERSION

@@ -1 +1 @@
-4.3.1.502
+4.3.1.503

src/lisflood/hydrological_modules/lakes.py

@@ -239,6 +239,14 @@ def dynamic_inloop(self, NoRoutingExecuted):
             # SI = S1/dtime - Qout1/2 + (Qin1 + Qin2)/2 = S1/dtime - Qout1/2 + LakeIn
             # Right hand part of the Modified Puls Method equation above
 
+            # Check LakeOutflowCC for negative values
+            if any(LakeStorageIndicator < 0):
+                # this can happen with large oscillations in the lake inflow
+                msg = "Negative or NaN outflow from lakes. " \
+                      "Consider increasing computation time step for routing (DtSecChannel) \n"
+                warnings.warn(LisfloodWarning(msg))
+                LakeStorageIndicator[LakeStorageIndicator < 0] = 0
+
             # Qout2
             self.var.LakeOutflowCC = np.square( -self.var.LakeFactor + np.sqrt(self.var.LakeFactorSqr + 2 * LakeStorageIndicator))
             # Qout2 average lake outflow in [m3/s] per timestep t (routing sub-step)
@@ -286,7 +294,7 @@ def dynamic_inloop(self, NoRoutingExecuted):
                     self.var.sumLakeInCC = self.var.LakeInflowCC * self.var.DtRouting
                     self.var.sumLakeOutCC = QLakeOutM3DtCC
                     # for timeseries output - in and outflow to the reservoir
-                    # is sumed up over the sub timesteps and stored in m/s
+                    # is summed up over the sub timesteps and stored in m/s
                     # set to zero at first timestep
                 else:
                     self.var.sumLakeInCC += self.var.LakeInflowCC * self.var.DtRouting

src/lisflood/hydrological_modules/mct.py

@@ -339,10 +339,11 @@ def MCTRouting_single(
     # q1m cannot be smaller than eps or it will cause instability
     if q1mm < 0:   # cmcheck <=0
         q1mm = 0
-        if ql < 0: ql = 0
+        ###if ql < 0: ql = 0
         # prevent water abstraction or open water evaporation from drying out the channel and keep extracting water
         # NOTE THIS GENERATES AN ERROR IN THE WATER BALANCE
         V11 = V00 + (q0mm + ql - q1mm) * dt
+        if V11 < 0: V11 = 0
 
     # q11 Outflow at O(t+dt)
     # q1m average outflow in time dt

src/lisflood/hydrological_modules/waterabstraction.py

@@ -10,6 +10,7 @@
 """
 from __future__ import absolute_import, print_function
 
+import os
 import warnings
 
 from pcraster import boolean, nominal, ifthen, defined, areamaximum, downstream, cover, lddrepair, ifthenelse, upstream, \
@@ -250,6 +251,14 @@ def initial(self):
             self.var.FractionAbstractedFromChannels = maskinfo.in_zero()
             self.var.areatotal_abstraction_SW_actual_irrigation_M3 = maskinfo.in_zero()
             self.var.areatotal_withdrawal_SW_actual_M3 = maskinfo.in_zero()
+
+            #### Prepare output file to store info on potential issue for Paddy Rice Water Abstraction
+            if option.get('repPaddyRiceDebug') is True:
+                self.debug_paddyrice_filename = os.path.join(settings.output_dir,f"debug_paddyrice.csv")
+                if os.path.exists(self.debug_paddyrice_filename):
+                    os.remove(self.debug_paddyrice_filename)
+
+
 
 
     def dynamic(self):
@@ -678,8 +687,33 @@ def dynamic(self):
 
             self.var.Theta1a.values[iveg] = self.var.W1a.values[iveg] / self.var.SoilDepth1a.values[ilanduse] 
             self.var.Theta1b.values[iveg] = self.var.W1b.values[iveg] / self.var.SoilDepth1b.values[ilanduse]
-
 
+            if option.get('repPaddyRiceDebug') is True:
+                # ************************************************************
+                # 20. Check if we abstracted more water for PaddyRice then available water, and write a txt file (known issue, to be fixed)
+                # (issue for Paddy Rice Water Abstraction)
+                # ************************************************************
+                unique_regions, index = np.unique(self.var.WUseRegionC, return_index=True)
+                areatotal_PaddyRiceWaterAbstractionFromSurfaceWaterM3 = np.take(np.bincount(self.var.WUseRegionC, weights=self.var.PaddyRiceWaterAbstractionFromSurfaceWaterM3), self.var.WUseRegionC)
+                for region, idx in zip(unique_regions, index):
+                    if (areatotal_PaddyRiceWaterAbstractionFromSurfaceWaterM3[idx]>0.0) and (self.var.areatotal_withdrawal_LakRes_actual_M3[idx] + \
+                                                                                             self.var.AreaTotalAvailableWaterFromChannelsM3[idx] - \
+                                                                                             areatotal_withdrawal_SW_required[idx] < 0.0):
+                        header = None
+                        if not os.path.exists(self.debug_paddyrice_filename):
+                            header = "Step, Region, areatotal_PaddyRiceWaterAbstractionFromSurfaceWaterM3, " \
+                                     "areatotal_withdrawal_LakRes_actual_M3, AreaTotalAvailableWaterFromChannelsM3, " \
+                                     "areatotal_withdrawal_SW_required, areatotal_withdrawal_SW_except_PaddyRice_required\n"
+                        message = f"{self.var.currentStep}, {region}, {areatotal_PaddyRiceWaterAbstractionFromSurfaceWaterM3[idx]}, " \
+                                  f"{self.var.areatotal_withdrawal_LakRes_actual_M3[idx]}, {self.var.AreaTotalAvailableWaterFromChannelsM3[idx]}, " \
+                                  f"{areatotal_withdrawal_SW_required[idx]}, {areatotal_withdrawal_SW_required[idx]-areatotal_PaddyRiceWaterAbstractionFromSurfaceWaterM3[idx]}\n"
+                        #print(message)
+                        # Open the file in write mode
+                        with open(self.debug_paddyrice_filename, 'a') as file:
+                            if header is not None:
+                                file.write(header)
+                            # Write the message to the file
+                            file.write(message)
 
 #from numba import njit
 #from builtins import max, min

src/lisfloodSettings_reference.xml

@@ -67,6 +67,8 @@ You can use builtin path variables in this template and reference to other paths
       <setoption choice="1"  name="writeNetcdfStack"/>
       <setoption choice="1"  name="writeNetcdf"/>
 
+      # option to Report PaddyRiceDebug information in CSV file for debugging Paddy Rice Water Abstraction issue on dry Channels
+		  <setoption choice="0"  name="repPaddyRiceDebug"/>
 
       #-----------------------------------------------------------
       # report time series