Generate topoindex_stream

F90subroutines/routing_interface.f90

@@ -128,7 +128,7 @@ END SUBROUTINE gethydrogrid
 
 SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, trip_bx, basin_bx, fac_bx, hierarchy_bx, topoind_bx, &
      &                rlen_bx, rdz_bx, rweight_bx, lshead_bx, nb_basin, basin_inbxid, basin_outlet, basin_outtp, &
-     &                basin_sz, basin_bxout, basin_bbout, basin_pts, basin_lshead, coast_pts, lontmp, lattmp)
+     &                basin_sz, basin_bxout, basin_bbout, basin_pts, basin_lshead, coast_pts, lontmp, lattmp, orog_bx)
   !
   USE ioipsl
   USE grid
@@ -149,6 +149,7 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, trip_bx, basin_bx,
   REAL, INTENT(inout)    :: rdz_bx(nbpt,ijdimmax,ijdimmax) !!
   REAL, INTENT(inout)    :: rweight_bx(nbpt,ijdimmax,ijdimmax) !! Weight of each river within the HTU
   REAL, INTENT(in)       :: lshead_bx(nbpt,ijdimmax,ijdimmax) !! Large scale heading for outflow points.
+  REAL, INTENT(in)       :: orog_bx(nbpt,ijdimmax,ijdimmax) !! 
   !
   !
   ! OUTPUT VARIABLES
@@ -182,7 +183,8 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, trip_bx, basin_bx,
           & basin_bx(ib,:,:), fac_bx(ib,:,:), hierarchy_bx(ib,:,:), topoind_bx(ib,:,:), &
           & rlen_bx(ib,:,:), rdz_bx(ib,:,:), rweight_bx(ib,:,:), lshead_bx(ib,:,:), diaglalo, nb_basin(ib), &
           & basin_inbxid(ib,:), basin_outlet(ib,:,:), basin_outtp(ib,:), basin_sz(ib,:), basin_bxout(ib,:), &
-          & basin_bbout(ib,:), basin_pts(ib,:,:,:), basin_lshead(ib,:), coast_pts(ib,:), lontmp(ib,:,:), lattmp(ib,:,:))
+          & basin_bbout(ib,:), basin_pts(ib,:,:,:), basin_lshead(ib,:), coast_pts(ib,:), lontmp(ib,:,:), lattmp(ib,:,:),&
+          & orog_bx(ib,:,:))
   ENDDO
   
 END SUBROUTINE findbasins

F90subroutines/routing_reg.f90

@@ -427,7 +427,7 @@ CONTAINS
 
   SUBROUTINE routing_reg_findbasins(nb_htu, nbv, ib, ijdimmax, nbi, nbj, trip, basin, fac, hierarchy, topoind, &
        & rlen, rdz, 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)
+       & basin_bxout, basin_bbout, basin_pts, basin_lshead, coast_pts, lontmp, lattmp, orog) !basin_topoindex_stream
     !
     IMPLICIT NONE
     !
@@ -440,6 +440,7 @@ CONTAINS
     REAL(r_std), INTENT(in) :: hierarchy(ijdimmax,ijdimmax) !!
     REAL(r_std), INTENT(in) :: fac(ijdimmax,ijdimmax) !!
     REAL(r_std), INTENT(in) :: lshead(ijdimmax,ijdimmax)
+    REAL(r_std), INTENT(in) :: orog(ijdimmax,ijdimmax)
     REAL(r_std), DIMENSION(:,:), INTENT(in) :: diaglalo !! Point (in Lat/Lon) where diagnostics will be printed.
     !
     ! Modified
@@ -460,6 +461,7 @@ CONTAINS
     INTEGER(i_std), INTENT(out) :: nb_basin               !! Number of sub-basins (unitless)
     INTEGER(i_std), INTENT(out) :: basin_inbxid(nb_htu)   !!
     REAL(r_std), INTENT(out)    :: basin_outlet(nb_htu,2) !!
+    INTEGER(i_std)              :: basin_outloc(nb_htu,2) !!
     REAL(r_std), INTENT(out)    :: basin_outtp(nb_htu)    !!
     INTEGER(i_std), INTENT(out) :: basin_sz(nb_htu)       !!
     INTEGER(i_std), INTENT(out) :: basin_bxout(nb_htu)    !!
@@ -467,6 +469,7 @@ CONTAINS
     INTEGER(i_std), INTENT(out) :: basin_bbout(nb_htu)    !!
     INTEGER(i_std), INTENT(out) :: basin_pts(nb_htu, nbv, 2) !!
     INTEGER(i_std), INTENT(out) :: coast_pts(nb_htu)      !! The coastal flow points (unitless)
