diff --git a/F90subroutines/routing_interface.f90 b/F90subroutines/routing_interface.f90
index a163c6437b718bd8fa140debd3b6a9f1e5ef3369..3d49719cba5b7ef1485c3f1216fc444c343cd8fd 100644
--- a/F90subroutines/routing_interface.f90
+++ b/F90subroutines/routing_interface.f90
@@ -772,7 +772,7 @@ SUBROUTINE correct_inflows(nbpt, nwbas, inflowmax, outflow_grid,&
   INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(inout) :: inflow_grid
 
   ! LOCAL
-  INTEGER(i_std) :: ig, nbas, ib, og, ob, inf, found
+  INTEGER(i_std) :: ig, nbas, ib, og, ob
 
   WRITE(numout,*) "Checking if the HTUs are in the inflows of their outflow"
 
diff --git a/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index e569b5b45f0176df41e5e53d8a812ccf765ba44f..2e3dada7f5c07876fe8a3a49659b9e2deff192da 100644
--- a/F90subroutines/routing_reg.f90
+++ b/F90subroutines/routing_reg.f90
@@ -2866,7 +2866,6 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b
     ! Update the information needed in the basin "totakeover"
     ! For the moment only area
     !
-    WRITE(numout, *) "XXXXX Arguments : inflow_number = ", inflow_number(ib,1:10)
     !
     basin_cg(ib, totakeover, 1) = (basin_cg(ib, totakeover, 1)*basin_area(ib, totakeover) &
     & + basin_cg(ib, tokill, 1)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))
diff --git a/HydroGrid.py b/HydroGrid.py
index da9dc9c341f855107268b44a38a76a87c4939eb1..0ebb1b0256ece24159b7ee16d570699d760ef996 100644
--- a/HydroGrid.py
+++ b/HydroGrid.py
@@ -93,9 +93,14 @@ class HydroGrid :
 class HydroData :
     def __init__(self, nf, box, index) :
         istr, iend, jstr, jend = box[:]
+        #
+        # Flow direction
+        #
         self.trip=gather(nf.variables["trip"][jstr:jend,istr:iend].astype(np.float32), index, 97)
         self.tripdesc=nf.variables["trip"].long_name
         #
+        # ID of basin
+        #
         self.basins=gather(nf.variables["basins"][jstr:jend,istr:iend], index, 999)
         self.basinsdesc=nf.variables["basins"].long_name
         att = getattrcontaining(nf, "basins", "max")
@@ -106,36 +111,33 @@ class HydroData :
             # This variable seems not to be used further
             self.basinsmax = part.domainmax(np.ma.max(ma.masked_where(self.basins < 1.e10, self.basins)))
         #
-        self.topoind=gather(nf.variables["topoind"][jstr:jend,istr:iend].astype(np.float32), index, 10)
-        self.topoinddesc=nf.variables["topoind"].long_name
-        att = getattrcontaining(nf, "topoind", "min")
-        if len(att) > 0 :
-            self.topoindmin=att[0]
-        else :
-            INFO("We need to scan full file to find minimum topoind over domain")
-            self.topoindmin=np.min(np.where(nf.variables["topoind"][:,:] <  1.e15))
-        #
-        #
-        self.topoindh=gather(nf.variables["topoind_h"][jstr:jend,istr:iend].astype(np.float32), index, 10)
-        self.topoindhdesc=nf.variables["topoind_h"].long_name
-        att = getattrcontaining(nf, "topoind_h", "min")
-        if len(att) > 0 :
-            self.topoindhmin=att[0]
-        else :
-            INFO("We need to scan full file to find minimum topoind_h over domain")
-            self.topoindhmin=np.min(np.where(nf.variables["topoind_h"][:,:] <  1.e15))
+        # Distance to ocean
         #
         self.disto=gather(nf.variables["disto"][jstr:jend,istr:iend].astype(np.float32), index, 0)
         self.distodesc=nf.variables["disto"].long_name
         #
+        # Flow accumulation
+        #
         self.fac=gather(nf.variables["fac"][jstr:jend,istr:iend].astype(np.float32), index, 0)
         self.facdesc=nf.variables["fac"].long_name
         #
+        # If Orography is present in the file then we can compute the topographic index.
+        # Else the topographic index needs to be present in the file.
+        #
         if "orog" in nf.variables.keys():
             self.orog = gather(nf.variables["orog"][jstr:jend,istr:iend].astype(np.float32), index, 0)
             self.orogdesc=nf.variables["orog"].long_name
         else:
             self.orog = gather(np.zeros((jend-jstr,iend-istr)).astype(np.float32), index)
+            self.topoind=gather(nf.variables["topoind"][jstr:jend,istr:iend].astype(np.float32), index, 10)
+            self.topoinddesc=nf.variables["topoind"].long_name
+            att = getattrcontaining(nf, "topoind", "min")
+            if len(att) > 0 :
+                self.topoindmin=att[0]
+            else :
+                INFO("We need to scan full file to find minimum topoind over domain")
+                self.topoindmin=np.min(np.where(nf.variables["topoind"][:,:] <  1.e15))
+        #
         #
         if "floodplains" in nf.variables.keys():
             self.floodplains = gather(nf.variables["floodplains"][jstr:jend,istr:iend].astype(np.float32), index, 0)
diff --git a/RoutingPreProc.py b/RoutingPreProc.py
index b3062b8f539df7390d2d9ee5a20ae44acddd7810..d9bd5976b1a4ab2391f35b9954e675cfb90ebc4a 100644
--- a/RoutingPreProc.py
+++ b/RoutingPreProc.py
@@ -63,7 +63,10 @@ w = RPP.compweights(wfile, part, modelgrid, hydrogrid)
 #
 #
 nbpt = len(w.index)
-nbvmax = part.domainmax(max(w.hpts))
+nbvmax = part.domainmax(w.maxhpts)
+#
+nbasmax = min(nbasmax,nbvmax)
+#
 INFO("nbpt : {0}".format(nbpt))
 INFO("nbvmax : {0}".format(nbvmax))
 INFO("nbasmax : {0}".format(nbasmax))