Index: wflow-py/wflow/reservoir_Sa.py
===================================================================
diff -u -r978510e581d205d904b07bc305a5d0550779c581 -rb0a0649fbc4d5ed289c09f5beffbd40a711e95ca
--- wflow-py/wflow/reservoir_Sa.py	(.../reservoir_Sa.py)	(revision 978510e581d205d904b07bc305a5d0550779c581)
+++ wflow-py/wflow/reservoir_Sa.py	(.../reservoir_Sa.py)	(revision b0a0649fbc4d5ed289c09f5beffbd40a711e95ca)
@@ -25,6 +25,38 @@
 import scipy
 import JarvisCoefficients
 
+def selectSaR(i):
+    """
+    not all functions are still in this file, the older functions can be found
+    (with the same numbering) in h:\My Documents\memo's\python scripts\wflow\
+    """
+    if i == 1:
+        name = 'agriZone_Jarvis'
+    elif i == 2:
+        name = 'agriZone_Ep'
+    elif i == 3:
+        name = 'agriZone_Ep_Sa'
+    elif i == 4:
+        name = 'agriZone_Ep_Sa_cropG'
+    elif i == 5:
+        name = 'agriZone_Ep_Sa_cropG_beta'
+    elif i == 6:        
+        name = 'agriZone_Ep_Sa_beta'
+    elif i == 7:
+        name = 'agriZone_Ep_Sa_beta_frost'
+    elif i == 8:
+        name = 'agriZone_Ep_Sa_beta_Fvar'
+    elif i == 9:
+        name = 'agriZone_hourlyEp_Sa_beta_Fvar'
+    elif i == 10:
+        name = 'agriZone_hourlyEp_Sa_beta_frost'
+    elif i == 11:
+        name = 'agriZone_hourlyEp_Sa_beta_frostSamax'
+    elif i == 12:
+        name = 'agriZone_Ep_Sa_beta_frostSamax'
+    elif i == 13:
+        name = 'agriZone_Ep_Sa_beta_frostSamax_surfTemp'
+    return name
 
 def agriZone_no_reservoir(self, k):
     """
@@ -42,7 +74,155 @@
     self.Fa_[k] = max(self.Pe_[k], 0)
     self.wbSa_[k] = self.Pe_[k] - self.Ea_[k] - self.Qa_[k] - self.Fa_[k] - self.Sa[k] + self.Sa_t[k]
       
+def agriZone_Jarvis(self,k):
+    """
+    - Potential evaporation is decreased by energy used for interception evaporation    
+    - Formula for evaporation based on Jarvis stress functions
+    - Outgoing fluxes are determined based on (value in previous timestep + inflow) 
+    and if this leads to negative storage, the outgoing fluxes are corrected to rato --> Eu is 
+    no longer taken into account for this correction
+    - Qa u is determined from overflow from Sa
+    - Code for ini-file: 1
+    """
+    self.Qa = max(self.Pe - (self.samax[k] - self.Sa_t[k]),0)    
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa)
+    self.SaN = min(self.Sa[k] / self.samax2, 1)
+    self.SuN = self.Su[k] / self.sumax[k]
     
+    JarvisCoefficients.calcEu(self,k,1)           #calculation of Ea based on Jarvis stress functions
+    self.Ea1 = self.Eu
+    
+    self.Fa1 = self.Fmin[k] + (self.Fmax[k] - self.Fmin[k]) * e ** (-self.decF[k] * self.SuN)
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa) - self.Fa1 - self.Ea1
+
+    self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0)
+    self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 > 0 , self.Fa1 + self.Ea1 , 1)) * self.Sa_diff
+    self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 > 0 , self.Fa1 + self.Ea1 , 1)) * self.Sa_diff
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa) - self.Ea - self.Fa   
+    self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k])    
+    self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) 
+
+    self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Fa - self.Sa[k] + self.Sa_t[k]
+    
+    self.Ea_[k] = self.Ea
+    self.Qa_[k] = self.Qa
+    self.Fa_[k] = self.Fa
+    
+
+def agriZone_Ep(self,k):
+    """
+    - Potential evaporation is decreased by energy used for interception evaporation    
+    - Formula for evaporation based on LP
+    - Outgoing fluxes are determined based on (value in previous timestep + inflow) 
+    and if this leads to negative storage, the outgoing fluxes are corrected to rato --> Eu is 
+    no longer taken into account for this correction
+    - Qa u is determined from overflow from Sa
+    - Code for ini-file: 2
+    """
+    JarvisCoefficients.calcEp(self,k)
+    self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0)  
+    
+    self.Qa = max(self.Pe - (self.samax[k] - self.Sa_t[k]),0)    
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa)
+    self.SaN = min(self.Sa[k] / self.samax2, 1)
+    self.SuN = self.Su[k] / self.sumax[k]
+    
+    self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax[k] * self.LP[k]),1)        
+    
+    self.Fa1 = self.Fmin[k] + (self.Fmax[k] - self.Fmin[k]) * e ** (-self.decF[k] * self.SuN)
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa) - self.Fa1 - self.Ea1
+
+    self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0)
+    self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 > 0 , self.Fa1 + self.Ea1 , 1)) * self.Sa_diff
+    self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 > 0 , self.Fa1 + self.Ea1 , 1)) * self.Sa_diff
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa) - self.Ea - self.Fa   
+    self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k])    
+    self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) 
+
+    self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Fa - self.Sa[k] + self.Sa_t[k]
+    
+    self.Ea_[k] = self.Ea
+    self.Qa_[k] = self.Qa
+    self.Fa_[k] = self.Fa
+
+def agriZone_Ep_Sa(self,k):
+    """
+    - Potential evaporation is decreased by energy used for interception evaporation    
+    - Formula for evaporation based on LP
+    - Outgoing fluxes are determined based on (value in previous timestep + inflow) 
+    and if this leads to negative storage, the outgoing fluxes are corrected to rato --> Eu is 
+    no longer taken into account for this correction
+    - Qa u is determined from overflow from Sa
+    - Fa is based on storage in Sa
+    - Code for ini-file: 3
+    """
+    JarvisCoefficients.calcEp(self,k)
+    self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0)   
+    
+    self.Qa = max(self.Pe - (self.samax[k] - self.Sa_t[k]),0)    
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa)
+    self.SaN = min(self.Sa[k] / self.samax2, 1)
+    self.SuN = self.Su[k] / self.sumax[k]
+    
+    self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax[k] * self.LP[k]),1)        
+    
+    self.Fa1 = ifthenelse(self.SaN > 0,self.Fmin[k] + (self.Fmax[k] - self.Fmin[k]) * e ** (-self.decF[k] * (1 - self.SaN)),0)
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa) - self.Fa1 - self.Ea1
+
+    self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0)
+    self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 > 0 , self.Fa1 + self.Ea1 , 1)) * self.Sa_diff
+    self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 > 0 , self.Fa1 + self.Ea1 , 1)) * self.Sa_diff
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa) - self.Ea - self.Fa   
+    self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k])    
+    self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) 
+
+    self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Fa - self.Sa[k] + self.Sa_t[k]
+    
+    self.Ea_[k] = self.Ea
+    self.Qa_[k] = self.Qa
+    self.Fa_[k] = self.Fa
+    
+    
+def agriZone_Ep_Sa_cropG(self,k):
+    """
+    - Potential evaporation is decreased by energy used for interception evaporation    
+    - Formula for evaporation based on LP
+    - Outgoing fluxes are determined based on (value in previous timestep + inflow) 
+    and if this leads to negative storage, the outgoing fluxes are corrected to rato --> Eu is 
+    no longer taken into account for this correction
+    - Qa u is determined from overflow from Sa
+    - Fa is based on storage in Sa
+    - Code for ini-file: 4
+    """
+    JarvisCoefficients.calcEp(self,k)
+    self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0)   
+    
+    self.samax2 = self.samax[k] * self.cropG   
+    self.Qaadd = max(self.Sa_t[k] - self.samax2,0)
+    
+    self.Qa = max(self.Pe - (self.samax2 - self.Sa_t[k]),0) + self.Qaadd
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa)
+    self.SaN = min(self.Sa[k] / self.samax2, 1)
+    self.SuN = self.Su[k] / self.sumax[k]
+    
+    self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1)        
+    
+    self.Fa1 = ifthenelse(self.SaN > 0,self.Fmin[k] + (self.Fmax[k] - self.Fmin[k]) * e ** (-self.decF[k] * (1 - self.SaN)),0)
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa) - self.Fa1 - self.Ea1
+
+    self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0)
+    self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 > 0 , self.Fa1 + self.Ea1 , 1)) * self.Sa_diff
+    self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 > 0 , self.Fa1 + self.Ea1 , 1)) * self.Sa_diff
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qa) - self.Ea - self.Fa   
+    self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k])    
+    self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) 
+
+    self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Fa - self.Sa[k] + self.Sa_t[k]
+    
+    self.Ea_[k] = self.Ea
+    self.Qa_[k] = self.Qa
+    self.Fa_[k] = self.Fa
+    
 def agriZone_Ep_Sa_cropG_beta(self,k):
     """
     - Potential evaporation is decreased by energy used for interception evaporation    
