Index: doc/plots/glacier-plot.py
===================================================================
diff -u
--- doc/plots/glacier-plot.py	(revision 0)
+++ doc/plots/glacier-plot.py	(revision 911009b83b9c498347ba719f6a3bbec960014691)
@@ -0,0 +1,22 @@
+# make plot of soil depth function for different values of M
+from pylab import *
+import wflow.wflow_lib as wl
+
+thetaS = 0.6
+Zi = arange(0,1000,5)
+Ksat = 100
+
+Snow = arange(8000,8700,10)
+
+#MM = [1541.0 * -7]
+
+ToIce = wl.sCurve(Snow,a=8300.,c=0.06)
+
+
+alot =Snow>8300
+#ToIce[alot] = Snow[alot] - 8300
+
+plot(Snow,ToIce)
+xlabel("Snow depth")
+ylabel("Fraction above 8300")
+show()
Index: doc/wflow_sbm.rst
===================================================================
diff -u -r531317c453fb999265e9b1dabd69dfd6532bd4e1 -r911009b83b9c498347ba719f6a3bbec960014691
--- doc/wflow_sbm.rst	(.../wflow_sbm.rst)	(revision 531317c453fb999265e9b1dabd69dfd6532bd4e1)
+++ doc/wflow_sbm.rst	(.../wflow_sbm.rst)	(revision 911009b83b9c498347ba719f6a3bbec960014691)
@@ -141,7 +141,7 @@
 
 Glaciers
 --------
-If snow modeling is enabled Glacier modelling can be enables by including the following three entries
+If snow modeling is enabled Glacier modelling can also be enabled by including the following three entries
 in the modelparameters section:
 
 ::
@@ -152,7 +152,7 @@
 
 
 *GlacierFrac* is a map that gives the fraction of each grid cell covered by a glacier as a number between zero and one.
-Furthermore two lokup tables must be defined: *G\_TT* and *G\_Cfmax*. If the air temperature,
+Furthermore two lookup tables must be defined: *G\_TT* and *G\_Cfmax*. If the air temperature,
 :math:`T_{a}`, is below  :math:`G\_TT (\approx0^{o}C)`
 precipitation occurs as snowfall, whereas it occurs as rainfall if
 :math:`T_{a}\geq G\_TT`.
@@ -165,14 +165,15 @@
 
 where :math:`Q_{m}` is the rate of glacier meltand $cfmax$ is the melting factor in :math:`mm/(^{o}C*day)`.
 
-Accumulated snow on top of the glacier is converted to ice (and will thus become part of the glacier mass) if
-the total snow depth > 8300 mm. If the snow pack is < 8300 mm the following equation is used to convert part of the snow pack
-to ice:
+Accumulated snow on top of the glacier is converted to ice (and will thus become part of the glacier store) if
+the total snow depth > 8300 mm. An S-curve is used to smooth the transition. A maximum of 8mm/day can be converted to ice.
 
-::
+.. plot:: plots/glacier-plot.py
 
-    ToIce = Snow/4150 * 2.0 * timestepsecs/basetimestep
 
+Becuase the glacier store (GlacierStore.map) cannot be initialized by running the model a couple of year a default
+initial state map should be supplied by placing a GlacierStore.map file in the staticmaps directory. This map gives
+the amount of water (in mm) within the Glaciers at each gridcell)
 
 The rainfall interception model
 -------------------------------
Index: wflow-py/wflow/wflow_lib.py
===================================================================
diff -u -rbda446ad8d5ec510a514a29575822e732ddc217d -r911009b83b9c498347ba719f6a3bbec960014691
--- wflow-py/wflow/wflow_lib.py	(.../wflow_lib.py)	(revision bda446ad8d5ec510a514a29575822e732ddc217d)
+++ wflow-py/wflow/wflow_lib.py	(.../wflow_lib.py)	(revision 911009b83b9c498347ba719f6a3bbec960014691)
@@ -784,7 +784,10 @@
     Output:
         - result
     """
-    s = 1.0/(b + exp(-c * (X-a)))
+    try:
+        s = 1.0/(b + exp(-c * (X-a)))
+    except:
+        s = 1.0 / (b + np.exp(-c * (X - a)))
     return s
 
 def sCurveSlope(X,a=0.0,b=1.0,c=1.0):
Index: wflow-py/wflow/wflow_sbm.py
===================================================================
diff -u -rf660893440d321b9250258a349112048127247df -r911009b83b9c498347ba719f6a3bbec960014691
--- wflow-py/wflow/wflow_sbm.py	(.../wflow_sbm.py)	(revision f660893440d321b9250258a349112048127247df)
+++ wflow-py/wflow/wflow_sbm.py	(.../wflow_sbm.py)	(revision 911009b83b9c498347ba719f6a3bbec960014691)
@@ -882,8 +882,10 @@
             self.SnowWater = self.ZeroMap
             self.TSoil = self.ZeroMap + 10.0
             self.CanopyStorage = self.ZeroMap
-            self.ReservoirVolume = self.ResMaxVolume * self.ResTargetFullFrac
-            self.GlacierStore = self.wf_readmap(os.path.join(self.Dir,"staticmaps","GlacierStore.map"), 55.0 * 1000)
+            if hasattr(self, 'ReserVoirLocs'):
+                self.ReservoirVolume = self.ResMaxVolume * self.ResTargetFullFrac
+            if hasattr(self, 'GlacierFrac'):
+                self.GlacierStore = self.wf_readmap(os.path.join(self.Dir,"staticmaps","GlacierStore.map"), 55.0 * 1000)
         else:
             self.logger.info("Setting initial conditions from state files")
             self.wf_resume(os.path.join(self.Dir,"instate"))
@@ -1035,10 +1037,12 @@
                 If below that a max amount of 2mm/day can be converted to glacier-ice
                 """
                 #TODO: document glacier module
-                self.Snow2Glacier = ifthenelse(self.Snow > 8300, self.Snow - 8300 , \
-                                               self.Snow/4150 * 2.0 * self.timestepsecs/self.basetimestep )
+                self.snowdist = sCurve(self.Snow,a=8300.,c=0.06)
+                self.Snow2Glacier = ifthenelse(self.Snow > 8300, self.snowdist * (self.Snow - 8300), self.ZeroMap)
 
                 self.Snow2Glacier = ifthenelse(self.GlacierFrac > 0.0, self.Snow2Glacier,self.ZeroMap)
+                # Max conversion to 8mm/day
+                self.Snow2Glacier = min(self.Snow2Glacier,8.0) * self.timestepsecs/self.basetimestep
 
                 self.Snow = self.Snow - (self.Snow2Glacier * self.GlacierFrac)