Skip to content

GitLab

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

Optimization of the fetch calculation

parent f36a0022
Hide whitespace changes
Inline Side-by-side
Showing with 92 additions and 33 deletions
F90subroutines/routing_interface.f90
View file @ cad03e7b
......@@ -333,7 +333,7 @@ SUBROUTINE areanorm(nbpt, nwbas, basin_count, basin_area, outflow_grid, basin_ar
!
END SUBROUTINE areanorm
SUBROUTINE fetch(nbpt, nwbas, nbcore, corepts, basin_count, basin_area, basin_id, basin_hierarchy, basin_fac, &
SUBROUTINE fetch(nbpt, nwbas, nbcore, nbhalo, fhalopts, corepts, basin_count, basin_area, basin_id, basin_hierarchy, basin_fac, &
& outflow_grid, outflow_basin, partial_sum, fetch_basin, outflow_uparea)
!
USE ioipsl
......@@ -347,6 +347,8 @@ SUBROUTINE fetch(nbpt, nwbas, nbcore, corepts, basin_count, basin_area, basin_id
INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
INTEGER(i_std), INTENT (in) :: nwbas !!
INTEGER(i_std), INTENT (in) :: nbcore
INTEGER(i_std), INTENT (in) :: nbhalo
INTEGER(i_std), DIMENSION(nbhalo), INTENT(in) :: fhalopts
INTEGER(i_std), DIMENSION(nbcore), INTENT(in) :: corepts
INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_area !!
......@@ -372,7 +374,8 @@ SUBROUTINE fetch(nbpt, nwbas, nbcore, corepts, basin_count, basin_area, basin_id
fcorepts(ic) = corepts(ic)+1
ENDDO
!
CALL routing_reg_fetch(nbpt, area, contfrac, nwbas, nbcore, fcorepts, basin_count, basin_area, basin_id, &
CALL routing_reg_fetch(nbpt, area, contfrac, nwbas, nbcore,nbhalo, fhalopts, &
& fcorepts, basin_count, basin_area, basin_id, &
& basin_hierarchy, basin_fac, outflow_grid, outflow_basin, partial_sum, fetch_basin, outflow_uparea)
!
......
F90subroutines/routing_reg.f90
View file @ cad03e7b
......@@ -2417,7 +2417,8 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
!! \n
!_ ================================================================================================================================
SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, basin_count, basin_area, basin_id, &
SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, nbhalo, &
& halopts, corepts, basin_count, basin_area, basin_id, &
& basin_hierarchy, basin_fac, outflow_grid, outflow_basin, partial_sum, fetch_basin, outflow_uparea)
!
IMPLICIT NONE
......@@ -2430,6 +2431,8 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, b
!
INTEGER(i_std) :: nwbas !!
INTEGER(i_std), INTENT (in) :: nbcore
INTEGER(i_std), INTENT (in) :: nbhalo
INTEGER(i_std), DIMENSION(nbcore), 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,nwbas), INTENT(in) :: basin_area !!
......@@ -2448,45 +2451,59 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, b
INTEGER(i_std) :: ig, ib, ic, ff(1), it, itt, igrif, ibasf, nboutflow !! Indices (unitless)
REAL(r_std), DIMENSION(nbpt*nbvmax) :: tmp_area !!
INTEGER(i_std), DIMENSION(nbpt*nbvmax,2) :: tmpindex !!
REAL(r_std) :: f
!
INTEGER(i_std), PARAMETER :: maxiterations=10000
INTEGER(i_std), DIMENSION(nbpt) :: in_core
LOGICAL :: inhalo
!
!_ ================================================================================================================================
!
! Fill array with 1 where nbpt is a core point
!
in_core(:) = 0
DO ig = 1, nbpt
CALL isin_halo(ig, nbhalo, nbpt, halopts, inhalo)
IF (.NOT. inhalo) in_core(ig) = 1
END DO
!
! 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 = outflow_grid(ig,ib)
ibasf = outflow_basin(ig,ib)
!
itt = 0
!
DO WHILE (igrif .GT. 0)
igrif = outflow_grid(ig,ib)
IF (igrif .GT. 0) THEN
IF (in_core(igrif) .EQ. 1) THEN
itt = 0
ibasf = outflow_basin(ig,ib)
f = partial_sum(ig,ib)
DO WHILE (igrif .GT. 0)
!
fetch_basin(igrif,ibasf) = MAX(fetch_basin(igrif,ibasf), partial_sum(ig, ib))
fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + f
!
it = outflow_grid(igrif, ibasf)
ibasf = outflow_basin(igrif, ibasf)
igrif = it
IF (igrif .GT. 0) THEN
IF (in_core(igrif) .EQ. 0) igrif = 0
END IF
itt = itt + 1
IF ( itt .GT. maxiterations) THEN
IF ( itt .GT. maxiterations+20) THEN
CALL ipslerr_p(3,'routing_reg_fetch','Problem in propagating partial sum','','')
ENDIF
IF ( itt .GT. maxiterations+20) THEN
CALL ipslerr_p(3,'routing_reg_fetch','Problem in propagating partial sum','','')
ENDIF
ENDIF
ENDDO
ENDIF
ENDDO
ENDDO
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 = corepts(ic)
......@@ -2502,12 +2519,16 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, b
!
DO WHILE (igrif .GT. 0 )
!
fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + basin_area(ig, ib)
!
it = outflow_grid(igrif, ibasf)
ibasf = outflow_basin(igrif, ibasf)
igrif = it
itt = itt + 1
IF (in_core(igrif) .EQ. 1) THEN
fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + basin_area(ig, ib)
!
it = outflow_grid(igrif, ibasf)
ibasf = outflow_basin(igrif, ibasf)
igrif = it
itt = itt + 1
ELSE
igrif = 0
END IF
IF ( itt .GT. maxiterations) THEN
WRITE(numout,&
"('Grid ',I5, ' and basin ',I5, ' did not converge after iteration ',I5)") ig, ib, itt
......@@ -2558,7 +2579,40 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, b
!
END SUBROUTINE routing_reg_fetch
!
!
!
!! ================================================================================================================================
!! SUBROUTINE : isin_halo
!!
!>\BRIEF This subroutine establishes the list of the points
! and indicates whether if they are in the core (=1) or in the halo (=0).
SUBROUTINE isin_halo(ig, nbhalo, nbpt, halopts, isinhalo)
IMPLICIT NONE
! INPUT
INTEGER(i_std), INTENT (in) :: ig
INTEGER(i_std), INTENT (in) :: nbhalo
INTEGER(i_std), DIMENSION(nbhalo), INTENT(in) :: halopts
INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
! OUTPUT
LOGICAL(i_std), INTENT(out) :: isinhalo
! LOCAL
LOGICAL(i_std) :: loop
INTEGER(i_std) :: nb
isinhalo = .FALSE.
loop = .TRUE.
DO nb=1,nbhalo
IF (loop) THEN
IF (halopts(nb) .EQ. ig) THEN
isinhalo = .TRUE.
loop = .FALSE.
END IF
END IF
END DO
END SUBROUTINE isin_halo
!! ================================================================================================================================
!! SUBROUTINE : routing_reg_truncate
!!
......
Interface.py
View file @ cad03e7b
......@@ -392,9 +392,11 @@ class HydroSuper :
maxdiff_sorted = prec*prec
iter_count = 0
#
fhalolist = [i+1 for i in range(self.nbpt) if i not in part.landcorelist]
#
while iter_count < part.size*3 and maxdiff_sorted > prec :
fetch_basin[:,:] = 0.0
outflow_uparea = routing_interface.fetch(part.landcorelist, self.basin_count, self.basin_area, self.basin_id, self.basin_hierarchy, \
outflow_uparea = routing_interface.fetch(fhalolist, part.landcorelist, self.basin_count, self.basin_area, self.basin_id, self.basin_hierarchy, \
self.basin_fac, self.outflow_grid, self.outflow_basin, partial_sum, fetch_basin)
partial_sum = np.copy(fetch_basin)
part.landsendtohalo(partial_sum, order='F')
......
Projections.py
View file @ cad03e7b
......@@ -8,7 +8,7 @@ import configparser
config=configparser.ConfigParser()
config.read("run.def")
#
EarthRadius=config.getfloat("OverAll", "EarthRadius")
EarthRadius=config.getfloat("OverAll", "EarthRadius", fallback=6370000.0)
rad_per_deg = np.pi/180.
deg_per_rad = 180./np.pi
#
......
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