diff --git a/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index 4c7b94445b23ff82c23891f4aeff41980e291c4b..215e5346e2ada071862cd71f45508fb910bba9c1 100644
--- a/F90subroutines/routing_reg.f90
+++ b/F90subroutines/routing_reg.f90
@@ -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
- IF ( dm1 .GT. 0 ) THEN
- IF ( angm1 <= 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, dm1, basin_count, basin_id, basin_hierarchy, basin_fac, &
- & basin_flowdir, nbcoastal, coastal_basin, bm1, bm1qual)
- ELSE
- bm1 = undef_int
- bm1qual = 0
- ENDIF
- ELSE
- bm1 = undef_int
- bm1qual = 0
- ENDIF
- !
- ! Decide between dp1 and dm1 which is the better solution !
- !
- dop = undef_int
- IF ( bp1 < undef_int .OR. bm1 < undef_int ) THEN
- IF ( bp1 < undef_int ) THEN
- dop = dp1
- bop = bp1
- ELSE IF ( bm1 < undef_int ) THEN
- dop = dm1
- bop = bm1
- ELSE
- IF ( bp1qual > bm1qual ) THEN
- dop = dp1
- bop = bp1
- ELSE IF ( bp1qual > bm1qual ) THEN
- dop = dm1
- bop = bm1
- ELSE
- IF ( angm1 < angp1 ) THEN
- dop = dm1
- bop = bm1
- ELSE
- dop = dp1
- bop = bp1
- ENDIF
- ENDIF
+ order_ref(:) = -1
+ inp = outflow_grid(sp,sb)
+ ! The first one is the location of the outflow grid
+ order_ref(1) = minloc(abs(neighbours_g(sp,:) - inp), 1)
+ ! Then we have to adjust in function of NbNeighb to add the location 2 by 2
+ nb_sym = (NbNeighb-1)/2
+ ind_nb = 2
+ DO nb=1,nb_sym
+ DO onb=-1,1,2
+ val = modulo(order_ref(1)+nb*onb,NbNeighb)
+ IF ( val == 0) val = NbNeighb
+ order_ref(ind_nb) = val
+ ind_nb = ind_nb + 1
+ ENDDO
+ ENDDO
+ ! If NbNeighb is odd
+ IF ( modulo(NbNeighb,2) == 0 ) THEN
+ val = modulo(order_ref(1)+4,NbNeighb)
+ IF ( val == 0) val = NbNeighb
+ order_ref(NbNeighb) = val
ENDIF
- ENDIF
- !
- ! Now that we know where we want the water to flow to we write the
- ! the information in the right fields.
- !
- IF ( dop > 0 .AND. dop < undef_int .AND. bop < undef_int ) THEN
!
- CALL routing_updateflow(sp, sb, dop, bop, nbpt, nwbas, nbxmax*8, outflow_grid, outflow_basin, &
- & inflow_number, inflow_grid, inflow_basin)
!
- ENDIF
- !
- ! If one of the outflow options is an ocean point then we should use it.
- ! We are a little more stringent with the angle.
- !
- IF ( outflow_basin(sp,sb) .EQ. undef_int .AND. &
- & outflow_grid(sp,sb) > 0 .AND. basin_flowdir(sp,sb) > 0) THEN
- IF ( nbocean >= 1 ) THEN
- IF ( (dp1 .LE. -1 .AND. angp1 < 45.0 ) .OR. &
- & (dm1 .LE. -1 .AND. angm1 < 45.0 ) ) THEN
- outflow_grid(sp,sb) = -2
- outflow_basin(sp,sb) = undef_int
- solved(sp,3) = solved(sp,3)+1
- IF ( debug .AND. routing_diagbox_g(sp,diaglalo) ) THEN
- WRITE(numout,*) "Linkup 3.0 - ", sp, sb, "becomes a coastal flow basin."
- ENDIF
- ENDIF
- ENDIF
- ENDIF
- !
- IF ( debug .AND. routing_diagbox_g(sp,diaglalo) ) THEN
- WRITE(numout,*) "Linkup 3.0 - In routing_reg_linkup working on sp & sb:", sp, sb
- WRITE(numout,*) "Linkup 3.0 - outflow_flowdir = ", basin_flowdir(sp,sb)
- WRITE(numout,*) "Linkup 3.0 - Tested dp1 & dm1 : ", dp1, dm1
- WRITE(numout,*) "Linkup 3.0 - Angles : ",angp1, angm1
- WRITE(numout,*) "Linkup 3.0 - Test output bp1 & bm1 = ", bp1, bm1
- WRITE(numout,*) "Linkup 3.0 - Selected bop & dop = ", bop, dop
- WRITE(numout,*) "Linkup 3.0 - outflow_grid =", outflow_grid(sp,sb)
- WRITE(numout,*) "Linkup 3.0 - outflow_basin = ", outflow_basin(sp,sb)
- ENDIF
- !
- IF ( outflow_basin(sp,sb) < undef_int ) solved(sp,3) = solved(sp,3)+1
- !
- ENDIF
- !
- ENDDO
- ENDDO
- !
- ! 4.0 If we still have not found the correct basin so we look into other directions
- ! East and West of the large scale flow direction we have. We work this with the
- ! the smallest angle to the large scale flow direction first. Remember the grid
- ! in the large scale flow direction was already tested (in 2.0).
- !
- nbtotest = INT(NbNeighb/2.0)
- !
- DO sp=1,nbpt
- !
- DO sb=1,basin_count(sp)
- !
- nbocean = 0
- wocean = zero
- !
- 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.
- !
- DO i=1,NbNeighb
- diffangle(i) = ABS(haversine_diffangle(headings_g(sp,i), basin_lshead(sp,sb)))
- ENDDO
- !
- ! The minimum was already done above. So we delete from the list
- !
- ff = MINLOC(diffangle)
- dls = neighbours_g(sp,ff(1))
- IF ( dls .LE. -1 ) THEN
- nbocean = nbocean + 1
- wocean = wocean + diffangle(ff(1))
- ENDIF
- diffangle(ff(1)) = undef_sechiba
- !
- DO it=1,nbtotest
- ff = MINLOC(diffangle)
- gridstotest(it) = neighbours_g(sp,ff(1))
- gridangle(it) = diffangle(ff(1))
+ ! If not found, now test the next points 2 by 2
!
- ! neighbours outside of the domain are considered as ocean points.
- !
- IF ( gridstotest(it) .LE. -1 ) THEN
- nbocean = nbocean + 1
- wocean = wocean + diffangle(ff(1))
- ENDIF
- diffangle(ff(1)) = undef_sechiba
- ENDDO
- !
- ! Now work through our list of points
- !
- DO it=1,nbtotest
- IF ( outflow_basin(sp,sb) .EQ. undef_int ) THEN
- dls = gridstotest(it)
- !
- IF ( dls .GT. 0 ) THEN
- dop = dls
- 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, bop, bopqual)
- ELSE
- bop = undef_int
- bopqual = 0
- ENDIF
- gridbasin(it) = bop
- !
- ! Now that we know where we want the water to flow to we write the
- ! the information into the right fields.
- !
- IF ( dls > 0 .AND. dls < undef_int .AND. bop < undef_int ) THEN
+ IF ( found == 0) THEN
+ nb = 1
+ DO WHILE (( found == 0 ) .AND. (nb .LE. nb_sym))
!
- CALL routing_updateflow(sp, sb, dls, bop, nbpt, nwbas, nbxmax*8, outflow_grid, outflow_basin, &
- & inflow_number, inflow_grid, inflow_basin)
- !
- ENDIF
- ENDIF
- ENDDO
- !
- IF ( debug .AND. routing_diagbox_g(sp,diaglalo) ) THEN
- WRITE(numout,*) "Linkup 4.0 - In routing_reg_linkup working on sp & sb:", sp, sb
- WRITE(numout,*) "Linkup 4.0 - Large scale heading =", basin_lshead(sp,sb)
- WRITE(numout,*) "Linkup 4.0 - Tested gridstotest : ", gridstotest(1:nbtotest)
- WRITE(numout,*) "Linkup 4.0 - Angles of tested : ", gridangle(1:nbtotest)
- WRITE(numout,*) "Linkup 4.0 - Tested outflow basin = ", gridbasin(1:nbtotest)
- WRITE(numout,*) "Linkup 4.0 - outflow_grid =", outflow_grid(sp,sb)
- WRITE(numout,*) "Linkup 4.0 - outflow_basin = ", outflow_basin(sp,sb)
- WRITE(numout,*) "Linkup 4.0 - Number of ocean options : ", nbocean, nbtotest
- IF ( basin_id(sp,sb) < invented_basins ) THEN
- WRITE(numout,*) "Linkup 4.0 - Hmin for basin :", hatoutflow(basin_id(sp,sb))
- WRITE(numout,*) "Linkup 4.0 - Hmin location :", minhloc(basin_id(sp,sb),1), minhloc(basin_id(sp,sb),2)
- ENDIF
- ENDIF
- !
- IF ( outflow_basin(sp,sb) < undef_int ) solved(sp,4) = solved(sp,4)+1
- !
- ENDIF
- !
- ! If we have not found anything, but we have an ocean point as neighbour, then go for it.
- ! This would be coastal flow then.
- !
- IF ( outflow_basin(sp,sb) .EQ. undef_int .AND. &
- & outflow_grid(sp,sb) > 0 .AND. basin_flowdir(sp,sb) > 0) THEN
- !
- IF ( basin_id(sp,sb) < invented_basins ) THEN
- crit = (basin_hierarchy(sp,sb)-hatoutflow(basin_id(sp,sb)))/(basin_hierarchy(sp,sb)+1)
- ELSE
- crit = zero
- ENDIF
- !
- ! If we have at least 2 ocean points within the examined points then ouflow points
- ! we put this basin into the ocean. If hierarchy is lower than the
- ! average at other outflow points, then we relax this condition.
- !
- IF ( crit > 0.1 ) THEN
- IF ( nbocean >= nbtotest/2.0 ) THEN
- !
- ! These ocean points should not be too far away from the large scale direction.
- ! In this phase we should also consider the size of the basin. If it small then we can
- ! neglect its and put it into the ocean.
- !
- IF ( wocean/nbocean < 90.0 .OR. (basin_area(sp,sb)/area_g(sp)*100 < 1.0) ) THEN
- outflow_grid(sp,sb) = -2
- outflow_basin(sp,sb) = undef_int
- solved(sp,4) = solved(sp,4)+1
- IF ( debug .AND. routing_diagbox_g(sp,diaglalo) ) THEN
- WRITE(numout,*) "Linkup 4.1 - In routing_reg_linkup working on sp & sb:", sp, sb
- WRITE(numout,*) "Linkup 4.1 - basin_id = ", basin_id(sp,sb)
- WRITE(numout,*) "Linkup 4.1 - hierarchy = ", basin_hierarchy(sp,sb)
- WRITE(numout,*) "Linkup 4.1 - crit =", crit
- WRITE(numout,*) "Linkup 4.1 - ocean direction = ", wocean/nbocean
- WRITE(numout,*) "Linkup 4.1 - other criteria : ", basin_area(sp,sb)/area_g(sp)*100
+ ! From the grid point + and - get bp and bpqual
+ ! First point : dpi1 -> bp1 et bp1qual
+ dp1 = neighbours_g(sp,order_ref(nb*2))
+ ! we know it's wrong but it will serve to make decision later
+ angp1 = routing_anglediff_g(sp, dp1, basin_flowdir(sp,sb))
+ ! Check if grid point
+ IF (dp1 .GT. 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
- ENDIF
- ENDIF
- ELSE
- !
- ! Here we are with a basin that has a hierarchy close to that of the outflow point. So
- ! we can be more relaxed about the criteria.
- !
- IF ( nbocean >= 1 ) THEN
- !
- ! These ocean points should not be too far away from the large scale direction.
- !
- IF ( wocean/nbocean < 180.0 .OR. (basin_area(sp,sb)/area_g(sp)*100 < 25.0) ) THEN
- outflow_grid(sp,sb) = -2
- outflow_basin(sp,sb) = undef_int
- solved(sp,4) = solved(sp,4)+1
- IF ( debug .AND. routing_diagbox_g(sp,diaglalo) ) THEN
- WRITE(numout,*) "Linkup 4.1 - ", sp, sb, "becomes a coastal flow basin. Relaxed condition"
+ !
+ ! Second point
+ ! dm1 -> bm1 et bm1qual
+ dm1 = neighbours_g(sp,order_ref(nb*2+1))
+ ! we know it's wrong but it will serve to make decision later
+ angm1 = routing_anglediff_g(sp, dm1, basin_flowdir(sp,sb))
+ IF (dm1 .GT. 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, dm1, basin_count, basin_id, basin_hierarchy, basin_fac, &
+ & basin_flowdir, nbcoastal, coastal_basin, bm1, bm1qual)
+ ELSE
+ bm1 = undef_int
+ bm1qual = 0
ENDIF
- ENDIF
- ENDIF
- ENDIF
- !
- ENDIF
- !
- !
- ENDDO
- ENDDO
- !
- ! 5.0
- !
- ! We probably have not yet solved all points. So we go through all points and basins again
- ! to see if we cannot find some local solution as a last resort.
- !
- IF ( debug .AND. routing_diagbox_g(sp,diaglalo) ) THEN
- WRITE(numout,*) "Linkup 5.0 - In routing_reg_linkup working on sp & sb:", sp, sb
- WRITE(numout,*) "Linkup 5.0 - outflow_basin :", outflow_basin(sp,sb)
- WRITE(numout,*) "Linkup 5.0 - outflow_grid :", outflow_grid(sp,sb)
- WRITE(numout,*) "Linkup 5.0 - basin_flowdir :", basin_flowdir(sp,sb)
- ENDIF
- !
- 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
- !
- ! Last resort, we look if a larger basin (i.e. we already have solved !) within this grid has
- ! a lower or the same hierarchy (within 1%).
- !
- sbint = undef_int
- !
- DO sba=1,basin_count(sp)
- IF ( sba .NE. sb ) THEN
- IF ( basin_id(sp,sb) == basin_id(sp,sba) ) THEN
!
- ! Only look for suitable hierrachy if we have a solution for the sub-basin : it has an outflow grid
- ! and basin or flows into the ocean.
+ ! Decide between dp1 and dm1 which is the better solution
!
- IF ( (outflow_grid(sp,sba) > 0 .AND. outflow_basin(sp,sba) < undef_int) .OR. &
- & outflow_grid(sp,sba) .EQ. -1 .OR. outflow_grid(sp,sba) .EQ. -2 ) THEN
- !
- ! calculation of the criterion for considering that the hierarchy are close enough.
- ! Watch out hierarchy can be zero. 1 is considered the
- ! minimum here.
+ dop = undef_int
+ IF ( bp1 < undef_int .OR. bm1 < undef_int ) THEN
+ IF ( bp1 < undef_int ) THEN
+ dop = dp1
+ bop = bp1
+ ELSE IF ( bm1 < undef_int ) THEN
+ dop = dm1
+ bop = bm1
+ ELSE
+ IF ( bp1qual > bm1qual ) THEN
+ dop = dp1
+ bop = bp1
+ ELSE IF ( bp1qual > bm1qual ) THEN
+ dop = dm1
+ bop = bm1
+ ELSE
+ ! ATTENTION !!
+ ! In the case both point are the same quality
+ ! we need of another factor of decision
+ ! for now we use the angles even if we know it's wrong
+ IF ( angm1 < angp1 ) THEN
+ dop = dm1
+ bop = bm1
+ ELSE
+ dop = dp1
+ bop = bp1
+ ENDIF
+ ENDIF
+ ENDIF
+ ENDIF
+
+ IF ( dop > 0 .AND. dop < undef_int .AND. bop < undef_int ) THEN
!
- ! Trung: If current sub-basin (sp,sb) is the closest
- ! point to the outlet of river, the value of
- ! basin_hierarchy(sp,sb) is always lower than
- ! basin_hierarchy(sp,sba). It means that crit
- ! is negative and statement IF (crit < 0.1) is
- ! not useful. => I changed it a little bit.
+ CALL routing_updateflow(sp, sb, dop, bop, nbpt, nwbas, nbxmax*8, outflow_grid, outflow_basin, &
+ & inflow_number, inflow_grid, inflow_basin)
!
- crit = (basin_hierarchy(sp,sb)-basin_hierarchy(sp,sba))/(basin_hierarchy(sp,sb)+1)
- IF ( crit < 0.01 .AND. crit >= zero ) THEN
- ! Trung: Put condition of flow accumulation here to
- ! avoid error of circulation in routing_reg_fetch
- IF (basin_fac(sp,sba) .GE. basin_fac(sp,sb)) THEN
- sbint = sba
- ENDIF
+ IF ( outflow_basin(sp,sb) == bop ) THEN
+ solved(sp,2) = solved(sp,2) + 1
+ found = 1
+ WRITE (numout,*) "Neighbours, output found at:" , nb, "level"
ENDIF
+ !
+ ELSE
+ nb = nb+1
ENDIF
- ENDIF
+ ENDDO
ENDIF
- ENDDO
!
- ! If this grid contains the lowest hierarchy of this basin ID, we can also merge without second thought.
+ ! Should we check the opposite point if odd NeighNb?
+ !IF ( modulo(NbNeighb,2) == 0 ) THEN
+ ! ! Do like first one
+ ! ! loc(NbNeighb)
+ ENDIF
!
- IF ( basin_id(sp,sb) < invented_basins ) THEN
- IF ( sp == minhloc(basin_id(sp,sb),1) ) THEN
- sbint = minhloc(basin_id(sp,sb),2)
+ ! If not found, and at least one of the neighbour is ocean -> coastal outflow
+ IF ( found == 0) THEN
+ WRITE (numout,*) "Is there a coastal outflow ? "
+ nbocean = 0
+ DO nb=1,NbNeighb
+ IF ( neighbours_g(sp,nb) .LE. -1 ) nbocean = nbocean +1
+ ENDDO
+ IF ( nbocean .LE. 1 ) THEN
+ outflow_grid(sp,sb) = -2
+ outflow_basin(sp,sb) = undef_int
+ found = 1
+ WRITE (numout,*) "There is more than one neighbour that is coastal outflow"
+ WRITE (numout,*) "Changed in coastal outflow"
+
+ solved(sp,3) = solved(sp,3)+1
ENDIF
ENDIF
- !
- ! Trung: Be careful when changing value of basin_hierarchy here !
- !
- IF ( sbint < undef_int ) THEN
- basin_hierarchy(sp,sb) = basin_hierarchy(sp,sbint)
- CALL routing_updateflow(sp, sb, sp, sbint, nbpt, nwbas, nbxmax*8, outflow_grid, outflow_basin, &
- & inflow_number, inflow_grid, inflow_basin)
- ELSE
- !
- ! Trung: Well, if you come to this step, it should be a special
- ! (or curious, I said) case. For example, when processing
- ! Sulina branch of river Danube with the forcing data
- ! from E2OFD (/bdd/ORCHIDEE_Forcing/BC/OOL/OL2/E2OFD/)
- ! in resolution of 0.25 degree. The grid box [45.125N, 29.375E]
- ! is a ocean point in the middle of 7 land points ("a
- ! hole"). But HydroSHEDS shows that Sulina branch flows
- ! through this grid cell (from [45.125N, 29.125E] to [45.125N, 29.625E]).
- ! Trung: If the "JUMP" does not work, it means a worse case.
- ! Or it means a worse case.
- ! The outlet of this river falls in the ocean point of
- ! forcing data. For example: with 0.25 degree forcing
- ! data of E2OFD, the grid point [ 45.625N ; 13.875E ]
- ! lost the outlet of basin ID 69666.
- ! We need to use this point (which should be the last point
- ! of this river with shortest hierarchy and highest flow accumulation)
- ! to become new outlet.
+
+ ! Equivalent to linkup 5
+ IF ( found == 0 ) THEN
+ WRITE (numout,*) "Eq. Linkup 5"
+ DO sba=1,basin_count(sp)
+ IF ( sba .NE. sb ) THEN
+ IF ( basin_id(sp,sb) == basin_id(sp,sba) ) THEN
+ !
+ ! Only look for suitable hierrachy if we have a solution for the sub-basin : it has an outflow grid
+ ! and basin or flows into the ocean.
+ !
+ IF ( (outflow_grid(sp,sba) > 0 .AND. outflow_basin(sp,sba) < undef_int) .OR. &
+ & outflow_grid(sp,sba) .EQ. -1 .OR. outflow_grid(sp,sba) .EQ. -2 ) THEN
+ !
+ ! calculation of the criterion for considering that the hierarchy are close enough.
+ ! Watch out hierarchy can be zero. 1 is considered the
+ ! minimum here.
+ !
+ ! Trung: If current sub-basin (sp,sb) is the closest
+ ! point to the outlet of river, the value of
+ ! basin_hierarchy(sp,sb) is always lower than
+ ! basin_hierarchy(sp,sba). It means that crit
+ ! is negative and statement IF (crit < 0.1) is
+ ! not useful. => I changed it a little bit.
+ !
+ crit = (basin_hierarchy(sp,sb)-basin_hierarchy(sp,sba))/(basin_hierarchy(sp,sb)+1)
+ IF ( crit < 0.01 .AND. crit >= zero ) THEN
+ ! Trung: Put condition of flow accumulation here to
+ ! avoid error of circulation in routing_reg_fetch
+ IF (basin_fac(sp,sba) .GE. basin_fac(sp,sb)) THEN
+ sbint = sba
+ ENDIF
+ ENDIF
+ ENDIF
+ ENDIF
+ ENDIF
+ ENDDO
!
- ! This is not good solution because there can be few
- ! outlet for one basin ID.
- ! But I need the model works now !!! => so, come back
- ! here later !
+ ! If this grid contains the lowest hierarchy of this basin ID, we can also merge without second thought.
!
- WRITE (numout,*) "Linkup : Made a NEW OUTLET at sp & sb: ",sp,sb
- ! Coastal flow or river flow is both ok here
- outflow_grid(sp,sb) = -2
- basin_hierarchy(sp,sb) = 0.00
- ENDIF
- !
- ! See if we could find the same basin in this grid with a lower hierarchy
- !
- IF ( debug .AND. routing_diagbox_g(sp,diaglalo) ) THEN
- IF ( basin_id(sp,sb) == basin_id(sp,sba) ) THEN
- WRITE(numout,*) "Linkup 5.0 - In routing_reg_linkup working on sp & sb :", sp,sb
- WRITE(numout,*) "Linkup 5.0 - Do we have an internal solution here :", outflow_basin(sp,sb)
- WRITE(numout,*) "Linkup 5.0 - ID = ", basin_id(sp,sb), basin_id(sp,sba)
- WRITE(numout,*) "Linkup 5.0 - H = ", basin_hierarchy(sp,sb), basin_hierarchy(sp,sba)
- WRITE(numout,*) "Linkup 5.0 - Area = ", basin_area(sp,sb), basin_area(sp,sba)
- WRITE(numout,*) "Linkup 5.0 - outflow_grid =", outflow_grid(sp,sb)
- WRITE(numout,*) "Linkup 5.0 - outflow_basin = ", outflow_basin(sp,sb)
+ IF ( basin_id(sp,sb) < invented_basins ) THEN
+ IF ( sp == minhloc(basin_id(sp,sb),1) ) THEN
+ sbint = minhloc(basin_id(sp,sb),2)
ENDIF
- ENDIF
- !
- IF ( outflow_basin(sp,sb) < undef_int ) solved(sp,5) = solved(sp,5)+1
- !
- ENDIF
- !
- ENDDO
- ENDDO
- !
- IF ( debug .AND. testbasinid > 0 ) THEN
- DO sp=1,nbpt
- DO sb=1,basin_count(sp)
- IF ( basin_id(sp,sb) == testbasinid .AND. routing_diagbox_g(sp,diaglalo) ) THEN
- WRITE(numout,*) "Linkup TEST : ", basin_id(sp,sb), "@", sp, sb
- WRITE(numout,*) "Linkup TEST H= ", basin_hierarchy(sp,sb)
- WRITE(numout,*) "Linkup TEST FAC= ", basin_fac(sp,sb)
- WRITE(numout,*) "Linkup TEST outflow_grid ", outflow_grid(sp,sb)
- WRITE(numout,*) "Linkup TEST outflow_basin ", outflow_basin(sp,sb)
- WRITE(numout,*) "Linkup TEST basin_lshead ", basin_lshead(sp,sb)
- WRITE(numout,*) "Linkup TEST basin_flowdir ", basin_flowdir(sp,sb)
- ENDIF
- ENDDO
- ENDDO
- ENDIF
- !
- ! 4.0 We hopefully have solved all the points but here we check it and print sole diagnostics.
- !
- ! Trung: I prefer to write the percent number here in real than integer
- ! because for large domain, 0.01% is really amount.
- !
- WRITE(numout,'(A16,I10,A8,I10,A14)') "Linkup : Out of ", SUM(basin_count), " basins ", SUM(solved), " were resolved."
- DO sb=1,5
- WRITE(numout,'(A9,7X,F10.4,A27,I2.2,A2)') "Linkup : ",REAL(SUM(solved(:,sb)), r_std)/SUM(basin_count)*100.0, &
- & " % of basins solved by case ", sb," ."
- ENDDO
- !
- DO sp=1,nbpt
- !
- IF ( debug ) THEN
- DO sb=1,5
- WRITE(numout,*) sp," Solutions found for case ", sb," : ",&
- & NINT(REAL(solved(sp,sb), r_std)/basin_count(sp)*100.0)," %"
- ENDDO
- IF ( SUM(solved(sp,:)) < basin_count(sp) ) THEN
- WRITE(numout,*) sp,"No solution found for ", basin_count(sp)-SUM(solved(sp,:))," basins."
- ENDIF
- ENDIF
- !
- DO sb=1,basin_count(sp)
- !
- ! OK this is it, we give up :-)
- !
- IF ( outflow_basin(sp,sb) .EQ. undef_int .AND. &
- & outflow_grid(sp,sb) > 0 .AND. basin_flowdir(sp,sb) > 0 ) THEN
- WRITE(numout,*) 'We could not find the basin into which we need to flow, Linkup FAILED'
- WRITE(numout,*) 'Linkup FAILED : Grid point ', sp, ' and basin ', sb
- WRITE(numout,*) 'Linkup FAILED : Coordinates : ', lalo_g(sp,2), lalo_g(sp,1)
- WRITE(numout,*) 'Linkup FAILED : Outflow direction :', basin_flowdir(sp,sb)
- WRITE(numout,*) 'Linkup FAILED : Large scale outflow direction :',basin_lshead(sp,sb)
- WRITE(numout,*) 'Linkup FAILED : Outlfow grid :', outflow_grid(sp,sb)
- WRITE(numout,*) 'Linkup FAILED : Outlfow basin :',outflow_basin(sp,sb)
- WRITE(numout,*) 'Linkup FAILED : basin ID:',basin_id(sp,sb)
- WRITE(numout,*) 'Linkup FAILED : Hierarchy :', basin_hierarchy(sp,sb)
- WRITE(numout,*) 'Linkup FAILED : Flow accumulation :', basin_fac(sp,sb)
- CALL ipslerr_p(3,'routing_reg_linkup', &
- "We could not find the basin into which we need to flow",'Try with debug=.TRUE.','stop routing_reg_linkup')
- ENDIF
- !
- ENDDO
- !
- ENDDO
- !
- ! Check for each outflow basin that it exists
- !
- DO sp=1,nbpt
- DO sb=1,basin_count(sp)
- !
- inp = outflow_grid(sp,sb)
- sbl = outflow_basin(sp,sb)
- IF ( inp .GT. 0 ) THEN
- IF ( basin_count(inp) .LT. sbl ) THEN
- WRITE(numout,*) 'Point :', sp, ' Basin :', sb
- WRITE(numout,*) 'Flows into point :', inp, ' basin :', sbl
- WRITE(numout,*) 'But it flows into now here as basin count = ', basin_count(inp)
+ ENDIF
!
- WRITE(numout,*) 'Input sub-basin :'
- WRITE(numout,*) 'Outflow direction :', basin_flowdir(sp,sb)
- WRITE(numout,*) 'Outlfow grid :', outflow_grid(sp,sb)
- WRITE(numout,*) 'Outlfow basin :',outflow_basin(sp,sb)
- WRITE(numout,*) 'Basin ID:',basin_id(sp,sb)
- WRITE(numout,*) 'Hierarchy :', basin_hierarchy(sp,sb)
- WRITE(numout,*) 'Output grid box:'
- WRITE(numout,*) 'Basin ID:',basin_id(inp,basin_count(inp))
- WRITE(numout,*) 'Hierarchy :', basin_hierarchy(inp,basin_count(inp))
+ ! Trung: Be careful when changing value of basin_hierarchy here !
!
- CALL ipslerr_p(3,'routing_reg_linkup','Problem with outflow','','')
+ IF (( sbint < undef_int ) .AND. (sbint .GT. 0) ) THEN
+ basin_hierarchy(sp,sb) = basin_hierarchy(sp,sbint)
+ CALL routing_updateflow(sp, sb, sp, sbint, nbpt, nwbas, nbxmax*8, outflow_grid, outflow_basin, &
+ & inflow_number, inflow_grid, inflow_basin)
+ WRITE (numout,*) "Lowest basin hierarchy in the grid file"
+
+ ELSE
+ !
+ ! Trung: Well, if you come to this step, it should be a special
+ ! (or curious, I said) case. For example, when processing
+ ! Sulina branch of river Danube with the forcing data
+ ! from E2OFD (/bdd/ORCHIDEE_Forcing/BC/OOL/OL2/E2OFD/)
+ ! in resolution of 0.25 degree. The grid box [45.125N, 29.375E]
+ ! is a ocean point in the middle of 7 land points ("a
+ ! hole"). But HydroSHEDS shows that Sulina branch flows
+ ! through this grid cell (from [45.125N, 29.125E] to [45.125N, 29.625E]).
+ ! Trung: If the "JUMP" does not work, it means a worse case.
+ ! Or it means a worse case.
+ ! The outlet of this river falls in the ocean point of
+ ! forcing data. For example: with 0.25 degree forcing
+ ! data of E2OFD, the grid point [ 45.625N ; 13.875E ]
+ ! lost the outlet of basin ID 69666.
+ ! We need to use this point (which should be the last point
+ ! of this river with shortest hierarchy and highest flow accumulation)
+ ! to become new outlet.
+ !
+ ! This is not good solution because there can be few
+ ! outlet for one basin ID.
+ ! But I need the model works now !!! => so, come back
+ ! here later !
+ !
+ WRITE (numout,*) "Linkup : Made a NEW OUTLET at sp & sb: ",sp,sb
+ ! Coastal flow or river flow is both ok here
+ outflow_grid(sp,sb) = -2
+ basin_hierarchy(sp,sb) = 0.00
+ ENDIF
+ IF ( outflow_basin(sp,sb) < undef_int ) solved(sp,4) = solved(sp,4)+1
ENDIF
ENDIF
ENDDO
@@ -2680,7 +2332,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
!
END SUBROUTINE routing_reg_linkup
!
- !
+ !
!! ================================================================================================================================
!! SUBROUTINE : routing_reg_fetch
!!
@@ -2943,7 +2595,7 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_la
INTEGER(i_std) :: nbtruncate !!
INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: indextrunc !!
!$OMP THREADPRIVATE(indextrunc)
- !
+ !
LOGICAL :: trunc_equa = .TRUE. !! Logical to choose if equality truncation is activated or not (true/false)
!
@@ -3019,7 +2671,7 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_la
tmp_area(multbas) = fetch_basin(ib,ij)
ENDIF
ENDDO
- !
+ !
! Test the new truncation priority
!
IF ( trunc_equa ) THEN
@@ -3604,12 +3256,12 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
SUBROUTINE allocategraph(nbpt, nbasmax)
INTEGER(i_std), INTENT(in) :: nbpt, nbasmax
-
+
INTEGER(i_std) :: ier
nbpt_save = nbpt
nbasmax_save = nbasmax
-
+
ALLOCATE (routing_area_glo(nbpt,nbasmax), stat=ier)
ALLOCATE (route_fetch_glo(nbpt,nbasmax), stat=ier)
ALLOCATE (routing_cg_glo(nbpt,nbasmax,2), stat=ier)
@@ -3622,6 +3274,6 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
ALLOCATE (route_outlet_glo(nbpt,nbasmax,2), stat=ier)
ALLOCATE (route_type_glo(nbpt,nbasmax), stat=ier)
ALLOCATE (origin_nbintobas_glo(nbpt), stat=ier)
-
+
END SUBROUTINE allocategraph
END MODULE routing_reg