Index: wflow/sphy/advanced_routing.py
===================================================================
diff -u
--- wflow/sphy/advanced_routing.py	(revision 0)
+++ wflow/sphy/advanced_routing.py	(revision 5bc2645dde878c05902af9d9e373036e8d9908b2)
@@ -0,0 +1,11 @@
+#-Function to rout the specific runoff
+def ROUT(self, pcr, rvolume, qold, qout, sres):
+    # Calculate the discharge Q (m3/d)
+    Q = pcr.accufractionflux(self.FlowDir, rvolume, self.QFRAC)
+    # Re-calculate Q, based on qold en kx, and assign Qout for cells being lake/reservoir
+    Q = pcr.ifthenelse(self.QFRAC==0, qout, (1-self.kx) * Q + (qold*24*3600) * self.kx)
+    # Only calculate inflow for lake/reservoir cells
+    Qin = pcr.ifthenelse(self.QFRAC==0, pcr.upstream(self.FlowDir, Q), 0)
+    sres = sres - qout + Qin
+    Q = Q / (24 * 3600)  #-only convert Q to m3/s
+    return sres, Q, Qin
Index: wflow/sphy/dynamic_veg.py
===================================================================
diff -u
--- wflow/sphy/dynamic_veg.py	(revision 0)
+++ wflow/sphy/dynamic_veg.py	(revision 5bc2645dde878c05902af9d9e373036e8d9908b2)
@@ -0,0 +1,18 @@
+#-Function that returns crop factor (Kc) and maximum storage (Smax)
+def Veg_function(pcr, ndvi, fpar_max, fpar_min, lai_max, ndvi_min, ndvi_max, kc_min, kc_max):
+    SR = (1 + ndvi)/(1 - ndvi)
+    SR_max = (1 + ndvi_max)/(1 - ndvi_max)
+    SR_min = (1 + ndvi_min)/(1 - ndvi_min)
+    FPAR = pcr.min((((SR - SR_min) * (fpar_max - fpar_min))/ (SR_max - SR_min)) + 0.001, 0.95)
+    LAI = lai_max * pcr.log10(1-FPAR)/pcr.log10(1-fpar_max)
+    Smax = 0.935 + 0.498*LAI - 0.00575*(LAI**2)            
+    Kc = kc_min + (kc_max - kc_min) * pcr.max(pcr.min((ndvi - ndvi_min)/(ndvi_max - ndvi_min), 1), 0)
+    return Kc, Smax  
+
+#-Function that returns the interception, precipitation throughfall, and remaining storage
+def Inter_function(pcr, S, Smax, Etr):
+    PreT = pcr.max(0, S - Smax)
+    S = S - PreT
+    Int = pcr.min(1.5 * Etr, S)
+    S = S - Int
+    return Int, PreT, S 
Index: wflow/sphy/lakes.py
===================================================================
diff -u
--- wflow/sphy/lakes.py	(revision 0)
+++ wflow/sphy/lakes.py	(revision 5bc2645dde878c05902af9d9e373036e8d9908b2)
@@ -0,0 +1,52 @@
+#-Function that updates the lake storage and lake level given a measured lake level. If no lake
+# level is measured, then the actual storage is not updated with a measured level. The function
+# returns the updated storage and lake level
+def UpdateLakeHStore(self, pcr, pcrm):
+    #-buffer actual storage
+    OldStorage = self.StorRES
+    #-Check if measured lake levels area available
+    try:
+        LakeLevel = pcr.readmap(pcrm.generateNameT(self.LLevel, self.counter))
+        Level = True
+    except:
+        Level = False
+    if Level:
+        #-update lake storage according to measured water level
+        self.StorRES = pcr.ifthenelse(self.UpdateLakeLevel, pcr.ifthenelse(pcr.defined(LakeLevel), pcr.ifthenelse(self.LakeSH_Func==1,\
+            self.LakeSH_exp_a * pcr.exp(self.LakeSH_exp_b * LakeLevel), pcr.ifthenelse(self.LakeSH_Func==2, self.LakeSH_pol_a1 \
+            * LakeLevel + self.LakeSH_pol_b, pcr.ifthenelse(self.LakeSH_Func==3, (self.LakeSH_pol_a2 * LakeLevel**2) + \
+            self.LakeSH_pol_a1 * LakeLevel + self.LakeSH_pol_b, (self.LakeSH_pol_a3 * LakeLevel**3) + (self.LakeSH_pol_a2 \
+            * LakeLevel**2) + (self.LakeSH_pol_a1 * LakeLevel + self.LakeSH_pol_b)))), self.StorRES), self.StorRES)
+        # prevent storage becoming negative for whatever reason
+        self.StorRES = pcr.max(self.StorRES, 0)
+        #-Update the lake level based on the storage for lakes where no levels are measured
+        LakeLevel = pcr.ifthenelse(self.UpdateLakeLevel, pcr.ifthenelse(pcr.defined(LakeLevel), LakeLevel, \
+            pcr.ifthenelse(self.LakeHS_Func==1, self.LakeHS_exp_a * pcr.exp(self.LakeHS_exp_b * self.StorRES), pcr.ifthenelse(self.LakeHS_Func==2, self.LakeHS_pol_a1 * \
+            self.StorRES + self.LakeHS_pol_b, pcr.ifthenelse(self.LakeHS_Func==3, (self.LakeHS_pol_a2 * self.StorRES**2) + \
+            self.LakeHS_pol_a1 * self.StorRES + self.LakeHS_pol_b, (self.LakeHS_pol_a3 * self.StorRES**3) + (self.LakeHS_pol_a2 *\
+            self.StorRES**2) + self.LakeHS_pol_a1 * self.StorRES + self.LakeHS_pol_b)))), pcr.ifthenelse(self.LakeHS_Func==1, \
+            self.LakeHS_exp_a * pcr.exp(self.LakeHS_exp_b * self.StorRES), pcr.ifthenelse(self.LakeHS_Func==2, self.LakeHS_pol_a1 * \
+            self.StorRES + self.LakeHS_pol_b, pcr.ifthenelse(self.LakeHS_Func==3, (self.LakeHS_pol_a2 * self.StorRES**2) + \
+            self.LakeHS_pol_a1 * self.StorRES + self.LakeHS_pol_b, (self.LakeHS_pol_a3 * self.StorRES**3) + (self.LakeHS_pol_a2 *\
+            self.StorRES**2) + self.LakeHS_pol_a1 * self.StorRES + self.LakeHS_pol_b))))
+
+    else:
+        # if no lake level map is available, then calculate the h based on storages
+        LakeLevel = pcr.ifthenelse(self.LakeHS_Func==1, self.LakeHS_exp_a * pcr.exp(self.LakeHS_exp_b * self.StorRES), \
+            pcr.ifthenelse(self.LakeHS_Func==2, self.LakeHS_pol_a1 * self.StorRES + self.LakeHS_pol_b, pcr.ifthenelse(\
+            self.LakeHS_Func==3, (self.LakeHS_pol_a2 * self.StorRES**2) + self.LakeHS_pol_a1 * self.StorRES + self.LakeHS_pol_b,\
+            (self.LakeHS_pol_a3 * self.StorRES**3) + (self.LakeHS_pol_a2 * self.StorRES**2) + self.LakeHS_pol_a1 * self.StorRES +\
+            self.LakeHS_pol_b)))
+    self.StorRES = pcr.ifthenelse(self.LakeID != 0, self.StorRES, OldStorage)
+    return LakeLevel, self.StorRES
+
+#-function that calculates the fraction of lake storage that is available for routing, and the lake outflow
+def QLake(self, pcr, LakeLevel):
+    Qout = pcr.ifthenelse(self.LakeQH_Func==1, self.LakeQH_exp_a * pcr.exp(self.LakeQH_exp_b * LakeLevel), pcr.ifthenelse(\
+        self.LakeQH_Func==2, self.LakeQH_pol_a1 * LakeLevel + self.LakeQH_pol_b, pcr.ifthenelse(self.LakeQH_Func==3, \
+        (self.LakeQH_pol_a2 * LakeLevel**2) + self.LakeQH_pol_a1 * LakeLevel + self.LakeQH_pol_b, (self.LakeQH_pol_a3 * \
+        LakeLevel**3) + (self.LakeQH_pol_a2 * LakeLevel**2) + self.LakeQH_pol_a1 * LakeLevel + self.LakeQH_pol_b)))
+    Qout = pcr.max(0, Qout)
+    Qout = Qout * 3600 * 24  #-convert to m3/d
+    Qout = pcr.cover(Qout, 0) #-for non-lake cells, Qout is zero
+    return Qout
Index: wflow/sphy/reservoirs.py
===================================================================
diff -u
--- wflow/sphy/reservoirs.py	(revision 0)
+++ wflow/sphy/reservoirs.py	(revision 5bc2645dde878c05902af9d9e373036e8d9908b2)
@@ -0,0 +1,29 @@
+#-Advanced reservoir     
+def QAdv(self, pcr):
+    DayNo = self.timecalc.julian(self)[0]
+    #-determine if it is flood or dry season
+    S1 = pcr.ifthenelse(self.ResFlStart < self.ResFlEnd, pcr.ifthenelse(DayNo>=self.ResFlStart, pcr.ifthenelse(DayNo<=self.ResFlEnd, pcr.boolean(1), pcr.boolean(0)), pcr.boolean(0)),\
+                        pcr.ifthenelse(DayNo>=self.ResFlEnd, pcr.ifthenelse(DayNo>=self.ResFlStart, pcr.boolean(1), pcr.boolean(0)), pcr.ifthenelse(DayNo<=self.ResFlEnd, \
+                        pcr.ifthenelse(DayNo<=self.ResFlStart, pcr.boolean(1), pcr.boolean(0)), pcr.boolean(0))))
+
+    S_avail = pcr.max(self.StorRES - self.ResPVOL, 0)
+    Q = pcr.max(pcr.ifthenelse(S1, self.ResMaxFl * S_avail / (self.ResEVOL - self.ResPVOL),\
+        self.ResDemFl * S_avail / (self.ResEVOL - self.ResPVOL)), self.StorRES - self.ResEVOL)
+    return Q
+
+#-Simple reservoir    
+def QSimple(self, pcr):
+    Q = pcr.max(self.ResKr * self.StorRES * (self.StorRES / self.ResSmax)**1.5, self.StorRES - self.ResSmax)
+    return Q
+    
+#-Calculates reservoir outflow and the fraction to release, depending on the type of reservoir (simple or advanced)    
+def QRes(self, pcr):
+    if self.ResSimple and self.ResAdvanced:
+        Qout = pcr.ifthenelse(self.ResFunc==1, QSimple(self, pcr), pcr.ifthenelse(self.ResFunc==2,\
+        QAdv(self, pcr), 0))
+    elif self.ResSimple:
+        Qout = pcr.ifthenelse(self.ResFunc==1, QSimple(self, pcr), 0)
+    else:
+        Qout = pcr.ifthenelse(self.ResFunc==2, QAdv(self, pcr), 0)
+    
+    return Qout
\ No newline at end of file