diff --git a/F90subroutines/routing_interface.f90 b/F90subroutines/routing_interface.f90
index d9cef5eff8b566d4b1bf63f2b487963d22968074..c55e738d640cbdd32a936e8449aeb04f3fbe9739 100644
--- a/F90subroutines/routing_interface.f90
+++ b/F90subroutines/routing_interface.f90
@@ -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)
diff --git a/Interface.py b/Interface.py
index 827855a7d40bd6e96a149f225c8d0afe272208cc..66a9767b63a1c660f04524fd759b1b8d56d68788 100644
--- a/Interface.py
+++ b/Interface.py
@@ -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
diff --git a/RoutingPreProc.py b/RoutingPreProc.py
index fa755077aafbe436c152e3a5271ff8a860809c3e..d1876ce0dd36f478088e08aded6c1b8d2ac997b2 100644
--- a/RoutingPreProc.py
+++ b/RoutingPreProc.py
@@ -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)
diff --git a/tests/Iberia_n48/run.def b/tests/Iberia_n48/run.def
index e252d30a3064c468aa45316579c5e5125ee09113..6cafae4b1e0d80aadcbe108c9bf3fbd953e31be9 100644
--- a/tests/Iberia_n48/run.def
+++ b/tests/Iberia_n48/run.def
@@ -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
#
diff --git a/tests/Iberia_n48_regular/run.def b/tests/Iberia_n48_regular/run.def
index 7011ae3e9322099c4f0664bcdfcb40c22afb0933..b0afac7840d2d0c3dab4bd0c8e24ad757175fb80 100644
--- a/tests/Iberia_n48_regular/run.def
+++ b/tests/Iberia_n48_regular/run.def
@@ -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
#
diff --git a/tests/Mallorca/run.def b/tests/Mallorca/run.def
index 7ddf7f6698cc16802ed6681db12298a6599b5f8d..981458b2511f0ee507e5cd69801ae568cb179186 100644
--- a/tests/Mallorca/run.def
+++ b/tests/Mallorca/run.def
@@ -17,6 +17,10 @@ Documentation = true
#
nbasmax = 35
#
+# Number of operation of simplification performed together
+#
+numop = 200
+#
# Output
#
GraphFile = MEDCORDEX_test_graph.nc
diff --git a/tests/Pantanal_n48_regular/BuildHTUs_Pantanal.pbs b/tests/Pantanal_n48_regular/BuildHTUs_Pantanal.pbs
new file mode 100644
index 0000000000000000000000000000000000000000..57016f0b37acc5bd3190f7a3c363c0f2cec679ac
--- /dev/null
+++ b/tests/Pantanal_n48_regular/BuildHTUs_Pantanal.pbs
@@ -0,0 +1,35 @@
+#!/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
diff --git a/tests/Pantanal_n48_regular/run.def b/tests/Pantanal_n48_regular/run.def
new file mode 100644
index 0000000000000000000000000000000000000000..18671c0534e42928b932fa9e534cbb28bd52609f
--- /dev/null
+++ b/tests/Pantanal_n48_regular/run.def
@@ -0,0 +1,32 @@
+[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