Index: examples/batches/run_topoflex.bat =================================================================== diff -u -r687a7cb59a549f3d4f41cd8b76eac92b1f8486d1 -rfacc3a7f729f6910dac37d4761992c70d238fcbd --- examples/batches/run_topoflex.bat (.../run_topoflex.bat) (revision 687a7cb59a549f3d4f41cd8b76eac92b1f8486d1) +++ examples/batches/run_topoflex.bat (.../run_topoflex.bat) (revision facc3a7f729f6910dac37d4761992c70d238fcbd) @@ -1,3 +1,3 @@ cd -python ../../wflow-py/Sandbox/wflow_topoflex/wflow_topoflex.py -C ../wflow_orientale_topoflex -R test_run_reservoirs_names -c wflow_orientale_topoflex_zero_res.ini -I -T 20 +python ../../wflow-py/wflow/wflow_topoflex.py -C ../wflow_orientale_topoflex -R test_run_reservoirs_names -c wflow_orientale_topoflex_zero_res.ini -I -T 20 pause \ No newline at end of file Index: examples/wflow_orientale_topoflex/wflow_orientale_topoflex_zero_res.ini =================================================================== diff -u -r687a7cb59a549f3d4f41cd8b76eac92b1f8486d1 -rfacc3a7f729f6910dac37d4761992c70d238fcbd --- examples/wflow_orientale_topoflex/wflow_orientale_topoflex_zero_res.ini (.../wflow_orientale_topoflex_zero_res.ini) (revision 687a7cb59a549f3d4f41cd8b76eac92b1f8486d1) +++ examples/wflow_orientale_topoflex/wflow_orientale_topoflex_zero_res.ini (.../wflow_orientale_topoflex_zero_res.ini) (revision facc3a7f729f6910dac37d4761992c70d238fcbd) @@ -19,7 +19,6 @@ maxTransitTime = 9 DistForcing = 3 maxGaugeId = 10 -Ks = 0.0004 #spinUp_time = 5 #NSEbreak = 0 @@ -28,58 +27,61 @@ timestepsecs = 3600 #selection of reservoir configuration - 'None' means reservoir should not be modelled for this class +selectSw = None, None selectSi = interception_overflow2, interception_overflow2 -selectSu= unsatZone_LP_beta, unsatZone_LP_beta -selectSus=None,None -selectSf=fastRunoff_lag2, fastRunoff_lag2 -selectSr=None,None +selectSa = None, None +selectSu = unsatZone_LP_beta, unsatZone_LP_beta +selectSf = fastRunoff_lag2, fastRunoff_lag2 +selectSfa = None, None +selectSr = None, None + # selection of Ss (lumped over entire catchment, so only one!) # selectSs = -#input time series -Pfile_1 = intss\1_P.tss -Efile_1 = intss\1_PET.tss -Tfile_1 = intss\1_T.tss -#Qfile_1 = D:/TEuser/Onderzoek/Promotie/modellen/OpenStreams/wflow/wflow_orientale_topoflex/intss/1_Qobs.tss -#Pfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_rain_10H.tss -#Efile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_Ep_10H.tss -#TDMfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_tempDMean_10H.tss -#RNfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_radNet_10H.tss -#RSfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_radSI_10H.tss -#SGfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_sgamma_10H.tss -#VPDfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_vpd_10H.tss -#Wfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_wind_10H.tss -#DSfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_dayS_10H.tss -#DEfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_dayE_10H.tss -#LAIfile_2 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_lai_10H.tss -#rst_lai_0 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_rstMin_laiEff_W4.tss -#rst_lai_1 = /u/euser/WFLOW/wflow/wflow_ourthe_testStructure/intss/2_rstMin_laiEff_HPPPA4.tss +# selection of routing routine (for combined fluxes of different classes, so only one!) +selectRout = routingQf_combined +#selectRout = routingQf_Qs_grid #wflow maps with percentages wflow_percent_0 = staticmaps/wflow_percentW4.map wflow_percent_1 = staticmaps/wflow_percentHPPPA4.map #constant model parameters - some are catchment dependent -Ks = 0.0004 +Ks = [0.0004, 0.0004] lamda = 2.45e6 -Cp = 1.01e-3 -rhoA = 1.29 +lamdaS = 2.83e6 +Cp = 1004 +rhoA = 1.205 rhoW = 1000 gamma = 0.066 JC_Topt = 301 #parameters for fluxes and storages -sumax = [140, 300] -beta = [0.1, 0.1] -D = [0, 0.10] -Kf = [0.1, 0.006] +sumax = [70,85] +samax = [0,0] +beta = [0.5, 0.5] +betaA = [0,0.22] +D = [0, 0.26] +Kf = [0.055, 0.0045] +Kfa = [0,0.055] Tf = [1, 3] -imax = [1.2, 2] -perc = [0, 0.000] +imax = [1.2, 2.2] +perc = [0, 0.008] cap = [0.09, 0] -LP = [0.5, 0.8] -Ce = [1, 1] +LP = [0.05,0.15] +redsu = [1,1] +Tfa = [0,0] +decF = [0,0] +Fmax = [0,0] +Fmin = [0,0] +Tt = [0.7, 0.7] +Tm = [1.4, 1.4] +Fm = [0.1, 0.1] +FrDur0 = [-1, -1, -1] +FrDur1 = [0, 0, 0] +dayDeg = [1, 1, 1] +ratFT = [5, 5, 5] #Jarvis stressfunctions JC_D05 = [1.5,1.5] Index: wflow-py/wflow/JarvisCoefficients.py =================================================================== diff -u -r49be3f4c4953a2a77f98e4ca62c6c3d118922167 -rfacc3a7f729f6910dac37d4761992c70d238fcbd --- wflow-py/wflow/JarvisCoefficients.py (.../JarvisCoefficients.py) (revision 49be3f4c4953a2a77f98e4ca62c6c3d118922167) +++ wflow-py/wflow/JarvisCoefficients.py (.../JarvisCoefficients.py) (revision facc3a7f729f6910dac37d4761992c70d238fcbd) @@ -35,6 +35,25 @@ # potential_evaporation(self,k) downscale_evaporation(self,k) +def calcEpSnow(self,k): + """ + REQUIRED INPUT: + daily input data + """ +# resistenceAeroD(self) +# potential_evaporation(self,k) + self.EpDaySnow = self.EpDay + downscale_evaporation_snow(self,k) + +def calcEpSnowHour(self,k): + """ + REQUIRED INPUT: + hourly input data + """ +# resistenceAeroD(self) +# potential_evaporation(self,k) + self.EpHour = self.PotEvaporation / self.lamda * self.lamdaS + def calcEu_laiFixed(self,k): """ REQUIRED INPUT: @@ -208,25 +227,49 @@ """ REQUIRED INPUT: daily evaporation (EpDay) - - hour of the day (x; derived from self.teller) + - hour of the day (x; derived from self.thestep) - start of the day (derived from global radiation) - end of the day (derived from global radiation) PARAMETERS: - """ - - x = self.teller - floor(self.teller/24) * 24 * scalar(self.TopoId) + + teller = numpy.max(pcr2numpy(self.thestep,nan)) + x = teller - floor(teller/24) * 24 * scalar(self.TopoId) DL = self.DE - self.DS + 1 P = 2 * pi / (2 * DL) # period SH = DL - 12 #horizontal shift of new signal aN = -1 * self.EpDay * P # nominator of the amplitude aDN = sin((P * (self.DE + SH)) * 180 / pi) - sin((P * (self.DS + SH)) * 180 / pi) # denominator of the amplitude ampl = aN / aDN # amplitude of new signal - self.EpHour = ifthenelse(pcrand(x >= self.DS, x <= self.DE), -1 * ampl * cos((P * (x + SH)) * 180 / pi), 0) + self.EpHour = max(ifthenelse(pcrand(x >= self.DS, x <= self.DE), -1 * ampl * cos((P * (x + SH)) * 180 / pi), 0), 0) + +def downscale_evaporation_snow(self,k): + """ + REQUIRED INPUT: + daily evaporation (EpDay) + - hour of the day (x; derived from self.teller) + - start of the day (derived from global radiation) + - end of the day (derived from global radiation) + PARAMETERS: + - + """ + teller = numpy.max(pcr2numpy(self.thestep,nan)) + x = teller - floor(teller/24) * 24 * scalar(self.TopoId) + DL = self.DE - self.DS + 1 + P = 2 * pi / (2 * DL) # period + SH = DL - 12 #horizontal shift of new signal + aN = -1 * self.EpDaySnow * P # nominator of the amplitude + aDN = sin((P * (self.DE + SH)) * 180 / pi) - sin((P * (self.DS + SH)) * 180 / pi) # denominator of the amplitude + ampl = aN / aDN # amplitude of new signal + self.EpHour = max(ifthenelse(pcrand(x >= self.DS, x <= self.DE), -1 * ampl * cos((P * (x + SH)) * 180 / pi), 0), 0) + + + \ No newline at end of file Index: wflow-py/wflow/reservoir_Sa.py =================================================================== diff -u -rc2a3a7569feb954f3646d329d41023cf0793b2ad -rfacc3a7f729f6910dac37d4761992c70d238fcbd --- wflow-py/wflow/reservoir_Sa.py (.../reservoir_Sa.py) (revision c2a3a7569feb954f3646d329d41023cf0793b2ad) +++ wflow-py/wflow/reservoir_Sa.py (.../reservoir_Sa.py) (revision facc3a7f729f6910dac37d4761992c70d238fcbd) @@ -38,9 +38,37 @@ name = 'agriZone_Ep_Sa' elif i == 4: name = 'agriZone_Ep_Sa_cropG' + elif i == 5: + name = 'agriZone_Ep_Sa_cropG_beta' + elif i == 6: + name = 'agriZone_Ep_Sa_beta' + elif i == 7: + name = 'agriZone_Ep_Sa_beta_frost' + elif i == 8: + name = 'agriZone_Ep_Sa_beta_Fvar' + elif i == 9: + name = 'agriZone_hourlyEp_Sa_beta_Fvar' + elif i == 10: + name = 'agriZone_hourlyEp_Sa_beta_frost' + elif i == 11: + name = 'agriZone_hourlyEp_Sa_beta_frostSamax' return name - +def agriZone_no_reservoir(self, k): + """ + This function is used when no unsaturated zone reservoir is used and only + passes fluxes from the upper reservoirs to the lower + self.Qa_[k] = 0. + self.Ea_[k] = 0. + self.Sa[k] = 0. + self.Fa_[k] = Pe + Storage in unsaturated zone = 0. + """ + self.Qa_[k] = 0. + self.Ea_[k] = 0. + self.Sa[k] = 0. + self.Fa_[k] = max(self.Pe_[k], 0) + self.wbSa_[k] = self.Pe_[k] - self.Ea_[k] - self.Qa_[k] - self.Fa_[k] - self.Sa[k] + self.Sa_t[k] def agriZone_Jarvis(self,k): """ @@ -60,8 +88,6 @@ JarvisCoefficients.calcEu(self,k,1) #calculation of Ea based on Jarvis stress functions self.Ea1 = self.Eu -# if self.teller == 45: -# pdb.set_trace() self.Fa1 = self.Fmin[k] + (self.Fmax[k] - self.Fmin[k]) * e ** (-self.decF[k] * self.SuN) self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa) - self.Fa1 - self.Ea1 @@ -170,9 +196,6 @@ self.samax2 = self.samax[k] * self.cropG self.Qaadd = max(self.Sa_t[k] - self.samax2,0) -# if self.teller == 66: -# pdb.set_trace() - self.Qa = max(self.Pe - (self.samax2 - self.Sa_t[k]),0) + self.Qaadd self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa) self.SaN = self.Sa[k] / self.samax2 @@ -195,3 +218,357 @@ self.Ea_[k] = self.Ea self.Qa_[k] = self.Qa self.Fa_[k] = self.Fa + +def agriZone_Ep_Sa_cropG_beta(self,k): + """ + - Potential evaporation is decreased by energy used for interception evaporation + - Formula for evaporation based on LP + - Outgoing fluxes are determined based on (value in previous timestep + inflow) + and if this leads to negative storage, the outgoing fluxes are corrected to rato --> Eu is + no longer taken into account for this correction + - Qa u is determined from overflow from Sa --> incorporation of beta function + - Fa is based on storage in Sa + - Code for ini-file: 5 + """ + + JarvisCoefficients.calcEp(self,k) + self.PotEvaporation = self.EpHour + + self.samax2 = self.samax[k] * self.cropG + self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0) + + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) + self.SaN = self.Sa[k] / self.samax2 + self.SuN = self.Su[k] / self.sumax[k] + + self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1) + self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k]) + self.Fa1 = ifthenelse(self.SaN > 0,self.Fmin[k] + (self.Fmax[k] - self.Fmin[k]) * e ** (-self.decF[k] * (1 - self.SaN)),0) + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Qa1 - self.Fa1 - self.Ea1 + + self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0) + self.Qa = self.Qa1 + (self.Qa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1, 1)) * self.Sa_diff + self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff + self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Ea - self.Fa - self.Qa + self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k]) + self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) + + self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Qaadd - self.Fa - self.Sa[k] + self.Sa_t[k] + + self.Ea_[k] = self.Ea + self.Qa_[k] = self.Qa + self.Qaadd + self.Fa_[k] = self.Fa + + +def agriZone_Ep_Sa_beta(self,k): + """ + - Potential evaporation is decreased by energy used for interception evaporation + - Formula for evaporation based on LP + - Outgoing fluxes are determined based on (value in previous timestep + inflow) + and if this leads to negative storage, the outgoing fluxes are corrected to rato --> Eu is + no longer taken into account for this correction + - Qa u is determined from overflow from Sa --> incorporation of beta function + - Fa is based on storage in Sa + - Code for ini-file: 6 + """ + + JarvisCoefficients.calcEp(self,k) + self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0) + + self.samax2 = self.samax[k] * scalar(self.TopoId) + self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0) + + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) + self.SaN = self.Sa[k] / self.samax2 + self.SuN = self.Su[k] / self.sumax[k] + + self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1) + self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k]) + self.Fa1 = ifthenelse(self.SaN > 0,self.Fmin[k] + (self.Fmax[k] - self.Fmin[k]) * e ** (-self.decF[k] * (1 - self.SaN)),0) + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Qa1 - self.Fa1 - self.Ea1 + + self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0) + self.Qa = self.Qa1 + (self.Qa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff + self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff + self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Ea - self.Fa - self.Qa + self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k]) + self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) + + self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Qaadd - self.Fa - self.Sa[k] + self.Sa_t[k] + + self.Ea_[k] = self.Ea + self.Qa_[k] = self.Qa + self.Qaadd + self.Fa_[k] = self.Fa + +def agriZone_Ep_Sa_beta_frost(self,k): + """ + - Potential evaporation is decreased by energy used for interception evaporation + - Formula for evaporation based on LP + - Outgoing fluxes are determined based on (value in previous timestep + inflow) + and if this leads to negative storage, the outgoing fluxes are corrected to rato --> Eu is + no longer taken into account for this correction + - Qa u is determined from overflow from Sa --> incorporation of beta function + - Fa is based on storage in Sa + - Fa is decreased in case of frozen soil + - Code for ini-file: 7 + """ + + JarvisCoefficients.calcEp(self,k) + self.PotEvaporation = self.EpHour + + self.samax2 = self.samax[k] * scalar(self.TopoId) + self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0) + self.FrDur[k] = min(self.FrDur[k] + (self.Tmean - 273.15) / 86400 * self.timestepsecs * self.dayDeg[k], 0) + + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) + self.SaN = self.Sa[k] / self.samax2 + self.SuN = self.Su[k] / self.sumax[k] + + self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1) + self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k]) + self.Ft = min(max(self.FrDur[k] / (self.FrDur1[k] - self.FrDur0[k]) - self.FrDur0[k] / (self.FrDur1[k] - self.FrDur0[k]), 0), 1) + self.Fa1 = self.Ft * ifthenelse(self.SaN > 0,self.Fmin[k] + (self.Fmax[k] - self.Fmin[k]) * e ** (-self.decF[k] * (1 - self.SaN)),0) + + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Qa1 - self.Fa1 - self.Ea1 + + self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0) + self.Qa = self.Qa1 + (self.Qa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0, self.Fa1 + self.Ea1 + self.Qa1, 1)) * self.Sa_diff + self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0, self.Fa1 + self.Ea1 + self.Qa1, 1)) * self.Sa_diff + self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0, self.Fa1 + self.Ea1 + self.Qa1, 1)) * self.Sa_diff + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Ea - self.Fa - self.Qa + self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k]) + self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) + + self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Qaadd - self.Fa - self.Sa[k] + self.Sa_t[k] + + self.Ea_[k] = self.Ea + self.Qa_[k] = self.Qa + self.Qaadd + self.Fa_[k] = self.Fa + self.Ft_[k] = self.Ft + +def agriZone_hourlyEp_Sa_beta_frost(self,k): + """ + - Potential evaporation is decreased by energy used for interception evaporation + - Formula for evaporation based on LP + - Outgoing fluxes are determined based on (value in previous timestep + inflow) + and if this leads to negative storage, the outgoing fluxes are corrected to rato --> Eu is + no longer taken into account for this correction + - Qa u is determined from overflow from Sa --> incorporation of beta function + - Fa is based on storage in Sa + - Fa is decreased in case of frozen soil + - Code for ini-file: 10 + """ + + #JarvisCoefficients.calcEp(self,k) + #self.PotEvaporation = self.EpHour + + self.samax2 = self.samax[k] * scalar(self.TopoId) + self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0) + self.FrDur[k] = min(self.FrDur[k] + (self.Temperature) / 86400 * self.timestepsecs * self.dayDeg[k], 0) + + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) + self.SaN = self.Sa[k] / self.samax2 + self.SuN = self.Su[k] / self.sumax[k] + + self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1) + self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k]) + self.Ft = min(max(self.FrDur[k] / (self.FrDur1[k] - self.FrDur0[k]) - self.FrDur0[k] / (self.FrDur1[k] - self.FrDur0[k]), 0), 1) + self.Fa1 = self.Ft * ifthenelse(self.SaN > 0,self.Fmin[k] + (self.Fmax[k] - self.Fmin[k]) * e ** (-self.decF[k] * (1 - self.SaN)),0) + + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Qa1 - self.Fa1 - self.Ea1 + + self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0) + self.Qa = self.Qa1 + (self.Qa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0, self.Fa1 + self.Ea1 + self.Qa1, 1)) * self.Sa_diff + self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0, self.Fa1 + self.Ea1 + self.Qa1, 1)) * self.Sa_diff + self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0, self.Fa1 + self.Ea1 + self.Qa1, 1)) * self.Sa_diff + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Ea - self.Fa - self.Qa + self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k]) + self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) + + self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Qaadd - self.Fa - self.Sa[k] + self.Sa_t[k] + + self.Ea_[k] = self.Ea + self.Qa_[k] = self.Qa + self.Qaadd + self.Fa_[k] = self.Fa + self.Ft_[k] = self.Ft + +def agriZone_hourlyEp_Sa_beta_frostSamax(self,k): + """ + - Potential evaporation is decreased by energy used for interception evaporation + - Formula for evaporation based on LP + - Outgoing fluxes are determined based on (value in previous timestep + inflow) + and if this leads to negative storage, the outgoing fluxes are corrected to rato --> Eu is + no longer taken into account for this correction + - Qa u is determined from overflow from Sa --> incorporation of beta function + - Fa is based on storage in Sa + - Fa is decreased in case of frozen soil + - Code for ini-file: 11 + """ + + #JarvisCoefficients.calcEp(self,k) + #self.PotEvaporation = self.EpHour + + self.FrDur[k] = min(self.FrDur[k] + (self.Temperature) * self.dayDeg[k], 0) + self.Ft = min(max(self.FrDur[k] / (self.FrDur1[k] - self.FrDur0[k]) - self.FrDur0[k] / (self.FrDur1[k] - self.FrDur0[k]), 0.1), 1) + + + self.samax2 = self.samax[k] * scalar(self.TopoId) * self.Ft + self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0) + + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) + self.SaN = self.Sa[k] / self.samax2 + self.SuN = self.Su[k] / self.sumax[k] + + self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1) + self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k]) + self.Fa1 = ifthenelse(self.SaN > 0,self.Fmin[k] + (self.Fmax[k] - self.Fmin[k]) * e ** (-self.decF[k] * (1 - self.SaN)),0) + + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Qa1 - self.Fa1 - self.Ea1 + + self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0) + self.Qa = self.Qa1 + (self.Qa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0, self.Fa1 + self.Ea1 + self.Qa1, 1)) * self.Sa_diff + self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0, self.Fa1 + self.Ea1 + self.Qa1, 1)) * self.Sa_diff + self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0, self.Fa1 + self.Ea1 + self.Qa1, 1)) * self.Sa_diff + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Ea - self.Fa - self.Qa + self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k]) + self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) + + self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Qaadd - self.Fa - self.Sa[k] + self.Sa_t[k] + + self.Ea_[k] = self.Ea + self.Qa_[k] = self.Qa + self.Qaadd + self.Fa_[k] = self.Fa + self.Ft_[k] = self.Ft + +def agriZone_Ep_Sa_beta_frostSamax(self,k): + """ + - Potential evaporation is decreased by energy used for interception evaporation + - Formula for evaporation based on LP + - Outgoing fluxes are determined based on (value in previous timestep + inflow) + and if this leads to negative storage, the outgoing fluxes are corrected to rato --> Eu is + no longer taken into account for this correction + - Qa u is determined from overflow from Sa --> incorporation of beta function + - Fa is based on storage in Sa + - Fa is decreased in case of frozen soil + - Code for ini-file: 12 + """ + + JarvisCoefficients.calcEp(self,k) + self.PotEvaporation = self.EpHour + + self.FrDur[k] = min(self.FrDur[k] + ifthenelse(self.Temperature > 0, self.ratFT[k] * self.dayDeg[k] * self.Temperature, self.Temperature) * self.dayDeg[k], 0) + self.Ft = min(max(self.FrDur[k] / (self.FrDur1[k] - self.FrDur0[k]) - self.FrDur0[k] / (self.FrDur1[k] - self.FrDur0[k]), 0.05), 1) + + + self.samax2 = self.samax[k] * scalar(self.TopoId) * self.Ft + self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0) + + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) + self.SaN = self.Sa[k] / self.samax2 + self.SuN = self.Su[k] / self.sumax[k] + + self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1) + self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k]) + + self.Fa1 = ifthenelse(self.SaN > 0,self.Fmin[k] + (self.Fmax[k] - self.Fmin[k]) * e ** (-self.decF[k] * (1 - self.SaN)),0) + + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Qa1 - self.Fa1 - self.Ea1 + + self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0) + self.Qa = self.Qa1 + (self.Qa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0, self.Fa1 + self.Ea1 + self.Qa1, 1)) * self.Sa_diff + self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0, self.Fa1 + self.Ea1 + self.Qa1, 1)) * self.Sa_diff + self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0, self.Fa1 + self.Ea1 + self.Qa1, 1)) * self.Sa_diff + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Ea - self.Fa - self.Qa + self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k]) + self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) + + self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Qaadd - self.Fa - self.Sa[k] + self.Sa_t[k] + + self.Ea_[k] = self.Ea + self.Qa_[k] = self.Qa + self.Qaadd + self.Fa_[k] = self.Fa + self.Ft_[k] = self.Ft + +def agriZone_Ep_Sa_beta_Fvar(self,k): + """ + - Potential evaporation is decreased by energy used for interception evaporation + - Formula for evaporation based on LP + - Outgoing fluxes are determined based on (value in previous timestep + inflow) + and if this leads to negative storage, the outgoing fluxes are corrected to rato --> Eu is + no longer taken into account for this correction + - Qa u is determined from overflow from Sa --> incorporation of beta function + - Fa is based on storage in Sa + - Code for ini-file: 8 + """ + + JarvisCoefficients.calcEp(self,k) + self.PotEvaporation = self.EpHour + + self.samax2 = self.samax[k] * scalar(self.TopoId) + self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0) + + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) + self.SaN = self.Sa[k] / self.samax2 + self.SuN = self.Su[k] / self.sumax[k] + + self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1) + self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k]) + self.Fa1 = self.cropG * ifthenelse(self.SaN > 0,self.Fmin[k] + (self.Fmax[k] - self.Fmin[k]) * e ** (-self.decF[k] * (1 - self.SaN)),0) + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Qa1 - self.Fa1 - self.Ea1 + + self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0) + self.Qa = self.Qa1 + (self.Qa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff + self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff + self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Ea - self.Fa - self.Qa + self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k]) + self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) + + self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Qaadd - self.Fa - self.Sa[k] + self.Sa_t[k] + + self.Ea_[k] = self.Ea + self.Qa_[k] = self.Qa + self.Qaadd + self.Fa_[k] = self.Fa + +def agriZone_hourlyEp_Sa_beta_Fvar(self,k): + """ + - Potential evaporation is decreased by energy used for interception evaporation + - Formula for evaporation based on LP + - Outgoing fluxes are determined based on (value in previous timestep + inflow) + and if this leads to negative storage, the outgoing fluxes are corrected to rato --> Eu is + no longer taken into account for this correction + - Qa u is determined from overflow from Sa --> incorporation of beta function + - Fa is based on storage in Sa + - Code for ini-file: 9 + """ + +# JarvisCoefficients.calcEp(self,k) +# self.PotEvaporation = self.EpHour + + self.samax2 = self.samax[k] * scalar(self.TopoId) + self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0) + + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) + self.SaN = self.Sa[k] / self.samax2 + self.SuN = self.Su[k] / self.sumax[k] + + self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1) + self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k]) + self.Fa1 = self.cropG * ifthenelse(self.SaN > 0,self.Fmin[k] + (self.Fmax[k] - self.Fmin[k]) * e ** (-self.decF[k] * (1 - self.SaN)),0) + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Qa1 - self.Fa1 - self.Ea1 + + self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0) + self.Qa = self.Qa1 + (self.Qa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff + self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff + self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Ea - self.Fa - self.Qa + self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k]) + self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) + + self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Qaadd - self.Fa - self.Sa[k] + self.Sa_t[k] + + self.Ea_[k] = self.Ea + self.Qa_[k] = self.Qa + self.Qaadd + self.Fa_[k] = self.Fa \ No newline at end of file Index: wflow-py/wflow/reservoir_Sf.py =================================================================== diff -u -rf95a643bab9f58880912cb18318fdb74dd6e10a9 -rfacc3a7f729f6910dac37d4761992c70d238fcbd --- wflow-py/wflow/reservoir_Sf.py (.../reservoir_Sf.py) (revision f95a643bab9f58880912cb18318fdb74dd6e10a9) +++ wflow-py/wflow/reservoir_Sf.py (.../reservoir_Sf.py) (revision facc3a7f729f6910dac37d4761992c70d238fcbd) @@ -32,14 +32,26 @@ """ This function is used when no unsaturated zone reservoir is used and only passes fluxes from the upper reservoirs to the lower - Qf = 0. + Qf = Qf_in. Storage in fast reservoir = 0. """ self.Qfin_[k] = self.Qu * (1 - self.D) self.Qf_[k] = self.Qfin_[k] self.Sf_[k] = 0. self.wbSf_[k] = self.Qfin_[k] - self.Qf_[k] - self.Sf[k] + self.Sf_t[k] + +def fastAgriRunoff_no_reservoir(self, k): + """ + This function is used when no unsaturated zone reservoir is used and only + passes fluxes from the upper reservoirs to the lower + Qfa = Qfa_in. + Storage in fast reservoir = 0. + """ + self.Qfa_[k] = self.Qa_[k] + self.Sfa[k] = 0. + self.wbSfa_[k] = self.Qfa_[k] - self.Qfa_[k] - self.Sfa[k] + self.Sfa_t[k] - sum(self.convQa[k]) + sum(self.convQa_t[k]) + def fastRunoff_lag2(self, k): """ - Lag is applied before inflow into the fast reservoir @@ -88,14 +100,106 @@ self.Qf_[k] = self.Qf # self.QuA_[k] = self.Qu +def fastRunoff_lag_forAgri_combined(self,k): + """ + - Lag is applied before inflow into the fast reservoir + - Lag formula is derived from Fenicia (2011) + - Outgoing fluxes are determined based on (value in previous timestep + inflow) + and if this leads to negative storage, the outgoing fluxes are corrected to rato + - not a semi analytical solution for Sf anymore + - When separate ditch/road fast reservoir is taken into account, this option should not be used + - Code for ini-file: 3 + """ + + if self.FR_L: + self.Qu = areatotal(self.Qu_[k] * self.percentArea ,nominal(self.TopoId)) + self.Qa = areatotal(self.Qa_[k] * self.percentArea ,nominal(self.TopoId)) + else: + self.Qu = self.Qu_[k] + self.Qa = self.Qa_[k] + self.Qfin = (1 - self.D[k]) * self.Qu + self.Qa + + + if self.D[k] < 1.00: + if self.convQu[k]: + self.QfinLag = self.convQu[k][-1] + self.Qf = self.Sf[k] * self.Kf[k] + self.Sf[k] = self.Sf[k] + self.QfinLag - self.Qf + + self.convQu[k].insert(0, 0 * scalar(self.TopoId)) #convolution Qu for following time steps + self.Tfmap = self.Tf[k] * scalar(self.TopoId) + del self.convQu[k][-1] + temp = [self.convQu[k][i] + (2/self.Tfmap-2/(self.Tfmap*(self.Tfmap+1))*(self.Tfmap-i))*self.Qfin for i in range(len(self.convQu[k]))] + self.convQu[k] = temp + + + else: + self.Qf = self.Sf[k] * self.Kf[k] + self.Sf[k] = self.Sf[k] + self.Qfin - self.Qf + + + else: + self.Qf = self.ZeroMap + self.Qfinput_[k] = self.ZeroMap + self.Qfin_[k] = self.ZeroMap + + self.wbSf_[k] = self.Qfin - self.Qf - self.Sf[k] + self.Sf_t[k] - sum(self.convQu[k]) + sum(self.convQu_t[k]) + + self.Qf_[k] = self.Qf + +def fastRunoff_lag_agriDitch(self,k): + """ + - Lag is applied before inflow into the fast reservoir + - Lag formula is derived from Fenicia (2011) + - Outgoing fluxes are determined based on (value in previous timestep + inflow) + and if this leads to negative storage, the outgoing fluxes are corrected to rato + - not a semi analytical solution for Sf anymore + - very fast responding reservoir to represent fast drainage via roads and ditches + - Code for ini-file: 4 + """ + + if self.FR_L: + self.Qa = areatotal(self.Qa_[k] * self.percentArea ,nominal(self.TopoId)) + else: + self.Qa = self.Qa_[k] + + self.Qfain = self.Qa + +# commented on 4 August 2015, as the output of this reservoir should not depent on the value of D +# if self.D[k] < 1.00: + if self.convQa[k]: + self.QfainLag = self.convQa[k][-1] + self.Qfa = self.Sfa[k] * self.Kfa[k] + self.Sfa[k] = self.Sfa[k] + self.QfainLag - self.Qfa + + self.convQa[k].insert(0, 0 * scalar(self.TopoId)) #convolution Qu for following time steps + self.Tfmap = self.Tfa[k] * scalar(self.TopoId) + del self.convQa[k][-1] + temp = [self.convQa[k][i] + (2/self.Tfmap-2/(self.Tfmap*(self.Tfmap+1))*(self.Tfmap-i))*self.Qfain for i in range(len(self.convQa[k]))] + self.convQa[k] = temp + + + else: + self.Qfa = self.Sfa[k] * self.Kfa[k] + self.Sfa[k] = self.Sfa[k] + self.Qfain - self.Qfa + +# commented on 4 August 2015, as the output of this reservoir should not depent on the value of D +# else: +# self.Qfa = self.ZeroMap +# self.Qfain_[k] = self.ZeroMap + + self.wbSfa_[k] = self.Qfain - self.Qfa - self.Sfa[k] + self.Sfa_t[k] - sum(self.convQa[k]) + sum(self.convQa_t[k]) + + self.Qfa_[k] = self.Qfa + def routingQf_combined(self): """ - Routing of fluxes from fast reservoir - Qf is devided over the reservoir numbers for the timesteps matching with the average travel time for a calcultation cell """ - + if nansum(pcr2numpy(self.Transit,NaN)) > 0: self.Qflag = self.trackQ[0] # first bucket is transferred to outlet self.trackQ.append(0 * scalar(self.TopoId)) # add new bucket for present time step @@ -115,6 +219,7 @@ self.Qtlag = self.Qflag_ / self.timestepsecs + self.Qs_ / 1000 * self.surfaceArea / self.timestepsecs self.QLagTot = areatotal(self.Qtlag, nominal(self.TopoId)) # catchment total runoff with looptijd + def routingQf_Qs_grid(self): """ @@ -132,6 +237,4 @@ # water balance of flux routing self.dSdt = self.Qstate-self.Qstate_t - self.WB_rout = (accuflux(self.TopoLdd, self.Qtotal - self.dSdt)-self.Qrout)/accuflux(self.TopoLdd, self.Qtotal) - - \ No newline at end of file + self.WB_rout = (accuflux(self.TopoLdd, self.Qtotal - self.dSdt)-self.Qrout)/accuflux(self.TopoLdd, self.Qtotal) \ No newline at end of file Index: wflow-py/wflow/reservoir_Si.py =================================================================== diff -u -rdc31062e54cbc5f26112e703b97ef3fd2d61ea11 -rfacc3a7f729f6910dac37d4761992c70d238fcbd --- wflow-py/wflow/reservoir_Si.py (.../reservoir_Si.py) (revision dc31062e54cbc5f26112e703b97ef3fd2d61ea11) +++ wflow-py/wflow/reservoir_Si.py (.../reservoir_Si.py) (revision facc3a7f729f6910dac37d4761992c70d238fcbd) @@ -50,15 +50,18 @@ self.Pe = max(self.Precipitation - (self.imax[k] - self.Si_t[k]),0) self.Si[k] = self.Si_t[k] + (self.Precipitation - self.Pe) - self.Ei = min(self.PotEvaporation, self.Si[k]) + self.Ei = ifthenelse(self.Sw[k] == 0, min(self.PotEvaporation, self.Si[k]), 0) #ifstatement added on 3-11-2015 for snow module self.Si[k] = self.Si[k] - self.Ei + self.wbSi_[k] = self.Precipitation - self.Ei - self.Pe - self.Si[k] + self.Si_t[k] + + self.Pe = self.Pe + self.Qw #added on 3-11-2015 for snow module + self.Ei = self.Ei + self.Ew + self.Ei_[k]=self.Ei self.Pe_[k]=self.Pe - self.PotEvaporation_ = self.PotEvaporation + self.PotEvaporation_ = self.PotEvaporation - self.wbSi_[k] = self.Precipitation - self.Ei - self.Pe - self.Si[k] + self.Si_t[k] - if self.URFR_L: self.Ei = areatotal(self.Ei * self.percentArea, nominal(self.TopoId)) self.Pe = areatotal(self.Pe * self.percentArea, nominal(self.TopoId)) @@ -75,20 +78,24 @@ """ JarvisCoefficients.calcEp(self,k) + self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0) self.Pe = max(self.Precipitation - (self.imax[k] - self.Si_t[k]),0) self.Si[k] = self.Si_t[k] + (self.Precipitation - self.Pe) - self.Ei = min(self.EpHour, self.Si[k]) + self.Ei = ifthenelse(self.Sw[k] == 0, min(self.PotEvaporation, self.Si[k]), 0) #ifstatement added on 3-11-2015 for snow module self.Si[k] = self.Si[k] - self.Ei + self.wbSi_[k] = self.Precipitation - self.Ei - self.Pe - self.Si[k] + self.Si_t[k] + + self.Pe = self.Pe + self.Qw #added on 3-11-2015 for snow module + self.Ei = self.Ei + self.Ew + self.Ei_[k]=self.Ei - self.Pe_[k]=self.Pe - self.Ep_[k]=self.EpHour - - self.wbSi_[k] = self.Precipitation - self.Ei - self.Pe - self.Si[k] + self.Si_t[k] + self.Pe_[k]=self.Pe + self.Ep_[k]=self.EpHour if self.URFR_L: self.Ei = areatotal(self.Ei * self.percentArea, nominal(self.TopoId)) self.Pe = areatotal(self.Pe * self.percentArea, nominal(self.TopoId)) self.PotEvaporation = areatotal(self.PotEvaporation * self.percentArea, nominal(self.TopoId)) - self.Si[k] = areatotal(self.Si[k] * self.percentArea, nominal(self.TopoId)) + self.Si[k] = areatotal(self.Si[k] * self.percentArea, nominal(self.TopoId)) \ No newline at end of file Index: wflow-py/wflow/reservoir_Ss.py =================================================================== diff -u -r3048a2b4b52a5328cd08bcfcb090d8b247dfdc20 -rfacc3a7f729f6910dac37d4761992c70d238fcbd --- wflow-py/wflow/reservoir_Ss.py (.../reservoir_Ss.py) (revision 3048a2b4b52a5328cd08bcfcb090d8b247dfdc20) +++ wflow-py/wflow/reservoir_Ss.py (.../reservoir_Ss.py) (revision facc3a7f729f6910dac37d4761992c70d238fcbd) @@ -50,8 +50,8 @@ self.QsinTemp = sum([a*b for a,b in zip(self.QsinClass ,self.percent)]) #fluxes are summed per cell according to percentage of class self.Qsin = areatotal(self.QsinTemp * self.percentArea ,nominal(self.TopoId)) #areatotal is taken, according to area percentage of cell - self.Qs = self.Ss * self.Ks - self.Ss = self.Ss * exp(-self.Ks) + self.Qsin + self.Qs = self.Ss * self.Ks[0] + self.Ss = self.Ss * exp(-self.Ks[0]) + self.Qsin self.wbSs = self.Qsin - self.Qs - self.Ss + self.Ss_t Index: wflow-py/wflow/reservoir_Su.py =================================================================== diff -u -r2021d28f559d05321903ec806b5dadb39f55b31a -rfacc3a7f729f6910dac37d4761992c70d238fcbd --- wflow-py/wflow/reservoir_Su.py (.../reservoir_Su.py) (revision 2021d28f559d05321903ec806b5dadb39f55b31a) +++ wflow-py/wflow/reservoir_Su.py (.../reservoir_Su.py) (revision facc3a7f729f6910dac37d4761992c70d238fcbd) @@ -78,7 +78,7 @@ self.Qu_[k] = max(self.Pe_[k], 0) self.Eu_[k] = 0. self.Perc_[k] = 0. - self.Su_[k] = 0. + self.Su[k] = 0. self.Cap_[k] = 0. self.wbSu_[k] = self.Pe - self.Eu - self.Qu - self.Perc + self.Cap - self.Su[k] + self.Su_t[k] @@ -225,13 +225,11 @@ JarvisCoefficients.calcEp(self,k) self.PotEvaporation = self.EpHour -# -# pdb.set_trace() + self.cropG_scal = pcr2numpy(self.cropG,NaN) if any(self.cropG_scal == 1): self.sumax2 = self.sumax[k] - elif any(self.cropG_scal > 0): -# pdb.set_trace() + elif any(self.cropG_scal > 0): self.sumax2 = self.sumax[k] * (1 - numpy.max(self.cropG_scal[self.cropG_scal >= 0]) * (1-self.redsu[k])) else: self.sumax2 = self.sumax[k] * self.redsu[k] @@ -591,4 +589,4 @@ self.Qu_[k] = self.Qa + self.Qaadd self.Fa_[k] = self.Fa self.Cap_[k] = self.Cap - self.Perc_[k] = self.Perc + self.Perc_[k] = self.Perc \ No newline at end of file Index: wflow-py/wflow/wflow_topoflex.py =================================================================== diff -u -r0c392f437b26afdbe80a0e4ffcc262cf5ad06540 -rfacc3a7f729f6910dac37d4761992c70d238fcbd --- wflow-py/wflow/wflow_topoflex.py (.../wflow_topoflex.py) (revision 0c392f437b26afdbe80a0e4ffcc262cf5ad06540) +++ wflow-py/wflow/wflow_topoflex.py (.../wflow_topoflex.py) (revision facc3a7f729f6910dac37d4761992c70d238fcbd) @@ -1,27 +1,29 @@ #!/usr/bin/python """ -Definition of the wflow_topoflex model. +Definition of the wflow_sceleton model. --------------------------------------- This simple model calculates soil temperature using air temperature as a forcing. Usage: -wflow_topoflex -C case -R Runid -c inifile +wflow_sceleton -C case -R Runid -c inifile -C: set the name of the case (directory) to run -R: set the name runId within the current case -c name of the config file (in the case directory) - +$Author: schelle $ +$Id: wflow_sceleton.py 898 2014-01-09 14:47:06Z schelle $ +$Rev: 898 $ """ -## TODO: add wflow prefix of put in seperate folder import reservoir_Si -# import reservoir_Sa +import reservoir_Sa +import reservoir_Sw import reservoir_Su import reservoir_Sf import reservoir_Ss @@ -33,11 +35,12 @@ import shutil, glob import getopt import time +import pdb from wflow.wf_DynamicFramework import * from wflow.wflow_adapt import * # import scipy -import copy +from copy import copy as copylist # TODO: see below """ @@ -104,15 +107,15 @@ # Meteo and other forcing # Meteo and other forcing - modelparameters.append( - self.ParamType(name="Temperature", stack='intss/T.tss', type="timeseries", default=10.0, verbose=False, - lookupmaps=['staticmaps/wflow_subcatch.map'])) - modelparameters.append( - self.ParamType(name="Precipitation", stack='intss/P.tss', type="timeseries", default=0.0, verbose=False, - lookupmaps=['staticmaps/wflow_subcatch.map'])) - modelparameters.append( - self.ParamType(name="PotEvaporation", stack='intss/PET.tss', type="timeseries", default=0.0, verbose=False, - lookupmaps=['staticmaps/wflow_subcatch.map'])) +# modelparameters.append( +# self.ParamType(name="Temperature", stack='intss/T.tss', type="timeseries", default=10.0, verbose=False, +# lookupmaps=['staticmaps/wflow_subcatch.map'])) +# modelparameters.append( +# self.ParamType(name="Precipitation", stack='intss/P.tss', type="timeseries", default=0.0, verbose=False, +# lookupmaps=['staticmaps/wflow_subcatch.map'])) +# modelparameters.append( +# self.ParamType(name="PotEvaporation", stack='intss/PET.tss', type="timeseries", default=0.0, verbose=False, +# lookupmaps=['staticmaps/wflow_subcatch.map'])) return modelparameters @@ -272,11 +275,21 @@ self.Classes = [x for x in range(len(self.NamesClasses))] # numbering of classes # selection of reservoir conceputalisatie - codes are described in reservoir files + self.selectSw = configget(self.config, "model", + "selectSw", "0, 0, 0").replace( + ' ', '').replace('[', '').replace( + ']', '').replace( + 'None', '').split(',') self.selectSi = configget(self.config, "model", "selectSi", "0, 0, 0").replace( ' ', '').replace('[', '').replace( ']', '').replace( 'None', '').split(',') + self.selectSa = configget(self.config, "model", + "selectSa", "0, 0, 0").replace( + ' ', '').replace('[', '').replace( + ']', '').replace( + 'None', '').split(',') self.selectSu = configget(self.config, "model", "selectSu", "0, 0, 0").replace( ' ', '').replace('[', '').replace( @@ -287,6 +300,11 @@ ' ', '').replace('[', '').replace( ']', '').replace( 'None', '').split(',') + self.selectSfa = configget(self.config, "model", + "selectSfa", "0, 0, 0").replace( + ' ', '').replace('[', '').replace( + ']', '').replace( + 'None', '').split(',') self.selectSs = configget(self.config, "model", "selectSs", "groundWaterCombined3") self.selectSr = configget(self.config, "model", "selectSr", "0, 0, 0").replace( @@ -347,19 +365,31 @@ self.samax = eval(str(configget(self.config, "model", "samax", "[0]"))) self.srmax = eval(str(configget(self.config, "model", "sumax", "[0]"))) self.beta = eval(str(configget(self.config, "model", "beta", "[0]"))) - self.famax = eval(str(configget(self.config, "model", "famax", "[0]"))) + self.Fmax = eval(str(configget(self.config,"model","Fmax","[0]"))) + self.Fmin = eval(str(configget(self.config,"model","Fmin","[0]"))) + self.decF = eval(str(configget(self.config,"model","decF","[0]"))) self.Ce = eval(str(configget(self.config, "model", "Ce", "[0]"))) self.Co = eval(str(configget(self.config, "model", "Co", "[0]"))) self.D = eval(str(configget(self.config, "model", "D", "[0]"))) self.Kf = eval(str(configget(self.config, "model", "Kf", "[0]"))) + self.Kfa = eval(str(configget(self.config, "model", "Kfa", "[0]"))) self.Tf = eval(str(configget(self.config, "model", "Tf", "[0]"))) + self.Tfa = eval(str(configget(self.config, "model", "Tfa", "[0]"))) self.imax = eval(str(configget(self.config, "model", "imax", "[0]"))) self.perc = eval(str(configget(self.config, "model", "perc", "[0]"))) self.cap = eval(str(configget(self.config, "model", "cap", "[0]"))) self.Kd = eval(str(configget(self.config, "model", "Kd", "[0]"))) self.Kr = eval(str(configget(self.config, "model", "Kr", "[0]"))) self.LP = eval(str(configget(self.config, "model", "LP", "[0]"))) self.Ks = eval(str(configget(self.config, "model", "Ks", "[0]"))) + self.dayDeg = eval(str(configget(self.config,"model","dayDeg","[0]"))) + self.FrDur0 = eval(str(configget(self.config,"model","FrDur0","[0]"))) + self.FrDur1 = eval(str(configget(self.config,"model","FrDur1","[0]"))) + self.ratFT = eval(str(configget(self.config,"model","ratFT","[0]"))) + self.Tt = eval(str(configget(self.config,"model","Tt","[0]"))) + self.Tm = eval(str(configget(self.config,"model","Tm","[0]"))) + self.Fm = eval(str(configget(self.config,"model","Fm","[0]"))) + # Jarvis stressfunctions self.JC_Topt = eval(str(configget(self.config, "model", "JC_Topt", "[0]"))) self.JC_D05 = eval(str(configget(self.config, "model", "JC_D05", "[0]"))) @@ -381,9 +411,8 @@ # initialise list for lag function self.convQu = [[0 * scalar(self.TopoId)] * self.Tf[i] for i in self.Classes] + self.convQa = [[0 * scalar(self.TopoId)] * self.Tfa[i] for i in self.Classes] - - self.wf_updateparameters() if self.scalarInput: self.gaugesMap = nominal(readmap(os.path.join(self.Dir, wflow_mgauges))) #: Map with locations of rainfall/evap/temp gauge(s). Only needed if the input to the model is not in maps @@ -435,9 +464,11 @@ # self.logger.info("Setting initial conditions to default (zero!)") self.logger.info("Setting initial conditions to preset values in main script!!") self.Si = [self.ZeroMap] * len(self.Classes) + self.Sw = [self.ZeroMap] * len(self.Classes) self.Su = [self.ZeroMap] * len(self.Classes) self.Sa = [self.ZeroMap] * len(self.Classes) self.Sf = [self.ZeroMap] * len(self.Classes) + self.Sfa = [self.ZeroMap] * len(self.Classes) self.Sr = [self.ZeroMap] * len(self.Classes) # self.Ss = [self.ZeroMap] * len(self.Classes) # for separate gw reservoir per class self.Ss = self.ZeroMap # for combined gw reservoir @@ -455,30 +486,35 @@ self.wbSi_ = [self.ZeroMap] * len(self.Classes) self.wbSu_ = [self.ZeroMap] * len(self.Classes) self.wbSa_ = [self.ZeroMap] * len(self.Classes) + self.wbSw_ = [self.ZeroMap] * len(self.Classes) self.wbSr_ = [self.ZeroMap] * len(self.Classes) self.wbSf_ = [self.ZeroMap] * len(self.Classes) + self.wbSfa_ = [self.ZeroMap] * len(self.Classes) self.wbSfrout = self.ZeroMap self.wbSs = self.ZeroMap self.Ei_ = [self.ZeroMap] * len(self.Classes) self.Pe_ = [self.ZeroMap] * len(self.Classes) self.Si_ = [self.ZeroMap] * len(self.Classes) self.Eu_ = [self.ZeroMap] * len(self.Classes) - self.Er_ = [self.ZeroMap] * len(self.Classes) + self.Ea_ = [self.ZeroMap] * len(self.Classes) + self.Ew_ = [self.ZeroMap] * len(self.Classes) self.Qu_ = [self.ZeroMap] * len(self.Classes) - self.Qd_ = [self.ZeroMap] * len(self.Classes) - self.Qo_ = [self.ZeroMap] * len(self.Classes) - self.Qr_ = [self.ZeroMap] * len(self.Classes) + self.Qw_ = [self.ZeroMap] * len(self.Classes) + self.Qa_ = [self.ZeroMap] * len(self.Classes) self.Cap_ = [self.ZeroMap] * len(self.Classes) self.Perc_ = [self.ZeroMap] * len(self.Classes) self.Fa_ = [self.ZeroMap] * len(self.Classes) self.Qf_ = [self.ZeroMap] * len(self.Classes) + self.Qfa_ = [self.ZeroMap] * len(self.Classes) # self.Qs_ = [self.ZeroMap] * len(self.Classes) # for separate gw reservoir per class self.Qs_ = self.ZeroMap # for combined gw reservoir self.Qflag_ = [self.ZeroMap] * len(self.Classes) self.Qfcub_ = [self.ZeroMap] * len(self.Classes) self.Ep_ = [self.ZeroMap] * len(self.Classes) self.EpD_ = [self.ZeroMap] * len(self.Classes) + self.FrDur = [self.ZeroMap] * len(self.Classes) + self.Ft_ = [self.ZeroMap] * len(self.Classes) self.JC_temp_ = [self.ZeroMap] * len(self.Classes) self.JC_vpd_ = [self.ZeroMap] * len(self.Classes) @@ -509,15 +545,18 @@ self.logger.debug("Step: "+str(int(self.thestep + self._d_firstTimeStep))+"/"+str(int(self._d_nrTimeSteps))) self.thestep = self.thestep + 1 - self.Si_t = copy.deepcopy(self.Si) - self.Su_t = copy.deepcopy(self.Su) - self.Sa_t = copy.deepcopy(self.Sa) - self.Sf_t = copy.deepcopy(self.Sf) - self.Sr_t = copy.deepcopy(self.Sr) + self.Si_t = copylist(self.Si) + self.Sw_t = copylist(self.Sw) + self.Su_t = copylist(self.Su) + self.Sa_t = copylist(self.Sa) + self.Sf_t = copylist(self.Sf) + self.Sfa_t = copylist(self.Sfa) + self.Sr_t = copylist(self.Sr) self.Ss_t = self.Ss - self.trackQ_t = copy.deepcopy(self.trackQ) # copylist(self.trackQ) - self.convQu_t = [copy.deepcopy(self.convQu[i]) for i in self.Classes] # copylist(self.convQu) - + self.trackQ_t = copylist(self.trackQ) # copylist(self.trackQ) + self.convQu_t = [copylist(self.convQu[i]) for i in self.Classes] # copylist(self.convQu) + self.convQa_t = [copylist(self.convQa[i]) for i in self.Classes] + # if self.scalarInput: # if self.InputSeries == 1: # self.Precipitation = timeinputscalar(self.precipTss, self.gaugesMap) @@ -551,35 +590,59 @@ self.PotEvaporation = areatotal(self.PotEvaporation * self.percentArea, nominal(self.TopoId)) self.Precipitation = areatotal(self.Precipitation * self.percentArea, nominal(self.TopoId)) self.Temperature = areaaverage(self.Temperature * self.percentArea, nominal(self.TopoId)) + + self.PrecipTotal = self.Precipitation + if self.selectSw[0] > 0: + self.Precipitation = ifthenelse(self.Temperature >= self.Tt[0], self.PrecipTotal,0) + self.PrecipitationSnow = ifthenelse(self.Temperature < self.Tt[0], self.PrecipTotal,0) + for k in self.Classes: # SNOW ================================================================================================= + if self.selectSw[k]: + eval_str = 'reservoir_Sw.{:s}(self, k)'.format(self.selectSw[k]) + else: + eval_str = 'reservoir_Sw.snow_no_reservoir(self, k)' + eval(eval_str) - # nu nog even niet gecodeerd - # INTERCEPTION ========================================================================================= if self.selectSi[k]: eval_str = 'reservoir_Si.{:s}(self, k)'.format(self.selectSi[k]) else: eval_str = 'reservoir_Si.interception_no_reservoir(self, k)' eval(eval_str) + + # AGRICULTURE ZONE ====================================================================================== + if self.selectSa[k]: + eval_str = 'reservoir_Sa.{:s}(self, k)'.format(self.selectSa[k]) + else: + eval_str = 'reservoir_Sa.agriZone_no_reservoir(self, k)' + eval(eval_str) + # UNSATURATED ZONE ====================================================================================== if self.selectSu[k]: eval_str = 'reservoir_Su.{:s}(self, k)'.format(self.selectSu[k]) else: - eval_str = 'reservoir_Si.unsatZone_no_reservoir(self, k)' + eval_str = 'reservoir_Su.unsatZone_no_reservoir(self, k)' eval(eval_str) # FAST RUNOFF RESERVOIR =================================================================================== if self.selectSf[k]: eval_str = 'reservoir_Sf.{:s}(self, k)'.format(self.selectSf[k]) else: - eval_str = 'reservoir_Si.fastRunoff_no_reservoir(self, k)' + eval_str = 'reservoir_Sf.fastRunoff_no_reservoir(self, k)' eval(eval_str) + #FAST AGRICULTURE DITCHES RUNOFF RESERVOIR =================================================================================== + if self.selectSfa[k]: + eval_str = 'reservoir_Sf.{:s}(self, k)'.format(self.selectSfa[k]) + else: + eval_str = 'reservoir_Sf.fastAgriRunoff_no_reservoir(self, k)' + eval(eval_str) + # RIPARIAN ZONE RESERVOIR ================================================================================== if self.selectSr[k]: eval_str = 'reservoir_Sr.{:s}(self, k)'.format(self.selectSr[k]) @@ -590,7 +653,8 @@ # TOTAL RUNOFF ============================================================================================= # self.Qfcub = (sum([x*y for x,y in zip(self.Qf_,self.percent)]) + sum([x*y for x,y in zip(self.Qo_,self.percent)]) + sum([x*y for x,y in zip(self.Qd_,self.percent)]) + sum([x*y for x,y in zip(self.Qr_,self.percent)]))/ 1000 * self.surfaceArea - self.Qftotal = (sum([x * y for x, y in zip(self.Qf_, self.percent)])) +# pdb.set_trace() + self.Qftotal = sum([x * y for x, y in zip(self.Qf_, self.percent)]) + sum([x*y for x,y in zip(self.Qfa_,self.percent)]) if self.selectSs: eval_str = 'reservoir_Ss.{:s}(self)'.format(self.selectSs) else: @@ -611,10 +675,12 @@ self.QtlagWB = (self.Qtlag / self.surfaceArea) * 1000 * self.timestepsecs self.convQuWB = [sum(self.convQu[i]) for i in self.Classes] self.convQuWB_t = [sum(self.convQu_t[i]) for i in self.Classes] + self.convQaWB = [sum(self.convQa[i]) for i in self.Classes] + self.convQaWB_t = [sum(self.convQa_t[i]) for i in self.Classes] self.trackQWB = (sum(self.trackQ) / self.surfaceArea) * 1000 self.trackQWB_t = (sum(self.trackQ_t) / self.surfaceArea) * 1000 self.WB = self.Precipitation - sum(multiply(self.Ei_, self.percent)) - sum( - multiply(self.Eu_, self.percent)) - sum(multiply(self.Er_, self.percent)) - self.QtlagWB - sum( + multiply(self.Eu_, self.percent)) - self.QtlagWB - sum( multiply(self.Si, self.percent)) + sum(multiply(self.Si_t, self.percent)) - sum( multiply(self.Su, self.percent)) + sum(multiply(self.Su_t, self.percent)) - sum( multiply(self.Sf, self.percent)) + sum(multiply(self.Sf_t, self.percent)) - sum( @@ -624,54 +690,53 @@ multiply(self.convQuWB, self.percent)) + sum(multiply(self.convQuWB_t, self.percent)) # #fuxes and states in m3/h - self.P = areatotal(self.Precipitation / 1000 * self.surfaceArea, nominal(self.TopoId)) - self.Ei = areatotal(sum(multiply(self.Ei_, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId)) - self.Eu = areatotal(sum(multiply(self.Eu_, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId)) - self.Er = areatotal(sum(multiply(self.Er_, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId)) + self.P = areatotal(self.PrecipTotal / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.Ei = areatotal(sum(multiply(self.Ei_,self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.Ea = areatotal(sum(multiply(self.Ea_,self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.Eu = areatotal(sum(multiply(self.Eu_,self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) +# self.Er = areatotal(sum(multiply(self.Er_,self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) +# self.QtotnoRout = areatotal(self.Qftotal + self.Qs_/ 1000 * self.surfaceArea,nominal(self.TopoId)) self.Qtot = self.QLagTot * self.timestepsecs - #self.QtotnoRout = areatotal(self.Qtotal, nominal(self.TopoId)) - self.SiWB = areatotal(sum(multiply(self.Si, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId)) - self.Si_WB = areatotal(sum(multiply(self.Si_t, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId)) - self.SuWB = areatotal(sum(multiply(self.Su, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId)) - self.Su_WB = areatotal(sum(multiply(self.Su_t, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId)) - self.SaWB = areatotal(sum(multiply(self.Sa, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId)) - self.Sa_WB = areatotal(sum(multiply(self.Sa_t, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId)) - self.SfWB = areatotal(sum(multiply(self.Sf, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId)) - self.Sf_WB = areatotal(sum(multiply(self.Sf_t, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId)) - self.SrWB = areatotal(sum(multiply(self.Sr, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId)) - self.Sr_WB = areatotal(sum(multiply(self.Sr_t, self.percent)) / 1000 * self.surfaceArea, nominal(self.TopoId)) - self.SsWB = areatotal(self.Ss / 1000 * self.surfaceArea, nominal(self.TopoId)) - self.Ss_WB = areatotal(self.Ss_t / 1000 * self.surfaceArea, nominal(self.TopoId)) - self.convQuWB = areatotal( - sum(multiply([sum(self.convQu[i]) for i in self.Classes], self.percent)) / 1000 * self.surfaceArea, - nominal(self.TopoId)) - self.convQu_WB = areatotal( - sum(multiply([sum(self.convQu_t[i]) for i in self.Classes], self.percent)) / 1000 * self.surfaceArea, - nominal(self.TopoId)) - self.trackQWB = areatotal(sum(self.trackQ) , nominal(self.TopoId)) - self.trackQ_WB = areatotal(sum(self.trackQ_t), nominal(self.TopoId)) - self.QstateWB = areatotal(sum(self.Qstate) * 3600, nominal(self.TopoId)) - self.Qstate_WB = areatotal(sum(self.Qstate_t) * 3600, nominal(self.TopoId)) + self.SiWB = areatotal(sum(multiply(self.Si,self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.Si_WB = areatotal(sum(multiply(self.Si_t,self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.SuWB = areatotal(sum(multiply(self.Su,self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.Su_WB = areatotal(sum(multiply(self.Su_t,self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.SaWB = areatotal(sum(multiply(self.Sa,self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.Sa_WB = areatotal(sum(multiply(self.Sa_t,self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.SfWB = areatotal(sum(multiply(self.Sf,self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.Sf_WB = areatotal(sum(multiply(self.Sf_t,self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.SfaWB = areatotal(sum(multiply(self.Sfa,self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.Sfa_WB = areatotal(sum(multiply(self.Sfa_t,self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.SrWB = areatotal(sum(multiply(self.Sr,self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.Sr_WB = areatotal(sum(multiply(self.Sr_t,self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.SwWB = areatotal(sum(multiply(self.Sw,self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.Sw_WB = areatotal(sum(multiply(self.Sw_t,self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.SsWB = areatotal(self.Ss / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.Ss_WB = areatotal(self.Ss_t / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.convQuWB = areatotal(sum(multiply([sum(self.convQu[i]) for i in self.Classes],self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.convQu_WB = areatotal(sum(multiply([sum(self.convQu_t[i]) for i in self.Classes],self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.convQaWB = areatotal(sum(multiply([sum(self.convQa[i]) for i in self.Classes],self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.convQa_WB = areatotal(sum(multiply([sum(self.convQa_t[i]) for i in self.Classes],self.percent)) / 1000 * self.surfaceArea,nominal(self.TopoId)) + self.trackQWB = areatotal(sum(self.trackQ),nominal(self.TopoId)) + self.trackQ_WB = areatotal(sum(self.trackQ_t),nominal(self.TopoId)) - # WBtot in m3/s - self.WBtot = ( - self.P - self.Ei - self.Eu - self.Er - self.Qtot - self.SiWB + self.Si_WB - self.SuWB + self.Su_WB - self.SaWB + self.Sa_WB - self.SfWB + self.Sf_WB - self.SrWB + self.Sr_WB - self.SsWB + self.Ss_WB - self.convQuWB + self.convQu_WB - self.trackQWB + self.trackQ_WB - self.QstateWB + self.Qstate_WB) / self.timestepsecs + #WBtot in m3/s + self.WBtot = (self.P - self.Ei - self.Ea - self.Eu - self.Qtot - self.SiWB + self.Si_WB - self.SuWB + self.Su_WB - self.SaWB + self.Sa_WB - self.SwWB + self.Sw_WB - self.SfWB + self.Sf_WB - self.SfaWB + self.Sfa_WB - self.SrWB + self.Sr_WB - self.SsWB + self.Ss_WB - self.convQuWB +self.convQu_WB - self.convQaWB +self.convQa_WB - self.trackQWB + self.trackQ_WB) / self.timestepsecs - # SUMMED FLUXES ====================================================================================== self.sumprecip = self.sumprecip + self.Precipitation # accumulated rainfall for water balance (m/h) self.sumevap = self.sumevap + sum(multiply(self.Ei_, self.percent)) + sum( multiply(self.Eu_, self.percent)) + sum( - multiply(self.Er_, self.percent)) # accumulated evaporation for water balance (m/h) + multiply(self.Ea_, self.percent)) + sum( + multiply(self.Ew_, self.percent)) # accumulated evaporation for water balance (m/h) try: self.sumpotevap = self.sumpotevap + self.PotEvaporation # accumulated potential evaporation (m/h) except: self.sumpotevap = self.EpHour self.sumrunoff = self.sumrunoff + self.Qtlag * 1000 * self.timestepsecs / self.surfaceArea # accumulated runoff for water balance (m/h) self.sumwb = self.sumwb + self.WB - self.sumE = sum(multiply(self.Ei_, self.percent)) + sum(multiply(self.Eu_, self.percent)) + sum( - multiply(self.Er_, self.percent)) + self.sumE = sum(multiply(self.Ei_, self.percent)) + sum(multiply(self.Eu_, self.percent)) # The main function is used to run the program from the command line