@@ -52,6 +232,7 @@
     no longer taken into account for this correction
     - Qa u is determined from overflow from Sa --> incorporation of beta function
     - Fa is based on storage in Sa
+    - Code for ini-file: 5
     """
 
     JarvisCoefficients.calcEp(self,k)
@@ -61,7 +242,7 @@
     self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0)
     
     self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd)
-    self.SaN = self.Sa[k] / self.samax2
+    self.SaN = min(self.Sa[k] / self.samax2, 1)
     self.SuN = self.Su[k] / self.sumax[k]
     
     self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1)        
@@ -93,16 +274,17 @@
     no longer taken into account for this correction
     - Qa u is determined from overflow from Sa --> incorporation of beta function
     - Fa is based on storage in Sa
+    - Code for ini-file: 6
     """
 
     JarvisCoefficients.calcEp(self,k)
     self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0)
     
-    self.samax2 = self.samax[k] * scalar(self.TopoId)
+    self.samax2 = self.samax[k] * scalar(self.catchArea)
     self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0)
     
     self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd)
-    self.SaN = max(self.Sa[k] / self.samax2, 1)
+    self.SaN = min(max(self.Sa[k] / self.samax2, 0), 1)
     self.SuN = self.Su[k] / self.sumax[k]
     
     self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1)        