+    REAL(r_std)   :: basin_topoindex_stream(nb_htu) !! INTENT(out)
     !
 !! LOCAL VARIABLES
     LOGICAL, PARAMETER :: debug=.FALSE.
@@ -906,9 +909,6 @@ CONTAINS
           tmpsz(jpp(1)) = -1
        ENDIF
     ENDDO
-    WRITE(*,*) "nb_basin, nbb", nb_basin, nbb
-    WRITE(*,*) "sortind", sortind(1:nb_basin)
-    WRITE(*,*) "tsz", tsz(sortind(1:nb_basin))
 
     basin_inbxid(1:nb_basin) = tbname(sortind(1:nb_basin))
     basin_outlet(1:nb_basin,1) = tlat(sortind(1:nb_basin))
@@ -918,6 +918,8 @@ CONTAINS
     basin_pts(1:nb_basin,:,:) = tpts(sortind(1:nb_basin),:,:)
     basin_bxout(1:nb_basin) = toutdir(sortind(1:nb_basin))
     basin_lshead(1:nb_basin) = toutlshead(sortind(1:nb_basin))
+    basin_outloc(1:nb_basin,1) = toutloc(sortind(1:nb_basin), 1)
+    basin_outloc(1:nb_basin,2) = toutloc(sortind(1:nb_basin), 2)
     !
     ! Correct toutbas after sorting
     !
@@ -1071,6 +1073,18 @@ CONTAINS
           ENDIF
        ENDDO
     ENDIF
+    
+    ! topoindex_stream à initialiser.
+    basin_topoindex_stream(:) = 0
+    DO ip=1,nb_basin
+        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, basin_topoindex_stream(ip))
+        IF ( basin_topoindex_stream(ip) .LT. 1 ) THEN
+            WRITE(numout,*) basin_topoindex_stream(ip)
+            CALL ipslerr_p(3,'routing_reg_findbasins','Error in the mainstream topoindex','','')
+        END IF 
+    END DO
     !
   END SUBROUTINE routing_reg_findbasins
   !
@@ -1822,6 +1836,247 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz,
    !
    END SUBROUTINE routing_reg_divbas_divide
    ! 
