Skip to content
Snippets Groups Projects
Commit e8810ece authored by POLCHER Jan's avatar POLCHER Jan :bicyclist_tone4:
Browse files

Added the function to dump the graph of HydroSuper, that is before the truncation.

parent b299412d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
F90subroutines/routing_interface.f90
+ 2
− 1
View file @ e8810ece
......@@ -154,7 +154,7 @@ SUBROUTINE findbasins(nbpt, nbvmax_in, nbxmax_in, nbi, nbj, trip_bx, basin_bx, f
!
INTEGER :: ib
!
diaglalo(1,:) = (/ 38.25, -7.25 /)
diaglalo(1,:) = (/ 39.6791, 2.6687 /)
!
IF ( nbvmax_in .NE. nbvmax .OR. nbxmax_in .NE. nbxmax ) THEN
WRITE(*,*) "nbvmax or nbxmax have changed !!"
......@@ -508,6 +508,7 @@ SUBROUTINE finalfetch(nbpt, nbasmax, nbcore, corepts, routing_area, basin_count,
!
USE defprec
USE ioipsl
USE constantes_var
!
!! INPUT VARIABLES
INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
......
Interface.py
+ 222
− 125
View file @ e8810ece
......@@ -71,6 +71,106 @@ def closeinterface(comm) :
#
#
#
def addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind) :
#
# Longitude
longitude = part.gather(procgrid.lon_full)
if part.rank == 0 :
lon=outnf.createVariable("lon", vtyp, ('y','x'), fill_value=NCFillValue)
lon.units="grid box centre degrees east"
lon.title="Longitude"
lon.axis="X"
lon[:,:] = longitude[:,:]
#
# Latitude
latitude = part.gather(procgrid.lat_full)
if part.rank == 0 :
lat=outnf.createVariable("lat", vtyp, ('y','x'), fill_value=NCFillValue)
lat.units="grid box centre degrees north"
lat.standard_name="grid latitude"
lat.title="Latitude"
lat.axis="Y"
lat[:] = latitude[:,:]
#
# Bounds of grid box
#
llonpoly=np.zeros((nbcorners,procgrid.nbland))
llatpoly=np.zeros((nbcorners,procgrid.nbland))
for i in range(procgrid.nbland) :
llonpoly[:,i] = [procgrid.polyll[i][ic][0] for ic in cornerind]
llatpoly[:,i] = [procgrid.polyll[i][ic][1] for ic in cornerind]
lon_bnd = procgrid.landscatter(llonpoly)
lat_bnd = procgrid.landscatter(llatpoly)
if part.rank == 0 :
lonbnd=outnf.createVariable("lon_bnd", vtyp, ('bnd','y','x'), fill_value=NCFillValue)
lonbnd.units="grid box corners degrees east"
lonbnd.title="Longitude of Corners"
latbnd=outnf.createVariable("lat_bnd", vtyp, ('bnd','y','x'), fill_value=NCFillValue)
latbnd.units="grid box corners degrees north"
latbnd.title="Latitude of Corners"
else :
lonbnd= np.zeros((1,1,1))
latbnd= np.zeros((1,1,1))
lonbnd[:,:,:] = part.gather(lon_bnd[:,:,:])
latbnd[:,:,:] = part.gather(lat_bnd[:,:,:])
#
# Land sea mask
#
if part.rank == 0 :
land=outnf.createVariable("land", vtyp, ('y','x'), fill_value=NCFillValue)
land.units="Land Sea mask"
land.standard_name="landsea mask"
land.title="Land"
land[:,:] = globalgrid.land[:,:]
# Area
areas = procgrid.landscatter(np.array(procgrid.area, dtype=np.float64))
areas[np.isnan(areas)] = NCFillValue
if part.rank == 0 :
area=outnf.createVariable("area", vtyp, ('y','x'), fill_value=NCFillValue)
area.units="m^2"
area.standard_name="grid area"
area.title="Area"
else :
area = np.zeros((1,1))
area[:,:] = part.gather(areas[:,:])
#
return
#
# Add environment to netCDF file
#
def addenvironment(outnf, procgrid, part, vtyp, NCFillValue, nbpt) :
#
nbpt_proc = np.arange(1,nbpt+1, dtype=vtyp)
proc = np.full(nbpt, part.rank, dtype=vtyp)
# Environment
# nbpt_proc
subpt = procgrid.landscatter(nbpt_proc[:], order='F')
subpt = subpt.astype(vtyp, copy=False)
subpt[np.isnan(subpt)] = NCFillValue
if part.rank == 0 :
subptgrid = outnf.createVariable("nbpt_proc", vtyp, ('y','x'), fill_value=NCFillValue)
subptgrid.title = "gridpoint reference inside each proc"
subptgrid.units = "-"
else :
subptgrid = np.zeros((1,1))
subptgrid[:,:] = part.gather(subpt)
#
# rank
procnum = procgrid.landscatter(proc[:], order='F')
procnum = procnum.astype(vtyp, copy=False)
procnum[np.isnan(procnum)] = NCFillValue
if part.rank == 0 :
procn = outnf.createVariable("proc_num", vtyp, ('y','x'), fill_value=NCFillValue)
procn.title = "rank"
procn.units = "-"
else :
procn = np.zeros((1,1))
procn[:,:] = part.gather(procnum)
#
return
#
#
#
def finalfetch(part, routing_area, basin_count, route_togrid, route_tobasin, fetch_in) :
#
fetch_out = np.zeros(routing_area.shape, dtype=np.float32, order='F')
......@@ -169,6 +269,10 @@ class HydroOverlap :
#
class HydroSuper :
def __init__(self, nbvmax, hydrodata, hydrooverlap) :
#
# Keep largest possible number of HTUs
#
self.nbhtuext = nbvmax
#
# Call findbasins
#
......@@ -191,6 +295,7 @@ class HydroSuper :
hydrooverlap.hierarchy_bx, hydrooverlap.fac_bx, hydrooverlap.topoind_bx, hydrodata.topoindmin, \
nb_basin, basin_inbxid, basin_outlet, basin_outtp, self.basin_sz, self.basin_pts, basin_bxout, \
basin_bbout, basin_lshead, coast_pts, hydrooverlap.nwbas)
self.nbpt = self.basin_count.shape[0]
return
#
def linkup(self, hydrodata) :
......@@ -249,12 +354,115 @@ class HydroSuper :
self.fetch_basin, self.largest_rivarea)
print("Rank :"+str(part.rank)+" Area of smallest large rivers : ", self.largest_rivarea, " Nb of Large rivers on proc : ",self.num_largest)
return
#
#
#
def dumpnetcdf(self, filename, globalgrid, procgrid, part) :
#
NCFillValue=1.0e20
vtyp=np.float64
cornerind=[0,2,4,6]
nbcorners = len(cornerind)
#
if part.rank == 0 :
outnf=Dataset(filename, 'w', format='NETCDF4_CLASSIC')
# Dimensions
outnf.createDimension('x', globalgrid.ni)
outnf.createDimension('y', globalgrid.nj)
outnf.createDimension('land', len(procgrid.area))
outnf.createDimension('htuext', self.nbhtuext)
outnf.createDimension('bnd', nbcorners)
else :
outnf = None
#
addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind)
addenvironment(outnf, procgrid, part, vtyp, NCFillValue, self.nbpt)
#
# Variables
# Once gathered on root_proc we transform them into float64. Careful, Integer variables do not have NaN !
#
# self.basin_id
bid = procgrid.landscatter(self.basin_id[:,:].astype(vtyp), order='F')
bid[np.isnan(bid)] = NCFillValue
if part.rank == 0 :
basinid = outnf.createVariable("HTUid", vtyp, ('htuext','y','x'), fill_value=NCFillValue)
basinid.title = "ID for each HTU"
basinid.units = "-"
else :
basinid = np.zeros((1,1,1))
basinid[:,:,:] = part.gather(bid)
#
#self.basin_count
bcnt = procgrid.landscatter(self.basin_count[:].astype(vtyp), order='F')
bcnt[np.isnan(bcnt)] = NCFillValue
if part.rank == 0 :
basincnt = outnf.createVariable("HTUcount", vtyp, ('y','x'), fill_value=NCFillValue)
basincnt.title = "HTU count"
basincnt.units = "-"
else :
basincnt = np.zeros((1,1))
basincnt[:,:] = part.gather(bcnt)
#
#self.basin_area
ba = procgrid.landscatter(self.basin_area[:,:].astype(vtyp), order='F')
ba[np.isnan(ba)] = NCFillValue
if part.rank == 0 :
basinarea = outnf.createVariable("HTUarea", vtyp, ('htuext','y','x'), fill_value=NCFillValue)
basinarea.title = "HTU area"
basinarea.units = "m^2"
else :
basinarea = np.zeros((1,1,1))
basinarea[:,:,:] = part.gather(ba)
#
#self.outflow_grid
og = procgrid.landscatter(self.outflow_grid[:,:].astype(vtyp), order='F')
og[np.isnan(og)] = NCFillValue
if part.rank == 0 :
outgrid = outnf.createVariable("HTUoutgrid", vtyp, ('htuext','y','x'), fill_value=NCFillValue)
outgrid.title = "HTU outflow grid"
outgrid.units = "-"
else :
outgrid = np.zeros((1,1,1))
outgrid[:,:,:] = part.gather(og)
#
#self.outflow_basin
ob = procgrid.landscatter(self.outflow_basin[:,:].astype(vtyp), order='F')
ob[np.isnan(ob)] = NCFillValue
if part.rank == 0 :
outbas = outnf.createVariable("HTUoutbasin", vtyp, ('htuext','y','x'), fill_value=NCFillValue)
outbas.title = "Outflow HTU of grid"
outbas.units = "-"
else :
outbas = np.zeros((1,1,1))
outbas[:,:,:] = part.gather(ob)
#
# Save the fetch of each basin
#
fe = procgrid.landscatter(self.fetch_basin[:,:].astype(vtyp), order='F')
fe[np.isnan(fe)] = NCFillValue
if part.rank == 0 :
fetch = outnf.createVariable("fetch", vtyp, ('htuext','y','x'), fill_value=NCFillValue)
fetch.title = "Fetch contributing to each HTU"
fetch.units = "m^2"
else :
fetch = np.zeros((1,1,1))
fetch[:,:,:] = part.gather(fe)
#
# Close file
#
if part.rank == 0 :
outnf.close()
#
return
#
#
#
class HydroGraph :
def __init__(self, nbasmax, hydrosuper, part) :
#
self.nbasmax = nbasmax
self.nbpt = hydrosuper.basin_count.shape[0]
#
self.routing_area, self.routing_cg, self.topo_resid, self.route_nbbasin, self.route_togrid, self.route_tobasin, self.route_nbintobas, \
self.global_basinid, self.route_outlet, self.route_type, self.origin_nbintobas, self.routing_fetch = \
routing_interface.truncate(nbasmax, hydrosuper.num_largest, part.landcorelist, hydrosuper.basin_count, \
......@@ -269,9 +477,6 @@ class HydroGraph :
self.num_largest = routing_interface.finalrivclass(part.landcorelist, self.route_togrid, self.route_tobasin, self.routing_fetch, \
hydrosuper.largest_rivarea)
#
nbpt = hydrosuper.basin_count.shape[0]
self.nbpt_proc = np.arange(1,nbpt+1)
self.proc = np.full(nbpt, part.rank)
return
#
def dumpnetcdf(self, filename, globalgrid, procgrid, part) :
......@@ -289,109 +494,11 @@ class HydroGraph :
outnf.createDimension('land', len(procgrid.area))
outnf.createDimension('htu', self.nbasmax)
outnf.createDimension('bnd', nbcorners)
#
# Coordinates
#
# Longitude
longitude = part.gather(procgrid.lon_full)
if part.rank == 0 :
lon=outnf.createVariable("lon", vtyp, ('y','x'), fill_value=NCFillValue)
lon.units="grid box centre degrees east"
lon.title="Longitude"
lon.axis="X"
lon[:,:] = longitude[:,:]
#
# Latitude
latitude = part.gather(procgrid.lat_full)
if part.rank == 0 :
lat=outnf.createVariable("lat", vtyp, ('y','x'), fill_value=NCFillValue)
lat.units="grid box centre degrees north"
lat.standard_name="grid latitude"
lat.title="Latitude"
lat.axis="Y"
lat[:] = latitude[:,:]
#
# Bounds of grid box
#
llonpoly=np.zeros((nbcorners,procgrid.nbland))
llatpoly=np.zeros((nbcorners,procgrid.nbland))
for i in range(procgrid.nbland) :
llonpoly[:,i] = [procgrid.polyll[i][ic][0] for ic in cornerind]
llatpoly[:,i] = [procgrid.polyll[i][ic][1] for ic in cornerind]
lon_bnd = procgrid.landscatter(llonpoly)
lat_bnd = procgrid.landscatter(llatpoly)
if part.rank == 0 :
lonbnd=outnf.createVariable("lon_bnd", vtyp, ('bnd','y','x'), fill_value=NCFillValue)
lonbnd.units="grid box corners degrees east"
lonbnd.title="Longitude of Corners"
latbnd=outnf.createVariable("lat_bnd", vtyp, ('bnd','y','x'), fill_value=NCFillValue)
latbnd.units="grid box corners degrees north"
latbnd.title="Latitude of Corners"
else :
lonbnd= np.zeros((1,1,1))
latbnd= np.zeros((1,1,1))
lonbnd[:,:,:] = part.gather(lon_bnd[:,:,:])
latbnd[:,:,:] = part.gather(lat_bnd[:,:,:])
#
# Land sea mask
#
if part.rank == 0 :
land=outnf.createVariable("land", vtyp, ('y','x'), fill_value=NCFillValue)
land.units="Land Sea mask"
land.standard_name="landsea mask"
land.title="Land"
land[:,:] = globalgrid.land[:,:]
# Area
areas = procgrid.landscatter(np.array(procgrid.area, dtype=np.float64))
areas[np.isnan(areas)] = NCFillValue
if part.rank == 0 :
area=outnf.createVariable("area", vtyp, ('y','x'), fill_value=NCFillValue)
area.units="m^2"
area.standard_name="grid area"
area.title="Area"
else :
area = np.zeros((1,1))
area[:,:] = part.gather(areas[:,:])
#
# Environment
# nbpt_proc
subpt = procgrid.landscatter(self.nbpt_proc[:], order='F')
subpt = subpt.astype(vtyp, copy=False)
subpt[np.isnan(subpt)] = NCFillValue
if part.rank == 0 :
subptgrid = outnf.createVariable("nbpt_proc", vtyp, ('y','x'), fill_value=NCFillValue)
subptgrid.title = "gridpoint reference inside each proc"
subptgrid.units = "-"
else :
subptgrid = np.zeros((1,1,1))
subptgrid[:,:] = part.gather(subpt)
#
# rank
procnum = procgrid.landscatter(self.proc[:], order='F')
procnum = procnum.astype(vtyp, copy=False)
procnum[np.isnan(procnum)] = NCFillValue
if part.rank == 0 :
procn = outnf.createVariable("proc_num", vtyp, ('y','x'), fill_value=NCFillValue)
procn.title = "rank"
procn.units = "-"
else :
procn = np.zeros((1,1,1))
procn[:,:] = part.gather(procnum)
#
# Variables
# Once gathered on root_proc we transform them into float64. Careful, Integer variables do not have NaN !
outnf = None
#
rarea = procgrid.landscatter(self.routing_area[:,:], order='F')
rarea = rarea.astype(vtyp, copy=False)
rarea[np.isnan(rarea)] = NCFillValue
if part.rank == 0 :
routingarea = outnf.createVariable("routingarea", vtyp, ('htu','y','x'), fill_value=NCFillValue)
routingarea.title = "Surface of basin"
routingarea.units = "m^2"
else :
routingarea = np.zeros((1,1,1))
routingarea[:,:,:] = part.gather(rarea)
addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind)
addenvironment(outnf, procgrid, part, vtyp, NCFillValue, self.nbpt)
#
# The field route_togrid is with indices on the local grid. That needs to be converted to the global grid.
#
......@@ -408,8 +515,7 @@ class HydroGraph :
grgrid = part.l2glandindex(self.route_togrid[:,:])
if itarget >+ 0 :
print(part.rank, itarget, " After route_togrid = ", self.route_togrid[itarget,:])
rgrid = procgrid.landscatter(grgrid, order='F')
rgrid = rgrid.astype(vtyp, copy=False)
rgrid = procgrid.landscatter(grgrid.astype(vtyp), order='F')
rgrid[rgrid >= RPP.IntFillValue] = NCFillValue
if part.rank == 0 :
routetogrid = outnf.createVariable("routetogrid", vtyp, ('htu','y','x'), fill_value=NCFillValue)
......@@ -419,7 +525,7 @@ class HydroGraph :
routetogrid = np.zeros((1,1,1))
routetogrid[:,:,:] = part.gather(rgrid)
#
rtobasin = procgrid.landscatter(self.route_tobasin[:,:], order='F')
rtobasin = procgrid.landscatter(self.route_tobasin[:,:].astype(vtyp), order='F')
rtobasin = rtobasin.astype(vtyp, copy=False)
rtobasin[rtobasin >= RPP.IntFillValue] = NCFillValue
if part.rank == 0 :
......@@ -430,8 +536,7 @@ class HydroGraph :
routetobasin = np.zeros((1,1,1))
routetobasin[:,:,:] = part.gather(rtobasin)
#
rid = procgrid.landscatter(self.global_basinid[:,:], order='F')
rid = rid.astype(vtyp, copy=False)
rid = procgrid.landscatter(self.global_basinid[:,:].astype(vtyp), order='F')
rid[rid >= RPP.IntFillValue] = NCFillValue
if part.rank == 0 :
basinid = outnf.createVariable("basinid", vtyp, ('htu','y','x'), fill_value=NCFillValue)
......@@ -441,8 +546,7 @@ class HydroGraph :
basinid = np.zeros((1,1,1))
basinid[:,:,:] = part.gather(rid)
#
rintobas = procgrid.landscatter(self.route_nbintobas[:])
rintobas = rintobas.astype(vtyp, copy=False)
rintobas = procgrid.landscatter(self.route_nbintobas[:].astype(vtyp))
rintobas[rintobas >= RPP.IntFillValue] = NCFillValue
if part.rank == 0 :
routenbintobas = outnf.createVariable("routenbintobas", vtyp, ('y','x'), fill_value=NCFillValue)
......@@ -452,8 +556,7 @@ class HydroGraph :
routenbintobas = np.zeros((1,1))
routenbintobas[:,:] = part.gather(rintobas)
#
onbintobas = procgrid.landscatter(self.origin_nbintobas[:])
onbintobas = onbintobas.astype(vtyp, copy=False)
onbintobas = procgrid.landscatter(self.origin_nbintobas[:].astype(vtyp))
onbintobas[onbintobas >= RPP.IntFillValue] = NCFillValue
if part.rank == 0 :
originnbintobas = outnf.createVariable("originnbintobas", vtyp, ('y','x'), fill_value=NCFillValue)
......@@ -463,8 +566,7 @@ class HydroGraph :
originnbintobas = np.zeros((1,1))
originnbintobas[:,:] = part.gather(onbintobas)
#
olat = procgrid.landscatter(self.route_outlet[:,:,0], order='F')
olat = olat.astype(vtyp, copy=False)
olat = procgrid.landscatter(self.route_outlet[:,:,0].astype(vtyp), order='F')
olat[np.isnan(olat)] = NCFillValue
if part.rank == 0 :
outletlat = outnf.createVariable("outletlat", vtyp, ('htu','y','x'), fill_value=NCFillValue)
......@@ -474,8 +576,7 @@ class HydroGraph :
outletlat = np.zeros((1,1,1))
outletlat[:,:,:] = part.gather(olat)
#
olon = procgrid.landscatter(self.route_outlet[:,:,1], order='F')
olon = olon.astype(vtyp, copy=False)
olon = procgrid.landscatter(self.route_outlet[:,:,1].astype(vtyp), order='F')
olon[np.isnan(olon)] = NCFillValue
if part.rank == 0 :
outletlon = outnf.createVariable("outletlon", vtyp, ('htu','y','x'), fill_value=NCFillValue)
......@@ -485,8 +586,7 @@ class HydroGraph :
outletlon = np.zeros((1,1,1))
outletlon[:,:,:] = part.gather(olon)
#
otype = procgrid.landscatter(self.route_type[:,:], order='F')
otype = otype.astype(vtyp, copy=False)
otype = procgrid.landscatter(self.route_type[:,:].astype(vtyp), order='F')
otype[np.isnan(otype)] = NCFillValue
if part.rank == 0 :
outlettype = outnf.createVariable("outlettype", vtyp, ('htu','y','x'), fill_value=NCFillValue)
......@@ -496,8 +596,7 @@ class HydroGraph :
outlettype = np.zeros((1,1,1))
outlettype[:,:,:] = part.gather(otype)
#
tind = procgrid.landscatter(self.topo_resid[:,:], order='F')
tind = tind.astype(vtyp, copy=False)
tind = procgrid.landscatter(self.topo_resid[:,:].astype(vtyp), order='F')
tind[np.isnan(tind)] = NCFillValue
if part.rank == 0 :
topoindex = outnf.createVariable("topoindex", vtyp, ('htu','y','x'), fill_value=NCFillValue)
......@@ -509,8 +608,7 @@ class HydroGraph :
#
# Save centre of gravity of HTU
#
cg = procgrid.landscatter(self.routing_cg[:,:,:], order='F')
cg = cg.astype(vtyp, copy=False)
cg = procgrid.landscatter(self.routing_cg[:,:,:].astype(vtyp), order='F')
cg[np.isnan(cg)] = NCFillValue
if part.rank == 0 :
CG_lon = outnf.createVariable("CG_lon", vtyp, ('htu','y','x'), fill_value=NCFillValue)
......@@ -527,8 +625,7 @@ class HydroGraph :
#
# Save the fetch of each basin
#
fe = procgrid.landscatter(self.routing_fetch[:,:], order='F')
fe = fe.astype(vtyp, copy=False)
fe = procgrid.landscatter(self.routing_fetch[:,:].astype(vtyp), order='F')
fe[np.isnan(fe)] = NCFillValue
if part.rank == 0 :
fetch = outnf.createVariable("fetch", vtyp, ('htu','y','x'), fill_value=NCFillValue)
......
RoutingPreProc.py
+ 1
− 0
View file @ e8810ece
......@@ -175,6 +175,7 @@ gc.collect()
print("=================== Compute fetch ====================")
hsuper.fetch(part)
hsuper.dumpnetcdf(OutGraphFile.replace(".nc","_HydroSuper.nc"), gg, modelgrid, part)
hgraph = IF.HydroGraph(nbasmax, hsuper, part)
......
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