Commit f8c324bd authored by Anthony's avatar Anthony
Browse files

Improvement to solve the inflow issue

parent e1d15ce4
......@@ -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"
......
......@@ -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','')
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment