From 3aaa4036c39d879b8196592a2232e3a5cf59d363 Mon Sep 17 00:00:00 2001
From: POLCHER Jan <jan.polcher@lmd.jussieu.fr>
Date: Fri, 29 May 2020 21:10:16 +0200
Subject: [PATCH] Improved the graph files so that they are compatible with the
IOIPSL routines used to read them in the model.
---
Interface.py | 38 +++++++++++++++++++-------------------
ModelGrid.py | 7 ++++++-
2 files changed, 25 insertions(+), 20 deletions(-)
diff --git a/Interface.py b/Interface.py
index 04e2a82..2e02691 100644
--- a/Interface.py
+++ b/Interface.py
@@ -88,7 +88,7 @@ def addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorne
lon.units="grid box centre degrees east"
lon.title="Longitude"
lon.axis="X"
- lon[:,:] = longitude[:,:]
+ lon[:,:] = globalgrid.lonmat[:,:]
#
# Latitude
latitude = part.gather(procgrid.lat_full)
@@ -98,7 +98,7 @@ def addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorne
lat.standard_name="grid latitude"
lat.title="Latitude"
lat.axis="Y"
- lat[:] = latitude[:,:]
+ lat[:,:] = globalgrid.latmat[:,:]
#
# Bounds of grid box
#
@@ -652,7 +652,7 @@ class HydroGraph :
outnf.createDimension('x', globalgrid.ni)
outnf.createDimension('y', globalgrid.nj)
outnf.createDimension('land', len(procgrid.area))
- outnf.createDimension('htu', self.nbasmax)
+ outnf.createDimension('z', self.nbasmax)
outnf.createDimension('bnd', nbcorners)
outnf.createDimension('inflow', self.max_inflow)
else :
@@ -688,18 +688,18 @@ class HydroGraph :
################
#
# 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")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z', '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")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z', '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")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z', '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")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_area", "area of basin", "m^2", self.routing_area[:,:], vtyp, "float")
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "basin_orog", "Mean orography", "m", self.routing_orog[:,:], vtyp, "float")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_orog", "Mean orography", "m", self.routing_orog[:,:], vtyp, "float")
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "basin_floodp", "area of floodplains", "m^2", self.routing_floodp[:,:], vtyp, "float")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_floodp", "area of floodplains", "m^2", self.routing_floodp[:,:], 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")
@@ -708,37 +708,37 @@ class HydroGraph :
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")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z','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")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z','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")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outlettype", "Type of outlet", "code", self.route_type[:,:], vtyp, "float")
#
# 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")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "topoindex", "Topographic index of the retention time", "m", self.topo_resid[:,:], vtyp, "float")
#
# Inflow number
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "route_innum", "Number of inflow", "-", self.route_innum[:,:], vtyp, "int")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "route_innum", "Number of inflow", "-", self.route_innum[:,:], vtyp, "int")
#
# 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")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow', 'z','y','x'), "route_ingrid", "Grid from which the water flows", "-", self.route_ingrid[:,:,:], vtyp, "int")
#
# 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")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow', 'z','y','x'), "route_inbasin", "Basin from which the water flows", "-", self.route_inbasin[:,:,:], vtyp, "int")
#
# Save centre of gravity of HTU
#
# 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, ('z','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")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z','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
#
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "fetch", "Fetch contributing to each HTU", "m^2", self.routing_fetch[:,:], vtyp, "float")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "fetch", "Fetch contributing to each HTU", "m^2", self.routing_fetch[:,:], vtyp, "float")
#
# Close the file
if part.rank == 0:
diff --git a/ModelGrid.py b/ModelGrid.py
index 81b7a4c..8a1233c 100644
--- a/ModelGrid.py
+++ b/ModelGrid.py
@@ -398,6 +398,8 @@ class GlobalGrid :
self.nj, self.ni = lon_full.shape
self.jgstart = 0
self.igstart = 0
+ self.lonmat = lon_full
+ self.latmat = lat_full
else :
dist=np.sqrt(np.power(lon_full-np.min(lonrange),2) + np.power(lat_full-np.min(latrange),2))
jmin,imin=np.unravel_index(np.argmin(dist, axis=None), dist.shape)
@@ -408,7 +410,10 @@ class GlobalGrid :
self.jgstart = jmin
self.ni = imax-imin+1
self.igstart = imin
- print("Shape : ", lon_full.shape, "N = ", self.nj, self.ni, " Start :", self.jgstart, self.igstart)
+ self.lonmat = lon_full[self.jgstart:self.jgstart+self.nj,self.igstart:self.igstart+self.ni]
+ self.latmat = lat_full[self.jgstart:self.jgstart+self.nj,self.igstart:self.igstart+self.ni]
+ INFO("Shape of full grid kept in GlobalGrid : "+str(self.lonmat.shape))
+
self.land = getland(geo, self.igstart, self.ni, self.jgstart, self.nj)
self.nbland = int(np.sum(self.land))
--
GitLab