diff --git a/VERSION b/VERSION
index 7ac554f1..58ca1887 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-4.3.1.502
\ No newline at end of file
+4.3.1.503
\ No newline at end of file
diff --git a/src/lisflood/hydrological_modules/lakes.py b/src/lisflood/hydrological_modules/lakes.py
index 234b592e..9134bfba 100755
--- a/src/lisflood/hydrological_modules/lakes.py
+++ b/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
diff --git a/src/lisflood/hydrological_modules/mct.py b/src/lisflood/hydrological_modules/mct.py
index 900d9b42..99f0ad22 100644
--- a/src/lisflood/hydrological_modules/mct.py
+++ b/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
diff --git a/src/lisflood/hydrological_modules/waterabstraction.py b/src/lisflood/hydrological_modules/waterabstraction.py
index 6ebcb2c4..350d13ca 100644
--- a/src/lisflood/hydrological_modules/waterabstraction.py
+++ b/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
diff --git a/src/lisfloodSettings_reference.xml b/src/lisfloodSettings_reference.xml
index 297c42d4..1cec854d 100644
--- a/src/lisfloodSettings_reference.xml
+++ b/src/lisfloodSettings_reference.xml
@@ -67,6 +67,8 @@ You can use builtin path variables in this template and reference to other paths
+ # option to Report PaddyRiceDebug information in CSV file for debugging Paddy Rice Water Abstraction issue on dry Channels
+
#-----------------------------------------------------------
# report time series