Index: examples/wflow_rhine_sbm/inmaps.nc
===================================================================
diff -u -r01a02941af2650cc7087fd1d02d95159062faad0 -r524657cbf82aea4dc0ef9e8de7d2a48a346a4993
Binary files differ
Index: wflow-py/Scripts/pcr2netcdf.py
===================================================================
diff -u -rf48b5321292915b892132463b750dbda2998e3cb -r524657cbf82aea4dc0ef9e8de7d2a48a346a4993
--- wflow-py/Scripts/pcr2netcdf.py	(.../pcr2netcdf.py)	(revision f48b5321292915b892132463b750dbda2998e3cb)
+++ wflow-py/Scripts/pcr2netcdf.py	(.../pcr2netcdf.py)	(revision 524657cbf82aea4dc0ef9e8de7d2a48a346a4993)
@@ -431,11 +431,10 @@
                 var.append(flst)
                 varname.append(flst)
             else:
-                mapstackname = flst
-                var =flst
-                varname =flst
+                mapstackname = a
+                var =a
+                varname =a
 
-
     # Use first timestep as clone-map
     logger = setlogger('pcr2netcdf.log','pcr2netcdf', thelevel = logging.DEBUG)
 
Index: wflow-py/wflow/wflow_lib.py
===================================================================
diff -u -r8b08f07d1577581c87b804290de74fd24c62340f -r524657cbf82aea4dc0ef9e8de7d2a48a346a4993
--- wflow-py/wflow/wflow_lib.py	(.../wflow_lib.py)	(revision 8b08f07d1577581c87b804290de74fd24c62340f)
+++ wflow-py/wflow/wflow_lib.py	(.../wflow_lib.py)	(revision 524657cbf82aea4dc0ef9e8de7d2a48a346a4993)
@@ -670,31 +670,46 @@
     return ifthen(result >=0,result)
 
 
-def hand(dem,ldd,threshold=50.0,stream=None):
+def derive_HAND(dem, ldd, accuThreshold, rivers=None, basin=None):
     """
-    Determines heigth above nearest drain.
-
+    Function derives Height-Above-Nearest-Drain.
+    See http://www.sciencedirect.com/science/article/pii/S003442570800120X
     Input:
-        ldd
-        dem
-        Optional:
-        threshold - to determine streams (number of upstreams cell needed to be a stream)
-        stream - if the threshols is not a stream map (boolean) should be supplied
-
-    Returns: Hand and the stream map
-
+        dem -- pcraster object float32, elevation data
+        ldd -- pcraster object direction, local drain directions
+        accuThreshold -- upstream amount of cells as threshold for river
+            delineation
+        rivers=None -- you can provide a rivers layer here. Pixels that are 
+                        identified as river should have a value > 0, other
+                        pixels a value of zero.
+        basin=None -- set a boolean pcraster map where areas with True are estimated using the nearest drain in ldd distance
+                        and areas with False by means of the nearest friction distance. Friction distance estimated using the 
+                        upstream area as weight (i.e. drains with a bigger upstream area have a lower friction)
+                        the spreadzone operator is used in this case.
+    Output:
+        hand -- pcraster bject float32, height, normalised to nearest stream
+        dist -- distance to nearest stream measured in cell lengths
+            according to D8 directions
     """
+    if rivers is None:
+        stream = ifthenelse(accuflux(ldd, 1) >= accuThreshold,
+                                boolean(1), boolean(0))
+    else:
+        stream = boolean(cover(rivers, 0))
+    
+    height_river = ifthenelse(stream, ordinal(dem*100), 0)
+    if basin is None:
+        up_elevation = scalar(subcatchment(ldd, height_river))
+    else:
+        drainage_surf = ifthen(rivers, accuflux(ldd, 1))
+        weight = 1./scalar(spreadzone(cover(ordinal(drainage_surf), 0), 0, 0))
+        up_elevation = ifthenelse(basin, scalar(subcatchment(ldd, height_river)), scalar(spreadzone(height_river, 0, weight)))
+        # replace areas outside of basin by a spread zone calculation.
+    hand = max(scalar(ordinal(dem*100))-up_elevation, 0)/100
+    dist = ldddist(ldd, stream, 1)
+        
+    return hand, dist
 
-    #setglobaloption("unittrue")
-    if threshold:
-        stream = boolean(ifthenelse(accuflux(ldd,1) > threshold, boolean(1), boolean(0)))
-
-    demdiff = dem - downstream(ldd,dem)
-    hand = ldddist(ldd,stream,max(demdiff,0))/celllength()
-
-    return hand, stream
-
-
     
 def sCurve(X,a=0.0,b=1.0,c=1.0):
     """
Index: wflow-py/wflow/wflow_sbm.py
===================================================================
diff -u -r3844888de0b2f066d4831fd900634056413c9acd -r524657cbf82aea4dc0ef9e8de7d2a48a346a4993
--- wflow-py/wflow/wflow_sbm.py	(.../wflow_sbm.py)	(revision 3844888de0b2f066d4831fd900634056413c9acd)
+++ wflow-py/wflow/wflow_sbm.py	(.../wflow_sbm.py)	(revision 524657cbf82aea4dc0ef9e8de7d2a48a346a4993)
@@ -288,16 +288,6 @@
             self.logger.info("Saving initial conditions for FEWS...")
             self.wf_suspend(self.Dir + "/outstate/")
 
-        #report(self.CumInwaterMM, self.SaveDir + "/outsum/CumInwaterMM.map")
-        #report(self.CumReinfilt, self.SaveDir + "/outsum/CumReinfilt.map")
-        #report(self.CumPrec, self.SaveDir + "/outsum/CumPrec.map")
-        #report(self.CumEvap, self.SaveDir + "/outsum/CumEvap.map")
-        #report(self.CumPotenTrans, self.SaveDir + "/outsum/CumPotenTrans.map")
-        #report(self.CumInt, self.SaveDir + "/outsum/CumInt.map")
-        #report(self.CumLeakage, self.SaveDir + "/outsum/CumLeakage.map")
-        #report(self.CumPotenEvap, self.SaveDir + "/outsum/CumPotenEvap.map")
-        #report(self.CumExfiltWater, self.SaveDir + "/outsum/CumExfiltWater.map")
-        #report(self.watbal, self.SaveDir + "/outsum/watbal.map")
 
     def parameters(self):
         """