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

Corrected a bug on the cmputation of the partial sums and classification of rivers.

parent 5c1f707a
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)
......@@ -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
......
......@@ -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)
......
......@@ -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
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment