diff --git a/DocumentationInterface b/DocumentationInterface
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0a431bb3357b8b4035f82d5e536d4a6d29583352 100644
--- a/DocumentationInterface
+++ b/DocumentationInterface
@@ -0,0 +1,264 @@
+initatmgrid(rank_bn,nbcore,polygons_in,centre_in,area_in,contfrac_in,neighbours_in,[nbpt,nbseg])
+
+Wrapper for ``initatmgrid``.
+
+Parameters
+----------
+rank_bn : input int
+nbcore : input int
+polygons_in : input rank-3 array('f') with bounds (nbpt,2 * nbseg + 1,2)
+centre_in : input rank-2 array('f') with bounds (nbpt,2)
+area_in : input rank-1 array('f') with bounds (nbpt)
+contfrac_in : input rank-1 array('f') with bounds (nbpt)
+neighbours_in : input rank-2 array('i') with bounds (nbpt,2 * nbseg)
+
+Other Parameters
+----------------
+nbpt : input int, optional
+ Default: shape(polygons_in,0)
+nbseg : input int, optional
+ Default: (shape(polygons_in,1)-1)/(2)
+====================================================================
+nbi,nbj,area_bx,trip_bx,basin_bx,topoind_bx,fac_bx,hierarchy_bx,lon_bx,lat_bx,lshead_bx,nwbas = gethydrogrid(nbxmax_in,sub_pts,sub_index,sub_area,max_basins,min_topoind,lon_rel,lat_rel,trip,basins,topoindex,fac,hierarchy,[nbpt,nbvmax_in])
+
+Wrapper for ``gethydrogrid``.
+
+Parameters
+----------
+nbxmax_in : input int
+sub_pts : input rank-1 array('i') with bounds (nbpt)
+sub_index : input rank-3 array('i') with bounds (nbpt,nbvmax_in,2)
+sub_area : input rank-2 array('f') with bounds (nbpt,nbvmax_in)
+max_basins : in/output rank-0 array(float,'f')
+min_topoind : input float
+lon_rel : input rank-2 array('f') with bounds (nbpt,nbvmax_in)
+lat_rel : input rank-2 array('f') with bounds (nbpt,nbvmax_in)
+trip : in/output rank-2 array('f') with bounds (nbpt,nbvmax_in)
+basins : in/output rank-2 array('f') with bounds (nbpt,nbvmax_in)
+topoindex : in/output rank-2 array('f') with bounds (nbpt,nbvmax_in)
+fac : in/output rank-2 array('f') with bounds (nbpt,nbvmax_in)
+hierarchy : in/output rank-2 array('f') with bounds (nbpt,nbvmax_in)
+
+Other Parameters
+----------------
+nbpt : input int, optional
+ Default: len(sub_pts)
+nbvmax_in : input int, optional
+ Default: shape(sub_index,1)
+
+Returns
+-------
+nbi : rank-1 array('i') with bounds (nbpt)
+nbj : rank-1 array('i') with bounds (nbpt)
+area_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+trip_bx : rank-3 array('i') with bounds (nbpt,nbxmax_in,nbxmax_in)
+basin_bx : rank-3 array('i') with bounds (nbpt,nbxmax_in,nbxmax_in)
+topoind_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+fac_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+hierarchy_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+lon_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+lat_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+lshead_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+nwbas : int
+====================================================================
+nb_basin,basin_inbxid,basin_outlet,basin_outtp,basin_sz,basin_bxout,basin_bbout,basin_pts,basin_lshead,coast_pts = findbasins(nbvmax_in,nbi,nbj,trip_bx,basin_bx,fac_bx,hierarchy_bx,topoind_bx,lshead_bx,lontmp,lattmp,[nbpt,nbxmax_in])
+
+Wrapper for ``findbasins``.
+
+Parameters
+----------
+nbvmax_in : input int
+nbi : input rank-1 array('i') with bounds (nbpt)
+nbj : input rank-1 array('i') with bounds (nbpt)
+trip_bx : in/output rank-3 array('i') with bounds (nbpt,nbxmax_in,nbxmax_in)
+basin_bx : in/output rank-3 array('i') with bounds (nbpt,nbxmax_in,nbxmax_in)
+fac_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+hierarchy_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+topoind_bx : in/output rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+lshead_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+lontmp : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+lattmp : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+
+Other Parameters
+----------------
+nbpt : input int, optional
+ Default: len(nbi)
+nbxmax_in : input int, optional
+ Default: shape(trip_bx,1)
+
+Returns
+-------
+nb_basin : rank-1 array('i') with bounds (nbpt)
+basin_inbxid : rank-2 array('i') with bounds (nbpt,nbvmax_in)
+basin_outlet : rank-3 array('f') with bounds (nbpt,nbvmax_in,2)
+basin_outtp : rank-2 array('f') with bounds (nbpt,nbvmax_in)
+basin_sz : rank-2 array('i') with bounds (nbpt,nbvmax_in)
+basin_bxout : rank-2 array('i') with bounds (nbpt,nbvmax_in)
+basin_bbout : rank-2 array('i') with bounds (nbpt,nbvmax_in)
+basin_pts : rank-4 array('i') with bounds (nbpt,nbvmax_in,nbvmax_in,2)
+basin_lshead : rank-2 array('f') with bounds (nbpt,nbvmax_in)
+coast_pts : rank-2 array('i') with bounds (nbpt,nbvmax_in)
+====================================================================
+basin_count,basin_notrun,basin_area,basin_cg,basin_hierarchy,basin_fac,basin_topoind,basin_id,basin_coor,basin_type,basin_flowdir,basin_lshead,outflow_grid,outflow_basin,nbcoastal,coastal_basin = globalize(area_bx,lon_bx,lat_bx,trip_bx,hierarchy_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,[nbpt,nbvmax_in,nbxmax_in])
+
+Wrapper for ``globalize``.
+
+Parameters
+----------
+area_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+lon_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+lat_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+trip_bx : input rank-3 array('i') with bounds (nbpt,nbxmax_in,nbxmax_in)
+hierarchy_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+fac_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+topoind_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+min_topoind : input float
+nb_basin : input rank-1 array('i') with bounds (nbpt)
+basin_inbxid : input rank-2 array('i') with bounds (nbpt,nbvmax_in)
+basin_outlet : input rank-3 array('f') with bounds (nbpt,nbvmax_in,2)
+basin_outtp : input rank-2 array('f') with bounds (nbpt,nbvmax_in)
+basin_sz : input rank-2 array('i') with bounds (nbpt,nbvmax_in)
+basin_pts : input rank-4 array('i') with bounds (nbpt,nbvmax_in,nbvmax_in,2)
+basin_bxout : input rank-2 array('i') with bounds (nbpt,nbvmax_in)
+basin_bbout : input rank-2 array('i') with bounds (nbpt,nbvmax_in)
+lshead : input rank-2 array('f') with bounds (nbpt,nbvmax_in)
+coast_pts : input rank-2 array('i') with bounds (nbpt,nbvmax_in)
+nwbas : input int
+
+Other Parameters
+----------------
+nbpt : input int, optional
+ Default: shape(area_bx,0)
+nbvmax_in : input int, optional
+ Default: shape(basin_inbxid,1)
+nbxmax_in : input int, optional
+ Default: shape(area_bx,1)
+
+Returns
+-------
+basin_count : rank-1 array('i') with bounds (nbpt)
+basin_notrun : rank-1 array('i') with bounds (nbpt)
+basin_area : rank-2 array('f') with bounds (nbpt,nwbas)
+basin_cg : rank-3 array('f') with bounds (nbpt,nwbas,2)
+basin_hierarchy : rank-2 array('f') with bounds (nbpt,nwbas)
+basin_fac : rank-2 array('f') with bounds (nbpt,nwbas)
+basin_topoind : rank-2 array('f') with bounds (nbpt,nwbas)
+basin_id : rank-2 array('i') with bounds (nbpt,nwbas)
+basin_coor : rank-3 array('f') with bounds (nbpt,nwbas,2)
+basin_type : rank-2 array('f') with bounds (nbpt,nwbas)
+basin_flowdir : rank-2 array('i') with bounds (nbpt,nwbas)
+basin_lshead : rank-2 array('f') with bounds (nbpt,nwbas)
+outflow_grid : rank-2 array('i') with bounds (nbpt,nwbas)
+outflow_basin : rank-2 array('i') with bounds (nbpt,nwbas)
+nbcoastal : rank-1 array('i') with bounds (nbpt)
+coastal_basin : rank-2 array('i') with bounds (nbpt,nwbas)
+====================================================================
+inflow_number,inflow_grid,inflow_basin = linkup(nbxmax_in,basin_count,basin_area,basin_id,basin_flowdir,basin_lshead,basin_hierarchy,basin_fac,outflow_grid,outflow_basin,nbcoastal,coastal_basin,invented_basins,[nbpt,nwbas])
+
+Wrapper for ``linkup``.
+
+Parameters
+----------
+nbxmax_in : input int
+basin_count : input rank-1 array('i') with bounds (nbpt)
+basin_area : input rank-2 array('f') with bounds (nbpt,nwbas)
+basin_id : input rank-2 array('i') with bounds (nbpt,nwbas)
+basin_flowdir : input rank-2 array('i') with bounds (nbpt,nwbas)
+basin_lshead : input rank-2 array('f') with bounds (nbpt,nwbas)
+basin_hierarchy : in/output rank-2 array('f') with bounds (nbpt,nwbas)
+basin_fac : input rank-2 array('f') with bounds (nbpt,nwbas)
+outflow_grid : in/output rank-2 array('i') with bounds (nbpt,nwbas)
+outflow_basin : in/output rank-2 array('i') with bounds (nbpt,nwbas)
+nbcoastal : in/output rank-1 array('i') with bounds (nbpt)
+coastal_basin : in/output rank-2 array('i') with bounds (nbpt,nwbas)
+invented_basins : input float
+
+Other Parameters
+----------------
+nbpt : input int, optional
+ Default: len(basin_count)
+nwbas : input int, optional
+ Default: shape(basin_area,1)
+
+Returns
+-------
+inflow_number : rank-2 array('i') with bounds (nbpt,8 * nbxmax_in)
+inflow_grid : rank-3 array('i') with bounds (nbpt,8 * nbxmax_in,8 * nbxmax_in)
+inflow_basin : rank-3 array('i') with bounds (nbpt,8 * nbxmax_in,8 * nbxmax_in)
+====================================================================
+outflow_uparea = fetch(corepts,basin_count,basin_area,basin_id,basin_hierarchy,basin_fac,outflow_grid,outflow_basin,partial_sum,fetch_basin,[nbpt,nwbas,nbcore])
+
+Wrapper for ``fetch``.
+
+Parameters
+----------
+corepts : input rank-1 array('i') with bounds (nbcore)
+basin_count : input rank-1 array('i') with bounds (nbpt)
+basin_area : in/output rank-2 array('f') with bounds (nbpt,nwbas)
+basin_id : input rank-2 array('i') with bounds (nbpt,nwbas)
+basin_hierarchy : input rank-2 array('f') with bounds (nbpt,nwbas)
+basin_fac : input rank-2 array('f') with bounds (nbpt,nwbas)
+outflow_grid : input rank-2 array('i') with bounds (nbpt,nwbas)
+outflow_basin : input rank-2 array('i') with bounds (nbpt,nwbas)
+partial_sum : input rank-2 array('f') with bounds (nbpt,nwbas)
+fetch_basin : in/output rank-2 array('f') with bounds (nbpt,nwbas)
+
+Other Parameters
+----------------
+nbpt : input int, optional
+ Default: len(basin_count)
+nwbas : input int, optional
+ Default: shape(basin_area,1)
+nbcore : input int, optional
+ Default: len(corepts)
+
+Returns
+-------
+outflow_uparea : rank-1 array('f') with bounds (nbpt*nwbas)
+====================================================================
+routing_area,routing_cg,topo_resid,route_togrid,route_tobasin,route_nbintobas,global_basinid,route_outlet,route_type,origin_nbintobas,routing_fetch = truncate(nbasmax,num_largest,basin_count,basin_notrun,basin_area,basin_cg,basin_topoind,fetch_basin,basin_id,basin_coor,basin_type,basin_flowdir,outflow_grid,outflow_basin,inflow_number,inflow_grid,inflow_basin,[nbpt,nbxmax_in,nwbas])
+
+Wrapper for ``truncate``.
+
+Parameters
+----------
+nbasmax : input int
+num_largest : input int
+basin_count : in/output rank-1 array('i') with bounds (nbpt)
+basin_notrun : in/output rank-1 array('i') with bounds (nbpt)
+basin_area : in/output rank-2 array('f') with bounds (nbpt,nwbas)
+basin_cg : in/output rank-3 array('f') with bounds (nbpt,nwbas,2)
+basin_topoind : in/output rank-2 array('f') with bounds (nbpt,nwbas)
+fetch_basin : in/output rank-2 array('f') with bounds (nbpt,nwbas)
+basin_id : in/output rank-2 array('i') with bounds (nbpt,nwbas)
+basin_coor : in/output rank-3 array('f') with bounds (nbpt,nwbas,2)
+basin_type : in/output rank-2 array('f') with bounds (nbpt,nwbas)
+basin_flowdir : in/output rank-2 array('i') with bounds (nbpt,nwbas)
+outflow_grid : in/output rank-2 array('i') with bounds (nbpt,nwbas)
+outflow_basin : in/output rank-2 array('i') with bounds (nbpt,nwbas)
+inflow_number : in/output rank-2 array('i') with bounds (nbpt,8 * nbxmax_in)
+inflow_grid : in/output rank-3 array('i') with bounds (nbpt,8 * nbxmax_in,8 * nbxmax_in)
+inflow_basin : in/output rank-3 array('i') with bounds (nbpt,8 * nbxmax_in,8 * nbxmax_in)
+
+Other Parameters
+----------------
+nbpt : input int, optional
+ Default: len(basin_count)
+nbxmax_in : input int, optional
+ Default: (shape(inflow_number,1))/(8)
+nwbas : input int, optional
+ Default: shape(basin_area,1)
+
+Returns
+-------
+routing_area : rank-2 array('f') with bounds (nbpt,nbasmax)
+routing_cg : rank-3 array('f') with bounds (nbpt,nbasmax,2)
+topo_resid : rank-2 array('f') with bounds (nbpt,nbasmax)
+route_togrid : rank-2 array('i') with bounds (nbpt,nbasmax)
+route_tobasin : rank-2 array('i') with bounds (nbpt,nbasmax)
+route_nbintobas : rank-1 array('i') with bounds (nbpt)
+global_basinid : rank-2 array('i') with bounds (nbpt,nbasmax)
+route_outlet : rank-3 array('f') with bounds (nbpt,nbasmax,2)
+route_type : rank-2 array('f') with bounds (nbpt,nbasmax)
+origin_nbintobas : rank-1 array('i') with bounds (nbpt)
+routing_fetch : rank-2 array('f') with bounds (nbpt,nbasmax)
diff --git a/F90subroutines/routing_interface.f90 b/F90subroutines/routing_interface.f90
index 67d52162c3098a9593e6a44870d8c9f3e96e4a99..59886d09cbac7e542237e3070123ed38b297c46d 100644
--- a/F90subroutines/routing_interface.f90
+++ b/F90subroutines/routing_interface.f90
@@ -286,7 +286,7 @@ SUBROUTINE linkup(nbpt, nbxmax_in, nwbas, basin_count, basin_area, basin_id, bas
!
END SUBROUTINE linkup
-SUBROUTINE areanorm(nbpt, nwbas, basin_count, basin_area, basin_area_norm)
+SUBROUTINE areanorm(nbpt, nwbas, basin_count, basin_area, outflow_grid, basin_area_norm)
!
USE grid
!
@@ -296,7 +296,8 @@ SUBROUTINE areanorm(nbpt, nwbas, basin_count, basin_area, basin_area_norm)
INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_area !!
!
- REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(out) :: basin_area_norm !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless)
+ REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(out) :: basin_area_norm !!
!
!! LOCAL VARIABLES
INTEGER(i_std) :: ib, ij
@@ -312,13 +313,22 @@ SUBROUTINE areanorm(nbpt, nwbas, basin_count, basin_area, basin_area_norm)
!
DO ij=1,basin_count(ib)
basin_area_norm(ib,ij) = basin_area(ib,ij)/totbasins*contarea
+ !
+ ! Simplify the outflow condition. Large rivers will be reset later in rivclassification.
+ ! We set all outflow points to coastal flow. This will be adjusted later once all procs have
+ ! exchanged their information. So we will have outflow_grid = -2.
+ !
+ IF ( outflow_grid(ib,ij) .EQ. -1) THEN
+ outflow_grid(ib,ij) = -2
+ ENDIF
+
ENDDO
!
ENDDO
!
END SUBROUTINE areanorm
-SUBROUTINE fetch(nbpt, nwbas, basin_count, basin_area, basin_id, basin_hierarchy, basin_fac, &
+SUBROUTINE fetch(nbpt, nwbas, nbcore, corepts, basin_count, basin_area, basin_id, basin_hierarchy, basin_fac, &
& outflow_grid, outflow_basin, partial_sum, fetch_basin, outflow_uparea)
!
USE ioipsl
@@ -331,22 +341,34 @@ SUBROUTINE fetch(nbpt, nwbas, basin_count, basin_area, basin_id, basin_hierarchy
!! INPUT VARIABLES
INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
INTEGER(i_std), INTENT (in) :: nwbas !!
+ INTEGER(i_std), INTENT (in) :: nbcore
+ INTEGER(i_std), DIMENSION(nbcore), INTENT(in) :: corepts
INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_area !!
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_id !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_hierarchy
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_fac
+ INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: outflow_grid !! Type of outflow on the grid box (unitless)
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: outflow_basin !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: partial_sum
!
!! IN-OUTPUT VARIABLES
- INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless)
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: fetch_basin !!
! Output
REAL(r_std), DIMENSION(nbpt*nwbas), INTENT(out) :: outflow_uparea
!
- CALL routing_reg_fetch(nbpt, area, contfrac, nwbas, basin_count, basin_area, basin_id, basin_hierarchy, basin_fac,&
- & outflow_grid, outflow_basin, partial_sum, fetch_basin, outflow_uparea)
+ !
+ !! Local
+ INTEGER(i_std) :: ic
+ INTEGER(i_std), DIMENSION(nbcore) :: fcorepts
+ !
+ !
+ DO ic=1,nbcore
+ fcorepts(ic) = corepts(ic)+1
+ ENDDO
+ !
+ CALL routing_reg_fetch(nbpt, area, contfrac, nwbas, nbcore, fcorepts, basin_count, basin_area, basin_id, &
+ & basin_hierarchy, basin_fac, outflow_grid, outflow_basin, partial_sum, fetch_basin, outflow_uparea)
!
END SUBROUTINE fetch
diff --git a/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index 6ed789ebc8fadeac5b4507d0dee950bb068d9315..72bd5df750dd8ba2e50146a83a50ec958e47f513 100644
--- a/F90subroutines/routing_reg.f90
+++ b/F90subroutines/routing_reg.f90
@@ -2690,8 +2690,8 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
!! \n
!_ ================================================================================================================================
-SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, basin_count, basin_area, basin_id, basin_hierarchy, basin_fac,&
- & outflow_grid, outflow_basin, partial_sum, fetch_basin, outflow_uparea)
+SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, basin_count, basin_area, basin_id, &
+ & basin_hierarchy, basin_fac, outflow_grid, outflow_basin, partial_sum, fetch_basin, outflow_uparea)
!
IMPLICIT NONE
!
@@ -2701,22 +2701,24 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, basin_count, basin
REAL(r_std), DIMENSION(nbpt), INTENT(in) :: gridarea !! The size of each grid box in X and Y (m)
REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of land in each grid box (unitless;0-1)
!
- INTEGER(i_std) :: nwbas !!
+ INTEGER(i_std) :: nwbas !!
+ INTEGER(i_std), INTENT (in) :: nbcore
+ INTEGER(i_std), DIMENSION(nbcore), INTENT(in) :: corepts
INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !!
- REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_area !!
+ REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_area !!
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_id !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_hierarchy
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_fac
+ INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: outflow_grid !! Type of outflow on the grid box (unitless)
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: outflow_basin !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: partial_sum
!
!! INPUT/OUTPUT VARIABLES
- INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless)
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: fetch_basin !!
REAL(r_std), DIMENSION(nbpt*nwbas), INTENT(out) :: outflow_uparea
!
!! LOCAL VARIABLES
- INTEGER(i_std) :: ib, ij, ff(1), it, itt, igrif, ibasf, nboutflow !! Indices (unitless)
+ INTEGER(i_std) :: ib, ij, ic, ff(1), it, itt, igrif, ibasf, nboutflow !! Indices (unitless)
REAL(r_std), DIMENSION(nbpt*nbvmax) :: tmp_area !!
INTEGER(i_std), DIMENSION(nbpt*nbvmax,2) :: tmpindex !!
!
@@ -2727,10 +2729,6 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, basin_count, basin
! Propagate first the partial sums from the halo into the domain
!
DO ib=1,nbpt
- !
- IF ( SUM(partial_sum(ib,:)) > 0 ) THEN
- WRITE(numout,*) ib, "Have some partial system to add in point", lalo(ib,2), lalo(ib,1)
- ENDIF
!
DO ij=1,basin_count(ib)
!
@@ -2763,48 +2761,39 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, basin_count, basin
!
! Add the areas contributed by the core region of the domain.
!
- DO ib=1,nbpt
+ DO ic=1,nbcore
+ ib = corepts(ic)
!
DO ij=1,basin_count(ib)
!
- IF (partial_sum(ib,ij) <= 0.01) THEN
- !
- fetch_basin(ib, ij) = fetch_basin(ib, ij) + basin_area(ib,ij)
- !
- igrif = outflow_grid(ib,ij)
- ibasf = outflow_basin(ib,ij)
+ fetch_basin(ib, ij) = fetch_basin(ib, ij) + basin_area(ib,ij)
+ !
+ igrif = outflow_grid(ib,ij)
+ ibasf = outflow_basin(ib,ij)
+ !
+ itt = 0
+ !
+ DO WHILE (igrif .GT. 0 )
!
- itt = 0
+ fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + basin_area(ib, ij)
!
- DO WHILE (igrif .GT. 0 )
- !
- IF (partial_sum(igrif,ibasf) <= 0.01) THEN
- ! In the core regions
- fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + basin_area(ib, ij)
- ELSE
- ! In the halo of the domain avoid double counting
- fetch_basin(igrif,ibasf) = MAX(fetch_basin(igrif,ibasf), partial_sum(ib, ij))
- ENDIF
+ it = outflow_grid(igrif, ibasf)
+ ibasf = outflow_basin(igrif, ibasf)
+ igrif = it
+ itt = itt + 1
+ IF ( itt .GT. maxiterations) THEN
+ WRITE(numout,&
+ "('Grid ',I5, ' and basin ',I5, ' did not converge after iteration ',I5)") ib, ij, itt
+ WRITE(numout, "('==> We flow through grid ',I5,' and basin ',I5)") igrif, ibasf
+ WRITE(numout,*) '==> Basin ID :', basin_id(igrif,ibasf), "Hierarchy :", basin_hierarchy(igrif,ibasf)
+ WRITE(numout,*) "==> Flow accumulation :", basin_fac(igrif,ibasf)
+ WRITE(numout,*) "==> Coordinates : ", lalo_g(igrif,2), lalo_g(igrif,1)
!
- it = outflow_grid(igrif, ibasf)
- ibasf = outflow_basin(igrif, ibasf)
- igrif = it
- itt = itt + 1
- IF ( itt .GT. maxiterations) THEN
- WRITE(numout,&
- "('Grid ',I5, ' and basin ',I5, ' did not converge after iteration ',I5)") ib, ij, itt
- WRITE(numout, "('==> We flow through grid ',I5,' and basin ',I5)") igrif, ibasf
- WRITE(numout,*) '==> Basin ID :', basin_id(igrif,ibasf), "Hierarchy :", basin_hierarchy(igrif,ibasf)
- WRITE(numout,*) "==> Flow accumulation :", basin_fac(igrif,ibasf)
- WRITE(numout,*) "==> Coordinates : ", lalo_g(igrif,2), lalo_g(igrif,1)
- !
- IF ( itt .GT. maxiterations+20) THEN
- CALL ipslerr_p(3,'routing_reg_fetch','Problem...','','')
- ENDIF
+ IF ( itt .GT. maxiterations+20) THEN
+ CALL ipslerr_p(3,'routing_reg_fetch','Problem...','','')
ENDIF
- ENDDO
- !
- ENDIF
+ ENDIF
+ ENDDO
!
ENDDO
!
@@ -2813,13 +2802,12 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, basin_count, basin
WRITE(numout,*) 'The smallest FETCH :', MINVAL(fetch_basin)
WRITE(numout,*) 'The largest FETCH :', MAXVAL(fetch_basin)
!
- ! We set all outflow points to coastal flow. This will be adjusted later once all procs have
- ! exchanged their information. So we will have outflow_grid = -2.
!
nboutflow = 0
outflow_uparea(:) = zero
!
- DO ib=1,nbpt
+ DO ic=1,nbcore
+ ib = corepts(ic)
!
DO ij=1,basin_count(ib)
!
@@ -2827,9 +2815,6 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, basin_count, basin
! We archive the upstream area of all outflow points so that we can sort them.
!
IF ( outflow_grid(ib,ij) .LT. 0) THEN
- IF ( outflow_grid(ib,ij) .EQ. -1) THEN
- outflow_grid(ib,ij) = -2
- ENDIF
nboutflow = nboutflow + 1
IF ( nboutflow <= nbpt*nwbas) THEN
outflow_uparea(nboutflow) = fetch_basin(ib,ij)
diff --git a/Interface.py b/Interface.py
index 69d5fe3f72d884b9a2f7aea28b57251a0c96ff46..d6b215d7f21ca2a0006925004cb6d2a3b3bb5535 100644
--- a/Interface.py
+++ b/Interface.py
@@ -162,30 +162,46 @@ class HydroSuper :
#
def fetch(self, part, modelgrid) :
#
- largest_pos = -50
+ largest_pos = 50
+ # Order of magnitude for the area precision in m^2.
+ prec = 100.0
#
- fetch_basin = np.zeros(self.basin_area.shape, dtype=np.float32, order='F')
+ fetch_basin = np.zeros(self.basin_area.shape, dtype=np.float64, order='F')
#
- self.basin_area = routing_interface.areanorm(self.basin_count, self.basin_area)
- partial_sum = np.zeros(self.basin_area.shape, dtype=np.float32, order='F')
+ self.basin_area = routing_interface.areanorm(self.basin_count, self.basin_area, self.outflow_grid)
+ partial_sum = np.zeros(self.basin_area.shape, dtype=np.float64, order='F')
#
- for i in range(part.size) :
+ old_sorted = np.zeros(largest_pos, dtype=np.float64, order='F')
+ #
+ maxdiff_sorted = prec*prec
+ iter_count = 0
+ #
+ while iter_count < part.size*3 and maxdiff_sorted > prec :
fetch_basin[:,:] = 0.0
- outflow_uparea = routing_interface.fetch(self.basin_count, self.basin_area, self.basin_id, self.basin_hierarchy, \
+ outflow_uparea = routing_interface.fetch(part.landcorelist, self.basin_count, self.basin_area, self.basin_id, self.basin_hierarchy, \
self.basin_fac, self.outflow_grid, self.outflow_basin, partial_sum, fetch_basin)
partial_sum = np.copy(fetch_basin)
part.landsendtohalo(partial_sum, order='F')
partial_sum = part.zerocore(partial_sum, order='F')
#
- yy=modelgrid.landscatter(np.sum(fetch_basin, axis=1)/np.sum(self.basin_area, axis=1))
+ # Find area the largest basins need at least to have.
+ #
+ xtmp = np.hstack(part.comm.allgather(outflow_uparea[np.where(outflow_uparea > 0.0)]))
+ # Precision in m^2 of the upstream areas when sorting.
+ sorted_outareas = (np.unique(np.rint(np.array(xtmp)/prec))*prec)[::-1]
+ # If mono-proc no need to iterate as fetch produces the full result.
+ if part.size == 1 :
+ maxdiff_sorted = 0.0
+ else :
+ maxdiff_sorted = np.max(np.abs(sorted_outareas[0:largest_pos]-old_sorted))
+ old_sorted[:] = sorted_outareas[0:largest_pos]
+ iter_count += 1
self.fetch_basin = np.copy(fetch_basin)
#
- # Find area the largest basins need at least to have.
+ # Upstream area of the smalest river we call largest rivers.
#
- xtmp = np.hstack(part.comm.allgather(outflow_uparea[np.where(outflow_uparea > 0.0)]))
- sorted_outareas = np.sort(np.array(xtmp))
- largest_rivarea = sorted_outareas[largest_pos]
+ largest_rivarea = sorted_outareas[largest_pos-1]
#
#
#
@@ -193,6 +209,7 @@ class HydroSuper :
print("Rank :"+str(part.rank)+" OUT of fetch =+=+=+=+=+=+=+=+=+= \n"+str(yy)+"\n =+=+=+=+=+=+=+=+=+=")
self.num_largest = routing_interface.rivclassification(self.basin_count, self.outflow_grid, self.outflow_basin, \
self.fetch_basin, largest_rivarea)
+ print("Rank :"+str(part.rank)+" Area of smallest large rivers : ", largest_rivarea, " Nb of Large rivers on proc : ",self.num_largest)
return
#
#
@@ -208,6 +225,7 @@ class HydroGraph :
hydrosuper.basin_outcoor, hydrosuper.basin_type, hydrosuper.basin_flowdir, \
hydrosuper.outflow_grid, hydrosuper.outflow_basin, \
hydrosuper.inflow_number,hydrosuper.inflow_grid,hydrosuper.inflow_basin)
+ print("Rank :"+str(part.rank)+" Out truncate diag")
yy=modelgrid.landscatter(np.sum(self.routing_fetch, axis=1)/np.sum(self.routing_area, axis=1))
print ("Rank :"+str(part.rank)+" OUT of truncate =+=+=+=+=+=+=+=+=+= \n"+str(yy)+"\n =+=+=+=+=+=+=+=+=+=")
return