Truncate : anticipating loop errors and handling HTUs flowing out of the domain

F90subroutines/routing_interface.f90

@@ -581,7 +581,7 @@ SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numo
   INTEGER(i_std), DIMENSION(nbpt,nbxmax_in,nbxmax_in), INTENT(inout) :: inflow_grid !!
 
   !! LOCAL
-  INTEGER(i_std) :: ib, op, tok, totak
+  INTEGER(i_std) :: ib, op, tok, totak, igrif, ibasf
 
 
   DO ib=1,nbpt
@@ -591,6 +591,22 @@ SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numo
          totak = totakeover(ib,op)
          IF (tok .GT. 0) THEN 
             WRITE(numout,*) "nbpt", ib, "tokill", tok, "totakover", totak
+            ! Test if tokill is downstream of totakeover (avoid loop)
+            igrif = outflow_grid(ib,totak)
+            ibasf = outflow_basin(ib,totak)
+            DO WHILE ( igrif >0 )
+              IF ((igrif .EQ. ib) .AND. (ibasf .EQ. tok)) THEN
+                !CALL ipslerr_p(3,'killbas','tokill is downstream totakeover','','')
+                igrif = 0
+                it = totak
+                totak = tok
+                tok = it
+              ELSE
+                it = outflow_grid(igrif,ibasf)
+                ibasf = outflow_basin(igrif,ibasf)            
+                igrif = it
+              END IF
+            END DO
             CALL routing_reg_killbas(nbpt, ib, tok, totak, nwbas, basin_count, basin_area, basin_cg, basin_topoind,&
                & fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, inflow_number,&
                & inflow_grid, inflow_basin)

F90subroutines/routing_reg.f90

@@ -1059,7 +1059,7 @@ CONTAINS
     !
 !! LOCAL VARIABLES
     REAL(r_std), PARAMETER :: flag=-9999. !!
-    LOGICAL, PARAMETER :: debug=.TRUE.
+    LOGICAL, PARAMETER :: debug=.FALSE.
     LOGICAL, PARAMETER :: checkgrid=.FALSE.
     CHARACTER(LEN=7) :: fmt !!
     CHARACTER(LEN=9) :: fmtr !!