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

The gathering of all outflow points from the atmopsheric grid and their...

The gathering of all outflow points from the atmopsheric grid and their catchments is now in a separate subrouting.
parent 098102aa
......@@ -41,6 +41,10 @@ MODULE routing_reg
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: route_fetch_glo !! Upstream area at each HTU
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: routing_floodcri_glo !! Upstream area at each HTU
!
! Number of neighbours on the HydroShed grid (regular Lat/Lon)
!
INTEGER(i_std), SAVE, DIMENSION(nbne,2) :: inc !! Index increments to follow rivers with trip.
!
! Options to compute the topographic index based on the information available in the Hydrological files.
! topo_option = 1 : Constant topographic index
! topo_option = 2 : Simple average without any weighting
......@@ -169,15 +173,20 @@ CONTAINS
!! LOCAL VARIABLES
INTEGER(i_std) :: ip, jp, ll(1), iloc, jloc !! Indices (unitless)
INTEGER(i_std) :: ipp, jpp
INTEGER(i_std), DIMENSION(8,2) :: inc
REAL(r_std) :: cenlon, cenlat, dlon, dlat, deslon, deslat, facti, factj
REAL(r_std) :: lonstr(ijdimmax*ijdimmax) !!
REAL(r_std) :: latstr(ijdimmax*ijdimmax) !!
!
!_ ================================================================================================================================
!
topoindexmin = min_topoind
hydro_meanlen = meanrlen
hydro_meandz = meanrdz
hydro_mindz = mindz
!
! The routing code (i=1, j=2)
! This codng could change with the hydrological file we read.
!
inc(1,1) = 0
inc(1,2) = -1
......@@ -196,11 +205,6 @@ CONTAINS
inc(8,1) = -1
inc(8,2) = -1
!
topoindexmin = min_topoind
hydro_meanlen = meanrlen
hydro_meandz = meanrdz
hydro_mindz = mindz
!
! Set everything to undef to locate easily empty points
!
trip_bx(:,:) = undef_int
......@@ -482,7 +486,6 @@ CONTAINS
CHARACTER(LEN=9) :: fmtr !!
!
LOGICAL :: newpoint !! (true/false)
INTEGER(i_std), DIMENSION(8,2) :: inc !!
INTEGER(i_std), DIMENSION(nbi,nbj) :: outtmp !!
!
INTEGER(i_std) :: nbb, p, cnt, ibas !! To count number of subbasin
......@@ -492,6 +495,8 @@ CONTAINS
INTEGER(i_std) :: tbname(nb_htu) !!
INTEGER(i_std) :: tsz(nb_htu) !!
INTEGER(i_std) :: tpts(nb_htu,nbv,2) !!
INTEGER(i_std) :: slsz
INTEGER(i_std) :: slbas(nb_htu)
INTEGER(i_std) :: toutdir(nb_htu) !!
INTEGER(i_std) :: toutbas(nb_htu) !!
INTEGER(i_std) :: toutloc(nb_htu,2) !!
......@@ -507,6 +512,7 @@ CONTAINS
!
INTEGER(i_std) :: rivpas(ijdimmax,ijdimmax) !! Number of passage through grid when following rivers
INTEGER(i_std) :: color(ijdimmax,ijdimmax)
REAL(r_std) :: fac_loc(ijdimmax,ijdimmax)
!
INTEGER(i_std) :: itrans !!
INTEGER(i_std) :: trans(nbvmax) !!
......@@ -524,7 +530,6 @@ CONTAINS
INTEGER(i_std) :: new_pts(nb_htu, nbv, 2) !!
INTEGER(i_std) :: oldorder, neworder !!
!
INTEGER(i_std), PARAMETER :: nbne=8
REAL(r_std) :: fac_out(nb_htu) !!
INTEGER(i_std) :: ff(1) !!
INTEGER(i_std) :: sortedtrifac(nbne) !
......@@ -549,25 +554,6 @@ CONTAINS
ENDDO
ENDIF
!
! The routing code (i=1, j=2)
!
inc(1,1) = 0
inc(1,2) = -1
inc(2,1) = 1
inc(2,2) = -1
inc(3,1) = 1
inc(3,2) = 0
inc(4,1) = 1
inc(4,2) = 1
inc(5,1) = 0
inc(5,2) = 1
inc(6,1) = -1
inc(6,2) = 1
inc(7,1) = -1
inc(7,2) = 0
inc(8,1) = -1
inc(8,2) = -1
!
! 1.0 Find number of outflow point (>= 97)
! So basically, we get nb_basin, basin_bxout, basin_inbxid
! (with nbb, toutdir, tbname)
......@@ -609,150 +595,59 @@ CONTAINS
! 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)
!
totsz = 0
tsz(:) = 0
tpts(:,:,:) = 0
CALL routing_reg_complocalfac(ib, nbi, nbj, nb_htu, nbv, nbb, basin, trip, outtmp, totsz, tsz, tpts, fac_loc)
!
DO ip=1,nbi
DO jp=1,nbj
IF ( basin(ip,jp) .LT. undef_int) THEN
! Think about basin .LT. undef_int here
!IF ( trip(ip,jp) .GT. 0 .AND. trip(ip,jp) .LT. undef_int ) THEN
IF ( trip(ip,jp) .GT. 0 .AND. trip(ip,jp) .LT. 109 ) THEN
p = trip(ip,jp)
cnt = 0
il = ip ; jl = jp
newpoint = .TRUE.
DO WHILE ( p .GT. 0 .AND. p .LT. 9 .AND. cnt .LT. nbi*nbj .AND. newpoint )
cnt = cnt + 1
ill = il + inc(p,1)
jll = jl + inc(p,2)
il = ill ; jl = jll
p = trip(il,jl)
IF ( outtmp(il,jl) .NE. -1 ) newpoint = .FALSE.
ENDDO
!
!
IF ( cnt .EQ. nbi*nbj) THEN
IF ( debug ) THEN
WRITE(numout,*) "Error: cnt .EQ. nbi*nbj, ib= ",ib
WRITE(numout,*) "Point: ", ip, jp
WRITE(numout,*) "Per: ", nbi, nbj
WRITE(numout,*) "Number of outlet: ", nbb, cnt
WRITE(numout,*) '+++++++++++++++++++ OUTLET IN FINDBASINS ++++++++++++++++++++'
WRITE(fmt,"('(',I3,'I6)')") nbi
DO je=1,nbj
WRITE(numout,fmt) outtmp(1:nbi,je)
ENDDO
!
WRITE(numout,*) '+++++++++++++++++++ TRIP IN FINDBASINS ++++++++++++++++++++'
DO je=1,nbj
WRITE(numout,fmt) trip(1:nbi,je)
ENDDO
!
WRITE(numout,*) '+++++++++++++++++++ BASIN IDs IN FINDBASINS ++++++++++++++++++++'
DO je=1,nbj
WRITE(numout,fmt) basin(1:nbi,je)
ENDDO
!
WRITE(numout,*) '+++++++++++++++++++ LONGITUDE IN FINDBASINS ++++++++++++++++++++'
WRITE(fmtr,"('(',I3,'F8.1)')") nbi
DO je=1,nbj
WRITE(numout,fmtr) lontmp(1:nbi,je)
ENDDO
!
WRITE(numout,*) '+++++++++++++++++++ LATITUDE IN FINDBASINS ++++++++++++++++++++'
WRITE(fmtr,"('(',I3,'F8.1)')") nbi
DO je=1,nbj
WRITE(numout,fmtr) lattmp(1:nbi,je)
ENDDO
!
WRITE(numout,*) '+++++++++++++++++++ HIERARCHY IN FINDBASINS ++++++++++++++++++++'
WRITE(fmtr,"('(',I3,'F8.1)')") nbi
DO je=1,nbj
WRITE(numout,fmtr) hierarchy(1:nbi,je)
ENDDO
CALL FLUSH(numout)
ENDIF
CALL ipslerr_p(3,'routing_reg_findbasins','We could not route point','','')
ENDIF
!
IF ( outtmp(il,jl) .NE. -1 ) THEN
ibas = outtmp(il,jl)
tsz(ibas) = tsz(ibas) + 1
IF ( tsz(ibas) .GT. nbvmax ) CALL ipslerr_p(3,'routing_reg_findbasins','nbvmax too small','second section','')
tpts(ibas, tsz(ibas), 1) = ip
tpts(ibas, tsz(ibas), 2) = jp
outtmp(ip,jp) = outtmp(il,jl)
ELSE
IF ( debug ) THEN
WRITE(numout,*) "Error: outtmp(il,jl) .EQ. -1 ", ib
WRITE(numout,*) "Point: ", ip, jp
WRITE(numout,*) "Per: ", nbi, nbj
WRITE(numout,*) "Number of outlet: ", nbb
WRITE(numout,*) '+++++++++++++++++++ OUTLET IN FINDBASINS ++++++++++++++++++++'
WRITE(fmt,"('(',I3,'I6)')") nbi
DO je=1,nbj
WRITE(numout,fmt) outtmp(1:nbi,je)
ENDDO
!
WRITE(numout,*) '+++++++++++++++++++ TRIP IN FINDBASINS ++++++++++++++++++++'
DO je=1,nbj
WRITE(numout,fmt) trip(1:nbi,je)
ENDDO
!
WRITE(numout,*) '+++++++++++++++++++ BASIN IDs IN FINDBASINS ++++++++++++++++++++'
DO je=1,nbj
WRITE(numout,fmt) basin(1:nbi,je)
ENDDO
!
WRITE(numout,*) '+++++++++++++++++++ LONGITUDE IN FINDBASINS ++++++++++++++++++++'
WRITE(fmtr,"('(',I3,'F8.1)')") nbi
DO je=1,nbj
WRITE(numout,fmtr) lontmp(1:nbi,je)
ENDDO
!
WRITE(numout,*) '+++++++++++++++++++ LATITUDE IN FINDBASINS ++++++++++++++++++++'
WRITE(fmtr,"('(',I3,'F8.1)')") nbi
DO je=1,nbj
WRITE(numout,fmtr) lattmp(1:nbi,je)
ENDDO
!
WRITE(numout,*) '+++++++++++++++++++ HIERARCHY IN FINDBASINS ++++++++++++++++++++'
WRITE(fmtr,"('(',I3,'F8.1)')") nbi
DO je=1,nbj
WRITE(numout,fmtr) hierarchy(1:nbi,je)
ENDDO
ENDIF
CALL FLUSH(numout)
WRITE (numout,*) ' outtmp(il,jl):', outtmp(il,jl)
WRITE (numout,*) ' trip(il,jl):', trip(il,jl)
CALL ipslerr_p(3,'routing_reg_findbasins','We could not find right outlet','','')
ENDIF
!
totsz = totsz + 1
ENDIF
! Testing this IF
ENDIF
ENDDO
slsz=0
DO ibas = 1, nbb
! 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. &
& tsz(ibas)/ REAL(totsz) < htufrac_min ) THEN
! This HTU could be deleted !
slsz = slsz+1
slbas(slsz) = ibas
ENDIF
ENDDO
!
IF ( debug .AND. routing_diagbox(ib, diaglalo) ) THEN
IF ( slsz > 0 ) THEN
color(:,:) = 0
DO jp = 1, slsz
ibas = slbas(jp)
! 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
DO ip=1,tsz(ibas)
color(tpts(ibas,ip,1),tpts(ibas,ip,2)) = outtmp(tpts(ibas,1,1),tpts(ibas,1,2))
ENDDO
ENDIF
ENDDO
ENDIF
!
!
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 outlet: ", nbb
WRITE(numout,*) '+++++++++++++++++++ OUTLET IN FINDBASINS ++++++++++++++++++++'
WRITE(fmt,"('(',I3,'I6)')") nbi
WRITE(numout,*) "Number of HTU smaller than ",htufrac_min, " = ", slsz
WRITE(numout,*) '+++++++++++++++++++ TRIP IN FINDBASINS ++++++++++++++++++++'
WRITE(fmt,"('(',I3,'I8)')") nbi
DO je=1,nbj
WRITE(numout,fmt) outtmp(1:nbi,je)
WRITE(numout,fmt) MIN(trip(1:nbi,je), 999)
ENDDO
!
WRITE(numout,*) '+++++++++++++++++++ TRIP IN FINDBASINS ++++++++++++++++++++'
WRITE(numout,*) '+++++++++++++++++++ LOCAL BASINS IN FINDBASINS ++++++++++++++++++++'
DO je=1,nbj
WRITE(numout,fmt) trip(1:nbi,je)
WRITE(numout,fmt) color(1:nbi,je)
ENDDO
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 global IN FINDBASINS ++++++++++++++++++++'
WRITE(fmtr,"('(',I3,'F8.1)')") nbi
DO je=1,nbj
WRITE(numout,fmtr) fac(1:nbi,je)
ENDDO
ENDIF
!
!!!!!!!!!!!!!!!!!!!!
! fac_lim : flow accumulation of the HTUs with the fourth largest fac
fac_out(:) = flag
......@@ -1120,8 +1015,160 @@ CONTAINS
!
END SUBROUTINE routing_reg_findbasins
!
!! ================================================================================================================================
!! SUBROUTINE : routing_reg_complocalfac
!!
!>\BRIEF This subroutine gathers all the points with the same outflow and computes the flow accumulation
!! over the entire catchment.
!!
!! DESCRIPTION (definitions, functional, design, flags) :
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S):
!!
!! REFERENCES : None
!!
!! FLOWCHART : None
!! \n
!_ ================================================================================================================================
SUBROUTINE routing_reg_complocalfac(ib, nbi, nbj, nb_htu, nbv, nbb, basin, trip, outtmp, totsz, tsz, tpts, fac_loc)
!
!! INPUT Variables
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) :: basin(:,:)
INTEGER(i_std), INTENT(in) :: trip(:,:)
!!
INTEGER(i_std), INTENT(inout) :: outtmp(:,:)
!! OUTPUT Variables
INTEGER(i_std), INTENT(out) :: totsz
INTEGER(i_std), INTENT(out) :: tsz(nb_htu)
INTEGER(i_std), INTENT(out) :: tpts(nb_htu,nbv,2)
REAL(r_std), INTENT(out) :: fac_loc(nbi,nbj)
!! LOCAL
INTEGER(i_std) :: ip, jp, il, jl, ill, jll
INTEGER(i_std) :: p, cnt, ibas, je
CHARACTER(LEN=7) :: fmt
CHARACTER(LEN=9) :: fmtr
!
! Initialisation
!
tsz(:) = 0
tpts(:,:,:) = 0
fac_loc(:,:) = zero
!
! Add first all outflow points
!
DO ip=1,nbi
DO jp=1,nbj
IF ( outtmp(ip,jp) .GE. 0 ) THEN
ibas = outtmp(ip,jp)
tsz(ibas) = 1
tpts(ibas, 1, 1) = ip
tpts(ibas, 1, 2) = jp
ENDIF
ENDDO
ENDDO
!
! Loops over all points of the hydro grid
!
DO ip=1,nbi
DO jp=1,nbj
IF ( outtmp(ip,jp) .LT. 0 ) THEN
IF ( trip(ip,jp) .GT. 0 .AND. trip(ip,jp) .LT. 9 ) THEN
p = trip(ip,jp)
cnt = 0
il = ip ; jl = jp
DO WHILE ( p .GT. 0 .AND. p .LT. 9 .AND. cnt .LT. nbi*nbj)
cnt = cnt + 1
ill = il + inc(p,1)
jll = jl + inc(p,2)
il = ill ; jl = jll
p = trip(il,jl)
fac_loc(il, jl) = fac_loc(il, jl) + 1
ENDDO
!
!
IF ( cnt .EQ. nbi*nbj) THEN
WRITE(numout,*) "Error: cnt .EQ. nbi*nbj on ib=", ib
WRITE(numout,*) "Point: ", ip, jp
WRITE(numout,*) "Per: ", nbi, nbj
WRITE(numout,*) "Number of outlet: ", nbb, cnt
WRITE(numout,*) '+++++++++++++++++++ OUTLET IN FINDBASINS ++++++++++++++++++++'
WRITE(fmt,"('(',I3,'I6)')") nbi
DO je=1,nbj
WRITE(numout,fmt) outtmp(1:nbi,je)
ENDDO
!
WRITE(numout,*) '+++++++++++++++++++ TRIP IN FINDBASINS ++++++++++++++++++++'
DO je=1,nbj
WRITE(numout,fmt) trip(1:nbi,je)
ENDDO
!
WRITE(numout,*) '+++++++++++++++++++ BASIN IDs IN FINDBASINS ++++++++++++++++++++'
DO je=1,nbj
WRITE(numout,fmt) basin(1:nbi,je)
ENDDO
!
CALL FLUSH(numout)
CALL ipslerr_p(3,'routing_reg_findbasins','We could not route point','','')
ENDIF
!
IF ( outtmp(il,jl) .NE. -1 ) THEN
ibas = outtmp(il,jl)
tsz(ibas) = tsz(ibas) + 1
IF ( tsz(ibas) .GT. nbvmax ) CALL ipslerr_p(3,'routing_reg_findbasins','nbvmax too small','second section','')
tpts(ibas, tsz(ibas), 1) = ip
tpts(ibas, tsz(ibas), 2) = jp
outtmp(ip,jp) = outtmp(il,jl)
ELSE
WRITE(numout,*) "Error: outtmp(il,jl) .EQ. -1 on ib=", ib
WRITE(numout,*) "Point: ", ip, jp
WRITE(numout,*) "End Point :" , il, jl
WRITE(numout,*) "Per: ", nbi, nbj
WRITE(numout,*) "Number of outlet: ", nbb
WRITE(numout,*) '+++++++++++++++++++ OUTLET IN FINDBASINS ++++++++++++++++++++'
WRITE(fmt,"('(',I3,'I6)')") nbi
DO je=1,nbj
WRITE(numout,fmt) outtmp(1:nbi,je)
ENDDO
!
WRITE(numout,*) '+++++++++++++++++++ TRIP IN FINDBASINS ++++++++++++++++++++'
DO je=1,nbj
WRITE(numout,fmt) MIN(trip(1:nbi,je), 999)
ENDDO
!
WRITE(numout,*) '+++++++++++++++++++ BASIN IDs IN FINDBASINS ++++++++++++++++++++'
DO je=1,nbj
WRITE(numout,fmt) basin(1:nbi,je)
ENDDO
!
WRITE (numout,*) ' outtmp(il,jl):', outtmp(il,jl)
WRITE (numout,*) ' trip(il,jl):', trip(il,jl)
CALL ipslerr_p(3,'routing_reg_findbasins','We could not find right outlet','','')
ENDIF
!
ENDIF
! Testing this IF
ENDIF
ENDDO
ENDDO
!
IF (SUM(tsz) == COUNT(trip .LT. undef_int)) THEN
totsz = SUM(tsz)
ELSE
CALL ipslerr_p(3,'routing_reg_complocalfac','We could not find an outlet for each pixel.','','')
ENDIF
!
END SUBROUTINE routing_reg_complocalfac
!
!
!! ================================================================================================================================
!! SUBROUTINE : routing_reg_divbas_search
!!
......@@ -1162,7 +1209,7 @@ CONTAINS
REAL(r_std), INTENT(in) :: headd
INTEGER(i_std), INTENT(in) :: tpts(:,:,:) !!
INTEGER(i_std), INTENT(in) :: basin(:,:) !!
REAL(r_std), INTENT(in) :: fac(:,:) !!
REAL(r_std), INTENT(in) :: fac(:,:) !! Global Flow Accumulation
REAL(r_std), INTENT(in) :: lon(:,:), lat(:,:) !!
INTEGER(i_std), INTENT(in) :: option !!
INTEGER(i_std), INTENT(in) :: trip(:,:) !!
......@@ -1183,8 +1230,8 @@ CONTAINS
CHARACTER(LEN=11) :: afmt !!
CHARACTER(LEN=13) :: afmtr !!
!
REAL(r_std), DIMENSION(nbi,nbj) :: factmp_loc !!
REAL(r_std), DIMENSION(nbi,nbj) :: factmp_glo !!
REAL(r_std), DIMENSION(nbi,nbj) :: factmp_loc !! Local Flow Accumulation
INTEGER(i_std), DIMENSION(nbi,nbj) :: triptmp !!
INTEGER(i_std), DIMENSION(nbi,nbj) :: triptemp !!
INTEGER(i_std), DIMENSION(nbi,nbj) :: basintmp !!
......@@ -1198,10 +1245,6 @@ CONTAINS
!
LOGICAL :: okpoint=.TRUE.
!
! Number of neighbours on the HydroShed grid (regular Lat/Lon)
!
INTEGER(i_std), PARAMETER :: nbne=8
INTEGER(i_std), DIMENSION(nbne,2) :: inc !!
!
! Explanation for dimension of below arrays:
!
......@@ -1251,25 +1294,6 @@ CONTAINS
!_ ================================================================================================================================
!
! The routing code (i=1, j=2)
!
inc(1,1) = 0
inc(1,2) = -1
inc(2,1) = 1
inc(2,2) = -1
inc(3,1) = 1
inc(3,2) = 0
inc(4,1) = 1
inc(4,2) = 1
inc(5,1) = 0
inc(5,2) = 1
inc(6,1) = -1
inc(6,2) = 1
inc(7,1) = -1
inc(7,2) = 0
inc(8,1) = -1
inc(8,2) = -1
!
! 0.0 Get the information of all the points belongs to the subbasin
! which need to divide. Flag for flow accumulation is -9999 (the value of
! flow accumulation should be always non negative).
......@@ -1547,11 +1571,6 @@ SUBROUTINE routing_reg_divbas_cut(nb_htu, nbv, nbi, nbj, ibas, iloc, jloc, tsz,
!
LOGICAL :: okpoint=.TRUE.
!
! Number of neighbours on the HydroShed grid (regular Lat/Lon)
!
INTEGER(i_std), PARAMETER :: nbne=8
INTEGER(i_std), DIMENSION(nbne,2) :: inc !!
!
! Explanation for dimension of below arrays:
!
! + Firstly, I suggest that you should read about the Pfafstetter Coding
......@@ -1600,24 +1619,7 @@ SUBROUTINE routing_reg_divbas_cut(nb_htu, nbv, nbi, nbj, ibas, iloc, jloc, tsz,
!_ ================================================================================================================================
!
! The routing code (i=1, j=2)
!
inc(1,1) = 0
inc(1,2) = -1
inc(2,1) = 1
inc(2,2) = -1
inc(3,1) = 1
inc(3,2) = 0
inc(4,1) = 1
inc(4,2) = 1
inc(5,1) = 0
inc(5,2) = 1
inc(6,1) = -1
inc(6,2) = 1
inc(7,1) = -1
inc(7,2) = 0
inc(8,1) = -1
inc(8,2) = -1
!
!
triptmp(:,:) = -1
triptemp(:,:) = -1
......@@ -1727,7 +1729,7 @@ SUBROUTINE routing_reg_divbas_cut(nb_htu, nbv, nbi, nbj, ibas, iloc, jloc, tsz,
ENDDO
!
IF ( checksz .NE. tsz) THEN
WRITE(numout,*) ' ater got lost while I tried to divide basins'
WRITE(numout,*) ' Water got lost while I tried to divide basins'
WRITE(numout,*) ' Number of new sub-basin : ', new_nb
ip = 2 !COUNT(divloc(:,1) .NE. 0)
WRITE(numout,*) ' Number of mainstrem sub-basin : ', oic
......@@ -1961,10 +1963,6 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz,
LOGICAL :: okpoint=.TRUE.
LOGICAL :: debug =.FALSE.
!
! Number of neighbours on the HydroShed grid (regular Lat/Lon)
!
INTEGER(i_std), PARAMETER :: nbne=8
INTEGER(i_std), DIMENSION(nbne,2) :: inc !!
!
INTEGER(i_std) :: main_loc(nbv,2), tri_loc(nbne,nbv,2)
REAL(r_std) :: main_len(nbv),main_orog(nbv),main_dz(nbv)
......@@ -1979,25 +1977,6 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz,
!_ ================================================================================================================================
!
! The routing code (i=1, j=2)
!
inc(1,1) = 0
inc(1,2) = -1
inc(2,1) = 1