Index: wflow-py/Sandbox/wflow_w3.py =================================================================== diff -u -r4c8b12536d7241d885afa05e95909a4ec0d755c3 -r12ea40dc08628f654753679e0972e87b7bb12f7a --- wflow-py/Sandbox/wflow_w3.py (.../wflow_w3.py) (revision 4c8b12536d7241d885afa05e95909a4ec0d755c3) +++ wflow-py/Sandbox/wflow_w3.py (.../wflow_w3.py) (revision 12ea40dc08628f654753679e0972e87b7bb12f7a) @@ -43,13 +43,14 @@ from wflow.wflow_adapt import * +# TODO: Make the script HRU independent (loop over the nr of HRU's) +# TODO: Add self.LAImax -#TODO: Make the script HRU independent (loop over the nr of HRU's) -#TODO: Add self.LAImax def usage(*args): sys.stdout = sys.stderr - for msg in args: print msg + for msg in args: + print msg print __doc__ sys.exit(0) @@ -59,53 +60,59 @@ define tanh for pcraster objects """ - return (exp(x)-exp(-x))/(exp(x) + exp(-x)) - -def interp_hand(z,hand,hand_perc): - - - z_lim = xarray.ufuncs.minimum( xarray.ufuncs.maximum(z,hand[0]), hand[hand_perc.size-1] ) # limit values within measured elevation range + return (exp(x) - exp(-x)) / (exp(x) + exp(-x)) - iLower = (hand.where(hand <= z_lim)) # find next lower elevation - PercLower = ((iLower*0+1.0).where(iLower==iLower.max(axis=0)) * hand_perc).max(axis=0, skipna = True) - zLower = iLower.where(iLower==iLower.max(axis=0)).max(axis=0, skipna = True) - iUpper = (hand.where(hand >= z_lim)) # find next higher elevation - PercUpper = ((iUpper*0+1.0).where(iUpper==iUpper.min(axis=0)) * hand_perc).max(axis=0, skipna = True) - zUpper = iUpper.where(iUpper==iUpper.min(axis=0)).max(axis=0, skipna = True) - - flim = PercLower + (PercUpper - PercLower) * xarray.ufuncs.fmax(0, xarray.ufuncs.fmin(1, (z_lim - zLower) / (zUpper - zLower))) +def interp_hand(z, hand, hand_perc): - pcr_flim = numpy2pcr(Scalar,flim.fillna(-999.0).values,-999.0) + z_lim = xarray.ufuncs.minimum( + xarray.ufuncs.maximum(z, hand[0]), hand[hand_perc.size - 1] + ) # limit values within measured elevation range + iLower = hand.where(hand <= z_lim) # find next lower elevation + PercLower = ( + (iLower * 0 + 1.0).where(iLower == iLower.max(axis=0)) * hand_perc + ).max(axis=0, skipna=True) + zLower = iLower.where(iLower == iLower.max(axis=0)).max(axis=0, skipna=True) + + iUpper = hand.where(hand >= z_lim) # find next higher elevation + PercUpper = ( + (iUpper * 0 + 1.0).where(iUpper == iUpper.min(axis=0)) * hand_perc + ).max(axis=0, skipna=True) + zUpper = iUpper.where(iUpper == iUpper.min(axis=0)).max(axis=0, skipna=True) + + flim = PercLower + (PercUpper - PercLower) * xarray.ufuncs.fmax( + 0, xarray.ufuncs.fmin(1, (z_lim - zLower) / (zUpper - zLower)) + ) + + pcr_flim = numpy2pcr(Scalar, flim.fillna(-999.0).values, -999.0) + return pcr_flim -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. @@ -116,13 +123,21 @@ 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 @@ -131,20 +146,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 @@ -155,189 +170,361 @@ 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_clone") - #self.Altitude=readmap(self.Dir + "/staticmaps/wflow_dem") - self.Altitude=readmap(self.Dir + "/staticmaps/wflow_clone") + 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 - # Add reading of parameters here + self.Fhru = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/fHRU.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.OpenWaterFrac = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/OpenWaterFrac.map"), 0.0, fail=True + ) + self.slope = ( + self.wf_readmap( + os.path.join(self.Dir, "staticmaps/slope.map"), 0.0, fail=True + ) + / 100.0 + ) + self.hveg = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/hveg.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.ER_coeff = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/ER_coeff.map"), 0.0, fail=True + ) + self.FsoilEmax = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/FsoilEmax.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.k_s = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/k_s.map"), 0.0, fail=True + ) + self.Lambda = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/lambda.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.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.fImp = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/fImp.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.fPotDeep = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/fPotDeep.map"), 0.0, fail=True + ) + self.porosity = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/porosity.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.theta_s = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/theta_s.map"), 0.0, fail=True + ) - self.Fhru = self.wf_readmap(os.path.join(self.Dir, "staticmaps/fHRU.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.OpenWaterFrac = self.wf_readmap(os.path.join(self.Dir, "staticmaps/OpenWaterFrac.map"),0.0,fail=True) - self.slope = self.wf_readmap(os.path.join(self.Dir, "staticmaps/slope.map"),0.0,fail=True)/100.0 - self.hveg = self.wf_readmap(os.path.join(self.Dir, "staticmaps/hveg.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.ER_coeff = self.wf_readmap(os.path.join(self.Dir, "staticmaps/ER_coeff.map"),0.0,fail=True) - self.FsoilEmax = self.wf_readmap(os.path.join(self.Dir, "staticmaps/FsoilEmax.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.k_s = self.wf_readmap(os.path.join(self.Dir, "staticmaps/k_s.map"),0.0,fail=True) - self.Lambda = self.wf_readmap(os.path.join(self.Dir, "staticmaps/lambda.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.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.fImp = self.wf_readmap(os.path.join(self.Dir, "staticmaps/fImp.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.fPotDeep = self.wf_readmap(os.path.join(self.Dir, "staticmaps/fPotDeep.map"),0.0,fail=True) - self.porosity = self.wf_readmap(os.path.join(self.Dir, "staticmaps/porosity.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.theta_s = self.wf_readmap(os.path.join(self.Dir, "staticmaps/theta_s.map"),0.0,fail=True) - - - self.alb_dry = self.wf_readmap(os.path.join(self.Dir, "staticmaps/alb_dry.map"),0.20,fail=False) - self.alb_wet = self.wf_readmap(os.path.join(self.Dir, "staticmaps/alb_wet.map"),0.15,fail=False) - self.alb_snow = self.wf_readmap(os.path.join(self.Dir, "staticmaps/alb_snow.map"),0.60,fail=False) - self.alb_water = self.wf_readmap(os.path.join(self.Dir, "staticmaps/alb_water.map"),0.05,fail=False) - self.Cg = self.wf_readmap(os.path.join(self.Dir, "staticmaps/Cg.map"),1.940,fail=False) - self.cGsmax = self.wf_readmap(os.path.join(self.Dir, "staticmaps/cGsmax.map"),0.020,fail=False) - self.d0 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/d0.map"),0.15,fail=False) - self.ds = self.wf_readmap(os.path.join(self.Dir, "staticmaps/ds.map"),0.85,fail=False) - self.dd = self.wf_readmap(os.path.join(self.Dir, "staticmaps/dd.map"),4.00,fail=False) - self.D50 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/D50.map"),700,fail=False) - self.ER_exp = self.wf_readmap(os.path.join(self.Dir, "staticmaps/ER_exp.map"),0.114,fail=False) - self.f_alb_Vc = self.wf_readmap(os.path.join(self.Dir, "staticmaps/f_alb_Vc.map"),0.4,fail=False) - self.Fgw_conn = self.wf_readmap(os.path.join(self.Dir, "staticmaps/Fgw_conn.map"),1,fail=False) - self.fvegref_G = self.wf_readmap(os.path.join(self.Dir, "staticmaps/fvegref_G.map"),0.15,fail=False) - self.FwaterE = self.wf_readmap(os.path.join(self.Dir, "staticmaps/FwaterE.map"),1,fail=False) - self.Gfrac_max = self.wf_readmap(os.path.join(self.Dir, "staticmaps/Gfrac_max.map"),0.15,fail=False) - self.InitLoss = self.wf_readmap(os.path.join(self.Dir, "staticmaps/InitLoss.map"),0,fail=False) - self.K_rout = self.wf_readmap(os.path.join(self.Dir, "staticmaps/K_rout.map"),0.5,fail=False) - self.Kr_coeff = self.wf_readmap(os.path.join(self.Dir, "staticmaps/Kr_coeff.map"),0.0741,fail=False) - self.LAIref = self.wf_readmap(os.path.join(self.Dir, "staticmaps/LAIref.map"),2.4,fail=False) - self.LUEmax = self.wf_readmap(os.path.join(self.Dir, "staticmaps/LUEmax.map"),0.0544,fail=False) - self.Pref_imp = self.wf_readmap(os.path.join(self.Dir, "staticmaps/Pref_imp.map"),10,fail=False) - self.R0 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/R0.map"),0.789,fail=False) - self.SLA = self.wf_readmap(os.path.join(self.Dir, "staticmaps/SLA.map"),5,fail=False) - self.slope_coeff = self.wf_readmap(os.path.join(self.Dir, "staticmaps/slope_coeff.map"),0.9518,fail=False) - self.snow_TTI = self.wf_readmap(os.path.join(self.Dir, "staticmaps/snow_TTI.map"),1,fail=False) - self.T24_snow = self.wf_readmap(os.path.join(self.Dir, "staticmaps/T24_snow.map"),18,fail=False) - self.Tmin = self.wf_readmap(os.path.join(self.Dir, "staticmaps/Tmin.map"),-10,fail=False) - self.Topt = self.wf_readmap(os.path.join(self.Dir, "staticmaps/Topt.map"),10,fail=False) - self.Tgrow = self.wf_readmap(os.path.join(self.Dir, "staticmaps/Tgrow.map"),200,fail=False) - self.Tsenc = self.wf_readmap(os.path.join(self.Dir, "staticmaps/Tsenc.map"),20,fail=False) - self.Ud0 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/Ud0.map"),6,fail=False) - self.Ug0 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/Ug0.map"),1,fail=False) - self.Us0 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/Us0.map"),6,fail=False) - self.Vc = self.wf_readmap(os.path.join(self.Dir, "staticmaps/Vc.map"),0.5,fail=False) - self.w0ref_alb = self.wf_readmap(os.path.join(self.Dir, "staticmaps/w0ref_alb.map"),0.3,fail=False) - - - ds_hand = xarray.open_dataset(os.path.join(self.Dir, "staticmaps/HAND.nc")) - - hand = ds_hand['HAND'] - self.HAND = xarray.concat([(hand[0]*0.0).expand_dims('z'),hand],dim='z') - - perc_HAND = ds_hand['percentile'] - self.perc_HAND = xarray.concat([(perc_HAND[0]*0.0).expand_dims('z'),perc_HAND],dim='z') - - psi_FC = self.wf_readmap(os.path.join(self.Dir, "staticmaps/psi_FC.map"),-3.3,fail=False) # m or hPa or 33 kPa - psi_FC0 = self.wf_readmap(os.path.join(self.Dir, "staticmaps/psi_FC0.map"),-0.5,fail=False) # m or hPa or 5 kPa - rapidly drainable theta for top soil - psi_ERRP = self.wf_readmap(os.path.join(self.Dir, "staticmaps/psi_ERRP.map"),-10,fail=False) # m or 100 kPa - assumed pressure at which soil moisture starts to limit soil evaporation (following D. Tran, 2015) - psi_d = self.wf_readmap(os.path.join(self.Dir, "staticmaps/psi_d.map"),-50,fail=False) # m assumed pressure at which soil moisture starts to limit soil water uptake - psi_PWP = self.wf_readmap(os.path.join(self.Dir, "staticmaps/psi_PWP.map"),-150,fail=False) # m - psi_res = self.wf_readmap(os.path.join(self.Dir, "staticmaps/psi_res.map"),-1e6,fail=False) # m - - - theta_FC = self.theta_s * min(1, (self.psi_s / psi_FC))**self.Lambda #fraction - theta_FC0 = self.theta_s * min(1, (self.psi_s / psi_FC0 ))**self.Lambda # fraction - theta_ERRP = self.theta_s *min(1, (self.psi_s / psi_ERRP ))**self.Lambda # fraction - theta_d = self.theta_s *min(1, (self.psi_s / psi_d ))**self.Lambda # fraction - theta_PWP = self.theta_s *min(1, (self.psi_s / psi_PWP ))**self.Lambda # fraction - theta_res = self.theta_s *min(1, (self.psi_s / psi_res ))**self.Lambda # fraction - - self.S0max = self.d0 * 1000 * (theta_FC0 - theta_res) # mm available storage for evaporation, note FC0 is used rather than theta_sat - self.Ssmax = self.ds * 1000 * (self.theta_s - theta_PWP); - self.Sdmax = self.dd * 1000 * (self.theta_s - theta_PWP); - self.K0sat = self.K0_scalar * self.k_s # mm/d - note that this is technically in fact not Ksat but K(theta_FC0) - self.Kssat = self.K0_scalar * (((self.ds + self.d0)/self.d0)**-self.Ksat_exp) * self.k_s - self.Kdsat = self.K0_scalar * (((self.dd + self.ds + self.d0)/self.d0)**-self.Ksat_exp) * self.k_s - self.w0limE = (theta_ERRP - theta_res) / (self.theta_s - theta_res) - self.wslimU = (theta_d - theta_PWP) /(self.theta_s - theta_PWP) - self.wdlimU = (theta_d - theta_PWP) /(self.theta_s - theta_PWP) + self.alb_dry = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/alb_dry.map"), 0.20, fail=False + ) + self.alb_wet = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/alb_wet.map"), 0.15, fail=False + ) + self.alb_snow = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/alb_snow.map"), 0.60, fail=False + ) + self.alb_water = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/alb_water.map"), 0.05, fail=False + ) + self.Cg = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/Cg.map"), 1.940, fail=False + ) + self.cGsmax = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/cGsmax.map"), 0.020, fail=False + ) + self.d0 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/d0.map"), 0.15, fail=False + ) + self.ds = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/ds.map"), 0.85, fail=False + ) + self.dd = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/dd.map"), 4.00, fail=False + ) + self.D50 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/D50.map"), 700, fail=False + ) + self.ER_exp = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/ER_exp.map"), 0.114, fail=False + ) + self.f_alb_Vc = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/f_alb_Vc.map"), 0.4, fail=False + ) + self.Fgw_conn = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/Fgw_conn.map"), 1, fail=False + ) + self.fvegref_G = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/fvegref_G.map"), 0.15, fail=False + ) + self.FwaterE = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/FwaterE.map"), 1, fail=False + ) + self.Gfrac_max = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/Gfrac_max.map"), 0.15, fail=False + ) + self.InitLoss = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/InitLoss.map"), 0, fail=False + ) + self.K_rout = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/K_rout.map"), 0.5, fail=False + ) + self.Kr_coeff = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/Kr_coeff.map"), 0.0741, fail=False + ) + self.LAIref = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/LAIref.map"), 2.4, fail=False + ) + self.LUEmax = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/LUEmax.map"), 0.0544, fail=False + ) + self.Pref_imp = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/Pref_imp.map"), 10, fail=False + ) + self.R0 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/R0.map"), 0.789, fail=False + ) + self.SLA = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/SLA.map"), 5, fail=False + ) + self.slope_coeff = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/slope_coeff.map"), 0.9518, fail=False + ) + self.snow_TTI = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/snow_TTI.map"), 1, fail=False + ) + self.T24_snow = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/T24_snow.map"), 18, fail=False + ) + self.Tmin = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/Tmin.map"), -10, fail=False + ) + self.Topt = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/Topt.map"), 10, fail=False + ) + self.Tgrow = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/Tgrow.map"), 200, fail=False + ) + self.Tsenc = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/Tsenc.map"), 20, fail=False + ) + self.Ud0 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/Ud0.map"), 6, fail=False + ) + self.Ug0 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/Ug0.map"), 1, fail=False + ) + self.Us0 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/Us0.map"), 6, fail=False + ) + self.Vc = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/Vc.map"), 0.5, fail=False + ) + self.w0ref_alb = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/w0ref_alb.map"), 0.3, fail=False + ) - 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)) + ds_hand = xarray.open_dataset(os.path.join(self.Dir, "staticmaps/HAND.nc")) - self.logger.info("Starting Dynamic run...") + hand = ds_hand["HAND"] + self.HAND = xarray.concat([(hand[0] * 0.0).expand_dims("z"), hand], dim="z") + perc_HAND = ds_hand["percentile"] + self.perc_HAND = xarray.concat( + [(perc_HAND[0] * 0.0).expand_dims("z"), perc_HAND], dim="z" + ) - def resume(self): - """ + psi_FC = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/psi_FC.map"), -3.3, fail=False + ) # m or hPa or 33 kPa + psi_FC0 = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/psi_FC0.map"), -0.5, fail=False + ) # m or hPa or 5 kPa - rapidly drainable theta for top soil + psi_ERRP = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/psi_ERRP.map"), -10, fail=False + ) # m or 100 kPa - assumed pressure at which soil moisture starts to limit soil evaporation (following D. Tran, 2015) + psi_d = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/psi_d.map"), -50, fail=False + ) # m assumed pressure at which soil moisture starts to limit soil water uptake + psi_PWP = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/psi_PWP.map"), -150, fail=False + ) # m + psi_res = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/psi_res.map"), -1e6, fail=False + ) # m + + theta_FC = ( + self.theta_s * min(1, (self.psi_s / psi_FC)) ** self.Lambda + ) # fraction + theta_FC0 = ( + self.theta_s * min(1, (self.psi_s / psi_FC0)) ** self.Lambda + ) # fraction + theta_ERRP = ( + self.theta_s * min(1, (self.psi_s / psi_ERRP)) ** self.Lambda + ) # fraction + theta_d = self.theta_s * min(1, (self.psi_s / psi_d)) ** self.Lambda # fraction + theta_PWP = ( + self.theta_s * min(1, (self.psi_s / psi_PWP)) ** self.Lambda + ) # fraction + theta_res = ( + self.theta_s * min(1, (self.psi_s / psi_res)) ** self.Lambda + ) # fraction + + self.S0max = ( + self.d0 * 1000 * (theta_FC0 - theta_res) + ) # mm available storage for evaporation, note FC0 is used rather than theta_sat + self.Ssmax = self.ds * 1000 * (self.theta_s - theta_PWP) + self.Sdmax = self.dd * 1000 * (self.theta_s - theta_PWP) + self.K0sat = ( + self.K0_scalar * self.k_s + ) # mm/d - note that this is technically in fact not Ksat but K(theta_FC0) + self.Kssat = ( + self.K0_scalar + * (((self.ds + self.d0) / self.d0) ** -self.Ksat_exp) + * self.k_s + ) + self.Kdsat = ( + self.K0_scalar + * (((self.dd + self.ds + self.d0) / self.d0) ** -self.Ksat_exp) + * self.k_s + ) + self.w0limE = (theta_ERRP - theta_res) / (self.theta_s - theta_res) + self.wslimU = (theta_d - theta_PWP) / (self.theta_s - theta_PWP) + self.wdlimU = (theta_d - theta_PWP) / (self.theta_s - theta_PWP) + + 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...") + + 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. """ - if self.reinit == 1: - self.logger.info("Setting initial conditions to default") - - self.Sg = cover(0) - self.Sr = cover(0) - self.Mleaf = 2./self.SLA - self.S0 = 0.2*self.w0limE*self.S0max - self.Ss = 0.2*self.wslimU*self.Ssmax - self.Sd = 0.2*self.wdlimU*self.Sdmax - self.FreeWater = cover(0) - self.DrySnow = cover(0) + if self.reinit == 1: + self.logger.info("Setting initial conditions to default") + self.Sg = cover(0) + self.Sr = cover(0) + self.Mleaf = 2. / self.SLA + self.S0 = 0.2 * self.w0limE * self.S0max + self.Ss = 0.2 * self.wslimU * self.Ssmax + self.Sd = 0.2 * self.wdlimU * self.Sdmax + self.FreeWater = cover(0) + self.DrySnow = cover(0) - else: - self.logger.info("Setting initial conditions from state files") - self.wf_resume(os.path.join(self.Dir,"instate")) - + else: + self.logger.info("Setting initial conditions from state files") + self.wf_resume(os.path.join(self.Dir, "instate")) + # for s in self.stateVariables(): + # exec "self." + s + " = cover(0)" - - - - #for s in self.stateVariables(): - # exec "self." + s + " = cover(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 @@ -347,480 +534,608 @@ """ 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 - self.WINDSPEED=cover(self.wf_readmap(self.WINDSPEED_mapstack, default=1.0), scalar(1.0)) - self.AIRPRESS=cover(self.wf_readmap(self.AIRPRESS_mapstack, default=980.0), scalar(980.0)) - # 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 + self.WINDSPEED = cover( + self.wf_readmap(self.WINDSPEED_mapstack, default=1.0), scalar(1.0) + ) + self.AIRPRESS = cover( + self.wf_readmap(self.AIRPRESS_mapstack, default=980.0), scalar(980.0) + ) + # 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, 980.0), scalar(980.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, 980.0), scalar(980.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 - pes = 610.8 * exp(17.27 * Ta/(237.3 + Ta)) + # 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 + pes = 610.8 * exp(17.27 * Ta / (237.3 + Ta)) # diagnostic equations - w0 = self.S0/self.S0max # (2.1) - ws = self.Ss/self.Ssmax # (2.1) - wd = self.Sd/self.Sdmax # (2.1) - - #Calculate vegetation parameters and cover fractions - self.LAI = self.SLA * self.Mleaf # (5.3) - fveg = max(1 - exp(- self.LAI / self.LAIref),0.000001) # (5.3) - fsoil = 1 - fveg - LUE = self.LUEmax * self.Vc * fveg - + w0 = self.S0 / self.S0max # (2.1) + ws = self.Ss / self.Ssmax # (2.1) + wd = self.Sd / self.Sdmax # (2.1) + + # Calculate vegetation parameters and cover fractions + self.LAI = self.SLA * self.Mleaf # (5.3) + fveg = max(1 - exp(-self.LAI / self.LAIref), 0.000001) # (5.3) + fsoil = 1 - fveg + LUE = self.LUEmax * self.Vc * fveg + # Calculate open water fraction - ChannelSurface = min(0,(0.007*self.Sr**0.75)) + ChannelSurface = min(0, (0.007 * self.Sr ** 0.75)) OpenWaterFrac = max(ChannelSurface, self.OpenWaterFrac) # Calculate snow cover fraction TotSnow = self.FreeWater + self.DrySnow - fsnow = min(1.0,0.05*TotSnow) # assumed; more analysis needed - + fsnow = min(1.0, 0.05 * TotSnow) # assumed; more analysis needed + # V5 'HANDometric' equations # requires self.porosity, self.HAND, self.perc_HAND - z_g = self.HAND[0] + pcr2numpy(self.Sg / (self.porosity * 1e3), np.nan) # groundwater table height in m AMSL (Sg=0 equates to drainage base) + z_g = self.HAND[0] + pcr2numpy( + self.Sg / (self.porosity * 1e3), np.nan + ) # groundwater table height in m AMSL (Sg=0 equates to drainage base) # saturated area (considers capillary rise, hence +0.3 m) - z = self.HAND[0] + pcr2numpy((self.Sg /(self.porosity * 1e3) + (-self.psi_s)), np.NaN) #bubbling pressure as indication of capillary fringe - fg = interp_hand(z,self.HAND,self.perc_HAND) / 100.0 + z = self.HAND[0] + pcr2numpy( + (self.Sg / (self.porosity * 1e3) + (-self.psi_s)), np.NaN + ) # bubbling pressure as indication of capillary fringe + fg = interp_hand(z, self.HAND, self.perc_HAND) / 100.0 # same for veg with access to gw - RD = 1.0 # assumed maximum depth of shallow root water uptake - z = self.HAND[0] + pcr2numpy((self.Sg / (self.porosity * 1e3)),np.nan) + RD - fUgShallow = (interp_hand(z,self.HAND,self.perc_HAND) /100.0) * (1.0 - self.fPotDeep ) - RD = 7.0 # assumed maximum depth of deep root water uptake - z = self.HAND[0] + pcr2numpy(self.Sg /(self.porosity * 1e3),np.nan) + RD - fUgDeep = interp_hand(z,self.HAND,self.perc_HAND) /100 * self.fPotDeep + RD = 1.0 # assumed maximum depth of shallow root water uptake + z = self.HAND[0] + pcr2numpy((self.Sg / (self.porosity * 1e3)), np.nan) + RD + fUgShallow = (interp_hand(z, self.HAND, self.perc_HAND) / 100.0) * ( + 1.0 - self.fPotDeep + ) + RD = 7.0 # assumed maximum depth of deep root water uptake + z = self.HAND[0] + pcr2numpy(self.Sg / (self.porosity * 1e3), np.nan) + RD + fUgDeep = interp_hand(z, self.HAND, self.perc_HAND) / 100 * self.fPotDeep fUg = fUgShallow + fUgDeep - - + # Spatialise these fractions (largely superfluous with 1 HRU) # Rewrite this part if > 1 HRU fw_local = ChannelSurface fwater = OpenWaterFrac - fsat = min(1,max( OpenWaterFrac, fg )) ## V5 - + fsat = min(1, max(OpenWaterFrac, fg)) ## V5 + # Aerodynamic conductance (3.7) - fh = ln(813/max(0.25,self.hveg)-5.45) # assume minimum roughness of 0.25 m + fh = ln(813 / max(0.25, self.hveg) - 5.45) # assume minimum roughness of 0.25 m # ADJUSTED FOR E2O WFEI DATA: uz at 1m screen height (see AWRA technical report) - ku1 = 0.359 / (fh*(fh+2.3)) - ga = max(0.001, ku1*self.u1) # minimum of 0.001 imposed to avoid issues + ku1 = 0.359 / (fh * (fh + 2.3)) + ga = max(0.001, ku1 * self.u1) # minimum of 0.001 imposed to avoid issues - - 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! self.Ept = self.E0 - + elif self.UseETPdata == 0: # CALCULATION OF PET # Conversions and coefficients (3.1) - fRH = pe/pes # relative air humidity -------------- check - cRE = 0.03449+4.27e-5*Ta - Caero = 0.176*(1+Ta/209.1)*(pair-0.417*pe)*(1-fRH) # removed fday as already daytime - keps = 1.4e-3*((Ta/187)**2+Ta/107+1)*(6.36*pair+pe)/pes - Rgeff = Rg /self.fday # this is original + fRH = ( + pe / pes + ) # relative air humidity -------------- check + cRE = 0.03449 + 4.27e-5 * Ta + Caero = ( + 0.176 * (1 + Ta / 209.1) * (pair - 0.417 * pe) * (1 - fRH) + ) # removed fday as already daytime + keps = 1.4e-3 * ((Ta / 187) ** 2 + Ta / 107 + 1) * (6.36 * pair + pe) / pes + Rgeff = Rg / self.fday # this is original - # albedo model + # albedo model alb_veg = self.f_alb_Vc * self.Vc - dryfrac = exp(-w0/self.w0ref_alb)*(1-fsat) - alb_soil = self.alb_wet + (self.alb_dry-self.alb_snow) * dryfrac - alb_ns = fveg*alb_veg+fsoil*alb_soil - alb = (1-fwater)*(1-fsnow)*alb_ns +fsnow*self.alb_snow + fwater*self.alb_water - - RSn = (1-alb)*Rgeff - + dryfrac = exp(-w0 / self.w0ref_alb) * (1 - fsat) + alb_soil = self.alb_wet + (self.alb_dry - self.alb_snow) * dryfrac + alb_ns = fveg * alb_veg + fsoil * alb_soil + alb = ( + (1 - fwater) * (1 - fsnow) * alb_ns + + fsnow * self.alb_snow + + fwater * self.alb_water + ) + + RSn = (1 - alb) * Rgeff + # long wave radiation balance (3.3 to 3.5) StefBolz = 5.67e-8 - Tkelv = Ta+273.16 - - RLin = self.LWdown # provided by E2O data (though not sure how good it is..) - RLout = 1 * StefBolz * Tkelv**4 #v0.5 # (3.4) - RLn = RLin-RLout - - self.fGR = self.Gfrac_max*(1-exp(-fsoil/self.fvegref_G)) - self.Rneff = max(1, (RSn+self.RLn)*(1-self.fGR) ) # original (assuming any condensation is already measured in rain and there is a minimum Rneff of 1 W/m2 (to prevent any zero issues) - - - - + Tkelv = Ta + 273.16 + + RLin = ( + self.LWdown + ) # provided by E2O data (though not sure how good it is..) + RLout = 1 * StefBolz * Tkelv ** 4 # v0.5 # (3.4) + RLn = RLin - RLout + + self.fGR = self.Gfrac_max * (1 - exp(-fsoil / self.fvegref_G)) + self.Rneff = max( + 1, (RSn + self.RLn) * (1 - self.fGR) + ) # original (assuming any condensation is already measured in rain and there is a minimum Rneff of 1 W/m2 (to prevent any zero issues) + # Potential evaporation (original) - kalpha = min(1.4, 1+Caero*ga/self.Rneff) # do not allow value higher as that implies a unlikely high rate of advection from nearby areas only likely to occur for wet canopy. - self.E0 = cRE*(1/(1+keps))*kalpha*self.Rneff*self.fday # for canopy - self.Ept = cRE*(1/(1+keps))*1.26*self.Rneff*self.fday # for open water - + kalpha = min( + 1.4, 1 + Caero * ga / self.Rneff + ) # do not allow value higher as that implies a unlikely high rate of advection from nearby areas only likely to occur for wet canopy. + self.E0 = ( + cRE * (1 / (1 + keps)) * kalpha * self.Rneff * self.fday + ) # for canopy + self.Ept = ( + cRE * (1 / (1 + keps)) * 1.26 * self.Rneff * self.fday + ) # for open water + # CALCULATION OF ET FLUXES AND ROOT WATER UPTAKE # Root water uptake constraint (4.4) - # For v5 no Uomax so temporarily bypassed here + # For v5 no Uomax so temporarily bypassed here 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 - + Gsmax = self.Gs_scalar * self.cGsmax * self.Vc + if self.UseETPdata == 1: fD = 1.0 - elif self.UseETPdata == 0: - 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 + elif self.UseETPdata == 0: + 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) # # Below seems to be in v5 U0 = scalar(0) - Us = max(0, min( (Usmax /(Usmax+Udmax+Ugmax))*Et, self.Ss-1e-2 ) ) - Ud = max(0, min( (Udmax /(Usmax+Udmax+Ugmax))*Et, self.Sd-1e-2 ) ) - Ug = max(0, min( (Ugmax /(Usmax+Udmax+Ugmax))*Et, self.Sd-1e-2 ) ) - - - Et = U0 + Us + Ud + Ug # to ensure mass balance + Us = max(0, min((Usmax / (Usmax + Udmax + Ugmax)) * Et, self.Ss - 1e-2)) + Ud = max(0, min((Udmax / (Usmax + Udmax + Ugmax)) * Et, self.Sd - 1e-2)) + Ug = max(0, min((Ugmax / (Usmax + Udmax + Ugmax)) * Et, self.Sd - 1e-2)) + Et = U0 + Us + Ud + Ug # to ensure mass balance + # Soil evaporation (4.5) - w0x = max(0, (self.S0 - U0)/self.S0max) # adjusted top soil water content - fsoilE = self.FsoilEmax*min(1,w0x/self.w0limE) - Es0 = (1-fsat)*fsoilE*(max(0,self.E0-Et)) - + w0x = max(0, (self.S0 - U0) / self.S0max) # adjusted top soil water content + fsoilE = self.FsoilEmax * min(1, w0x / self.w0limE) + Es0 = (1 - fsat) * fsoilE * (max(0, self.E0 - Et)) + # Groundwater evaporation (4.6) - Eg0 = max(0,fsat-fwater)*self.FsoilEmax*max(0,self.E0-Et) + Eg0 = max(0, fsat - fwater) * self.FsoilEmax * max(0, self.E0 - Et) Es = Es0 + Eg0 - + # Open water evaporation (4.7) # uses Priestley-Taylor - Erl = fw_local*self.FwaterE*self.Ept # from local river channels - Err = (fwater-fw_local)*self.FwaterE*self.Ept # from remaining open water + Erl = fw_local * self.FwaterE * self.Ept # from local river channels + Err = (fwater - fw_local) * self.FwaterE * self.Ept # from remaining open water 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(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(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-(self.snow_TT-self.snow_TTI/2))/self.snow_TTI,1)) - SnowFrac = 1 - RainFrac #fraction of precipitation which falls as snow + Temperature = T24 # Dimmie, let op: tijdelijke regel!! + RainFrac = max( + 0, + min((Temperature - (self.snow_TT - self.snow_TTI / 2)) / self.snow_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 = self.snow_Cfmax*max(0,Temperature-self.snow_TT) # Potential snow melt, based on temperature - PotRefreezing = self.snow_Cfmax*self.snow_Cfr*max(self.snow_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 + SnowFall = SnowFrac * Pn # snowfall depth + RainFall = RainFrac * Pn # rainfall depth + PotSnowMelt = self.snow_Cfmax * max( + 0, Temperature - self.snow_TT + ) # Potential snow melt, based on temperature + PotRefreezing = ( + self.snow_Cfmax * self.snow_Cfr * max(self.snow_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.FreeWater * self.snow_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 ## v5 ## - 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 ## v5 ## + 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 # combined runoff + QR = Rhof + Rsof + Rmelt # combined runoff I = Ps - Rhof - Rsof # SOIL WATER BALANCES (2.1 & 2.4) # Soil hydrology from v5 (Viney et al., 2015) http://www.bom.gov.au/water/landscape/static/publications/Viney_et_al_2015_AWRA_L_5.0_model_description.pdf - Kr_0s = self.K0sat/self.Kssat - Rh_0s = pcr_tanh(self.slope_coeff*self.slope*w0)*pcr_tanh(self.Kr_coeff*(Kr_0s-1.0)*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.0) * w0 + ) # general 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*((S0/self.S0max)**2) - IF0 = Rh_0s*Km*((S0/self.S0max)**2) + S0 = (-B + ((B ** 2 - 4 * A * C) ** 0.5)) / (2 * A) + D0 = (1 - Rh_0s) * Km * ((S0 / self.S0max) ** 2) + IF0 = Rh_0s * Km * ((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 = ifthenelse(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 = ifthenelse(imap, self.S0max, S0) # enforce mass balance (there can be small numerical errors in quadratic equation) - 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 root zone water balance (Ss) (2.4) - 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 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 = (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 = (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 = ifthenelse(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 = ifthenelse(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 root zone water balance (Sd) (2.4) # general case - A = self.Kdsat/(self.Sdmax**2) + A = self.Kdsat / (self.Sdmax ** 2) B = 1.0 - 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 = (Sd+Ds)<=Ud - Ud = ifthenelse(imap,(self.Sd+Ds),Ud) - Sd = ifthenelse(imap,0,Sd) - Dd = ifthenelse(imap,0,Dd) + imap = (Sd + Ds) <= Ud + Ud = ifthenelse(imap, (self.Sd + Ds), Ud) + Sd = ifthenelse(imap, 0, Sd) + Dd = ifthenelse(imap, 0, Dd) # 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 = ifthenelse(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 = ifthenelse(imap, self.Sdmax, Sd) # enforce mass balance (for numerical & rounding errors - Sd = max(0, min(Sd,self.Sdmax)) - massbal = self.Sd + Ds - Ud - Dd - IFd - Sd + 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 # add to runoff - + IFs = IFs + IFd # add up to interflow + QR = QR + IF0 + IFs # add to runoff + # CATCHMENT WATER BALANCE # Groundwater store water balance (Sg) (2.5) NetGf = Dd - Eg0 - Ug self.Sg = self.Sg + NetGf - Sg_fd = max(self.Sg,0) - Qg = min(Sg_fd, (1-exp(-self.K_gw))*Sg_fd) + Sg_fd = max(self.Sg, 0) + Qg = min(Sg_fd, (1 - exp(-self.K_gw)) * Sg_fd) self.Sg = self.Sg - Qg # Surface water store water balance (Sr) (2.7) - self.Sr = max(0, self.Sr + QR - Erl + 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 - Erl + Qg) + self.Qtot = max(0, min(self.Sr, (1 - exp(-self.K_rout)) * self.Sr)) + self.Sr = self.Sr - self.Qtot # VEGETATION ADJUSTMENT (5.7-5.8) - #how to deal with self.LAImax (not set now)? - #self.LAImax = 8.0 + # how to deal with self.LAImax (not set now)? + # self.LAImax = 8.0 - 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) - + 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) + # VEGETATION ADJUSTMENT (5.4-5.6) - dMleaf = -ln(1-fveq)*self.LAIref/self.SLA-self.Mleaf - Mleafnet = scalar(dMleaf>0)*(dMleaf/self.Tgrow) +scalar(dMleaf<0)*dMleaf/self.Tsenc + dMleaf = -ln(1 - fveq) * self.LAIref / self.SLA - self.Mleaf + 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.LAI/self.LAIref) #(5.3) + fveg = 1 - exp(-self.LAI / self.LAIref) # (5.3) # in case this is desired as output: - self.w0 = self.S0/self.S0max #(2.1) - self.TotSnow = self.DrySnow + self.FreeWater + self.w0 = self.S0 / self.S0max # (2.1) + 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_w3" # "D:/trambaue/_Projects/GLOFFIS/201501/GLOFFIS_SA/Modules/openstreams_w3ra/" + caseName = ( + "openstreams_w3" + ) # "D:/trambaue/_Projects/GLOFFIS/201501/GLOFFIS_SA/Modules/openstreams_w3ra/" runId = "run_default" - configfile="wflow_w3.ini" - _lastTimeStep = 0 - _firstTimeStep = 0 - timestepsecs=86400 + configfile = "wflow_w3.ini" + _lastTimeStep = 0 + _firstTimeStep = 0 + timestepsecs = 86400 - 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:') - + opts, args = getopt.getopt(argv, "C:S:T:c:s:R:") + for o, a in opts: - if o == '-C': caseName = a - if o == '-R': runId = a - if o == '-c': configfile = a - if o == '-s': timestepsecs = int(a) - - starttime = dt.datetime(1990,01,01) + if o == "-C": + caseName = a + if o == "-R": + runId = a + if o == "-c": + configfile = a + if o == "-s": + timestepsecs = int(a) + starttime = dt.datetime(1990, 01, 01) + 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_w3",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_w3", + 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, 'run', '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 == '-T': - configset(myModel.config, 'run', 'endtime', a, overwrite=True) - if o == '-S': - configset(myModel.config, 'run', 'starttime', 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, "run", "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 == "-T": + configset(myModel.config, "run", "endtime", a, overwrite=True) + if o == "-S": + configset(myModel.config, "run", "starttime", a, overwrite=True) dynModelFw.setupFramework() dynModelFw._runInitial() dynModelFw._runResume() - #dynModelFw._runDynamic(0,0) + # dynModelFw._runDynamic(0,0) dynModelFw._runDynamic(_firstTimeStep, _lastTimeStep) dynModelFw._runSuspend() dynModelFw._wf_shutdown() - + if __name__ == "__main__": - main() \ No newline at end of file + main()