Index: examples/wflow_orientale_topoflex_distrib/wflow_orientale_topoflex_distrib_accutravel.ini =================================================================== diff -u --- examples/wflow_orientale_topoflex_distrib/wflow_orientale_topoflex_distrib_accutravel.ini (revision 0) +++ examples/wflow_orientale_topoflex_distrib/wflow_orientale_topoflex_distrib_accutravel.ini (revision 08e6dbe2c7a682c7095170bd1bf89d3674fefb7c) @@ -0,0 +1,222 @@ + +[framework] +outputformat = 1 + +# Model parameters and settings +[model] +AnnualDischarge=2290 +# Alpha for wiver-width estimation 5 for mountain stream 60 for the river rhine +Alpha=120 +ModelSnow=1 +ScalarInput=1 +InputSeries=1 +InterpolationMethod=inv +Tslice=1 +UpdMaxDist=300000.0 +P_style = 1 +L_IRURFR = 0 +L_URFR = 0 +L_FR = 0 +maxTransitTime = 9 +DistForcing = 3 +maxGaugeId = 10 +Ks = 0.0004 +spinUp_time = 5 +NSEbreak = 0 + +#W=wetland(0) H=hillslope(1) P=plateau(2) WD = drained wetland(3) +classes = ['W','H'] +timestepsecs = 3600 + +#selection of reservoir configuration - 'None' means reservoir should not be modelled for this class +selectSi = interception_overflow2, interception_overflow2 +selectSu= unsatZone_LP_beta, unsatZone_LP_beta +selectSus=None,None +selectSf=fastRunoff_lag2, fastRunoff_lag2 +selectSr=None,None + +# selection of Ss (lumped over entire catchment, so only one!) +# selectSs = + +#input time series +Pfile_1 = intss\1_P.tss +Efile_1 = intss\1_PET.tss +Tfile_1 = intss\1_T.tss +#Qfile_1 = D:/TEuser/Onderzoek/Promotie/modellen/OpenStreams/wflow/wflow_orientale_topoflex/intss/1_Qobs.tss +#Pfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_rain_10H.tss +#Efile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_Ep_10H.tss +#TDMfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_tempDMean_10H.tss +#RNfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_radNet_10H.tss +#RSfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_radSI_10H.tss +#SGfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_sgamma_10H.tss +#VPDfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_vpd_10H.tss +#Wfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_wind_10H.tss +#DSfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_dayS_10H.tss +#DEfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_dayE_10H.tss +#LAIfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_lai_10H.tss +#rst_lai_0 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_rstMin_laiEff_W4.tss +#rst_lai_1 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_rstMin_laiEff_HPPPA4.tss + +#wflow maps with percentages +wflow_percent_0 = staticmaps/wflow_percentW4.map +wflow_percent_1 = staticmaps/wflow_percentHPPPA4.map + +#constant model parameters - some are catchment dependent +Ks = 0.0004 +lamda = 2.45e6 +Cp = 1.01e-3 +rhoA = 1.29 +rhoW = 1000 +gamma = 0.066 +JC_Topt = 301 + +#parameters for fluxes and storages +sumax = [140, 300] +beta = [0.1, 0.1] +D = [0, 0.10] +Kf = [0.1, 0.006] +Tf = [1, 3] +imax = [1.2, 2] +perc = [0, 0.000] +cap = [0.09, 0] +LP = [0.5, 0.8] +Ce = [1, 1] + +#Jarvis stressfunctions +JC_D05 = [1.5,1.5] +JC_cd1 = [3,3] +JC_cd2 = [0.1,0.1] +JC_cr = [100,100] +JC_cuz = [0.07,0.07] +SuFC = [0.98,0.98] +SuWP = [0.1,0.1] +JC_rstmin = [150,250] + +#parameters not used for this configuration +samax = [0,0] +famax = [0,0] +sumin = [0,0] +susmax1 = [0,0] +susmax2 = [0,0] +susmax3 = [0,0] +srmax = [0,0] +Co = [1,1] +Kd = [0,0] +Kr =[0,0] + + +[layout] +# if set to zero the cell-size is given in lat/long (the default) +sizeinmetres = 0 + +[outputmaps] +self.Qstate=Qstate +self.Qrout=Qrout +self.QLagTot=QLagTot +self.WB_rout=WB_rout +self.dSdt=dSdt +#self.Pe=pe +#self.Ei=Ei +#self.Si=si +#self.Qfin=Qfin +#self.QfinLag=QfinLag +#self.RunoffLaged=runlag +#self.Su=Su +#self.Qu=Qu +#self.evaporation=act +#self.Su_diff=sudiff + + +# List all timeseries in tss format to save in this section. Timeseries are +# produced per subcatchment. + +[outputtss_0] +samplemap=staticmaps/wflow_mgauges.map +#states +self.Si[0]=SiW.tss +self.Si[1]=SiH.tss +self.Sf[1]=SfH.tss +self.Sf[0]=SfW.tss +self.Su[1]=SuH.tss +self.Su[0]=SuW.tss +#fluxen +self.Precipitation=Prec.tss +self.Qu_[0]=QuW.tss +self.Qu_[1]=QuH.tss +self.Ei_[0]=EiW.tss +self.Ei_[1]=EiH.tss +self.Eu_[0]=EuW.tss +self.Eu_[1]=EuH.tss +self.Pe_[0]=peW.tss +self.Pe_[1]=peH.tss +self.Perc_[1]=PercH.tss +self.Cap_[0]=CapW.tss +self.Qf_[1] = QfH.tss +self.Qf_[0] = QfW.tss +self.Qftotal = Qftotal.tss +self.Qtlag = Qtlag.tss + +[outputtss_1] +samplemap = staticmaps/wflow_gauges.map +#states +self.Ss = Ss.tss +#fluxen +self.Qs = Qs.tss +self.QLagTot = runLag.tss +self.Qrout = Qrout.tss +self.WBtot = WB.tss + +#[outputtss] +##stages +#self.Si[0]=SiW +#self.Si[1]=SiH +#self.Sf[1]=SfH +#self.Sf[0]=SfW +#self.Su[1]=SuH +#self.Su[0]=SuW +##fluxen +#self.Precipitation=Prec +#self.Qu_[0]=QuW +#self.Qu_[1]=QuH +#self.Ei_[0]=EiW +#self.Ei_[1]=EiH +#self.Eu_[0]=EuW +#self.Eu_[1]=EuH +#self.Pe_[0]=peW +#self.Pe_[1]=peH +#self.Perc_[1]=PercH +#self.Cap_[0]=CapW +#self.Qf_[1] = QfH +#self.Qf_[0] = QfW +#self.Qfcub = Qfcub +#self.Qtlag = Qtlag +##waterbalance + + + +[modelparameters] +# Format: +# name=stack,type,default,verbose[lookupmap_1],[lookupmap_2],lookupmap_n] +# example: +# RootingDepth=monthlyclim/ROOT,monthyclim,100,1 + +# - name - Name of the parameter (internal variable) +# - stack - Name of the mapstack (representation on disk or in mem) relative to case +# - type - Type of parameter (default = static) +# - default - Default value if map/tbl is not present +# - set to 1 to be verbose if a map is missing +# - lookupmap - maps to be used in the lookuptable in the case the type is statictbl + +#Possible types are:: +# - staticmap: Read at startup from map +# - statictbl: Read at startup from tbl +# - tbl: Read each timestep from tbl and at startup +# - timeseries: read map for each timestep +# - monthlyclim: read a map corresponding to the current month (12 maps in total) +# - dailyclim: read a map corresponding to the current day of the year +# - hourlyclim: read a map corresponding to the current hour of the day (24 in total) (not implemented yet) +# - tss: read a tss file and link to lookupmap (only one allowed) a map using timeinputscalar +Precipitation=intss/1_P.tss,tss,0.0,1,staticmaps/wflow_mgauges.map +Temperature=intss/1_T.tss,tss,10.5,1,staticmaps/wflow_mgauges.map +PotEvaporation=intss/1_PET.tss,tss,0.0,1,staticmaps/wflow_mgauges.map + Index: wflow-py/Sandbox/wflow_topoflex/reservoir_Sf.py =================================================================== diff -u -rd9373dad57e311dd6b2fefe99d78b6e91d750b74 -r08e6dbe2c7a682c7095170bd1bf89d3674fefb7c --- wflow-py/Sandbox/wflow_topoflex/reservoir_Sf.py (.../reservoir_Sf.py) (revision d9373dad57e311dd6b2fefe99d78b6e91d750b74) +++ wflow-py/Sandbox/wflow_topoflex/reservoir_Sf.py (.../reservoir_Sf.py) (revision 08e6dbe2c7a682c7095170bd1bf89d3674fefb7c) @@ -95,19 +95,19 @@ - Qf is devided over the reservoir numbers for the timesteps matching with the average travel time for a calcultation cell """ - if sum(pcr2numpy(self.Transit,NaN)) > 0: + if nansum(pcr2numpy(self.Transit,NaN)) > 0: self.Qflag = self.trackQ[0] # first bucket is transferred to outlet self.trackQ.append(0 * scalar(self.TopoId)) # add new bucket for present time step del self.trackQ[0] # remove first bucket (transferred to outlet) temp = [self.trackQ[i] + ifthenelse(rounddown(self.Transit) == i*scalar(self.TopoId), - (self.Transit - i*scalar(self.TopoId)) * self.Qfcub, + (self.Transit - i*scalar(self.TopoId)) * self.Qftotal / 1000 * self.surfaceArea, ifthenelse(roundup(self.Transit) == i*scalar(self.TopoId), - (i*scalar(self.TopoId) - self.Transit) * self.Qfcub, 0)) for i in range(len(self.trackQ))] + (i*scalar(self.TopoId) - self.Transit) * self.Qftotal / 1000 * self.surfaceArea, 0)) for i in range(len(self.trackQ))] self.trackQ = temp else: - self.Qflag = self.Qfcub + self.Qflag = self.Qftotal - self.wbSfrout = self.Qfcub - self.Qflag - sum(self.trackQ) + sum (self.trackQ_t) + self.wbSfrout = self.Qftotal - self.Qflag - sum(self.trackQ) + sum (self.trackQ_t) self.Qflag_ = self.Qflag \ No newline at end of file Index: wflow-py/Sandbox/wflow_topoflex/wflow_topoflex.py =================================================================== diff -u -rd9373dad57e311dd6b2fefe99d78b6e91d750b74 -r08e6dbe2c7a682c7095170bd1bf89d3674fefb7c --- wflow-py/Sandbox/wflow_topoflex/wflow_topoflex.py (.../wflow_topoflex.py) (revision d9373dad57e311dd6b2fefe99d78b6e91d750b74) +++ wflow-py/Sandbox/wflow_topoflex/wflow_topoflex.py (.../wflow_topoflex.py) (revision 08e6dbe2c7a682c7095170bd1bf89d3674fefb7c) @@ -29,12 +29,12 @@ import reservoir_Sus import JarvisCoefficients -import pdb import numpy import os import os.path import shutil, glob import getopt +import time from wflow.wf_DynamicFramework import * from wflow.wflow_adapt import * @@ -136,7 +136,7 @@ :var TSoil: Temperature of the soil [oC] """ - states = ['Si', 'Su', 'Sus', 'Sf', 'Ss'] + states = ['Si', 'Su', 'Sus', 'Sf', 'Ss', 'Qstate'] return states @@ -177,6 +177,7 @@ [report(self.Sf[i], self.SaveDir + "/outmaps/Sf" + self.NamesClasses[i] + ".map") for i in self.Classes] [report(self.Sr[i], self.SaveDir + "/outmaps/Sr" + self.NamesClasses[i] + ".map") for i in self.Classes] report(self.Ss, self.SaveDir + "/outmaps/Ss.map") + report(self.Qstate, self.SaveDir + "/outmaps/Qstate.map") [report(self.percent[i], self.SaveDir + "/outmaps/percent" + self.NamesClasses[i] + ".map") for i in self.Classes] @@ -207,8 +208,6 @@ #: in this model but always good to keep in mind. setglobaloption("unittrue") - self.teller = 0 - self.thestep = scalar(0) #: files to be used in case of timesries (scalar) input to the model # files for forcing data @@ -319,6 +318,8 @@ "model", "wflow_surfaceArea", "staticmaps/wflow_surfaceArea.map") wflow_transit = configget(self.config, "model", "wflow_transit", "staticmaps/wflow_transit.map") + wflow_velocity = configget(self.config, + "model", "wflow_velocity", "staticmaps/wflow_velocity.map") wflow_percent = [configget(self.config, "model", "wflow_percent_" + str(self.Classes[i]), "staticmaps/wflow_percent" + str(self.Classes[i]) + ".map") for i in self.Classes] @@ -341,6 +342,7 @@ self.totalArea = areatotal(self.surfaceArea, nominal(self.TopoId)) self.percentArea = self.surfaceArea / self.totalArea self.Transit = scalar(readmap(os.path.join(self.Dir, wflow_transit))) #: Map with surface area per cell + self.velocity = scalar(readmap(os.path.join(self.Dir, wflow_velocity))) #: Map with surface area per cell self.gaugesR = nominal(readmap(os.path.join(self.Dir, wflow_gauges))) self.percent = [] for i in self.Classes: @@ -448,6 +450,7 @@ self.Sr = [self.ZeroMap] * len(self.Classes) # self.Ss = [self.ZeroMap] * len(self.Classes) # for separate gw reservoir per class self.Ss = self.ZeroMap # for combined gw reservoir + self.Qstate = self.ZeroMap # for combined gw reservoir # set initial storage values self.Sa = [x + y for (x, y) in zip(self.Su, [130 * scalar(self.TopoId)] * len(self.Classes))] @@ -508,12 +511,12 @@ This is where all the time dependent functions are executed. Time dependent output should also be saved here. + :type self: object """ + # TODO: change rainfall .tss files into grids self.wf_updateparameters() # read the temperature map for each step (see parameters()) - self.teller = self.teller + 1 - - # self.logger.debug("Step: "+str(int(self.thestep + self._d_firstTimeStep))+"/"+str(int(self._d_nrTimeSteps))) + self.logger.debug("Step: "+str(int(self.thestep + self._d_firstTimeStep))+"/"+str(int(self._d_nrTimeSteps))) self.thestep = self.thestep + 1 self.Si_t = copylist(self.Si) @@ -602,18 +605,29 @@ # TOTAL RUNOFF ============================================================================================= # self.Qfcub = (sum([x*y for x,y in zip(self.Qf_,self.percent)]) + sum([x*y for x,y in zip(self.Qo_,self.percent)]) + sum([x*y for x,y in zip(self.Qd_,self.percent)]) + sum([x*y for x,y in zip(self.Qr_,self.percent)]))/ 1000 * self.surfaceArea - self.Qfcub = (sum([x * y for x, y in zip(self.Qf_, self.percent)])) / 1000 * self.surfaceArea + self.Qftotal = (sum([x * y for x, y in zip(self.Qf_, self.percent)])) + time_0 = time.time() reservoir_Sf.routingQf_combined(self) - + time_1 = time.time() + self.logger.info('Routing with Tanja buckets took {:f} seconds'.format(time_1-time_0)) if self.selectSs: eval_str = 'reservoir_Ss.{:s}(self)'.format(self.selectSs) else: eval_str = 'reservoir_Ss.groundWater_no_reservoir(self)' eval(eval_str) - # for separate gw reservoir per class - # self.Qtlag = self.Qflag_/ self.timestepsecs + sum([x*y for x,y in zip(self.Qs_,self.percent)])/ 1000 * self.surfaceArea/ self.timestepsecs - # for combinzed gw reservoir + # ROUTING + self.Qtot = self.Qftotal + self.Qs_ # total local discharge in mm/hour + self.Qtotal = self.Qtot / 1000 * self.surfaceArea / self.timestepsecs # total local discharge in m3/s + self.Qstate_t = self.Qstate + time_0 = time.time() + self.Qrout = accutraveltimeflux(self.TopoLdd, self.Qstate + self.Qtotal, self.velocity) + self.Qstate = accutraveltimestate(self.TopoLdd, self.Qstate + self.Qtotal, self.velocity) + time_1 = time.time() + self.logger.info('Routing with accutraveltimeflux took {:f} seconds'.format(time_1-time_0)) + # water balance of flux routing + self.dSdt = self.Qstate-self.Qstate_t + self.WB_rout = (accuflux(self.TopoLdd, self.Qtotal - self.dSdt)-self.Qrout)/accuflux(self.TopoLdd, self.Qtotal) self.Qtlag = self.Qflag_ / self.timestepsecs + self.Qs_ / 1000 * self.surfaceArea / self.timestepsecs self.QLagTot = areatotal(self.Qtlag, nominal(self.TopoId)) # catchment total runoff with looptijd @@ -640,8 +654,8 @@ self.Ei = areatotal(sum(multiply(self.Ei_, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId)) self.Eu = areatotal(sum(multiply(self.Eu_, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId)) self.Er = areatotal(sum(multiply(self.Er_, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId)) - self.QtotnoRout = areatotal(self.Qfcub + self.Qs_ / 1000 * self.surfaceArea, nominal(self.TopoId)) - self.Qtot = self.QLagTot * self.timestepsecs + #self.Qtot = self.QLagTot * self.timestepsecs + self.QtotnoRout = areatotal(self.Qtotal, nominal(self.TopoId)) self.SiWB = areatotal(sum(multiply(self.Si, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId)) self.Si_WB = areatotal(sum(multiply(self.Si_t, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId)) self.SuWB = areatotal(sum(multiply(self.Su, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId))