diff --git a/F90subroutines/routing_interface.f90 b/F90subroutines/routing_interface.f90
index 89f67adda14e1604b5878baeda2225bd027b4ab2..1f5b8492bc0bff10f85a116a7c8753f86b6980c2 100644
--- a/F90subroutines/routing_interface.f90
+++ b/F90subroutines/routing_interface.f90
@@ -45,7 +45,7 @@ END SUBROUTINE closeinterface
 
 SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area,  max_basins, min_topoind, &
      & lon_rel, lat_rel, trip, basins, topoindex, fac, hierarchy, nbi, nbj, area_bx, trip_bx, basin_bx, &
-     & topoind_bx, fac_bx, hierarchy_bx, lon_bx, lat_bx, lshead_bx, nwbas)
+     & topoind_bx, fac_bx, hierarchy_bx, lon_bx, lat_bx, lshead_bx)
   !
   USE ioipsl
   USE grid
@@ -83,7 +83,6 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area
   REAL, INTENT(out) :: topoind_bx(nbpt,nbxmax_in,nbxmax_in) !! Topographic index of the residence time for each of the smaller boxes (m)
   INTEGER, INTENT(out) :: trip_bx(nbpt,nbxmax_in,nbxmax_in) !! The trip field for each of the smaller boxes (unitless)
   INTEGER, INTENT(out) :: basin_bx(nbpt,nbxmax_in,nbxmax_in) !!
-  INTEGER, INTENT(out) :: nwbas
   !
   INTEGER :: ii, ib
   REAL :: resolution(nbpt,2)
@@ -99,8 +98,6 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area
      resolution(ii,2) = SQRT(area(ii))
   ENDDO
 
-  nwbas = MAXVAL(sub_pts)
-
   DO ib=1,nbpt
      CALL routing_reg_getgrid(nbpt, ib, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
           & lon_rel, lat_rel, lalo, resolution, contfrac, trip, basins, topoindex, fac, hierarchy, &
diff --git a/Interface.py b/Interface.py
index 9ef45cb585002eb3ba888beedfabb70f706145c4..471ca5dbc17f8a00b23c4f096c16dfe37bce7f63 100644
--- a/Interface.py
+++ b/Interface.py
@@ -259,12 +259,13 @@ class HydroOverlap :
         print("GETHYDROGRID : nbvmax = ", nbvmax)
         print("GETHYDROGRID : nbxmax = ", nbxmax)
         self.nbi, self.nbj, self.area_bx, self.trip_bx, self.basin_bx, self.topoind_bx, self.fac_bx, self.hierarchy_bx, \
-            self.lon_bx, self.lat_bx, self.lshead_bx, self.nwbas = \
+            self.lon_bx, self.lat_bx, self.lshead_bx = \
                     routing_interface.gethydrogrid(nbxmax, sub_pts, sub_index, sub_area, \
                     hydrodata.basinsmax, hydrodata.topoindmin, sub_lon, sub_lat, trip_tmp, basins_tmp, topoind_tmp, fac_tmp, hierarchy_tmp)
         #
         # Plot some diagnostics for the hydrology grid within the atmospheric meshes.
         #
+        self.nwbas = nbvmax
         # Clean-up these arrays so that they are easy to use in Python.
         self.lon_bx[self.lon_bx > 360.]=np.nan
         self.lat_bx[self.lat_bx > 90.]=np.nan