diff --git a/Interface.py b/Interface.py
index be42d37af04d9f8035b85a81b100c1b9bf3e0a84..14ce51df7e0792dca92f9441ec18f92778e2554b 100644
--- a/Interface.py
+++ b/Interface.py
@@ -296,14 +296,17 @@ class HydroSuper :
         #
         # nb_htu can be adjusted with self.nwbas
         # nb_htu can be lowered with a larger maxpercent (routing_reg.f90)
-        nb_htu = 600 
+        nb_htu = 600
+        nb_htu = nbvmax
         nbv = nbvmax
         #
         # Call findbasins
         #
         nb_basin, basin_inbxid, basin_outlet, basin_outtp, self.basin_sz, basin_bxout, basin_bbout, self.basin_pts, basin_lshead, coast_pts = \
-                    routing_interface.findbasins(nbpt = self.nbpt, nb_htu = nb_htu, nbv = nbv, nbi = hydrooverlap.nbi, nbj = hydrooverlap.nbj, trip_bx = hydrooverlap.trip_bx, \
-                                                 basin_bx = hydrooverlap.basin_bx, fac_bx = hydrooverlap.fac_bx, hierarchy_bx = hydrooverlap.hierarchy_bx, \
+                    routing_interface.findbasins(nbpt = self.nbpt, nb_htu = nb_htu, nbv = nbv, nbi = hydrooverlap.nbi, nbj = hydrooverlap.nbj, \
+                                                 trip_bx = hydrooverlap.trip_bx, \
+                                                 basin_bx = hydrooverlap.basin_bx, fac_bx = hydrooverlap.fac_bx, \
+                                                 hierarchy_bx = hydrooverlap.hierarchy_bx, \
                                                  topoind_bx = hydrooverlap.topoind_bx, lshead_bx = hydrooverlap.lshead_bx, \
                                                  lontmp = hydrooverlap.lon_bx, lattmp = hydrooverlap.lat_bx)
         #
@@ -322,10 +325,12 @@ class HydroSuper :
             self.basin_orog, self.basin_floodp, self.basin_fac, self.basin_topoind, \
             self.basin_id, self.basin_outcoor, self.basin_type, self.basin_flowdir, \
             self.basin_lshead, self.outflow_grid, self.outflow_basin, self.nbcoastal, self.coastal_basin = \
-                    routing_interface.globalize(nbpt = self.nbpt, nb_htu = nb_htu,nbv = nbv, area_bx = hydrooverlap.area_bx, lon_bx = lon_bx_tmp, lat_bx = lat_bx_tmp, trip_bx = hydrooverlap.trip_bx, \
+                    routing_interface.globalize(nbpt = self.nbpt, nb_htu = nb_htu,nbv = nbv, area_bx = hydrooverlap.area_bx, lon_bx = lon_bx_tmp, \
+                                                lat_bx = lat_bx_tmp, trip_bx = hydrooverlap.trip_bx, \
                                                 hierarchy_bx = hydrooverlap.hierarchy_bx, orog_bx = hydrooverlap.orog_bx, floodp_bx =  hydrooverlap.floodp_bx,\
                                                 fac_bx = hydrooverlap.fac_bx, topoind_bx = hydrooverlap.topoind_bx, min_topoind = hydrodata.topoindmin, \
-                                                nb_basin = nb_basin, basin_inbxid = basin_inbxid, basin_outlet = basin_outlet, basin_outtp = basin_outtp, basin_sz = self.basin_sz, basin_pts = self.basin_pts, basin_bxout = basin_bxout, \
+                                                nb_basin = nb_basin, basin_inbxid = basin_inbxid, basin_outlet = basin_outlet, basin_outtp = basin_outtp, \
+                                                basin_sz = self.basin_sz, basin_pts = self.basin_pts, basin_bxout = basin_bxout, \
                                                 basin_bbout = basin_bbout, lshead = basin_lshead, coast_pts = coast_pts, nwbas = self.nwbas)
 
         # Memory management
@@ -395,13 +400,15 @@ class HydroSuper :
 
     def check_fetch(self):
 
-        routing_interface.checkfetch(nbpt = self.nbpt, nwbas = self.nwbas, fetch_basin = self.fetch_basin, outflow_grid = self.outflow_grid, outflow_basin = self.outflow_basin, basin_count = self.basin_count)
+        routing_interface.checkfetch(nbpt = self.nbpt, nwbas = self.nwbas, fetch_basin = self.fetch_basin, outflow_grid = self.outflow_grid, \
+                                     outflow_basin = self.outflow_basin, basin_count = self.basin_count)
 
         return
 
     def check_routing(self):
 
-        routing_interface.checkrouting(nbpt = self.nbpt, nwbas = self.nwbas, outflow_grid = self.outflow_grid, outflow_basin = self.outflow_basin, basin_count = self.basin_count)
+        routing_interface.checkrouting(nbpt = self.nbpt, nwbas = self.nwbas, outflow_grid = self.outflow_grid, outflow_basin = self.outflow_basin, \
+                                       basin_count = self.basin_count)
 
         return
     #
@@ -549,13 +556,19 @@ class HydroGraph :
         self.routing_area, self.routing_orog, self.routing_floodp, self.routing_cg, self.topo_resid, self.route_nbbasin,\
             self.route_togrid, self.route_tobasin, self.route_nbintobas, self.global_basinid, \
             self.route_outlet, self.route_type, self.origin_nbintobas, self.routing_fetch = \
-                                    routing_interface.finish_truncate(nbpt = self.nbpt, nbxmax_in = nbxmax_in, nbasmax = nbasmax, nwbas = nwbas, num_largest = hydrosuper.num_largest, gridarea = modelgrid.area, cfrac = modelgrid.contfrac, basin_count = hydrosuper.basin_count, \
-                                                               basin_notrun = hydrosuper.basin_notrun, basin_area = hydrosuper.basin_area, \
-                                                               basin_orog = hydrosuper.basin_orog, basin_floodp = hydrosuper.basin_floodp, basin_cg = hydrosuper.basin_cg, \
-                                                               basin_topoind = hydrosuper.basin_topoind, fetch_basin = hydrosuper.fetch_basin, basin_id = hydrosuper.basin_id, \
-                                                               basin_coor = hydrosuper.basin_outcoor, basin_type = hydrosuper.basin_type, basin_flowdir = hydrosuper.basin_flowdir, \
-                                                               outflow_grid = hydrosuper.outflow_grid, outflow_basin = hydrosuper.outflow_basin, \
-                                                               inflow_number = hydrosuper.inflow_number, inflow_grid = hydrosuper.inflow_grid, inflow_basin = hydrosuper.inflow_basin)
+                                    routing_interface.finish_truncate(nbpt = self.nbpt, nbxmax_in = nbxmax_in, nbasmax = nbasmax, nwbas = nwbas, \
+                                                                      num_largest = hydrosuper.num_largest, gridarea = modelgrid.area, \
+                                                                      cfrac = modelgrid.contfrac, basin_count = hydrosuper.basin_count, \
+                                                                      basin_notrun = hydrosuper.basin_notrun, basin_area = hydrosuper.basin_area, \
+                                                                      basin_orog = hydrosuper.basin_orog, basin_floodp = hydrosuper.basin_floodp, \
+                                                                      basin_cg = hydrosuper.basin_cg, \
+                                                                      basin_topoind = hydrosuper.basin_topoind, fetch_basin = hydrosuper.fetch_basin, \
+                                                                      basin_id = hydrosuper.basin_id, \
+                                                                      basin_coor = hydrosuper.basin_outcoor, basin_type = hydrosuper.basin_type, \
+                                                                      basin_flowdir = hydrosuper.basin_flowdir, \
+                                                                      outflow_grid = hydrosuper.outflow_grid, outflow_basin = hydrosuper.outflow_basin, \
+                                                                      inflow_number = hydrosuper.inflow_number, inflow_grid = hydrosuper.inflow_grid, \
+                                                                      inflow_basin = hydrosuper.inflow_basin)
 
         #
         self.routing_fetch = finalfetch(part, self.routing_area, self.route_nbbasin, self.route_togrid, self.route_tobasin, self.routing_fetch)
@@ -566,7 +579,12 @@ class HydroGraph :
         # Inflows
         self.max_inflow = part.domainmax(np.max(hydrosuper.inflow_number))
         gingrid = part.l2glandindex( hydrosuper.inflow_grid[:,:,:self.max_inflow])
-        self.route_innum, self.route_ingrid, self.route_inbasin = routing_interface.finish_inflows(nbpt = self.nbpt, nbxmax_in = nbxmax_in, nbasmax = nbasmax, inf_max = self.max_inflow, basin_count = hydrosuper.basin_count, inflow_number = hydrosuper.inflow_number, inflow_grid = gingrid, inflow_basin = hydrosuper.inflow_basin[:,:,:self.max_inflow])
+        self.route_innum, self.route_ingrid, self.route_inbasin = routing_interface.finish_inflows(nbpt = self.nbpt, nbxmax_in = nbxmax_in, \
+                                                                                                   nbasmax = nbasmax, inf_max = self.max_inflow, \
+                                                                                                   basin_count = hydrosuper.basin_count, \
+                                                                                                   inflow_number = hydrosuper.inflow_number, \
+                                                                                                   inflow_grid = gingrid, \
+                                                                                                   inflow_basin = hydrosuper.inflow_basin[:,:,:self.max_inflow])
 
         return
     #