Index: wflow-py/wflow/wflow_routing.py
===================================================================
diff -u -r241799f88889c8c189027b512b229bbd4430aa34 -rca8480f9ef36c09ba55cb0b5a1fa263dec2b14c1
--- wflow-py/wflow/wflow_routing.py	(.../wflow_routing.py)	(revision 241799f88889c8c189027b512b229bbd4430aa34)
+++ wflow-py/wflow/wflow_routing.py	(.../wflow_routing.py)	(revision ca8480f9ef36c09ba55cb0b5a1fa263dec2b14c1)
@@ -619,10 +619,17 @@
             self.OutflowDwn = upstream(self.TopoLddOrg,cover(self.Outflow,scalar(0.0)))
             self.Inflow = cover(self.OutflowDwn,self.Inflow)
 
+        # Check if we do not try to abstract more runoff then present
+        self.InflowKinWaveCell = upstream(self.TopoLdd, self.SurfaceRunoff)
+        # The extraction should be equal to the discharge upstream cell.
+        # You should not make the abstraction depended on the downstream cell, because they are correlated.
+        # During a stationary sum they will get equal to each other.
+        MaxExtract = self.InflowKinWaveCell + self.Inwater
+        self.SurfaceWaterSupply = ifthenelse (self.Inflow < 0.0 , min(MaxExtract,-1.0 * self.Inflow), self.ZeroMap)
+        self.OldSurfaceRunoff=self.SurfaceRunoff
+        self.OldInwater=self.Inwater
+        self.Inwater = self.Inwater + ifthenelse(self.SurfaceWaterSupply> 0, -1.0 * self.SurfaceWaterSupply,self.Inflow)
 
-        self.Inwater = self.Inwater + self.Inflow  # Add abstractions/inflows in m^3/sec
-
-
         ##########################################################################
         # Runoff calculation via Kinematic wave ##################################
         ##########################################################################
@@ -634,8 +641,41 @@
         self.SurfaceRunoffMM = self.SurfaceRunoff * self.QMMConv  # SurfaceRunoffMM (mm) from SurfaceRunoff (m3/s)
         self.updateRunOff()
         self.InflowKinWaveCell = upstream(self.TopoLdd, self.SurfaceRunoff)
-        self.MassBalKinWave = (-self.KinWaveVolume + self.OldKinWaveVolume) / self.timestepsecs + self.InflowKinWaveCell + self.Inwater - self.SurfaceRunoff
 
+        maxit = 10  # max number of iteration in abstraction calculations
+        nrit = 0
+        if float(mapminimum(spatial(self.Inflow))) < 0.0:
+            while True:
+                nrit += 1
+                oldsup = self.SurfaceWaterSupply
+                ##########################################################################
+                # Iterate to make a better estimation for the supply #####################
+                # (Runoff calculation via Kinematic wave) ################################
+                ##########################################################################
+                MaxExtract = self.InflowKinWaveCell + self.OldInwater
+                self.SurfaceWaterSupply = ifthenelse(self.Inflow < 0.0, min(MaxExtract, -1.0 * self.Inflow),
+                                                     self.ZeroMap)
+                self.Inwater = self.OldInwater + ifthenelse(self.SurfaceWaterSupply > 0, -1.0 * self.SurfaceWaterSupply,
+                                                            self.Inflow)
+                # per distance along stream
+                q = self.Inwater / self.DCL
+                # discharge (m3/s)
+                self.SurfaceRunoff = kinematic(self.TopoLdd, self.OldSurfaceRunoff, q, self.Alpha, self.Beta,
+                                               self.Tslice,
+                                               self.timestepsecs, self.DCL)  # m3/s
+                self.SurfaceRunoffMM = self.SurfaceRunoff * self.QMMConv  # SurfaceRunoffMM (mm) from SurfaceRunoff (m3/s)
+                self.updateRunOff()
+                self.InflowKinWaveCell = upstream(self.TopoLdd, self.OldSurfaceRunoff)
+
+                deltasup = float(mapmaximum(abs(oldsup - self.SurfaceWaterSupply)))
+
+                if deltasup < 0.01 or nrit >= maxit:
+                    break
+
+        self.MassBalKinWave = (-self.KinWaveVolume + self.OldKinWaveVolume) / self.timestepsecs + \
+                                self.InflowKinWaveCell + self.Inwater - self.SurfaceRunoff
+
+
         Runoff = self.SurfaceRunoff
 
         # Updating