Index: wflow-py/wflow/wflow_hbv.py
===================================================================
diff -u -r70d5816d74c2254eb2f2c848d5b65c19a8ff4bf7 -r3e39e84af48f1bcb5ec0d243748147be223674f2
--- wflow-py/wflow/wflow_hbv.py (.../wflow_hbv.py) (revision 70d5816d74c2254eb2f2c848d5b65c19a8ff4bf7)
+++ wflow-py/wflow/wflow_hbv.py (.../wflow_hbv.py) (revision 3e39e84af48f1bcb5ec0d243748147be223674f2)
@@ -16,7 +16,7 @@
# You should have received a copy of the GNU General Public License
# along with this program. If not, see .
-#TODO: split off routing
+# TODO: split off routing
"""
Run the wflow_hbv hydrological model..
@@ -91,16 +91,15 @@
from wflow.wflow_adapt import *
from wflow.wflow_adapt import *
-#import scipy
-#import pcrut
+# import scipy
+# import pcrut
-
wflow = "wflow_hbv"
#: columns used in updating
-updateCols = [] #: columns used in updating
+updateCols = [] #: columns used in updating
""" Column used in updating """
@@ -111,45 +110,44 @@
- *args: command line arguments given
"""
sys.stdout = sys.stderr
- for msg in args: print(msg)
+ for msg in args:
+ print(msg)
print(__doc__)
sys.exit(0)
+
class WflowModel(DynamicModel):
- """
+ """
The user defined model class.
"""
+ def __init__(self, cloneMap, Dir, RunDir, configfile):
+ DynamicModel.__init__(self)
+ self.caseName = os.path.abspath(Dir)
+ 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)
- def __init__(self, cloneMap,Dir,RunDir,configfile):
- DynamicModel.__init__(self)
- self.caseName = os.path.abspath(Dir)
- 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)
-
-
- def updateRunOff(self):
- """
+ def updateRunOff(self):
+ """
Updates the kinematic wave reservoir
"""
- self.WaterLevel=(self.Alpha*pow(self.SurfaceRunoff,self.Beta))/self.Bw
- # wetted perimeter (m)
- P=self.Bw+(2*self.WaterLevel)
- # Alpha
- self.Alpha=self.AlpTerm*pow(P,self.AlpPow)
- self.OldKinWaveVolume = self.KinWaveVolume
- self.KinWaveVolume = self.WaterLevel * self.Bw * self.DCL
+ self.WaterLevel = (self.Alpha * pow(self.SurfaceRunoff, self.Beta)) / self.Bw
+ # wetted perimeter (m)
+ P = self.Bw + (2 * self.WaterLevel)
+ # Alpha
+ self.Alpha = self.AlpTerm * pow(P, self.AlpPow)
+ self.OldKinWaveVolume = self.KinWaveVolume
+ self.KinWaveVolume = self.WaterLevel * self.Bw * self.DCL
-
- def stateVariables(self):
- """
+ def stateVariables(self):
+ """
returns a list of state variables that are essential to the model.
This list is essential for the resume and suspend functions to work.
@@ -165,93 +163,140 @@
:var self.InterceptionStorage: Amount of water on the Canopy [mm]
"""
- states = ['FreeWater', 'SoilMoisture',
- 'UpperZoneStorage',
- 'LowerZoneStorage',
- 'InterceptionStorage',
- 'SurfaceRunoff',
- 'WaterLevel',
- 'DrySnow']
+ states = [
+ "FreeWater",
+ "SoilMoisture",
+ "UpperZoneStorage",
+ "LowerZoneStorage",
+ "InterceptionStorage",
+ "SurfaceRunoff",
+ "WaterLevel",
+ "DrySnow",
+ ]
- if hasattr(self,'ReserVoirSimpleLocs'):
- states.append('ReservoirVolume')
+ if hasattr(self, "ReserVoirSimpleLocs"):
+ states.append("ReservoirVolume")
- if hasattr(self,'ReserVoirComplexLocs'):
- states.append('ReservoirWaterLevel')
+ if hasattr(self, "ReserVoirComplexLocs"):
+ states.append("ReservoirWaterLevel")
- return states
+ return states
-
# The following are made to better connect to deltashell/openmi
- def supplyCurrentTime(self):
- """
+ def supplyCurrentTime(self):
+ """
gets the current time in seconds after the start of the run
Ouput:
- time in seconds since the start of the model run
"""
- return self.currentTimeStep() * int(configget(self.config,'model','timestepsecs','86400'))
+ return self.currentTimeStep() * int(
+ configget(self.config, "model", "timestepsecs", "86400")
+ )
- def parameters(self):
- """
+ def parameters(self):
+ """
Define all model parameters here that the framework should handle for the model
See wf_updateparameters and the parameters section of the ini file
If you use this make sure to all wf_updateparameters at the start of the dynamic section
and at the start/end of the initial section
"""
- modelparameters = []
+ 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))
- # Meteo and other forcing
+ # Meteo and other forcing
- 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.Seepage_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Seepage",
- "/inmaps/SE") # timeseries for rainfall "/inmaps/SE" # 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="PotEvaporation",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="Seepage",stack=self.Seepage_mapstack,type="timeseries",default=0.0,verbose=False,lookupmaps=[]))
+ 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.Seepage_mapstack = self.Dir + configget(
+ self.config, "inputmapstacks", "Seepage", "/inmaps/SE"
+ ) # timeseries for rainfall "/inmaps/SE" # 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="PotEvaporation",
+ 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="Seepage",
+ stack=self.Seepage_mapstack,
+ type="timeseries",
+ default=0.0,
+ verbose=False,
+ lookupmaps=[],
+ )
+ )
+ return modelparameters
-
- return modelparameters
-
-
- def suspend(self):
- """
+ def suspend(self):
+ """
Suspends the model to disk. All variables needed to restart the model
are saved to disk as pcraster maps. Use resume() to re-read them
"""
+ self.logger.info("Saving initial conditions...")
+ self.wf_suspend(os.path.join(self.SaveDir, "outstate"))
- self.logger.info("Saving initial conditions...")
- self.wf_suspend(os.path.join(self.SaveDir,"outstate"))
+ if self.OverWriteInit:
+ self.logger.info("Saving initial conditions over start conditions...")
+ self.wf_suspend(os.path.join(self.SaveDir, "instate"))
- if self.OverWriteInit:
- self.logger.info("Saving initial conditions over start conditions...")
- self.wf_suspend(os.path.join(self.SaveDir,"instate"))
+ if self.fewsrun:
+ self.logger.info("Saving initial conditions for FEWS...")
+ self.wf_suspend(os.path.join(self.Dir, "outstate"))
+ def initial(self):
- if self.fewsrun:
- self.logger.info("Saving initial conditions for FEWS...")
- self.wf_suspend(os.path.join(self.Dir, "outstate"))
-
-
-
- def initial(self):
-
- """
+ """
Initial part of the model, executed only once. Reads all static model
information (parameters) and sets-up the variables used in modelling.
@@ -295,406 +340,720 @@
"""
- global statistics
- global multpars
- global updateCols
+ global statistics
+ global multpars
+ global updateCols
- setglobaloption("unittrue")
+ setglobaloption("unittrue")
+ self.thestep = scalar(0)
- self.thestep = scalar(0)
+ #: files to be used in case of timesries (scalar) input to the model
- #: files to be used in case of timesries (scalar) input to the model
+ #: name of the tss file with precipitation data ("../intss/P.tss")
+ self.precipTss = "../intss/P.tss"
+ self.evapTss = (
+ "../intss/PET.tss"
+ ) #: name of the tss file with potential evap data ("../intss/PET.tss")
+ self.tempTss = (
+ "../intss/T.tss"
+ ) #: name of the tss file with temperature data ("../intss/T.tss")
+ self.inflowTss = (
+ "../intss/Inflow.tss"
+ ) #: NOT TESTED name of the tss file with inflow data ("../intss/Inflow.tss")
+ self.SeepageTss = (
+ "../intss/Seepage.tss"
+ ) #: NOT TESTED name of the tss file with seepage data ("../intss/Seepage.tss")"
- #: name of the tss file with precipitation data ("../intss/P.tss")
- self.precipTss = "../intss/P.tss"
- self.evapTss="../intss/PET.tss" #: name of the tss file with potential evap data ("../intss/PET.tss")
- self.tempTss="../intss/T.tss" #: name of the tss file with temperature data ("../intss/T.tss")
- self.inflowTss="../intss/Inflow.tss" #: NOT TESTED name of the tss file with inflow data ("../intss/Inflow.tss")
- self.SeepageTss="../intss/Seepage.tss" #: NOT TESTED name of the tss file with seepage data ("../intss/Seepage.tss")"
+ self.logger.info("running for " + str(self.nrTimeSteps()) + " timesteps")
-
- self.logger.info("running for " + str(self.nrTimeSteps()) + " timesteps")
-
-
# Set and get defaults from ConfigFile here ###################################
- self.scalarInput = int(configget(self.config,"model","ScalarInput","0"))
- self.Tslice = int(configget(self.config,"model","Tslice","1"))
- self.interpolMethod = configget(self.config,"model","InterpolationMethod","inv")
- self.reinit = int(configget(self.config,"run","reinit","0"))
- self.fewsrun = int(configget(self.config,"run","fewsrun","0"))
- self.OverWriteInit = int(configget(self.config,"model","OverWriteInit","0"))
- self.updating = int(configget(self.config,"model","updating","0"))
- self.updateFile = configget(self.config,"model","updateFile","no_set")
+ self.scalarInput = int(configget(self.config, "model", "ScalarInput", "0"))
+ self.Tslice = int(configget(self.config, "model", "Tslice", "1"))
+ self.interpolMethod = configget(
+ self.config, "model", "InterpolationMethod", "inv"
+ )
+ self.reinit = int(configget(self.config, "run", "reinit", "0"))
+ self.fewsrun = int(configget(self.config, "run", "fewsrun", "0"))
+ self.OverWriteInit = int(configget(self.config, "model", "OverWriteInit", "0"))
+ self.updating = int(configget(self.config, "model", "updating", "0"))
+ self.updateFile = configget(self.config, "model", "updateFile", "no_set")
- self.sCatch = int(configget(self.config,"model","sCatch","0"))
- self.intbl = configget(self.config,"model","intbl","intbl")
- self.P_style = int(configget(self.config,"model","P_style","1"))
- self.PET_style = int(configget(self.config,"model","PET_style","1"))
- self.TEMP_style = int(configget(self.config,"model","TEMP_style","1"))
+ self.sCatch = int(configget(self.config, "model", "sCatch", "0"))
+ self.intbl = configget(self.config, "model", "intbl", "intbl")
+ self.P_style = int(configget(self.config, "model", "P_style", "1"))
+ self.PET_style = int(configget(self.config, "model", "PET_style", "1"))
+ self.TEMP_style = int(configget(self.config, "model", "TEMP_style", "1"))
+ 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"))
+ 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.SetKquickFlow = int(configget(self.config, "model", "SetKquickFlow", "0"))
+ self.MassWasting = int(configget(self.config, "model", "MassWasting", "0"))
+ self.SubCatchFlowOnly = int(
+ configget(self.config, "model", "SubCatchFlowOnly", "0")
+ )
- 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"))
- 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.SetKquickFlow=int(configget(self.config,'model','SetKquickFlow','0'))
- self.MassWasting = int(configget(self.config,"model","MassWasting","0"))
- self.SubCatchFlowOnly = int(configget(self.config, 'model', 'SubCatchFlowOnly', '0'))
+ # 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_mgauges = configget(
+ self.config, "model", "wflow_mgauges", "staticmaps/wflow_mgauges.map"
+ )
+ wflow_riverwidth = configget(
+ self.config, "model", "wflow_riverwidth", "staticmaps/wflow_riverwidth.map"
+ )
- # 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_mgauges = configget(self.config,"model","wflow_mgauges","staticmaps/wflow_mgauges.map")
- wflow_riverwidth = configget(self.config,"model","wflow_riverwidth","staticmaps/wflow_riverwidth.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 = ifthen(subcatch > 0, subcatch)
+ if self.sCatch > 0:
+ subcatch = ifthen(subcatch == sCatch, subcatch)
+ self.Altitude = self.wf_readmap(
+ os.path.join(self.Dir, wflow_dem), 0.0, fail=True
+ ) * scalar(
+ defined(subcatch)
+ ) #: The digital elevation map (DEM)
+ self.TopoLdd = self.wf_readmap(
+ os.path.join(self.Dir, wflow_ldd), 0.0, fail=True
+ ) #: The local drinage definition map (ldd)
+ self.TopoId = ordinal(
+ self.wf_readmap(os.path.join(self.Dir, wflow_subcatch), 0.0, fail=True)
+ ) #: Map define the area over which the calculations are done (mask)
+ self.River = cover(
+ boolean(
+ self.wf_readmap(os.path.join(self.Dir, wflow_river), 0.0, fail=True)
+ ),
+ 0,
+ ) #: river network map. Fro those cell that belong to a river a specific width is used in the kinematic wave caulations
+ self.RiverLength = self.wf_readmap(
+ os.path.join(self.Dir, wflow_riverlength), 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
+ )
+ # 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 = self.wf_readmap(
+ os.path.join(self.Dir, wflow_landuse), 0.0, fail=True
+ ) #: Map with lan-use/cover classes
+ self.LandUse = cover(self.LandUse, nominal(ordinal(subcatch) > 0))
+ self.Soil = self.wf_readmap(
+ os.path.join(self.Dir, wflow_soil), 0.0, fail=True
+ ) #: Map with soil classes
+ self.Soil = cover(self.Soil, nominal(ordinal(subcatch) > 0))
+ self.OutputLoc = self.wf_readmap(
+ os.path.join(self.Dir, wflow_gauges), 0.0, fail=True
+ ) #: Map with locations of output gauge(s)
+ self.InflowLoc = nominal(
+ self.wf_readmap(os.path.join(self.Dir, wflow_inflow), 0.0)
+ ) #: Map with location of abstractions/inflows.
+ self.SeepageLoc = self.wf_readmap(
+ os.path.join(self.Dir, wflow_inflow), 0.0
+ ) #: Seapage from external model (if configured)
+ RiverWidth = self.wf_readmap(os.path.join(self.Dir, wflow_riverwidth), 0.0)
- # 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 = ifthen(subcatch > 0, subcatch)
- if self.sCatch > 0:
- subcatch = ifthen(subcatch == sCatch,subcatch)
+ # Temperature correction per cell to add
+ self.TempCor = self.wf_readmap(
+ os.path.join(
+ self.Dir,
+ configget(
+ self.config,
+ "model",
+ "TemperatureCorrectionMap",
+ "staticmap/swflow_tempcor.map",
+ ),
+ ),
+ 0.0,
+ )
- self.Altitude=self.wf_readmap(os.path.join(self.Dir,wflow_dem),0.0,fail=True) * scalar(defined(subcatch)) #: The digital elevation map (DEM)
- self.TopoLdd=self.wf_readmap(os.path.join(self.Dir, wflow_ldd),0.0,fail=True) #: The local drinage definition map (ldd)
- self.TopoId=ordinal(self.wf_readmap(os.path.join(self.Dir, wflow_subcatch),0.0,fail=True) ) #: Map define the area over which the calculations are done (mask)
- self.River=cover(boolean(self.wf_readmap(os.path.join(self.Dir, wflow_river),0.0,fail=True)),0) #: river network map. Fro those cell that belong to a river a specific width is used in the kinematic wave caulations
- self.RiverLength=self.wf_readmap(os.path.join(self.Dir, wflow_riverlength),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)
+ if self.scalarInput:
+ self.gaugesMap = self.wf_readmap(
+ os.path.join(self.Dir, wflow_mgauges), 0.0, fail=True
+ ) #: Map with locations of rainfall/evap/temp gauge(s). Only needed if the input to the model is not in maps
+ self.OutputId = self.wf_readmap(
+ os.path.join(self.Dir, wflow_subcatch), 0.0, fail=True
+ ) # location of subcatchment
- # 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=self.wf_readmap(os.path.join(self.Dir , wflow_landuse),0.0,fail=True)#: Map with lan-use/cover classes
- self.LandUse=cover(self.LandUse,nominal(ordinal(subcatch) > 0))
- self.Soil=self.wf_readmap(os.path.join(self.Dir , wflow_soil),0.0,fail=True)#: Map with soil classes
- self.Soil=cover(self.Soil,nominal(ordinal(subcatch) > 0))
- self.OutputLoc=self.wf_readmap(os.path.join(self.Dir , wflow_gauges),0.0,fail=True) #: Map with locations of output gauge(s)
- self.InflowLoc=nominal(self.wf_readmap(os.path.join(self.Dir , wflow_inflow),0.0)) #: Map with location of abstractions/inflows.
- self.SeepageLoc=self.wf_readmap(os.path.join(self.Dir , wflow_inflow),0.0) #: Seapage from external model (if configured)
- RiverWidth=self.wf_readmap(os.path.join(self.Dir, wflow_riverwidth),0.0)
+ self.ZeroMap = 0.0 * scalar(defined(self.Altitude)) # map with only zero's
+ # 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.Seepage_mapstack = self.Dir + configget(
+ self.config, "inputmapstacks", "Seepage", "/inmaps/SE"
+ ) # timeseries for rainfall "/inmaps/SE" # in/outflow locations (abstractions)
+ # For in memory override:
+ self.P = self.ZeroMap
+ self.PET = self.ZeroMap
+ self.TEMP = self.ZeroMap
+ # Set static initial values here #########################################
- # Temperature correction per cell to add
- self.TempCor=self.wf_readmap(os.path.join(self.Dir , configget(self.config,"model","TemperatureCorrectionMap","staticmap/swflow_tempcor.map")),0.0)
+ self.Latitude = ycoordinate(boolean(self.Altitude))
+ self.Longitude = xcoordinate(boolean(self.Altitude))
+ self.logger.info("Linking parameters to landuse, catchment and soil...")
- if self.scalarInput:
- self.gaugesMap=self.wf_readmap(os.path.join(self.Dir , wflow_mgauges),0.0,fail=True) #: Map with locations of rainfall/evap/temp gauge(s). Only needed if the input to the model is not in maps
- self.OutputId=self.wf_readmap(os.path.join(self.Dir , wflow_subcatch),0.0,fail=True) # location of subcatchment
+ self.Beta = scalar(0.6) # For sheetflow
+ # self.M=lookupscalar(self.Dir + "/" + modelEnv['intbl'] + "/M.tbl" ,self.LandUse,subcatch,self.Soil) # Decay parameter in Topog_sbm
+ self.N = lookupscalar(
+ self.Dir + "/" + self.intbl + "/N.tbl", self.LandUse, subcatch, self.Soil
+ ) # Manning overland flow
+ """ *Parameter:* Manning's N for all non-river cells """
+ self.NRiver = lookupscalar(
+ self.Dir + "/" + self.intbl + "/N_River.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ ) # Manning river
+ """ Manning's N for all cells that are marked as a river """
- self.ZeroMap=0.0*scalar(defined(self.Altitude)) #map with only zero's
+ self.wf_updateparameters()
- # 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.Seepage_mapstack=self.Dir + configget(self.config,"inputmapstacks","Seepage","/inmaps/SE") # timeseries for rainfall "/inmaps/SE" # in/outflow locations (abstractions)
- # For in memory override:
- self.P = self.ZeroMap
- self.PET = self.ZeroMap
- self.TEMP = self.ZeroMap
- # Set static initial values here #########################################
+ self.ReserVoirLocs = self.ZeroMap
+ 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)
+ )
+ else:
+ self.nrresSimple = 0
- self.Latitude = ycoordinate(boolean(self.Altitude))
- self.Longitude = xcoordinate(boolean(self.Altitude))
+ 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, res_area * 0.0 + 1.0
+ )
- self.logger.info("Linking parameters to landuse, catchment and soil...")
+ # read files
+ self.sh = {}
+ res_ids = ifthen(self.ResStorFunc == 2, self.ReserVoirComplexLocs)
+ 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.hq = {}
+ res_ids = ifthen(self.ResOutflowFunc == 1, self.ReserVoirComplexLocs)
+ 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.Beta = scalar(0.6) # For sheetflow
- #self.M=lookupscalar(self.Dir + "/" + modelEnv['intbl'] + "/M.tbl" ,self.LandUse,subcatch,self.Soil) # Decay parameter in Topog_sbm
- self.N=lookupscalar(self.Dir + "/" + self.intbl + "/N.tbl",self.LandUse,subcatch,self.Soil) # Manning overland flow
- """ *Parameter:* Manning's N for all non-river cells """
- self.NRiver=lookupscalar(self.Dir + "/" + self.intbl + "/N_River.tbl",self.LandUse,subcatch,self.Soil) # Manning river
- """ Manning's N for all cells that are marked as a river """
+ else:
+ self.nrresComplex = 0
- self.wf_updateparameters()
+ 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.ReserVoirLocs = self.ZeroMap
+ # HBV Soil params
+ self.FC = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/FC.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 260.0,
+ )
+ self.BetaSeepage = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/BetaSeepage.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 1.8,
+ ) # exponent in soil runoff generation equation
+ self.LP = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/LP.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 0.53000,
+ ) # fraction of Fieldcapacity below which actual evaporation=potential evaporation (LP)
+ self.K4 = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/K4.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 0.02307,
+ ) # Recession constant baseflow #K4=0.07; BASEFLOW:LINEARRESERVOIR
+ if self.SetKquickFlow:
+ self.KQuickFlow = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/KQuickFlow.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 0.09880,
+ ) # recession rate at flow HQ #KHQ=0.2; OUTFLOWUPPERZONE_NONLINEARRESERVOIR
+ self.SUZ = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/SUZ.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 100.0,
+ ) # Level over wich K0 is used
+ self.K0 = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/K0.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 0.3,
+ ) # K0
+ else:
+ self.KHQ = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/KHQ.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 0.09880,
+ ) # recession rate at flow HQ #KHQ=0.2; OUTFLOWUPPERZONE_NONLINEARRESERVOIR
+ self.HQ = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/HQ.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 3.27000,
+ ) # high flow rate HQ for which recession rate of upper reservoir is known #HQ=3.76;
+ self.AlphaNL = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/AlphaNL.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 1.1,
+ ) # measure of non-linearity of upper reservoir #Alpha=1.6;
- 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))
- else:
- self.nrresSimple = 0
+ self.PERC = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/PERC.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 0.4000,
+ ) # percolation from Upper to Lowerzone (mm/day)
+ self.CFR = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/CFR.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 0.05000,
+ ) # refreezing efficiency constant in refreezing of freewater in snow
+ # self.FoCfmax=self.readtblDefault(self.Dir + "/" + modelEnv['intbl'] + "/FoCfmax.tbl",self.LandUse,subcatch,self.Soil, 0.6000) # correcton factor for snow melt/refreezing in forested and non-forested areas
+ self.Pcorr = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/Pcorr.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 1.0,
+ ) # correction factor for precipitation
+ self.RFCF = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/RFCF.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 1.0,
+ ) # correction factor for rainfall
+ self.SFCF = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/SFCF.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 1.0,
+ ) # correction factor for snowfall
+ self.Cflux = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/Cflux.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 2.0,
+ ) # maximum capillary rise from runoff response routine to soil moisture routine
+ self.ICF = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/ICF.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 2.0,
+ ) # maximum interception storage (in forested AND non-forested areas)
+ self.CEVPF = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/CEVPF.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 1.0,
+ ) # correction factor for potential evaporation (1.15 in in forested areas )
+ self.EPF = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/EPF.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 0.0,
+ ) # exponent of correction factor for evaporation on days with precipitation
+ self.ECORR = self.readtblDefault(
+ self.Dir + "/" + self.intbl + "/ECORR.tbl",
+ self.LandUse,
+ subcatch,
+ self.Soil,
+ 1.0,
+ ) # evap correction
+ # Soil Moisture parameters
+ self.ECALT = self.ZeroMap + 0.00000 # evaporation lapse per 100m
+ # self.Ecorr=self.ZeroMap+1 # correction factor for evaporation
+ # 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,
+ )
+ # 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,
+ )
+ # 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,
+ )
- 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, res_area*0.0 + 1.0)
+ # Determine real slope and cell length
+ 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),
+ )
+ Terrain_angle = scalar(atan(self.Slope))
+ temp = (
+ catchmenttotal(cover(1.0), self.TopoLdd)
+ * self.reallength
+ * 0.001
+ * 0.001
+ * self.reallength
+ )
+ self.QMMConvUp = cover(self.timestepsecs * 0.001) / temp
- #read files
- self.sh = {}
- res_ids = ifthen(self.ResStorFunc == 2, self.ReserVoirComplexLocs)
- 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.hq = {}
- res_ids = ifthen(self.ResOutflowFunc == 1, self.ReserVoirComplexLocs)
- 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)
+ # Multiply parameters with a factor (for calibration etc) -P option in command line
+ self.wf_multparameters()
+ self.N = ifthenelse(self.River, self.NRiver, self.N)
+ # 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
+ # "Noah J. Finnegan et al 2005 Controls on the channel width of rivers:
+ # Implications for modeling fluvial incision of bedrock"
+ 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)
+ )
+ # Use supplied riverwidth if possible, else calulate
+ RiverWidth = ifthenelse(RiverWidth <= 0.0, W, RiverWidth)
- else:
- self.nrresComplex = 0
+ self.SnowWater = self.ZeroMap
+ # Which columns/gauges to use/ignore in kinematic wave updating
+ self.UpdateMap = self.ZeroMap
- 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))
+ if self.updating:
+ _tmp = pcr2numpy(self.OutputLoc, 0.0)
+ gaugear = _tmp
+ touse = numpy.zeros(gaugear.shape, dtype="int")
+ for thecol in updateCols:
+ idx = (gaugear == thecol).nonzero()
+ touse[idx] = thecol
+ self.UpdateMap = numpy2pcr(Nominal, touse, 0.0)
+ # 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,
+ )
+ # self.DistToUpdPt = ldddist(self.TopoLdd,boolean(cover(self.OutputId,0.0)),1)
+ # * self.reallength/celllength()
- #HBV Soil params
- self.FC=self.readtblDefault(self.Dir + "/" + self.intbl + "/FC.tbl",self.LandUse,subcatch,self.Soil,260.0)
- self.BetaSeepage= self.readtblDefault(self.Dir + "/" + self.intbl + "/BetaSeepage.tbl",self.LandUse,subcatch,self.Soil,1.8) # exponent in soil runoff generation equation
- self.LP= self.readtblDefault(self.Dir + "/" + self.intbl + "/LP.tbl",self.LandUse,subcatch,self.Soil, 0.53000) # fraction of Fieldcapacity below which actual evaporation=potential evaporation (LP)
- self.K4= self.readtblDefault(self.Dir + "/" + self.intbl + "/K4.tbl",self.LandUse,subcatch,self.Soil, 0.02307) # Recession constant baseflow #K4=0.07; BASEFLOW:LINEARRESERVOIR
- if self.SetKquickFlow:
- self.KQuickFlow= self.readtblDefault(self.Dir + "/" + self.intbl + "/KQuickFlow.tbl",self.LandUse,subcatch,self.Soil, 0.09880) # recession rate at flow HQ #KHQ=0.2; OUTFLOWUPPERZONE_NONLINEARRESERVOIR
- self.SUZ= self.readtblDefault(self.Dir + "/" + self.intbl + "/SUZ.tbl",self.LandUse,subcatch,self.Soil, 100.0) # Level over wich K0 is used
- self.K0= self.readtblDefault(self.Dir + "/" + self.intbl + "/K0.tbl",self.LandUse,subcatch,self.Soil, 0.3) # K0
- else:
- self.KHQ= self.readtblDefault(self.Dir + "/" + self.intbl + "/KHQ.tbl",self.LandUse,subcatch,self.Soil, 0.09880) # recession rate at flow HQ #KHQ=0.2; OUTFLOWUPPERZONE_NONLINEARRESERVOIR
- self.HQ= self.readtblDefault(self.Dir + "/" + self.intbl + "/HQ.tbl",self.LandUse,subcatch,self.Soil, 3.27000) # high flow rate HQ for which recession rate of upper reservoir is known #HQ=3.76;
- self.AlphaNL= self.readtblDefault(self.Dir + "/" + self.intbl + "/AlphaNL.tbl",self.LandUse,subcatch,self.Soil, 1.1) # measure of non-linearity of upper reservoir #Alpha=1.6;
+ # Initializing of variables
+ self.logger.info("Initializing of model variables..")
+ self.TopoLdd = lddmask(self.TopoLdd, boolean(self.TopoId))
+ catchmentcells = maptotal(scalar(self.TopoId))
- self.PERC= self.readtblDefault(self.Dir + "/" + self.intbl + "/PERC.tbl",self.LandUse,subcatch,self.Soil, 0.4000) # percolation from Upper to Lowerzone (mm/day)
- self.CFR=self.readtblDefault(self.Dir + "/" + self.intbl + "/CFR.tbl",self.LandUse,subcatch,self.Soil, 0.05000) # refreezing efficiency constant in refreezing of freewater in snow
- #self.FoCfmax=self.readtblDefault(self.Dir + "/" + modelEnv['intbl'] + "/FoCfmax.tbl",self.LandUse,subcatch,self.Soil, 0.6000) # correcton factor for snow melt/refreezing in forested and non-forested areas
- self.Pcorr=self.readtblDefault(self.Dir + "/" + self.intbl + "/Pcorr.tbl",self.LandUse,subcatch,self.Soil, 1.0) # correction factor for precipitation
- self.RFCF=self.readtblDefault(self.Dir + "/" + self.intbl + "/RFCF.tbl",self.LandUse,subcatch,self.Soil,1.0) # correction factor for rainfall
- self.SFCF=self.readtblDefault(self.Dir + "/" + self.intbl + "/SFCF.tbl",self.LandUse,subcatch,self.Soil, 1.0) # correction factor for snowfall
- self.Cflux= self.readtblDefault(self.Dir + "/" + self.intbl + "/Cflux.tbl",self.LandUse,subcatch,self.Soil, 2.0) # maximum capillary rise from runoff response routine to soil moisture routine
- self.ICF= self.readtblDefault(self.Dir + "/" + self.intbl + "/ICF.tbl",self.LandUse,subcatch,self.Soil, 2.0) # maximum interception storage (in forested AND non-forested areas)
- self.CEVPF= self.readtblDefault(self.Dir + "/" + self.intbl + "/CEVPF.tbl",self.LandUse,subcatch,self.Soil, 1.0) # correction factor for potential evaporation (1.15 in in forested areas )
- self.EPF= self.readtblDefault(self.Dir + "/" + self.intbl + "/EPF.tbl",self.LandUse,subcatch,self.Soil, 0.0) # exponent of correction factor for evaporation on days with precipitation
- self.ECORR= self.readtblDefault(self.Dir + "/" + self.intbl + "/ECORR.tbl",self.LandUse,subcatch,self.Soil, 1.0) # evap correction
- # Soil Moisture parameters
- self.ECALT= self.ZeroMap+0.00000 # evaporation lapse per 100m
- #self.Ecorr=self.ZeroMap+1 # correction factor for evaporation
+ # Limit lateral flow per subcatchment (make pits at all subcatch boundaries)
+ # 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))
+ # Used to seperate output per LandUse/management classes
+ # OutZones = self.LandUse
+ # report(self.reallength,"rl.map")
+ # report(catchmentcells,"kk.map")
+ self.QMMConv = self.timestepsecs / (
+ self.reallength * self.reallength * 0.001
+ ) # m3/s --> mm
+ self.ToCubic = (
+ self.reallength * self.reallength * 0.001
+ ) / self.timestepsecs # m3/s
+ 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 (weigthted for upstream area)
+ self.sumlevel = self.ZeroMap #: accumulated level for water balance
+ self.sumpotevap = self.ZeroMap # accumulated runoff for water balance
+ self.sumsoilevap = self.ZeroMap
+ self.sumtemp = self.ZeroMap # accumulated runoff for water balance
+ self.ForecQ_qmec = (
+ self.ZeroMap
+ ) # Extra inflow to kinematic wave reservoir for forcing in m^/sec
+ self.KinWaveVolume = self.ZeroMap
+ self.OldKinWaveVolume = self.ZeroMap
+ self.Qvolume = self.ZeroMap
+ self.Q = self.ZeroMap
+ self.suminflow = self.ZeroMap
+ # cntd
+ self.FieldCapacity = self.FC #: total water holding capacity of the soil
+ self.Treshold = (
+ self.LP * self.FieldCapacity
+ ) # Threshold soilwaterstorage above which AE=PE
+ # CatSurface=maptotal(scalar(ifthen(scalar(self.TopoId)>scalar(0.0),scalar(1.0)))) # catchment surface (in km2)
- # 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)
- # 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)
- #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.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)
+ )
- # Determine real slope and cell length
- 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))
- Terrain_angle=scalar(atan(self.Slope))
- temp = catchmenttotal(cover(1.0), self.TopoLdd) * self.reallength * 0.001 * 0.001 * self.reallength
- self.QMMConvUp = cover(self.timestepsecs * 0.001)/temp
+ # Set DCL to riverlength if that is longer that the basic length calculated from grid
+ drainlength = detdrainlength(self.TopoLdd, self.xl, self.yl)
+ self.DCL = max(drainlength, self.RiverLength) # m
+ # Multiply with Factor (taken from upscaling operation, defaults to 1.0 if no map is supplied
+ self.DCL = self.DCL * max(1.0, self.RiverLengthFac)
- # Multiply parameters with a factor (for calibration etc) -P option in command line
+ # water depth (m)
+ # set width for kinematic wave to cell width for all cells
+ self.Bw = detdrainwidth(self.TopoLdd, self.xl, self.yl)
+ # However, in the main river we have real flow so set the width to the
+ # width of the river
- self.wf_multparameters()
- self.N=ifthenelse(self.River, self.NRiver, self.N)
+ self.Bw = ifthenelse(self.River, RiverWidth, self.Bw)
+ # term for Alpha
+ self.AlpTerm = pow((self.N / (sqrt(self.Slope))), self.Beta)
+ # power for Alpha
+ self.AlpPow = (2.0 / 3.0) * self.Beta
+ # initial approximation for Alpha
- # 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
- # "Noah J. Finnegan et al 2005 Controls on the channel width of rivers:
- # Implications for modeling fluvial incision of bedrock"
+ # calculate catchmentsize
+ self.upsize = catchmenttotal(self.xl * self.yl, self.TopoLdd)
+ self.csize = areamaximum(self.upsize, self.TopoId)
- 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)
- # Use supplied riverwidth if possible, else calulate
- RiverWidth = ifthenelse(RiverWidth <=0.0, W, RiverWidth)
+ self.logger.info("End of initial section.")
- self.SnowWater = self.ZeroMap
-
-
- # Which columns/gauges to use/ignore in kinematic wave updating
- self.UpdateMap = self.ZeroMap
-
- if self.updating:
- _tmp =pcr2numpy(self.OutputLoc,0.0)
- gaugear= _tmp
- touse = numpy.zeros(gaugear.shape,dtype='int')
-
- for thecol in updateCols:
- idx = (gaugear == thecol).nonzero()
- touse[idx] = thecol
-
- self.UpdateMap = numpy2pcr(Nominal,touse,0.0)
- # 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)
- #self.DistToUpdPt = ldddist(self.TopoLdd,boolean(cover(self.OutputId,0.0)),1)
- #* self.reallength/celllength()
-
-
- # Initializing of variables
- self.logger.info("Initializing of model variables..")
- self.TopoLdd=lddmask(self.TopoLdd,boolean(self.TopoId))
- catchmentcells=maptotal(scalar(self.TopoId))
-
-
- # Limit lateral flow per subcatchment (make pits at all subcatch boundaries)
- # 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))
-
- # Used to seperate output per LandUse/management classes
- #OutZones = self.LandUse
- #report(self.reallength,"rl.map")
- #report(catchmentcells,"kk.map")
- self.QMMConv = self.timestepsecs/(self.reallength * self.reallength * 0.001) #m3/s --> mm
- self.ToCubic = (self.reallength * self.reallength * 0.001) / self.timestepsecs # m3/s
- 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 (weigthted for upstream area)
- self.sumlevel=self.ZeroMap #: accumulated level for water balance
- self.sumpotevap=self.ZeroMap #accumulated runoff for water balance
- self.sumsoilevap=self.ZeroMap
- self.sumtemp=self.ZeroMap #accumulated runoff for water balance
- self.ForecQ_qmec=self.ZeroMap # Extra inflow to kinematic wave reservoir for forcing in m^/sec
- self.KinWaveVolume=self.ZeroMap
- self.OldKinWaveVolume=self.ZeroMap
- self.Qvolume=self.ZeroMap
- self.Q=self.ZeroMap
- self.suminflow=self.ZeroMap
- # cntd
- self.FieldCapacity=self.FC #: total water holding capacity of the soil
- self.Treshold=self.LP*self.FieldCapacity # Threshold soilwaterstorage above which AE=PE
- #CatSurface=maptotal(scalar(ifthen(scalar(self.TopoId)>scalar(0.0),scalar(1.0)))) # catchment surface (in km2)
-
-
- 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))
-
-
-
- # Set DCL to riverlength if that is longer that the basic length calculated from grid
- drainlength = detdrainlength(self.TopoLdd,self.xl,self.yl)
-
- self.DCL=max(drainlength,self.RiverLength) # m
- # Multiply with Factor (taken from upscaling operation, defaults to 1.0 if no map is supplied
- self.DCL = self.DCL * max(1.0,self.RiverLengthFac)
-
- # water depth (m)
- # set width for kinematic wave to cell width for all cells
- self.Bw=detdrainwidth(self.TopoLdd,self.xl,self.yl)
- # However, in the main river we have real flow so set the width to the
- # width of the river
-
- self.Bw=ifthenelse(self.River, RiverWidth, self.Bw)
-
- # term for Alpha
- self.AlpTerm=pow((self.N/(sqrt(self.Slope))),self.Beta)
- # power for Alpha
- self.AlpPow=(2.0/3.0)*self.Beta
- # initial approximation for Alpha
-
- # calculate catchmentsize
- self.upsize=catchmenttotal(self.xl * self.yl,self.TopoLdd)
- self.csize=areamaximum(self.upsize,self.TopoId)
-
-
- self.logger.info("End of initial section.")
-
-
- def default_summarymaps(self):
- """
+ 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
Example:
"""
- lst = ['self.Cfmax','self.csize','self.upsize','self.TTI','self.TT','self.WHC',
- 'self.Slope','self.N','self.xl','self.yl','self.reallength','self.DCL','self.Bw',]
+ lst = [
+ "self.Cfmax",
+ "self.csize",
+ "self.upsize",
+ "self.TTI",
+ "self.TT",
+ "self.WHC",
+ "self.Slope",
+ "self.N",
+ "self.xl",
+ "self.yl",
+ "self.reallength",
+ "self.DCL",
+ "self.Bw",
+ ]
- return lst
+ return lst
- def resume(self):
- """ read initial state maps (they are output of a previous call to suspend()) """
+ def resume(self):
+ """ read initial state maps (they are output of a previous call to suspend()) """
- if self.reinit == 1:
- self.logger.info("Setting initial conditions to default (zero!)")
- self.FreeWater = cover(0.0) #: Water on surface (state variable [mm])
- self.SoilMoisture = self.FC #: Soil moisture (state variable [mm])
- self.UpperZoneStorage = 0.2 * self.FC #: Storage in Upper Zone (state variable [mm])
- self.LowerZoneStorage = 1.0/(3.0 * self.K4) #: Storage in Uppe Zone (state variable [mm])
- self.InterceptionStorage = cover(0.0) #: Interception Storage (state variable [mm])
- self.SurfaceRunoff = cover(0.0) #: Discharge in kinimatic wave (state variable [m^3/s])
- self.WaterLevel = cover(0.0) #: Water level in kinimatic wave (state variable [m])
- self.DrySnow=cover(0.0) #: Snow amount (state variable [mm])
- if hasattr(self, 'ReserVoirSimpleLocs'):
- self.ReservoirVolume = self.ResMaxVolume * self.ResTargetFullFrac
- if hasattr(self, 'ReserVoirComplexLocs'):
- self.ReservoirWaterLevel = cover(0.0)
- else:
- self.wf_resume(os.path.join(self.Dir, "instate"))
+ if self.reinit == 1:
+ self.logger.info("Setting initial conditions to default (zero!)")
+ self.FreeWater = cover(0.0) #: Water on surface (state variable [mm])
+ self.SoilMoisture = self.FC #: Soil moisture (state variable [mm])
+ self.UpperZoneStorage = (
+ 0.2 * self.FC
+ ) #: Storage in Upper Zone (state variable [mm])
+ self.LowerZoneStorage = 1.0 / (
+ 3.0 * self.K4
+ ) #: Storage in Uppe Zone (state variable [mm])
+ self.InterceptionStorage = cover(
+ 0.0
+ ) #: Interception Storage (state variable [mm])
+ self.SurfaceRunoff = cover(
+ 0.0
+ ) #: Discharge in kinimatic wave (state variable [m^3/s])
+ self.WaterLevel = cover(
+ 0.0
+ ) #: Water level in kinimatic wave (state variable [m])
+ self.DrySnow = cover(0.0) #: Snow amount (state variable [mm])
+ if hasattr(self, "ReserVoirSimpleLocs"):
+ self.ReservoirVolume = self.ResMaxVolume * self.ResTargetFullFrac
+ if hasattr(self, "ReserVoirComplexLocs"):
+ self.ReservoirWaterLevel = cover(0.0)
+ else:
+ self.wf_resume(os.path.join(self.Dir, "instate"))
- P=self.Bw+(2.0*self.WaterLevel)
- self.Alpha=self.AlpTerm*pow(P,self.AlpPow)
+ P = self.Bw + (2.0 * self.WaterLevel)
+ self.Alpha = self.AlpTerm * pow(P, self.AlpPow)
- self.OldSurfaceRunoff = self.SurfaceRunoff
+ self.OldSurfaceRunoff = self.SurfaceRunoff
- self.SurfaceRunoffMM=self.SurfaceRunoff * self.QMMConv
+ self.SurfaceRunoffMM = self.SurfaceRunoff * self.QMMConv
# Determine initial kinematic wave volume
- self.KinWaveVolume = self.WaterLevel * self.Bw * self.DCL
- self.OldKinWaveVolume = self.KinWaveVolume
- self.initstorage=self.FreeWater + self.DrySnow + self.SoilMoisture + self.UpperZoneStorage + self.LowerZoneStorage \
- + self.InterceptionStorage
+ self.KinWaveVolume = self.WaterLevel * self.Bw * self.DCL
+ self.OldKinWaveVolume = self.KinWaveVolume
+ self.initstorage = (
+ self.FreeWater
+ + self.DrySnow
+ + self.SoilMoisture
+ + self.UpperZoneStorage
+ + self.LowerZoneStorage
+ + self.InterceptionStorage
+ )
- if not self.SetKquickFlow:
- self.KQuickFlow=(self.KHQ**(1.0+self.AlphaNL))*(self.HQ**-self.AlphaNL) # recession rate of the upper reservoir, KHQ*UHQ=HQ=kquickflow*(UHQ**alpha)
+ if not self.SetKquickFlow:
+ self.KQuickFlow = (self.KHQ ** (1.0 + self.AlphaNL)) * (
+ self.HQ ** -self.AlphaNL
+ ) # recession rate of the upper reservoir, KHQ*UHQ=HQ=kquickflow*(UHQ**alpha)
+ def dynamic(self):
-
- def dynamic(self):
-
- """
+ """
Below a list of variables that can be save to disk as maps or as
timeseries (see ini file for syntax):
@@ -730,264 +1089,416 @@
:var self.ToCubic: Mutiplier to convert mm to m^3/s for fluxes
"""
- self.wf_updateparameters() # read forcing an dynamic parameters
- self.Precipitation = max(0.0,self.Precipitation) * self.Pcorr
+ self.wf_updateparameters() # read forcing an dynamic parameters
+ self.Precipitation = max(0.0, self.Precipitation) * self.Pcorr
- #self.Precipitation=cover(self.wf_readmap(self.P_mapstack,0.0),0.0) * self.Pcorr
- #self.PotEvaporation=cover(self.wf_readmap(self.PET_mapstack,0.0),0.0)
- #self.Inflow=cover(self.wf_readmap(self.Inflow_mapstack,0.0,verbose=False),0.0)
- # These ar ALWAYS 0 at present!!!
- #self.Inflow=pcrut.readmapSave(self.Inflow_mapstack,0.0)
- if self.ExternalQbase:
- self.Seepage = cover(self.wf_readmap(self.Seepage_mapstack,0.0),0.0)
- else:
- self.Seepage=cover(0.0)
- self.Temperature=cover(self.wf_readmap(self.TEMP_mapstack,10.0),10.0)
- self.Temperature = self.Temperature + self.TempCor
+ # self.Precipitation=cover(self.wf_readmap(self.P_mapstack,0.0),0.0) * self.Pcorr
+ # self.PotEvaporation=cover(self.wf_readmap(self.PET_mapstack,0.0),0.0)
+ # self.Inflow=cover(self.wf_readmap(self.Inflow_mapstack,0.0,verbose=False),0.0)
+ # These ar ALWAYS 0 at present!!!
+ # self.Inflow=pcrut.readmapSave(self.Inflow_mapstack,0.0)
+ if self.ExternalQbase:
+ self.Seepage = cover(self.wf_readmap(self.Seepage_mapstack, 0.0), 0.0)
+ else:
+ self.Seepage = cover(0.0)
+ self.Temperature = cover(self.wf_readmap(self.TEMP_mapstack, 10.0), 10.0)
+ self.Temperature = self.Temperature + self.TempCor
- # Multiply input parameters with a factor (for calibration etc) -p option in command line (no also in ini)
+ # Multiply input parameters with a factor (for calibration etc) -p option in command line (no also in ini)
- self.wf_multparameters()
+ self.wf_multparameters()
- RainFrac=ifthenelse(1.0*self.TTI == 0.0,ifthenelse(self.Temperature <= self.TT,scalar(0.0),scalar(1.0)),min((self.Temperature-(self.TT-self.TTI/2.0))/self.TTI,scalar(1.0)))
- RainFrac=max(RainFrac,scalar(0.0)) #fraction of precipitation which falls as rain
- SnowFrac=1.0-RainFrac #fraction of self.Precipitation which falls as snow
+ RainFrac = ifthenelse(
+ 1.0 * self.TTI == 0.0,
+ ifthenelse(self.Temperature <= self.TT, scalar(0.0), scalar(1.0)),
+ min(
+ (self.Temperature - (self.TT - self.TTI / 2.0)) / self.TTI, scalar(1.0)
+ ),
+ )
+ RainFrac = max(
+ RainFrac, scalar(0.0)
+ ) # fraction of precipitation which falls as rain
+ SnowFrac = 1.0 - RainFrac # fraction of self.Precipitation which falls as snow
- self.Precipitation = self.SFCF * SnowFrac * self.Precipitation + self.RFCF * RainFrac * self.Precipitation # different correction for rainfall and snowfall
+ self.Precipitation = (
+ self.SFCF * SnowFrac * self.Precipitation
+ + self.RFCF * RainFrac * self.Precipitation
+ ) # different correction for rainfall and snowfall
- #Water onto the canopy
- Interception=min(self.Precipitation,self.ICF-self.InterceptionStorage)#: Interception in mm/timestep
- self.InterceptionStorage=self.InterceptionStorage+Interception #: Current interception storage
- self.Precipitation=self.Precipitation-Interception
+ # Water onto the canopy
+ Interception = min(
+ self.Precipitation, self.ICF - self.InterceptionStorage
+ ) #: Interception in mm/timestep
+ self.InterceptionStorage = (
+ self.InterceptionStorage + Interception
+ ) #: Current interception storage
+ self.Precipitation = self.Precipitation - Interception
+ self.PotEvaporation = (
+ exp(-self.EPF * self.Precipitation) * self.ECORR * self.PotEvaporation
+ ) # correction for potential evaporation on wet days
+ self.PotEvaporation = self.CEVPF * self.PotEvaporation # Correct per landuse
- self.PotEvaporation=exp(-self.EPF*self.Precipitation)*self.ECORR * self.PotEvaporation # correction for potential evaporation on wet days
- self.PotEvaporation=self.CEVPF*self.PotEvaporation # Correct per landuse
+ self.IntEvap = min(
+ self.InterceptionStorage, self.PotEvaporation
+ ) #: Evaporation from interception storage
+ self.InterceptionStorage = self.InterceptionStorage - self.IntEvap
- self.IntEvap=min(self.InterceptionStorage,self.PotEvaporation) #: Evaporation from interception storage
- self.InterceptionStorage=self.InterceptionStorage-self.IntEvap
+ # I nthe origal HBV code
+ RestEvap = max(0.0, self.PotEvaporation - self.IntEvap)
- # I nthe origal HBV code
- RestEvap = max(0.0,self.PotEvaporation-self.IntEvap)
+ if hasattr(self, "ReserVoirComplexLocs"):
+ self.ReserVoirPotEvap = self.PotEvaporation
+ self.ReserVoirPrecip = self.Precipitation
- if hasattr(self, 'ReserVoirComplexLocs'):
- self.ReserVoirPotEvap = self.PotEvaporation
- self.ReserVoirPrecip = self.Precipitation
+ self.PotEvaporation = self.filter_P_PET * self.PotEvaporation
+ self.Precipitation = self.filter_P_PET * self.Precipitation
- self.PotEvaporation = self.filter_P_PET * self.PotEvaporation
- self.Precipitation = self.filter_P_PET * self.Precipitation
+ SnowFall = SnowFrac * self.Precipitation #: snowfall depth
+ RainFall = RainFrac * self.Precipitation #: rainfall depth
+ PotSnowMelt = ifthenelse(
+ self.Temperature > self.TT,
+ self.Cfmax * (self.Temperature - self.TT),
+ scalar(0.0),
+ ) # Potential snow melt, based on temperature
+ PotRefreezing = ifthenelse(
+ self.Temperature < self.TT,
+ self.Cfmax * self.CFR * (self.TT - self.Temperature),
+ 0.0,
+ ) # Potential refreezing, based on temperature
+ Refreezing = ifthenelse(
+ self.Temperature < self.TT, min(PotRefreezing, self.FreeWater), 0.0
+ ) # actual refreezing
+ self.SnowMelt = min(PotSnowMelt, self.DrySnow) # actual snow melt
+ self.DrySnow = (
+ self.DrySnow + SnowFall + Refreezing - self.SnowMelt
+ ) # dry snow content
+ self.FreeWater = self.FreeWater - Refreezing # free water content in snow
+ MaxFreeWater = self.DrySnow * self.WHC
+ self.FreeWater = self.FreeWater + self.SnowMelt + RainFall
+ InSoil = max(
+ self.FreeWater - MaxFreeWater, 0.0
+ ) # abundant water in snow pack which goes into soil
+ self.FreeWater = self.FreeWater - InSoil
+ RainAndSnowmelt = RainFall + self.SnowMelt
- SnowFall=SnowFrac*self.Precipitation #: snowfall depth
- RainFall=RainFrac*self.Precipitation #: rainfall depth
- PotSnowMelt=ifthenelse(self.Temperature > self.TT,self.Cfmax*(self.Temperature-self.TT),scalar(0.0)) #Potential snow melt, based on temperature
- PotRefreezing=ifthenelse(self.Temperature < self.TT, self.Cfmax*self.CFR*(self.TT-self.Temperature),0.0) #Potential refreezing, based on temperature
+ self.SnowCover = ifthenelse(self.DrySnow > 0, scalar(1), scalar(0))
+ self.NrCell = areatotal(self.SnowCover, self.TopoId)
- Refreezing=ifthenelse(self.Temperature < self.TT,min(PotRefreezing,self.FreeWater),0.0) #actual refreezing
- self.SnowMelt=min(PotSnowMelt,self.DrySnow) #actual snow melt
- self.DrySnow=self.DrySnow+SnowFall+Refreezing-self.SnowMelt #dry snow content
- self.FreeWater=self.FreeWater-Refreezing #free water content in snow
- MaxFreeWater=self.DrySnow*self.WHC
- self.FreeWater=self.FreeWater+self.SnowMelt+RainFall
- InSoil = max(self.FreeWater-MaxFreeWater,0.0) #abundant water in snow pack which goes into soil
- self.FreeWater=self.FreeWater-InSoil
- RainAndSnowmelt = RainFall + self.SnowMelt
+ # first part of precipitation is intercepted
+ # Interception=min(InSoil,self.ICF-self.InterceptionStorage)#: Interception in mm/timestep
+ # self.InterceptionStorage=self.InterceptionStorage+Interception #: Current interception storage
+ # NetInSoil=InSoil-Interception
+ NetInSoil = InSoil
- self.SnowCover = ifthenelse(self.DrySnow >0, scalar(1), scalar(0))
- self.NrCell= areatotal(self.SnowCover,self.TopoId)
+ self.SoilMoisture = self.SoilMoisture + NetInSoil
+ DirectRunoff = max(
+ self.SoilMoisture - self.FieldCapacity, 0.0
+ ) # if soil is filled to capacity: abundant water runs of directly
+ self.SoilMoisture = self.SoilMoisture - DirectRunoff
+ NetInSoil = NetInSoil - DirectRunoff # net water which infiltrates into soil
- #first part of precipitation is intercepted
- #Interception=min(InSoil,self.ICF-self.InterceptionStorage)#: Interception in mm/timestep
- #self.InterceptionStorage=self.InterceptionStorage+Interception #: Current interception storage
- #NetInSoil=InSoil-Interception
- NetInSoil=InSoil
+ MaxSnowPack = 10000.0
+ if self.MassWasting:
+ # Masswasting of snow
+ # 5.67 = tan 80 graden
+ SnowFluxFrac = min(0.5, self.Slope / 5.67) * min(
+ 1.0, self.DrySnow / MaxSnowPack
+ )
+ MaxFlux = SnowFluxFrac * self.DrySnow
+ self.DrySnow = accucapacitystate(self.TopoLdd, self.DrySnow, MaxFlux)
+ self.FreeWater = accucapacitystate(
+ self.TopoLdd, self.FreeWater, SnowFluxFrac * self.FreeWater
+ )
+ else:
+ SnowFluxFrac = self.ZeroMap
+ MaxFlux = self.ZeroMap
- self.SoilMoisture=self.SoilMoisture+NetInSoil
- DirectRunoff=max(self.SoilMoisture-self.FieldCapacity,0.0) #if soil is filled to capacity: abundant water runs of directly
- self.SoilMoisture=self.SoilMoisture-DirectRunoff
- NetInSoil=NetInSoil-DirectRunoff #net water which infiltrates into soil
+ # IntEvap=min(self.InterceptionStorage,self.PotEvaporation) #: Evaporation from interception storage
+ # self.InterceptionStorage=self.InterceptionStorage-IntEvap
- MaxSnowPack = 10000.0
- if self.MassWasting:
- # Masswasting of snow
- # 5.67 = tan 80 graden
- SnowFluxFrac = min(0.5,self.Slope/5.67) * min(1.0,self.DrySnow/MaxSnowPack)
- MaxFlux = SnowFluxFrac * self.DrySnow
- self.DrySnow = accucapacitystate(self.TopoLdd,self.DrySnow, MaxFlux)
- self.FreeWater = accucapacitystate(self.TopoLdd,self.FreeWater,SnowFluxFrac * self.FreeWater )
- else:
- SnowFluxFrac = self.ZeroMap
- MaxFlux= self.ZeroMap
+ # I nthe origal HBV code
+ # RestEvap = max(0.0,self.PotEvaporation-IntEvap)
+ self.SoilEvap = ifthenelse(
+ self.SoilMoisture > self.Treshold,
+ min(self.SoilMoisture, RestEvap),
+ min(
+ self.SoilMoisture,
+ min(
+ RestEvap, self.PotEvaporation * (self.SoilMoisture / self.Treshold)
+ ),
+ ),
+ )
+ #: soil evapotranspiration
+ self.SoilMoisture = (
+ self.SoilMoisture - self.SoilEvap
+ ) # evaporation from soil moisture storage
- #IntEvap=min(self.InterceptionStorage,self.PotEvaporation) #: Evaporation from interception storage
- #self.InterceptionStorage=self.InterceptionStorage-IntEvap
+ self.ActEvap = (
+ self.IntEvap + self.SoilEvap
+ ) #: Sum of evaporation components (IntEvap+SoilEvap)
+ self.HBVSeepage = (
+ (min(self.SoilMoisture / self.FieldCapacity, 1)) ** self.BetaSeepage
+ ) * NetInSoil # runoff water from soil
+ self.SoilMoisture = self.SoilMoisture - self.HBVSeepage
- # I nthe origal HBV code
- #RestEvap = max(0.0,self.PotEvaporation-IntEvap)
+ Backtosoil = min(
+ self.FieldCapacity - self.SoilMoisture, DirectRunoff
+ ) # correction for extremely wet periods: soil is filled to capacity
+ self.DirectRunoff = DirectRunoff - Backtosoil
+ self.SoilMoisture = self.SoilMoisture + Backtosoil
+ self.InUpperZone = (
+ self.DirectRunoff + self.HBVSeepage
+ ) # total water available for runoff
- self.SoilEvap=ifthenelse(self.SoilMoisture > self.Treshold,min(self.SoilMoisture,RestEvap),\
- min(self.SoilMoisture,min(RestEvap,self.PotEvaporation*(self.SoilMoisture/self.Treshold))))
- #: soil evapotranspiration
- self.SoilMoisture=self.SoilMoisture-self.SoilEvap #evaporation from soil moisture storage
+ # Steps is always 1 at the moment
+ # calculations for Upper zone
+ self.UpperZoneStorage = (
+ self.UpperZoneStorage + self.InUpperZone
+ ) # incoming water from soil
+ self.Percolation = min(
+ self.PERC, self.UpperZoneStorage - self.InUpperZone / 2
+ ) # Percolation
+ self.UpperZoneStorage = self.UpperZoneStorage - self.Percolation
+ self.CapFlux = self.Cflux * (
+ ((self.FieldCapacity - self.SoilMoisture) / self.FieldCapacity)
+ ) #: Capillary flux flowing back to soil
+ self.CapFlux = min(self.UpperZoneStorage, self.CapFlux)
+ self.CapFlux = min(self.FieldCapacity - self.SoilMoisture, self.CapFlux)
+ self.UpperZoneStorage = self.UpperZoneStorage - self.CapFlux
+ self.SoilMoisture = self.SoilMoisture + self.CapFlux
+ if not self.SetKquickFlow:
+ self.QuickFlow = min(
+ ifthenelse(
+ self.Percolation < self.PERC,
+ 0,
+ self.KQuickFlow
+ * (
+ (
+ self.UpperZoneStorage
+ - min(self.InUpperZone / 2, self.UpperZoneStorage)
+ )
+ ** (1.0 + self.AlphaNL)
+ ),
+ ),
+ self.UpperZoneStorage,
+ )
+ self.UpperZoneStorage = max(
+ ifthenelse(
+ self.Percolation < self.PERC,
+ self.UpperZoneStorage,
+ self.UpperZoneStorage - self.QuickFlow,
+ ),
+ 0,
+ )
+ # QuickFlow_temp = max(0,self.KQuickFlow*(self.UpperZoneStorage**(1.0+self.AlphaNL)))
+ # self.QuickFlow = min(QuickFlow_temp,self.UpperZoneStorage)
+ self.RealQuickFlow = self.ZeroMap
+ else:
+ self.QuickFlow = self.KQuickFlow * self.UpperZoneStorage
+ self.RealQuickFlow = max(0, self.K0 * (self.UpperZoneStorage - self.SUZ))
+ self.UpperZoneStorage = (
+ self.UpperZoneStorage - self.QuickFlow - self.RealQuickFlow
+ )
+ """Quickflow volume in mm/timestep"""
+ # self.UpperZoneStorage=self.UpperZoneStorage-self.QuickFlow-self.RealQuickFlow
- self.ActEvap=self.IntEvap+self.SoilEvap #: Sum of evaporation components (IntEvap+SoilEvap)
- self.HBVSeepage=((min(self.SoilMoisture/self.FieldCapacity,1))**self.BetaSeepage)*NetInSoil #runoff water from soil
- self.SoilMoisture=self.SoilMoisture-self.HBVSeepage
+ # calculations for Lower zone
+ self.LowerZoneStorage = self.LowerZoneStorage + self.Percolation
+ self.BaseFlow = min(
+ self.LowerZoneStorage, self.K4 * self.LowerZoneStorage
+ ) #: Baseflow in mm/timestep
+ self.LowerZoneStorage = self.LowerZoneStorage - self.BaseFlow
+ # Direct runoff generation
+ if self.ExternalQbase:
+ DirectRunoffStorage = self.QuickFlow + self.Seepage + self.RealQuickFlow
+ else:
+ DirectRunoffStorage = self.QuickFlow + self.BaseFlow + self.RealQuickFlow
- Backtosoil=min(self.FieldCapacity-self.SoilMoisture,DirectRunoff) #correction for extremely wet periods: soil is filled to capacity
- self.DirectRunoff=DirectRunoff-Backtosoil
- self.SoilMoisture=self.SoilMoisture+Backtosoil
- self.InUpperZone=self.DirectRunoff+self.HBVSeepage # total water available for runoff
+ self.InSoil = InSoil
+ self.RainAndSnowmelt = RainAndSnowmelt
+ self.NetInSoil = NetInSoil
+ self.InwaterMM = max(0.0, DirectRunoffStorage)
+ self.Inwater = self.InwaterMM * self.ToCubic
- # Steps is always 1 at the moment
- # calculations for Upper zone
- self.UpperZoneStorage=self.UpperZoneStorage+self.InUpperZone #incoming water from soil
- self.Percolation=min(self.PERC,self.UpperZoneStorage-self.InUpperZone/2) #Percolation
- self.UpperZoneStorage=self.UpperZoneStorage-self.Percolation
- self.CapFlux=self.Cflux*(((self.FieldCapacity-self.SoilMoisture)/self.FieldCapacity)) #: Capillary flux flowing back to soil
- self.CapFlux=min(self.UpperZoneStorage,self.CapFlux)
- self.CapFlux=min(self.FieldCapacity-self.SoilMoisture,self.CapFlux)
- self.UpperZoneStorage=self.UpperZoneStorage-self.CapFlux
- self.SoilMoisture=self.SoilMoisture+self.CapFlux
+ # only run the reservoir module if needed
- if not self.SetKquickFlow:
- self.QuickFlow=min(ifthenelse(self.Percolation 0:
+ self.ReservoirVolume, self.Outflow, self.ResPercFull, self.DemandRelease = simplereservoir(
+ self.ReservoirVolume,
+ self.SurfaceRunoff,
+ self.ResMaxVolume,
+ self.ResTargetFullFrac,
+ self.ResMaxRelease,
+ self.ResDemand,
+ self.ResTargetMinFrac,
+ self.ReserVoirSimpleLocs,
+ timestepsecs=self.timestepsecs,
+ )
+ self.OutflowDwn = upstream(
+ self.TopoLddOrg, cover(self.Outflow, scalar(0.0))
+ )
+ self.Inflow = self.OutflowDwn + cover(self.Inflow, self.ZeroMap)
+ # else:
+ # self.Inflow= cover(self.Inflow,self.ZeroMap)
- # calculations for Lower zone
- self.LowerZoneStorage=self.LowerZoneStorage+self.Percolation
- self.BaseFlow=min(self.LowerZoneStorage,self.K4*self.LowerZoneStorage) #: Baseflow in mm/timestep
- self.LowerZoneStorage=self.LowerZoneStorage-self.BaseFlow
- # Direct runoff generation
- if self.ExternalQbase:
- DirectRunoffStorage=self.QuickFlow+self.Seepage+self.RealQuickFlow
- else:
- DirectRunoffStorage=self.QuickFlow+self.BaseFlow+self.RealQuickFlow
+ elif self.nrresComplex > 0:
+ self.ReservoirWaterLevel, self.Outflow, self.ReservoirPrecipitation, self.ReservoirEvaporation, self.ReservoirVolume = 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.Outflow, scalar(0.0))
+ )
+ self.Inflow = self.OutflowDwn + cover(self.Inflow, self.ZeroMap)
+ else:
+ self.Inflow = cover(self.Inflow, self.ZeroMap)
- self.InSoil = InSoil
- self.RainAndSnowmelt = RainAndSnowmelt
- self.NetInSoil = NetInSoil
- self.InwaterMM=max(0.0,DirectRunoffStorage)
- self.Inwater=self.InwaterMM * self.ToCubic
+ self.QuickFlowCubic = (self.QuickFlow + self.RealQuickFlow) * self.ToCubic
+ self.BaseFlowCubic = self.BaseFlow * self.ToCubic
- #only run the reservoir module if needed
+ self.SurfaceWaterSupply = ifthenelse(
+ self.Inflow < 0.0,
+ max(-1.0 * self.Inwater, self.SurfaceRunoff),
+ self.ZeroMap,
+ )
+ self.Inwater = self.Inwater + ifthenelse(
+ self.SurfaceWaterSupply > 0, -1.0 * self.SurfaceWaterSupply, self.Inflow
+ )
- if self.nrresSimple > 0:
- self.ReservoirVolume, self.Outflow, self.ResPercFull,\
- self.DemandRelease = simplereservoir(self.ReservoirVolume, self.SurfaceRunoff,\
- self.ResMaxVolume, self.ResTargetFullFrac,
- self.ResMaxRelease, self.ResDemand,
- self.ResTargetMinFrac, self.ReserVoirSimpleLocs,
- timestepsecs=self.timestepsecs)
- self.OutflowDwn = upstream(self.TopoLddOrg,cover(self.Outflow,scalar(0.0)))
- self.Inflow = self.OutflowDwn + cover(self.Inflow,self.ZeroMap)
- #else:
- # self.Inflow= cover(self.Inflow,self.ZeroMap)
+ ##########################################################################
+ # Runoff calculation via Kinematic wave ##################################
+ ##########################################################################
+ # per distance along stream
+ q = self.Inwater / self.DCL + self.ForecQ_qmec / self.DCL
+ self.OldSurfaceRunoff = self.SurfaceRunoff
- elif self.nrresComplex > 0:
- self.ReservoirWaterLevel, self.Outflow, self.ReservoirPrecipitation, self.ReservoirEvaporation,\
- self.ReservoirVolume = 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.Outflow,scalar(0.0)))
- self.Inflow = self.OutflowDwn + cover(self.Inflow,self.ZeroMap)
- else:
- self.Inflow= cover(self.Inflow,self.ZeroMap)
+ self.SurfaceRunoff = kinematic(
+ self.TopoLdd,
+ self.SurfaceRunoff,
+ 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.QuickFlowCubic = (self.QuickFlow + self.RealQuickFlow) * self.ToCubic
- self.BaseFlowCubic = self.BaseFlow * self.ToCubic
+ self.updateRunOff()
+ InflowKinWaveCell = upstream(self.TopoLdd, self.SurfaceRunoff)
+ self.MassBalKinWave = (
+ (self.KinWaveVolume - self.OldKinWaveVolume) / self.timestepsecs
+ + InflowKinWaveCell
+ + self.Inwater
+ - self.SurfaceRunoff
+ )
+ Runoff = self.SurfaceRunoff
- self.SurfaceWaterSupply = ifthenelse (self.Inflow < 0.0 , max(-1.0 * self.Inwater,self.SurfaceRunoff), self.ZeroMap)
- self.Inwater = self.Inwater + ifthenelse(self.SurfaceWaterSupply> 0, -1.0 * self.SurfaceWaterSupply,self.Inflow)
+ # Updating
+ # --------
+ # Assume a tss file with as many columns as outpulocs. Start updating for each non-missing value and start with the
+ # first column (nr 1). Assumes that outputloc and columns match!
+ if self.updating:
+ QM = timeinputscalar(self.updateFile, self.UpdateMap) * self.QMMConv
- ##########################################################################
- # Runoff calculation via Kinematic wave ##################################
- ##########################################################################
- # per distance along stream
- q=self.Inwater/self.DCL + self.ForecQ_qmec/self.DCL
- self.OldSurfaceRunoff=self.SurfaceRunoff
+ # Now update the state. Just add to the Ustore
+ # self.UStoreDepth = result
+ # 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(QM, self.UpdateMap) / areamaximum(
+ self.SurfaceRunoffMM, self.UpdateMap
+ )
- self.SurfaceRunoff = kinematic(self.TopoLdd, self.SurfaceRunoff,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)
+ 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),
+ )
+ # update/nudge self.UStoreDepth for the whole upstream area,
+ # not sure how much this helps or worsens things
+ UpdSoil = True
+ if UpdSoil:
+ toadd = (self.UpperZoneStorage * UpRatioSoil) - self.UpperZoneStorage
+ self.UpperZoneStorage = self.UpperZoneStorage + toadd
- self.updateRunOff()
- InflowKinWaveCell=upstream(self.TopoLdd,self.SurfaceRunoff)
- self.MassBalKinWave = (self.KinWaveVolume - self.OldKinWaveVolume)/self.timestepsecs + InflowKinWaveCell + self.Inwater - self.SurfaceRunoff
- Runoff=self.SurfaceRunoff
+ # Update the kinematic wave reservoir up to a maximum upstream distance
+ # TODO: add (much smaller) downstream updating also?
+ MM = (1.0 - self.UpRatioKyn) / self.UpdMaxDist
+ self.UpRatioKyn = MM * self.DistToUpdPt + self.UpRatioKyn
- # Updating
- # --------
- # Assume a tss file with as many columns as outpulocs. Start updating for each non-missing value and start with the
- # first column (nr 1). Assumes that outputloc and columns match!
+ self.SurfaceRunoff = self.SurfaceRunoff * self.UpRatioKyn
+ self.SurfaceRunoffMM = (
+ self.SurfaceRunoff * self.QMMConv
+ ) # SurfaceRunoffMM (mm) from SurfaceRunoff (m3/s)
+ self.updateRunOff()
- if self.updating:
- QM = timeinputscalar(self.updateFile, self.UpdateMap) * self.QMMConv
+ Runoff = self.SurfaceRunoff
- # Now update the state. Just add to the Ustore
- # self.UStoreDepth = result
- # 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(QM, self.UpdateMap)/areamaximum(self.SurfaceRunoffMM, self.UpdateMap)
+ self.QCatchmentMM = self.SurfaceRunoff * self.QMMConvUp
+ # self.RunoffCoeff = self.QCatchmentMM/catchmenttotal(self.Precipitation, self.TopoLdd)/catchmenttotal(cover(1.0), self.TopoLdd)
- 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.sumprecip = (
+ self.sumprecip + self.Precipitation
+ ) # accumulated rainfall for water balance
+ self.sumevap = (
+ self.sumevap + self.ActEvap
+ ) # accumulated evaporation for water balance
+ self.sumsoilevap = self.sumsoilevap + self.SoilEvap
+ self.sumpotevap = self.sumpotevap + self.PotEvaporation
+ self.sumtemp = self.sumtemp + self.Temperature
+ self.sumrunoff = (
+ self.sumrunoff + self.InwaterMM
+ ) # accumulated Cell runoff for water balance
+ self.sumlevel = self.sumlevel + self.WaterLevel
+ self.suminflow = self.suminflow + self.Inflow
+ self.storage = (
+ self.FreeWater
+ + self.DrySnow
+ + self.SoilMoisture
+ + self.UpperZoneStorage
+ + self.LowerZoneStorage
+ )
+ # + self.InterceptionStorage
+ self.watbal = (
+ (self.initstorage - self.storage)
+ + self.sumprecip
+ - self.sumsoilevap
+ - self.sumrunoff
+ )
- # update/nudge self.UStoreDepth for the whole upstream area,
- # not sure how much this helps or worsens things
- UpdSoil = True
- if UpdSoil:
- toadd = (self.UpperZoneStorage * UpRatioSoil) - self.UpperZoneStorage
- self.UpperZoneStorage = self.UpperZoneStorage + toadd
- # Update the kinematic wave reservoir up to a maximum upstream distance
- # TODO: add (much smaller) downstream updating also?
- 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.updateRunOff()
-
- Runoff=self.SurfaceRunoff
-
-
- self.QCatchmentMM = self.SurfaceRunoff * self.QMMConvUp
- #self.RunoffCoeff = self.QCatchmentMM/catchmenttotal(self.Precipitation, self.TopoLdd)/catchmenttotal(cover(1.0), self.TopoLdd)
-
- self.sumprecip=self.sumprecip + self.Precipitation #accumulated rainfall for water balance
- self.sumevap=self.sumevap + self.ActEvap #accumulated evaporation for water balance
- self.sumsoilevap = self.sumsoilevap + self.SoilEvap
- self.sumpotevap=self.sumpotevap + self.PotEvaporation
- self.sumtemp=self.sumtemp + self.Temperature
- self.sumrunoff=self.sumrunoff + self.InwaterMM #accumulated Cell runoff for water balance
- self.sumlevel=self.sumlevel + self.WaterLevel
- self.suminflow=self.suminflow + self.Inflow
- self.storage=self.FreeWater + self.DrySnow + self.SoilMoisture + self.UpperZoneStorage + self.LowerZoneStorage \
- #+ self.InterceptionStorage
- self.watbal=(self.initstorage - self.storage)+self.sumprecip-self.sumsoilevap-self.sumrunoff
-
-
-
-
# The main function is used to run the program from the command line
+
def main(argv=None):
"""
Perform command line execution of the model.
@@ -996,15 +1507,15 @@
global updateCols
caseName = "default_hbv"
runId = "run_default"
- configfile="wflow_hbv.ini"
- LogFileName="wflow.log"
+ configfile = "wflow_hbv.ini"
+ LogFileName = "wflow.log"
_lastTimeStep = 0
_firstTimeStep = 0
- fewsrun=False
- runinfoFile="runinfo.xml"
- timestepsecs=86400
- wflow_cloneMap = 'wflow_subcatch.map'
- NoOverWrite=1
+ fewsrun = False
+ runinfoFile = "runinfo.xml"
+ timestepsecs = 86400
+ wflow_cloneMap = "wflow_subcatch.map"
+ NoOverWrite = 1
loglevel = logging.DEBUG
if argv is None:
@@ -1016,86 +1527,115 @@
## Main model starts here
########################################################################
try:
- opts, args = getopt.getopt(argv, 'c:QXS:F:hC:Ii:T:R:u:s:P:p:Xx:U:fl:L:')
+ opts, args = getopt.getopt(argv, "c:QXS:F:hC:Ii:T:R:u:s:P:p:Xx:U:fl:L:")
except getopt.error as msg:
pcrut.usage(msg)
for o, a in opts:
- if o == '-F':
+ if o == "-F":
runinfoFile = a
fewsrun = True
- if o == '-C': caseName = a
- if o == '-R': runId = a
- if o == '-L': LogFileName = a
- if o == '-l': exec("loglevel = logging." + a)
- if o == '-c': configfile = a
- if o == '-s': timestepsecs = int(a)
- if o == '-h': usage()
- if o == '-f': NoOverWrite = 0
+ if o == "-C":
+ caseName = a
+ if o == "-R":
+ runId = a
+ if o == "-L":
+ LogFileName = a
+ if o == "-l":
+ exec("loglevel = logging." + a)
+ if o == "-c":
+ configfile = a
+ if o == "-s":
+ timestepsecs = int(a)
+ if o == "-h":
+ usage()
+ if o == "-f":
+ NoOverWrite = 0
-
-
if fewsrun:
- ts = getTimeStepsfromRuninfo(runinfoFile,timestepsecs)
+ ts = getTimeStepsfromRuninfo(runinfoFile, timestepsecs)
starttime = getStartTimefromRuninfo(runinfoFile)
- if (ts):
- _lastTimeStep = ts# * 86400/timestepsecs
+ if ts:
+ _lastTimeStep = ts # * 86400/timestepsecs
_firstTimeStep = 1
else:
print("Failed to get timesteps from runinfo file: " + runinfoFile)
sys.exit(2)
else:
- starttime = dt.datetime(1990,1,1)
+ starttime = dt.datetime(1990, 1, 1)
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,logfname=LogFileName,level=loglevel,doSetupFramework=False)
+ myModel = WflowModel(wflow_cloneMap, caseName, runId, configfile)
+ dynModelFw = wf_DynamicFramework(
+ myModel, _lastTimeStep, firstTimestep=_firstTimeStep, datetimestart=starttime
+ )
+ dynModelFw.createRunId(
+ NoOverWrite=NoOverWrite,
+ logfname=LogFileName,
+ level=loglevel,
+ 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':
- exec("zz =" + a)
+ 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":
+ exec("zz =" + a)
updateCols = zz
- if o == '-T':
- configset(myModel.config, 'run', 'endtime', a, overwrite=True)
- if o == '-S':
- configset(myModel.config, 'run', 'starttime', a, 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.logger.info("Command line: " + str(argv))
dynModelFw._runInitial()
dynModelFw._runResume()
- #dynModelFw._runDynamic(0,0)
+ # dynModelFw._runDynamic(0,0)
dynModelFw._runDynamic(_firstTimeStep, _lastTimeStep)
dynModelFw._runSuspend()
dynModelFw._wf_shutdown()
-
os.chdir("../../")