@@ -124,7 +306,52 @@
     self.Qa_[k] = self.Qa + self.Qaadd
     self.Fa_[k] = self.Fa
     
-   
+def agriZone_Ep_Sa_beta_frost(self,k):
+    """
+    - Potential evaporation is decreased by energy used for interception evaporation    
+    - Formula for evaporation based on LP
+    - Outgoing fluxes are determined based on (value in previous timestep + inflow) 
+    and if this leads to negative storage, the outgoing fluxes are corrected to rato --> Eu is 
+    no longer taken into account for this correction
+    - Qa u is determined from overflow from Sa --> incorporation of beta function
+    - Fa is based on storage in Sa
+    - Fa is decreased in case of frozen soil
+    - Code for ini-file: 7
+    """
+
+    JarvisCoefficients.calcEp(self,k)
+    self.PotEvaporation = self.EpHour    
+    
+    self.samax2 = self.samax[k] * scalar(self.catchArea)
+    self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0)
+    self.FrDur[k] = min(self.FrDur[k] + (self.Tmean - 273.15) / 86400 * self.timestepsecs * self.dayDeg[k], 0)
+    
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd)
+    self.SaN = min(self.Sa[k] / self.samax2, 1)
+    self.SuN = self.Su[k] / self.sumax[k]
+    
+    self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1)        
+    self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k])
+    self.Ft = min(max(self.FrDur[k] / (self.FrDur1[k] - self.FrDur0[k]) - self.FrDur0[k] / (self.FrDur1[k] - self.FrDur0[k]), 0), 1)
+    self.Fa1 = self.Ft * ifthenelse(self.SaN > 0,self.Fmin[k] + (self.Fmax[k] - self.Fmin[k]) * e ** (-self.decF[k] * (1 - self.SaN)),0)
+    
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Qa1 - self.Fa1 - self.Ea1
+
+    self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0)
+    self.Qa = self.Qa1 + (self.Qa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0, self.Fa1 + self.Ea1 + self.Qa1, 1)) * self.Sa_diff    
+    self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0, self.Fa1 + self.Ea1 + self.Qa1, 1)) * self.Sa_diff
+    self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0, self.Fa1 + self.Ea1 + self.Qa1, 1)) * self.Sa_diff
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Ea - self.Fa - self.Qa
+    self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k])    
+    self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) 
+
+    self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Qaadd - self.Fa - self.Sa[k] + self.Sa_t[k]
+    
+    self.Ea_[k] = self.Ea
+    self.Qa_[k] = self.Qa + self.Qaadd
+    self.Fa_[k] = self.Fa
+    self.Ft_[k] = self.Ft
+    
 def agriZone_hourlyEp_Sa_beta_frost(self,k):
     """
     - Potential evaporation is decreased by energy used for interception evaporation    
@@ -135,14 +362,18 @@
     - Qa u is determined from overflow from Sa --> incorporation of beta function
     - Fa is based on storage in Sa
     - Fa is decreased in case of frozen soil
+    - Code for ini-file: 10
     """
