Index: wflow-py/wflow/wflow_sphy.py =================================================================== diff -u -ra0638f1d565ab6a51c556ca0ee26269f16993557 -r12ea40dc08628f654753679e0972e87b7bb12f7a --- wflow-py/wflow/wflow_sphy.py (.../wflow_sphy.py) (revision a0638f1d565ab6a51c556ca0ee26269f16993557) +++ wflow-py/wflow/wflow_sphy.py (.../wflow_sphy.py) (revision 12ea40dc08628f654753679e0972e87b7bb12f7a) @@ -16,10 +16,9 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -#TODO: split off routing +# TODO: split off routing - import numpy import sys import os @@ -33,17 +32,18 @@ from wflow.wflow_adapt import * from wflow.wflow_adapt import * import pcraster as pcr -#from wflow.wflow_lib import reporting -#import scipy -#import pcrut +# from wflow.wflow_lib import reporting +# import scipy +# import pcrut + wflow = "wflow_sphy" #: columns used in updating -updateCols = [] #: columns used in updating +updateCols = [] #: columns used in updating """ Column used in updating """ @@ -54,45 +54,44 @@ - *args: command line arguments given """ sys.stdout = sys.stderr - for msg in args: print msg + for msg in args: + print msg print __doc__ sys.exit(0) + class WflowModel(DynamicModel): - """ + """ The user defined model class. """ + def __init__(self, cloneMap, Dir, RunDir, configfile): + DynamicModel.__init__(self) + self.caseName = os.path.abspath(Dir) + self.clonemappath = os.path.join(os.path.abspath(Dir), "staticmaps", cloneMap) + setclone(self.clonemappath) + self.runId = RunDir + self.Dir = os.path.abspath(Dir) + self.configfile = configfile + self.SaveDir = os.path.join(self.Dir, self.runId) - def __init__(self, cloneMap,Dir,RunDir,configfile): - DynamicModel.__init__(self) - self.caseName = os.path.abspath(Dir) - self.clonemappath = os.path.join(os.path.abspath(Dir),"staticmaps",cloneMap) - setclone(self.clonemappath) - self.runId = RunDir - self.Dir = os.path.abspath(Dir) - self.configfile = configfile - self.SaveDir = os.path.join(self.Dir,self.runId) - - - def updateRunOff(self): # - this may not be required - """ + def updateRunOff(self): # - this may not be required + """ Updates the kinematic wave reservoir """ - self.WaterLevel=(self.Alpha*pow(self.SurfaceRunoff,self.Beta))/self.Bw - # wetted perimeter (m) - P=self.Bw+(2*self.WaterLevel) - # Alpha - self.Alpha=self.AlpTerm*pow(P,self.AlpPow) - self.OldKinWaveVolume = self.KinWaveVolume - self.KinWaveVolume = self.WaterLevel * self.Bw * self.DCL + self.WaterLevel = (self.Alpha * pow(self.SurfaceRunoff, self.Beta)) / self.Bw + # wetted perimeter (m) + P = self.Bw + (2 * self.WaterLevel) + # Alpha + self.Alpha = self.AlpTerm * pow(P, self.AlpPow) + self.OldKinWaveVolume = self.KinWaveVolume + self.KinWaveVolume = self.WaterLevel * self.Bw * self.DCL - - def stateVariables(self): - """ + def stateVariables(self): + """ returns a list of state variables that are essential to the model. This list is essential for the resume and suspend functions to work. @@ -108,1363 +107,2129 @@ :var self.InterceptionStorage: Amount of water on the Canopy [mm] """ - states = ['RootWater','SubWater', 'CapRise', 'RootDrain','SubDrain','GwRecharge','Gw','H_gw','SnowWatStore'] ## -> complete list with required states + states = [ + "RootWater", + "SubWater", + "CapRise", + "RootDrain", + "SubDrain", + "GwRecharge", + "Gw", + "H_gw", + "SnowWatStore", + ] ## -> complete list with required states - # if hasattr(self,'ReserVoirSimpleLocs'): ## -> add states that may be required if certain modules are used - # states.append('ReservoirVolume') + # if hasattr(self,'ReserVoirSimpleLocs'): ## -> add states that may be required if certain modules are used + # states.append('ReservoirVolume') - # if hasattr(self,'ReserVoirComplexLocs'): - # states.append('ReservoirWaterLevel') + # if hasattr(self,'ReserVoirComplexLocs'): + # states.append('ReservoirWaterLevel') - return states + return states - # The following are made to better connect to deltashell/openmi - def supplyCurrentTime(self): # - this may not be required - """ + def supplyCurrentTime(self): # - this may not be required + """ gets the current time in seconds after the start of the run Ouput: - time in seconds since the start of the model run """ - return self.currentTimeStep() * int(configget(self.config,'model','timestepsecs','86400')) + return self.currentTimeStep() * int( + configget(self.config, "model", "timestepsecs", "86400") + ) - def parameters(self): - """ + def parameters(self): + """ Define all model parameters here that the framework should handle for the model See wf_updateparameters and the parameters section of the ini file If you use this make sure to all wf_updateparameters at the start of the dynamic section and at the start/end of the initial section """ - modelparameters = [] + modelparameters = [] - #Static model parameters e.g. - #modelparameters.append(self.ParamType(name="RunoffGeneratingGWPerc",stack="intbl/RunoffGeneratingGWPerc.tbl",type="static",default=0.1)) + # Static model parameters e.g. + # modelparameters.append(self.ParamType(name="RunoffGeneratingGWPerc",stack="intbl/RunoffGeneratingGWPerc.tbl",type="static",default=0.1)) - # Meteo and other forcing + # Meteo and other forcing - self.Prec_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Prec","/inmaps/Prec") # timeseries for rainfall - self.Tair_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Tair","/inmaps/Tair") # timeseries for rainfall "/inmaps/TEMP" # global radiation - self.Tmax_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Tmax","/inmaps/Tmax") # timeseries for rainfall "/inmaps/TEMP" # global radiation - self.Tmin_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Tmin","/inmaps/Tmin") # timeseries for rainfall "/inmaps/TEMP" # global radiation - - - + self.Prec_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "Prec", "/inmaps/Prec" + ) # timeseries for rainfall + self.Tair_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "Tair", "/inmaps/Tair" + ) # timeseries for rainfall "/inmaps/TEMP" # global radiation + self.Tmax_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "Tmax", "/inmaps/Tmax" + ) # timeseries for rainfall "/inmaps/TEMP" # global radiation + self.Tmin_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "Tmin", "/inmaps/Tmin" + ) # timeseries for rainfall "/inmaps/TEMP" # global radiation - # Meteo and other forcing - modelparameters.append(self.ParamType(name="Prec",stack=self.Prec_mapstack,type="timeseries",default=0.0,verbose=True,lookupmaps=[])) - modelparameters.append(self.ParamType(name="Tair",stack=self.Tair_mapstack,type="timeseries",default=0.0,verbose=True,lookupmaps=[])) - modelparameters.append(self.ParamType(name="Tmax",stack=self.Tmax_mapstack,type="timeseries",default=0.0,verbose=True,lookupmaps=[])) - modelparameters.append(self.ParamType(name="Tmin",stack=self.Tmin_mapstack,type="timeseries",default=0.0,verbose=True,lookupmaps=[])) - - + # Meteo and other forcing + modelparameters.append( + self.ParamType( + name="Prec", + stack=self.Prec_mapstack, + type="timeseries", + default=0.0, + verbose=True, + lookupmaps=[], + ) + ) + modelparameters.append( + self.ParamType( + name="Tair", + stack=self.Tair_mapstack, + type="timeseries", + default=0.0, + verbose=True, + lookupmaps=[], + ) + ) + modelparameters.append( + self.ParamType( + name="Tmax", + stack=self.Tmax_mapstack, + type="timeseries", + default=0.0, + verbose=True, + lookupmaps=[], + ) + ) + modelparameters.append( + self.ParamType( + name="Tmin", + stack=self.Tmin_mapstack, + type="timeseries", + default=0.0, + verbose=True, + lookupmaps=[], + ) + ) + return modelparameters + def suspend(self): + """ + Suspends the model to disk. All variables needed to restart the model + are saved to disk as pcraster maps. Use resume() to re-read them + """ + self.logger.info("Saving initial conditions...") + self.wf_suspend(os.path.join(self.SaveDir, "outstate")) - return modelparameters + if self.OverWriteInit: + self.logger.info("Saving initial conditions over start conditions...") + self.wf_suspend(os.path.join(self.SaveDir, "instate")) + if self.fewsrun: + self.logger.info("Saving initial conditions for FEWS...") + self.wf_suspend(os.path.join(self.Dir, "outstate")) - def suspend(self): - """ - Suspends the model to disk. All variables needed to restart the model - are saved to disk as pcraster maps. Use resume() to re-read them - """ + def initial(self): + global statistics + global multpars + global updateCols - self.logger.info("Saving initial conditions...") - self.wf_suspend(os.path.join(self.SaveDir,"outstate")) + setglobaloption("unittrue") - if self.OverWriteInit: - self.logger.info("Saving initial conditions over start conditions...") - self.wf_suspend(os.path.join(self.SaveDir,"instate")) + self.thestep = scalar(0) + #: files to be used in case of timesries (scalar) input to the model - if self.fewsrun: - self.logger.info("Saving initial conditions for FEWS...") - self.wf_suspend(os.path.join(self.Dir, "outstate")) + # #: name of the tss file with precipitation data ("../intss/P.tss") + # self.precipTss = "../intss/P.tss" + # self.evapTss="../intss/PET.tss" #: name of the tss file with potential evap data ("../intss/PET.tss") + # self.tempTss="../intss/T.tss" #: name of the tss file with temperature data ("../intss/T.tss") + # self.inflowTss="../intss/Inflow.tss" #: NOT TESTED name of the tss file with inflow data ("../intss/Inflow.tss") + # self.SeepageTss="../intss/Seepage.tss" #: NOT TESTED name of the tss file with seepage data ("../intss/Seepage.tss")" + self.logger.info("running for " + str(self.nrTimeSteps()) + " timesteps") + # Set and get defaults from ConfigFile here ################################### + # self.scalarInput = int(configget(self.config,"model","ScalarInput","0")) + # self.Tslice = int(configget(self.config,"model","Tslice","1")) + # self.interpolMethod = configget(self.config,"model","InterpolationMethod","inv") + self.reinit = int(configget(self.config, "run", "reinit", "0")) + self.fewsrun = int(configget(self.config, "run", "fewsrun", "0")) + self.OverWriteInit = int(configget(self.config, "model", "OverWriteInit", "0")) + # self.updating = int(configget(self.config,"model","updating","0")) + # self.updateFile = configget(self.config,"model","updateFile","no_set") - def initial(self): + # self.sCatch = int(configget(self.config,"model","sCatch","0")) + # self.intbl = configget(self.config,"model","intbl","intbl") + # self.P_style = int(configget(self.config,"model","P_style","1")) + # self.PET_style = int(configget(self.config,"model","PET_style","1")) + # self.TEMP_style = int(configget(self.config,"model","TEMP_style","1")) - global statistics - global multpars - global updateCols + # self.modelSnow = int(configget(self.config,"model","ModelSnow","1")) + # sizeinmetres = int(configget(self.config,"layout","sizeinmetres","0")) + # alf = float(configget(self.config,"model","Alpha","60")) + # Qmax = float(configget(self.config,"model","AnnualDischarge","300")) + # self.UpdMaxDist =float(configget(self.config,"model","UpdMaxDist","100")) + # self.MaxUpdMult =float(configget(self.config,"model","MaxUpdMult","1.3")) + # self.MinUpdMult =float(configget(self.config,"model","MinUpdMult","0.7")) + # self.UpFrac =float(configget(self.config,"model","UpFrac","0.8")) + # self.ExternalQbase=int(configget(self.config,'model','ExternalQbase','0')) + # self.SetKquickFlow=int(configget(self.config,'model','SetKquickFlow','0')) + # self.MassWasting = int(configget(self.config,"model","MassWasting","0")) + # self.SubCatchFlowOnly = int(configget(self.config, 'model', 'SubCatchFlowOnly', '0')) - setglobaloption("unittrue") + # Print model info + print "The Spatial Processes in HYdrology (SPHY) model is " "developed and owned by FutureWater, Wageningen, The Netherlands" + print "Version 2.1" + print " " + # Read the modules to be used + self.GlacFLAG = int(configget(self.config, "MODULES", "GlacFLAG", "0")) + self.SnowFLAG = int(configget(self.config, "MODULES", "SnowFLAG", "0")) + self.RoutFLAG = int(configget(self.config, "MODULES", "RoutFLAG", "0")) + self.ResFLAG = int(configget(self.config, "MODULES", "ResFLAG", "0")) + self.LakeFLAG = int(configget(self.config, "MODULES", "LakeFLAG", "0")) + self.DynVegFLAG = int(configget(self.config, "MODULES", "DynVegFLAG", "0")) + self.GroundFLAG = int(configget(self.config, "MODULES", "GroundFLAG", "0")) - self.thestep = scalar(0) + # import the required modules + import datetime, calendar + from wflow.sphy import reporting as reporting + from wflow.sphy import timecalc as timecalc + from wflow.sphy import ET as ET + from wflow.sphy import rootzone as rootzone + from wflow.sphy import subzone as subzone - #: files to be used in case of timesries (scalar) input to the model + # from wflow.wflow_lib import * + from math import pi - # #: name of the tss file with precipitation data ("../intss/P.tss") - # self.precipTss = "../intss/P.tss" - # self.evapTss="../intss/PET.tss" #: name of the tss file with potential evap data ("../intss/PET.tss") - # self.tempTss="../intss/T.tss" #: name of the tss file with temperature data ("../intss/T.tss") - # self.inflowTss="../intss/Inflow.tss" #: NOT TESTED name of the tss file with inflow data ("../intss/Inflow.tss") - # self.SeepageTss="../intss/Seepage.tss" #: NOT TESTED name of the tss file with seepage data ("../intss/Seepage.tss")" + # -standard python modules + self.datetime = datetime + self.calendar = calendar + self.pi = pi + # -FW defined modules + self.reporting = reporting + self.timecalc = timecalc + self.ET = ET + self.rootzone = rootzone + self.subzone = subzone + del datetime, calendar, pi, reporting, timecalc, ET, rootzone, subzone + # -import additional modules if required + if self.GlacFLAG == 1: + self.SnowFLAG = 1 + self.GroundFLAG = 1 + import wflow.sphy.glacier as glacier # glacier melting processes + self.glacier = glacier + del glacier + if self.SnowFLAG == 1: + import wflow.sphy.snow as snow # snow melt processes - self.logger.info("running for " + str(self.nrTimeSteps()) + " timesteps") + self.snow = snow + del snow + if self.RoutFLAG == 1: + import wflow.sphy.routing as routing # simple routing scheme + self.routing = routing + del routing + if self.LakeFLAG == 1: + import lakes # import lake module - # Set and get defaults from ConfigFile here ################################### - # self.scalarInput = int(configget(self.config,"model","ScalarInput","0")) - # self.Tslice = int(configget(self.config,"model","Tslice","1")) - # self.interpolMethod = configget(self.config,"model","InterpolationMethod","inv") - self.reinit = int(configget(self.config,"run","reinit","0")) - self.fewsrun = int(configget(self.config,"run","fewsrun","0")) - self.OverWriteInit = int(configget(self.config,"model","OverWriteInit","0")) - # self.updating = int(configget(self.config,"model","updating","0")) - # self.updateFile = configget(self.config,"model","updateFile","no_set") + self.lakes = lakes + del lakes + if self.ResFLAG == 1: + import reservoirs # import reservoir module - # self.sCatch = int(configget(self.config,"model","sCatch","0")) - # self.intbl = configget(self.config,"model","intbl","intbl") - # self.P_style = int(configget(self.config,"model","P_style","1")) - # self.PET_style = int(configget(self.config,"model","PET_style","1")) - # self.TEMP_style = int(configget(self.config,"model","TEMP_style","1")) + self.reservoirs = reservoirs + del reservoirs + if self.LakeFLAG == 1 or self.ResFLAG == 1: + import advanced_routing # overwrite the simple routing scheme + self.routing = advanced_routing + del advanced_routing + self.RoutFLAG = 0 + if self.DynVegFLAG == 1: + import dynamic_veg # dynamic crop growth using ndvi or kc time-series - # self.modelSnow = int(configget(self.config,"model","ModelSnow","1")) - # sizeinmetres = int(configget(self.config,"layout","sizeinmetres","0")) - # alf = float(configget(self.config,"model","Alpha","60")) - # Qmax = float(configget(self.config,"model","AnnualDischarge","300")) - # self.UpdMaxDist =float(configget(self.config,"model","UpdMaxDist","100")) - # self.MaxUpdMult =float(configget(self.config,"model","MaxUpdMult","1.3")) - # self.MinUpdMult =float(configget(self.config,"model","MinUpdMult","0.7")) - # self.UpFrac =float(configget(self.config,"model","UpFrac","0.8")) - # self.ExternalQbase=int(configget(self.config,'model','ExternalQbase','0')) - # self.SetKquickFlow=int(configget(self.config,'model','SetKquickFlow','0')) - # self.MassWasting = int(configget(self.config,"model","MassWasting","0")) - # self.SubCatchFlowOnly = int(configget(self.config, 'model', 'SubCatchFlowOnly', '0')) + self.dynamic_veg = dynamic_veg + del dynamic_veg + if self.GroundFLAG == 1: + import wflow.sphy.groundwater as groundwater # groundwater storage as third storage layer. This is used instead of a fixed bottomflux - # Print model info - print 'The Spatial Processes in HYdrology (SPHY) model is ' \ - 'developed and owned by FutureWater, Wageningen, The Netherlands' - print 'Version 2.1' - print ' ' + self.groundwater = groundwater + del groundwater - # Read the modules to be used - self.GlacFLAG = int(configget(self.config,'MODULES','GlacFLAG','0')) - self.SnowFLAG = int(configget(self.config,'MODULES','SnowFLAG','0')) - self.RoutFLAG = int(configget(self.config,'MODULES','RoutFLAG','0')) - self.ResFLAG = int(configget(self.config,'MODULES','ResFLAG','0')) - self.LakeFLAG = int(configget(self.config,'MODULES','LakeFLAG','0')) - self.DynVegFLAG = int(configget(self.config,'MODULES','DynVegFLAG','0')) - self.GroundFLAG = int(configget(self.config,'MODULES','GroundFLAG','0')) - - # import the required modules - import datetime, calendar - from wflow.sphy import reporting as reporting - from wflow.sphy import timecalc as timecalc - from wflow.sphy import ET as ET - from wflow.sphy import rootzone as rootzone - from wflow.sphy import subzone as subzone + # -read the input and output directories from the configuration file + # self.inpath = config.get('DIRS', 'inputdir') + # self.inpathforcingT = config.get('DIRS','inputforcingdirT') + # self.inpathforcingP = config.get('DIRS','inputforcingdirP') + # self.outpath = config.get('DIRS', 'outputdir') - #from wflow.wflow_lib import * - from math import pi - #-standard python modules - self.datetime = datetime - self.calendar = calendar - self.pi = pi - #-FW defined modules - self.reporting = reporting - self.timecalc = timecalc - self.ET = ET - self.rootzone = rootzone - self.subzone = subzone - del datetime, calendar, pi, reporting, timecalc, ET, rootzone, subzone - #-import additional modules if required - if self.GlacFLAG == 1: - self.SnowFLAG = 1 - self.GroundFLAG = 1 - import wflow.sphy.glacier as glacier # glacier melting processes - self.glacier = glacier - del glacier - if self.SnowFLAG == 1: - import wflow.sphy.snow as snow # snow melt processes - self.snow = snow - del snow - if self.RoutFLAG == 1: - import wflow.sphy.routing as routing # simple routing scheme - self.routing = routing - del routing - if self.LakeFLAG == 1: - import lakes # import lake module - self.lakes = lakes - del lakes - if self.ResFLAG == 1: - import reservoirs # import reservoir module - self.reservoirs = reservoirs - del reservoirs - if self.LakeFLAG == 1 or self.ResFLAG == 1: - import advanced_routing # overwrite the simple routing scheme - self.routing = advanced_routing - del advanced_routing - self.RoutFLAG = 0 - if self.DynVegFLAG == 1: - import dynamic_veg # dynamic crop growth using ndvi or kc time-series - self.dynamic_veg = dynamic_veg - del dynamic_veg - if self.GroundFLAG == 1: - import wflow.sphy.groundwater as groundwater # groundwater storage as third storage layer. This is used instead of a fixed bottomflux - self.groundwater = groundwater - del groundwater - - #-read the input and output directories from the configuration file - # self.inpath = config.get('DIRS', 'inputdir') - # self.inpathforcingT = config.get('DIRS','inputforcingdirT') - # self.inpathforcingP = config.get('DIRS','inputforcingdirP') - # self.outpath = config.get('DIRS', 'outputdir') - - # self.starttime = configget(self.config,"run","starttime","0") - # ds = dt.datetime.strptime(self.starttime, '%Y-%m-%d %H:%M:%S %Z') - # self.endtime = configget(self.config,"run","endtime","0") - # de = dt.datetime.strptime(self.endtime, '%Y-%m-%d %H:%M:%S %Z') - - # #-set the timing criteria - # sy = ds.year - # sm = ds.month - # sd = ds.day - # ey = de.year - # em = de.month - # ed = de.day - # self.startdate = self.datetime.datetime(sy,sm,sd) - # self.enddate = self.datetime.datetime(ey,em,ed) - - # #-get start date of first forcing file in forcing directory - # syF = config.getint('TIMING', 'startyear_F') - # smF = config.getint('TIMING', 'startmonth_F') - # sdF = config.getint('TIMING', 'startday_F') - # self.startdateF = self.datetime.datetime(syF, smF, sdF) - - #-set the global options - setglobaloption('radians') - #-set the 2000 julian date number - self.julian_date_2000 = 2451545 - #-set the option to calculate the fluxes in mm for the upstream area - self.mm_rep_FLAG = int(configget(self.config,'REPORTING','mm_rep_FLAG','1')) - - # #-setting clone map - # clonemap = self.inpath + config.get('GENERAL','mask') ##->check - # setclone(clonemap) - # self.clone = readmap(clonemap) - - #-read general maps - self.DEM = readmap(os.path.join(self.Dir,"staticmaps", configget(self.config, 'GENERAL','dem','dem.map'))) #-> This has to be implemented for all readmap functions - self.Slope = readmap(os.path.join(self.Dir,"staticmaps", configget(self.config, 'GENERAL','Slope','slope.map'))) - self.Locations = readmap(os.path.join(self.Dir,"staticmaps", configget(self.config, 'GENERAL','locations','outlets.map'))) - - #-read soil maps - #self.Soil = readmap(self.inpath + config.get('SOIL','Soil')) - self.RootFieldMap = readmap(os.path.join(self.Dir,"staticmaps", configget(self.config,'SOIL','RootFieldMap','root_field.map'))) - self.RootSatMap = readmap(os.path.join(self.Dir,"staticmaps", configget(self.config,'SOIL','RootSatMap','root_sat.map'))) - self.RootDryMap = readmap(os.path.join(self.Dir,"staticmaps", configget(self.config,'SOIL','RootDryMap','root_dry.map'))) - self.RootWiltMap = readmap(os.path.join(self.Dir,"staticmaps", configget(self.config,'SOIL','RootWiltMap','root_wilt.map'))) - self.RootKsat = readmap(os.path.join(self.Dir,"staticmaps", configget(self.config,'SOIL','RootKsat','root_ksat.map'))) - self.SubSatMap = readmap(os.path.join(self.Dir,"staticmaps", configget(self.config,'SOIL','SubSatMap', 'sub_sat.map'))) - self.SubFieldMap = readmap(os.path.join(self.Dir,"staticmaps", configget(self.config,'SOIL','SubFieldMap','sub_field.map'))) - self.SubKsat = readmap(os.path.join(self.Dir,"staticmaps", configget(self.config,'SOIL','SubKsat','sub_ksat.map'))) - self.RootDrainVel = self.RootKsat * self.Slope - - #-Read and set the soil parameters - pars = ['CapRiseMax','RootDepthFlat','SubDepthFlat'] - for i in pars: - try: - #setattr(self, i, readmap(self.inpath + config.get('SOILPARS',i))) - setattr(self, i, readmap(os.path.join(self.Dir,"staticmaps", configget(self.config,'SOILPARS',i,i)))) - except: - #setattr(self, i, config.getfloat('SOILPARS',i)) - setattr(self, i, float(configget(self.config,'SOILPARS',i,i))) - if self.GroundFLAG == 0: # if groundwater module is not used, read seepage and gwl_base - self.SeepStatFLAG = config.getint('SOILPARS','SeepStatic') - if self.SeepStatFLAG == 0: # set the seepage map series - self.Seepmaps = self.inpath + config.get('SOILPARS', 'SeePage') - else: #-set a static map or value for seepage - try: - self.SeePage = readmap(self.inpath + config.get('SOILPARS','SeePage')) - except: - self.SeePage = config.getfloat('SOILPARS','SeePage') - try: - self.GWL_base = readmap(self.inpath + config.get('SOILPARS','GWL_base')) - except: - self.GWL_base = config.getfloat('SOILPARS','GWL_base') - - self.SubDrainVel = self.SubKsat * self.Slope - else: # if groundwater module is used, then read the groundwater soil parameters - pars = ['GwDepth','GwSat','deltaGw','BaseThresh','alphaGw','YieldGw'] - for i in pars: - try: - setattr(self, i, readmap(self.inpath + config.get('GROUNDW_PARS',i))) - except: - setattr(self, i, float(configget(self.config,'GROUNDW_PARS',i,i))) - - #-calculate soil properties - self.RootField = self.RootFieldMap * self.RootDepthFlat - self.RootSat = self.RootSatMap * self.RootDepthFlat - self.RootDry = self.RootDryMap * self.RootDepthFlat - self.RootWilt = self.RootWiltMap * self.RootDepthFlat - self.SubSat = self.SubSatMap * self.SubDepthFlat - self.SubField = self.SubFieldMap * self.SubDepthFlat - self.RootTT = (self.RootSat - self.RootField) / self.RootKsat - self.SubTT = (self.SubSat - self.SubField) / self.SubKsat - # soil max and soil min for scaling of gwl if groundwater module is not used - if self.GroundFLAG == 0: - self.SoilMax = self.RootSat + self.SubSat - self.SoilMin = self.RootDry + self.SubField - - #-read the crop coefficient table if the dynamic vegetation module is not used - if self.DynVegFLAG == 0: - self.KcStatFLAG = int(configget(self.config,'LANDUSE', 'KCstatic','kc.tbl')) - if self.KcStatFLAG == 1: - #-read land use map and kc table - self.LandUse = readmap(os.path.join(self.Dir,"staticmaps",configget(self.config,'LANDUSE','LandUse','landuse.map'))) - self.kc_table = os.path.join(self.Dir,"staticmaps",configget(self.config,'LANDUSE','CropFac','kc.tbl')) - self.Kc = lookupscalar(self.kc_table, self.LandUse) - else: - #-set the kc map series - self.Kcmaps = self.inpath + config.get('LANDUSE', 'KC') - #-Use the dynamic vegetation module - else: - #-set the ndvi map series to be read - self.ndvi = self.inpath + config.get('DYNVEG', 'NDVI') - #-read the vegetation parameters - pars = ['NDVImax','NDVImin','NDVIbase','KCmax','KCmin','LAImax','FPARmax','FPARmin'] - for i in pars: - try: - setattr(self, i, readmap(self.inpath + config.get('DYNVEG', i))) - except: - setattr(self, i, config.getfloat('DYNVEG', i)) - - #-read and set glacier maps and parameters if glacier module is used - if self.GlacFLAG == 1: - # self.GlacFracCI = readmap(self.inpath + config.get('GLACIER','GlacFracCI')) - # self.GlacFracDB = readmap(self.inpath + config.get('GLACIER','GlacFracDB')) - self.GlacFracCI = readmap(os.path.join(self.Dir,"staticmaps", configget(self.config,'GLACIER','GlacFracCI','glacier_clean.map'))) - self.GlacFracDB = readmap(os.path.join(self.Dir,"staticmaps", configget(self.config,'GLACIER','GlacFracDB','glacier_debris.map'))) - pars = ['DDFG','DDFDG','GlacF'] - for i in pars: - try: - setattr(self, i, readmap(self.inpath + config.get('GLACIER',i))) - except: - #setattr(self, i, config.getfloat('GLACIER',i)) - setattr(self, i, float(configget(self.config,'GLACIER',i,i))) - - #-read and set snow maps and parameters if snow modules are used - if self.SnowFLAG == 1: - pars = ['Tcrit','SnowSC','DDFS'] - for i in pars: - try: - setattr(self, i, readmap(self.inpath + config.get('SNOW',i))) - except: -# setattr(self, i, float(configget(self.config,'SNOW',i,i))) - setattr(self, i, float(configget(self.config,'SNOW',i,i))) - - #-read and set climate forcing and the calculation of etref - - self.Prec_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Prec","/inmaps/Prec") # timeseries for rainfall - self.Tair_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Tair","/inmaps/Tair") # timeseries for rainfall "/inmaps/TEMP" # global radiation - self.Tmax_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Tmax","/inmaps/Tmax") # timeseries for rainfall "/inmaps/TEMP" # global radiation - self.Tmin_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Tmin","/inmaps/Tmin") # timeseries for rainfall "/inmaps/TEMP" # global radiation - - - - - # self.Prec = self.inpathforcingP + config.get('CLIMATE','Prec') - # self.Tair = self.inpathforcingT + config.get('CLIMATE','Tair') - #self.ETREF_FLAG = config.getint('ETREF','ETREF_FLAG') ##-> for now should be zero. - self.ETREF_FLAG = int(configget(self.config,'ETREF','ETREF_FLAG',0)) ##-> for now should be zero. - #-determine the use of a given etref time-series or calculate etref using Hargreaves - if self.ETREF_FLAG == 1: - self.ETref = self.inpath + config.get('ETREF','ETref') - else: - #self.Lat = readmap(self.inpath + config.get('ETREF','Lat')) - self.Lat = readmap(os.path.join(self.Dir,"staticmaps",configget(self.config,'ETREF','Lat','latitude.map'))) - self.Tmax_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Tmax","/inmaps/Tmax") # timeseries for rainfall "/inmaps/TEMP" # global radiation - self.Tmin_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Tmin","/inmaps/Tmin") # timeseries for rainfall "/inmaps/TEMP" # global radiation - #self.Gsc = config.getfloat('ETREF', 'Gsc') - self.Gsc = float(configget(self.config,'ETREF','Gsc',0.0820)) - from wflow.sphy import hargreaves - self.Hargreaves = hargreaves - del hargreaves - - #-read and set routing maps and parameters - if self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1: - #self.FlowDir = readmap(self.inpath + config.get('ROUTING','flowdir')) - self.FlowDir = readmap(os.path.join(self.Dir,"staticmaps",configget(self.config,'ROUTING','flowdir','ldd.map'))) - try: - self.kx = readmap(self.inpath + config.get('ROUTING','kx')) - except: - #self.kx = config.getfloat('ROUTING','kx') - self.kx = float(configget(self.config,'ROUTING','kx',1)) - - setglobaloption('matrixtable') - #-read lake maps and parameters if lake module is used - if self.LakeFLAG == 1: - # nominal map with lake IDs - self.LakeID = cover(readmap(self.inpath + config.get('LAKE','LakeId')), 0) - # lookup table with function for each lake (exp, 1-order poly, 2-order poly, 3-order poly) - LakeFunc_Tab = self.inpath + config.get('LAKE', 'LakeFunc') - # lookup table with Qh-coeficients for each lake - LakeQH_Tab = self.inpath + config.get('LAKE', 'LakeQH') - # lookup table with Sh-coeficients for each lake - LakeSH_Tab = self.inpath + config.get('LAKE', 'LakeSH') - # lookup table with hS-coeficients for each lake - LakeHS_Tab = self.inpath + config.get('LAKE', 'LakeHS') - # create lake coefficient maps - self.LakeQH_Func = lookupnominal(LakeFunc_Tab, 1, self.LakeID) - self.LakeSH_Func = lookupnominal(LakeFunc_Tab, 2, self.LakeID) - self.LakeHS_Func = lookupnominal(LakeFunc_Tab, 3, self.LakeID) - # Read QH coefficients - self.LakeQH_exp_a = lookupscalar(LakeQH_Tab, 1, self.LakeID) - self.LakeQH_exp_b = lookupscalar(LakeQH_Tab, 2, self.LakeID) - self.LakeQH_pol_b = lookupscalar(LakeQH_Tab, 3, self.LakeID) - self.LakeQH_pol_a1 = lookupscalar(LakeQH_Tab, 4, self.LakeID) - self.LakeQH_pol_a2 = lookupscalar(LakeQH_Tab, 5, self.LakeID) - self.LakeQH_pol_a3 = lookupscalar(LakeQH_Tab, 6, self.LakeID) - # Read SH coefficients - self.LakeSH_exp_a = lookupscalar(LakeSH_Tab, 1, self.LakeID) - self.LakeSH_exp_b = lookupscalar(LakeSH_Tab, 2, self.LakeID) - self.LakeSH_pol_b = lookupscalar(LakeSH_Tab, 3, self.LakeID) - self.LakeSH_pol_a1 = lookupscalar(LakeSH_Tab, 4, self.LakeID) - self.LakeSH_pol_a2 = lookupscalar(LakeSH_Tab, 5, self.LakeID) - self.LakeSH_pol_a3 = lookupscalar(LakeSH_Tab, 6, self.LakeID) - # Read HS coefficients - self.LakeHS_exp_a = lookupscalar(LakeHS_Tab, 1, self.LakeID) - self.LakeHS_exp_b = lookupscalar(LakeHS_Tab, 2, self.LakeID) - self.LakeHS_pol_b = lookupscalar(LakeHS_Tab, 3, self.LakeID) - self.LakeHS_pol_a1 = lookupscalar(LakeHS_Tab, 4, self.LakeID) - self.LakeHS_pol_a2 = lookupscalar(LakeHS_Tab, 5, self.LakeID) - self.LakeHS_pol_a3 = lookupscalar(LakeHS_Tab, 6, self.LakeID) - #-read water level maps and parameters if available - try: - self.UpdateLakeLevel = readmap(self.inpath + config.get('LAKE','updatelakelevel')) - self.LLevel = self.inpath + config.get('LAKE','LakeFile') - print 'measured lake levels will be used to update lake storage' - except: - pass - - #-read reservior maps and parameters if reservoir module is used - if self.ResFLAG == 1: - # nominal map with reservoir IDs - self.ResID = cover(readmap(self.inpath + config.get('RESERVOIR','ResId')), 0) - # lookup table with operational scheme to use (simple or advanced) - ResFunc_Tab = self.inpath + config.get('RESERVOIR', 'ResFuncStor') - # Reservoir function - self.ResFunc = cover(lookupscalar(ResFunc_Tab, 1, self.ResID), 0) - try: - # lookup table with coefficients for simple reservoirs - ResSimple_Tab = self.inpath + config.get('RESERVOIR', 'ResSimple') - # Read coefficients for simple reservoirs - self.ResKr = lookupscalar(ResSimple_Tab, 1, self.ResID) - self.ResSmax = lookupscalar(ResSimple_Tab, 2, self.ResID) * 10**6 # convert to m3 - self.ResSimple = True - except: - self.ResSimple = False - try: - # lookup table with coefficients for advanced reservoirs - ResAdvanced_Tab = self.inpath + config.get('RESERVOIR', 'ResAdv') - # Read coefficients for advanced reservoirs - self.ResEVOL = lookupscalar(ResAdvanced_Tab, 1, self.ResID) * 10**6 # convert to m3 - self.ResPVOL = lookupscalar(ResAdvanced_Tab, 2, self.ResID) * 10**6 # convert to m3 - self.ResMaxFl = lookupscalar(ResAdvanced_Tab, 3, self.ResID) * 10**6 # convert to m3/d - self.ResDemFl = lookupscalar(ResAdvanced_Tab, 4, self.ResID) * 10**6 # convert to m3/d - self.ResFlStart = lookupscalar(ResAdvanced_Tab, 5, self.ResID) - self.ResFlEnd = lookupscalar(ResAdvanced_Tab, 6, self.ResID) - self.ResAdvanced = True - except: - self.ResAdvanced = False - - #-below is original initial part from sphy - - - #-get the correct forcing file number, depending on the start date of your simulation - #-and the start date of the first forcing file in your forcing directory. - #self.counter = (self.startdate - self.startdateF).days - # #-initial date - # self.curdate = self.startdate - #self.curdate = configget(self.config,"run","starttime","0") - #print self.curdate - #-initial soil properties - #-initial rootwater content + # self.starttime = configget(self.config,"run","starttime","0") + # ds = dt.datetime.strptime(self.starttime, '%Y-%m-%d %H:%M:%S %Z') + # self.endtime = configget(self.config,"run","endtime","0") + # de = dt.datetime.strptime(self.endtime, '%Y-%m-%d %H:%M:%S %Z') + # #-set the timing criteria + # sy = ds.year + # sm = ds.month + # sd = ds.day + # ey = de.year + # em = de.month + # ed = de.day + # self.startdate = self.datetime.datetime(sy,sm,sd) + # self.enddate = self.datetime.datetime(ey,em,ed) - self.RootWater = self.RootField - self.SubWater = self.SubField - - # if not config.get('SOIL_INIT','RootWater'): - # self.RootWater = self.RootField - # else: - # try: - # self.RootWater = config.getfloat('SOIL_INIT','RootWater') - # except: - # self.RootWater = readmap(self.inpath + config.get('SOIL_INIT','RootWater')) - # #-initial water content in subsoil - # if not config.get('SOIL_INIT','SubWater'): - # self.SubWater = self.SubField - # else: - # try: - # self.SubWater = config.getfloat('SOIL_INIT','SubWater') - # except: - # self.SubWater = readmap(self.inpath + config.get('SOIL_INIT','SubWater')) - #-initial water storage in rootzone + subsoil - self.SoilWater = self.RootWater + self.SubWater - #-initial capillary rise - self.CapRise = configget(self.config,'SOIL_INIT','CapRise',3) - # try: - # self.CapRise = config.getfloat('SOIL_INIT','CapRise') - # except: - # self.CapRise = readmap(self.inpath + config.get('SOIL_INIT','CapRise')) - #-initial drainage from rootzone - self.RootDrain = configget(self.config,'SOIL_INIT','RootDrain',3) - # try: - # self.RootDrain = config.getfloat('SOIL_INIT','RootDrain') - # except: - # self.RootDrain = readmap(self.inpath + config.get('SOIL_INIT','RootDrain')) - - if self.DynVegFLAG == 1: - #-initial canopy storage - self.Scanopy = 0 - #-initial ndvi if first map is not provided - self.ndviOld = scalar((self.NDVImax + self.NDVImin)/2) - elif self.KcStatFLAG == 0: - #-set initial kc value to one, if kc map is not available for first timestep - self.KcOld = scalar(1) - - #-initial groundwater properties - if self.GroundFLAG == 1: - self.GwRecharge = float(configget(self.config,'GROUNDW_INIT','GwRecharge',0)) - self.BaseR = float(configget(self.config,'GROUNDW_INIT','BaseR',1)) - self.Gw = float(configget(self.config,'GROUNDW_INIT','Gw',1)) - self.H_gw = float(configget(self.config,'GROUNDW_INIT','H_gw',1)) - # #-initial groundwater recharge - # try: - # self.GwRecharge = config.getfloat('GROUNDW_INIT','GwRecharge') - # except: - # self.GwRecharge = readmap(self.inpath + config.get('GROUNDW_INIT','GwRecharge')) - # #-initial baseflow - # try: - # self.BaseR = config.getfloat('GROUNDW_INIT','BaseR') - # except: - # self.BaseR = readmap(self.inpath + config.get('GROUNDW_INIT','BaseR')) - # #-initial groundwater storage - # try: - # self.Gw = config.getfloat('GROUNDW_INIT','Gw') - # except: - # self.Gw = readmap(self.inpath + config.get('GROUNDW_INIT','Gw')) - # #-initial groundwater level - # try: - # self.H_gw = config.getfloat('GROUNDW_INIT','H_gw') - # except: - # self.H_gw = readmap(self.inpath + config.get('GROUNDW_INIT','H_gw')) - # self.H_gw = max((self.RootDepthFlat + self.SubDepthFlat + self.GwDepth)/1000 - self.H_gw, 0) - # else: - # #-initial drainage from subsoil - # try: - # self.SubDrain = config.getfloat('SOIL_INIT','SubDrain') - # except: - # self.SubDrain = readmap(self.inpath + config.get('SOIL_INIT','SubDrain')) - # #-initial seepage value if seepage map series is used - # if self.SeepStatFLAG == 0: - # self.SeepOld = scalar(0) - - #-initial snow properties - if self.SnowFLAG == 1: - try: - #self.SnowStore = config.getfloat('SNOW_INIT','SnowIni') - self.SnowStore = float(configget(self.config,'SNOW_INIT','SnowIni',0)) - except: - self.SnowStore = readmap(self.inpath + config.get('SNOW_INIT','SnowIni')) - #-initial water stored in snowpack - try: - self.SnowWatStore = float(configget(self.config,'SNOW_INIT','SnowWatStore',0)) - #self.SnowWatStore = config.getfloat('SNOW_INIT','SnowWatStore') - except: - self.SnowWatStore = readmap(self.inpath + config.get('SNOW_INIT','SnowWatStore')) - self.TotalSnowStore = self.SnowStore + self.SnowWatStore - - #-initial glacier properties - if self.GlacFLAG == 1: + # #-get start date of first forcing file in forcing directory + # syF = config.getint('TIMING', 'startyear_F') + # smF = config.getint('TIMING', 'startmonth_F') + # sdF = config.getint('TIMING', 'startday_F') + # self.startdateF = self.datetime.datetime(syF, smF, sdF) + + # -set the global options + setglobaloption("radians") + # -set the 2000 julian date number + self.julian_date_2000 = 2451545 + # -set the option to calculate the fluxes in mm for the upstream area + self.mm_rep_FLAG = int(configget(self.config, "REPORTING", "mm_rep_FLAG", "1")) + + # #-setting clone map + # clonemap = self.inpath + config.get('GENERAL','mask') ##->check + # setclone(clonemap) + # self.clone = readmap(clonemap) + + # -read general maps + self.DEM = readmap( + os.path.join( + self.Dir, + "staticmaps", + configget(self.config, "GENERAL", "dem", "dem.map"), + ) + ) # -> This has to be implemented for all readmap functions + self.Slope = readmap( + os.path.join( + self.Dir, + "staticmaps", + configget(self.config, "GENERAL", "Slope", "slope.map"), + ) + ) + self.Locations = readmap( + os.path.join( + self.Dir, + "staticmaps", + configget(self.config, "GENERAL", "locations", "outlets.map"), + ) + ) + + # -read soil maps + # self.Soil = readmap(self.inpath + config.get('SOIL','Soil')) + self.RootFieldMap = readmap( + os.path.join( + self.Dir, + "staticmaps", + configget(self.config, "SOIL", "RootFieldMap", "root_field.map"), + ) + ) + self.RootSatMap = readmap( + os.path.join( + self.Dir, + "staticmaps", + configget(self.config, "SOIL", "RootSatMap", "root_sat.map"), + ) + ) + self.RootDryMap = readmap( + os.path.join( + self.Dir, + "staticmaps", + configget(self.config, "SOIL", "RootDryMap", "root_dry.map"), + ) + ) + self.RootWiltMap = readmap( + os.path.join( + self.Dir, + "staticmaps", + configget(self.config, "SOIL", "RootWiltMap", "root_wilt.map"), + ) + ) + self.RootKsat = readmap( + os.path.join( + self.Dir, + "staticmaps", + configget(self.config, "SOIL", "RootKsat", "root_ksat.map"), + ) + ) + self.SubSatMap = readmap( + os.path.join( + self.Dir, + "staticmaps", + configget(self.config, "SOIL", "SubSatMap", "sub_sat.map"), + ) + ) + self.SubFieldMap = readmap( + os.path.join( + self.Dir, + "staticmaps", + configget(self.config, "SOIL", "SubFieldMap", "sub_field.map"), + ) + ) + self.SubKsat = readmap( + os.path.join( + self.Dir, + "staticmaps", + configget(self.config, "SOIL", "SubKsat", "sub_ksat.map"), + ) + ) + self.RootDrainVel = self.RootKsat * self.Slope + + # -Read and set the soil parameters + pars = ["CapRiseMax", "RootDepthFlat", "SubDepthFlat"] + for i in pars: + try: + # setattr(self, i, readmap(self.inpath + config.get('SOILPARS',i))) + setattr( + self, + i, + readmap( + os.path.join( + self.Dir, + "staticmaps", + configget(self.config, "SOILPARS", i, i), + ) + ), + ) + except: + # setattr(self, i, config.getfloat('SOILPARS',i)) + setattr(self, i, float(configget(self.config, "SOILPARS", i, i))) + if ( + self.GroundFLAG == 0 + ): # if groundwater module is not used, read seepage and gwl_base + self.SeepStatFLAG = config.getint("SOILPARS", "SeepStatic") + if self.SeepStatFLAG == 0: # set the seepage map series + self.Seepmaps = self.inpath + config.get("SOILPARS", "SeePage") + else: # -set a static map or value for seepage + try: + self.SeePage = readmap( + self.inpath + config.get("SOILPARS", "SeePage") + ) + except: + self.SeePage = config.getfloat("SOILPARS", "SeePage") + try: + self.GWL_base = readmap( + self.inpath + config.get("SOILPARS", "GWL_base") + ) + except: + self.GWL_base = config.getfloat("SOILPARS", "GWL_base") + + self.SubDrainVel = self.SubKsat * self.Slope + else: # if groundwater module is used, then read the groundwater soil parameters + pars = ["GwDepth", "GwSat", "deltaGw", "BaseThresh", "alphaGw", "YieldGw"] + for i in pars: + try: + setattr( + self, i, readmap(self.inpath + config.get("GROUNDW_PARS", i)) + ) + except: + setattr( + self, i, float(configget(self.config, "GROUNDW_PARS", i, i)) + ) + + # -calculate soil properties + self.RootField = self.RootFieldMap * self.RootDepthFlat + self.RootSat = self.RootSatMap * self.RootDepthFlat + self.RootDry = self.RootDryMap * self.RootDepthFlat + self.RootWilt = self.RootWiltMap * self.RootDepthFlat + self.SubSat = self.SubSatMap * self.SubDepthFlat + self.SubField = self.SubFieldMap * self.SubDepthFlat + self.RootTT = (self.RootSat - self.RootField) / self.RootKsat + self.SubTT = (self.SubSat - self.SubField) / self.SubKsat + # soil max and soil min for scaling of gwl if groundwater module is not used + if self.GroundFLAG == 0: + self.SoilMax = self.RootSat + self.SubSat + self.SoilMin = self.RootDry + self.SubField + + # -read the crop coefficient table if the dynamic vegetation module is not used + if self.DynVegFLAG == 0: + self.KcStatFLAG = int( + configget(self.config, "LANDUSE", "KCstatic", "kc.tbl") + ) + if self.KcStatFLAG == 1: + # -read land use map and kc table + self.LandUse = readmap( + os.path.join( + self.Dir, + "staticmaps", + configget(self.config, "LANDUSE", "LandUse", "landuse.map"), + ) + ) + self.kc_table = os.path.join( + self.Dir, + "staticmaps", + configget(self.config, "LANDUSE", "CropFac", "kc.tbl"), + ) + self.Kc = lookupscalar(self.kc_table, self.LandUse) + else: + # -set the kc map series + self.Kcmaps = self.inpath + config.get("LANDUSE", "KC") + # -Use the dynamic vegetation module + else: + # -set the ndvi map series to be read + self.ndvi = self.inpath + config.get("DYNVEG", "NDVI") + # -read the vegetation parameters + pars = [ + "NDVImax", + "NDVImin", + "NDVIbase", + "KCmax", + "KCmin", + "LAImax", + "FPARmax", + "FPARmin", + ] + for i in pars: + try: + setattr(self, i, readmap(self.inpath + config.get("DYNVEG", i))) + except: + setattr(self, i, config.getfloat("DYNVEG", i)) + + # -read and set glacier maps and parameters if glacier module is used + if self.GlacFLAG == 1: + # self.GlacFracCI = readmap(self.inpath + config.get('GLACIER','GlacFracCI')) + # self.GlacFracDB = readmap(self.inpath + config.get('GLACIER','GlacFracDB')) + self.GlacFracCI = readmap( + os.path.join( + self.Dir, + "staticmaps", + configget( + self.config, "GLACIER", "GlacFracCI", "glacier_clean.map" + ), + ) + ) + self.GlacFracDB = readmap( + os.path.join( + self.Dir, + "staticmaps", + configget( + self.config, "GLACIER", "GlacFracDB", "glacier_debris.map" + ), + ) + ) + pars = ["DDFG", "DDFDG", "GlacF"] + for i in pars: + try: + setattr(self, i, readmap(self.inpath + config.get("GLACIER", i))) + except: + # setattr(self, i, config.getfloat('GLACIER',i)) + setattr(self, i, float(configget(self.config, "GLACIER", i, i))) + + # -read and set snow maps and parameters if snow modules are used + if self.SnowFLAG == 1: + pars = ["Tcrit", "SnowSC", "DDFS"] + for i in pars: + try: + setattr(self, i, readmap(self.inpath + config.get("SNOW", i))) + except: + # setattr(self, i, float(configget(self.config,'SNOW',i,i))) + setattr(self, i, float(configget(self.config, "SNOW", i, i))) + + # -read and set climate forcing and the calculation of etref + + self.Prec_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "Prec", "/inmaps/Prec" + ) # timeseries for rainfall + self.Tair_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "Tair", "/inmaps/Tair" + ) # timeseries for rainfall "/inmaps/TEMP" # global radiation + self.Tmax_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "Tmax", "/inmaps/Tmax" + ) # timeseries for rainfall "/inmaps/TEMP" # global radiation + self.Tmin_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "Tmin", "/inmaps/Tmin" + ) # timeseries for rainfall "/inmaps/TEMP" # global radiation + + # self.Prec = self.inpathforcingP + config.get('CLIMATE','Prec') + # self.Tair = self.inpathforcingT + config.get('CLIMATE','Tair') + # self.ETREF_FLAG = config.getint('ETREF','ETREF_FLAG') ##-> for now should be zero. + self.ETREF_FLAG = int( + configget(self.config, "ETREF", "ETREF_FLAG", 0) + ) ##-> for now should be zero. + # -determine the use of a given etref time-series or calculate etref using Hargreaves + if self.ETREF_FLAG == 1: + self.ETref = self.inpath + config.get("ETREF", "ETref") + else: + # self.Lat = readmap(self.inpath + config.get('ETREF','Lat')) + self.Lat = readmap( + os.path.join( + self.Dir, + "staticmaps", + configget(self.config, "ETREF", "Lat", "latitude.map"), + ) + ) + self.Tmax_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "Tmax", "/inmaps/Tmax" + ) # timeseries for rainfall "/inmaps/TEMP" # global radiation + self.Tmin_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "Tmin", "/inmaps/Tmin" + ) # timeseries for rainfall "/inmaps/TEMP" # global radiation + # self.Gsc = config.getfloat('ETREF', 'Gsc') + self.Gsc = float(configget(self.config, "ETREF", "Gsc", 0.0820)) + from wflow.sphy import hargreaves + + self.Hargreaves = hargreaves + del hargreaves + + # -read and set routing maps and parameters + if self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1: + # self.FlowDir = readmap(self.inpath + config.get('ROUTING','flowdir')) + self.FlowDir = readmap( + os.path.join( + self.Dir, + "staticmaps", + configget(self.config, "ROUTING", "flowdir", "ldd.map"), + ) + ) + try: + self.kx = readmap(self.inpath + config.get("ROUTING", "kx")) + except: + # self.kx = config.getfloat('ROUTING','kx') + self.kx = float(configget(self.config, "ROUTING", "kx", 1)) + + setglobaloption("matrixtable") + # -read lake maps and parameters if lake module is used + if self.LakeFLAG == 1: + # nominal map with lake IDs + self.LakeID = cover(readmap(self.inpath + config.get("LAKE", "LakeId")), 0) + # lookup table with function for each lake (exp, 1-order poly, 2-order poly, 3-order poly) + LakeFunc_Tab = self.inpath + config.get("LAKE", "LakeFunc") + # lookup table with Qh-coeficients for each lake + LakeQH_Tab = self.inpath + config.get("LAKE", "LakeQH") + # lookup table with Sh-coeficients for each lake + LakeSH_Tab = self.inpath + config.get("LAKE", "LakeSH") + # lookup table with hS-coeficients for each lake + LakeHS_Tab = self.inpath + config.get("LAKE", "LakeHS") + # create lake coefficient maps + self.LakeQH_Func = lookupnominal(LakeFunc_Tab, 1, self.LakeID) + self.LakeSH_Func = lookupnominal(LakeFunc_Tab, 2, self.LakeID) + self.LakeHS_Func = lookupnominal(LakeFunc_Tab, 3, self.LakeID) + # Read QH coefficients + self.LakeQH_exp_a = lookupscalar(LakeQH_Tab, 1, self.LakeID) + self.LakeQH_exp_b = lookupscalar(LakeQH_Tab, 2, self.LakeID) + self.LakeQH_pol_b = lookupscalar(LakeQH_Tab, 3, self.LakeID) + self.LakeQH_pol_a1 = lookupscalar(LakeQH_Tab, 4, self.LakeID) + self.LakeQH_pol_a2 = lookupscalar(LakeQH_Tab, 5, self.LakeID) + self.LakeQH_pol_a3 = lookupscalar(LakeQH_Tab, 6, self.LakeID) + # Read SH coefficients + self.LakeSH_exp_a = lookupscalar(LakeSH_Tab, 1, self.LakeID) + self.LakeSH_exp_b = lookupscalar(LakeSH_Tab, 2, self.LakeID) + self.LakeSH_pol_b = lookupscalar(LakeSH_Tab, 3, self.LakeID) + self.LakeSH_pol_a1 = lookupscalar(LakeSH_Tab, 4, self.LakeID) + self.LakeSH_pol_a2 = lookupscalar(LakeSH_Tab, 5, self.LakeID) + self.LakeSH_pol_a3 = lookupscalar(LakeSH_Tab, 6, self.LakeID) + # Read HS coefficients + self.LakeHS_exp_a = lookupscalar(LakeHS_Tab, 1, self.LakeID) + self.LakeHS_exp_b = lookupscalar(LakeHS_Tab, 2, self.LakeID) + self.LakeHS_pol_b = lookupscalar(LakeHS_Tab, 3, self.LakeID) + self.LakeHS_pol_a1 = lookupscalar(LakeHS_Tab, 4, self.LakeID) + self.LakeHS_pol_a2 = lookupscalar(LakeHS_Tab, 5, self.LakeID) + self.LakeHS_pol_a3 = lookupscalar(LakeHS_Tab, 6, self.LakeID) + # -read water level maps and parameters if available + try: + self.UpdateLakeLevel = readmap( + self.inpath + config.get("LAKE", "updatelakelevel") + ) + self.LLevel = self.inpath + config.get("LAKE", "LakeFile") + print "measured lake levels will be used to update lake storage" + except: + pass + + # -read reservior maps and parameters if reservoir module is used + if self.ResFLAG == 1: + # nominal map with reservoir IDs + self.ResID = cover( + readmap(self.inpath + config.get("RESERVOIR", "ResId")), 0 + ) + # lookup table with operational scheme to use (simple or advanced) + ResFunc_Tab = self.inpath + config.get("RESERVOIR", "ResFuncStor") + # Reservoir function + self.ResFunc = cover(lookupscalar(ResFunc_Tab, 1, self.ResID), 0) + try: + # lookup table with coefficients for simple reservoirs + ResSimple_Tab = self.inpath + config.get("RESERVOIR", "ResSimple") + # Read coefficients for simple reservoirs + self.ResKr = lookupscalar(ResSimple_Tab, 1, self.ResID) + self.ResSmax = ( + lookupscalar(ResSimple_Tab, 2, self.ResID) * 10 ** 6 + ) # convert to m3 + self.ResSimple = True + except: + self.ResSimple = False + try: + # lookup table with coefficients for advanced reservoirs + ResAdvanced_Tab = self.inpath + config.get("RESERVOIR", "ResAdv") + # Read coefficients for advanced reservoirs + self.ResEVOL = ( + lookupscalar(ResAdvanced_Tab, 1, self.ResID) * 10 ** 6 + ) # convert to m3 + self.ResPVOL = ( + lookupscalar(ResAdvanced_Tab, 2, self.ResID) * 10 ** 6 + ) # convert to m3 + self.ResMaxFl = ( + lookupscalar(ResAdvanced_Tab, 3, self.ResID) * 10 ** 6 + ) # convert to m3/d + self.ResDemFl = ( + lookupscalar(ResAdvanced_Tab, 4, self.ResID) * 10 ** 6 + ) # convert to m3/d + self.ResFlStart = lookupscalar(ResAdvanced_Tab, 5, self.ResID) + self.ResFlEnd = lookupscalar(ResAdvanced_Tab, 6, self.ResID) + self.ResAdvanced = True + except: + self.ResAdvanced = False + + # -below is original initial part from sphy + + # -get the correct forcing file number, depending on the start date of your simulation + # -and the start date of the first forcing file in your forcing directory. + # self.counter = (self.startdate - self.startdateF).days + # #-initial date + # self.curdate = self.startdate + # self.curdate = configget(self.config,"run","starttime","0") + # print self.curdate + # -initial soil properties + # -initial rootwater content + + self.RootWater = self.RootField + self.SubWater = self.SubField + + # if not config.get('SOIL_INIT','RootWater'): + # self.RootWater = self.RootField + # else: # try: - # self.GlacFrac = config.getfloat('GLACIER_INIT','GlacFrac') - # except: - # self.GlacFrac = readmap(self.inpath + config.get('GLACIER_INIT','GlacFrac')) - self.GlacFrac = readmap(os.path.join(self.Dir,"staticmaps",configget(self.config,'GLACIER_INIT','GlacFrac','glacierfraction.map'))) - print self.GlacFrac - #-initial routed total runoff and of individual components - if self.RoutFLAG == 1 or self.LakeFLAG==1 or self.ResFLAG==1: - #-initial routed total runoff - try: - self.QRAold = config.getfloat('ROUT_INIT','QRA_init') - except: - try: - self.QRAold = readmap(self.inpath + config.get('ROUT_INIT','QRA_init')) - except: - self.QRAold = 0 - #-initial routed runoff for the individual components - pars = ['RainRA','SnowRA','GlacRA','BaseRA'] - self.RainRAold = 0 - self.SnowRAold = 0 - self.GlacRAold = 0 - self.BaseRAold = 0 - self.RainRA_FLAG = True - self.SnowRA_FLAG = True - self.GlacRA_FLAG = True - self.BaseRA_FLAG = True - # for i in pars: - # try: - # setattr(self, i + 'old', readmap(self.inpath + config.get('ROUT_INIT', i + '_init'))) - # setattr(self, i + '_FLAG', True) - # except: - # try: - # #setattr(self, i + 'old', config.getfloat('ROUT_INIT', i + '_init')) - # setattr(self, i + '_FLAG', True) - # print RainRA_init - # except: - # setattr(self, i + '_FLAG', False) - - #-initial storage in lakes and reservoirs - if self.LakeFLAG == 1 or self.ResFLAG == 1: - #-Read initial storages from table/reservoir file - if self.LakeFLAG == 1: - LakeStor_Tab = self.inpath + config.get('LAKE', 'LakeStor') - self.StorRES = cover(lookupscalar(LakeStor_Tab, 1, self.LakeID), 0) * 10**6 # convert to m3 - #-Qfrac for lake cells should be zero, else 1 - self.QFRAC = ifthenelse(self.LakeID != 0, scalar(0), 1) - if self.ResFLAG == 1: - ResStor_Tab = self.inpath + config.get('RESERVOIR', 'ResFuncStor') - ResStor = cover(lookupscalar(ResStor_Tab, 2, self.ResID), 0) * 10**6 # convert to m3 - try: - self.StorRES = self.StorRES + ResStor - #-Qfrac for reservoir cells should be zero, else 1 - self.QFRAC = ifthenelse(self.ResID != 0, scalar(0), self.QFRAC) - except: - self.StorRES = ResStor - #-Qfrac for reservoir cells should be zero, else 1 - self.QFRAC = ifthenelse(self.ResID != 0, scalar(0), 1) - - #-initial storage in lakes/reservoirs of individual flow components - pars = ['RainRA','SnowRA','GlacRA','BaseRA'] - for i in pars: - column = pars.index(i) # identify column to be read from lake or reservoir table - try: #-try to sum the storages read from the lake and reservoir tables if both thse modules are used - setattr(self, i + 'stor', (cover(lookupscalar(LakeStor_Tab, column + 2, self.LakeID), 0) + \ - cover(lookupscalar(ResStor_Tab, column + 3, self.ResID), 0)) * 10**6) - if eval('self.' + i + '_FLAG'): - setattr(self, i + '_FLAG', True) - else: - setattr(self, i + '_FLAG', False) - except: - try: #-try to read the storages from the lake table - setattr(self, i + 'stor', cover(lookupscalar(LakeStor_Tab, column + 2, self.LakeID), 0) * 10**6) - if eval('self.' + i + '_FLAG'): - setattr(self, i + '_FLAG', True) - else: - setattr(self, i + '_FLAG', False) - except: #-try to read the storages from the reservoir table - try: - setattr(self, i + 'stor', cover(lookupscalar(ResStor_Tab, column + 3, self.ResID), 0) * 10**6) - if eval('self.' + i + '_FLAG'): - setattr(self, i + '_FLAG', True) - else: - setattr(self, i + '_FLAG', False) - except: - setattr(self, i + '_FLAG', False) - - #-Initial values for reporting and setting of time-series - #-set time-series reporting for mm flux from upstream area for prec and eta - if self.mm_rep_FLAG == 1 and (self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1): - self.PrecSubBasinTSS = TimeoutputTimeseries("PrecSubBasinTSS", self, self.Locations, noHeader=False) - self.ETaSubBasinTSS = TimeoutputTimeseries("ETaSubBasinTSS", self, self.Locations, noHeader=False) - if self.GlacFLAG == 1: - pars = ['wbal','GWL','TotPrec','TotPrecF','TotPrecEF','TotIntF','TotRain','TotRainF','TotETpot','TotETpotF','TotETact','TotETactF','TotSnow','TotSnowF','TotSnowMelt','TotSnowMeltF','TotGlacMelt','TotGlacMeltF','TotRootRF','TotRootDF','TotRootPF',\ - 'TotSubPF','TotCapRF','TotGlacPercF','TotGwRechargeF','TotRainRF','TotBaseRF','TotSnowRF','TotGlacRF','TotRF','RainRAtot','SnowRAtot','GlacRAtot','BaseRAtot','QallRAtot'] - #-set time-series reporting for mm fluxes from upstream area if glacier and routing/reservoir modules are used - if self.mm_rep_FLAG == 1 and (self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1): - self.GMeltSubBasinTSS = TimeoutputTimeseries("GMeltSubBasinTSS", self, self.Locations, noHeader=False) - self.QSNOWSubBasinTSS = TimeoutputTimeseries("QSNOWSubBasinTSS", self, self.Locations, noHeader=False) - self.QRAINSubBasinTSS = TimeoutputTimeseries("QRAINSubBasinTSS", self, self.Locations, noHeader=False) - self.QGLACSubBasinTSS = TimeoutputTimeseries("QGLACSubBasinTSS", self, self.Locations, noHeader=False) - self.QBASFSubBasinTSS = TimeoutputTimeseries("QBASFSubBasinTSS", self, self.Locations, noHeader=False) - self.QTOTSubBasinTSS = TimeoutputTimeseries("QTOTSubBasinTSS", self, self.Locations, noHeader=False) - elif self.SnowFLAG == 1: - if self.GroundFLAG == 1: - pars = ['wbal','GWL','TotPrec','TotPrecF','TotPrecEF','TotIntF','TotRain','TotRainF','TotETpot','TotETpotF','TotETact','TotETactF','TotSnow','TotSnowF','TotSnowMelt','TotSnowMeltF','TotRootRF','TotRootDF','TotRootPF',\ - 'TotSubPF','TotCapRF','TotGwRechargeF','TotRainRF','TotBaseRF','TotSnowRF','TotRF','RainRAtot','SnowRAtot','BaseRAtot','QallRAtot'] - #-set time-series reporting for mm fluxes from upstream area if snow, groundwater and routing/reservoir modules are used - if self.mm_rep_FLAG == 1 and (self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1): - self.QSNOWSubBasinTSS = TimeoutputTimeseries("QSNOWSubBasinTSS", self, self.Locations, noHeader=False) - self.QRAINSubBasinTSS = TimeoutputTimeseries("QRAINSubBasinTSS", self, self.Locations, noHeader=False) - self.QBASFSubBasinTSS = TimeoutputTimeseries("QBASFSubBasinTSS", self, self.Locations, noHeader=False) - self.QTOTSubBasinTSS = TimeoutputTimeseries("QTOTSubBasinTSS", self, self.Locations, noHeader=False) - else: - pars = ['wbal','GWL','TotPrec','TotPrecF','TotPrecEF','TotIntF','TotRain','TotRainF','TotETpot','TotETpotF','TotETact','TotETactF','TotSnow','TotSnowF','TotSnowMelt','TotSnowMeltF','TotRootRF','TotRootDF','TotRootPF',\ - 'TotSubDF','TotCapRF','TotSeepF','TotRainRF','TotSnowRF','TotRF','RainRAtot','SnowRAtot','BaseRAtot','QallRAtot'] - #-set time-series reporting for mm fluxes from upstream area if snow and routing/reservoir modules are used - if self.mm_rep_FLAG == 1 and (self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1): - self.SeepSubBasinTSS = TimeoutputTimeseries("SeepSubBasinTSS", self, self.Locations, noHeader=False) - self.QSNOWSubBasinTSS = TimeoutputTimeseries("QSNOWSubBasinTSS", self, self.Locations, noHeader=False) - self.QRAINSubBasinTSS = TimeoutputTimeseries("QRAINSubBasinTSS", self, self.Locations, noHeader=False) - self.QBASFSubBasinTSS = TimeoutputTimeseries("QBASFSubBasinTSS", self, self.Locations, noHeader=False) - self.QTOTSubBasinTSS = TimeoutputTimeseries("QTOTSubBasinTSS", self, self.Locations, noHeader=False) - else: - if self.GroundFLAG == 1: - pars = ['wbal','GWL','TotPrec','TotPrecF','TotPrecEF','TotIntF','TotRain','TotRainF','TotETpot','TotETpotF','TotETact','TotETactF','TotRootRF','TotRootDF','TotRootPF',\ - 'TotSubPF','TotCapRF','TotGwRechargeF','TotRainRF','TotBaseRF','TotRF','RainRAtot','BaseRAtot','QallRAtot'] - #-set time-series reporting for mm fluxes from upstream area if groundwater and routing/reservoir modules are used - if self.mm_rep_FLAG == 1 and (self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1): - self.QRAINSubBasinTSS = TimeoutputTimeseries("QRAINSubBasinTSS", self, self.Locations, noHeader=False) - self.QBASFSubBasinTSS = TimeoutputTimeseries("QBASFSubBasinTSS", self, self.Locations, noHeader=False) - self.QTOTSubBasinTSS = TimeoutputTimeseries("QTOTSubBasinTSS", self, self.Locations, noHeader=False) - else: - pars = ['wbal','GWL','TotPrec','TotPrecF','TotPrecEF','TotIntF','TotRain','TotRainF','TotETpot','TotETpotF','TotETact','TotETactF','TotRootRF','TotRootDF','TotRootPF',\ - 'TotSubDF','TotCapRF','TotSeepF','TotRainRF','TotRF','RainRAtot','BaseRAtot','QallRAtot'] - #-set time-series reporting for mm fluxes from upstream area if routing/reservoir modules are used - if self.mm_rep_FLAG == 1 and (self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1): - self.SeepSubBasinTSS = TimeoutputTimeseries("SeepSubBasinTSS", self, self.Locations, noHeader=False) - self.QRAINSubBasinTSS = TimeoutputTimeseries("QRAINSubBasinTSS", self, self.Locations, noHeader=False) - self.QBASFSubBasinTSS = TimeoutputTimeseries("QBASFSubBasinTSS", self, self.Locations, noHeader=False) - self.QTOTSubBasinTSS = TimeoutputTimeseries("QTOTSubBasinTSS", self, self.Locations, noHeader=False) - #-remove routing output from reported list of parameters if these modules are not used - if self.RoutFLAG == 0 and self.ResFLAG == 0 and self.LakeFLAG == 0: - rpars = ['RainRAtot','SnowRAtot','GlacRAtot','BaseRAtot','QallRAtot'] - for i in rpars: - try: - j = pars.index(i) - del pars[j] - except: - pass - #-set reporting options and read initial values - for i in pars: - mapoutops = configget(self.config,'REPORTING',i+'_mapoutput',i+'_mapoutput') - #mapoutops = config.get('REPORTING', i+'_mapoutput') - TSoutops = configget(self.config,'REPORTING',i+'_TSoutput',i+'_TSoutput') - #TSoutops = config.get('REPORTING', i+'_TSoutput') - if mapoutops == 'NONE' and TSoutops == 'NONE': - print i + ' will NOT be reported' - else: - print i + ' will be reported' - fname = configget(self.config,'REPORTING', i+'_fname', i+'_fname') - #fname = config.get('REPORTING', i+'_fname') - setattr(self, i+'_fname', fname) - # try: - # setattr(self, i, readmap(self.inpath + config.get('INITTOT', i))) - # except: - # try: - # setattr(self, i, config.getfloat('INITTOT', i)) - # except: - # setattr(self, i, 0.) - setattr(self, i, 0.) # use this instead of the commented part above, because it is more logical to always zero as initial condition for reporting - if mapoutops != 'NONE': - mapoutops = mapoutops.split(",") - for j in mapoutops: - if j == 'D': - setattr(self, i+'_Day', eval('self.'+i)) - setattr(self, i+'_Day_map', 1) - elif j == 'M': - setattr(self, i+'_Month', eval('self.'+i)) - setattr(self, i+'_Month_map', 1) - elif j == 'Y': - setattr(self, i+'_Year', eval('self.'+i)) - setattr(self, i+'_Year_map', 1) - else: - setattr(self, i+'_Final', eval('self.'+i)) - setattr(self, i+'_Final_map', 1) - if TSoutops != 'NONE': - TSoutops = TSoutops.split(",") - for j in TSoutops: - if j == 'D': - setattr(self, i+'_Day', eval('self.'+i)) - setattr(self, i+'_DayTS', eval('TimeoutputTimeseries("'+fname+'DTS'+'", self, self.Locations, noHeader=False)')) - elif j == 'M': - setattr(self, i+'_Month', eval('self.'+i)) - setattr(self, i+'_MonthTS', eval('TimeoutputTimeseries("'+fname+'MTS'+'", self, self.Locations, noHeader=False)')) - elif j == 'Y': - setattr(self, i+'_Year', eval('self.'+i)) - setattr(self, i+'_YearTS', eval('TimeoutputTimeseries("'+fname+'YTS'+'", self, self.Locations, noHeader=False)')) - - #-set reporting of water balances for lakes - if self.LakeFLAG == 1 and config.getint('REPORTING', 'Lake_wbal') ==1: - self.LakeInTSS = pcrm.TimeoutputTimeseries("LakeInTSS", self, self.LakeID, noHeader=True) - self.LakeOutTSS = pcrm.TimeoutputTimeseries("LakeOutTSS", self, self.LakeID, noHeader=True) - self.LakeStorTSS = pcrm.TimeoutputTimeseries("LakeStorTSS", self, self.LakeID, noHeader=True) - if self.RainRA_FLAG==1: #-set reporting of water balances for individual components - self.LakeRainInTSS = pcrm.TimeoutputTimeseries("LakeRainInTSS", self, self.LakeID, noHeader=True) - self.LakeRainOutTSS = pcrm.TimeoutputTimeseries("LakeRainOutTSS", self, self.LakeID, noHeader=True) - self.LakeRainStorTSS = pcrm.TimeoutputTimeseries("LakeRainStorTSS", self, self.LakeID, noHeader=True) - if self.SnowRA_FLAG==1: - self.LakeSnowInTSS = pcrm.TimeoutputTimeseries("LakeSnowInTSS", self, self.LakeID, noHeader=True) - self.LakeSnowOutTSS = pcrm.TimeoutputTimeseries("LakeSnowOutTSS", self, self.LakeID, noHeader=True) - self.LakeSnowStorTSS = pcrm.TimeoutputTimeseries("LakeSnowStorTSS", self, self.LakeID, noHeader=True) - if self.GlacRA_FLAG==1: - self.LakeGlacInTSS = pcrm.TimeoutputTimeseries("LakeGlacInTSS", self, self.LakeID, noHeader=True) - self.LakeGlacOutTSS = pcrm.TimeoutputTimeseries("LakeGlacOutTSS", self, self.LakeID, noHeader=True) - self.LakeGlacStorTSS = pcrm.TimeoutputTimeseries("LakeGlacStorTSS", self, self.LakeID, noHeader=True) - if self.BaseRA_FLAG==1: - self.LakeBaseInTSS = pcrm.TimeoutputTimeseries("LakeBaseInTSS", self, self.LakeID, noHeader=True) - self.LakeBaseOutTSS = pcrm.TimeoutputTimeseries("LakeBaseOutTSS", self, self.LakeID, noHeader=True) - self.LakeBaseStorTSS = pcrm.TimeoutputTimeseries("LakeBaseStorTSS", self, self.LakeID, noHeader=True) - #-set reporting of water balances for reservoirs - if self.ResFLAG == 1 and config.getint('REPORTING', 'Res_wbal') == 1: - self.ResInTSS = pcrm.TimeoutputTimeseries("ResInTSS", self, self.ResID, noHeader=True) - self.ResOutTSS = pcrm.TimeoutputTimeseries("ResOutTSS", self, self.ResID, noHeader=True) - self.ResStorTSS = pcrm.TimeoutputTimeseries("ResStorTSS", self, self.ResID, noHeader=True) - if self.RainRA_FLAG==1: #-set reporting of water balances for individual components - self.ResRainInTSS = pcrm.TimeoutputTimeseries("ResRainInTSS", self, self.ResID, noHeader=True) - self.ResRainOutTSS = pcrm.TimeoutputTimeseries("ResRainOutTSS", self, self.ResID, noHeader=True) - self.ResRainStorTSS = pcrm.TimeoutputTimeseries("ResRainStorTSS", self, self.ResID, noHeader=True) - if self.SnowRA_FLAG==1: - self.ResSnowInTSS = pcrm.TimeoutputTimeseries("ResSnowInTSS", self, self.ResID, noHeader=True) - self.ResSnowOutTSS = pcrm.TimeoutputTimeseries("ResSnowOutTSS", self, self.ResID, noHeader=True) - self.ResSnowStorTSS = pcrm.TimeoutputTimeseries("ResSnowStorTSS", self, self.ResID, noHeader=True) - if self.GlacRA_FLAG==1: - self.ResGlacInTSS = pcrm.TimeoutputTimeseries("ResGlacInTSS", self, self.ResID, noHeader=True) - self.ResGlacOutTSS = pcrm.TimeoutputTimeseries("ResGlacOutTSS", self, self.ResID, noHeader=True) - self.ResGlacStorTSS = pcrm.TimeoutputTimeseries("ResGlacStorTSS", self, self.ResID, noHeader=True) - if self.BaseRA_FLAG==1: - self.ResBaseInTSS = pcrm.TimeoutputTimeseries("ResBaseInTSS", self, self.ResID, noHeader=True) - self.ResBaseOutTSS = pcrm.TimeoutputTimeseries("ResBaseOutTSS", self, self.ResID, noHeader=True) - self.ResBaseStorTSS = pcrm.TimeoutputTimeseries("ResBaseStorTSS", self, self.ResID, noHeader=True) - - - - - - - - - # if self.scalarInput: - # self.gaugesMap=self.wf_readmap(os.path.join(self.Dir , wflow_mgauges),0.0,fail=True) #: Map with locations of rainfall/evap/temp gauge(s). Only needed if the input to the model is not in maps - # self.OutputId=self.wf_readmap(os.path.join(self.Dir , wflow_subcatch),0.0,fail=True) # location of subcatchment - - self.ZeroMap=0.0*scalar(defined(self.DEM)) #map with only zero's - - - # For in memory override: - #self.Prec, self.Tair, self.Tmax, self.Tmin = self.ZeroMap - - - # Set static initial values here ######################################### - self.Latitude = ycoordinate(boolean(self.ZeroMap)) - self.Longitude = xcoordinate(boolean(self.ZeroMap)) - - # self.logger.info("Linking parameters to landuse, catchment and soil...") - - # self.Beta = scalar(0.6) # For sheetflow - # #self.M=lookupscalar(self.Dir + "/" + modelEnv['intbl'] + "/M.tbl" ,self.LandUse,subcatch,self.Soil) # Decay parameter in Topog_sbm - # self.N=lookupscalar(self.Dir + "/" + self.intbl + "/N.tbl",self.LandUse,subcatch,self.Soil) # Manning overland flow - # """ *Parameter:* Manning's N for all non-river cells """ - # self.NRiver=lookupscalar(self.Dir + "/" + self.intbl + "/N_River.tbl",self.LandUse,subcatch,self.Soil) # Manning river - # """ Manning's N for all cells that are marked as a river """ - - self.wf_updateparameters() - - - # Multiply parameters with a factor (for calibration etc) -P option in command line - - self.wf_multparameters() + # self.RootWater = config.getfloat('SOIL_INIT','RootWater') + # except: + # self.RootWater = readmap(self.inpath + config.get('SOIL_INIT','RootWater')) + # #-initial water content in subsoil + # if not config.get('SOIL_INIT','SubWater'): + # self.SubWater = self.SubField + # else: + # try: + # self.SubWater = config.getfloat('SOIL_INIT','SubWater') + # except: + # self.SubWater = readmap(self.inpath + config.get('SOIL_INIT','SubWater')) + # -initial water storage in rootzone + subsoil + self.SoilWater = self.RootWater + self.SubWater + # -initial capillary rise + self.CapRise = configget(self.config, "SOIL_INIT", "CapRise", 3) + # try: + # self.CapRise = config.getfloat('SOIL_INIT','CapRise') + # except: + # self.CapRise = readmap(self.inpath + config.get('SOIL_INIT','CapRise')) + # -initial drainage from rootzone + self.RootDrain = configget(self.config, "SOIL_INIT", "RootDrain", 3) + # try: + # self.RootDrain = config.getfloat('SOIL_INIT','RootDrain') + # except: + # self.RootDrain = readmap(self.inpath + config.get('SOIL_INIT','RootDrain')) + if self.DynVegFLAG == 1: + # -initial canopy storage + self.Scanopy = 0 + # -initial ndvi if first map is not provided + self.ndviOld = scalar((self.NDVImax + self.NDVImin) / 2) + elif self.KcStatFLAG == 0: + # -set initial kc value to one, if kc map is not available for first timestep + self.KcOld = scalar(1) + # -initial groundwater properties + if self.GroundFLAG == 1: + self.GwRecharge = float( + configget(self.config, "GROUNDW_INIT", "GwRecharge", 0) + ) + self.BaseR = float(configget(self.config, "GROUNDW_INIT", "BaseR", 1)) + self.Gw = float(configget(self.config, "GROUNDW_INIT", "Gw", 1)) + self.H_gw = float(configget(self.config, "GROUNDW_INIT", "H_gw", 1)) + # #-initial groundwater recharge + # try: + # self.GwRecharge = config.getfloat('GROUNDW_INIT','GwRecharge') + # except: + # self.GwRecharge = readmap(self.inpath + config.get('GROUNDW_INIT','GwRecharge')) + # #-initial baseflow + # try: + # self.BaseR = config.getfloat('GROUNDW_INIT','BaseR') + # except: + # self.BaseR = readmap(self.inpath + config.get('GROUNDW_INIT','BaseR')) + # #-initial groundwater storage + # try: + # self.Gw = config.getfloat('GROUNDW_INIT','Gw') + # except: + # self.Gw = readmap(self.inpath + config.get('GROUNDW_INIT','Gw')) + # #-initial groundwater level + # try: + # self.H_gw = config.getfloat('GROUNDW_INIT','H_gw') + # except: + # self.H_gw = readmap(self.inpath + config.get('GROUNDW_INIT','H_gw')) + # self.H_gw = max((self.RootDepthFlat + self.SubDepthFlat + self.GwDepth)/1000 - self.H_gw, 0) + # else: + # #-initial drainage from subsoil + # try: + # self.SubDrain = config.getfloat('SOIL_INIT','SubDrain') + # except: + # self.SubDrain = readmap(self.inpath + config.get('SOIL_INIT','SubDrain')) + # #-initial seepage value if seepage map series is used + # if self.SeepStatFLAG == 0: + # self.SeepOld = scalar(0) + # -initial snow properties + if self.SnowFLAG == 1: + try: + # self.SnowStore = config.getfloat('SNOW_INIT','SnowIni') + self.SnowStore = float( + configget(self.config, "SNOW_INIT", "SnowIni", 0) + ) + except: + self.SnowStore = readmap( + self.inpath + config.get("SNOW_INIT", "SnowIni") + ) + # -initial water stored in snowpack + try: + self.SnowWatStore = float( + configget(self.config, "SNOW_INIT", "SnowWatStore", 0) + ) + # self.SnowWatStore = config.getfloat('SNOW_INIT','SnowWatStore') + except: + self.SnowWatStore = readmap( + self.inpath + config.get("SNOW_INIT", "SnowWatStore") + ) + self.TotalSnowStore = self.SnowStore + self.SnowWatStore - def default_summarymaps(self): ##-maybe not needed. check later - """ + # -initial glacier properties + if self.GlacFLAG == 1: + # try: + # self.GlacFrac = config.getfloat('GLACIER_INIT','GlacFrac') + # except: + # self.GlacFrac = readmap(self.inpath + config.get('GLACIER_INIT','GlacFrac')) + self.GlacFrac = readmap( + os.path.join( + self.Dir, + "staticmaps", + configget( + self.config, "GLACIER_INIT", "GlacFrac", "glacierfraction.map" + ), + ) + ) + print self.GlacFrac + # -initial routed total runoff and of individual components + if self.RoutFLAG == 1 or self.LakeFLAG == 1 or self.ResFLAG == 1: + # -initial routed total runoff + try: + self.QRAold = config.getfloat("ROUT_INIT", "QRA_init") + except: + try: + self.QRAold = readmap( + self.inpath + config.get("ROUT_INIT", "QRA_init") + ) + except: + self.QRAold = 0 + # -initial routed runoff for the individual components + pars = ["RainRA", "SnowRA", "GlacRA", "BaseRA"] + self.RainRAold = 0 + self.SnowRAold = 0 + self.GlacRAold = 0 + self.BaseRAold = 0 + self.RainRA_FLAG = True + self.SnowRA_FLAG = True + self.GlacRA_FLAG = True + self.BaseRA_FLAG = True + # for i in pars: + # try: + # setattr(self, i + 'old', readmap(self.inpath + config.get('ROUT_INIT', i + '_init'))) + # setattr(self, i + '_FLAG', True) + # except: + # try: + # #setattr(self, i + 'old', config.getfloat('ROUT_INIT', i + '_init')) + # setattr(self, i + '_FLAG', True) + # print RainRA_init + # except: + # setattr(self, i + '_FLAG', False) + + # -initial storage in lakes and reservoirs + if self.LakeFLAG == 1 or self.ResFLAG == 1: + # -Read initial storages from table/reservoir file + if self.LakeFLAG == 1: + LakeStor_Tab = self.inpath + config.get("LAKE", "LakeStor") + self.StorRES = ( + cover(lookupscalar(LakeStor_Tab, 1, self.LakeID), 0) * 10 ** 6 + ) # convert to m3 + # -Qfrac for lake cells should be zero, else 1 + self.QFRAC = ifthenelse(self.LakeID != 0, scalar(0), 1) + if self.ResFLAG == 1: + ResStor_Tab = self.inpath + config.get("RESERVOIR", "ResFuncStor") + ResStor = ( + cover(lookupscalar(ResStor_Tab, 2, self.ResID), 0) * 10 ** 6 + ) # convert to m3 + try: + self.StorRES = self.StorRES + ResStor + # -Qfrac for reservoir cells should be zero, else 1 + self.QFRAC = ifthenelse(self.ResID != 0, scalar(0), self.QFRAC) + except: + self.StorRES = ResStor + # -Qfrac for reservoir cells should be zero, else 1 + self.QFRAC = ifthenelse(self.ResID != 0, scalar(0), 1) + + # -initial storage in lakes/reservoirs of individual flow components + pars = ["RainRA", "SnowRA", "GlacRA", "BaseRA"] + for i in pars: + column = pars.index( + i + ) # identify column to be read from lake or reservoir table + try: # -try to sum the storages read from the lake and reservoir tables if both thse modules are used + setattr( + self, + i + "stor", + ( + cover( + lookupscalar(LakeStor_Tab, column + 2, self.LakeID), 0 + ) + + cover( + lookupscalar(ResStor_Tab, column + 3, self.ResID), 0 + ) + ) + * 10 ** 6, + ) + if eval("self." + i + "_FLAG"): + setattr(self, i + "_FLAG", True) + else: + setattr(self, i + "_FLAG", False) + except: + try: # -try to read the storages from the lake table + setattr( + self, + i + "stor", + cover( + lookupscalar(LakeStor_Tab, column + 2, self.LakeID), 0 + ) + * 10 ** 6, + ) + if eval("self." + i + "_FLAG"): + setattr(self, i + "_FLAG", True) + else: + setattr(self, i + "_FLAG", False) + except: # -try to read the storages from the reservoir table + try: + setattr( + self, + i + "stor", + cover( + lookupscalar(ResStor_Tab, column + 3, self.ResID), 0 + ) + * 10 ** 6, + ) + if eval("self." + i + "_FLAG"): + setattr(self, i + "_FLAG", True) + else: + setattr(self, i + "_FLAG", False) + except: + setattr(self, i + "_FLAG", False) + + # -Initial values for reporting and setting of time-series + # -set time-series reporting for mm flux from upstream area for prec and eta + if self.mm_rep_FLAG == 1 and ( + self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1 + ): + self.PrecSubBasinTSS = TimeoutputTimeseries( + "PrecSubBasinTSS", self, self.Locations, noHeader=False + ) + self.ETaSubBasinTSS = TimeoutputTimeseries( + "ETaSubBasinTSS", self, self.Locations, noHeader=False + ) + if self.GlacFLAG == 1: + pars = [ + "wbal", + "GWL", + "TotPrec", + "TotPrecF", + "TotPrecEF", + "TotIntF", + "TotRain", + "TotRainF", + "TotETpot", + "TotETpotF", + "TotETact", + "TotETactF", + "TotSnow", + "TotSnowF", + "TotSnowMelt", + "TotSnowMeltF", + "TotGlacMelt", + "TotGlacMeltF", + "TotRootRF", + "TotRootDF", + "TotRootPF", + "TotSubPF", + "TotCapRF", + "TotGlacPercF", + "TotGwRechargeF", + "TotRainRF", + "TotBaseRF", + "TotSnowRF", + "TotGlacRF", + "TotRF", + "RainRAtot", + "SnowRAtot", + "GlacRAtot", + "BaseRAtot", + "QallRAtot", + ] + # -set time-series reporting for mm fluxes from upstream area if glacier and routing/reservoir modules are used + if self.mm_rep_FLAG == 1 and ( + self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1 + ): + self.GMeltSubBasinTSS = TimeoutputTimeseries( + "GMeltSubBasinTSS", self, self.Locations, noHeader=False + ) + self.QSNOWSubBasinTSS = TimeoutputTimeseries( + "QSNOWSubBasinTSS", self, self.Locations, noHeader=False + ) + self.QRAINSubBasinTSS = TimeoutputTimeseries( + "QRAINSubBasinTSS", self, self.Locations, noHeader=False + ) + self.QGLACSubBasinTSS = TimeoutputTimeseries( + "QGLACSubBasinTSS", self, self.Locations, noHeader=False + ) + self.QBASFSubBasinTSS = TimeoutputTimeseries( + "QBASFSubBasinTSS", self, self.Locations, noHeader=False + ) + self.QTOTSubBasinTSS = TimeoutputTimeseries( + "QTOTSubBasinTSS", self, self.Locations, noHeader=False + ) + elif self.SnowFLAG == 1: + if self.GroundFLAG == 1: + pars = [ + "wbal", + "GWL", + "TotPrec", + "TotPrecF", + "TotPrecEF", + "TotIntF", + "TotRain", + "TotRainF", + "TotETpot", + "TotETpotF", + "TotETact", + "TotETactF", + "TotSnow", + "TotSnowF", + "TotSnowMelt", + "TotSnowMeltF", + "TotRootRF", + "TotRootDF", + "TotRootPF", + "TotSubPF", + "TotCapRF", + "TotGwRechargeF", + "TotRainRF", + "TotBaseRF", + "TotSnowRF", + "TotRF", + "RainRAtot", + "SnowRAtot", + "BaseRAtot", + "QallRAtot", + ] + # -set time-series reporting for mm fluxes from upstream area if snow, groundwater and routing/reservoir modules are used + if self.mm_rep_FLAG == 1 and ( + self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1 + ): + self.QSNOWSubBasinTSS = TimeoutputTimeseries( + "QSNOWSubBasinTSS", self, self.Locations, noHeader=False + ) + self.QRAINSubBasinTSS = TimeoutputTimeseries( + "QRAINSubBasinTSS", self, self.Locations, noHeader=False + ) + self.QBASFSubBasinTSS = TimeoutputTimeseries( + "QBASFSubBasinTSS", self, self.Locations, noHeader=False + ) + self.QTOTSubBasinTSS = TimeoutputTimeseries( + "QTOTSubBasinTSS", self, self.Locations, noHeader=False + ) + else: + pars = [ + "wbal", + "GWL", + "TotPrec", + "TotPrecF", + "TotPrecEF", + "TotIntF", + "TotRain", + "TotRainF", + "TotETpot", + "TotETpotF", + "TotETact", + "TotETactF", + "TotSnow", + "TotSnowF", + "TotSnowMelt", + "TotSnowMeltF", + "TotRootRF", + "TotRootDF", + "TotRootPF", + "TotSubDF", + "TotCapRF", + "TotSeepF", + "TotRainRF", + "TotSnowRF", + "TotRF", + "RainRAtot", + "SnowRAtot", + "BaseRAtot", + "QallRAtot", + ] + # -set time-series reporting for mm fluxes from upstream area if snow and routing/reservoir modules are used + if self.mm_rep_FLAG == 1 and ( + self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1 + ): + self.SeepSubBasinTSS = TimeoutputTimeseries( + "SeepSubBasinTSS", self, self.Locations, noHeader=False + ) + self.QSNOWSubBasinTSS = TimeoutputTimeseries( + "QSNOWSubBasinTSS", self, self.Locations, noHeader=False + ) + self.QRAINSubBasinTSS = TimeoutputTimeseries( + "QRAINSubBasinTSS", self, self.Locations, noHeader=False + ) + self.QBASFSubBasinTSS = TimeoutputTimeseries( + "QBASFSubBasinTSS", self, self.Locations, noHeader=False + ) + self.QTOTSubBasinTSS = TimeoutputTimeseries( + "QTOTSubBasinTSS", self, self.Locations, noHeader=False + ) + else: + if self.GroundFLAG == 1: + pars = [ + "wbal", + "GWL", + "TotPrec", + "TotPrecF", + "TotPrecEF", + "TotIntF", + "TotRain", + "TotRainF", + "TotETpot", + "TotETpotF", + "TotETact", + "TotETactF", + "TotRootRF", + "TotRootDF", + "TotRootPF", + "TotSubPF", + "TotCapRF", + "TotGwRechargeF", + "TotRainRF", + "TotBaseRF", + "TotRF", + "RainRAtot", + "BaseRAtot", + "QallRAtot", + ] + # -set time-series reporting for mm fluxes from upstream area if groundwater and routing/reservoir modules are used + if self.mm_rep_FLAG == 1 and ( + self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1 + ): + self.QRAINSubBasinTSS = TimeoutputTimeseries( + "QRAINSubBasinTSS", self, self.Locations, noHeader=False + ) + self.QBASFSubBasinTSS = TimeoutputTimeseries( + "QBASFSubBasinTSS", self, self.Locations, noHeader=False + ) + self.QTOTSubBasinTSS = TimeoutputTimeseries( + "QTOTSubBasinTSS", self, self.Locations, noHeader=False + ) + else: + pars = [ + "wbal", + "GWL", + "TotPrec", + "TotPrecF", + "TotPrecEF", + "TotIntF", + "TotRain", + "TotRainF", + "TotETpot", + "TotETpotF", + "TotETact", + "TotETactF", + "TotRootRF", + "TotRootDF", + "TotRootPF", + "TotSubDF", + "TotCapRF", + "TotSeepF", + "TotRainRF", + "TotRF", + "RainRAtot", + "BaseRAtot", + "QallRAtot", + ] + # -set time-series reporting for mm fluxes from upstream area if routing/reservoir modules are used + if self.mm_rep_FLAG == 1 and ( + self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1 + ): + self.SeepSubBasinTSS = TimeoutputTimeseries( + "SeepSubBasinTSS", self, self.Locations, noHeader=False + ) + self.QRAINSubBasinTSS = TimeoutputTimeseries( + "QRAINSubBasinTSS", self, self.Locations, noHeader=False + ) + self.QBASFSubBasinTSS = TimeoutputTimeseries( + "QBASFSubBasinTSS", self, self.Locations, noHeader=False + ) + self.QTOTSubBasinTSS = TimeoutputTimeseries( + "QTOTSubBasinTSS", self, self.Locations, noHeader=False + ) + # -remove routing output from reported list of parameters if these modules are not used + if self.RoutFLAG == 0 and self.ResFLAG == 0 and self.LakeFLAG == 0: + rpars = ["RainRAtot", "SnowRAtot", "GlacRAtot", "BaseRAtot", "QallRAtot"] + for i in rpars: + try: + j = pars.index(i) + del pars[j] + except: + pass + # -set reporting options and read initial values + for i in pars: + mapoutops = configget( + self.config, "REPORTING", i + "_mapoutput", i + "_mapoutput" + ) + # mapoutops = config.get('REPORTING', i+'_mapoutput') + TSoutops = configget( + self.config, "REPORTING", i + "_TSoutput", i + "_TSoutput" + ) + # TSoutops = config.get('REPORTING', i+'_TSoutput') + if mapoutops == "NONE" and TSoutops == "NONE": + print i + " will NOT be reported" + else: + print i + " will be reported" + fname = configget(self.config, "REPORTING", i + "_fname", i + "_fname") + # fname = config.get('REPORTING', i+'_fname') + setattr(self, i + "_fname", fname) + # try: + # setattr(self, i, readmap(self.inpath + config.get('INITTOT', i))) + # except: + # try: + # setattr(self, i, config.getfloat('INITTOT', i)) + # except: + # setattr(self, i, 0.) + setattr( + self, i, 0. + ) # use this instead of the commented part above, because it is more logical to always zero as initial condition for reporting + if mapoutops != "NONE": + mapoutops = mapoutops.split(",") + for j in mapoutops: + if j == "D": + setattr(self, i + "_Day", eval("self." + i)) + setattr(self, i + "_Day_map", 1) + elif j == "M": + setattr(self, i + "_Month", eval("self." + i)) + setattr(self, i + "_Month_map", 1) + elif j == "Y": + setattr(self, i + "_Year", eval("self." + i)) + setattr(self, i + "_Year_map", 1) + else: + setattr(self, i + "_Final", eval("self." + i)) + setattr(self, i + "_Final_map", 1) + if TSoutops != "NONE": + TSoutops = TSoutops.split(",") + for j in TSoutops: + if j == "D": + setattr(self, i + "_Day", eval("self." + i)) + setattr( + self, + i + "_DayTS", + eval( + 'TimeoutputTimeseries("' + + fname + + "DTS" + + '", self, self.Locations, noHeader=False)' + ), + ) + elif j == "M": + setattr(self, i + "_Month", eval("self." + i)) + setattr( + self, + i + "_MonthTS", + eval( + 'TimeoutputTimeseries("' + + fname + + "MTS" + + '", self, self.Locations, noHeader=False)' + ), + ) + elif j == "Y": + setattr(self, i + "_Year", eval("self." + i)) + setattr( + self, + i + "_YearTS", + eval( + 'TimeoutputTimeseries("' + + fname + + "YTS" + + '", self, self.Locations, noHeader=False)' + ), + ) + + # -set reporting of water balances for lakes + if self.LakeFLAG == 1 and config.getint("REPORTING", "Lake_wbal") == 1: + self.LakeInTSS = pcrm.TimeoutputTimeseries( + "LakeInTSS", self, self.LakeID, noHeader=True + ) + self.LakeOutTSS = pcrm.TimeoutputTimeseries( + "LakeOutTSS", self, self.LakeID, noHeader=True + ) + self.LakeStorTSS = pcrm.TimeoutputTimeseries( + "LakeStorTSS", self, self.LakeID, noHeader=True + ) + if ( + self.RainRA_FLAG == 1 + ): # -set reporting of water balances for individual components + self.LakeRainInTSS = pcrm.TimeoutputTimeseries( + "LakeRainInTSS", self, self.LakeID, noHeader=True + ) + self.LakeRainOutTSS = pcrm.TimeoutputTimeseries( + "LakeRainOutTSS", self, self.LakeID, noHeader=True + ) + self.LakeRainStorTSS = pcrm.TimeoutputTimeseries( + "LakeRainStorTSS", self, self.LakeID, noHeader=True + ) + if self.SnowRA_FLAG == 1: + self.LakeSnowInTSS = pcrm.TimeoutputTimeseries( + "LakeSnowInTSS", self, self.LakeID, noHeader=True + ) + self.LakeSnowOutTSS = pcrm.TimeoutputTimeseries( + "LakeSnowOutTSS", self, self.LakeID, noHeader=True + ) + self.LakeSnowStorTSS = pcrm.TimeoutputTimeseries( + "LakeSnowStorTSS", self, self.LakeID, noHeader=True + ) + if self.GlacRA_FLAG == 1: + self.LakeGlacInTSS = pcrm.TimeoutputTimeseries( + "LakeGlacInTSS", self, self.LakeID, noHeader=True + ) + self.LakeGlacOutTSS = pcrm.TimeoutputTimeseries( + "LakeGlacOutTSS", self, self.LakeID, noHeader=True + ) + self.LakeGlacStorTSS = pcrm.TimeoutputTimeseries( + "LakeGlacStorTSS", self, self.LakeID, noHeader=True + ) + if self.BaseRA_FLAG == 1: + self.LakeBaseInTSS = pcrm.TimeoutputTimeseries( + "LakeBaseInTSS", self, self.LakeID, noHeader=True + ) + self.LakeBaseOutTSS = pcrm.TimeoutputTimeseries( + "LakeBaseOutTSS", self, self.LakeID, noHeader=True + ) + self.LakeBaseStorTSS = pcrm.TimeoutputTimeseries( + "LakeBaseStorTSS", self, self.LakeID, noHeader=True + ) + # -set reporting of water balances for reservoirs + if self.ResFLAG == 1 and config.getint("REPORTING", "Res_wbal") == 1: + self.ResInTSS = pcrm.TimeoutputTimeseries( + "ResInTSS", self, self.ResID, noHeader=True + ) + self.ResOutTSS = pcrm.TimeoutputTimeseries( + "ResOutTSS", self, self.ResID, noHeader=True + ) + self.ResStorTSS = pcrm.TimeoutputTimeseries( + "ResStorTSS", self, self.ResID, noHeader=True + ) + if ( + self.RainRA_FLAG == 1 + ): # -set reporting of water balances for individual components + self.ResRainInTSS = pcrm.TimeoutputTimeseries( + "ResRainInTSS", self, self.ResID, noHeader=True + ) + self.ResRainOutTSS = pcrm.TimeoutputTimeseries( + "ResRainOutTSS", self, self.ResID, noHeader=True + ) + self.ResRainStorTSS = pcrm.TimeoutputTimeseries( + "ResRainStorTSS", self, self.ResID, noHeader=True + ) + if self.SnowRA_FLAG == 1: + self.ResSnowInTSS = pcrm.TimeoutputTimeseries( + "ResSnowInTSS", self, self.ResID, noHeader=True + ) + self.ResSnowOutTSS = pcrm.TimeoutputTimeseries( + "ResSnowOutTSS", self, self.ResID, noHeader=True + ) + self.ResSnowStorTSS = pcrm.TimeoutputTimeseries( + "ResSnowStorTSS", self, self.ResID, noHeader=True + ) + if self.GlacRA_FLAG == 1: + self.ResGlacInTSS = pcrm.TimeoutputTimeseries( + "ResGlacInTSS", self, self.ResID, noHeader=True + ) + self.ResGlacOutTSS = pcrm.TimeoutputTimeseries( + "ResGlacOutTSS", self, self.ResID, noHeader=True + ) + self.ResGlacStorTSS = pcrm.TimeoutputTimeseries( + "ResGlacStorTSS", self, self.ResID, noHeader=True + ) + if self.BaseRA_FLAG == 1: + self.ResBaseInTSS = pcrm.TimeoutputTimeseries( + "ResBaseInTSS", self, self.ResID, noHeader=True + ) + self.ResBaseOutTSS = pcrm.TimeoutputTimeseries( + "ResBaseOutTSS", self, self.ResID, noHeader=True + ) + self.ResBaseStorTSS = pcrm.TimeoutputTimeseries( + "ResBaseStorTSS", self, self.ResID, noHeader=True + ) + + # if self.scalarInput: + # self.gaugesMap=self.wf_readmap(os.path.join(self.Dir , wflow_mgauges),0.0,fail=True) #: Map with locations of rainfall/evap/temp gauge(s). Only needed if the input to the model is not in maps + # self.OutputId=self.wf_readmap(os.path.join(self.Dir , wflow_subcatch),0.0,fail=True) # location of subcatchment + + self.ZeroMap = 0.0 * scalar(defined(self.DEM)) # map with only zero's + + # For in memory override: + # self.Prec, self.Tair, self.Tmax, self.Tmin = self.ZeroMap + + # Set static initial values here ######################################### + self.Latitude = ycoordinate(boolean(self.ZeroMap)) + self.Longitude = xcoordinate(boolean(self.ZeroMap)) + + # self.logger.info("Linking parameters to landuse, catchment and soil...") + + # self.Beta = scalar(0.6) # For sheetflow + # #self.M=lookupscalar(self.Dir + "/" + modelEnv['intbl'] + "/M.tbl" ,self.LandUse,subcatch,self.Soil) # Decay parameter in Topog_sbm + # self.N=lookupscalar(self.Dir + "/" + self.intbl + "/N.tbl",self.LandUse,subcatch,self.Soil) # Manning overland flow + # """ *Parameter:* Manning's N for all non-river cells """ + # self.NRiver=lookupscalar(self.Dir + "/" + self.intbl + "/N_River.tbl",self.LandUse,subcatch,self.Soil) # Manning river + # """ Manning's N for all cells that are marked as a river """ + + self.wf_updateparameters() + + # Multiply parameters with a factor (for calibration etc) -P option in command line + + self.wf_multparameters() + + def default_summarymaps(self): ##-maybe not needed. check later + """ Returns a list of default summary-maps at the end of a run. This is model specific. You can also add them to the [summary]section of the ini file but stuff you think is crucial to the model should be listed here Example: """ - #lst = ['self.Cfmax','self.csize','self.upsize','self.TTI','self.TT','self.WHC', - # 'self.Slope','self.N','self.xl','self.yl','self.reallength','self.DCL','self.Bw',] - lst = ['self.GlacFrac'] + # lst = ['self.Cfmax','self.csize','self.upsize','self.TTI','self.TT','self.WHC', + # 'self.Slope','self.N','self.xl','self.yl','self.reallength','self.DCL','self.Bw',] + lst = ["self.GlacFrac"] - return lst + return lst - def resume(self): - """ read initial state maps (they are output of a previous call to suspend()) """ + def resume(self): + """ read initial state maps (they are output of a previous call to suspend()) """ - if self.reinit == 1: #-to be defined for sphy model state variables!!! - self.logger.info("Setting initial conditions to default (zero!)") - self.RootWater = RootFieldMap - self.SubWater = SubFieldMap - self.CapRise = 3 - self.RootDrain = 3 - self.SubDrain = 3 - self.GwRecharge = 2 - self.BaseR = 1 - self.Gw = 1500 - self.H_gw = 3 - self.SnowIni = cover(0.0) - self.SnowWatStore = cover(0.0) - self.QRA_init = cover(0.0) - self.RainRA_init = cover(0.0) - self.BaseRA_init = cover(0.0) - self.SnowRA_init = cover(0.0) - self.GlacRA_init = cover(0.0) - #self.FreeWater = cover(0.0) #: Water on surface (state variable [mm]) - #self.SoilMoisture = self.FC #: Soil moisture (state variable [mm]) - #self.UpperZoneStorage = 0.2 * self.FC #: Storage in Upper Zone (state variable [mm]) - #self.LowerZoneStorage = 1.0/(3.0 * self.K4) #: Storage in Uppe Zone (state variable [mm]) - #self.InterceptionStorage = cover(0.0) #: Interception Storage (state variable [mm]) - #self.SurfaceRunoff = cover(0.0) #: Discharge in kinimatic wave (state variable [m^3/s]) - #self.WaterLevel = cover(0.0) #: Water level in kinimatic wave (state variable [m]) - #self.DrySnow=cover(0.0) #: Snow amount (state variable [mm]) - #if hasattr(self, 'ReserVoirSimpleLocs'): - # self.ReservoirVolume = self.ResMaxVolume * self.ResTargetFullFrac - #if hasattr(self, 'ReserVoirComplexLocs'): - # self.ReservoirWaterLevel = cover(0.0) - else: - self.wf_resume(os.path.join(self.Dir, "instate")) - + if self.reinit == 1: # -to be defined for sphy model state variables!!! + self.logger.info("Setting initial conditions to default (zero!)") + self.RootWater = RootFieldMap + self.SubWater = SubFieldMap + self.CapRise = 3 + self.RootDrain = 3 + self.SubDrain = 3 + self.GwRecharge = 2 + self.BaseR = 1 + self.Gw = 1500 + self.H_gw = 3 + self.SnowIni = cover(0.0) + self.SnowWatStore = cover(0.0) + self.QRA_init = cover(0.0) + self.RainRA_init = cover(0.0) + self.BaseRA_init = cover(0.0) + self.SnowRA_init = cover(0.0) + self.GlacRA_init = cover(0.0) + # self.FreeWater = cover(0.0) #: Water on surface (state variable [mm]) + # self.SoilMoisture = self.FC #: Soil moisture (state variable [mm]) + # self.UpperZoneStorage = 0.2 * self.FC #: Storage in Upper Zone (state variable [mm]) + # self.LowerZoneStorage = 1.0/(3.0 * self.K4) #: Storage in Uppe Zone (state variable [mm]) + # self.InterceptionStorage = cover(0.0) #: Interception Storage (state variable [mm]) + # self.SurfaceRunoff = cover(0.0) #: Discharge in kinimatic wave (state variable [m^3/s]) + # self.WaterLevel = cover(0.0) #: Water level in kinimatic wave (state variable [m]) + # self.DrySnow=cover(0.0) #: Snow amount (state variable [mm]) + # if hasattr(self, 'ReserVoirSimpleLocs'): + # self.ReservoirVolume = self.ResMaxVolume * self.ResTargetFullFrac + # if hasattr(self, 'ReserVoirComplexLocs'): + # self.ReservoirWaterLevel = cover(0.0) + else: + self.wf_resume(os.path.join(self.Dir, "instate")) + def dynamic(self): - def dynamic(self): - - self.wf_updateparameters() # read forcing an dynamic parameters - #self.counter+=1 - - #print str(self.curdate.day)+'-'+str(self.curdate.month)+'-'+str(self.curdate.year)+' t = '+str(self.counter) + self.wf_updateparameters() # read forcing an dynamic parameters + # self.counter+=1 - # Snow and glacier fraction settings - if self.GlacFLAG == 0: - self.GlacFrac = 0 - if self.SnowFLAG == 0: - self.SnowStore = scalar(0) - SnowFrac = ifthenelse(self.SnowStore > 0, scalar(1 - self.GlacFrac), 0) - RainFrac = ifthenelse(self.SnowStore == 0, scalar(1 - self.GlacFrac), 0) - - #-Read the precipitation time-series - Precip = self.Prec - #-Report Precip - self.reporting.reporting(self, pcr, 'TotPrec', Precip) - self.reporting.reporting(self, pcr, 'TotPrecF', Precip * (1-self.GlacFrac)) + # print str(self.curdate.day)+'-'+str(self.curdate.month)+'-'+str(self.curdate.year)+' t = '+str(self.counter) - #-Temperature and reference evapotranspiration - Temp = self.Tair - if self.ETREF_FLAG == 0: - TempMax = self.Tmax - TempMin = self.Tmin - ETref = self.Hargreaves.Hargreaves(pcr, self.Hargreaves.extrarad(self, pcr), Temp, TempMax, TempMin) - else: - ETref = readmap(generateNameT(self.ETref, self.counter)) + # Snow and glacier fraction settings + if self.GlacFLAG == 0: + self.GlacFrac = 0 + if self.SnowFLAG == 0: + self.SnowStore = scalar(0) + SnowFrac = ifthenelse(self.SnowStore > 0, scalar(1 - self.GlacFrac), 0) + RainFrac = ifthenelse(self.SnowStore == 0, scalar(1 - self.GlacFrac), 0) - #-Interception and effective precipitation - #-Update canopy storage - if self.DynVegFLAG == 1: - #-try to read the ndvi map series. If not available, then use ndvi old - try: - ndvi = readmap(pcrm.generateNameT(self.ndvi, self.counter)) - self.ndviOld = ndvi - except: - ndvi = self.ndviOld - #-fill missing ndvi values with ndvi base - ndvi = ifthenelse(defined(ndvi) == 1, ndvi, self.NDVIbase) - #-calculate the vegetation parameters - vegoutput = self.dynamic_veg.Veg_function(pcr, ndvi, self.FPARmax, self.FPARmin, self.LAImax, self.NDVImin, self.NDVImax, self.KCmin, self.KCmax) - #-Kc - self.Kc = vegoutput[0] - #-Update canopy storage - self.Scanopy = self.Scanopy + Precip - #-interception and effective precipitation - intercep = self.dynamic_veg.Inter_function(pcr, self.Scanopy, vegoutput[1], ETref) - #-interception - Int = intercep[0] - #-report interception corrected for fraction - self.reporting.reporting(self, pcr, 'TotIntF', Int * (1-self.GlacFrac)) - #-effective precipitation - Precip = intercep[1] - #-Report effective precipitation corrected for fraction - self.reporting.reporting(self, pcr, 'TotPrecEF', Precip * (1-self.GlacFrac)) - #-canopy storage - self.Scanopy = intercep[2] - elif self.KcStatFLAG == 0: - #-Try to read the KC map series - try: - self.Kc = readmap(pcrm.generateNameT(self.Kcmaps, self.counter)) - self.KcOld = self.Kc - except: - self.Kc = self.KcOld - #-report mm effective precipitation for sub-basin averages - if self.mm_rep_FLAG == 1 and (self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1): - self.PrecSubBasinTSS.sample(catchmenttotal(Precip * (1-self.GlacFrac), self.FlowDir) / catchmenttotal(1, self.FlowDir)) + # -Read the precipitation time-series + Precip = self.Prec + # -Report Precip + self.reporting.reporting(self, pcr, "TotPrec", Precip) + self.reporting.reporting(self, pcr, "TotPrecF", Precip * (1 - self.GlacFrac)) - # Snow and rain - if self.SnowFLAG == 1: - #-Snow and rain differentiation - Snow = ifthenelse(Temp >= self.Tcrit, 0, Precip) - Rain = ifthenelse(Temp < self.Tcrit, 0, Precip) - #-Report Snow - self.reporting.reporting(self, pcr, 'TotSnow', Snow) - self.reporting.reporting(self, pcr, 'TotSnowF', Snow * (1-self.GlacFrac)) - #-Snow melt - PotSnowMelt = self.snow.PotSnowMelt(pcr, Temp, self.DDFS) - ActSnowMelt = self.snow.ActSnowMelt(pcr, self.SnowStore, PotSnowMelt) - #-Report snow melt - self.reporting.reporting(self, pcr, 'TotSnowMelt', ActSnowMelt) - self.reporting.reporting(self, pcr, 'TotSnowMeltF', ActSnowMelt * SnowFrac) - #-Update snow store - self.SnowStore = self.snow.SnowStoreUpdate(pcr, self.SnowStore, Snow, ActSnowMelt, Temp, self.SnowWatStore) - #-Caclulate the maximum amount of water that can be stored in snowwatstore - MaxSnowWatStore = self.snow.MaxSnowWatStorage(self.SnowSC, self.SnowStore) - OldSnowWatStore = self.SnowWatStore - #-Calculate the actual amount of water stored in snowwatstore - self.SnowWatStore = self.snow.SnowWatStorage(pcr, Temp, MaxSnowWatStore, self.SnowWatStore, ActSnowMelt, Rain) - #-Changes in total water storage in snow (SnowStore and SnowWatStore) - OldTotalSnowStore = self.TotalSnowStore - self.TotalSnowStore = self.snow.TotSnowStorage(self.SnowStore, self.SnowWatStore, SnowFrac, RainFrac) - #-Snow runoff - SnowR = self.snow.SnowR(pcr, self.SnowWatStore, MaxSnowWatStore, ActSnowMelt, Rain, OldSnowWatStore, SnowFrac) - #-Report Snow runoff - self.reporting.reporting(self, pcr, 'TotSnowRF', SnowR) - else: - Rain = Precip - SnowR = 0 - OldTotalSnowStore = 0 - self.TotalSnowStore = 0 - #-Report Rain - self.reporting.reporting(self, pcr, 'TotRain', Rain) - self.reporting.reporting(self, pcr, 'TotRainF', Rain * (1-self.GlacFrac)) + # -Temperature and reference evapotranspiration + Temp = self.Tair + if self.ETREF_FLAG == 0: + TempMax = self.Tmax + TempMin = self.Tmin + ETref = self.Hargreaves.Hargreaves( + pcr, self.Hargreaves.extrarad(self, pcr), Temp, TempMax, TempMin + ) + else: + ETref = readmap(generateNameT(self.ETref, self.counter)) - #-Glacier calculations - if self.GlacFLAG == 1: - #-Glacier melt from clean ice glaciers - GlacCIMelt = self.glacier.GlacCDMelt(pcr, Temp, self.DDFG, self.GlacFracCI) - #-Glacier melt from debris covered glaciers - GlacDCMelt = self.glacier.GlacCDMelt(pcr, Temp, self.DDFDG, self.GlacFracDB) - #-Total melt from glaciers - GlacMelt = self.glacier.GMelt(GlacCIMelt, GlacDCMelt) - self.GlacMelt = GlacMelt - #-Report glacier melt - self.reporting.reporting(self, pcr, 'TotGlacMelt', GlacMelt) - self.reporting.reporting(self, pcr, 'TotGlacMeltF', GlacMelt * self.GlacFrac) - if self.mm_rep_FLAG == 1 and (self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1): - self.GMeltSubBasinTSS.sample(catchmenttotal(GlacMelt * self.GlacFrac, self.FlowDir) / catchmenttotal(1, self.FlowDir)) - #-Glacier runoff - GlacR = self.glacier.GlacR(self.GlacF, GlacMelt, self.GlacFrac) - #-Report glacier runoff - self.reporting.reporting(self, pcr, 'TotGlacRF', GlacR) - #-Glacier percolation to groundwater - GlacPerc = self.glacier.GPerc(self.GlacF, GlacMelt, self.GlacFrac) - #-Report glacier percolation to groundwater - self.reporting.reporting(self, pcr, 'TotGlacPercF', GlacPerc) - else: - GlacR = 0 - GlacMelt = 0 - GlacPerc = 0 + # -Interception and effective precipitation + # -Update canopy storage + if self.DynVegFLAG == 1: + # -try to read the ndvi map series. If not available, then use ndvi old + try: + ndvi = readmap(pcrm.generateNameT(self.ndvi, self.counter)) + self.ndviOld = ndvi + except: + ndvi = self.ndviOld + # -fill missing ndvi values with ndvi base + ndvi = ifthenelse(defined(ndvi) == 1, ndvi, self.NDVIbase) + # -calculate the vegetation parameters + vegoutput = self.dynamic_veg.Veg_function( + pcr, + ndvi, + self.FPARmax, + self.FPARmin, + self.LAImax, + self.NDVImin, + self.NDVImax, + self.KCmin, + self.KCmax, + ) + # -Kc + self.Kc = vegoutput[0] + # -Update canopy storage + self.Scanopy = self.Scanopy + Precip + # -interception and effective precipitation + intercep = self.dynamic_veg.Inter_function( + pcr, self.Scanopy, vegoutput[1], ETref + ) + # -interception + Int = intercep[0] + # -report interception corrected for fraction + self.reporting.reporting(self, pcr, "TotIntF", Int * (1 - self.GlacFrac)) + # -effective precipitation + Precip = intercep[1] + # -Report effective precipitation corrected for fraction + self.reporting.reporting( + self, pcr, "TotPrecEF", Precip * (1 - self.GlacFrac) + ) + # -canopy storage + self.Scanopy = intercep[2] + elif self.KcStatFLAG == 0: + # -Try to read the KC map series + try: + self.Kc = readmap(pcrm.generateNameT(self.Kcmaps, self.counter)) + self.KcOld = self.Kc + except: + self.Kc = self.KcOld + # -report mm effective precipitation for sub-basin averages + if self.mm_rep_FLAG == 1 and ( + self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1 + ): + self.PrecSubBasinTSS.sample( + catchmenttotal(Precip * (1 - self.GlacFrac), self.FlowDir) + / catchmenttotal(1, self.FlowDir) + ) - #-Potential evapotranspiration (THIS SHOULD STILL BE IMPROVED WITH DYNAMIC VEGETATION MODULE) - ETpot = self.ET.ETpot(ETref, self.Kc) - #-Report ETpot - self.reporting.reporting(self, pcr, 'TotETpot', ETpot) - self.reporting.reporting(self, pcr, 'TotETpotF', ETpot * RainFrac) - - #-Rootzone calculations - self.RootWater = self.RootWater + ifthenelse(RainFrac > 0, Rain, 0) + self.CapRise - #-Rootzone runoff - RootRunoff = self.rootzone.RootRunoff(pcr, RainFrac, self.RootWater, self.RootSat) - self.RootWater = self.RootWater - RootRunoff - #-Actual evapotranspiration - etreddry = max(min((self.RootWater - self.RootDry) / (self.RootWilt - self.RootDry), 1), 0) - ETact = self.ET.ETact(pcr, ETpot, self.RootWater, self.RootSat, etreddry, RainFrac) - #-Report the actual evapotranspiration - self.reporting.reporting(self, pcr, 'TotETact', ETact) - #-Actual evapotranspiration, corrected for rain fraction - ActETact = ETact * RainFrac - #-Report the actual evapotranspiration, corrected for rain fraction - self.reporting.reporting(self, pcr, 'TotETactF', ActETact) - if self.mm_rep_FLAG == 1 and (self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1): - self.ETaSubBasinTSS.sample(catchmenttotal(ActETact, self.FlowDir) / catchmenttotal(1, self.FlowDir)) - #-Update rootwater content - self.RootWater = max(self.RootWater - ETact, 0) - #-Rootwater drainage - self.RootDrain = self.rootzone.RootDrainage(pcr, self.RootWater, self.RootDrain, self.RootField, self.RootSat, self.RootDrainVel, self.RootTT) - #-Update rootwater content - self.RootWater = self.RootWater - self.RootDrain - #-Rootwater percolation - rootperc = self.rootzone.RootPercolation(pcr, self.RootWater, self.SubWater, self.RootField, self.RootTT, self.SubSat) - #-Report rootzone percolation, corrected for fraction - self.reporting.reporting(self, pcr, 'TotRootPF', rootperc * (1 - self.GlacFrac)) - #-Update rootwater content - self.RootWater = self.RootWater - rootperc + # Snow and rain + if self.SnowFLAG == 1: + # -Snow and rain differentiation + Snow = ifthenelse(Temp >= self.Tcrit, 0, Precip) + Rain = ifthenelse(Temp < self.Tcrit, 0, Precip) + # -Report Snow + self.reporting.reporting(self, pcr, "TotSnow", Snow) + self.reporting.reporting(self, pcr, "TotSnowF", Snow * (1 - self.GlacFrac)) + # -Snow melt + PotSnowMelt = self.snow.PotSnowMelt(pcr, Temp, self.DDFS) + ActSnowMelt = self.snow.ActSnowMelt(pcr, self.SnowStore, PotSnowMelt) + # -Report snow melt + self.reporting.reporting(self, pcr, "TotSnowMelt", ActSnowMelt) + self.reporting.reporting(self, pcr, "TotSnowMeltF", ActSnowMelt * SnowFrac) + # -Update snow store + self.SnowStore = self.snow.SnowStoreUpdate( + pcr, self.SnowStore, Snow, ActSnowMelt, Temp, self.SnowWatStore + ) + # -Caclulate the maximum amount of water that can be stored in snowwatstore + MaxSnowWatStore = self.snow.MaxSnowWatStorage(self.SnowSC, self.SnowStore) + OldSnowWatStore = self.SnowWatStore + # -Calculate the actual amount of water stored in snowwatstore + self.SnowWatStore = self.snow.SnowWatStorage( + pcr, Temp, MaxSnowWatStore, self.SnowWatStore, ActSnowMelt, Rain + ) + # -Changes in total water storage in snow (SnowStore and SnowWatStore) + OldTotalSnowStore = self.TotalSnowStore + self.TotalSnowStore = self.snow.TotSnowStorage( + self.SnowStore, self.SnowWatStore, SnowFrac, RainFrac + ) + # -Snow runoff + SnowR = self.snow.SnowR( + pcr, + self.SnowWatStore, + MaxSnowWatStore, + ActSnowMelt, + Rain, + OldSnowWatStore, + SnowFrac, + ) + # -Report Snow runoff + self.reporting.reporting(self, pcr, "TotSnowRF", SnowR) + else: + Rain = Precip + SnowR = 0 + OldTotalSnowStore = 0 + self.TotalSnowStore = 0 + # -Report Rain + self.reporting.reporting(self, pcr, "TotRain", Rain) + self.reporting.reporting(self, pcr, "TotRainF", Rain * (1 - self.GlacFrac)) - #-Sub soil calculations - self.SubWater = self.SubWater + rootperc - if self.GroundFLAG == 0: - if self.SeepStatFLAG == 0: - try: - self.SeePage = readmap(pcrm.generateNameT(self.Seepmaps, self.counter)) - self.SeepOld = self.SeePage - except: - self.SeePage = self.SeepOld - #-Report seepage - self.reporting.reporting(self, pcr, 'TotSeepF', scalar(self.SeePage)) - self.SubWater = min(max(self.SubWater - self.SeePage, 0), self.SubSat) - if self.mm_rep_FLAG == 1 and (self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1): - self.SeepSubBasinTSS.sample(catchmenttotal(self.SeePage, self.FlowDir) / catchmenttotal(1, self.FlowDir)) - #-Capillary rise - self.CapRise = self.subzone.CapilRise(pcr, self.SubField, self.SubWater, self.CapRiseMax, self.RootWater, self.RootSat, self.RootField) - #-Report capillary rise, corrected for fraction - self.reporting.reporting(self, pcr, 'TotCapRF', self.CapRise * (1-self.GlacFrac)) - #-Update sub soil water content - self.SubWater = self.SubWater - self.CapRise - if self.GroundFLAG == 1: # sub percolation will be calculated instead of subdrainage - subperc = self.subzone.SubPercolation(pcr, self.SubWater, self.SubField, self.SubTT, self.Gw, self.GwSat) - ActSubPerc = subperc * (1-self.GlacFrac) - #-Report the subzone percolation, corrected for the fraction - self.reporting.reporting(self, pcr, 'TotSubPF', ActSubPerc) - #-Update sub soil water content - self.SubWater = self.SubWater - subperc - else: # sub drainage will be calculated instead of sub percolation - self.SubDrain = self.subzone.SubDrainage(pcr, self.SubWater, self.SubField, self.SubSat, self.SubDrainVel, self.SubDrain, self.SubTT) - #-Report drainage from subzone - self.reporting.reporting(self, pcr, 'TotSubDF', self.SubDrain) - #-Update sub soil water content - self.SubWater = self.SubWater - self.SubDrain - - #-Changes in soil water storage - OldSoilWater = self.SoilWater - self.SoilWater = (self.RootWater + self.SubWater) * (1-self.GlacFrac) + # -Glacier calculations + if self.GlacFLAG == 1: + # -Glacier melt from clean ice glaciers + GlacCIMelt = self.glacier.GlacCDMelt(pcr, Temp, self.DDFG, self.GlacFracCI) + # -Glacier melt from debris covered glaciers + GlacDCMelt = self.glacier.GlacCDMelt(pcr, Temp, self.DDFDG, self.GlacFracDB) + # -Total melt from glaciers + GlacMelt = self.glacier.GMelt(GlacCIMelt, GlacDCMelt) + self.GlacMelt = GlacMelt + # -Report glacier melt + self.reporting.reporting(self, pcr, "TotGlacMelt", GlacMelt) + self.reporting.reporting( + self, pcr, "TotGlacMeltF", GlacMelt * self.GlacFrac + ) + if self.mm_rep_FLAG == 1 and ( + self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1 + ): + self.GMeltSubBasinTSS.sample( + catchmenttotal(GlacMelt * self.GlacFrac, self.FlowDir) + / catchmenttotal(1, self.FlowDir) + ) + # -Glacier runoff + GlacR = self.glacier.GlacR(self.GlacF, GlacMelt, self.GlacFrac) + # -Report glacier runoff + self.reporting.reporting(self, pcr, "TotGlacRF", GlacR) + # -Glacier percolation to groundwater + GlacPerc = self.glacier.GPerc(self.GlacF, GlacMelt, self.GlacFrac) + # -Report glacier percolation to groundwater + self.reporting.reporting(self, pcr, "TotGlacPercF", GlacPerc) + else: + GlacR = 0 + GlacMelt = 0 + GlacPerc = 0 - #-Rootzone runoff - RootR = RootRunoff * RainFrac - #-Report rootzone runoff, corrected for fraction - self.reporting.reporting(self, pcr, 'TotRootRF', RootR) - #-Rootzone drainage - RootD = self.RootDrain * (1-self.GlacFrac) - #-Report rootzone drainage, corrected for fraction - self.reporting.reporting(self, pcr, 'TotRootDF', RootD) - #-Rain runoff - RainR = RootR + RootD - #-Report rain runoff - self.reporting.reporting(self, pcr, 'TotRainRF', RainR) + # -Potential evapotranspiration (THIS SHOULD STILL BE IMPROVED WITH DYNAMIC VEGETATION MODULE) + ETpot = self.ET.ETpot(ETref, self.Kc) + # -Report ETpot + self.reporting.reporting(self, pcr, "TotETpot", ETpot) + self.reporting.reporting(self, pcr, "TotETpotF", ETpot * RainFrac) - #-Groundwater calculations - if self.GroundFLAG == 1: - GwOld = self.Gw - #-Groundwater recharge - self.GwRecharge = self.groundwater.GroundWaterRecharge(pcr, self.deltaGw, self.GwRecharge, ActSubPerc, GlacPerc) - #-Report groundwater recharge - self.reporting.reporting(self, pcr, 'TotGwRechargeF', self.GwRecharge) - #-Update groundwater storage - self.Gw = self.Gw + self.GwRecharge - #-Baseflow - self.BaseR = self.groundwater.BaseFlow(pcr, self.Gw, self.BaseR, self.GwRecharge, self.BaseThresh, self.alphaGw) - #-Report Baseflow - self.reporting.reporting(self, pcr, 'TotBaseRF', self.BaseR) - #-Update groundwater storage - self.Gw = self.Gw - self.BaseR - #-Calculate groundwater level - self.H_gw = self.groundwater.HLevel(pcr, self.H_gw, self.alphaGw, self.GwRecharge, self.YieldGw) - #-Report groundwater - self.reporting.reporting(self, pcr, 'GWL', ((self.SubDepthFlat + self.RootDepthFlat + self.GwDepth)/1000 - self.H_gw)*-1) + # -Rootzone calculations + self.RootWater = ( + self.RootWater + ifthenelse(RainFrac > 0, Rain, 0) + self.CapRise + ) + # -Rootzone runoff + RootRunoff = self.rootzone.RootRunoff( + pcr, RainFrac, self.RootWater, self.RootSat + ) + self.RootWater = self.RootWater - RootRunoff + # -Actual evapotranspiration + etreddry = max( + min((self.RootWater - self.RootDry) / (self.RootWilt - self.RootDry), 1), 0 + ) + ETact = self.ET.ETact( + pcr, ETpot, self.RootWater, self.RootSat, etreddry, RainFrac + ) + # -Report the actual evapotranspiration + self.reporting.reporting(self, pcr, "TotETact", ETact) + # -Actual evapotranspiration, corrected for rain fraction + ActETact = ETact * RainFrac + # -Report the actual evapotranspiration, corrected for rain fraction + self.reporting.reporting(self, pcr, "TotETactF", ActETact) + if self.mm_rep_FLAG == 1 and ( + self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1 + ): + self.ETaSubBasinTSS.sample( + catchmenttotal(ActETact, self.FlowDir) / catchmenttotal(1, self.FlowDir) + ) + # -Update rootwater content + self.RootWater = max(self.RootWater - ETact, 0) + # -Rootwater drainage + self.RootDrain = self.rootzone.RootDrainage( + pcr, + self.RootWater, + self.RootDrain, + self.RootField, + self.RootSat, + self.RootDrainVel, + self.RootTT, + ) + # -Update rootwater content + self.RootWater = self.RootWater - self.RootDrain + # -Rootwater percolation + rootperc = self.rootzone.RootPercolation( + pcr, self.RootWater, self.SubWater, self.RootField, self.RootTT, self.SubSat + ) + # -Report rootzone percolation, corrected for fraction + self.reporting.reporting(self, pcr, "TotRootPF", rootperc * (1 - self.GlacFrac)) + # -Update rootwater content + self.RootWater = self.RootWater - rootperc - else: - #-Use drainage from subsoil as baseflow - self.BaseR = self.SubDrain - #-Groundwater level as scaled between min and max measured gwl - SoilAct = self.RootWater + self.SubWater; - SoilRel = (SoilAct - self.SoilMin) / (self.SoilMax - self.SoilMin) # scale between 0 (dry) and 1 (wet) - GWL = self.GWL_base - (SoilRel-0.5) * self.GWL_base - #-Report groundwater - self.reporting.reporting(self, pcr, 'GWL', GWL) + # -Sub soil calculations + self.SubWater = self.SubWater + rootperc + if self.GroundFLAG == 0: + if self.SeepStatFLAG == 0: + try: + self.SeePage = readmap( + pcrm.generateNameT(self.Seepmaps, self.counter) + ) + self.SeepOld = self.SeePage + except: + self.SeePage = self.SeepOld + # -Report seepage + self.reporting.reporting(self, pcr, "TotSeepF", scalar(self.SeePage)) + self.SubWater = min(max(self.SubWater - self.SeePage, 0), self.SubSat) + if self.mm_rep_FLAG == 1 and ( + self.RoutFLAG == 1 or self.ResFLAG == 1 or self.LakeFLAG == 1 + ): + self.SeepSubBasinTSS.sample( + catchmenttotal(self.SeePage, self.FlowDir) + / catchmenttotal(1, self.FlowDir) + ) + # -Capillary rise + self.CapRise = self.subzone.CapilRise( + pcr, + self.SubField, + self.SubWater, + self.CapRiseMax, + self.RootWater, + self.RootSat, + self.RootField, + ) + # -Report capillary rise, corrected for fraction + self.reporting.reporting( + self, pcr, "TotCapRF", self.CapRise * (1 - self.GlacFrac) + ) + # -Update sub soil water content + self.SubWater = self.SubWater - self.CapRise + if ( + self.GroundFLAG == 1 + ): # sub percolation will be calculated instead of subdrainage + subperc = self.subzone.SubPercolation( + pcr, self.SubWater, self.SubField, self.SubTT, self.Gw, self.GwSat + ) + ActSubPerc = subperc * (1 - self.GlacFrac) + # -Report the subzone percolation, corrected for the fraction + self.reporting.reporting(self, pcr, "TotSubPF", ActSubPerc) + # -Update sub soil water content + self.SubWater = self.SubWater - subperc + else: # sub drainage will be calculated instead of sub percolation + self.SubDrain = self.subzone.SubDrainage( + pcr, + self.SubWater, + self.SubField, + self.SubSat, + self.SubDrainVel, + self.SubDrain, + self.SubTT, + ) + # -Report drainage from subzone + self.reporting.reporting(self, pcr, "TotSubDF", self.SubDrain) + # -Update sub soil water content + self.SubWater = self.SubWater - self.SubDrain - #-Report Total runoff - self.reporting.reporting(self, pcr, 'TotRF', self.BaseR + RainR + SnowR + GlacR) + # -Changes in soil water storage + OldSoilWater = self.SoilWater + self.SoilWater = (self.RootWater + self.SubWater) * (1 - self.GlacFrac) - #-Water balance - if self.GroundFLAG == 1: - waterbalance = Precip * (1-self.GlacFrac) + GlacMelt * self.GlacFrac - ActETact - GlacR - SnowR - RainR -\ - self.BaseR - (self.SoilWater-OldSoilWater) - (self.TotalSnowStore-OldTotalSnowStore) - (self.Gw-GwOld) - elif self.GroundFLAG == 0: - waterbalance = Precip - ActETact - self.SeePage - SnowR - RainR - self.BaseR - (self.SoilWater-OldSoilWater) - (self.TotalSnowStore-OldTotalSnowStore) - self.reporting.reporting(self, pcr, 'wbal', waterbalance) + # -Rootzone runoff + RootR = RootRunoff * RainFrac + # -Report rootzone runoff, corrected for fraction + self.reporting.reporting(self, pcr, "TotRootRF", RootR) + # -Rootzone drainage + RootD = self.RootDrain * (1 - self.GlacFrac) + # -Report rootzone drainage, corrected for fraction + self.reporting.reporting(self, pcr, "TotRootDF", RootD) + # -Rain runoff + RainR = RootR + RootD + # -Report rain runoff + self.reporting.reporting(self, pcr, "TotRainRF", RainR) - #-Routing for lake and/or reservoir modules - if self.LakeFLAG == 1 or self.ResFLAG == 1: - #-Update storage in lakes/reservoirs (m3) with specific runoff - self.StorRES = self.StorRES + ifthenelse(self.QFRAC==0, 0.001 * cellarea() * (self.BaseR + RainR + GlacR + SnowR), 0) - OldStorage = self.StorRES - #-Calculate lake/reservoir outflow volumes - if self.LakeFLAG ==1 and self.ResFLAG ==1: - tempvar = self.lakes.UpdateLakeHStore(self, pcr, pcrm) - LakeLevel = tempvar[0] - self.StorRES = tempvar[1] - LakeQ = self.lakes.QLake(self, pcr, LakeLevel) - ResQ = self.reservoirs.QRes(self, pcr) - Qout = ifthenelse(self.ResID != 0, ResQ, ifthenelse(self.LakeID!=0, LakeQ, 0)) - elif self.LakeFLAG ==1: - tempvar = self.lakes.UpdateLakeHStore(self, pcr, pcrm) - LakeLevel = tempvar[0] - self.StorRES = tempvar[1] - Qout = self.lakes.QLake(self, pcr, LakeLevel) - else: - Qout = self.reservoirs.QRes(self, pcr) + # -Groundwater calculations + if self.GroundFLAG == 1: + GwOld = self.Gw + # -Groundwater recharge + self.GwRecharge = self.groundwater.GroundWaterRecharge( + pcr, self.deltaGw, self.GwRecharge, ActSubPerc, GlacPerc + ) + # -Report groundwater recharge + self.reporting.reporting(self, pcr, "TotGwRechargeF", self.GwRecharge) + # -Update groundwater storage + self.Gw = self.Gw + self.GwRecharge + # -Baseflow + self.BaseR = self.groundwater.BaseFlow( + pcr, self.Gw, self.BaseR, self.GwRecharge, self.BaseThresh, self.alphaGw + ) + # -Report Baseflow + self.reporting.reporting(self, pcr, "TotBaseRF", self.BaseR) + # -Update groundwater storage + self.Gw = self.Gw - self.BaseR + # -Calculate groundwater level + self.H_gw = self.groundwater.HLevel( + pcr, self.H_gw, self.alphaGw, self.GwRecharge, self.YieldGw + ) + # -Report groundwater + self.reporting.reporting( + self, + pcr, + "GWL", + ( + (self.SubDepthFlat + self.RootDepthFlat + self.GwDepth) / 1000 + - self.H_gw + ) + * -1, + ) - #-Calculate volume available for routing (=outflow lakes/reservoir + cell specific runoff) - RunoffVolume = upstream(self.FlowDir, Qout) + ifthenelse(self.QFRAC==0, 0, 0.001 * cellarea() * (self.BaseR + RainR + GlacR + SnowR)) - #-Routing of total flow - tempvar = self.routing.ROUT(self, pcr, RunoffVolume, self.QRAold, Qout, self.StorRES) - self.StorRES = tempvar[0] - Q = tempvar[1] - Qin = tempvar[2] - self.QRAold = Q - self.reporting.reporting(self, pcr, 'QallRAtot', Q) - #-report flux in mm - if self.mm_rep_FLAG == 1: - self.QTOTSubBasinTSS.sample(((Q * 3600 * 24) / catchmenttotal(cellarea(), self.FlowDir)) * 1000) - #-report lake and reservoir waterbalance - if self.LakeFLAG == 1 and config.getint('REPORTING', 'Lake_wbal') ==1: - self.LakeInTSS.sample(Qin) - self.LakeOutTSS.sample(Qout) - self.LakeStorTSS.sample(self.StorRES) - if self.ResFLAG == 1 and config.getint('REPORTING', 'Res_wbal') ==1: - self.ResInTSS.sample(Qin) - self.ResOutTSS.sample(Qout) - self.ResStorTSS.sample(self.StorRES) + else: + # -Use drainage from subsoil as baseflow + self.BaseR = self.SubDrain + # -Groundwater level as scaled between min and max measured gwl + SoilAct = self.RootWater + self.SubWater + SoilRel = (SoilAct - self.SoilMin) / ( + self.SoilMax - self.SoilMin + ) # scale between 0 (dry) and 1 (wet) + GWL = self.GWL_base - (SoilRel - 0.5) * self.GWL_base + # -Report groundwater + self.reporting.reporting(self, pcr, "GWL", GWL) - #-Routing of individual contributers - #-Snow routing - if self.SnowRA_FLAG == 1 and self.SnowFLAG == 1: - self.SnowRAstor = self.SnowRAstor + ifthenelse(self.QFRAC==0, SnowR * 0.001 * cellarea(), 0) - cQfrac = cover(self.SnowRAstor / OldStorage, 0) - cQout = cQfrac * Qout - cRunoffVolume = upstream(self.FlowDir, cQout) + ifthenelse(self.QFRAC==0, 0, 0.001 * cellarea() * SnowR) - tempvar = self.routing.ROUT(self, pcr, cRunoffVolume, self.SnowRAold, cQout, self.SnowRAstor) - self.SnowRAstor = tempvar[0] - SnowRA = tempvar[1] - cQin = tempvar[2] - self.SnowRAold = SnowRA - self.reporting.reporting(self, pcr, 'SnowRAtot', SnowRA) - if self.mm_rep_FLAG == 1: - self.QSNOWSubBasinTSS.sample(((SnowRA * 3600 * 24) / catchmenttotal(cellarea(), self.FlowDir)) *1000) - #-report lake and reservoir waterbalance - if self.LakeFLAG == 1 and config.getint('REPORTING', 'Lake_wbal') ==1: - self.LakeSnowInTSS.sample(cQin) - self.LakeSnowOutTSS.sample(cQout) - self.LakeSnowStorTSS.sample(self.SnowRAstor) - if self.ResFLAG == 1 and config.getint('REPORTING', 'Res_wbal') ==1: - self.ResSnowInTSS.sample(cQin) - self.ResSnowOutTSS.sample(cQout) - self.ResSnowStorTSS.sample(self.SnowRAstor) - #-Rain routing - if self.RainRA_FLAG == 1: - self.RainRAstor = self.RainRAstor + ifthenelse(self.QFRAC==0, RainR * 0.001 * cellarea(), 0) - cQfrac = cover(self.RainRAstor / OldStorage, 0) - cQout = cQfrac * Qout - cRunoffVolume = upstream(self.FlowDir, cQout) + ifthenelse(self.QFRAC==0, 0, 0.001 * cellarea() * RainR) - tempvar = self.routing.ROUT(self, pcr, cRunoffVolume, self.RainRAold, cQout, self.RainRAstor) - self.RainRAstor = tempvar[0] - RainRA = tempvar[1] - cQin = tempvar[2] - self.RainRAold = RainRA - self.reporting.reporting(self, pcr, 'RainRAtot', RainRA) - if self.mm_rep_FLAG == 1: - self.QRAINSubBasinTSS.sample(((RainRA * 3600 * 24) / catchmenttotal(cellarea(), self.FlowDir)) *1000) - #-report lake and reservoir waterbalance - if self.LakeFLAG == 1 and config.getint('REPORTING', 'Lake_wbal') ==1: - self.LakeRainInTSS.sample(cQin) - self.LakeRainOutTSS.sample(cQout) - self.LakeRainStorTSS.sample(self.RainRAstor) - if self.ResFLAG == 1 and config.getint('REPORTING', 'Res_wbal') ==1: - self.ResRainInTSS.sample(cQin) - self.ResRainOutTSS.sample(cQout) - self.ResRainStorTSS.sample(self.RainRAstor) - #-Glacier routing - if self.GlacRA_FLAG == 1 and self.GlacFLAG == 1: - self.GlacRAstor = self.GlacRAstor + ifthenelse(self.QFRAC==0, GlacR * 0.001 * cellarea(), 0) - cQfrac = cover(self.GlacRAstor / OldStorage, 0) - cQout = cQfrac * Qout - cRunoffVolume = upstream(self.FlowDir, cQout) + ifthenelse(self.QFRAC==0, 0, 0.001 * cellarea() * GlacR) - tempvar = self.routing.ROUT(self, pcr, cRunoffVolume, self.GlacRAold, cQout, self.GlacRAstor) - self.GlacRAstor = tempvar[0] - GlacRA = tempvar[1] - cQin = tempvar[2] - self.GlacRAold = GlacRA - self.reporting.reporting(self, pcr, 'GlacRAtot', GlacRA) - if self.mm_rep_FLAG == 1: - self.QGLACSubBasinTSS.sample(((GlacRA * 3600 * 24) / catchmenttotal(cellarea(), self.FlowDir)) *1000) - #-report lake and reservoir waterbalance - if self.LakeFLAG == 1 and config.getint('REPORTING', 'Lake_wbal') ==1: - self.LakeGlacInTSS.sample(cQin) - self.LakeGlacOutTSS.sample(cQout) - self.LakeGlacStorTSS.sample(self.GlacRAstor) - if self.ResFLAG == 1 and config.getint('REPORTING', 'Res_wbal') ==1: - self.ResGlacInTSS.sample(cQin) - self.ResGlacOutTSS.sample(cQout) - self.ResGlacStorTSS.sample(self.GlacRAstor) - #-Baseflow routing - if self.BaseRA_FLAG == 1: - self.BaseRAstor = self.BaseRAstor + ifthenelse(self.QFRAC==0, self.BaseR * 0.001 * cellarea(), 0) - cQfrac = cover(self.BaseRAstor / OldStorage, 0) - cQout = cQfrac * Qout - cRunoffVolume = upstream(self.FlowDir, cQout) + ifthenelse(self.QFRAC==0, 0, 0.001 * cellarea() * self.BaseR) - tempvar = self.routing.ROUT(self, pcr, cRunoffVolume, self.BaseRAold, cQout, self.BaseRAstor) - self.BaseRAstor = tempvar[0] - BaseRA = tempvar[1] - cQin = tempvar[2] - self.BaseRAold = BaseRA - self.reporting.reporting(self, pcr, 'BaseRAtot', BaseRA) - if self.mm_rep_FLAG == 1: - self.QBASFSubBasinTSS.sample(((BaseRA * 3600 * 24) / catchmenttotal(cellarea(), self.FlowDir)) *1000) - #-report lake and reservoir waterbalance - if self.LakeFLAG == 1 and config.getint('REPORTING', 'Lake_wbal') ==1: - self.LakeBaseInTSS.sample(cQin) - self.LakeBaseOutTSS.sample(cQout) - self.LakeBaseStorTSS.sample(self.BaseRAstor) - if self.ResFLAG == 1 and config.getint('REPORTING', 'Res_wbal') ==1: - self.ResBaseInTSS.sample(cQin) - self.ResBaseOutTSS.sample(cQout) - self.ResBaseStorTSS.sample(self.BaseRAstor) - - #-Normal routing module - elif self.RoutFLAG == 1: - #-Rout total runoff - Q = self.routing.ROUT(pcr, self.BaseR + RainR + GlacR + SnowR, self.QRAold, self.FlowDir, self.kx) - self.QRAold = Q - self.reporting.reporting(self, pcr, 'QallRAtot', Q) - if self.mm_rep_FLAG == 1: - self.QTOTSubBasinTSS.sample(((Q * 3600 * 24) / catchmenttotal(cellarea(), self.FlowDir)) * 1000) - #-Snow routing - if self.SnowRA_FLAG == 1 and self.SnowFLAG == 1: - SnowRA = self.routing.ROUT(pcr, SnowR, self.SnowRAold, self.FlowDir, self.kx) - self.SnowRAold = SnowRA - self.reporting.reporting(self, pcr, 'SnowRAtot', SnowRA) - if self.mm_rep_FLAG == 1: - self.QSNOWSubBasinTSS.sample(((SnowRA * 3600 * 24) / catchmenttotal(cellarea(), self.FlowDir)) *1000) - #-Rain routing - if self.RainRA_FLAG == 1: - RainRA = self.routing.ROUT(pcr, RainR, self.RainRAold, self.FlowDir, self.kx) - self.RainRAold = RainRA - self.reporting.reporting(self, pcr, 'RainRAtot', RainRA) - if self.mm_rep_FLAG == 1: - self.QRAINSubBasinTSS.sample(((RainRA * 3600 * 24) / catchmenttotal(cellarea(), self.FlowDir)) *1000) - #-Glacier routing - if self.GlacRA_FLAG == 1 and self.GlacFLAG == 1: - GlacRA = self.routing.ROUT(pcr, GlacR, self.GlacRAold, self.FlowDir, self.kx) - self.GlacRAold = GlacRA - self.reporting.reporting(self, pcr, 'GlacRAtot', GlacRA) - if self.mm_rep_FLAG == 1: - self.QGLACSubBasinTSS.sample(((GlacRA * 3600 * 24) / catchmenttotal(cellarea(), self.FlowDir)) *1000) - #-Baseflow routing - if self.BaseRA_FLAG == 1: - BaseRA = self.routing.ROUT(pcr, self.BaseR, self.BaseRAold, self.FlowDir, self.kx) - self.BaseRAold = BaseRA - self.reporting.reporting(self, pcr, 'BaseRAtot', BaseRA) - if self.mm_rep_FLAG == 1: - self.QBASFSubBasinTSS.sample(((BaseRA * 3600 * 24) / catchmenttotal(cellarea(), self.FlowDir)) *1000) + # -Report Total runoff + self.reporting.reporting(self, pcr, "TotRF", self.BaseR + RainR + SnowR + GlacR) + # -Water balance + if self.GroundFLAG == 1: + waterbalance = ( + Precip * (1 - self.GlacFrac) + + GlacMelt * self.GlacFrac + - ActETact + - GlacR + - SnowR + - RainR + - self.BaseR + - (self.SoilWater - OldSoilWater) + - (self.TotalSnowStore - OldTotalSnowStore) + - (self.Gw - GwOld) + ) + elif self.GroundFLAG == 0: + waterbalance = ( + Precip + - ActETact + - self.SeePage + - SnowR + - RainR + - self.BaseR + - (self.SoilWater - OldSoilWater) + - (self.TotalSnowStore - OldTotalSnowStore) + ) + self.reporting.reporting(self, pcr, "wbal", waterbalance) + # -Routing for lake and/or reservoir modules + if self.LakeFLAG == 1 or self.ResFLAG == 1: + # -Update storage in lakes/reservoirs (m3) with specific runoff + self.StorRES = self.StorRES + ifthenelse( + self.QFRAC == 0, + 0.001 * cellarea() * (self.BaseR + RainR + GlacR + SnowR), + 0, + ) + OldStorage = self.StorRES + # -Calculate lake/reservoir outflow volumes + if self.LakeFLAG == 1 and self.ResFLAG == 1: + tempvar = self.lakes.UpdateLakeHStore(self, pcr, pcrm) + LakeLevel = tempvar[0] + self.StorRES = tempvar[1] + LakeQ = self.lakes.QLake(self, pcr, LakeLevel) + ResQ = self.reservoirs.QRes(self, pcr) + Qout = ifthenelse( + self.ResID != 0, ResQ, ifthenelse(self.LakeID != 0, LakeQ, 0) + ) + elif self.LakeFLAG == 1: + tempvar = self.lakes.UpdateLakeHStore(self, pcr, pcrm) + LakeLevel = tempvar[0] + self.StorRES = tempvar[1] + Qout = self.lakes.QLake(self, pcr, LakeLevel) + else: + Qout = self.reservoirs.QRes(self, pcr) + # -Calculate volume available for routing (=outflow lakes/reservoir + cell specific runoff) + RunoffVolume = upstream(self.FlowDir, Qout) + ifthenelse( + self.QFRAC == 0, + 0, + 0.001 * cellarea() * (self.BaseR + RainR + GlacR + SnowR), + ) + # -Routing of total flow + tempvar = self.routing.ROUT( + self, pcr, RunoffVolume, self.QRAold, Qout, self.StorRES + ) + self.StorRES = tempvar[0] + Q = tempvar[1] + Qin = tempvar[2] + self.QRAold = Q + self.reporting.reporting(self, pcr, "QallRAtot", Q) + # -report flux in mm + if self.mm_rep_FLAG == 1: + self.QTOTSubBasinTSS.sample( + ((Q * 3600 * 24) / catchmenttotal(cellarea(), self.FlowDir)) * 1000 + ) + # -report lake and reservoir waterbalance + if self.LakeFLAG == 1 and config.getint("REPORTING", "Lake_wbal") == 1: + self.LakeInTSS.sample(Qin) + self.LakeOutTSS.sample(Qout) + self.LakeStorTSS.sample(self.StorRES) + if self.ResFLAG == 1 and config.getint("REPORTING", "Res_wbal") == 1: + self.ResInTSS.sample(Qin) + self.ResOutTSS.sample(Qout) + self.ResStorTSS.sample(self.StorRES) + # -Routing of individual contributers + # -Snow routing + if self.SnowRA_FLAG == 1 and self.SnowFLAG == 1: + self.SnowRAstor = self.SnowRAstor + ifthenelse( + self.QFRAC == 0, SnowR * 0.001 * cellarea(), 0 + ) + cQfrac = cover(self.SnowRAstor / OldStorage, 0) + cQout = cQfrac * Qout + cRunoffVolume = upstream(self.FlowDir, cQout) + ifthenelse( + self.QFRAC == 0, 0, 0.001 * cellarea() * SnowR + ) + tempvar = self.routing.ROUT( + self, pcr, cRunoffVolume, self.SnowRAold, cQout, self.SnowRAstor + ) + self.SnowRAstor = tempvar[0] + SnowRA = tempvar[1] + cQin = tempvar[2] + self.SnowRAold = SnowRA + self.reporting.reporting(self, pcr, "SnowRAtot", SnowRA) + if self.mm_rep_FLAG == 1: + self.QSNOWSubBasinTSS.sample( + ( + (SnowRA * 3600 * 24) + / catchmenttotal(cellarea(), self.FlowDir) + ) + * 1000 + ) + # -report lake and reservoir waterbalance + if self.LakeFLAG == 1 and config.getint("REPORTING", "Lake_wbal") == 1: + self.LakeSnowInTSS.sample(cQin) + self.LakeSnowOutTSS.sample(cQout) + self.LakeSnowStorTSS.sample(self.SnowRAstor) + if self.ResFLAG == 1 and config.getint("REPORTING", "Res_wbal") == 1: + self.ResSnowInTSS.sample(cQin) + self.ResSnowOutTSS.sample(cQout) + self.ResSnowStorTSS.sample(self.SnowRAstor) + # -Rain routing + if self.RainRA_FLAG == 1: + self.RainRAstor = self.RainRAstor + ifthenelse( + self.QFRAC == 0, RainR * 0.001 * cellarea(), 0 + ) + cQfrac = cover(self.RainRAstor / OldStorage, 0) + cQout = cQfrac * Qout + cRunoffVolume = upstream(self.FlowDir, cQout) + ifthenelse( + self.QFRAC == 0, 0, 0.001 * cellarea() * RainR + ) + tempvar = self.routing.ROUT( + self, pcr, cRunoffVolume, self.RainRAold, cQout, self.RainRAstor + ) + self.RainRAstor = tempvar[0] + RainRA = tempvar[1] + cQin = tempvar[2] + self.RainRAold = RainRA + self.reporting.reporting(self, pcr, "RainRAtot", RainRA) + if self.mm_rep_FLAG == 1: + self.QRAINSubBasinTSS.sample( + ( + (RainRA * 3600 * 24) + / catchmenttotal(cellarea(), self.FlowDir) + ) + * 1000 + ) + # -report lake and reservoir waterbalance + if self.LakeFLAG == 1 and config.getint("REPORTING", "Lake_wbal") == 1: + self.LakeRainInTSS.sample(cQin) + self.LakeRainOutTSS.sample(cQout) + self.LakeRainStorTSS.sample(self.RainRAstor) + if self.ResFLAG == 1 and config.getint("REPORTING", "Res_wbal") == 1: + self.ResRainInTSS.sample(cQin) + self.ResRainOutTSS.sample(cQout) + self.ResRainStorTSS.sample(self.RainRAstor) + # -Glacier routing + if self.GlacRA_FLAG == 1 and self.GlacFLAG == 1: + self.GlacRAstor = self.GlacRAstor + ifthenelse( + self.QFRAC == 0, GlacR * 0.001 * cellarea(), 0 + ) + cQfrac = cover(self.GlacRAstor / OldStorage, 0) + cQout = cQfrac * Qout + cRunoffVolume = upstream(self.FlowDir, cQout) + ifthenelse( + self.QFRAC == 0, 0, 0.001 * cellarea() * GlacR + ) + tempvar = self.routing.ROUT( + self, pcr, cRunoffVolume, self.GlacRAold, cQout, self.GlacRAstor + ) + self.GlacRAstor = tempvar[0] + GlacRA = tempvar[1] + cQin = tempvar[2] + self.GlacRAold = GlacRA + self.reporting.reporting(self, pcr, "GlacRAtot", GlacRA) + if self.mm_rep_FLAG == 1: + self.QGLACSubBasinTSS.sample( + ( + (GlacRA * 3600 * 24) + / catchmenttotal(cellarea(), self.FlowDir) + ) + * 1000 + ) + # -report lake and reservoir waterbalance + if self.LakeFLAG == 1 and config.getint("REPORTING", "Lake_wbal") == 1: + self.LakeGlacInTSS.sample(cQin) + self.LakeGlacOutTSS.sample(cQout) + self.LakeGlacStorTSS.sample(self.GlacRAstor) + if self.ResFLAG == 1 and config.getint("REPORTING", "Res_wbal") == 1: + self.ResGlacInTSS.sample(cQin) + self.ResGlacOutTSS.sample(cQout) + self.ResGlacStorTSS.sample(self.GlacRAstor) + # -Baseflow routing + if self.BaseRA_FLAG == 1: + self.BaseRAstor = self.BaseRAstor + ifthenelse( + self.QFRAC == 0, self.BaseR * 0.001 * cellarea(), 0 + ) + cQfrac = cover(self.BaseRAstor / OldStorage, 0) + cQout = cQfrac * Qout + cRunoffVolume = upstream(self.FlowDir, cQout) + ifthenelse( + self.QFRAC == 0, 0, 0.001 * cellarea() * self.BaseR + ) + tempvar = self.routing.ROUT( + self, pcr, cRunoffVolume, self.BaseRAold, cQout, self.BaseRAstor + ) + self.BaseRAstor = tempvar[0] + BaseRA = tempvar[1] + cQin = tempvar[2] + self.BaseRAold = BaseRA + self.reporting.reporting(self, pcr, "BaseRAtot", BaseRA) + if self.mm_rep_FLAG == 1: + self.QBASFSubBasinTSS.sample( + ( + (BaseRA * 3600 * 24) + / catchmenttotal(cellarea(), self.FlowDir) + ) + * 1000 + ) + # -report lake and reservoir waterbalance + if self.LakeFLAG == 1 and config.getint("REPORTING", "Lake_wbal") == 1: + self.LakeBaseInTSS.sample(cQin) + self.LakeBaseOutTSS.sample(cQout) + self.LakeBaseStorTSS.sample(self.BaseRAstor) + if self.ResFLAG == 1 and config.getint("REPORTING", "Res_wbal") == 1: + self.ResBaseInTSS.sample(cQin) + self.ResBaseOutTSS.sample(cQout) + self.ResBaseStorTSS.sample(self.BaseRAstor) + # -Normal routing module + elif self.RoutFLAG == 1: + # -Rout total runoff + Q = self.routing.ROUT( + pcr, + self.BaseR + RainR + GlacR + SnowR, + self.QRAold, + self.FlowDir, + self.kx, + ) + self.QRAold = Q + self.reporting.reporting(self, pcr, "QallRAtot", Q) + if self.mm_rep_FLAG == 1: + self.QTOTSubBasinTSS.sample( + ((Q * 3600 * 24) / catchmenttotal(cellarea(), self.FlowDir)) * 1000 + ) + # -Snow routing + if self.SnowRA_FLAG == 1 and self.SnowFLAG == 1: + SnowRA = self.routing.ROUT( + pcr, SnowR, self.SnowRAold, self.FlowDir, self.kx + ) + self.SnowRAold = SnowRA + self.reporting.reporting(self, pcr, "SnowRAtot", SnowRA) + if self.mm_rep_FLAG == 1: + self.QSNOWSubBasinTSS.sample( + ( + (SnowRA * 3600 * 24) + / catchmenttotal(cellarea(), self.FlowDir) + ) + * 1000 + ) + # -Rain routing + if self.RainRA_FLAG == 1: + RainRA = self.routing.ROUT( + pcr, RainR, self.RainRAold, self.FlowDir, self.kx + ) + self.RainRAold = RainRA + self.reporting.reporting(self, pcr, "RainRAtot", RainRA) + if self.mm_rep_FLAG == 1: + self.QRAINSubBasinTSS.sample( + ( + (RainRA * 3600 * 24) + / catchmenttotal(cellarea(), self.FlowDir) + ) + * 1000 + ) + # -Glacier routing + if self.GlacRA_FLAG == 1 and self.GlacFLAG == 1: + GlacRA = self.routing.ROUT( + pcr, GlacR, self.GlacRAold, self.FlowDir, self.kx + ) + self.GlacRAold = GlacRA + self.reporting.reporting(self, pcr, "GlacRAtot", GlacRA) + if self.mm_rep_FLAG == 1: + self.QGLACSubBasinTSS.sample( + ( + (GlacRA * 3600 * 24) + / catchmenttotal(cellarea(), self.FlowDir) + ) + * 1000 + ) + # -Baseflow routing + if self.BaseRA_FLAG == 1: + BaseRA = self.routing.ROUT( + pcr, self.BaseR, self.BaseRAold, self.FlowDir, self.kx + ) + self.BaseRAold = BaseRA + self.reporting.reporting(self, pcr, "BaseRAtot", BaseRA) + if self.mm_rep_FLAG == 1: + self.QBASFSubBasinTSS.sample( + ( + (BaseRA * 3600 * 24) + / catchmenttotal(cellarea(), self.FlowDir) + ) + * 1000 + ) # The main function is used to run the program from the command line + def main(argv=None): """ Perform command line execution of the model. @@ -1473,15 +2238,15 @@ global updateCols caseName = "wflow_ganga_sphy" runId = "run_default" - configfile="wflow_sphy.ini" - LogFileName="wflow.log" + configfile = "wflow_sphy.ini" + LogFileName = "wflow.log" _lastTimeStep = 0 _firstTimeStep = 0 - fewsrun=False - runinfoFile="runinfo.xml" - timestepsecs=86400 - wflow_cloneMap = 'wflow_subcatch.map' - NoOverWrite=1 + fewsrun = False + runinfoFile = "runinfo.xml" + timestepsecs = 86400 + wflow_cloneMap = "wflow_subcatch.map" + NoOverWrite = 1 loglevel = logging.DEBUG if argv is None: @@ -1493,86 +2258,111 @@ ## Main model starts here ######################################################################## try: - opts, args = getopt.getopt(argv, 'c:QXS:F:hC:Ii:T:R:u:s:P:p:Xx:U:fl:L:') + opts, args = getopt.getopt(argv, "c:QXS:F:hC:Ii:T:R:u:s:P:p:Xx:U:fl:L:") except getopt.error, msg: pcrut.usage(msg) for o, a in opts: - if o == '-F': + if o == "-F": runinfoFile = a fewsrun = True - if o == '-C': caseName = a - if o == '-R': runId = a - if o == '-L': LogFileName = a - if o == '-l': exec "loglevel = logging." + a - if o == '-c': configfile = a - if o == '-s': timestepsecs = int(a) - if o == '-h': usage() - if o == '-f': NoOverWrite = 0 + if o == "-C": + caseName = a + if o == "-R": + runId = a + if o == "-L": + LogFileName = a + if o == "-l": + exec "loglevel = logging." + a + if o == "-c": + configfile = a + if o == "-s": + timestepsecs = int(a) + if o == "-h": + usage() + if o == "-f": + NoOverWrite = 0 - - if fewsrun: - ts = getTimeStepsfromRuninfo(runinfoFile,timestepsecs) + ts = getTimeStepsfromRuninfo(runinfoFile, timestepsecs) starttime = getStartTimefromRuninfo(runinfoFile) - if (ts): - _lastTimeStep = ts# * 86400/timestepsecs + if ts: + _lastTimeStep = ts # * 86400/timestepsecs _firstTimeStep = 1 else: print "Failed to get timesteps from runinfo file: " + runinfoFile sys.exit(2) else: - starttime = dt.datetime(1990,01,01) + starttime = dt.datetime(1990, 01, 01) if _lastTimeStep < _firstTimeStep: - print "The starttimestep (" + str(_firstTimeStep) +") is smaller than the last timestep (" + str(_lastTimeStep) + ")" + print "The starttimestep (" + str( + _firstTimeStep + ) + ") is smaller than the last timestep (" + str(_lastTimeStep) + ")" usage() - myModel = WflowModel(wflow_cloneMap, caseName,runId,configfile) - dynModelFw = wf_DynamicFramework(myModel, _lastTimeStep,firstTimestep=_firstTimeStep,datetimestart=starttime) - dynModelFw.createRunId(NoOverWrite=NoOverWrite,logfname=LogFileName,level=loglevel,doSetupFramework=False) + myModel = WflowModel(wflow_cloneMap, caseName, runId, configfile) + dynModelFw = wf_DynamicFramework( + myModel, _lastTimeStep, firstTimestep=_firstTimeStep, datetimestart=starttime + ) + dynModelFw.createRunId( + NoOverWrite=NoOverWrite, + logfname=LogFileName, + level=loglevel, + doSetupFramework=False, + ) - for o, a in opts: - if o == '-P': - left = a.split('=')[0] - right = a.split('=')[1] - configset(myModel.config,'variable_change_once',left,right,overwrite=True) - if o == '-p': - left = a.split('=')[0] - right = a.split('=')[1] - configset(myModel.config,'variable_change_timestep',left,right,overwrite=True) - if o == '-X': configset(myModel.config,'model','OverWriteInit','1',overwrite=True) - if o == '-I': configset(myModel.config,'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',"0",overwrite=True) - if o == '-Q': configset(myModel.config,'model','ExternalQbase','1',overwrite=True) - if o == '-U': - configset(myModel.config,'model','updateFile',a,overwrite=True) - configset(myModel.config,'model','updating',"1",overwrite=True) - if o == '-u': - exec "zz =" + a + if o == "-P": + left = a.split("=")[0] + right = a.split("=")[1] + configset( + myModel.config, "variable_change_once", left, right, overwrite=True + ) + if o == "-p": + left = a.split("=")[0] + right = a.split("=")[1] + configset( + myModel.config, "variable_change_timestep", left, right, overwrite=True + ) + if o == "-X": + configset(myModel.config, "model", "OverWriteInit", "1", overwrite=True) + if o == "-I": + configset(myModel.config, "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", "0", overwrite=True) + if o == "-Q": + configset(myModel.config, "model", "ExternalQbase", "1", overwrite=True) + if o == "-U": + configset(myModel.config, "model", "updateFile", a, overwrite=True) + configset(myModel.config, "model", "updating", "1", overwrite=True) + if o == "-u": + exec "zz =" + a updateCols = zz - if o == '-T': - configset(myModel.config, 'run', 'endtime', a, overwrite=True) - if o == '-S': - configset(myModel.config, 'run', 'starttime', a, overwrite=True) + if o == "-T": + configset(myModel.config, "run", "endtime", a, overwrite=True) + if o == "-S": + configset(myModel.config, "run", "starttime", a, overwrite=True) dynModelFw.setupFramework() dynModelFw.logger.info("Command line: " + str(argv)) dynModelFw._runInitial() dynModelFw._runResume() - #dynModelFw._runDynamic(0,0) + # dynModelFw._runDynamic(0,0) dynModelFw._runDynamic(_firstTimeStep, _lastTimeStep) dynModelFw._runSuspend() dynModelFw._wf_shutdown() - os.chdir("../../")