diff --git a/F90subroutines/routing_interface.f90 b/F90subroutines/routing_interface.f90
index d2ac0285d8dc97f225f4d399693ad16fc9a0ba82..bba1e001e0af1eceea9b9a20be7e9f184070a643 100644
--- a/F90subroutines/routing_interface.f90
+++ b/F90subroutines/routing_interface.f90
@@ -113,7 +113,7 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area
 END SUBROUTINE gethydrogrid
 
 
-SUBROUTINE findbasins(nbpt, nbvmax_in, nbxmax_in, nbi, nbj, trip_bx, basin_bx, fac_bx, hierarchy_bx, topoind_bx, lshead_bx, &
+SUBROUTINE findbasins(nbpt, nb_htu, nbv, nbxmax_in, nbi, nbj, trip_bx, basin_bx, fac_bx, hierarchy_bx, topoind_bx, lshead_bx, &
      &                nb_basin, basin_inbxid, basin_outlet, basin_outtp, basin_sz, basin_bxout, basin_bbout, basin_pts, &
      &                basin_lshead, coast_pts, lontmp, lattmp)
   !
@@ -124,7 +124,7 @@ SUBROUTINE findbasins(nbpt, nbvmax_in, nbxmax_in, nbi, nbj, trip_bx, basin_bx, f
   !
   ! Arguments
   !
-  INTEGER, INTENT(in)    :: nbvmax_in, nbxmax_in
+  INTEGER, INTENT(in)    :: nb_htu, nbv, nbxmax_in
   INTEGER, INTENT(in)    :: nbpt !! Domain size (unitless)
   INTEGER, INTENT(in)    :: nbi(nbpt), nbj(nbpt) !! Number of point in x and y within the grid (unitless)
   INTEGER, INTENT(inout) :: trip_bx(nbpt,nbxmax_in,nbxmax_in) !! The trip field for each of the smaller boxes (unitless)
@@ -137,15 +137,15 @@ SUBROUTINE findbasins(nbpt, nbvmax_in, nbxmax_in, nbi, nbj, trip_bx, basin_bx, f
   !
   ! OUTPUT VARIABLES
   INTEGER, INTENT(out) :: nb_basin(nbpt)               !! Number of sub-basins (unitless)
-  INTEGER, INTENT(out) :: basin_inbxid(nbpt,nbvmax_in)   !!
-  REAL, INTENT(out)    :: basin_outlet(nbpt,nbvmax_in,2) !!
-  REAL, INTENT(out)    :: basin_outtp(nbpt,nbvmax_in)    !!
-  INTEGER, INTENT(out) :: basin_sz(nbpt,nbvmax_in)       !!
-  INTEGER, INTENT(out) :: basin_bxout(nbpt,nbvmax_in)    !!
-  INTEGER, INTENT(out) :: basin_bbout(nbpt,nbvmax_in)    !!
-  INTEGER, INTENT(out) :: basin_pts(nbpt,nbvmax_in, nbvmax_in, 2) !!
-  REAL, INTENT(out)    :: basin_lshead(nbpt,nbvmax_in)   !!
-  INTEGER, INTENT(out) :: coast_pts(nbpt,nbvmax_in)      !! The coastal flow points (unitless)
+  INTEGER, INTENT(out) :: basin_inbxid(nbpt,nb_htu)   !!
+  REAL, INTENT(out)    :: basin_outlet(nbpt,nb_htu,2) !!
+  REAL, INTENT(out)    :: basin_outtp(nbpt,nb_htu)    !!
+  INTEGER, INTENT(out) :: basin_sz(nbpt,nb_htu)       !!
+  INTEGER, INTENT(out) :: basin_bxout(nbpt,nb_htu)    !!
+  INTEGER, INTENT(out) :: basin_bbout(nbpt,nb_htu)    !!
+  INTEGER, INTENT(out) :: basin_pts(nbpt,nb_htu, nbv, 2) !!
+  REAL, INTENT(out)    :: basin_lshead(nbpt,nb_htu)   !!
+  INTEGER, INTENT(out) :: coast_pts(nbpt,nb_htu)      !! The coastal flow points (unitless)
   !
   ! For debug and get coordinate of river outlet
   !
@@ -156,15 +156,16 @@ SUBROUTINE findbasins(nbpt, nbvmax_in, nbxmax_in, nbi, nbj, trip_bx, basin_bx, f
   !
   INTEGER :: ib
   !
-  diaglalo(1,:) = (/ 39.6791, 2.6687 /)
+  !diaglalo(1,:) = (/ 39.6791, 2.6687 /)
   !
-  IF ( nbvmax_in .NE. nbvmax .OR. nbxmax_in .NE. nbxmax ) THEN
-     WRITE(*,*) "nbvmax or nbxmax have changed !!"
+  IF ( nbxmax_in .NE. nbxmax ) THEN
+     WRITE(*,*) "nbxmax has changed !!"
      STOP
   ENDIF
 
   DO ib=1,nbpt
-     CALL routing_reg_findbasins(ib, nbi(ib), nbj(ib), trip_bx(ib,:,:), basin_bx(ib,:,:), fac_bx(ib,:,:), hierarchy_bx(ib,:,:), &
+     CALL routing_reg_findbasins(nb_htu, nbv, ib, nbi(ib), nbj(ib), trip_bx(ib,:,:), &
+          & basin_bx(ib,:,:), fac_bx(ib,:,:), hierarchy_bx(ib,:,:), &
           & topoind_bx(ib,:,:), lshead_bx(ib,:,:), diaglalo, &
           & nb_basin(ib), basin_inbxid(ib,:), basin_outlet(ib,:,:), basin_outtp(ib,:), basin_sz(ib,:), basin_bxout(ib,:), &
           & basin_bbout(ib,:), basin_pts(ib,:,:,:), basin_lshead(ib,:), coast_pts(ib,:), lontmp(ib,:,:), lattmp(ib,:,:))
@@ -172,7 +173,7 @@ SUBROUTINE findbasins(nbpt, nbvmax_in, nbxmax_in, nbi, nbj, trip_bx, basin_bx, f
 
 END SUBROUTINE findbasins
 
-SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_bx, orog_bx, floodp_bx, &
+SUBROUTINE globalize(nbpt, nb_htu, nbv, nbxmax_in, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_bx, orog_bx, floodp_bx, &
           & fac_bx, topoind_bx, &
           & min_topoind, nb_basin, basin_inbxid, basin_outlet, basin_outtp, basin_sz, basin_pts, basin_bxout, &
           & basin_bbout, lshead, coast_pts, nwbas, basin_count, basin_notrun, basin_area, basin_cg, basin_hierarchy, &
@@ -189,7 +190,7 @@ SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, lon_bx, lat_bx, trip_b
   !
   !! INPUT VARIABLES
   INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
-  INTEGER(i_std), INTENT(in)  :: nbvmax_in, nbxmax_in
+  INTEGER(i_std), INTENT(in)  :: nbv,nb_htu, nbxmax_in
   REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in)    :: area_bx      !! Area of each small box in the grid box (m^2)
   REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in)    :: lon_bx       !! Longitude of each small box in the grid box
   REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in)    :: lat_bx       !! Latitude of each small box in the grid box
@@ -201,14 +202,14 @@ SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, lon_bx, lat_bx, trip_b
   REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in)    :: topoind_bx   !! Topographic index of the residence time for each of the smaller boxes (m)
   REAL(r_std), INTENT (in)                                         :: min_topoind  !! The current minimum of topographic index (m)
   INTEGER(i_std), INTENT (in), DIMENSION(nbpt)                     :: nb_basin     !! Number of sub-basins (unitless)
-  INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nbvmax_in)     :: basin_inbxid, basin_sz !! ID of basin, number of points in the basin
-  REAL(r_std), INTENT (in), DIMENSION(nbpt,nbvmax_in,2)      :: basin_outlet !! Outlet coordinate of each subgrid basin
-  REAL(r_std), INTENT (in), DIMENSION(nbpt,nbvmax_in)        :: basin_outtp  !!
-  INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nbvmax_in,nbvmax_in,2) :: basin_pts  !! Points in each basin
-  INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nbvmax_in)     :: basin_bxout  !! outflow direction
-  INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nbvmax_in)     :: basin_bbout  !! outflow sub-basin
-  REAL(r_std), INTENT (in), DIMENSION(nbpt,nbvmax_in)        :: lshead       !! Large scale heading of outflow.
-  INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nbvmax_in)     :: coast_pts !! The coastal flow points (unitless)
+  INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nb_htu)     :: basin_inbxid, basin_sz !! ID of basin, number of points in the basin
+  REAL(r_std), INTENT (in), DIMENSION(nbpt,nb_htu,2)      :: basin_outlet !! Outlet coordinate of each subgrid basin
+  REAL(r_std), INTENT (in), DIMENSION(nbpt,nb_htu)        :: basin_outtp  !!
+  INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nb_htu,nbv,2) :: basin_pts  !! Points in each basin
+  INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nb_htu)     :: basin_bxout  !! outflow direction
+  INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nb_htu)     :: basin_bbout  !! outflow sub-basin
+  REAL(r_std), INTENT (in), DIMENSION(nbpt,nb_htu)        :: lshead       !! Large scale heading of outflow.
+  INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nb_htu)     :: coast_pts !! The coastal flow points (unitless)
   !! Output VARIABLES
   INTEGER(i_std), INTENT (in) :: nwbas !!
   INTEGER(i_std), INTENT (out), DIMENSION(nbpt)       :: basin_count         !!
@@ -232,13 +233,14 @@ SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, lon_bx, lat_bx, trip_b
   !!
   INTEGER(i_std) :: ib
   !!
-  IF ( nbvmax_in .NE. nbvmax .OR. nbxmax_in .NE. nbxmax ) THEN
-     WRITE(*,*) "GLOBALIZE : nbvmax or nbxmax have changed !!"
+  IF ( nbxmax_in .NE. nbxmax ) THEN
+     WRITE(*,*) "GLOBALIZE : nbxmax has changed !!"
      STOP
   ENDIF
   !!
   DO ib=1,nbpt
-     CALL routing_reg_globalize(nbpt, ib, neighbours, area_bx(ib,:,:), lon_bx(ib,:,:), lat_bx(ib,:,:), trip_bx(ib,:,:), &
+     CALL routing_reg_globalize(nbpt,nb_htu, nbv, ib, neighbours, area_bx(ib,:,:),& 
+          & lon_bx(ib,:,:), lat_bx(ib,:,:), trip_bx(ib,:,:), &
           & hierarchy_bx(ib,:,:), orog_bx(ib,:,:), floodp_bx(ib,:,:), fac_bx(ib,:,:), &
           & topoind_bx(ib,:,:), min_topoind, nb_basin(ib), basin_inbxid(ib,:), basin_outlet(ib,:,:), basin_outtp(ib,:), &
           & basin_sz(ib,:), basin_pts(ib,:,:,:), basin_bxout(ib,:), basin_bbout(ib,:), lshead(ib,:), coast_pts(ib,:), nwbas, &
diff --git a/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index d005d2104527b7bf6e7ee6e5193a82a6017a6f0c..4c061942a47f2772df237228d43dd48f67fb4274 100644
--- a/F90subroutines/routing_reg.f90
+++ b/F90subroutines/routing_reg.f90
@@ -363,13 +363,14 @@ CONTAINS
 !! \n
 !_ ================================================================================================================================
 
-  SUBROUTINE routing_reg_findbasins(ib, nbi, nbj, trip, basin, fac, hierarchy, topoind, lshead, diaglalo, &
+  SUBROUTINE routing_reg_findbasins(nb_htu, nbv, ib, nbi, nbj, trip, basin, fac, hierarchy, topoind, lshead, diaglalo, &
        & nb_basin, basin_inbxid, basin_outlet,  basin_outtp, basin_sz, basin_bxout, basin_bbout, basin_pts, &
        & basin_lshead, coast_pts, lontmp, lattmp)
     !
     IMPLICIT NONE
     !
 !! INPUT VARIABLES
+    INTEGER(i_std), INTENT(in) :: nb_htu, nbv
     INTEGER(i_std), INTENT(in) :: ib !!
     INTEGER(i_std), INTENT(in) :: nbi !! Number of point in x within the grid (unitless)
     INTEGER(i_std), INTENT(in) :: nbj !! Number of point in y within the grid (unitless)
@@ -391,15 +392,15 @@ CONTAINS
     !
 !! OUTPUT VARIABLES
     INTEGER(i_std), INTENT(out) :: nb_basin               !! Number of sub-basins (unitless)
-    INTEGER(i_std), INTENT(out) :: basin_inbxid(nbvmax)   !!
-    REAL(r_std), INTENT(out)    :: basin_outlet(nbvmax,2) !!
-    REAL(r_std), INTENT(out)    :: basin_outtp(nbvmax)    !!
-    INTEGER(i_std), INTENT(out) :: basin_sz(nbvmax)       !!
-    INTEGER(i_std), INTENT(out) :: basin_bxout(nbvmax)    !!
-    REAL(r_std), INTENT(out)    :: basin_lshead(nbvmax)   !!
-    INTEGER(i_std), INTENT(out) :: basin_bbout(nbvmax)    !!
-    INTEGER(i_std), INTENT(out) :: basin_pts(nbvmax, nbvmax, 2) !!
-    INTEGER(i_std), INTENT(out) :: coast_pts(nbvmax)      !! The coastal flow points (unitless)
+    INTEGER(i_std), INTENT(out) :: basin_inbxid(nb_htu)   !!
+    REAL(r_std), INTENT(out)    :: basin_outlet(nb_htu,2) !!
+    REAL(r_std), INTENT(out)    :: basin_outtp(nb_htu)    !!
+    INTEGER(i_std), INTENT(out) :: basin_sz(nb_htu)       !!
+    INTEGER(i_std), INTENT(out) :: basin_bxout(nb_htu)    !!
+    REAL(r_std), INTENT(out)    :: basin_lshead(nb_htu)   !!
+    INTEGER(i_std), INTENT(out) :: basin_bbout(nb_htu)    !!
+    INTEGER(i_std), INTENT(out) :: basin_pts(nb_htu, nbv, 2) !!
+    INTEGER(i_std), INTENT(out) :: coast_pts(nb_htu)      !! The coastal flow points (unitless)
     !
 !! LOCAL VARIABLES
     LOGICAL, PARAMETER :: debug=.TRUE.
@@ -414,16 +415,16 @@ CONTAINS
     INTEGER(i_std) :: il, jl, ill, jll, ij, iz!!
     INTEGER(i_std) :: totsz, checksz !!
     !
-    INTEGER(i_std) :: tbname(nbvmax) !!
-    INTEGER(i_std) :: tsz(nbvmax) !!
-    INTEGER(i_std) :: tpts(nbvmax,nbvmax,2) !!
-    INTEGER(i_std) :: toutdir(nbvmax) !!
-    INTEGER(i_std) :: toutbas(nbvmax) !!
-    INTEGER(i_std) :: toutloc(nbvmax,2) !!
-    REAL(r_std) :: tlon(nbvmax) !!
-    REAL(r_std) :: tlat(nbvmax) !!
-    REAL(r_std) :: touttp(nbvmax) !!
-    REAL(r_std) :: toutlshead(nbvmax) !!
+    INTEGER(i_std) :: tbname(nb_htu) !!
+    INTEGER(i_std) :: tsz(nb_htu) !!
+    INTEGER(i_std) :: tpts(nb_htu,nbv,2) !!
+    INTEGER(i_std) :: toutdir(nb_htu) !!
+    INTEGER(i_std) :: toutbas(nb_htu) !!
+    INTEGER(i_std) :: toutloc(nb_htu,2) !!
+    REAL(r_std) :: tlon(nb_htu) !!
+    REAL(r_std) :: tlat(nb_htu) !!
+    REAL(r_std) :: touttp(nb_htu) !!
+    REAL(r_std) :: toutlshead(nb_htu) !!
     !
     INTEGER(i_std) :: tmpsz(nbvmax) !!
     INTEGER(i_std) :: ip, jp, jpp(1), ipb !!
@@ -435,15 +436,15 @@ CONTAINS
     !
     INTEGER(i_std) :: new_nb !! Number of sub-basins (unitless)
     INTEGER(i_std) :: mp !! Number of dividing sub-basin in mainstream (unitless)
-    INTEGER(i_std) :: new_bname(nbvmax) !!
-    INTEGER(i_std) :: new_outdir(nbvmax) !!
-    REAL(r_std)    :: new_heading(nbvmax)!!
-    INTEGER(i_std) :: new_outbas(nbvmax) !!
-    INTEGER(i_std) :: new_outloc(nbvmax,2) !!
-    REAL(r_std)    :: new_lon(nbvmax) !!
-    REAL(r_std)    :: new_lat(nbvmax) !!
-    INTEGER(i_std) :: new_sz(nbvmax) !!
-    INTEGER(i_std) :: new_pts(nbvmax, nbvmax, 2) !!
+    INTEGER(i_std) :: new_bname(nb_htu) !!
+    INTEGER(i_std) :: new_outdir(nb_htu) !!
+    REAL(r_std)    :: new_heading(nb_htu)!!
+    INTEGER(i_std) :: new_outbas(nb_htu) !!
+    INTEGER(i_std) :: new_outloc(nb_htu,2) !!
+    REAL(r_std)    :: new_lon(nb_htu) !!
+    REAL(r_std)    :: new_lat(nb_htu) !!
+    INTEGER(i_std) :: new_sz(nb_htu) !!
+    INTEGER(i_std) :: new_pts(nb_htu, nbv, 2) !!
     INTEGER(i_std) :: oldorder, neworder !!
     !
     INTEGER(i_std) :: option !! Option to calculate topoindex
@@ -694,7 +695,8 @@ CONTAINS
           !
           ibas = trans(ipb)
           !
-          CALL routing_reg_divbas(nbi, nbj, ibas, toutloc(ibas,1), toutloc(ibas,2), tsz(ibas), toutbas(ibas), toutdir(ibas), &
+          CALL routing_reg_divbas(nb_htu, nbv, nbi, nbj, ibas, toutloc(ibas,1), toutloc(ibas,2),&
+             & tsz(ibas), toutbas(ibas), toutdir(ibas), &
              & toutlshead(ibas), tpts, trip, basin, fac, lontmp, lattmp, &
              & new_nb, mp, new_bname, new_outdir, new_heading, new_outbas, new_lat, new_lon, new_outloc, new_sz, new_pts)
           !
@@ -1031,13 +1033,14 @@ CONTAINS
 !! \n
 !_ ================================================================================================================================
 
-  SUBROUTINE routing_reg_divbas(nbi, nbj, ibas, iloc, jloc, tsz, tout, toutd, headd, &
+  SUBROUTINE routing_reg_divbas(nb_htu, nbv, nbi, nbj, ibas, iloc, jloc, tsz, tout, toutd, headd, &
                 & tpts, trip, basin, fac, lon, lat,      &
                 & new_nb, oic, new_bname, new_outdir, new_head, new_outbas, new_lat, new_lon, new_outloc, new_sz, new_pts)
     !
     IMPLICIT NONE
     !
 !! INPUT VARIABLES
+    INTEGER(i_std), INTENT(in) :: nb_htu, nbv
     INTEGER(i_std), INTENT(in) :: nbi !! Number of point in x within the grid (unitless)
     INTEGER(i_std), INTENT(in) :: nbj !! Number of point in y within the grid (unitless)
     INTEGER(i_std), INTENT(in) :: ibas !! Order of the basin will be divided
@@ -1058,15 +1061,15 @@ CONTAINS
 !! OUTPUT VARIABLES
     INTEGER(i_std), INTENT(out) :: new_nb !! Number of sub-basins (unitless)
     INTEGER(i_std), INTENT(out) :: oic
-    INTEGER(i_std), INTENT(out) :: new_bname(nbvmax) !!
-    INTEGER(i_std), INTENT(out) :: new_outdir(nbvmax) !!
-    REAL(r_std), INTENT(out)    :: new_head(nbvmax) !!
-    INTEGER(i_std), INTENT(out) :: new_outbas(nbvmax) !!
-    INTEGER(i_std), INTENT(out) :: new_outloc(nbvmax,2) !!
-    REAL(r_std), INTENT(out) :: new_lon(nbvmax) !!
-    REAL(r_std), INTENT(out) :: new_lat(nbvmax) !!
-    INTEGER(i_std), INTENT(out) :: new_sz(nbvmax) !!
-    INTEGER(i_std), INTENT(out) :: new_pts(nbvmax, nbvmax, 2) !!
+    INTEGER(i_std), INTENT(out) :: new_bname(nb_htu) !!
+    INTEGER(i_std), INTENT(out) :: new_outdir(nb_htu) !!
+    REAL(r_std), INTENT(out)    :: new_head(nb_htu) !!
+    INTEGER(i_std), INTENT(out) :: new_outbas(nb_htu) !!
+    INTEGER(i_std), INTENT(out) :: new_outloc(nb_htu,2) !!
+    REAL(r_std), INTENT(out) :: new_lon(nb_htu) !!
+    REAL(r_std), INTENT(out) :: new_lat(nb_htu) !!
+    INTEGER(i_std), INTENT(out) :: new_sz(nb_htu) !!
+    INTEGER(i_std), INTENT(out) :: new_pts(nb_htu, nbv, 2) !!
     !
 !! LOCAL VARIABLES
     REAL(r_std), PARAMETER :: flag=-9999. !!
@@ -1116,13 +1119,13 @@ CONTAINS
     !    river, and get the tributaries. To store location and flow accumulation
     !    of points belongs to main stream and tributaries:
     !
-    INTEGER(i_std) :: main_loc(nbvmax,2), tri_loc(nbne,nbvmax,2)
-    REAL(r_std) :: main_fac(nbvmax), tri_fac(nbne,nbvmax)
+    INTEGER(i_std) :: main_loc(nbv,2), tri_loc(nbne,nbv,2)
+    REAL(r_std) :: main_fac(nbv), tri_fac(nbne,nbv)
     REAL(r_std) :: tmptri_fac(nbne)  ! Flow accumulation of all neighbour points
     INTEGER(i_std) :: sortedtrifac(nbne) ! And sort of tmptri_fac(nbne)
-    REAL(r_std) :: alltri_fac(nbvmax) ! Sort flow accumulation of all tributaries
-    INTEGER(i_std) :: alltri_loc1(nbvmax), alltri_loc2(nbvmax) ! And their location
-    INTEGER(i_std) :: sortedallfac(nbvmax) ! Sort alltri_fac(nbvmax)
+    REAL(r_std) :: alltri_fac(nbv) ! Sort flow accumulation of all tributaries
+    INTEGER(i_std) :: alltri_loc1(nbv), alltri_loc2(nbv) ! And their location
+    INTEGER(i_std) :: sortedallfac(nbv) ! Sort alltri_fac(nbv)
     !
     !     *  Note: again, there are few dimensions should be exactly [nbne-1]
     !     or [nbne-2]. Because when you follow upstream the mainstream, there are
@@ -1646,7 +1649,7 @@ CONTAINS
 !! \n
 !_ ================================================================================================================================
 
-SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_bx, &
+SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, neighbours, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_bx, &
        & orog_bx, floodp_bx, fac_bx, topoind_bx, &
        & min_topoind, nb_basin, basin_inbxid, basin_outlet, basin_outtp, basin_sz, basin_pts, basin_bxout, &
        & basin_bbout, lshead, coast_pts, nwbas, basin_count, basin_notrun, basin_area, basin_cg, basin_hierarchy, &
@@ -1657,7 +1660,7 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx,
     IMPLICIT NONE
     !
 !! INPUT VARIABLES
-    INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
+    INTEGER(i_std), INTENT (in) :: nbpt, nb_htu, nbv !! Domain size (unitless)
     INTEGER(i_std), INTENT (in) :: ib !! Current basin (unitless)
     INTEGER(i_std), INTENT(in)  :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
                                                                          !! (1=North and then clockwise)
@@ -1673,14 +1676,14 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx,
     REAL(r_std), DIMENSION(nbxmax,nbxmax)    :: topoind_bx   !! Topographic index of the residence time for each of the smaller boxes (m)
     REAL(r_std)                              :: min_topoind  !! The current minimum of topographic index (m)
     INTEGER(i_std)                           :: nb_basin     !! Number of sub-basins (unitless)
-    INTEGER(i_std), DIMENSION(nbvmax)        :: basin_inbxid, basin_sz !! ID of basin, number of points in the basin
-    REAL(r_std), DIMENSION(nbvmax,2)         :: basin_outlet !! Outlet coordinate of each subgrid basin (lat,lon)
+    INTEGER(i_std), DIMENSION(nb_htu)        :: basin_inbxid, basin_sz !! ID of basin, number of points in the basin
+    REAL(r_std), DIMENSION(nb_htu,2)         :: basin_outlet !! Outlet coordinate of each subgrid basin (lat,lon)
     REAL(r_std), DIMENSION(nbvmax)           :: basin_outtp  !!
-    INTEGER(i_std), DIMENSION(nbvmax,nbvmax,2) :: basin_pts  !! Points in each basin
-    INTEGER(i_std), DIMENSION(nbvmax)        :: basin_bxout  !! outflow direction
-    INTEGER(i_std), DIMENSION(nbvmax)        :: basin_bbout  !! outflow sub-basin
-    REAL(r_std), DIMENSION(nbvmax)           :: lshead       !! Large scale heading of outflow.
-    INTEGER(i_std)                           :: coast_pts(nbvmax) !! The coastal flow points (unitless)
+    INTEGER(i_std), DIMENSION(nb_htu,nbv,2) :: basin_pts  !! Points in each basin
+    INTEGER(i_std), DIMENSION(nb_htu)        :: basin_bxout  !! outflow direction
+    INTEGER(i_std), DIMENSION(nb_htu)        :: basin_bbout  !! outflow sub-basin
+    REAL(r_std), DIMENSION(nb_htu)           :: lshead       !! Large scale heading of outflow.
+    INTEGER(i_std)                           :: coast_pts(nb_htu) !! The coastal flow points (unitless)
     ! global maps
     INTEGER(i_std) :: nwbas !!
     INTEGER(i_std), DIMENSION(nbpt) :: basin_count         !!
diff --git a/Interface.py b/Interface.py
index adf4d7452ef7e884d37b70f2f2dc256810ffa21e..be42d37af04d9f8035b85a81b100c1b9bf3e0a84 100644
--- a/Interface.py
+++ b/Interface.py
@@ -286,21 +286,31 @@ class HydroOverlap :
 #
 #
 class HydroSuper :
-    def __init__(self, nbvmax, hydrodata, hydrooverlap, nbasmax) :
+    def __init__(self, nbvmax, hydrodata, hydrooverlap, nbasmax, part) :
         #
         # Keep largest possible number of HTUs
         #
         self.nbasmax = nbasmax
         self.nbhtuext = nbvmax
-        self.nwbas = hydrooverlap.nwbas
+        self.nbpt = hydrooverlap.nbi.shape[0]
+        #
+        # nb_htu can be adjusted with self.nwbas
+        # nb_htu can be lowered with a larger maxpercent (routing_reg.f90)
+        nb_htu = 600 
+        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(nbvmax, hydrooverlap.nbi, hydrooverlap.nbj, hydrooverlap.trip_bx, \
-                                                 hydrooverlap.basin_bx, hydrooverlap.fac_bx, hydrooverlap.hierarchy_bx, \
-                                                 hydrooverlap.topoind_bx, hydrooverlap.lshead_bx, \
-                                                 hydrooverlap.lon_bx, hydrooverlap.lat_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)
+        #
+        # Adjust nwbas to the maximum found over the domain
+        #
+        self.nwbas = part.domainmax(np.max(nb_basin))
+        print("Maximum number of basin created : {0}".format(self.nwbas))
         #
         # Call Globalize
         #
@@ -312,14 +322,15 @@ 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(hydrooverlap.area_bx, lon_bx_tmp, lat_bx_tmp, hydrooverlap.trip_bx, \
-                                                hydrooverlap.hierarchy_bx, hydrooverlap.orog_bx, hydrooverlap.floodp_bx,\
-                                                hydrooverlap.fac_bx, hydrooverlap.topoind_bx, hydrodata.topoindmin, \
-                                                nb_basin, basin_inbxid, basin_outlet, basin_outtp, self.basin_sz, self.basin_pts, basin_bxout, \
-                                                basin_bbout, basin_lshead, coast_pts, hydrooverlap.nwbas)
-
-        self.nbpt = self.basin_count.shape[0]
+                    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, \
+                                                basin_bbout = basin_bbout, lshead = basin_lshead, coast_pts = coast_pts, nwbas = self.nwbas)
 
+        # Memory management
+        del basin_bbout; del basin_lshead; del coast_pts; del basin_bxout; del self.basin_pts;
+        del basin_outtp; del basin_outlet; del basin_inbxid; del nb_basin
         return
     #
     def linkup(self, hydrodata) :
diff --git a/RoutingPreProc.py b/RoutingPreProc.py
index c97d29534646caa3132cf6faa898a5c831fe5de0..52cf319b3dbb9c85a185551c4ba2ca0f60a95389 100644
--- a/RoutingPreProc.py
+++ b/RoutingPreProc.py
@@ -83,11 +83,11 @@ hoverlap = IF.HydroOverlap(nbpt, nbvmax, w.hpts, w.index, w.area, w.lon, w.lat,
 #
 # Do some memory management and synchronize procs.
 #
-comm.Barrier()
 del w
 gc.collect()
+comm.Barrier()
 
-hsuper = IF.HydroSuper(nbvmax, hydrodata, hoverlap, nbasmax)
+hsuper = IF.HydroSuper(nbvmax, hydrodata, hoverlap, nbasmax, part)
 #
 INFO("=================== LINKUP ====================")
 hsuper.linkup(hydrodata)
@@ -95,10 +95,10 @@ hsuper.linkup(hydrodata)
 #
 # Do some memory management and synchronize procs.
 #
- 
-comm.Barrier()
+
 del hoverlap
-gc.collect()
+gc.collect() 
+comm.Barrier()
 
 INFO("=================== Compute fetch ====================")
 t = time.time()