Commit 4934a56c authored by POLCHER Jan's avatar POLCHER Jan 🚴🏾
Browse files

Added a function to go from local land index to global land index.

parent 02825e89
......@@ -23,6 +23,7 @@ xcuserdata/
*.mod
*.pdone
*.bak
*.jnl
# Backup files
*~.nib
......
......@@ -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 :
......
......@@ -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
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