diff --git a/Interface.py b/Interface.py
index 04e2a82e991fb636b3e242da120027190ea79d80..2e026913c8631b4b6472699f4ccabf7026b800d4 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 81b7a4c35bb01720219b9f2d0fac9a2af5184c1f..8a1233ce27e92ddd39237bcccc893556d03bcd19 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))