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/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index 57b0d017df7ef9353c87c95238a17e0ba8758222..0fb14293419f4207a94abc08b25b991c0fe9a37e 100644
--- a/F90subroutines/routing_reg.f90
+++ b/F90subroutines/routing_reg.f90
@@ -2291,13 +2291,20 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
                          dop = minhloc(basin_id(sp,sb),1)
                          bop = minhloc(basin_id(sp,sb),2)
                          IF ((dop  < undef_int ) .AND. (dop .GT. 0) .AND. (bop < undef_int) .AND. (bop .GT. 0)) THEN
-                            !basin_hierarchy(sp,sb) = basin_hierarchy(dop,bop)
+
                             CALL routing_updateflow(sp, sb, dop, bop, nbpt,nwbas, nbxmax*8, outflow_grid, outflow_basin, &
                                  &   inflow_number, inflow_grid, inflow_basin)
-                         found = 1
-                         WRITE(numout,*) "Lowest basin hierarchy in the neighbours grid points"
-                         WRITE(numout,*) "hierarch of (sp,sb)", basin_hierarchy(sp,sb) 
-                         WRITE(numout,*) "Lowest basinid hierarch", basin_hierarchy(dop,bop)
+                            ! It is possible that the lowest hierarchy is in two grid cells
+                            ! In this case, if we have an error for this routing_updateflow
+                            ! We go to next step and make it a coastal flow 
+                            IF (outflow_basin(sp, sb) < undef_int) THEN
+                               found = 1
+                               WRITE(numout,*) "Lowest basin hierarchy in the neighbours grid points"
+                            ELSE
+                               WRITE(numout,*) "Lowest hierarchy may be in two different grid points"
+                            END IF
+                            WRITE(numout,*) "hierarch of (sp,sb)", basin_hierarchy(sp,sb) 
+                            WRITE(numout,*) "Lowest basinid hierarch", basin_hierarchy(dop,bop)
                          END IF
                       END IF
                    END DO
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