Commit 20b3f417 authored by POLCHER Jan's avatar POLCHER Jan 🚴🏾
Browse files

Now the result with 1 or 2 processors is about the same.

The fact that performing the sums in different order still induces differences but they are small.
parent 77cb5d70
......@@ -216,7 +216,7 @@ 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])
routing_area,routing_cg,topo_resid,route_nbbasin,route_togrid,route_tobasin,route_nbintobas,global_basinid,route_outlet,route_type,origin_nbintobas,routing_fetch = truncate(nbasmax,num_largest,corepts,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,nbcore])
Wrapper for ``truncate``.
......@@ -224,6 +224,7 @@ Parameters
----------
nbasmax : input int
num_largest : input int
corepts : input rank-1 array('i') with bounds (nbcore)
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)
......@@ -248,12 +249,15 @@ nbxmax_in : input int, optional
Default: (shape(inflow_number,1))/(8)
nwbas : input int, optional
Default: shape(basin_area,1)
nbcore : input int, optional
Default: len(corepts)
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_nbbasin : rank-1 array('i') with bounds (nbpt)
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)
......
......@@ -300,26 +300,26 @@ SUBROUTINE areanorm(nbpt, nwbas, basin_count, basin_area, outflow_grid, basin_ar
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(out) :: basin_area_norm !!
!
!! LOCAL VARIABLES
INTEGER(i_std) :: ib, ij
INTEGER(i_std) :: ib, ig
REAL(r_std) :: contarea !!
REAL(r_std) :: totbasins !!
!
! Normalize the area of all basins
!
DO ib=1,nbpt
DO ig=1,nbpt
!
totbasins = SUM(basin_area(ib,1:basin_count(ib)))
contarea = area(ib)*contfrac(ib)
totbasins = SUM(basin_area(ig,1:basin_count(ig)))
contarea = area(ig)*contfrac(ig)
!
DO ij=1,basin_count(ib)
basin_area_norm(ib,ij) = basin_area(ib,ij)/totbasins*contarea
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(ib,ij) .EQ. -1) THEN
outflow_grid(ib,ij) = -2
IF ( outflow_grid(ig,ib) .EQ. -1) THEN
outflow_grid(ig,ib) = -2
ENDIF
ENDDO
......@@ -373,13 +373,16 @@ SUBROUTINE fetch(nbpt, nwbas, nbcore, corepts, basin_count, basin_area, basin_id
!
END SUBROUTINE fetch
SUBROUTINE rivclassification(nbpt, nwbas, basin_count, outflow_grid, outflow_basin, fetch_basin, largest_rivarea, num_largest)
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(in) :: outflow_basin !!
......@@ -388,18 +391,24 @@ SUBROUTINE rivclassification(nbpt, nwbas, basin_count, outflow_grid, outflow_bas
INTEGER(i_std), INTENT (out) :: num_largest
!
!! LOCAL VARIABLES
INTEGER(i_std) :: ib, ij
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
!
DO ib=1,nbpt
DO i=1,nbcore
ig = fcorepts(i)
!
DO ij=1,basin_count(ib)
DO ib=1,basin_count(ig)
!
IF (outflow_grid(ib,ij) .LT. 0 .AND. fetch_basin(ib,ij) .GE. largest_rivarea) THEN
IF (outflow_grid(ig,ib) .LT. 0 .AND. fetch_basin(ig,ib) .GE. largest_rivarea) THEN
num_largest = num_largest + 1
outflow_grid(ib,ij) = -1
outflow_grid(ig,ib) = -1
ENDIF
!
ENDDO
......@@ -407,10 +416,10 @@ SUBROUTINE rivclassification(nbpt, nwbas, basin_count, outflow_grid, outflow_bas
!
END SUBROUTINE rivclassification
SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, 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, routing_area, routing_cg, topo_resid, route_togrid, route_tobasin, route_nbintobas, &
& global_basinid, route_outlet, route_type, origin_nbintobas, routing_fetch)
SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, nbcore, corepts, 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, routing_area, routing_cg, topo_resid, route_nbbasin, route_togrid, &
& route_tobasin, route_nbintobas, global_basinid, route_outlet, route_type, origin_nbintobas, routing_fetch)
!
USE ioipsl
USE grid
......@@ -424,6 +433,8 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, basin_count, b
INTEGER(i_std), INTENT (in) :: nbxmax_in, nbasmax
INTEGER(i_std), INTENT (in) :: nwbas !!
INTEGER(i_std), INTENT (in) :: num_largest
INTEGER(i_std), INTENT (in) :: nbcore
INTEGER(i_std), DIMENSION(nbcore), INTENT(in) :: corepts
!
INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_count !!
INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_notrun !!
......@@ -450,6 +461,7 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, basin_count, b
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 (m)
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)
......@@ -458,18 +470,26 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, basin_count, b
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 (-)
!
INTEGER(i_std) :: ic
INTEGER(i_std), DIMENSION(nbcore) :: fcorepts
!
IF ( nbxmax_in .NE. nbxmax ) THEN
WRITE(*,*) "TRUNCATE : nbxmax has changed !!"
STOP
ENDIF
CALL routing_reg_truncate(nbpt, nbasmax, area, contfrac, nwbas, 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)
!
DO ic=1,nbcore
fcorepts(ic) = corepts(ic)+1
ENDDO
CALL routing_reg_truncate(nbpt, nbasmax, area, contfrac, nwbas, num_largest, nbcore, fcorepts, 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)
routing_area(:,:) = routing_area_glo(:,:)
routing_cg(:,:,:) = routing_cg_glo(:,:,:)
topo_resid(:,:) = topo_resid_glo(:,:)
topo_resid(:,:) = topo_resid_glo(:,:)
route_nbbasin(:) = route_count_glo(:)
route_togrid(:,:) = route_togrid_glo(:,:)
route_tobasin(:,:) = route_tobasin_glo(:,:)
......@@ -482,3 +502,180 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, basin_count, b
origin_nbintobas(:) = origin_nbintobas_glo(:)
END SUBROUTINE truncate
SUBROUTINE finalfetch(nbpt, nbasmax, nbcore, corepts, routing_area, basin_count, route_togrid, route_tobasin, partial_sum, &
& fetch_basin, outflow_uparea)
!
USE defprec
USE ioipsl
!
!! INPUT VARIABLES
INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
INTEGER(i_std), INTENT(in) :: nbasmax !!
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,nbasmax), INTENT(in) :: routing_area !! Surface of basin (m^2)
INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(in) :: route_togrid !! Grid into which the basin flows (unitless)
INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(in) :: route_tobasin !! Basin in to which the water goes (unitless)
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(in) :: partial_sum
!
!! INPUT/OUTPUT VARIABLES
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(inout) :: fetch_basin !!
! Out
REAL(r_std), DIMENSION(nbpt*nbasmax), INTENT(out) :: outflow_uparea
!
! LOCAL
INTEGER(i_std) :: ic, ib, ig, igrif, ibasf, itt, ff(1)
INTEGER(i_std), PARAMETER :: maxiterations=10000
INTEGER(i_std), DIMENSION(nbcore) :: fcorepts
!
!
DO ic=1,nbcore
fcorepts(ic) = corepts(ic)+1
ENDDO
!
! Propagate first the partial sums from the halo into the domain
!
DO ig=1,nbpt
!
DO ib=1,basin_count(ig)
IF (partial_sum(ig,ib) > EPSILON(partial_sum)) THEN
!
fetch_basin(ig,ib) = MAX(fetch_basin(ig,ib), partial_sum(ig,ib))
!
igrif = route_togrid(ig,ib)
ibasf = route_tobasin(ig,ib)
!
itt = 0
!
DO WHILE ( ibasf .LE. nbasmax )
!
fetch_basin(igrif,ibasf) = MAX(fetch_basin(igrif,ibasf), partial_sum(ig,ib))
!
it = route_togrid(igrif,ibasf)
ibasf = route_tobasin(igrif,ibasf)
igrif = it
itt = itt + 1
IF ( itt .GT. maxiterations) THEN
IF ( itt .GT. maxiterations+20) THEN
CALL ipslerr_p(3,'routing_interface','Problem in propagating partial sum','','')
ENDIF
ENDIF
ENDDO
ENDIF
ENDDO
ENDDO
!
! Add the areas contributed by the core region of the domain.
!
DO ic=1,nbcore
ig = fcorepts(ic)
!
DO ib=1,basin_count(ig)
!
fetch_basin(ig,ib) = fetch_basin(ig,ib) + routing_area(ig,ib)
!
igrif = route_togrid(ig,ib)
ibasf = route_tobasin(ig,ib)
!
itt = 0
!
DO WHILE ( ibasf .LE. nbasmax )
!
fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + routing_area(ig,ib)
!
it = route_togrid(igrif,ibasf)
ibasf = route_tobasin(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)") ig, ib, itt
WRITE(numout, "('==> We flow through grid ',I5,' and basin ',I5)") igrif, ibasf
!
IF ( itt .GT. maxiterations+20) THEN
CALL ipslerr_p(3,'routing_reg_fetch','Problem with computing final fetch','','')
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
!
! Get the largest upstream areas in this core
!
nboutflow = 0
outflow_uparea(:) = zero
!
DO ic=1,nbcore
ig = fcorepts(ic)
!
DO ib=1,basin_count(ig)
!
! We do not need any more the river flow flag as we are going to reset it.
! We archive the upstream area of all outflow points so that we can sort them.
!
IF ( route_tobasin(ig,ib) .GT. nbasmax) THEN
nboutflow = nboutflow + 1
IF ( nboutflow <= nbpt*nbasmax) THEN
outflow_uparea(nboutflow) = fetch_basin(ig,ib)
ELSE
! Not enought place so we replace a smaller one.
IF (fetch_basin(ig,ib) > MINVAL(outflow_uparea)) THEN
ff = MINLOC(outflow_uparea)
outflow_uparea(ff(1)) = fetch_basin(ig,ib)
ELSE
! Ignore basin
ENDIF
ENDIF
ENDIF
!
ENDDO
ENDDO
!
END SUBROUTINE finalfetch
SUBROUTINE finalrivclass(nbpt, nbasmax, nbcore, corepts, route_togrid, route_tobasin, routing_fetch, largest_rivarea, num_largest)
!
USE routing_reg
!
INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
INTEGER(i_std), INTENT (in) :: nbasmax !!
INTEGER(i_std), INTENT (in) :: nbcore
INTEGER(i_std), DIMENSION(nbcore), INTENT(in) :: corepts
INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(inout) :: route_togrid !! Grid into which the basin flows (unitless)
INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(inout) :: route_tobasin !! Basin in to which the water goes (unitless)
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(in) :: routing_fetch !! Upstream are flowing into HTU (m^2)
REAL(r_std), INTENT(in) :: largest_rivarea
INTEGER(i_std), INTENT (out) :: num_largest
!
! LOCAL
INTEGER(i_std) :: i, ib, ig
INTEGER(i_std), DIMENSION(nbcore) :: fcorepts
!
!
DO i=1,nbcore
fcorepts(i) = corepts(i)+1
ENDDO
!
! Redo the the distinction between river outflow and coastal flow. We can not
! take into account the return flow points.
!
num_largest = 0
DO i = 1, nbcore
ig = fcorepts(i)
DO ib = 1, nbasmax
IF (route_tobasin_glo(ig,ib) .GT. nbasmax) THEN
IF (routing_fetch(ig,ib) .GE. largest_rivarea ) THEN
num_largest = num_largest + 1
route_tobasin_glo(ig,ib) = nbasmax + 3
route_tobasin(ig,ib) = nbasmax + 3
ENDIF
ENDIF
ENDDO
ENDDO
!
WRITE(numout,*) 'NUMBER OF RIVERS :', COUNT(route_tobasin_glo .GE. nbasmax + 3)
!
END SUBROUTINE finalrivclass
......@@ -14,17 +14,18 @@ MODULE routing_reg
INTEGER(i_std), SAVE :: nbpt_save
INTEGER(i_std), SAVE :: nbasmax_save
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: routing_area_glo !! Surface of basin (m^2)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: routing_cg_glo !! Centre of gravity of HTU (Lat, Lon)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: topo_resid_glo !! Topographic index of the retention time (m)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: routing_area_glo !! Surface of basin (m^2)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: routing_cg_glo !! Centre of gravity of HTU (Lat, Lon)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: topo_resid_glo !! Topographic index of the retention time (m)
INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: route_count_glo !! Number of basins per grid
INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: route_togrid_glo !! Grid into which the basin flows (unitless)
INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: route_tobasin_glo !! Basin in to which the water goes (unitless)
INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: route_nbintobas_glo !! Number of basin into current one (unitless)
INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: origin_nbintobas_glo !! Number of sub-grid basin into current one before truncation (unitless)
INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: route_nbintobas_glo !! Number of basin into current one (unitless)
INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: origin_nbintobas_glo !! Number of sub-grid basin into current one before truncation (unitless)
INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: global_basinid_glo !! ID of basin (unitless)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: route_outlet_glo !! Coordinate of outlet (-)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: route_type_glo !! Coordinate of outlet (-)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: route_fetch_glo !! Upstream area at each HTU
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: route_outlet_glo !! Coordinate of outlet (-)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: route_type_glo !! Coordinate of outlet (-)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: route_fetch_glo !! Upstream area at each HTU
CONTAINS
!! ================================================================================================================================
......@@ -1608,7 +1609,7 @@ CONTAINS
!
ip = COUNT(triptemp(1:nbi,1:nbj) .EQ. -1)
IF ( ip .NE. nbi*nbj ) THEN
WRITE(*,*) 'Unprocessed points :', nbi*nbj-ip
WRITE(numout,*) 'Unprocessed points :', nbi*nbj-ip
ENDIF
!
END SUBROUTINE routing_reg_divbas
......@@ -2718,7 +2719,7 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, b
REAL(r_std), DIMENSION(nbpt*nwbas), INTENT(out) :: outflow_uparea
!
!! LOCAL VARIABLES
INTEGER(i_std) :: ib, ij, ic, ff(1), it, itt, igrif, ibasf, nboutflow !! Indices (unitless)
INTEGER(i_std) :: ig, ib, 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 !!
!
......@@ -2728,22 +2729,22 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, b
!
! Propagate first the partial sums from the halo into the domain
!
DO ib=1,nbpt
DO ig=1,nbpt
!
DO ij=1,basin_count(ib)
DO ib=1,basin_count(ig)
!
IF (partial_sum(ib,ij) > EPSILON(partial_sum)) THEN
IF (partial_sum(ig,ib) > EPSILON(partial_sum)) THEN
!
fetch_basin(ib, ij) = MAX(fetch_basin(ib, ij), partial_sum(ib,ij))
fetch_basin(ig, ib) = MAX(fetch_basin(ig, ib), partial_sum(ig,ib))
!
igrif = outflow_grid(ib,ij)
ibasf = outflow_basin(ib,ij)
igrif = outflow_grid(ig,ib)
ibasf = outflow_basin(ig,ib)
!
itt = 0
!
DO WHILE (igrif .GT. 0)
!
fetch_basin(igrif,ibasf) = MAX(fetch_basin(igrif,ibasf), partial_sum(ib, ij))
fetch_basin(igrif,ibasf) = MAX(fetch_basin(igrif,ibasf), partial_sum(ig, ib))
!
it = outflow_grid(igrif, ibasf)
ibasf = outflow_basin(igrif, ibasf)
......@@ -2762,20 +2763,20 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, b
! Add the areas contributed by the core region of the domain.
!
DO ic=1,nbcore
ib = corepts(ic)
ig = corepts(ic)
!
DO ij=1,basin_count(ib)
DO ib=1,basin_count(ig)
!
fetch_basin(ib, ij) = fetch_basin(ib, ij) + basin_area(ib,ij)
fetch_basin(ig,ib) = fetch_basin(ig,ib) + basin_area(ig,ib)
!
igrif = outflow_grid(ib,ij)
ibasf = outflow_basin(ib,ij)
igrif = outflow_grid(ig,ib)
ibasf = outflow_basin(ig,ib)
!
itt = 0
!
DO WHILE (igrif .GT. 0 )
!
fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + basin_area(ib, ij)
fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + basin_area(ig, ib)
!
it = outflow_grid(igrif, ibasf)
ibasf = outflow_basin(igrif, ibasf)
......@@ -2783,7 +2784,7 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, b
itt = itt + 1
IF ( itt .GT. maxiterations) THEN
WRITE(numout,&
"('Grid ',I5, ' and basin ',I5, ' did not converge after iteration ',I5)") ib, ij, itt
"('Grid ',I5, ' and basin ',I5, ' did not converge after iteration ',I5)") ig, ib, 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)
......@@ -2807,22 +2808,22 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, b
outflow_uparea(:) = zero
!
DO ic=1,nbcore
ib = corepts(ic)
ig = corepts(ic)
!
DO ij=1,basin_count(ib)
DO ib=1,basin_count(ig)
!
! We do not need any more the river flow flag as we are going to reset it.
! 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(ig,ib) .LT. 0) THEN
nboutflow = nboutflow + 1
IF ( nboutflow <= nbpt*nwbas) THEN
outflow_uparea(nboutflow) = fetch_basin(ib,ij)
outflow_uparea(nboutflow) = fetch_basin(ig,ib)
ELSE
! Not enought place so we replace a smaller one.
IF (fetch_basin(ib,ij) > MINVAL(outflow_uparea)) THEN
IF (fetch_basin(ig,ib) > MINVAL(outflow_uparea)) THEN
ff = MINLOC(outflow_uparea)
outflow_uparea(ff(1)) = fetch_basin(ib,ij)
outflow_uparea(ff(1)) = fetch_basin(ig,ib)
ELSE
! Ignore basin
ENDIF
......@@ -2874,9 +2875,9 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, b
!! \n
!_ ================================================================================================================================
SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, 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)
SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_largest, nbcore, corepts, 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)
!
IMPLICIT NONE
!
......@@ -2892,7 +2893,9 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_la
REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of land in each grid box (unitless;0-1)
!
INTEGER(i_std), INTENT(in) :: nwbas !!
INTEGER(i_std), INTENT (in) :: num_largest
INTEGER(i_std), INTENT (in) :: num_largest
INTEGER(i_std), INTENT (in) :: nbcore
INTEGER(i_std), DIMENSION(nbcore), INTENT(in) :: corepts
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 !!
......@@ -2913,7 +2916,7 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_la
INTEGER(i_std), DIMENSION(nbpt,nbxmax*8,nbxmax*8), INTENT(inout) :: inflow_grid !!
!
!! LOCAL VARIABLES
INTEGER(i_std) :: ib, ij, ibf, ijf, igrif, ibasf, cnt, pold, bold, ff(2) !! Indices (unitless)
INTEGER(i_std) :: ib, ij, ic, ibf, ijf, igrif, ibasf, cnt, pold, bold, ff(2) !! Indices (unitless)
INTEGER(i_std) :: ii, kbas, sbas, ik, iter, ibt, obj !! Indices (unitless)
REAL(r_std), DIMENSION(nbpt,nbasmax) :: floflo !!
REAL(r_std), DIMENSION(nbpt) :: gridarea_tmp !!
......@@ -2942,7 +2945,7 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_la
!
! Read option for truncation method
!
trunc_equa = .TRUE.
trunc_equa = .FALSE.
!!$ CALL getin('ROUTING_TRUNCEQUA', trunc_equa)
!
IF ( .NOT. ALLOCATED(indextrunc)) THEN
......@@ -2955,10 +2958,12 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_la
!
DO iter = 1, MAXVAL(basin_count) - nbasmax +3
!
! Get the points over which we wish to truncate
! Get the points over which we wish to truncate. We only work with points of the
! core area of the domain.
!
nbtruncate = 0
DO ib = 1, nbpt
DO ic = 1, nbcore
ib = corepts(ic)
IF ( basin_count(ib) .GT. nbasmax ) THEN
nbtruncate = nbtruncate + 1
indextrunc(nbtruncate) = ib
......@@ -3029,6 +3034,7 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_la
ENDIF
ENDIF
!
!
! 2) If we have basins flowing into the same grid but different basins then we put them
! together. Obviously we first work with the grid which has most streams running into it
! and putting the smallest in the largests catchments.
......@@ -3173,7 +3179,6 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_la
!
ENDDO
!