#!/usr/bin/python # Wflow is Free software, see below: # # Copyright (c) J. Schellekens/Deltares 2005-2013 # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see . #TODO: split off routing import numpy import sys import os import os.path import shutil, glob import getopt import datetime as dt from wflow.wf_DynamicFramework import * from wflow.wf_DynamicFramework import * from wflow.wflow_adapt import * from wflow.wflow_adapt import * import pcraster as pcr #from wflow.wflow_lib import reporting #import scipy #import pcrut wflow = "wflow_sphy" #: columns used in updating updateCols = [] #: columns used in updating """ Column used in updating """ def usage(*args): """ Print usage information - *args: command line arguments given """ sys.stdout = sys.stderr 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 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 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. This function is specific for each model and **must** be present. :var self.SurfaceRunoff: Surface runoff in the kin-wave resrvoir [m^3/s] :var self.WaterLevel: Water level in the kin-wave resrvoir [m] :var self.DrySnow: Snow pack [mm] :var self.FreeWater: Available free water [mm] :var self.UpperZoneStorage: Water in the upper zone [mm] :var self.LowerZoneStorage: Water in the lower zone [mm] :var self.SoilMoisture: Soil moisture [mm] :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 # 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') return states # The following are made to better connect to deltashell/openmi 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')) 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 = [] #Static model parameters e.g. #modelparameters.append(self.ParamType(name="RunoffGeneratingGWPerc",stack="intbl/RunoffGeneratingGWPerc.tbl",type="static",default=0.1)) # 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 # 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")) if self.OverWriteInit: self.logger.info("Saving initial conditions over start conditions...") self.wf_suspend(os.path.join(self.SaveDir,"instate")) if self.fewsrun: self.logger.info("Saving initial conditions for FEWS...") self.wf_suspend(os.path.join(self.Dir, "outstate")) def initial(self): global statistics global multpars global updateCols setglobaloption("unittrue") self.thestep = scalar(0) #: files to be used in case of timesries (scalar) input to the model # #: name of the tss file with precipitation data ("../intss/P.tss") # self.precipTss = "../intss/P.tss" # self.evapTss="../intss/PET.tss" #: name of the tss file with potential evap data ("../intss/PET.tss") # self.tempTss="../intss/T.tss" #: name of the tss file with temperature data ("../intss/T.tss") # self.inflowTss="../intss/Inflow.tss" #: NOT TESTED name of the tss file with inflow data ("../intss/Inflow.tss") # self.SeepageTss="../intss/Seepage.tss" #: NOT TESTED name of the tss file with seepage data ("../intss/Seepage.tss")" self.logger.info("running for " + str(self.nrTimeSteps()) + " timesteps") # Set and get defaults from ConfigFile here ################################### # self.scalarInput = int(configget(self.config,"model","ScalarInput","0")) # self.Tslice = int(configget(self.config,"model","Tslice","1")) # self.interpolMethod = configget(self.config,"model","InterpolationMethod","inv") self.reinit = int(configget(self.config,"run","reinit","0")) self.fewsrun = int(configget(self.config,"run","fewsrun","0")) self.OverWriteInit = int(configget(self.config,"model","OverWriteInit","0")) # self.updating = int(configget(self.config,"model","updating","0")) # self.updateFile = configget(self.config,"model","updateFile","no_set") # self.sCatch = int(configget(self.config,"model","sCatch","0")) # self.intbl = configget(self.config,"model","intbl","intbl") # self.P_style = int(configget(self.config,"model","P_style","1")) # self.PET_style = int(configget(self.config,"model","PET_style","1")) # self.TEMP_style = int(configget(self.config,"model","TEMP_style","1")) # self.modelSnow = int(configget(self.config,"model","ModelSnow","1")) # sizeinmetres = int(configget(self.config,"layout","sizeinmetres","0")) # alf = float(configget(self.config,"model","Alpha","60")) # Qmax = float(configget(self.config,"model","AnnualDischarge","300")) # self.UpdMaxDist =float(configget(self.config,"model","UpdMaxDist","100")) # self.MaxUpdMult =float(configget(self.config,"model","MaxUpdMult","1.3")) # self.MinUpdMult =float(configget(self.config,"model","MinUpdMult","0.7")) # self.UpFrac =float(configget(self.config,"model","UpFrac","0.8")) # self.ExternalQbase=int(configget(self.config,'model','ExternalQbase','0')) # self.SetKquickFlow=int(configget(self.config,'model','SetKquickFlow','0')) # self.MassWasting = int(configget(self.config,"model","MassWasting","0")) # self.SubCatchFlowOnly = int(configget(self.config, 'model', 'SubCatchFlowOnly', '0')) # 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')) # import the required modules import datetime, calendar import wflow.reporting as reporting import wflow.timecalc as timecalc import wflow.ET as ET import wflow.rootzone as rootzone import wflow.subzone as subzone #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 glacier # glacier melting processes self.glacier = glacier del glacier if self.SnowFLAG == 1: import snow # snow melt processes self.snow = snow del snow if self.RoutFLAG == 1: import 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 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)) 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.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: # 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'] return lst 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")) 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) # 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)) #-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)) #-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)) # 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)) #-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 #-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 #-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) #-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) #-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) 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) #-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. """ global multpars global updateCols caseName = "wflow_ganga_sphy" runId = "run_default" 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 loglevel = logging.DEBUG if argv is None: argv = sys.argv[1:] if len(argv) == 0: usage() return ## 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:') except getopt.error, msg: pcrut.usage(msg) for o, a in opts: 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 fewsrun: ts = getTimeStepsfromRuninfo(runinfoFile,timestepsecs) starttime = getStartTimefromRuninfo(runinfoFile) 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) if _lastTimeStep < _firstTimeStep: 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) 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 updateCols = zz 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(_firstTimeStep, _lastTimeStep) dynModelFw._runSuspend() dynModelFw._wf_shutdown() os.chdir("../../") if __name__ == "__main__": main()