Index: wflow-py/wflow/reservoir_Sw.py
===================================================================
diff -u -r368a19849555f5d349fdb349010530b3e7207d95 -r978510e581d205d904b07bc305a5d0550779c581
--- wflow-py/wflow/reservoir_Sw.py	(.../reservoir_Sw.py)	(revision 368a19849555f5d349fdb349010530b3e7207d95)
+++ wflow-py/wflow/reservoir_Sw.py	(.../reservoir_Sw.py)	(revision 978510e581d205d904b07bc305a5d0550779c581)
@@ -15,6 +15,7 @@
 """
 
 import numpy
+import pdb
 from copy import copy as copylist
 
 try:
@@ -24,15 +25,7 @@
 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
@@ -52,11 +45,8 @@
 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
@@ -77,23 +67,50 @@
     
     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
+    """
+   
+    JarvisCoefficients.calcEpSnow(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
 
-def snowHour(self,k):
+    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 
-    - for hourly input  data
-    - Code for ini-file: 2
+    - snow melt based on degree day factor and minimum surface temperature
+    - meltfactor increases with temperature
     """
-#    JarvisCoefficients.calcEpSnowHour(self,k)
-#    self.PotEvaporation = self.EpHour
-#    self.PotEvaporation = ifthenelse(self.EpHour > 0, self.EpHour, 0)
+   
+    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.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.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
 
@@ -104,11 +121,66 @@
     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()
+
+def snow_rain_Tsurf(self,k):
+    """
+    - snow melt based on degree day factor and minimum surface temperature
+    - meltfactor increases with temperature
+    """
+   
+    JarvisCoefficients.calcEpSnow(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.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
+    """
+    JarvisCoefficients.calcEpSnow(self,k)
+    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
+