diff --git a/VERSION b/VERSION
index c3ef9b53..9286fc2e 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-4.3.1.505
\ No newline at end of file
+4.3.1.506
\ No newline at end of file
diff --git a/src/lisflood/Lisflood_initial.py b/src/lisflood/Lisflood_initial.py
index 346c1dff..caa48104 100644
--- a/src/lisflood/Lisflood_initial.py
+++ b/src/lisflood/Lisflood_initial.py
@@ -138,6 +138,7 @@ def __init__(self):
self.surface_routing_module = surface_routing(self)
self.reservoir_module = Reservoir(self) # get_reservoir(option['reservoirHanazaki'])
self.lakes_module = lakes(self)
+
self.polder_module = polder(self)
self.waterabstraction_module = waterabstraction(self)
self.indicatorcalc_module = indicatorcalc(self)
@@ -189,10 +190,12 @@ def __init__(self):
self.snow_module.initial()
self.frost_module.initial()
+
self.leafarea_module.initial()
- self.soilloop_module.initial()
+ self.soilloop_module.initial()
self.soil_module.initial()
+
self.routing_module.initial()
self.groundwater_module.initial()
@@ -201,6 +204,8 @@ def __init__(self):
self.inflow_module.initial()
self.surface_routing_module.initial()
+ # At this point LddChan and LddKinematic do not have any structure reservoirs/lakes MCT confluence
+
self.reservoir_module.initial()
self.lakes_module.initial()
self.polder_module.initial()
@@ -209,18 +214,30 @@ def __init__(self):
self.structures_module.initial()
# Structures such as reservoirs and lakes are modelled by interrupting the channel flow paths
+ # At this point LddKinematic and LddChan have pits upstream of (structures) reservoirs and lakes
# ----------------------------------------------------------------------
# ----------------------------------------------------------------------
self.routing_module.initialSecond()
# CHANNEL INITIAL SPLIT UP IN SECOND CHANNEL
+
self.surface_routing_module.initialSecond()
+ # MCT confluence must be in the LddKinematic when I get here
+ self.routing_module.initialKinematicWave()
+
if option.get('MCTRouting'):
self.routing_module.initialMCT()
# initialising Muskingum-Cunge-Todini routing for channel
+ # self.routing_module.initialKinematicWave()
+ # this cannot be here because I need river_router in the MCT initialization
+
+
+
+
+
self.evapowater_module.initial()
self.riceirrigation_module.initial()
self.waterabstraction_module.initial()
@@ -397,4 +414,4 @@ def defsoil(self, name_1, name_2=None, name_3=None, coords=None):
def deffraction(self, variable):
"""Weighted sum over the soil fractions of each pixel"""
ax_veg = variable.dims.index("vegetation")
- return _vegSum(ax_veg, variable.values, self.SoilFraction.values)
+ return _vegSum(ax_veg, variable.values, self.SoilFraction.values)
\ No newline at end of file
diff --git a/src/lisflood/global_modules/default_options.py b/src/lisflood/global_modules/default_options.py
index 94181f80..aa2a9968 100644
--- a/src/lisflood/global_modules/default_options.py
+++ b/src/lisflood/global_modules/default_options.py
@@ -1288,6 +1288,7 @@
'repwateruseGauges': False,
'repwateruseSites': False,
'riceIrrigation': False,
+ 'simulateCalibrationPoints': False,
'simulateLakes': False,
'simulatePF': False,
'simulatePolders': False,
diff --git a/src/lisflood/hydrological_modules/mct.py b/src/lisflood/hydrological_modules/mct.py
index 99f0ad22..2a39f339 100644
--- a/src/lisflood/hydrological_modules/mct.py
+++ b/src/lisflood/hydrological_modules/mct.py
@@ -20,11 +20,13 @@ def __init__(
dt, # computation time step for routing [s]
river_router, # class
mapping_mct, # MCT pixels mapping
+ CalibPointsIds, # calibration points pixels ids
):
# Process flow direction matrix: downstream and upstream lookups, and routing orders
flow_dir = decodeFlowMatrix(rebuildFlowMatrix(compressed_encoded_ldd, land_mask))
self.downstream_lookup, self.upstream_lookup = streamLookups(flow_dir, land_mask)
+ # Inside streamLookups each land_mask pixel is assigned a unique id from 0 to numpix-1, numbered by row
self.num_upstream_pixels = (self.upstream_lookup != -1).sum(1).astype(int) # astype for cython import in windows (to avoid 'long long' buffer dtype mismatch)
# Routing order: decompose domain into batches; within each batch, pixels can be routed in parallel
self._setMCTRoutingOrders()
@@ -37,7 +39,7 @@ def __init__(
self.dt = dt
self.river_router = river_router
self.mapping_mct = mapping_mct
-
+ self.CalibPointsIds = CalibPointsIds
def _setMCTRoutingOrders(self):
"""Compute the MCT wave routing order. Pixels are grouped in sets with the same order.
@@ -103,6 +105,7 @@ def routing(
PrevCm0, # Courant number in input: at time t; in output: at time t+dt
PrevDm0, # Reynolds number in input: at time t; in output: at time t+dt
ChanM3, # V11 as output
+ self.CalibPointsIds,# inflow points used by the calibration suite
)
@@ -130,6 +133,7 @@ def mct_routing(
PrevCm0, # Courant number in input: at time t; in output: at time t+dt
PrevDm0, # Reynolds number in input: at time t; in output: at time t+dt
ChanM3, # V11 as output
+ CalibPointsIds, # inflow points used by the calibration suite
):
"""This function implements Muskingum-Cunge-Todini routing method
MCT routing is calculated on MCT pixels only but gets inflow from both Kinematic/Split and MCT upstream pixels.
@@ -167,11 +171,22 @@ def mct_routing(
q00 = 0.0
q0m = 0.0
q01 = 0.0
+ ql = SideflowChanMCT[kinpix] # ← ql defined BEFORE the loop
for ups_ix in range(num_upstream_pixels[kinpix]):
ups_pix = upstream_pixels[ups_ix] # upstream pixel id
- q00 += ChanQ_0[ups_pix] # Inflow (x) to the pixel at previous step t (instant)
- q0m += ChanQAvgDt[ups_pix] # Average inflow (x) to the pixel at previous step t (average)
- q01 += ChanQ[ups_pix] # Inflow (x) at current step t+dt (instant)
+ #####################################################################################################
+ # This is necessary for EFAS6/GloFAs5 calibration
+ if np.any(CalibPointsIds == ups_pix):
+ # this upstream pixel is a calibration point - add to sideflow
+ ql += ChanQAvgDt[ups_pix] # avoid +=
+ # Sideflow during step dt including contribution from calibration pixel
+ else:
+ # not a calibration point - go as usual
+ q00 += ChanQ_0[ups_pix] # Inflow (x) to the pixel at previous step t (instant)
+ q0m += ChanQAvgDt[ups_pix] # Average inflow (x) to the pixel at previous step t (average)
+ q01 += ChanQ[ups_pix] # Inflow (x) at current step t+dt (instant)
+
+ #####################################################################################################
# get outflow from the pixel at previous step t
q10 = ChanQ_0[kinpix] # Outflow (x+dx) from the pixel at previous step t (instant)
@@ -181,7 +196,7 @@ def mct_routing(
Cm0 = PrevCm0[kinpix] # Courant number at the end of previous step t
Dm0 = PrevDm0[kinpix] # Reynolds number at the end of previous step t
- ql = SideflowChanMCT[kinpix] # Sideflow during step dt
+ # ql = SideflowChanMCT[kinpix] # Sideflow during step dt
# static data
xpix = ChanLength[kinpix] # Channel length
@@ -239,7 +254,7 @@ def MCTRouting_single(
# Calc O' first guess for the outflow at time t+dt
# O'(t+dt)=O(t)+(I(t+dt)-I(t))
- q11 = q10 + (q01 - q00)
+ q11 = q10 + (q01 - q00) # + ql
# check for negative and zero discharge values
# zero outflow is not allowed
@@ -250,6 +265,11 @@ def MCTRouting_single(
# qm0 = (I(t)+O(t))/2
# qm0 = (q00 + q10) / 2.
+ # ql correction for hydraulic state computation
+ # Only positive ql contributes (ignore abstraction),
+ # /2 represents uniform distribution along reach
+ # ql_correction = max(ql, 0.0) / 2.0
+
# Calc O(t+dt)=q11 at time t+dt using MCT equations
for i in range(2): # repeat 2 times for accuracy
@@ -273,7 +293,7 @@ def MCTRouting_single(
# Calc reference discharge time t+dt
# Q(t+dt)=(I(t+dt)+O'(t+dt))/2
- qm1 = (q01 + q11) / 2.0
+ qm1 = (q01 + q11) / 2.0 # + ql_correction
# cm
if qm1 <= eps : # cmcheck ==0 #tpk
qm1 = eps #tpk
@@ -291,10 +311,24 @@ def MCTRouting_single(
Cm1 = ck1 * dt / xpix / Beta1
# Calc MCT parameters
+ # Guard Diffusivity
+ # Dm1 = min(max(Dm1, 0.0), 1.0)
den = 1 + Cm1 + Dm1
+
c1 = (-1 + Cm1 + Dm1) / den
- c2 = (1 + Cm0 - Dm0) / den * (Cm1 / Cm0)
- c3 = (1 - Cm0 + Dm0) / den * (Cm1 / Cm0)
+
+ # Guard against Cm0 near zero (first timestep or after dry conditions)
+ if Cm0 > eps:
+ cm_ratio = Cm1 / Cm0
+ if cm_ratio > 3.0:
+ cm_ratio = 3.0
+ elif cm_ratio < 0.1:
+ cm_ratio = 0.1
+ else:
+ cm_ratio = 1.0
+
+ c2 = (1 + Cm0 - Dm0) / den * cm_ratio
+ c3 = (1 - Cm0 + Dm0) / den * cm_ratio
c4 = (2 * Cm1) / den
# cmcheck
@@ -304,9 +338,23 @@ def MCTRouting_single(
# Mass balance equation that takes into consideration the lateral flow
q11 = c1 * q01 + c2 * q00 + c3 * q10 + c4 * ql
+ # Diagnostic - now only fires when the CLAMPED ratio is at the ceiling
+ # or other suspicious conditions
+ # if cm_ratio >= 3.0 - 1e-9: # ratio was clamped
+ # print("CLAMPED Cm0=", Cm0, "Cm1=", Cm1, "raw_ratio=", Cm1/max(Cm0, 0.001), "q01=", q01)
+
if q11 < 0: # cmcheck <=0 #tpk
q11 = 0 #tpk
+ # NEW: mass-balance upper bound
+ # if ql > 0:
+ # ql_abs = ql
+ # else:
+ # ql_abs = -ql
+ max_q11 = q01 + q00 + abs(ql) + V00 / dt
+ if q11 > max_q11:
+ q11 = max_q11
+
#### end of for loop
# # cmcheck
@@ -326,6 +374,10 @@ def MCTRouting_single(
else:
V11 = (1 - Dm1) * dt / (2 * Cm1) * q01 + (1 + Dm1) * dt / (2 * Cm1) * q11
# V11 = k1 * (x1 * q01 + (1. - x1) * q11) # MUST be the same as above!
+ # transition = min(1.0, q11 / (q11 + eps)) # Sigmoid-like transition
+ # V11_formula1 = V00 + (q00 + q01 - q10 - q11) * dt / 2
+ # V11_formula2 = (1 - Dm1) * dt / (2 * Cm1) * q01 + (1 + Dm1) * dt / (2 * Cm1) * q11
+ # V11 = transition * V11_formula2 + (1 - transition) * V11_formula1
if V11 < 0 : #tpk
V11 = 0 #tpk
@@ -334,6 +386,14 @@ def MCTRouting_single(
# calc average discharge outflow q1m for MCT channels during routing sub step dt
# Calculate average outflow using water balance for MCT channel grid cell over sub-routing step
q1mm = q0mm + ql + (V00 - V11) / dt
+ # Ensure q1mm is consistent with q11 (instantaneous outflow)
+ # q1mm should be within reasonable bounds of q11
+ # if q11 > 0:
+ # ratio = q1mm / q11
+ # if ratio > 10:
+ # q1mm = q11
+ # elif ratio < 0.1:
+ # q1mm = q11
# cmcheck
# q1m cannot be smaller than eps or it will cause instability
diff --git a/src/lisflood/hydrological_modules/routing.py b/src/lisflood/hydrological_modules/routing.py
index feb2fc7d..886ef50c 100644
--- a/src/lisflood/hydrological_modules/routing.py
+++ b/src/lisflood/hydrological_modules/routing.py
@@ -50,7 +50,8 @@ class routing(HydroModule):
'ChanBottomWMult', 'ChanDepthTMult', 'ChanSMult'],
'SplitRouting': ['CrossSection2AreaInitValue', 'PrevSideflowInitValue', 'CalChanMan2'],
'dynamicWave': ['ChannelsDynamic'],
- 'MCTRouting': ['ChannelsMCT', 'ChanGradMaxMCT', 'PrevCmMCTInitValue', 'PrevDmMCTInitValue', 'CalChanMan3']}
+ 'MCTRouting': ['ChannelsMCT', 'ChanGradMaxMCT', 'PrevCmMCTInitValue', 'PrevDmMCTInitValue', 'CalChanMan3'],
+ 'simulateCalibrationPoints': ['CalibrationPoints']}
module_name = 'Routing'
def __init__(self, routing_variable):
@@ -140,6 +141,7 @@ def initial(self):
self.var.LddToChan = lddrepair(ifthenelse(self.var.IsChannelPcr, 5, self.var.Ldd)) #pcr
self.var.LddToChanNp=compressArray(self.var.LddToChan) #np
# Routing of runoff (incl. groundwater) to the river channel
+ # LDD for routing runoff (incl. groundwater) to the channel
if option['dynamicWave']:
pass
@@ -174,6 +176,45 @@ def initial(self):
self.var.LddKinematic = self.var.LddChan
self.var.LddKinematicNp = compressArray(self.var.LddKinematic) # np
+
+ # At this point, LddChan and LddKinematic do not have sinks at reservoirs/lakes or MCT headwater and MCT confluences
+ # LddMCT does not exist yet
+
+ # ************************************************************
+ # ***** MCT DRAINAGE NETWORK GEOMETRY - LDD *****************
+ # ************************************************************
+
+ # This is done here to be able to add MCT headwater pixels to structures
+ if option['MCTRouting']:
+
+ self.var.IsChannelMCTPcr = boolean(loadmap('ChannelsMCT', pcr=True)) # pcr
+ # load mask of MCT river grid cells
+ self.var.IsChannelMCT = np.bool8(compressArray(self.var.IsChannelMCTPcr)) # bool
+
+ # even if MCT is active, it should be deactivated if there is no MCT cell in the domain
+ if self.var.IsChannelMCT.sum() == 0:
+ warnings.warn(LisfloodWarning('There are no MCT grid cell. MCT routing is deactivated'))
+ option['MCTRouting'] = False
+ # rebuild lists of reported files with MCTRouting = False
+ settings.build_reportedmaps_dicts()
+
+ if option['MCTRouting'] and not option['InitLisflood']:
+
+ self.var.IsChannelMCTPcr = boolean(decompress(self.var.IsChannelMCT)) # pcr
+ # Identify channel pixels where Muskingum-Cunge-Todini is used
+
+ self.var.mctmask = np.bool8(pcr2numpy(self.var.IsChannelMCTPcr,0))
+ # Create a mask with cells using MCT
+
+ self.var.IsChannelKinematicPcr = (self.var.IsChannelPcr == 1) & (self.var.IsChannelMCTPcr == 0) #pcr
+ self.var.IsChannelKinematic = np.bool8(compressArray(self.var.IsChannelKinematicPcr)) #np
+ # Identify channel pixels where Kinematic wave is used instead of MCT
+
+
+ # ************************************************************
+ # ***** MCT DRAINAGE NETWORK GEOMETRY - LDD *****************
+ # ************************************************************
+
self.var.AtLastPoint = boolean(pit(self.var.Ldd)) #pcr
# Assign True to each of the grid cells where there are outlet points
# Function 'pit' assigns a unique number starting from 1 to pit cells (ldd=5) in the Ldd
@@ -248,6 +289,10 @@ def initial(self):
TotalCrossSectionAreaHalfBankFull = BankFullPerc * self.var.TotalCrossSectionAreaBankFull
# set BankFullPerc to 0.5 for half bankfull
+ # Channel volume initialization for MCT cells
+ # TotalCrossSectionAreaHalfBankFull = np.where(self.var.IsChannelKinematic, TotalCrossSectionAreaHalfBankFull, 0.01 * self.var.TotalCrossSectionAreaBankFull)
+ # set initial volume in MCT cells to 1% of bankfull
+
TotalCrossSectionAreaInitValue = loadmap('TotalCrossSectionAreaInitValue')
self.var.TotalCrossSectionArea = np.where(TotalCrossSectionAreaInitValue == -9999, TotalCrossSectionAreaHalfBankFull, TotalCrossSectionAreaInitValue)
# Total cross-sectional area [m2]: if initial value in binding equals -9999 the value at half bankfull is used,
@@ -459,7 +504,7 @@ def initialSecond(self):
# Over bankful discharge starts at QLimit
# lower discharge limit for second line of routing
# set to multiple of average discharge (map from prerun)
- # QSplitMult =2 is around 90 to 95% of Q
+ # QSplitMult = 2 is around 90 to 95% of Q
self.var.QLimit = loadmap_base('AvgDis') * loadmap('QSplitMult')
# Using loadmap_base function as we don't want to cache avgdis in the calibration
@@ -486,7 +531,15 @@ def initialSecond(self):
# Total (virtual) outflow from river channel when second routing line is active (= using increased Manning coeff 2)
self.var.ChanQKin = (self.var.ChanM3Kin * self.var.InvChanLength * self.var.InvChannelAlpha) ** (self.var.InvBeta)
# (Real) outflow from main channel when second line of routing is active (= using riverbed Manning coeff 2)
-
+
+ def initialKinematicWave(self):
+ """ Initialization of the parallel kinematic wave router for Kinematic routing and SplitRouting:
+ main channel-only routing if self.var.ChannelAlpha2 is None; else split-routing(main channel + floodplains).
+ Initialization uses LDD for kinematic routing (LddKinematic)
+ """
+ settings = LisSettings.instance()
+ option = settings.options
+ flags = settings.flags
# ************************************************************
# ***** INITIALISE PARALLEL KINEMATIC WAVE ROUTER ************
@@ -505,7 +558,7 @@ def initialSecond(self):
if option['InitLisflood'] and option['repMBTs']:
# Calculate initial water storage in rivers (no lakes no reservoirs)
# self.var.StorageStepINIT= self.var.ChanM3Kin
- self.var.StorageStepINIT = self.var.ChanM3
+ self.var.StorageStepINIT = self.var.ChanM3.copy()
# Initial water volume in river channels
self.var.DischargeM3StructuresIni = maskinfo.in_zero()
if option['simulateReservoirs']:
@@ -515,7 +568,7 @@ def initialSecond(self):
self.var.StorageStepINIT = np.take(np.bincount(self.var.Catchments, weights=self.var.StorageStepINIT), self.var.Catchments)
if not option['InitLisflood'] and option['repMBTs']:
- self.var.StorageStepINIT = self.var.ChanM3
+ self.var.StorageStepINIT = self.var.ChanM3.copy()
# DisStructure = np.where(self.var.IsUpsOfStructureKinematicC, self.var.ChanQ * self.var.DtRouting, 0)
DisStructure = np.where(self.var.IsUpsOfStructureChanC, self.var.ChanQ * self.var.DtRouting, 0)
if not(option['SplitRouting']):
@@ -554,28 +607,13 @@ def initialMCT(self):
# ***** INITIALISATION FOR MCT ROUTING ************
# ************************************************************
- # even if MCT is active, it should be deactivated if there is no MCT cell in the domain
- if option['MCTRouting']:
- self.var.IsChannelMCTPcr = boolean(loadmap('ChannelsMCT', pcr=True)) #pcr
- self.var.IsChannelMCT = np.bool8(compressArray(self.var.IsChannelMCTPcr)) #bool
- if self.var.IsChannelMCT.sum()==0:
- warnings.warn(LisfloodWarning('There are no MCT grid cell. MCT routing is deactivated'))
- option['MCTRouting'] = False
- # rebuild lists of reported files with MCTRouting = False
- settings.build_reportedmaps_dicts()
if option['MCTRouting'] and not option['InitLisflood']:
maskinfo = MaskInfo.instance()
- self.var.IsChannelMCTPcr = boolean(decompress(self.var.IsChannelMCT)) # pcr
- # Identify channel pixels where Muskingum-Cunge-Todini is used
-
- self.var.mctmask = np.bool8(pcr2numpy(self.var.IsChannelMCTPcr,0))
- # Create a mask with cells using MCT
-
- self.var.IsChannelKinematicPcr = (self.var.IsChannelPcr == 1) & (self.var.IsChannelMCTPcr == 0) #pcr
- self.var.IsChannelKinematic = np.bool8(compressArray(self.var.IsChannelKinematicPcr)) #np
- # Identify channel pixels where Kinematic wave is used instead of MCT
+ # self.var.IsChannelKinematicPcr = (self.var.IsChannelPcr == 1) & (self.var.IsChannelMCTPcr == 0) #pcr
+ # self.var.IsChannelKinematic = np.bool8(compressArray(self.var.IsChannelKinematicPcr)) #np
+ # # Identify channel pixels where Kinematic wave is used instead of MCT
self.var.LddMCT = lddmask(self.var.LddChan, self.var.IsChannelMCTPcr) #pcr
# Ldd for MCT routing
@@ -591,7 +629,6 @@ def initialMCT(self):
self.var.ChanGrad[MCT_slope_mask] = ChanGradMaxMCT
# set max channel slope for MCT pixels
- # cmcheck
# This could become a calibration parameter if we want to use MCT+SplitRouting
self.var.ChanManMCT = (self.var.ChanMan / self.var.CalChanMan) * loadmap('CalChanMan3')
# Mannings coefficient for MCT pixels (same as second line of split routing)
@@ -603,12 +640,29 @@ def initialMCT(self):
self.var.PrevDm0 = np.where(PrevDmMCT == -9999, maskinfo.in_zero(), PrevDmMCT) #np
# Reynolds number (Dm) for MCT at previous time step t0
+ # ************************************************************
+ # ***** CALIBRATION POINTS ********
+ # ************************************************************
+ CalibPoints = maskinfo.in_zero()
+ if option['simulateCalibrationPoints']:
+ CalibPoints = loadmap('CalibrationPoints') # 1D array size all catchment pixels
+ # read location of calibration points
+
+ inAr = np.arange(maskinfo.info.mapC[0], dtype="int32") # np
+ # Assign a number to each non-missing pixel as cell id, by row starting from 0
+ CalibPointsIds = inAr[CalibPoints > 0]
+ # pixel id of calibration points
+
+
# ************************************************************
# ***** INITIALISE MUSKINGUM-CUNGE-TODINI WAVE ROUTER ********
# ************************************************************
mct_ldd = self.compress_mct(compressArray(self.var.LddMCT))
# Compress LddMCT to array with MCT pixels only
+ # mct_CalInflowPoints = self.compress_mct(self.var.CalInflowPoints)
+ # # Compress CalInflowPoints to array with MCT pixels only
+
mapping_mct = self.compress_mct(range(len(self.var.ChanLength)))
# create mapping from global domain pixels index to MCT pixels index
@@ -622,7 +676,8 @@ def initialMCT(self):
self.var.ChanSdXdY, # Riverbed side slope
self.var.DtRouting, # computation time step for routing [s]
self.river_router, # class
- mapping_mct # MCT pixels mapping
+ mapping_mct, # MCT pixels mapping
+ CalibPointsIds, # id of calibrationn points in full LDD
)
@@ -702,7 +757,7 @@ def dynamic(self, NoRoutingExecuted):
else:
self.var.AddedTRUN += np.take(np.bincount(self.var.Catchments, weights=self.var.ToChanM3RunoffDt.copy()),self.var.Catchments)
if option['inflow']:
- self.var.AddedTRUN += np.take(np.bincount(self.var.Catchments, weights=self.var.QInDt),self.var.Catchments)
+ self.var.AddedTRUN += np.take(np.bincount(self.var.Catchments, weights=self.var.QInDt),self.var.Catchments)
if option['openwaterevapo']:
self.var.AddedTRUN -= np.take(np.bincount(self.var.Catchments, weights=self.var.EvaAddM3Dt.copy()),self.var.Catchments)
if option['wateruse']:
@@ -782,9 +837,6 @@ def dynamic(self, NoRoutingExecuted):
# First, Kinematic/Split routing is solved on all pixels (including MCT pixels) then results are updated
# for the MCT pixels.
- # Sideflow contribution to MCT grid cells expressed in [m3/s]
- SideflowChanMCT = np.where(self.var.IsChannelMCT, SideflowChanM3 * self.var.InvDtRouting, 0) #Ql
-
# Grab outflow at the end of the previous routing step t for all pixels) - current state of the MCT pixel
ChanQ_0 = self.var.ChanQ.copy() # Outflow (x+dx) at time t (end of previous routing step) (instant) -> used to calc q00
ChanM3_0 = self.var.ChanM3.copy() # Channel storage at time t (end of previous routing step) (instant) V00
@@ -795,6 +847,10 @@ def dynamic(self, NoRoutingExecuted):
self.var.ChanM3 = ChanM3
self.var.ChanQAvgDt = ChanQAvgDt # -> used to calc q0m
+ # Sideflow contribution to MCT grid cells expressed in [m3/s]
+ SideflowChanMCT = np.where(self.var.IsChannelMCT, SideflowChanM3 * self.var.InvDtRouting, 0) #Ql
+ # SideflowChanMCTM3 = np.where(self.var.IsChannelMCT, SideflowChanM3, 0)
+
# Solve MCT routing and update current state at MCT pixels
self.mct_river_router.routing(
ChanQ_0, # -> used to calc q00
@@ -807,7 +863,32 @@ def dynamic(self, NoRoutingExecuted):
self.var.PrevDm0, # Reynolds number in input: at time t; in output: at time t+dt
self.var.ChanM3 # Channel storage volume. In input: at time t V00; in output: at time t+dt V11
)
+ problematic = np.where(self.var.PrevDm0[self.var.IsChannelMCT] > 2)[0]
+ # Case 1: ChanQAvgDt too small relative to ChanQ (flood arrival)
+ mct_avg_too_small = (
+ self.var.IsChannelMCT &
+ (self.var.ChanQ > 0.1) &
+ (self.var.ChanQAvgDt < 0.1 * self.var.ChanQ)
+ )
+ self.var.ChanQAvgDt = np.where(mct_avg_too_small, self.var.ChanQ, self.var.ChanQAvgDt)
+
+ # Case 2: ChanQAvgDt too large relative to ChanQ (flood recession / cold start)
+ mct_avg_too_large = (
+ self.var.IsChannelMCT &
+ (self.var.ChanQAvgDt > 0.1) &
+ (self.var.ChanQ < 0.1 * self.var.ChanQAvgDt)
+ )
+ self.var.ChanQAvgDt = np.where(mct_avg_too_large, self.var.ChanQ, self.var.ChanQAvgDt)
+ # if 95 <= self.var.currentStep <= 110:
+ # mct_pixels = self.var.IsChannelMCT
+ # cm0 = self.var.PrevCm0[mct_pixels] # this is now Cm1 after routing
+ # # We can't get cm_ratio directly but we can see Cm distribution
+ # print(f"Step {self.var.currentStep}, substep {NoRoutingExecuted}:")
+ # print(f" Cm: min={cm0.min():.2f} mean={cm0.mean():.2f} max={cm0.max():.2f}")
+ # print(f" Dm: min={self.var.PrevDm0[mct_pixels].min():.2f} "
+ # f"mean={self.var.PrevDm0[mct_pixels].mean():.2f} "
+ # f"max={self.var.PrevDm0[mct_pixels].max():.2f}")
##################################################################3
if flags['debug']:
# checking Courant number for potential instability in MCT
@@ -824,11 +905,19 @@ def dynamic(self, NoRoutingExecuted):
bad = too_large | too_small
- if np.any(bad):
- warnings.warn(LisfloodWarning("WARNING! At least one ChanQ is >> or << ChanQAvgDt. Consider increasing DtRouting step or using kinematic routing"))
- # # list 'bad' cells
- # bad_indices = np.where(dismask)[0][bad]
- # print("Bad indices:", bad_indices)
+ # if np.any(bad):
+ # warnings.warn(LisfloodWarning("WARNING! At least one ChanQ is >> or << ChanQAvgDt. Consider increasing DtRouting step or using kinematic routing"))
+ # # # list 'bad' cells
+ # # bad_indices = np.where(dismask)[0][bad]
+ # # print("Bad indices:", bad_indices)
+ # bad_indices = np.where(dismask)[0][bad]
+ # print("Step:", NoRoutingExecuted)
+ # print("Bad indices:", bad_indices)
+ # print("ChanQ at bad:", self.var.ChanQ[bad_indices])
+ # print("ChanQAvgDt at bad:", self.var.ChanQAvgDt[bad_indices])
+ # print("PrevCm0 at bad:", self.var.PrevCm0[bad_indices])
+ # print("PrevDm0 at bad:", self.var.PrevDm0[bad_indices])
+ # print("IsChannelMCT at bad:", self.var.IsChannelMCT[bad_indices])
##################################################################3
else:
diff --git a/tests/data/LF_ETRS89_UseCase/settings/mct_cold.xml b/tests/data/LF_ETRS89_UseCase/settings/mct_cold.xml
index e4b4879f..91caf0ce 100644
--- a/tests/data/LF_ETRS89_UseCase/settings/mct_cold.xml
+++ b/tests/data/LF_ETRS89_UseCase/settings/mct_cold.xml
@@ -47,6 +47,7 @@ You can use builtin path variables in this template and reference to other paths
+
@@ -1681,6 +1682,13 @@ You can use builtin path variables in this template and reference to other paths
+
+
+ location of calibration points
+ OPTIONAL: nominal map with locations of calibration points
+
+
+
OPTIONAL: observed or simulated input hydrographs as
@@ -5476,6 +5484,13 @@ LFBINDING: MORE LOW-LEVEL CONTROL OVER MODEL IN- AND OUTPUT
+
+
+ location of calibration points
+ OPTIONAL: nominal map with locations of calibration points
+
+
+
Observed or simulated input hydrographs as
diff --git a/tests/data/LF_ETRS89_UseCase/settings/mct_warm.xml b/tests/data/LF_ETRS89_UseCase/settings/mct_warm.xml
index 639730fc..90d26839 100644
--- a/tests/data/LF_ETRS89_UseCase/settings/mct_warm.xml
+++ b/tests/data/LF_ETRS89_UseCase/settings/mct_warm.xml
@@ -47,6 +47,7 @@ You can use builtin path variables in this template and reference to other paths
+
@@ -1650,6 +1651,13 @@ You can use builtin path variables in this template and reference to other paths
+
+
+ location of calibration points
+ OPTIONAL: nominal map with locations of calibration points
+
+
+
OPTIONAL: observed or simulated input hydrographs as
@@ -5426,6 +5434,13 @@ LFBINDING: MORE LOW-LEVEL CONTROL OVER MODEL IN- AND OUTPUT
inflow hydrographs
+
+
+
+ location of calibration points
+ OPTIONAL: nominal map with locations of calibration points
+
+
diff --git a/tests/test_mct_warmstart_slow.py b/tests/test_mct_warmstart_slow.py
index 88d26fd5..baa6b3c8 100644
--- a/tests/test_mct_warmstart_slow.py
+++ b/tests/test_mct_warmstart_slow.py
@@ -77,6 +77,42 @@ def test_mct_warmstart_6h(self):
report_steps = '38220..38830'
self.run_warmstart_by_dtsec('mct_all', dt_sec, dt_sec_channel, step_end, step_start, calendar_day_start,report_steps=report_steps)
+ def test_mct_only_warmstart_daily_withcalibpoints(self):
+ calendar_day_start = '02/01/1990 06:00'
+ step_start = '02/01/2016 06:00'
+ step_end = '31/12/2016 06:00'
+ dt_sec = 86400
+ dt_sec_channel = 3600
+ report_steps = '9496..9861'
+ self.run_warmstart_by_dtsec('mct_calibpoints', dt_sec, dt_sec_channel, step_end, step_start, calendar_day_start,report_steps=report_steps)
+
+ def test_mct_warmstart_daily_withcalibpoints(self):
+ calendar_day_start = '02/01/1990 06:00'
+ step_start = '02/01/2016 06:00'
+ step_end = '31/12/2016 06:00'
+ dt_sec = 86400
+ dt_sec_channel = 3600
+ report_steps = '9496..9861'
+ self.run_warmstart_by_dtsec('mct_all_calibpoints', dt_sec, dt_sec_channel, step_end, step_start, calendar_day_start,report_steps=report_steps)
+
+ def test_mct_only_warmstart_6h_withcalibpoints(self):
+ calendar_day_start = '02/01/1990 06:00'
+ step_start = '01/03/2016 06:00'
+ step_end = '31/07/2016 06:00'
+ dt_sec = 21600
+ dt_sec_channel = 3600
+ report_steps = '38220..38830'
+ self.run_warmstart_by_dtsec('mct_calibpoints', dt_sec, dt_sec_channel, step_end, step_start, calendar_day_start,report_steps=report_steps)
+
+ def test_mct_warmstart_6h_withcalibpoints(self):
+ calendar_day_start = '02/01/1990 06:00'
+ step_start = '01/03/2016 06:00'
+ step_end = '31/07/2016 06:00'
+ dt_sec = 21600
+ dt_sec_channel = 3600
+ report_steps = '38220..38830'
+ self.run_warmstart_by_dtsec('mct_all_calibpoints', dt_sec, dt_sec_channel, step_end, step_start, calendar_day_start,report_steps=report_steps)
+
def run_warmstart_by_dtsec(self, mct_case, dt_sec, dt_sec_channel, step_end, step_start, calendar_day_start,report_steps='1..9999'):
mk_path_out(os.path.join(self.case_dir, 'out'))
@@ -101,6 +137,25 @@ def run_warmstart_by_dtsec(self, mct_case, dt_sec, dt_sec_channel, step_end, ste
]
opts_to_unset = ['repMBTs',
'simulateReservoirs']
+ elif mct_case == 'mct_calibpoints':
+ opts_to_set = ['repStateMaps','repDischargeMaps',
+ 'TransLoss',
+ 'simulateCalibrationPoints']
+ opts_to_unset = ['repMBTs', 'simulateReservoirs', 'simulateLakes']
+ elif mct_case == 'mct_all_calibpoints':
+ opts_to_set = ['repStateMaps',
+ 'repDischargeMaps',
+ 'wateruse',
+ 'drainedIrrigation',
+ 'riceIrrigation',
+ 'openwaterevapo',
+ 'simulateLakes',
+ 'TransLoss',
+ 'simulateCalibrationPoints'
+ ]
+ opts_to_unset = ['repMBTs',
+ 'simulateReservoirs']
+
settings_longrun = setoptions(self.settings_files['cold'],
opts_to_set=opts_to_set,