Index: codes/icosagcm/devel/Python/dynamico/meshes.py
===================================================================
--- codes/icosagcm/devel/Python/dynamico/meshes.py	(revision 868)
+++ codes/icosagcm/devel/Python/dynamico/meshes.py	(revision 869)
@@ -22,5 +22,5 @@
 INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
 
-radian=180/math.pi
+radian=180/math.pi # convert from radians to degrees
 
 #------------------- Hybrid mass-based coordinate -------------
@@ -131,4 +131,5 @@
         Riv2                           = zeros((nx*ny,4)  )
         wee                            = zeros((2*nx*ny,4))
+        primal_bounds_lon, primal_bounds_lat, dual_bounds_lon, dual_bounds_lat = zeros( (nx*ny,4), 4)
 
         primal_deg,dual_deg                 = indices( nx*ny,   2)
@@ -146,17 +147,17 @@
                 ij=index(x,y)-1
                 lon_i[ij], lat_i[ij] = x, y
-                put(index(x,y)-1, primal_deg,(primal_edge,primal_vertex,primal_ne),
-                    ((indexu(x,y),index(x,y),1), 
-                    (indexv(x,y),index(x-1,y),1),
-                    (indexu(x-1,y),index(x-1,y-1),-1),
-                    (indexv(x,y-1),index(x,y-1),-1) ))
+                put(ij, primal_deg,(primal_edge,primal_vertex,primal_ne, primal_bounds_lon, primal_bounds_lat),
+                    ((indexu(x,y),   index(x,y),      1, x+.5, y+.5), 
+                     (indexv(x,y),   index(x-1,y),    1, x-.5, y+.5),
+                     (indexu(x-1,y), index(x-1,y-1), -1, x-.5, y-.5),
+                     (indexv(x,y-1), index(x,y-1),   -1, x+.5, y-.5) ))
                 # dual cells
                 ij=index(x,y)-1
                 lon_v[ij], lat_v[ij] = x+.5, y+.5
-                put(ij,dual_deg,(dual_edge,dual_vertex,dual_ne,Riv2),
-                   ((indexv(x+1,y),index(x,y),1,.25),
-                    (indexu(x,y+1),index(x+1,y),-1,.25),
-                    (indexv(x,y),index(x+1,y+1),-1,.25),
-                    (indexu(x,y),index(x,y+1),1,.25) ))
+                put(ij,dual_deg,(dual_edge,dual_vertex,dual_ne,Riv2, dual_bounds_lon, dual_bounds_lat),
+                    ((indexv(x+1,y),index(x,y),      1, .25, x,    y),
+                     (indexu(x,y+1),index(x+1,y),   -1, .25, x+1., y),
+                     (indexv(x,y),  index(x+1,y+1), -1, .25, x+1., y+1.),
+                     (indexu(x,y),  index(x,y+1),    1, .25, x,    y+1.) ))
                 # edges :
                 # left and right are adjacent primal cells
@@ -193,5 +194,7 @@
         self.set_members(locals(), 'primal_num',  'dual_num', 'edge_num',
                          'primal_deg', 'primal_edge', 'primal_ne', 'primal_vertex',
+                         'primal_bounds_lon', 'primal_bounds_lat',
                          'dual_deg', 'dual_edge', 'dual_ne', 'dual_vertex',
+                         'dual_bounds_lon', 'dual_bounds_lat',
                          'left','right','down','up',
                          'trisk_deg','trisk','Riv2','wee',
@@ -253,9 +256,13 @@
                      [("primal_edge", "i4", self.primal_edge),
                       ("primal_ne", "i4", self.primal_ne),
-                      ("primal_vertex", "i4", self.primal_vertex)] )
+                      ("primal_vertex", "i4", self.primal_vertex),
+                      ("primal_bounds_lon", "f8", self.primal_bounds_lon),
+                      ("primal_bounds_lat", "f8", self.primal_bounds_lat)] )
 
         create_vars( ("dual_cell","dual_edge_or_vertex"), 
                      [("dual_edge", "i4", self.dual_edge),
                       ("dual_vertex","i4",self.dual_vertex),
+                      ("dual_bounds_lon", "f8", self.dual_bounds_lon),
+                      ("dual_bounds_lat", "f8", self.dual_bounds_lat),
                       ("dual_ne","i4",self.dual_ne),
                       ("Riv2","f8",self.Riv2)] )
