Index: wflow-py/wflow/wflow_sbm.py =================================================================== diff -u -r9dd1a78b1a2e6aa0e67e970d6235e4923b8bfc42 -r21a737cc033ad3a927e53eb8d3cd820574c155e9 --- wflow-py/wflow/wflow_sbm.py (.../wflow_sbm.py) (revision 9dd1a78b1a2e6aa0e67e970d6235e4923b8bfc42) +++ wflow-py/wflow/wflow_sbm.py (.../wflow_sbm.py) (revision 21a737cc033ad3a927e53eb8d3cd820574c155e9) @@ -100,6 +100,7 @@ import os.path import shutil, glob import getopt + try: from wflow.wf_DynamicFramework import * except ImportError: @@ -108,12 +109,12 @@ from wflow.wflow_funcs import * except ImportError: from wflow_funcs import * - + try: - from wflow.wflow_adapt import * + from wflow.wflow_adapt import * except ImportError: - from wflow_adapt import * - + from wflow_adapt import * + import scipy import ConfigParser @@ -127,15 +128,16 @@ # Dictionary with parameters and multipliers (used in calibration) multpars = {} multdynapars = {} + + def usage(*args): sys.stdout = sys.stderr for msg in args: print msg print __doc__ sys.exit(0) - -def actEvap_SBM(RootingDepth,WTable, UStoreDepth,FirstZoneDepth, PotTrans,smoothpar): +def actEvap_SBM(RootingDepth, WTable, UStoreDepth, FirstZoneDepth, PotTrans, smoothpar): """ Actual evaporation function: @@ -153,33 +155,33 @@ - ActEvap, FirstZoneDepth, UStoreDepth ActEvapUStore """ - + # Step 1 from saturated zone, use rootingDepth as a limiting factor #rootsinWater = WTable < RootingDepth #ActEvapSat = ifthenelse(rootsinWater,min(PotTrans,FirstZoneDepth),0.0) # new method: # use sCurve to determine if the roots are wet.At the moment this ise set # to be a 0-1 curve - wetroots = sCurve(WTable,a=RootingDepth,c=smoothpar) - ActEvapSat = min(PotTrans * wetroots,FirstZoneDepth) - + wetroots = sCurve(WTable, a=RootingDepth, c=smoothpar) + ActEvapSat = min(PotTrans * wetroots, FirstZoneDepth) + FirstZoneDepth = FirstZoneDepth - ActEvapSat RestPotEvap = PotTrans - ActEvapSat - + # now try unsat store - AvailCap = min(1.0,max (0.0,(WTable - RootingDepth)/(RootingDepth + 1.0))) - + AvailCap = min(1.0, max(0.0, (WTable - RootingDepth) / (RootingDepth + 1.0))) + #AvailCap = max(0.0,ifthenelse(WTable < RootingDepth, WTable/RootingDepth, RootingDepth/WTable)) - MaxExtr = AvailCap * UStoreDepth - ActEvapUStore = min(MaxExtr,RestPotEvap,UStoreDepth) + MaxExtr = AvailCap * UStoreDepth + ActEvapUStore = min(MaxExtr, RestPotEvap, UStoreDepth) UStoreDepth = UStoreDepth - ActEvapUStore - + ActEvap = ActEvapSat + ActEvapUStore - return ActEvap, FirstZoneDepth, UStoreDepth, ActEvapUStore + return ActEvap, FirstZoneDepth, UStoreDepth, ActEvapUStore -def SnowPackHBV (Snow,SnowWater,Precipitation,Temperature,TTI,TT,Cfmax,WHC): +def SnowPackHBV(Snow, SnowWater, Precipitation, Temperature, TTI, TT, Cfmax, WHC): """ HBV Type snowpack modelling using a Temperature degree factor. All correction factors (RFCF and SFCF) are set to 1. The refreezing efficiency factor is set to 0.05. @@ -191,39 +193,38 @@ :returns: Snow,SnowMelt,Precipitation """ - - RFCF= 1.0 # correction factor for rainfall - 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 - - SnowFall=SnowFrac*Precipitation #snowfall depth - RainFall=RainFrac*Precipitation #rainfall depth - PotSnowMelt=ifthenelse(Temperature > TT,Cfmax*(Temperature-TT),scalar(0.0)) #Potential snow melt, based on temperature - PotRefreezing=ifthenelse(Temperature < TT, Cfmax*CFR*(TT-Temperature),0.0) #Potential refreezing, based on temperature - Refreezing=ifthenelse(Temperature < TT,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 - MaxSnowWater=Snow*WHC # Max water in the snow - 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 + RFCF = 1.0 # correction factor for rainfall + 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 + SnowFall = SnowFrac * Precipitation #snowfall depth + RainFall = RainFrac * Precipitation #rainfall depth + PotSnowMelt = ifthenelse(Temperature > TT, Cfmax * (Temperature - TT), + scalar(0.0)) #Potential snow melt, based on temperature + PotRefreezing = ifthenelse(Temperature < TT, Cfmax * CFR * (TT - Temperature), + 0.0) #Potential refreezing, based on temperature + Refreezing = ifthenelse(Temperature < TT, 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 + MaxSnowWater = Snow * WHC # Max water in the snow + 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 + class WflowModel(DynamicModel): - """ + """ .. versionchanged:: 0.91 - Calculation of GWScale moved to resume() to allow fitting. @@ -232,32 +233,29 @@ """ - def __init__(self, cloneMap,Dir,RunDir,configfile): - DynamicModel.__init__(self) - setclone(Dir + "/staticmaps/" + cloneMap) - self.runId=RunDir - self.caseName=Dir - self.Dir = Dir + "/" - self.configfile = configfile + def __init__(self, cloneMap, Dir, RunDir, configfile): + DynamicModel.__init__(self) + setclone(Dir + "/staticmaps/" + cloneMap) + self.runId = RunDir + self.caseName = Dir + self.Dir = Dir + "/" + self.configfile = configfile - - - - def updateRunOff(self): - """ + def updateRunOff(self): + """ Updates the kinematic wave reservoir. Should be run after updates to Q """ - 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. @@ -274,50 +272,49 @@ :var self.FirstZoneDepth: Water in the saturated store [mm] :var self.CanopyStorage: Amount of water on the Canopy [mm] """ - states = ['SurfaceRunoff', 'WaterLevel', - 'FirstZoneDepth', - 'Snow', - 'TSoil', + states = ['SurfaceRunoff', 'WaterLevel', + 'FirstZoneDepth', + 'Snow', + 'TSoil', 'UStoreDepth', - 'SnowWater','CanopyStorage'] - - return states - + 'SnowWater', 'CanopyStorage'] - - def supplyCurrentTime(self): - """ + return states + + + def supplyCurrentTime(self): + """ gets the current time in seconds after the start of the run """ - return self.currentTimeStep() * self.timestepsecs - - - def suspend(self): - - self.logger.info("Saving initial conditions...") - self.wf_suspend(self.SaveDir + "/outstate/") - - if self.OverWriteInit: + return self.currentTimeStep() * self.timestepsecs + + + def suspend(self): + + self.logger.info("Saving initial conditions...") + self.wf_suspend(self.SaveDir + "/outstate/") + + if self.OverWriteInit: self.logger.info("Saving initial conditions over start conditions...") self.wf_suspend(self.SaveDir + "/instate/") - if self.fewsrun: - self.logger.info("Saving initial conditions for FEWS...") - self.wf_suspend(self.Dir + "/outstate/") + if self.fewsrun: + self.logger.info("Saving initial conditions for FEWS...") + self.wf_suspend(self.Dir + "/outstate/") - report(self.CumInwaterMM,self.SaveDir + "/outsum/CumInwaterMM.map") - report(self.CumReinfilt,self.SaveDir + "/outsum/CumReinfilt.map") - report(self.CumPrec,self.SaveDir + "/outsum/CumPrec.map") - report(self.CumEvap,self.SaveDir + "/outsum/CumEvap.map") - report(self.CumPotenTrans,self.SaveDir + "/outsum/CumPotenTrans.map") - report(self.CumInt,self.SaveDir + "/outsum/CumInt.map") - report(self.CumLeakage,self.SaveDir + "/outsum/CumLeakage.map") - report(self.CumPotenEvap,self.SaveDir + "/outsum/CumPotenEvap.map") - report(self.CumExfiltWater,self.SaveDir + "/outsum/CumExfiltWater.map") - report(self.watbal,self.SaveDir + "/outsum/watbal.map") - - def initial(self): - """ + report(self.CumInwaterMM, self.SaveDir + "/outsum/CumInwaterMM.map") + report(self.CumReinfilt, self.SaveDir + "/outsum/CumReinfilt.map") + report(self.CumPrec, self.SaveDir + "/outsum/CumPrec.map") + report(self.CumEvap, self.SaveDir + "/outsum/CumEvap.map") + report(self.CumPotenTrans, self.SaveDir + "/outsum/CumPotenTrans.map") + report(self.CumInt, self.SaveDir + "/outsum/CumInt.map") + report(self.CumLeakage, self.SaveDir + "/outsum/CumLeakage.map") + report(self.CumPotenEvap, self.SaveDir + "/outsum/CumPotenEvap.map") + report(self.CumExfiltWater, self.SaveDir + "/outsum/CumExfiltWater.map") + report(self.watbal, self.SaveDir + "/outsum/watbal.map") + + def initial(self): + """ Initial part of the model, executed only once. Reads all static data from disk @@ -365,426 +362,456 @@ :var w_soil.tbl: Soil temperature smooth factor. Given for daily timesteps. (0.1125) [-] Wigmosta, M. S., L. J. Lane, J. D. Tagestad, and A. M. Coleman (2009). """ - global statistics - global multpars - global updateCols + global statistics + global multpars + global updateCols - self.thestep = scalar(0) - self.basetimestep = 86400 - self.SSSF=False - setglobaloption("unittrue") + self.thestep = scalar(0) + self.basetimestep = 86400 + self.SSSF = False + setglobaloption("unittrue") - - self.precipTss="/intss/P.tss" - self.evapTss="/intss/PET.tss" - self.tempTss="/intss/T.tss" - self.inflowTss="/intss/Inflow.tss" + self.precipTss = "/intss/P.tss" + self.evapTss = "/intss/PET.tss" + self.tempTss = "/intss/T.tss" + self.inflowTss = "/intss/Inflow.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,"model","reinit","0")) - self.fewsrun = int(configget(self.config,"model","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.timestepsecs = int(configget(self.config,"model","timestepsecs","86400")) - 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 - 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')) - - + # 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, "model", "reinit", "0")) + self.fewsrun = int(configget(self.config, "model", "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") - # 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(readmap(self.Dir + wflow_subcatch)) # Determines the area of calculations (all cells > 0) - subcatch = ifthen(subcatch > 0, subcatch) + self.sCatch = int(configget(self.config, "model", "sCatch", "0")) + self.intbl = configget(self.config, "model", "intbl", "intbl") + self.timestepsecs = int(configget(self.config, "model", "timestepsecs", "86400")) + 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 + Qmax = float(configget(self.config, "model", "AnnualDischarge", "300")) + self.UpdMaxDist = float(configget(self.config, "model", "UpdMaxDist", "100")) - - self.Altitude=readmap(self.Dir + wflow_dem)# * scalar(defined(subcatch)) # DEM - self.TopoLdd=readmap(self.Dir + wflow_ldd) # Local - self.TopoId=readmap(self.Dir + wflow_subcatch) # area map - self.River=cover(boolean(readmap(self.Dir + wflow_river)),0) - - self.RiverLength=cover(pcrut.readmapSave(self.Dir + wflow_riverlength,0.0),0.0) - # Factor to multiply riverlength with (defaults to 1.0) - self.RiverLengthFac=pcrut.readmapSave(self.Dir + wflow_riverlength_fact,1.0) + 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")) - # 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=readmap(self.Dir + wflow_landuse) - self.LandUse=cover(self.LandUse,nominal(ordinal(subcatch) > 0)) - self.Soil=readmap(self.Dir + wflow_soil) - self.Soil=cover(self.Soil,nominal(ordinal(subcatch) > 0)) - self.OutputLoc=readmap(self.Dir + wflow_gauges) # location of output gauge(s) - self.InflowLoc=pcrut.readmapSave(self.Dir + wflow_inflow,0.0) # location abstractions/inflows. - self.RiverWidth=pcrut.readmapSave(self.Dir + wflow_riverwidth,0.0) - - # Experimental - self.RunoffGenSigmaFunction = int(configget(self.config,'model','RunoffGenSigmaFunction','0')) - self.RunoffGeneratingGWPerc = float(configget(self.config,'defaultfortbl','RunoffGeneratingGWPerc','0.1')) - - if self.scalarInput: - self.gaugesMap=readmap(self.Dir + wflow_mgauges) # location of rainfall/evap/temp gauge(s) - self.OutputId=readmap(self.Dir + wflow_subcatch) # location of subcatchment - # Temperature correction poer cell to add + #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.TempCor=pcrut.readmapSave(self.Dir + configget(self.config,"model","TemperatureCorrectionMap","staticmaps/wflow_tempcor.map"),0.0) - - self.ZeroMap=0.0*scalar(subcatch) #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) - # Set static initial values here ######################################### - self.SoilAlbedo = 0.1 # Not used at the moment - self.pi = 3.1416 - self.e = 2.7183 - self.SScale = 100.0 - - self.Latitude = ycoordinate(boolean(self.Altitude)) - self.Longitude = xcoordinate(boolean(self.Altitude)) - - self.logger.info("Linking parameters to landuse, catchment and soil...") - self.RunoffGeneratingGWPerc=self.readtblDefault(self.Dir + "/" + self.intbl + "/RunoffGeneratingGWPerc.tbl",self.LandUse,subcatch,self.Soil,self.RunoffGeneratingGWPerc) - - self.Cmax=self.readtblDefault(self.Dir + "/" + self.intbl + "/MaxCanopyStorage.tbl",self.LandUse,subcatch,self.Soil,1.0) - self.EoverR=self.readtblDefault(self.Dir + "/" + self.intbl + "/EoverR.tbl",self.LandUse,subcatch,self.Soil,0.1) - # self.Albedo=lookupscalar(self.Dir + "\intbl\Albedo.tbl",self.LandUse) # Not used anymore - self.CanopyGapFraction=self.readtblDefault(self.Dir + "/" + self.intbl + "/CanopyGapFraction.tbl",self.LandUse,subcatch,self.Soil,0.1) - 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 + # 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") - # 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) # - # 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 - # 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) - # thickness of the soil - self.FirstZoneThickness = self.readtblDefault(self.Dir + "/" + self.intbl + "/FirstZoneCapacity.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.FirstZoneMinCapacity = self.readtblDefault(self.Dir + "/" + self.intbl + "/FirstZoneMinCapacity.tbl",self.LandUse,subcatch,self.Soil,500.0) + # 2: Input base maps ######################################################## + subcatch = ordinal(readmap(self.Dir + wflow_subcatch)) # Determines the area of calculations (all cells > 0) + subcatch = ifthen(subcatch > 0, subcatch) - # FirstZoneKsatVer = $2\inmaps\FirstZoneKsatVer.map - self.FirstZoneKsatVer=self.readtblDefault(self.Dir + "/" + self.intbl + "/FirstZoneKsatVer.tbl",self.LandUse,subcatch,self.Soil,3000.0) * self.timestepsecs/self.basetimestep - 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.NRiver=self.readtblDefault(self.Dir + "/" + self.intbl + "/N_River.tbl",self.LandUse,subcatch,self.Soil,0.036) # Manning river - self.WaterFrac=self.readtblDefault(self.Dir + "/" + self.intbl + "/WaterFrac.tbl",self.LandUse,subcatch,self.Soil,0.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) - # 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) - # 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.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 - - # 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)) - self.Slope=max(0.001,self.Slope*celllength()/self.reallength) - Terrain_angle=scalar(atan(self.Slope)) + self.Altitude = readmap(self.Dir + wflow_dem) # * scalar(defined(subcatch)) # DEM + self.TopoLdd = readmap(self.Dir + wflow_ldd) # Local + self.TopoId = readmap(self.Dir + wflow_subcatch) # area map + self.River = cover(boolean(readmap(self.Dir + wflow_river)), 0) - - # Multiply parameters with a factor (for calibration etc) -P option in command line - for k, v in multpars.iteritems(): - if self.sCatch > 0: - estr = k + "= ifthenelse(self.TopoId == self.sCatch," + k + "*" + str(v)+ "," + k + ")" - else: - estr = k + "=" + k + "*" + str(v) - self.logger.info("Parameter multiplication: " + estr) - exec estr - - self.N=ifthenelse(self.River, self.NRiver, self.N) - + self.RiverLength = cover(pcrut.readmapSave(self.Dir + wflow_riverlength, 0.0), 0.0) + # Factor to multiply riverlength with (defaults to 1.0) + self.RiverLengthFac = pcrut.readmapSave(self.Dir + wflow_riverlength_fact, 1.0) - # 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 - self.RiverWidth = ifthenelse(self.RiverWidth <=0.0, W, self.RiverWidth) - - + # 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 = readmap(self.Dir + wflow_landuse) + self.LandUse = cover(self.LandUse, nominal(ordinal(subcatch) > 0)) + self.Soil = readmap(self.Dir + wflow_soil) + self.Soil = cover(self.Soil, nominal(ordinal(subcatch) > 0)) + self.OutputLoc = readmap(self.Dir + wflow_gauges) # location of output gauge(s) + self.InflowLoc = pcrut.readmapSave(self.Dir + wflow_inflow, 0.0) # location abstractions/inflows. + self.RiverWidth = pcrut.readmapSave(self.Dir + wflow_riverwidth, 0.0) - # 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 - WIMax = areamaximum(WI, self.TopoId) * WIMaxScale - self.FirstZoneThickness = max(min(self.FirstZoneThickness,(WI/WIMax) * self.FirstZoneThickness), self.FirstZoneMinCapacity) - - self.FirstZoneCapacity = self.FirstZoneThickness * (self.thetaS -self.thetaR) - - # limit roots to top 99% of first zone - self.RootingDepth = min(self.FirstZoneThickness * 0.99,self.RootingDepth) + # Experimental + self.RunoffGenSigmaFunction = int(configget(self.config, 'model', 'RunoffGenSigmaFunction', '0')) + self.RunoffGeneratingGWPerc = float(configget(self.config, 'defaultfortbl', 'RunoffGeneratingGWPerc', '0.1')) - # subgrid runoff generation - self.DemMax=readmap(self.Dir + "/staticmaps/wflow_demmax") - self.DrainageBase=readmap(self.Dir + "/staticmaps/wflow_demmin") - self.CC = min(100.0,-log(1.0/0.1 - 1)/min(-0.1,self.DrainageBase - self.Altitude)) - - - #self.GWScale = (self.DemMax-self.DrainageBase)/self.FirstZoneThickness / self.RunoffGeneratingGWPerc - # 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') + if self.scalarInput: + self.gaugesMap = readmap(self.Dir + wflow_mgauges) # location of rainfall/evap/temp gauge(s) + self.OutputId = readmap(self.Dir + wflow_subcatch) # location of subcatchment + # Temperature correction poer cell to add - for thecol in updateCols: - idx = (gaugear == thecol).nonzero() - touse[idx] = thecol - - self.UpdateMap = numpy2pcr(Nominal,touse,0.0) - # Calulate 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)) - - # Used to seperate output per LandUse/management classes - OutZones = self.LandUse - - 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.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.sumleakage=self.ZeroMap - self.CumReinfilt=self.ZeroMap - self.sumoutflow=self.ZeroMap - self.sumsnowmelt=self.ZeroMap - self.CumRad=self.ZeroMap - self.SnowMelt=self.ZeroMap - self.CumPrec=self.ZeroMap - self.CumInwaterMM=self.ZeroMap - self.CumInfiltExcess=self.ZeroMap - self.CumExfiltWater=self.ZeroMap - self.CumSurfaceWater=self.ZeroMap - self.CumEvap=self.ZeroMap - self.CumPotenEvap=self.ZeroMap - self.CumPotenTrans=self.ZeroMap - self.CumInt=self.ZeroMap - self.CumRad=self.ZeroMap - self.CumLeakage=self.ZeroMap - self.CumPrecPol=self.ZeroMap - self.FirstZoneFlux=self.ZeroMap - self.FreeWaterDepth=self.ZeroMap - self.SumCellWatBal=self.ZeroMap - self.PathInfiltExceeded=self.ZeroMap - self.SoilInfiltExceeded=self.ZeroMap - self.CumOutFlow=self.ZeroMap - self.CumCellInFlow=self.ZeroMap - self.CumIF=self.ZeroMap - self.CumActInfilt=self.ZeroMap - 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) - - # Multiply with Factor (taken from upscaling operation, defaults to 1.0 if no map is supplied - self.DCL = drainlength * max(1.0,self.RiverLengthFac) - - self.DCL=max(self.DCL,self.RiverLength) # m - - # 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, self.RiverWidth, self.Bw) - - # Add rivers to the WaterFrac, but check with waterfrac map - self.RiverFrac = min(1.0,ifthenelse(self.River,(self.RiverWidth*self.DCL)/(self.xl*self.yl),0)) - self.WaterFrac = self.WaterFrac - ifthenelse((self.RiverFrac + self.WaterFrac) > 1.0, self.RiverFrac + self.WaterFrac - 1.0, 0.0) - - - # 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") - self.AlpTerm=pow((self.N/(sqrt(self.Slope* riverslopecor))),self.Beta) - # power for Alpha - self.AlpPow=(2.0/3.0)*self.Beta - # initial approximation for Alpha - - #self.initstorage=areaaverage(self.FirstZoneDepth,self.TopoId)+areaaverage(self.UStoreDepth,self.TopoId)#+areaaverage(self.Snow,self.TopoId) - - # calculate catchmentsize - self.upsize=catchmenttotal(self.xl * self.yl,self.TopoLdd) - self.csize=areamaximum(self.upsize,self.TopoId) - # Save some summary maps - self.logger.info("Saving summary maps...") - if self.modelSnow: - report(self.Cfmax,self.Dir + "/" + self.runId + "/outsum/Cfmax.map") - report(self.TTI,self.Dir + "/" + self.runId + "/outsum/TTI.map") - report(self.TT,self.Dir + "/" + self.runId + "/outsum/TT.map") - report(self.WHC,self.Dir + "/" + self.runId + "/outsum/WHC.map") - - report(self.Cmax,self.Dir + "/" + self.runId + "/outsum/Cmax.map") - report(self.csize,self.Dir + "/" + self.runId + "/outsum/CatchmentSize.map") - report(self.upsize,self.Dir + "/" + self.runId + "/outsum/UpstreamSize.map") - report(self.EoverR,self.Dir + "/" + self.runId + "/outsum/EoverR.map") - report(self.RootingDepth,self.Dir + "/" + self.runId + "/outsum/RootingDepth.map") - report(self.CanopyGapFraction,self.Dir + "/" + self.runId + "/outsum/CanopyGapFraction.map") - report(self.InfiltCapSoil,self.Dir + "/" + self.runId + "/outsum/InfiltCapSoil.map") - report(self.InfiltCapPath,self.Dir + "/" + self.runId + "/outsum/InfiltCapPath.map") - report(self.PathFrac,self.Dir + "/" + self.runId + "/outsum/PathFrac.map") - report(self.thetaR,self.Dir + "/" + self.runId + "/outsum/thetaR.map") - report(self.thetaS,self.Dir + "/" + self.runId + "/outsum/thetaS.map") - report(self.FirstZoneMinCapacity,self.Dir + "/" + self.runId + "/outsum/FirstZoneMinCapacity.map") - report(self.FirstZoneKsatVer,self.Dir + "/" + self.runId + "/outsum/FirstZoneKsatVer.map") - report(self.M,self.Dir + "/" + self.runId + "/outsum/M.map") - report(self.FirstZoneCapacity,self.Dir + "/" + self.runId + "/outsum/FirstZoneCapacity.map") - report(Terrain_angle,self.Dir + "/" + self.runId + "/outsum/angle.map") - report(self.Slope,self.Dir + "/" + self.runId + "/outsum/slope.map") - report(WI,self.Dir + "/" + self.runId + "/outsum/WI.map") - report(self.CC,self.Dir + "/" + self.runId + "/outsum/CC.map") - report(self.N,self.Dir + "/" + self.runId + "/outsum/N.map") - report(self.RiverFrac,self.Dir + "/" + self.runId + "/outsum/RiverFrac.map") - - report(self.xl,self.Dir + "/" + self.runId + "/outsum/xl.map") - report(self.yl,self.Dir + "/" + self.runId + "/outsum/yl.map") - report(self.reallength,self.Dir + "/" + self.runId + "/outsum/rl.map") - report(self.DCL,self.Dir + "/" + self.runId + "/outsum/DCL.map") - report(self.Bw,self.Dir + "/" + self.runId + "/outsum/Bw.map") - report(ifthen(self.River,self.Bw),self.Dir + "/" + self.runId + "/outsum/RiverWidth.map") - if self.updating: - report(self.DistToUpdPt,self.Dir + "/" + self.runId + "/outsum/DistToUpdPt.map") - + self.TempCor = pcrut.readmapSave( + self.Dir + configget(self.config, "model", "TemperatureCorrectionMap", "staticmaps/wflow_tempcor.map"), 0.0) - self.SaveDir = self.Dir + "/" + self.runId + "/" - self.logger.info("End of initial section") + self.ZeroMap = 0.0 * scalar(subcatch) #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) - def resume(self): - if self.reinit == 1: - self.logger.info("Setting initial conditions to default") - self.FirstZoneDepth = self.FirstZoneCapacity * 0.85 - self.UStoreDepth = self.FirstZoneCapacity * 0.0 - self.WaterLevel = self.ZeroMap - self.SurfaceRunoff = self.ZeroMap - self.Snow = self.ZeroMap - self.SnowWater = self.ZeroMap - self.TSoil = self.ZeroMap + 10.0 - self.CanopyStorage = self.ZeroMap - - else: - self.logger.info("Setting initial conditions from state files") - self.wf_resume(self.Dir + "/instate/") - + # Set static initial values here ######################################### + self.SoilAlbedo = 0.1 # Not used at the moment + self.pi = 3.1416 + self.e = 2.7183 + self.SScale = 100.0 - - P=self.Bw+(2.0*self.WaterLevel) - self.Alpha=self.AlpTerm*pow(P,self.AlpPow) - self.OldSurfaceRunoff = self.SurfaceRunoff - - self.SurfaceRunoffMM=self.SurfaceRunoff * self.QMMConv + self.Latitude = ycoordinate(boolean(self.Altitude)) + self.Longitude = xcoordinate(boolean(self.Altitude)) + + self.logger.info("Linking parameters to landuse, catchment and soil...") + self.RunoffGeneratingGWPerc = self.readtblDefault(self.Dir + "/" + self.intbl + "/RunoffGeneratingGWPerc.tbl", + self.LandUse, subcatch, self.Soil, + self.RunoffGeneratingGWPerc) + + self.Cmax = self.readtblDefault(self.Dir + "/" + self.intbl + "/MaxCanopyStorage.tbl", self.LandUse, subcatch, + self.Soil, 1.0) + self.EoverR = self.readtblDefault(self.Dir + "/" + self.intbl + "/EoverR.tbl", self.LandUse, subcatch, + self.Soil, 0.1) + # self.Albedo=lookupscalar(self.Dir + "\intbl\Albedo.tbl",self.LandUse) # Not used anymore + self.CanopyGapFraction = self.readtblDefault(self.Dir + "/" + self.intbl + "/CanopyGapFraction.tbl", + self.LandUse, subcatch, self.Soil, 0.1) + 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 + + # 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) # + # 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 + # 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) + # thickness of the soil + self.FirstZoneThickness = self.readtblDefault(self.Dir + "/" + self.intbl + "/FirstZoneCapacity.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.FirstZoneMinCapacity = self.readtblDefault(self.Dir + "/" + self.intbl + "/FirstZoneMinCapacity.tbl", + self.LandUse, subcatch, self.Soil, 500.0) + + # FirstZoneKsatVer = $2\inmaps\FirstZoneKsatVer.map + self.FirstZoneKsatVer = self.readtblDefault(self.Dir + "/" + self.intbl + "/FirstZoneKsatVer.tbl", self.LandUse, + subcatch, self.Soil, 3000.0) * self.timestepsecs / self.basetimestep + 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.NRiver = self.readtblDefault(self.Dir + "/" + self.intbl + "/N_River.tbl", self.LandUse, subcatch, + self.Soil, 0.036) # Manning river + self.WaterFrac = self.readtblDefault(self.Dir + "/" + self.intbl + "/WaterFrac.tbl", self.LandUse, subcatch, + self.Soil, 0.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) + # 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) + # 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.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 + + # 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)) + self.Slope = max(0.001, self.Slope * celllength() / self.reallength) + Terrain_angle = scalar(atan(self.Slope)) + + + # Multiply parameters with a factor (for calibration etc) -P option in command line + for k, v in multpars.iteritems(): + if self.sCatch > 0: + estr = k + "= ifthenelse(self.TopoId == self.sCatch," + k + "*" + str(v) + "," + k + ")" + else: + estr = k + "=" + k + "*" + str(v) + self.logger.info("Parameter multiplication: " + estr) + exec estr + + 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 + self.RiverWidth = ifthenelse(self.RiverWidth <= 0.0, W, self.RiverWidth) + + + + # 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 + WIMax = areamaximum(WI, self.TopoId) * WIMaxScale + self.FirstZoneThickness = max(min(self.FirstZoneThickness, (WI / WIMax) * self.FirstZoneThickness), + self.FirstZoneMinCapacity) + + self.FirstZoneCapacity = self.FirstZoneThickness * (self.thetaS - self.thetaR) + + # limit roots to top 99% of first zone + self.RootingDepth = min(self.FirstZoneThickness * 0.99, self.RootingDepth) + + # subgrid runoff generation + self.DemMax = readmap(self.Dir + "/staticmaps/wflow_demmax") + self.DrainageBase = readmap(self.Dir + "/staticmaps/wflow_demmin") + self.CC = min(100.0, -log(1.0 / 0.1 - 1) / min(-0.1, self.DrainageBase - self.Altitude)) + + + #self.GWScale = (self.DemMax-self.DrainageBase)/self.FirstZoneThickness / self.RunoffGeneratingGWPerc + # 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') + + 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)) + + # Used to seperate output per LandUse/management classes + OutZones = self.LandUse + + 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.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.sumleakage = self.ZeroMap + self.CumReinfilt = self.ZeroMap + self.sumoutflow = self.ZeroMap + self.sumsnowmelt = self.ZeroMap + self.CumRad = self.ZeroMap + self.SnowMelt = self.ZeroMap + self.CumPrec = self.ZeroMap + self.CumInwaterMM = self.ZeroMap + self.CumInfiltExcess = self.ZeroMap + self.CumExfiltWater = self.ZeroMap + self.CumSurfaceWater = self.ZeroMap + self.CumEvap = self.ZeroMap + self.CumPotenEvap = self.ZeroMap + self.CumPotenTrans = self.ZeroMap + self.CumInt = self.ZeroMap + self.CumRad = self.ZeroMap + self.CumLeakage = self.ZeroMap + self.CumPrecPol = self.ZeroMap + self.FirstZoneFlux = self.ZeroMap + self.FreeWaterDepth = self.ZeroMap + self.SumCellWatBal = self.ZeroMap + self.PathInfiltExceeded = self.ZeroMap + self.SoilInfiltExceeded = self.ZeroMap + self.CumOutFlow = self.ZeroMap + self.CumCellInFlow = self.ZeroMap + self.CumIF = self.ZeroMap + self.CumActInfilt = self.ZeroMap + 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) + + # Multiply with Factor (taken from upscaling operation, defaults to 1.0 if no map is supplied + self.DCL = drainlength * max(1.0, self.RiverLengthFac) + + self.DCL = max(self.DCL, self.RiverLength) # m + + # 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, self.RiverWidth, self.Bw) + + # Add rivers to the WaterFrac, but check with waterfrac map + self.RiverFrac = min(1.0, ifthenelse(self.River, (self.RiverWidth * self.DCL) / (self.xl * self.yl), 0)) + self.WaterFrac = self.WaterFrac - ifthenelse((self.RiverFrac + self.WaterFrac) > 1.0, + self.RiverFrac + self.WaterFrac - 1.0, 0.0) + + + # 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") + self.AlpTerm = pow((self.N / (sqrt(self.Slope * riverslopecor))), self.Beta) + # power for Alpha + self.AlpPow = (2.0 / 3.0) * self.Beta + # initial approximation for Alpha + + #self.initstorage=areaaverage(self.FirstZoneDepth,self.TopoId)+areaaverage(self.UStoreDepth,self.TopoId)#+areaaverage(self.Snow,self.TopoId) + + # calculate catchmentsize + self.upsize = catchmenttotal(self.xl * self.yl, self.TopoLdd) + self.csize = areamaximum(self.upsize, self.TopoId) + # Save some summary maps + self.logger.info("Saving summary maps...") + if self.modelSnow: + report(self.Cfmax, self.Dir + "/" + self.runId + "/outsum/Cfmax.map") + report(self.TTI, self.Dir + "/" + self.runId + "/outsum/TTI.map") + report(self.TT, self.Dir + "/" + self.runId + "/outsum/TT.map") + report(self.WHC, self.Dir + "/" + self.runId + "/outsum/WHC.map") + + report(self.Cmax, self.Dir + "/" + self.runId + "/outsum/Cmax.map") + report(self.csize, self.Dir + "/" + self.runId + "/outsum/CatchmentSize.map") + report(self.upsize, self.Dir + "/" + self.runId + "/outsum/UpstreamSize.map") + report(self.EoverR, self.Dir + "/" + self.runId + "/outsum/EoverR.map") + report(self.RootingDepth, self.Dir + "/" + self.runId + "/outsum/RootingDepth.map") + report(self.CanopyGapFraction, self.Dir + "/" + self.runId + "/outsum/CanopyGapFraction.map") + report(self.InfiltCapSoil, self.Dir + "/" + self.runId + "/outsum/InfiltCapSoil.map") + report(self.InfiltCapPath, self.Dir + "/" + self.runId + "/outsum/InfiltCapPath.map") + report(self.PathFrac, self.Dir + "/" + self.runId + "/outsum/PathFrac.map") + report(self.thetaR, self.Dir + "/" + self.runId + "/outsum/thetaR.map") + report(self.thetaS, self.Dir + "/" + self.runId + "/outsum/thetaS.map") + report(self.FirstZoneMinCapacity, self.Dir + "/" + self.runId + "/outsum/FirstZoneMinCapacity.map") + report(self.FirstZoneKsatVer, self.Dir + "/" + self.runId + "/outsum/FirstZoneKsatVer.map") + report(self.M, self.Dir + "/" + self.runId + "/outsum/M.map") + report(self.FirstZoneCapacity, self.Dir + "/" + self.runId + "/outsum/FirstZoneCapacity.map") + report(Terrain_angle, self.Dir + "/" + self.runId + "/outsum/angle.map") + report(self.Slope, self.Dir + "/" + self.runId + "/outsum/slope.map") + report(WI, self.Dir + "/" + self.runId + "/outsum/WI.map") + report(self.CC, self.Dir + "/" + self.runId + "/outsum/CC.map") + report(self.N, self.Dir + "/" + self.runId + "/outsum/N.map") + report(self.RiverFrac, self.Dir + "/" + self.runId + "/outsum/RiverFrac.map") + + report(self.xl, self.Dir + "/" + self.runId + "/outsum/xl.map") + report(self.yl, self.Dir + "/" + self.runId + "/outsum/yl.map") + report(self.reallength, self.Dir + "/" + self.runId + "/outsum/rl.map") + report(self.DCL, self.Dir + "/" + self.runId + "/outsum/DCL.map") + report(self.Bw, self.Dir + "/" + self.runId + "/outsum/Bw.map") + report(ifthen(self.River, self.Bw), self.Dir + "/" + self.runId + "/outsum/RiverWidth.map") + if self.updating: + report(self.DistToUpdPt, self.Dir + "/" + self.runId + "/outsum/DistToUpdPt.map") + + self.SaveDir = self.Dir + "/" + self.runId + "/" + self.logger.info("End of initial section") + + + def resume(self): + + if self.reinit == 1: + self.logger.info("Setting initial conditions to default") + self.FirstZoneDepth = self.FirstZoneCapacity * 0.85 + self.UStoreDepth = self.FirstZoneCapacity * 0.0 + self.WaterLevel = self.ZeroMap + self.SurfaceRunoff = self.ZeroMap + self.Snow = self.ZeroMap + self.SnowWater = self.ZeroMap + self.TSoil = self.ZeroMap + 10.0 + self.CanopyStorage = self.ZeroMap + + else: + self.logger.info("Setting initial conditions from state files") + self.wf_resume(self.Dir + "/instate/") + + P = self.Bw + (2.0 * self.WaterLevel) + self.Alpha = self.AlpTerm * pow(P, self.AlpPow) + self.OldSurfaceRunoff = self.SurfaceRunoff + + self.SurfaceRunoffMM = self.SurfaceRunoff * self.QMMConv # Determine initial kinematic wave volume - self.KinWaveVolume = self.WaterLevel * self.Bw * self.DCL - self.OldKinWaveVolume = self.KinWaveVolume - - self.SurfaceRunoffMM=self.SurfaceRunoff * self.QMMConv - self.InitialStorage = self.FirstZoneDepth + self.UStoreDepth + self.CanopyStorage - self.CellStorage = self.FirstZoneDepth + self.UStoreDepth - - # Determine actual water depth - self.zi = max(0.0,self.FirstZoneThickness - self.FirstZoneDepth/(self.thetaS -self.thetaR)) + self.KinWaveVolume = self.WaterLevel * self.Bw * self.DCL + self.OldKinWaveVolume = self.KinWaveVolume + + self.SurfaceRunoffMM = self.SurfaceRunoff * self.QMMConv + self.InitialStorage = self.FirstZoneDepth + self.UStoreDepth + self.CanopyStorage + self.CellStorage = self.FirstZoneDepth + self.UStoreDepth + + # Determine actual water depth + self.zi = max(0.0, self.FirstZoneThickness - self.FirstZoneDepth / (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.FirstZoneThickness / self.RunoffGeneratingGWPerc + 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.FirstZoneThickness / self.RunoffGeneratingGWPerc - - - - - def dynamic(self): - """ + + def dynamic(self): + """ Stuf that is done for each timestep @@ -827,438 +854,450 @@ :var self.ToCubic: Mutiplier to convert mm to m^3/s for fluxes """ - self.logger.debug("Step: "+str(int(self.thestep + self._d_firstTimeStep))+"/"+str(int(self._d_nrTimeSteps))) - self.thestep = self.thestep + 1 - - if self.scalarInput: + self.logger.debug( + "Step: " + str(int(self.thestep + self._d_firstTimeStep)) + "/" + str(int(self._d_nrTimeSteps))) + self.thestep = self.thestep + 1 + + if self.scalarInput: # gaugesmap not yet finished. Should be a map with cells that # hold the gauges with an unique id - self.Precipitation = timeinputscalar(self.caseName + self.precipTss,self.gaugesMap) + self.Precipitation = timeinputscalar(self.caseName + self.precipTss, self.gaugesMap) if (os.path.exists(self.caseName + self.inflowTss)): - self.Inflow = cover(timeinputscalar(self.caseName + self.inflowTss,nominal(self.InflowLoc)),0) + self.Inflow = cover(timeinputscalar(self.caseName + self.inflowTss, nominal(self.InflowLoc)), 0) else: self.Inflow = self.ZeroMap - self.Precipitation = pcrut.interpolategauges(self.Precipitation,self.interpolMethod) - self.PotenEvap=timeinputscalar(self.caseName + self.evapTss,self.gaugesMap) - self.PotenEvap = pcrut.interpolategauges(self.PotenEvap,self.interpolMethod) + self.Precipitation = pcrut.interpolategauges(self.Precipitation, self.interpolMethod) + self.PotenEvap = timeinputscalar(self.caseName + self.evapTss, self.gaugesMap) + self.PotenEvap = pcrut.interpolategauges(self.PotenEvap, self.interpolMethod) if self.modelSnow: - self.Temperature=timeinputscalar(self.caseName + self.tempTss,self.gaugesMap) - self.Temperature = pcrut.interpolategauges(self.Temperature,self.interpolMethod) + self.Temperature = timeinputscalar(self.caseName + self.tempTss, self.gaugesMap) + self.Temperature = pcrut.interpolategauges(self.Temperature, self.interpolMethod) self.Temperature = self.Temperature + self.TempCor - else: - self.Precipitation=cover(self.wf_readmap(self.P_mapstack,0.0),0) - self.PotenEvap=cover(self.wf_readmap(self.PET_mapstack,0.0),0) + else: + self.Precipitation = cover(self.wf_readmap(self.P_mapstack, 0.0), 0) + self.PotenEvap = cover(self.wf_readmap(self.PET_mapstack, 0.0), 0) #self.Inflow=cover(self.wf_readmap(self.Inflow),0) if (os.path.exists(self.caseName + self.inflowTss)): - self.Inflow = cover(timeinputscalar(self.caseName + self.inflowTss,nominal(self.InflowLoc)),0) + self.Inflow = cover(timeinputscalar(self.caseName + self.inflowTss, nominal(self.InflowLoc)), 0) else: - self.Inflow=pcrut.readmapSave(self.Inflow_mapstack,0.0) + self.Inflow = pcrut.readmapSave(self.Inflow_mapstack, 0.0) #self.Inflow=spatial(scalar(0.0)) if self.modelSnow: - self.Temperature=cover(self.wf_readmap(self.TEMP_mapstack,10.0),10.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 - for k, v in multdynapars.iteritems(): - estr = k + "=" + k + "*" + str(v) - self.logger.debug("Dynamic Parameter multiplication: " + estr) - exec estr - - self.PotEvap = self.PotenEvap # - #TODO: Snow modelling if enabled _ need to be moved as it breaks the scalar input - """ + + + # Multiply input parameters with a factor (for calibration etc) -p option in command line + for k, v in multdynapars.iteritems(): + estr = k + "=" + k + "*" + str(v) + self.logger.debug("Dynamic Parameter multiplication: " + estr) + exec estr + + self.PotEvap = self.PotenEvap # + #TODO: Snow modelling if enabled _ need to be moved as it breaks the scalar input + """ .. todo:: Snow modelling if enabled _ needs to be moved as it breaks the scalar input """ - if self.modelSnow: + 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.Precipitation = SnowPackHBV(self.Snow,self.SnowWater,self.Precipitation,self.Temperature,self.TTI,self.TT,self.Cfmax,self.WHC) - - ########################################################################## - # Interception according to a modified Gash model - # TODOs: add sub-daily interception! - ########################################################################## + self.Snow, self.SnowWater, self.SnowMelt, self.Precipitation = SnowPackHBV(self.Snow, self.SnowWater, + self.Precipitation, + self.Temperature, self.TTI, + self.TT, self.Cfmax, self.WHC) - if self.timestepsecs >= (23 * 3600): - ThroughFall, Interception, StemFlow, self.CanopyStorage = rainfall_interception_gash(self.Cmax,self.EoverR,self.CanopyGapFraction, self.Precipitation,self.CanopyStorage) - PotTrans = cover(max(0.0,self.PotEvap - Interception),0.0) # now in mm - else: - NetInterception, ThroughFall, StemFlow, LeftOver, Interception, self.CanopyStorage = rainfall_interception_modrut(self.Precipitation,self.PotEvap,self.CanopyStorage,self.CanopyGapFraction,self.Cmax) - PotTrans = cover(max(0.0,LeftOver),0.0) # now in mm + ########################################################################## + # Interception according to a modified Gash model + # TODOs: add sub-daily interception! + ########################################################################## - ########################################################################## - # Start with the soil calculations ###################################### - ########################################################################## + if self.timestepsecs >= (23 * 3600): + ThroughFall, Interception, StemFlow, self.CanopyStorage = rainfall_interception_gash(self.Cmax, self.EoverR, + self.CanopyGapFraction, + self.Precipitation, + self.CanopyStorage) + PotTrans = cover(max(0.0, self.PotEvap - Interception), 0.0) # now in mm + else: + NetInterception, ThroughFall, StemFlow, LeftOver, Interception, self.CanopyStorage = rainfall_interception_modrut( + self.Precipitation, self.PotEvap, self.CanopyStorage, self.CanopyGapFraction, self.Cmax) + PotTrans = cover(max(0.0, LeftOver), 0.0) # now in mm + ########################################################################## + # Start with the soil calculations ###################################### + ########################################################################## - self.ExfiltWater=self.ZeroMap - FreeWaterDepth=self.ZeroMap - - ########################################################################## - # Determine infiltration into Unsaturated store...######################## - ########################################################################## - # Add precipitation surplus FreeWater storage... - FreeWaterDepth= ThroughFall + StemFlow - UStoreCapacity = self.FirstZoneCapacity - self.FirstZoneDepth - self.UStoreDepth - - # Runoff onto water boddies and river network - self.RunoffOpenWater = self.RiverFrac * self.WaterFrac * FreeWaterDepth - #self.RunoffOpenWater = self.ZeroMap - FreeWaterDepth = FreeWaterDepth - self.RunoffOpenWater - - if self.RunoffGenSigmaFunction: - self.AbsoluteGW=self.DemMax-(self.zi*self.GWScale) - self.SubCellFrac = sCurve(self.AbsoluteGW,c=self.CC,a=self.Altitude+1.0) - self.SubCellRunoff = self.SubCellFrac * FreeWaterDepth - self.SubCellGWRunoff = min(self.SubCellFrac * self.FirstZoneDepth, self.SubCellFrac * self.Slope * self.FirstZoneKsatVer * exp(-self.f * self.zi)) - self.FirstZoneDepth=self.FirstZoneDepth-self.SubCellGWRunoff - FreeWaterDepth = FreeWaterDepth - self.SubCellRunoff - else: - self.AbsoluteGW=self.DemMax-(self.zi*self.GWScale) - self.SubCellFrac = spatial(scalar(0.0)) - self.SubCellGWRunoff = spatial(scalar(0.0)) - self.SubCellRunoff = spatial(scalar(0.0)) - - #----->> - # First determine if the soil infiltration capacity can deal with the - # amount of water - # split between infiltration in undisturbed soil and compacted areas (paths) - SoilInf = FreeWaterDepth * (1- self.PathFrac) - PathInf = FreeWaterDepth * self.PathFrac - if self.modelSnow: - # soilInfRedu = ifthenelse(self.TSoil < 0.0 , self.cf_soil, 1.0) - bb = 1.0/(1.0 - self.cf_soil) - soilInfRedu = sCurve(self.TSoil,a=self.ZeroMap,b=bb,c=8.0) + self.cf_soil - else: - soilInfRedu = 1.0 - MaxInfiltSoil= min(self.InfiltCapSoil*soilInfRedu,SoilInf) - - self.SoilInfiltExceeded=self.SoilInfiltExceeded + scalar(self.InfiltCapSoil*soilInfRedu < SoilInf) - InfiltSoil = min(MaxInfiltSoil, UStoreCapacity) - self.UStoreDepth = self.UStoreDepth + InfiltSoil - UStoreCapacity = UStoreCapacity - InfiltSoil - FreeWaterDepth = FreeWaterDepth - InfiltSoil - # <------- - MaxInfiltPath= min(self.InfiltCapPath*soilInfRedu,PathInf) - #self.PathInfiltExceeded=self.PathInfiltExceeded + ifthenelse(self.InfiltCapPath < FreeWaterDepth, scalar(1), scalar(0)) - self.PathInfiltExceeded=self.PathInfiltExceeded + scalar(self.InfiltCapPath*soilInfRedu < PathInf) - InfiltPath = min(MaxInfiltPath, UStoreCapacity) - self.UStoreDepth = self.UStoreDepth + InfiltPath - UStoreCapacity = UStoreCapacity - InfiltPath - FreeWaterDepth = FreeWaterDepth - InfiltPath - - self.ActInfilt = InfiltPath + InfiltSoil + self.ExfiltWater = self.ZeroMap + FreeWaterDepth = self.ZeroMap - self.InfiltExcess = ifthenelse (UStoreCapacity > 0.0, FreeWaterDepth, 0.0) - self.CumInfiltExcess=self.CumInfiltExcess+self.InfiltExcess - - self.ActEvap, self.FirstZoneDepth, self.UStoreDepth, self.ActEvapUStore = actEvap_SBM(self.RootingDepth,self.zi,self.UStoreDepth,self.FirstZoneDepth, PotTrans,self.rootdistpar) - #self.ActEvap = self.ZeroMap - #self.ActEvapUStore = self.ZeroMap - ########################################################################## - # Transfer of water from unsaturated to saturated store...################ - ########################################################################## - self.zi = max(0.0,self.FirstZoneThickness - self.FirstZoneDepth/(self.thetaS -self.thetaR)) # Determine actual water depth - Ksat = self.FirstZoneKsatVer * exp(-self.f * self.zi) - self.DeepKsat = self.FirstZoneKsatVer * exp(-self.f * self.FirstZoneThickness) - - # Determine saturation deficit. NB, as noted by Vertessy and Elsenbeer 1997 - # this deficit does NOT take into account the water in the unsaturated zone - SaturationDeficit = self.FirstZoneCapacity - self.FirstZoneDepth - + ########################################################################## + # Determine infiltration into Unsaturated store...######################## + ########################################################################## + # Add precipitation surplus FreeWater storage... + FreeWaterDepth = ThroughFall + StemFlow + UStoreCapacity = self.FirstZoneCapacity - self.FirstZoneDepth - self.UStoreDepth - # now the actual tranfer to the saturated store.. - self.Transfer = min(self.UStoreDepth,ifthenelse (SaturationDeficit <= 0.00001, 0.0, Ksat * self.UStoreDepth/(SaturationDeficit+1))) - # Determine Ksat at base - #DeepTransfer = min(self.UStoreDepth,ifthenelse (SaturationDeficit <= 0.00001, 0.0, DeepKsat * self.UStoreDepth/(SaturationDeficit+1))) - ActLeakage = 0.0 - # Now add leakage. to deeper groundwater - #ActLeakage = cover(max(0,min(self.MaxLeakage,ActLeakage)),0) - - # Now look if there is Seeapage - - #ActLeakage = ifthenelse(self.Seepage > 0.0, -1.0 * Seepage, ActLeakage) - self.FirstZoneDepth = self.FirstZoneDepth + self.Transfer - ActLeakage - self.UStoreDepth = self.UStoreDepth - self.Transfer + # Runoff onto water boddies and river network + self.RunoffOpenWater = self.RiverFrac * self.WaterFrac * FreeWaterDepth + #self.RunoffOpenWater = self.ZeroMap + FreeWaterDepth = FreeWaterDepth - self.RunoffOpenWater - # Determine % saturated - #Sat = ifthenelse(self.FirstZoneDepth >= (self.FirstZoneCapacity*0.999), scalar(1.0), scalar(0.0)) - self.Sat = max(self.SubCellFrac,scalar(self.FirstZoneDepth >= (self.FirstZoneCapacity*0.999))) - #PercSat = areaaverage(scalar(Sat),self.TopoId) * 100 + if self.RunoffGenSigmaFunction: + self.AbsoluteGW = self.DemMax - (self.zi * self.GWScale) + self.SubCellFrac = sCurve(self.AbsoluteGW, c=self.CC, a=self.Altitude + 1.0) + self.SubCellRunoff = self.SubCellFrac * FreeWaterDepth + self.SubCellGWRunoff = min(self.SubCellFrac * self.FirstZoneDepth, + self.SubCellFrac * self.Slope * self.FirstZoneKsatVer * exp(-self.f * self.zi)) + self.FirstZoneDepth = self.FirstZoneDepth - self.SubCellGWRunoff + FreeWaterDepth = FreeWaterDepth - self.SubCellRunoff + else: + self.AbsoluteGW = self.DemMax - (self.zi * self.GWScale) + self.SubCellFrac = spatial(scalar(0.0)) + self.SubCellGWRunoff = spatial(scalar(0.0)) + self.SubCellRunoff = spatial(scalar(0.0)) + #----->> + # First determine if the soil infiltration capacity can deal with the + # amount of water + # split between infiltration in undisturbed soil and compacted areas (paths) - ########################################################################## - # Horizontal (downstream) transport of water ############################# - ########################################################################## - - if self.waterdem: - waterDem = self.Altitude - (self.zi * 0.001) - waterLdd = lddcreate(waterDem,1E35,1E35,1E35,1E35) - #waterLdd = lddcreate(waterDem,1,1,1,1) - waterSlope=max(0.00001,slope(waterDem)*celllength()/self.reallength) - - self.zi = max(0.0,self.FirstZoneThickness - self.FirstZoneDepth/(self.thetaS -self.thetaR)) # Determine actual water depth + SoilInf = FreeWaterDepth * (1 - self.PathFrac) + PathInf = FreeWaterDepth * self.PathFrac + if self.modelSnow: + # soilInfRedu = ifthenelse(self.TSoil < 0.0 , self.cf_soil, 1.0) + bb = 1.0 / (1.0 - self.cf_soil) + soilInfRedu = sCurve(self.TSoil, a=self.ZeroMap, b=bb, c=8.0) + self.cf_soil + else: + soilInfRedu = 1.0 + MaxInfiltSoil = min(self.InfiltCapSoil * soilInfRedu, SoilInf) - if self.waterdem: - MaxHor = max(0.0,min(self.FirstZoneKsatVer * waterSlope * exp(-SaturationDeficit/self.M),self.FirstZoneDepth)) - self.FirstZoneFlux = accucapacityflux (waterLdd, self.FirstZoneDepth, MaxHor) - self.FirstZoneDepth = accucapacitystate (waterLdd, self.FirstZoneDepth, MaxHor) - else: - # - #MaxHor = max(0,min(self.FirstZoneKsatVer * self.Slope * exp(-SaturationDeficit/self.M),self.FirstZoneDepth*(self.thetaS-self.thetaR))) * timestepsecs/basetimestep - MaxHor = max(0.0,min(self.FirstZoneKsatVer * self.Slope * exp(-SaturationDeficit/self.M),self.FirstZoneDepth)) - self.FirstZoneFlux = accucapacityflux (self.TopoLdd, self.FirstZoneDepth, MaxHor) - self.FirstZoneDepth = accucapacitystate (self.TopoLdd, self.FirstZoneDepth, MaxHor) - + self.SoilInfiltExceeded = self.SoilInfiltExceeded + scalar(self.InfiltCapSoil * soilInfRedu < SoilInf) + InfiltSoil = min(MaxInfiltSoil, UStoreCapacity) + self.UStoreDepth = self.UStoreDepth + InfiltSoil + UStoreCapacity = UStoreCapacity - InfiltSoil + FreeWaterDepth = FreeWaterDepth - InfiltSoil + # <------- + MaxInfiltPath = min(self.InfiltCapPath * soilInfRedu, PathInf) + #self.PathInfiltExceeded=self.PathInfiltExceeded + ifthenelse(self.InfiltCapPath < FreeWaterDepth, scalar(1), scalar(0)) + self.PathInfiltExceeded = self.PathInfiltExceeded + scalar(self.InfiltCapPath * soilInfRedu < PathInf) + InfiltPath = min(MaxInfiltPath, UStoreCapacity) + self.UStoreDepth = self.UStoreDepth + InfiltPath + UStoreCapacity = UStoreCapacity - InfiltPath + FreeWaterDepth = FreeWaterDepth - InfiltPath + self.ActInfilt = InfiltPath + InfiltSoil + self.InfiltExcess = ifthenelse(UStoreCapacity > 0.0, FreeWaterDepth, 0.0) + self.CumInfiltExcess = self.CumInfiltExcess + self.InfiltExcess - ########################################################################## - # Determine returnflow from first zone ########################## - ########################################################################## - self.ExfiltWaterFrac = sCurve(self.FirstZoneDepth,a=self.FirstZoneCapacity,c=5.0) - self.ExfiltWater=self.ExfiltWaterFrac * (self.FirstZoneDepth - self.FirstZoneCapacity) - #self.ExfiltWater=ifthenelse (self.FirstZoneDepth - self.FirstZoneCapacity > 0 , self.FirstZoneDepth - self.FirstZoneCapacity , 0.0) - self.FirstZoneDepth=self.FirstZoneDepth - self.ExfiltWater - - - # Re-determine UStoreCapacity - UStoreCapacity = self.FirstZoneCapacity - self.FirstZoneDepth - self.UStoreDepth - #Determine capilary rise - self.zi = max(0.0,self.FirstZoneThickness - self.FirstZoneDepth/(self.thetaS -self.thetaR)) # Determine actual water depth - Ksat = self.FirstZoneKsatVer * exp(-self.f * self.zi) + self.ActEvap, self.FirstZoneDepth, self.UStoreDepth, self.ActEvapUStore = actEvap_SBM(self.RootingDepth, + self.zi, self.UStoreDepth, + self.FirstZoneDepth, + PotTrans, + self.rootdistpar) + #self.ActEvap = self.ZeroMap + #self.ActEvapUStore = self.ZeroMap + ########################################################################## + # Transfer of water from unsaturated to saturated store...################ + ########################################################################## + self.zi = max(0.0, self.FirstZoneThickness - self.FirstZoneDepth / ( + self.thetaS - self.thetaR)) # Determine actual water depth + Ksat = self.FirstZoneKsatVer * exp(-self.f * self.zi) + self.DeepKsat = self.FirstZoneKsatVer * exp(-self.f * self.FirstZoneThickness) - MaxCapFlux = max(0.0,min(Ksat,self.ActEvapUStore,UStoreCapacity,self.FirstZoneDepth)) - # 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.RootingDepth, self.CapScale/(self.CapScale + self.zi -self.RootingDepth), 0.0) - self.CapFlux = MaxCapFlux * CapFluxScale - - - self.UStoreDepth = self.UStoreDepth + self.CapFlux - self.FirstZoneDepth = self.FirstZoneDepth - self.CapFlux - - # org SurfaceWater = self.SurfaceRunoff * self.DCL * self.QMMConv # SurfaceWater (mm) from SurfaceRunoff (m3/s) - SurfaceWater = self.SurfaceRunoff * self.QMMConv # SurfaceWater (mm) from SurfaceRunoff (m3/s) - self.CumSurfaceWater = self.CumSurfaceWater + SurfaceWater + # Determine saturation deficit. NB, as noted by Vertessy and Elsenbeer 1997 + # this deficit does NOT take into account the water in the unsaturated zone + SaturationDeficit = self.FirstZoneCapacity - self.FirstZoneDepth - # Estimate water that may re-infiltrate - if self.reInfilt: - Reinfilt = max(0,min(SurfaceWater,min(self.InfiltCapSoil,UStoreCapacity))) - self.CumReinfilt=self.CumReinfilt + Reinfilt + + # now the actual tranfer to the saturated store.. + self.Transfer = min(self.UStoreDepth, ifthenelse(SaturationDeficit <= 0.00001, 0.0, + Ksat * self.UStoreDepth / (SaturationDeficit + 1))) + # Determine Ksat at base + #DeepTransfer = min(self.UStoreDepth,ifthenelse (SaturationDeficit <= 0.00001, 0.0, DeepKsat * self.UStoreDepth/(SaturationDeficit+1))) + ActLeakage = 0.0 + # Now add leakage. to deeper groundwater + #ActLeakage = cover(max(0,min(self.MaxLeakage,ActLeakage)),0) + + # Now look if there is Seeapage + + #ActLeakage = ifthenelse(self.Seepage > 0.0, -1.0 * Seepage, ActLeakage) + self.FirstZoneDepth = self.FirstZoneDepth + self.Transfer - ActLeakage + self.UStoreDepth = self.UStoreDepth - self.Transfer + + # Determine % saturated + #Sat = ifthenelse(self.FirstZoneDepth >= (self.FirstZoneCapacity*0.999), scalar(1.0), scalar(0.0)) + self.Sat = max(self.SubCellFrac, scalar(self.FirstZoneDepth >= (self.FirstZoneCapacity * 0.999))) + #PercSat = areaaverage(scalar(Sat),self.TopoId) * 100 + + + ########################################################################## + # Horizontal (downstream) transport of water ############################# + ########################################################################## + + if self.waterdem: + waterDem = self.Altitude - (self.zi * 0.001) + waterLdd = lddcreate(waterDem, 1E35, 1E35, 1E35, 1E35) + #waterLdd = lddcreate(waterDem,1,1,1,1) + waterSlope = max(0.00001, slope(waterDem) * celllength() / self.reallength) + + self.zi = max(0.0, self.FirstZoneThickness - self.FirstZoneDepth / ( + self.thetaS - self.thetaR)) # Determine actual water depth + + if self.waterdem: + MaxHor = max(0.0, min(self.FirstZoneKsatVer * waterSlope * exp(-SaturationDeficit / self.M), + self.FirstZoneDepth)) + self.FirstZoneFlux = accucapacityflux(waterLdd, self.FirstZoneDepth, MaxHor) + self.FirstZoneDepth = accucapacitystate(waterLdd, self.FirstZoneDepth, MaxHor) + else: + # + #MaxHor = max(0,min(self.FirstZoneKsatVer * self.Slope * exp(-SaturationDeficit/self.M),self.FirstZoneDepth*(self.thetaS-self.thetaR))) * timestepsecs/basetimestep + MaxHor = max(0.0, min(self.FirstZoneKsatVer * self.Slope * exp(-SaturationDeficit / self.M), + self.FirstZoneDepth)) + self.FirstZoneFlux = accucapacityflux(self.TopoLdd, self.FirstZoneDepth, MaxHor) + self.FirstZoneDepth = accucapacitystate(self.TopoLdd, self.FirstZoneDepth, MaxHor) + + + + + ########################################################################## + # Determine returnflow from first zone ########################## + ########################################################################## + self.ExfiltWaterFrac = sCurve(self.FirstZoneDepth, a=self.FirstZoneCapacity, c=5.0) + self.ExfiltWater = self.ExfiltWaterFrac * (self.FirstZoneDepth - self.FirstZoneCapacity) + #self.ExfiltWater=ifthenelse (self.FirstZoneDepth - self.FirstZoneCapacity > 0 , self.FirstZoneDepth - self.FirstZoneCapacity , 0.0) + self.FirstZoneDepth = self.FirstZoneDepth - self.ExfiltWater + + + # Re-determine UStoreCapacity + UStoreCapacity = self.FirstZoneCapacity - self.FirstZoneDepth - self.UStoreDepth + #Determine capilary rise + self.zi = max(0.0, self.FirstZoneThickness - self.FirstZoneDepth / ( + self.thetaS - self.thetaR)) # Determine actual water depth + Ksat = self.FirstZoneKsatVer * exp(-self.f * self.zi) + + MaxCapFlux = max(0.0, min(Ksat, self.ActEvapUStore, UStoreCapacity, self.FirstZoneDepth)) + # 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.RootingDepth, + self.CapScale / (self.CapScale + self.zi - self.RootingDepth), 0.0) + self.CapFlux = MaxCapFlux * CapFluxScale + + self.UStoreDepth = self.UStoreDepth + self.CapFlux + self.FirstZoneDepth = self.FirstZoneDepth - self.CapFlux + + # org SurfaceWater = self.SurfaceRunoff * self.DCL * self.QMMConv # SurfaceWater (mm) from SurfaceRunoff (m3/s) + SurfaceWater = self.SurfaceRunoff * self.QMMConv # SurfaceWater (mm) from SurfaceRunoff (m3/s) + self.CumSurfaceWater = self.CumSurfaceWater + SurfaceWater + + # Estimate water that may re-infiltrate + if self.reInfilt: + Reinfilt = max(0, min(SurfaceWater, min(self.InfiltCapSoil, UStoreCapacity))) + self.CumReinfilt = self.CumReinfilt + Reinfilt self.UStoreDepth = self.UStoreDepth + Reinfilt - else: + else: Reinfilt = self.ZeroMap - - - self.InwaterMM=max(0.0,self.ExfiltWater + FreeWaterDepth + self.SubCellRunoff + self.SubCellGWRunoff + self.RunoffOpenWater - Reinfilt) - self.Inwater=self.InwaterMM * self.ToCubic # m3/s - - 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.FreeWaterDepthCubic=FreeWaterDepth * self.ToCubic - self.ReinfiltCubic=-1.0 * Reinfilt * self.ToCubic - self.Inwater=self.Inwater + self.Inflow # Add abstractions/inflows in m^3/sec - - ########################################################################## - # 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.SurfaceRunoffMM=self.SurfaceRunoff*self.QMMConv # SurfaceRunoffMM (mm) from SurfaceRunoff (m3/s) - self.updateRunOff() - self.InflowKinWaveCell=upstream(self.TopoLdd,self.SurfaceRunoff) - self.MassBalKinWave = (self.KinWaveVolume - self.OldKinWaveVolume)/self.timestepsecs + self.InflowKinWaveCell + self.Inwater - self.SurfaceRunoff + self.InwaterMM = max(0.0, + self.ExfiltWater + FreeWaterDepth + self.SubCellRunoff + self.SubCellGWRunoff + self.RunoffOpenWater - Reinfilt) + self.Inwater = self.InwaterMM * self.ToCubic # m3/s - - - Runoff=self.SurfaceRunoff + 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.FreeWaterDepthCubic = FreeWaterDepth * self.ToCubic + self.ReinfiltCubic = -1.0 * Reinfilt * self.ToCubic + self.Inwater = self.Inwater + self.Inflow # Add abstractions/inflows in m^3/sec - # 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! + ########################################################################## + # Runoff calculation via Kinematic wave ################################## + ########################################################################## + # per distance along stream + q = self.Inwater / self.DCL + # discharge (m3/s) - if self.updating: - QM = timeinputscalar(self.updateFile, self.UpdateMap) * self.QMMConv - - # 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) - - 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.UStoreDepth * UpRatioSoil) - self.UStoreDepth - self.UStoreDepth = self.UStoreDepth + 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.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.updateRunOff() - - Runoff=self.SurfaceRunoff - - ########################################################################## - # water balance ########################################### - ########################################################################## + self.InflowKinWaveCell = upstream(self.TopoLdd, self.SurfaceRunoff) + self.MassBalKinWave = ( + self.KinWaveVolume - self.OldKinWaveVolume) / self.timestepsecs + self.InflowKinWaveCell + self.Inwater - self.SurfaceRunoff - # Single cell based water budget - CellStorage = self.UStoreDepth+self.FirstZoneDepth + self.CanopyStorage - DeltaStorage = CellStorage - self.InitialStorage - OutFlow = self.FirstZoneFlux - CellInFlow = upstream(self.TopoLdd,scalar(self.FirstZoneFlux)); - #CellWatBal = ActInfilt - self.ActEvap - self.ExfiltWater - ActLeakage + Reinfilt + IF - OutFlow + (OldCellStorage - CellStorage) - #SumCellWatBal = SumCellWatBal + CellWatBal; + Runoff = self.SurfaceRunoff - self.CumOutFlow = self.CumOutFlow + OutFlow - self.CumActInfilt = self.CumActInfilt + self.ActInfilt - self.CumCellInFlow = self.CumCellInFlow + CellInFlow - self.CumPrec=self.CumPrec+self.Precipitation - self.CumEvap=self.CumEvap+self.ActEvap - self.CumPotenTrans=self.CumPotenTrans+PotTrans - self.CumPotenEvap=self.CumPotenEvap+self.PotenEvap - if self.timestepsecs >= (23 * 3600): - self.CumInt=self.CumInt+Interception - else: - self.CumInt=self.CumInt+NetInterception - - self.CumLeakage=self.CumLeakage+ActLeakage - self.CumInwaterMM=self.CumInwaterMM+self.InwaterMM - self.CumExfiltWater=self.CumExfiltWater+self.ExfiltWater - # Water budget - #self.watbal = self.CumPrec- self.CumEvap - self.CumInt - self.CumInwaterMM - DeltaStorage - self.CumOutFlow + self.CumIF - #self.watbal = self.CumActInfilt - self.CumEvap - self.CumExfiltWater - DeltaStorage - self.CumOutFlow + self.CumIF - self.watbal = self.CumPrec + self.CumCellInFlow - self.CumOutFlow- self.CumEvap - self.CumLeakage - self.CumInwaterMM - self.CumInt - DeltaStorage + self.CumReinfilt + # 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 + # 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) -def main(argv=None): - + 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.UStoreDepth * UpRatioSoil) - self.UStoreDepth + self.UStoreDepth = self.UStoreDepth + 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 + + ########################################################################## + # water balance ########################################### + ########################################################################## + + # Single cell based water budget + CellStorage = self.UStoreDepth + self.FirstZoneDepth + self.CanopyStorage + DeltaStorage = CellStorage - self.InitialStorage + OutFlow = self.FirstZoneFlux + CellInFlow = upstream(self.TopoLdd, scalar(self.FirstZoneFlux)); + #CellWatBal = ActInfilt - self.ActEvap - self.ExfiltWater - ActLeakage + Reinfilt + IF - OutFlow + (OldCellStorage - CellStorage) + #SumCellWatBal = SumCellWatBal + CellWatBal; + + self.CumOutFlow = self.CumOutFlow + OutFlow + self.CumActInfilt = self.CumActInfilt + self.ActInfilt + self.CumCellInFlow = self.CumCellInFlow + CellInFlow + self.CumPrec = self.CumPrec + self.Precipitation + self.CumEvap = self.CumEvap + self.ActEvap + self.CumPotenTrans = self.CumPotenTrans + PotTrans + self.CumPotenEvap = self.CumPotenEvap + self.PotenEvap + if self.timestepsecs >= (23 * 3600): + self.CumInt = self.CumInt + Interception + else: + self.CumInt = self.CumInt + NetInterception + + self.CumLeakage = self.CumLeakage + ActLeakage + self.CumInwaterMM = self.CumInwaterMM + self.InwaterMM + self.CumExfiltWater = self.CumExfiltWater + self.ExfiltWater + # Water budget + #self.watbal = self.CumPrec- self.CumEvap - self.CumInt - self.CumInwaterMM - DeltaStorage - self.CumOutFlow + self.CumIF + #self.watbal = self.CumActInfilt - self.CumEvap - self.CumExfiltWater - DeltaStorage - self.CumOutFlow + self.CumIF + self.watbal = self.CumPrec + self.CumCellInFlow - self.CumOutFlow - self.CumEvap - self.CumLeakage - self.CumInwaterMM - self.CumInt - DeltaStorage + self.CumReinfilt + + +def main(argv=None): """ Perform command line execution of the model. - """ + """ caseName = "default_sbm" - global multpars + global multpars runId = "run_default" - configfile="wflow_sbm.ini" + configfile = "wflow_sbm.ini" _lastTimeStep = 0 _firstTimeStep = 1 - LogFileName="wflow.log" - fewsrun=False - runinfoFile="runinfo.xml" - timestepsecs=86400 + LogFileName = "wflow.log" + fewsrun = False + runinfoFile = "runinfo.xml" + timestepsecs = 86400 wflow_cloneMap = 'wflow_subcatch.map' - _NoOverWrite=1 + _NoOverWrite = 1 global updateCols loglevel = logging.DEBUG if argv is None: argv = sys.argv[1:] if len(argv) == 0: usage() - return - + return + ## Main model starts here ######################################################################## try: opts, args = getopt.getopt(argv, 'XF:L: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 == '-F': + if o == '-F': runinfoFile = a fewsrun = True - if o == '-P': - exec ("multpars =" + a,globals(), globals()) - if o == '-p': + if o == '-P': + exec ("multpars =" + a, globals(), globals()) + if o == '-p': exec "multdynapars =" + a - exec ("multdynapars =" + a,globals(), globals()) + exec ("multdynapars =" + a, globals(), globals()) 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 == '-T': _lastTimeStep=int(a) - if o == '-S': _firstTimeStep=int(a) + if o == '-T': _lastTimeStep = int(a) + if o == '-S': _firstTimeStep = int(a) if o == '-h': usage() if o == '-f': _NoOverWrite = 0 if o == '-l': exec "loglevel = logging." + a - - - if fewsrun: - ts = getTimeStepsfromRuninfo(runinfoFile,timestepsecs) + if fewsrun: + ts = getTimeStepsfromRuninfo(runinfoFile, timestepsecs) if (ts): - _lastTimeStep = ts# * 86400/timestepsecs - _firstTimeStep = 1 + _lastTimeStep = ts # * 86400/timestepsecs + _firstTimeStep = 1 else: print "Failed to get timesteps from runinfo file: " + runinfoFile exit(2) - + 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() print LogFileName - myModel = WflowModel(wflow_cloneMap, caseName,runId,configfile) - dynModelFw = wf_DynamicFramework(myModel, _lastTimeStep,firstTimestep=_firstTimeStep) - dynModelFw.createRunId(NoOverWrite=_NoOverWrite,level=loglevel,logfname=LogFileName) - + myModel = WflowModel(wflow_cloneMap, caseName, runId, configfile) + dynModelFw = wf_DynamicFramework(myModel, _lastTimeStep, firstTimestep=_firstTimeStep) + dynModelFw.createRunId(NoOverWrite=_NoOverWrite, level=loglevel, logfname=LogFileName) + for o, a in opts: - if o == '-X': configset(myModel.config,'model','OverWriteInit','1',overwrite=True) - if o == '-I': configset(myModel.config,'model','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',"1",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 == '-X': configset(myModel.config, 'model', 'OverWriteInit', '1', overwrite=True) + if o == '-I': configset(myModel.config, 'model', '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', "1", 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 == '-E': configset(myModel.config,'model','reInfilt','1',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 == '-W': configset(myModel.config, 'model', 'waterdem', '1', overwrite=True) - dynModelFw._runInitial() dynModelFw._runResume() - dynModelFw._runDynamic(_firstTimeStep,_lastTimeStep) + dynModelFw._runDynamic(_firstTimeStep, _lastTimeStep) dynModelFw._runSuspend() - - fp = open(caseName + "/" + runId + "/runinfo/configofrun.ini",'wb') - myModel.config.write(fp ) + + fp = open(caseName + "/" + runId + "/runinfo/configofrun.ini", 'wb') + myModel.config.write(fp) dynModelFw._wf_shutdown() - - if __name__ == "__main__": main()