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

Implementation of the simplificatino of river branches in...

Implementation of the simplificatino of river branches in routing_reg_finbasins. It allows to avoid small outflows from the atmopsheric grid and reduces the work on truncate later on.
parent 428a0435
...@@ -30,7 +30,7 @@ SUBROUTINE initatmgrid(rankmpi, nbcore, nbpt, nbseg, polygons_in, centre_in, are ...@@ -30,7 +30,7 @@ SUBROUTINE initatmgrid(rankmpi, nbcore, nbpt, nbseg, polygons_in, centre_in, are
CALL grid_init(nbpt, nbseg, polygons_in, centre_in, area_in, contfrac_in, neighbours_loc) CALL grid_init(nbpt, nbseg, polygons_in, centre_in, area_in, contfrac_in, neighbours_loc)
! !
numout=100+rankmpi numout=100+rankmpi
WRITE(outfname,'(A15,I4.4,A6,I4.4,A4)') "Out_RoutingReg_",rankmpi,"_from_",nbcore,".txt" WRITE(outfname,'(A15,I4.4,A6,I4.4,A4)') "Out_RoutingReg_",INT(rankmpi),"_from_",INT(nbcore),".txt"
OPEN(unit=numout, file=outfname) OPEN(unit=numout, file=outfname)
! !
END SUBROUTINE initatmgrid END SUBROUTINE initatmgrid
...@@ -126,7 +126,7 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, maxpercent_in, sub_pts, sub_i ...@@ -126,7 +126,7 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, maxpercent_in, sub_pts, sub_i
END SUBROUTINE gethydrogrid END SUBROUTINE gethydrogrid
SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, trip_bx, basin_bx, fac_bx, hierarchy_bx, topoind_bx, & SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, nbasmax, 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, & & 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_sz, basin_bxout, basin_bbout, basin_pts, basin_lshead, coast_pts, lontmp, lattmp, orog_bx, &
& area_bx, basin_topoindex_stream, basin_rlen_stream, basin_dz_stream) & area_bx, basin_topoindex_stream, basin_rlen_stream, basin_dz_stream)
...@@ -141,10 +141,11 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, trip_bx, basin_bx, ...@@ -141,10 +141,11 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, trip_bx, basin_bx,
INTEGER, INTENT(in) :: nb_htu, nbv, ijdimmax INTEGER, INTENT(in) :: nb_htu, nbv, ijdimmax
INTEGER, INTENT(in) :: nbpt !! Domain size (unitless) 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(in) :: nbi(nbpt), nbj(nbpt) !! Number of point in x and y within the grid (unitless)
INTEGER, INTENT(in) :: nbasmax
INTEGER, INTENT(inout) :: trip_bx(nbpt,ijdimmax,ijdimmax) !! The trip field for each of the smaller boxes (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) !! INTEGER, INTENT(inout) :: basin_bx(nbpt,ijdimmax,ijdimmax) !!
REAL, INTENT(in) :: fac_bx(nbpt,ijdimmax,ijdimmax) !! Flow accumulation REAL, INTENT(inout) :: 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) :: 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) :: 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) :: rlen_bx(nbpt,ijdimmax,ijdimmax) !!
REAL, INTENT(inout) :: rdz_bx(nbpt,ijdimmax,ijdimmax) !! REAL, INTENT(inout) :: rdz_bx(nbpt,ijdimmax,ijdimmax) !!
...@@ -184,13 +185,15 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, trip_bx, basin_bx, ...@@ -184,13 +185,15 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, trip_bx, basin_bx,
IF ( debug) WRITE(numout,*) "Memory Mgt findbasin : nbvmax, nb_htu, nbv = ", nbvmax, nb_htu, nbv IF ( debug) WRITE(numout,*) "Memory Mgt findbasin : nbvmax, nb_htu, nbv = ", nbvmax, nb_htu, nbv
DO ib=1,nbpt DO ib=1,nbpt
CALL routing_reg_findbasins(nb_htu, nbv, ib, ijdimmax, nbi(ib), nbj(ib), trip_bx(ib,:,:), & CALL routing_reg_findbasins(nb_htu, nbv, ib, ijdimmax, nbi(ib), nbj(ib), nbasmax, trip_bx(ib,:,:), &
& basin_bx(ib,:,:), fac_bx(ib,:,:), hierarchy_bx(ib,:,:), topoind_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), & & 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_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,:,:),& & basin_bbout(ib,:), basin_pts(ib,:,:,:), basin_lshead(ib,:), coast_pts(ib,:), lontmp(ib,:,:), lattmp(ib,:,:),&
& orog_bx(ib,:,:), area_bx(ib,:,:),basin_topoindex_stream(ib,:), basin_rlen_stream(ib,:), basin_dz_stream(ib,:)) & orog_bx(ib,:,:), area_bx(ib,:,:),basin_topoindex_stream(ib,:), basin_rlen_stream(ib,:), basin_dz_stream(ib,:))
ENDDO ENDDO
WRITE(numout,*) "XXXX range of nb_basin after findbasin : ", MINVAL(nb_basin), sum(nb_basin)/nbpt, MAXVAL(nb_basin)
END SUBROUTINE findbasins END SUBROUTINE findbasins
......
...@@ -430,7 +430,7 @@ CONTAINS ...@@ -430,7 +430,7 @@ CONTAINS
!! \n !! \n
!_ ================================================================================================================================ !_ ================================================================================================================================
SUBROUTINE routing_reg_findbasins(nb_htu, nbv, ib, ijdimmax, nbi, nbj, trip, basin, fac, hierarchy, topoind, & SUBROUTINE routing_reg_findbasins(nb_htu, nbv, ib, ijdimmax, nbi, nbj, nbasmax, trip, basin, fac, hierarchy, topoind, &
& rlen, rdz, rweight, lshead, diaglalo, nb_basin, basin_inbxid, basin_outlet, basin_outtp, basin_sz, & & rlen, rdz, rweight, lshead, diaglalo, nb_basin, basin_inbxid, basin_outlet, basin_outtp, basin_sz, &
& basin_bxout, basin_bbout, basin_pts, basin_lshead, coast_pts, lontmp, lattmp, orog, area, & & basin_bxout, basin_bbout, basin_pts, basin_lshead, coast_pts, lontmp, lattmp, orog, area, &
& basin_topoindex_stream, basin_rlen_stream, basin_dz_stream ) & basin_topoindex_stream, basin_rlen_stream, basin_dz_stream )
...@@ -443,8 +443,7 @@ CONTAINS ...@@ -443,8 +443,7 @@ CONTAINS
INTEGER(i_std), INTENT(in) :: ijdimmax INTEGER(i_std), INTENT(in) :: ijdimmax
INTEGER(i_std), INTENT(in) :: nbi !! Number of point in x within the grid (unitless) INTEGER(i_std), INTENT(in) :: nbi !! Number of point in x within the grid (unitless)
INTEGER(i_std), INTENT(in) :: nbj !! Number of point in y within the grid (unitless) INTEGER(i_std), INTENT(in) :: nbj !! Number of point in y within the grid (unitless)
REAL(r_std), INTENT(in) :: hierarchy(ijdimmax,ijdimmax) !! INTEGER(i_std), INTENT(in) :: nbasmax
REAL(r_std), INTENT(in) :: fac(ijdimmax,ijdimmax) !!
REAL(r_std), INTENT(in) :: lshead(ijdimmax,ijdimmax) REAL(r_std), INTENT(in) :: lshead(ijdimmax,ijdimmax)
REAL(r_std), INTENT(in) :: orog(ijdimmax,ijdimmax) REAL(r_std), INTENT(in) :: orog(ijdimmax,ijdimmax)
REAL(r_std), INTENT(in) :: area(ijdimmax,ijdimmax) REAL(r_std), INTENT(in) :: area(ijdimmax,ijdimmax)
...@@ -453,6 +452,8 @@ CONTAINS ...@@ -453,6 +452,8 @@ CONTAINS
! Modified ! Modified
! !
INTEGER(i_std), INTENT(inout) :: trip(ijdimmax,ijdimmax) !! The trip field (unitless) INTEGER(i_std), INTENT(inout) :: trip(ijdimmax,ijdimmax) !! The trip field (unitless)
REAL(r_std), INTENT(inout) :: fac(ijdimmax,ijdimmax) !!
REAL(r_std), INTENT(inout) :: hierarchy(ijdimmax,ijdimmax) !!
INTEGER(i_std), INTENT(inout) :: basin(ijdimmax,ijdimmax) !! INTEGER(i_std), INTENT(inout) :: basin(ijdimmax,ijdimmax) !!
REAL(r_std), INTENT(inout) :: topoind(ijdimmax,ijdimmax) !! Topographic index of the residence time (km) REAL(r_std), INTENT(inout) :: topoind(ijdimmax,ijdimmax) !! Topographic index of the residence time (km)
REAL(r_std), INTENT(inout) :: rlen(ijdimmax,ijdimmax) !! River length (m) REAL(r_std), INTENT(inout) :: rlen(ijdimmax,ijdimmax) !! River length (m)
...@@ -495,7 +496,7 @@ CONTAINS ...@@ -495,7 +496,7 @@ CONTAINS
INTEGER(i_std) :: tbname(nb_htu) !! INTEGER(i_std) :: tbname(nb_htu) !!
INTEGER(i_std) :: tsz(nb_htu) !! INTEGER(i_std) :: tsz(nb_htu) !!
INTEGER(i_std) :: tpts(nb_htu,nbv,2) !! INTEGER(i_std) :: tpts(nb_htu,nbv,2) !!
INTEGER(i_std) :: slsz INTEGER(i_std) :: slsz, slsz_old
INTEGER(i_std) :: slbas(nb_htu) INTEGER(i_std) :: slbas(nb_htu)
INTEGER(i_std) :: toutdir(nb_htu) !! INTEGER(i_std) :: toutdir(nb_htu) !!
INTEGER(i_std) :: toutbas(nb_htu) !! INTEGER(i_std) :: toutbas(nb_htu) !!
...@@ -538,9 +539,14 @@ CONTAINS ...@@ -538,9 +539,14 @@ CONTAINS
INTEGER(i_std) :: option, il_trib, jl_trib, changes !! INTEGER(i_std) :: option, il_trib, jl_trib, changes !!
INTEGER(i_std) :: nbdiv, il_upstr, jl_upstr INTEGER(i_std) :: nbdiv, il_upstr, jl_upstr
REAL(r_std) :: fac_glo_trib, fac_loc_trib, fac_loc_main, fac_loc_upstr REAL(r_std) :: fac_glo_trib, fac_loc_trib, fac_loc_main, fac_loc_upstr
REAL(r_std), PARAMETER :: htufrac_min=0.01 REAL(r_std) :: htufrac_min=0.1
! !
!_ ================================================================================================================================ !_ ================================================================================================================================
!
! We get some estimate of the smallest HTU we will try to maintain. Smaller than a fraction
! htufrac_min of the atmopsheric grid will be merged with others at this stage.
!
htufrac_min = 1.0/MAX(REAL(nbasmax),4.0)
! !
IF ( debug .AND. routing_diagbox(ib, diaglalo) ) THEN IF ( debug .AND. routing_diagbox(ib, diaglalo) ) THEN
WRITE(numout,*) "Point: ", ib WRITE(numout,*) "Point: ", ib
...@@ -554,50 +560,16 @@ CONTAINS ...@@ -554,50 +560,16 @@ CONTAINS
ENDDO ENDDO
ENDIF ENDIF
! !
! 1.0 Find number of outflow point (>= 97) ! 1.0 Follow the flow of the water to get size and gather point for each
! So basically, we get nb_basin, basin_bxout, basin_inbxid ! sub-basin. The lists of all the basins with an outflow of the gird are produced.
! (with nbb, toutdir, tbname) ! So we get basin_sz, basin_pts (with tsz, tpts)
!
nbb = 0
tbname(:) = undef_int
toutdir(:) = 0
toutloc(:,:) = 0
outtmp(:,:) = -1
!
toutbas(:) = undef_int
touttp(:) = undef_int
toutlshead(:) = undef_sechiba
!
DO ip=1,nbi
DO jp=1,nbj
!
!
IF ( basin(ip,jp) .LT. undef_int) THEN
! Because, there is maybe difference between undef_int and undef in routing.nc
!IF ( trip(ip,jp) .GE. 97 .AND. trip(ip,jp) .LT. undef_int ) THEN
IF ( trip(ip,jp) .GE. 97 .AND. trip(ip,jp) .LT. 109 ) THEN
nbb = nbb + 1
IF ( nbb .GT. nbvmax ) CALL ipslerr_p(3,'routing_reg_findbasins','nbvmax too small','first section','')
tbname(nbb) = basin(ip,jp)
toutdir(nbb) = trip(ip,jp)-100
toutlshead(nbb) = lshead(ip,jp)
touttp(nbb) = trip(ip,jp)
toutloc(nbb,1) = ip
toutloc(nbb,2) = jp
tlat(nbb) = lattmp(ip,jp)
tlon(nbb) = lontmp(ip,jp)
outtmp(ip,jp) = nbb
ENDIF
ENDIF
ENDDO
ENDDO
!
! 2.0 Follow the flow of the water to get size and gather point for each
! sub-basin. So we get basin_sz, basin_pts (with tsz, tpts)
! !
CALL routing_reg_complocalfac(ib, nbi, nbj, nb_htu, nbv, nbb, basin, trip, outtmp, totsz, tsz, tpts, fac_loc) CALL routing_reg_localbas(ib, nbi, nbj, nb_htu, nbv, basin, trip, lshead, lattmp, lontmp,&
& nbb, tbname, toutdir, toutbas, toutlshead, touttp, toutloc, tlat, tlon, &
& outtmp, totsz, tsz, tpts, fac_loc)
! !
slsz=0 slsz=0
slsz_old=-1
DO ibas = 1, nbb DO ibas = 1, nbb
! We know the first element is the outflow point of the HTU ! We know the first element is the outflow point of the HTU
IF ( fac(tpts(ibas,1,1),tpts(ibas,1,2)) == fac_loc(tpts(ibas,1,1),tpts(ibas,1,2)) .AND. & IF ( fac(tpts(ibas,1,1),tpts(ibas,1,2)) == fac_loc(tpts(ibas,1,1),tpts(ibas,1,2)) .AND. &
...@@ -608,25 +580,50 @@ CONTAINS ...@@ -608,25 +580,50 @@ CONTAINS
ENDIF ENDIF
ENDDO ENDDO
! !
IF ( slsz > 0 ) THEN ! 1.1 If we have basins smaller than htufrac_min we try to redirect them into another sub-basin
color(:,:) = 0 ! of the same river. This will minimize the number of outflow of this grid for a given river.
DO jp = 1, slsz !
ibas = slbas(jp) DO WHILE ( slsz > 0 .AND. slsz .NE. slsz_old)
!
IF ( debug ) THEN
WRITE(numout,*) "XXXX ",ib," Going into simplify with : nbb=",nbb," slsz=", slsz
WRITE(numout,*) '+++++++++++++++++++ fac_loc IN FINDBASINS ++++++++++++++++++++'
WRITE(fmtr,"('(',I3,'F8.1)')") nbi
DO je=1,nbj
WRITE(numout,fmtr) fac_loc(1:nbi,je)
ENDDO
WRITE(numout,*) '+++++++++++++++++++ fac IN FINDBASINS ++++++++++++++++++++'
WRITE(fmtr,"('(',I3,'F8.1)')") nbi
DO je=1,nbj
WRITE(numout,fmtr) fac(1:nbi,je)
ENDDO
ENDIF
!
CALL routing_reg_simplifyborders(ib, nbi, nbj, nb_htu, nbv, nbb, basin, trip, fac, hierarchy, &
& slsz, slbas, tsz, tpts, fac_loc)
!
CALL routing_reg_localbas(ib, nbi, nbj, nb_htu, nbv, basin, trip, lshead, lattmp, lontmp, &
& nbb, tbname, toutdir, toutbas, toutlshead, touttp, toutloc, tlat, tlon, &
& outtmp, totsz, tsz, tpts, fac_loc)
!
slsz_old=slsz
slsz=0
DO ibas = 1, nbb
! We know the first element is the outflow point of the HTU ! We know the first element is the outflow point of the HTU
IF ( fac(tpts(ibas,1,1),tpts(ibas,1,2)) == fac_loc(tpts(ibas,1,1),tpts(ibas,1,2)) ) THEN IF ( fac(tpts(ibas,1,1),tpts(ibas,1,2)) == fac_loc(tpts(ibas,1,1),tpts(ibas,1,2)) .AND. &
DO ip=1,tsz(ibas) & tsz(ibas)/ REAL(totsz) < htufrac_min ) THEN
color(tpts(ibas,ip,1),tpts(ibas,ip,2)) = outtmp(tpts(ibas,1,1),tpts(ibas,1,2)) ! This HTU could be deleted !
ENDDO slsz = slsz+1
slbas(slsz) = ibas
ENDIF ENDIF
ENDDO ENDDO
ENDIF ENDDO
! !
! !
IF ( debug .AND. slsz > 0 ) THEN IF ( debug .AND. slsz > 0 ) THEN
WRITE(numout,*) "=== To check grid before dividing sub-basin ==="
WRITE(numout,*) "Number of points: ", nbi, nbj WRITE(numout,*) "Number of points: ", nbi, nbj
WRITE(numout,*) "Number of outlet: ", nbb WRITE(numout,*) "Number of outlet: ", nbb
WRITE(numout,*) "Number of HTU smaller than ",htufrac_min, " = ", slsz WRITE(numout,*) "Number of HTU smaller than ",htufrac_min, " slsz = ", slsz
WRITE(numout,*) '+++++++++++++++++++ TRIP IN FINDBASINS ++++++++++++++++++++' WRITE(numout,*) '+++++++++++++++++++ TRIP IN FINDBASINS ++++++++++++++++++++'
WRITE(fmt,"('(',I3,'I8)')") nbi WRITE(fmt,"('(',I3,'I8)')") nbi
DO je=1,nbj DO je=1,nbj
...@@ -634,7 +631,7 @@ CONTAINS ...@@ -634,7 +631,7 @@ CONTAINS
ENDDO ENDDO
WRITE(numout,*) '+++++++++++++++++++ LOCAL BASINS IN FINDBASINS ++++++++++++++++++++' WRITE(numout,*) '+++++++++++++++++++ LOCAL BASINS IN FINDBASINS ++++++++++++++++++++'
DO je=1,nbj DO je=1,nbj
WRITE(numout,fmt) color(1:nbi,je) WRITE(numout,fmt) outtmp(1:nbi,je)
ENDDO ENDDO
WRITE(numout,*) '+++++++++++++++++++ fac_loc IN FINDBASINS ++++++++++++++++++++' WRITE(numout,*) '+++++++++++++++++++ fac_loc IN FINDBASINS ++++++++++++++++++++'
WRITE(fmtr,"('(',I3,'F8.1)')") nbi WRITE(fmtr,"('(',I3,'F8.1)')") nbi
...@@ -1016,10 +1013,108 @@ CONTAINS ...@@ -1016,10 +1013,108 @@ CONTAINS
END SUBROUTINE routing_reg_findbasins END SUBROUTINE routing_reg_findbasins
! !
!! ================================================================================================================================ !! ================================================================================================================================
!! SUBROUTINE : routing_reg_complocalfac !! SUBROUTINE : routing_reg_simplifyborders
!!
!>\BRIEF This subroutine will simplify the borders of the atmospheric grid by joining
!! rivers whith different branches flowing into the same direction.
!!
!! DESCRIPTION (definitions, functional, design, flags) :
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S):
!!
!! REFERENCES : None
!!
!! FLOWCHART : None
!! \n
!_ ================================================================================================================================
SUBROUTINE routing_reg_simplifyborders(ib, nbi, nbj, nb_htu, nbv, nbb, basin, trip, fac, hierarchy, &
& slsz, slbas, tsz, tpts, fac_loc)
!
! INPUT
!
INTEGER(i_std), INTENT(in) :: ib !! Grid box on which we are currently working
INTEGER(i_std), INTENT(in) :: nbi !! Number of point in x within the grid (unitless)
INTEGER(i_std), INTENT(in) :: nbj !! Number of point in y within the grid (unitless)
INTEGER(i_std), INTENT(in) :: nb_htu
INTEGER(i_std), INTENT(in) :: nbv
INTEGER(i_std), INTENT(in) :: nbb !! Number of outflow points
!
INTEGER(i_std), INTENT(in) :: slsz
INTEGER(i_std), INTENT(in) :: slbas(nb_htu)
INTEGER(i_std), INTENT(in) :: tsz(nb_htu)
INTEGER(i_std), INTENT(in) :: tpts(nb_htu,nbv,2)
!
! INOUT
!
INTEGER(i_std), INTENT(in) :: basin(:,:)
INTEGER(i_std), INTENT(inout) :: trip(:,:)
REAL(r_std), INTENT(inout) :: fac(:,:)
REAL(r_std), INTENT(inout) :: hierarchy(:,:)
REAL(r_std), INTENT(inout) :: fac_loc(:,:)
!
! LOCAL
!
INTEGER(i_std) :: is, ir, ibas, ip, jp, ipp, jpp
INTEGER(i_std) :: nbopt, opt(nbne,2), optdir(nbne)
REAL(r_std) :: optsz(nbne)
INTEGER(i_std) :: ff(1)
!
! Go through all te small basins to see if we can merge them.
!
DO is=1,slsz
ibas = slbas(is)
! We know the first point in the list tpts is the outflow point.
ip=tpts(ibas,1,1)
jp=tpts(ibas,1,2)
!
nbopt=0
optsz(:)=zero
DO ir=1,nbne
ipp = ip+inc(ir,1)
jpp = jp+inc(ir,2)
IF ( ipp .GT. 0 .AND. ipp .LE. nbi .AND. jpp .GT. 0 .AND. jpp .LE. nbj ) THEN
! Three conditions to merge outflow points !
! 1) The target point has to be an outflow as well
! 2) It has to belong to the same basin
! 3) Its flow accumulation has to be larger or equal.
IF ( trip(ipp,jpp) .GE. 97 .AND. basin(ip,jp) == basin(ipp,jpp) .AND. &
& fac_loc(ip,jp) .LE. fac_loc(ipp,jpp) ) THEN
nbopt = nbopt+1
opt(nbopt,1)=ipp
opt(nbopt,2)=jpp
optdir(nbopt)=ir
optsz(nbopt)=fac_loc(ipp,jpp)
ENDIF
ENDIF
ENDDO
!
! Exploite options found
!
IF ( nbopt > 0 ) THEN
ff=MAXLOC(optsz)
ir=ff(1)
! Modify the fields which need to be adjusted to eliminate the outflow point
trip(ip,jp) = optdir(ir)
! Add the flow accumulation of the redirected basin +1 for the extra pixel gathered.
fac(opt(ir,1),opt(ir,2)) = fac(opt(ir,1),opt(ir,2)) + fac(ip,jp) + 1
! As we use the minimum hierarchy of the HTU no need to change its value here.
! But that could be different.
! So the lines are here as a reminder.
!!hierarchy() = hierarchy()
ENDIF
ENDDO
!
END SUBROUTINE routing_reg_simplifyborders
!
!
!! ================================================================================================================================
!! SUBROUTINE : routing_reg_localbas
!! !!
!>\BRIEF This subroutine gathers all the points with the same outflow and computes the flow accumulation !>\BRIEF This subroutine gathers all the points with the same outflow and computes the flow accumulation
!! over the entire catchment. !! over the entire catchment. It essentially determines the local basins relative to the atmospheric
!! grid.
!! !!
!! DESCRIPTION (definitions, functional, design, flags) : !! DESCRIPTION (definitions, functional, design, flags) :
!! !!
...@@ -1033,7 +1128,9 @@ CONTAINS ...@@ -1033,7 +1128,9 @@ CONTAINS
!! \n !! \n
!_ ================================================================================================================================ !_ ================================================================================================================================
SUBROUTINE routing_reg_complocalfac(ib, nbi, nbj, nb_htu, nbv, nbb, basin, trip, outtmp, totsz, tsz, tpts, fac_loc) SUBROUTINE routing_reg_localbas(ib, nbi, nbj, nb_htu, nbv, basin, trip, lshead, lattmp, lontmp, &
& nbb, tbname, toutdir, toutbas, toutlshead, touttp, toutloc, tlat, tlon, &
& outtmp, totsz, tsz, tpts, fac_loc)
! !
!! INPUT Variables !! INPUT Variables
INTEGER(i_std), INTENT(in) :: ib !! Grid box on which we are currently working INTEGER(i_std), INTENT(in) :: ib !! Grid box on which we are currently working
...@@ -1041,15 +1138,26 @@ CONTAINS ...@@ -1041,15 +1138,26 @@ CONTAINS
INTEGER(i_std), INTENT(in) :: nbj !! Number of point in y within the grid (unitless) INTEGER(i_std), INTENT(in) :: nbj !! Number of point in y within the grid (unitless)
INTEGER(i_std), INTENT(in) :: nb_htu INTEGER(i_std), INTENT(in) :: nb_htu
INTEGER(i_std), INTENT(in) :: nbv INTEGER(i_std), INTENT(in) :: nbv
INTEGER(i_std), INTENT(in) :: nbb !! Number of outflow points
INTEGER(i_std), INTENT(in) :: basin(:,:) INTEGER(i_std), INTENT(in) :: basin(:,:)
INTEGER(i_std), INTENT(in) :: trip(:,:) INTEGER(i_std), INTENT(in) :: trip(:,:)
!! REAL(r_std), INTENT(in) :: lshead(:,:)
INTEGER(i_std), INTENT(inout) :: outtmp(:,:) REAL(r_std), INTENT(in) :: lattmp(:,:) !! Latitude
REAL(r_std), INTENT(in) :: lontmp(:,:) !! Longitude
!
!! OUTPUT Variables !! OUTPUT Variables
INTEGER(i_std), INTENT(out) :: nbb !! Number of outflow points
INTEGER(i_std), INTENT(out) :: tbname(nb_htu) !!
INTEGER(i_std), INTENT(out) :: totsz INTEGER(i_std), INTENT(out) :: totsz
INTEGER(i_std), INTENT(out) :: tsz(nb_htu) INTEGER(i_std), INTENT(out) :: tsz(nb_htu)
INTEGER(i_std), INTENT(out) :: tpts(nb_htu,nbv,2) INTEGER(i_std), INTENT(out) :: tpts(nb_htu,nbv,2)
INTEGER(i_std), INTENT(out) :: toutdir(nb_htu) !!
INTEGER(i_std), INTENT(out) :: toutbas(nb_htu) !!
REAL(r_std), INTENT(out) :: toutlshead(nb_htu) !!
REAL(r_std), INTENT(out) :: touttp(nb_htu) !!
INTEGER(i_std), INTENT(out) :: toutloc(nb_htu,2) !!
REAL(r_std), INTENT(out) :: tlon(nb_htu) !!
REAL(r_std), INTENT(out) :: tlat(nb_htu) !!
INTEGER(i_std), INTENT(out) :: outtmp(:,:)
REAL(r_std), INTENT(out) :: fac_loc(nbi,nbj) REAL(r_std), INTENT(out) :: fac_loc(nbi,nbj)
!! LOCAL !! LOCAL
INTEGER(i_std) :: ip, jp, il, jl, ill, jll INTEGER(i_std) :: ip, jp, il, jl, ill, jll
...@@ -1059,6 +1167,37 @@ CONTAINS ...@@ -1059,6 +1167,37 @@ CONTAINS
! !
! Initialisation ! Initialisation
! !
nbb = 0
tbname(:) = undef_int
toutdir(:) = 0
toutloc(:,:) = 0
outtmp(:,:) = -1
!
toutbas(:) = undef_int
touttp(:) = undef_int
toutlshead(:) = undef_sechiba
!
DO ip=1,nbi
DO jp=1,nbj
!
IF ( basin(ip,jp) .LT. undef_int) THEN
IF ( trip(ip,jp) .GE. 97 .AND. trip(ip,jp) .LT. 109 ) THEN
nbb = nbb + 1
IF ( nbb .GT. nbvmax ) CALL ipslerr_p(3,'routing_reg_findbasins','nbvmax too small','first section','')
tbname(nbb) = basin(ip,jp)
toutdir(nbb) = trip(ip,jp)-100
toutlshead(nbb) = lshead(ip,jp)
touttp(nbb) = trip(ip,jp)
toutloc(nbb,1) = ip
toutloc(nbb,2) = jp
tlat(nbb) = lattmp(ip,jp)
tlon(nbb) = lontmp(ip,jp)
outtmp(ip,jp) = nbb
ENDIF
ENDIF
ENDDO
ENDDO
!
tsz(:) = 0 tsz(:) = 0
tpts(:,:,:) = 0 tpts(:,:,:) = 0
fac_loc(:,:) = zero fac_loc(:,:) = zero
...@@ -1117,13 +1256,13 @@ CONTAINS ...@@ -1117,13 +1256,13 @@ CONTAINS
ENDDO ENDDO
! !
CALL FLUSH(numout) CALL FLUSH(numout)
CALL ipslerr_p(3,'routing_reg_findbasins','We could not route point','','') CALL ipslerr_p(3,'routing_reg_localbas','We could not route point','','')
ENDIF ENDIF
! !
IF ( outtmp(il,jl) .NE. -1 ) THEN IF ( outtmp(il,jl) .NE. -1 ) THEN
ibas = outtmp(il,jl) ibas = outtmp(il,jl)
tsz(ibas) = tsz(ibas) + 1 tsz(ibas) = tsz(ibas) + 1
IF ( tsz(ibas) .GT. nbvmax ) CALL ipslerr_p(3,'routing_reg_findbasins','nbvmax too small','second section','') IF ( tsz(ibas) .GT. nbvmax ) CALL ipslerr_p(3,'routing_reg_localbas','nbvmax too small','second section','')
tpts(ibas, tsz(ibas), 1) = ip tpts(ibas, tsz(ibas), 1) = ip
tpts(ibas, tsz(ibas), 2) = jp tpts(ibas, tsz(ibas), 2) = jp
outtmp(ip,jp) = outtmp(il,jl) outtmp(ip,jp) = outtmp(il,jl)
...@@ -1151,7 +1290,7 @@ CONTAINS ...@@ -1151,7 +1290,7 @@ CONTAINS
! !
WRITE (numout,*) ' outtmp(il,jl):', outtmp(il,jl) WRITE (numout,*) ' outtmp(il,jl):', outtmp(il,jl)
WRITE (numout,*) ' trip(il,jl):', trip(il,jl) WRITE (numout,*) ' trip(il,jl):', trip(il,jl)
CALL ipslerr_p(3,'routing_reg_findbasins','We could not find right outlet','','') CALL ipslerr_p(3,'routing_reg_localbas','We could not find right outlet','','')
ENDIF ENDIF
! !
ENDIF ENDIF
...@@ -1163,10 +1302,10 @@ CONTAINS ...@@ -1163,10 +1302,10 @@ CONTAINS
IF (SUM(tsz) == COUNT(trip .LT. undef_int)) THEN IF (SUM(tsz) == COUNT(trip .LT. undef_int)) THEN
totsz = SUM(tsz)