Skip to content

GitLab

Commit d6b18e5a authored by POLCHER Jan's avatar POLCHER Jan 🚴🏾
Browse files

Improves the treatment of "undefined HTUs". There was a lack of consistency... 

Improves the treatment of "undefined HTUs". There was a lack of consistency between the Python and F90 part. This is now clarified and the routing-graph.nc file can more easily be processed.
parent 3a5c94dc
Hide whitespace changes
Inline Side-by-side
Showing with 85 additions and 42 deletions
F90subroutines/routing_interface.f90
View file @ d6b18e5a
......@@ -440,7 +440,8 @@ SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridare
& inflow_number, inflow_grid, inflow_basin, floodcri, routing_area, routing_orog_mean, &
& routing_orog_min, routing_orog_max, routing_floodp, routing_beta,&
& routing_cg, topo_resid, route_nbbasin, route_togrid, route_tobasin, route_nbintobas, &
& global_basinid, route_outlet, route_type, origin_nbintobas, routing_fetch,routing_floodcri)
& global_basinid, route_outlet, route_type, origin_nbintobas, routing_fetch,routing_floodcri, &
& rfillval, ifillval)
!
USE ioipsl
USE grid
......@@ -476,6 +477,9 @@ SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridare
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless)
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_basin !!
!
REAL(r_std), INTENT(in) :: rfillval
INTEGER(i_std), INTENT(in) :: ifillval
!
!
!
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: inflow_number !!
......@@ -505,33 +509,56 @@ SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridare
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: routing_floodcri
!
!
INTEGER(i_std) :: ib, ij
!
CALL routing_reg_end_truncate(nbpt, nbasmax, gridarea, cfrac, gridcenters, nwbas, inflowmax, num_largest, &
& basin_count, basin_notrun, basin_area, basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp,&
& basin_cg, basin_topoind, basin_beta_fp, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, &
& outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin, floodcri)
routing_area(:,:) = routing_area_glo(:,:)
routing_orog_mean(:,:) = routing_orog_mean_glo(:,:)
routing_orog_min(:,:) = routing_orog_min_glo(:,:)
routing_orog_max(:,:) = routing_orog_max_glo(:,:)
routing_floodp(:,:) = routing_floodp_glo(:,:)
routing_cg(:,:,:) = routing_cg_glo(:,:,:)
topo_resid(:,:) = topo_resid_glo(:,:)
route_nbbasin(:) = route_count_glo(:)
route_togrid(:,:) = route_togrid_glo(:,:)
route_tobasin(:,:) = route_tobasin_glo(:,:)
routing_floodcri(:,:) = routing_floodcri_glo(:,:)
routing_fetch(:,:) = route_fetch_glo(:,:)
routing_beta(:,:) = routing_beta_glo(:,:)
route_nbintobas(:) = route_nbintobas_glo(:)
global_basinid(:,:) = global_basinid_glo(:,:)
route_outlet(:,:,:) = route_outlet_glo(:,:,:)
route_type(:,:) = route_type_glo(:,:)
origin_nbintobas(:) = origin_nbintobas_glo(:)
routing_area(:,:) = rfillval
routing_orog_mean(:,:) = rfillval
routing_orog_min(:,:) = rfillval
routing_orog_max(:,:) = rfillval
routing_floodp(:,:) = rfillval
routing_cg(:,:,:) = rfillval
topo_resid(:,:) = rfillval
route_nbbasin(:) = ifillval
route_togrid(:,:) = ifillval
route_tobasin(:,:) = ifillval
routing_fetch(:,:) = rfillval
routing_beta(:,:) = rfillval
route_nbintobas(:) = ifillval
global_basinid(:,:) = ifillval
route_outlet(:,:,:) = rfillval
route_type(:,:) = rfillval
origin_nbintobas(:) = ifillval
DO ij=1,nbpt
DO ib=1,basin_count(ij)
!
routing_area(ij,ib) = routing_area_glo(ij,ib)
routing_orog_mean(ij,ib) = routing_orog_mean_glo(ij,ib)
routing_orog_min(ij,ib) = routing_orog_min_glo(ij,ib)
routing_orog_max(ij,ib) = routing_orog_max_glo(ij,ib)
routing_floodp(ij,ib) = routing_floodp_glo(ij,ib)
routing_cg(ij,ib,:) = routing_cg_glo(ij,ib,:)
topo_resid(ij,ib) = topo_resid_glo(ij,ib)
route_nbbasin(:) = route_count_glo(:)
route_togrid(ij,ib) = route_togrid_glo(ij,ib)
route_tobasin(ij,ib) = route_tobasin_glo(ij,ib)
routing_floodcri(ij,ib) = routing_floodcri_glo(ij,ib)
routing_fetch(ij,ib) = route_fetch_glo(ij,ib)
routing_beta(ij,ib) = routing_beta_glo(ij,ib)
!
route_nbintobas(ij) = route_nbintobas_glo(ij)
global_basinid(ij,ib) = global_basinid_glo(ij,ib)
route_outlet(ij,ib,:) = route_outlet_glo(ij,ib,:)
route_type(ij,ib) = route_type_glo(ij,ib)
origin_nbintobas(ij) = origin_nbintobas_glo(ij)
ENDDO
ENDDO
!
END SUBROUTINE finish_truncate
......
F90subroutines/routing_reg.f90
View file @ d6b18e5a
......@@ -12,6 +12,7 @@ MODULE routing_reg
INTEGER(i_std), SAVE :: nbpt_save
INTEGER(i_std), SAVE :: nbasmax_save
REAL(r_std), SAVE :: topoindexmin = 0.0
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: routing_area_glo !! Surface of basin (m^2)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: routing_orog_mean_glo !! Mean topography (m)
......@@ -153,6 +154,8 @@ CONTAINS
inc(8,1) = -1
inc(8,2) = -1
!
topoindexmin = min_topoind
!
! Set everything to undef to locate easily empty points
!
trip_bx(:,:) = undef_int
......@@ -201,6 +204,11 @@ CONTAINS
ll = MINLOC(ABS(latstr(1:nbj) - lat_rel(ib, ip)))
jloc = ll(1)
!
IF ( topoindex(ib, ip) < 0 .OR. topoindex(ib, ip) > 10000 ) THEN
WRITE(numout,*) "There is an issue with Topoindex : ", topoindex(ib, ip), min_topoind
WRITE(numout,*) "Check other variables : ", trip(ib, ip), basins(ib, ip), sub_area(ib, ip)
STOP
ENDIF
trip_bx(iloc, jloc) = NINT(trip(ib, ip))
basin_bx(iloc, jloc) = NINT(basins(ib, ip))
area_bx(iloc, jloc) = sub_area(ib, ip)
......@@ -2700,22 +2708,32 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
! For non existing basins the route_tobasin variable is put to zero. This will allow us
! to pick up the number of basin afterwards.
!
route_togrid_glo(ib,:) = ib
route_tobasin_glo(ib,:) = 0
routing_area_glo(ib,:) = zero
route_fetch_glo(ib,:) = zero
route_count_glo(ib) = basin_count(ib)
route_nbintobas_glo(ib) = basin_count(ib)
origin_nbintobas_glo(ib) = basin_notrun(ib)
routing_floodp_glo(ib,:) = 0
routing_orog_mean_glo(ib,:) = 0
routing_orog_min_glo(ib,:) = 0
routing_orog_max_glo(ib,:) = 0
routing_floodcri_glo(ib,:) = 0
!
routing_floodp_glo(ib,:) = 0
routing_beta_glo(ib,:) = 0
routing_cg_glo(ib,:,:) = 0.0
!
topo_resid_glo(ib,:) = topoindexmin
global_basinid_glo(ib,:) = zero
route_outlet_glo(ib,:,:) = zero
!
route_fetch_glo(ib,:) = zero
route_type_glo(ib,:) = zero
!
route_togrid_glo(ib,:) = ib
route_tobasin_glo(ib,:) = 0
routing_floodcri_glo(ib,:) = 0
!
route_count_glo(ib) = basin_count(ib)
route_nbintobas_glo(ib) = basin_count(ib)
origin_nbintobas_glo(ib) = basin_notrun(ib)
! If the land point doesn't overlap with any hydrogrid point
! (cf 0.5° trip on coast or some particular grid points )
IF ( basin_count(ib) .EQ. 0) THEN
IF ( basin_count(ib) .EQ. 1) THEN
basin_area(ib,1) = contfrac(ib)*gridarea(ib)
basin_orog_min(ib,1) = 0
basin_orog_max(ib,1) = 0
......@@ -2726,7 +2744,7 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
basin_cg(ib,1,2) = gridcenters(ib,2)
! To define in function of the minimum topoindex on the map
! (to revise)
basin_topoind(ib,1) = 1
basin_topoind(ib,1) = topoindexmin
basin_id(ib,1) = 0
basin_coor(ib,1,1) = basin_cg(ib,1,1)
basin_coor(ib,1,2) = basin_cg(ib,1,2)
......@@ -2753,20 +2771,19 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
routing_floodp_glo(ib,ij) = basin_floodp(ib,ij)/basin_area(ib,ij)
routing_beta_glo(ib,ij) = basin_beta_fp(ib,ij)
routing_cg_glo(ib,ij,:) = basin_cg(ib,ij,:)
!
topo_resid_glo(ib,ij) = basin_topoind(ib,ij)
global_basinid_glo(ib,ij) = basin_id(ib,ij)
route_outlet_glo(ib,ij,1) = basin_coor(ib,ij,1)
route_outlet_glo(ib,ij,2) = basin_coor(ib,ij,2)
!
route_fetch_glo(ib,ij) = fetch_basin(ib,ij)
!
route_type_glo(ib,ij) = basin_type(ib,ij)
!
route_togrid_glo(ib,ij) = outflow_grid(ib,ij)
route_tobasin_glo(ib,ij) = outflow_basin(ib,ij)
routing_floodcri_glo(ib,ij) = floodcri(ib,ij)*2000
!
!
ENDDO
ENDDO
!
......
GraphHydro.py
View file @ d6b18e5a
......@@ -24,6 +24,8 @@ import routing_interface
import NcHelp as NH
import RPPtools as RPP
#
NCFillValue=1.0e20
#
# Add time constants and possibly other prameters to the graph file
#
def addparameters(outnf, version, part, tcst, vtyp, NCFillValue) :
......@@ -84,10 +86,8 @@ class HydroGraph :
basin_flowdir = hydrosuper.basin_flowdir, \
outflow_grid = hydrosuper.outflow_grid, outflow_basin = hydrosuper.outflow_basin, \
inflow_number = hydrosuper.inflow_number, inflow_grid = hydrosuper.inflow_grid, \
inflow_basin = hydrosuper.inflow_basin, floodcri = hydrosuper.floodcri)
#
#
inflow_basin = hydrosuper.inflow_basin, floodcri = hydrosuper.floodcri, \
rfillval = NCFillValue, ifillval = RPP.IntFillValue)
#
self.num_largest = routing_interface.finalrivclass(part.landcorelist, self.route_togrid, self.route_tobasin, self.routing_fetch, \
hydrosuper.largest_rivarea)
......@@ -99,7 +99,7 @@ class HydroGraph :
routing_interface.finish_inflows(nbpt = self.nbpt, nwbas = nwbas, nbasmax = nbasmax, inf_max = self.max_inflow, \
basin_count = hydrosuper.basin_count, inflow_number = hydrosuper.inflow_number, \
inflow_grid = gingrid, inflow_basin = hydrosuper.inflow_basin[:,:,:self.max_inflow])
#
return
#
#
......@@ -165,7 +165,6 @@ class HydroGraph :
#
def dumpnetcdf(self, filename, version, globalgrid, procgrid, part, tcst, locations) :
#
NCFillValue=1.0e20
vtyp=np.float64
cornerind=[0,2,4,6]
nbcorners = len(cornerind)
......
Interface.py
View file @ d6b18e5a
......@@ -212,7 +212,7 @@ class HydroSuper :
nb_basin = nb_basin, basin_inbxid = basin_inbxid, basin_outlet = basin_outlet, basin_outtp = basin_outtp, \
basin_sz = self.basin_sz, basin_pts = self.basin_pts, basin_bxout = basin_bxout, \
basin_bbout = basin_bbout, lshead = basin_lshead, coast_pts = coast_pts, nwbas = self.nwbas)
#
# Memory management
del basin_bbout; del basin_lshead; del coast_pts; del basin_bxout; del self.basin_pts;
del basin_outtp; del basin_outlet; del basin_inbxid; del nb_basin
......
RoutingPreProc.py
View file @ d6b18e5a
......@@ -132,7 +132,7 @@ TR(hsuper, part, comm, modelgrid, numop = numop)
print("Rank : {0} - Basin_count After Truncate : ".format(part.rank), hsuper.basin_count)
hsuper.get_floodcri()
INFO("=================== HydroGraph ====================")
INFO("=================== GraphHydro ====================")
hgraph = GR.HydroGraph(nbasmax, hsuper, part, modelgrid)
INFO("=================== Locate gauging stations and other structures ====================")
......
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