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