Correction topoindex calculation + add topo_option 7. topo_option is set by default to 4.

F90subroutines/routing_reg.f90

@@ -58,8 +58,10 @@ MODULE routing_reg
   !                   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 
+  ! topo_option = 7 : Calculate rlen/rdz for each pixel (from difference of orog and fac)
+  !                   then average the value for the HTU
   !
-  INTEGER(i_std), SAVE :: topo_option = 3 !! Option to calculate topoindex
+  INTEGER(i_std), SAVE :: topo_option = 4 !! Option to calculate topoindex
   !
   ! Options to compute the properties when merging HTUs in routing_reg_killbas
   ! kill_option = 1 : The old way of simpling building the weighted average.
@@ -440,10 +442,10 @@ CONTAINS
 !! \n
 !_ ================================================================================================================================
 
-  SUBROUTINE routing_reg_findbasins(nb_htu, nbv, ib, ijdimmax, nbi, nbj, nbasmax, trip, basin, fac, hierarchy, topoind, &
-       & rlen, rdz, rweight, lshead, diaglalo, nb_basin, basin_inbxid, basin_outlet,  basin_outtp, basin_sz, &
+  SUBROUTINE routing_reg_findbasins(nb_htu, nbv, ib, ijdimmax, nbi, nbj, nbasmax, trip, basin, fac, hierarchy, topoind_bx, &
+       & rlen_bx, rdz_bx, rweight, lshead, diaglalo, nb_basin, basin_inbxid, basin_outlet,  basin_outtp, basin_sz, &
        & basin_bxout, basin_bbout, basin_pts, basin_lshead, coast_pts, lontmp, lattmp, orog, area, &
-       & basin_topoindex_stream, basin_rlen_stream, basin_dz_stream )
+       & basin_topoindex_stream, basin_rlen_stream, basin_dz_stream)
     !
     IMPLICIT NONE
     !
@@ -465,9 +467,9 @@ CONTAINS
     REAL(r_std), INTENT(inout)    :: fac(ijdimmax,ijdimmax) !!
     REAL(r_std), INTENT(inout)    :: hierarchy(ijdimmax,ijdimmax) !!
     INTEGER(i_std), INTENT(inout) :: basin(ijdimmax,ijdimmax) !!
-    REAL(r_std), INTENT(inout) :: topoind(ijdimmax,ijdimmax) !! Topographic index of the residence time (km)
-    REAL(r_std), INTENT(inout) :: rlen(ijdimmax,ijdimmax) !! River length (m)
-    REAL(r_std), INTENT(inout) :: rdz(ijdimmax,ijdimmax) !! Elevation change (m)
+    REAL(r_std), INTENT(inout) :: topoind_bx(ijdimmax,ijdimmax) !! Topographic index of the residence time (km)
+    REAL(r_std), INTENT(inout) :: rlen_bx(ijdimmax,ijdimmax) !! River length (m)
+    REAL(r_std), INTENT(inout) :: rdz_bx(ijdimmax,ijdimmax) !! Elevation change (m)
     REAL(r_std), INTENT(out)   :: rweight(ijdimmax,ijdimmax) !!
     !
     ! For debug and get coordinate of river outlet
@@ -515,6 +517,10 @@ CONTAINS
     REAL(r_std) :: tlat(nb_htu) !!
     REAL(r_std) :: touttp(nb_htu) !!
     REAL(r_std) :: toutlshead(nb_htu) !!
+    REAL(r_std) :: topoind(ijdimmax,ijdimmax) !!
+    REAL(r_std) :: rlen(ijdimmax,ijdimmax) !! River length (m)
+    REAL(r_std) :: rdz(ijdimmax,ijdimmax) !! Elevation change (m)
+
     !
     INTEGER(i_std) :: tmpsz(nbvmax) !!
     INTEGER(i_std) :: ip, jp, jpp(1), ipb !!
@@ -560,6 +566,8 @@ CONTAINS
     !
     htufrac_min = 1.0/MAX(REAL(nbasmax),4.0)
     htufrac_min = 0.05
+    rlen(:,:) = 0
+    rdz(:,:) = 0
     !
     IF ( debug .AND. routing_diagbox(ib, diaglalo) ) THEN
        WRITE(numout,*) "Point: ", ib
