Docker-in-Docker (DinD) capabilities of public runners deactivated. More info

Commit 125ee7f0 authored by POLCHER Jan's avatar POLCHER Jan 🚴🏾
Browse files

Now the sum of the topoindex is done correctly along the graph and then...

Now the sum of the topoindex is done correctly along the graph and then averaged over all branches. This is done for the first HTUs which are selected in routing_reg_findbasin.
parent e85c014a
......@@ -41,7 +41,7 @@ MODULE routing_reg
! values are then averaged.
! topo_option = 4 : As above but the average topoindex is computed for each river.
!
INTEGER(i_std), SAVE :: topo_option = 4 !! Option to calculate topoindex
INTEGER(i_std), SAVE :: topo_option = 3 !! Option to calculate topoindex
!
CONTAINS
!! ================================================================================================================================
......@@ -453,6 +453,10 @@ CONTAINS
INTEGER(i_std) :: ie, je !!
INTEGER(i_std) :: sortind(nbvmax) !!
!
INTEGER(i_std) :: rivpas(ijdimmax,ijdimmax) !! Number of passage through grid when following rivers
INTEGER(i_std) :: rivlen(ijdimmax,ijdimmax) !! number of grids along the river
INTEGER(i_std) :: color(ijdimmax,ijdimmax)
!
INTEGER(i_std) :: itrans !!
INTEGER(i_std) :: trans(nbvmax) !!
!
......@@ -976,10 +980,19 @@ CONTAINS
! 7.0 If we calculate topoindex with option 3
! we need to accumulate the topoind_bx along the flow
!
DO ij=1,nb_basin
DO iz=1,basin_sz(ij)
ip = basin_pts(ij,iz,1)
jp = basin_pts(ij,iz,2)
color(ip,jp) = ij
ENDDO
ENDDO
!
! Work only for topo_option number 3 or 4
!
IF ( topo_option .EQ. 3 .OR. topo_option .EQ. 4 ) THEN
IF ( topo_option .EQ. 3 ) THEN
rivlen(:,:) = 1
rivpas(:,:) = 0
! Loop for all sub-basins and each point belongs to that sub-basin
DO ij=1,nb_basin
DO iz=1,basin_sz(ij)
......@@ -991,7 +1004,8 @@ CONTAINS
p = trip(ip,jp)
cnt = 0
il = ip ; jl = jp
DO WHILE ( p .GT. 0 .AND. p .LT. 9 .AND. cnt .LT. nbi*nbj )
! Only follow the river witin the current HTU
DO WHILE ( ij == color(il,jl) .AND. 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)
......@@ -999,18 +1013,38 @@ CONTAINS
! this point still belongs to current sub-basin. Here I
! assume that the sub-routine routing_reg_divbas worked well.
! So we don't need to check.
topoind(basin_pts(ij,iz,1),basin_pts(ij,iz,2))=topoind(basin_pts(ij,iz,1),basin_pts(ij,iz,2))+topoind(ill,jll)
topoind(ip,jp)=topoind(ip,jp)+topoind(ill,jll)
rivpas(ill,jll) = rivpas(ill,jll)+1
rivlen(ip,jp) = rivlen(ip,jp)+1
!
il = ill ; jl = jll
p = trip(il,jl)
ENDDO
IF ( topo_option .EQ. 4 ) THEN
! Divide by the number of points of the river to get the mean value.
topoind(basin_pts(ij,iz,1),basin_pts(ij,iz,2))=topoind(basin_pts(ij,iz,1),basin_pts(ij,iz,2))/cnt
ENDIF
ENDIF
!
ENDDO
! Second pass through all points of the HTU to keep only the topoindex of full rivers.
! These are those throug which we did not pass when following the rivers.
cnt = 0
DO iz=1,basin_sz(ij)
ip = basin_pts(ij,iz,1)
jp = basin_pts(ij,iz,2)
IF ( rivpas(ip,jp) == 0 ) THEN
! Found a river head
topoind(ip,jp) = topoind(ip,jp) * rivlen(ip,jp)
cnt = cnt + rivlen(ip,jp)
ELSE
topoind(ip,jp) = 0.0
ENDIF
ENDDO
! Finish the weight computation of the rivers followed
IF ( cnt > 0 ) THEN
DO iz=1,basin_sz(ij)
ip = basin_pts(ij,iz,1)
jp = basin_pts(ij,iz,2)
topoind(ip,jp) = topoind(ip,jp)/cnt
ENDDO
ENDIF
ENDDO
ENDIF
!
......@@ -1604,9 +1638,9 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
!_ ================================================================================================================================
!
!
IF ( topo_option .LT. 1 .OR. topo_option .GT. 4 ) THEN
IF ( topo_option .LT. 1 .OR. topo_option .GT. 3 ) THEN
WRITE(numout,*) ' You chose wrong topo_option for calculating topoindex: ', topo_option
WRITE(numout,*) ' It should be: 1, 2, 4 or 4'
WRITE(numout,*) ' It should be: 1, 2, or 3'
STOP 'routing_reg_globalize'
ENDIF
!
......@@ -1764,33 +1798,11 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
ENDDO
basin_topoind(ib,ij) = basin_topoind(ib,ij)/REAL(basin_sz(ij),r_std)
ENDIF
! Minus outlet hierarchy
IF ( topo_option .EQ. 3 .OR. topo_option .EQ. 4 ) THEN
! All topoindex with HTU have already been weighted and those not at the root of the river set to zero.
IF ( topo_option .EQ. 3 ) THEN
DO iz=1,basin_sz(ij)
area_frac = area_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))/REAL(basin_area(ib,ij))
basin_topoind(ib,ij) = basin_topoind(ib,ij) + topoind_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))*area_frac
basin_topoind(ib,ij) = basin_topoind(ib,ij) + topoind_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
ENDDO
!
! Trung: In routing_reg_globalize:
! Trung: Checking negative value of basin_topoind
!
IF ( debug ) THEN
IF ( basin_topoind(ib,ij) .LT. 0. ) THEN
WRITE (numout,*) 'In routing_reg_globalize: We got negative value in ',ib,ij
WRITE (numout,*) 'Start checking ... '
WRITE (numout,*) 'basin_area ', REAL(basin_area(ib,ij))
WRITE (numout,*) ' ==================> '
DO iz=1,basin_sz(ij)
area_frac = area_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))/REAL(basin_area(ib,ij))
WRITE (numout,*) 'iz = ',iz
WRITE (numout,*) 'topoind_bx :', topoind_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
WRITE (numout,*) 'area_bx ', area_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
WRITE (numout,*) 'area_frac ', area_frac
WRITE (numout,*) 'topoind_bx * area_frac ', topoind_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))*area_frac
ENDDO
ENDIF
ENDIF
!
ENDIF
!
! To make sure that it has the lowest number if this is an outflow point we reset basin_hierarchy
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment