Index: wflow-py/Sandbox/wflow_topoflex/reservoir_Sus.py =================================================================== diff -u -r5f72a83356683168a37f127593b6f33c4eccd3b5 -r36b469410a4a5af24c3902cbebb595e9da83c7bb --- wflow-py/Sandbox/wflow_topoflex/reservoir_Sus.py (.../reservoir_Sus.py) (revision 5f72a83356683168a37f127593b6f33c4eccd3b5) +++ wflow-py/Sandbox/wflow_topoflex/reservoir_Sus.py (.../reservoir_Sus.py) (revision 36b469410a4a5af24c3902cbebb595e9da83c7bb) @@ -10,27 +10,28 @@ """ - import numpy import pdb 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 + def selectSusR(i): if i == 1: - name = 'unsatSatZone_GWout' - elif i ==2: - name = 'unsatSatZone_noGWout' + name = "unsatSatZone_GWout" + elif i == 2: + name = "unsatSatZone_noGWout" elif i == 3: - name = 'unsatSatZone_noGWout_VSA' + name = "unsatSatZone_noGWout_VSA" return name -def unsatSatZone_noGWout_VSA(self,k): + +def unsatSatZone_noGWout_VSA(self, k): """ - Potential evaporation is decreased by energy used for interception evaporation - Evaporation rate constrained by moisture stress via LP @@ -41,31 +42,50 @@ area is estimated, linear relation is derived from HAND distribution in wetland) - Code for ini-file: 3 """ - -# pdb.set_trace() - A = min(self.Co[k] * (1 + ((self.Sus_t[k]-self.susmax2[k])*(1-self.Co[k])) / ((self.susmax3[k]-self.susmax2[k])*self.Co[k])),1) - self.Qo = ifthenelse (self.Sus_t[k] >= self.susmax2[k], self.Pe * A,0) + + # pdb.set_trace() + A = min( + self.Co[k] + * ( + 1 + + ((self.Sus_t[k] - self.susmax2[k]) * (1 - self.Co[k])) + / ((self.susmax3[k] - self.susmax2[k]) * self.Co[k]) + ), + 1, + ) + self.Qo = ifthenelse(self.Sus_t[k] >= self.susmax2[k], self.Pe * A, 0) self.Sus[k] = self.Sus_t[k] + (self.Pe - self.Qo) - self.Eu1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sus[k] / (self.susmax2[k] * self.LP[k]),1) + self.Eu1 = max((self.PotEvaporation - self.Ei), 0) * min( + self.Sus[k] / (self.susmax2[k] * self.LP[k]), 1 + ) self.Qd1 = self.Kd[k] * max((self.Sus[k] - self.susmax1[k]), 0) - + self.Sus[k] = self.Sus_t[k] + (self.Pe - self.Qo) - self.Eu1 - self.Qd1 - + self.Su_diff = ifthenelse(self.Sus[k] < 0, self.Sus[k], 0) - self.Eu = self.Eu1 + (self.Eu1 / ifthenelse(self.Qd1 + self.Eu1 > 0 , self.Qd1 + self.Eu1, 1)) * self.Su_diff - self.Qd = self.Qd1 + (self.Qd1/ifthenelse(self.Qd1 + self.Eu1 > 0 , self.Qd1 + self.Eu1, 1)) * self.Su_diff + self.Eu = ( + self.Eu1 + + (self.Eu1 / ifthenelse(self.Qd1 + self.Eu1 > 0, self.Qd1 + self.Eu1, 1)) + * self.Su_diff + ) + self.Qd = ( + self.Qd1 + + (self.Qd1 / ifthenelse(self.Qd1 + self.Eu1 > 0, self.Qd1 + self.Eu1, 1)) + * self.Su_diff + ) self.Sus[k] = self.Sus_t[k] + self.Pe - self.Qo - self.Eu - self.Qd - self.Sus[k] = ifthenelse(self.Sus[k] < 0, 0 , self.Sus[k]) - self.Su_diff2 = ifthen(self.Sus[k] < 0, self.Sus[k]) - - self.wbSus_[k] = self.Pe - self.Qo - self.Qd - self.Eu - self.Sus[k] + self.Sus_t[k] - + self.Sus[k] = ifthenelse(self.Sus[k] < 0, 0, self.Sus[k]) + self.Su_diff2 = ifthen(self.Sus[k] < 0, self.Sus[k]) + + self.wbSus_[k] = self.Pe - self.Qo - self.Qd - self.Eu - self.Sus[k] + self.Sus_t[k] + self.Eu_[k] = self.Eu self.Qo_[k] = self.Qo self.Qd_[k] = self.Qd self.Su_diff_[k] = self.Su_diff -def unsatSatZone_noGWout(self,k): + +def unsatSatZone_noGWout(self, k): """ - Potential evaporation is decreased by energy used for interception evaporation - Evaporation rate constrained by moisture stress via LP @@ -74,29 +94,40 @@ - the saturated part of the reservoir has only drainage outflow - Code for ini-file: 2 """ - - self.Qo = max (self.Pe - (self.susmax2[k] - self.Sus_t[k]), 0) + + self.Qo = max(self.Pe - (self.susmax2[k] - self.Sus_t[k]), 0) self.Sus[k] = self.Sus_t[k] + self.Pe - self.Qo - self.Eu1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sus[k] / (self.susmax2[k] * self.LP[k]),1) + self.Eu1 = max((self.PotEvaporation - self.Ei), 0) * min( + self.Sus[k] / (self.susmax2[k] * self.LP[k]), 1 + ) self.Qd1 = self.Kd[k] * max((self.Sus[k] - self.susmax1[k]), 0) - + self.Sus[k] = self.Sus_t[k] + self.Pe - self.Qo - self.Eu1 - self.Qd1 - + self.Su_diff = ifthenelse(self.Sus[k] < 0, self.Sus[k], 0) - self.Eu = self.Eu1 + (self.Eu1 / ifthenelse(self.Qd1 + self.Eu1 > 0 , self.Qd1 + self.Eu1, 1)) * self.Su_diff - self.Qd = self.Qd1 + (self.Qd1/ifthenelse(self.Qd1 + self.Eu1 > 0 , self.Qd1 + self.Eu1, 1)) * self.Su_diff + self.Eu = ( + self.Eu1 + + (self.Eu1 / ifthenelse(self.Qd1 + self.Eu1 > 0, self.Qd1 + self.Eu1, 1)) + * self.Su_diff + ) + self.Qd = ( + self.Qd1 + + (self.Qd1 / ifthenelse(self.Qd1 + self.Eu1 > 0, self.Qd1 + self.Eu1, 1)) + * self.Su_diff + ) self.Sus[k] = self.Sus_t[k] + self.Pe - self.Qo - self.Eu - self.Qd - self.Sus[k] = ifthenelse(self.Sus[k] < 0, 0 , self.Sus[k]) - self.Su_diff2 = ifthen(self.Sus[k] < 0, self.Sus[k]) - - self.wbSus_[k] = self.Pe - self.Qo - self.Qd - self.Eu - self.Sus[k] + self.Sus_t[k] - + self.Sus[k] = ifthenelse(self.Sus[k] < 0, 0, self.Sus[k]) + self.Su_diff2 = ifthen(self.Sus[k] < 0, self.Sus[k]) + + self.wbSus_[k] = self.Pe - self.Qo - self.Qd - self.Eu - self.Sus[k] + self.Sus_t[k] + self.Eu_[k] = self.Eu self.Qo_[k] = self.Qo self.Qd_[k] = self.Qd self.Su_diff_[k] = self.Su_diff -def unsatSatZone_GWout(self,k): + +def unsatSatZone_GWout(self, k): """ - Potential evaporation is decreased by energy used for interception evaporation - Evaporation rate constrained by moisture stress via LP @@ -106,28 +137,40 @@ - NIET UITGEBREID GETEST VOOR SLUITEN WB, NU OOK CONFLICT MET Qs - Code for ini-file: 1 """ - - self.Qo = max (self.Pe - (self.susmax2[k] - self.Sus_t[k]), 0) + + self.Qo = max(self.Pe - (self.susmax2[k] - self.Sus_t[k]), 0) self.Sus[k] = self.Sus_t[k] + self.Pe - self.Qo - self.Eu1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sus[k] / (self.susmax2[k] * self.LP[k]),1) + self.Eu1 = max((self.PotEvaporation - self.Ei), 0) * min( + self.Sus[k] / (self.susmax2[k] * self.LP[k]), 1 + ) self.Qd1 = self.Kd[k] * max((self.Sus[k] - self.susmax1[k]), 0) - + self.Sus[k] = self.Sus_t[k] + self.Pe - self.Qo - self.Eu1 - self.Qd1 - + self.Su_diff = ifthenelse(self.Sus[k] < 0, self.Sus[k], 0) - self.Eu = self.Eu1 + (self.Eu1 / ifthenelse(self.Qd1 + self.Eu1 > 0 , self.Qd1 + self.Eu1, 1)) * self.Su_diff - self.Qd = self.Qd1 + (self.Qd1/ifthenelse(self.Qd1 + self.Eu1 > 0 , self.Qd1 + self.Eu1, 1)) * self.Su_diff + self.Eu = ( + self.Eu1 + + (self.Eu1 / ifthenelse(self.Qd1 + self.Eu1 > 0, self.Qd1 + self.Eu1, 1)) + * self.Su_diff + ) + self.Qd = ( + self.Qd1 + + (self.Qd1 / ifthenelse(self.Qd1 + self.Eu1 > 0, self.Qd1 + self.Eu1, 1)) + * self.Su_diff + ) self.Sus[k] = self.Sus_t[k] + self.Pe - self.Qo - self.Eu - self.Qd - self.Sus[k] = ifthenelse(self.Sus[k] < 0, 0 , self.Sus[k]) - self.Su_diff2 = ifthen(self.Sus[k] < 0, self.Sus[k]) - - self.Qus = self.Ks[k] * self.Sus[k] + self.Sus[k] = ifthenelse(self.Sus[k] < 0, 0, self.Sus[k]) + self.Su_diff2 = ifthen(self.Sus[k] < 0, self.Sus[k]) + + self.Qus = self.Ks[k] * self.Sus[k] self.Sus[k] = self.Sus[k] - self.Qus - - self.wbSus_[k] = self.Pe - self.Qo - self.Qd - self.Eu - self.Qus - self.Sus + self.Sus_t - + + self.wbSus_[k] = ( + self.Pe - self.Qo - self.Qd - self.Eu - self.Qus - self.Sus + self.Sus_t + ) + self.Eu_[k] = self.Eu self.Qo_[k] = self.Qo self.Qd_[k] = self.Qd self.Qus_[k] = self.Qus - self.Su_diff_[k] = self.Su_diff \ No newline at end of file + self.Su_diff_[k] = self.Su_diff