From eedd88162120cdace73c029c6ed2fd57b393ba5e Mon Sep 17 00:00:00 2001
From: Anthony <anthony.schrapffer@polytechnique.fr>
Date: Wed, 15 Apr 2020 13:56:04 +0200
Subject: [PATCH] Add the inflows information in dumpnetcdf and simplify the
dumpnetcdf by using functions
---
F90subroutines/routing_interface.f90 | 48 +++++++
Interface.py | 191 +++++++++++----------------
2 files changed, 128 insertions(+), 111 deletions(-)
diff --git a/F90subroutines/routing_interface.f90 b/F90subroutines/routing_interface.f90
index d9cef5e..c55e738 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 827855a..66a9767 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
--
GitLab