Index: wflow-py/Sandbox/wflow_w3ra_v2.py =================================================================== diff -u -r3cfc57a80ec1f8649e6f630c9e04cb8d42054cf1 -r3e39e84af48f1bcb5ec0d243748147be223674f2 --- wflow-py/Sandbox/wflow_w3ra_v2.py (.../wflow_w3ra_v2.py) (revision 3cfc57a80ec1f8649e6f630c9e04cb8d42054cf1) +++ wflow-py/Sandbox/wflow_w3ra_v2.py (.../wflow_w3ra_v2.py) (revision 3e39e84af48f1bcb5ec0d243748147be223674f2) @@ -38,13 +38,15 @@ from wflow.wf_DynamicFramework import * from wflow.wflow_adapt import * - -#TODO: Make the script HRU independent (loop over the nr of HRU's) -#TODO: +# TODO: Make the script HRU independent (loop over the nr of HRU's) +# TODO: + + def usage(*args): sys.stdout = sys.stderr - for msg in args: print(msg) + for msg in args: + print(msg) print(__doc__) sys.exit(0) @@ -54,34 +56,32 @@ define tanh for pcraster objects """ - return (exp(x)-exp(-x))/(exp(x) + exp(-x)) + return (exp(x) - exp(-x)) / (exp(x) + exp(-x)) -class WflowModel(DynamicModel): - """ +class WflowModel(DynamicModel): + """ The user defined model class. T """ - - def __init__(self, cloneMap,Dir,RunDir,configfile): - """ + + def __init__(self, cloneMap, Dir, RunDir, configfile): + """ *Required* The init function **must** contain what is shown below. Other functionality may be added by you if needed. """ - DynamicModel.__init__(self) - setclone(Dir + "/staticmaps/" + cloneMap) - self.runId=RunDir - self.caseName=Dir - self.Dir = Dir - self.configfile = configfile - self.SaveDir = self.Dir + "/" + self.runId + "/" - + DynamicModel.__init__(self) + setclone(Dir + "/staticmaps/" + cloneMap) + self.runId = RunDir + self.caseName = Dir + self.Dir = Dir + self.configfile = configfile + self.SaveDir = self.Dir + "/" + self.runId + "/" - - def stateVariables(self): - """ + def stateVariables(self): + """ *Required* Returns a list of state variables that are essential to the model. @@ -92,13 +92,22 @@ this function must return and empty array (states = []) """ - states = ['S0','Ss','Sd','Mleaf','FreeWater','DrySnow','Sg','Sr','OpenWaterFrac'] - - return states - + states = [ + "S0", + "Ss", + "Sd", + "Mleaf", + "FreeWater", + "DrySnow", + "Sg", + "Sr", + "OpenWaterFrac", + ] - def suspend(self): - """ + return states + + def suspend(self): + """ *Required* Suspends the model to disk. All variables needed to restart the model @@ -107,20 +116,20 @@ This function is required. """ - - self.logger.info("Saving initial conditions...") - #: It is advised to use the wf_suspend() function - #: here which will suspend the variables that are given by stateVariables - #: function. - self.wf_suspend(self.SaveDir + "/outstate/") - if self.fewsrun: - self.logger.info("Saving initial conditions for FEWS...") - self.wf_suspend(self.Dir + "/outstate/") - - def initial(self): - - """ + self.logger.info("Saving initial conditions...") + #: It is advised to use the wf_suspend() function + #: here which will suspend the variables that are given by stateVariables + #: function. + self.wf_suspend(self.SaveDir + "/outstate/") + + if self.fewsrun: + self.logger.info("Saving initial conditions for FEWS...") + self.wf_suspend(self.Dir + "/outstate/") + + def initial(self): + + """ *Required* Initial part of the model, executed only once. It reads all static model @@ -131,147 +140,255 @@ structure used in the other models. """ - #: pcraster option to calculate with units or cells. Not really an issue - #: in this model but always good to keep in mind. - setglobaloption("unittrue") - setglobaloption("radians") # Needed as W3RA was originally written in matlab + #: pcraster option to calculate with units or cells. Not really an issue + #: in this model but always good to keep in mind. + setglobaloption("unittrue") + setglobaloption("radians") # Needed as W3RA was originally written in matlab + self.timestepsecs = int( + configget(self.config, "model", "timestepsecs", "86400") + ) + self.UseETPdata = int( + configget(self.config, "model", "UseETPdata", "1") + ) # 1: Use ETP data, 0: Compute ETP from meteorological variables + self.logger.debug("use DATA: " + str(self.UseETPdata)) + self.basetimestep = 86400 + self.SaveMapDir = self.Dir + "/" + self.runId + "/outmaps" - self.timestepsecs = int(configget(self.config,'model','timestepsecs','86400')) - self.UseETPdata = int(configget(self.config,'model','UseETPdata','1')) # 1: Use ETP data, 0: Compute ETP from meteorological variables - self.logger.debug('use DATA: ' + str(self.UseETPdata)) - self.basetimestep=86400 - self.SaveMapDir = self.Dir + "/" + self.runId + "/outmaps" + # Define here the W3RA mapstacks (best to read these via netcdf) - # Define here the W3RA mapstacks (best to read these via netcdf) + self.TMAX_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "TMAX", "/inmaps/TMAX" + ) + self.TMIN_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "TMIN", "/inmaps/TMIN" + ) + self.TDAY_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "TDAY", "/inmaps/TDAY" + ) + self.EPOT_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "EPOT", "/inmaps/EPOT" + ) + self.PRECIP_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "PRECIP", "/inmaps/PRECIP" + ) + self.RAD_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "RAD", "/inmaps/RAD" + ) + # self.WINDSPEED_mapstack=self.Dir + configget(self.config,"inputmapstacks","WINDSPEED","/inmaps/ClimatologyMapFiles/WINDS/WNDSPEED") + # self.AIRPRESS_mapstack=self.Dir + configget(self.config,"inputmapstacks","AIRPRESS","/inmaps/ClimatologyMapFiles/AIRPRESS/AIRPRESS") + self.ALBEDO_mapstack = self.Dir + configget( + self.config, + "inputmapstacks", + "ALBEDO", + "/inmaps/ClimatologyMapFiles/ALBEDO/ALBEDO", + ) + self.WINDSPEED_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "WINDSPEED", "/inmaps/WIND" + ) + self.AIRPRESS_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "AIRPRESS", "/inmaps/PRES" + ) - self.TMAX_mapstack=self.Dir + configget(self.config,"inputmapstacks","TMAX","/inmaps/TMAX") - self.TMIN_mapstack=self.Dir + configget(self.config,"inputmapstacks","TMIN","/inmaps/TMIN") - self.TDAY_mapstack=self.Dir + configget(self.config,"inputmapstacks","TDAY","/inmaps/TDAY") - self.EPOT_mapstack=self.Dir + configget(self.config,"inputmapstacks","EPOT","/inmaps/EPOT") - self.PRECIP_mapstack=self.Dir + configget(self.config,"inputmapstacks","PRECIP","/inmaps/PRECIP") - self.RAD_mapstack=self.Dir + configget(self.config,"inputmapstacks","RAD","/inmaps/RAD") - #self.WINDSPEED_mapstack=self.Dir + configget(self.config,"inputmapstacks","WINDSPEED","/inmaps/ClimatologyMapFiles/WINDS/WNDSPEED") - #self.AIRPRESS_mapstack=self.Dir + configget(self.config,"inputmapstacks","AIRPRESS","/inmaps/ClimatologyMapFiles/AIRPRESS/AIRPRESS") - self.ALBEDO_mapstack=self.Dir + configget(self.config,"inputmapstacks","ALBEDO","/inmaps/ClimatologyMapFiles/ALBEDO/ALBEDO") - self.WINDSPEED_mapstack=self.Dir + configget(self.config,"inputmapstacks","WINDSPEED","/inmaps/WIND") - self.AIRPRESS_mapstack=self.Dir + configget(self.config,"inputmapstacks","AIRPRESS","/inmaps/PRES") + self.Altitude = readmap(self.Dir + "/staticmaps/wflow_dem") - self.Altitude=readmap(self.Dir + "/staticmaps/wflow_dem") + self.fewsrun = int(configget(self.config, "model", "fewsrun", "0")) - self.fewsrun = int(configget(self.config, "model", "fewsrun", "0")) + self.latitude = ycoordinate(boolean(self.Altitude)) - self.latitude = ycoordinate(boolean(self.Altitude)) + # Add reading of parameters here + self.ER_coef = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/ER_coef.map"), 0.0, fail=True + ) + self.flmp = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/flmp.map"), 0.0, fail=True + ) + self.fPotDeep = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/fPotDeep.map"), 0.0, fail=True + ) + self.FsoilEmax = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/FsoilEmax.map"), 0.0, fail=True + ) + self.Gs_scalar = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/Gs_scalar.map"), 0.0, fail=True + ) + self.hand_percentile0 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile0.map"), 0.0, fail=True + ) + self.hand_percentile1 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile1.map"), 0.0, fail=True + ) + self.hand_percentile2 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile2.map"), 0.0, fail=True + ) + self.hand_percentile3 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile3.map"), 0.0, fail=True + ) + self.hand_percentile4 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile4.map"), 0.0, fail=True + ) + self.hand_percentile5 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile5.map"), 0.0, fail=True + ) + self.hand_percentile6 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile6.map"), 0.0, fail=True + ) + self.hand_percentile7 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile7.map"), 0.0, fail=True + ) + self.hand_percentile8 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile8.map"), 0.0, fail=True + ) + self.hand_percentile9 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile9.map"), 0.0, fail=True + ) + self.hand_percentile10 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile10.map"), 0.0, fail=True + ) + self.hand_percentile11 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile11.map"), 0.0, fail=True + ) + self.hand_percentile12 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile12.map"), 0.0, fail=True + ) + self.hand_percentile13 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile13.map"), 0.0, fail=True + ) + self.hand_percentile14 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile14.map"), 0.0, fail=True + ) + self.hand_percentile15 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile15.map"), 0.0, fail=True + ) + self.hand_percentile16 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile16.map"), 0.0, fail=True + ) + self.hand_percentile17 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile17.map"), 0.0, fail=True + ) + self.hand_percentile18 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile18.map"), 0.0, fail=True + ) + self.hand_percentile19 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hand_percentile19.map"), 0.0, fail=True + ) + self.hveg = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hveg.map"), 0.0, fail=True + ) + self.K_gw = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/K_gw.map"), 0.0, fail=True + ) + self.k_s = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/k_s.map"), 0.0, fail=True + ) + self.K0_scalar = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/K0_scalar.map"), 0.0, fail=True + ) + self.Ksat_exp = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/Ksat_exp.map"), 0.0, fail=True + ) + self.lambd = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/lambda.map"), 0.0, fail=True + ) + self.OpenWaterFrac = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/OpenWaterFrac.map"), 0.0, fail=True + ) + self.porosity = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/porosity.map"), 0.0, fail=True + ) + self.Pref = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/pref.map"), 0.0, fail=True + ) + self.psi_s = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/psi_s.map"), 0.0, fail=True + ) + self.S_sls = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/S_sls.map"), 0.0, fail=True + ) + self.slope = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/slope.map"), 0.0, fail=True + ) + self.snow_Cfmax = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/snow_Cfmax.map"), 0.0, fail=True + ) + self.snow_Cfr = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/snow_Cfr.map"), 0.0, fail=True + ) + self.snow_TT = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/snow_TT.map"), 0.0, fail=True + ) + self.snow_WHC = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/snow_WHC.map"), 0.0, fail=True + ) + self.T_offset = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/T_offset.map"), 0.0, fail=True + ) + self.theta_s = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/theta_s.map"), 0.0, fail=True + ) + # self.K_rout = self.wf_readmap(os.path.join(self.Dir, "staticmaps/k_rout.map"),0.0,fail=True) + # self.Sgref = self.wf_readmap(os.path.join(self.Dir, "staticmaps/sgref.map"),0.0,fail=True) + # self.alb_dry = self.wf_readmap(os.path.join(self.Dir, "staticmaps/alb_dry.map"),0.0,fail=True) + # self.alb_wet = self.wf_readmap(os.path.join(self.Dir, "staticmaps/alb_wet.map"),0.0,fail=True) + # self.beta = self.wf_readmap(os.path.join(self.Dir, "staticmaps/beta.map"),0.0,fail=True) + # self.cGsmax = self.wf_readmap(os.path.join(self.Dir, "staticmaps/cgsmax.map"),0.0,fail=True) + # self.ER_frac_ref = self.wf_readmap(os.path.join(self.Dir, "staticmaps/er_frac_ref.map"),0.0,fail=True) + # self.Fhru = self.wf_readmap(os.path.join(self.Dir, "staticmaps/fHRU.map"),0.0,fail=True) + # self.FdrainFC = self.wf_readmap(os.path.join(self.Dir, "staticmaps/fdrainfc.map"),0.0,fail=True) + # self.Fgw_conn = self.wf_readmap(os.path.join(self.Dir, "staticmaps/fgw_conn.map"),0.0,fail=True) + # self.SLA = self.wf_readmap(os.path.join(self.Dir, "staticmaps/sla.map"),0.0,fail=True) + # self.LAIref = self.wf_readmap(os.path.join(self.Dir, "staticmaps/lairef.map"),0.0,fail=True) + # self.fvegref_G = self.wf_readmap(os.path.join(self.Dir, "staticmaps/fvegref_g.map"),0.0,fail=True) + # self.FwaterE = self.wf_readmap(os.path.join(self.Dir, "staticmaps/fwatere.map"),0.0,fail=True) + # self.Gfrac_max = self.wf_readmap(os.path.join(self.Dir, "staticmaps/gfrac_max.map"),0.0,fail=True) + # self.InitLoss = self.wf_readmap(os.path.join(self.Dir, "staticmaps/initloss.map"),0.0,fail=True) + # self.LAImax = self.wf_readmap(os.path.join(self.Dir, "staticmaps/laimax.map"),0.0,fail=True) + # self.S0FC = self.wf_readmap(os.path.join(self.Dir, "staticmaps/s0fc.map"),0.0,fail=True) + # self.SsFC = self.wf_readmap(os.path.join(self.Dir, "staticmaps/ssfc.map"),0.0,fail=True) + # self.SdFC = self.wf_readmap(os.path.join(self.Dir, "staticmaps/sdfc.map"),0.0,fail=True) + # self.Vc = self.wf_readmap(os.path.join(self.Dir, "staticmaps/vc.map"),0.0,fail=True) + # self.w0ref_alb = self.wf_readmap(os.path.join(self.Dir, "staticmaps/w0ref_alb.map"),0.0,fail=True) + # self.Us0 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/us0.map"),0.0,fail=True) + # self.Ud0 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/ud0.map"),0.0,fail=True) + # self.wslimU = self.wf_readmap(os.path.join(self.Dir, "staticmaps/wslimu.map"),0.0,fail=True) + # self.wdlimU = self.wf_readmap(os.path.join(self.Dir, "staticmaps/wdlimu.map"),0.0,fail=True) + # self.w0limE = self.wf_readmap(os.path.join(self.Dir, "staticmaps/w0lime.map"),0.0,fail=True) + # self.Tgrow = self.wf_readmap(os.path.join(self.Dir, "staticmaps/tgrow.map"),0.0,fail=True) + # self.Tsenc = self.wf_readmap(os.path.join(self.Dir, "staticmaps/tsenc.map"),0.0,fail=True) - # Add reading of parameters here - self.ER_coef = self.wf_readmap(os.path.join(self.Dir, "staticmaps/ER_coef.map"),0.0,fail=True) - self.flmp = self.wf_readmap(os.path.join(self.Dir, "staticmaps/flmp.map"),0.0,fail=True) - self.fPotDeep = self.wf_readmap(os.path.join(self.Dir, "staticmaps/fPotDeep.map"),0.0,fail=True) - self.FsoilEmax = self.wf_readmap(os.path.join(self.Dir, "staticmaps/FsoilEmax.map"),0.0,fail=True) - self.Gs_scalar = self.wf_readmap(os.path.join(self.Dir, "staticmaps/Gs_scalar.map"),0.0,fail=True) - self.hand_percentile0 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile0.map"),0.0,fail=True) - self.hand_percentile1 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile1.map"),0.0,fail=True) - self.hand_percentile2 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile2.map"),0.0,fail=True) - self.hand_percentile3 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile3.map"),0.0,fail=True) - self.hand_percentile4 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile4.map"),0.0,fail=True) - self.hand_percentile5 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile5.map"),0.0,fail=True) - self.hand_percentile6 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile6.map"),0.0,fail=True) - self.hand_percentile7 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile7.map"),0.0,fail=True) - self.hand_percentile8 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile8.map"),0.0,fail=True) - self.hand_percentile9 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile9.map"),0.0,fail=True) - self.hand_percentile10 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile10.map"),0.0,fail=True) - self.hand_percentile11 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile11.map"),0.0,fail=True) - self.hand_percentile12 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile12.map"),0.0,fail=True) - self.hand_percentile13 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile13.map"),0.0,fail=True) - self.hand_percentile14 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile14.map"),0.0,fail=True) - self.hand_percentile15 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile15.map"),0.0,fail=True) - self.hand_percentile16 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile16.map"),0.0,fail=True) - self.hand_percentile17 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile17.map"),0.0,fail=True) - self.hand_percentile18 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile18.map"),0.0,fail=True) - self.hand_percentile19 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hand_percentile19.map"),0.0,fail=True) - self.hveg = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hveg.map"),0.0,fail=True) - self.K_gw = self.wf_readmap(os.path.join(self.Dir, "staticmaps/K_gw.map"),0.0,fail=True) - self.k_s = self.wf_readmap(os.path.join(self.Dir, "staticmaps/k_s.map"),0.0,fail=True) - self.K0_scalar = self.wf_readmap(os.path.join(self.Dir, "staticmaps/K0_scalar.map"),0.0,fail=True) - self.Ksat_exp = self.wf_readmap(os.path.join(self.Dir, "staticmaps/Ksat_exp.map"),0.0,fail=True) - self.lambd = self.wf_readmap(os.path.join(self.Dir, "staticmaps/lambda.map"),0.0,fail=True) - self.OpenWaterFrac = self.wf_readmap(os.path.join(self.Dir, "staticmaps/OpenWaterFrac.map"),0.0,fail=True) - self.porosity = self.wf_readmap(os.path.join(self.Dir, "staticmaps/porosity.map"),0.0,fail=True) - self.Pref = self.wf_readmap(os.path.join(self.Dir, "staticmaps/pref.map"),0.0,fail=True) - self.psi_s = self.wf_readmap(os.path.join(self.Dir, "staticmaps/psi_s.map"),0.0,fail=True) - self.S_sls = self.wf_readmap(os.path.join(self.Dir, "staticmaps/S_sls.map"),0.0,fail=True) - self.slope = self.wf_readmap(os.path.join(self.Dir, "staticmaps/slope.map"),0.0,fail=True) - self.snow_Cfmax = self.wf_readmap(os.path.join(self.Dir, "staticmaps/snow_Cfmax.map"),0.0,fail=True) - self.snow_Cfr = self.wf_readmap(os.path.join(self.Dir, "staticmaps/snow_Cfr.map"),0.0,fail=True) - self.snow_TT = self.wf_readmap(os.path.join(self.Dir, "staticmaps/snow_TT.map"),0.0,fail=True) - self.snow_WHC = self.wf_readmap(os.path.join(self.Dir, "staticmaps/snow_WHC.map"),0.0,fail=True) - self.T_offset = self.wf_readmap(os.path.join(self.Dir, "staticmaps/T_offset.map"),0.0,fail=True) - self.theta_s = self.wf_readmap(os.path.join(self.Dir, "staticmaps/theta_s.map"),0.0,fail=True) - #self.K_rout = self.wf_readmap(os.path.join(self.Dir, "staticmaps/k_rout.map"),0.0,fail=True) - #self.Sgref = self.wf_readmap(os.path.join(self.Dir, "staticmaps/sgref.map"),0.0,fail=True) - #self.alb_dry = self.wf_readmap(os.path.join(self.Dir, "staticmaps/alb_dry.map"),0.0,fail=True) - #self.alb_wet = self.wf_readmap(os.path.join(self.Dir, "staticmaps/alb_wet.map"),0.0,fail=True) - #self.beta = self.wf_readmap(os.path.join(self.Dir, "staticmaps/beta.map"),0.0,fail=True) - #self.cGsmax = self.wf_readmap(os.path.join(self.Dir, "staticmaps/cgsmax.map"),0.0,fail=True) - #self.ER_frac_ref = self.wf_readmap(os.path.join(self.Dir, "staticmaps/er_frac_ref.map"),0.0,fail=True) - #self.Fhru = self.wf_readmap(os.path.join(self.Dir, "staticmaps/fHRU.map"),0.0,fail=True) - #self.FdrainFC = self.wf_readmap(os.path.join(self.Dir, "staticmaps/fdrainfc.map"),0.0,fail=True) - #self.Fgw_conn = self.wf_readmap(os.path.join(self.Dir, "staticmaps/fgw_conn.map"),0.0,fail=True) - #self.SLA = self.wf_readmap(os.path.join(self.Dir, "staticmaps/sla.map"),0.0,fail=True) - #self.LAIref = self.wf_readmap(os.path.join(self.Dir, "staticmaps/lairef.map"),0.0,fail=True) - #self.fvegref_G = self.wf_readmap(os.path.join(self.Dir, "staticmaps/fvegref_g.map"),0.0,fail=True) - #self.FwaterE = self.wf_readmap(os.path.join(self.Dir, "staticmaps/fwatere.map"),0.0,fail=True) - #self.Gfrac_max = self.wf_readmap(os.path.join(self.Dir, "staticmaps/gfrac_max.map"),0.0,fail=True) - #self.InitLoss = self.wf_readmap(os.path.join(self.Dir, "staticmaps/initloss.map"),0.0,fail=True) - #self.LAImax = self.wf_readmap(os.path.join(self.Dir, "staticmaps/laimax.map"),0.0,fail=True) - #self.S0FC = self.wf_readmap(os.path.join(self.Dir, "staticmaps/s0fc.map"),0.0,fail=True) - #self.SsFC = self.wf_readmap(os.path.join(self.Dir, "staticmaps/ssfc.map"),0.0,fail=True) - #self.SdFC = self.wf_readmap(os.path.join(self.Dir, "staticmaps/sdfc.map"),0.0,fail=True) - #self.Vc = self.wf_readmap(os.path.join(self.Dir, "staticmaps/vc.map"),0.0,fail=True) - #self.w0ref_alb = self.wf_readmap(os.path.join(self.Dir, "staticmaps/w0ref_alb.map"),0.0,fail=True) - #self.Us0 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/us0.map"),0.0,fail=True) - #self.Ud0 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/ud0.map"),0.0,fail=True) - #self.wslimU = self.wf_readmap(os.path.join(self.Dir, "staticmaps/wslimu.map"),0.0,fail=True) - #self.wdlimU = self.wf_readmap(os.path.join(self.Dir, "staticmaps/wdlimu.map"),0.0,fail=True) - #self.w0limE = self.wf_readmap(os.path.join(self.Dir, "staticmaps/w0lime.map"),0.0,fail=True) - #self.Tgrow = self.wf_readmap(os.path.join(self.Dir, "staticmaps/tgrow.map"),0.0,fail=True) - #self.Tsenc = self.wf_readmap(os.path.join(self.Dir, "staticmaps/tsenc.map"),0.0,fail=True) + self.wf_multparameters() + # Static, for the computation of Aerodynamic conductance (3.7) + self.fh = ln(813. / max(self.hveg, 0.25) - 5.45) + self.ku1 = 0.305 / (self.fh * (self.fh + 2.3)) - self.wf_multparameters() - # Static, for the computation of Aerodynamic conductance (3.7) - self.fh = ln(813./max(self.hveg,0.25)-5.45) - self.ku1 = 0.305/(self.fh*(self.fh+2.3)) + self.logger.info("Starting Dynamic run...") - self.logger.info("Starting Dynamic run...") - - - def resume(self): - """ + def resume(self): + """ *Required* This function is required. Read initial state maps (they are output of a previous call to suspend()). The implementation shown here is the most basic setup needed. """ - self.logger.info("Reading initial conditions...") - #: It is advised to use the wf_resume() function - #: here which pick up the variable save by a call to wf_suspend() - try: - self.wf_resume(self.Dir + "/instate/") - except: - self.logger.warning("Cannot load initial states, setting to default") - for s in self.stateVariables(): - exec("self." + s + " = cover(1.0)") + self.logger.info("Reading initial conditions...") + #: It is advised to use the wf_resume() function + #: here which pick up the variable save by a call to wf_suspend() + try: + self.wf_resume(self.Dir + "/instate/") + except: + self.logger.warning("Cannot load initial states, setting to default") + for s in self.stateVariables(): + exec("self." + s + " = cover(1.0)") - - def default_summarymaps(self): - """ + def default_summarymaps(self): + """ *Optional* Return a default list of variables to report as summary maps in the outsum dir. """ - return [] + return [] - def parameters(self): + def parameters(self): """ Define all model parameters here that the framework should handle for the model See wf_updateparameters and the parameters section of the ini file @@ -281,462 +398,578 @@ """ modelparameters = [] - #Static model parameters e.g. - #modelparameters.append(self.ParamType(name="RunoffGeneratingGWPerc",stack="intbl/RunoffGeneratingGWPerc.tbl",type="static",default=0.1)) + # Static model parameters e.g. + # modelparameters.append(self.ParamType(name="RunoffGeneratingGWPerc",stack="intbl/RunoffGeneratingGWPerc.tbl",type="static",default=0.1)) # 3: Input time series ################################################### - #self.P_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Precipitation", + # 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", + # 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", + # 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", + # self.Inflow_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Inflow", # "/inmaps/IF") # timeseries for rainfall "/inmaps/IF" # in/outflow locations (abstractions) # Meteo and other forcing - #modelparameters.append(self.ParamType(name="Precipitation",stack=self.P_mapstack,type="timeseries",default=0.0,verbose=True,lookupmaps=[])) - #modelparameters.append(self.ParamType(name="PotenEvap",stack=self.PET_mapstack,type="timeseries",default=0.0,verbose=True,lookupmaps=[])) - #modelparameters.append(self.ParamType(name="Temperature",stack=self.TEMP_mapstack,type="timeseries",default=10.0,verbose=True,lookupmaps=[])) - #modelparameters.append(self.ParamType(name="Inflow",stack=self.Inflow_mapstack,type="timeseries",default=0.0,verbose=False,lookupmaps=[])) + # modelparameters.append(self.ParamType(name="Precipitation",stack=self.P_mapstack,type="timeseries",default=0.0,verbose=True,lookupmaps=[])) + # modelparameters.append(self.ParamType(name="PotenEvap",stack=self.PET_mapstack,type="timeseries",default=0.0,verbose=True,lookupmaps=[])) + # modelparameters.append(self.ParamType(name="Temperature",stack=self.TEMP_mapstack,type="timeseries",default=10.0,verbose=True,lookupmaps=[])) + # modelparameters.append(self.ParamType(name="Inflow",stack=self.Inflow_mapstack,type="timeseries",default=0.0,verbose=False,lookupmaps=[])) return modelparameters - - - def dynamic(self): + def dynamic(self): """ *Required* This is where all the time dependent functions are executed. Time dependent output should also be saved here. """ - #print 'useETPdata' , self.UseETPdata - #Put the W3RA here. Stuff from W3RA_timestep_model.m - #read meteo from file + # print 'useETPdata' , self.UseETPdata + # Put the W3RA here. Stuff from W3RA_timestep_model.m + # read meteo from file self.logger.debug("Running for: " + str(self.currentdatetime)) - self.PRECIP=cover(self.wf_readmap(self.PRECIP_mapstack, 0.0), scalar(0.0)) # mm - + self.PRECIP = cover( + self.wf_readmap(self.PRECIP_mapstack, 0.0), scalar(0.0) + ) # mm + if self.UseETPdata == 1: - self.TDAY=cover(self.wf_readmap(self.TDAY_mapstack, 10.0), scalar(10.0)) # T in degC - self.EPOT=cover(self.wf_readmap(self.EPOT_mapstack, 0.0), scalar(0.0)) # mm - # print "Using climatology for wind, air pressure and albedo." + self.TDAY = cover( + self.wf_readmap(self.TDAY_mapstack, 10.0), scalar(10.0) + ) # T in degC + self.EPOT = cover( + self.wf_readmap(self.EPOT_mapstack, 0.0), scalar(0.0) + ) # mm + # print "Using climatology for wind, air pressure and albedo." elif self.UseETPdata == 0: - self.TMIN=cover(self.wf_readmap(self.TMIN_mapstack, 10.0), scalar(10.0)) # T in degC - self.TMAX=cover(self.wf_readmap(self.TMAX_mapstack, 10.0), scalar(10.0)) # T in degC - self.RAD=cover(self.wf_readmap(self.RAD_mapstack, 10.0), scalar(10.0))# W m-2 s-1 - self.WINDSPEED=cover(self.wf_readmap(self.WINDSPEED_mapstack, 10.0), scalar(10.0))# ms-1 - self.AIRPRESS=cover(self.wf_readmap(self.AIRPRESS_mapstack, 10.0), scalar(10.0))# Pa - self.ALBEDO=cover(self.wf_readmapClimatology(self.ALBEDO_mapstack, default=0.1), scalar(0.1)) + self.TMIN = cover( + self.wf_readmap(self.TMIN_mapstack, 10.0), scalar(10.0) + ) # T in degC + self.TMAX = cover( + self.wf_readmap(self.TMAX_mapstack, 10.0), scalar(10.0) + ) # T in degC + self.RAD = cover( + self.wf_readmap(self.RAD_mapstack, 10.0), scalar(10.0) + ) # W m-2 s-1 + self.WINDSPEED = cover( + self.wf_readmap(self.WINDSPEED_mapstack, 10.0), scalar(10.0) + ) # ms-1 + self.AIRPRESS = cover( + self.wf_readmap(self.AIRPRESS_mapstack, 10.0), scalar(10.0) + ) # Pa + self.ALBEDO = cover( + self.wf_readmapClimatology(self.ALBEDO_mapstack, default=0.1), + scalar(0.1), + ) - self.wf_multparameters() - doy=self.currentdatetime.timetuple().tm_yday + doy = self.currentdatetime.timetuple().tm_yday - #conversion daylength + # conversion daylength setglobaloption("radians") - m = scalar(1)-tan((self.latitude*scalar(pi)/scalar(180)))*tan(((scalar(23.439)*scalar(pi)/scalar(180))*cos(scalar(2)*scalar(pi)*(doy+scalar(9))/scalar(365.25)))) - self.fday = min(max(scalar(0.02),scalar(acos(scalar(1)-min(max(scalar(0),m),scalar(2))))/scalar(pi)),scalar(1)) #fraction daylength - + m = scalar(1) - tan((self.latitude * scalar(pi) / scalar(180))) * tan( + ( + (scalar(23.439) * scalar(pi) / scalar(180)) + * cos(scalar(2) * scalar(pi) * (doy + scalar(9)) / scalar(365.25)) + ) + ) + self.fday = min( + max( + scalar(0.02), + scalar(acos(scalar(1) - min(max(scalar(0), m), scalar(2)))) + / scalar(pi), + ), + scalar(1), + ) # fraction daylength # Assign forcing and estimate effective meteorological variables - Pg = self.PRECIP # mm - + Pg = self.PRECIP # mm + if self.UseETPdata == 1: Ta = self.TDAY # T in degC T24 = self.TDAY # T in degC elif self.UseETPdata == 0: - Rg = max(self.RAD,scalar(0.0001)) # already in W m-2 s-1; set minimum of 0.01 to avoid numerical problems - Ta = self.TMIN+scalar(0.75)*(self.TMAX-self.TMIN) # T in degC - T24 = self.TMIN+scalar(0.5)*(self.TMAX-self.TMIN) # T in degC - pex = min(scalar(17.27)*(self.TMIN)/(scalar(237.3)+self.TMIN),scalar(10)) # T in degC - pe = min(scalar(610.8)*(exp(pex)),scalar(10000.0)) # Mean actual vapour pressure, from dewpoint temperature + Rg = max( + self.RAD, scalar(0.0001) + ) # already in W m-2 s-1; set minimum of 0.01 to avoid numerical problems + Ta = self.TMIN + scalar(0.75) * (self.TMAX - self.TMIN) # T in degC + T24 = self.TMIN + scalar(0.5) * (self.TMAX - self.TMIN) # T in degC + pex = min( + scalar(17.27) * (self.TMIN) / (scalar(237.3) + self.TMIN), scalar(10) + ) # T in degC + pe = min( + scalar(610.8) * (exp(pex)), scalar(10000.0) + ) # Mean actual vapour pressure, from dewpoint temperature # windspeed is at 1m - #u2 = scalar(WindFactor)*self.WINDSPEED*(scalar(1)-(scalar(1)-self.fday)*scalar(0.25))/self.fday - self.u1 = self.WINDSPEED*(scalar(1)-(scalar(1)-self.fday)*scalar(0.25))/self.fday - pair = self.AIRPRESS # already in Pa + # u2 = scalar(WindFactor)*self.WINDSPEED*(scalar(1)-(scalar(1)-self.fday)*scalar(0.25))/self.fday + self.u1 = ( + self.WINDSPEED + * (scalar(1) - (scalar(1) - self.fday) * scalar(0.25)) + / self.fday + ) + pair = self.AIRPRESS # already in Pa - # diagnostic equations - self.LAI = self.SLA*self.Mleaf # (5.3) - fveg = max(1 - exp(-self.LAI/self.LAIref),0.000001) # (5.3) + self.LAI = self.SLA * self.Mleaf # (5.3) + fveg = max(1 - exp(-self.LAI / self.LAIref), 0.000001) # (5.3) # Vc = max(0,EVI-0.07)/fveg fsoil = 1 - fveg - w0 = self.S0/self.S0max # (2.) - ws = self.Ss/self.Ssmax # (2.1) - wd = self.Sd/self.Sdmax # (2.1) + w0 = self.S0 / self.S0max # (2.) + ws = self.Ss / self.Ssmax # (2.1) + wd = self.Sd / self.Sdmax # (2.1) + TotSnow = self.FreeWater + self.DrySnow + # wSnow = self.FreeWater/(TotSnow+1e-5) - TotSnow = self.FreeWater+self.DrySnow - #wSnow = self.FreeWater/(TotSnow+1e-5) - # Spatialise catchment fractions - #Sgfree = max(self.Sg,0.0) + # Sgfree = max(self.Sg,0.0) # JS: Not sure if this is translated properly.... - #for i=1:par.Nhru - ChannelSurface = min(0,(0.007*self.Sr**0.75)) + # for i=1:par.Nhru + ChannelSurface = min(0, (0.007 * self.Sr ** 0.75)) OpenWaterFrac = max(ChannelSurface, self.OpenWaterFrac) - + # !! HANDometric functions go here !! - #fsat = min(1.0,max(min(0.005,0.007*self.Sr**0.75),Sgfree/self.Sgref)) - #Sghru = self.Sg + # fsat = min(1.0,max(min(0.005,0.007*self.Sr**0.75),Sgfree/self.Sgref)) + # Sghru = self.Sg # CALCULATION OF PET # Conversions and coefficients (3.1) - pesx = min((scalar(17.27)*Ta/(scalar(237.3)+Ta)),scalar(10)) - pes = min(scalar((scalar(610.8))*exp(pesx)),scalar(10000)) # saturated vapour pressure - # fRH = pe/pes # relative air humidity -------------- check - cRE = 0.03449+4.27e-5*Ta - # Caero = self.fday*0.176*(1+Ta/209.1)*(pair-0.417*pe)*(1-fRH) -------------- check + pesx = min((scalar(17.27) * Ta / (scalar(237.3) + Ta)), scalar(10)) + pes = min( + scalar((scalar(610.8)) * exp(pesx)), scalar(10000) + ) # saturated vapour pressure + # fRH = pe/pes # relative air humidity -------------- check + cRE = 0.03449 + 4.27e-5 * Ta + # Caero = self.fday*0.176*(1+Ta/209.1)*(pair-0.417*pe)*(1-fRH) -------------- check # keps = 1.4e-3*((Ta/187)**2+Ta/107+1)*(6.36*pair+pe)/pes - ga = max(0.001, self.ku1*self.u1 ) + ga = max(0.001, self.ku1 * self.u1) - if self.UseETPdata == 1: - self.E0 = max(self.EPOT,0) - keps = 0.655E-3 * pair / pes # See Appendix A3 (http://www.clw.csiro.au/publications/waterforahealthycountry/2010/wfhc-aus-water-resources-assessment-system.pdf) -------------------------------- check! - + self.E0 = max(self.EPOT, 0) + keps = ( + 0.655E-3 * pair / pes + ) # See Appendix A3 (http://www.clw.csiro.au/publications/waterforahealthycountry/2010/wfhc-aus-water-resources-assessment-system.pdf) -------------------------------- check! + elif self.UseETPdata == 0: # Aerodynamic conductance (3.7) - ns_alb = self.ALBEDO # available as parameter climatology - Rgeff = Rg/self.fday + ns_alb = self.ALBEDO # available as parameter climatology + Rgeff = Rg / self.fday # shortwave radiation balance (3.2) - #alb_veg = 0.452*Vc - #alb_soil = alb_wet+(alb_dry-alb_wet)*exp(-w0/w0ref_alb) + # alb_veg = 0.452*Vc + # alb_soil = alb_wet+(alb_dry-alb_wet)*exp(-w0/w0ref_alb) # new equations for snow albedo - alb_snow = self.SNOWALBEDO # available as parameter climatology - fsnow = min(1.0,0.05*TotSnow) # assumed; ideally some lit research needed - #alb = fveg*alb_veg+(fsoil-fsnow)*alb_soil +fsnow*alb_snow - #alb = albedo - alb = (1-fsnow)*ns_alb +fsnow*alb_snow - RSn = (1-alb)*Rgeff + alb_snow = self.SNOWALBEDO # available as parameter climatology + fsnow = min( + 1.0, 0.05 * TotSnow + ) # assumed; ideally some lit research needed + # alb = fveg*alb_veg+(fsoil-fsnow)*alb_soil +fsnow*alb_snow + # alb = albedo + alb = (1 - fsnow) * ns_alb + fsnow * alb_snow + RSn = (1 - alb) * Rgeff # long wave radiation balance (3.3 to 3.5) StefBolz = 5.67e-8 - Tkelv = Ta+273.16 -# self.RLin = (0.65*(pe/Tkelv)**0.14)*StefBolz*Tkelv**4 # (3.3) - self.RLin - self.LWDOWN # available from meteo forcing - RLout = StefBolz*Tkelv**4.0 # (3.4) - self.RLn = self.RLin-RLout - - self.fGR = self.Gfrac_max*(1-exp(-fsoil/self.fvegref_G)) - self.Rneff = max(1, (RSn+self.RLn)*(1-self.fGR) ) - - fRH = pe/pes # relative air humidity - Caero = self.fday*0.176*(1+Ta/209.1)*(pair-0.417*pe)*(1-fRH) # -------------- check - keps = 1.4e-3*((Ta/187)**2+Ta/107+1)*(6.36*pair+pe)/pes - + Tkelv = Ta + 273.16 + # self.RLin = (0.65*(pe/Tkelv)**0.14)*StefBolz*Tkelv**4 # (3.3) + self.RLin - self.LWDOWN # available from meteo forcing + RLout = StefBolz * Tkelv ** 4.0 # (3.4) + self.RLn = self.RLin - RLout + + self.fGR = self.Gfrac_max * (1 - exp(-fsoil / self.fvegref_G)) + self.Rneff = max(1, (RSn + self.RLn) * (1 - self.fGR)) + + fRH = pe / pes # relative air humidity + Caero = ( + self.fday * 0.176 * (1 + Ta / 209.1) * (pair - 0.417 * pe) * (1 - fRH) + ) # -------------- check + keps = 1.4e-3 * ((Ta / 187) ** 2 + Ta / 107 + 1) * (6.36 * pair + pe) / pes + # Potential evaporation - kalpha = 1+Caero*ga/self.Rneff - self.E0 = cRE*(1/(1+keps))*kalpha*self.Rneff*self.fday - self.E0 = max(self.E0,0) - self.Eeq = cRE*(1/(1+keps))*1*self.Rneff*self.fday - self.Ept = cRE*(1/(1+keps))*1.26*self.Rneff*self.fday - + kalpha = 1 + Caero * ga / self.Rneff + self.E0 = cRE * (1 / (1 + keps)) * kalpha * self.Rneff * self.fday + self.E0 = max(self.E0, 0) + self.Eeq = cRE * (1 / (1 + keps)) * 1 * self.Rneff * self.fday + self.Ept = cRE * (1 / (1 + keps)) * 1.26 * self.Rneff * self.fday + # CALCULATION OF ET FLUXES AND ROOT WATER UPTAKE # Root water uptake constraint (4.4) U0max = scalar(0) - Usmax = max(0, self.Us0*min(1,ws/self.wslimU)) ##0-waarden omdat ws1 bevat 0-waarden (zie regel 116) - Udmax = max(0, self.Ud0*min(1,wd/self.wdlimU)) ##0-waarden omdat wd1 bevat 0-waarden (zie regel 118) - Ugmax = max(0, self.Ug0*max(0,fUg-fsat)) - Umax = max(Usmax, max(Udmax,Ugmax)) + Usmax = max( + 0, self.Us0 * min(1, ws / self.wslimU) + ) ##0-waarden omdat ws1 bevat 0-waarden (zie regel 116) + Udmax = max( + 0, self.Ud0 * min(1, wd / self.wdlimU) + ) ##0-waarden omdat wd1 bevat 0-waarden (zie regel 118) + Ugmax = max(0, self.Ug0 * max(0, fUg - fsat)) + Umax = max(Usmax, max(Udmax, Ugmax)) # Maximum transpiration (4.3) - Gsmax = self.Gs_scalar*self.cGsmax*self.Vc - VPD = max(0,pes-pe) - fD = self.Cg/(1+VPD/self.D50) - gs = fveg*fD*Gsmax - ft = 1/(1+(keps/(1+keps))*ga/gs) - Etmax = ft*self.E0 + Gsmax = self.Gs_scalar * self.cGsmax * self.Vc + VPD = max(0, pes - pe) + fD = self.Cg / (1 + VPD / self.D50) + gs = fveg * fD * Gsmax + ft = 1 / (1 + (keps / (1 + keps)) * ga / gs) + Etmax = ft * self.E0 # Actual transpiration (4.1) Et = min(Umax, Etmax) - + # # Root water uptake distribution (2.3) - U0 = max(min((U0max/(U0max + Usmax + Udmax + Ugmax))*Et, self.S0-1e-2),0) - Us = max(min((Usmax/(U0max + Usmax + Udmax + Ugmax))*Et, self.Ss-1e-2),0) - Ud = max(min((Udmax/(U0max + Usmax + Udmax + Ugmax))*Et, self.Sd-1e-2),0) - Ug = max(min((Ugmax/(U0max + Usmax + Udmax + Ugmax))*Et, self.Sd-1e-2),0) - Et = U0 + Us + Ud + Ug # to ensure mass balance + U0 = max(min((U0max / (U0max + Usmax + Udmax + Ugmax)) * Et, self.S0 - 1e-2), 0) + Us = max(min((Usmax / (U0max + Usmax + Udmax + Ugmax)) * Et, self.Ss - 1e-2), 0) + Ud = max(min((Udmax / (U0max + Usmax + Udmax + Ugmax)) * Et, self.Sd - 1e-2), 0) + Ug = max(min((Ugmax / (U0max + Usmax + Udmax + Ugmax)) * Et, self.Sd - 1e-2), 0) + Et = U0 + Us + Ud + Ug # to ensure mass balance # Soil evaporation (4.5) - w0x = max(0, (self.S0 - U0)/self.S0max ) # (2.1) - fsoilE = self.FsoilEmax*min(1,w0x/self.w0limE) - Es0 = (1-fsat)*fsoilE*(max(0,self.Eeq-Et)) + w0x = max(0, (self.S0 - U0) / self.S0max) # (2.1) + fsoilE = self.FsoilEmax * min(1, w0x / self.w0limE) + Es0 = (1 - fsat) * fsoilE * (max(0, self.Eeq - Et)) # Groundwater evaporation (4.6) - Eg0 = max(0,fsat-fwater)*self.FsoilEmax*max(0,self.Eeq-Et) + Eg0 = max(0, fsat - fwater) * self.FsoilEmax * max(0, self.Eeq - Et) Es = Es0 + Eg0 # Open water evaporation (4.7) - Erl = fw_local*self.FwaterE*self.Ept - Err = (fwater-fw_local)*self.FwaterE*self.Ept + Erl = fw_local * self.FwaterE * self.Ept + Err = (fwater - fw_local) * self.FwaterE * self.Ept Er = Erl + Err # Rainfall interception evaporation (4.2) - Sveg = self.S_sls*self.LAI - fER = fveg*self.ER_coeff*max(0.05,self.hveg)**self.ER_exp - Pwet = max(0, (scalar(Sveg>0 & fER>0 & (fER/fveg)<1)*-ln(1-fER/fveg)*Sveg/fER)) - Ei = scalar(self.T24>0)*(scalar(Pg=Pwet)*(fveg*Pwet+fER*(Pg-Pwet))) + Sveg = self.S_sls * self.LAI + fER = fveg * self.ER_coeff * max(0.05, self.hveg) ** self.ER_exp + Pwet = max( + 0, + ( + scalar(Sveg > 0 & fER > 0 & (fER / fveg) < 1) + * -ln(1 - fER / fveg) + * Sveg + / fER + ), + ) + Ei = scalar(self.T24 > 0) * ( + scalar(Pg < Pwet) * fveg * Pg + + scalar(Pg >= Pwet) * (fveg * Pwet + fER * (Pg - Pwet)) + ) - Edry=Et+Es+Er - self.EACT=Edry+Ei # for output only + Edry = Et + Es + Er + self.EACT = Edry + Ei # for output only # HBV snow routine # Matlab: function [FreeWater,DrySnow,InSoil]=snow_submodel(Precipitation,Temperature,FreeWater,DrySnow) # derived from HBV-96 shared by Jaap Schellekens (Deltares) in May 2011 # original in PCraster, adapted to Matlab by Albert van Dijk # HBV snow routine - Pn = Pg-Ei + Pn = Pg - Ei # Snow routine parameters # parameters # Partitioning into fractions rain and snow - Temperature = T24 # Dimmie, let op: tijdelijke regel!! - RainFrac = max(0,min((Temperature-(TT-TTI/2))/TTI,1)) - SnowFrac = 1 - RainFrac #fraction of precipitation which falls as snow + Temperature = T24 # Dimmie, let op: tijdelijke regel!! + RainFrac = max(0, min((Temperature - (TT - TTI / 2)) / TTI, 1)) + SnowFrac = 1 - RainFrac # fraction of precipitation which falls as snow # Snowfall/melt calculations - SnowFall = SnowFrac*Pn # snowfall depth - RainFall = RainFrac*Pn # rainfall depth - PotSnowMelt = Cfmax*max(0,Temperature-TT) # Potential snow melt, based on temperature - PotRefreezing = Cfmax*CFR*max(TT-Temperature,0) # Potential refreezing, based on temperature - Refreezing = min(PotRefreezing,self.FreeWater) # actual refreezing - SnowMelt = min(PotSnowMelt,self.DrySnow) # actual snow melt - self.DrySnow = self.DrySnow + SnowFall + Refreezing -SnowMelt # dry snow content - self.FreeWater = self.FreeWater - Refreezing # free water content in snow - MaxFreeWater = self.DrySnow*WHC + SnowFall = SnowFrac * Pn # snowfall depth + RainFall = RainFrac * Pn # rainfall depth + PotSnowMelt = Cfmax * max( + 0, Temperature - TT + ) # Potential snow melt, based on temperature + PotRefreezing = ( + Cfmax * CFR * max(TT - Temperature, 0) + ) # Potential refreezing, based on temperature + Refreezing = min(PotRefreezing, self.FreeWater) # actual refreezing + SnowMelt = min(PotSnowMelt, self.DrySnow) # actual snow melt + self.DrySnow = ( + self.DrySnow + SnowFall + Refreezing - SnowMelt + ) # dry snow content + self.FreeWater = self.FreeWater - Refreezing # free water content in snow + MaxFreeWater = self.DrySnow * WHC self.FreeWater = self.FreeWater + SnowMelt + RainFall - InSoil = max(self.FreeWater-MaxFreeWater,0) # abundant water in snow pack which goes into soil + InSoil = max( + self.FreeWater - MaxFreeWater, 0 + ) # abundant water in snow pack which goes into soil self.FreeWater = self.FreeWater - InSoil # End of Snow Module - Rmelt = scalar(Temperature<0)*InSoil # runs off if soil still frozen - Ps = scalar(Temperature>=0)*InSoil - + Rmelt = scalar(Temperature < 0) * InSoil # runs off if soil still frozen + Ps = scalar(Temperature >= 0) * InSoil + # CALCULATION OF WATER BALANCES # surface water fluxes (2.2) Rsof = fsat * Ps - Pi = max(0, Ps-self.InitLoss) - Rhof_soil = max(0,1-fsat-self.fImp)*(Pi - self.Pref*pcr_tanh(Pi/self.Pref)) # CHECK IF THIS GOES OK IN PYTHON - Rhof_imp = self.fImp*(Pi - self.Pref_imp*pcr_tanh(Pi/self.Pref_imp)) # CHECK IF THIS GOES OK IN PYTHON + Pi = max(0, Ps - self.InitLoss) + Rhof_soil = max(0, 1 - fsat - self.fImp) * ( + Pi - self.Pref * pcr_tanh(Pi / self.Pref) + ) # CHECK IF THIS GOES OK IN PYTHON + Rhof_imp = self.fImp * ( + Pi - self.Pref_imp * pcr_tanh(Pi / self.Pref_imp) + ) # CHECK IF THIS GOES OK IN PYTHON Rhof = Rhof_soil + Rhof_imp QR = Rhof + Rsof + Rmelt - I = Ps - Rhof - Rsof + I = Ps - Rhof - Rsof # SOIL WATER BALANCES (2.1 & 2.4) # Topsoil water balance (S0) # top soil layer hydrology - Kr_0s = self.K0sat/self.Kssat - Rh_0s = pcr_tanh(self.slope_coeff*self.slope*w0)*pcr_tanh(self.Kr_coeff*(Kr_0s-1)*w0) + Kr_0s = self.K0sat / self.Kssat + Rh_0s = pcr_tanh(self.slope_coeff * self.slope * w0) * pcr_tanh( + self.Kr_coeff * (Kr_0s - 1) * w0 + ) # general default case - Km = (self.K0sat*self.Kssat)**0.5 - A = Km/(self.S0max**2) + Km = (self.K0sat * self.Kssat) ** 0.5 + A = Km / (self.S0max ** 2) B = 1 - C = -(self.S0+I+Es) - S0 = (-B + ((B**2-4*A*C)**0.5))/(2*A) - D0 = (1-Rh_0s)*Km*((self.S0/self.S0max)**2) - IF0 = Rh_0s*Km*((self.S0/self.S0max)**2) + C = -(self.S0 + I + Es) + S0 = (-B + ((B ** 2 - 4 * A * C) ** 0.5)) / (2 * A) + D0 = (1 - Rh_0s) * Km * ((self.S0 / self.S0max) ** 2) + IF0 = Rh_0s * Km * ((self.S0 / self.S0max) ** 2) # depletion case - imap = (self.S0+I)<=Es - Es = ifthenelse(imap,(self.S0+I),Es) - S0 = ifthenelse(imap,0,S0) - D0 = ifthenelse(imap,0,D0) - IF0 = ifthenelse(imap,0,IF0) + imap = (self.S0 + I) <= Es + Es = ifthenelse(imap, (self.S0 + I), Es) + S0 = ifthenelse(imap, 0, S0) + D0 = ifthenelse(imap, 0, D0) + IF0 = ifthenelse(imap, 0, IF0) # saturation case - imap = (self.S0max-self.S0-self.K0sat)<=(I-Es) - D0 = ifthenelse(imap,(1-Rh_0s)*self.K0sat,D0) - IF0 = ifthenelse(imap,Rh_0s*self.K0sat+(self.S0-self.S0max-self.K0sat+I-Es),IF0) - S0 = ifthenselse(imap,self.S0max,S0) + imap = (self.S0max - self.S0 - self.K0sat) <= (I - Es) + D0 = ifthenelse(imap, (1 - Rh_0s) * self.K0sat, D0) + IF0 = ifthenelse( + imap, Rh_0s * self.K0sat + (self.S0 - self.S0max - self.K0sat + I - Es), IF0 + ) + S0 = ifthenselse(imap, self.S0max, S0) # enforce mass balance (for numerical & rounding errors - S0 = max(0, min(S0,self.S0max)) + S0 = max(0, min(S0, self.S0max)) massbal = self.S0 + I - Es - D0 - IF0 - S0 - D0 = D0 + (1-Rh_0s)*massbal - IF0 = IF0 + Rh_0s*massbal + D0 = D0 + (1 - Rh_0s) * massbal + IF0 = IF0 + Rh_0s * massbal self.S0 = S0 # Update state - + # shallow soil layer hydrology - Kr_sd = self.Kssat/self.Kdsat - Rh_sd = pcr_tanh(self.slope_coeff*self.slope*ws)*pcr_tanh(self.Kr_coeff*(Kr_sd-1)*ws) + Kr_sd = self.Kssat / self.Kdsat + Rh_sd = pcr_tanh(self.slope_coeff * self.slope * ws) * pcr_tanh( + self.Kr_coeff * (Kr_sd - 1) * ws + ) # general default case - Km = (self.Kssat*self.Kdsat)**0.5 - A = Km/(self.Ssmax**2) + Km = (self.Kssat * self.Kdsat) ** 0.5 + A = Km / (self.Ssmax ** 2) B = 1 - C = -(self.Ss+D0-Us) - Ss = (-B + ((B**2-4*A*C)**0.5))/(2*A) - Ds = (1-Rh_sd)*Km*((Ss/self.Ssmax)**2) - IFs = Rh_sd*Km*((Ss/self.Ssmax)**2) + C = -(self.Ss + D0 - Us) + Ss = (-B + ((B ** 2 - 4 * A * C) ** 0.5)) / (2 * A) + Ds = (1 - Rh_sd) * Km * ((Ss / self.Ssmax) ** 2) + IFs = Rh_sd * Km * ((Ss / self.Ssmax) ** 2) # depletion case - imap = (self.Ss+D0)<=Us - Us = ifthenelse(imap,(self.Ss+D0),Us) - Ss = ifthenelse(imap,0,Ss) - Ds = ifthenelse(imap,0,Ds) - IFs = ifthenelse(imap,0,IFs) + imap = (self.Ss + D0) <= Us + Us = ifthenelse(imap, (self.Ss + D0), Us) + Ss = ifthenelse(imap, 0, Ss) + Ds = ifthenelse(imap, 0, Ds) + IFs = ifthenelse(imap, 0, IFs) # saturation case - imap = (self.Ssmax-self.Ss-self.Kssat)<=(D0-Us) - Ds = ifthenelse(imap,(1-Rh_sd)*self.Kssat,Ds) - IFs = ifthenelse(imap,Rh_sd*self.Kssat+(self.Ss-self.Ssmax-self.Kssat+D0-Us),IFs) - Ss = ifthenselse(imap,self.Ssmax,Ss) + imap = (self.Ssmax - self.Ss - self.Kssat) <= (D0 - Us) + Ds = ifthenelse(imap, (1 - Rh_sd) * self.Kssat, Ds) + IFs = ifthenelse( + imap, + Rh_sd * self.Kssat + (self.Ss - self.Ssmax - self.Kssat + D0 - Us), + IFs, + ) + Ss = ifthenselse(imap, self.Ssmax, Ss) # enforce mass balance (for numerical & rounding errors - Ss = max(0, min(Ss,self.Ssmax)) + Ss = max(0, min(Ss, self.Ssmax)) massbal = self.Ss + D0 - Us - Ds - IFs - Ss - Ds = Ds + (1-Rh_sd)*massbal - IFs = IFs + Rh_sd*massbal + Ds = Ds + (1 - Rh_sd) * massbal + IFs = IFs + Rh_sd * massbal self.Ss = Ss # Update state # Deep soil layer hydrology - A = self.Kdsat/(self.Sdmax**2) + A = self.Kdsat / (self.Sdmax ** 2) B = 1 - C = -(self.Sd+Ds-Ud) - Sd = (-B + ((B**2-4*A*C)**0.5))/(2*A) - Dd = self.Kdsat*((Sd/self.Sdmax)**2) - IFd = 0*Dd; + C = -(self.Sd + Ds - Ud) + Sd = (-B + ((B ** 2 - 4 * A * C) ** 0.5)) / (2 * A) + Dd = self.Kdsat * ((Sd / self.Sdmax) ** 2) + IFd = 0 * Dd # depletion case - imap = (self.Sd+Ds)<=Ud - Ud = ifthenelse(imap,(self.Sd+Ds),Ud) - Sd = ifthenelse(imap,0,Sd) - Dd = ifthenelse(imap,0,Dd) - IFd = ifthenelse(imap,0,IFd) + imap = (self.Sd + Ds) <= Ud + Ud = ifthenelse(imap, (self.Sd + Ds), Ud) + Sd = ifthenelse(imap, 0, Sd) + Dd = ifthenelse(imap, 0, Dd) + IFd = ifthenelse(imap, 0, IFd) # saturation case - imap = (self.Sdmax-self.Sd-self.Kdsat)<=(Ds-Ud) - Dd = ifthenelse(imap,self.Kdsat,Dd) - IFd = ifthenelse(imap,self.Sd-self.Sdmax-self.Kdsat+Ds-Ud,IFd) - Sd = ifthenselse(imap,self.Sdmax,Sd) + imap = (self.Sdmax - self.Sd - self.Kdsat) <= (Ds - Ud) + Dd = ifthenelse(imap, self.Kdsat, Dd) + IFd = ifthenelse(imap, self.Sd - self.Sdmax - self.Kdsat + Ds - Ud, IFd) + Sd = ifthenselse(imap, self.Sdmax, Sd) # enforce mass balance (for numerical & rounding errors - Sd = max(0, min(Sd,self.Sdmax)) + Sd = max(0, min(Sd, self.Sdmax)) massbal = self.Sd + Ds - Ud - Dd - IFd - Sd Dd = Dd + massbal self.Sd = Sd # Update state - IFs = IFs + IFd ; # add up to interflow - QR = QR + IF0 + IFs; - + IFs = IFs + IFd + # add up to interflow + QR = QR + IF0 + IFs + # CATCHMENT WATER BALANCE # Groundwater store water balance (Sg) (2.5) NetGf = Dd - Eg - Ug self.Sg = self.Sg + NetGf - Sgfree = max(self.Sg,0) - Qg = min(Sgfree, (1-exp(-self.K_gw))*Sgfree) + Sgfree = max(self.Sg, 0) + Qg = min(Sgfree, (1 - exp(-self.K_gw)) * Sgfree) self.Sg = self.Sg - Qg # Surface water store water balance (Sr) (2.7) - self.Sr = max(0, self.Sr + QR - Er + Qg ) - self.Qtot = max(0, min(self.Sr, (1-exp(-self.K_rout))*self.Sr) ) - self.Sr = self.Sr - self.Qtot + self.Sr = max(0, self.Sr + QR - Er + Qg) + self.Qtot = max(0, min(self.Sr, (1 - exp(-self.K_rout)) * self.Sr)) + self.Sr = self.Sr - self.Qtot # VEGETATION ADJUSTMENT (5) - fvmax = 1-exp(-max(self.LAImax,0.002778)/self.LAIref) - fveq = (1/max((self.E0/Umax)-1,1e-3))*(keps/(1+keps))*(ga/(fd*Gsmax)) - fveq = min(fveq,fvmax) - dMleaf = -ln(1-fveq)*self.LAIref/self.SLA-self.Mleaf + fvmax = 1 - exp(-max(self.LAImax, 0.002778) / self.LAIref) + fveq = ( + (1 / max((self.E0 / Umax) - 1, 1e-3)) + * (keps / (1 + keps)) + * (ga / (fd * Gsmax)) + ) + fveq = min(fveq, fvmax) + dMleaf = -ln(1 - fveq) * self.LAIref / self.SLA - self.Mleaf - #Mleafnet1 = dMleaf1 * (dMleaf1/self.Tgrow1) + dMleaf1 * dMleaf1/self.Tsenc1 - #Mleafnet2 = dMleaf2 * (dMleaf1/self.Tgrow2) + dMleaf2 * dMleaf2/self.Tsenc2 - Mleafnet = scalar(dMleaf>0)*(dMleaf/self.Tgrow) +scalar(dMleaf<0)*dMleaf/self.Tsenc + # Mleafnet1 = dMleaf1 * (dMleaf1/self.Tgrow1) + dMleaf1 * dMleaf1/self.Tsenc1 + # Mleafnet2 = dMleaf2 * (dMleaf1/self.Tgrow2) + dMleaf2 * dMleaf2/self.Tsenc2 + Mleafnet = ( + scalar(dMleaf > 0) * (dMleaf / self.Tgrow) + + scalar(dMleaf < 0) * dMleaf / self.Tsenc + ) - self.Mleaf = self.Mleaf + Mleafnet - self.LAI = self.SLA*self.Mleaf # (5.3) + self.LAI = self.SLA * self.Mleaf # (5.3) - fveg = 1 - exp(-self.LAI1/self.LAIref1) #(5.3) + fveg = 1 - exp(-self.LAI1 / self.LAIref1) # (5.3) # in case this is desired as output: - self.w0 = self.S0/self.S0max - self.TotSnow = self.DrySnow + self.FreeWater + self.w0 = self.S0 / self.S0max + self.TotSnow = self.DrySnow + self.FreeWater + # The main function is used to run the program from the command line -def main(argv=None): + +def main(argv=None): """ *Optional* Perform command line execution of the model. This example uses the getopt module to parse the command line options. The user can set the caseName, the runDir, the timestep and the configfile. - """ + """ global multpars - caseName = "../openstreams_w3ra" # "D:/trambaue/_Projects/GLOFFIS/201501/GLOFFIS_SA/Modules/openstreams_w3ra/" + caseName = ( + "../openstreams_w3ra" + ) # "D:/trambaue/_Projects/GLOFFIS/201501/GLOFFIS_SA/Modules/openstreams_w3ra/" runId = "run_default" - configfile="wflow_W3RA.ini" - _lastTimeStep = 15 - _firstTimeStep = 0 - timestepsecs=86400 + configfile = "wflow_W3RA.ini" + _lastTimeStep = 15 + _firstTimeStep = 0 + timestepsecs = 86400 fewsrun = False - wflow_cloneMap = 'wflow_subcatch.map' + wflow_cloneMap = "wflow_subcatch.map" runinfoFile = "runinfo.xml" - _NoOverWrite=False + _NoOverWrite = False loglevel = logging.DEBUG LogFileName = "wflow.log" - - # This allows us to use the model both on the command line and to call + # This allows us to use the model both on the command line and to call # the model usinge main function from another python script. - + if argv is None: argv = sys.argv[1:] if len(argv) == 0: usage() - return + return - opts, args = getopt.getopt(argv, 'C:S:T:c:s:R:F:') - + opts, args = getopt.getopt(argv, "C:S:T:c:s:R:F:") + for o, a in opts: - if o == '-F': + if o == "-F": runinfoFile = a fewsrun = True - if o == '-C': caseName = a - if o == '-R': runId = a - if o == '-c': configfile = a - if o == '-s': timestepsecs = int(a) - if o == '-T': _lastTimeStep=int(a) - if o == '-S': _firstTimeStep=int(a) - - if (len(opts) <=1): + 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 len(opts) <= 1: usage() if fewsrun: ts = getTimeStepsfromRuninfo(runinfoFile, timestepsecs) starttime = getStartTimefromRuninfo(runinfoFile) - if (ts): + if ts: _lastTimeStep = ts _firstTimeStep = 1 else: print("Failed to get timesteps from runinfo file: " + runinfoFile) exit(2) else: - starttime = dt.datetime(1990,1,1) + starttime = dt.datetime(1990, 1, 1) if _lastTimeStep < _firstTimeStep: - print("The starttimestep (" + str(_firstTimeStep) + ") is smaller than the last timestep (" + str( - _lastTimeStep) + ")") + print( + "The starttimestep (" + + str(_firstTimeStep) + + ") is smaller than the last timestep (" + + str(_lastTimeStep) + + ")" + ) usage() - myModel = WflowModel(wflow_cloneMap, caseName,runId,configfile) - dynModelFw = wf_DynamicFramework(myModel, _lastTimeStep,firstTimestep=_firstTimeStep,datetimestart=starttime) - dynModelFw.createRunId(NoOverWrite=_NoOverWrite, level=loglevel, logfname=LogFileName,model="wflow_W3RA",doSetupFramework=False) + myModel = WflowModel(wflow_cloneMap, caseName, runId, configfile) + dynModelFw = wf_DynamicFramework( + myModel, _lastTimeStep, firstTimestep=_firstTimeStep, datetimestart=starttime + ) + dynModelFw.createRunId( + NoOverWrite=_NoOverWrite, + level=loglevel, + logfname=LogFileName, + model="wflow_W3RA", + 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 == "-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) dynModelFw.setupFramework() dynModelFw._runInitial() dynModelFw._runResume() - dynModelFw._runDynamic(0,0) + dynModelFw._runDynamic(0, 0) dynModelFw._runSuspend() dynModelFw._wf_shutdown() - + if __name__ == "__main__": - main() \ No newline at end of file + main()