Index: wflow/wflow_lib.py
===================================================================
diff -u -rc427e9b948a7fbad40a1c828b68355cb277dc5e0 -ra9498adee6baab0a0abaa331041be8948510167b
--- wflow/wflow_lib.py	(.../wflow_lib.py)	(revision c427e9b948a7fbad40a1c828b68355cb277dc5e0)
+++ wflow/wflow_lib.py	(.../wflow_lib.py)	(revision a9498adee6baab0a0abaa331041be8948510167b)
@@ -133,11 +133,15 @@
     inflow = pcr.ifthen(pcr.boolean(ReserVoirLocs), inflow)
 
     prec_av = pcr.cover(
-        pcr.ifthen(pcr.boolean(ReserVoirLocs), pcr.areaaverage(precip, ReservoirSimpleAreas)),
+        pcr.ifthen(
+            pcr.boolean(ReserVoirLocs), pcr.areaaverage(precip, ReservoirSimpleAreas)
+        ),
         pcr.scalar(0.0),
     )
     pet_av = pcr.cover(
-        pcr.ifthen(pcr.boolean(ReserVoirLocs), pcr.areaaverage(pet, ReservoirSimpleAreas)),
+        pcr.ifthen(
+            pcr.boolean(ReserVoirLocs), pcr.areaaverage(pet, ReservoirSimpleAreas)
+        ),
         pcr.scalar(0.0),
     )
 
@@ -230,8 +234,12 @@
 
     inflow = pcr.ifthen(pcr.boolean(ReserVoirLocs), inflow)
 
-    prec_av = pcr.ifthen(pcr.boolean(ReserVoirLocs), pcr.areaaverage(precip, ReservoirComplexAreas))
-    pet_av = pcr.ifthen(pcr.boolean(ReserVoirLocs), pcr.areaaverage(pet, ReservoirComplexAreas))
+    prec_av = pcr.ifthen(
+        pcr.boolean(ReserVoirLocs), pcr.areaaverage(precip, ReservoirComplexAreas)
+    )
+    pet_av = pcr.ifthen(
+        pcr.boolean(ReserVoirLocs), pcr.areaaverage(pet, ReservoirComplexAreas)
+    )
 
     np_reslocs = pcr.pcr2numpy(ReserVoirLocs, 0.0)
     np_linkedreslocs = pcr.pcr2numpy(LinkedReserVoirLocs, 0.0)
@@ -274,7 +282,7 @@
         np_outflow = pcr.pcr2numpy(outflow, np.nan)
         np_outflow_linked = np_reslocs * 0.0
 
