diff --git a/.gitignore b/.gitignore
index 1559e986d01fb016e518648d86e6a06b1dc4a123..aca73777d794e08652ab8c25b007386b073e9881 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 1b5f33ad6d89be844dd9f2016e53df92b04e01d5..ede10559d96a3166f263ba02a5847259d4e0f1ad 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 332c569149d1e0c3ef4f373f2f544c38f04b6443..09a8b0db4c6c80ee6540de5c9acb63b995d6f316 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