+!! ================================================================================================================================
+!! SUBROUTINE : mainstream_topoindex
+!!
+!>\BRIEF This subroutine calculate the topoindex over the main stream
+!!
+!! DESCRIPTION (definitions, functional, design, flags) : None
+!!
+!! RECENT CHANGE(S): None
+!!
+!! MAIN OUTPUT VARIABLE(S):
+!! topoindex_stream
+!!
+!! REFERENCES : None
+!!
+!! FLOWCHART : None
+!! \n
+!_ ================================================================================================================================
+  !
+  SUBROUTINE mainstream_topoindex(nb_htu, nbv, nbi, nbj, ibas, iloc, jloc, tsz, tout, toutd, headd, &
+      & tpts, trip, basin, fac, lon, lat, orog, rlen, rdz, topoindex_stream)
+   !
+   IMPLICIT NONE
+   !
+   !! INPUT VARIABLES
+   INTEGER(i_std), INTENT(in) :: nb_htu, nbv
+   INTEGER(i_std), INTENT(in) :: nbi !! Number of point in x within the grid (unitless)
+   INTEGER(i_std), INTENT(in) :: nbj !! Number of point in y within the grid (unitless)
+   INTEGER(i_std), INTENT(in) :: ibas !! Order of the basin will be divided
+   INTEGER(i_std), INTENT(in) :: iloc, jloc !! Outlet location
+   !
+   INTEGER(i_std), INTENT(in) :: tsz !! Size of sub-basin
+   INTEGER(i_std), INTENT(in) :: tout !! Outlet type
+   INTEGER(i_std), INTENT(in) :: toutd !! Outlet direction
+   REAL(r_std), INTENT(in)    :: headd
+   INTEGER(i_std), INTENT(in) :: tpts(:,:,:) !!
+   INTEGER(i_std), INTENT(in) :: basin(:,:) !!
+   REAL(r_std), INTENT(in)    :: fac(:,:) !!
+   REAL(r_std), INTENT(in)    :: orog(:,:), rlen(:,:), rdz(:,:) !!
+   REAL(r_std), INTENT(in)    :: lon(:,:), lat(:,:) !!
+   INTEGER(i_std), INTENT(in) :: trip(:,:) !!
+
+   !
+   !! OUTPUT VARIABLES
+   REAL(r_std), INTENT(out) :: topoindex_stream !!
+   !
+   !! LOCAL VARIABLES
+   REAL(r_std), PARAMETER :: flag=-9999. !!
+   LOGICAL, PARAMETER :: debug=.FALSE.
+   LOGICAL, PARAMETER :: checkgrid=.FALSE.
+   CHARACTER(LEN=7) :: fmt !!
+   CHARACTER(LEN=9) :: fmtr !!
+   CHARACTER(LEN=11) :: afmt !!
+   CHARACTER(LEN=13) :: afmtr !!
+   !
+   REAL(r_std), DIMENSION(nbi,nbj) :: factmp_glo !!
+   REAL(r_std) :: distance, dorog !!   
+   INTEGER(i_std), DIMENSION(nbi,nbj) :: triptmp, triptemp !!
+   !
+   INTEGER(i_std) :: il, jl, ill, jll, jp !!
+   INTEGER(i_std) :: ie, je !!
+   INTEGER(i_std) :: p, cnt, k, l, ic!!
+   INTEGER(i_std) :: ip, isz !!
+   INTEGER(i_std) :: ff(1) !!
+   !
+   LOGICAL :: okpoint=.TRUE.
+   !
+   ! Number of neighbours on the HydroShed grid (regular Lat/Lon)
+   !
+   INTEGER(i_std), PARAMETER :: nbne=8
+   INTEGER(i_std), DIMENSION(nbne,2) :: inc !!
+   !
+   INTEGER(i_std) :: main_loc(nbv,2), tri_loc(nbne,nbv,2)
+   INTEGER(r_std) :: main_len(nbv), main_orog(nbv), main_dz(nbv)
+   REAL(r_std) :: main_fac(nbv), tri_fac(nbne,nbv)
+   REAL(r_std) :: tmptri_fac(nbne)  ! Flow accumulation of all neighbour points
+   INTEGER(i_std) :: sortedtrifac(nbne) ! And sort of tmptri_fac(nbne)
+   REAL(r_std) :: alltri_fac(nbv) ! Sort flow accumulation of all tributaries
+
+   !_ ================================================================================================================================
+   !
+   ! The routing code (i=1, j=2)
+   !
+   inc(1,1) = 0
+   inc(1,2) = -1
+   inc(2,1) = 1
+   inc(2,2) = -1
+   inc(3,1) = 1
+   inc(3,2) = 0
+   inc(4,1) = 1
+   inc(4,2) = 1
+   inc(5,1) = 0
+   inc(5,2) = 1
+   inc(6,1) = -1
+   inc(6,2) = 1
+   inc(7,1) = -1
+   inc(7,2) = 0
+   inc(8,1) = -1
+   inc(8,2) = -1
+   !
+   ! 0.0 Get the information of all the points belongs to the subbasin
+   ! which need to divide. Flag for flow accumulation is -9999 (the value of
+   ! flow accumulation should be always non negative).
+   !
+   factmp_glo(:,:) = flag
+   triptmp(:,:) = -1
+   triptemp(:,:) = -1
+   !
+   DO isz = 1, tsz
+      il = tpts(ibas, isz, 1)
+      jl = tpts(ibas, isz, 2)
+      !
+      factmp_glo(il,jl) = fac(il,jl)
+      triptmp(il,jl) = trip(il,jl)
+      triptemp(il,jl) = trip(il,jl)
+   ENDDO
+
+   ! 1.0 If the first version of this subroutine doesn't work well,
+   ! we need to write the validation step of input arrays here (for
+   ! example, the input sub-basin should only has 1 outlet).
+   ! (October 2016: it still works well :))
+   !
+   IF ( triptmp(iloc,jloc) .LT. -4 .OR. triptmp(iloc,jloc) .GE. 109 ) THEN
+      CALL ipslerr_p(3,'routing_reg_divbas','The TRIP of outlet is wrong','','')
+   ENDIF
+   IF ( debug ) THEN
+      WRITE(numout,*) 'Basin number: ', ibas
+      WRITE(numout,*) 'Basin size: ', tsz
+      WRITE(numout,*) 'Basin out: ', tout, toutd
+   ENDIF
+   !
+   ! 2.0 
+   !
+   main_loc(:,:) = 0
+   main_fac(:) = flag
+   main_len(:) = 0
+   main_orog(:) = 0
+   main_dz(:) = 0
+   tri_loc(:,:,:) = 0
+   tri_fac(:,:) = flag
+   !
+   ! Start!
+   !
+   cnt = 0
+   okpoint = .TRUE.
+   il = iloc ; jl = jloc
+   !
+   DO WHILE ( il .GT. 0 .AND. il .LE. nbi .AND. jl .GT. 0 .AND. jl .LE. nbj &
+      & .AND. cnt .LT. nbi*nbj .AND. okpoint )
+      !
+      ! Count the length of the mainstream
+      !
+      cnt = cnt + 1
+      main_loc(cnt,1) = il
+      main_loc(cnt,2) = jl
+      main_len(cnt) =  rlen(il,jl) ! Vaut toujours 0 ? 
+      main_orog(cnt) = orog(il,jl)
+      main_dz(cnt) = rdz(il,jl)
+      main_fac(cnt) = factmp_glo(il,jl)
+      !
+      ! Look for [nbne] neighbour points to find the tributaries
+      ! October 2016: so far, nbne = 8 (for 8 directions of TRIP)
+      !
+      l = 0
+      DO k = 1, nbne
+         ill = il + inc(k,1)
+         jll = jl + inc(k,2)
+         ! If this point is still in the grid box
+         IF ( ill .GT. 0 .AND. ill .LE. nbi .AND. jll .GT. 0 .AND. jll .LE. nbj ) THEN
+            ! If this point still belongs to the current sub-basin
+            IF ( triptmp(ill,jll) .GT. 0 .AND. triptmp(ill,jll) .LE. nbne ) THEN
+               ie = ill + inc(INT(triptmp(ill,jll)),1)
+               je = jll + inc(INT(triptmp(ill,jll)),2)
+               ! If this neighbour point contributes water for current point
+               IF ( ie .EQ. il .AND. je .EQ. jl ) THEN
+                  ! Store all points here
+                  l = l + 1
+                  tri_fac(l,cnt) = factmp_glo(ill,jll)
+                  tri_loc(l,cnt,1) = ill
+                  tri_loc(l,cnt,2) = jll
+                  ! Mark the processed points here
+                  triptmp(ill,jll) = -2
+               ENDIF
+            ENDIF
+         ENDIF
+      ENDDO
+      !
+      ! Remettre la partie suivante
+      !
+      IF ( l > 0 ) THEN
+         !
+         ! Sort the FAC of all neighbour points
+         tmptri_fac(:) = flag
+         sortedtrifac(:) = 0
+         !
+         tmptri_fac(:) = tri_fac(:,cnt)
+         !tmptrifac -> local
+         DO k = 1, l
+            ff = MAXLOC(tmptri_fac)
+            sortedtrifac(k) = ff(1)
+            tmptri_fac(ff(1)) = flag
+         ENDDO
+         ! Continue trace upstream: the neighbour point with highest flow
+         ! accumulation will belongs to main stream. We move to this point
+         ! (change value of il, jl).
+         IF ( tri_fac(sortedtrifac(1),cnt) .NE. flag ) THEN
+            il = tri_loc(sortedtrifac(1),cnt,1)
+            jl = tri_loc(sortedtrifac(1),cnt,2)
+            !
+            tri_fac(sortedtrifac(1),cnt) = flag
+            !
+         ELSE 
+            okpoint = .FALSE.
+         END IF
+      ELSE
+         okpoint = .FALSE.
+      END IF
+   ENDDO
+
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!
+   ! Main river output
+   il = main_loc(1,1); jl = main_loc(1,2)
+   ! Main river - most upstream pixel 
+   ill = main_loc(cnt,1); jll = main_loc(cnt,2)
+
+   IF (cnt .GT. 1) THEN
+      dorog = MAXVAL(main_orog(:cnt)) - MINVAL(main_orog(:cnt))
+      dorog = MAX(0.1, dorog)
+   ELSE
+      ! if the main river has only one pixel we use directly rdz 
+      dorog = MAX(0.1, rdz(il,jl))
+   END IF
+  
+   ! For the distance : sum rlen along the main river 
+   DO l = 1,cnt
+      distance = SUM(main_len(:cnt))
+   END DO 
+
+   topoindex_stream = SQRT(distance **3 / dorog)
+!
+END SUBROUTINE mainstream_topoindex
+
 !! ================================================================================================================================
 !! SUBROUTINE : routing_reg_globalize
 !!

Interface.py

@@ -179,7 +179,7 @@ class HydroSuper :
                                              hierarchy_bx = hydrooverlap.hierarchy_bx, topoind_bx = hydrooverlap.topoind_bx, \
                                              rlen_bx = hydrooverlap.rlen_bx, rdz_bx = hydrooverlap.rdz_bx, \
                                              rweight_bx = hydrooverlap.rweight_bx, lshead_bx = hydrooverlap.lshead_bx, \
-                                             lontmp = hydrooverlap.lon_bx, lattmp = hydrooverlap.lat_bx)
+                                             lontmp = hydrooverlap.lon_bx, lattmp = hydrooverlap.lat_bx, orog_bx = hydrooverlap.orog_bx)
         #
         # Adjust nwbas to the maximum found over the domain
         #