Index: wflow-py/wflow/reservoir_Su.py =================================================================== diff -u -rb0a0649fbc4d5ed289c09f5beffbd40a711e95ca -r36b469410a4a5af24c3902cbebb595e9da83c7bb --- wflow-py/wflow/reservoir_Su.py (.../reservoir_Su.py) (revision b0a0649fbc4d5ed289c09f5beffbd40a711e95ca) +++ wflow-py/wflow/reservoir_Su.py (.../reservoir_Su.py) (revision 36b469410a4a5af24c3902cbebb595e9da83c7bb) @@ -12,59 +12,61 @@ from copy import copy as copylist try: - from wflow.wf_DynamicFramework import * + from wflow.wf_DynamicFramework import * except ImportError: - from wf_DynamicFramework import * + from wf_DynamicFramework import * 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' + name = "unsatZone_LP_beta" elif i == 2: - name = 'unsatZone_LP' + name = "unsatZone_LP" elif i == 3: - name = 'unsatZone_SF' + name = "unsatZone_SF" elif i == 4: - name = 'unsatZone_EisEp' + name = "unsatZone_EisEp" elif i == 5: - name = 'unsatZone_drain' + name = "unsatZone_drain" elif i == 6: - name = 'unsatZone_drain_split' + 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' + 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' + name = "unsatZone_LP_beta_Jarvis" + elif i == 13: + name = "unsatZone_LP_beta_Ep" elif i == 14: - name = 'unsatZone_withAgri_Ep' + name = "unsatZone_withAgri_Ep" elif i == 15: - name = 'unsatZone_withAgri_Jarvis' + name = "unsatZone_withAgri_Jarvis" elif i == 16: - name = 'unsatZone_forAgri_Jarvis' + name = "unsatZone_forAgri_Jarvis" elif i == 17: - name = 'unsatZone_forAgri_Ep' + name = "unsatZone_forAgri_Ep" elif i == 18: - name = 'unsatZone_forAgri_Jarvis_cropG' + name = "unsatZone_forAgri_Jarvis_cropG" elif i == 19: - name = 'unsatZone_forAgri_Ep_cropG' - elif i == 20: - name = 'unsatZone_LP_beta_Ep_cropG' + 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 @@ -80,9 +82,12 @@ self.Perc_[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] + self.wbSu_[k] = ( + self.Pe - self.Eu - self.Qu - self.Perc + self.Cap - self.Su[k] + self.Su_t[k] + ) -def unsatZone_LP_beta(self,k): + +def unsatZone_LP_beta(self, k): """ - Potential evaporation is decreased by energy used for interception evaporation - Formula for evaporation linear until LP, from than with potential rate @@ -91,39 +96,95 @@ - 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.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.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[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.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.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 - - -def unsatZone_LP_beta_Jarvis(self,k): + + +def unsatZone_LP_beta_Jarvis(self, k): """ - Potential evaporation is decreased by energy used for interception evaporation - Formula for evaporation based on Jarvis stress functions @@ -133,38 +194,67 @@ - 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) + 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] - - JarvisCoefficients.calcEu(self,k,1) #calculation of Eu based on Jarvis stress functions - + self.SiN = self.Si[k] / self.imax[k] + + JarvisCoefficients.calcEu( + self, k, 1 + ) # calculation of Eu based on Jarvis stress functions + 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.Eu - self.Perc1 self.Su_diff = ifthenelse(self.Su[k] < 0, self.Su[k], 0) - self.Qu = self.Qu1 + (self.Qu1/ifthenelse(self.Qu1 + self.Perc1 > 0 , self.Qu1 + self.Perc1 , 1)) * self.Su_diff - self.Perc = ifthenelse (self.Perc1 > 0, self.Perc1 + (self.Perc1/ifthenelse(self.Qu1 + self.Perc1 > 0 , self.Qu1 + 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.Qu = ( + self.Qu1 + + (self.Qu1 / ifthenelse(self.Qu1 + self.Perc1 > 0, self.Qu1 + self.Perc1, 1)) + * self.Su_diff + ) + self.Perc = ifthenelse( + self.Perc1 > 0, + self.Perc1 + + (self.Perc1 / ifthenelse(self.Qu1 + self.Perc1 > 0, self.Qu1 + 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.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 - - -def unsatZone_LP_beta_Ep(self,k): + + +def unsatZone_LP_beta_Ep(self, k): """ - Potential evaporation is calculated with formula in 'JarvisCoefficients', but without using the Jarvis stress functions @@ -175,44 +265,99 @@ - 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) + + # 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.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[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.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.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.Su_diff_[k] = self.Su_diff + # self.Quadd_[k] = self.Quadd self.Epot_[k] = self.PotEvaporation -def unsatZone_LP_beta_Ep_Ei(self,k): + +def unsatZone_LP_beta_Ep_Ei(self, k): """ - Potential evaporation is calculated with formula in 'JarvisCoefficients', but without using the Jarvis stress functions @@ -223,45 +368,103 @@ - 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) + + # 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 = 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.SiN = self.Si[k] / self.imax[k] + + 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 - self.Su[k] = self.Su_t[k] + (self.Pe - self.Quadd) - self.Qu1 - self.Eu1 - self.Perc1 + 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.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.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.Su_diff_[k] = self.Su_diff + # self.Quadd_[k] = self.Quadd self.Epot_[k] = self.PotEvaporation -def unsatZone_LP_beta_Ep_percD(self,k): + +def unsatZone_LP_beta_Ep_percD(self, k): """ - Potential evaporation is calculated with formula in 'JarvisCoefficients', but without using the Jarvis stress functions @@ -272,42 +475,97 @@ - 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) + 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.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[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.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.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): + + +def unsatZone_LP_beta_Ep_percD(self, k): """ - Potential evaporation is calculated with formula in 'JarvisCoefficients', but without using the Jarvis stress functions @@ -318,43 +576,97 @@ - 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) + 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.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[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.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.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): +def unsatZone_LP_beta_Ep_percDvar(self, k): """ - Potential evaporation is calculated with formula in 'JarvisCoefficients', but without using the Jarvis stress functions @@ -365,47 +677,120 @@ - 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) + + 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.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.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 = 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[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.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.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): + +def unsatZone_LP_beta_Ep_cropG(self, k): """ - Potential evaporation is calculated with formula in 'JarvisCoefficients', but without using the Jarvis stress functions @@ -417,50 +802,107 @@ - root zone storage for crop land is decreased in autumn and winter - Code for ini-file: 20 """ - - JarvisCoefficients.calcEp(self,k) - self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0) - self.cropG_scal = pcr2numpy(self.cropG,NaN) + JarvisCoefficients.calcEp(self, k) + self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0), 0) + + self.cropG_scal = pcr2numpy(self.cropG, NaN) if any(self.cropG_scal == 1): self.sumax2 = self.sumax[k] - 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])) + 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] - - self.Su[k] = ifthenelse(self.Su_t[k] + self.Pe > self.sumax2, self.sumax2, self.Su_t[k] + self.Pe) - self.Quadd = ifthenelse(self.Su_t[k] + self.Pe > self.sumax2, self.Su_t[k] + self.Pe - self.sumax2, 0) + + self.Su[k] = ifthenelse( + self.Su_t[k] + self.Pe > self.sumax2, self.sumax2, self.Su_t[k] + self.Pe + ) + self.Quadd = ifthenelse( + self.Su_t[k] + self.Pe > self.sumax2, self.Su_t[k] + self.Pe - self.sumax2, 0 + ) self.SuN = self.Su[k] / self.sumax2 - self.SiN = self.Si[k] / self.imax[k] - - self.Eu1 = max((self.PotEvaporation - self.Ei),0) * min(self.Su[k] / (self.sumax2 * self.LP[k]),1) - + self.SiN = self.Si[k] / self.imax[k] + + self.Eu1 = max((self.PotEvaporation - self.Ei), 0) * min( + self.Su[k] / (self.sumax2 * 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[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.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.sumax2), 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.Cap = min(self.cap[k] * (1 - self.Su[k] / self.sumax2), 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 -def unsatZone_forAgri_Jarvis(self,k): + +def unsatZone_forAgri_Jarvis(self, k): """ - Potential evaporation is decreased by energy used for interception evaporation - Formula for evaporation based on Jarvis stress functions @@ -471,37 +913,90 @@ - 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) + 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] - - JarvisCoefficients.calcEu(self,k,2) #calculation of Eu based on Jarvis stress functions + self.SiN = self.Si[k] / self.imax[k] + + JarvisCoefficients.calcEu( + self, k, 2 + ) # calculation of Eu based on Jarvis stress functions self.Eu1 = self.Eu - + self.Qu1 = (self.Fa - self.Quadd) * (1 - (1 - self.SuN) ** self.beta[k]) self.Perc1 = self.perc[k] * self.SuN 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.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.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(self,k): + +def unsatZone_forAgri_Ep(self, k): """ - Potential evaporation is decreased by energy used for interception evaporation - Formula for evaporation based on beta/LP @@ -512,40 +1007,96 @@ - 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) + 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.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.Qu1 = (self.Fa - self.Quadd) * (1 - (1 - self.SuN) ** self.beta[k]) self.Perc1 = self.perc[k] * self.SuN 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.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.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_percD(self,k): + +def unsatZone_forAgri_Ep_percD(self, k): """ - Potential evaporation is decreased by energy used for interception evaporation - Formula for evaporation based on beta/LP @@ -556,41 +1107,97 @@ - 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) + 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.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.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.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): + +def unsatZone_forAgri_Ep_percDvar(self, k): """ - Potential evaporation is decreased by energy used for interception evaporation - Formula for evaporation based on beta/LP @@ -601,45 +1208,119 @@ - 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) + 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.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.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.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): + +def unsatZone_forAgri_hourlyEp(self, k): """ - Potential evaporation is decreased by energy used for interception evaporation - Formula for evaporation based on beta/LP @@ -650,38 +1331,93 @@ - 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) - self.Quadd = ifthenelse(self.Su_t[k] + self.Fa > self.sumax[k], self.Su_t[k] + self.Fa - self.sumax[k], 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 31 mrt 2016 - + 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 31 mrt 2016 + self.Qu1 = (self.Fa - self.Quadd) * (1 - (1 - self.SuN) ** self.beta[k]) self.Perc1 = self.perc[k] * self.SuN 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.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.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_Jarvis_cropG(self,k): +def unsatZone_forAgri_Jarvis_cropG(self, k): """ - Potential evaporation is decreased by energy used for interception evaporation - Formula for evaporation based on Jarvis stress functions @@ -692,43 +1428,94 @@ - inflow is infiltration from agriculture reservoir - Code for ini-file: 18 """ - self.cropG_scal = pcr2numpy(self.cropG,NaN) + self.cropG_scal = pcr2numpy(self.cropG, NaN) if any(self.cropG_scal == 1): self.sumax2 = self.sumax[k] else: self.sumax2 = self.sumax[k] * self.redsu[k] - - self.Su[k] = ifthenelse(self.Su_t[k] + self.Fa > self.sumax2, self.sumax2, self.Su_t[k] + self.Fa) - self.Quadd = ifthenelse(self.Su_t[k] + self.Fa > self.sumax2, self.Su_t[k] + self.Fa - self.sumax2, 0) + + self.Su[k] = ifthenelse( + self.Su_t[k] + self.Fa > self.sumax2, self.sumax2, self.Su_t[k] + self.Fa + ) + self.Quadd = ifthenelse( + self.Su_t[k] + self.Fa > self.sumax2, self.Su_t[k] + self.Fa - self.sumax2, 0 + ) self.SuN = self.Su[k] / self.sumax2 - self.SiN = self.Si[k] / self.imax[k] - - JarvisCoefficients.calcEu(self,k,2) #calculation of Eu based on Jarvis stress functions + self.SiN = self.Si[k] / self.imax[k] + + JarvisCoefficients.calcEu( + self, k, 2 + ) # calculation of Eu based on Jarvis stress functions self.Eu1 = self.Eu - + self.Qu1 = (self.Fa - self.Quadd) * (1 - (1 - self.SuN) ** self.beta[k]) self.Perc1 = self.perc[k] * self.SuN 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.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.sumax2), 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.Cap = min(self.cap[k] * (1 - self.Su[k] / self.sumax2), 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_cropG(self,k): + +def unsatZone_forAgri_Ep_cropG(self, k): """ - Potential evaporation is decreased by energy used for interception evaporation - Formula for evaporation based on beta/LP @@ -739,46 +1526,97 @@ - inflow is infiltration from agriculture reservoir - Code for ini-file: 19 """ - - JarvisCoefficients.calcEp(self,k) - self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0) - - self.cropG_scal = pcr2numpy(self.cropG,NaN) + + JarvisCoefficients.calcEp(self, k) + self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0), 0) + + self.cropG_scal = pcr2numpy(self.cropG, NaN) if any(self.cropG_scal == 1): self.sumax2 = self.sumax[k] else: self.sumax2 = self.sumax[k] * self.redsu[k] - - self.Su[k] = ifthenelse(self.Su_t[k] + self.Fa > self.sumax2, self.sumax2, self.Su_t[k] + self.Fa) - self.Quadd = ifthenelse(self.Su_t[k] + self.Fa > self.sumax2, self.Su_t[k] + self.Fa - self.sumax2, 0) + + self.Su[k] = ifthenelse( + self.Su_t[k] + self.Fa > self.sumax2, self.sumax2, self.Su_t[k] + self.Fa + ) + self.Quadd = ifthenelse( + self.Su_t[k] + self.Fa > self.sumax2, self.Su_t[k] + self.Fa - self.sumax2, 0 + ) self.SuN = self.Su[k] / self.sumax2 - self.SiN = self.Si[k] / self.imax[k] - - self.Eu1 = max((self.PotEvaporation - self.Ei),0) * min(self.Su[k] / (self.sumax2 * self.LP[k]),1) - + self.SiN = self.Si[k] / self.imax[k] + + self.Eu1 = max((self.PotEvaporation - self.Ei), 0) * min( + self.Su[k] / (self.sumax2 * self.LP[k]), 1 + ) + self.Qu1 = (self.Fa - self.Quadd) * (1 - (1 - self.SuN) ** self.beta[k]) self.Perc1 = self.perc[k] * self.SuN 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.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.sumax2), 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.Cap = min(self.cap[k] * (1 - self.Su[k] / self.sumax2), 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_withAgri(self,k): + +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 @@ -787,44 +1625,92 @@ - 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.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.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.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.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.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.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): + + +def unsatZone_withAgri_Ep(self, k): """ - Potential evaporation is calculated with formula in 'JarvisCoefficients', but without using the Jarvis stress functions @@ -836,46 +1722,95 @@ - Code for ini-file: 14 """ - JarvisCoefficients.calcEp(self,k) - self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0) + 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.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.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.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.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.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.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): + + +def unsatZone_withAgri_Jarvis(self, k): """ - Potential evaporation is calculated with formula in 'JarvisCoefficients', but without using the Jarvis stress functions @@ -887,37 +1822,64 @@ - 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.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 - + + 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.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.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.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.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 + self.Perc_[k] = self.Perc