Merge branch 'master' of gitlab.in2p3.fr:ipsl/lmd/intro/routingpp

F90subroutines/routing_reg.f90

@@ -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
 !! ================================================================================================================================
@@ -1018,7 +1021,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)
@@ -1039,22 +1042,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)
@@ -1924,9 +1931,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
     !
@@ -2087,8 +2094,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)
@@ -2100,11 +2108,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)
@@ -2139,7 +2147,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