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

Merge gitlab.in2p3.fr:jan.polcher/routingpp

Conflicts:
	RoutingPreProc.py
parents 3acee208 8216ff4e
......@@ -496,6 +496,54 @@ SUBROUTINE finish_truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridare
END SUBROUTINE finish_truncate
SUBROUTINE finish_inflows(nbpt,nbxmax_in, nbasmax, inf_max, basin_count, inflow_number, inflow_grid, inflow_basin,&
& route_innum, route_ingrid, route_inbasin)
!
USE ioipsl
USE grid
USE routing_tools
USE routing_reg
!
! Arguments
!
INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
INTEGER(i_std), INTENT (in) :: nbxmax_in, nbasmax
INTEGER(i_std), INTENT (in) :: inf_max !!
!
INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_count !!
INTEGER(i_std), DIMENSION(nbpt,nbxmax_in), INTENT(in) :: inflow_number !!
INTEGER(i_std), DIMENSION(nbpt,nbxmax_in,inf_max), INTENT(in) :: inflow_basin !!
INTEGER(i_std), DIMENSION(nbpt,nbxmax_in,inf_max), INTENT(in) :: inflow_grid !!
!
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: route_innum
REAL(r_std), DIMENSION(nbpt,nbasmax,inf_max), INTENT(out) :: route_ingrid
REAL(r_std), DIMENSION(nbpt,nbasmax,inf_max), INTENT(out) :: route_inbasin
!
!
!
! inflow_number -> route_innum
route_innum(:,:) = 0
DO ig=1,nbpt
nbas = basin_count(ig)
route_innum(ig,:nbas) = inflow_number(ig,:nbas)
END DO
! inflow_grid -> route_ingrid
! inflow_basin -> route_inbasin
route_ingrid(:,:,:) = 0
route_inbasin(:,:,:) = 0
DO ig=1,nbpt
nbas = basin_count(ig)
DO ib=1,nbas
nin = route_innum(ig,ib)
route_ingrid(ig,ib,:nin) = inflow_grid(ig,ib,:nin)
route_inbasin(ig,ib,:nin) = inflow_basin(ig,ib,:nin)
END DO
END DO
END SUBROUTINE finish_inflows
SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numops, basin_count, basin_area, &
& basin_cg, basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, &
& inflow_number, inflow_grid, inflow_basin)
......
......@@ -529,7 +529,35 @@ class HydroGraph :
self.num_largest = routing_interface.finalrivclass(part.landcorelist, self.route_togrid, self.route_tobasin, self.routing_fetch, \
hydrosuper.largest_rivarea)
#
# Inflows
self.max_inflow = part.domainmax(np.max(hydrosuper.inflow_number))
gingrid = part.l2glandindex( hydrosuper.inflow_grid[:,:,:self.max_inflow])
self.route_innum, self.route_ingrid, self.route_inbasin = routing_interface.finish_inflows(nbpt = self.nbpt, nbxmax_in = nbxmax_in, 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
#
#
#
def add_variable(self,outnf, procgrid, NCFillValue, part, coord, name, title, units, data, vtyp, orig_type = "float"):
var = procgrid.landscatter(data.astype(vtyp), order='F')
if orig_type == "float":
var[np.isnan(var)] = NCFillValue
elif orig_type == "int":
var[var>=RPP.IntFillValue] = NCFillValue
if part.rank == 0:
ncvar = outnf.createVariable(name, vtyp, coord, fill_value=NCFillValue)
ncvar.title = title
ncvar.units = units
else :
ncvar = np.zeros((1,1,1))
ncvar[:] = part.gather(var)
#
#
#
def dumpnetcdf(self, filename, globalgrid, procgrid, part) :
#
......@@ -546,14 +574,23 @@ class HydroGraph :
outnf.createDimension('land', len(procgrid.area))
outnf.createDimension('htu', self.nbasmax)
outnf.createDimension('bnd', nbcorners)
outnf.createDimension('inflow', self.max_inflow)
else :
outnf = None
#
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.
# land grid index -> to facilitate the analyses of the routing
#
nbpt_loc = np.zeros((self.nbpt,1)).astype(np.int32)
nbpt_loc[:,0] = np.arange(1, self.nbpt+1)
nbpt_glo = part.l2glandindex(nbpt_loc)
self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "nbpt_glo", "Grid point Global", "-", nbpt_glo[:,0], vtyp)
#
################
#
# TEST: l2glandindex
itarget=-1
for il in range(procgrid.nbland) :
lo = procgrid.lon_full[procgrid.indP[il][0],procgrid.indP[il][1]]
......@@ -564,130 +601,62 @@ class HydroGraph :
if itarget >+ 0 :
print(part.rank, itarget, " Before route_togrid = ", self.route_togrid[itarget,:])
# Conversion
grgrid = part.l2glandindex(self.route_togrid[:,:])
if itarget >+ 0 :
print(part.rank, itarget, " After route_togrid = ", self.route_togrid[itarget,:])
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)
routetogrid.title = "Grid into which the basin flows"
routetogrid.units = "-"
else :
routetogrid = np.zeros((1,1,1))
routetogrid[:,:,:] = part.gather(rgrid)
#
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 :
routetobasin = outnf.createVariable("routetobasin", vtyp, ('htu','y','x'), fill_value=NCFillValue)
routetobasin.title = "Basin in to which the water goes"
routetobasin.units = "-"
else :
routetobasin = np.zeros((1,1,1))
routetobasin[:,:,:] = part.gather(rtobasin)
#
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)
basinid.title = "ID of basin"
basinid.units = "-"
else :
basinid = np.zeros((1,1,1))
basinid[:,:,:] = part.gather(rid)
#
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)
routenbintobas.title = "Number of basin into current one"
routenbintobas.units = "-"
else :
routenbintobas = np.zeros((1,1))
routenbintobas[:,:] = part.gather(rintobas)
#
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)
originnbintobas.title = "Number of sub-grid basin into current one before truncation"
originnbintobas.units = "-"
else :
originnbintobas = np.zeros((1,1))
originnbintobas[:,:] = part.gather(onbintobas)
################
#
# The field route_togrid is with indices on the local grid. That needs to be converted to the global grid.
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "routetogrid", "Grid into which the basin flows", "-", grgrid, vtyp, "int")
# route to basin
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "routetobasin", "Basin in to which the water goes", "-", self.route_tobasin[:,:], vtyp, "int")
# basin id
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "basinid", "ID of basin", "-", self.global_basinid[:,:], vtyp, "int")
#
# basin area
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "basin_area", "area of basin", "m^2", self.routing_area[:,:], vtyp, "float")
# route number into basin
self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "routenbintobas", "Number of basin into current one", "-", self.route_nbintobas[:], vtyp, "int")
#
# original number into basin
self.add_variable(outnf, procgrid, NCFillValue, part, ( 'y','x'), "originnbintobas", "Number of sub-grid basin into current one before truncation", "-", self.origin_nbintobas[:], vtyp, "int")
#
# latitude of outlet
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "outletlat", "Latitude of Outlet", "degrees north", self.route_outlet[:,:,0], vtyp, "float")
# longitude of outlet
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "outletlon", "Longitude of Outlet", "degrees east", self.route_outlet[:,:,1], vtyp, "float")
# type of outlet
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "outlettype", "Type of outlet", "code", self.route_type[:,:], vtyp, "float")
#
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)
outletlat.title = "Latitude of Outlet"
outletlat.title = "degrees north"
else :
outletlat = np.zeros((1,1,1))
outletlat[:,:,:] = part.gather(olat)
# topographic index
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "topoindex", "Topographic index of the retention time", "m", self.topo_resid[:,:], vtyp, "float")
#
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)
outletlon.title = "Longitude of outlet"
outletlon.units = "degrees east"
else :
outletlon = np.zeros((1,1,1))
outletlon[:,:,:] = part.gather(olon)
# Inflow number
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "route_innum", "Number of inflow", "-", self.route_innum[:,:], vtyp, "int")
#
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)
outlettype.title = "Type of outlet"
outlettype.units = "code"
else :
outlettype = np.zeros((1,1,1))
outlettype[:,:,:] = part.gather(otype)
# Inflow grid
#gingrid = part.l2glandindex(self.inflow_grid[:,:,:inflow_size])
self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow', 'htu','y','x'), "route_ingrid", "Grid from which the water flows", "-", self.route_ingrid[:,:,:], vtyp, "int")
#
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)
topoindex.title = "Topographic index of the retention time"
topoindex.units = "m"
else :
topoindex = np.zeros((1,1,1))
topoindex[:,:,:] = part.gather(tind)
# Inflow basin
self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow', 'htu','y','x'), "route_inbasin", "Basin from which the water flows", "-", self.route_inbasin[:,:,:], vtyp, "int")
#
# Save centre of gravity of HTU
#
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)
CG_lon.title = "Longitude of centre of gravity of HTU"
CG_lon.units = "degrees east"
CG_lat = outnf.createVariable("CG_lat", vtyp, ('htu','y','x'), fill_value=NCFillValue)
CG_lat.title = "Latitude of centre of gravity of HTU"
CG_lat.units = "degrees north"
else :
CG_lon = np.zeros((1,1,1))
CG_lat = np.zeros((1,1,1))
CG_lon[:,:,:] = part.gather(cg[1,:,:,:])
CG_lat[:,:,:] = part.gather(cg[0,:,:,:])
# Check if it works
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "CG_lon", "Longitude of centre of gravity of HTU", "degrees east", self.routing_cg[:,:,1], vtyp, "float")
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "CG_lat", "Latitude of centre of gravity of HTU", "degrees north", self.routing_cg[:,:,0], vtyp, "float")
#
# Save the fetch of each basin
#
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)
fetch.title = "Fetch contributing to each HTU"
fetch.units = "m^2"
else :
fetch = np.zeros((1,1,1))
fetch[:,:,:] = part.gather(fe)
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "fetch", "Fetch contributing to each HTU", "m^2", self.routing_fetch[:,:], vtyp, "float")
#
if part.rank == 0 :
# Close the file
if part.rank == 0:
outnf.close()
#
return
......
......@@ -15,9 +15,10 @@ import time
# Gert the information from the configuration file.
#
import configparser
config=configparser.ConfigParser({"DiagLon":"0.0, 0.0", "DiagLat":"0.0, 0.0"})
config=configparser.ConfigParser({"DiagLon":"0.0, 0.0", "DiagLat":"0.0, 0.0", "numop" : 100})
config.read("run.def")
nbasmax=int(config.getfloat("OverAll", "nbasmax"))
numop=int(config.getfloat("OverAll", "numop"))
OutGraphFile=config.get("OverAll","GraphFile")
lonint=np.array(config.get("OverAll", "DiagLon").split(","),dtype=float)
latint=np.array(config.get("OverAll", "DiagLat").split(","),dtype=float)
......@@ -140,7 +141,7 @@ comm.Barrier()
hsuper.dumpnetcdf(OutGraphFile.replace(".nc","_HydroSuper.nc"), gg, modelgrid, part)
print("Rank : {0} - Basin_count Before Truncate : ".format(part.rank), hsuper.basin_count)
hs = TR(hsuper, part, comm, modelgrid, numop = 200)
hs = TR(hsuper, part, comm, modelgrid, numop = numop)
print("Rank : {0} - Basin_count After Truncate : ".format(part.rank), hsuper.basin_count)
hgraph = IF.HydroGraph(nbasmax, hsuper, part, modelgrid)
......
......@@ -15,7 +15,11 @@ Documentation = true
#
# Configuration for the graph to be generated
#
nbasmax = 100
nbasmax = 80
#
# Number of operation of simplification performed together
#
numop = 200
#
# Output
#
......
......@@ -15,7 +15,11 @@ Documentation = true
#
# Configuration for the graph to be generated
#
nbasmax = 50
nbasmax = 35
#
# Number of operation of simplification performed together
#
numop = 100
#
# Output
#
......
......@@ -17,6 +17,10 @@ Documentation = true
#
nbasmax = 35
#
# Number of operation of simplification performed together
#
numop = 200
#
# Output
#
GraphFile = MEDCORDEX_test_graph.nc
......
#!/bin/bash
#
#PBS -N BuildHTUs_P
#
#PBS -j oe
#PBS -l nodes=1:ppn=48
#PBS -l walltime=128:00:00
#PBS -l mem=252gb
#PBS -l vmem=252gb
#
cd ${PBS_O_WORKDIR}
export NSLOTS=$(($PBS_NUM_NODES*$PBS_NUM_PPN))
#
# Set the right Python 3 Anaconda environment
#
source ../../Environment
#
# Clean-up. Weights are kept for future runs.
#
/bin/rm -f DocumentationInterface *.nc *.txt
#
# Run the Python code to generate the HTUs and write them into a netCDF file.
#
mpirun -n ${NSLOTS} python ../../RoutingPreProc.py
if [ $? -gt 0 ] ; then
echo "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
echo "X Run on Iberian Peninsula failed X"
echo "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
else
echo "======================================="
echo "= Run on Iberian Peninsula successful ="
echo "======================================="
fi
#
ls -l
[OverAll]
#
#
#
EarthRadius = 6370000.
#
ModelGridFile = /bdd/ORCHIDEE_Forcing/BC/OOL/OL2/WFDEI_CRU/WFDEI_CRU_1994.nc
WEST_EAST = -68.75, -53.25
SOUTH_NORTH = -24.5, -14.25
HydroFile = /homedata/aschrapffer/RoutingNc_AS/SouthAmerica_FP/routing.nc
#
# FORTRAN interface parameters
#
Documentation = true
#
# Configuration for the graph to be generated
#
nbasmax = 50
#
# Number of operation of simplification performed together
#
numop = 200
#
# Output
#
GraphFile = WFDEI_CRU_Pantanal_test_graph.nc
#
# Diagnostics
# You need to provide an interval in longitude and Latitude.
#
DiagLon = 3.9, 3.9
DiagLat = 40.0, 40.0
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