#!/usr/bin/python # Wflow is Free software, see below: # # Copyright (c) J. Schellekens 2005-2011 # # 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 . """ Run the wflow_cqf hydrological model.. wflow_cqf is a model develope oftare the cqflow model and is developed specifically for the rio chiquito catchment usage :: wflow_cqf [-h][-v level][-F runinfofile][-L logfile][-C casename][-R runId] [-c configfile][-T last_step][-S first_step][-s seconds][-W][-E][-N][-U discharge] [-P parameter multiplication][-X][-f][-I][-i tbl_dir][-x subcatchId][-u updatecols] [-p inputparameter multiplication] -F: if set wflow is expected to be run by FEWS. It will determine the timesteps from the runinfo.xml file and save the output initial conditions to an alternate location. The runinfo.xml file should be located in the inmaps directory of the case. Also set fewsrun=1 in the .ini file! -X: save state at the end of the run over the initial conditions at the start -f: Force overwrite of existing results -T: Set last timestep -S: Set the start timestep (default = 1) -N: No lateral flow, use runoff response function to generate fast runoff -s: Set the model timesteps in seconds -I: re-initialize the initial model conditions with default -i: Set input table directory (default is intbl) -x: run for subcatchment only (e.g. -x 1) -C: set the name of the case (directory) to run -R: set the name runId within the current case -L: set the logfile -E: Switch on reinfiltration of overland flow -c: name of wflow the configuration file (default: Casename/wflow_cqf.ini). -h: print usage information -W: If set, this flag indicates that an ldd is created for the water level for each timestep. If not the water is assumed to flow according to the DEM. Wflow will run a lot slower with this option. Most of the time (shallow soil, steep topography) you do not need this option. Also, if you need it you migth actually need another model. -U: The argument to this option should be a .tss file with measured discharge in [m^3/s] which the program will use to update the internal state to match the measured flow. The number of columns in this file should match the number of gauges. -u: list of gauges/columns to use in update. Format: -u [1 , 4 ,13] The above example uses column 1, 4 and 13 Note that this also sets the order in which the updating takes place! In general specify downstream gauges first. -P: set parameter multiply dictionary (e.g: -P {'self.FirstZoneDepth' : 1.2} to increase self.FirstZoneDepth by 20%, multiply with 1.2) -p: set input parameter (dynamic, e.g. precip) multiply dictionary (e.g: -p {'self.Precipitation' : 1.2} to increase Precipitation by 20%, multiply with 1.2) -v: set verbosity level $Author: schelle $ $Id: wflow_cqf.py 909 2014-01-16 15:23:32Z schelle $ $Rev: 909 $ """ import numpy #import pcrut import os import os.path import shutil, glob import getopt from wflow.wf_DynamicFramework import * from wflow.wflow_funcs import * from wflow.wflow_adapt import * import scipy import ConfigParser wflow = "wflow_cqf: " wflowVersion = "$Revision: 909 $ $Date: 2014-01-16 16:23:32 +0100 (Thu, 16 Jan 2014) $" updateCols = [] # Dictionary with parameters and multipliers (used in calibration) multpars = {} multdynapars = {} def usage(*args): sys.stdout = sys.stderr for msg in args: print msg print __doc__ sys.exit(0) def actEvap_SBM(RootingDepth,WTable, UStoreDepth,FirstZoneDepth, PotTrans,smoothpar): """ Actual evaporation function: - first try to get demand from the saturated zone, using the rootingdepth as a limiting factor - secondly try to get the remaining water from the unsaturated store Input: - RootingDepth,WTable, UStoreDepth,FirstZoneDepth, PotTrans Output: - ActEvap, FirstZoneDepth, UStoreDepth ActEvapUStore .. todo:: add option to take length of roots in saturated zone into account """ # Step 1 from saturated zone, use rootingDepth as a limiting factor #rootsinWater = WTable < RootingDepth #ActEvapSat = ifthenelse(rootsinWater,min(PotTrans,FirstZoneDepth),0.0) # new method: # use sCurve to determine if the roots are wet.At the moment this ise set # to be a 0-1 curve wetroots = sCurve(WTable,a=RootingDepth,c=smoothpar) ActEvapSat = min(PotTrans * wetroots,FirstZoneDepth) FirstZoneDepth = FirstZoneDepth - ActEvapSat RestPotEvap = PotTrans - ActEvapSat # now try unsat store AvailCap = min(1.0,max (0.0,(WTable - RootingDepth)/(RootingDepth + 1.0))) #AvailCap = max(0.0,ifthenelse(WTable < RootingDepth, WTable/RootingDepth, RootingDepth/WTable)) MaxExtr = AvailCap * UStoreDepth ActEvapUStore = min(MaxExtr,RestPotEvap,UStoreDepth) UStoreDepth = UStoreDepth - ActEvapUStore ActEvap = ActEvapSat + ActEvapUStore return ActEvap, FirstZoneDepth, UStoreDepth, ActEvapUStore class WflowModel(DynamicModel): def __init__(self, cloneMap,Dir,RunDir,configfile): DynamicModel.__init__(self) setclone(Dir + "/staticmaps/" + cloneMap) self.runId=RunDir self.caseName=Dir self.Dir = Dir + "/" self.configfile = configfile def updateRunOff(self): """ Updates the kinematic wave reservoir. Should be run after updates to Q """ self.WaterLevel=(self.Alpha*pow(self.SurfaceRunoff,self.Beta))/self.Bw # wetted perimeter (m) P=self.Bw+(2*self.WaterLevel) # Alpha self.Alpha=self.AlpTerm*pow(P,self.AlpPow) self.OldKinWaveVolume = self.KinWaveVolume self.KinWaveVolume = self.WaterLevel * self.Bw * self.DCL 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. - CanopyStorage is any needed for subdaily steps """ states = ['SurfaceRunoff', 'WaterLevel', 'FirstZoneDepth', 'UStoreDepth', 'CanopyStorage'] return states def supplyCurrentTime(self): """ gets the current time in seconds after the start of the run """ return self.currentTimeStep() * self.timestepsecs def readtblDefault(self,pathtotbl,landuse,subcatch,soil, default): """ First check if a prepared map of the same name is present in the staticmaps directory. next try to read a tbl file to match a landuse, catchment and soil map. Returns the default value if the tbl file is not found. Input: - pathtotbl,landuse,subcatch,soil, default Output: - map constructed from tbl file or map with default value """ mapname = os.path.dirname(pathtotbl) + "/../staticmaps/" + os.path.splitext(os.path.basename(pathtotbl))[0]+".map" if os.path.exists(mapname): self.logger.info("reading map parameter file: " + mapname) rest = cover(readmap(mapname),default) else: if os.path.isfile(pathtotbl): rest=lookupscalar(pathtotbl,landuse,subcatch,soil) # self.logger.info("Creating map from table: " + pathtotbl) else: self.logger.warn("tbl file not found (" + pathtotbl + ") returning default value: " + str(default)) rest = scalar(default) return rest def suspend(self): self.logger.info("Saving initial conditions...") self.wf_suspend(self.SaveDir + "/outstate/") if self.OverWriteInit: self.logger.info("Saving initial conditions over start conditions...") self.wf_suspend(self.SaveDir + "/instate/") if self.fewsrun: self.logger.info("Saving initial conditions for FEWS...") self.wf_suspend(self.Dir + "/outstate/") report(self.CumInwaterMM,self.SaveDir + "/outsum/CumInwaterMM.map") report(self.CumReinfilt,self.SaveDir + "/outsum/CumReinfilt.map") report(self.CumPrec,self.SaveDir + "/outsum/CumPrec.map") report(self.CumEvap,self.SaveDir + "/outsum/CumEvap.map") report(self.CumInt,self.SaveDir + "/outsum/CumInt.map") report(self.CumLeakage,self.SaveDir + "/outsum/CumLeakage.map") report(self.CumPotenEvap,self.SaveDir + "/outsum/CumPotenEvap.map") report(self.CumExfiltWater,self.SaveDir + "/outsum/CumExfiltWater.map") report(self.watbal,self.SaveDir + "/outsum/watbal.map") def initial(self): """Initial part of the model, executed only once """ global statistics global multpars self.thestep = scalar(0) self.basetimestep = 86400 self.SSSF=False setglobaloption("unittrue") intbl = "intbl" self.precipTss="/intss/P.tss" self.evapTss="/intss/PET.tss" self.tempTss="/intss/T.tss" self.inflowTss="/intss/Inflow.tss" self.SeepageTss="/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,"model","reinit","0")) self.fewsrun = int(configget(self.config,"model","fewsrun","0")) self.OverWriteInit = int(configget(self.config,"model","OverWriteInit","0")) self.updating = int(configget(self.config,"model","updating","0")) self.updateFile = configget(self.config,"model","updateFile","no_set") self.sCatch = int(configget(self.config,"model","sCatch","0")) self.intbl = configget(self.config,"model","intbl","intbl") self.timestepsecs = int(configget(self.config,"model","timestepsecs","86400")) self.modelSnow = int(configget(self.config,"model","ModelSnow","1")) sizeinmetres = int(configget(self.config,"layout","sizeinmetres","0")) alf = float(configget(self.config,"model","Alpha","60")) #TODO: make this into a list for all gauges or a map Qmax = float(configget(self.config,"model","AnnualDischarge","300")) self.UpdMaxDist =float(configget(self.config,"model","UpdMaxDist","100")) self.ExternalQbase=int(configget(self.config,'model','ExternalQbase','0')) self.waterdem=int(configget(self.config,'model','waterdem','0')) WIMaxScale=float(configget(self.config,'model','WIMaxScale','0.8')) self.reInfilt=int(configget(self.config,'model','reInfilt','0')) # static maps to use (normally default) wflow_subcatch = configget(self.config,"model","wflow_subcatch","staticmaps/wflow_subcatch.map") wflow_dem = configget(self.config,"model","wflow_dem","staticmaps/wflow_dem.map") wflow_ldd = configget(self.config,"model","wflow_ldd","staticmaps/wflow_ldd.map") wflow_river = configget(self.config,"model","wflow_river","staticmaps/wflow_river.map") wflow_riverlength = configget(self.config,"model","wflow_riverlength","staticmaps/wflow_riverlength.map") wflow_riverlength_fact = configget(self.config,"model","wflow_riverlength_fact","staticmaps/wflow_riverlength_fact.map") wflow_landuse = configget(self.config,"model","wflow_landuse","staticmaps/wflow_landuse.map") wflow_soil = configget(self.config,"model","wflow_soil","staticmaps/wflow_soil.map") wflow_gauges = configget(self.config,"model","wflow_gauges","staticmaps/wflow_gauges.map") wflow_inflow = configget(self.config,"model","wflow_inflow","staticmaps/wflow_inflow.map") wflow_mgauges = configget(self.config,"model","wflow_mgauges","staticmaps/wflow_mgauges.map") # 2: Input base maps ######################################################## subcatch=ordinal(readmap(self.Dir + wflow_subcatch)) # Determines the area of calculations (all cells > 0) subcatch = ifthen(subcatch > 0, subcatch) if self.sCatch > 0: subcatch = ifthen(subcatch == sCatch,subcatch) self.Altitude=readmap(self.Dir + wflow_dem)# * scalar(defined(subcatch)) # DEM self.TopoLdd=readmap(self.Dir + wflow_ldd) # Local self.TopoId=ordinal(readmap(self.Dir + wflow_subcatch)) # area map self.River=cover(boolean(readmap(self.Dir + wflow_river)),0) self.RiverLength=pcrut.readmapSave(self.Dir + wflow_riverlength,0.0) # Factor to multiply riverlength with (defaults to 1.0) self.RiverLengthFac=pcrut.readmapSave(self.Dir + wflow_riverlength_fact,1.0) # read landuse and soilmap and make sure there are no missing points related to the # subcatchment map. Currently sets the lu and soil type type to 1 self.LandUse=readmap(self.Dir + wflow_landuse) self.LandUse=cover(self.LandUse,nominal(ordinal(subcatch) > 0)) self.Soil=readmap(self.Dir + wflow_soil) self.Soil=cover(self.Soil,nominal(ordinal(subcatch) > 0)) self.OutputLoc=ordinal(readmap(self.Dir + wflow_gauges)) # location of output gauge(s) self.InflowLoc=pcrut.readmapSave(self.Dir + wflow_inflow,0.0) # location abstractions/inflows. self.SeepageLoc=pcrut.readmapSave(self.Dir + wflow_inflow,0.0) # location abstractions/inflows. # Experimental self.RunoffGenSigmaFunction = int(configget(self.config,'model','RunoffGenSigmaFunction','0')) self.RunoffGeneratingGWPerc = float(configget(self.config,'defaultfortbl','RunoffGeneratingGWPerc','0.1')) self.RunoffGeneratingThickness = float(configget(self.config,'defaultfortbl','RunoffGeneratingThickness','0.0')) if self.scalarInput: self.gaugesMap=readmap(self.Dir + wflow_mgauges) # location of rainfall/evap/temp gauge(s) self.OutputId=ordinal(readmap(self.Dir + wflow_subcatch)) # location of subcatchment # Temperature correction poer cell to add self.TempCor=pcrut.readmapSave(self.Dir + configget(self.config,"model","TemperatureCorrectionMap","staticmaps/wflow_tempcor.map"),0.0) self.ZeroMap=0.0*scalar(subcatch) #map with only zero's # 3: Input time series ################################################### self.P_mapstack=self.Dir + configget(self.config,"inputmapstacks","Precipitation","/inmaps/P") # timeseries for rainfall self.PET_mapstack=self.Dir + configget(self.config,"inputmapstacks","EvapoTranspiration","/inmaps/PET") # timeseries for rainfall"/inmaps/PET" # potential evapotranspiration self.TEMP_mapstack=self.Dir + configget(self.config,"inputmapstacks","Temperature","/inmaps/TEMP") # timeseries for rainfall "/inmaps/TEMP" # global radiation self.Inflow_mapstack=self.Dir + configget(self.config,"inputmapstacks","Inflow","/inmaps/IF") # timeseries for rainfall "/inmaps/IF" # in/outflow locations (abstractions) self.RH_mapstack=self.Dir + configget(self.config,"inputmapstacks","RH","/inmaps/RH") # timeseries for rainfall "/inmaps/SE" # in/outflow locations (abstractions) self.WindSpeed_mapstack=self.Dir + configget(self.config,"inputmapstacks","WindSpeed","/inmaps/wins") # timeseries for rainfall "/inmaps/SE" # in/outflow locations (abstractions) self.RAD_mapstack=self.Dir + configget(self.config,"inputmapstacks","Radiation","/inmaps/RAD") # timeseries for rainfall "/inmaps/SE" # in/outflow locations (abstractions) self.Seepage_mapstack=self.Dir + configget(self.config,"inputmapstacks","Seepage","/inmaps/SE") # timeseries for rainfall "/inmaps/SE" # in/outflow locations (abstractions) self.HP_mapstack=self.Dir + configget(self.config,"inputmapstacks","HP","/inmaps/HP") # timeseries for rainfall "/inmaps/SE" # in/outflow locations (abstractions) self.WaterCatch_mapstack=self.Dir + configget(self.config,"inputmapstacks","WaterCatch","/inmaps/WC") # timeseries for rainfall "/inmaps/SE" # in/outflow locations (abstractions) # 3: Input time series ################################################### # Set static initial values here ######################################### self.SoilAlbedo = 0.1 # Not used at the moment self.pi = 3.1416 self.e = 2.7183 self.SScale = 100.0 self.Latitude = ycoordinate(boolean(self.Altitude)) self.Longitude = xcoordinate(boolean(self.Altitude)) self.logger.info("Linking parameters to landuse, catchment and soil...") self.RunoffGeneratingGWPerc=self.readtblDefault(self.Dir + "/" + self.intbl + "/RunoffGeneratingGWPerc.tbl",self.LandUse,subcatch,self.Soil,self.RunoffGeneratingGWPerc) self.RunoffGeneratingThickness=self.readtblDefault(self.Dir + "/" + self.intbl + "/RunoffGeneratingThickness.tbl",self.LandUse,subcatch,self.Soil,self.RunoffGeneratingThickness) self.Cmax=self.readtblDefault(self.Dir + "/" + self.intbl + "/MaxCanopyStorage.tbl",self.LandUse,subcatch,self.Soil,1.0) self.EoverR=self.readtblDefault(self.Dir + "/" + self.intbl + "/EoverR.tbl",self.LandUse,subcatch,self.Soil,0.1) # self.Albedo=lookupscalar(self.Dir + "\intbl\Albedo.tbl",self.LandUse) # Not used anymore self.CanopyGapFraction=self.readtblDefault(self.Dir + "/" + self.intbl + "/CanopyGapFraction.tbl",self.LandUse,subcatch,self.Soil,0.1) self.RootingDepth=self.readtblDefault(self.Dir + "/" + self.intbl + "/RootingDepth.tbl",self.LandUse,subcatch,self.Soil,750.0) #rooting depth #: rootdistpar determien how roots are linked to water table.The number shoudl be negative. A high number means that all roots are wet if #: the water table is above the lowest part of the roots. A lower number smooths this. self.rootdistpar=self.readtblDefault(self.Dir + "/" + self.intbl + "/rootdistpar.tbl",self.LandUse,subcatch,self.Soil,-80000.0) #rrootdistpar # Soil parameters # infiltration capacity if the soil [mm/day] self.InfiltCapSoil=self.readtblDefault(self.Dir + "/" + self.intbl + "/InfiltCapSoil.tbl",self.LandUse,subcatch,self.Soil,100.0) * self.timestepsecs/self.basetimestep self.CapScale=self.readtblDefault(self.Dir + "/" + self.intbl + "/CapScale.tbl",self.LandUse,subcatch,self.Soil,100.0) # # infiltration capacity of the compacted self.InfiltCapPath=self.readtblDefault(self.Dir + "/" + self.intbl + "/InfiltCapPath.tbl",self.LandUse,subcatch,self.Soil,10.0) * self.timestepsecs/self.basetimestep self.MaxLeakage=self.readtblDefault(self.Dir + "/" + self.intbl + "/MaxLeakage.tbl",self.LandUse,subcatch,self.Soil,0.0) * self.timestepsecs/self.basetimestep # areas (paths) in [mm/day] # Fraction area with compacted soil (Paths etc.) self.PathFrac=self.readtblDefault(self.Dir + "/" + self.intbl + "/PathFrac.tbl",self.LandUse,subcatch,self.Soil,0.01) # thickness of the soil self.FirstZoneThickness = self.readtblDefault(self.Dir + "/" + self.intbl + "/FirstZoneCapacity.tbl",self.LandUse,subcatch,self.Soil,2000.0) self.thetaR = self.readtblDefault(self.Dir + "/" + self.intbl + "/thetaR.tbl",self.LandUse,subcatch,self.Soil,0.01) self.thetaS = self.readtblDefault(self.Dir + "/" + self.intbl + "/thetaS.tbl",self.LandUse,subcatch,self.Soil,0.6) # minimum thickness of soild self.FirstZoneMinCapacity = self.readtblDefault(self.Dir + "/" + self.intbl + "/FirstZoneMinCapacity.tbl",self.LandUse,subcatch,self.Soil,500.0) # FirstZoneKsatVer = $2\inmaps\FirstZoneKsatVer.map self.FirstZoneKsatVer=self.readtblDefault(self.Dir + "/" + self.intbl + "/FirstZoneKsatVer.tbl",self.LandUse,subcatch,self.Soil,3000.0) * self.timestepsecs/self.basetimestep self.Beta = scalar(0.6) # For sheetflow self.M=self.readtblDefault(self.Dir + "/" + self.intbl + "/M.tbl" ,self.LandUse,subcatch,self.Soil,300.0) # Decay parameter in Topog_cqf self.N=self.readtblDefault(self.Dir + "/" + self.intbl + "/N.tbl",self.LandUse,subcatch,self.Soil,0.072) # Manning overland flow self.NRiver=self.readtblDefault(self.Dir + "/" + self.intbl + "/N_River.tbl",self.LandUse,subcatch,self.Soil,0.036) # Manning river self.WaterFrac=self.readtblDefault(self.Dir + "/" + self.intbl + "/WaterFrac.tbl",self.LandUse,subcatch,self.Soil,0.0) # Fraction Open water #cqflow specific stuff self.Albedo=self.readtblDefault(self.Dir + "/" + self.intbl + "/Albedo.tbl",self.LandUse,subcatch,self.Soil,0.18) # self.LeafAreaIndex=self.readtblDefault(self.Dir + "/" + self.intbl + "/LeafAreaIndex.tbl",self.LandUse,subcatch,self.Soil,4.0) # self.WindSpeedHeigth=self.readtblDefault(self.Dir + "/" + self.intbl + "/WindSpeedHeigth.tbl",self.LandUse,subcatch,self.Soil,4.0) # self.VegetationHeigth=self.readtblDefault(self.Dir + "/" + self.intbl + "/VegetationHeigth.tbl",self.LandUse,subcatch,self.Soil,4.0) # self.xl,self.yl,self.reallength = pcrut.detRealCellLength(self.ZeroMap,sizeinmetres) self.Slope= slope(self.Altitude) #self.Slope=ifthen(boolean(self.TopoId),max(0.001,self.Slope*celllength()/self.reallength)) self.Slope=max(0.001,self.Slope*celllength()/self.reallength) Terrain_angle=scalar(atan(self.Slope)) # Multiply parameters with a factor (for calibration etc) -P option in command line for k, v in multpars.iteritems(): estr = k + "=" + k + "*" + str(v) self.logger.info("Parameter multiplication: " + estr) exec estr self.N=ifthenelse(self.River, self.NRiver, self.N) # Determine river width from DEM, upstream area and yearly average discharge # Scale yearly average Q at outlet with upstream are to get Q over whole catchment # Alf ranges from 5 to > 60. 5 for hardrock. large values for sediments # "Noah J. Finnegan et al 2005 Controls on the channel width of rivers: # Implications for modeling fluvial incision of bedrock" upstr = catchmenttotal(1, self.TopoLdd) Qscale = upstr/mapmaximum(upstr) * Qmax W = (alf * (alf + 2.0)**(0.6666666667))**(0.375) * Qscale**(0.375) * (max(0.0001,windowaverage(self.Slope,celllength() * 4.0)))**(-0.1875) * self.N **(0.375) RiverWidth = W # soil thickness based on topographical index (see Environmental modelling: finding simplicity in complexity) # 1: calculate wetness index # 2: Scale the capacity (now actually a max capacity) based on the index, also apply a minmum capacity WI = ln(accuflux(self.TopoLdd,1)/self.Slope) # Topographical wetnesss. Scale WI by zone/subcatchment assuming these ara also geological units WIMax = areamaximum(WI, self.TopoId) * WIMaxScale self.FirstZoneThickness = max(min(self.FirstZoneThickness,(WI/WIMax) * self.FirstZoneThickness), self.FirstZoneMinCapacity) self.FirstZoneCapacity = self.FirstZoneThickness * (self.thetaS -self.thetaR) # limit roots to top 99% of first zone self.RootingDepth = min(self.FirstZoneThickness * 0.99,self.RootingDepth) # subgrid runoff generation self.DemMax=readmap(self.Dir + "/staticmaps/wflow_demmax") self.DrainageBase=readmap(self.Dir + "/staticmaps/wflow_demmin") self.CC = min(100.0,-log(1.0/0.1 - 1)/min(-0.1,self.DrainageBase - self.Altitude)) #if maptotal(self.RunoffGeneratingThickness <= 0.0): self.GWScale = (self.DemMax-self.DrainageBase)/self.FirstZoneThickness / self.RunoffGeneratingGWPerc #else: # self.GWScale = (self.DemMax-self.DrainageBase)/min(self.RunoffGeneratingThickness, self.FirstZoneThickness) # Which columns/gauges to use/ignore in updating self.UpdateMap = self.ZeroMap if self.updating: touse = numpy.zeros(gaugear.shape,dtype='int') for thecol in updateCols: idx = (gaugear == thecol).nonzero() touse[idx] = thecol self.UpdateMap = numpy2pcr(Nominal,touse,0.0) # Calulate distance to updating points (upstream) annd use to scale the correction # ldddist returns zero for cell at the gauges so add 1.0 tp result self.DistToUpdPt = cover(min(ldddist(self.TopoLdd,boolean(cover(self.UpdateMap,0)),1) * self.reallength/celllength(),self.UpdMaxDist),self.UpdMaxDist) # Initializing of variables self.logger.info("Initializing of model variables..") self.TopoLdd=lddmask(self.TopoLdd,boolean(self.TopoId)) catchmentcells=maptotal(scalar(self.TopoId)) # Used to seperate output per LandUse/management classes OutZones = self.LandUse self.QMMConv = self.timestepsecs/(self.reallength * self.reallength * 0.001) #m3/s --> mm self.ToCubic = (self.reallength * self.reallength * 0.001) / self.timestepsecs # m3/s self.KinWaveVolume=self.ZeroMap self.OldKinWaveVolume=self.ZeroMap self.sumprecip=self.ZeroMap #accumulated rainfall for water balance self.sumevap=self.ZeroMap #accumulated evaporation for water balance self.sumrunoff=self.ZeroMap #accumulated runoff for water balance self.sumint=self.ZeroMap #accumulated interception for water balance self.sumleakage=self.ZeroMap self.CumReinfilt=self.ZeroMap self.sumoutflow=self.ZeroMap self.sumsnowmelt=self.ZeroMap self.CumRad=self.ZeroMap self.SnowMelt=self.ZeroMap self.CumPrec=self.ZeroMap self.CumInwaterMM=self.ZeroMap self.CumInfiltExcess=self.ZeroMap self.CumExfiltWater=self.ZeroMap self.CumSurfaceWater=self.ZeroMap self.CumEvap=self.ZeroMap self.CumPotenEvap=self.ZeroMap self.CumInt=self.ZeroMap self.CumRad=self.ZeroMap self.CumLeakage=self.ZeroMap self.CumPrecPol=self.ZeroMap self.FirstZoneFlux=self.ZeroMap self.FreeWaterDepth=self.ZeroMap self.SumCellWatBal=self.ZeroMap self.PathInfiltExceeded=self.ZeroMap self.SoilInfiltExceeded=self.ZeroMap self.CumOutFlow=self.ZeroMap self.CumCellInFlow=self.ZeroMap self.CumIF=self.ZeroMap self.CumSeepage=self.ZeroMap self.CumActInfilt=self.ZeroMap self.Aspect=scalar(aspect(self.Altitude))# aspect [deg] self.Aspect = ifthenelse(self.Aspect <= 0.0 , scalar(0.001),self.Aspect) # On Flat areas the Aspect function fails, fill in with average... self.Aspect = ifthenelse (defined(self.Aspect), self.Aspect, areaaverage(self.Aspect,self.TopoId)) # Set DCL to riverlength if that is longer that the basic length calculated from grid drainlength = detdrainlength(self.TopoLdd,self.xl,self.yl) self.DCL=max(drainlength,self.RiverLength) # m # Multiply with Factor (taken from upscaling operation, defaults to 1.0 if no map is supplied self.DCL = self.DCL * max(1.0,self.RiverLengthFac) # water depth (m) # set width for kinematic wave to cell width for all cells self.Bw=detdrainwidth(self.TopoLdd,self.xl,self.yl) # However, in the main river we have real flow so set the width to the # width of the river self.Bw=ifthenelse(self.River, RiverWidth, self.Bw) # Add rivers to the WaterFrac, but check with waterfrac map self.RiverFrac = min(1.0,ifthenelse(self.River,(RiverWidth*self.DCL)/(self.xl*self.yl),0)) self.WaterFrac = self.WaterFrac - ifthenelse((self.RiverFrac + self.WaterFrac) > 1.0, self.RiverFrac + self.WaterFrac - 1.0, 0.0) # term for Alpha self.AlpTerm=pow((self.N/(sqrt(self.Slope))),self.Beta) # power for Alpha self.AlpPow=(2.0/3.0)*self.Beta # initial approximation for Alpha #self.initstorage=areaaverage(self.FirstZoneDepth,self.TopoId)+areaaverage(self.UStoreDepth,self.TopoId)#+areaaverage(self.Snow,self.TopoId) # calculate catchmentsize self.upsize=catchmenttotal(self.xl * self.yl,self.TopoLdd) self.csize=areamaximum(self.upsize,self.TopoId) # Save some summary maps self.logger.info("Saving summary maps...") if self.modelSnow: report(self.Cfmax,self.Dir + "/" + self.runId + "/outsum/Cfmax.map") report(self.TTI,self.Dir + "/" + self.runId + "/outsum/TTI.map") report(self.TT,self.Dir + "/" + self.runId + "/outsum/TT.map") report(self.WHC,self.Dir + "/" + self.runId + "/outsum/WHC.map") report(self.Cmax,self.Dir + "/" + self.runId + "/outsum/Cmax.map") report(self.csize,self.Dir + "/" + self.runId + "/outsum/CatchmentSize.map") report(self.upsize,self.Dir + "/" + self.runId + "/outsum/UpstreamSize.map") report(self.EoverR,self.Dir + "/" + self.runId + "/outsum/EoverR.map") report(self.RootingDepth,self.Dir + "/" + self.runId + "/outsum/RootingDepth.map") report(self.CanopyGapFraction,self.Dir + "/" + self.runId + "/outsum/CanopyGapFraction.map") report(self.InfiltCapSoil,self.Dir + "/" + self.runId + "/outsum/InfiltCapSoil.map") report(self.InfiltCapPath,self.Dir + "/" + self.runId + "/outsum/InfiltCapPath.map") report(self.PathFrac,self.Dir + "/" + self.runId + "/outsum/PathFrac.map") report(self.thetaR,self.Dir + "/" + self.runId + "/outsum/thetaR.map") report(self.thetaS,self.Dir + "/" + self.runId + "/outsum/thetaS.map") report(self.FirstZoneMinCapacity,self.Dir + "/" + self.runId + "/outsum/FirstZoneMinCapacity.map") report(self.FirstZoneKsatVer,self.Dir + "/" + self.runId + "/outsum/FirstZoneKsatVer.map") report(self.M,self.Dir + "/" + self.runId + "/outsum/M.map") report(self.FirstZoneCapacity,self.Dir + "/" + self.runId + "/outsum/FirstZoneCapacity.map") report(Terrain_angle,self.Dir + "/" + self.runId + "/outsum/angle.map") report(self.Slope,self.Dir + "/" + self.runId + "/outsum/slope.map") report(WI,self.Dir + "/" + self.runId + "/outsum/WI.map") report(self.CC,self.Dir + "/" + self.runId + "/outsum/CC.map") report(self.N,self.Dir + "/" + self.runId + "/outsum/N.map") report(self.RiverFrac,self.Dir + "/" + self.runId + "/outsum/RiverFrac.map") report(self.xl,self.Dir + "/" + self.runId + "/outsum/xl.map") report(self.yl,self.Dir + "/" + self.runId + "/outsum/yl.map") report(self.reallength,self.Dir + "/" + self.runId + "/outsum/rl.map") report(self.DCL,self.Dir + "/" + self.runId + "/outsum/DCL.map") report(self.Bw,self.Dir + "/" + self.runId + "/outsum/Bw.map") report(ifthen(self.River,self.Bw),self.Dir + "/" + self.runId + "/outsum/RiverWidth.map") if self.updating: report(self.DistToUpdPt,self.Dir + "/" + self.runId + "/outsum/DistToUpdPt.map") self.SaveDir = self.Dir + "/" + self.runId + "/" self.logger.info("Starting Dynamic run...") def resume(self): if self.reinit == 1: self.logger.info("Setting initial conditions to default") self.FirstZoneDepth = self.FirstZoneCapacity * 0.85 self.UStoreDepth = self.FirstZoneCapacity * 0.0 self.WaterLevel = self.ZeroMap self.SurfaceRunoff = self.ZeroMap self.Snow = self.ZeroMap self.SnowWater = self.ZeroMap self.TSoil = self.ZeroMap + 10.0 self.CanopyStorage = self.ZeroMap else: self.logger.info("Setting initial conditions from state files") self.wf_resume(self.Dir + "/instate/") P=self.Bw+(2.0*self.WaterLevel) self.Alpha=self.AlpTerm*pow(P,self.AlpPow) self.OldSurfaceRunoff = self.SurfaceRunoff self.SurfaceRunoffMM=self.SurfaceRunoff * self.QMMConv # Determine initial kinematic wave volume self.KinWaveVolume = self.WaterLevel * self.Bw * self.DCL self.OldKinWaveVolume = self.KinWaveVolume self.SurfaceRunoffMM=self.SurfaceRunoff * self.QMMConv self.InitialStorage = self.FirstZoneDepth + self.UStoreDepth self.CellStorage = self.FirstZoneDepth + self.UStoreDepth # Determine actual water depth self.zi = max(0.0,self.FirstZoneThickness - self.FirstZoneDepth/(self.thetaS -self.thetaR)) # TOPOG_cqf type soil stuff self.f = (self.thetaS -self.thetaR)/self.M def dynamic(self): """ Stuf that is done for each timestep Below a list of variables that can be save to disk as maps or as timeseries (see ini file for syntax): :var self.SurfaceRunoff: Surface runoff in the kinematic wave [m^3/s] :var self.ActEvap: Actual EvapoTranspiration [mm] :var self.WaterLevel: Water level in the kinematic wave [m] (above the bottom) :var self.ActInfilt: Actual infiltration into the unsaturated zone [mm] :var self.CanopyStorage: actual canopystorage (only for subdaily timesteps) [mm] :var self.FirstZoneDepth: Amount of water in the saturated store [mm] :var self.UStoreDepth: Amount of water in the unsaturated store [mm] :var self.TSoil: Top soil temperature [oC] :var self.FirstZoneDepth: amount of available water in the saturated part of the soil [mm] :var self.UStoreDepth: amount of available water in the unsaturated zone [mm] :var self.Transfer: downward flux from unsaturated to saturated zone [mm] :var self.CapFlux: capilary flux from saturated to unsaturated zone [mm] Static variables: :var self.Altitude: The altitude of each cell [m] :var self.Bw: Width of the river [m] :var self.River: booolean map indicating the presence of a river [-] :var self.DLC: length of the river within a cell [m] :var self.ToCubic: Mutiplier to convert mm to m^3/s for fluxes """ self.logger.debug("Step: "+str(int(self.thestep + self._d_firstTimeStep))+"/"+str(int(self._d_nrTimeSteps))) self.thestep = self.thestep + 1 self.Precipitation=cover(self.wf_readmap(self.P_mapstack,0.0),0) self.HP=cover(self.wf_readmap(self.HP_mapstack,0.0),0) #self.PotenEvap=cover(self.wf_readmap(self.PET_mapstack,0.0),0) self.Radiation=cover(self.wf_readmap(self.RAD_mapstack,0.0),0) #Inflow=cover(self.readmap(self.Inflow),0) self.Inflow=pcrut.readmapSave(self.Inflow_mapstack,0.0) self.Seepage=pcrut.readmapSave(self.Seepage_mapstack,0.0) #Inflow=spatial(scalar(0.0)) self.Temperature=cover(self.wf_readmap(self.TEMP_mapstack,0.0),0) self.RH=cover(self.wf_readmap(self.RH_mapstack,0.0),0) self.WindSpeed=cover(self.wf_readmap(self.WindSpeed_mapstack,0.0),0) self.WaterCatch=cover(self.wf_readmap(self.WaterCatch_mapstack,0.0),0) for k, v in multdynapars.iteritems(): estr = k + "=" + k + "*" + str(v) self.logger.debug("Dynamic Parameter multiplication: " + estr) exec estr #PotEvap = self.PotenEvap # self.CanopyStorage = self.CanopyStorage + self.WaterCatch ShortWave = self.Radiation # Change this to measure later!!!! ########################################################################## # Calculate Penman Montieth evaporation ################################## ########################################################################## # Estimate Soil Heat Flux from radiation and leaf area index G = self.Radiation * (1- self.Albedo) * 1.0/(self.LeafAreaIndex+12) # Determine Esat and Delta according to Calder 1990 Esat = 6.1078*exp(17.2694*self.Temperature/(self.Temperature+237.3)) Delta = Esat*17.2694*237.3/sqr(self.Temperature+237.3) # Determine Eact using relative humidity Eact = self.RH * Esat / 100 # Determine specific heat of air Lambda = 4185.5 * (751.78 - (0.5655 * (self.Temperature + 273.15))) # Now determine Gamma p = 900.0 # pressure in mb cp = 1005.0 # J/(kgK) Gamma = (cp * p)/(0.622 * Lambda) # mbC # density of dry air in kg/m^3 rho = 1.201 * (290 * ( p - 0.378 * Eact)/(1000 * (self.Temperature + 273.15))); #rho = hml_rho(p,Eact,Temperature); # At present set A to Radiation * (1- Albedo) and split according to # wetted part of the Canopy WetPart = min(1.0,self.CanopyStorage/self.Cmax); Atrans = (self.Radiation -G) * (1 - self.Albedo) * (1 - WetPart); A = (self.Radiation - G) * (1 - self.Albedo); Acanopy = (self.Radiation - G) * (1 - self.Albedo); Aint = Acanopy - Atrans; # Potential in mm temp var, needed for check AtransMM = Atrans/(Delta+Gamma)/Lambda; AintMM = Aint/(Delta+Gamma)/Lambda; # Determine Ra as a function of Windspeed and canopy parameters # Calculates ra (aerodynamic resistance) according to Arnouds function # CQ Specific function! z = self.WindSpeedHeigth; Zom = self.VegetationHeigth * 0.123; Zoh = 0.25 * Zom; d = 0.66 * self.VegetationHeigth; Ra = 4.72 * ln((z-d)/Zom) * ln((z-d)/Zoh)/(1 + 0.54 * self.WindSpeed); # Now the actual formula, this is for Interception, rs is zero VPD = (Esat - Eact); n = Delta + Gamma; # for interception Rs = 0; Rs = self.ZeroMap tmp = Rs/Ra nn = Delta + (Gamma * (1 + tmp)) t = (Delta * A) + (rho * cp * VPD / Ra) #t = (Delta * Aint) + (rho * cp * VPD / Ra); EA = (t / n) PotEvap = EA/Lambda * self.timestepsecs # now in mm # Now the actual formula, this is for Transpiration # Determine Rs seperate for Pasture and Forest (Hard coded, should be paramiterized in files later) # Reference equations InVPD = ifthenelse(VPD < 0.01, 0.01, VPD) InWave = ifthenelse(ShortWave < 5.0, 5.0, ShortWave) # CQ Specific function! PasRs=exp(1.05 * ln(InVPD) - 0.651 * ln(InWave) + 5.89) # CQ Specific function! ForRs=exp(0.867 * ln(InVPD) - 0.000831 * InWave + 2.81) Rs = max(0.5,min(1000,ifthenelse (scalar(self.LandUse) > 1.0, ForRs,PasRs))) # No transpiration at nigth, this is of no use to the trees. Rs = ifthenelse (InWave < 10.0, 500, Rs) ########################################################################## # Interception according to a modified Rutter model with hourly timesteps# ########################################################################## p = self.CanopyGapFraction; pt = 0.1 * p # Amount of P that falls on the canopy Pfrac = (1 - p -pt) * self.Precipitation # S cannot be larger than Cmax, no gravity drainage bolow that DD = ifthenelse (self.CanopyStorage > self.Cmax , self.Cmax - self.CanopyStorage , 0.0) self.CanopyStorage = self.CanopyStorage - DD # Add the precipitation that falls on the canopy to the store self.CanopyStorage = self.CanopyStorage + Pfrac # Now do the Evap, make sure the store does not get negative dC = -1 * min(self.CanopyStorage, PotEvap) self.CanopyStorage = self.CanopyStorage + dC LeftOver = PotEvap +dC; # Amount of evap not used # Now drain the canopy storage again if needed... D = ifthenelse (self.CanopyStorage > self.Cmax , self.CanopyStorage - self.Cmax , 0.0) self.CanopyStorage = self.CanopyStorage - D # Calculate throughfall ThroughFall = DD + D + p * self.Precipitation StemFlow = self.Precipitation * pt # Calculate interception, this is NET Interception NetInterception = self.Precipitation - ThroughFall - StemFlow Interception = -dC # Determine Evnergy left over for transpiration Atrans = ifthenelse (self.CanopyStorage > 0.001 , ifthenelse(PotEvap > 0 , LeftOver/PotEvap * A , A) , A) t = (Delta * Atrans) + (rho * cp * VPD/ Ra) #t = (Delta * A) + (rho * cp * VPD/ Ra); tmp = Rs/Ra n = Delta + (Gamma * (1 + tmp)) EA = (t / n) PotTrans = EA/Lambda * self.timestepsecs # now in mm RestPotEvap= PotTrans # TODOL bring timeseries export also to the framework # sample timeseries # Do runoff always #self.runTss.sample(Runoff) #self.levTss.sample(self.WaterLevel) ########################################################################## # Start with the soil calculations ###################################### ########################################################################## self.ExfiltWater=self.ZeroMap FreeWaterDepth=self.ZeroMap ########################################################################## # Determine infiltration into Unsaturated store...######################## ########################################################################## # Add precipitation surplus FreeWater storage... FreeWaterDepth= ThroughFall + StemFlow UStoreCapacity = self.FirstZoneCapacity - self.FirstZoneDepth - self.UStoreDepth # Runoff onto water boddies and river network self.RunoffOpenWater = self.RiverFrac * self.WaterFrac * FreeWaterDepth #self.RunoffOpenWater = self.ZeroMap FreeWaterDepth = FreeWaterDepth - self.RunoffOpenWater if self.RunoffGenSigmaFunction: self.AbsoluteGW=self.DemMax-(self.zi*self.GWScale) self.SubCellFrac = sCurve(self.AbsoluteGW,c=self.CC,a=self.Altitude+1.0) self.SubCellRunoff = self.SubCellFrac * FreeWaterDepth self.SubCellGWRunoff = min(self.SubCellFrac * self.FirstZoneDepth, self.SubCellFrac * self.Slope * self.FirstZoneKsatVer * exp(-self.f * self.zi) ) self.FirstZoneDepth=self.FirstZoneDepth-self.SubCellGWRunoff FreeWaterDepth = FreeWaterDepth - self.SubCellRunoff else: self.AbsoluteGW=self.DemMax-(self.zi*self.GWScale) self.SubCellFrac = spatial(scalar(0.0)) self.SubCellGWRunoff = spatial(scalar(0.0)) self.SubCellRunoff = spatial(scalar(0.0)) #----->> # First determine if the soil infiltration capacity can deal with the # amount of water # split between infiltration in undisturbed soil and compacted areas (paths) SoilInf = FreeWaterDepth * (1- self.PathFrac) PathInf = FreeWaterDepth * self.PathFrac if self.modelSnow: soilInfRedu = ifthenelse(self.TSoil < 0.0 , self.cf_soil, 1.0) else: soilInfRedu = 1.0 MaxInfiltSoil= min(self.InfiltCapSoil*soilInfRedu,SoilInf) self.SoilInfiltExceeded=self.SoilInfiltExceeded + scalar(self.InfiltCapSoil*soilInfRedu < SoilInf) InfiltSoil = min(MaxInfiltSoil, UStoreCapacity) self.UStoreDepth = self.UStoreDepth + InfiltSoil UStoreCapacity = UStoreCapacity - InfiltSoil FreeWaterDepth = FreeWaterDepth - InfiltSoil # <------- MaxInfiltPath= min(self.InfiltCapPath*soilInfRedu,PathInf) #self.PathInfiltExceeded=self.PathInfiltExceeded + ifthenelse(self.InfiltCapPath < FreeWaterDepth, scalar(1), scalar(0)) self.PathInfiltExceeded=self.PathInfiltExceeded + scalar(self.InfiltCapPath*soilInfRedu < PathInf) InfiltPath = min(MaxInfiltPath, UStoreCapacity) self.UStoreDepth = self.UStoreDepth + InfiltPath UStoreCapacity = UStoreCapacity - InfiltPath FreeWaterDepth = FreeWaterDepth - InfiltPath self.ActInfilt = InfiltPath + InfiltSoil self.InfiltExcess = ifthenelse (UStoreCapacity > 0.0, FreeWaterDepth, 0.0) self.CumInfiltExcess=self.CumInfiltExcess+self.InfiltExcess self.ActEvap, self.FirstZoneDepth, self.UStoreDepth, self.ActEvapUStore = actEvap_SBM(self.RootingDepth,self.zi,self.UStoreDepth,self.FirstZoneDepth, PotTrans,self.rootdistpar) #self.ActEvap = self.ZeroMap #self.ActEvapUStore = self.ZeroMap ########################################################################## # Transfer of water from unsaturated to saturated store...################ ########################################################################## self.zi = max(0.0,self.FirstZoneThickness - self.FirstZoneDepth/(self.thetaS -self.thetaR)) # Determine actual water depth Ksat = self.FirstZoneKsatVer * exp(-self.f * self.zi) self.DeepKsat = self.FirstZoneKsatVer * exp(-self.f * self.FirstZoneThickness) # Determine saturation deficit. NB, as noted by Vertessy and Elsenbeer 1997 # this deficit does NOT take into account the water in the unsaturated zone SaturationDeficit = self.FirstZoneCapacity - self.FirstZoneDepth # now the actual tranfer to the saturated store.. self.Transfer = min(self.UStoreDepth,ifthenelse (SaturationDeficit <= 0.00001, 0.0, Ksat * self.UStoreDepth/(SaturationDeficit+1))) # Determine Ksat at base #DeepTransfer = min(self.UStoreDepth,ifthenelse (SaturationDeficit <= 0.00001, 0.0, DeepKsat * self.UStoreDepth/(SaturationDeficit+1))) # Now add leakage # Limit to MaxLeakage/day. Leakage percentage gets bigger if the # storm is bigger (macropores start kicking in... ActLeakage = max(0,min(self.MaxLeakage,self.Transfer * exp(0.01 * self.Transfer)/e)); self.Transfer = self.Transfer - ActLeakage; # Now add leakage. to deeper groundwater #ActLeakage = cover(max(0,min(self.MaxLeakage* timestepsecs/basetimestep,ActLeakage)),0) # Now look if there is Seeapage #ActLeakage = ifthenelse(self.Seepage > 0.0, -1.0 * Seepage, ActLeakage) self.FirstZoneDepth = self.FirstZoneDepth + self.Transfer - ActLeakage self.UStoreDepth = self.UStoreDepth - self.Transfer # Determine % saturated #Sat = ifthenelse(self.FirstZoneDepth >= (self.FirstZoneCapacity*0.999), scalar(1.0), scalar(0.0)) self.Sat = max(self.SubCellFrac,scalar(self.FirstZoneDepth >= (self.FirstZoneCapacity*0.999))) #PercSat = areaaverage(scalar(Sat),self.TopoId) * 100 ########################################################################## # Horizontal (downstream) transport of water ############################# ########################################################################## if self.waterdem: waterDem = self.Altitude - (self.zi * 0.001) waterLdd = lddcreate(waterDem,1E35,1E35,1E35,1E35) #waterLdd = lddcreate(waterDem,1,1,1,1) waterSlope=max(0.00001,slope(waterDem)*celllength()/self.reallength) self.zi = max(0.0,self.FirstZoneThickness - self.FirstZoneDepth/(self.thetaS -self.thetaR)) # Determine actual water depth if self.waterdem: MaxHor = max(0.0,min(self.FirstZoneKsatVer * waterSlope * exp(-SaturationDeficit/self.M),self.FirstZoneDepth)) self.FirstZoneFlux = accucapacityflux (waterLdd, self.FirstZoneDepth, MaxHor) self.FirstZoneDepth = accucapacitystate (waterLdd, self.FirstZoneDepth, MaxHor) else: # #MaxHor = max(0,min(self.FirstZoneKsatVer * self.Slope * exp(-SaturationDeficit/self.M),self.FirstZoneDepth*(self.thetaS-self.thetaR))) * timestepsecs/basetimestep MaxHor = max(0.0,min(self.FirstZoneKsatVer * self.Slope * exp(-SaturationDeficit/self.M),self.FirstZoneDepth)) self.FirstZoneFlux = accucapacityflux (self.TopoLdd, self.FirstZoneDepth, MaxHor) self.FirstZoneDepth = accucapacitystate (self.TopoLdd, self.FirstZoneDepth, MaxHor) ########################################################################## # Determine returnflow from first zone ########################## ########################################################################## self.ExfiltWaterFrac = sCurve(self.FirstZoneDepth,a=self.FirstZoneCapacity,c=5.0) self.ExfiltWater=self.ExfiltWaterFrac * (self.FirstZoneDepth - self.FirstZoneCapacity) #self.ExfiltWater=ifthenelse (self.FirstZoneDepth - self.FirstZoneCapacity > 0 , self.FirstZoneDepth - self.FirstZoneCapacity , 0.0) self.FirstZoneDepth=self.FirstZoneDepth - self.ExfiltWater # Re-determine UStoreCapacity UStoreCapacity = self.FirstZoneCapacity - self.FirstZoneDepth - self.UStoreDepth #Determine capilary rise self.zi = max(0.0,self.FirstZoneThickness - self.FirstZoneDepth/(self.thetaS -self.thetaR)) # Determine actual water depth Ksat = self.FirstZoneKsatVer * exp(-self.f * self.zi) MaxCapFlux = max(0.0,min(Ksat,self.ActEvapUStore,UStoreCapacity,self.FirstZoneDepth)) # No capilary flux is roots are in water, max flux if very near to water, lower flux if distance is large CapFluxScale = ifthenelse(self.zi > self.RootingDepth, self.CapScale/(self.CapScale + self.zi -self.RootingDepth), 0.0) self.CapFlux = MaxCapFlux * CapFluxScale self.UStoreDepth = self.UStoreDepth + self.CapFlux self.FirstZoneDepth = self.FirstZoneDepth - self.CapFlux # org SurfaceWater = self.SurfaceRunoff * self.DCL * self.QMMConv # SurfaceWater (mm) from SurfaceRunoff (m3/s) SurfaceWater = self.SurfaceRunoff * self.QMMConv # SurfaceWater (mm) from SurfaceRunoff (m3/s) self.CumSurfaceWater = self.CumSurfaceWater + SurfaceWater # Estimate water that may re-infiltrate if self.reInfilt: Reinfilt = max(0,min(SurfaceWater,min(self.InfiltCapSoil,UStoreCapacity))) self.CumReinfilt=self.CumReinfilt + Reinfilt self.UStoreDepth = self.UStoreDepth + Reinfilt else: Reinfilt = self.ZeroMap self.InwaterMM=max(0.0,self.ExfiltWater + FreeWaterDepth + self.SubCellRunoff + self.SubCellGWRunoff + self.RunoffOpenWater - Reinfilt) self.Inwater=self.InwaterMM * self.ToCubic # m3/s self.ExfiltWaterCubic=self.ExfiltWater * self.ToCubic self.SubCellGWRunoffCubic = self.SubCellGWRunoff * self.ToCubic self.SubCellRunoffCubic = self.SubCellRunoff * self.ToCubic self.InfiltExcessCubic = self.InfiltExcess * self.ToCubic self.FreeWaterDepthCubic=FreeWaterDepth * self.ToCubic self.ReinfiltCubic=-1.0 * Reinfilt * self.ToCubic self.Inwater=self.Inwater + self.Inflow # Add abstractions/inflows in m^3/sec ########################################################################## # Runoff calculation via Kinematic wave ################################## ########################################################################## # per distance along stream q=self.Inwater/self.DCL # discharge (m3/s) self.SurfaceRunoff = kinematic(self.TopoLdd, self.SurfaceRunoff,q,self.Alpha, self.Beta,self.Tslice,self.timestepsecs,self.DCL) # m3/s self.SurfaceRunoffMM=self.SurfaceRunoff*self.QMMConv # SurfaceRunoffMM (mm) from SurfaceRunoff (m3/s) self.updateRunOff() self.InflowKinWaveCell=upstream(self.TopoLdd,self.SurfaceRunoff) self.MassBalKinWave = (self.KinWaveVolume - self.OldKinWaveVolume)/self.timestepsecs + self.InflowKinWaveCell + self.Inwater - self.SurfaceRunoff Runoff=self.SurfaceRunoff # Updating # -------- # Assume a tss file with as many columns as outpulocs. Start updating for each non-missing value and start with the # first column (nr 1). Assumes that outputloc and columns match! if self.updating: QM = timeinputscalar(updateFile, self.UpdateMap) * self.QMMConv # Now update the state. Just add to the Ustore # self.UStoreDepth = result # No determine multiplication ratio for each gauge influence area. # For missing gauges 1.0 is assumed (no change). # UpDiff = areamaximum(QM, self.UpdateMap) - areamaximum(self.SurfaceRunoffMM, self.UpdateMap) UpRatio = areamaximum(QM, self.UpdateMap)/areamaximum(self.SurfaceRunoffMM, self.UpdateMap) UpRatio = cover(areaaverage(UpRatio,self.TopoId),1.0) # Now split between Soil and Kyn wave UpRatioKyn = min(MaxUpdMult,max(MinUpdMult,(UpRatio - 1.0) * UpFrac + 1.0)) UpRatioSoil = min(MaxUpdMult,max(MinUpdMult,(UpRatio - 1.0) * (1.0 - UpFrac) + 1.0)) # update/nudge self.UStoreDepth for the whole upstream area, # not sure how much this helps or worsens things if UpdSoil: toadd = min((self.UStoreDepth * UpRatioSoil) - self.UStoreDepth,StorageDeficit * 0.95) self.UStoreDepth = self.UStoreDepth + toadd # Update the kinematic wave reservoir up to a maximum upstream distance # TODO: add (much smaller) downstream updating also? MM = (1.0 - UpRatioKyn)/self.UpdMaxDist UpRatioKyn = MM * self.DistToUpdPt + UpRatioKyn self.SurfaceRunoff = self.SurfaceRunoff * UpRatioKyn self.SurfaceRunoffMM=self.SurfaceRunoff*self.QMMConv # SurfaceRunoffMM (mm) from SurfaceRunoff (m3/s) self.updateRunOff() Runoff=self.SurfaceRunoff ########################################################################## # water balance ########################################### ########################################################################## # Single cell based water budget CellStorage = self.UStoreDepth+self.FirstZoneDepth DeltaStorage = CellStorage - self.InitialStorage OutFlow = self.FirstZoneFlux CellInFlow = upstream(self.TopoLdd,scalar(self.FirstZoneFlux)); #CellWatBal = ActInfilt - self.ActEvap - self.ExfiltWater - ActLeakage + Reinfilt + IF - OutFlow + (OldCellStorage - CellStorage) #SumCellWatBal = SumCellWatBal + CellWatBal; self.CumOutFlow = self.CumOutFlow + OutFlow self.CumActInfilt = self.CumActInfilt + self.ActInfilt self.CumCellInFlow = self.CumCellInFlow + CellInFlow self.CumPrec=self.CumPrec+self.Precipitation self.CumEvap=self.CumEvap+self.ActEvap self.CumPotenEvap=self.CumPotenEvap+PotTrans self.CumInt=self.CumInt+Interception self.CumLeakage=self.CumLeakage+ActLeakage self.CumInwaterMM=self.CumInwaterMM+self.InwaterMM self.CumExfiltWater=self.CumExfiltWater+self.ExfiltWater # Water budget #self.watbal = self.CumPrec- self.CumEvap - self.CumInt - self.CumInwaterMM - DeltaStorage - self.CumOutFlow + self.CumIF #self.watbal = self.CumActInfilt - self.CumEvap - self.CumExfiltWater - DeltaStorage - self.CumOutFlow + self.CumIF self.watbal = self.CumPrec + self.CumCellInFlow - self.CumOutFlow- self.CumEvap - self.CumLeakage - self.CumInwaterMM - self.CumInt - DeltaStorage + self.CumReinfilt def main(): """ Perform command line execution of the model. """ caseName = "default_cqf" global multpars runId = "run_default" configfile="wflow_cqf.ini" _lastTimeStep = 0 _firstTimeStep = 1 fewsrun=False runinfoFile="runinfo.xml" timestepsecs=86400 wflow_cloneMap = 'wflow_subcatch.map' NoOverWrite=1 ## Main model starts here ######################################################################## try: opts, args = getopt.getopt(sys.argv[1:], 'XF:L:hC:Ii:v:S:T:WNR:u:s:EP:p:Xx:U:fOc:') except getopt.error, msg: pcrut.usage(msg) for o, a in opts: if o == '-F': runinfoFile = a fewsrun = True if o == '-P': exec ("multpars =" + a,globals(), globals()) if o == '-p': exec "multdynapars =" + a exec ("multdynapars =" + a,globals(), globals()) if o == '-C': caseName = a if o == '-R': runId = a if o == '-c': configfile = a if o == '-s': timestepsecs = int(a) if o == '-T': _lastTimeStep=int(a) if o == '-S': _firstTimeStep=int(a) if o == '-h': usage() if o == '-f': NoOverWrite = 0 if fewsrun: ts = getTimeStepsfromRuninfo(runinfoFile) if (ts): _lastTimeStep = ts * 86400/timestepsecs _firstTimeStep = 1 else: print "Failed to get timesteps from runinfo file: " + runinfoFile exit(2) if _lastTimeStep < _firstTimeStep: usage() myModel = WflowModel(wflow_cloneMap, caseName,runId,configfile) dynModelFw = wf_DynamicFramework(myModel, _lastTimeStep,firstTimestep=_firstTimeStep) dynModelFw.createRunId(NoOverWrite=NoOverWrite) for o, a in opts: if o == '-X': configset(myModel.config,'model','OverWriteInit','1',overwrite=True) if o == '-I': configset(myModel.config,'model','reinit','1',overwrite=True) if o == '-i': configset(myModel.config,'model','intbl',a,overwrite=True) if o == '-s': configset(myModel.config,'model','timestepsecs',a,overwrite=True) if o == '-x': configset(myModel.config,'model','sCatch',a,overwrite=True) if o == '-c': configset(myModel.config,'model','configfile', a,overwrite=True) if o == '-M': configset(myModel.config,'model','MassWasting',"1",overwrite=True) if o == '-N': configset(myModel.config,'model','nolateral','1',overwrite=True) if o == '-Q': configset(myModel.config,'model','ExternalQbase','1',overwrite=True) if o == '-U': configset(myModel.config,'model','updateFile',a,overwrite=True) configset(myModel.config,'model','updating',"1",overwrite=True) if o == '-u': print a exec "updateCols =" + a if o == '-E': configset(myModel.config,'model','reInfilt','1',overwrite=True) if o == '-R': runId = a if o == '-W': configset(myModel.config,'model','waterdem','1',overwrite=True) dynModelFw._runInitial() dynModelFw._runResume() dynModelFw._runDynamic(_firstTimeStep,_lastTimeStep) dynModelFw._runSuspend() fp = open(caseName + "/" + runId + "/runinfo/configofrun.ini",'wb') myModel.config.write(fp ) if __name__ == "__main__": main()