-        with np.errstate(invalid='ignore'):
+        with np.errstate(invalid="ignore"):
             if np_outflow[np_outflow < 0] is not None:
                 np_outflow_linked[
                     np.in1d(np_reslocs, np_linkedreslocs[np_outflow < 0]).reshape(
@@ -301,7 +309,7 @@
         )
 
         np_outflow_nz = np_outflow * 0.0
-        with np.errstate(invalid='ignore'):
+        with np.errstate(invalid="ignore"):
             np_outflow_nz[np_outflow > 0] = np_outflow[np_outflow > 0]
         _outflow.append(np_outflow_nz)
 
@@ -528,7 +536,9 @@
     """
     strorder = pcr.streamorder(ldd)
     strorder = pcr.ifthen(strorder >= pcr.ordinal(order), strorder)
-    dist = pcr.max(pcr.celllength(), pcr.ifthen(pcr.boolean(strorder), pcr.downstreamdist(ldd)))
+    dist = pcr.max(
+        pcr.celllength(), pcr.ifthen(pcr.boolean(strorder), pcr.downstreamdist(ldd))
+    )
 
     return pcr.catchmenttotal(pcr.cover(dist, 0), ldd), dist, strorder
 
@@ -552,11 +562,18 @@
 
     strorder = pcr.streamorder(ldd)
     strorder = pcr.ifthen(strorder >= order, strorder)
-    dist = pcr.cover(pcr.max(pcr.celllength(), pcr.ifthen(pcr.boolean(strorder), pcr.downstreamdist(ldd))), 0)
+    dist = pcr.cover(
+        pcr.max(
+            pcr.celllength(), pcr.ifthen(pcr.boolean(strorder), pcr.downstreamdist(ldd))
+        ),
+        0,
+    )
     totdist = pcr.max(
         pcr.ifthen(
             pcr.boolean(strorder),
-            pcr.windowtotal(pcr.ifthen(pcr.boolean(strorder), dist), pcr.celllength() * factor),
+            pcr.windowtotal(
+                pcr.ifthen(pcr.boolean(strorder), dist), pcr.celllength() * factor
+            ),
         ),
         dist,
     )
@@ -648,7 +665,9 @@
         - outlet map (single point in the map)
     """
     largest = pcr.mapmaximum(pcr.catchmenttotal(pcr.spatial(pcr.scalar(1.0)), ldd))
-    outlet = pcr.ifthen(pcr.catchmenttotal(1.0, ldd) == largest, pcr.spatial(pcr.scalar(1.0)))
+    outlet = pcr.ifthen(
+        pcr.catchmenttotal(1.0, ldd) == largest, pcr.spatial(pcr.scalar(1.0))
+    )
 
     return outlet
 
@@ -793,15 +812,20 @@
     elif assign_existing:
         # unaccounted areas are added to largest nearest draining basin
         if up_area is None:
-            up_area = pcr.ifthen(pcr.boolean(pcr.cover(stream_ge, 0)), pcr.accuflux(ldd, 1))
+            up_area = pcr.ifthen(
+                pcr.boolean(pcr.cover(stream_ge, 0)), pcr.accuflux(ldd, 1)
+            )
         riverid = pcr.ifthen(pcr.boolean(pcr.cover(stream_ge, 0)), subcatch)
 
         friction = 1.0 / pcr.scalar(
             pcr.spreadzone(pcr.cover(pcr.ordinal(up_area), 0), 0, 0)
         )  # *(pcr.scalar(ldd)*0+1)
         delta = pcr.ifthen(
             pcr.scalar(ldd) >= 0,
-            pcr.ifthen(pcr.cover(subcatch, 0) == 0, pcr.spreadzone(pcr.cover(riverid, 0), 0, friction)),
+            pcr.ifthen(
+                pcr.cover(subcatch, 0) == 0,
+                pcr.spreadzone(pcr.cover(riverid, 0), 0, friction),
+            ),
         )
         subcatch = pcr.ifthenelse(pcr.boolean(pcr.cover(subcatch, 0)), subcatch, delta)
 
@@ -836,7 +860,13 @@
     pts = pcr.ifthen((pcr.scalar(sttd) - pcr.scalar(stt)) > 0.0, sttd)
     dif = pcr.upstream(
         ldd,
-        pcr.cover(pcr.ifthen(large, pcr.uniqueid(pcr.boolean(pcr.ifthen(stt == pcr.ordinal(oorder), pts)))), 0),
+        pcr.cover(
+            pcr.ifthen(
+                large,
+                pcr.uniqueid(pcr.boolean(pcr.ifthen(stt == pcr.ordinal(oorder), pts))),
+            ),
+            0,
+        ),
     )
     dif = pcr.cover(pcr.scalar(outl), dif)  # Add catchment outlet
     dif = pcr.ordinal(pcr.uniqueid(pcr.boolean(dif)))
@@ -880,20 +910,31 @@
     if fill:
         for order in range(oorder, maxorder):
             m_pts = pcr.ifthen((pcr.scalar(sttd) - pcr.scalar(order)) > 0.0, sttd)
-            m_dif = pcr.uniqueid(pcr.boolean(pcr.ifthen(stt == pcr.ordinal(order), m_pts)))
+            m_dif = pcr.uniqueid(
+                pcr.boolean(pcr.ifthen(stt == pcr.ordinal(order), m_pts))
+            )
             dif = pcr.uniqueid(pcr.boolean(pcr.cover(m_dif, dif)))
 
         for myorder in range(oorder - 1, stoporder, -1):
             sc = pcr.subcatchment(ldd, pcr.nominal(dif))
             m_pts = pcr.ifthen((pcr.scalar(sttd) - pcr.scalar(stt)) > 0.0, sttd)
-            m_dif = pcr.uniqueid(pcr.boolean(pcr.ifthen(stt == pcr.ordinal(myorder - 1), m_pts)))
-            dif = pcr.uniqueid(pcr.boolean(pcr.cover(pcr.ifthen(pcr.scalar(sc) == 0, m_dif), dif)))
+            m_dif = pcr.uniqueid(
+                pcr.boolean(pcr.ifthen(stt == pcr.ordinal(myorder - 1), m_pts))
+            )
+            dif = pcr.uniqueid(
+                pcr.boolean(pcr.cover(pcr.ifthen(pcr.scalar(sc) == 0, m_dif), dif))
+            )
 
         if fillcomplete:
             sc = pcr.subcatchment(ldd, pcr.nominal(dif))
             cs, m_dif, stt = subcatch_order_a(ldd, stoporder)
             dif = pcr.uniqueid(
-                pcr.boolean(pcr.cover(pcr.ifthen(pcr.scalar(sc) == 0, pcr.ordinal(m_dif)), pcr.ordinal(dif)))
+                pcr.boolean(
+                    pcr.cover(
+                        pcr.ifthen(pcr.scalar(sc) == 0, pcr.ordinal(m_dif)),
+                        pcr.ordinal(dif),
+                    )
+                )
             )
 
     scsize = pcr.catchmenttotal(1, ldd)
@@ -997,7 +1038,7 @@
     # Loop over points and "burn in" map
     for n in range(0, xcor.size):
         if Verbose:
-            print (n)
+            print(n)
         diffx = x - xcor[n]
         diffy = y - ycor[n]
         col_ = np.absolute(diffx) <= (cell_length * tolerance)  # cellsize
@@ -1030,7 +1071,9 @@
         pcr.ifthenelse(
             draindir == 8,
             yl,
-            pcr.ifthenelse(draindir == 4, xl, pcr.ifthenelse(draindir == 6, xl, slantlength)),
+            pcr.ifthenelse(
+                draindir == 4, xl, pcr.ifthenelse(draindir == 6, xl, slantlength)
+            ),
         ),
     )
 
