Skip to content

GitLab

Commit f1cd2338 authored by Anthony Schrapffer's avatar Anthony Schrapffer
Browse files

Correction of Fetch issue

parent 98541ecb
Hide whitespace changes
Inline Side-by-side
Showing with 56 additions and 46 deletions
F90subroutines/routing_interface.f90
View file @ f1cd2338
......@@ -843,8 +843,8 @@ SUBROUTINE checkfetch(nbpt, nwbas, fetch_basin, outflow_grid, outflow_basin, bas
END SUBROUTINE checkfetch
SUBROUTINE finalfetch(nbpt, nbasmax, nbcore, corepts, routing_area, basin_count, route_togrid, route_tobasin, partial_sum, &
SUBROUTINE finalfetch(nbpt, nbasmax, nbcore, nbhalo, in_core, halopts, corepts, &
& routing_area, basin_count, route_togrid, route_tobasin, partial_sum, &
& fetch_basin, outflow_uparea)
!
USE defprec
......@@ -855,6 +855,9 @@ SUBROUTINE finalfetch(nbpt, nbasmax, nbcore, corepts, routing_area, basin_count,
INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
INTEGER(i_std), INTENT(in) :: nbasmax !!
INTEGER(i_std), INTENT (in) :: nbcore
INTEGER(i_std), INTENT (in) :: nbhalo
INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: in_core
INTEGER(i_std), DIMENSION(nbhalo), INTENT(in) :: halopts
INTEGER(i_std), DIMENSION(nbcore), INTENT(in) :: corepts
INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(in) :: routing_area !! Surface of basin (m^2)
......@@ -878,38 +881,40 @@ SUBROUTINE finalfetch(nbpt, nbasmax, nbcore, corepts, routing_area, basin_count,
ENDDO
!
! Propagate first the partial sums from the halo into the domain
! We propagate the HTUs in the halo flowing into the core domain
! until the flow flows out of the domain
!
DO ig=1,nbpt
!
DO ic=1,nbhalo
ig = halopts(ic)
DO ib=1,basin_count(ig)
IF (partial_sum(ig,ib) > EPSILON(partial_sum)) THEN
!
fetch_basin(ig,ib) = MAX(fetch_basin(ig,ib), partial_sum(ig,ib))
!
igrif = route_togrid(ig,ib)
ibasf = route_tobasin(ig,ib)
!
itt = 0
!
DO WHILE ( ibasf .LE. nbasmax )
!
fetch_basin(igrif,ibasf) = MAX(fetch_basin(igrif,ibasf), partial_sum(ig,ib))
!
it = route_togrid(igrif,ibasf)
ibasf = route_tobasin(igrif,ibasf)
igrif = it
itt = itt + 1
IF ( itt .GT. maxiterations) THEN
IF ( itt .GT. maxiterations+20) THEN
CALL ipslerr_p(3,'routing_interface','Problem in propagating partial sum','','')
ENDIF
ENDIF
ENDDO
ENDIF
ENDDO
ENDDO
igrif = route_togrid(ig,ib)
ibasf = route_tobasin(ig,ib)
IF (ibasf .LE. nbasmax) THEN
IF (in_core(igrif) .EQ. 1) THEN
itt = 0
f = partial_sum(ig,ib)
DO WHILE (ibasf .LE. nbasmax)
!
fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + f
!
it = route_togrid(igrif, ibasf)
ibasf = route_tobasin(igrif, ibasf)
igrif = it
IF (ibasf .LE. nbasmax) THEN
IF (in_core(igrif) .EQ. 0) ibasf = nbasmax+1
END IF
itt = itt + 1
IF ( itt .GT. maxiterations) THEN
CALL ipslerr_p(3,'routing_reg_fetch','Problem in propagating partial sum','','')
ENDIF
ENDDO
END IF
END IF
END DO
END DO
!
! Add the areas contributed by the core region of the domain.
! and propagate them till the flow flows out of the core domain.
!
DO ic=1,nbcore
ig = fcorepts(ic)
......@@ -923,25 +928,25 @@ SUBROUTINE finalfetch(nbpt, nbasmax, nbcore, corepts, routing_area, basin_count,
!
itt = 0
!
DO WHILE ( ibasf .LE. nbasmax )
DO WHILE (ibasf .LE. nbasmax)
!
fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + routing_area(ig,ib)
!
it = route_togrid(igrif,ibasf)
ibasf = route_tobasin(igrif,ibasf)
igrif = it
itt = itt + 1
IF (in_core(igrif) .EQ. 1) THEN
fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + routing_area(ig, ib)
!
it = route_togrid(igrif, ibasf)
ibasf = route_tobasin(igrif, ibasf)
igrif = it
itt = itt + 1
ELSE
ibasf = nbasmax + 1
END IF
IF ( itt .GT. maxiterations) THEN
WRITE(numout,&
"('Grid ',I5, ' and basin ',I5, ' did not converge after iteration ',I5)") ig, ib, itt
WRITE(numout, "('==> We flow through grid ',I5,' and basin ',I5)") igrif, ibasf
!
IF ( itt .GT. maxiterations+20) THEN
CALL ipslerr_p(3,'routing_reg_fetch','Problem with computing final fetch','','')
ENDIF
CALL ipslerr_p(3,'routing_reg_fetch','Problem...','','')
ENDIF
ENDDO
!
ENDDO
!
ENDDO
!
! Get the largest upstream areas in this core
......
Interface.py
View file @ f1cd2338
......@@ -184,13 +184,17 @@ def finalfetch(part, routing_area, basin_count, route_togrid, route_tobasin, fet
fetch_out = np.zeros(routing_area.shape, dtype=np.float32, order='F')
partial_sum = np.zeros(routing_area.shape, dtype=np.float32, order='F')
old_sorted = np.zeros(largest_pos, dtype=np.float32, order='F')
#
nbpt = routing_area.shape[0]
fhalolist = [i+1 for i in range(nbpt) if i not in part.landcorelist]
in_core = np.array([1 if i in part.landcorelist else 0 for i in range(nbpt)])
#
maxdiff_sorted = prec*prec
iter_count = 0
#
while iter_count < part.size*3 and maxdiff_sorted > prec :
fetch_out[:,:] = 0.0
outflow_uparea = routing_interface.finalfetch(part.landcorelist, routing_area, basin_count, route_togrid, \
outflow_uparea = routing_interface.finalfetch(in_core,fhalolist, part.landcorelist, routing_area, basin_count, route_togrid, \
route_tobasin, partial_sum, fetch_out)
partial_sum = np.copy(fetch_out)
part.landsendtohalo(partial_sum, order='F')
......@@ -647,7 +651,8 @@ class HydroGraph :
inflow_basin = hydrosuper.inflow_basin)
#
self.routing_fetch = finalfetch(part, self.routing_area, self.route_nbbasin, self.route_togrid, self.route_tobasin, self.routing_fetch)
# This step is no more necessary
#self.routing_fetch = finalfetch(part, self.routing_area, self.route_nbbasin, self.route_togrid, self.route_tobasin, self.routing_fetch)
#
self.num_largest = routing_interface.finalrivclass(part.landcorelist, self.route_togrid, self.route_tobasin, self.routing_fetch, \
hydrosuper.largest_rivarea)
......
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