Some clean-up in the options to compute topoindex per HTU. A 4th option was...

Some clean-up in the options to compute topoindex per HTU. A 4th option was introduced to take the mean of topoindex along the river within the HTU.

F90subroutines/routing_reg.f90

@@ -33,7 +33,16 @@ MODULE routing_reg
   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
   REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)    :: routing_floodcri_glo !! Upstream area at each HTU
-
+  !
+  ! 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
+  ! topo_option = 3 : The sum over all rivers of the HTU are computed in _findbasin. These
+  !                   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
+  !
 CONTAINS
 !! ================================================================================================================================
 !! SUBROUTINE : routing_reg_getgrid
@@ -460,8 +469,6 @@ CONTAINS
     INTEGER(i_std) :: new_pts(nb_htu, nbv, 2) !!
     INTEGER(i_std) :: oldorder, neworder !!
     !
-    INTEGER(i_std) :: option !! Option to calculate topoindex
-
 !_ ================================================================================================================================
     !
     IF ( debug .AND. routing_diagbox(ib, diaglalo) ) THEN
@@ -969,31 +976,10 @@ CONTAINS
     ! 7.0 If we calculate topoindex with option 3
     !     we need to accumulate the topoind_bx along the flow
     !
-    ! Get option for calculating topoindex
-    !
-    !Config Key = TOPOINDEX
-    !Config Desc = Options to calculate topoindex
-    !Config Def = 3
-    !Config Help = This flag allows to calculate topoindex with different options
-    !Config Units = [-]
-    !
-    option = 3
-!!$    CALL getin('TOPOINDEX', option)
-!!$    !
-!!$    ! Checking option: until now (May 2016) there is only 3 acceptable options
-!!$    ! = 1 : topoindex .EQ. 1000. everywhere
-!!$    ! = 2 : simple average as the method in the older routing scheme
-!!$    ! = 3 : accumulate the "topoind" from routing.nc along the flow (as Agnes DUCHARNE's suggestion)
-!!$    !
-!!$    IF ( option .LT. 1 .OR. option .GT. 3 ) THEN
-!!$       WRITE(numout,*) ' You chose wrong option for calculating topoindex: ', option
-!!$       WRITE(numout,*) ' It should be: 1, 2 or 3'
-!!$       STOP 'routing_reg_globalize'
-!!$    ENDIF
-    !
-    ! Work only using option number 3
-    !
-    IF ( option .EQ. 3 ) THEN
+    !
+    ! Work only for topo_option number 3 or 4
+    !
+    IF ( topo_option .EQ. 3 .OR. topo_option .EQ. 4 ) THEN
        ! Loop for all sub-basins and each point belongs to that sub-basin
        DO ij=1,nb_basin
           DO iz=1,basin_sz(ij)
@@ -1018,6 +1004,10 @@ CONTAINS
                    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
@@ -1598,7 +1588,6 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
     !
     ! To calculate topoindex in new method
     !
-    INTEGER(i_std) :: option !! Option to calculate topoindex
     REAL(r_std) :: area_frac !! Area fraction for calculate topoindex
     REAL(r_std) :: temp !! local variable
     !
@@ -1614,21 +1603,15 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
     !
 !_ ================================================================================================================================
     !
-    ! Get option for calculating topoindex
-    !
-    option = 3
-    !
-    ! Checking option: until now (May 2016) there is only 3 acceptable options
-    ! = 1 : topoindex .EQ. 1000. everywhere
-    ! = 2 : simple average as the method in the older routing scheme
-    ! = 3 : accumulate the "topoind" from routing.nc along the flow (as Agnes DUCHARNE's suggestion)
     !
-    IF ( option .LT. 1 .OR. option .GT. 3 ) THEN
-       WRITE(numout,*) ' You chose wrong option for calculating topoindex: ', option
-       WRITE(numout,*) ' It should be: 1, 2 or 3'
+    IF ( topo_option .LT. 1 .OR. topo_option .GT. 4 ) THEN
+       WRITE(numout,*) ' You chose wrong topo_option for calculating topoindex: ', topo_option
+       WRITE(numout,*) ' It should be: 1, 2, 4 or 4'
        STOP 'routing_reg_globalize'
     ENDIF
     !
+    basin_topoind(ib,:) = undef_sechiba
+    !
     DO ij=1, nb_basin
        !
        ! Count the basins and keep their ID
@@ -1771,18 +1754,18 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
        !
        basin_topoind(ib,ij) = zero
        ! Constant value (1000.)
-       IF ( option .EQ. 1 ) THEN
+       IF ( topo_option .EQ. 1 ) THEN
           basin_topoind(ib,ij) = 1000.
        ENDIF
        ! Simple average (as the old method)
-       IF ( option .EQ. 2 ) THEN
+       IF ( topo_option .EQ. 2 ) THEN
           DO iz=1,basin_sz(ij)
              basin_topoind(ib,ij) = basin_topoind(ib,ij) + topoind_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
           ENDDO
           basin_topoind(ib,ij) = basin_topoind(ib,ij)/REAL(basin_sz(ij),r_std)
        ENDIF
        ! Minus outlet hierarchy
-       IF ( option .EQ. 3 ) THEN
+       IF ( topo_option .EQ. 3 .OR. topo_option .EQ. 4 ) 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
@@ -1813,14 +1796,7 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
        ! To make sure that it has the lowest number if this is an outflow point we reset basin_hierarchy
        !
        IF (basin_bxout(ij) .LT. 0 .AND. basin_bxout(ij) .NE. -4) THEN
-          IF ( option .NE. 1 ) THEN
-             basin_hierarchy(ib,ij) = min_topoind
-             !basin_fac(ib,ij) = min_topoind
-          ELSE
              basin_hierarchy(ib,ij) = 0.00
-             !basin_fac(ib,ij) = 0.00
-          ENDIF
-          basin_topoind(ib,ij) = min_topoind
        ENDIF
        !
        !