Dsum is the current reference configuration until we solve the HTU decomposssition issues.

F90subroutines/routing_reg.f90

@@ -47,7 +47,8 @@ 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 : The average of the dz have to be larger than min(dz) and smaller than orog_max-orog_min
+  ! topo_option = 5 : We only keel the river segment length and elevation change of the grid point
+  !                   with the largest flow accumulation.
   !
   INTEGER(i_std), SAVE :: topo_option = 4 !! Option to calculate topoindex
   !
@@ -58,12 +59,13 @@ MODULE routing_reg
   !
   ! Option to limit the dz when computing topoindex
   !
-  INTEGER(i_std), SAVE :: dzrang_limiter = 1
+  INTEGER(i_std), SAVE :: dzrang_limiter = 0
   ! Tested combinations :
   ! Asum : topo_option = 3, kill_option = 1, dzrang_limiter = 0
   ! Bsum : topo_option = 4, kill_option = 1, dzrang_limiter = 0
   ! 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
   !
 CONTAINS
 !! ================================================================================================================================
@@ -1031,7 +1033,7 @@ CONTAINS
     !
     ! Work only for topo_option number 3 or 4
     !
-    IF ( topo_option .EQ. 3 .OR. topo_option .EQ. 4) THEN
+    IF ( topo_option .GE. 3 ) THEN
        rivpas(:,:) = 0
        ! Loop for all sub-basins and each point belongs to that sub-basin
        DO ij=1,nb_basin
@@ -1667,6 +1669,7 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
     !
     INTEGER(i_std)   :: ij, iz !! Indices (unitless)
     INTEGER(i_std)   :: ip, jp, cnt
+    REAL(r_std)      :: fac
     CHARACTER(LEN=4) :: hierar_method = 'OUTP' !!
     !
     ! To calculate topoindex in new method
@@ -1689,7 +1692,7 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
     !
     IF ( topo_option .LT. 1 .OR. topo_option .GT. 5 ) THEN
        WRITE(numout,*) ' You chose wrong topo_option for calculating topoindex: ', topo_option
-       WRITE(numout,*) ' It should be: 1, 2, 3 or 4'
+       WRITE(numout,*) ' It should be: 1, 2, 3, 4 or 5'
        STOP 'routing_reg_globalize'
     ENDIF
     !
@@ -1859,9 +1862,11 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
           ENDDO
        ENDIF
        !
-       ! Compute mean length of rivers in HTU
+       ! Compute mean length of rivers in HTU in all cases for topo_option <= 4.
+       ! For topo_option == 5 we use the dominant river only.
        !
        cnt = 0
+       IF ( topo_option .LE. 4 ) THEN
           DO iz=1,basin_sz(ij)
              ip = basin_pts(ij,iz,1)
              jp = basin_pts(ij,iz,2)
@@ -1872,6 +1877,25 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
                 basin_rdz(ib,ij) = basin_rdz(ib,ij)+rdz_bx(ip,jp)
              ENDIF
           ENDDO
+       ELSE IF ( topo_option .EQ. 5 ) THEN
+          cnt = 1
+          fac = 0.0
+          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
+             IF ( rlen_bx(ip,jp) .GT. 0 .AND. fac_bx(ip,jp) > fac ) THEN
+                basin_rlen(ib,ij) = rlen_bx(ip,jp)
+                basin_rdz(ib,ij) = rdz_bx(ip,jp)
+             ENDIF
+          ENDDO
+       ELSE
+          WRITE(numout,*) ' We should not be here for topoindex = ', topo_option
+          STOP 'routing_reg_globalize'
+       ENDIF
+       !
+       !
+       !
        IF ( cnt > 0 ) THEN
           basin_rlen(ib,ij) = basin_rlen(ib,ij)/cnt
           IF ( dzrang_limiter == 0 ) THEN
@@ -1887,7 +1911,7 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
        ENDIF
        !
        !
-       IF ( topo_option .EQ. 4 ) THEN
+       IF ( topo_option .EQ. 4 .OR. topo_option .EQ. 5 ) THEN
           basin_topoind(ib,ij) = SQRT(basin_rlen(ib,ij)**3./(basin_rdz(ib,ij)*1.0e6))
        ENDIF
        !