routing_interface.f90

SUBROUTINE initatmgrid(rankmpi, nbcore, nbpt, nbseg, polygons_in, centre_in, area_in, contfrac_in, neighbours_in)
  !
  USE grid
  USE ioipsl
  !
  IMPLICIT NONE
  !
  !! INPUT VARIABLES
  INTEGER, INTENT(in) :: rankmpi   !! rank of processor
  INTEGER, INTENT(in) :: nbcore
  INTEGER, INTENT(in) :: nbpt   !! Atmospheric Domain size (unitless)
  INTEGER, INTENT(in) :: nbseg  !! Number of segments in the polygone
  REAL, INTENT(in) :: polygons_in(nbpt,2*nbseg+1,2)
  REAL, INTENT(in) :: centre_in(nbpt,2)
  REAL, INTENT(in) :: contfrac_in(nbpt)
  REAL, INTENT(in) :: area_in(nbpt)
  INTEGER, INTENT(in) :: neighbours_in(nbpt,nbseg*2)
  !
  !! Local
  INTEGER :: i
  INTEGER :: neighbours_loc(nbpt,nbseg*2)
  CHARACTER(LEN=33) :: outfname
  !
  DO i=1,nbpt
     ! Move from C to F indexing.
     neighbours_loc(i,:) = neighbours_in(i,:)+1
  ENDDO
  !
  CALL grid_init(nbpt, nbseg, polygons_in, centre_in, area_in, contfrac_in, neighbours_loc)
  !
  numout=100+rankmpi
  WRITE(outfname,'(A15,I4.4,A6,I4.4,A4)') "Out_RoutingReg_",rankmpi,"_from_",nbcore,".txt"
  OPEN(unit=numout, file=outfname)
  !
END SUBROUTINE initatmgrid

SUBROUTINE closeinterface()
  !
  USE ioipsl
  !
  CLOSE(unit=numout)
  !
END SUBROUTINE closeinterface

SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, maxpercent_in, sub_pts, sub_index, sub_area,  max_basins, &
     & min_topoind, meanrlen, meanrdz, minrdz, lon_rel, lat_rel, trip, basins, topoindex, rlen, rdz, fac, hierarchy, &
     & orog, floodp, nbi, nbj, area_bx, trip_bx, basin_bx, topoind_bx, rlen_bx, rdz_bx, rweight_bx, fac_bx, &
     & hierarchy_bx, orog_bx, floodp_bx, lon_bx, lat_bx, lshead_bx)
  !
  USE ioipsl
  USE grid
  USE routing_tools
  USE routing_reg
  !
  IMPLICIT NONE
  !
  INTEGER, INTENT(in) :: nbvmax_in, ijdimmax, maxpercent_in
  !
  INTEGER, INTENT(in) :: nbpt !! Domain size (unitless)
  INTEGER, INTENT(in) :: sub_pts(nbpt) !! Number of high resolution points on this grid (unitless)
  INTEGER, INTENT(in) :: sub_index(nbpt,nbvmax_in,2) !! Indices of the points we need on the fine grid (unitless)
  REAL, INTENT(inout) :: max_basins !! The current maximum of basins
  REAL, INTENT(in) :: min_topoind !! The current minimum of topographic index (km)
  REAL, INTENT(in) :: meanrlen, meanrdz, minrdz !! Mean length and altitude change in the original grid (m)
  REAL, INTENT(in) :: sub_area(nbpt,nbvmax_in) !! Area on the fine grid (m^2)
  !
  REAL, INTENT(in) :: lon_rel(nbpt,nbvmax_in) !!
  REAL, INTENT(in) :: lat_rel(nbpt,nbvmax_in) !! coordinates of the fine grid
  !
  REAL, INTENT(inout) :: trip(nbpt,nbvmax_in) !! The trip field (unitless)
  REAL, INTENT(inout) :: basins(nbpt,nbvmax_in) !! data on the fine grid
  REAL, INTENT(inout) :: topoindex(nbpt,nbvmax_in) !! Topographic index of the residence time (m)
  REAL, INTENT(inout) :: rlen(nbpt,nbvmax_in) !!
  REAL, INTENT(inout) :: rdz(nbpt,nbvmax_in) !!
  REAL, INTENT(inout) :: hierarchy(nbpt,nbvmax_in) !! data on the fine grid
  REAL, INTENT(inout) :: fac(nbpt,nbvmax_in) !! data on the fine grid
  REAL, INTENT(inout) :: orog(nbpt,nbvmax_in) !! data on the fine grid
  REAL, INTENT(inout) :: floodp(nbpt,nbvmax_in) !! data on the fine grid
  !
  !! OUTPUT VARIABLES
  INTEGER, INTENT(out) :: nbi(nbpt), nbj(nbpt) !! Number of point in x and y within the grid (unitless)
  REAL, INTENT(out) :: area_bx(nbpt,ijdimmax,ijdimmax) !! Area of each small box in the grid box (m^2)
  REAL, INTENT(out) :: hierarchy_bx(nbpt,ijdimmax,ijdimmax) !! Level in the basin of the point
  REAL, INTENT(out) :: fac_bx(nbpt,ijdimmax,ijdimmax) !! Flow accumulation
  REAL, INTENT(out) :: lon_bx(nbpt,ijdimmax,ijdimmax) !!
  REAL, INTENT(out) :: lat_bx(nbpt,ijdimmax,ijdimmax) !!
  REAL, INTENT(out) :: lshead_bx(nbpt,ijdimmax,ijdimmax) !! Large scale heading for outflow points.
  REAL, INTENT(out) :: topoind_bx(nbpt,ijdimmax,ijdimmax) !! Topographic index of the residence time for each of the smaller boxes (m)
  REAL, INTENT(out) :: rlen_bx(nbpt,ijdimmax,ijdimmax) !!
  REAL, INTENT(out) :: rdz_bx(nbpt,ijdimmax,ijdimmax) !!
  REAL, INTENT(out) :: rweight_bx(nbpt,ijdimmax,ijdimmax) !!
  INTEGER, INTENT(out) :: trip_bx(nbpt,ijdimmax,ijdimmax) !! The trip field for each of the smaller boxes (unitless)
  INTEGER, INTENT(out) :: basin_bx(nbpt,ijdimmax,ijdimmax) !!
  REAL, INTENT(out) :: orog_bx(nbpt,ijdimmax,ijdimmax) !!  Orography (m)
  REAL, INTENT(out) :: floodp_bx(nbpt,ijdimmax,ijdimmax) !! floodplains (m^2)
  !
  INTEGER :: ii, ib
  REAL :: resolution(nbpt,2)
  LOGICAL        :: debug = .FALSE.
  !
  ! nbvmax is still used to dimension wome variables in routing_reg.f90.
  ! It is transfered here but should be argument to the various subroutines.
  !
  nbvmax = nbvmax_in
  maxpercent = maxpercent_in
  !
  IF ( debug ) WRITE(numout,*) "Memory Mgt getgrid : nbvmax, ijdimmax = ", nbvmax, ijdimmax
  !
  DO ii=1,nbpt
     resolution(ii,1) = SQRT(area(ii))
     resolution(ii,2) = SQRT(area(ii))
  ENDDO

  DO ib=1,nbpt
     CALL routing_reg_getgrid(nbpt, ib, ijdimmax, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
          & meanrlen, meanrdz, minrdz, lon_rel, lat_rel, lalo, resolution, contfrac, trip, basins, topoindex, &
          & rlen, rdz, fac, hierarchy, orog, floodp, nbi(ib), nbj(ib), area_bx(ib,:,:), trip_bx(ib,:,:), &
          & basin_bx(ib,:,:), topoind_bx(ib,:,:), rlen_bx(ib,:,:), rdz_bx(ib,:,:), fac_bx(ib,:,:), &
          & hierarchy_bx(ib,:,:), orog_bx(ib,:,:), floodp_bx(ib,:,:), lon_bx(ib,:,:), lat_bx(ib,:,:), &
          & lshead_bx(ib,:,:))
  ENDDO
  !
  rweight_bx(:,:,:) = zero
  !
END SUBROUTINE gethydrogrid


SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, trip_bx, basin_bx, fac_bx, hierarchy_bx, topoind_bx, &
     &                rlen_bx, rdz_bx, rweight_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, orog_bx, &
     & basin_topoindex_stream, basin_rlen_stream, basin_dz_stream)
  !
  USE ioipsl
  USE grid
  USE routing_tools
  USE routing_reg
  !
  ! Arguments
  !
  INTEGER, INTENT(in)    :: nb_htu, nbv, ijdimmax
  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,ijdimmax,ijdimmax) !! The trip field for each of the smaller boxes (unitless)
  INTEGER, INTENT(inout) :: basin_bx(nbpt,ijdimmax,ijdimmax) !!
  REAL, INTENT(in)       :: fac_bx(nbpt,ijdimmax,ijdimmax) !! Flow accumulation
  REAL, INTENT(in)       :: hierarchy_bx(nbpt,ijdimmax,ijdimmax) !! Level in the basin of the point
  REAL, INTENT(inout)    :: topoind_bx(nbpt,ijdimmax,ijdimmax) !! Topographic index of the residence time for each of the smaller boxes (m)
  REAL, INTENT(inout)    :: rlen_bx(nbpt,ijdimmax,ijdimmax) !!
  REAL, INTENT(inout)    :: rdz_bx(nbpt,ijdimmax,ijdimmax) !!
  REAL, INTENT(inout)    :: rweight_bx(nbpt,ijdimmax,ijdimmax) !! Weight of each river within the HTU
  REAL, INTENT(in)       :: lshead_bx(nbpt,ijdimmax,ijdimmax) !! Large scale heading for outflow points.
  REAL, INTENT(in)       :: orog_bx(nbpt,ijdimmax,ijdimmax) !! 
  !
  !
  ! OUTPUT VARIABLES
  INTEGER, INTENT(out) :: nb_basin(nbpt)              !! Number of sub-basins (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)
  REAL, INTENT(out)    :: basin_topoindex_stream(nbpt,nb_htu)   !!
  REAL, INTENT(out)    :: basin_rlen_stream(nbpt,nb_htu)   !!
  REAL, INTENT(out)    :: basin_dz_stream(nbpt,nb_htu)   !!
  !
  ! For debug and get coordinate of river outlet
  !
  REAL, INTENT(in) :: lontmp(nbpt,ijdimmax,ijdimmax) !! Longitude
  REAL, INTENT(in) :: lattmp(nbpt,ijdimmax,ijdimmax) !! Latitude
  !
  ! Local
  !
  INTEGER :: ib
  LOGICAL        :: debug = .FALSE.
  !
  !diaglalo(1,:) = (/ 39.6791, 2.6687 /)
  !
  IF ( debug) WRITE(numout,*) "Memory Mgt findbasin : nbvmax, nb_htu, nbv = ", nbvmax, nb_htu, nbv

  DO ib=1,nbpt
     CALL routing_reg_findbasins(nb_htu, nbv, ib, ijdimmax, nbi(ib), nbj(ib), trip_bx(ib,:,:), &
          & basin_bx(ib,:,:), fac_bx(ib,:,:), hierarchy_bx(ib,:,:), topoind_bx(ib,:,:), &
          & rlen_bx(ib,:,:), rdz_bx(ib,:,:), rweight_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,:,:),&
          & orog_bx(ib,:,:), basin_topoindex_stream(ib,:), basin_rlen_stream(ib,:), basin_dz_stream(ib,:))
  ENDDO
  
END SUBROUTINE findbasins

SUBROUTINE globalize(nbpt, nb_htu, nbv, ijdimmax, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_bx, orog_bx, floodp_bx, &
          & fac_bx, topoind_bx, rlen_bx, rdz_bx, rweight_bx, &
          & min_topoind, nb_basin, basin_inbxid, basin_outlet, basin_outtp, basin_sz, basin_pts, basin_bxout, &
          & basin_bbout, lshead, coast_pts, nwbas, basin_rlen_stream, basin_dz_stream, basin_count, &
          & basin_notrun, basin_area, basin_cg, basin_hierarchy, &
          & basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_fac, basin_topoind, &
          & basin_rlen, basin_rdz, basin_id, basin_coor, basin_type, basin_flowdir, basin_lshead, &
          & basin_beta_fp, outflow_grid, outflow_basin, nbcoastal, coastal_basin)
  !
  USE ioipsl
  USE grid
  USE routing_tools
  USE routing_reg
  !
  ! Arguments
  !
  !! INPUT VARIABLES
  INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
  INTEGER(i_std), INTENT(in)  :: nbv, nb_htu, ijdimmax
  REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax)    :: area_bx      !! Area of each small box in the grid box (m^2)
  REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax)    :: lon_bx       !! Longitude of each small box in the grid box
  REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax)    :: lat_bx       !! Latitude of each small box in the grid box
  INTEGER(i_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax) :: trip_bx      !! The trip field for each of the smaller boxes (unitless)
  REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax)    :: hierarchy_bx !! Level in the basin of the point
  REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax)    :: orog_bx      !! Orography (m)
  REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax)    :: floodp_bx    !! Floodplains area (m^2)
  REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax)    :: fac_bx       !! Level in the basin of the point
  REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax)    :: topoind_bx   !! Topographic index of the residence time for each of the smaller boxes (km)
  REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax)    :: rlen_bx      !!
  REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax)    :: rdz_bx       !!
  REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax)    :: rweight_bx   !!
  REAL(r_std), INTENT (in)                                       :: min_topoind  !! The current minimum of topographic index (km)
  INTEGER(i_std), INTENT (in), DIMENSION(nbpt)                   :: nb_basin     !! Number of sub-basins (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)
  REAL(r_std), INTENT (in), DIMENSION(nbpt,nb_htu)        :: basin_rlen_stream, basin_dz_stream
  !! Output VARIABLES
  INTEGER(i_std), INTENT (in) :: nwbas !!
  INTEGER(i_std), INTENT (out), DIMENSION(nbpt)       :: basin_count         !!
  INTEGER(i_std), INTENT (out), DIMENSION(nbpt)       :: basin_notrun        !!
  INTEGER(i_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_id      !!
  REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas,2)  :: basin_coor     !!
  REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas)    :: basin_type       !!
  INTEGER(i_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_flowdir !! Water flow directions in the basin (unitless)
  REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas)    :: basin_area       !!
  REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas,2)  :: basin_cg         !! Centre of gravity of the HTU
  REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas)    :: basin_hierarchy  !!
  REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas)    :: basin_orog_mean  !!
  REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas)    :: basin_orog_min  !!
  REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas)    :: basin_orog_max  !!
  REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas)    :: basin_floodp  !!
  REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas)    :: basin_fac  !!
  REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas)    :: basin_topoind    !! Topographic index of the residence time for a basin (km)
  REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas)    :: basin_rlen
  REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas)    :: basin_rdz
  REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas)    :: basin_lshead     !! Large scale heading out of the grid box.
  REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas)    :: basin_beta_fp  !! Beta parameter for floodplains
  INTEGER(i_std), INTENT (out), DIMENSION(nbpt,nwbas) :: outflow_grid  !! Type of outflow on the grid box (unitless)
  INTEGER(i_std), INTENT (out), DIMENSION(nbpt,nwbas) :: outflow_basin !!
  INTEGER(i_std), INTENT (out), DIMENSION(nbpt)       :: nbcoastal           !!
  INTEGER(i_std), INTENT (out), DIMENSION(nbpt,nwbas) :: coastal_basin !!
  !!
  INTEGER(i_std) :: ib
  !!
  DO ib=1,nbpt
     CALL routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, 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,:,:), rlen_bx(ib,:,:), rdz_bx(ib,:,:), rweight_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, basin_count, basin_notrun, basin_area, basin_cg, basin_hierarchy,&
          & basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_fac,&
          & basin_topoind, basin_rlen, basin_rdz, basin_id, basin_coor, basin_type, &
          & basin_flowdir, basin_lshead, basin_beta_fp, outflow_grid, outflow_basin, nbcoastal, coastal_basin)
     ! Once all possible operation on topoindex have been performed, we keep
     ! rlen_stream and rdz_stream as rlen and rdz
     basin_rlen(ib,:) = basin_rlen_stream(ib,:nwbas)
     basin_rdz(ib,:) = basin_dz_stream(ib,:nwbas)
  ENDDO
  !
