Index: wflow-py/wflow/reservoir_Sf.py
===================================================================
diff -u -r6f7894006c07223c46d4d7022d97e446214717bc -rf95a643bab9f58880912cb18318fdb74dd6e10a9
--- wflow-py/wflow/reservoir_Sf.py	(.../reservoir_Sf.py)	(revision 6f7894006c07223c46d4d7022d97e446214717bc)
+++ wflow-py/wflow/reservoir_Sf.py	(.../reservoir_Sf.py)	(revision f95a643bab9f58880912cb18318fdb74dd6e10a9)
@@ -95,6 +95,7 @@
     - Qf is devided over the reservoir numbers for the timesteps matching with the average
     travel time for a calcultation cell
     """
+        
     if nansum(pcr2numpy(self.Transit,NaN)) > 0:
         self.Qflag = self.trackQ[0]  # first bucket is transferred to outlet
         self.trackQ.append(0 * scalar(self.TopoId))  # add new bucket for present time step
@@ -111,3 +112,26 @@
     self.wbSfrout = self.Qftotal - self.Qflag - sum(self.trackQ) + sum (self.trackQ_t)
     
     self.Qflag_ = self.Qflag
+    
+    self.Qtlag = self.Qflag_ / self.timestepsecs + self.Qs_ / 1000 * self.surfaceArea / self.timestepsecs
+    self.QLagTot = areatotal(self.Qtlag, nominal(self.TopoId))  # catchment total runoff with looptijd
+   
+def routingQf_Qs_grid(self):
+    """
+    - Routing of both Qf and Qs
+    - based on a velocity map
+    """
+    self.Qtot = self.Qftotal + self.Qs_  # total local discharge in mm/hour
+    self.Qtotal = self.Qtot / 1000 * self.surfaceArea / self.timestepsecs  # total local discharge in m3/s
+    self.Qstate_t = self.Qstate
+    self.Qrout = accutraveltimeflux(self.TopoLdd, self.Qstate + self.Qtotal, self.velocity)
+    self.Qstate = accutraveltimestate(self.TopoLdd, self.Qstate + self.Qtotal, self.velocity)
+    
+    self.Qtlag = self.Qrout
+    self.QLagTot = self.Qrout
+    
+    # water balance of flux routing
+    self.dSdt = self.Qstate-self.Qstate_t
+    self.WB_rout = (accuflux(self.TopoLdd, self.Qtotal - self.dSdt)-self.Qrout)/accuflux(self.TopoLdd, self.Qtotal)
+    
+    
\ No newline at end of file