Index: wflow-py/wflow/wflow_sbm.py =================================================================== diff -u -ra1856bef241b96a2ac269325164777e9c34cfa42 -r12ea40dc08628f654753679e0972e87b7bb12f7a --- wflow-py/wflow/wflow_sbm.py (.../wflow_sbm.py) (revision a1856bef241b96a2ac269325164777e9c34cfa42) +++ wflow-py/wflow/wflow_sbm.py (.../wflow_sbm.py) (revision 12ea40dc08628f654753679e0972e87b7bb12f7a) @@ -85,7 +85,8 @@ """ import numpy -#import pcrut + +# import pcrut import sys import os import os.path @@ -103,11 +104,11 @@ updateCols = [] - def usage(*args): sys.stdout = sys.stderr """Way""" - for msg in args: print msg + for msg in args: + print msg print __doc__ sys.exit(0) @@ -126,7 +127,24 @@ return RestPotEvap, FirstZoneDepth, ActEvapSat -def actEvap_unsat_SBM(RootingDepth, WTable, UStoreDepth, zi_layer, UStoreLayerThickness, sumLayer, RestPotEvap, maskLayer, ZeroMap, layerIndex, sumActEvapUStore,c, L, thetaS, thetaR,ust=0): +def actEvap_unsat_SBM( + RootingDepth, + WTable, + UStoreDepth, + zi_layer, + UStoreLayerThickness, + sumLayer, + RestPotEvap, + maskLayer, + ZeroMap, + layerIndex, + sumActEvapUStore, + c, + L, + thetaS, + thetaR, + ust=0, +): """ Actual evaporation function: @@ -146,84 +164,119 @@ - ActEvap, FirstZoneDepth, UStoreDepth ActEvapUStore """ - - #AvailCap is fraction of unsat zone containing roots + # AvailCap is fraction of unsat zone containing roots if ust >= 1: AvailCap = UStoreDepth * 0.99 else: - AvailCap = ifthenelse(layerIndex 0, UStoreDepth / L, 0) vwc = ifthenelse(vwc > 0, vwc, 0.0000001) - head = hb / (((vwc)/(thetaS - thetaR))**(1/par_lambda)) # Note that in the original formula, thetaR is extracted from vwc, but thetaR is not part of the numerical vwc calculation + head = hb / ( + ((vwc) / (thetaS - thetaR)) ** (1 / par_lambda) + ) # Note that in the original formula, thetaR is extracted from vwc, but thetaR is not part of the numerical vwc calculation head = ifthenelse(head <= hb, 1, head) head = cover(head, 0) - + # Transform h to a reduction coefficient value according to Feddes et al. (1978). - alpha = ifthenelse(head <= h1, 0, ifthenelse(head >= h4, 0, ifthenelse(head < h2, (head-h1)/(h2-h1), ifthenelse(head > h3, 1-(head - h3)/(h4 - h3), 1)))) + alpha = ifthenelse( + head <= h1, + 0, + ifthenelse( + head >= h4, + 0, + ifthenelse( + head < h2, + (head - h1) / (h2 - h1), + ifthenelse(head > h3, 1 - (head - h3) / (h4 - h3), 1), + ), + ), + ) - ActEvapUStore = (ifthenelse(layerIndex>zi_layer, ZeroMap,min(MaxExtr, RestPotEvap, UStoreDepth))) * alpha - + ActEvapUStore = ( + ifthenelse( + layerIndex > zi_layer, ZeroMap, min(MaxExtr, RestPotEvap, UStoreDepth) + ) + ) * alpha + UStoreDepth = ifthenelse( + layerIndex > zi_layer, maskLayer, UStoreDepth - ActEvapUStore + ) - UStoreDepth = ifthenelse(layerIndex>zi_layer, maskLayer, UStoreDepth - ActEvapUStore) - - RestPotEvap = RestPotEvap - ActEvapUStore sumActEvapUStore = ActEvapUStore + sumActEvapUStore - - - return UStoreDepth, sumActEvapUStore, RestPotEvap, ActEvapUStore -def soilevap_SBM(CanopyGapFraction,PotTransSoil,SoilWaterCapacity,SatWaterDepth,UStoreLayerDepth,zi,thetaS,thetaR,UStoreLayerThickness): +def soilevap_SBM( + CanopyGapFraction, + PotTransSoil, + SoilWaterCapacity, + SatWaterDepth, + UStoreLayerDepth, + zi, + thetaS, + thetaR, + UStoreLayerThickness, +): # Split between bare soil and vegetation - #potsoilevap = (1.0 - CanopyGapFraction) * PotTransSoil + # potsoilevap = (1.0 - CanopyGapFraction) * PotTransSoil - #PotTrans = CanopyGapFraction * PotTransSoil + # PotTrans = CanopyGapFraction * PotTransSoil SaturationDeficit = SoilWaterCapacity - SatWaterDepth # Linear reduction of soil moisture evaporation based on deficit - soilevap = ifthenelse(len(UStoreLayerThickness)==1, PotTransSoil * min(1.0, SaturationDeficit/SoilWaterCapacity), - PotTransSoil * min(1.0, ifthenelse(zi <= UStoreLayerThickness[0], UStoreLayerDepth[0]/(UStoreLayerThickness[0]*(thetaS-thetaR)), - zi/((zi+1.0)*(thetaS-thetaR))))) + soilevap = ifthenelse( + len(UStoreLayerThickness) == 1, + PotTransSoil * min(1.0, SaturationDeficit / SoilWaterCapacity), + PotTransSoil + * min( + 1.0, + ifthenelse( + zi <= UStoreLayerThickness[0], + UStoreLayerDepth[0] / (UStoreLayerThickness[0] * (thetaS - thetaR)), + zi / ((zi + 1.0) * (thetaS - thetaR)), + ), + ), + ) return soilevap def sum_UstoreLayerDepth(UStoreLayerThickness, ZeroMap, UStoreLayerDepth): sum_UstoreLayerDepth = ZeroMap for n in arange(0, len(UStoreLayerThickness)): - sum_UstoreLayerDepth = sum_UstoreLayerDepth + cover(UStoreLayerDepth[n],ZeroMap) + sum_UstoreLayerDepth = sum_UstoreLayerDepth + cover( + UStoreLayerDepth[n], ZeroMap + ) return sum_UstoreLayerDepth - def SnowPackHBV(Snow, SnowWater, Precipitation, Temperature, TTI, TT, TTM, Cfmax, WHC): """ HBV Type snowpack modelling using a Temperature degree factor. All correction @@ -245,32 +298,45 @@ CFR = 0.05000 # refreeing efficiency constant in refreezing of freewater in snow SFCF = 1.0 # correction factor for snowfall - RainFrac = ifthenelse(1.0 * TTI == 0.0, ifthenelse(Temperature <= TT, scalar(0.0), scalar(1.0)), - min((Temperature - (TT - TTI / 2)) / TTI, scalar(1.0))); - RainFrac = max(RainFrac, scalar(0.0)) #fraction of precipitation which falls as rain - SnowFrac = 1 - RainFrac #fraction of precipitation which falls as snow - Precipitation = SFCF * SnowFrac * Precipitation + RFCF * RainFrac * Precipitation # different correction for rainfall and snowfall + RainFrac = ifthenelse( + 1.0 * TTI == 0.0, + ifthenelse(Temperature <= TT, scalar(0.0), scalar(1.0)), + min((Temperature - (TT - TTI / 2)) / TTI, scalar(1.0)), + ) + RainFrac = max( + RainFrac, scalar(0.0) + ) # fraction of precipitation which falls as rain + SnowFrac = 1 - RainFrac # fraction of precipitation which falls as snow + Precipitation = ( + SFCF * SnowFrac * Precipitation + RFCF * RainFrac * Precipitation + ) # different correction for rainfall and snowfall - SnowFall = SnowFrac * Precipitation #snowfall depth - RainFall = RainFrac * Precipitation #rainfall depth - PotSnowMelt = ifthenelse(Temperature > TTM, Cfmax * (Temperature - TTM), - scalar(0.0)) #Potential snow melt, based on temperature - PotRefreezing = ifthenelse(Temperature < TTM, Cfmax * CFR * (TTM - Temperature), - 0.0) #Potential refreezing, based on temperature - Refreezing = ifthenelse(Temperature < TTM, min(PotRefreezing, SnowWater), 0.0) #actual refreezing + SnowFall = SnowFrac * Precipitation # snowfall depth + RainFall = RainFrac * Precipitation # rainfall depth + PotSnowMelt = ifthenelse( + Temperature > TTM, Cfmax * (Temperature - TTM), scalar(0.0) + ) # Potential snow melt, based on temperature + PotRefreezing = ifthenelse( + Temperature < TTM, Cfmax * CFR * (TTM - Temperature), 0.0 + ) # Potential refreezing, based on temperature + Refreezing = ifthenelse( + Temperature < TTM, min(PotRefreezing, SnowWater), 0.0 + ) # actual refreezing # No landuse correction here - SnowMelt = min(PotSnowMelt, Snow) #actual snow melt - Snow = Snow + SnowFall + Refreezing - SnowMelt #dry snow content - SnowWater = SnowWater - Refreezing #free water content in snow + SnowMelt = min(PotSnowMelt, Snow) # actual snow melt + Snow = Snow + SnowFall + Refreezing - SnowMelt # dry snow content + SnowWater = SnowWater - Refreezing # free water content in snow MaxSnowWater = Snow * WHC # Max water in the snow - SnowWater = SnowWater + SnowMelt + RainFall # Add all water and potentially supersaturate the snowpack + SnowWater = ( + SnowWater + SnowMelt + RainFall + ) # Add all water and potentially supersaturate the snowpack RainFall = max(SnowWater - MaxSnowWater, 0.0) # rain + surpluss snowwater SnowWater = SnowWater - RainFall - return Snow, SnowWater, SnowMelt, RainFall,SnowFall + return Snow, SnowWater, SnowMelt, RainFall, SnowFall -def GlacierMelt(GlacierStore, Snow, Temperature, TT, Cfmax): +def GlacierMelt(GlacierStore, Snow, Temperature, TT, Cfmax): """ Glacier modelling using a Temperature degree factor. Melting only occurs if the snow cover > 10 mm @@ -283,17 +349,18 @@ :returns: GlacierStore,GlacierMelt, """ + PotMelt = ifthenelse( + Temperature > TT, Cfmax * (Temperature - TT), scalar(0.0) + ) # Potential snow melt, based on temperature - PotMelt = ifthenelse(Temperature > TT, Cfmax * (Temperature - TT), - scalar(0.0)) # Potential snow melt, based on temperature - - GlacierMelt = ifthenelse(Snow > 10.0,min(PotMelt, GlacierStore),cover(0.0)) # actual Glacier melt + GlacierMelt = ifthenelse( + Snow > 10.0, min(PotMelt, GlacierStore), cover(0.0) + ) # actual Glacier melt GlacierStore = GlacierStore - GlacierMelt # dry snow content + return GlacierStore, GlacierMelt - return GlacierStore, GlacierMelt - class WflowModel(DynamicModel): """ .. versionchanged:: 0.91 @@ -309,14 +376,14 @@ def __init__(self, cloneMap, Dir, RunDir, configfile): DynamicModel.__init__(self) - self.UStoreLayerDepth= [] + self.UStoreLayerDepth = [] self.caseName = os.path.abspath(Dir) - self.clonemappath = os.path.join(os.path.abspath(Dir),"staticmaps",cloneMap) + self.clonemappath = os.path.join(os.path.abspath(Dir), "staticmaps", cloneMap) setclone(self.clonemappath) self.runId = RunDir self.Dir = os.path.abspath(Dir) self.configfile = configfile - self.SaveDir = os.path.join(self.Dir,self.runId) + self.SaveDir = os.path.join(self.Dir, self.runId) def irrigationdemand(self, pottrans, acttrans, irareas): """ @@ -330,10 +397,10 @@ :return: demand """ - Et_diff = areaaverage(pottrans-acttrans, nominal(irareas)) + Et_diff = areaaverage(pottrans - acttrans, nominal(irareas)) # Now determine demand in m^3/s for each area - sqmarea = areatotal(self.reallength * self.reallength,nominal(irareas)) - m3sec = Et_diff * sqmarea/1000.0/self.timestepsecs + sqmarea = areatotal(self.reallength * self.reallength, nominal(irareas)) + m3sec = Et_diff * sqmarea / 1000.0 / self.timestepsecs return Et_diff, m3sec @@ -369,22 +436,28 @@ :var self.ReservoirVolume: Volume of each reservoir [m^3] :var self.GlacierStore: Thickness of the Glacier in a gridcell [mm] """ - states = ['SurfaceRunoff', 'WaterLevel', - 'SatWaterDepth','Snow', - 'TSoil','UStoreLayerDepth','SnowWater', - 'CanopyStorage'] - if hasattr(self, 'GlacierFrac'): - states.append('GlacierStore') - - if hasattr(self,'ReserVoirSimpleLocs'): - states.append('ReservoirVolume') + states = [ + "SurfaceRunoff", + "WaterLevel", + "SatWaterDepth", + "Snow", + "TSoil", + "UStoreLayerDepth", + "SnowWater", + "CanopyStorage", + ] + if hasattr(self, "GlacierFrac"): + states.append("GlacierStore") - if hasattr(self,'ReserVoirComplexLocs'): - states.append('ReservoirWaterLevel') + if hasattr(self, "ReserVoirSimpleLocs"): + states.append("ReservoirVolume") - if hasattr(self,'nrpaddyirri'): + if hasattr(self, "ReserVoirComplexLocs"): + states.append("ReservoirWaterLevel") + + if hasattr(self, "nrpaddyirri"): if self.nrpaddyirri > 0: - states.append('PondingDepth') + states.append("PondingDepth") return states def supplyCurrentTime(self): @@ -396,14 +469,12 @@ def suspend(self): self.logger.info("Saving initial conditions...") - self.wf_suspend(os.path.join(self.SaveDir,"outstate")) + self.wf_suspend(os.path.join(self.SaveDir, "outstate")) if self.OverWriteInit: self.logger.info("Saving initial conditions over start conditions...") self.wf_suspend(self.SaveDir + "/instate/") - - def parameters(self): """ Define all model parameters here that the framework should handle for the model @@ -413,49 +484,138 @@ """ modelparameters = [] - #Static model parameters e.g. - #modelparameters.append(self.ParamType(name="RunoffGeneratingGWPerc",stack="intbl/RunoffGeneratingGWPerc.tbl",type="static",default=0.1)) + # Static model parameters e.g. + # modelparameters.append(self.ParamType(name="RunoffGeneratingGWPerc",stack="intbl/RunoffGeneratingGWPerc.tbl",type="static",default=0.1)) # 3: Input time series ################################################### - self.P_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Precipitation", - "/inmaps/P") # timeseries for rainfall - self.PET_mapstack = self.Dir + configget(self.config, "inputmapstacks", "EvapoTranspiration", - "/inmaps/PET") # timeseries for rainfall"/inmaps/PET" # potential evapotranspiration - self.TEMP_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Temperature", - "/inmaps/TEMP") # timeseries for rainfall "/inmaps/TEMP" # global radiation - self.Inflow_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Inflow", - "/inmaps/IF") # timeseries for rainfall "/inmaps/IF" # in/outflow locations (abstractions) + self.P_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "Precipitation", "/inmaps/P" + ) # timeseries for rainfall + self.PET_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "EvapoTranspiration", "/inmaps/PET" + ) # timeseries for rainfall"/inmaps/PET" # potential evapotranspiration + self.TEMP_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "Temperature", "/inmaps/TEMP" + ) # timeseries for rainfall "/inmaps/TEMP" # global radiation + self.Inflow_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "Inflow", "/inmaps/IF" + ) # timeseries for rainfall "/inmaps/IF" # in/outflow locations (abstractions) # Meteo and other forcing - modelparameters.append(self.ParamType(name="Precipitation",stack=self.P_mapstack,type="timeseries",default=0.0,verbose=True,lookupmaps=[])) - modelparameters.append(self.ParamType(name="PotenEvap",stack=self.PET_mapstack,type="timeseries",default=0.0,verbose=True,lookupmaps=[])) - modelparameters.append(self.ParamType(name="Temperature",stack=self.TEMP_mapstack,type="timeseries",default=10.0,verbose=True,lookupmaps=[])) - modelparameters.append(self.ParamType(name="Inflow",stack=self.Inflow_mapstack,type="timeseries",default=0.0,verbose=False,lookupmaps=[])) - + modelparameters.append( + self.ParamType( + name="Precipitation", + stack=self.P_mapstack, + type="timeseries", + default=0.0, + verbose=True, + lookupmaps=[], + ) + ) + modelparameters.append( + self.ParamType( + name="PotenEvap", + stack=self.PET_mapstack, + type="timeseries", + default=0.0, + verbose=True, + lookupmaps=[], + ) + ) + modelparameters.append( + self.ParamType( + name="Temperature", + stack=self.TEMP_mapstack, + type="timeseries", + default=10.0, + verbose=True, + lookupmaps=[], + ) + ) + modelparameters.append( + self.ParamType( + name="Inflow", + stack=self.Inflow_mapstack, + type="timeseries", + default=0.0, + verbose=False, + lookupmaps=[], + ) + ) - modelparameters.append(self.ParamType(name="IrrigationAreas", stack='staticmaps/wflow_irrigationareas.map', - type="staticmap", default=0.0, verbose=False, lookupmaps=[])) - modelparameters.append(self.ParamType(name="IrrigationSurfaceIntakes", stack='staticmaps/wflow_irrisurfaceintake.map', - type="staticmap", default=0.0, verbose=False, lookupmaps=[])) - modelparameters.append(self.ParamType(name="IrrigationPaddyAreas", stack='staticmaps/wflow_irrigationpaddyareas.map', - type="staticmap", default=0.0, verbose=False, lookupmaps=[])) modelparameters.append( - self.ParamType(name="IrrigationSurfaceReturn", stack='staticmaps/wflow_irrisurfacereturns.map', - type="staticmap", default=0.0, verbose=False, lookupmaps=[])) + self.ParamType( + name="IrrigationAreas", + stack="staticmaps/wflow_irrigationareas.map", + type="staticmap", + default=0.0, + verbose=False, + lookupmaps=[], + ) + ) + modelparameters.append( + self.ParamType( + name="IrrigationSurfaceIntakes", + stack="staticmaps/wflow_irrisurfaceintake.map", + type="staticmap", + default=0.0, + verbose=False, + lookupmaps=[], + ) + ) + modelparameters.append( + self.ParamType( + name="IrrigationPaddyAreas", + stack="staticmaps/wflow_irrigationpaddyareas.map", + type="staticmap", + default=0.0, + verbose=False, + lookupmaps=[], + ) + ) + modelparameters.append( + self.ParamType( + name="IrrigationSurfaceReturn", + stack="staticmaps/wflow_irrisurfacereturns.map", + type="staticmap", + default=0.0, + verbose=False, + lookupmaps=[], + ) + ) + modelparameters.append( + self.ParamType( + name="h_max", + stack="staticmaps/wflow_hmax.map", + type="staticmap", + default=0.0, + verbose=False, + lookupmaps=[], + ) + ) + modelparameters.append( + self.ParamType( + name="h_min", + stack="staticmaps/wflow_hmin.map", + type="staticmap", + default=0.0, + verbose=False, + lookupmaps=[], + ) + ) + modelparameters.append( + self.ParamType( + name="h_p", + stack="staticmaps/wflow_hp.map", + type="staticmap", + default=0.0, + verbose=False, + lookupmaps=[], + ) + ) - - modelparameters.append(self.ParamType(name="h_max", stack='staticmaps/wflow_hmax.map', - type="staticmap", default=0.0, verbose=False, lookupmaps=[])) - modelparameters.append(self.ParamType(name="h_min", stack='staticmaps/wflow_hmin.map', - type="staticmap", default=0.0, verbose=False, lookupmaps=[])) - modelparameters.append(self.ParamType(name="h_p", stack='staticmaps/wflow_hp.map', - type="staticmap", default=0.0, verbose=False, lookupmaps=[])) - - return modelparameters - - def initial(self): """ Initial part of the model, executed only once. Reads all static data from disk @@ -509,17 +669,13 @@ global multpars global updateCols - self.thestep = scalar(0) self.basetimestep = 86400 self.SSSF = False setglobaloption("unittrue") - - self.logger.info("running for " + str(self.nrTimeSteps()) + " timesteps") - # Set and get defaults from ConfigFile here ################################### self.Tslice = int(configget(self.config, "model", "Tslice", "1")) @@ -528,19 +684,24 @@ self.updating = int(configget(self.config, "model", "updating", "0")) self.updateFile = configget(self.config, "model", "updateFile", "no_set") self.LateralMethod = int(configget(self.config, "model", "lateralmethod", "1")) - self.TransferMethod = int(configget(self.config, "model", "transfermethod", "1")) + self.TransferMethod = int( + configget(self.config, "model", "transfermethod", "1") + ) self.maxitsupply = int(configget(self.config, "model", "maxitsupply", "5")) self.UST = int(configget(self.config, "model", "Whole_UST_Avail", "0")) self.NRiverMethod = int(configget(self.config, "model", "nrivermethod", "1")) - if self.LateralMethod == 1: - self.logger.info("Applying the original topog_sbm lateral transfer formulation") + self.logger.info( + "Applying the original topog_sbm lateral transfer formulation" + ) elif self.LateralMethod == 2: self.logger.warn("Using alternate wflow lateral transfer formulation") if self.TransferMethod == 1: - self.logger.info("Applying the original topog_sbm vertical transfer formulation") + self.logger.info( + "Applying the original topog_sbm vertical transfer formulation" + ) elif self.TransferMethod == 2: self.logger.warn("Using alternate wflow vertical transfer formulation") @@ -550,78 +711,141 @@ self.modelSnow = int(configget(self.config, "model", "ModelSnow", "1")) sizeinmetres = int(configget(self.config, "layout", "sizeinmetres", "0")) alf = float(configget(self.config, "model", "Alpha", "60")) - #TODO: make this into a list for all gauges or a map + # TODO: make this into a list for all gauges or a map Qmax = float(configget(self.config, "model", "AnnualDischarge", "300")) self.UpdMaxDist = float(configget(self.config, "model", "UpdMaxDist", "100")) self.MaxUpdMult = float(configget(self.config, "model", "MaxUpdMult", "1.3")) self.MinUpdMult = float(configget(self.config, "model", "MinUpdMult", "0.7")) self.UpFrac = float(configget(self.config, "model", "UpFrac", "0.8")) - #self.ExternalQbase=int(configget(self.config,'model','ExternalQbase','0')) - self.waterdem = int(configget(self.config, 'model', 'waterdem', '0')) - WIMaxScale = float(configget(self.config, 'model', 'WIMaxScale', '0.8')) - self.reInfilt = int(configget(self.config, 'model', 'reInfilt', '0')) - self.MassWasting = int(configget(self.config,"model","MassWasting","0")) + # self.ExternalQbase=int(configget(self.config,'model','ExternalQbase','0')) + self.waterdem = int(configget(self.config, "model", "waterdem", "0")) + WIMaxScale = float(configget(self.config, "model", "WIMaxScale", "0.8")) + self.reInfilt = int(configget(self.config, "model", "reInfilt", "0")) + self.MassWasting = int(configget(self.config, "model", "MassWasting", "0")) - self.nrLayers = int(configget(self.config,"model","nrLayers",'1')) + self.nrLayers = int(configget(self.config, "model", "nrLayers", "1")) - - - # static maps to use (normally default) - wflow_subcatch = configget(self.config, "model", "wflow_subcatch", "staticmaps/wflow_subcatch.map") - wflow_dem = configget(self.config, "model", "wflow_dem", "staticmaps/wflow_dem.map") - wflow_ldd = configget(self.config, "model", "wflow_ldd", "staticmaps/wflow_ldd.map") - wflow_river = configget(self.config, "model", "wflow_river", "staticmaps/wflow_river.map") - wflow_riverlength = configget(self.config, "model", "wflow_riverlength", "staticmaps/wflow_riverlength.map") - wflow_riverlength_fact = configget(self.config, "model", "wflow_riverlength_fact", - "staticmaps/wflow_riverlength_fact.map") - wflow_landuse = configget(self.config, "model", "wflow_landuse", "staticmaps/wflow_landuse.map") - wflow_soil = configget(self.config, "model", "wflow_soil", "staticmaps/wflow_soil.map") - wflow_gauges = configget(self.config, "model", "wflow_gauges", "staticmaps/wflow_gauges.map") - wflow_inflow = configget(self.config, "model", "wflow_inflow", "staticmaps/wflow_inflow.map") - wflow_riverwidth = configget(self.config, "model", "wflow_riverwidth", "staticmaps/wflow_riverwidth.map") - wflow_streamorder = configget(self.config, "model", "wflow_streamorder", "staticmaps/wflow_streamorder.map") + wflow_subcatch = configget( + self.config, "model", "wflow_subcatch", "staticmaps/wflow_subcatch.map" + ) + wflow_dem = configget( + self.config, "model", "wflow_dem", "staticmaps/wflow_dem.map" + ) + wflow_ldd = configget( + self.config, "model", "wflow_ldd", "staticmaps/wflow_ldd.map" + ) + wflow_river = configget( + self.config, "model", "wflow_river", "staticmaps/wflow_river.map" + ) + wflow_riverlength = configget( + self.config, + "model", + "wflow_riverlength", + "staticmaps/wflow_riverlength.map", + ) + wflow_riverlength_fact = configget( + self.config, + "model", + "wflow_riverlength_fact", + "staticmaps/wflow_riverlength_fact.map", + ) + wflow_landuse = configget( + self.config, "model", "wflow_landuse", "staticmaps/wflow_landuse.map" + ) + wflow_soil = configget( + self.config, "model", "wflow_soil", "staticmaps/wflow_soil.map" + ) + wflow_gauges = configget( + self.config, "model", "wflow_gauges", "staticmaps/wflow_gauges.map" + ) + wflow_inflow = configget( + self.config, "model", "wflow_inflow", "staticmaps/wflow_inflow.map" + ) + wflow_riverwidth = configget( + self.config, "model", "wflow_riverwidth", "staticmaps/wflow_riverwidth.map" + ) + wflow_streamorder = configget( + self.config, + "model", + "wflow_streamorder", + "staticmaps/wflow_streamorder.map", + ) - # 2: Input base maps ######################################################## - subcatch = ordinal(self.wf_readmap(os.path.join(self.Dir,wflow_subcatch),0.0,fail=True)) # Determines the area of calculations (all cells > 0) + subcatch = ordinal( + self.wf_readmap(os.path.join(self.Dir, wflow_subcatch), 0.0, fail=True) + ) # Determines the area of calculations (all cells > 0) subcatch = ifthen(subcatch > 0, subcatch) - self.Altitude = self.wf_readmap(os.path.join(self.Dir,wflow_dem),0.0,fail=True) # * scalar(defined(subcatch)) # DEM - self.TopoLdd = ldd(self.wf_readmap(os.path.join(self.Dir,wflow_ldd),0.0,fail=True)) # Local - self.TopoId = ordinal(self.wf_readmap(os.path.join(self.Dir,wflow_subcatch),0.0,fail=True)) # area map - self.River = cover(boolean(self.wf_readmap(os.path.join(self.Dir,wflow_river),0.0,fail=True)), 0) + self.Altitude = self.wf_readmap( + os.path.join(self.Dir, wflow_dem), 0.0, fail=True + ) # * scalar(defined(subcatch)) # DEM + self.TopoLdd = ldd( + self.wf_readmap(os.path.join(self.Dir, wflow_ldd), 0.0, fail=True) + ) # Local + self.TopoId = ordinal( + self.wf_readmap(os.path.join(self.Dir, wflow_subcatch), 0.0, fail=True) + ) # area map + self.River = cover( + boolean( + self.wf_readmap(os.path.join(self.Dir, wflow_river), 0.0, fail=True) + ), + 0, + ) - self.RiverLength = cover(self.wf_readmap(os.path.join(self.Dir,wflow_riverlength), 0.0), 0.0) + self.RiverLength = cover( + self.wf_readmap(os.path.join(self.Dir, wflow_riverlength), 0.0), 0.0 + ) # Factor to multiply riverlength with (defaults to 1.0) - self.RiverLengthFac = self.wf_readmap(os.path.join(self.Dir,wflow_riverlength_fact), 1.0) + self.RiverLengthFac = self.wf_readmap( + os.path.join(self.Dir, wflow_riverlength_fact), 1.0 + ) # read landuse and soilmap and make sure there are no missing points related to the # subcatchment map. Currently sets the lu and soil type type to 1 - self.LandUse = ordinal(self.wf_readmap(os.path.join(self.Dir,wflow_landuse),0.0,fail=True)) + self.LandUse = ordinal( + self.wf_readmap(os.path.join(self.Dir, wflow_landuse), 0.0, fail=True) + ) self.LandUse = cover(self.LandUse, ordinal(subcatch > 0)) - self.Soil = ordinal(self.wf_readmap(os.path.join(self.Dir,wflow_soil),0.0,fail=True)) + self.Soil = ordinal( + self.wf_readmap(os.path.join(self.Dir, wflow_soil), 0.0, fail=True) + ) self.Soil = cover(self.Soil, ordinal(subcatch > 0)) - self.OutputLoc = ordinal(self.wf_readmap(os.path.join(self.Dir,wflow_gauges),0.0,fail=True) ) # location of output gauge(s) - self.InflowLoc = ordinal(self.wf_readmap(os.path.join(self.Dir,wflow_inflow), 0.0) ) # location abstractions/inflows. - self.RiverWidth = self.wf_readmap(os.path.join(self.Dir,wflow_riverwidth), 0.0) + self.OutputLoc = ordinal( + self.wf_readmap(os.path.join(self.Dir, wflow_gauges), 0.0, fail=True) + ) # location of output gauge(s) + self.InflowLoc = ordinal( + self.wf_readmap(os.path.join(self.Dir, wflow_inflow), 0.0) + ) # location abstractions/inflows. + self.RiverWidth = self.wf_readmap(os.path.join(self.Dir, wflow_riverwidth), 0.0) # Experimental - self.RunoffGenSigmaFunction = int(configget(self.config, 'model', 'RunoffGenSigmaFunction', '0')) - self.SubCatchFlowOnly = int(configget(self.config, 'model', 'SubCatchFlowOnly', '0')) - self.OutputId = ordinal(self.wf_readmap(os.path.join(self.Dir,wflow_subcatch),0.0,fail=True)) # location of subcatchment + self.RunoffGenSigmaFunction = int( + configget(self.config, "model", "RunoffGenSigmaFunction", "0") + ) + self.SubCatchFlowOnly = int( + configget(self.config, "model", "SubCatchFlowOnly", "0") + ) + self.OutputId = ordinal( + self.wf_readmap(os.path.join(self.Dir, wflow_subcatch), 0.0, fail=True) + ) # location of subcatchment # Temperature correction poer cell to add - - self.TempCor = self.wf_readmap( - self.Dir + configget(self.config, "model", "TemperatureCorrectionMap", "staticmaps/wflow_tempcor.map"), 0.0) + self.Dir + + configget( + self.config, + "model", + "TemperatureCorrectionMap", + "staticmaps/wflow_tempcor.map", + ), + 0.0, + ) - self.ZeroMap = 0.0 * scalar(subcatch) #map with only zero's + self.ZeroMap = 0.0 * scalar(subcatch) # map with only zero's - - # Set static initial values here ######################################### self.pi = 3.1416 self.e = 2.7183 @@ -633,85 +857,210 @@ self.logger.info("Linking parameters to landuse, catchment and soil...") self.wf_updateparameters() - self.RunoffGeneratingGWPerc = self.readtblDefault(self.Dir + "/" + self.intbl + "/RunoffGeneratingGWPerc.tbl", - self.LandUse, subcatch, self.Soil, - 0.1) + self.RunoffGeneratingGWPerc = self.readtblDefault( + self.Dir + "/" + self.intbl + "/RunoffGeneratingGWPerc.tbl", + self.LandUse, + subcatch, + self.Soil, + 0.1, + ) - if hasattr(self,"LAI"): + if hasattr(self, "LAI"): # Sl must also be defined - if not hasattr(self,"Sl"): - logging.error("Sl (specific leaf storage) not defined! Needed becausee LAI is defined.") + if not hasattr(self, "Sl"): + logging.error( + "Sl (specific leaf storage) not defined! Needed becausee LAI is defined." + ) logging.error("Please add it to the modelparameters section. e.g.:") - logging.error("Sl=inmaps/clim/LCtoSpecificLeafStorage.tbl,tbl,0.5,1,inmaps/clim/LC.map") - if not hasattr(self,"Kext"): - logging.error("Kext (canopy extinction coefficient) not defined! Needed becausee LAI is defined.") + logging.error( + "Sl=inmaps/clim/LCtoSpecificLeafStorage.tbl,tbl,0.5,1,inmaps/clim/LC.map" + ) + if not hasattr(self, "Kext"): + logging.error( + "Kext (canopy extinction coefficient) not defined! Needed becausee LAI is defined." + ) logging.error("Please add it to the modelparameters section. e.g.:") - logging.error("Kext=inmaps/clim/LCtoExtinctionCoefficient.tbl,tbl,0.5,1,inmaps/clim/LC.map") - if not hasattr(self,"Swood"): - logging.error("Swood wood (branches, trunks) canopy storage not defined! Needed becausee LAI is defined.") + logging.error( + "Kext=inmaps/clim/LCtoExtinctionCoefficient.tbl,tbl,0.5,1,inmaps/clim/LC.map" + ) + if not hasattr(self, "Swood"): + logging.error( + "Swood wood (branches, trunks) canopy storage not defined! Needed becausee LAI is defined." + ) logging.error("Please add it to the modelparameters section. e.g.:") - logging.error("Swood=inmaps/clim/LCtoBranchTrunkStorage.tbl,tbl,0.5,1,inmaps/clim/LC.map") + logging.error( + "Swood=inmaps/clim/LCtoBranchTrunkStorage.tbl,tbl,0.5,1,inmaps/clim/LC.map" + ) self.Cmax = self.Sl * self.LAI + self.Swood self.CanopyGapFraction = exp(-self.Kext * self.LAI) # TODO: Add MAXLAI and CWf lookup else: - self.Cmax = self.readtblDefault(self.Dir + "/" + self.intbl + "/MaxCanopyStorage.tbl", self.LandUse, subcatch, - self.Soil, 1.0) - self.CanopyGapFraction = self.readtblDefault(self.Dir + "/" + self.intbl + "/CanopyGapFraction.tbl", - self.LandUse, subcatch, self.Soil, 0.1) - self.EoverR = self.readtblDefault(self.Dir + "/" + self.intbl + "/EoverR.tbl", self.LandUse, subcatch, - self.Soil, 0.1) + self.Cmax = self.readtblDefault( + self.Dir + "/" + self.intbl + "/MaxCanopyStorage.tbl", + self.LandUse, + subcatch, + self.Soil, + 1.0, + ) + self.CanopyGapFraction = self.readtblDefault( + self.Dir + "/" + self.intbl + "/CanopyGapFraction.tbl", + self.LandUse, + subcatch, + self.Soil, + 0.1, + ) + self.EoverR = self.readtblDefault( + self.Dir + "/" + self.intbl + "/EoverR.tbl", + self.LandUse, + subcatch, + self.Soil, + 0.1, + ) - if not hasattr(self,'DemandReturnFlowFraction'): + if not hasattr(self, "DemandReturnFlowFraction"): self.DemandReturnFlowFraction = self.ZeroMap - self.RootingDepth = self.readtblDefault(self.Dir + "/" + self.intbl + "/RootingDepth.tbl", self.LandUse, - subcatch, self.Soil, 750.0) #rooting depth + self.RootingDepth = self.readtblDefault( + self.Dir + "/" + self.intbl + "/RootingDepth.tbl", + self.LandUse, + subcatch, + self.Soil, + 750.0, + ) # rooting depth #: rootdistpar determien how roots are linked to water table. - self.rootdistpar = self.readtblDefault(self.Dir + "/" + self.intbl + "/rootdistpar.tbl", self.LandUse, subcatch, - self.Soil, -8000) #rrootdistpar + self.rootdistpar = self.readtblDefault( + self.Dir + "/" + self.intbl + "/rootdistpar.tbl", + self.LandUse, + subcatch, + self.Soil, + -8000, + ) # rrootdistpar # Soil parameters # infiltration capacity if the soil [mm/day] - self.InfiltCapSoil = self.readtblDefault(self.Dir + "/" + self.intbl + "/InfiltCapSoil.tbl", self.LandUse, - subcatch, self.Soil, 100.0) * self.timestepsecs / self.basetimestep - self.CapScale = self.readtblDefault(self.Dir + "/" + self.intbl + "/CapScale.tbl", self.LandUse, subcatch, - self.Soil, 100.0) # + self.InfiltCapSoil = ( + self.readtblDefault( + self.Dir + "/" + self.intbl + "/InfiltCapSoil.tbl", + self.LandUse, + subcatch, + self.Soil, + 100.0, + ) + * self.timestepsecs + / self.basetimestep + ) + self.CapScale = self.readtblDefault( + self.Dir + "/" + self.intbl + "/CapScale.tbl", + self.LandUse, + subcatch, + self.Soil, + 100.0, + ) # # infiltration capacity of the compacted - self.InfiltCapPath = self.readtblDefault(self.Dir + "/" + self.intbl + "/InfiltCapPath.tbl", self.LandUse, - subcatch, self.Soil, 10.0) * self.timestepsecs / self.basetimestep - self.MaxLeakage = self.readtblDefault(self.Dir + "/" + self.intbl + "/MaxLeakage.tbl", self.LandUse, subcatch, - self.Soil, 0.0) * self.timestepsecs / self.basetimestep - self.MaxPercolation = self.readtblDefault(self.Dir + "/" + self.intbl + "/MaxPercolation.tbl", self.LandUse, subcatch, - self.Soil, 0.0) * self.timestepsecs / self.basetimestep + self.InfiltCapPath = ( + self.readtblDefault( + self.Dir + "/" + self.intbl + "/InfiltCapPath.tbl", + self.LandUse, + subcatch, + self.Soil, + 10.0, + ) + * self.timestepsecs + / self.basetimestep + ) + self.MaxLeakage = ( + self.readtblDefault( + self.Dir + "/" + self.intbl + "/MaxLeakage.tbl", + self.LandUse, + subcatch, + self.Soil, + 0.0, + ) + * self.timestepsecs + / self.basetimestep + ) + self.MaxPercolation = ( + self.readtblDefault( + self.Dir + "/" + self.intbl + "/MaxPercolation.tbl", + self.LandUse, + subcatch, + self.Soil, + 0.0, + ) + * self.timestepsecs + / self.basetimestep + ) - # areas (paths) in [mm/day] # Fraction area with compacted soil (Paths etc.) - self.PathFrac = self.readtblDefault(self.Dir + "/" + self.intbl + "/PathFrac.tbl", self.LandUse, subcatch, - self.Soil, 0.01) + self.PathFrac = self.readtblDefault( + self.Dir + "/" + self.intbl + "/PathFrac.tbl", + self.LandUse, + subcatch, + self.Soil, + 0.01, + ) # thickness of the soil - self.SoilThickness = self.readtblDefault(self.Dir + "/" + self.intbl + "/SoilThickness.tbl", - self.LandUse, subcatch, self.Soil, 2000.0) - self.thetaR = self.readtblDefault(self.Dir + "/" + self.intbl + "/thetaR.tbl", self.LandUse, subcatch, - self.Soil, 0.01) - self.thetaS = self.readtblDefault(self.Dir + "/" + self.intbl + "/thetaS.tbl", self.LandUse, subcatch, - self.Soil, 0.6) + self.SoilThickness = self.readtblDefault( + self.Dir + "/" + self.intbl + "/SoilThickness.tbl", + self.LandUse, + subcatch, + self.Soil, + 2000.0, + ) + self.thetaR = self.readtblDefault( + self.Dir + "/" + self.intbl + "/thetaR.tbl", + self.LandUse, + subcatch, + self.Soil, + 0.01, + ) + self.thetaS = self.readtblDefault( + self.Dir + "/" + self.intbl + "/thetaS.tbl", + self.LandUse, + subcatch, + self.Soil, + 0.6, + ) # minimum thickness of soild - self.SoilMinThickness = self.readtblDefault(self.Dir + "/" + self.intbl + "/SoilMinThickness.tbl", - self.LandUse, subcatch, self.Soil, 500.0) + self.SoilMinThickness = self.readtblDefault( + self.Dir + "/" + self.intbl + "/SoilMinThickness.tbl", + self.LandUse, + subcatch, + self.Soil, + 500.0, + ) # KsatVer = $2\inmaps\KsatVer.map - self.KsatVer = self.readtblDefault(self.Dir + "/" + self.intbl + "/KsatVer.tbl", self.LandUse, - subcatch, self.Soil, 3000.0) * self.timestepsecs / self.basetimestep - self.MporeFrac = self.readtblDefault(self.Dir + "/" + self.intbl + "/MporeFrac.tbl", self.LandUse, - subcatch, self.Soil, 0.0) + self.KsatVer = ( + self.readtblDefault( + self.Dir + "/" + self.intbl + "/KsatVer.tbl", + self.LandUse, + subcatch, + self.Soil, + 3000.0, + ) + * self.timestepsecs + / self.basetimestep + ) + self.MporeFrac = self.readtblDefault( + self.Dir + "/" + self.intbl + "/MporeFrac.tbl", + self.LandUse, + subcatch, + self.Soil, + 0.0, + ) - self.KsatHorFrac = self.readtblDefault(self.Dir + "/" + self.intbl + "/KsatHorFrac.tbl", self.LandUse, - subcatch, self.Soil, 1.0) + self.KsatHorFrac = self.readtblDefault( + self.Dir + "/" + self.intbl + "/KsatHorFrac.tbl", + self.LandUse, + subcatch, + self.Soil, + 1.0, + ) # Check if we have irrigation areas tt = pcr2numpy(self.IrrigationAreas, 0.0) @@ -722,117 +1071,217 @@ self.Beta = scalar(0.6) # For sheetflow - self.M = self.readtblDefault(self.Dir + "/" + self.intbl + "/M.tbl", self.LandUse, subcatch, self.Soil, - 300.0) # Decay parameter in Topog_sbm - self.N = self.readtblDefault(self.Dir + "/" + self.intbl + "/N.tbl", self.LandUse, subcatch, self.Soil, - 0.072) # Manning overland flow + self.M = self.readtblDefault( + self.Dir + "/" + self.intbl + "/M.tbl", + self.LandUse, + subcatch, + self.Soil, + 300.0, + ) # Decay parameter in Topog_sbm + self.N = self.readtblDefault( + self.Dir + "/" + self.intbl + "/N.tbl", + self.LandUse, + subcatch, + self.Soil, + 0.072, + ) # Manning overland flow if self.NRiverMethod == 1: - self.NRiver = self.readtblDefault(self.Dir + "/" + self.intbl + "/N_River.tbl", self.LandUse, subcatch, - self.Soil, 0.036) # Manning river + self.NRiver = self.readtblDefault( + self.Dir + "/" + self.intbl + "/N_River.tbl", + self.LandUse, + subcatch, + self.Soil, + 0.036, + ) # Manning river if self.NRiverMethod == 2: - self.NRiver = self.readtblFlexDefault(self.Dir + "/" + self.intbl + "/N_River.tbl", - 0.036, wflow_streamorder) - - self.WaterFrac = self.readtblDefault(self.Dir + "/" + self.intbl + "/WaterFrac.tbl", self.LandUse, subcatch, - self.Soil, 0.0) # Fraction Open water - self.et_RefToPot = self.readtblDefault(self.Dir + "/" + self.intbl + "/et_reftopot.tbl", self.LandUse, subcatch, - self.Soil, 1.0) # Fraction Open water + self.NRiver = self.readtblFlexDefault( + self.Dir + "/" + self.intbl + "/N_River.tbl", 0.036, wflow_streamorder + ) + self.WaterFrac = self.readtblDefault( + self.Dir + "/" + self.intbl + "/WaterFrac.tbl", + self.LandUse, + subcatch, + self.Soil, + 0.0, + ) # Fraction Open water + self.et_RefToPot = self.readtblDefault( + self.Dir + "/" + self.intbl + "/et_reftopot.tbl", + self.LandUse, + subcatch, + self.Soil, + 1.0, + ) # Fraction Open water if self.modelSnow: # HBV Snow parameters # critical temperature for snowmelt and refreezing: TTI= 1.000 - self.TTI = self.readtblDefault(self.Dir + "/" + self.intbl + "/TTI.tbl", self.LandUse, subcatch, self.Soil, - 1.0) + self.TTI = self.readtblDefault( + self.Dir + "/" + self.intbl + "/TTI.tbl", + self.LandUse, + subcatch, + self.Soil, + 1.0, + ) # TT = -1.41934 # defines interval in which precipitation falls as rainfall and snowfall - self.TT = self.readtblDefault(self.Dir + "/" + self.intbl + "/TT.tbl", self.LandUse, subcatch, self.Soil, - -1.41934) - self.TTM = self.readtblDefault(self.Dir + "/" + self.intbl + "/TTM.tbl", self.LandUse, subcatch, self.Soil, - -1.41934) - #Cfmax = 3.75653 # meltconstant in temperature-index - self.Cfmax = self.readtblDefault(self.Dir + "/" + self.intbl + "/Cfmax.tbl", self.LandUse, subcatch, - self.Soil, 3.75653) + self.TT = self.readtblDefault( + self.Dir + "/" + self.intbl + "/TT.tbl", + self.LandUse, + subcatch, + self.Soil, + -1.41934, + ) + self.TTM = self.readtblDefault( + self.Dir + "/" + self.intbl + "/TTM.tbl", + self.LandUse, + subcatch, + self.Soil, + -1.41934, + ) + # Cfmax = 3.75653 # meltconstant in temperature-index + self.Cfmax = self.readtblDefault( + self.Dir + "/" + self.intbl + "/Cfmax.tbl", + self.LandUse, + subcatch, + self.Soil, + 3.75653, + ) # WHC= 0.10000 # fraction of Snowvolume that can store water - self.WHC = self.readtblDefault(self.Dir + "/" + self.intbl + "/WHC.tbl", self.LandUse, subcatch, self.Soil, - 0.1) + self.WHC = self.readtblDefault( + self.Dir + "/" + self.intbl + "/WHC.tbl", + self.LandUse, + subcatch, + self.Soil, + 0.1, + ) # Wigmosta, M. S., L. J. Lane, J. D. Tagestad, and A. M. Coleman (2009). - self.w_soil = self.readtblDefault(self.Dir + "/" + self.intbl + "/w_soil.tbl", self.LandUse, subcatch, - self.Soil, 0.9 * 3.0 / 24.0) * self.timestepsecs / self.basetimestep + self.w_soil = ( + self.readtblDefault( + self.Dir + "/" + self.intbl + "/w_soil.tbl", + self.LandUse, + subcatch, + self.Soil, + 0.9 * 3.0 / 24.0, + ) + * self.timestepsecs + / self.basetimestep + ) - self.cf_soil = min(0.99, - self.readtblDefault(self.Dir + "/" + self.intbl + "/cf_soil.tbl", self.LandUse, subcatch, - self.Soil, 0.038)) # Ksat reduction factor fro frozen soi + self.cf_soil = min( + 0.99, + self.readtblDefault( + self.Dir + "/" + self.intbl + "/cf_soil.tbl", + self.LandUse, + subcatch, + self.Soil, + 0.038, + ), + ) # Ksat reduction factor fro frozen soi # We are modelling gletchers # Determine real slope and cell length - self.xl, self.yl, self.reallength = pcrut.detRealCellLength(self.ZeroMap, sizeinmetres) + self.xl, self.yl, self.reallength = pcrut.detRealCellLength( + self.ZeroMap, sizeinmetres + ) self.Slope = slope(self.Altitude) - #self.Slope=ifthen(boolean(self.TopoId),max(0.001,self.Slope*celllength()/self.reallength)) + # self.Slope=ifthen(boolean(self.TopoId),max(0.001,self.Slope*celllength()/self.reallength)) self.Slope = max(0.00001, self.Slope * celllength() / self.reallength) Terrain_angle = scalar(atan(self.Slope)) - - self.N = ifthenelse(self.River, self.NRiver, self.N) - if (hasattr(self,'ReserVoirSimpleLocs') or hasattr(self,'ReserVoirComplexLocs')): + if hasattr(self, "ReserVoirSimpleLocs") or hasattr( + self, "ReserVoirComplexLocs" + ): self.ReserVoirLocs = self.ZeroMap self.filter_P_PET = self.ZeroMap + 1.0 - if hasattr(self,'ReserVoirSimpleLocs'): + if hasattr(self, "ReserVoirSimpleLocs"): # Check if we have simple and or complex reservoirs tt_simple = pcr2numpy(self.ReserVoirSimpleLocs, 0.0) self.nrresSimple = tt_simple.max() - self.ReserVoirLocs = self.ReserVoirLocs + cover(scalar(self.ReserVoirSimpleLocs)) + self.ReserVoirLocs = self.ReserVoirLocs + cover( + scalar(self.ReserVoirSimpleLocs) + ) areamap = self.reallength * self.reallength res_area = areatotal(spatial(areamap), self.ReservoirSimpleAreas) resarea_pnt = ifthen(boolean(self.ReserVoirSimpleLocs), res_area) - self.ResSimpleArea = ifthenelse(cover(self.ResSimpleArea,scalar(0.0)) > 0, self.ResSimpleArea, cover(resarea_pnt,scalar(0.0))) - self.filter_P_PET = ifthenelse(boolean(cover(res_area,scalar(0.0))), res_area*0.0, self.filter_P_PET) + self.ResSimpleArea = ifthenelse( + cover(self.ResSimpleArea, scalar(0.0)) > 0, + self.ResSimpleArea, + cover(resarea_pnt, scalar(0.0)), + ) + self.filter_P_PET = ifthenelse( + boolean(cover(res_area, scalar(0.0))), res_area * 0.0, self.filter_P_PET + ) else: self.nrresSimple = 0 - - - if hasattr(self, 'ReserVoirComplexLocs'): + + if hasattr(self, "ReserVoirComplexLocs"): tt_complex = pcr2numpy(self.ReserVoirComplexLocs, 0.0) self.nrresComplex = tt_complex.max() - self.ReserVoirLocs = self.ReserVoirLocs + cover(scalar(self.ReserVoirComplexLocs)) - res_area = cover(scalar(self.ReservoirComplexAreas),0.0) - self.filter_P_PET = ifthenelse(res_area > 0, res_area*0.0, self.filter_P_PET) - - #read files + self.ReserVoirLocs = self.ReserVoirLocs + cover( + scalar(self.ReserVoirComplexLocs) + ) + res_area = cover(scalar(self.ReservoirComplexAreas), 0.0) + self.filter_P_PET = ifthenelse( + res_area > 0, res_area * 0.0, self.filter_P_PET + ) + + # read files self.sh = {} res_ids = ifthen(self.ResStorFunc == 2, self.ReserVoirComplexLocs) - np_res_ids = pcr2numpy(res_ids,0) + np_res_ids = pcr2numpy(res_ids, 0) np_res_ids_u = np.unique(np_res_ids[nonzero(np_res_ids)]) if np.size(np_res_ids_u) > 0: for item in nditer(np_res_ids_u): - self.sh[int(item)] = loadtxt(self.Dir + "/" + self.intbl + "/Reservoir_SH_" + str(item) + ".tbl") + self.sh[int(item)] = loadtxt( + self.Dir + + "/" + + self.intbl + + "/Reservoir_SH_" + + str(item) + + ".tbl" + ) self.hq = {} res_ids = ifthen(self.ResOutflowFunc == 1, self.ReserVoirComplexLocs) - np_res_ids = pcr2numpy(res_ids,0) + np_res_ids = pcr2numpy(res_ids, 0) np_res_ids_u = np.unique(np_res_ids[nonzero(np_res_ids)]) if size(np_res_ids_u) > 0: for item in nditer(np_res_ids_u): - self.hq[int(item)] = loadtxt(self.Dir + "/" + self.intbl + "/Reservoir_HQ_" + str(item) + ".tbl", skiprows=3) - + self.hq[int(item)] = loadtxt( + self.Dir + + "/" + + self.intbl + + "/Reservoir_HQ_" + + str(item) + + ".tbl", + skiprows=3, + ) + else: self.nrresComplex = 0 - - if (self.nrresSimple + self.nrresComplex) > 0: - self.ReserVoirLocs =ordinal(self.ReserVoirLocs) - self.logger.info("A total of " + str(self.nrresSimple) + " simple reservoirs and " + str(self.nrresComplex) + " complex reservoirs found.") + if (self.nrresSimple + self.nrresComplex) > 0: + self.ReserVoirLocs = ordinal(self.ReserVoirLocs) + self.logger.info( + "A total of " + + str(self.nrresSimple) + + " simple reservoirs and " + + str(self.nrresComplex) + + " complex reservoirs found." + ) self.ReserVoirDownstreamLocs = downstream(self.TopoLdd, self.ReserVoirLocs) self.TopoLddOrg = self.TopoLdd - self.TopoLdd = lddrepair(cover(ifthen(boolean(self.ReserVoirLocs), ldd(5)), self.TopoLdd)) + self.TopoLdd = lddrepair( + cover(ifthen(boolean(self.ReserVoirLocs), ldd(5)), self.TopoLdd) + ) tt_filter = pcr2numpy(self.filter_P_PET, 1.0) self.filterResArea = tt_filter.min() - # Determine river width from DEM, upstream area and yearly average discharge # Scale yearly average Q at outlet with upstream are to get Q over whole catchment # Alf ranges from 5 to > 60. 5 for hardrock. large values for sediments @@ -841,93 +1290,156 @@ upstr = catchmenttotal(1, self.TopoLdd) Qscale = upstr / mapmaximum(upstr) * Qmax - W = (alf * (alf + 2.0) ** (0.6666666667)) ** (0.375) * Qscale ** (0.375) * ( - max(0.0001, windowaverage(self.Slope, celllength() * 4.0))) ** (-0.1875) * self.N ** (0.375) + W = ( + (alf * (alf + 2.0) ** (0.6666666667)) ** (0.375) + * Qscale ** (0.375) + * (max(0.0001, windowaverage(self.Slope, celllength() * 4.0))) ** (-0.1875) + * self.N ** (0.375) + ) # Use supplied riverwidth if possible, else calulate self.RiverWidth = ifthenelse(self.RiverWidth <= 0.0, W, self.RiverWidth) # Only allow reinfiltration in river cells by default - if not hasattr(self,'MaxReinfilt'): - self.MaxReinfilt = ifthenelse(self.River, self.ZeroMap + 999.0, self.ZeroMap) + if not hasattr(self, "MaxReinfilt"): + self.MaxReinfilt = ifthenelse( + self.River, self.ZeroMap + 999.0, self.ZeroMap + ) # soil thickness based on topographical index (see Environmental modelling: finding simplicity in complexity) # 1: calculate wetness index # 2: Scale the capacity (now actually a max capacity) based on the index, also apply a minmum capacity - WI = ln(accuflux(self.TopoLdd, - 1) / self.Slope) # Topographical wetnesss. Scale WI by zone/subcatchment assuming these ara also geological units + WI = ln( + accuflux(self.TopoLdd, 1) / self.Slope + ) # Topographical wetnesss. Scale WI by zone/subcatchment assuming these ara also geological units WIMax = areamaximum(WI, self.TopoId) * WIMaxScale - self.SoilThickness = max(min(self.SoilThickness, (WI / WIMax) * self.SoilThickness), - self.SoilMinThickness) + self.SoilThickness = max( + min(self.SoilThickness, (WI / WIMax) * self.SoilThickness), + self.SoilMinThickness, + ) self.SoilWaterCapacity = self.SoilThickness * (self.thetaS - self.thetaR) # determine number of layers based on total soil thickness # assign thickness, unsaturated water store and transfer to these layers (initializing) - UStoreLayerThickness = configget(self.config,"model","UStoreLayerThickness",'0') - if UStoreLayerThickness != '0': - self.USatLayers = len(UStoreLayerThickness.split(',')) + UStoreLayerThickness = configget( + self.config, "model", "UStoreLayerThickness", "0" + ) + if UStoreLayerThickness != "0": + self.USatLayers = len(UStoreLayerThickness.split(",")) self.maxLayers = self.USatLayers + 1 else: UStoreLayerThickness = self.SoilThickness self.USatLayers = 1 self.maxLayers = self.USatLayers - - self.UStoreLayerThickness= [] - self.UStoreLayerDepth= [] + self.UStoreLayerThickness = [] + self.UStoreLayerDepth = [] self.T = [] self.maskLayer = [] self.SumThickness = self.ZeroMap self.nrLayersMap = self.ZeroMap - - for n in arange(0,self.maxLayers): + for n in arange(0, self.maxLayers): self.SumLayer = self.SumThickness if self.USatLayers > 1 and n < self.USatLayers: - UstoreThick_temp = float(UStoreLayerThickness.split(',')[n])+self.ZeroMap - UstoreThick = min(UstoreThick_temp,max(self.SoilThickness-self.SumLayer,0.0)) + UstoreThick_temp = ( + float(UStoreLayerThickness.split(",")[n]) + self.ZeroMap + ) + UstoreThick = min( + UstoreThick_temp, max(self.SoilThickness - self.SumLayer, 0.0) + ) else: UstoreThick_temp = mapmaximum(self.SoilThickness) - self.SumLayer - UstoreThick = min(UstoreThick_temp,max(self.SoilThickness-self.SumLayer,0.0)) - + UstoreThick = min( + UstoreThick_temp, max(self.SoilThickness - self.SumLayer, 0.0) + ) + self.SumThickness = UstoreThick_temp + self.SumThickness - self.nrLayersMap = ifthenelse((self.SoilThickness>=self.SumThickness) | (self.SoilThickness-self.SumLayer>self.ZeroMap) , self.nrLayersMap + 1 ,self.nrLayersMap) + self.nrLayersMap = ifthenelse( + (self.SoilThickness >= self.SumThickness) + | (self.SoilThickness - self.SumLayer > self.ZeroMap), + self.nrLayersMap + 1, + self.nrLayersMap, + ) - self.UStoreLayerThickness.append(ifthenelse((self.SumThickness<=self.SoilThickness) | (self.SoilThickness-self.SumLayer>self.ZeroMap),UstoreThick,0.0)) - self.UStoreLayerDepth.append(ifthen((self.SumThickness<=self.SoilThickness) | (self.SoilThickness-self.SumLayer>self.ZeroMap), self.SoilThickness*0.0)) - self.T.append(ifthen((self.SumThickness<=self.SoilThickness) | (self.SoilThickness-self.SumLayer>self.ZeroMap), self.SoilThickness*0.0)) - self.maskLayer.append(ifthen((self.SumThickness<=self.SoilThickness) | (self.SoilThickness-self.SumLayer>self.ZeroMap), self.SoilThickness*0.0)) + self.UStoreLayerThickness.append( + ifthenelse( + (self.SumThickness <= self.SoilThickness) + | (self.SoilThickness - self.SumLayer > self.ZeroMap), + UstoreThick, + 0.0, + ) + ) + self.UStoreLayerDepth.append( + ifthen( + (self.SumThickness <= self.SoilThickness) + | (self.SoilThickness - self.SumLayer > self.ZeroMap), + self.SoilThickness * 0.0, + ) + ) + self.T.append( + ifthen( + (self.SumThickness <= self.SoilThickness) + | (self.SoilThickness - self.SumLayer > self.ZeroMap), + self.SoilThickness * 0.0, + ) + ) + self.maskLayer.append( + ifthen( + (self.SumThickness <= self.SoilThickness) + | (self.SoilThickness - self.SumLayer > self.ZeroMap), + self.SoilThickness * 0.0, + ) + ) self.KsatVerFrac = [] self.c = [] - for n in arange(0,len(self.UStoreLayerThickness)): - self.KsatVerFrac.append(self.readtblLayersDefault(self.Dir + "/" + self.intbl + "/KsatVerFrac.tbl", self.LandUse, - subcatch, self.Soil, n, 1.0)) - self.c.append(self.readtblLayersDefault(self.Dir + "/" + self.intbl + "/c.tbl", self.LandUse, - subcatch, self.Soil, n, 10.0)) + for n in arange(0, len(self.UStoreLayerThickness)): + self.KsatVerFrac.append( + self.readtblLayersDefault( + self.Dir + "/" + self.intbl + "/KsatVerFrac.tbl", + self.LandUse, + subcatch, + self.Soil, + n, + 1.0, + ) + ) + self.c.append( + self.readtblLayersDefault( + self.Dir + "/" + self.intbl + "/c.tbl", + self.LandUse, + subcatch, + self.Soil, + n, + 10.0, + ) + ) - # limit roots to top 99% of first zone self.RootingDepth = min(self.SoilThickness * 0.99, self.RootingDepth) # subgrid runoff generation, determine CC (sharpness of S-Curve) for upper # en lower part and take average self.DemMax = readmap(self.Dir + "/staticmaps/wflow_demmax") self.DrainageBase = readmap(self.Dir + "/staticmaps/wflow_demmin") - self.CClow = min(100.0, - ln(1.0 / 0.1 - 1) / min(-0.1, self.DrainageBase - self.Altitude)) - self.CCup = min(100.0, - ln(1.0 / 0.1 - 1) / min(-0.1, self.Altitude - self.DemMax)) + self.CClow = min( + 100.0, -ln(1.0 / 0.1 - 1) / min(-0.1, self.DrainageBase - self.Altitude) + ) + self.CCup = min( + 100.0, -ln(1.0 / 0.1 - 1) / min(-0.1, self.Altitude - self.DemMax) + ) self.CC = (self.CClow + self.CCup) * 0.5 - # Which columns/gauges to use/ignore in updating self.UpdateMap = self.ZeroMap if self.updating: _tmp = pcr2numpy(self.OutputLoc, 0.0) gaugear = _tmp - touse = numpy.zeros(gaugear.shape, dtype='int') + touse = numpy.zeros(gaugear.shape, dtype="int") for thecol in updateCols: idx = (gaugear == thecol).nonzero() @@ -937,10 +1449,15 @@ # Calculate distance to updating points (upstream) annd use to scale the correction # ldddist returns zero for cell at the gauges so add 1.0 tp result self.DistToUpdPt = cover( - min(ldddist(self.TopoLdd, boolean(cover(self.UpdateMap, 0)), 1) * self.reallength / celllength(), - self.UpdMaxDist), self.UpdMaxDist) + min( + ldddist(self.TopoLdd, boolean(cover(self.UpdateMap, 0)), 1) + * self.reallength + / celllength(), + self.UpdMaxDist, + ), + self.UpdMaxDist, + ) - # Initializing of variables self.logger.info("Initializing of model variables..") self.TopoLdd = lddmask(self.TopoLdd, boolean(self.TopoId)) @@ -950,24 +1467,34 @@ # This is very handy for Ribasim etc... if self.SubCatchFlowOnly > 0: self.logger.info("Creating subcatchment-only drainage network (ldd)") - ds = downstream(self.TopoLdd,self.TopoId) - usid = ifthenelse(ds != self.TopoId,self.TopoId,0) - self.TopoLdd = lddrepair(ifthenelse(boolean(usid),ldd(5),self.TopoLdd)) + ds = downstream(self.TopoLdd, self.TopoId) + usid = ifthenelse(ds != self.TopoId, self.TopoId, 0) + self.TopoLdd = lddrepair(ifthenelse(boolean(usid), ldd(5), self.TopoLdd)) # Used to seperate output per LandUse/management classes OutZones = self.LandUse - self.QMMConv = self.timestepsecs / (self.reallength * self.reallength * 0.001) #m3/s --> actial mm of water over the cell - #self.QMMConvUp = 1000.0 * self.timestepsecs / ( catchmenttotal(cover(1.0), self.TopoLdd) * self.reallength * self.reallength) #m3/s --> mm over upstreams - temp = catchmenttotal(cover(1.0), self.TopoLdd) * self.reallength * 0.001 * 0.001 * self.reallength - self.QMMConvUp = cover(self.timestepsecs * 0.001)/temp - self.ToCubic = (self.reallength * self.reallength * 0.001) / self.timestepsecs # m3/s + self.QMMConv = self.timestepsecs / ( + self.reallength * self.reallength * 0.001 + ) # m3/s --> actial mm of water over the cell + # self.QMMConvUp = 1000.0 * self.timestepsecs / ( catchmenttotal(cover(1.0), self.TopoLdd) * self.reallength * self.reallength) #m3/s --> mm over upstreams + temp = ( + catchmenttotal(cover(1.0), self.TopoLdd) + * self.reallength + * 0.001 + * 0.001 + * self.reallength + ) + self.QMMConvUp = cover(self.timestepsecs * 0.001) / temp + self.ToCubic = ( + self.reallength * self.reallength * 0.001 + ) / self.timestepsecs # m3/s self.KinWaveVolume = self.ZeroMap self.OldKinWaveVolume = self.ZeroMap - self.sumprecip = self.ZeroMap #accumulated rainfall for water balance - self.sumevap = self.ZeroMap #accumulated evaporation for water balance - self.sumrunoff = self.ZeroMap #accumulated runoff for water balance - self.sumint = self.ZeroMap #accumulated interception for water balance + self.sumprecip = self.ZeroMap # accumulated rainfall for water balance + self.sumevap = self.ZeroMap # accumulated evaporation for water balance + self.sumrunoff = self.ZeroMap # accumulated runoff for water balance + self.sumint = self.ZeroMap # accumulated interception for water balance self.sumleakage = self.ZeroMap self.CumReinfilt = self.ZeroMap self.sumoutflow = self.ZeroMap @@ -999,7 +1526,9 @@ self.Aspect = scalar(aspect(self.Altitude)) # aspect [deg] self.Aspect = ifthenelse(self.Aspect <= 0.0, scalar(0.001), self.Aspect) # On Flat areas the Aspect function fails, fill in with average... - self.Aspect = ifthenelse(defined(self.Aspect), self.Aspect, areaaverage(self.Aspect, self.TopoId)) + self.Aspect = ifthenelse( + defined(self.Aspect), self.Aspect, areaaverage(self.Aspect, self.TopoId) + ) # Set DCL to riverlength if that is longer that the basic length calculated from grid drainlength = detdrainlength(self.TopoLdd, self.xl, self.yl) @@ -1017,14 +1546,19 @@ self.Bw = ifthenelse(self.River, self.RiverWidth, self.Bw) # Add rivers to the WaterFrac, but check with waterfrac map and correct - self.RiverFrac = min(1.0, ifthenelse(self.River, (self.RiverWidth * self.DCL) / (self.xl * self.yl), 0)) - self.WaterFrac = min(1.0,self.WaterFrac + self.RiverFrac) + self.RiverFrac = min( + 1.0, + ifthenelse( + self.River, (self.RiverWidth * self.DCL) / (self.xl * self.yl), 0 + ), + ) + self.WaterFrac = min(1.0, self.WaterFrac + self.RiverFrac) # term for Alpha # Correct slope for extra length of the river in a gridcel riverslopecor = drainlength / self.DCL - #report(riverslopecor,"cor.map") - #report(self.Slope * riverslopecor,"slope.map") + # report(riverslopecor,"cor.map") + # report(self.Slope * riverslopecor,"slope.map") self.AlpTerm = pow((self.N / (sqrt(self.Slope * riverslopecor))), self.Beta) # power for Alpha self.AlpPow = (2.0 / 3.0) * self.Beta @@ -1036,50 +1570,58 @@ # Save some summary maps self.logger.info("Saving summary maps...") - #self.IF = self.ZeroMap + # self.IF = self.ZeroMap self.logger.info("End of initial section") - - def default_summarymaps(self): - """ + """ Returns a list of default summary-maps at the end of a run. This is model specific. You can also add them to the [summary]section of the ini file but stuff you think is crucial to the model should be listed here """ - lst = ['self.RiverWidth', - 'self.Cmax', 'self.csize', 'self.upsize', - 'self.EoverR', 'self.RootingDepth', - 'self.CanopyGapFraction', 'self.InfiltCapSoil', - 'self.InfiltCapPath', - 'self.PathFrac', - 'self.thetaR', - 'self.thetaS', - 'self.SoilMinThickness', 'self.SoilThickness', 'self.nrLayersMap', - 'self.KsatVer', - 'self.M', - 'self.SoilWaterCapacity', - 'self.et_RefToPot', - 'self.Slope', - 'self.CC', - 'self.N', - 'self.RiverFrac', - 'self.WaterFrac', - 'self.xl', 'self.yl', 'self.reallength', - 'self.DCL', - 'self.Bw', - 'self.PathInfiltExceeded','self.SoilInfiltExceeded'] + lst = [ + "self.RiverWidth", + "self.Cmax", + "self.csize", + "self.upsize", + "self.EoverR", + "self.RootingDepth", + "self.CanopyGapFraction", + "self.InfiltCapSoil", + "self.InfiltCapPath", + "self.PathFrac", + "self.thetaR", + "self.thetaS", + "self.SoilMinThickness", + "self.SoilThickness", + "self.nrLayersMap", + "self.KsatVer", + "self.M", + "self.SoilWaterCapacity", + "self.et_RefToPot", + "self.Slope", + "self.CC", + "self.N", + "self.RiverFrac", + "self.WaterFrac", + "self.xl", + "self.yl", + "self.reallength", + "self.DCL", + "self.Bw", + "self.PathInfiltExceeded", + "self.SoilInfiltExceeded", + ] - return lst + return lst - def resume(self): if self.reinit == 1: self.logger.info("Setting initial conditions to default") self.SatWaterDepth = self.SoilWaterCapacity * 0.85 - #for n in arange(0,self.nrLayers): + # for n in arange(0,self.nrLayers): # self.UStoreLayerDepth[n] = self.ZeroMap # TODO: move UStoreLayerDepth from initial to here @@ -1089,21 +1631,22 @@ self.SnowWater = self.ZeroMap self.TSoil = self.ZeroMap + 10.0 self.CanopyStorage = self.ZeroMap - if hasattr(self, 'ReserVoirSimpleLocs'): + if hasattr(self, "ReserVoirSimpleLocs"): self.ReservoirVolume = self.ResMaxVolume * self.ResTargetFullFrac - if hasattr(self, 'ReserVoirComplexLocs'): + if hasattr(self, "ReserVoirComplexLocs"): self.ReservoirWaterLevel = cover(0.0) - if hasattr(self, 'GlacierFrac'): - self.GlacierStore = self.wf_readmap(os.path.join(self.Dir,"staticmaps","GlacierStore.map"), 55.0 * 1000) + if hasattr(self, "GlacierFrac"): + self.GlacierStore = self.wf_readmap( + os.path.join(self.Dir, "staticmaps", "GlacierStore.map"), + 55.0 * 1000, + ) if self.nrpaddyirri > 0: self.PondingDepth = self.ZeroMap else: self.logger.info("Setting initial conditions from state files") - self.wf_resume(os.path.join(self.Dir,"instate")) + self.wf_resume(os.path.join(self.Dir, "instate")) - - P = self.Bw + (2.0 * self.WaterLevel) self.Alpha = self.AlpTerm * pow(P, self.AlpPow) self.OldSurfaceRunoff = self.SurfaceRunoff @@ -1114,19 +1657,30 @@ self.OldKinWaveVolume = self.KinWaveVolume self.QCatchmentMM = self.SurfaceRunoff * self.QMMConvUp - self.InitialStorage = self.SatWaterDepth + sum_list_cover(self.UStoreLayerDepth,self.ZeroMap) + self.CanopyStorage - self.CellStorage = self.SatWaterDepth + sum_list_cover(self.UStoreLayerDepth,self.ZeroMap) + self.InitialStorage = ( + self.SatWaterDepth + + sum_list_cover(self.UStoreLayerDepth, self.ZeroMap) + + self.CanopyStorage + ) + self.CellStorage = self.SatWaterDepth + sum_list_cover( + self.UStoreLayerDepth, self.ZeroMap + ) # Determine actual water depth - self.zi = max(0.0, self.SoilThickness - self.SatWaterDepth / (self.thetaS - self.thetaR)) + self.zi = max( + 0.0, self.SoilThickness - self.SatWaterDepth / (self.thetaS - self.thetaR) + ) # TOPOG_SBM type soil stuff self.f = (self.thetaS - self.thetaR) / self.M # NOTE:: This line used to be in the initial section. As a result # simulations will now be different as it used to be before # the rescaling of the FirstZoneThickness - self.GWScale = (self.DemMax - self.DrainageBase) / self.SoilThickness / self.RunoffGeneratingGWPerc + self.GWScale = ( + (self.DemMax - self.DrainageBase) + / self.SoilThickness + / self.RunoffGeneratingGWPerc + ) - def dynamic(self): """ Stuf that is done for each timestep of the model @@ -1183,32 +1737,34 @@ :var self.ToCubic: Mutiplier to convert mm to m^3/s for fluxes """ - # Read forcing data and dynamic parameters self.wf_updateparameters() - self.Precipitation = max(0.0,self.Precipitation) + self.Precipitation = max(0.0, self.Precipitation) # NB This may interfere with lintul link - if hasattr(self,"LAI"): + if hasattr(self, "LAI"): # Sl must also be defined ##TODO: add MAXLAI and CWf self.Cmax = self.Sl * self.LAI + self.Swood self.CanopyGapFraction = exp(-self.Kext * self.LAI) self.Ewet = (1 - exp(-self.Kext * self.LAI)) * self.PotenEvap - self.EoverR = ifthenelse(self.Precipitation > 0.0, \ - min(0.25,cover(self.Ewet/max(0.0001,self.Precipitation),0.0)), 0.0) - if hasattr(self,'MAXLAI') and hasattr(self,'CWf'): + self.EoverR = ifthenelse( + self.Precipitation > 0.0, + min(0.25, cover(self.Ewet / max(0.0001, self.Precipitation), 0.0)), + 0.0, + ) + if hasattr(self, "MAXLAI") and hasattr(self, "CWf"): # Adjust rootinggdepth - self.ActRootingDepth = self.CWf * (self.RootingDepth * self.LAI/max(0.001,self.MAXLAI))\ - + ((1- self.CWf) * self.RootingDepth) + self.ActRootingDepth = self.CWf * ( + self.RootingDepth * self.LAI / max(0.001, self.MAXLAI) + ) + ((1 - self.CWf) * self.RootingDepth) else: self.ActRootingDepth = self.RootingDepth else: self.ActRootingDepth = self.RootingDepth - - #Apply forcing data corrections + # Apply forcing data corrections self.PotenEvap = self.PotenEvap * self.et_RefToPot if self.modelSnow: self.Temperature = self.Temperature + self.TempCor @@ -1218,11 +1774,13 @@ if (self.nrresSimple + self.nrresComplex) > 0 and self.filterResArea == 0: self.ReserVoirPotEvap = self.PotenEvap self.ReserVoirPrecip = self.Precipitation - + self.PotenEvap = self.filter_P_PET * self.PotenEvap self.Precipitation = self.filter_P_PET * self.Precipitation - self.OrgStorage = sum_list_cover(self.UStoreLayerDepth,self.ZeroMap) + self.SatWaterDepth + self.OrgStorage = ( + sum_list_cover(self.UStoreLayerDepth, self.ZeroMap) + self.SatWaterDepth + ) self.OldCanopyStorage = self.CanopyStorage if self.nrpaddyirri > 0: self.OldPondingDepth = self.PondingDepth @@ -1231,95 +1789,154 @@ if self.modelSnow: self.TSoil = self.TSoil + self.w_soil * (self.Temperature - self.TSoil) # return Snow,SnowWater,SnowMelt,RainFall - self.Snow, self.SnowWater, self.SnowMelt, self.PrecipitationPlusMelt,self.SnowFall = SnowPackHBV(self.Snow, self.SnowWater, - self.Precipitation, - self.Temperature, self.TTI, - self.TT, self.TTM, self.Cfmax, self.WHC) + self.Snow, self.SnowWater, self.SnowMelt, self.PrecipitationPlusMelt, self.SnowFall = SnowPackHBV( + self.Snow, + self.SnowWater, + self.Precipitation, + self.Temperature, + self.TTI, + self.TT, + self.TTM, + self.Cfmax, + self.WHC, + ) MaxSnowPack = 10000.0 if self.MassWasting: # Masswasting of dry snow # 5.67 = tan 80 graden - SnowFluxFrac = min(0.5,self.Slope/5.67) * min(1.0,self.Snow/MaxSnowPack) + SnowFluxFrac = min(0.5, self.Slope / 5.67) * min( + 1.0, self.Snow / MaxSnowPack + ) MaxFlux = SnowFluxFrac * self.Snow - self.Snow = accucapacitystate(self.TopoLdd,self.Snow, MaxFlux) + self.Snow = accucapacitystate(self.TopoLdd, self.Snow, MaxFlux) else: SnowFluxFrac = self.ZeroMap - MaxFlux= self.ZeroMap + MaxFlux = self.ZeroMap - self.SnowCover = ifthenelse(self.Snow >0, scalar(1), scalar(0)) - self.NrCell= areatotal(self.SnowCover,self.TopoId) + self.SnowCover = ifthenelse(self.Snow > 0, scalar(1), scalar(0)) + self.NrCell = areatotal(self.SnowCover, self.TopoId) - if hasattr(self,'GlacierFrac'): + if hasattr(self, "GlacierFrac"): """ Run Glacier module and add the snowpack on-top of it. Snow becomes ice when pressure is about 830 k/m^2, e.g 8300 mm If below that a max amount of 2mm/day can be converted to glacier-ice """ - #TODO: document glacier module - self.snowdist = sCurve(self.Snow,a=8300.,c=0.06) - self.Snow2Glacier = ifthenelse(self.Snow > 8300, self.snowdist * (self.Snow - 8300), self.ZeroMap) + # TODO: document glacier module + self.snowdist = sCurve(self.Snow, a=8300., c=0.06) + self.Snow2Glacier = ifthenelse( + self.Snow > 8300, self.snowdist * (self.Snow - 8300), self.ZeroMap + ) - self.Snow2Glacier = ifthenelse(self.GlacierFrac > 0.0, self.Snow2Glacier,self.ZeroMap) + self.Snow2Glacier = ifthenelse( + self.GlacierFrac > 0.0, self.Snow2Glacier, self.ZeroMap + ) # Max conversion to 8mm/day - self.Snow2Glacier = min(self.Snow2Glacier,8.0) * self.timestepsecs/self.basetimestep + self.Snow2Glacier = ( + min(self.Snow2Glacier, 8.0) * self.timestepsecs / self.basetimestep + ) self.Snow = self.Snow - (self.Snow2Glacier * self.GlacierFrac) - self.GlacierStore, self.GlacierMelt = GlacierMelt(self.GlacierStore + self.Snow2Glacier,self.Snow,self.Temperature,\ - self.G_TT, self.G_Cfmax) + self.GlacierStore, self.GlacierMelt = GlacierMelt( + self.GlacierStore + self.Snow2Glacier, + self.Snow, + self.Temperature, + self.G_TT, + self.G_Cfmax, + ) # Convert to mm per grid cell and add to snowmelt self.GlacierMelt = self.GlacierMelt * self.GlacierFrac - self.PrecipitationPlusMelt = self.PrecipitationPlusMelt + self.GlacierMelt + self.PrecipitationPlusMelt = ( + self.PrecipitationPlusMelt + self.GlacierMelt + ) else: self.PrecipitationPlusMelt = self.Precipitation ########################################################################## # Interception according to a modified Gash model ########################################################################## if self.timestepsecs >= (23 * 3600): - self.ThroughFall, self.Interception, self.StemFlow, self.CanopyStorage = rainfall_interception_gash(self.Cmax, self.EoverR, - self.CanopyGapFraction, - self.PrecipitationPlusMelt, - self.CanopyStorage,maxevap=self.PotEvap) + self.ThroughFall, self.Interception, self.StemFlow, self.CanopyStorage = rainfall_interception_gash( + self.Cmax, + self.EoverR, + self.CanopyGapFraction, + self.PrecipitationPlusMelt, + self.CanopyStorage, + maxevap=self.PotEvap, + ) - self.PotTransSoil = cover(max(0.0, self.PotEvap - self.Interception), 0.0) # now in mm + self.PotTransSoil = cover( + max(0.0, self.PotEvap - self.Interception), 0.0 + ) # now in mm else: NetInterception, self.ThroughFall, self.StemFlow, LeftOver, Interception, self.CanopyStorage = rainfall_interception_modrut( - self.PrecipitationPlusMelt, self.PotEvap, self.CanopyStorage, self.CanopyGapFraction, self.Cmax) + self.PrecipitationPlusMelt, + self.PotEvap, + self.CanopyStorage, + self.CanopyGapFraction, + self.Cmax, + ) self.PotTransSoil = cover(max(0.0, LeftOver), 0.0) # now in mm - self.Interception=NetInterception + self.Interception = NetInterception # Start with the soil calculations # -------------------------------- # Code to be able to force zi from the outside # - self.SatWaterDepth = (self.thetaS - self.thetaR) * (self.SoilThickness - self.zi) + self.SatWaterDepth = (self.thetaS - self.thetaR) * ( + self.SoilThickness - self.zi + ) - self.AvailableForInfiltration = self.ThroughFall + self.StemFlow + self.IRSupplymm + self.AvailableForInfiltration = ( + self.ThroughFall + self.StemFlow + self.IRSupplymm + ) self.oldIRSupplymm = self.IRSupplymm - UStoreCapacity = self.SoilWaterCapacity - self.SatWaterDepth - sum_list_cover(self.UStoreLayerDepth,self.ZeroMap) + UStoreCapacity = ( + self.SoilWaterCapacity + - self.SatWaterDepth + - sum_list_cover(self.UStoreLayerDepth, self.ZeroMap) + ) # Runoff from water bodies and river network - self.RunoffOpenWater = min(1.0,self.RiverFrac + self.WaterFrac) * self.AvailableForInfiltration - #self.RunoffOpenWater = self.ZeroMap - self.AvailableForInfiltration = self.AvailableForInfiltration - self.RunoffOpenWater + self.RunoffOpenWater = ( + min(1.0, self.RiverFrac + self.WaterFrac) * self.AvailableForInfiltration + ) + # self.RunoffOpenWater = self.ZeroMap + self.AvailableForInfiltration = ( + self.AvailableForInfiltration - self.RunoffOpenWater + ) - if self.RunoffGenSigmaFunction: self.AbsoluteGW = self.DemMax - (self.zi * self.GWScale) # Determine saturated fraction of cell self.SubCellFrac = sCurve(self.AbsoluteGW, c=self.CC, a=self.Altitude + 1.0) # Make sure total of SubCellFRac + WaterFRac + RiverFrac <=1 to avoid double counting - Frac_correction = ifthenelse((self.SubCellFrac + self.RiverFrac + self.WaterFrac) > 1.0, - self.SubCellFrac + self.RiverFrac + self.WaterFrac - 1.0, 0.0) - self.SubCellRunoff = (self.SubCellFrac - Frac_correction) * self.AvailableForInfiltration - self.SubCellGWRunoff = min(self.SubCellFrac * self.SatWaterDepth, - max(0.0,self.SubCellFrac * self.Slope * self.KsatVer * \ - self.KsatHorFrac* exp(-self.f * self.zi))) + Frac_correction = ifthenelse( + (self.SubCellFrac + self.RiverFrac + self.WaterFrac) > 1.0, + self.SubCellFrac + self.RiverFrac + self.WaterFrac - 1.0, + 0.0, + ) + self.SubCellRunoff = ( + self.SubCellFrac - Frac_correction + ) * self.AvailableForInfiltration + self.SubCellGWRunoff = min( + self.SubCellFrac * self.SatWaterDepth, + max( + 0.0, + self.SubCellFrac + * self.Slope + * self.KsatVer + * self.KsatHorFrac + * exp(-self.f * self.zi), + ), + ) self.SatWaterDepth = self.SatWaterDepth - self.SubCellGWRunoff - self.AvailableForInfiltration = self.AvailableForInfiltration - self.SubCellRunoff + self.AvailableForInfiltration = ( + self.AvailableForInfiltration - self.SubCellRunoff + ) else: self.AbsoluteGW = self.DemMax - (self.zi * self.GWScale) self.SubCellFrac = spatial(scalar(0.0)) @@ -1337,19 +1954,22 @@ else: soilInfRedu = 1.0 MaxInfiltSoil = min(self.InfiltCapSoil * soilInfRedu, SoilInf) - self.SoilInfiltExceeded = self.SoilInfiltExceeded + scalar(self.InfiltCapSoil * soilInfRedu < SoilInf) + self.SoilInfiltExceeded = self.SoilInfiltExceeded + scalar( + self.InfiltCapSoil * soilInfRedu < SoilInf + ) MaxInfiltPath = min(self.InfiltCapPath * soilInfRedu, PathInf) - self.PathInfiltExceeded = self.PathInfiltExceeded + scalar(self.InfiltCapPath * soilInfRedu < PathInf) + self.PathInfiltExceeded = self.PathInfiltExceeded + scalar( + self.InfiltCapPath * soilInfRedu < PathInf + ) - InfiltSoilPath = min(MaxInfiltPath+MaxInfiltSoil,max(0.0,UStoreCapacity)) + InfiltSoilPath = min(MaxInfiltPath + MaxInfiltSoil, max(0.0, UStoreCapacity)) self.In = InfiltSoilPath - self.ActInfilt = InfiltSoilPath # JS Ad this to be compatible with rest + self.ActInfilt = InfiltSoilPath # JS Ad this to be compatible with rest self.SumThickness = self.ZeroMap self.ZiLayer = self.ZeroMap - # Limit rootingdepth (if set externally) self.ActRootingDepth = min(self.SoilThickness * 0.99, self.ActRootingDepth) @@ -1361,104 +1981,221 @@ # Determine Open Water EVAP. Later subtract this from water that # enters the Kinematic wave self.RestEvap = self.potsoilopenwaterevap - #self.RestEvap = (self.PotTrans - self.Transpiration) + self.potsoilopenwaterevap - self.ActEvapOpenWater = min(self.WaterLevel * 1000.0 * self.WaterFrac ,self.WaterFrac * self.RestEvap) + # self.RestEvap = (self.PotTrans - self.Transpiration) + self.potsoilopenwaterevap + self.ActEvapOpenWater = min( + self.WaterLevel * 1000.0 * self.WaterFrac, self.WaterFrac * self.RestEvap + ) self.RestEvap = self.RestEvap - self.ActEvapOpenWater self.RE = self.RestEvap self.ActEvapPond = self.ZeroMap if self.nrpaddyirri > 0: - self.ActEvapPond = min(self.PondingDepth,self.RestEvap) - self.PondingDepth = self.PondingDepth - self.ActEvapPond + self.ActEvapPond = min(self.PondingDepth, self.RestEvap) + self.PondingDepth = self.PondingDepth - self.ActEvapPond self.RestEvap = self.RestEvap - self.ActEvapPond - # Go from top to bottom layer self.zi_t = self.zi - for n in arange(0,len(self.UStoreLayerThickness)): + for n in arange(0, len(self.UStoreLayerThickness)): # Find layer with zi level - self.ZiLayer = ifthenelse(self.zi > self.SumThickness, min(self.ZeroMap + n,self.nrLayersMap-1), self.ZiLayer) + self.ZiLayer = ifthenelse( + self.zi > self.SumThickness, + min(self.ZeroMap + n, self.nrLayersMap - 1), + self.ZiLayer, + ) self.SumThickness = self.UStoreLayerThickness[n] + self.SumThickness - self.SaturationDeficit = self.SoilWaterCapacity - self.SatWaterDepth + # self.RestPotEvap, self.SatWaterDepth, self.ActEvapSat = actEvap_sat_SBM(self.ActRootingDepth, self.zi, self.SatWaterDepth, self.PotTrans, self.rootdistpar) - #self.RestPotEvap, self.SatWaterDepth, self.ActEvapSat = actEvap_sat_SBM(self.ActRootingDepth, self.zi, self.SatWaterDepth, self.PotTrans, self.rootdistpar) - - self.ActEvapUStore = self.ZeroMap - self.SumThickness = self.ZeroMap l_Thickness = [] self.storage = [] l_T = [] - for n in arange(0,len(self.UStoreLayerThickness)): + for n in arange(0, len(self.UStoreLayerThickness)): l_T.append(self.SumThickness) self.SumLayer = self.SumThickness self.SumThickness = self.UStoreLayerThickness[n] + self.SumThickness l_Thickness.append(self.SumThickness) # Height of unsat zone in layer n - self.L = ifthenelse(self.ZiLayer == n, ifthenelse(self.ZeroMap + n > 0, self.zi - l_Thickness[n-1], self.zi), self.UStoreLayerThickness[n] ) + self.L = ifthenelse( + self.ZiLayer == n, + ifthenelse(self.ZeroMap + n > 0, self.zi - l_Thickness[n - 1], self.zi), + self.UStoreLayerThickness[n], + ) # Depth for calculation of vertical fluxes (bottom layer or zi) - self.z = ifthenelse(self.ZiLayer == n, self.zi , self.SumThickness) - self.storage.append(self.L*(self.thetaS-self.thetaR)) + self.z = ifthenelse(self.ZiLayer == n, self.zi, self.SumThickness) + self.storage.append(self.L * (self.thetaS - self.thetaR)) # First layer is treated differently than layers below first layer if n == 0: - DownWard = InfiltSoilPath#MaxInfiltPath+MaxInfiltSoil + DownWard = InfiltSoilPath # MaxInfiltPath+MaxInfiltSoil self.UStoreLayerDepth[n] = self.UStoreLayerDepth[n] + DownWard - self.soilevap = soilevap_SBM(self.CanopyGapFraction,self.RestEvap,self.SoilWaterCapacity,self.SatWaterDepth,self.UStoreLayerDepth,self.zi,self.thetaS,self.thetaR,self.UStoreLayerThickness) - #assume soil evaporation is from first soil layer + self.soilevap = soilevap_SBM( + self.CanopyGapFraction, + self.RestEvap, + self.SoilWaterCapacity, + self.SatWaterDepth, + self.UStoreLayerDepth, + self.zi, + self.thetaS, + self.thetaR, + self.UStoreLayerThickness, + ) + # assume soil evaporation is from first soil layer if self.nrpaddyirri > 0: - self.soilevap = ifthenelse(self.PondingDepth > 0.0, 0.0,min(self.soilevap, self.UStoreLayerDepth[0])) + self.soilevap = ifthenelse( + self.PondingDepth > 0.0, + 0.0, + min(self.soilevap, self.UStoreLayerDepth[0]), + ) else: self.soilevap = min(self.soilevap, self.UStoreLayerDepth[n]) self.UStoreLayerDepth[n] = self.UStoreLayerDepth[n] - self.soilevap - self.PotTrans = self.PotTransSoil - self.soilevap - self.ActEvapOpenWater - self.RestPotEvap, self.SatWaterDepth, self.ActEvapSat = actEvap_sat_SBM(self.ActRootingDepth, self.zi, self.SatWaterDepth, self.PotTrans, self.rootdistpar) - self.UStoreLayerDepth[n], self.ActEvapUStore, self.RestPotEvap, self.ET = actEvap_unsat_SBM(self.ActRootingDepth, self.zi, self.UStoreLayerDepth[n], - self.ZiLayer, self.UStoreLayerThickness[n], - self.SumLayer, self.RestPotEvap, self.maskLayer[n], self.ZeroMap, self.ZeroMap+n, self.ActEvapUStore, self.c[n], self.L, self.thetaS, self.thetaR, self.UST) + self.PotTrans = ( + self.PotTransSoil - self.soilevap - self.ActEvapOpenWater + ) + self.RestPotEvap, self.SatWaterDepth, self.ActEvapSat = actEvap_sat_SBM( + self.ActRootingDepth, + self.zi, + self.SatWaterDepth, + self.PotTrans, + self.rootdistpar, + ) + self.UStoreLayerDepth[ + n + ], self.ActEvapUStore, self.RestPotEvap, self.ET = actEvap_unsat_SBM( + self.ActRootingDepth, + self.zi, + self.UStoreLayerDepth[n], + self.ZiLayer, + self.UStoreLayerThickness[n], + self.SumLayer, + self.RestPotEvap, + self.maskLayer[n], + self.ZeroMap, + self.ZeroMap + n, + self.ActEvapUStore, + self.c[n], + self.L, + self.thetaS, + self.thetaR, + self.UST, + ) - if len(self.UStoreLayerThickness) > 1: - st = self.KsatVerFrac[n]*self.KsatVer * exp(-self.f*self.z) * \ - min(((self.UStoreLayerDepth[n]/(self.L*(self.thetaS-self.thetaR)))**self.c[n]),1.0) - self.T[n] = ifthenelse(self.SaturationDeficit <= 0.00001, 0.0, min(self.UStoreLayerDepth[n],st)) - self.T[n] = ifthenelse(self.ZiLayer==n,self.maskLayer[n],self.T[n]) + st = ( + self.KsatVerFrac[n] + * self.KsatVer + * exp(-self.f * self.z) + * min( + ( + ( + self.UStoreLayerDepth[n] + / (self.L * (self.thetaS - self.thetaR)) + ) + ** self.c[n] + ), + 1.0, + ) + ) + self.T[n] = ifthenelse( + self.SaturationDeficit <= 0.00001, + 0.0, + min(self.UStoreLayerDepth[n], st), + ) + self.T[n] = ifthenelse( + self.ZiLayer == n, self.maskLayer[n], self.T[n] + ) self.UStoreLayerDepth[n] = self.UStoreLayerDepth[n] - self.T[n] else: - self.UStoreLayerDepth[n] = ifthenelse(self.ZiLayer 0 or hasattr(self,"IrriDemandExternal"): - if not hasattr(self,"IrriDemandExternal"): # if not given - self.IrriDemand, self.IrriDemandm3 = self.irrigationdemand(self.PotTrans,self.Transpiration,self.IrrigationAreas) - IRDemand = idtoid(self.IrrigationAreas, self.IrrigationSurfaceIntakes, self.IrriDemandm3) * -1.0 + if self.nrirri > 0 or hasattr(self, "IrriDemandExternal"): + if not hasattr(self, "IrriDemandExternal"): # if not given + self.IrriDemand, self.IrriDemandm3 = self.irrigationdemand( + self.PotTrans, self.Transpiration, self.IrrigationAreas + ) + IRDemand = ( + idtoid( + self.IrrigationAreas, + self.IrrigationSurfaceIntakes, + self.IrriDemandm3, + ) + * -1.0 + ) else: IRDemand = self.IrriDemandExternal # loop over irrigation areas and assign Q to linked river extraction points - self.Inflow = cover(IRDemand,self.Inflow) + self.Inflow = cover(IRDemand, self.Inflow) - ########################################################################## # Transfer of water from unsaturated to saturated store...################ ########################################################################## @@ -1468,240 +2205,422 @@ # Optional Macrco-Pore transfer (not yet implemented for # layers > 1) self.MporeTransfer = self.ActInfilt * self.MporeFrac self.SatWaterDepth = self.SatWaterDepth + self.MporeTransfer - #self.UStoreLayerDepth = self.UStoreLayerDepth - self.MporeTransfer + # self.UStoreLayerDepth = self.UStoreLayerDepth - self.MporeTransfer self.SaturationDeficit = self.SoilWaterCapacity - self.SatWaterDepth Ksat = self.ZeroMap - for n in arange(0,len(self.UStoreLayerThickness)): - Ksat = Ksat + ifthenelse(self.ZiLayer==n,self.KsatVerFrac[n]*self.KsatVer * exp(-self.f*self.zi),0.0) + for n in arange(0, len(self.UStoreLayerThickness)): + Ksat = Ksat + ifthenelse( + self.ZiLayer == n, + self.KsatVerFrac[n] * self.KsatVer * exp(-self.f * self.zi), + 0.0, + ) self.DeepKsat = self.KsatVer * exp(-self.f * self.SoilThickness) # now the actual transfer to the saturated store from layers with zi self.Transfer = self.ZeroMap - for n in arange(0,len(self.UStoreLayerThickness)): + for n in arange(0, len(self.UStoreLayerThickness)): if self.TransferMethod == 1: - self.L = ifthen(self.ZiLayer == n, ifthenelse(self.ZeroMap + n > 0, self.zi - l_Thickness[n-1], self.zi)) - self.Transfer = self.Transfer + ifthenelse(self.ZiLayer==n,min(cover(self.UStoreLayerDepth[n],0.0), - ifthenelse(self.SaturationDeficit <= 0.00001, 0.0,self.KsatVerFrac[n]*self.KsatVer * exp(-self.f*self.zi) * (min(cover(self.UStoreLayerDepth[n],0.0),(self.L+0.0001)*(self.thetaS-self.thetaR))) / (self.SaturationDeficit + 1))),0.0) + self.L = ifthen( + self.ZiLayer == n, + ifthenelse( + self.ZeroMap + n > 0, self.zi - l_Thickness[n - 1], self.zi + ), + ) + self.Transfer = self.Transfer + ifthenelse( + self.ZiLayer == n, + min( + cover(self.UStoreLayerDepth[n], 0.0), + ifthenelse( + self.SaturationDeficit <= 0.00001, + 0.0, + self.KsatVerFrac[n] + * self.KsatVer + * exp(-self.f * self.zi) + * ( + min( + cover(self.UStoreLayerDepth[n], 0.0), + (self.L + 0.0001) * (self.thetaS - self.thetaR), + ) + ) + / (self.SaturationDeficit + 1), + ), + ), + 0.0, + ) if self.TransferMethod == 2: - self.L = ifthen(self.ZiLayer == n, ifthenelse(self.ZeroMap + n > 0, self.zi - l_Thickness[n-1], self.zi)) - st = ifthen(self.ZiLayer==n, self.KsatVer * exp(-self.f*self.zi) * min((self.UStoreLayerDepth[n] /((self.L+0.0001)*(self.thetaS-self.thetaR))),1.0)**self.c[n]) - self.Transfer = self.Transfer + ifthenelse(self.ZiLayer==n, min(self.UStoreLayerDepth[n], - ifthenelse(self.SaturationDeficit <= 0.00001, 0.0, st)),0.0) + self.L = ifthen( + self.ZiLayer == n, + ifthenelse( + self.ZeroMap + n > 0, self.zi - l_Thickness[n - 1], self.zi + ), + ) + st = ifthen( + self.ZiLayer == n, + self.KsatVer + * exp(-self.f * self.zi) + * min( + ( + self.UStoreLayerDepth[n] + / ((self.L + 0.0001) * (self.thetaS - self.thetaR)) + ), + 1.0, + ) + ** self.c[n], + ) + self.Transfer = self.Transfer + ifthenelse( + self.ZiLayer == n, + min( + self.UStoreLayerDepth[n], + ifthenelse(self.SaturationDeficit <= 0.00001, 0.0, st), + ), + 0.0, + ) - #check soil moisture + # check soil moisture self.ToExtra = self.ZeroMap - for n in arange(len(self.UStoreLayerThickness)-1, -1, -1): - #self.UStoreLayerDepth[n] = ifthenelse(self.ZiLayer<=n, self.UStoreLayerDepth[n] + self.ToExtra,self.UStoreLayerDepth[n]) - diff = ifthenelse(self.ZiLayer == n, max(0.0,(cover(self.UStoreLayerDepth[n],0.0) - self.Transfer)-self.storage[n]), max(self.ZeroMap,cover(self.UStoreLayerDepth[n],0.0) - \ - ifthenelse(self.zi <= l_T[n],0.0, self.storage[n]))) - self.ToExtra = ifthenelse(diff>0,diff,self.ZeroMap) - self.UStoreLayerDepth[n] = self.UStoreLayerDepth[n]-diff + for n in arange(len(self.UStoreLayerThickness) - 1, -1, -1): + # self.UStoreLayerDepth[n] = ifthenelse(self.ZiLayer<=n, self.UStoreLayerDepth[n] + self.ToExtra,self.UStoreLayerDepth[n]) + diff = ifthenelse( + self.ZiLayer == n, + max( + 0.0, + (cover(self.UStoreLayerDepth[n], 0.0) - self.Transfer) + - self.storage[n], + ), + max( + self.ZeroMap, + cover(self.UStoreLayerDepth[n], 0.0) + - ifthenelse(self.zi <= l_T[n], 0.0, self.storage[n]), + ), + ) + self.ToExtra = ifthenelse(diff > 0, diff, self.ZeroMap) + self.UStoreLayerDepth[n] = self.UStoreLayerDepth[n] - diff - if n>0: - self.UStoreLayerDepth[n-1] = self.UStoreLayerDepth[n-1] + self.ToExtra + if n > 0: + self.UStoreLayerDepth[n - 1] = ( + self.UStoreLayerDepth[n - 1] + self.ToExtra + ) - #self.UStoreLayerDepth[n] = ifthenelse(self.ZiLayer<=n, self.UStoreLayerDepth[n]-diff,self.UStoreLayerDepth[n]) + # self.UStoreLayerDepth[n] = ifthenelse(self.ZiLayer<=n, self.UStoreLayerDepth[n]-diff,self.UStoreLayerDepth[n]) SatFlow = self.ToExtra - UStoreCapacity = self.SoilWaterCapacity - self.SatWaterDepth - sum_list_cover(self.UStoreLayerDepth,self.ZeroMap) + UStoreCapacity = ( + self.SoilWaterCapacity + - self.SatWaterDepth + - sum_list_cover(self.UStoreLayerDepth, self.ZeroMap) + ) - MaxCapFlux = max(0.0, min(Ksat, self.ActEvapUStore, UStoreCapacity, self.SatWaterDepth)) + MaxCapFlux = max( + 0.0, min(Ksat, self.ActEvapUStore, UStoreCapacity, self.SatWaterDepth) + ) - # No capilary flux is roots are in water, max flux if very near to water, lower flux if distance is large - CapFluxScale = ifthenelse(self.zi > self.ActRootingDepth, - self.CapScale / (self.CapScale + self.zi - self.ActRootingDepth) * self.timestepsecs/self.basetimestep, 0.0) + CapFluxScale = ifthenelse( + self.zi > self.ActRootingDepth, + self.CapScale + / (self.CapScale + self.zi - self.ActRootingDepth) + * self.timestepsecs + / self.basetimestep, + 0.0, + ) self.CapFlux = MaxCapFlux * CapFluxScale ToAdd = self.CapFlux sumLayer = self.ZeroMap - #Now add capflux to the layers one by one (from bottom to top) - for n in arange(len(self.UStoreLayerThickness)-1, -1, -1): + # Now add capflux to the layers one by one (from bottom to top) + for n in arange(len(self.UStoreLayerThickness) - 1, -1, -1): - L = ifthenelse(self.ZiLayer == n, ifthenelse(self.ZeroMap + n > 0, self.zi - l_Thickness[n-1], self.zi), self.UStoreLayerThickness[n] ) - thisLayer = ifthenelse(self.ZiLayer <= n,min(ToAdd,max(L*(self.thetaS-self.thetaR)-self.UStoreLayerDepth[n],0.0)), 0.0) - self.UStoreLayerDepth[n] = ifthenelse(self.ZiLayer <= n, self.UStoreLayerDepth[n] + thisLayer,self.UStoreLayerDepth[n] ) - ToAdd = ToAdd - cover(thisLayer,0.0) - sumLayer = sumLayer + cover(thisLayer,0.0) + L = ifthenelse( + self.ZiLayer == n, + ifthenelse(self.ZeroMap + n > 0, self.zi - l_Thickness[n - 1], self.zi), + self.UStoreLayerThickness[n], + ) + thisLayer = ifthenelse( + self.ZiLayer <= n, + min( + ToAdd, + max( + L * (self.thetaS - self.thetaR) - self.UStoreLayerDepth[n], 0.0 + ), + ), + 0.0, + ) + self.UStoreLayerDepth[n] = ifthenelse( + self.ZiLayer <= n, + self.UStoreLayerDepth[n] + thisLayer, + self.UStoreLayerDepth[n], + ) + ToAdd = ToAdd - cover(thisLayer, 0.0) + sumLayer = sumLayer + cover(thisLayer, 0.0) # Determine Ksat at base - self.DeepTransfer = min(self.SatWaterDepth,self.DeepKsat) - #ActLeakage = 0.0 + self.DeepTransfer = min(self.SatWaterDepth, self.DeepKsat) + # ActLeakage = 0.0 # Now add leakage. to deeper groundwater - self.ActLeakage = cover(max(0.0,min(self.MaxLeakage,self.DeepTransfer)),0) - self.Percolation = cover(max(0.0,min(self.MaxPercolation,self.DeepTransfer)),0) + self.ActLeakage = cover(max(0.0, min(self.MaxLeakage, self.DeepTransfer)), 0) + self.Percolation = cover( + max(0.0, min(self.MaxPercolation, self.DeepTransfer)), 0 + ) + # self.ActLeakage = ifthenelse(self.Seepage > 0.0, -1.0 * self.Seepage, self.ActLeakage) + self.SatWaterDepth = ( + self.SatWaterDepth + + self.Transfer + - sumLayer + - self.ActLeakage + - self.Percolation + ) - #self.ActLeakage = ifthenelse(self.Seepage > 0.0, -1.0 * self.Seepage, self.ActLeakage) - self.SatWaterDepth = self.SatWaterDepth + self.Transfer - sumLayer - self.ActLeakage - self.Percolation + for n in arange(0, len(self.UStoreLayerThickness)): + self.UStoreLayerDepth[n] = ifthenelse( + self.ZiLayer == n, + self.UStoreLayerDepth[n] - self.Transfer, + self.UStoreLayerDepth[n], + ) - - for n in arange(0,len(self.UStoreLayerThickness)): - self.UStoreLayerDepth[n] = ifthenelse(self.ZiLayer==n,self.UStoreLayerDepth[n] - self.Transfer, self.UStoreLayerDepth[n]) - # Determine % saturated taking into account subcell fraction - self.Sat = max(self.SubCellFrac, scalar(self.SatWaterDepth >= (self.SoilWaterCapacity * 0.999))) + self.Sat = max( + self.SubCellFrac, + scalar(self.SatWaterDepth >= (self.SoilWaterCapacity * 0.999)), + ) ########################################################################## # Horizontal (downstream) transport of water ############################# ########################################################################## - self.zi = max(0.0, self.SoilThickness - self.SatWaterDepth / ( - self.thetaS - self.thetaR)) # Determine actual water depth + self.zi = max( + 0.0, self.SoilThickness - self.SatWaterDepth / (self.thetaS - self.thetaR) + ) # Determine actual water depth - # Re-Determine saturation deficit. NB, as noted by Vertessy and Elsenbeer 1997 # this deficit does NOT take into account the water in the unsaturated zone self.SaturationDeficit = self.SoilWaterCapacity - self.SatWaterDepth - #self.logger.debug("Waterdem set to Altitude....") + # self.logger.debug("Waterdem set to Altitude....") self.WaterDem = self.Altitude - (self.zi * 0.001) - self.waterSlope = max(0.000001, slope(self.WaterDem) * celllength() / self.reallength) + self.waterSlope = max( + 0.000001, slope(self.WaterDem) * celllength() / self.reallength + ) if self.waterdem: self.waterLdd = lddcreate(self.WaterDem, 1E35, 1E35, 1E35, 1E35) - #waterLdd = lddcreate(waterDem,1,1,1,1) + # waterLdd = lddcreate(waterDem,1,1,1,1) + # TODO: We should make a couple of itterations here... - #TODO: We should make a couple of itterations here... - if self.waterdem: if self.LateralMethod == 1: - Lateral = self.KsatHorFrac * self.KsatVer * self.waterSlope * exp(-self.SaturationDeficit / self.M) + Lateral = ( + self.KsatHorFrac + * self.KsatVer + * self.waterSlope + * exp(-self.SaturationDeficit / self.M) + ) elif self.LateralMethod == 2: - #Lateral = Ksat * self.waterSlope - Lateral = self.KsatHorFrac * self.KsatVer * (exp(-self.f*self.zi)-exp(-self.f*self.SoilThickness))*(1/self.f)/(self.SoilThickness-self.zi)*self.waterSlope + # Lateral = Ksat * self.waterSlope + Lateral = ( + self.KsatHorFrac + * self.KsatVer + * (exp(-self.f * self.zi) - exp(-self.f * self.SoilThickness)) + * (1 / self.f) + / (self.SoilThickness - self.zi) + * self.waterSlope + ) MaxHor = max(0.0, min(Lateral, self.SatWaterDepth)) - self.SatWaterFlux = accucapacityflux(self.waterLdd, self.SatWaterDepth, MaxHor) - self.SatWaterDepth = accucapacitystate(self.waterLdd, self.SatWaterDepth, MaxHor) + self.SatWaterFlux = accucapacityflux( + self.waterLdd, self.SatWaterDepth, MaxHor + ) + self.SatWaterDepth = accucapacitystate( + self.waterLdd, self.SatWaterDepth, MaxHor + ) else: # - #MaxHor = max(0,min(self.KsatVer * self.Slope * exp(-SaturationDeficit/self.M),self.SatWaterDepth*(self.thetaS-self.thetaR))) * timestepsecs/basetimestep - #MaxHor = max(0.0, min(self.KsatVer * self.Slope * exp(-selield' object does not support item assignmentf.SaturationDeficit / self.M), + # MaxHor = max(0,min(self.KsatVer * self.Slope * exp(-SaturationDeficit/self.M),self.SatWaterDepth*(self.thetaS-self.thetaR))) * timestepsecs/basetimestep + # MaxHor = max(0.0, min(self.KsatVer * self.Slope * exp(-selield' object does not support item assignmentf.SaturationDeficit / self.M), # self.SatWaterDepth)) if self.LateralMethod == 1: - Lateral = self.KsatHorFrac * self.KsatVer * self.waterSlope * exp(-self.SaturationDeficit / self.M) + Lateral = ( + self.KsatHorFrac + * self.KsatVer + * self.waterSlope + * exp(-self.SaturationDeficit / self.M) + ) elif self.LateralMethod == 2: - #Lateral = Ksat * self.waterSlope - Lateral = self.KsatHorFrac * self.KsatVer * (exp(-self.f*self.zi)-exp(-self.f*self.SoilThickness))*(1/self.f)/(self.SoilThickness-self.zi+1.0)*self.waterSlope + # Lateral = Ksat * self.waterSlope + Lateral = ( + self.KsatHorFrac + * self.KsatVer + * (exp(-self.f * self.zi) - exp(-self.f * self.SoilThickness)) + * (1 / self.f) + / (self.SoilThickness - self.zi + 1.0) + * self.waterSlope + ) - MaxHor = max(0.0, min(Lateral, self.SatWaterDepth)) - #MaxHor = self.ZeroMap - self.SatWaterFlux = accucapacityflux(self.TopoLdd, self.SatWaterDepth, MaxHor) - self.SatWaterDepth = accucapacitystate(self.TopoLdd, self.SatWaterDepth, MaxHor) + # MaxHor = self.ZeroMap + self.SatWaterFlux = accucapacityflux( + self.TopoLdd, self.SatWaterDepth, MaxHor + ) + self.SatWaterDepth = accucapacitystate( + self.TopoLdd, self.SatWaterDepth, MaxHor + ) ########################################################################## # Determine returnflow from first zone ########################## ########################################################################## - self.ExfiltWaterFrac = sCurve(self.SatWaterDepth, a=self.SoilWaterCapacity, c=5.0) - self.ExfiltWater = self.ExfiltWaterFrac * (self.SatWaterDepth - self.SoilWaterCapacity) - #self.ExfiltWater=ifthenelse (self.SatWaterDepth - self.SoilWaterCapacity > 0 , self.SatWaterDepth - self.SoilWaterCapacity , 0.0) + self.ExfiltWaterFrac = sCurve( + self.SatWaterDepth, a=self.SoilWaterCapacity, c=5.0 + ) + self.ExfiltWater = self.ExfiltWaterFrac * ( + self.SatWaterDepth - self.SoilWaterCapacity + ) + # self.ExfiltWater=ifthenelse (self.SatWaterDepth - self.SoilWaterCapacity > 0 , self.SatWaterDepth - self.SoilWaterCapacity , 0.0) self.SatWaterDepth = self.SatWaterDepth - self.ExfiltWater # Re-determine UStoreCapacity - self.zi = max(0.0, self.SoilThickness - self.SatWaterDepth / ( - self.thetaS - self.thetaR)) # Determine actual water depth + self.zi = max( + 0.0, self.SoilThickness - self.SatWaterDepth / (self.thetaS - self.thetaR) + ) # Determine actual water depth self.SumThickness = self.ZeroMap self.ZiLayer = self.ZeroMap - for n in arange(0,len(self.UStoreLayerThickness)): + for n in arange(0, len(self.UStoreLayerThickness)): # Find layer with zi level - self.ZiLayer = ifthenelse(self.zi > self.SumThickness, min(self.ZeroMap + n,self.nrLayersMap-1), self.ZiLayer) + self.ZiLayer = ifthenelse( + self.zi > self.SumThickness, + min(self.ZeroMap + n, self.nrLayersMap - 1), + self.ZiLayer, + ) self.SumThickness = self.UStoreLayerThickness[n] + self.SumThickness self.SumThickness = self.ZeroMap l_Thickness = [] self.storage = [] - self.L =[] + self.L = [] l_T = [] - #redistribute soil moisture (balance) + # redistribute soil moisture (balance) for n in arange(len(self.UStoreLayerThickness)): self.SumLayer = self.SumThickness l_T.append(self.SumThickness) self.SumThickness = self.UStoreLayerThickness[n] + self.SumThickness l_Thickness.append(self.SumThickness) # Height of unsat zone in layer n - self.L.append(ifthenelse(self.ZiLayer == n, ifthenelse(self.ZeroMap + n > 0, self.zi - l_Thickness[n-1], self.zi), self.UStoreLayerThickness[n] )) + self.L.append( + ifthenelse( + self.ZiLayer == n, + ifthenelse( + self.ZeroMap + n > 0, self.zi - l_Thickness[n - 1], self.zi + ), + self.UStoreLayerThickness[n], + ) + ) - self.storage.append(self.L[n]*(self.thetaS-self.thetaR)) + self.storage.append(self.L[n] * (self.thetaS - self.thetaR)) - - self.ExfiltFromUstore = self.ZeroMap + for n in arange(len(self.UStoreLayerThickness) - 1, -1, -1): + diff = max( + self.ZeroMap, + cover(self.UStoreLayerDepth[n], 0.0) + - ifthenelse(self.zi <= l_T[n], 0.0, self.storage[n]), + ) + self.ExfiltFromUstore = ifthenelse(diff > 0, diff, self.ZeroMap) + self.UStoreLayerDepth[n] = self.UStoreLayerDepth[n] - diff - for n in arange(len(self.UStoreLayerThickness)-1, -1, -1): - diff = max(self.ZeroMap,cover(self.UStoreLayerDepth[n],0.0) - ifthenelse(self.zi <= l_T[n],0.0, self.storage[n])) - self.ExfiltFromUstore = ifthenelse(diff>0,diff,self.ZeroMap) - self.UStoreLayerDepth[n] = self.UStoreLayerDepth[n]-diff + if n > 0: + self.UStoreLayerDepth[n - 1] = ( + self.UStoreLayerDepth[n - 1] + self.ExfiltFromUstore + ) - if n>0: - self.UStoreLayerDepth[n-1] = self.UStoreLayerDepth[n-1] + self.ExfiltFromUstore - # Re-determine UStoreCapacityield' object does not support item assignment - UStoreCapacity = self.SoilWaterCapacity - self.SatWaterDepth - sum_list_cover(self.UStoreLayerDepth,self.ZeroMap) + UStoreCapacity = ( + self.SoilWaterCapacity + - self.SatWaterDepth + - sum_list_cover(self.UStoreLayerDepth, self.ZeroMap) + ) + # self.AvailableForInfiltration = self.AvailableForInfiltration - InfiltSoilPath - SatFlow #MaxInfiltPath+MaxInfiltSoil + SatFlow - #self.AvailableForInfiltration = self.AvailableForInfiltration - InfiltSoilPath - SatFlow #MaxInfiltPath+MaxInfiltSoil + SatFlow + self.ActInfilt = ( + InfiltSoilPath - SatFlow + ) # MaxInfiltPath+MaxInfiltSoil - SatFlow - self.ActInfilt = InfiltSoilPath - SatFlow#MaxInfiltPath+MaxInfiltSoil - SatFlow + self.InfiltExcess = self.AvailableForInfiltration - InfiltSoilPath + SatFlow - self.InfiltExcess = self.AvailableForInfiltration - InfiltSoilPath + SatFlow + self.ExcessWater = self.AvailableForInfiltration - InfiltSoilPath + SatFlow - self.ExcessWater = self.AvailableForInfiltration - InfiltSoilPath + SatFlow - self.CumInfiltExcess = self.CumInfiltExcess + self.InfiltExcess + # self.ExfiltFromUstore = ifthenelse(self.zi == 0.0,self.ExfiltFromUstore,self.ZeroMap) - - #self.ExfiltFromUstore = ifthenelse(self.zi == 0.0,self.ExfiltFromUstore,self.ZeroMap) - self.ExfiltWater = self.ExfiltWater + self.ExfiltFromUstore self.inund = self.ExfiltWater + self.ExcessWater ponding_add = self.ZeroMap if self.nrpaddyirri > 0: - ponding_add = cover(min(ifthen(self.h_p > 0,self.inund),self.h_p-self.PondingDepth),0.0) + ponding_add = cover( + min(ifthen(self.h_p > 0, self.inund), self.h_p - self.PondingDepth), 0.0 + ) self.PondingDepth = self.PondingDepth + ponding_add - irr_depth = ifthenelse(self.PondingDepth < self.h_min, self.h_max - self.PondingDepth, 0.0) * self.CRPST - sqmarea = areatotal(self.reallength * self.reallength, self.IrrigationPaddyAreas) - self.IrriDemandm3 = cover((irr_depth/1000.0)*sqmarea,0) - IRDemand = idtoid(self.IrrigationPaddyAreas, self.IrrigationSurfaceIntakes, self.IrriDemandm3) * (-1.0 / self.timestepsecs) + irr_depth = ( + ifthenelse( + self.PondingDepth < self.h_min, self.h_max - self.PondingDepth, 0.0 + ) + * self.CRPST + ) + sqmarea = areatotal( + self.reallength * self.reallength, self.IrrigationPaddyAreas + ) + self.IrriDemandm3 = cover((irr_depth / 1000.0) * sqmarea, 0) + IRDemand = idtoid( + self.IrrigationPaddyAreas, + self.IrrigationSurfaceIntakes, + self.IrriDemandm3, + ) * (-1.0 / self.timestepsecs) - self.IRDemand= IRDemand - self.Inflow = cover(IRDemand,self.Inflow) + self.IRDemand = IRDemand + self.Inflow = cover(IRDemand, self.Inflow) self.irr_depth = irr_depth + UStoreCapacity = ( + self.SoilWaterCapacity + - self.SatWaterDepth + - sum_list_cover(self.UStoreLayerDepth, self.ZeroMap) + ) + self.UStoreDepth = sum_list_cover(self.UStoreLayerDepth, self.ZeroMap) - UStoreCapacity = self.SoilWaterCapacity - self.SatWaterDepth - sum_list_cover(self.UStoreLayerDepth,self.ZeroMap) - self.UStoreDepth = sum_list_cover(self.UStoreLayerDepth,self.ZeroMap) + Ksat = self.KsatVer * exp(-self.f * self.zi) - Ksat = self.KsatVer * exp(-self.f*self.zi) - - - SurfaceWater = self.WaterLevel/1000.0 # SurfaceWater (mm) + SurfaceWater = self.WaterLevel / 1000.0 # SurfaceWater (mm) self.CumSurfaceWater = self.CumSurfaceWater + SurfaceWater - # Estimate water that may re-infiltrate # - Never more that 70% of the available water # - self.MaxReinFilt: a map with reinfilt locations (usually the river mak) can be supplied) # - take into account that the river may not cover the whole cell if self.reInfilt: - self.reinfiltwater = min(self.MaxReinfilt,max(0, min(SurfaceWater * self.RiverWidth/self.reallength * 0.7, - min(self.InfiltCapSoil * (1.0 - self.PathFrac), UStoreCapacity)))) + self.reinfiltwater = min( + self.MaxReinfilt, + max( + 0, + min( + SurfaceWater * self.RiverWidth / self.reallength * 0.7, + min(self.InfiltCapSoil * (1.0 - self.PathFrac), UStoreCapacity), + ), + ), + ) self.CumReinfilt = self.CumReinfilt + self.reinfiltwater # TODO: This still has to be reworked fro the differnt layers self.UStoreDepth = self.UStoreDepth + self.reinfiltwater @@ -1710,63 +2629,116 @@ # The Max here may lead to watbal error. However, if inwaterMMM becomes < 0, the kinematic wave becomes very slow...... if self.reInfilt: - self.InwaterMM = self.ExfiltWater + self.ExcessWater + self.SubCellRunoff + \ - self.SubCellGWRunoff + self.RunoffOpenWater - \ - self.reinfiltwater - self.ActEvapOpenWater - ponding_add + self.InwaterMM = ( + self.ExfiltWater + + self.ExcessWater + + self.SubCellRunoff + + self.SubCellGWRunoff + + self.RunoffOpenWater + - self.reinfiltwater + - self.ActEvapOpenWater + - ponding_add + ) else: - self.InwaterMM = max(0.0,self.ExfiltWater + self.ExcessWater + self.SubCellRunoff + \ - self.SubCellGWRunoff + self.RunoffOpenWater - \ - self.reinfiltwater - self.ActEvapOpenWater - ponding_add) + self.InwaterMM = max( + 0.0, + self.ExfiltWater + + self.ExcessWater + + self.SubCellRunoff + + self.SubCellGWRunoff + + self.RunoffOpenWater + - self.reinfiltwater + - self.ActEvapOpenWater + - ponding_add, + ) self.Inwater = self.InwaterMM * self.ToCubic # m3/s - - #only run the reservoir module if needed + # only run the reservoir module if needed if self.nrresSimple > 0: - self.ReservoirVolume, self.OutflowSR, self.ResPercFull, self.ResPrecipSR, self.ResEvapSR,\ - self.DemandRelease = simplereservoir(self.ReservoirVolume, self.SurfaceRunoff,self.ResSimpleArea,\ - self.ResMaxVolume, self.ResTargetFullFrac, - self.ResMaxRelease, self.ResDemand, - self.ResTargetMinFrac, self.ReserVoirSimpleLocs, - self.ReserVoirPrecip, self.ReserVoirPotEvap, - self.ReservoirSimpleAreas, timestepsecs=self.timestepsecs) - self.OutflowDwn = upstream(self.TopoLddOrg, cover(self.OutflowSR,scalar(0.0))) - self.Inflow = self.OutflowDwn + cover(self.Inflow,self.ZeroMap) + self.ReservoirVolume, self.OutflowSR, self.ResPercFull, self.ResPrecipSR, self.ResEvapSR, self.DemandRelease = simplereservoir( + self.ReservoirVolume, + self.SurfaceRunoff, + self.ResSimpleArea, + self.ResMaxVolume, + self.ResTargetFullFrac, + self.ResMaxRelease, + self.ResDemand, + self.ResTargetMinFrac, + self.ReserVoirSimpleLocs, + self.ReserVoirPrecip, + self.ReserVoirPotEvap, + self.ReservoirSimpleAreas, + timestepsecs=self.timestepsecs, + ) + self.OutflowDwn = upstream( + self.TopoLddOrg, cover(self.OutflowSR, scalar(0.0)) + ) + self.Inflow = self.OutflowDwn + cover(self.Inflow, self.ZeroMap) elif self.nrresComplex > 0: - self.ReservoirWaterLevel, self.OutflowCR, self.ResPrecipCR, self.ResEvapCR,\ - self.ReservoirVolumeCR = complexreservoir(self.ReservoirWaterLevel, self.ReserVoirComplexLocs, self.LinkedReservoirLocs, self.ResArea,\ - self.ResThreshold, self.ResStorFunc, self.ResOutflowFunc, self.sh, self.hq, self.Res_b, - self.Res_e, self.SurfaceRunoff,self.ReserVoirPrecip, self.ReserVoirPotEvap, - self.ReservoirComplexAreas, self.wf_supplyJulianDOY(), timestepsecs=self.timestepsecs) - self.OutflowDwn = upstream(self.TopoLddOrg, cover(self.OutflowCR,scalar(0.0))) - self.Inflow = self.OutflowDwn + cover(self.Inflow,self.ZeroMap) + self.ReservoirWaterLevel, self.OutflowCR, self.ResPrecipCR, self.ResEvapCR, self.ReservoirVolumeCR = complexreservoir( + self.ReservoirWaterLevel, + self.ReserVoirComplexLocs, + self.LinkedReservoirLocs, + self.ResArea, + self.ResThreshold, + self.ResStorFunc, + self.ResOutflowFunc, + self.sh, + self.hq, + self.Res_b, + self.Res_e, + self.SurfaceRunoff, + self.ReserVoirPrecip, + self.ReserVoirPotEvap, + self.ReservoirComplexAreas, + self.wf_supplyJulianDOY(), + timestepsecs=self.timestepsecs, + ) + self.OutflowDwn = upstream( + self.TopoLddOrg, cover(self.OutflowCR, scalar(0.0)) + ) + self.Inflow = self.OutflowDwn + cover(self.Inflow, self.ZeroMap) else: - self.Inflow= cover(self.Inflow,self.ZeroMap) - + self.Inflow = cover(self.Inflow, self.ZeroMap) + self.ExfiltWaterCubic = self.ExfiltWater * self.ToCubic self.SubCellGWRunoffCubic = self.SubCellGWRunoff * self.ToCubic self.SubCellRunoffCubic = self.SubCellRunoff * self.ToCubic self.InfiltExcessCubic = self.InfiltExcess * self.ToCubic self.ReinfiltCubic = -1.0 * self.reinfiltwater * self.ToCubic - #self.Inwater = self.Inwater + self.Inflow # Add abstractions/inflows in m^3/sec + # self.Inwater = self.Inwater + self.Inflow # Add abstractions/inflows in m^3/sec # Check if we do not try to abstract more runoff then present - self.InflowKinWaveCell = upstream(self.TopoLdd, self.SurfaceRunoff) #NG The extraction should be equal to the discharge upstream cell. You should not make the abstraction depended on the downstream cell, because they are correlated. During a stationary sum they will get equal to each other. - MaxExtract = self.InflowKinWaveCell + self.Inwater #NG + self.InflowKinWaveCell = upstream( + self.TopoLdd, self.SurfaceRunoff + ) # NG The extraction should be equal to the discharge upstream cell. You should not make the abstraction depended on the downstream cell, because they are correlated. During a stationary sum they will get equal to each other. + MaxExtract = self.InflowKinWaveCell + self.Inwater # NG # MaxExtract = self.SurfaceRunoff + self.Inwater - self.SurfaceWaterSupply = ifthenelse (self.Inflow < 0.0 , min(MaxExtract,-1.0 * self.Inflow), self.ZeroMap) - self.OldSurfaceRunoff=self.SurfaceRunoff #NG Store for iteration - self.OldInwater=self.Inwater - self.Inwater = self.Inwater + ifthenelse(self.SurfaceWaterSupply> 0, -1.0 * self.SurfaceWaterSupply,self.Inflow) + self.SurfaceWaterSupply = ifthenelse( + self.Inflow < 0.0, min(MaxExtract, -1.0 * self.Inflow), self.ZeroMap + ) + self.OldSurfaceRunoff = self.SurfaceRunoff # NG Store for iteration + self.OldInwater = self.Inwater + self.Inwater = self.Inwater + ifthenelse( + self.SurfaceWaterSupply > 0, -1.0 * self.SurfaceWaterSupply, self.Inflow + ) - ########################################################################## # Runoff calculation via Kinematic wave ################################## ########################################################################## # per distance along stream q = self.Inwater / self.DCL # discharge (m3/s) - self.SurfaceRunoff = kinematic(self.TopoLdd, self.SurfaceRunoff, q, self.Alpha, self.Beta, self.Tslice, - self.timestepsecs, self.DCL) # m3/s + self.SurfaceRunoff = kinematic( + self.TopoLdd, + self.SurfaceRunoff, + q, + self.Alpha, + self.Beta, + self.Tslice, + self.timestepsecs, + self.DCL, + ) # m3/s # If inflow is negative we have abstractions. Check if demand can be met (by looking # at the flow in the upstream cell) and iterate if needed @@ -1782,20 +2754,44 @@ # (Runoff calculation via Kinematic wave) ################################ ########################################################################## MaxExtract = self.InflowKinWaveCell + self.OldInwater - self.SurfaceWaterSupply = ifthenelse (self.Inflow < 0.0 , min(MaxExtract,-1.0 * self.Inflow),\ - self.ZeroMap) + self.SurfaceWaterSupply = ifthenelse( + self.Inflow < 0.0, min(MaxExtract, -1.0 * self.Inflow), self.ZeroMap + ) # Fraction of demand that is not used but flows back into the river get fracttion and move to return locations - self.DemandReturnFlow = cover(idtoid(self.IrrigationSurfaceIntakes,self.IrrigationSurfaceReturn, - self.DemandReturnFlowFraction * self.SurfaceWaterSupply),0.0) + self.DemandReturnFlow = cover( + idtoid( + self.IrrigationSurfaceIntakes, + self.IrrigationSurfaceReturn, + self.DemandReturnFlowFraction * self.SurfaceWaterSupply, + ), + 0.0, + ) - self.Inwater = self.OldInwater + ifthenelse(self.SurfaceWaterSupply> 0, -1.0 * self.SurfaceWaterSupply,\ - self.Inflow) + self.DemandReturnFlow + self.Inwater = ( + self.OldInwater + + ifthenelse( + self.SurfaceWaterSupply > 0, + -1.0 * self.SurfaceWaterSupply, + self.Inflow, + ) + + self.DemandReturnFlow + ) # per distance along stream q = self.Inwater / self.DCL # discharge (m3/s) - self.SurfaceRunoff = kinematic(self.TopoLdd, self.OldSurfaceRunoff, q, self.Alpha, self.Beta, self.Tslice, - self.timestepsecs, self.DCL) # m3/s - self.SurfaceRunoffMM = self.SurfaceRunoff * self.QMMConv # SurfaceRunoffMM (mm) from SurfaceRunoff (m3/s) + self.SurfaceRunoff = kinematic( + self.TopoLdd, + self.OldSurfaceRunoff, + q, + self.Alpha, + self.Beta, + self.Tslice, + self.timestepsecs, + self.DCL, + ) # m3/s + self.SurfaceRunoffMM = ( + self.SurfaceRunoff * self.QMMConv + ) # SurfaceRunoffMM (mm) from SurfaceRunoff (m3/s) self.InflowKinWaveCell = upstream(self.TopoLdd, self.OldSurfaceRunoff) deltasup = float(mapmaximum(abs(oldsup - self.SurfaceWaterSupply))) @@ -1806,30 +2802,51 @@ self.InflowKinWaveCell = upstream(self.TopoLdd, self.SurfaceRunoff) self.updateRunOff() else: - self.SurfaceRunoffMM = self.SurfaceRunoff * self.QMMConv # SurfaceRunoffMM (mm) from SurfaceRunoff (m3/s) + self.SurfaceRunoffMM = ( + self.SurfaceRunoff * self.QMMConv + ) # SurfaceRunoffMM (mm) from SurfaceRunoff (m3/s) self.updateRunOff() # Now add the supply that is linked to irrigation areas to extra precip if self.nrirri > 0: # loop over irrigation areas and spread-out the supply over the area - IRSupplymm = idtoid(self.IrrigationSurfaceIntakes, self.IrrigationAreas, - self.SurfaceWaterSupply * (1 - self.DemandReturnFlowFraction)) - sqmarea = areatotal(self.reallength * self.reallength, nominal(self.IrrigationAreas)) + IRSupplymm = idtoid( + self.IrrigationSurfaceIntakes, + self.IrrigationAreas, + self.SurfaceWaterSupply * (1 - self.DemandReturnFlowFraction), + ) + sqmarea = areatotal( + self.reallength * self.reallength, nominal(self.IrrigationAreas) + ) - self.IRSupplymm = cover(IRSupplymm/ (sqmarea / 1000.0 / self.timestepsecs),0.0) + self.IRSupplymm = cover( + IRSupplymm / (sqmarea / 1000.0 / self.timestepsecs), 0.0 + ) if self.nrpaddyirri > 0: # loop over irrigation areas and spread-out the supply over the area - IRSupplymm = idtoid(self.IrrigationSurfaceIntakes, ifthen(self.IrriDemandm3 > 0,self.IrrigationPaddyAreas), self.SurfaceWaterSupply) - sqmarea = areatotal(self.reallength * self.reallength, nominal(ifthen(self.IrriDemandm3 > 0,self.IrrigationPaddyAreas))) + IRSupplymm = idtoid( + self.IrrigationSurfaceIntakes, + ifthen(self.IrriDemandm3 > 0, self.IrrigationPaddyAreas), + self.SurfaceWaterSupply, + ) + sqmarea = areatotal( + self.reallength * self.reallength, + nominal(ifthen(self.IrriDemandm3 > 0, self.IrrigationPaddyAreas)), + ) - self.IRSupplymm = cover(((IRSupplymm * self.timestepsecs * 1000) / sqmarea),0.0) + self.IRSupplymm = cover( + ((IRSupplymm * self.timestepsecs * 1000) / sqmarea), 0.0 + ) + self.MassBalKinWave = ( + (-self.KinWaveVolume + self.OldKinWaveVolume) / self.timestepsecs + + self.InflowKinWaveCell + + self.Inwater + - self.SurfaceRunoff + ) - self.MassBalKinWave = (-self.KinWaveVolume + self.OldKinWaveVolume) / self.timestepsecs +\ - self.InflowKinWaveCell + self.Inwater - self.SurfaceRunoff - Runoff = self.SurfaceRunoff # Updating @@ -1845,12 +2862,20 @@ # No determine multiplication ratio for each gauge influence area. # For missing gauges 1.0 is assumed (no change). # UpDiff = areamaximum(QM, self.UpdateMap) - areamaximum(self.SurfaceRunoffMM, self.UpdateMap) - UpRatio = areamaximum(self.QM, self.UpdateMap) / areamaximum(self.SurfaceRunoffMM, self.UpdateMap) + UpRatio = areamaximum(self.QM, self.UpdateMap) / areamaximum( + self.SurfaceRunoffMM, self.UpdateMap + ) UpRatio = cover(areaaverage(UpRatio, self.TopoId), 1.0) # Now split between Soil and Kyn wave - self.UpRatioKyn = min(self.MaxUpdMult, max(self.MinUpdMult, (UpRatio - 1.0) * self.UpFrac + 1.0)) - UpRatioSoil = min(self.MaxUpdMult, max(self.MinUpdMult, (UpRatio - 1.0) * (1.0 - self.UpFrac) + 1.0)) + self.UpRatioKyn = min( + self.MaxUpdMult, + max(self.MinUpdMult, (UpRatio - 1.0) * self.UpFrac + 1.0), + ) + UpRatioSoil = min( + self.MaxUpdMult, + max(self.MinUpdMult, (UpRatio - 1.0) * (1.0 - self.UpFrac) + 1.0), + ) # update/nudge self.UStoreDepth for the whole upstream area, # not sure how much this helps or worsens things @@ -1864,59 +2889,104 @@ MM = (1.0 - self.UpRatioKyn) / self.UpdMaxDist self.UpRatioKyn = MM * self.DistToUpdPt + self.UpRatioKyn self.SurfaceRunoff = self.SurfaceRunoff * self.UpRatioKyn - self.SurfaceRunoffMM = self.SurfaceRunoff * self.QMMConv # SurfaceRunoffMM (mm) from SurfaceRunoff (m3/s) + self.SurfaceRunoffMM = ( + self.SurfaceRunoff * self.QMMConv + ) # SurfaceRunoffMM (mm) from SurfaceRunoff (m3/s) self.updateRunOff() Runoff = self.SurfaceRunoff - # Determine Soil moisture profile # self.vwc, self.vwcRoot: volumetric water content [m3/m3] per soil layer and root zone (including thetaR and saturated store) # self.vwc_perc, self.vwc_percRoot: volumetric water content [%] per soil layer and root zone (including thetaR and saturated store) # self.RootStore_sat: root water storage [mm] in saturated store (excluding thetaR) - # self.RootStore_unsat: root water storage [mm] in unsaturated store (excluding thetaR) - # self.RootStore: total root water storage [mm] (excluding thetaR) - - self.RootStore_sat = max(0.0,self.ActRootingDepth-self.zi)*(self.thetaS-self.thetaR) - + # self.RootStore_unsat: root water storage [mm] in unsaturated store (excluding thetaR) + # self.RootStore: total root water storage [mm] (excluding thetaR) + + self.RootStore_sat = max(0.0, self.ActRootingDepth - self.zi) * ( + self.thetaS - self.thetaR + ) + self.RootStore_unsat = self.ZeroMap self.SumThickness = self.ZeroMap self.vwc = [] self.vwc_perc = [] - + for n in arange(len(self.UStoreLayerThickness)): - - fracRoot = ifthenelse(self.ZiLayer > n, min(1.0,max(0.0,(min(self.ActRootingDepth,self.zi)-self.SumThickness)/self.UStoreLayerThickness[n])), - min(1.0,max(0.0, (self.ActRootingDepth-self.SumThickness) /(self.zi + 1 - self.SumThickness)))) - + + fracRoot = ifthenelse( + self.ZiLayer > n, + min( + 1.0, + max( + 0.0, + (min(self.ActRootingDepth, self.zi) - self.SumThickness) + / self.UStoreLayerThickness[n], + ), + ), + min( + 1.0, + max( + 0.0, + (self.ActRootingDepth - self.SumThickness) + / (self.zi + 1 - self.SumThickness), + ), + ), + ) + self.SumThickness = self.UStoreLayerThickness[n] + self.SumThickness - self.vwc.append(ifthenelse(self.ZiLayer > n, self.UStoreLayerDepth[n]/self.UStoreLayerThickness[n] + self.thetaR, - (((self.UStoreLayerDepth[n] + (self.thetaS-self.thetaR) * min(self.UStoreLayerThickness[n],(self.SumThickness-self.zi)))/self.UStoreLayerThickness[n]) + self.thetaR))) - - self.vwc_perc.append((self.vwc[n]/self.thetaS) * 100.0) - - self.RootStore_unsat = self.RootStore_unsat + cover(fracRoot*self.UStoreLayerDepth[n],0.0) - + self.vwc.append( + ifthenelse( + self.ZiLayer > n, + self.UStoreLayerDepth[n] / self.UStoreLayerThickness[n] + + self.thetaR, + ( + ( + ( + self.UStoreLayerDepth[n] + + (self.thetaS - self.thetaR) + * min( + self.UStoreLayerThickness[n], + (self.SumThickness - self.zi), + ) + ) + / self.UStoreLayerThickness[n] + ) + + self.thetaR + ), + ) + ) + self.vwc_perc.append((self.vwc[n] / self.thetaS) * 100.0) + + self.RootStore_unsat = self.RootStore_unsat + cover( + fracRoot * self.UStoreLayerDepth[n], 0.0 + ) + self.RootStore = self.RootStore_sat + self.RootStore_unsat - self.vwcRoot = self.RootStore/self.ActRootingDepth + self.thetaR - self.vwc_percRoot = (self.vwcRoot/self.thetaS) * 100.0 - + self.vwcRoot = self.RootStore / self.ActRootingDepth + self.thetaR + self.vwc_percRoot = (self.vwcRoot / self.thetaS) * 100.0 # 2: ########################################################################## # water balance ########################################### ########################################################################## self.QCatchmentMM = self.SurfaceRunoff * self.QMMConvUp - self.RunoffCoeff = self.QCatchmentMM/catchmenttotal(self.PrecipitationPlusMelt, self.TopoLdd)/catchmenttotal(cover(1.0), self.TopoLdd) - #self.AA = catchmenttotal(self.PrecipitationPlusMelt, self.TopoLdd) - #self.BB = catchmenttotal(cover(1.0), self.TopoLdd) + self.RunoffCoeff = ( + self.QCatchmentMM + / catchmenttotal(self.PrecipitationPlusMelt, self.TopoLdd) + / catchmenttotal(cover(1.0), self.TopoLdd) + ) + # self.AA = catchmenttotal(self.PrecipitationPlusMelt, self.TopoLdd) + # self.BB = catchmenttotal(cover(1.0), self.TopoLdd) # Single cell based water budget. snow not included yet. - self.CellStorage = sum_list_cover(self.UStoreLayerDepth,self.ZeroMap) + self.SatWaterDepth + self.CellStorage = ( + sum_list_cover(self.UStoreLayerDepth, self.ZeroMap) + self.SatWaterDepth + ) - self.sumUstore = sum_list_cover(self.UStoreLayerDepth,self.ZeroMap) + self.sumUstore = sum_list_cover(self.UStoreLayerDepth, self.ZeroMap) self.DeltaStorage = self.CellStorage - self.OrgStorage OutFlow = self.SatWaterFlux @@ -1940,18 +3010,38 @@ self.CumInwaterMM = self.CumInwaterMM + self.InwaterMM self.CumExfiltWater = self.CumExfiltWater + self.ExfiltWater - self.SoilWatbal = self.ActInfilt + self.reinfiltwater + CellInFlow - self.Transpiration - self.soilevap -\ - self.ExfiltWater - self.SubCellGWRunoff - self.DeltaStorage -\ - self.SatWaterFlux - self.InterceptionWatBal = self.PrecipitationPlusMelt - self.Interception -self.StemFlow - self.ThroughFall -\ - (self.OldCanopyStorage - self.CanopyStorage) - self.SurfaceWatbal = self.PrecipitationPlusMelt + self.oldIRSupplymm - self.Interception - \ - self.ExcessWater - self.RunoffOpenWater - self.SubCellRunoff - \ - self.ActInfilt -\ - (self.OldCanopyStorage - self.CanopyStorage) + self.SoilWatbal = ( + self.ActInfilt + + self.reinfiltwater + + CellInFlow + - self.Transpiration + - self.soilevap + - self.ExfiltWater + - self.SubCellGWRunoff + - self.DeltaStorage + - self.SatWaterFlux + ) + self.InterceptionWatBal = ( + self.PrecipitationPlusMelt + - self.Interception + - self.StemFlow + - self.ThroughFall + - (self.OldCanopyStorage - self.CanopyStorage) + ) + self.SurfaceWatbal = ( + self.PrecipitationPlusMelt + + self.oldIRSupplymm + - self.Interception + - self.ExcessWater + - self.RunoffOpenWater + - self.SubCellRunoff + - self.ActInfilt + - (self.OldCanopyStorage - self.CanopyStorage) + ) self.watbal = self.SoilWatbal + self.SurfaceWatbal + def main(argv=None): """ Perform command line execution of the model. @@ -1966,7 +3056,7 @@ runinfoFile = "runinfo.xml" timestepsecs = 86400 - wflow_cloneMap = 'wflow_subcatch.map' + wflow_cloneMap = "wflow_subcatch.map" _NoOverWrite = 1 global updateCols loglevel = logging.DEBUG @@ -1980,69 +3070,99 @@ ## Process command-line options # ######################################################################## try: - opts, args = getopt.getopt(argv, 'XL:hC:Ii:v:S:T:WR:u:s:EP:p:Xx:U:fOc:l:') + opts, args = getopt.getopt(argv, "XL:hC:Ii:v:S:T:WR:u:s:EP:p:Xx:U:fOc:l:") except getopt.error, msg: pcrut.usage(msg) for o, a in opts: - if o == '-C': caseName = a - if o == '-R': runId = a - if o == '-c': configfile = a - if o == '-L': LogFileName = a - if o == '-s': timestepsecs = int(a) - if o == '-h': usage() - if o == '-f': _NoOverWrite = 0 - if o == '-l': exec "loglevel = logging." + a + if o == "-C": + caseName = a + if o == "-R": + runId = a + if o == "-c": + configfile = a + if o == "-L": + LogFileName = a + if o == "-s": + timestepsecs = int(a) + if o == "-h": + usage() + if o == "-f": + _NoOverWrite = 0 + if o == "-l": + exec "loglevel = logging." + a + starttime = dt.datetime(1990, 01, 01) - starttime = dt.datetime(1990,01,01) - if _lastTimeStep < _firstTimeStep: - print "The starttimestep (" + str(_firstTimeStep) + ") is smaller than the last timestep (" + str( - _lastTimeStep) + ")" + print "The starttimestep (" + str( + _firstTimeStep + ) + ") is smaller than the last timestep (" + str(_lastTimeStep) + ")" usage() myModel = WflowModel(wflow_cloneMap, caseName, runId, configfile) - dynModelFw = wf_DynamicFramework(myModel, _lastTimeStep, firstTimestep=_firstTimeStep,datetimestart=starttime) - dynModelFw.createRunId(NoOverWrite=_NoOverWrite, level=loglevel, logfname=LogFileName,model="wflow_sbm",doSetupFramework=False) + dynModelFw = wf_DynamicFramework( + myModel, _lastTimeStep, firstTimestep=_firstTimeStep, datetimestart=starttime + ) + dynModelFw.createRunId( + NoOverWrite=_NoOverWrite, + level=loglevel, + logfname=LogFileName, + model="wflow_sbm", + doSetupFramework=False, + ) for o, a in opts: - if o == '-P': - left = a.split('=')[0] - right = a.split('=')[1] - configset(myModel.config,'variable_change_once',left,right,overwrite=True) - if o == '-p': - left = a.split('=')[0] - right = a.split('=')[1] - configset(myModel.config,'variable_change_timestep',left,right,overwrite=True) - if o == '-X': configset(myModel.config, 'model', 'OverWriteInit', '1', overwrite=True) - if o == '-I': configset(myModel.config, 'run', 'reinit', '1', overwrite=True) - if o == '-i': configset(myModel.config, 'model', 'intbl', a, overwrite=True) - if o == '-s': configset(myModel.config, 'model', 'timestepsecs', a, overwrite=True) - if o == '-x': configset(myModel.config, 'model', 'sCatch', a, overwrite=True) - if o == '-c': configset(myModel.config, 'model', 'configfile', a, overwrite=True) - if o == '-M': configset(myModel.config, 'model', 'MassWasting', "0", overwrite=True) - if o == '-Q': configset(myModel.config, 'model', 'ExternalQbase', '1', overwrite=True) - if o == '-U': - configset(myModel.config, 'model', 'updateFile', a, overwrite=True) - configset(myModel.config, 'model', 'updating', "1", overwrite=True) - if o == '-u': + if o == "-P": + left = a.split("=")[0] + right = a.split("=")[1] + configset( + myModel.config, "variable_change_once", left, right, overwrite=True + ) + if o == "-p": + left = a.split("=")[0] + right = a.split("=")[1] + configset( + myModel.config, "variable_change_timestep", left, right, overwrite=True + ) + if o == "-X": + configset(myModel.config, "model", "OverWriteInit", "1", overwrite=True) + if o == "-I": + configset(myModel.config, "run", "reinit", "1", overwrite=True) + if o == "-i": + configset(myModel.config, "model", "intbl", a, overwrite=True) + if o == "-s": + configset(myModel.config, "model", "timestepsecs", a, overwrite=True) + if o == "-x": + configset(myModel.config, "model", "sCatch", a, overwrite=True) + if o == "-c": + configset(myModel.config, "model", "configfile", a, overwrite=True) + if o == "-M": + configset(myModel.config, "model", "MassWasting", "0", overwrite=True) + if o == "-Q": + configset(myModel.config, "model", "ExternalQbase", "1", overwrite=True) + if o == "-U": + configset(myModel.config, "model", "updateFile", a, overwrite=True) + configset(myModel.config, "model", "updating", "1", overwrite=True) + if o == "-u": zz = [] exec "zz =" + a updateCols = zz - if o == '-E': configset(myModel.config, 'model', 'reInfilt', '1', overwrite=True) - if o == '-R': runId = a - if o == '-W': configset(myModel.config, 'model', 'waterdem', '1', overwrite=True) - if o == '-T': - configset(myModel.config, 'run', 'endtime', a, overwrite=True) - if o == '-S': - configset(myModel.config, 'run', 'starttime', a, overwrite=True) + if o == "-E": + configset(myModel.config, "model", "reInfilt", "1", overwrite=True) + if o == "-R": + runId = a + if o == "-W": + configset(myModel.config, "model", "waterdem", "1", overwrite=True) + if o == "-T": + configset(myModel.config, "run", "endtime", a, overwrite=True) + if o == "-S": + configset(myModel.config, "run", "starttime", a, overwrite=True) - dynModelFw.setupFramework() dynModelFw._runInitial() dynModelFw._runResume() - #dynModelFw._runDynamic(0, 0) + # dynModelFw._runDynamic(0, 0) dynModelFw._runDynamic(_firstTimeStep, _lastTimeStep) dynModelFw._runSuspend() dynModelFw._wf_shutdown()