Index: wflow-py/wflow/reservoir_Sa.py
===================================================================
diff -u -rb0a0649fbc4d5ed289c09f5beffbd40a711e95ca -rf02904810f997f1e5a4ff7cc186c1ea3075d01ef
--- wflow-py/wflow/reservoir_Sa.py	(.../reservoir_Sa.py)	(revision b0a0649fbc4d5ed289c09f5beffbd40a711e95ca)
+++ wflow-py/wflow/reservoir_Sa.py	(.../reservoir_Sa.py)	(revision f02904810f997f1e5a4ff7cc186c1ea3075d01ef)
@@ -1,625 +1,666 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Wed Feb 04 14:52:30 2015
-
-@author: teuser
-"""
-
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Apr 03 16:31:35 2014
-
-@author: TEuser
-
-List all function versions
-"""
-
-import numpy
-import pdb
-from copy import copy as copylist
-
-try:
-    from  wflow.wf_DynamicFramework import *
-except ImportError:
-    from  wf_DynamicFramework import *
-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):
-    """
-    This function is used when no unsaturated zone reservoir is used and only
-    passes fluxes from the upper reservoirs to the lower
-    self.Qa_[k] = 0.
-    self.Ea_[k] = 0.
-    self.Sa[k] = 0.
-    self.Fa_[k] = Pe
-    Storage in unsaturated zone = 0.
-    """
-    self.Qa_[k] = 0.
-    self.Ea_[k] = 0.
-    self.Sa[k] = 0.
-    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    
-    - 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: 5
-    """
-
-    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.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 = 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_Ep_Sa_beta(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: 6
-    """
-
-    JarvisCoefficients.calcEp(self,k)
-    self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0)
-    
-    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(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)        
-    self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k])
-    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.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_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    
-    - 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: 10
-    """
-
-    #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.Temperature) / 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_frostSamax(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: 11
-    """
-
-    #JarvisCoefficients.calcEp(self,k)
-    #self.PotEvaporation = self.EpHour    
-    
-    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.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[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 = 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_Ep_Sa_beta_frostSamax(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: 12
-    """
-
-    JarvisCoefficients.calcEp(self,k)
-    self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0)   
-     
-    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.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[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 = 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_Ep_Sa_beta_frostSamax_surfTemp(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: 13
-    """
-
-    JarvisCoefficients.calcEp(self,k)
-    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.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[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 = 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_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
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Feb 04 14:52:30 2015
+
+@author: teuser
+"""
+
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Apr 03 16:31:35 2014
+
+@author: TEuser
+
+List all function versions
+"""
+
+import numpy
+import pdb
+from copy import copy as copylist
+
+try:
+    from  wflow.wf_DynamicFramework import *
+except ImportError:
+    from  wf_DynamicFramework import *
+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):
+    """
+    This function is used when no unsaturated zone reservoir is used and only
+    passes fluxes from the upper reservoirs to the lower
+    self.Qa_[k] = 0.
+    self.Ea_[k] = 0.
+    self.Sa[k] = 0.
+    self.Fa_[k] = Pe
+    Storage in unsaturated zone = 0.
+    """
+    self.Qa_[k] = 0.
+    self.Ea_[k] = 0.
+    self.Sa[k] = 0.
+    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    
+    - 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: 5
+    """
+
+    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.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 = 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_Ep_Sa_beta(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: 6
+    """
+
+    JarvisCoefficients.calcEp(self,k)
+    self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0)
+    
+    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(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)        
+    self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k])
+    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.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(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:
+    """
+
+    #JarvisCoefficients.calcEp(self,k)
+    #self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0)
+    
+    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(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)        
+    self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k])
+    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.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_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    
+    - 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: 10
+    """
+
+    #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.Temperature) / 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_frostSamax(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: 11
+    """
+
+    #JarvisCoefficients.calcEp(self,k)
+    #self.PotEvaporation = self.EpHour    
+    
+    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.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[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 = 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_Ep_Sa_beta_frostSamax(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: 12
+    """
+
+    JarvisCoefficients.calcEp(self,k)
+    self.PotEvaporation = cover(ifthenelse(self.EpHour >= 0, self.EpHour, 0),0)   
+     
+    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.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[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 = 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_Ep_Sa_beta_frostSamax_surfTemp(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: 13
+    """
+
+    JarvisCoefficients.calcEp(self,k)
+    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.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[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 = 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_Ep_Sa_beta_Fvar(self,k):
+    """
+    - Potential evaporation is decreased by energy used for interception evaporation    
+    - Formula for evaporation based on LP
+    - Outgoing fluxes are determined based on (value in previous timestep + inflow) 
+    and if this leads to negative storage, the outgoing fluxes are corrected to rato --> Eu is 
+    no longer taken into account for this correction
+    - Qa u is determined from overflow from Sa --> incorporation of beta function
+    - Fa is based on storage in Sa
+    - Code for ini-file: 8
+    """
+
+    JarvisCoefficients.calcEp(self,k)
+    self.PotEvaporation = self.EpHour    
+    
+    self.samax2 = self.samax[k] * scalar(self.catchArea)
+    self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0)
+    
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd)
+    self.SaN = min(self.Sa[k] / self.samax2, 1)
+    self.SuN = self.Su[k] / self.sumax[k]
+    
+    self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1)        
+    self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k])
+    self.Fa1 = self.cropG * ifthenelse(self.SaN > 0,self.Fmin[k] + (self.Fmax[k] - self.Fmin[k]) * e ** (-self.decF[k] * (1 - self.SaN)),0)
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Qa1 - self.Fa1 - self.Ea1
+
+    self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0)
+    self.Qa = self.Qa1 + (self.Qa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff    
+    self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff
+    self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Ea - self.Fa - self.Qa
+    self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k])    
+    self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) 
+
+    self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Qaadd - self.Fa - self.Sa[k] + self.Sa_t[k]
+    
+    self.Ea_[k] = self.Ea
+    self.Qa_[k] = self.Qa + self.Qaadd
+    self.Fa_[k] = self.Fa
+    
+def agriZone_hourlyEp_Sa_beta_Fvar(self,k):
+    """
+    - Potential evaporation is decreased by energy used for interception evaporation    
+    - Formula for evaporation based on LP
+    - Outgoing fluxes are determined based on (value in previous timestep + inflow) 
+    and if this leads to negative storage, the outgoing fluxes are corrected to rato --> Eu is 
+    no longer taken into account for this correction
+    - Qa u is determined from overflow from Sa --> incorporation of beta function
+    - Fa is based on storage in Sa
+    - Code for ini-file: 9
+    """
+
+#    JarvisCoefficients.calcEp(self,k)
+#    self.PotEvaporation = self.EpHour    
+    
+    self.samax2 = self.samax[k] * scalar(self.catchArea)
+    self.Qaadd = max(self.Sa_t[k] + self.Pe - self.samax2,0)
+    
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd)
+    self.SaN = min(self.Sa[k] / self.samax2, 1)
+    self.SuN = self.Su[k] / self.sumax[k]
+    
+    self.Ea1 = max((self.PotEvaporation - self.Ei),0) * min(self.Sa[k] / (self.samax2 * self.LP[k]),1)        
+    self.Qa1 = (self.Pe - self.Qaadd) * (1 - (1 - self.SaN) ** self.beta[k])
+    self.Fa1 = self.cropG * ifthenelse(self.SaN > 0,self.Fmin[k] + (self.Fmax[k] - self.Fmin[k]) * e ** (-self.decF[k] * (1 - self.SaN)),0)
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Qa1 - self.Fa1 - self.Ea1
+
+    self.Sa_diff = ifthenelse(self.Sa[k] < 0, self.Sa[k], 0)
+    self.Qa = self.Qa1 + (self.Qa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff    
+    self.Fa = self.Fa1 + (self.Fa1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff
+    self.Ea = self.Ea1 + (self.Ea1/ifthenelse(self.Fa1 + self.Ea1 + self.Qa1 > 0 , self.Fa1 + self.Ea1 + self.Qa1 , 1)) * self.Sa_diff
+    self.Sa[k] = self.Sa_t[k] + (self.Pe - self.Qaadd) - self.Ea - self.Fa - self.Qa
+    self.Sa[k] = ifthenelse(self.Sa[k] < 0, 0 , self.Sa[k])    
+    self.Sa_diff2 = ifthen(self.Sa[k] < 0, self.Sa[k]) 
+
+    self.wbSa_[k] = self.Pe - self.Ea - self.Qa - self.Qaadd - self.Fa - self.Sa[k] + self.Sa_t[k]
+    
+    self.Ea_[k] = self.Ea
+    self.Qa_[k] = self.Qa + self.Qaadd
     self.Fa_[k] = self.Fa
\ No newline at end of file
Index: wflow-py/wflow/reservoir_Sw.py
===================================================================
diff -u -rb0a0649fbc4d5ed289c09f5beffbd40a711e95ca -rf02904810f997f1e5a4ff7cc186c1ea3075d01ef
--- wflow-py/wflow/reservoir_Sw.py	(.../reservoir_Sw.py)	(revision b0a0649fbc4d5ed289c09f5beffbd40a711e95ca)
+++ wflow-py/wflow/reservoir_Sw.py	(.../reservoir_Sw.py)	(revision f02904810f997f1e5a4ff7cc186c1ea3075d01ef)
@@ -1,272 +1,279 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Wed Feb 04 14:52:30 2015
-
-@author: teuser
-"""
-
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Apr 03 16:31:35 2014
-
-@author: TEuser
-
-List all function versions
-"""
-
-import numpy
-import pdb
-from copy import copy as copylist
-
-try:
-    from  wflow.wf_DynamicFramework import *
-except ImportError:
-    from  wf_DynamicFramework import *
-import scipy
-import JarvisCoefficients
-
-def selectSwR(i):
-
-    if i == 1:
-        name = 'snow'
-    if i == 2:
-	name = 'snowHour'
-
-    return name
-
-def snow_no_reservoir(self, k):
-    """
-    This function is used when no snow zone reservoir is used and only
-    passes fluxes from the upper reservoirs to the lower
-    Qw = Psnow
-    Ew = 0.
-    Storage in snow zone = 0.
-
-    !!!still needs a final check!!!    
-    
-    """
-    self.Qw_[k] = max(self.PrecipitationSnow, 0)
-    self.Ew_[k] = 0.
-    self.Sw[k] = 0.
-    self.wbSw_[k] = self.Precipitation - self.Ew_[k] - self.Qw_[k] - self.Sw[k] + self.Sw_t[k]  
-      
-def snow(self,k):
-    """
-    - snow melt based on degree day factor and 
-    - 
-    - Code for ini-file: 1
-    """
-    JarvisCoefficients.calcEpSnow(self,k)
-    self.PotEvaporation = self.EpHour
-    self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0)
-    
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow
-    
-    self.Ew1 = max(min(self.PotEvaporation,self.Sw[k]),0)
-    self.Qw1 = max(self.Fm[k] * (self.Temperature - self.Tm[k]), 0)
-    
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1
-
-    self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0)
-    self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
-    self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw
-    self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k])    
-    self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) 
-
-    self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k]
-    
-    self.Ew_[k] = self.Ew
-    self.Qw_[k] = self.Qw
-    
-def snow_rain(self,k):
-    """
-    - snow melt based on degree day factor and minimum surface temperature
-    - meltfactor increases with temperature
-    - 
-    - Code for ini-file: 6
-    """
-   
-    JarvisCoefficients.calcEpSnow(self,k)
-    #self.PotEvaporation = self.EpHour
-    self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0)
-    
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow
-    
-    self.Fm2 = max(self.Fm[k] * self.Precipitation, self.Fm[k])
-    self.Ew1 = max(min(self.PotEvaporation, self.Sw[k]),0)
-    self.Qw1 = max(min(self.Fm2 * (self.Temperature - self.Tm[k]), self.Sw[k]), 0)
-    
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1
-
-    self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0)
-    self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
-    self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw
-    self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k])    
-    self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) 
-
-    self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k]
-    
-    self.Ew_[k] = self.Ew
-    self.Qw_[k] = self.Qw
-
-def snow_rain_hourlyEp(self,k):
-    """
-    - snow melt based on degree day factor and minimum surface temperature
-    - meltfactor increases with temperature
-    - 
-    - Code for ini-file: 6
-    """
-   
-    JarvisCoefficients.calcEpSnowHour(self,k)
-    self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0)
-    
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow
-    
-    self.Fm2 = max(self.Fm[k] * self.Precipitation, self.Fm[k])
-    self.Ew1 = max(min(self.PotEvaporation, self.Sw[k]),0)
-    self.Qw1 = max(min(self.Fm2 * (self.Temperature - self.Tm[k]), self.Sw[k]), 0)
-    
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1
-
-    self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0)
-    self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
-    self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw
-    self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k])    
-    self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) 
-
-    self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k]
-    
-    self.Ew_[k] = self.Ew
-    self.Qw_[k] = self.Qw
-
-def snow_rain_Tsurf(self,k):
-    """
-    - snow melt based on degree day factor and minimum surface temperature
-    - meltfactor increases with temperature
-    - 
-    - Code for ini-file: 3
-    """
-   
-    JarvisCoefficients.calcEpSnow(self,k)
-    #self.PotEvaporation = self.EpHour
-    self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0)
-    
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow
-    
-    self.Fm2 = max(self.Fm[k] * self.Precipitation, self.Fm[k])
-    self.Ew1 = max(min(self.PotEvaporation, self.Sw[k]),0)
-    self.Qw1 = max(min(self.Fm2 * (self.TempSurf - self.Tm[k]), self.Sw[k]), 0)
-    
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1
-
-    self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0)
-    self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
-    self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw
-    self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k])    
-    self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) 
-
-    self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k]
-    
-    self.Ew_[k] = self.Ew
-    self.Qw_[k] = self.Qw
-
-def snow_rain_TsurfAir(self,k):
-    """
-    - snow melt based on degree day factor and minimum surface temperature
-    - meltfactor increases with temperature
-    - 
-    - Code for ini-file: 4
-    """
-    JarvisCoefficients.calcEpSnow(self,k)
-    #self.PotEvaporation = self.EpHour
-    self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0)
-    
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow
-    self.Temp = (self.TempSurf + self.Temperature) / 2
-    
-    self.Fm2 = max(self.Fm[k] * self.Precipitation, self.Fm[k])
-    self.Ew1 = max(min(self.PotEvaporation, self.Sw[k]),0)
-    self.Qw1 = max(min(self.Fm2 * (self.Temp - self.Tm[k]), self.Sw[k]), 0)
-    
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1
-
-    self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0)
-    self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
-    self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw
-    self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k])    
-    self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) 
-
-    self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k]
-    
-    self.Ew_[k] = self.Ew
-    self.Qw_[k] = self.Qw
-
-def snow_rain_Tsurf_noEw(self,k):
-    """
-    - snow melt based on degree day factor and minimum surface temperature
-    - meltfactor increases with temperature
-    - 
-    - Code for ini-file: 5
-    """
-    JarvisCoefficients.calcEpSnow(self,k)
-    #self.PotEvaporation = self.EpHour
-    self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0)
-    
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow
-    
-    self.Fm2 = max(self.Fm[k] * self.Precipitation, self.Fm[k])
-    self.Ew1 = 0
-    self.Qw1 = max(min(self.Fm2 * (self.TempSurf - self.Tm[k]), self.Sw[k]), 0)
-    
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1
-
-    self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0)
-    self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
-    self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw
-    self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k])    
-    self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) 
-
-    self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k]
-    
-    self.Ew_[k] = self.Ew
-    self.Qw_[k] = self.Qw
-
-def snowHour(self,k):
-    """
-    - snow melt based on degree day factor and 
-    - for hourly input  data
-    - Code for ini-file: 2
-    """
-#    JarvisCoefficients.calcEpSnowHour(self,k)
-#    self.PotEvaporation = self.EpHour
-#    self.PotEvaporation = ifthenelse(self.EpHour > 0, self.EpHour, 0)
-    
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow
-    
-    self.Ew1 = max(min(self.PotEvaporation,self.Sw[k]),0)
-#    self.Ew1 = 0
-    self.Qw1 = max(self.Fm[k] * (self.Temperature - self.Tm[k]), 0)
-    
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1
-
-    self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0)
-    self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
-    self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
-    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw
-    self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k])    
-    self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) 
-
-#    if any(pcr2numpy(self.Sw[k],nan) > 0):    
-#        pdb.set_trace()
-    self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k]
-    
-    self.Ew_[k] = self.Ew
-    self.Qw_[k] = self.Qw
-#    if any(pcr2numpy(self.Qw,nan) > 0):    
-#        pdb.set_trace()
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Feb 04 14:52:30 2015
+
+@author: teuser
+"""
+
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Apr 03 16:31:35 2014
+
+@author: TEuser
+
+List all function versions
+"""
+
+import numpy
+import pdb
+from copy import copy as copylist
+
+try:
+    from  wflow.wf_DynamicFramework import *
+except ImportError:
+    from  wf_DynamicFramework import *
+import scipy
+import JarvisCoefficients
+
+def selectSwR(i):
+
+    if i == 1:
+        name = 'snow'
+    if i == 2:
+	name = 'snowHour'
+
+    return name
+
+def snow_no_reservoir(self, k):
+    """
+    This function is used when no snow zone reservoir is used and only
+    passes fluxes from the upper reservoirs to the lower
+    Qw = Psnow
+    Ew = 0.
+    Storage in snow zone = 0.
+
+    !!!still needs a final check!!!    
+    
+    """
+    try:
+    	JarvisCoefficients.calcEpSnow(self,k)
+    except:
+    	JarvisCoefficients.calcEpSnowHour(self,k)
+    self.PotEvaporation = self.EpHour
+    self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0)
+    
+    self.Qw_[k] = max(self.PrecipitationSnow, 0)
+    self.Ew_[k] = 0.
+    self.Sw[k] = 0.
+    self.wbSw_[k] = self.Precipitation - self.Ew_[k] - self.Qw_[k] - self.Sw[k] + self.Sw_t[k]  
+      
+def snow(self,k):
+    """
+    - snow melt based on degree day factor and 
+    - 
+    - Code for ini-file: 1
+    """
+    JarvisCoefficients.calcEpSnow(self,k)
+    self.PotEvaporation = self.EpHour
+    self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0)
+    
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow
+    
+    self.Ew1 = max(min(self.PotEvaporation,self.Sw[k]),0)
+    self.Qw1 = max(self.Fm[k] * (self.Temperature - self.Tm[k]), 0)
+    
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1
+
+    self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0)
+    self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
+    self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw
+    self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k])    
+    self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) 
+
+    self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k]
+    
+    self.Ew_[k] = self.Ew
+    self.Qw_[k] = self.Qw
+    
+def snow_rain(self,k):
+    """
+    - snow melt based on degree day factor and minimum surface temperature
+    - meltfactor increases with temperature
+    - 
+    - Code for ini-file: 6
+    """
+   
+    JarvisCoefficients.calcEpSnow(self,k)
+    #self.PotEvaporation = self.EpHour
+    self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0)
+    
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow
+    
+    self.Fm2 = max(self.Fm[k] * self.Precipitation, self.Fm[k])
+    self.Ew1 = max(min(self.PotEvaporation, self.Sw[k]),0)
+    self.Qw1 = max(min(self.Fm2 * (self.Temperature - self.Tm[k]), self.Sw[k]), 0)
+    
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1
+
+    self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0)
+    self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
+    self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw
+    self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k])    
+    self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) 
+
+    self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k]
+    
+    self.Ew_[k] = self.Ew
+    self.Qw_[k] = self.Qw
+
+def snow_rain_hourlyEp(self,k):
+    """
+    - snow melt based on degree day factor and minimum surface temperature
+    - meltfactor increases with temperature
+    - 
+    - Code for ini-file: 6
+    """
+   
+    JarvisCoefficients.calcEpSnowHour(self,k)
+    self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0)
+    
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow
+    
+    self.Fm2 = max(self.Fm[k] * self.Precipitation, self.Fm[k])
+    self.Ew1 = max(min(self.PotEvaporation, self.Sw[k]),0)
+    self.Qw1 = max(min(self.Fm2 * (self.Temperature - self.Tm[k]), self.Sw[k]), 0)
+    
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1
+
+    self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0)
+    self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
+    self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw
+    self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k])    
+    self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) 
+
+    self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k]
+    
+    self.Ew_[k] = self.Ew
+    self.Qw_[k] = self.Qw
+
+def snow_rain_Tsurf(self,k):
+    """
+    - snow melt based on degree day factor and minimum surface temperature
+    - meltfactor increases with temperature
+    - 
+    - Code for ini-file: 3
+    """
+   
+    JarvisCoefficients.calcEpSnow(self,k)
+    #self.PotEvaporation = self.EpHour
+    self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0)
+    
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow
+    
+    self.Fm2 = max(self.Fm[k] * self.Precipitation, self.Fm[k])
+    self.Ew1 = max(min(self.PotEvaporation, self.Sw[k]),0)
+    self.Qw1 = max(min(self.Fm2 * (self.TempSurf - self.Tm[k]), self.Sw[k]), 0)
+    
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1
+
+    self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0)
+    self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
+    self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw
+    self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k])    
+    self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) 
+
+    self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k]
+    
+    self.Ew_[k] = self.Ew
+    self.Qw_[k] = self.Qw
+
+def snow_rain_TsurfAir(self,k):
+    """
+    - snow melt based on degree day factor and minimum surface temperature
+    - meltfactor increases with temperature
+    - 
+    - Code for ini-file: 4
+    """
+    JarvisCoefficients.calcEpSnow(self,k)
+    #self.PotEvaporation = self.EpHour
+    self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0)
+    
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow
+    self.Temp = (self.TempSurf + self.Temperature) / 2
+    
+    self.Fm2 = max(self.Fm[k] * self.Precipitation, self.Fm[k])
+    self.Ew1 = max(min(self.PotEvaporation, self.Sw[k]),0)
+    self.Qw1 = max(min(self.Fm2 * (self.Temp - self.Tm[k]), self.Sw[k]), 0)
+    
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1
+
+    self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0)
+    self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
+    self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw
+    self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k])    
+    self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) 
+
+    self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k]
+    
+    self.Ew_[k] = self.Ew
+    self.Qw_[k] = self.Qw
+
+def snow_rain_Tsurf_noEw(self,k):
+    """
+    - snow melt based on degree day factor and minimum surface temperature
+    - meltfactor increases with temperature
+    - 
+    - Code for ini-file: 5
+    """
+    JarvisCoefficients.calcEpSnow(self,k)
+    #self.PotEvaporation = self.EpHour
+    self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0)
+    
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow
+    
+    self.Fm2 = max(self.Fm[k] * self.Precipitation, self.Fm[k])
+    self.Ew1 = 0
+    self.Qw1 = max(min(self.Fm2 * (self.TempSurf - self.Tm[k]), self.Sw[k]), 0)
+    
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1
+
+    self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0)
+    self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
+    self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw
+    self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k])    
+    self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) 
+
+    self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k]
+    
+    self.Ew_[k] = self.Ew
+    self.Qw_[k] = self.Qw
+
+def snowHour(self,k):
+    """
+    - snow melt based on degree day factor and 
+    - for hourly input  data
+    - Code for ini-file: 2
+    """
+#    JarvisCoefficients.calcEpSnowHour(self,k)
+#    self.PotEvaporation = self.EpHour
+#    self.PotEvaporation = ifthenelse(self.EpHour > 0, self.EpHour, 0)
+    
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow
+    
+    self.Ew1 = max(min(self.PotEvaporation,self.Sw[k]),0)
+#    self.Ew1 = 0
+    self.Qw1 = max(self.Fm[k] * (self.Temperature - self.Tm[k]), 0)
+    
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1
+
+    self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0)
+    self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
+    self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff
+    self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw
+    self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k])    
+    self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) 
+
+#    if any(pcr2numpy(self.Sw[k],nan) > 0):    
+#        pdb.set_trace()
+    self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k]
+    
+    self.Ew_[k] = self.Ew
+    self.Qw_[k] = self.Qw
+#    if any(pcr2numpy(self.Qw,nan) > 0):    
+#        pdb.set_trace()