END SUBROUTINE globalize

SUBROUTINE linkup(nbpt, ijdimmax, nwbas, inflowmax, basin_count, basin_area, basin_id, basin_flowdir, &
     & basin_lshead, basin_hierarchy, basin_fac, outflow_grid, outflow_basin, inflow_number, inflow_grid, &
     & inflow_basin, nbcoastal, coastal_basin, invented_basins)
  !
  USE ioipsl
  USE grid
  USE routing_tools
  USE routing_reg
  !
  ! Arguments
  !
  !! INPUT VARIABLES
  INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
  INTEGER(i_std), INTENT (in) :: ijdimmax, inflowmax
  REAL(r_std), INTENT(in) :: invented_basins !!
  !
  INTEGER(i_std), INTENT (in) :: nwbas !!
  INTEGER(i_std), INTENT (in), DIMENSION(nbpt) :: basin_count !!
  INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nwbas) :: basin_id !!
  INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nwbas) :: basin_flowdir !! Water flow directions in the basin (unitless)
  REAL(r_std), INTENT (in), DIMENSION(nbpt,nwbas) :: basin_lshead  !! Large scale flow direction out of the basin.
  REAL(r_std), INTENT (in), DIMENSION(nbpt,nwbas) :: basin_area !!
  REAL(r_std), INTENT (inout), DIMENSION(nbpt,nwbas) :: basin_hierarchy !!
  REAL(r_std), INTENT (in), DIMENSION(nbpt,nwbas) :: basin_fac !!
  INTEGER(i_std), INTENT(inout), DIMENSION(nbpt,nwbas) :: outflow_grid !! Type of outflow on the grid box (unitless)
  INTEGER(i_std), INTENT(inout), DIMENSION(nbpt,nwbas) :: outflow_basin !!
  !
  INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nwbas) :: inflow_number !!
  INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nwbas,inflowmax) :: inflow_basin !!
  INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nwbas,inflowmax) :: inflow_grid !!
  !
  INTEGER(i_std), INTENT(inout), DIMENSION(nbpt) :: nbcoastal !!
  INTEGER(i_std), INTENT(inout), DIMENSION(nbpt,nwbas) :: coastal_basin !!
  !
  !
  LOGICAL        :: debug = .FALSE.
  !
  IF ( debug ) WRITE(numout,*) "Memory Mgt Linkup : nbvmax, ijdimmax, nwbas, inflowmax = ", nbvmax, ijdimmax, nwbas, inflowmax

  CALL routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basin_count, basin_area, basin_id, basin_flowdir, &
       & basin_lshead, basin_hierarchy, basin_fac, diaglalo, outflow_grid, outflow_basin, inflow_number, inflow_grid, &
       & inflow_basin, nbcoastal, coastal_basin, invented_basins)
  !
