diff --git a/Diags/plotgraphs.py b/Diags/plotgraphs.py
new file mode 100644
index 0000000000000000000000000000000000000000..c3d3af55f9e9154f4d6783f71d49359e4a7fd1b7
--- /dev/null
+++ b/Diags/plotgraphs.py
@@ -0,0 +1,115 @@
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+from mpl_toolkits.basemap import Basemap, cm
+import numpy as np
+from netCDF4 import Dataset
+
+NCfile="../MEDCORDEX_test_graph.nc"
+GraphicFile = "test.png"
+title="Test"
+#
+#####################################################################
+#
+class LandGrid :
+ def __init__(self, land) :
+ nj,ni = land.shape
+ self.indP=[]
+ self.nbland = 0
+ for i in range(ni) :
+ for j in range(nj) :
+ if (land[j,i] > 0 ) :
+ self.nbland = self.nbland+1
+ self.indP.append([j,i])
+ return
+
+ def land2ij(self, il) :
+ j,i = self.indP[il]
+ return j,i
+#
+#####################################################################
+#
+#
+nc = Dataset(NCfile, 'r')
+
+lon_full = nc.variables["lon"][:,:]
+lat_full = nc.variables["lat"][:,:]
+lon_bnd = nc.variables["lon_bnd"][:,:,:]
+lat_bnd = nc.variables["lat_bnd"][:,:,:]
+land = nc.variables["land"][:,:]
+routetogrid = nc.variables["routetogrid"][:,:,:]
+routetohtu = nc.variables["routetobasin"][:,:,:]
+cglon = nc.variables["CG_lon"][:,:,:]
+cglat = nc.variables["CG_lat"][:,:,:]
+
+nc.close()
+
+nhtu, nj, ni = routetogrid.shape
+landindex = LandGrid(land)
+print "Nb land points = ", landindex.nbland, 'Nb HTU : ', nhtu
+
+#
+# Get some average size of the grid boxes.
+#
+lon_res = np.nanmean(np.abs(np.gradient(lon_bnd, axis=0)))
+lat_res = np.nanmean(np.abs(np.gradient(lat_bnd, axis=0)))
+#
+#
+#
+fig = plt.figure(1)
+ax = fig.add_subplot(111)
+ax.set_title(title)
+ax.set_xlabel('Longitude',labelpad=20)
+ax.set_ylabel('Latitude',labelpad=20)
+ax.axis([0, 10, 0, 10])
+#
+m = Basemap(projection='cyl', resolution='i', \
+ llcrnrlon=np.nanmin(lon_bnd)-lon_res, llcrnrlat=np.nanmin(lat_bnd)-lat_res, \
+ urcrnrlon=np.nanmax(lon_bnd)+lon_res, urcrnrlat=np.nanmax(lat_bnd)+lat_res)
+
+m.drawcoastlines(linewidth=0.5)
+m.drawcountries(linewidth=0.5)
+m.drawparallels(np.arange(-180.,180.,0.1),labels=[1,0,0,0]) # draw parallels
+m.drawmeridians(np.arange(-80.,80.,0.1),labels=[0,0,0,1]) # draw meridians
+#
+cmap = 'RdBu'
+cmap = 'Pastel1'
+cmap = 'prism'
+#
+#
+#
+for i in range(ni) :
+ for j in range(nj) :
+ if land[j,i] > 0 :
+ lo = np.append(lon_bnd[:,j,i],lon_bnd[0,j,i])
+ la = np.append(lat_bnd[:,j,i],lat_bnd[0,j,i])
+ xl,yl=m(lo, la)
+ plt.plot(xl,yl, color="m")
+#
+# Draw the graph
+#
+for il in range(landindex.nbland) :
+ j,i = landindex.land2ij(il)
+ print il,"== routetogrid ==",routetogrid[:,j,i]
+ print il,"== routetohtu ==",routetohtu[:,j,i]
+ for ih in range(nhtu) :
+ lo = cglon[ih,j,i]
+ la = cglat[ih,j,i]
+ xl,yl=m(lo, la)
+ if int(routetohtu[ih,j,i]-1) < nhtu :
+ plt.plot(xl,yl, marker=11)
+ else :
+ plt.plot(xl,yl, marker="o")
+ if int(routetohtu[ih,j,i]-1) < nhtu :
+ ile = int(routetogrid[ih,j,i]-1)
+ ie,je = landindex.land2ij(ile)
+ ihe = int(routetohtu[ih,j,i]-1)
+ loe = cglon[ihe,je,ie]
+ lae = cglat[ihe,je,ie]
+ xle,yle=m(loe, lae)
+ plt.plot([xl,xle],[yl,yle], color="b")
+#
+#
+#
+plt.savefig(GraphicFile, dpi=450)
+plt.close()
diff --git a/Interface.py b/Interface.py
index 6ba11149b971a6197f8cfdc38637c484d02d595a..1b5f33ad6d89be844dd9f2016e53df92b04e01d5 100644
--- a/Interface.py
+++ b/Interface.py
@@ -177,6 +177,7 @@ class HydroGraph :
NCFillValue=1.0e20
vtyp=np.float64
cornerind=[0,2,4,6]
+ nbcorners = len(cornerind)
#
if part.rank == 0 :
outnf=Dataset(filename, 'w', format='NETCDF4_CLASSIC')
@@ -185,7 +186,7 @@ class HydroGraph :
outnf.createDimension('y', globalgrid.nj)
outnf.createDimension('land', len(procgrid.area))
outnf.createDimension('htu', self.nbasmax)
- outnf.createDimension('bnd', 4)
+ outnf.createDimension('bnd', nbcorners)
#
# Coordinates
#
@@ -197,11 +198,6 @@ class HydroGraph :
lon.title="Longitude"
lon.axis="X"
lon[:,:] = longitude[:,:]
- # # Longitude bounds
- # lonbnd=outnf.createVariable("lon_bnd", vtyp, ('bnd','y','x'), fill_value=NCFillValue)
- # lonbnd.units="grid box corners degrees east"
- # lonbnd.title="Longitude Corners"
- # lonbnd[:,:] = np.array(procgrid.polyll)[:,cornerind,0]
#
# Latitude
latitude = part.gather(procgrid.lat_full)
@@ -212,12 +208,28 @@ class HydroGraph :
lat.title="Latitude"
lat.axis="Y"
lat[:] = latitude[:,:]
- # # Latitude bounds
- # latbnd=outnf.createVariable("lat_bnd", vtyp, ('bnd','land'), fill_value=NCFillValue)
- # latbnd.units="grid box corners degrees north"
- # latbnd.title="Latitude Corners"
- # latbnd[:,:] = np.array(procgrid.polyll)[:,cornerind,1]
#
+ # Bounds of grid box
+ #
+ llonpoly=np.zeros((nbcorners,procgrid.nbland))
+ llatpoly=np.zeros((nbcorners,procgrid.nbland))
+ for i in range(procgrid.nbland) :
+ llonpoly[:,i] = [procgrid.polyll[i][ic][0] for ic in cornerind]
+ llatpoly[:,i] = [procgrid.polyll[i][ic][1] for ic in cornerind]
+ lon_bnd = procgrid.landscatter(llonpoly)
+ lat_bnd = procgrid.landscatter(llatpoly)
+ if part.rank == 0 :
+ lonbnd=outnf.createVariable("lon_bnd", vtyp, ('bnd','y','x'), fill_value=NCFillValue)
+ lonbnd.units="grid box corners degrees east"
+ lonbnd.title="Longitude of Corners"
+ latbnd=outnf.createVariable("lat_bnd", vtyp, ('bnd','y','x'), fill_value=NCFillValue)
+ latbnd.units="grid box corners degrees north"
+ latbnd.title="Latitude of Corners"
+ else :
+ lonbnd= np.zeros((1,1,1))
+ latbnd= np.zeros((1,1,1))
+ lonbnd[:,:,:] = part.gather(lon_bnd[:,:,:])
+ latbnd[:,:,:] = part.gather(lat_bnd[:,:,:])
#
# Land sea mask
#
diff --git a/SampleScripts/BuildHTUs.pbs b/SampleScripts/BuildHTUs.pbs
index 87481b26f3fedb466071d5f103b984960263c4bc..1f2c909fe96b9615e571242ff66faeee0a6c17f8 100644
--- a/SampleScripts/BuildHTUs.pbs
+++ b/SampleScripts/BuildHTUs.pbs
@@ -5,8 +5,8 @@
#PBS -j oe
#PBS -l nodes=1:ppn=42
#PBS -l walltime=48:00:00
-#PBS -l mem=120gb
-#PBS -l vmem=120gb
+#PBS -l mem=160gb
+#PBS -l vmem=160gb
#
cd ${PBS_O_WORKDIR}/..
diff --git a/diagplots.py b/diagplots.py
index 946db565f19c49c295b432a39c11a35e4849775e..39d7b7749c79b6bca8d58cf930e56f05ba22444d 100644
--- a/diagplots.py
+++ b/diagplots.py
@@ -5,6 +5,11 @@ import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, cm
import numpy as np
+import getargs
+log_master, log_world = getargs.getLogger()
+INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
+INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
+
inc=np.zeros((8,2), dtype=np.int8)
inc[0,0] = 1; inc[0,1] =0
inc[1,0] = 1; inc[1,1] =1
@@ -204,7 +209,7 @@ def SUPERMESHPLOT(filename, loninterval, latinterval, hoverlap, hsuper, atmpolys
lat = hoverlap.lat_bx
basin_sz = hsuper.basin_sz
basin_pts = hsuper.basin_pts
- outlet = hsuper.basin_coor
+ outlet = hsuper.basin_outcoor
#
if np.nanmin(lon) <= max(loninterval) and np.nanmax(lon) > min(loninterval) and \
np.nanmin(lat) <= max(latinterval) and np.nanmax(lat) > min(latinterval) :
@@ -265,6 +270,7 @@ def SUPERMESHPLOT(filename, loninterval, latinterval, hoverlap, hsuper, atmpolys
if ibo == ib :
# Same grid box
if hsuper.outflow_basin[ib,iz] < undef_int :
+ INFO("hsuper.outflow_basin[ib,iz] = "+str(hsuper.outflow_basin[ib,iz]))
ERROR("We cannot be here as we flow into an HTU of the same grid but we do not know the HTU index")
else :
izo = hsuper.outflow_basin[ib,iz]-1
diff --git a/run.def b/run.def
index 8c7083b3605a8bd7e7295b935443e37b8deb0bed..1cca8b5a374462067cfd83f6b8f1b49e42f6b090 100644
--- a/run.def
+++ b/run.def
@@ -24,5 +24,5 @@ GraphFile = MEDCORDEX_test_graph.nc
# Diagnostics
# You need to provide an interval in longitude and Latitude.
#
-DiagLon = 3.9, 3.9
-DiagLat = 40.0, 40.0
+DiagLon = 2.3, 3.5
+DiagLat = 39.0, 40.1