From 4934a56caab024b3d2448088580778bb808eb5ca Mon Sep 17 00:00:00 2001
From: Jan Polcher <jan.polcher@lmd.jussieu.fr>
Date: Sat, 15 Jun 2019 19:42:39 +0200
Subject: [PATCH] Added a function to go from local land index to global land
index.
---
.gitignore | 1 +
Interface.py | 5 ++++-
Partition.py | 36 +++++++++++++++++++++++++++++++-----
3 files changed, 36 insertions(+), 6 deletions(-)
diff --git a/.gitignore b/.gitignore
index 1559e98..aca7377 100644
--- a/.gitignore
+++ b/.gitignore
@@ -23,6 +23,7 @@ xcuserdata/
*.mod
*.pdone
*.bak
+*.jnl
# Backup files
*~.nib
diff --git a/Interface.py b/Interface.py
index 1b5f33a..ede1055 100644
--- a/Interface.py
+++ b/Interface.py
@@ -265,7 +265,10 @@ class HydroGraph :
routingarea = np.zeros((1,1,1))
routingarea[:,:,:] = part.gather(rarea)
#
- rgrid = procgrid.landscatter(self.route_togrid[:,:], order='F')
+ # The field route_togrid is with indices on the local grid. That needs to be converted to the global grid.
+ #
+ grgrid = part.l2glandindex(self.route_togrid[:,:])
+ rgrid = procgrid.landscatter(grgrid, order='F')
rgrid = rgrid.astype(vtyp, copy=False)
rgrid[rgrid >= RPP.IntFillValue] = NCFillValue
if part.rank == 0 :
diff --git a/Partition.py b/Partition.py
index 332c569..09a8b0d 100644
--- a/Partition.py
+++ b/Partition.py
@@ -158,17 +158,28 @@ def haloreceivelist(halosource_map, rank) :
#
# Get 1D indices for the land points
#
-def landindexmap(land) :
- nj,ni=land.shape
+def landindexmap(istart, ni, jstart, nj, land) :
+ gnj,gni=land.shape
+ gindland=np.zeros((gnj,gni), dtype=np.int32)
+ gindland[:,:]=-1
+ n=0
+ for i in range(gni) :
+ for j in range(gnj) :
+ if (land[j,i] > 0 ) :
+ gindland[j,i] = n
+ n += 1
+ #
indland=np.zeros((nj,ni), dtype=np.int32)
indland[:,:]=-1
+ local2global=[]
n=0
for i in range(ni) :
for j in range(nj) :
- if (land[j,i] > 0 ) :
+ if (land[jstart+j,istart+i] > 0 ) :
indland[j,i] = n
+ local2global.append(gindland[jstart+j,istart+i])
n += 1
- return indland
+ return n, indland, np.array(local2global, dtype=np.int32)
#
# Convert indices to local
#
@@ -230,7 +241,7 @@ class partition :
self.ihstart = part[rank]["ihstart"]
self.jhstart = part[rank]["jhstart"]
#
- landindmap = landindexmap(land[self.jhstart:self.jhstart+self.njh,self.ihstart:self.ihstart+self.nih])
+ self.nbland, landindmap, self.l2glandind = landindexmap(self.ihstart, self.nih, self.jhstart, self.njh, land)
#
if wunit != "None" :
wunit.write("Offsets with halo (j-i) : "+str(self.jhstart)+"-"+str(self.ihstart)+'\n')
@@ -372,3 +383,18 @@ class partition :
xf = modelgrid.landscatter(xl, order=order)
xout = self.gather(xf)
return xout
+ #
+ # Convert local index of land points to global index
+ #
+ def l2glandindex(self, x) :
+ nl,nh = x.shape
+ y = np.zeros(x.shape, dtype=x.dtype)
+ if nl == self.nbland :
+ for i in range(nl) :
+ for j in range(nh) :
+ # Land indices are in FORTRAN !!
+ y[i,j] = self.l2glandind[x[i,j]-1]+1
+ else :
+ ERROR("The first dimension does not have the length of the number of land points")
+ sys.exit()
+ return y
--
GitLab