Some corrections for the topoindex_stream (+ units in km & default value)

F90subroutines/routing_reg.f90

@@ -993,7 +993,23 @@ CONTAINS
        WRITE(numout,*) 'Basin sized :', basin_sz(1:nb_basin)
     ENDIF
     !
-    ! 7.0 If we calculate topoindex with option 3
+    ! 7.0 Calculate the stream topoindex before 8.0
+    !       because rlen/rdz may be changed
+    !
+    basin_topoindex_stream(:) = 0
+    DO ip=1,nb_basin
+        ! Another more efficient options would be 
+        !     to calculate topoindex_stream from disto
+        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
+    !
+    ! 8.0 If we calculate topoindex with option 3
     !     we need to accumulate the topoind_bx along the flow
     !
     DO ij=1,nb_basin
@@ -1072,18 +1088,6 @@ 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
   !
@@ -2068,8 +2072,8 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz,
       dorog = MAX(0.1, rdz(il,jl))
       distance = main_len(1)
    END IF
-  
-   topoindex_stream = SQRT(distance **3 / dorog)
+   ! topoindex is in km while distance and dorog are in meters
+   topoindex_stream = SQRT(distance **3 / dorog) / 1000
 !
 END SUBROUTINE mainstream_topoindex
 
@@ -3343,6 +3347,9 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
          ! To define in function of the minimum topoindex on the map
          ! (to revise)
          basin_topoind(ib,1) = SQRT(hydro_meanlen**3./(hydro_meandz*1.0e6))
+         ! When isolated artificial basin -> topoindex stream has no use
+         ! This is just not to keep it at 0
+         basin_topoindex_stream(ib,1) = SQRT(hydro_meanlen**3./(hydro_meandz*1.0e6))
          basin_rlen(ib,1) = hydro_meanlen
          basin_rdz(ib,1) = hydro_meandz
          basin_id(ib,1) = 0