@@ -1060,7 +1103,9 @@
         pcr.ifthenelse(
             draindir == 8,
             xl,
-            pcr.ifthenelse(draindir == 4, yl, pcr.ifthenelse(draindir == 6, yl, slantwidth)),
+            pcr.ifthenelse(
+                draindir == 4, yl, pcr.ifthenelse(draindir == 6, yl, slantwidth)
+            ),
         ),
     )
     return drainwidth
@@ -1076,7 +1121,9 @@
 
     result = pcr.ordinal(pcr.cover(-1))
     for l, u, c in zip(lower, upper, classes):
-        result = pcr.cover(pcr.ifthen(inmap >= l, pcr.ifthen(inmap < u, pcr.ordinal(c))), result)
+        result = pcr.cover(
+            pcr.ifthen(inmap >= l, pcr.ifthen(inmap < u, pcr.ordinal(c))), result
+        )
 
     return pcr.ifthen(result >= 0, result)
 
@@ -1103,7 +1150,9 @@
             according to D8 directions
     """
     if rivers is None:
-        stream = pcr.ifthenelse(pcr.accuflux(ldd, 1) >= accuThreshold, pcr.boolean(1), pcr.boolean(0))
+        stream = pcr.ifthenelse(
+            pcr.accuflux(ldd, 1) >= accuThreshold, pcr.boolean(1), pcr.boolean(0)
+        )
     else:
         stream = pcr.boolean(pcr.cover(rivers, 0))
 
@@ -1112,7 +1161,9 @@
         up_elevation = pcr.scalar(pcr.subcatchment(ldd, height_river))
     else:
         drainage_surf = pcr.ifthen(rivers, pcr.accuflux(ldd, 1))
-        weight = 1.0 / pcr.scalar(pcr.spreadzone(pcr.cover(pcr.ordinal(drainage_surf), 0), 0, 0))
+        weight = 1.0 / pcr.scalar(
+            pcr.spreadzone(pcr.cover(pcr.ordinal(drainage_surf), 0), 0, 0)
+        )
         up_elevation = pcr.ifthenelse(
             basin,
             pcr.scalar(pcr.subcatchment(ldd, height_river)),