Index: wflow-py/wflow/reservoir_Sa.py =================================================================== diff -u -r978510e581d205d904b07bc305a5d0550779c581 -rb0a0649fbc4d5ed289c09f5beffbd40a711e95ca --- wflow-py/wflow/reservoir_Sa.py (.../reservoir_Sa.py) (revision 978510e581d205d904b07bc305a5d0550779c581) +++ wflow-py/wflow/reservoir_Sa.py (.../reservoir_Sa.py) (revision b0a0649fbc4d5ed289c09f5beffbd40a711e95ca) @@ -25,6 +25,38 @@ import scipy import JarvisCoefficients +def selectSaR(i): + """ + not all functions are still in this file, the older functions can be found + (with the same numbering) in h:\My Documents\memo's\python scripts\wflow\ + """ + if i == 1: + name = 'agriZone_Jarvis' + elif i == 2: + name = 'agriZone_Ep' + elif i == 3: + 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' + elif i == 12: + name = 'agriZone_Ep_Sa_beta_frostSamax' + elif i == 13: + name = 'agriZone_Ep_Sa_beta_frostSamax_surfTemp' + return name def agriZone_no_reservoir(self, k): """ @@ -42,7 +74,155 @@ 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): + """ + - Potential evaporation is decreased by energy used for interception evaporation + - Formula for evaporation based on Jarvis stress functions + - 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 + - Code for ini-file: 1 + """ + self.Qa = max(self.Pe - (self.samax[k] - self.Sa_t[k]),0) + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa) + self.SaN = min(self.Sa[k] / self.samax2, 1) + self.SuN = self.Su[k] / self.sumax[k] + JarvisCoefficients.calcEu(self,k,1) #calculation of Ea based on Jarvis stress functions + self.Ea1 = self.Eu + + 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 + + self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0) + self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 > 0 , self.Fa1 + self.Ea1 , 1)) * self.Sa_diff + self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 > 0 , self.Fa1 + self.Ea1 , 1)) * self.Sa_diff + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa) - self.Ea - self.Fa + 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.Fa - self.Sa[k] + self.Sa_t[k] + + self.Ea_[k] = self.Ea + self.Qa_[k] = self.Qa + self.Fa_[k] = self.Fa + + +def agriZone_Ep(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 + - Code for ini-file: 2 + """ + JarvisCoefficients.calcEp(self,k) + self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0) + + self.Qa = max(self.Pe - (self.samax[k] - self.Sa_t[k]),0) + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa) + self.SaN = min(self.Sa[k] / self.samax2, 1) + self.SuN = self.Su[k] / self.sumax[k] + + self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax[k] * self.LP[k]),1) + + 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 + + self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0) + self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 > 0 , self.Fa1 + self.Ea1 , 1)) * self.Sa_diff + self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 > 0 , self.Fa1 + self.Ea1 , 1)) * self.Sa_diff + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa) - self.Ea - self.Fa + 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.Fa - self.Sa[k] + self.Sa_t[k] + + self.Ea_[k] = self.Ea + self.Qa_[k] = self.Qa + self.Fa_[k] = self.Fa + +def agriZone_Ep_Sa(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 + - Fa is based on storage in Sa + - Code for ini-file: 3 + """ + JarvisCoefficients.calcEp(self,k) + self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0) + + self.Qa = max(self.Pe - (self.samax[k] - self.Sa_t[k]),0) + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa) + self.SaN = min(self.Sa[k] / self.samax2, 1) + self.SuN = self.Su[k] / self.sumax[k] + + self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax[k] * self.LP[k]),1) + + 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.Qa) - self.Fa1 - self.Ea1 + + self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0) + self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 > 0 , self.Fa1 + self.Ea1 , 1)) * self.Sa_diff + self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 > 0 , self.Fa1 + self.Ea1 , 1)) * self.Sa_diff + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa) - self.Ea - self.Fa + 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.Fa - self.Sa[k] + self.Sa_t[k] + + self.Ea_[k] = self.Ea + self.Qa_[k] = self.Qa + self.Fa_[k] = self.Fa + + +def agriZone_Ep_Sa_cropG(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 + - Fa is based on storage in Sa + - Code for ini-file: 4 + """ + JarvisCoefficients.calcEp(self,k) + self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0) + + self.samax2 = self.samax[k] * self.cropG + self.Qaadd = max(self.Sa_t[k] - self.samax2,0) + + 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 = min(self.Sa[k] / self.samax2, 1) + 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.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.Qa) - self.Fa1 - self.Ea1 + + self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0) + self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 > 0 , self.Fa1 + self.Ea1 , 1)) * self.Sa_diff + self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 > 0 , self.Fa1 + self.Ea1 , 1)) * self.Sa_diff + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa) - self.Ea - self.Fa + 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.Fa - self.Sa[k] + self.Sa_t[k] + + 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 @@ -52,6 +232,7 @@ 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) @@ -61,7 +242,7 @@ 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.SaN = min(self.Sa[k] / self.samax2, 1) 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) @@ -93,16 +274,17 @@ 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.samax2 = self.samax[k] * scalar(self.catchArea) 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 = max(self.Sa[k] / self.samax2, 1) + self.SaN = min(max(self.Sa[k] / self.samax2, 0), 1) 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) @@ -124,7 +306,52 @@ 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.catchArea) + 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 = min(self.Sa[k] / self.samax2, 1) + 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 @@ -135,14 +362,18 @@ - 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.samax2 = self.samax[k] * scalar(self.catchArea) 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.SaN = min(self.Sa[k] / self.samax2, 1) 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) @@ -176,12 +407,18 @@ 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] + ifthenelse(self.Temperature > 0, self.ratFT[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]), self.samin[k]), 1) + 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.samax2 = self.samax[k] * scalar(self.catchArea) * 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) @@ -219,6 +456,7 @@ - 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) @@ -227,12 +465,12 @@ self.FrDur[k] = min(self.FrDur[k] + ifthenelse(self.Temperature > 0, self.ratFT[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]), self.samin[k]), 1) - self.samax2 = self.samax[k] * scalar(self.TopoId) * self.Ft + self.samax2 = self.samax[k] * scalar(self.catchArea) * 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 = min(self.Sa[k] / self.samax2, 1) - self.SuN = self.Su[k] / self.sumax + 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]) @@ -266,19 +504,20 @@ - 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: 13 """ JarvisCoefficients.calcEp(self,k) - self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0) + self.PotEvaporation = self.EpHour self.FrDur[k] = min(self.FrDur[k] + ifthenelse(self.TempSurf > 0, self.ratFT[k] * self.TempSurf, self.TempSurf) * 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]), self.samin[k]), 1) - self.samax2 = self.samax[k] * scalar(self.TopoId) * self.Ft + self.samax2 = self.samax[k] * scalar(self.catchArea) * 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.SaN = min(self.Sa[k] / self.samax2, 1) 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) @@ -303,3 +542,84 @@ 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.catchArea) + 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 = min(self.Sa[k] / self.samax2, 1) + 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.catchArea) + 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 = min(self.Sa[k] / self.samax2, 1) + 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 -r978510e581d205d904b07bc305a5d0550779c581 -rb0a0649fbc4d5ed289c09f5beffbd40a711e95ca --- wflow-py/wflow/reservoir_Sf.py (.../reservoir_Sf.py) (revision 978510e581d205d904b07bc305a5d0550779c581) +++ wflow-py/wflow/reservoir_Sf.py (.../reservoir_Sf.py) (revision b0a0649fbc4d5ed289c09f5beffbd40a711e95ca) @@ -17,6 +17,17 @@ from wf_DynamicFramework import * import scipy +def selectSfR(i): + """ + not all functions are still in this file, the older functions can be found + (with the same numbering) in h:\My Documents\memo's\python scripts\wflow\ + """ + if i == 1: + name = 'fastRunoff_lag' + if i == 2: + name = 'fastRunoff_lag2' + return name + def fastRunoff_no_reservoir(self, k): """ This function is used when no unsaturated zone reservoir is used and only @@ -48,6 +59,7 @@ - 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 + - Code for ini-file: 2 """ if self.FR_L: @@ -63,16 +75,21 @@ 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) + self.convQu[k].insert(0, 0 * scalar(self.catchArea)) #convolution Qu for following time steps + self.Tfmap = self.Tf[k] * scalar(self.catchArea) 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 +# self.Qfin_[k] = self.Qfin +# self.Qfinput_[k] = self.QfinLag + else: self.Qf = self.Sf[k] * self.Kf[k] self.Sf[k] = self.Sf[k] + self.Qfin - self.Qf +# self.Qfin_[k] = self.Qfin +# self.Qfinput_[k] = self.Qfin else: self.Qf = self.ZeroMap self.Qfinput_[k] = self.ZeroMap @@ -81,6 +98,7 @@ 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 +# self.QuA_[k] = self.Qu def fastRunoff_lag_forAgri_combined(self,k): @@ -91,6 +109,7 @@ 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: @@ -109,8 +128,8 @@ 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) + self.convQu[k].insert(0, 0 * scalar(self.catchArea)) #convolution Qu for following time steps + self.Tfmap = self.Tf[k] * scalar(self.catchArea) 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 @@ -138,6 +157,7 @@ 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: @@ -147,25 +167,74 @@ 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) + self.convQa[k].insert(0, 0 * scalar(self.catchArea)) #convolution Qu for following time steps + self.Tfmap = self.Tfa[k] * scalar(self.catchArea) 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 fastRunoff_lag_agriDitch_reInfilt(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: 5 + """ + + 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 + + if self.convQa[k]: + self.QfainLag = self.convQa[k][-1] + self.Qfa = self.Sfa[k] * self.Kfa[k] + self.Percfa = ifthenelse(self.Ft_[k] == 1, self.Sfa[k] * self.perc[k] * 200, 0) + self.Sfa[k] = self.Sfa[k] + self.QfainLag - self.Qfa - self.Percfa + + self.convQa[k].insert(0, 0 * scalar(self.catchArea)) #convolution Qu for following time steps + self.Tfmap = self.Tfa[k] * scalar(self.catchArea) + 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.Percfa = ifthenelse(self.Ft_[k] == 1, self.Sfa[k] * self.perc[k] * 200, 0) + self.Sfa[k] = self.Sfa[k] + self.Qfain - self.Qfa - self.Percfa + + self.wbSfa_[k] = self.Qfain - self.Qfa - self.Sfa[k] + self.Sfa_t[k] - sum(self.convQa[k]) + sum(self.convQa_t[k]) - self.Percfa + + self.Qfa_[k] = self.Qfa + self.Percfa_[k] = self.Percfa + def routingQf_combined(self): """ - Routing of fluxes from fast reservoir @@ -175,12 +244,12 @@ 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 + self.trackQ.append(0 * scalar(self.catchArea)) # add new bucket for present time step del self.trackQ[0] # remove first bucket (transferred to outlet) - temp = [self.trackQ[i] + ifthenelse(rounddown(self.Transit) == i*scalar(self.TopoId), - (self.Transit - i*scalar(self.TopoId)) * self.Qftotal / 1000 * self.surfaceArea, - ifthenelse(roundup(self.Transit) == i*scalar(self.TopoId), - (i*scalar(self.TopoId) - self.Transit) * self.Qftotal / 1000 * self.surfaceArea, 0)) for i in range(len(self.trackQ))] + temp = [self.trackQ[i] + ifthenelse(rounddown(self.Transit) == i*scalar(self.catchArea), + (self.Transit - i*scalar(self.catchArea)) * self.Qftotal / 1000 * self.surfaceArea, + ifthenelse(roundup(self.Transit) == i*scalar(self.catchArea), + (i*scalar(self.catchArea) - self.Transit) * self.Qftotal / 1000 * self.surfaceArea, 0)) for i in range(len(self.trackQ))] self.trackQ = temp else: @@ -201,7 +270,9 @@ """ self.Qtot = self.Qftotal + self.Qs_ # total local discharge in mm/hour self.Qtotal = self.Qtot / 1000 * self.surfaceArea / self.timestepsecs # total local discharge in m3/s + self.QtotNoRout = accuflux(self.TopoLdd, self.Qtotal) self.Qstate_t = self.Qstate + self.Qtest = self.Qstate + self.Qtotal self.Qrout = accutraveltimeflux(self.TopoLdd, self.Qstate + self.Qtotal, self.velocity) self.Qstate = accutraveltimestate(self.TopoLdd, self.Qstate + self.Qtotal, self.velocity) @@ -210,4 +281,24 @@ # 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) + +def routingQf_Qs_grid_mm(self): + """ + - Routing of both Qf and Qs + - based on a velocity map + """ + self.Qtot = self.Qftotal + self.Qs_ # total local discharge in mm/hour + self.Qtotal = self.Qtot + self.QtotNoRout = accuflux(self.TopoLdd, self.Qtotal) + self.Qstate_t = self.Qstate + self.Qtest = self.Qstate + self.Qtotal + self.Qrout = accutraveltimeflux(self.TopoLdd, self.Qstate + self.Qtotal, self.velocity) + self.Qstate = accutraveltimestate(self.TopoLdd, self.Qstate + self.Qtotal, self.velocity) + + self.Qtlag = self.Qrout + self.QLagTot = self.Qrout + + # 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 Index: wflow-py/wflow/reservoir_Si.py =================================================================== diff -u -r978510e581d205d904b07bc305a5d0550779c581 -rb0a0649fbc4d5ed289c09f5beffbd40a711e95ca --- wflow-py/wflow/reservoir_Si.py (.../reservoir_Si.py) (revision 978510e581d205d904b07bc305a5d0550779c581) +++ wflow-py/wflow/reservoir_Si.py (.../reservoir_Si.py) (revision b0a0649fbc4d5ed289c09f5beffbd40a711e95ca) @@ -17,6 +17,19 @@ import JarvisCoefficients +def selectSiR(i): + """ + not all functions are still in this file, the older functions can be found + (with the same numbering) in h:\My Documents\memo's\python scripts\wflow\ + """ + if i == 1: + name = 'interception_overflow' + if i == 2: + name = 'interception_overflow2' + if i == 3: + name = 'interception_overflow_Ep' + return name + def interception_no_reservoir(self, k): """ Effective rainfall = rainfall @@ -32,6 +45,7 @@ """ - Effective rainfall is all that does not fit into the interception reservoir - Outgoing fluxes are determined separately + - Code for ini-file: 2 """ self.Pe = max(self.Precipitation - (self.imax[k] - self.Si_t[k]),0) @@ -60,6 +74,7 @@ - Effective rainfall is all that does not fit into the interception reservoir - Outgoing fluxes are determined separately - this version cannot be run with Su averaged (for the current code) + - Code for ini-file: 3 """ JarvisCoefficients.calcEp(self,k) Index: wflow-py/wflow/reservoir_Su.py =================================================================== diff -u -r978510e581d205d904b07bc305a5d0550779c581 -rb0a0649fbc4d5ed289c09f5beffbd40a711e95ca --- wflow-py/wflow/reservoir_Su.py (.../reservoir_Su.py) (revision 978510e581d205d904b07bc305a5d0550779c581) +++ wflow-py/wflow/reservoir_Su.py (.../reservoir_Su.py) (revision b0a0649fbc4d5ed289c09f5beffbd40a711e95ca) @@ -18,6 +18,53 @@ import scipy import JarvisCoefficients +def selectSuR(i): + """ + not all functions are still in this file, the older functions can be found + (with the same numbering) in h:\My Documents\memo's\python scripts\wflow\ + """ + if i == 1: + name = 'unsatZone_LP_beta' + elif i == 2: + name = 'unsatZone_LP' + elif i == 3: + name = 'unsatZone_SF' + elif i == 4: + name = 'unsatZone_EisEp' + elif i == 5: + name = 'unsatZone_drain' + elif i == 6: + name = 'unsatZone_drain_split' + elif i == 7: + name = 'unsatZone_drain_Sumin' + elif i == 8: + name = 'unsatZone_drain_DZ' + elif i == 9: + name = 'unsatZone_drain_split_DZ' + elif i == 10: + name = 'unsatZone_withAgri' + elif i == 11: + name = 'unsatZone_withAgri_oldBeta' + elif i == 12: + name = 'unsatZone_LP_beta_Jarvis' + elif i == 13: + name = 'unsatZone_LP_beta_Ep' + elif i == 14: + name = 'unsatZone_withAgri_Ep' + elif i == 15: + name = 'unsatZone_withAgri_Jarvis' + elif i == 16: + name = 'unsatZone_forAgri_Jarvis' + elif i == 17: + name = 'unsatZone_forAgri_Ep' + elif i == 18: + name = 'unsatZone_forAgri_Jarvis_cropG' + elif i == 19: + name = 'unsatZone_forAgri_Ep_cropG' + elif i == 20: + name = 'unsatZone_LP_beta_Ep_cropG' + return name + def unsatZone_no_reservoir(self, k): """ This function is used when no unsaturated zone reservoir is used and only @@ -41,8 +88,9 @@ - Formula for evaporation linear until LP, from than with potential rate - 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 + - Qu is determined with a beta function (same as in HBV?) + - Code for ini-file: 1 """ - self.Su[k] = ifthenelse(self.Su_t[k] + self.Pe > self.sumax[k], self.sumax[k], self.Su_t[k] + self.Pe) self.Quadd = ifthenelse(self.Su_t[k] + self.Pe > self.sumax[k], self.Su_t[k] + self.Pe - self.sumax[k], 0) self.SuN = self.Su[k] / self.sumax[k] @@ -71,7 +119,9 @@ self.Qu_[k] = self.Qu + self.Quadd self.Cap_[k] = self.Cap self.Perc_[k] = self.Perc - +# self.Su_diff_[k] = self.Su_diff +# self.Quadd_[k] = self.Quadd + def unsatZone_LP_beta_Jarvis(self,k): """ @@ -81,6 +131,7 @@ and if this leads to negative storage, the outgoing fluxes are corrected to rato --> Eu is no longer taken into account for this correction - Qu is determined with a beta function (same as in HBV?) + - Code for ini-file: 12 """ self.Su[k] = ifthenelse(self.Su_t[k] + self.Pe > self.sumax[k], self.sumax[k], self.Su_t[k] + self.Pe) self.Quadd = ifthenelse(self.Su_t[k] + self.Pe > self.sumax[k], self.Su_t[k] + self.Pe - self.sumax[k], 0) @@ -109,6 +160,8 @@ self.Qu_[k] = self.Qu + self.Quadd self.Cap_[k] = self.Cap self.Perc_[k] = self.Perc +# self.Su_diff_[k] = self.Su_diff +# self.Quadd_[k] = self.Quadd def unsatZone_LP_beta_Ep(self,k): @@ -120,17 +173,68 @@ - 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 - Qu is determined with a beta function (same as in HBV?) + - Code for ini-file: 13 """ + + #pdb.set_trace() + JarvisCoefficients.calcEp(self,k) + self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0) + + self.Su[k] = ifthenelse(self.Su_t[k] + self.Pe > self.sumax[k], self.sumax[k], self.Su_t[k] + self.Pe) + self.Quadd = ifthenelse(self.Su_t[k] + self.Pe > self.sumax[k], self.Su_t[k] + self.Pe - self.sumax[k], 0) + self.SuN = self.Su[k] / self.sumax[k] + self.SiN = self.Si[k] / self.imax[k] + + self.Eu1 = max((self.PotEvaporation - self.Ei),0) * min(self.Su[k] / (self.sumax[k] * self.LP[k]),1) + + self.Qu1 = (self.Pe - self.Quadd) * (1 - (1 - self.SuN) ** self.beta[k]) + self.Perc1 = self.perc[k] * self.SuN + self.Su[k] = self.Su_t[k] + (self.Pe - self.Quadd) - self.Qu1 - self.Eu1 - self.Perc1 + self.Su_diff = ifthenelse(self.Su[k] < 0, self.Su[k], 0) + self.Eu = self.Eu1 + (self.Eu1 / ifthenelse(self.Qu1 + self.Eu1 +self.Perc1 > 0 , self.Qu1 + self.Eu1 + self.Perc1 , 1)) * self.Su_diff + self.Qu = self.Qu1 + (self.Qu1/ifthenelse(self.Qu1 + self.Eu1 + self.Perc1 > 0 , self.Qu1 + self.Eu1 + self.Perc1 , 1)) * self.Su_diff + self.Perc = ifthenelse (self.Perc1 > 0, self.Perc1 + (self.Perc1/ifthenelse(self.Qu1 + self.Eu1 + self.Perc1 > 0 , self.Qu1 + self.Eu1 + self.Perc1 , 1)) * self.Su_diff, self.Perc1) + self.Su[k] = self.Su_t[k] + (self.Pe - self.Quadd) - self.Eu - self.Qu - self.Perc + self.Su[k] = ifthenelse(self.Su[k] < 0, 0 , self.Su[k]) + self.Su_diff2 = ifthen(self.Su[k] < 0, self.Su[k]) + + self.Cap = min(self.cap[k] * (1 - self.Su[k] / self.sumax[k]), self.Ss) + self.Su[k] = self.Su[k] + self.Cap + + self.wbSu_[k] = self.Pe - self.Eu - self.Qu - self.Quadd - self.Perc + self. Cap - self.Su[k] + self.Su_t[k] + + self.Eu_[k] = self.Eu + self.Qu_[k] = self.Qu + self.Quadd + self.Cap_[k] = self.Cap + self.Perc_[k] = self.Perc +# self.Su_diff_[k] = self.Su_diff +# self.Quadd_[k] = self.Quadd + self.Epot_[k] = self.PotEvaporation + +def unsatZone_LP_beta_Ep_Ei(self,k): + """ + - Potential evaporation is calculated with formula in 'JarvisCoefficients', but without + using the Jarvis stress functions + - Potential evaporation is decreased by energy used for interception evaporation + - Formula for evaporation linear until LP, from than with potential rate + - 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 + - Qu is determined with a beta function (same as in HBV?) + - Code for ini-file: 13 + """ + + #pdb.set_trace() JarvisCoefficients.calcEp(self,k) self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0) - self.Su[k] = ifthenelse(self.Su_t[k] + self.Pe > self.sumax, self.sumax, self.Su_t[k] + self.Pe) - self.Quadd = ifthenelse(self.Su_t[k] + self.Pe > self.sumax, self.Su_t[k] + self.Pe - self.sumax, 0) - self.SuN = self.Su[k] / self.sumax + self.Su[k] = ifthenelse(self.Su_t[k] + self.Pe > self.sumax[k], self.sumax[k], self.Su_t[k] + self.Pe) + self.Quadd = ifthenelse(self.Su_t[k] + self.Pe > self.sumax[k], self.Su_t[k] + self.Pe - self.sumax[k], 0) + self.SuN = self.Su[k] / self.sumax[k] self.SiN = self.Si[k] / self.imax[k] - self.Eu1 = max((self.PotEvaporation - self.Ei),0) * min(self.Su[k] / (self.sumax * self.LP[k]),1) + self.Eu1 = ifthenelse(self.SiN == 1, 0, max((self.PotEvaporation),0) * min(self.Su[k] / (self.sumax[k] * self.LP[k]),1)) +# self.Eu1 = max((self.PotEvaporation),0) * min(self.Su[k] / (self.sumax[k] * self.LP[k]),1) self.Qu1 = (self.Pe - self.Quadd) * (1 - (1 - self.SuN) ** self.beta[k]) self.Perc1 = self.perc[k] * self.SuN @@ -144,7 +248,7 @@ self.Su[k] = ifthenelse(self.Su[k] < 0, 0 , self.Su[k]) self.Su_diff2 = ifthen(self.Su[k] < 0, self.Su[k]) - self.Cap = min(self.cap[k] * (1 - self.Su[k] / self.sumax), self.Ss) + self.Cap = min(self.cap[k] * (1 - self.Su[k] / self.sumax[k]), self.Ss) self.Su[k] = self.Su[k] + self.Cap self.wbSu_[k] = self.Pe - self.Eu - self.Qu - self.Quadd - self.Perc + self. Cap - self.Su[k] + self.Su_t[k] @@ -153,7 +257,154 @@ self.Qu_[k] = self.Qu + self.Quadd self.Cap_[k] = self.Cap self.Perc_[k] = self.Perc +# self.Su_diff_[k] = self.Su_diff +# self.Quadd_[k] = self.Quadd + self.Epot_[k] = self.PotEvaporation +def unsatZone_LP_beta_Ep_percD(self,k): + """ + - Potential evaporation is calculated with formula in 'JarvisCoefficients', but without + using the Jarvis stress functions + - Potential evaporation is decreased by energy used for interception evaporation + - Formula for evaporation linear until LP, from than with potential rate + - 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 + - Qu is determined with a beta function (same as in HBV?) + - Code for ini-file: 13 + """ + + JarvisCoefficients.calcEp(self,k) + self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0) + + self.Su[k] = ifthenelse(self.Su_t[k] + self.Pe > self.sumax[k], self.sumax[k], self.Su_t[k] + self.Pe) + self.Quadd = ifthenelse(self.Su_t[k] + self.Pe > self.sumax[k], self.Su_t[k] + self.Pe - self.sumax[k], 0) + self.SuN = self.Su[k] / self.sumax[k] + self.SiN = self.Si[k] / self.imax[k] + + self.Eu1 = max((self.PotEvaporation - self.Ei),0) * min(self.Su[k] / (self.sumax[k] * self.LP[k]),1) + + self.percDeep = self.percD[-1] + self.Qu1 = (self.Pe - self.Quadd) * (1 - (1 - self.SuN) ** self.beta[k]) + self.Perc1 = self.perc[k] * self.SuN + self.percDeep + self.Su[k] = self.Su_t[k] + (self.Pe - self.Quadd) - self.Qu1 - self.Eu1 - self.Perc1 + + self.Su_diff = ifthenelse(self.Su[k] < 0, self.Su[k], 0) + self.Eu = self.Eu1 + (self.Eu1 / ifthenelse(self.Qu1 + self.Eu1 +self.Perc1 > 0 , self.Qu1 + self.Eu1 + self.Perc1 , 1)) * self.Su_diff + self.Qu = self.Qu1 + (self.Qu1/ifthenelse(self.Qu1 + self.Eu1 + self.Perc1 > 0 , self.Qu1 + self.Eu1 + self.Perc1 , 1)) * self.Su_diff + self.Perc = ifthenelse (self.Perc1 > 0, self.Perc1 + (self.Perc1/ifthenelse(self.Qu1 + self.Eu1 + self.Perc1 > 0 , self.Qu1 + self.Eu1 + self.Perc1 , 1)) * self.Su_diff, self.Perc1) + self.Su[k] = self.Su_t[k] + (self.Pe - self.Quadd) - self.Eu - self.Qu - self.Perc + self.Su[k] = ifthenelse(self.Su[k] < 0, 0 , self.Su[k]) + self.Su_diff2 = ifthen(self.Su[k] < 0, self.Su[k]) + + self.Cap = min(self.cap[k] * (1 - self.Su[k] / self.sumax[k]), self.Ss) + self.Su[k] = self.Su[k] + self.Cap + + self.wbSu_[k] = self.Pe - self.Eu - self.Qu - self.Quadd - self.Perc + self. Cap - self.Su[k] + self.Su_t[k] + + self.Eu_[k] = self.Eu + self.Qu_[k] = self.Qu + self.Quadd + self.Cap_[k] = self.Cap + self.Perc_[k] = self.Perc + self.Epot_[k] = self.PotEvaporation + +def unsatZone_LP_beta_Ep_percD(self,k): + """ + - Potential evaporation is calculated with formula in 'JarvisCoefficients', but without + using the Jarvis stress functions + - Potential evaporation is decreased by energy used for interception evaporation + - Formula for evaporation linear until LP, from than with potential rate + - 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 + - Qu is determined with a beta function (same as in HBV?) + - Code for ini-file: 13 + """ + + JarvisCoefficients.calcEp(self,k) + self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0) + + self.Su[k] = ifthenelse(self.Su_t[k] + self.Pe > self.sumax[k], self.sumax[k], self.Su_t[k] + self.Pe) + self.Quadd = ifthenelse(self.Su_t[k] + self.Pe > self.sumax[k], self.Su_t[k] + self.Pe - self.sumax[k], 0) + self.SuN = self.Su[k] / self.sumax[k] + self.SiN = self.Si[k] / self.imax[k] + + self.Eu1 = max((self.PotEvaporation - self.Ei),0) * min(self.Su[k] / (self.sumax[k] * self.LP[k]),1) + + self.percDeep = self.percD[-1] + self.Qu1 = (self.Pe - self.Quadd) * (1 - (1 - self.SuN) ** self.beta[k]) + self.Perc1 = self.perc[k] * self.SuN + self.percDeep + self.Su[k] = self.Su_t[k] + (self.Pe - self.Quadd) - self.Qu1 - self.Eu1 - self.Perc1 + + self.Su_diff = ifthenelse(self.Su[k] < 0, self.Su[k], 0) + self.Eu = self.Eu1 + (self.Eu1 / ifthenelse(self.Qu1 + self.Eu1 +self.Perc1 > 0 , self.Qu1 + self.Eu1 + self.Perc1 , 1)) * self.Su_diff + self.Qu = self.Qu1 + (self.Qu1/ifthenelse(self.Qu1 + self.Eu1 + self.Perc1 > 0 , self.Qu1 + self.Eu1 + self.Perc1 , 1)) * self.Su_diff + self.Perc = ifthenelse (self.Perc1 > 0, self.Perc1 + (self.Perc1/ifthenelse(self.Qu1 + self.Eu1 + self.Perc1 > 0 , self.Qu1 + self.Eu1 + self.Perc1 , 1)) * self.Su_diff, self.Perc1) + self.Su[k] = self.Su_t[k] + (self.Pe - self.Quadd) - self.Eu - self.Qu - self.Perc + self.Su[k] = ifthenelse(self.Su[k] < 0, 0 , self.Su[k]) + self.Su_diff2 = ifthen(self.Su[k] < 0, self.Su[k]) + + self.Cap = min(self.cap[k] * (1 - self.Su[k] / self.sumax[k]), self.Ss) + self.Su[k] = self.Su[k] + self.Cap + + self.wbSu_[k] = self.Pe - self.Eu - self.Qu - self.Quadd - self.Perc + self. Cap - self.Su[k] + self.Su_t[k] + + self.Eu_[k] = self.Eu + self.Qu_[k] = self.Qu + self.Quadd + self.Cap_[k] = self.Cap + self.Perc_[k] = self.Perc + self.Epot_[k] = self.PotEvaporation + + +def unsatZone_LP_beta_Ep_percDvar(self,k): + """ + - Potential evaporation is calculated with formula in 'JarvisCoefficients', but without + using the Jarvis stress functions + - Potential evaporation is decreased by energy used for interception evaporation + - Formula for evaporation linear until LP, from than with potential rate + - 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 + - Qu is determined with a beta function (same as in HBV?) + - Code for ini-file: 13 + """ + + JarvisCoefficients.calcEp(self,k) + self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0) + + self.Su[k] = ifthenelse(self.Su_t[k] + self.Pe > self.sumax[k], self.sumax[k], self.Su_t[k] + self.Pe) + self.Quadd = ifthenelse(self.Su_t[k] + self.Pe > self.sumax[k], self.Su_t[k] + self.Pe - self.sumax[k], 0) + self.SuN = self.Su[k] / self.sumax[k] + self.SiN = self.Si[k] / self.imax[k] + + self.drought = ifthenelse(self.SuN < self.LP[k], self.TopoId, ifthenelse(pcrand(self.SuN < 0.8, self.drought), self.TopoId, boolean(scalar(self.TopoId) * 0))) + self.stijg = max(min(scalar(ifthenelse(self.drought == 1, self.stijg + self.Su_t[k] - self.Su_t2[k], 0)), self.sumax[k] * 100), 0) + + self.Eu1 = max((self.PotEvaporation - self.Ei),0) * min(self.Su[k] / (self.sumax[k] * self.LP[k]),1) + + self.Qu1 = (self.Pe - self.Quadd) * (1 - (1 - self.SuN) ** self.beta[k]) +# self.percDeep = max(10 * (1 - self.Ss / 30) * self.perc[k], 0) + self.percDeep = 0.8 * self.stijg * self.perc[k] + self.Perc1 = self.perc[k] * self.SuN + self.percDeep + self.Su[k] = self.Su_t[k] + (self.Pe - self.Quadd) - self.Qu1 - self.Eu1 - self.Perc1 + + self.Su_diff = ifthenelse(self.Su[k] < 0, self.Su[k], 0) + self.Eu = self.Eu1 + (self.Eu1 / ifthenelse(self.Qu1 + self.Eu1 +self.Perc1 > 0 , self.Qu1 + self.Eu1 + self.Perc1 , 1)) * self.Su_diff + self.Qu = self.Qu1 + (self.Qu1/ifthenelse(self.Qu1 + self.Eu1 + self.Perc1 > 0 , self.Qu1 + self.Eu1 + self.Perc1 , 1)) * self.Su_diff + self.Perc = ifthenelse (self.Perc1 > 0, self.Perc1 + (self.Perc1/ifthenelse(self.Qu1 + self.Eu1 + self.Perc1 > 0 , self.Qu1 + self.Eu1 + self.Perc1 , 1)) * self.Su_diff, self.Perc1) + self.Su[k] = self.Su_t[k] + (self.Pe - self.Quadd) - self.Eu - self.Qu - self.Perc + self.Su[k] = ifthenelse(self.Su[k] < 0, 0 , self.Su[k]) + self.Su_diff2 = ifthen(self.Su[k] < 0, self.Su[k]) + + self.Cap = min(self.cap[k] * (1 - self.Su[k] / self.sumax[k]), self.Ss) + self.Su[k] = self.Su[k] + self.Cap + + self.wbSu_[k] = self.Pe - self.Eu - self.Qu - self.Quadd - self.Perc + self. Cap - self.Su[k] + self.Su_t[k] + + self.Eu_[k] = self.Eu + self.Qu_[k] = self.Qu + self.Quadd + self.Cap_[k] = self.Cap + self.Perc_[k] = self.Perc + self.Epot_[k] = self.PotEvaporation + self.percDeep_[k] = self.percDeep + def unsatZone_LP_beta_Ep_cropG(self,k): """ - Potential evaporation is calculated with formula in 'JarvisCoefficients', but without @@ -164,6 +415,7 @@ and if this leads to negative storage, the outgoing fluxes are corrected to rato - Qu is determined with a beta function (same as in HBV?) - root zone storage for crop land is decreased in autumn and winter + - Code for ini-file: 20 """ JarvisCoefficients.calcEp(self,k) @@ -205,6 +457,8 @@ self.Qu_[k] = self.Qu + self.Quadd self.Cap_[k] = self.Cap self.Perc_[k] = self.Perc +# self.Su_diff_[k] = self.Su_diff +# self.Quadd_[k] = self.Quadd def unsatZone_forAgri_Jarvis(self,k): """ @@ -215,6 +469,7 @@ no longer taken into account for this correction - Qu is determined with a beta function (same as in HBV?) - inflow is infiltration from agriculture reservoir + - Code for ini-file: 16 """ self.Su[k] = ifthenelse(self.Su_t[k] + self.Fa > self.sumax[k], self.sumax[k], self.Su_t[k] + self.Fa) self.Quadd = ifthenelse(self.Su_t[k] + self.Fa > self.sumax[k], self.Su_t[k] + self.Fa - self.sumax[k], 0) @@ -255,6 +510,7 @@ no longer taken into account for this correction - Qu is determined with a beta function (same as in HBV?) - inflow is infiltration from agriculture reservoir + - Code for ini-file: 17 """ JarvisCoefficients.calcEp(self,k) @@ -265,7 +521,7 @@ self.SuN = self.Su[k] / self.sumax[k] self.SiN = self.Si[k] / self.imax[k] - self.Eu1 = ifthenelse(self.Ft_[k] == 1, max((self.PotEvaporation - self.Ei - self.Ea),0) * min(self.Su[k] / (self.sumax * self.LP[k]),1), 0) # no transpiration in case of frozen soil. Added on 22 feb 2016 + self.Eu1 = ifthenelse(self.Ft_[k] == 1, max((self.PotEvaporation - self.Ei - self.Ea),0) * min(self.Su[k] / (self.sumax[k] * self.LP[k]),1), 0) # no transpiration in case of frozen soil. Added on 22 feb 2016 self.Qu1 = (self.Fa - self.Quadd) * (1 - (1 - self.SuN) ** self.beta[k]) self.Perc1 = self.perc[k] * self.SuN @@ -289,6 +545,100 @@ self.Cap_[k] = self.Cap self.Perc_[k] = self.Perc +def unsatZone_forAgri_Ep_percD(self,k): + """ + - Potential evaporation is decreased by energy used for interception evaporation + - Formula for evaporation based on beta/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 + - Qu is determined with a beta function (same as in HBV?) + - inflow is infiltration from agriculture reservoir + - Code for ini-file: 17 + """ + + JarvisCoefficients.calcEp(self,k) + self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0) + + self.Su[k] = ifthenelse(self.Su_t[k] + self.Fa > self.sumax[k], self.sumax[k], self.Su_t[k] + self.Fa) + self.Quadd = ifthenelse(self.Su_t[k] + self.Fa > self.sumax[k], self.Su_t[k] + self.Fa - self.sumax[k], 0) + self.SuN = self.Su[k] / self.sumax[k] + self.SiN = self.Si[k] / self.imax[k] + + self.Eu1 = ifthenelse(self.Ft_[k] == 1, max((self.PotEvaporation - self.Ei - self.Ea),0) * min(self.Su[k] / (self.sumax[k] * self.LP[k]),1), 0) # no transpiration in case of frozen soil. Added on 22 feb 2016 + + self.percDeep = self.percD[-1] + self.Qu1 = (self.Fa - self.Quadd) * (1 - (1 - self.SuN) ** self.beta[k]) + self.Perc1 = self.perc[k] * self.SuN + self.percDeep + self.Su[k] = self.Su_t[k] + (self.Fa - self.Quadd) - self.Qu1 - self.Eu - self.Perc1 + + self.Su_diff = ifthenelse(self.Su[k] < 0, self.Su[k], 0) + self.Eu = self.Eu1 + (self.Eu1 / ifthenelse(self.Qu1 + self.Eu1 + self.Perc1 > 0 , self.Qu1 + self.Eu1 + self.Perc1 , 1)) * self.Su_diff + self.Qu = self.Qu1 + (self.Qu1 / ifthenelse(self.Qu1 + self.Eu1 + self.Perc1 > 0 , self.Qu1 + self.Eu1 + self.Perc1 , 1)) * self.Su_diff + self.Perc = ifthenelse (self.Perc1 > 0, self.Perc1 + (self.Perc1 / ifthenelse(self.Qu1 + self.Eu1 + self.Perc1 > 0 , self.Qu1 + self.Eu1 + self.Perc1 , 1)) * self.Su_diff, self.Perc1) + self.Su[k] = self.Su_t[k] + (self.Fa - self.Quadd) - self.Eu - self.Qu - self.Perc + self.Su[k] = ifthenelse(self.Su[k] < 0, 0 , self.Su[k]) + self.Su_diff2 = ifthen(self.Su[k] < 0, self.Su[k]) + + self.Cap = min(self.cap[k] * (1 - self.Su[k] / self.sumax[k]), self.Ss) + self.Su[k] = self.Su[k] + self.Cap + + self.wbSu_[k] = self.Fa - self.Eu - self.Qu - self.Quadd - self.Perc + self. Cap - self.Su[k] + self.Su_t[k] + + self.Eu_[k] = self.Eu + self.Qu_[k] = self.Qu + self.Quadd + self.Cap_[k] = self.Cap + self.Perc_[k] = self.Perc + +def unsatZone_forAgri_Ep_percDvar(self,k): + """ + - Potential evaporation is decreased by energy used for interception evaporation + - Formula for evaporation based on beta/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 + - Qu is determined with a beta function (same as in HBV?) + - inflow is infiltration from agriculture reservoir + - Code for ini-file: 17 + """ + + JarvisCoefficients.calcEp(self,k) + self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0) + + self.Su[k] = ifthenelse(self.Su_t[k] + self.Fa > self.sumax[k], self.sumax[k], self.Su_t[k] + self.Fa) + self.Quadd = ifthenelse(self.Su_t[k] + self.Fa > self.sumax[k], self.Su_t[k] + self.Fa - self.sumax[k], 0) + self.SuN = self.Su[k] / self.sumax[k] + self.SiN = self.Si[k] / self.imax[k] + + self.drought = ifthenelse(self.SuN < self.LP[k], self.TopoId, ifthenelse(pcrand(self.SuN < 0.8, self.drought), self.TopoId, boolean(scalar(self.TopoId) * 0))) + self.stijg = max(min(scalar(ifthenelse(self.drought == 1, self.stijg + self.Su_t[k] - self.Su_t2[k], 0)), self.sumax[k] * 100), 0) +# self.stijg = max(self.Su_t[k] - self.Su_t2[k], 0) + + self.Eu1 = ifthenelse(self.Ft_[k] == 1, max((self.PotEvaporation - self.Ei - self.Ea),0) * min(self.Su[k] / (self.sumax[k] * self.LP[k]),1), 0) # no transpiration in case of frozen soil. Added on 22 feb 2016 + + self.percDeep = 0.8 * self.stijg * self.perc[k] + self.Qu1 = (self.Fa - self.Quadd) * (1 - (1 - self.SuN) ** self.beta[k]) + self.Perc1 = self.perc[k] * self.SuN + self.percDeep + self.Su[k] = self.Su_t[k] + (self.Fa - self.Quadd) - self.Qu1 - self.Eu - self.Perc1 + + self.Su_diff = ifthenelse(self.Su[k] < 0, self.Su[k], 0) + self.Eu = self.Eu1 + (self.Eu1 / ifthenelse(self.Qu1 + self.Eu1 + self.Perc1 > 0 , self.Qu1 + self.Eu1 + self.Perc1 , 1)) * self.Su_diff + self.Qu = self.Qu1 + (self.Qu1 / ifthenelse(self.Qu1 + self.Eu1 + self.Perc1 > 0 , self.Qu1 + self.Eu1 + self.Perc1 , 1)) * self.Su_diff + self.Perc = ifthenelse (self.Perc1 > 0, self.Perc1 + (self.Perc1 / ifthenelse(self.Qu1 + self.Eu1 + self.Perc1 > 0 , self.Qu1 + self.Eu1 + self.Perc1 , 1)) * self.Su_diff, self.Perc1) + self.Su[k] = self.Su_t[k] + (self.Fa - self.Quadd) - self.Eu - self.Qu - self.Perc + self.Su[k] = ifthenelse(self.Su[k] < 0, 0 , self.Su[k]) + self.Su_diff2 = ifthen(self.Su[k] < 0, self.Su[k]) + + self.Cap = min(self.cap[k] * (1 - self.Su[k] / self.sumax[k]), self.Ss) + self.Su[k] = self.Su[k] + self.Cap + + self.wbSu_[k] = self.Fa - self.Eu - self.Qu - self.Quadd - self.Perc + self. Cap - self.Su[k] + self.Su_t[k] + + self.Eu_[k] = self.Eu + self.Qu_[k] = self.Qu + self.Quadd + self.Cap_[k] = self.Cap + self.Perc_[k] = self.Perc + def unsatZone_forAgri_hourlyEp(self,k): """ - Potential evaporation is decreased by energy used for interception evaporation @@ -298,6 +648,7 @@ no longer taken into account for this correction - Qu is determined with a beta function (same as in HBV?) - inflow is infiltration from agriculture reservoir + - Code for ini-file: 25 """ self.Su[k] = ifthenelse(self.Su_t[k] + self.Fa > self.sumax[k], self.sumax[k], self.Su_t[k] + self.Fa) @@ -339,6 +690,7 @@ no longer taken into account for this correction - Qu is determined with a beta function (same as in HBV?) - inflow is infiltration from agriculture reservoir + - Code for ini-file: 18 """ self.cropG_scal = pcr2numpy(self.cropG,NaN) if any(self.cropG_scal == 1): @@ -385,6 +737,7 @@ no longer taken into account for this correction - Qu is determined with a beta function (same as in HBV?) - inflow is infiltration from agriculture reservoir + - Code for ini-file: 19 """ JarvisCoefficients.calcEp(self,k) @@ -425,4 +778,146 @@ self.Cap_[k] = self.Cap self.Perc_[k] = self.Perc +def unsatZone_withAgri(self,k): + """ + - Potential evaporation is decreased by energy used for interception evaporation + - Formula for evaporation linear until LP, from than with potential rate + - 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 + - Qu is determined with a beta function (same as in HBV?) + - Code for ini-file: 10 + """ + self.Sa[k] = ifthenelse(self.Sa_t[k] + self.Pe > self.samax[k], self.samax[k], self.Sa_t[k] + self.Pe) + self.Qaadd = ifthenelse(self.Sa_t[k] + self.Pe > self.samax[k], self.Sa_t[k] + self.Pe - self.samax[k], 0) + self.SaN = self.Sa[k] / self.samax[k] + + self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax[k] * self.LP[k]),1) + self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k]) + self.Fa1 = self.famax[k] * (self.sumax[k] - self.Su[k]) / self.sumax[k] + + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Qa1 - self.Ea1 - self.Fa1 + + self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0) + self.Ea = self.Ea1 + (self.Ea1 / ifthenelse(self.Qa1 + self.Ea1 +self.Fa1 > 0 , self.Qa1 + self.Ea1 + self.Fa1 , 1)) * self.Sa_diff + self.Qa = self.Qa1 + (self.Qa1/ifthenelse(self.Qa1 + self.Ea1 + self.Fa1 > 0 , self.Qa1 + self.Ea1 + self.Fa1 , 1)) * self.Sa_diff + self.Fa = ifthenelse (self.Fa1 > 0, self.Fa1 + (self.Fa1/ifthenelse(self.Qa1 + self.Ea1 + self.Fa1 > 0 , self.Qa1 + self.Ea1 + self.Fa1 , 1)) * self.Sa_diff, self.Fa1) + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Ea - self.Qa - self.Fa + self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k]) + self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) + + self.Capa = min(self.cap[k] * (1 - self.Sa[k] / self.samax[k]), self.Su[k]) + self.Sa[k] = self.Sa[k] + self.Capa + + self.Su[k] = self.Su_t[k] + self.Fa - self.Capa + self.Perc = self.perc[k] * (self.Su[k] / self.sumax[k]) + self.Su[k] = self.Su[k] - self.Perc + + self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Qaadd - self.Fa + self. Capa - self.Sa[k] + self.Sa_t[k] + self.wbSu_[k] = self.Fa - self.Perc - self. Capa - self.Su[k] + self.Su_t[k] + + self.Eu_[k] = self.Ea + self.Qu_[k] = self.Qa + self.Qaadd + self.Fa_[k] = self.Fa + self.Cap_[k] = self.Cap + self.Perc_[k] = self.Perc +# self.Su_diff_[k] = self.Su_diff +# self.Quadd_[k] = self.Qaadd + +def unsatZone_withAgri_Ep(self,k): + """ + - Potential evaporation is calculated with formula in 'JarvisCoefficients', but without + using the Jarvis stress functions + - Potential evaporation is decreased by energy used for interception evaporation + - Formula for evaporation linear until LP, from than with potential rate + - 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 + - Qu is determined with a beta function (same as in HBV?) + - Code for ini-file: 14 + """ + + JarvisCoefficients.calcEp(self,k) + self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0) + + self.Sa[k] = ifthenelse(self.Sa_t[k] + self.Pe > self.samax[k], self.samax[k], self.Sa_t[k] + self.Pe) + self.Qaadd = ifthenelse(self.Sa_t[k] + self.Pe > self.samax[k], self.Sa_t[k] + self.Pe - self.samax[k], 0) + self.SaN = self.Sa[k] / self.samax[k] + + self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax[k] * self.LP[k]),1) + self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k]) + self.Fa1 = self.famax[k] * (self.sumax[k] - self.Su[k]) / self.sumax[k] + + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Qa1 - self.Ea1 - self.Fa1 + + self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0) + self.Ea = self.Ea1 + (self.Ea1 / ifthenelse(self.Qa1 + self.Ea1 +self.Fa1 > 0 , self.Qa1 + self.Ea1 + self.Fa1 , 1)) * self.Sa_diff + self.Qa = self.Qa1 + (self.Qa1/ifthenelse(self.Qa1 + self.Ea1 + self.Fa1 > 0 , self.Qa1 + self.Ea1 + self.Fa1 , 1)) * self.Sa_diff + self.Fa = ifthenelse (self.Fa1 > 0, self.Fa1 + (self.Fa1/ifthenelse(self.Qa1 + self.Ea1 + self.Fa1 > 0 , self.Qa1 + self.Ea1 + self.Fa1 , 1)) * self.Sa_diff, self.Fa1) + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Ea - self.Qa - self.Fa + self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k]) + self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) + + self.Capa = min(self.cap[k] * (1 - self.Sa[k] / self.samax[k]), self.Su[k]) + self.Sa[k] = self.Sa[k] + self.Capa + + self.Su[k] = self.Su_t[k] + self.Fa - self.Capa + self.Perc = self.perc[k] * (self.Su[k] / self.sumax[k]) + self.Su[k] = self.Su[k] - self.Perc + + self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Qaadd - self.Fa + self. Capa - self.Sa[k] + self.Sa_t[k] + self.wbSu_[k] = self.Fa - self.Perc - self. Capa - self.Su[k] + self.Su_t[k] + + self.Eu_[k] = self.Ea + self.Qu_[k] = self.Qa + self.Qaadd + self.Fa_[k] = self.Fa + self.Cap_[k] = self.Cap + self.Perc_[k] = self.Perc +# self.Su_diff_[k] = self.Su_diff +# self.Quadd_[k] = self.Qaadd + +def unsatZone_withAgri_Jarvis(self,k): + """ + - Potential evaporation is calculated with formula in 'JarvisCoefficients', but without + using the Jarvis stress functions + - Potential evaporation is decreased by energy used for interception evaporation + - Formula for evaporation linear until LP, from than with potential rate + - 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 + - Qu is determined with a beta function (same as in HBV?) + - Code for ini-file: 15 + """ + + self.Sa[k] = ifthenelse(self.Sa_t[k] + self.Pe > self.samax[k], self.samax[k], self.Sa_t[k] + self.Pe) + self.Qaadd = ifthenelse(self.Sa_t[k] + self.Pe > self.samax[k], self.Sa_t[k] + self.Pe - self.samax[k], 0) + self.SaN = self.Sa[k] / self.samax[k] + + JarvisCoefficients.calcEu(self,k,1) #calculation of Eu based on Jarvis stress functions + self.Ea = self.Eu + + self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k]) + self.Fa1 = self.famax[k] * (self.sumax[k] - self.Su[k]) / self.sumax[k] + + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Qa1 - self.Ea - self.Fa1 + + self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0) + self.Qa = self.Qa1 + (self.Qa1/ifthenelse(self.Qa1 + self.Fa1 > 0 , self.Qa1 + self.Fa1 , 1)) * self.Sa_diff + self.Fa = ifthenelse (self.Fa1 > 0, self.Fa1 + (self.Fa1/ifthenelse(self.Qa1 + self.Fa1 > 0 , self.Qa1 + self.Fa1 , 1)) * self.Sa_diff, self.Fa1) + self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Ea - self.Qa - self.Fa + self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k]) + self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) + + self.Capa = min(self.cap[k] * (1 - self.Sa[k] / self.samax[k]), self.Su[k]) + self.Sa[k] = self.Sa[k] + self.Capa + + self.Su[k] = self.Su_t[k] + self.Fa - self.Capa + self.Perc = self.perc[k] * (self.Su[k] / self.sumax[k]) + self.Su[k] = self.Su[k] - self.Perc + + self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Qaadd - self.Fa + self. Capa - self.Sa[k] + self.Sa_t[k] + self.wbSu_[k] = self.Fa - self.Perc - self. Capa - self.Su[k] + self.Su_t[k] + + self.Eu_[k] = self.Ea + self.Qu_[k] = self.Qa + self.Qaadd + self.Fa_[k] = self.Fa + self.Cap_[k] = self.Cap + self.Perc_[k] = self.Perc \ No newline at end of file Index: wflow-py/wflow/reservoir_Sw.py =================================================================== diff -u -r978510e581d205d904b07bc305a5d0550779c581 -rb0a0649fbc4d5ed289c09f5beffbd40a711e95ca --- wflow-py/wflow/reservoir_Sw.py (.../reservoir_Sw.py) (revision 978510e581d205d904b07bc305a5d0550779c581) +++ wflow-py/wflow/reservoir_Sw.py (.../reservoir_Sw.py) (revision b0a0649fbc4d5ed289c09f5beffbd40a711e95ca) @@ -25,7 +25,15 @@ import scipy import JarvisCoefficients +def selectSwR(i): + if i == 1: + name = 'snow' + if i == 2: + name = 'snowHour' + + return name + def snow_no_reservoir(self, k): """ This function is used when no snow zone reservoir is used and only @@ -45,8 +53,11 @@ def snow(self,k): """ - snow melt based on degree day factor and + - + - Code for ini-file: 1 """ JarvisCoefficients.calcEpSnow(self,k) + self.PotEvaporation = self.EpHour self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0) self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow @@ -72,9 +83,12 @@ """ - snow melt based on degree day factor and minimum surface temperature - meltfactor increases with temperature + - + - Code for ini-file: 6 """ JarvisCoefficients.calcEpSnow(self,k) + #self.PotEvaporation = self.EpHour self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0) self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow @@ -101,6 +115,8 @@ """ - snow melt based on degree day factor and minimum surface temperature - meltfactor increases with temperature + - + - Code for ini-file: 6 """ JarvisCoefficients.calcEpSnowHour(self,k) @@ -130,9 +146,12 @@ """ - snow melt based on degree day factor and minimum surface temperature - meltfactor increases with temperature + - + - Code for ini-file: 3 """ JarvisCoefficients.calcEpSnow(self,k) + #self.PotEvaporation = self.EpHour self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0) self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow @@ -159,8 +178,11 @@ """ - snow melt based on degree day factor and minimum surface temperature - meltfactor increases with temperature + - + - Code for ini-file: 4 """ JarvisCoefficients.calcEpSnow(self,k) + #self.PotEvaporation = self.EpHour self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0) self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow @@ -184,3 +206,67 @@ self.Ew_[k] = self.Ew self.Qw_[k] = self.Qw +def snow_rain_Tsurf_noEw(self,k): + """ + - snow melt based on degree day factor and minimum surface temperature + - meltfactor increases with temperature + - + - Code for ini-file: 5 + """ + JarvisCoefficients.calcEpSnow(self,k) + #self.PotEvaporation = self.EpHour + self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0) + + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow + + self.Fm2 = max(self.Fm[k] * self.Precipitation, self.Fm[k]) + self.Ew1 = 0 + self.Qw1 = max(min(self.Fm2 * (self.TempSurf - self.Tm[k]), self.Sw[k]), 0) + + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1 + + self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0) + self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff + self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw + self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k]) + self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) + + self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k] + + self.Ew_[k] = self.Ew + self.Qw_[k] = self.Qw + +def snowHour(self,k): + """ + - snow melt based on degree day factor and + - for hourly input data + - Code for ini-file: 2 + """ +# JarvisCoefficients.calcEpSnowHour(self,k) +# self.PotEvaporation = self.EpHour +# self.PotEvaporation = ifthenelse(self.EpHour > 0, self.EpHour, 0) + + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow + + self.Ew1 = max(min(self.PotEvaporation,self.Sw[k]),0) +# self.Ew1 = 0 + self.Qw1 = max(self.Fm[k] * (self.Temperature - self.Tm[k]), 0) + + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1 + + self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0) + self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff + self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw + self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k]) + self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) + +# if any(pcr2numpy(self.Sw[k],nan) > 0): +# pdb.set_trace() + self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k] + + self.Ew_[k] = self.Ew + self.Qw_[k] = self.Qw +# if any(pcr2numpy(self.Qw,nan) > 0): +# pdb.set_trace() Index: wflow-py/wflow/wflow_topoflex.py =================================================================== diff -u -r08acee7cae427d82254c9bb8eb62a8eb88353dda -rb0a0649fbc4d5ed289c09f5beffbd40a711e95ca --- wflow-py/wflow/wflow_topoflex.py (.../wflow_topoflex.py) (revision 08acee7cae427d82254c9bb8eb62a8eb88353dda) +++ wflow-py/wflow/wflow_topoflex.py (.../wflow_topoflex.py) (revision b0a0649fbc4d5ed289c09f5beffbd40a711e95ca) @@ -352,6 +352,8 @@ # static maps to use (normally default) wflow_subcatch = configget(self.config, "model", "wflow_subcatch", "staticmaps/wflow_subcatch.map") + wflow_catchArea = configget(self.config, + "model", "wflow_subcatch", "staticmaps/wflow_catchmentAreas.map") wflow_dem = configget(self.config, "model", "wflow_dem", "staticmaps/wflow_dem.map") wflow_ldd = configget(self.config, @@ -387,6 +389,7 @@ self.TopoLdd = readmap(os.path.join(self.Dir, wflow_ldd)) #: The local drinage definition map (ldd) self.TopoId = readmap( os.path.join(self.Dir, wflow_subcatch)) #: Map define the area over which the calculations are done (mask) + self.catchArea = scalar(ifthen(self.TopoId > 0, scalar(1))) self.LandUse=ordinal(self.wf_readmap(os.path.join(self.Dir , wflow_landuse),0.0,fail=True))#: Map with lan-use/cover classes self.LandUse=cover(self.LandUse,ordinal(ordinal(subcatch) > 0)) self.Soil=ordinal(self.wf_readmap(os.path.join(self.Dir , wflow_soil),0.0,fail=True))#: Map with soil classes @@ -407,9 +410,9 @@ self.Tf = eval(str(configget(self.config, "model", "Tf", "[0]"))) self.Tfa = eval(str(configget(self.config, "model", "Tfa", "[0]"))) - # MODEL PARAMETERS - BASED ON TABLES - self.imax = [self.readtblDefault2(self.Dir + "/" + self.intbl + "/imax" + self.NamesClasses[i] + ".tbl",self.LandUse,subcatch,self.Soil,1.5) for i in self.Classes] - self.sumax = [self.readtblDefault2(self.Dir + "/" + self.intbl + "/sumax" + self.NamesClasses[i] + ".tbl",self.LandUse,subcatch,self.Soil,50) for i in self.Classes] + # MODEL PARAMETERS - BASED ON TABLES + self.imax = [self.readtblDefault2(self.Dir + "/" + self.intbl + "/imax" + self.NamesClasses[i] + ".tbl",self.LandUse,subcatch,self.Soil,1.5) for i in self.Classes] + self.sumax = [self.readtblDefault2(self.Dir + "/" + self.intbl + "/sumax" + self.NamesClasses[i] + ".tbl",self.LandUse,subcatch,self.Soil,70) for i in self.Classes] self.samax = [self.readtblDefault2(self.Dir + "/" + self.intbl + "/samax" + self.NamesClasses[i] + ".tbl",self.LandUse,subcatch,self.Soil,50) for i in self.Classes] self.samin = [self.readtblDefault2(self.Dir + "/" + self.intbl + "/samin" + self.NamesClasses[i] + ".tbl",self.LandUse,subcatch,self.Soil,0.1) for i in self.Classes] self.beta = [self.readtblDefault2(self.Dir + "/" + self.intbl + "/beta" + self.NamesClasses[i] + ".tbl",self.LandUse,subcatch,self.Soil,0.4) for i in self.Classes] @@ -436,11 +439,11 @@ self.lamdaS = eval(str(configget(self.config, "model", "lamdaS", "[0]"))) # initialise list for routing - self.trackQ = [0 * scalar(self.TopoId)] * int(self.maxTransit) + self.trackQ = [0 * scalar(self.catchArea)] * int(self.maxTransit) # 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.convQu = [[0 * scalar(self.catchArea)] * self.Tf[i] for i in self.Classes] + self.convQa = [[0 * scalar(self.catchArea)] * self.Tfa[i] for i in self.Classes] if self.scalarInput: self.gaugesMap = nominal(readmap(os.path.join(self.Dir, @@ -471,7 +474,7 @@ self.sumrunoff = self.ZeroMap # accumulated runoff for water balance (weigthted for upstream area) self.sumpotevap = self.ZeroMap # accumulated runoff for water balance self.sumtemp = self.ZeroMap # accumulated runoff for water balance - self.Q = self.ZeroMap + self.Q = self.ZeroMap self.sumwb = self.ZeroMap # Define timeseries outputs There seems to be a bug and the .tss files are @@ -499,13 +502,14 @@ self.Sf = [self.ZeroMap] * len(self.Classes) self.Sfa = [self.ZeroMap] * len(self.Classes) self.Ss = self.ZeroMap # for combined gw reservoir - self.Qstate = self.ZeroMap # for combined gw reservoir - self.Qstate_t = self.ZeroMap + self.Qstate = self.catchArea * 0 # for combined gw reservoir + self.Qstate_t = self.catchArea * 0 # set initial storage values - self.Sa = [0.05 * self.samax[i] * scalar(self.TopoId) for i in self.Classes] - self.Su = [self.sumax[i] * scalar(self.TopoId) for i in self.Classes] - self.Ss = self.Ss + 30 * scalar(self.TopoId) # for combined gw reservoir +# pdb.set_trace() + self.Sa = [0.05 * self.samax[i] * scalar(self.catchArea) for i in self.Classes] + self.Su = [self.sumax[i] * scalar(self.catchArea) for i in self.Classes] + self.Ss = self.Ss + 30 * scalar(self.catchArea) # for combined gw reservoir else: # self.wf_resume(self.Dir + "/instate/") @@ -606,9 +610,12 @@ # TODO: change rainfall .tss files into grids self.wf_updateparameters() # read the temperature map for each step (see parameters()) - self.logger.debug("Step: "+str(int(self.thestep + self._d_firstTimeStep))+"/"+str(int(self._d_nrTimeSteps))) + # self.logger.debug("Step: "+str(int(self.thestep + self._d_firstTimeStep))+"/"+str(int(self._d_nrTimeSteps))) self.thestep = self.thestep + 1 + #if self.thestep == 26: +# pdb.set_trace() + self.Si_t = copylist(self.Si) self.Sw_t = copylist(self.Sw) self.Su_t = copylist(self.Su) @@ -633,6 +640,9 @@ self.EpDay2 = self.EpDay self.EpDaySnow2 = self.EpDaySnow + #if self.thestep >= 45: + #pdb.set_trace() + for k in self.Classes: # SNOW ================================================================================================= @@ -684,10 +694,10 @@ 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)]) # SLOW RUNOFF RESERVOIR =========================================================================== - if self.selectSs: - eval_str = 'reservoir_Ss.{:s}(self)'.format(self.selectSs) - else: - eval_str = 'reservoir_Ss.groundWater_no_reservoir(self)' + if self.selectSs: + eval_str = 'reservoir_Ss.{:s}(self)'.format(self.selectSs) + else: + eval_str = 'reservoir_Ss.groundWater_no_reservoir(self)' eval(eval_str) # ROUTING @@ -745,9 +755,13 @@ 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.QstateWB = areatotal(sum(self.Qstate) * 0.0405, nominal(self.TopoId)) +# self.Qstate_WB = areatotal(sum(self.Qstate_t) * 0.0405, nominal(self.TopoId)) +# self.QstateWB = areatotal(self.Qstate, nominal(self.TopoId)) +# self.Qstate_WB = areatotal(self.Qstate_t, nominal(self.TopoId)) +# #WBtot in m3/s - self.WBtot = (self.P - self.Ei + self.EwiCorr - self.Ew - 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.SsWB + self.Ss_WB - self.convQuWB +self.convQu_WB - self.convQaWB +self.convQa_WB - self.trackQWB + self.trackQ_WB) / self.timestepsecs + self.WBtot = (self.P - self.Ei + self.EwiCorr - self.Ew - 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.SsWB + self.Ss_WB - self.convQuWB +self.convQu_WB - self.convQaWB +self.convQa_WB - self.trackQWB + self.trackQ_WB - self.QstateWB + self.Qstate_WB) / self.timestepsecs # SUMMED FLUXES ====================================================================================== self.sumprecip = self.sumprecip + self.Precipitation # accumulated rainfall for water balance (m/h)