From f8c324bd7d156b3624a4013c0b2028fe78aa2b83 Mon Sep 17 00:00:00 2001
From: Anthony <anthony.schrapffer@polytechnique.fr>
Date: Tue, 5 May 2020 16:04:21 +0200
Subject: [PATCH] Improvement to solve the inflow issue
---
F90subroutines/routing_reg.f90 | 74 +++++++++++++++++++-------------
F90subroutines/routing_tools.f90 | 13 +++---
2 files changed, 52 insertions(+), 35 deletions(-)
diff --git a/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index 1ed64bc..a01c768 100644
--- a/F90subroutines/routing_reg.f90
+++ b/F90subroutines/routing_reg.f90
@@ -2034,6 +2034,8 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
CALL FLUSH(numout)
!
inflow_number(:,:) = 0
+ inflow_basin(:,:,:) = 0
+ inflow_grid(:,:,:) = 0
!
! Variable to keep track of how many points were solved at each stage
!
@@ -2074,7 +2076,14 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
WRITE(numout,*) "*****"
WRITE(numout,*) "Linkup 0 - sp, sb = ", sp, sb
WRITE(numout,*) "Linkup 0 - outflow_grid = ", outflow_grid(sp,sb)
- IF ( outflow_grid(sp,sb) == 0 ) WRITE(numout,*) "Linkup 1.0 -- Flow out of Halo zone"
+ !
+ found = 0
+ IF ( outflow_grid(sp,sb) == 0 ) THEN
+ found = 1
+ WRITE(numout,*) "Linkup 1.0 -- Flow out of Halo zone"
+ ELSE
+ found = 0
+ END IF
!
IF ( outflow_grid(sp,sb) < 0 ) THEN
IF ( outflow_grid(sp,sb) .EQ. -4 ) THEN
@@ -2083,28 +2092,34 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
CALL routing_updateflow(sp, sb, sp, bop, nbpt, nwbas, inflowmax, outflow_grid, outflow_basin, &
& inflow_number, inflow_grid, inflow_basin)
IF ( outflow_basin(sp,sb) == bop ) THEN
+ found = 1
solved(sp,1) = solved(sp,1) + 1
+ WRITE(numout,*) sp,sb,"flow in the same grid, basin:",bop
+ ELSE
+ WRITE(*,*) sp, sb, "flows in the same grid but has an error !"
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.
+ found = 1
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.
+ found = 1
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.
+ found = 1
solved(sp,1) = solved(sp,1) + 1
WRITE(numout,*) sp,sb,"is a river outflow"
ENDIF
- ELSE IF ( outflow_grid(sp,sb) .GT. 0 ) THEN
- found = 0
+ END IF
+ IF ( found .EQ. 0 ) THEN
!
! Deal with the first one as usual
!
@@ -2307,8 +2322,10 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
basin_hierarchy(sp,sb) = basin_hierarchy(sp,sbint)
CALL routing_updateflow(sp, sb, sp, sbint, nbpt, nwbas, inflowmax, outflow_grid, outflow_basin, &
& inflow_number, inflow_grid,inflow_basin)
- found = 1
- WRITE (numout,*) "Lowest basin hierarchy in the grid file"
+ IF (outflow_basin(sp,sb) == sbint) THEN
+ found = 1
+ WRITE (numout,*) "Lowest basin hierarchy in the grid file"
+ END IF
ENDIF
ENDIF
!
@@ -2328,7 +2345,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
! It is possible that the lowest hierarchy is in two grid cells
! In this case, if we have an error for this routing_updateflow
! We go to next step and make it a coastal flow
- IF (outflow_basin(sp, sb) < undef_int) THEN
+ IF (outflow_basin(sp, sb) .EQ. bop) THEN
found = 1
WRITE(numout,*) "Lowest basin hierarchy in the neighbours grid points"
ELSE
@@ -2877,8 +2894,6 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b
igrif = it
ENDDO
!
- !
- !
! Redirect the flows which went into the basin to be killed before we change everything
!
DO inf = 1, inflow_number(ib, tokill)
@@ -2899,29 +2914,30 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b
inb = outflow_basin(ib, tokill)
doshift = .FALSE.
!
- DO inf = 1, inflow_number(ing, inb)
- IF ( doshift ) THEN
- inflow_grid(ing, inb, inf-1) = inflow_grid(ing, inb, inf)
- inflow_basin(ing, inb, inf-1) = inflow_basin(ing, inb, inf)
- ENDIF
+ ! check if tokill outflow has inflow
+ IF (inflow_number(ing,inb) .EQ. 0 ) THEN
+ WRITE(numout,*) 'tokill basin outflow has no inflow'
+ CALL ipslerr_p(3,'routing_reg_killbas','Inflows incoherence - no inflow','','')
+ END IF
+ ! Search tokill in the inflow of its outflow
+ inf = 1
+ DO WHILE ((inf .LE. inflow_number(ing,inb)) .AND. (.NOT. doshift))
IF ( inflow_grid(ing, inb, inf) .EQ. ib .AND. inflow_basin(ing, inb, inf) .EQ. tokill) THEN
doshift = .TRUE.
- ENDIF
- ENDDO
- !
- ! This is only to allow for the last check
- !
- inf = inflow_number(ing, inb)
- IF ( inflow_grid(ing, inb, inf) .EQ. ib .AND. inflow_basin(ing, inb, inf) .EQ. tokill) THEN
- doshift = .TRUE.
- ENDIF
- !
- IF ( .NOT. doshift ) THEN
- WRITE(numout,*) 'Strange we did not find the basin to kill in the downstream basin'
- CALL ipslerr_p(3,'routing_reg_killbas','Basin not found','','')
- ENDIF
- inflow_number(ing, inb) = inflow_number(ing, inb) - 1
+ ELSE
+ inf = inf +1
+ END IF
+ END DO
+
+ IF (inf .GT. inflow_number(ing,inb)) THEN
+ WRITE(numout,*) 'tokill basin is not in the inflows of its outflow'
+ CALL ipslerr_p(3,'routing_reg_killbas','Inflows incoherence','','')
+ END IF
+ ! Inflows are ok -> change them
+ inflow_grid(ing, inb, inf:inflow_number(ing,inb)-1) = inflow_grid(ing, inb, inf+1:inflow_number(ing,inb))
+ inflow_basin(ing, inb, inf:inflow_number(ing,inb)-1) = inflow_basin(ing, inb, inf+1:inflow_number(ing,inb))
+ inflow_number(ing,inb) = inflow_number(ing,inb)-1
ENDIF
!
! Now remove from the arrays the information of basin "tokill"
diff --git a/F90subroutines/routing_tools.f90 b/F90subroutines/routing_tools.f90
index dd65654..cf6f35b 100644
--- a/F90subroutines/routing_tools.f90
+++ b/F90subroutines/routing_tools.f90
@@ -605,9 +605,9 @@ CONTAINS
!
INTEGER(i_std), DIMENSION(pt,bas), INTENT(inout) :: out_grid !! Grid to which we flow
INTEGER(i_std), DIMENSION(pt,bas), INTENT(inout) :: out_basin !! Basin in the out_grid
- INTEGER(i_std), DIMENSION(pt,basmax), INTENT(inout) :: in_number !! Number of basins flowing into this grid
- INTEGER(i_std), DIMENSION(pt,basmax,basmax), INTENT(inout) :: in_grid !! Source grid for this grid
- INTEGER(i_std), DIMENSION(pt,basmax,basmax), INTENT(inout) :: in_basin !! Source basin for this grid
+ INTEGER(i_std), DIMENSION(pt,bas), INTENT(inout) :: in_number !! Number of basins flowing into this grid
+ INTEGER(i_std), DIMENSION(pt,bas,basmax), INTENT(inout) :: in_grid !! Source grid for this grid
+ INTEGER(i_std), DIMENSION(pt,bas,basmax), INTENT(inout) :: in_basin !! Source basin for this grid
!
! LOCAL
!
@@ -631,14 +631,15 @@ CONTAINS
ENDIF
!
IF ( transfer ) THEN
- out_grid(sp,sb) = ep
- out_basin(sp,sb) = eb
- !
in_number(ep,eb) = in_number(ep,eb) + 1
+
IF ( in_number(ep,eb) .LE. basmax ) THEN
+ out_grid(sp,sb) = ep
+ out_basin(sp,sb) = eb
in_grid(ep,eb,in_number(ep,eb)) = sp
in_basin(ep,eb,in_number(ep,eb)) = sb
ELSE
+ in_number(ep,eb) = in_number(ep,eb) - 1
WRITE(numout,*) "routing_updateflow : Too many basins flowing into this point ", ep, eb
CALL ipslerr_p(3,'routing_updateflow', &
"Too many basins flowing into this point ",'Increase nbxmax or its derived quantities','')
--
GitLab