diff --git a/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index 9ffbac79d6fe6757c76ab1e5cce91913e050f1c7..2a21806882b161f665aba7398420cce234873177 100644
--- a/F90subroutines/routing_reg.f90
+++ b/F90subroutines/routing_reg.f90
@@ -2054,28 +2054,30 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
IF ( outflow_basin(sp,sb) == bop ) THEN
solved(sp,1) = solved(sp,1) + 1
ENDIF
+ WRITE(numout,*) sp,sb,"flow in the same grid, basin:",bop
!
ELSE IF ( outflow_grid(sp,sb) .EQ. -3 ) THEN
! Return flow
! Nothing to do but just remember it is done.
solved(sp,1) = solved(sp,1) + 1
+ WRITE(numout,*) sp,sb,"is a return flow"
ELSE IF ( outflow_grid(sp,sb) .EQ. -2 ) THEN
! Coastal flow
! Nothing to do but just remember it is done.
solved(sp,1) = solved(sp,1) + 1
+ WRITE(numout,*) sp,sb,"is a coastal flow"
ELSE IF ( outflow_grid(sp,sb) .EQ. -1 ) THEN
! River flow
! Nothing to do but just remember it is done.
solved(sp,1) = solved(sp,1) + 1
+ WRITE(numout,*) sp,sb,"is a river outflow"
ENDIF
- WRITE(numout,*) sp,sb,"is an outflow ", outflow_grid(sp,sb)
ELSE IF ( outflow_grid(sp,sb) .GT. 0 ) THEN
found = 0
!
! Deal with the first one as usual
!
inp = outflow_grid(sp,sb)
- IF ( inp == 0.) WRITE(numout,*) "ERROR 1"
!
CALL routing_reg_bestsubbasin(sp, sb, basin_id(sp,sb), basin_hierarchy(sp,sb), basin_fac(sp,sb), &
& basin_flowdir(sp,sb), invented_basins, &
@@ -2223,27 +2225,10 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
! ! loc(NbNeighb)
ENDIF
!
- ! If not found, and at least one of the neighbour is ocean -> coastal outflow
- IF ( found == 0) THEN
- WRITE (numout,*) "Is there a coastal outflow ? "
- nbocean = 0
- DO nb=1,NbNeighb
- IF ( neighbours_g(sp,nb) .LE. -1 ) nbocean = nbocean +1
- ENDDO
- IF ( nbocean .LE. 1 ) THEN
- outflow_grid(sp,sb) = -2
- outflow_basin(sp,sb) = undef_int
- found = 1
- WRITE (numout,*) "There is more than one neighbour that is coastal outflow"
- WRITE (numout,*) "Changed in coastal outflow"
-
- solved(sp,3) = solved(sp,3)+1
- ENDIF
- ENDIF
-
- ! Equivalent to linkup 5
+ ! Looking for a solution in the grid -> HTU with a similar hierarchy or lowest hierarchy
+ !
IF ( found == 0 ) THEN
- WRITE (numout,*) "Eq. Linkup 5"
+ WRITE (numout,*) "Looking for a solution in the grid"
DO sba=1,basin_count(sp)
IF ( sba .NE. sb ) THEN
IF ( basin_id(sp,sb) == basin_id(sp,sba) ) THEN
@@ -2285,47 +2270,71 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
sbint = minhloc(basin_id(sp,sb),2)
ENDIF
ENDIF
- !
- ! Trung: Be careful when changing value of basin_hierarchy here !
- !
- IF (( sbint < undef_int ) .AND. (sbint .GT. 0) ) THEN
+
+ IF (( sbint < undef_int ) .AND. (sbint .GT. 0) .AND. (sbint .NE. sb)) THEN
basin_hierarchy(sp,sb) = basin_hierarchy(sp,sbint)
CALL routing_updateflow(sp, sb, sp, sbint, nbpt, nwbas, nbxmax*8, outflow_grid, outflow_basin, &
- & inflow_number, inflow_grid, inflow_basin)
+ & inflow_number, inflow_grid,inflow_basin)
+ found = 1
WRITE (numout,*) "Lowest basin hierarchy in the grid file"
-
- ELSE
- !
- ! Trung: Well, if you come to this step, it should be a special
- ! (or curious, I said) case. For example, when processing
- ! Sulina branch of river Danube with the forcing data
- ! from E2OFD (/bdd/ORCHIDEE_Forcing/BC/OOL/OL2/E2OFD/)
- ! in resolution of 0.25 degree. The grid box [45.125N, 29.375E]
- ! is a ocean point in the middle of 7 land points ("a
- ! hole"). But HydroSHEDS shows that Sulina branch flows
- ! through this grid cell (from [45.125N, 29.125E] to [45.125N, 29.625E]).
- ! Trung: If the "JUMP" does not work, it means a worse case.
- ! Or it means a worse case.
- ! The outlet of this river falls in the ocean point of
- ! forcing data. For example: with 0.25 degree forcing
- ! data of E2OFD, the grid point [ 45.625N ; 13.875E ]
- ! lost the outlet of basin ID 69666.
- ! We need to use this point (which should be the last point
- ! of this river with shortest hierarchy and highest flow accumulation)
- ! to become new outlet.
- !
- ! This is not good solution because there can be few
- ! outlet for one basin ID.
- ! But I need the model works now !!! => so, come back
- ! here later !
- !
- WRITE (numout,*) "Linkup : Made a NEW OUTLET at sp & sb: ",sp,sb
- ! Coastal flow or river flow is both ok here
- outflow_grid(sp,sb) = -2
- basin_hierarchy(sp,sb) = 0.00
ENDIF
- IF ( outflow_basin(sp,sb) < undef_int ) solved(sp,4) = solved(sp,4)+1
ENDIF
+ !
+ ! Looking if the minimun hierarchy is in a neighbour grid point
+ IF (found .EQ. 0) THEN
+ dop = undef_int
+ bop = undef_int
+ IF ( basin_id(sp,sb) < invented_basins ) THEN
+ DO nb=1,NbNeighb
+ IF (( neighbours_g(sp,order_ref(nb)) == minhloc(basin_id(sp,sb),1)) .AND. (found == 0)) THEN
+ dop = minhloc(basin_id(sp,sb),1)
+ bop = minhloc(basin_id(sp,sb),2)
+ IF ((dop < undef_int ) .AND. (dop .GT. 0) .AND. (bop < undef_int) .AND. (bop .GT. 0)) THEN
+ !basin_hierarchy(sp,sb) = basin_hierarchy(dop,bop)
+ CALL routing_updateflow(sp, sb, dop, bop, nbpt,nwbas, nbxmax*8, outflow_grid, outflow_basin, &
+ & inflow_number, inflow_grid, inflow_basin)
+ found = 1
+ WRITE(numout,*) "Lowest basin hierarchy in the neighbours grid points"
+ WRITE(numout,*) "hierarch of (sp,sb)", basin_hierarchy(sp,sb)
+ WRITE(numout,*) "Lowest basinid hierarch", basin_hierarchy(dop,bop)
+ END IF
+ END IF
+ END DO
+ END IF
+ ENDIF
+ !
+ ! Last option : coastal outflow
+ IF (found .EQ. 0) THEN
+ !
+ ! Trung: Well, if you come to this step, it should be a special
+ ! (or curious, I said) case. For example, when processing
+ ! Sulina branch of river Danube with the forcing data
+ ! from E2OFD (/bdd/ORCHIDEE_Forcing/BC/OOL/OL2/E2OFD/)
+ ! in resolution of 0.25 degree. The grid box [45.125N, 29.375E]
+ ! is a ocean point in the middle of 7 land points ("a
+ ! hole"). But HydroSHEDS shows that Sulina branch flows
+ ! through this grid cell (from [45.125N, 29.125E] to [45.125N, 29.625E]).
+ ! Trung: If the "JUMP" does not work, it means a worse case.
+ ! Or it means a worse case.
+ ! The outlet of this river falls in the ocean point of
+ ! forcing data. For example: with 0.25 degree forcing
+ ! data of E2OFD, the grid point [ 45.625N ; 13.875E ]
+ ! lost the outlet of basin ID 69666.
+ ! We need to use this point (which should be the last point
+ ! of this river with shortest hierarchy and highest flow accumulation)
+ ! to become new outlet.
+ !
+ ! This is not good solution because there can be few
+ ! outlet for one basin ID.
+ ! But I need the model works now !!! => so, come back
+ ! here later !
+ !
+ WRITE (numout,*) "Linkup : Made a NEW OUTLET at sp & sb: ",sp,sb
+ ! Coastal flow or river flow is both ok here
+ outflow_grid(sp,sb) = -2
+ basin_hierarchy(sp,sb) = 0.00
+ ENDIF
+ IF ( outflow_basin(sp,sb) < undef_int ) solved(sp,4) = solved(sp,4)+1
ENDIF
ENDDO
ENDDO