+
+    #JarvisCoefficients.calcEp(self,k)
+    #self.PotEvaporation = self.EpHour    
     
-    self.samax2 = self.samax[k] * scalar(self.TopoId)
+    self.samax2 = self.samax[k] * scalar(self.catchArea)
     self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0)
     self.FrDur[k] = min(self.FrDur[k] + (self.Temperature) / 86400 * self.timestepsecs * self.dayDeg[k], 0)
     
     self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd)
-    self.SaN = self.Sa[k] / self.samax2
+    self.SaN = min(self.Sa[k] / self.samax2, 1)
     self.SuN = self.Su[k] / self.sumax[k]
     
     self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1)        
@@ -176,12 +407,18 @@
     no longer taken into account for this correction
     - Qa u is determined from overflow from Sa --> incorporation of beta function
     - Fa is based on storage in Sa
-    """ 
+    - Fa is decreased in case of frozen soil
+    - Code for ini-file: 11
+    """
+
+    #JarvisCoefficients.calcEp(self,k)
+    #self.PotEvaporation = self.EpHour    
     
-    self.FrDur[k] = min(self.FrDur[k] + ifthenelse(self.Temperature > 0, self.ratFT[k] * self.Temperature, self.Temperature) * self.dayDeg[k], 0)
-    self.Ft = min(max(self.FrDur[k] / (self.FrDur1[k] - self.FrDur0[k]) - self.FrDur0[k] / (self.FrDur1[k] - self.FrDur0[k]), self.samin[k]), 1)
+    self.FrDur[k] = min(self.FrDur[k] + (self.Temperature) * self.dayDeg[k], 0)    
+    self.Ft = min(max(self.FrDur[k] / (self.FrDur1[k] - self.FrDur0[k]) - self.FrDur0[k] / (self.FrDur1[k] - self.FrDur0[k]), 0.1), 1)
     
-    self.samax2 = self.samax[k] * scalar(self.TopoId) * self.Ft
+    
+    self.samax2 = self.samax[k] * scalar(self.catchArea) * self.Ft
     self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0)
     
     self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd)
@@ -219,6 +456,7 @@
     - Qa u is determined from overflow from Sa --> incorporation of beta function
     - Fa is based on storage in Sa
     - Fa is decreased in case of frozen soil
+    - Code for ini-file: 12
     """
 
     JarvisCoefficients.calcEp(self,k)
@@ -227,12 +465,12 @@
     self.FrDur[k] = min(self.FrDur[k] + ifthenelse(self.Temperature > 0, self.ratFT[k] * self.Temperature, self.Temperature) * self.dayDeg[k], 0)
     self.Ft = min(max(self.FrDur[k] / (self.FrDur1[k] - self.FrDur0[k]) - self.FrDur0[k] / (self.FrDur1[k] - self.FrDur0[k]), self.samin[k]), 1)
     
-    self.samax2 = self.samax[k] * scalar(self.TopoId) * self.Ft
+    self.samax2 = self.samax[k] * scalar(self.catchArea) * self.Ft
     self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0)
     
     self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd)
     self.SaN = min(self.Sa[k] / self.samax2, 1)
-    self.SuN = self.Su[k] / self.sumax
+    self.SuN = self.Su[k] / self.sumax[k]
     
     self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1)        
     self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k])
@@ -266,19 +504,20 @@
     - Qa u is determined from overflow from Sa --> incorporation of beta function
     - Fa is based on storage in Sa
     - Fa is decreased in case of frozen soil
+    - Code for ini-file: 13
     """
 
     JarvisCoefficients.calcEp(self,k)
