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

Add an option to compute the mean topoindex with all points along the rivers of the HTU.

parent 694db137
......@@ -47,8 +47,10 @@ MODULE routing_reg
! values are then averaged over the HTU.
! topo_option = 4 : The length and elevation change of all rivers within the HTU are computed.
! These values are then used to compute the topoindex of the HTU.
! topo_option = 5 : We only keel the river segment length and elevation change of the grid point
! topo_option = 5 : We only keep the river segment length and elevation change of the grid point
! with the largest flow accumulation.
! topo_option = 6 : Keep and average of all sums of topoindex which are build for all pixels
! within the HTU
!
INTEGER(i_std), SAVE :: topo_option = 4 !! Option to calculate topoindex
!
......@@ -66,6 +68,7 @@ MODULE routing_reg
! Csum : topo_option = 4, kill_option = 2, dzrang_limiter = 0
! Dsum : topo_option = 4, kill_option = 2, dzrang_limiter = 1
! Esum : topo_option = 5, kill_option = 2, dzrang_limiter = 1
! Xsum : topo_option = 6, kill_option = 2, dzrang_limiter = 0
!
CONTAINS
!! ================================================================================================================================
......@@ -1046,7 +1049,7 @@ CONTAINS
p = trip(ip,jp)
cnt = 0
il = ip ; jl = jp
! Only follow the river witin the current HTU
! Only follow the river within 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)
......@@ -1067,22 +1070,26 @@ CONTAINS
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.
! Second pass through all points of the HTU to keep only the topoindex of full river segments.
! These are those through 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
cnt = cnt + rweight(ip,jp)
IF ( topo_option .EQ. 6 ) THEN
cnt = cnt + rweight(ip,jp)
ELSE
topoind(ip,jp) = zero
rlen(ip,jp) = zero
rdz(ip,jp) = zero
rweight(ip,jp) = zero
IF ( rivpas(ip,jp) == 0 ) THEN
cnt = cnt + rweight(ip,jp)
ELSE
topoind(ip,jp) = zero
rlen(ip,jp) = zero
rdz(ip,jp) = zero
rweight(ip,jp) = zero
ENDIF
ENDIF
ENDDO
! Finish the weight computation of the rivers followed
! Finish the weight computation of the rivers followed with a normalisation
IF ( cnt > 0 ) THEN
DO iz=1,basin_sz(ij)
ip = basin_pts(ij,iz,1)
......@@ -1690,9 +1697,9 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
!_ ================================================================================================================================
!
!
IF ( topo_option .LT. 1 .OR. topo_option .GT. 5 ) THEN
IF ( topo_option .LT. 1 .OR. topo_option .GT. 6 ) THEN
WRITE(numout,*) ' You chose wrong topo_option for calculating topoindex: ', topo_option
WRITE(numout,*) ' It should be: 1, 2, 3, 4 or 5'
WRITE(numout,*) ' It should be: 1, 2, 3, 4, 5 or 6'
STOP 'routing_reg_globalize'
ENDIF
!
......@@ -1853,8 +1860,9 @@ 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
! Weight averaged topoindex of rivers within HTU. Those not at the root of the river set to zero.
IF ( topo_option .EQ. 3 ) THEN
! Weight averaged topoindex of rivers within HTU. Those not at the root of the river set to zero in topo_options=3.
! For topo_options=6 all points within the HTU are counted.
IF ( topo_option .EQ. 3 .OR. topo_option .EQ. 6 ) THEN
DO iz=1,basin_sz(ij)
ip = basin_pts(ij,iz,1)
jp = basin_pts(ij,iz,2)
......@@ -1866,11 +1874,11 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
! For topo_option == 5 we use the dominant river only.
!
cnt = 0
IF ( topo_option .LE. 4 ) THEN
IF ( topo_option .LE. 4 .OR. topo_option .EQ. 6 ) THEN
DO iz=1,basin_sz(ij)
ip = basin_pts(ij,iz,1)
jp = basin_pts(ij,iz,2)
! Both rlen and rdz are zero outside of the first point in a river of the HTU
! Both rlen and rdz are zero outside of the first point in a river of the HTU (all cases except 6)
IF ( rlen_bx(ip,jp) .GT. 0 ) THEN
cnt = cnt + 1
basin_rlen(ib,ij) = basin_rlen(ib,ij)+rlen_bx(ip,jp)
......@@ -1905,7 +1913,7 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
basin_rdz(ib,ij) = MAX(MIN(basin_rdz(ib,ij)/cnt, basin_orog_max(ib,ij)-basin_orog_min(ib,ij)), hydro_mindz)
ENDIF
ELSE
! Mean values of oringal hydro date are used to fill-in missing points.
! Mean values of original hydro data are used to fill-in missing points.
basin_rlen(ib,ij) = hydro_meanlen
basin_rdz(ib,ij) = hydro_meandz
ENDIF
......
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