Index: wflow-py/UnitTests/TestFrameWork.py =================================================================== diff -u -r652aea29305298e160136ae309f1f80c94569bec -rc7476868c506252d099ab57c9909e317b43aeeee --- wflow-py/UnitTests/TestFrameWork.py (.../TestFrameWork.py) (revision 652aea29305298e160136ae309f1f80c94569bec) +++ wflow-py/UnitTests/TestFrameWork.py (.../TestFrameWork.py) (revision c7476868c506252d099ab57c9909e317b43aeeee) @@ -37,7 +37,7 @@ dynModelFw._runSuspend() # saves the state variables dynModelFw._wf_shutdown() - # nore read the csv results acn check of they match the first run + # Now read the csv results acn check of they match the first run # Sum should be approx c 4.569673676 my_data = wf.genfromtxt(os.path.join(caseName,runId,"tes.csv"), delimiter=',') Index: wflow-py/wflow/pcrut.py =================================================================== diff -u -r3d8b5fe212103dd76367b768bb968449e5022d39 -rc7476868c506252d099ab57c9909e317b43aeeee --- wflow-py/wflow/pcrut.py (.../pcrut.py) (revision 3d8b5fe212103dd76367b768bb968449e5022d39) +++ wflow-py/wflow/pcrut.py (.../pcrut.py) (revision c7476868c506252d099ab57c9909e317b43aeeee) @@ -253,107 +253,3 @@ return rmat - - -# def correctrad(Day,Hour,Lat,Lon,Slope,Aspect,Altitude): -# """ Determines radiation over a DEM assuming clear sky""" -# -# Sc = 1367.0 # Solar constant (Gates, 1980) [W/m2] -# Trans = 0.6 # Transmissivity tau (Gates, 1980) -# pi = 3.1416 -# AtmPcor = pow(((288.0-0.0065*Altitude)/288.0),5.256) -# Lat = Lat * pi/180 -# ########################################################################## -# # Calculate Solar Angle and correct radiation ############################ -# ########################################################################## -# # Solar geometry -# # ---------------------------- -# # SolDec :declination sun per day between +23 & -23 [deg] -# # HourAng :hour angle [-] of sun during day -# # SolAlt :solar altitude [deg], height of sun above horizon -# # SolDec = -23.4*cos(360*(Day+10)/365); -# # Now added a new function that should work on all latitudes! -# theta =(Day-1)*2 * pi/365 # day expressed in radians -# -# SolDec =0.006918-0.399912 * cos(theta)+0.070257 * sin(theta) - 0.006758 * cos(2*theta)+0.000907 * sin(2*theta) - 0.002697 * cos(3*theta)+0.001480 * sin(3*theta) -# -# HourAng = 180/pi * 15*(Hour-12.01) -# SolAlt = scalar(asin(scalar(sin(Lat)*sin(SolDec)+cos(Lat)*cos(SolDec)*cos(HourAng)))) -# -# # Solar azimuth -# # ---------------------------- -# # SolAzi :angle solar beams to N-S axes earth [deg] -# SolAzi = scalar(acos((sin(SolDec)*cos(Lat)-cos(SolDec)* sin(Lat)*cos(HourAng))/cos(SolAlt))) -# SolAzi = ifthenelse(Hour <= 12, SolAzi, 360 - SolAzi) -# -# # Surface azimuth -# # ---------------------------- -# # cosIncident :cosine of angle of incident; angle solar beams to angle surface -# cosIncident = sin(SolAlt)*cos(Slope)+cos(SolAlt)*sin(Slope)*cos(SolAzi-Aspect) -# -# -# # Critical angle sun -# # ---------------------------- -# # HoriAng :tan maximum angle over DEM in direction sun, 0 if neg -# # CritSun :tan of maximum angle in direction solar beams -# # Shade :cell in sun 1, in shade 0 -# # NOTE: for a changing DEM in time use following 3 statements and put a # -# # for the 4th CritSun statement -# HoriAng = cover(horizontan(Altitude,directional(SolAzi)),0) -# #HoriAng = horizontan(Altitude,directional(SolAzi)) -# HoriAng = ifthenelse(HoriAng < 0, scalar(0), HoriAng) -# CritSun = ifthenelse(SolAlt > 90, scalar(0), scalar(atan(HoriAng))) -# Shade = SolAlt > CritSun -# -# # Radiation outer atmosphere -# # ---------------------------- -# OpCorr = Trans**((sqrt(1229+(614*sin(SolAlt))**2) -614*sin(SolAlt))*AtmPcor) # correction for air masses [-] -# Sout = Sc*(1+0.034*cos(360*Day/365)) # radiation outer atmosphere [W/m2] -# Snor = Sout*OpCorr # rad on surface normal to the beam [W/m2] -# -# # Radiation at DEM -# # ---------------------------- -# # Sdir :direct sunlight on a horizontal surface [W/m2] if no shade -# # Sdiff :diffuse light [W/m2] for shade and no shade -# # Stot :total incomming light Sdir+Sdiff [W/m2] at Hour -# # Radiation :avg of Stot(Hour) and Stot(Hour-HourStep) -# # NOTE: PradM only valid for HourStep & DayStep = 1 -# SdirCor = ifthenelse(Snor*cosIncident*scalar(Shade)<0,0.0,Snor*cosIncident*scalar(Shade)) -# Sdir = ifthenelse(Snor*cosIncident<0,0.0,Snor*cosIncident) -# Sdiff = ifthenelse(Sout*(0.271-0.294*OpCorr)*sin(SolAlt)<0, 0.0, Sout*(0.271-0.294*OpCorr)*sin(SolAlt)) -# AtmosDiffFrac = ifthenelse(Sdir > 0, Sdiff/Sdir, 1) -# -# # Stot = cover(Sdir+Sdiff,windowaverage(Sdir+Sdiff,3)); # Rad [W/m2] -# Stot = Sdir + Sdiff # Rad [W/m2] -# StotCor = SdirCor + Sdiff # Rad [W/m2] -# -# -# return scalar(SolAlt) -# -# -# def GenRadMaps(SaveDir,SaveName,Lat,Lon,Slope,Aspect,Altitude,debug): -# """ Generates daily radiation maps for a whole year. -# It does so by running correctrad for a whole year with hourly -# steps and averaging this per day.""" -# -# -# for Day in range(1,366): -# avgrad = 0.0 * Altitude -# for Hour in range(6,19): -# print "day: " + str(Day) + " Hour: " + str(Hour) -# crad = correctrad(Day,Hour,Lat,Lon,Slope,Aspect,Altitude) -# if debug: -# nr = "%0.3d" % Hour -# report(crad,SaveDir + "/" + str(Day) + "/" + SaveName + "00000." + nr) -# avgrad=avgrad + crad -# -# nr = "%0.3d" % Day -# report(avgrad/13.0,SaveDir + "/" + SaveName + "00000." + nr) -# - - - - - - - \ No newline at end of file Index: wflow-py/wflow/wflow_hbv.py =================================================================== diff -u -r8e380ffbfaedb9df9fccee3c03ebb53e346d1a71 -rc7476868c506252d099ab57c9909e317b43aeeee --- wflow-py/wflow/wflow_hbv.py (.../wflow_hbv.py) (revision 8e380ffbfaedb9df9fccee3c03ebb53e346d1a71) +++ wflow-py/wflow/wflow_hbv.py (.../wflow_hbv.py) (revision c7476868c506252d099ab57c9909e317b43aeeee) @@ -480,6 +480,9 @@ self.Slope= slope(self.Altitude) self.Slope=ifthen(boolean(self.TopoId),max(0.001,self.Slope*celllength()/self.reallength)) Terrain_angle=scalar(atan(self.Slope)) + temp = catchmenttotal(cover(1.0), self.TopoLdd) * self.reallength * 0.001 * 0.001 * self.reallength + self.QMMConvUp = cover(self.timestepsecs * 0.001)/temp + # Multiply parameters with a factor (for calibration etc) -P option in command line for k, v in multpars.iteritems(): @@ -865,7 +868,9 @@ Runoff=self.SurfaceRunoff - + self.QCatchmentMM = self.SurfaceRunoff * self.QMMConvUp + self.RunoffCoeff = self.QCatchmentMM/catchmenttotal(self.Precipitation, self.TopoLdd)/catchmenttotal(cover(1.0), self.TopoLdd) + self.sumprecip=self.sumprecip + self.Precipitation #accumulated rainfall for water balance self.sumevap=self.sumevap + ActEvap #accumulated evaporation for water balance self.sumpotevap=self.sumpotevap + self.PotEvaporation Index: wflow-py/wflow/wflow_sbm.py =================================================================== diff -u -re2a41c70e2142cd8157fdcc4264039e194e0522d -rc7476868c506252d099ab57c9909e317b43aeeee --- wflow-py/wflow/wflow_sbm.py (.../wflow_sbm.py) (revision e2a41c70e2142cd8157fdcc4264039e194e0522d) +++ wflow-py/wflow/wflow_sbm.py (.../wflow_sbm.py) (revision c7476868c506252d099ab57c9909e317b43aeeee) @@ -723,6 +723,7 @@ #self.QMMConvUp = 1000.0 * self.timestepsecs / ( catchmenttotal(cover(1.0), self.TopoLdd) * self.reallength * self.reallength) #m3/s --> mm over upstreams temp = catchmenttotal(cover(1.0), self.TopoLdd) * self.reallength * 0.001 * 0.001 * self.reallength self.QMMConvUp = cover(self.timestepsecs * 0.001)/temp + self.ToCubic = (self.reallength * self.reallength * 0.001) / self.timestepsecs # m3/s self.KinWaveVolume = self.ZeroMap self.OldKinWaveVolume = self.ZeroMap @@ -793,8 +794,6 @@ 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) @@ -815,17 +814,11 @@ 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 - - """ lst = ['self.RiverWidth', - 'self.Cmax', - 'self.csize', - 'self.upsize', - 'self.EoverR', - 'self.RootingDepth', - 'self.CanopyGapFraction', - 'self.InfiltCapSoil', + 'self.Cmax', 'self.csize', 'self.upsize', + 'self.EoverR', 'self.RootingDepth', + 'self.CanopyGapFraction', 'self.InfiltCapSoil', 'self.InfiltCapPath', 'self.PathFrac', 'self.thetaR', @@ -840,11 +833,10 @@ 'self.N', 'self.RiverFrac', 'self.WaterFrac', - 'self.xl', - 'self.yl', - 'self.reallength', + 'self.xl', 'self.yl', 'self.reallength', 'self.DCL', - 'self.Bw'] + 'self.Bw', + 'self.PathInfiltExceeded','self.SoilInfiltExceeded'] return lst @@ -897,13 +889,17 @@ def dynamic(self): """ - Stuf that is done for each timestep + Stuf that is done for each timestep of the model + ------------------------------------------------ Below a list of variables that can be save to disk as maps or as timeseries (see ini file for syntax): - *Dynamic variables* + Dynamic variables + ~~~~~~~~~~~~~~~~~ + All these can be saved per timestep if needed (see the config file [outputmaps] section). + :var self.Precipitation: Gross precipitation [mm] :var self.Temperature: Air temperature [oC] :var self.PotenEvap: Potential evapotranspiration [mm] @@ -933,9 +929,11 @@ :var self.Transfer: downward flux from unsaturated to saturated zone [mm] :var self.CapFlux: capilary flux from saturated to unsaturated zone [mm] :var self.CanopyStorage: Amount of water on the Canopy [mm] + :var self.RunoffCoeff: Runoff coefficient (Q/P) for each cell taking into accoutn the whole upstream area [-] - *Static variables* + Static variables + ~~~~~~~~~~~~~~~~ :var self.Altitude: The altitude of each cell [m] :var self.Bw: Width of the river [m] @@ -1005,15 +1003,11 @@ self.Interception=NetInterception - ########################################################################## - # Start with the soil calculations ###################################### - ########################################################################## - ########################################################################## - # Determine infiltration into Unsaturated store...######################## - ########################################################################## - # Add precipitation surplus FreeWater storage... + # Start with the soil calculations + # -------------------------------- + self.AvailableForInfiltration = ThroughFall + StemFlow UStoreCapacity = self.FirstZoneCapacity - self.FirstZoneDepth - self.UStoreDepth @@ -1060,7 +1054,6 @@ self.AvailableForInfiltration = self.AvailableForInfiltration - 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 @@ -1253,7 +1246,6 @@ self.MassBalKinWave = (self.KinWaveVolume - self.OldKinWaveVolume) / self.timestepsecs + self.InflowKinWaveCell + self.Inwater - self.SurfaceRunoff Runoff = self.SurfaceRunoff - self.QCatchmentMM = self.SurfaceRunoff * self.QMMConvUp # Updating # -------- @@ -1294,6 +1286,10 @@ # water balance ########################################### ########################################################################## + + self.QCatchmentMM = self.SurfaceRunoff * self.QMMConvUp + self.RunoffCoeff = self.QCatchmentMM/catchmenttotal(self.PrecipitationPlusMelt, self.TopoLdd)/catchmenttotal(cover(1.0), self.TopoLdd) + # Single cell based water budget. snow not included yet. CellStorage = self.CanopyStorage CellStorage = self.UStoreDepth + self.FirstZoneDepth + self.LowerZoneStorage