END SUBROUTINE linkup

SUBROUTINE areanorm(nbpt, nwbas, basin_count, basin_area, outflow_grid, basin_area_norm)
  !
  USE grid
  !
  !! INPUT VARIABLES
  INTEGER(i_std), INTENT (in)                     :: nbpt !! Domain size (unitless)
  INTEGER(i_std), INTENT (in)                     :: nwbas !!
  INTEGER(i_std), DIMENSION(nbpt), INTENT(in)     :: basin_count !!
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in)  :: basin_area !!
  !
  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, ig
  REAL(r_std) :: contarea !!
  REAL(r_std) :: totbasins !!
  !
  ! Normalize the area of all basins
  !
  DO ig=1,nbpt
     !
     totbasins = SUM(basin_area(ig,1:basin_count(ig)))
     contarea = area(ig)*contfrac(ig)
     !
     DO ib=1,basin_count(ig)
        basin_area_norm(ig,ib) = basin_area(ig,ib)/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(ig,ib) .EQ. -1) THEN
           outflow_grid(ig,ib) = -2
        ENDIF

     ENDDO
     !
  ENDDO
  !
END SUBROUTINE areanorm

SUBROUTINE fetch(nbpt, nwbas, nbcore, nbhalo, fhalopts, corepts, basin_count, basin_area, basin_id, basin_hierarchy, basin_fac, &
     & outflow_grid, outflow_basin, partial_sum, fetch_basin, outflow_uparea)
  !
  USE ioipsl
  USE grid
  USE routing_tools
  USE routing_reg
  !
  ! Arguments
  !
  !! 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), INTENT (in) :: nbhalo
  INTEGER(i_std), DIMENSION(nbhalo), INTENT(in)        :: fhalopts
  INTEGER(i_std), DIMENSION(nbcore), INTENT(in)        :: corepts
  INTEGER(i_std), DIMENSION(nbpt), INTENT(in)          :: basin_count !!
  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
  !
  !! IN-OUTPUT VARIABLES
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: fetch_basin !!
  ! Output
  REAL(r_std), DIMENSION(nbpt*nwbas), INTENT(out) :: 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,nbhalo, fhalopts, &
       & fcorepts, basin_count, basin_area, basin_id, &
       & basin_hierarchy, basin_fac, outflow_grid, outflow_basin, partial_sum, fetch_basin, outflow_uparea)

  !
END SUBROUTINE fetch

SUBROUTINE rivclassification(nbpt, nwbas, nbcore, corepts, basin_count, outflow_grid, outflow_basin, fetch_basin, &
     & largest_rivarea, num_largest)
  !
  USE defprec
  !
  !! 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 !!
  INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid
  INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_basin !!
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in)       :: fetch_basin !!
  REAL(r_std), INTENT(in) :: largest_rivarea
  INTEGER(i_std), INTENT (out) :: num_largest
  !
  !! LOCAL VARIABLES
  INTEGER(i_std) :: i, ib, ig
  INTEGER(i_std), DIMENSION(nbcore)               :: fcorepts
  !
  !
  DO i=1,nbcore
     fcorepts(i) = corepts(i)+1
  ENDDO
  !
  num_largest = 0
  !
  ! Just a dummy use to avoid f2py warnings
  outflow_basin(1,1) = outflow_basin(1,1)
  !
  DO i=1,nbcore
     ig = fcorepts(i)
     !
     DO ib=1,basin_count(ig)
        !
        IF (outflow_grid(ig,ib) .LT. 0 .AND. outflow_grid(ig,ib) .GT. -3  .AND. fetch_basin(ig,ib) .GE. largest_rivarea) THEN
           num_largest = num_largest + 1
           outflow_grid(ig,ib) = -1
        ENDIF
        !
     ENDDO
  ENDDO
  !
END SUBROUTINE rivclassification


SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridarea, cfrac, gridcenters, basin_count, &
     & basin_notrun, basin_area, basin_orog_mean, basin_orog_min, basin_orog_max, &
     & basin_floodp, basin_cg, basin_topoind, basin_topoindex_stream, basin_rlen, &
     & basin_rdz, basin_beta_fp, fetch_basin, &
     & basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, &
     & inflow_number, inflow_grid, inflow_basin, floodcri, routing_area, routing_orog_mean, &
     & routing_orog_min, routing_orog_max, routing_floodp, routing_beta, &
     & routing_cg, topo_resid, topo_resid_stream, topo_rlen, topo_rdz, & 
     & route_nbbasin, route_togrid, route_tobasin, route_nbintobas, &
     & global_basinid, route_outlet, route_type, origin_nbintobas, routing_fetch,routing_floodcri, &
     & rfillval, ifillval)
  !
  USE ioipsl
  USE grid
  USE routing_tools
  USE routing_reg
  !
  ! Arguments
  !
  !! INPUT VARIABLES
  INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
  INTEGER(i_std), INTENT (in) :: inflowmax, nbasmax
  INTEGER(i_std), INTENT (in) :: nwbas !!
  INTEGER(i_std), INTENT (in) :: num_largest
  REAL(r_std), DIMENSION(nbpt), INTENT(in)          :: gridarea
  REAL(r_std), DIMENSION(nbpt), INTENT(in)          :: cfrac
  REAL(r_std), DIMENSION(nbpt,2), INTENT(in)          :: gridcenters
  !
  INTEGER(i_std), DIMENSION(nbpt), INTENT(inout)       :: basin_count !!
  INTEGER(i_std), DIMENSION(nbpt), INTENT(inout)       :: basin_notrun !!
  INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_id !!
  REAL(r_std), DIMENSION(nbpt,nwbas,2), INTENT(inout)  :: basin_coor !!
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_type !!
  INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_flowdir !! Water flow directions in the basin (unitless)
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_area !!
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_orog_mean !!
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_orog_min !!
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_orog_max !!
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_floodp !!
  REAL(r_std), DIMENSION(nbpt,nwbas,2), INTENT(inout)  :: basin_cg
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_topoind !! Topographic index of the residence time for a basin (km)
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_topoindex_stream !! Topographic index of the residence time for a basin (km) - stream res.
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_rlen
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_rdz
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_beta_fp
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: fetch_basin !!
  INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless)
  INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_basin !!
  !
  REAL(r_std), INTENT(in)    :: rfillval
  INTEGER(i_std), INTENT(in) :: ifillval
  !
  !
  !
  INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout)          :: inflow_number !!
  INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(inout) :: inflow_basin !!
  INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(inout) :: inflow_grid !!
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: floodcri !!
  !
  ! Output variables
  !
  REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: routing_area !! Surface of basin (m^2)
  REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: routing_orog_mean !! Mean Orography (m)
  REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: routing_orog_min !! Mean Orography (m)
  REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: routing_orog_max !! Mean Orography (m)
  REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: routing_floodp !! Surface of Floodplains (m^2)
  REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: routing_beta
  REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: routing_fetch !! Upstream are flowing into HTU (m^2)
  REAL(r_std), DIMENSION(nbpt,nbasmax,2), INTENT(out)  :: routing_cg !! Centre of gravity of HTU (Lat, Lon)
  REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: topo_resid !! Topographic index of the retention time (km)
  REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: topo_resid_stream !! Topographic index of the retention time (km)
  REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: topo_rlen !!
  REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: topo_rdz  !!
  INTEGER(i_std), DIMENSION(nbpt), INTENT(out)         :: route_nbbasin !!
  INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: route_togrid !! Grid into which the basin flows (unitless)
  INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: route_tobasin !! Basin in to which the water goes (unitless)
  INTEGER(i_std), DIMENSION(nbpt), INTENT(out)         :: route_nbintobas !! Number of basin into current one (unitless)
  INTEGER(i_std), DIMENSION(nbpt), INTENT(out)         :: origin_nbintobas !! Number of sub-grid basin into current one before truncation (unitless)
  INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: global_basinid !! ID of basin (unitless)
  REAL(r_std), DIMENSION(nbpt,nbasmax,2), INTENT(out)  :: route_outlet !! Coordinate of outlet (-)
  REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: route_type !! Coordinate of outlet (-)
  REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: routing_floodcri
  !
  !
  INTEGER(i_std) :: ib, ij
  !
  CALL routing_reg_end_truncate(nbpt, nbasmax, gridarea, cfrac, gridcenters, nwbas, inflowmax, num_largest, &
       & basin_count, basin_notrun, basin_area, basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp,&
       & basin_cg, basin_topoind, basin_topoindex_stream, basin_rlen, basin_rdz, & 
       & basin_beta_fp, fetch_basin, basin_id, basin_coor, &
       & basin_type, basin_flowdir, outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin, floodcri)

  routing_area(:,:) = rfillval
  routing_orog_mean(:,:) = rfillval
  routing_orog_min(:,:) = rfillval
  routing_orog_max(:,:) = rfillval
  routing_floodp(:,:) = rfillval
  routing_cg(:,:,:) = rfillval
  topo_resid(:,:) = rfillval
  topo_resid_stream(:,:) = rfillval
  topo_rlen(:,:) = rfillval
  topo_rdz(:,:) = rfillval
  route_nbbasin(:) = ifillval
  route_togrid(:,:) = ifillval
  route_tobasin(:,:) = ifillval
  routing_fetch(:,:) = rfillval
  routing_beta(:,:) = rfillval
  route_nbintobas(:) = ifillval
  global_basinid(:,:) = ifillval
  route_outlet(:,:,:) = rfillval
  route_type(:,:) = rfillval
  origin_nbintobas(:) = ifillval

  DO ij=1,nbpt
     DO ib=1,basin_count(ij)
        !
        routing_area(ij,ib) = routing_area_glo(ij,ib)
        routing_orog_mean(ij,ib) = routing_orog_mean_glo(ij,ib)
        routing_orog_min(ij,ib) = routing_orog_min_glo(ij,ib)
        routing_orog_max(ij,ib) = routing_orog_max_glo(ij,ib)
        routing_floodp(ij,ib) = routing_floodp_glo(ij,ib)
        routing_cg(ij,ib,:) = routing_cg_glo(ij,ib,:)
        !
        topo_resid(ij,ib) = topo_resid_glo(ij,ib)
        topo_resid_stream(ij,ib) = topo_resid_stream_glo(ij,ib)
        topo_rlen(ij,ib) = topo_rlen_glo(ij,ib)
        topo_rdz(ij,ib) = topo_rdz_glo(ij,ib)
        !
        route_nbbasin(:) = route_count_glo(:)
        route_togrid(ij,ib) = route_togrid_glo(ij,ib)
        route_tobasin(ij,ib) = route_tobasin_glo(ij,ib)
        routing_floodcri(ij,ib) = routing_floodcri_glo(ij,ib)
        routing_fetch(ij,ib) = route_fetch_glo(ij,ib)
        routing_beta(ij,ib) = routing_beta_glo(ij,ib)
        !
        route_nbintobas(ij) = route_nbintobas_glo(ij)
        global_basinid(ij,ib) = global_basinid_glo(ij,ib)
        route_outlet(ij,ib,:) = route_outlet_glo(ij,ib,:)
        route_type(ij,ib) = route_type_glo(ij,ib)
        origin_nbintobas(ij) = origin_nbintobas_glo(ij)
     ENDDO
  ENDDO
  !
