Skip to content

GitLab

Commit adbeb681 authored by Anthony's avatar Anthony
Browse files

Parallelization of the truncate + subroutine to check fetch coherence and detect loop

parent 328dc536
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 693 additions and 629 deletions
F90subroutines/routing_interface.f90
View file @ adbeb681
......@@ -417,7 +417,7 @@ SUBROUTINE rivclassification(nbpt, nwbas, nbcore, corepts, basin_count, outflow_
END SUBROUTINE rivclassification
SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridarea, cfrac, basin_count, basin_notrun, basin_area, &
SUBROUTINE finish_truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridarea, cfrac, 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)
......@@ -474,7 +474,7 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridarea, cfra
!
!
CALL routing_reg_truncate(nbpt, nbasmax, gridarea, cfrac, nwbas, num_largest, basin_count, &
CALL routing_reg_end_truncate(nbpt, nbasmax, gridarea, cfrac, 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)
......@@ -493,7 +493,188 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridarea, cfra
route_type(:,:) = route_type_glo(:,:)
origin_nbintobas(:) = origin_nbintobas_glo(:)
END SUBROUTINE truncate
END SUBROUTINE finish_truncate
SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numops, basin_count, 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)
!
!
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) :: nbxmax_in, 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,2), INTENT(inout) :: basin_cg
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_topoind !! Topographic index of the residence time for a basin (m)
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 !!
!
! Changed nwbas to nbxmax*4 and *8 to save the memory
!
INTEGER(i_std), DIMENSION(nbpt,nbxmax_in), INTENT(inout) :: inflow_number !!
INTEGER(i_std), DIMENSION(nbpt,nbxmax_in,nbxmax_in), INTENT(inout) :: inflow_basin !!
INTEGER(i_std), DIMENSION(nbpt,nbxmax_in,nbxmax_in), INTENT(inout) :: inflow_grid !!
!! LOCAL
INTEGER(i_std) :: ib, op, tok, totak
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
WRITE(numout,*) "nbpt", ib, "tokill", tok, "totakover", totak
CALL routing_reg_killbas(nbpt, ib, tok, totak, nwbas, basin_count, 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)
END IF
END IF
END DO
END DO
! REPRENDRE LA FIN DE TRUNCATE
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
flag(:,:) = 0
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) .AND. (flag(igrif,ibasf) .NE. -99))
flag(igrif, ibasf) = -99
it = outflow_grid(igrif,ibasf)
ibasf = outflow_basin(igrif,ibasf)
igrif = it
END DO
ELSE IF (flag(igrif,ibasf) .EQ. -1) THEN
test = test + 1
flag(ig,ib) = -99
igrif = outflow_grid(ig,ib)
ibasf = outflow_basin(ig,ib)
DO WHILE (flag(igrif,ibasf) .EQ. -1)
flag(igrif, ibasf) = -99
it = outflow_grid(igrif,ibasf)
ibasf = outflow_basin(igrif,ibasf)
igrif = it
END DO
END IF
END IF
END DO
END DO
WRITE(numout,*) "**** Fetch has",test, "loop errors"
END SUBROUTINE checkrouting
SUBROUTINE checkfetch(nbpt, nwbas, fetch_basin, 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 !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: fetch_basin !!
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) :: ig, ib, it, bt, igrif, ibasf, basnum, test
WRITE(numout,*) "Checking Fetch coherence"
test = 0
DO ig=1,nbpt
basnum = basin_count(ig)
DO ib=1,basnum
igrif = outflow_grid(ig,ib)
ibasf = outflow_basin(ig,ib)
IF (igrif .GT. 0) THEN
IF (fetch_basin(ig,ib) .GT. fetch_basin(igrif,ibasf)) THEN
WRITE(numout,*) "****ERROR Fetch, nbpt:",ig
WRITE(numout,*) igrif, ibasf
WRITE(numout,*) it, bt
test = 1
END IF
END IF
END DO
END DO
IF (test .EQ. 0) WRITE(numout,*) "**** Fetch has no error"
END SUBROUTINE checkfetch
SUBROUTINE finalfetch(nbpt, nbasmax, nbcore, corepts, routing_area, basin_count, route_togrid, route_tobasin, partial_sum, &
& fetch_basin, outflow_uparea)
......
F90subroutines/routing_reg.f90
View file @ adbeb681
......@@ -2554,7 +2554,7 @@ 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, &
SUBROUTINE routing_reg_end_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)
!
......@@ -2623,352 +2623,6 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_la
!
! Read option for truncation method
!
trunc_equa = .FALSE.
!!$ CALL getin('ROUTING_TRUNCEQUA', trunc_equa)
!
IF ( .NOT. ALLOCATED(indextrunc)) THEN
ALLOCATE(indextrunc(nbpt))
ENDIF
!
! We have to go through the grid as least as often as we have to reduce the number of basins
! For good measure we add 3 more passages.
!
!
DO iter = 1, MAXVAL(basin_count) - nbasmax +3
!
! 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
IF ( basin_count(ib) .GT. nbasmax ) THEN
nbtruncate = nbtruncate + 1
indextrunc(nbtruncate) = ib
ENDIF
ENDDO
!
! Go through the basins which need to be truncated.
!
DO ibt=1,nbtruncate
!
ib = indextrunc(ibt)
!
! Check if we have basin which flows into a basin in the same grid
! kbas = basin we will have to kill
! sbas = basin which takes over kbas
!
kbas = 0
sbas = 0
!
! 1) Merge two basins which flow into the ocean as coastal or return flow
! (outflow_grid = -2 or -3). Well obviously only if we have more than 1 and
! have not found anything yet!
!
IF ( (COUNT(outflow_grid(ib,1:basin_count(ib)) .EQ. -2) .GT. 1 .OR.&
& COUNT(outflow_grid(ib,1:basin_count(ib)) .EQ. -3) .GT. 1) .AND.&
& kbas*sbas .EQ. 0) THEN
!
multbas = 0
multbas_sz(:) = 0
!
IF ( COUNT(outflow_grid(ib,1:basin_count(ib)) .EQ. -2) .GT. 1 ) THEN
obj = -2
ELSE
obj = -3
ENDIF
!
! First we get the list of all basins which go out as coastal or return flow (obj)
!
DO ij=1,basin_count(ib)
IF ( outflow_grid(ib,ij) .EQ. obj ) THEN
multbas = multbas + 1
multbas_sz(multbas) = ij
tmp_area(multbas) = fetch_basin(ib,ij)
ENDIF
ENDDO
!
! Test the new truncation priority
!
IF ( trunc_equa ) THEN
!
! Take the smallest to be transfered to the second smallest
!
iml = MINLOC(tmp_area(1:multbas), MASK = tmp_area(1:multbas) .GT. zero)
kbas = multbas_sz(iml(1))
tmp_area(iml(1)) = -1
iml = MINLOC(tmp_area(1:multbas), MASK = tmp_area(1:multbas) .GT. zero)
sbas = multbas_sz(iml(1))
!
ELSE
!
! Take the smallest to be transfered to the largest
!
iml = MAXLOC(tmp_area(1:multbas), MASK = tmp_area(1:multbas) .GT. zero)
sbas = multbas_sz(iml(1))
iml = MINLOC(tmp_area(1:multbas), MASK = tmp_area(1:multbas) .GT. zero)
kbas = multbas_sz(iml(1))
!
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.
!
IF ( kbas*sbas .EQ. 0) THEN
!
tmp_ids(1:basin_count(ib)) = outflow_grid(ib,1:basin_count(ib))
multbas = 0
multbas_sz(:) = 0
!
! First obtain the list of basins which flow into the same basin
!
DO ij=1,basin_count(ib)
IF ( outflow_grid(ib,ij) .GT. 0 .AND.&
& COUNT(tmp_ids(1:basin_count(ib)) .EQ. outflow_grid(ib,ij)) .GT. 1) THEN
multbas = multbas + 1
DO ii=1,basin_count(ib)
IF ( tmp_ids(ii) .EQ. outflow_grid(ib,ij)) THEN
multbas_sz(multbas) = multbas_sz(multbas) + 1
multbas_list(multbas,multbas_sz(multbas)) = ii
tmp_ids(ii) = -99
ENDIF
ENDDO
ELSE
tmp_ids(ij) = -99
ENDIF
ENDDO
!
! Did we come up with any basins to deal with this way ?
!
IF ( multbas .GT. 0 ) THEN
!
iml = MAXLOC(multbas_sz(1:multbas))
ik = iml(1)
!
! Take the smallest and largest of these basins !
!
DO ii=1,multbas_sz(ik)
tmp_area(ii) = fetch_basin(ib, multbas_list(ik,ii))
ENDDO
!
IF ( trunc_equa ) THEN
iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
kbas = multbas_list(ik,iml(1))
tmp_area(iml(1)) = -1
iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
sbas = multbas_list(ik,iml(1))
ELSE
iml = MAXLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
sbas = multbas_list(ik,iml(1))
iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
kbas = multbas_list(ik,iml(1))
ENDIF
!
ENDIF
!
ENDIF
!
! 3) If we have twice the same basin we put them together even if they flow into different
! directions. If one of them goes to the ocean it takes the advantage.
!
IF ( kbas*sbas .EQ. 0) THEN
!
tmp_ids(1:basin_count(ib)) = basin_id(ib,1:basin_count(ib))
multbas = 0
multbas_sz(:) = 0
!
! First obtain the list of basins which have sub-basins in this grid box.
! (these are identified by their IDs)
!
DO ij=1,basin_count(ib)
IF ( COUNT(tmp_ids(1:basin_count(ib)) .EQ. basin_id(ib,ij)) .GT. 1) THEN
multbas = multbas + 1
DO ii=1,basin_count(ib)
IF ( tmp_ids(ii) .EQ. basin_id(ib,ij)) THEN
multbas_sz(multbas) = multbas_sz(multbas) + 1
multbas_list(multbas,multbas_sz(multbas)) = ii
tmp_ids(ii) = -99
ENDIF
ENDDO
ELSE
tmp_ids(ij) = -99
ENDIF
ENDDO
!
! We are going to work on the basin with the largest number of sub-basins.
! (IF we have a basin which has subbasins !)
!
IF ( multbas .GT. 0 ) THEN
!
iml = MAXLOC(multbas_sz(1:multbas))
ik = iml(1)
!
! If one of the basins goes to the ocean then it is going to have the priority
!
tmp_area(:) = zero
IF ( COUNT(outflow_grid(ib,multbas_list(ik,1:multbas_sz(ik))) .LT. 0) .GT. 0) THEN
DO ii=1,multbas_sz(ik)
IF ( outflow_grid(ib,multbas_list(ik,ii)) .LT. 0 .AND. sbas .EQ. 0 ) THEN
sbas = multbas_list(ik,ii)
ELSE
tmp_area(ii) = fetch_basin(ib, multbas_list(ik,ii))
ENDIF
ENDDO
! take the smallest of the subbasins
iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
kbas = multbas_list(ik,iml(1))
ELSE
!
! Else we take simply the largest and smallest
!
DO ii=1,multbas_sz(ik)
tmp_area(ii) = fetch_basin(ib, multbas_list(ik,ii))
ENDDO
!
IF ( trunc_equa ) THEN
iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
kbas = multbas_list(ik,iml(1))
tmp_area(iml(1)) = -1
iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
sbas = multbas_list(ik,iml(1))
ELSE
iml = MAXLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
sbas = multbas_list(ik,iml(1))
iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
kbas = multbas_list(ik,iml(1))
ENDIF
!
ENDIF
!
!
ENDIF
ENDIF
!
! Then we call routing_reg_killbas to clean up the basins in this grid
!
IF ( kbas .GT. 0 .AND. sbas .GT. 0 ) THEN
CALL routing_reg_killbas(nbpt, ib, kbas, sbas, nwbas, basin_count, 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)
ENDIF
!
ENDDO
!
ENDDO
!
! We only merge a basin which flows into a basin of the same grid if it's
! necessary.
!
IF ( COUNT(basin_count .GT. nbasmax) .GT. 0 ) THEN
!
! Get the points over which we wish to truncate
!
nbtruncate = 0
DO ib = 1, nbpt
IF ( basin_count(ib) .GT. nbasmax ) THEN
nbtruncate = nbtruncate + 1
indextrunc(nbtruncate) = ib
ENDIF
ENDDO
!
! Go through the basins which need to be truncated.
!
DO ibt=1,nbtruncate
!
ib = indextrunc(ibt)
!
! Check if we have basin which flows into a basin in the same grid
! kbas = basin we will have to kill
! sbas = basin which takes over kbas
!
kbas = 0
sbas = 0
!
! 4) Can we find a basin which flows into a basin of the same grid ?
!
DO ij=1,basin_count(ib)
DO ii=1,basin_count(ib)
IF ( outflow_grid(ib,ii) .EQ. ib .AND. outflow_basin(ib, ii) .EQ. ij .AND. kbas*sbas .NE. 0) THEN
kbas = ii
sbas = ij
ENDIF
ENDDO
ENDDO
!
! Then we call routing_reg_killbas to clean up the basins in this grid
!
IF ( kbas .GT. 0 .AND. sbas .GT. 0 ) THEN
CALL routing_reg_killbas(nbpt, ib, kbas, sbas, nwbas, basin_count, 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)
ENDIF
!
ENDDO
!
ENDIF
!
! If there are any grids left with too many basins we need to take out the big hammer !
! We will only do it if this represents less than 5% of all points.
!
IF ( COUNT(basin_count .GT. nbasmax) .GT. 0 ) THEN
!
!
IF ( COUNT(basin_count .GT. nbasmax)/nbpt*100 .GT. 5 ) THEN
WRITE(numout,*) 'We have ', COUNT(basin_count .GT. nbasmax)/nbpt*100, '% of all points which do not yet'
WRITE(numout,*) 'have the right trunctaction. That is too much to apply a brutal method'
DO ib = 1, nbpt
IF ( basin_count(ib) .GT. nbasmax ) THEN
!
WRITE(numout,*) 'We did not find a basin which could be supressed. We will'
WRITE(numout,*) 'not be able to reduce the truncation in grid ', ib
DO ij=1,basin_count(ib)
WRITE(numout,*) 'grid, basin nb and id :', ib, ij, basin_id(ib,ij)
WRITE(numout,*) 'Outflow grid and basin ->', outflow_grid(ib,ij), outflow_basin(ib, ij)
ENDDO
ENDIF
ENDDO
CALL ipslerr_p(3,'routing_reg_truncate','No basin found which could be supressed.','','')
ELSE
!
!
DO ib = 1,nbpt
DO WHILE ( basin_count(ib) .GT. nbasmax )
!
! Here we simply put the smallest basins into the largest ones. It is really a brute force
! method but it will only be applied if everything has failed.
!
DO ii = 1,basin_count(ib)
tmp_area(ii) = fetch_basin(ib, ii)
ENDDO
!
IF ( trunc_equa ) THEN
iml = MINLOC(tmp_area(1:basin_count(ib)), MASK = tmp_area(1:basin_count(ib)) .GT. 0.)
kbas = iml(1)
tmp_area(iml(1)) = -1
iml = MINLOC(tmp_area(1:basin_count(ib)), MASK = tmp_area(1:basin_count(ib)) .GT. 0.)
sbas =iml(1)
ELSE
iml = MAXLOC(tmp_area(1:basin_count(ib)), MASK = tmp_area(1:basin_count(ib)) .GT. 0.)
sbas =iml(1)
iml = MINLOC(tmp_area(1:basin_count(ib)), MASK = tmp_area(1:basin_count(ib)) .GT. 0.)
kbas = iml(1)
ENDIF
!
IF ( kbas .GT. 0 .AND. sbas .GT. 0 ) THEN
CALL routing_reg_killbas(nbpt, ib, kbas, sbas, nwbas, basin_count, 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)
ENDIF
ENDDO