Index: wflow-py/wflow/reservoir_Ss.py
===================================================================
diff -u -rfacc3a7f729f6910dac37d4761992c70d238fcbd -r9aacf70d29eac315b5782685f8b1ec1bbc91c984
--- wflow-py/wflow/reservoir_Ss.py	(.../reservoir_Ss.py)	(revision facc3a7f729f6910dac37d4761992c70d238fcbd)
+++ wflow-py/wflow/reservoir_Ss.py	(.../reservoir_Ss.py)	(revision 9aacf70d29eac315b5782685f8b1ec1bbc91c984)
@@ -51,13 +51,18 @@
     self.Qsin = areatotal(self.QsinTemp * self.percentArea ,nominal(self.TopoId))  #areatotal is taken, according to area percentage of cell
  
     self.Qs = self.Ss * self.Ks[0]    
-    self.Ss = self.Ss * exp(-self.Ks[0]) + self.Qsin
+    # self.Ss = self.Ss * exp(-self.Ks[0]) + self.Qsin
+
+    # add a gain/lossterm based on constant value
+    self.Gain = ifthenelse(self.Closure < 0.0 , - min(self.Ss, -self.Closure) , self.Closure) # if negative, limited by storage in Ss
+    self.Ss = self.Ss * exp(-self.Ks[0]) + self.Qsin + self.Gain
     
-    self.wbSs = self.Qsin - self.Qs - self.Ss + self.Ss_t    
+    self.wbSs = self.Qsin - self.Qs - self.Ss + self.Ss_t + self.Gain
     
     self.Qs_ = self.Qs
     self.Recharge_ = self.Recharge
     self.QsinClass_ = self.QsinClass
     self.QsinTemp_ = self.QsinTemp
     self.Qsin_ = self.Qsin
+    self.Gain_ = self.Gain
     
Index: wflow-py/wflow/wflow_topoflex.py
===================================================================
diff -u -rad5e049f1419c25bc30c20d2d0e1bcc56dcb4381 -r9aacf70d29eac315b5782685f8b1ec1bbc91c984
--- wflow-py/wflow/wflow_topoflex.py	(.../wflow_topoflex.py)	(revision ad5e049f1419c25bc30c20d2d0e1bcc56dcb4381)
+++ wflow-py/wflow/wflow_topoflex.py	(.../wflow_topoflex.py)	(revision 9aacf70d29eac315b5782685f8b1ec1bbc91c984)
@@ -375,6 +375,8 @@
                                    "model", "wflow_subcatch", "staticmaps/wflow_catchmentAreas.map")
         wflow_dem = configget(self.config,
                               "model", "wflow_dem", "staticmaps/wflow_dem.map")
+        wflow_maxSlope = configget(self.config,
+                              "model", "wflow_maxSlope", "staticmaps/wflow_maxSlope.map")
         wflow_ldd = configget(self.config,
                               "model", "wflow_ldd", "staticmaps/wflow_ldd.map")
         wflow_landuse  = configget(self.config,
@@ -409,6 +411,7 @@
 
         self.Altitude = readmap(os.path.join(self.Dir, wflow_dem)) * scalar(
             defined(subcatch))  #: The digital elevation map (DEM)
+        self.maxSlope = self.wf_readmap(os.path.join(self.Dir, wflow_maxSlope),0.0)
         self.TopoLdd = readmap(os.path.join(self.Dir, wflow_ldd))  #: The local drinage definition map (ldd)
         self.TopoId = readmap(
             os.path.join(self.Dir, wflow_subcatch))  #: Map define the area over which the calculations are done (mask)
@@ -464,6 +467,7 @@
         self.Tm = [self.readtblDefault2(self.Dir + "/" + self.intbl + "/Tm" + self.NamesClasses[i] + ".tbl",self.LandUse,subcatch,self.Soil,2) for i in self.Classes]
         self.Fm = [self.readtblDefault2(self.Dir + "/" + self.intbl + "/Fm" + self.NamesClasses[i] + ".tbl",self.LandUse,subcatch,self.Soil,0.2) for i in self.Classes]
         self.ECORR= self.readtblDefault2(self.Dir + "/" + self.intbl + "/ECORR.tbl",self.LandUse,subcatch,self.Soil, 1.0)
+        self.Closure = self.readtblDefault2(self.Dir + "/" + self.intbl + "/Closure.tbl",self.LandUse,subcatch,self.Soil, 0.0)
 
         #kinematic wave parameters
         self.Beta = scalar(0.6) # For sheetflow
@@ -478,7 +482,7 @@
         self.lamdaS = eval(str(configget(self.config, "model", "lamdaS", "[0]")))
 
         # initialise list for routing
-        self.trackQ = [0 * scalar(self.catchArea)] * int(self.maxTransit)
+        self.trackQ = [0 * scalar(self.catchArea)] * int(self.maxTransit) # list * scalar ---> list wordt zoveel x gekopieerd als scalar.
 
         # initialise list for lag function
         self.convQu = [[0 * scalar(self.catchArea)] * self.Tf[i] for i in self.Classes]
@@ -496,6 +500,7 @@
         self.xl,self.yl,self.reallength = pcrut.detRealCellLength(self.ZeroMap,sizeinmetres)
         self.Slope= slope(self.Altitude)
         self.Slope=ifthen(boolean(self.TopoId),max(0.001,self.Slope*celllength()/self.reallength))
+        self.Slope=ifthenelse(self.maxSlope>0.0, self.maxSlope, self.Slope)
         Terrain_angle=scalar(atan(self.Slope))
         temp = catchmenttotal(cover(1.0), self.TopoLdd) * self.reallength * 0.001 * 0.001 *  self.reallength
         self.QMMConvUp = cover(self.timestepsecs * 0.001)/temp