END SUBROUTINE finish_truncate


SUBROUTINE finish_inflows(nbpt, nwbas, nbasmax, inf_max, basin_count, inflow_number, inflow_grid, inflow_basin,&
       & route_innum, route_ingrid, route_inbasin)
  !
  USE ioipsl
  USE grid
  USE routing_tools
  USE routing_reg
  !
  ! Arguments
  !
  INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
  INTEGER(i_std), INTENT (in) :: nwbas, nbasmax
  INTEGER(i_std), INTENT (in) :: inf_max !!
  !
  INTEGER(i_std), DIMENSION(nbpt), INTENT(inout)             :: basin_count !!
  INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in)          :: inflow_number !!
  INTEGER(i_std), DIMENSION(nbpt,nwbas,inf_max), INTENT(in)  :: inflow_basin !!
  INTEGER(i_std), DIMENSION(nbpt,nwbas,inf_max), INTENT(in)  :: inflow_grid !!
  !
  REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)            :: route_innum
  REAL(r_std), DIMENSION(nbpt,nbasmax,inf_max), INTENT(out)    :: route_ingrid
  REAL(r_std), DIMENSION(nbpt,nbasmax,inf_max), INTENT(out)    :: route_inbasin
  !
  !
  !
  ! inflow_number -> route_innum
  route_innum(:,:) = 0
  DO ig=1,nbpt
    nbas = basin_count(ig)
    route_innum(ig,:nbas) = inflow_number(ig,:nbas)
  END DO

  ! inflow_grid -> route_ingrid
  ! inflow_basin -> route_inbasin
  route_ingrid(:,:,:) = 0
  route_inbasin(:,:,:) = 0
  DO ig=1,nbpt
    nbas = basin_count(ig)
    DO ib=1,nbas
      nin = route_innum(ig,ib)
      route_ingrid(ig,ib,:nin) = inflow_grid(ig,ib,:nin)
      route_inbasin(ig,ib,:nin) = inflow_basin(ig,ib,:nin)
    END DO
  END DO

END SUBROUTINE finish_inflows


