Commit 5b030bc4 authored by Anthony Schrapffer's avatar Anthony Schrapffer
Browse files

Routing_reg with new linkup version

parent 9cb5c400
......@@ -26,7 +26,7 @@ MODULE routing_reg
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
!! ================================================================================================================================
!! SUBROUTINE : routing_reg_getgrid
......@@ -68,7 +68,7 @@ CONTAINS
!! REFERENCES : None
!!
!! FLOWCHART : None
!! \n
!! \n
!_ ================================================================================================================================
SUBROUTINE routing_reg_getgrid(nbpt, ib, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
......@@ -248,7 +248,7 @@ CONTAINS
ELSE
!
! It can also be a border point if the neighbour is not defined
!
!
IF ( basin_bx(ipp,jpp) > undef_int-1 ) THEN
IF (routing_nextgrid(ib,trip_bx(ip,jp)) < -1 ) THEN
trip_bx(ip,jp) = 98
......@@ -359,7 +359,7 @@ CONTAINS
IMPLICIT NONE
!
!! INPUT VARIABLES
INTEGER(i_std), INTENT(in) :: ib !!
INTEGER(i_std), INTENT(in) :: ib !!
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)
REAL(r_std), INTENT(in) :: hierarchy(:,:) !!
......@@ -391,7 +391,7 @@ CONTAINS
INTEGER(i_std), INTENT(out) :: coast_pts(nbvmax) !! The coastal flow points (unitless)
!
!! LOCAL VARIABLES
LOGICAL, PARAMETER :: debug=.FALSE.
LOGICAL, PARAMETER :: debug=.TRUE.
CHARACTER(LEN=7) :: fmt !!
CHARACTER(LEN=9) :: fmtr !!
!
......@@ -656,7 +656,7 @@ CONTAINS
ENDIF
!
! 4.0 Use the level 1 of the Pfafstetter coding system to divide the large
! sub-grid basin. The large sub-grid basin here is defined by the area
! sub-grid basin. The large sub-grid basin here is defined by the area
! larger than 9% of grid box area. This threshold should be discussed later.
!
! Make a loop through all sub-grid basins to find if there is big sub-grid basin.
......@@ -669,7 +669,7 @@ CONTAINS
!
DO WHILE ( COUNT(tsz(:)/REAL(totsz) >= maxpercent/REAL(100) ) > 0 .AND. COUNT(tsz(:) >= 4 ) > 0 )
!
cnt = 0
cnt = 0
DO ibas = 1, nbb
! Only take the subbasin is bigger than 4 points
IF (( tsz(ibas)/REAL(totsz) .GE. maxpercent/REAL(100) ) .AND. ( tsz(ibas) .GE. 4 )) THEN
......@@ -679,15 +679,15 @@ CONTAINS
ENDDO
!
! cnt should be greater or equal 1 here
DO ipb = 1, cnt
DO ipb = 1, cnt
!
ibas = trans(ipb)
!
CALL routing_reg_divbas(nbi, nbj, ibas, toutloc(ibas,1), toutloc(ibas,2), tsz(ibas), toutbas(ibas), toutdir(ibas), &
& toutlshead(ibas), tpts, trip, basin, fac, lontmp, lattmp, &
CALL routing_reg_divbas(nbi, nbj, ibas, toutloc(ibas,1), toutloc(ibas,2), tsz(ibas), toutbas(ibas), toutdir(ibas), &
& toutlshead(ibas), tpts, trip, basin, fac, lontmp, lattmp, &
& new_nb, mp, new_bname, new_outdir, new_heading, new_outbas, new_lat, new_lon, new_outloc, new_sz, new_pts)
!
! If the dividing step is ok
! If the dividing step is ok
!
IF ( new_nb .NE. 0 ) THEN
!
......@@ -713,7 +713,7 @@ CONTAINS
!
tsz(ibas) = new_sz(1)
tpts(ibas,:,:) = new_pts(1,:,:)
DO p = 1, tsz(ibas)
DO p = 1, tsz(ibas)
outtmp(tpts(ibas,p,1),tpts(ibas,p,2)) = ibas
ENDDO
!
......@@ -730,7 +730,7 @@ CONTAINS
!
IF ( new_outbas(ip) .EQ. 1 .AND. ip .NE. mp ) THEN
toutbas(nbb+ip-1) = ibas
ELSE IF ( ip .NE. mp ) THEN
ELSE IF ( ip .NE. mp ) THEN
toutbas(nbb+ip-1) = new_outbas(ip)+nbb-1
ELSE
toutbas(nbb+ip-1) = new_outbas(ip)
......@@ -744,7 +744,7 @@ CONTAINS
!
tsz(nbb+ip-1) = new_sz(ip)
tpts(nbb+ip-1,:,:) = new_pts(ip,:,:)
DO p = 1, tsz(nbb+ip-1)
DO p = 1, tsz(nbb+ip-1)
outtmp(tpts(nbb+ip-1,p,1),tpts(nbb+ip-1,p,2)) = nbb+ip-1
ENDDO
ENDDO
......@@ -760,7 +760,7 @@ CONTAINS
!
WRITE(numout,*) 'Can not divide sub-basins (routing_reg_findbasins): ',ibas
CALL ipslerr_p(3,'routing_reg_findbasins','new_nb = 0','Something is wrong with routing_reg_divbas','')
!
!
ENDIF
ENDDO
!
......@@ -880,7 +880,7 @@ CONTAINS
ENDIF
ENDDO
IF ( neworder .NE. 0 ) THEN
basin_bbout(ip) = neworder
basin_bbout(ip) = neworder
ELSE
WRITE(numout,*) 'ip, sortind(ip) :', ip, sortind(ip)
WRITE(numout,*) 'toutbas(sortind(ip)) :', toutbas(sortind(ip))
......@@ -895,10 +895,10 @@ CONTAINS
WRITE(numout,*) 'tbname(sortind(1:nb_basin)) :', tbname(sortind(1:nb_basin))
WRITE(numout,*) 'toutbas(1:nbb) :', toutbas(1:nbb)
!
CALL ipslerr_p(3,'routing_reg_findbasins','We got problem when sorting toutbas','','')
CALL ipslerr_p(3,'routing_reg_findbasins','We got problem when sorting toutbas','','')
ENDIF
ENDIF
ENDDO
ENDDO
!
IF ( debug .AND. routing_diagbox(ib, diaglalo) ) THEN
WRITE(numout,*) "Grid box: ", ib
......@@ -1029,7 +1029,7 @@ CONTAINS
!! INPUT VARIABLES
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) :: ibas !! Order of the basin will be divided
INTEGER(i_std), INTENT(in) :: ibas !! Order of the basin will be divided
INTEGER(i_std), INTENT(in) :: iloc, jloc !! Outlet location
!
INTEGER(i_std), INTENT(in) :: tsz !! Size of sub-basin
......@@ -1059,7 +1059,7 @@ CONTAINS
!
!! LOCAL VARIABLES
REAL(r_std), PARAMETER :: flag=-9999. !!
LOGICAL, PARAMETER :: debug=.FALSE.
LOGICAL, PARAMETER :: debug=.TRUE.
LOGICAL, PARAMETER :: checkgrid=.FALSE.
CHARACTER(LEN=7) :: fmt !!
CHARACTER(LEN=9) :: fmtr !!
......@@ -1107,19 +1107,19 @@ CONTAINS
!
INTEGER(i_std) :: main_loc(nbvmax,2), tri_loc(nbne,nbvmax,2)
REAL(r_std) :: main_fac(nbvmax), tri_fac(nbne,nbvmax)
REAL(r_std) :: tmptri_fac(nbne) ! Flow accumulation of all neighbour points
REAL(r_std) :: tmptri_fac(nbne) ! Flow accumulation of all neighbour points
INTEGER(i_std) :: sortedtrifac(nbne) ! And sort of tmptri_fac(nbne)
REAL(r_std) :: alltri_fac(nbvmax) ! Sort flow accumulation of all tributaries
INTEGER(i_std) :: alltri_loc1(nbvmax), alltri_loc2(nbvmax) ! And their location
INTEGER(i_std) :: sortedallfac(nbvmax) ! Sort alltri_fac(nbvmax)
!
INTEGER(i_std) :: sortedallfac(nbvmax) ! Sort alltri_fac(nbvmax)
!
! * Note: again, there are few dimensions should be exactly [nbne-1]
! or [nbne-2]. Because when you follow upstream the mainstream, there are
! always at least ONE neighbour points will not belong to TRIBUTARY.
! But to make your life easier, I put all "nbne" here.
!
! - Identify the 4 greatest tributaries:
!
!
REAL(r_std) :: tmpmain_fac(4) ! Flow accumulations of 4 greatest tributaries
INTEGER(i_std) :: tmpmain_loc(4,2), sortedmainfac(4) ! And their locations and sorted values
INTEGER(i_std) :: numtri(4), usetri_loc(4,4,2) ! Counting and store location when processing these 4 points
......@@ -1128,7 +1128,7 @@ CONTAINS
! stem: there can be 5 inter-basins. So we need 9 dividing points:
!
INTEGER(i_std) :: divloc(9,2) ! Location of dividing points (maximum is 9)
!
!
!_ ================================================================================================================================
......@@ -1178,7 +1178,7 @@ CONTAINS
! Print out information of grid box
!
IF ( checkgrid ) THEN
!
!
WRITE(numout,*) " Routing_reg_divbas: Grid before dividing "
WRITE(numout,*) " Size: ", nbi, nbj
WRITE(fmt,"('(',I3,'I6)')") nbi
......@@ -1225,7 +1225,7 @@ CONTAINS
WRITE(numout,*) 'Basin out: ', tout, toutd
ENDIF
!
! 2.0 From the outlet of the subbasin, trace back to find the main stream
! 2.0 From the outlet of the subbasin, trace back to find the main stream
! and identify the 4 greatest tributaries. If there are less than 4
! tributraries, take them all.
!
......@@ -1269,7 +1269,7 @@ CONTAINS
IF ( ie .EQ. il .AND. je .EQ. jl ) THEN
! Store all points here
l = l + 1
tri_fac(l,cnt) = factmp(ill,jll)
tri_fac(l,cnt) = factmp(ill,jll)
tri_loc(l,cnt,1) = ill
tri_loc(l,cnt,2) = jll
! Mark the processed points here
......@@ -1298,7 +1298,7 @@ CONTAINS
! accumulation will belongs to main stream. We move to this point
! (change value of il, jl).
IF ( tri_fac(sortedtrifac(1),cnt) .NE. flag ) THEN
il = tri_loc(sortedtrifac(1),cnt,1)
il = tri_loc(sortedtrifac(1),cnt,1)
jl = tri_loc(sortedtrifac(1),cnt,2)
!
IF ( checkgrid ) THEN
......@@ -1353,7 +1353,7 @@ CONTAINS
ENDDO
IF ( debug ) WRITE (numout,*) 'Length of mainstream: ', cnt
!
! 3.0 Sort the flow accumulation of the tributaries and find the
! 3.0 Sort the flow accumulation of the tributaries and find the
! divided location.
!
alltri_fac(:) = flag
......@@ -1366,7 +1366,7 @@ CONTAINS
! Actually, we should only DO l = 1, nbne-2
! Because, there are usually 2 neighbour points belong to main stream
! and maximum (nbne-2) neighbour points can be tributaries.
! DO l = 1, nbne-2
! DO l = 1, nbne-2
DO l = 1, nbne
IF ( tri_fac(l,k) .NE. flag ) THEN
ip = ip + 1
......@@ -1383,7 +1383,7 @@ CONTAINS
IF ( ip .GT. 0 ) THEN
! Original output
oic = 0
!
!
sortedallfac(:) = 0
DO k = 1, ip
ff = MAXLOC(alltri_fac)
......@@ -1414,14 +1414,14 @@ CONTAINS
tmpmain_fac(l) = main_fac(ik)
!
numtri(l) = 1
usetri_loc(numtri(l),l,1) = tri_loc(il,ik,1)
usetri_loc(numtri(l),l,1) = tri_loc(il,ik,1)
usetri_loc(numtri(l),l,2) = tri_loc(il,ik,2)
!
ELSE
! If this point in main stream have more than 1 tributary
! belongs to top 4 greatest tributaries
numtri(l) = numtri(l) + 1
usetri_loc(numtri(l),l,1) = tri_loc(il,ik,1)
usetri_loc(numtri(l),l,1) = tri_loc(il,ik,1)
usetri_loc(numtri(l),l,2) = tri_loc(il,ik,2)
!
ENDIF
......@@ -1455,12 +1455,12 @@ CONTAINS
!
new_outbas(ic) = tout
! Number of mainstream sub-basin
oic = ic
oic = ic
!
DO k = l, 1, -1
!
!DO il = 1, nbne-2
! Maximum 4 tributaries:
! Maximum 4 tributaries:
DO il = 1, 4
IF ( usetri_loc(il,sortedmainfac(k),1) .NE. 0 ) THEN
ic = ic + 1
......@@ -1475,7 +1475,7 @@ CONTAINS
ENDDO
! Save number of dividing points
l = ic
!
!
! If we don't have any tributary
ELSE
oic = 2
......@@ -1502,7 +1502,7 @@ CONTAINS
k = ik
ELSE
IF ( ik .LE. (l - oic) ) THEN
k = ik + oic
k = ik + oic
ELSE
k = ik - (l - oic)
ENDIF
......@@ -1510,8 +1510,8 @@ CONTAINS
! I'm so stupid careful here with this IF
IF ( divloc(k,1) .NE. 0 .AND. divloc(k,2) .NE. 0 ) THEN
!
new_sz(k) = 1
new_pts(k,new_sz(k),1) = divloc(k,1)
new_sz(k) = 1
new_pts(k,new_sz(k),1) = divloc(k,1)
new_pts(k,new_sz(k),2) = divloc(k,2)
!
new_bname(k) = basin(divloc(k,1),divloc(k,2))
......@@ -1535,7 +1535,7 @@ CONTAINS
!
!
triptemp(divloc(k,1), divloc(k,2)) = -1
!
!
DO isz = 1, tsz
il = tpts(ibas, isz, 1)
jl = tpts(ibas, isz, 2)
......@@ -1552,16 +1552,16 @@ CONTAINS
p = trip(ie,je)
IF ( ie .EQ. divloc(k,1) .AND. je .EQ. divloc(k,2) ) THEN
okpoint = .FALSE.
new_sz(k) = new_sz(k) + 1
new_pts(k,new_sz(k),1) = il
new_sz(k) = new_sz(k) + 1
new_pts(k,new_sz(k),1) = il
new_pts(k,new_sz(k),2) = jl
triptemp(tpts(ibas, isz, 1), tpts(ibas, isz, 2)) = -1
ENDIF
ENDDO
ENDIF
ENDDO
ENDDO
ENDIF
ENDDO
ENDDO
!
! 4.0 Make simple verification
!
......@@ -1570,7 +1570,7 @@ CONTAINS
DO ip=1,new_nb
checksz = checksz + new_sz(ip)
ENDDO
!
!
IF ( checksz .NE. tsz) THEN
WRITE(numout,*) ' Water got lost while I tried to divide basins'
WRITE(numout,*) ' Number of new sub-basin : ', new_nb
......@@ -1631,7 +1631,7 @@ CONTAINS
!!
!! REFERENCES : None
!!
!! FLOWCHART : None
!! FLOWCHART : None
!! \n
!_ ================================================================================================================================
......@@ -1660,7 +1660,7 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx,
INTEGER(i_std) :: nb_basin !! Number of sub-basins (unitless)
INTEGER(i_std), DIMENSION(nbvmax) :: basin_inbxid, basin_sz !! ID of basin, number of points in the basin
REAL(r_std), DIMENSION(nbvmax,2) :: basin_outlet !! Outlet coordinate of each subgrid basin (lat,lon)
REAL(r_std), DIMENSION(nbvmax) :: basin_outtp !!
REAL(r_std), DIMENSION(nbvmax) :: basin_outtp !!
INTEGER(i_std), DIMENSION(nbvmax,nbvmax,2) :: basin_pts !! Points in each basin
INTEGER(i_std), DIMENSION(nbvmax) :: basin_bxout !! outflow direction
INTEGER(i_std), DIMENSION(nbvmax) :: basin_bbout !! outflow sub-basin
......@@ -1681,7 +1681,7 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx,
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_topoind !! Topographic index of the residence time for a basin (m)
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_lshead !! Large scale heading out of the grid box.
INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_grid !! Type of outflow on the grid box (unitless)
INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_basin !!
INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_basin !!
INTEGER(i_std), DIMENSION(nbpt) :: nbcoastal !!
INTEGER(i_std), DIMENSION(nbpt,nwbas) :: coastal_basin !!
!
......@@ -1917,7 +1917,7 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx,
!! For outflow_basin we have the following conventions :
!! outflow_basin = basin id in the case of the basin not flowing out of the current grid (i.e. outflow_grid=-4).
!! Else outflow_basin = undef_int and the objective of this routine is to find what that basin is. This work
!! essentially occurs in the routine routing_reg_bestsubbasin. But for that we need to find the grid box where
!! essentially occurs in the routine routing_reg_bestsubbasin. But for that we need to find the grid box where
!! we will look for the right basin. This operation is performed here in 5 successive steps. The order correspond
!! to the solution we would prefer.
!! 1.0 : We will look in the grid provided by outflow_grid for the right basin. This neighbour has been obtained
......@@ -1985,12 +1985,16 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
REAL(r_std), ALLOCATABLE, DIMENSION(:) :: hatoutflow
INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: minhloc
!
INTEGER(i_std) , DIMENSION(NbNeighb) :: order_ref
INTEGER(i_std) :: nb_sym, nb, val, ind_nb, onb
INTEGER(i_std) :: found
!
!! PARAMETERS
LOGICAL, PARAMETER :: debug = .FALSE. !! (true/false)
LOGICAL, PARAMETER :: debug = .TRUE. !! (true/false)
!
!_ ================================================================================================================================
!
!
!
testbasinid = 195
!
......@@ -2036,6 +2040,11 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
! At this point any of the outflows is designated by a negative values in
! outflow_grid
!
WRITE(numout,*) "*****"
WRITE(numout,*) "Linkup 0 - sp, sb = ", sp, sb
WRITE(numout,*) "Linkup 0 - outflow_grid = ", outflow_grid(sp,sb)
IF ( outflow_grid(sp,sb) == 0 ) WRITE(numout,*) "Linkup 1.0 -- Flow out of Halo zone"
!
IF ( outflow_grid(sp,sb) < 0 ) THEN
IF ( outflow_grid(sp,sb) .EQ. -4 ) THEN
! Flow into a basin of the same grid
......@@ -2059,620 +2068,263 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
! Nothing to do but just remember it is done.
solved(sp,1) = solved(sp,1) + 1
ENDIF
WRITE(numout,*) sp,sb,"is an outflow ", outflow_grid(sp,sb)
ELSE IF ( outflow_grid(sp,sb) .GT. 0 ) THEN
found = 0
!
! Deal with the first one as usual
!
inp = outflow_grid(sp,sb)
IF ( inp == 0.) WRITE(numout,*) "ERROR 1"
!
CALL routing_reg_bestsubbasin(sp, sb, basin_id(sp,sb), basin_hierarchy(sp,sb), basin_fac(sp,sb), &
& basin_flowdir(sp,sb), invented_basins, &
& nbpt, nwbas, inp, basin_count, basin_id, basin_hierarchy, basin_fac, &
& basin_flowdir, nbcoastal, coastal_basin, bop, bopqual)
!
! Trung: Because I think it's for nothing, I commented this line:
!ang = routing_anglediff_g(sp, inp, basin_flowdir(sp,sb))
!
! If the basin is suitable (bop < undef_int) then take it
!
IF ( bop .LT. undef_int ) THEN
!
CALL routing_updateflow(sp, sb, inp, bop, nbpt, nwbas, nbxmax*8, outflow_grid, outflow_basin, &
& inflow_number, inflow_grid, inflow_basin)
IF ( outflow_basin(sp,sb) == bop ) THEN
solved(sp,1) = solved(sp,1) + 1
found = 1
WRITE(numout,*) "Solution found in the original outflow_grid"
ENDIF
!
ENDIF
!
ENDIF
!
IF ( debug .AND. routing_diagbox_g(sp,diaglalo) ) THEN
IF ( outflow_grid(sp,sb) == 0 ) WRITE(numout,*) "Linkup 1.0 -- Flow out of Halo zone"
WRITE(numout,*) "Linkup 1.0 - In routing_reg_linkup working on sp & sb:",sp, sb
WRITE(numout,*) "Linkup 1.0 - Coords : ", lalo_g(sp,2), lalo_g(sp,1)
WRITE(numout,*) "Linkup 1.0 - outflow_flowdir = ", basin_flowdir(sp,sb)
WRITE(numout,*) "Linkup 1.0 - Large scale heading =", basin_lshead(sp,sb)
WRITE(numout,*) "Linkup 1.0 - Hierarchy =", basin_hierarchy(sp,sb)
WRITE(numout,*) "Linkup 1.0 - Basin % of grid =", basin_area(sp,sb)/area_g(sp)*100
WRITE(numout,*) "Linkup 1.0 - outflow_grid =", outflow_grid(sp,sb)
WRITE(numout,*) "Linkup 1.0 - ID = ", basin_id(sp,sb)
IF ( outflow_grid(sp,sb) > 0 ) THEN
WRITE(numout,*) "Linkup 1.0 - Coords outflow: ", lalo_g(outflow_grid(sp,sb),2), lalo_g(outflow_grid(sp,sb),1)
ENDIF
WRITE(numout,*) "Linkup 1.0 - outflow_basin = ", outflow_basin(sp,sb)
IF ( outflow_grid(sp,sb) > 0 .AND. outflow_basin(sp,sb) < undef_int ) THEN
WRITE(numout,*) "Linkup 1.0 - Outflow Hierarchy =", basin_hierarchy(outflow_grid(sp,sb),outflow_basin(sp,sb))
ENDIF
ENDIF
ENDDO
ENDDO
!
! 2.0 : Do we have a valid outflow basin ? If not we need to test other options.
! The reason for not finding the outflow basin could be that the outflow_grid correspons to a
! small scale feature not corresponding to the orientation of the grid. So we use the large
! scale heading of the flow out of the grid in order to test another neighbour.
!
DO sp=1,nbpt
!
DO sb=1,basin_count(sp)
!
IF ( outflow_basin(sp,sb) .EQ. undef_int .AND. &
& outflow_grid(sp,sb) > 0 .AND. basin_flowdir(sp,sb) > 0) THEN
!
DO i=1,NbNeighb
diffangle(i) = ABS(haversine_diffangle(headings_g(sp,i), basin_lshead(sp,sb)))
ENDDO
ff = MINLOC(diffangle)
dls = neighbours_g(sp,ff(1))
!
IF ( dls .LE. -1 ) THEN
! This flows out of the domain so put the basin into coastal flow
outflow_grid(sp,sb) = -2
outflow_basin(sp,sb) = undef_int
solved(sp,2) = solved(sp,2)+1
IF ( debug .AND. routing_diagbox_g(sp,diaglalo) ) THEN
WRITE(numout,*) "Linkup 2.0 - ", sp, sb, "becomes a coastal flow basin."
ENDIF
ELSE
IF ( dls > 0 ) THEN
CALL routing_reg_bestsubbasin(sp, sb, basin_id(sp,sb), basin_hierarchy(sp,sb), &
& basin_fac(sp,sb), basin_flowdir(sp,sb), invented_basins, &
& nbpt, nwbas, dls, basin_count, basin_id, basin_hierarchy, &
& basin_fac, basin_flowdir, nbcoastal, coastal_basin, bls, blsqual)
ELSE
bls = undef_int
blsqual = 0
ENDIF
ENDIF
!
IF ( dls > 0 .AND. bls < undef_int ) THEN
IF ( found == 0 ) THEN
WRITE (numout,*) "Establishing the list of neighbours"
! Organize the location of the neighbours to visit by order of priority
! first the outflow grid then 2 by 2 till the opposite side (by +1/-1 - +2/-2 ...)
! if NbNeighb is odd we have to had the opposite neighbour
! ex : if NbNeighb = 8 and outflow grid is at the bottom right
!
CALL routing_updateflow(sp, sb, dls, bls, nbpt, nwbas, nbxmax*8, outflow_grid, outflow_basin, &
& inflow_number, inflow_grid, inflow_basin)
! 8 | 7 | 5
! ---|---|---
! 6 | x | 3
! ---|---|---
! 4 | 2 | 1
!
ENDIF
!
IF ( debug .AND. routing_diagbox_g(sp,diaglalo) ) THEN
WRITE(numout,*) "Linkup 2.0 - sp, sb = ", sp, sb
WRITE(numout,*) "Linkup 2.0 - Large scale heading =", basin_lshead(sp,sb)
WRITE(numout,*) "Linkup 2.0 - Large scale test =", dls, bls
WRITE(numout,*) "Linkup 2.0 - outflow_grid =", outflow_grid(sp,sb)
WRITE(numout,*) "Linkup 2.0 - outflow_basin = ", outflow_basin(sp,sb)
ENDIF
!
IF ( outflow_basin(sp,sb) < undef_int ) solved(sp,2) = solved(sp,2)+1
!
ENDIF
!
ENDDO
ENDDO
!
!
! 3.0 Look within the general direction of the indicated flow to see if we find a matching basin.
!
! In case we did not find the correct basin we start to look
! around. We find two options for the outflow direction (dp1 & dm1) and the corresponding
! basin index (bp1 & bm1).
! These options are only acceptable if the angle to the original direction is less or equal
! to 45 degrees.
!
DO sp=1,nbpt
!
DO sb=1,basin_count(sp)
!
nbocean = 0
!
IF ( outflow_basin(sp,sb) .EQ. undef_int .AND. &
& outflow_grid(sp,sb) > 0 .AND. basin_flowdir(sp,sb) > 0) THEN
!
! Try to neighboring directions on the model grid.
!
dp1 = routing_nextgrid_g(sp, basin_flowdir(sp,sb), +1)
angp1 = routing_anglediff_g(sp, dp1, basin_flowdir(sp,sb))
dm1 = routing_nextgrid_g(sp, basin_flowdir(sp,sb), -1)
angm1 = routing_anglediff_g(sp, dm1, basin_flowdir(sp,sb))
!
! How many are ocean points ?
!
IF ( dp1 .LE. -1 ) nbocean = nbocean + 1
IF ( dm1 .LE. -1 ) nbocean = nbocean + 1
!
! Compute bp1 and bm1
!
IF ( dp1 .GT. 0 ) THEN
IF ( angp1 <= 45.0 ) THEN
CALL routing_reg_bestsubbasin(sp, sb, basin_id(sp,sb), basin_hierarchy(sp,sb), basin_fac(sp,sb), &
& basin_flowdir(sp,sb), invented_basins, &
& nbpt, nwbas, dp1, basin_count, basin_id, basin_hierarchy, basin_fac, &
& basin_flowdir, nbcoastal, coastal_basin, bp1, bp1qual)
ELSE
bp1 = undef_int
bp1qual = 0
ENDIF
ELSE
bp1 = undef_int
bp1qual = 0
ENDIF