Commit eedd8816 authored by Anthony's avatar Anthony
Browse files

Add the inflows information in dumpnetcdf and simplify the dumpnetcdf by using functions

parent ae16979e
...@@ -496,6 +496,54 @@ SUBROUTINE finish_truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridare ...@@ -496,6 +496,54 @@ SUBROUTINE finish_truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridare
END SUBROUTINE finish_truncate 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, & 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, & & basin_cg, basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, &
& inflow_number, inflow_grid, inflow_basin) & inflow_number, inflow_grid, inflow_basin)
......
...@@ -529,7 +529,35 @@ class HydroGraph : ...@@ -529,7 +529,35 @@ class HydroGraph :
self.num_largest = routing_interface.finalrivclass(part.landcorelist, self.route_togrid, self.route_tobasin, self.routing_fetch, \ self.num_largest = routing_interface.finalrivclass(part.landcorelist, self.route_togrid, self.route_tobasin, self.routing_fetch, \
hydrosuper.largest_rivarea) 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 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) : def dumpnetcdf(self, filename, globalgrid, procgrid, part) :
# #
...@@ -546,14 +574,23 @@ class HydroGraph : ...@@ -546,14 +574,23 @@ class HydroGraph :
outnf.createDimension('land', len(procgrid.area)) outnf.createDimension('land', len(procgrid.area))
outnf.createDimension('htu', self.nbasmax) outnf.createDimension('htu', self.nbasmax)
outnf.createDimension('bnd', nbcorners) outnf.createDimension('bnd', nbcorners)
outnf.createDimension('inflow', self.max_inflow)
else : else :
outnf = None outnf = None
# #
addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind) addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind)
addenvironment(outnf, procgrid, part, vtyp, NCFillValue, self.nbpt) 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 itarget=-1
for il in range(procgrid.nbland) : for il in range(procgrid.nbland) :
lo = procgrid.lon_full[procgrid.indP[il][0],procgrid.indP[il][1]] lo = procgrid.lon_full[procgrid.indP[il][0],procgrid.indP[il][1]]
...@@ -564,130 +601,62 @@ class HydroGraph : ...@@ -564,130 +601,62 @@ class HydroGraph :
if itarget >+ 0 : if itarget >+ 0 :
print(part.rank, itarget, " Before route_togrid = ", self.route_togrid[itarget,:]) print(part.rank, itarget, " Before route_togrid = ", self.route_togrid[itarget,:])
# Conversion
grgrid = part.l2glandindex(self.route_togrid[:,:]) grgrid = part.l2glandindex(self.route_togrid[:,:])
if itarget >+ 0 : if itarget >+ 0 :
print(part.rank, itarget, " After route_togrid = ", self.route_togrid[itarget,:]) 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 : # The field route_togrid is with indices on the local grid. That needs to be converted to the global grid.
routetogrid = outnf.createVariable("routetogrid", vtyp, ('htu','y','x'), fill_value=NCFillValue) self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "routetogrid", "Grid into which the basin flows", "-", grgrid, vtyp, "int")
routetogrid.title = "Grid into which the basin flows" # route to basin
routetogrid.units = "-" self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "routetobasin", "Basin in to which the water goes", "-", self.route_tobasin[:,:], vtyp, "int")
else : # basin id
routetogrid = np.zeros((1,1,1)) self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "basinid", "ID of basin", "-", self.global_basinid[:,:], vtyp, "int")
routetogrid[:,:,:] = part.gather(rgrid) #
# # basin area
rtobasin = procgrid.landscatter(self.route_tobasin[:,:].astype(vtyp), order='F') self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "basin_area", "area of basin", "m^2", self.routing_area[:,:], vtyp, "float")
rtobasin = rtobasin.astype(vtyp, copy=False) # route number into basin
rtobasin[rtobasin >= RPP.IntFillValue] = NCFillValue self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "routenbintobas", "Number of basin into current one", "-", self.route_nbintobas[:], vtyp, "int")
if part.rank == 0 : #
routetobasin = outnf.createVariable("routetobasin", vtyp, ('htu','y','x'), fill_value=NCFillValue) # original number into basin
routetobasin.title = "Basin in to which the water goes" 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")
routetobasin.units = "-" #
else : # latitude of outlet
routetobasin = np.zeros((1,1,1)) self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "outletlat", "Latitude of Outlet", "degrees north", self.route_outlet[:,:,0], vtyp, "float")
routetobasin[:,:,:] = part.gather(rtobasin) # 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")
rid = procgrid.landscatter(self.global_basinid[:,:].astype(vtyp), order='F') # type of outlet
rid[rid >= RPP.IntFillValue] = NCFillValue self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "outlettype", "Type of outlet", "code", self.route_type[:,:], vtyp, "float")
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)
# #
olat = procgrid.landscatter(self.route_outlet[:,:,0].astype(vtyp), order='F') # topographic index
olat[np.isnan(olat)] = NCFillValue self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "topoindex", "Topographic index of the retention time", "m", self.topo_resid[:,:], vtyp, "float")
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)
# #
olon = procgrid.landscatter(self.route_outlet[:,:,1].astype(vtyp), order='F')
olon[np.isnan(olon)] = NCFillValue # Inflow number
if part.rank == 0 : self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "route_innum", "Number of inflow", "-", self.route_innum[:,:], vtyp, "int")
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)
# #
otype = procgrid.landscatter(self.route_type[:,:].astype(vtyp), order='F') # Inflow grid
otype[np.isnan(otype)] = NCFillValue #gingrid = part.l2glandindex(self.inflow_grid[:,:,:inflow_size])
if part.rank == 0 : self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow', 'htu','y','x'), "route_ingrid", "Grid from which the water flows", "-", self.route_ingrid[:,:,:], vtyp, "int")
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)
# #
tind = procgrid.landscatter(self.topo_resid[:,:].astype(vtyp), order='F') # Inflow basin
tind[np.isnan(tind)] = NCFillValue self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow', 'htu','y','x'), "route_inbasin", "Basin from which the water flows", "-", self.route_inbasin[:,:,:], vtyp, "int")
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)
# #
# Save centre of gravity of HTU # Save centre of gravity of HTU
# #
cg = procgrid.landscatter(self.routing_cg[:,:,:].astype(vtyp), order='F') # Check if it works
cg[np.isnan(cg)] = NCFillValue 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")
if part.rank == 0 :
CG_lon = outnf.createVariable("CG_lon", vtyp, ('htu','y','x'), fill_value=NCFillValue) 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")
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,:,:,:])
# #
# Save the fetch of each basin # Save the fetch of each basin
# #
fe = procgrid.landscatter(self.routing_fetch[:,:].astype(vtyp), order='F') self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "fetch", "Fetch contributing to each HTU", "m^2", self.routing_fetch[:,:], vtyp, "float")
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)
# #
if part.rank == 0 : # Close the file
if part.rank == 0:
outnf.close() outnf.close()
# #
return return
......
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