SUBROUTINE killbas(nbpt, inflowmax, nbasmax, nwbas, ops, tokill, totakeover, numops, basin_count, basin_area, &
     & basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_cg, &
     & basin_topoind, basin_topoindex_stream, basin_rlen, basin_rdz, basin_beta_fp, &
     & fetch_basin, basin_id, basin_coor, basin_type, &
     & basin_flowdir, outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin)
  !
  !
  USE ioipsl
  USE grid
  USE routing_tools
  USE routing_reg
  !
  !! INPUT VARIABLES
  INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
  INTEGER(i_std), INTENT (in) :: inflowmax, nbasmax
  INTEGER(i_std), INTENT (in) :: nwbas !!
  INTEGER(i_std), INTENT (in) :: ops !!
  INTEGER(i_std), DIMENSION(nbpt, ops), INTENT (in) :: tokill, totakeover
  INTEGER(i_std), DIMENSION(nbpt), INTENT (in) :: numops
  !
  INTEGER(i_std), DIMENSION(nbpt), INTENT(inout)       :: basin_count !!
  INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_id !!
  REAL(r_std), DIMENSION(nbpt,nwbas,2), INTENT(inout)  :: basin_coor !!
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_type !!
  INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_flowdir !! Water flow directions in the basin (unitless)
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_area !!
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_orog_mean !!
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_orog_min !!
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_orog_max !!

  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_floodp !!
  REAL(r_std), DIMENSION(nbpt,nwbas,2), INTENT(inout)  :: basin_cg
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_topoind !! Topographic index of the residence time for a basin (km)
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_topoindex_stream !! Topographic index of the residence time for a basin (km)
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_rlen
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_rdz
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_beta_fp
  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: fetch_basin !!
  INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless)
  INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_basin !!
  !
  !
  !
  INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout)          :: inflow_number !!
  INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(inout) :: inflow_basin !!
  INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(inout) :: inflow_grid !!

  !! LOCAL
  INTEGER(i_std) :: ib, op, tok, totak, igrif, ibasf
  LOGICAL        :: debug = .FALSE.
  !
  !
  !
  DO ib=1,nbpt
    DO op=1,numops(ib)
      IF (basin_count(ib) > nbasmax) THEN
         tok = tokill(ib,op)
         totak = totakeover(ib,op)
         IF (tok .GT. 0) THEN
            IF ( debug ) THEN
               WRITE(numout,*) " "
               WRITE(numout,*) "numop", op
               WRITE(numout,*) "nbpt", ib, "tokill", tok, "totakover", totak, "basin_count", basin_count(ib)
               WRITE(numout,*) "outflow tokill", outflow_grid(ib,tok), outflow_basin(ib,tok)
               WRITE(numout,*) "inflows totakeover", inflow_number(ib,totak)
               DO it=1,inflow_number(ib,totak)
                 WRITE(numout,*) "->", inflow_grid(ib,totak,it), inflow_basin(ib,totak,it)
               END DO
            ENDIF
            ! Test if tokill is downstream of totakeover (avoid loop)
            igrif = outflow_grid(ib,totak)
            ibasf = outflow_basin(ib,totak)
            DO WHILE ( igrif >0 )
              IF ((igrif .EQ. ib) .AND. (ibasf .EQ. tok)) THEN
                !CALL ipslerr_p(3,'killbas','tokill is downstream totakeover','','')
                igrif = 0
                it = totak
                totak = tok
                tok = it
              ELSE
                it = outflow_grid(igrif,ibasf)
                ibasf = outflow_basin(igrif,ibasf)
                igrif = it
              END IF
            END DO
            CALL routing_reg_killbas(nbpt, ib, tok, totak, nwbas, inflowmax, basin_count, basin_area, &
               & basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_cg,&
               & basin_topoind, basin_topoindex_stream, basin_rlen, &
               & basin_rdz, basin_beta_fp, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, &
               & outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin)
         END IF
      END IF
    END DO
  END DO

END SUBROUTINE killbas


SUBROUTINE checkrouting(nbpt, nwbas, outflow_grid, outflow_basin, basin_count)
  !
  USE ioipsl
  USE grid
  USE routing_tools
  USE routing_reg
  !
  !! INPUT VARIABLES
  INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
  INTEGER(i_std), INTENT (in) :: nwbas !!

  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 !!

  INTEGER(i_std), DIMENSION(nbpt),       INTENT(in)       :: basin_count !!
  ! Local
  INTEGER(i_std), DIMENSION(nbpt,nwbas)       :: flag !!
  INTEGER(i_std) :: ig, ib, it, igrif, ibasf, basnum, test
  LOGICAL        :: debug = .FALSE.
  !
  flag(:,:) = 0
  !
  IF ( debug ) WRITE(numout,*) "Checking routing"
  test = 0

  DO ig=1,nbpt
    basnum = basin_count(ig)
    DO ib=1,basnum
      IF (flag(ig,ib) .EQ. 0) THEN
        flag(ig,ib) = -1
        igrif = outflow_grid(ig,ib)
        ibasf = outflow_basin(ig,ib)
      !
        DO WHILE ((igrif .GT. 0) .AND. (flag(igrif,ibasf) .NE. -1) .AND. (flag(igrif,ibasf) .NE. -99))
          flag(igrif, ibasf) = -1
          it = outflow_grid(igrif,ibasf)
          ibasf = outflow_basin(igrif,ibasf)
          igrif = it
        END DO


        IF ((flag(igrif,ibasf) .EQ. -99) .OR. (igrif .LE. 0)) THEN
          flag(ig,ib) = -99
          igrif = outflow_grid(ig,ib)
          ibasf = outflow_basin(ig,ib)

          DO WHILE ((igrif .GT. 0)</