-    self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0)      
+    self.PotEvaporation = self.EpHour    
      
     self.FrDur[k] = min(self.FrDur[k] + ifthenelse(self.TempSurf > 0, self.ratFT[k] * self.TempSurf, self.TempSurf) * self.dayDeg[k], 0)
     self.Ft = min(max(self.FrDur[k] / (self.FrDur1[k] - self.FrDur0[k]) - self.FrDur0[k] / (self.FrDur1[k] - self.FrDur0[k]), self.samin[k]), 1)
     
-    self.samax2 = self.samax[k] * scalar(self.TopoId) * self.Ft
+    self.samax2 = self.samax[k] * scalar(self.catchArea) * self.Ft
     self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0)
     
     self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd)
-    self.SaN = self.Sa[k] / self.samax2
+    self.SaN = min(self.Sa[k] / self.samax2, 1)
     self.SuN = self.Su[k] / self.sumax[k]
     
     self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1)        
@@ -303,3 +542,84 @@
     self.Fa_[k] = self.Fa
     self.Ft_[k] = self.Ft
 
+def agriZone_Ep_Sa_beta_Fvar(self,k):
+    """
+    - Potential evaporation is decreased by energy used for interception evaporation    
+    - Formula for evaporation based on LP
+    - Outgoing fluxes are determined based on (value in previous timestep + inflow) 
+    and if this leads to negative storage, the outgoing fluxes are corrected to rato --> Eu is 
+    no longer taken into account for this correction
+    - Qa u is determined from overflow from Sa --> incorporation of beta function
+    - Fa is based on storage in Sa
+    - Code for ini-file: 8
+    """
+
+    JarvisCoefficients.calcEp(self,k)
+    self.PotEvaporation = self.EpHour    
+    
+    self.samax2 = self.samax[k] * scalar(self.catchArea)
+    self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0)
+    
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd)
+    self.SaN = min(self.Sa[k] / self.samax2, 1)
+    self.SuN = self.Su[k] / self.sumax[k]
+    
+    self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1)        
+    self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k])
+    self.Fa1 = self.cropG * ifthenelse(self.SaN > 0,self.Fmin[k] + (self.Fmax[k] - self.Fmin[k]) * e ** (-self.decF[k] * (1 - self.SaN)),0)
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Qa1 - self.Fa1 - self.Ea1
+
+    self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0)
+    self.Qa = self.Qa1 + (self.Qa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff    
+    self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff
+    self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Ea - self.Fa - self.Qa
+    self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k])    
+    self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) 
+
+    self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Qaadd - self.Fa - self.Sa[k] + self.Sa_t[k]
+    
+    self.Ea_[k] = self.Ea
+    self.Qa_[k] = self.Qa + self.Qaadd
+    self.Fa_[k] = self.Fa
+    
+def agriZone_hourlyEp_Sa_beta_Fvar(self,k):
+    """
+    - Potential evaporation is decreased by energy used for interception evaporation    
+    - Formula for evaporation based on LP
+    - Outgoing fluxes are determined based on (value in previous timestep + inflow) 
+    and if this leads to negative storage, the outgoing fluxes are corrected to rato --> Eu is 
+    no longer taken into account for this correction
+    - Qa u is determined from overflow from Sa --> incorporation of beta function
+    - Fa is based on storage in Sa
+    - Code for ini-file: 9
+    """
+
+#    JarvisCoefficients.calcEp(self,k)
+#    self.PotEvaporation = self.EpHour    
+    
+    self.samax2 = self.samax[k] * scalar(self.catchArea)
+    self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0)
+    
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd)
+    self.SaN = min(self.Sa[k] / self.samax2, 1)
+    self.SuN = self.Su[k] / self.sumax[k]
+    
+    self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1)        
+    self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k])
+    self.Fa1 = self.cropG * ifthenelse(self.SaN > 0,self.Fmin[k] + (self.Fmax[k] - self.Fmin[k]) * e ** (-self.decF[k] * (1 - self.SaN)),0)
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Qa1 - self.Fa1 - self.Ea1
+
+    self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0)
+    self.Qa = self.Qa1 + (self.Qa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff    
+    self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff
+    self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Ea - self.Fa - self.Qa
+    self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k])    
+    self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) 
+
+    self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Qaadd - self.Fa - self.Sa[k] + self.Sa_t[k]
+    
+    self.Ea_[k] = self.Ea
+    self.Qa_[k] = self.Qa + self.Qaadd
+    self.Fa_[k] = self.Fa
\ No newline at end of file