@@ -966,7 +974,7 @@ CONTAINS
         ! WRITE(numout,*) "POINT", ib, ip
         CALL mainstream_topoindex(nb_htu, nbv, nbi, nbj, ip, basin_outloc(ip,1), basin_outloc(ip,2),&
                 & basin_sz(ip), basin_bbout(ip), basin_bxout(ip), basin_lshead(ip), &
-                & basin_pts, trip, basin, fac, lontmp, lattmp, orog, rlen, rdz, area, & 
+                & basin_pts, trip, basin, fac, lontmp, lattmp, orog, rlen_bx, rdz_bx, area, & 
                 & basin_topoindex_stream(ip), basin_rlen_stream(ip), basin_dz_stream(ip))
         IF ( basin_topoindex_stream(ip) .LE. 0 ) THEN 
             ! We may have very small value (due to large slope)
@@ -1017,9 +1025,9 @@ 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(ip,jp) = topoind(ip,jp)+topoind(ill,jll)
-                   rlen(ip,jp) = rlen(ip,jp)+rlen(ill,jll)
-                   rdz(ip,jp) = rdz(ip,jp)+rdz(ill,jll)
+                   topoind(ip,jp) = topoind(ip,jp)+topoind_bx(ill,jll)
+                   rlen(ip,jp) = rlen(ip,jp)+rlen_bx(ill,jll)
+                   rdz(ip,jp) = rdz(ip,jp)+rdz_bx(ill,jll)
                    rivpas(ill,jll) = rivpas(ill,jll)+1
                    rweight(ip,jp) = rweight(ip,jp)+1.0
                    !
@@ -1057,6 +1065,12 @@ CONTAINS
              ENDDO
           ENDIF
        ENDDO
+       ! rlen and rdz are used to make the calculation
+       ! then the results are stored in rlen_bx and rdz_bx
+       ! thus, the difference between the pixel and the outflow 
+       ! replace the local values (between the pixel and the next pixel
+       rlen_bx(:,:) = rlen(:,:)
+       rdz_bx(:,:) = rdz(:,:)
     ENDIF
     !
   END SUBROUTINE routing_reg_findbasins
@@ -2478,7 +2492,7 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
     lowlim_std = hydro_dzprecision
     uplim_std = hydro_dzprecision*200
     !
-    IF ( topo_option .LT. 1 .OR. topo_option .GT. 6 ) THEN
+    IF ( topo_option .LT. 1 .OR. topo_option .GT. 7 ) THEN
        WRITE(numout,*) ' You chose wrong topo_option for calculating topoindex: ', topo_option
        WRITE(numout,*) ' It should be: 1, 2, 3, 4, 5 or 6'
        STOP 'routing_reg_globalize'
@@ -2666,6 +2680,8 @@ 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
+          basin_rlen(ib,ij) = basin_rlen(ib,ij)/cnt
+          basin_rdz(ib,ij) = basin_rdz(ib,ij)/cnt
        ELSE IF ( topo_option .EQ. 5 ) THEN
           cnt = 1
           fac = 0.0
@@ -2678,6 +2694,22 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
                 basin_rdz(ib,ij) = rdz_bx(ip,jp)
              ENDIF
           ENDDO
+       ELSE IF (topo_option .EQ. 7) THEN
+           cnt = basin_sz(ij)
+           IF (cnt .EQ. 1) THEN
+              ! Dans ce cas mettre la valeur de mainstream_topoindex
+              basin_rlen(ib,ij) = rlen_bx(basin_pts(ij,1,1),basin_pts(ij,1,2)) 
+              basin_rdz(ib,ij) = rdz_bx(basin_pts(ij,1,1),basin_pts(ij,1,2))
+           ELSE
+              DO iz=1,basin_sz(ij)
+                 ip = basin_pts(ij,iz,1)
+                 jp = basin_pts(ij,iz,2)
+                 basin_rlen(ib,ij) = basin_rlen(ib,ij) + (hierarchy_bx(ip,jp) - basin_hierarchy(ib,ij)) / basin_sz(ij)
+                 basin_rdz(ib,ij) = basin_rdz(ib,ij) + (orog_bx(ip,jp)-basin_orog_min(ib,ij)) / basin_sz(ij)
+              END DO
+           END IF
+           basin_rdz(ib,ij) = MAX(basin_rdz(ib,ij),0.1)
+           basin_rlen(ib,ij) = MAX(hydro_meanlen,basin_rlen(ib,ij))
        ELSE
           WRITE(numout,*) ' We should not be here for topoindex = ', topo_option
           STOP 'routing_reg_globalize'
@@ -2700,8 +2732,7 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
        ENDIF
        !
        !
-       IF ( topo_option .EQ. 4 .OR. topo_option .EQ. 5 ) THEN
-          
+       IF ( topo_option .EQ. 4 .OR. topo_option .EQ. 5  .OR. topo_option .EQ. 7 ) THEN
           basin_topoind(ib,ij) = SQRT(basin_rlen(ib,ij)**3./(basin_rdz(ib,ij)*1.0e6))
        ENDIF
        !