@@ -347,5 +354,7 @@
         'lat_e':'latEdge','lon_e':'lonEdge','angle_e':'angleEdge',
         'wee':'weightsOnEdge','Riv2':'kiteAreasOnVertex',
-        'primal_ne':'signOnCell', 'dual_ne':'signOnVertex'}
+        'primal_ne':'signOnCell', 'dual_ne':'signOnVertex',
+        'primal_bounds_lon':None, 'primal_bounds_lat':None,
+        'dual_bounds_lon':None, 'dual_bounds_lat':None }
     def __init__(self,gridfilename):
         try:
@@ -377,5 +386,5 @@
         # general case
         name=self.translate[name]
-        return parallel.PArray(pdim, self.nc.variables[name])
+        return None if name is None else parallel.PArray(pdim, self.nc.variables[name])
     def normalize(self, mesh):
         (trisk_deg, trisk, Ai,
@@ -397,5 +406,17 @@
         # this is now done by apply_map
         # mesh.le, mesh.de, mesh.Av, mesh.Ai = radius*le, radius*de, r2*Av, r2*Ai
-        
+    def get_bounds(self, mesh):
+        def bounds(deg, vertex, lon, lat):
+            num, maxdeg = deg.size, deg.max()
+            vertex_lon, vertex_lat = np.zeros((num, maxdeg)), np.zeros((num, maxdeg))
+            for cell in range(num):
+                for ivertex in range(deg[cell]):
+                    vtx = vertex[cell,ivertex] - 1
+                    vertex_lon[cell, ivertex] = lon[vtx]
+                    vertex_lat[cell, ivertex] = lat[vtx]
+            return vertex_lon, vertex_lat
+        mesh.primal_bounds_lon, mesh.primal_bounds_lat = bounds(mesh.primal_deg, mesh.primal_vertex, mesh.lon_v, mesh.lat_v)
+        mesh.dual_bounds_lon, mesh.dual_bounds_lat = bounds(mesh.dual_deg, mesh.dual_vertex, mesh.lon_i, mesh.lat_i)
+
 def compute_ne(num,deg,edges,left,right):
     ne = np.zeros_like(edges)
@@ -450,4 +471,6 @@
         self.lon_v, self.lat_v = mapping.map(lon_v, lat_v)
         self.lon_e, self.lat_e = mapping.map(lon_e, lat_e)
+        self.primal_bounds_lon, self.primal_bounds_lat = mapping.map( self.primal_bounds_lon, self.primal_bounds_lat)
+        self.dual_bounds_lon, self.dual_bounds_lat = mapping.map( self.dual_bounds_lon, self.dual_bounds_lat)
         self.ref_lon_i, self.ref_lat_i = lon_i, lat_i
         self.ref_lon_e, self.ref_lat_e = lon_e, lat_e
@@ -500,4 +523,5 @@
         primal_deg, primal_vertex, primal_edge, primal_ne, lon_i, lat_i, Ai = get_pvars(
             dim_primal, 'primal_deg', 'primal_vertex', 'primal_edge', 'primal_ne', 'lon_i', 'lat_i', 'Ai' )
+        primal_bounds_lon, primal_bounds_lat = get_pvars(dim_primal, 'primal_bounds_lon', 'primal_bounds_lat')
         edge_deg, trisk_deg, left_right, down_up, trisk, le, de, wee, lon_e, lat_e, angle_e = get_pvars(
              dim_edge, 'edge_deg', 'trisk_deg', 'edge_left_right', 'edge_down_up', 
@@ -505,4 +529,5 @@
         dual_deg, dual_vertex, dual_edge, dual_ne, Riv2, lon_v, lat_v, Av = get_pvars(
             dim_dual, 'dual_deg', 'dual_vertex', 'dual_edge', 'dual_ne', 'Riv2', 'lon_v', 'lat_v', 'Av')
+        dual_bounds_lon, dual_bounds_lat = get_pvars(dim_dual, 'dual_bounds_lon', 'dual_bounds_lat')
         if gridfile.meshtype == 'curvilinear':
             self.primal_i, self.primal_j = get_pvars(dim_primal, 'primal_i','primal_j')
@@ -521,6 +546,8 @@
         self.edge2edge   = Stencil_glob(trisk_deg, trisk, wee)
         self.set_members(locals(), 'dim_primal', 'dim_dual', 'dim_edge',
-                         'primal_deg', 'primal_vertex', 'primal_edge', 
+                         'primal_deg', 'primal_vertex', 'primal_edge',
+                         'primal_bounds_lon', 'primal_bounds_lat',
                          'trisk_deg', 'trisk', 'dual_deg', 'dual_edge',
+                         'dual_bounds_lon', 'dual_bounds_lat',
                          'lon_v', 'lat_v', 'Av', 'lon_i', 'lat_i', 'Ai', 
                          'le', 'de', 'lon_e', 'lat_e', 'angle_e')
@@ -640,4 +667,6 @@
         # read from mesh file : coordinates, lengths and areas in reference domain
         lon_i, lat_i, Ai = pmesh.get(get_all_cells, 'lon_i', 'lat_i', 'Ai')
+
+
         lon_v, lat_v, Av = pmesh.get(get_all_duals, 'lon_v', 'lat_v', 'Av')
         le, de, lon_e, lat_e, angle_e = pmesh.get(get_all_edges, 'le', 'de', 'lon_e', 'lat_e', 'angle_e')
@@ -657,4 +686,13 @@
                          'le', 'de', 'le_de', 'lon_e', 'lat_e', 'angle_e')
 
+        # obtain bounds of primal and dual cells (for XIOS)
+        if pmesh.primal_bounds_lon is not None:
+            INFO('Reading cell bounds from mesh file.')
+            self.primal_bounds_lon, self.primal_bounds_lat = pmesh.get(get_all_cells, 'primal_bounds_lon', 'primal_bounds_lat')
+            self.dual_bounds_lon, self.dual_bounds_lat = pmesh.get(get_all_duals, 'dual_bounds_lon', 'dual_bounds_lat')
+        else:
+            INFO('Mesh file does not contain information abound cell bounds.')
+            pmesh.gridfile.get_bounds(self)
+
         # normalize data read from file depending on file format
         pmesh.gridfile.normalize(self)
@@ -662,4 +700,5 @@
         self.apply_map(mapping)
         # now latitudes and longitudes refer to the physical domain
+
         self.fv = mapping.coriolis(self.lon_v,self.lat_v) # Coriolis parameter
         # propagate info to Fortran side
Index: codes/icosagcm/devel/Python/test/py/dump_partition.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/dump_partition.py	(revision 868)
+++ codes/icosagcm/devel/Python/test/py/dump_partition.py	(revision 869)
@@ -107,8 +107,11 @@
     def create_vars(self,f,dimname, info):
         for vname, vtype, vdata in info:
+            vdata = np.asarray(vdata)
+            print('creating variable %s with shape %s.' %(vname, vdata.shape))
             var = f.createVariable(vname,vtype,dimname)
             var[:] = vdata
     
     def meshwrite(self,name):
+        mesh = self._mesh
         """ Writing Mesh information to disk for Fortran segment"""
     
@@ -124,4 +127,12 @@
     
         #----------------DEFINING DIMENSIONS--------------------
+        max_primal_deg = mesh.primal_deg.max()
+        max_dual_deg = mesh.dual_deg.max()
+        max_trisk_deg = mesh.trisk_deg.max()
+        def crop(dim, *data) : return [ x[:,0:dim] for x in data ]
+        mesh.primal_edge, mesh.primal_ne, mesh.primal_vertex = crop(max_primal_deg, mesh.primal_edge, mesh.primal_ne, mesh.primal_vertex)
+        mesh.dual_edge, mesh.dual_ne = crop(max_dual_deg, mesh.dual_edge, mesh.dual_ne)
+        mesh.trisk, mesh.wee = crop(max_trisk_deg, mesh.trisk, mesh.wee)
+
         for dimname, dimsize in [("primal_cell", self._mesh.primal_deg.size),
                                     ("dual_cell", self._mesh.dual_deg.size),
@@ -144,8 +155,10 @@
     
         self.create_vars(f,("primal_cell","primal_edge_or_vertex"),
-                    [("primal_vertex","i4",self._mesh.primal_vertex),
-                    ("primal_edge","i4",self._mesh.primal_edge),
-                    ("primal_ne","i4",self._mesh.primal_ne) ])
-    
+                         [("primal_vertex","i4",self._mesh.primal_vertex),
+                          ("primal_bounds_lon","f8",self._mesh.primal_bounds_lon),
+                          ("primal_bounds_lat","f8",self._mesh.primal_bounds_lat),
+                          ("primal_edge","i4",self._mesh.primal_edge),
+                          ("primal_ne","i4",self._mesh.primal_ne) ])
+        
         self.create_vars(f,"dual_cell",
                     [("dual_deg","i4",self._mesh.dual_deg),
@@ -161,7 +174,9 @@
         self.create_vars(f,("dual_cell","dual_edge_or_vertex"),
                     [("dual_edge","i4",self._mesh.dual_edge),
-                    ("dual_vertex","i4",self._mesh.dual_vertex),
-                    ("dual_ne","i4",self._mesh.dual_ne),
-                    ("Riv2","f8",self._mesh.Riv2)])
+                     ("dual_vertex","i4",self._mesh.dual_vertex),
+                     ("dual_bounds_lon","f8",self._mesh.dual_bounds_lon),
+                     ("dual_bounds_lat","f8",self._mesh.dual_bounds_lat),
+                     ("dual_ne","i4",self._mesh.dual_ne),
+                     ("Riv2","f8",self._mesh.Riv2)])
     
         self.create_vars(f,"edge",
Index: codes/icosagcm/devel/Python/test/py/write_Cartesian_mesh.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/write_Cartesian_mesh.py	(revision 868)
+++ codes/icosagcm/devel/Python/test/py/write_Cartesian_mesh.py	(revision 869)
@@ -1,15 +1,11 @@
 from dynamico import getargs
-getargs.add("--nx", type=int,
-                    default=64, choices=None,
-                    help="number of x points")
-getargs.add("--ny", type=int,
-                    default=64, choices=None,
-                    help="number of y points")
-getargs.add("--Lx", type=float,
-                    default=8., choices=None,
-                    help="Lx")
-getargs.add("--Ly", type=float,
-                    default=8., choices=None,
-                    help="Ly")
+getargs.add("--nx", type=int, default=64, 
+                    help="number of points along the x axis")
+getargs.add("--ny", type=int, default=0,
+                    help="number of points along the y axis, defaults to nx")
+getargs.add("--Lx", type=float, default=0., 
+                    help="domain size along x, defaults to nx")
+getargs.add("--Ly", type=float, default=0., 
+                    help="domain size along y, defaults to ny")
 
 args = getargs.parse()
@@ -22,4 +18,7 @@
 
 nx, ny, Lx, Ly, llm, nqdyn = args.nx, args.ny,args.Lx, args.Ly, 1, 1
+if ny==0 : ny=nx
+if Lx==0. : Lx=nx
+if Ly==0. : Ly=ny
 
 dx,dy=Lx/nx,Ly/ny
