Index: codes/icosagcm/devel/Python/dynamico/meshes.py
===================================================================
--- codes/icosagcm/devel/Python/dynamico/meshes.py	(revision 758)
+++ codes/icosagcm/devel/Python/dynamico/meshes.py	(revision 759)
@@ -182,5 +182,6 @@
                          'trisk_deg','trisk','Riv2','wee',
                          'Aiv','lon_i','lat_i','lon_v','lat_v',
-                         'le_de','le','de','lon_e','lat_e','angle_e')
+                         'le_de','le','de','lon_e','lat_e','angle_e',
+                         'primal_i','primal_j')
     def ncwrite(self, name):
         """The method writes Cartesian mesh on the disc.
@@ -209,5 +210,5 @@
         
         f.description = "Cartesian_mesh"
-        f.nx, f.ny = self.nx, self.ny
+        f.meshtype, f.nx, f.ny = 'curvilinear', self.nx, self.ny
 
         def create_vars(dimname, info):
@@ -218,7 +219,10 @@
         create_vars("primal_cell", 
                     [("primal_deg","i4",self.primal_deg), 
+                     ("primal_i","i4",self.primal_i),
+                     ("primal_j","i4",self.primal_j),
                      ("Ai","f8",self.Aiv),
                      ("lon_i","f8",self.lon_i),
                      ("lat_i","f8",self.lat_i)] )
+
         create_vars("dual_cell", 
                     [("dual_deg","i4",self.dual_deg), 
@@ -293,5 +297,8 @@
     start_index=1 # indexing starts at 1 as in Fortran
     def __init__(self,gridfilename):
-        self.nc = cdf.Dataset(gridfilename, "r")
+        nc = cdf.Dataset(gridfilename, "r")
+        self.nc, self.meshtype = nc, nc.meshtype
+        if self.meshtype == 'curvilinear':
+            self.nx, self.ny = nc.nx, nc.ny
         print('#NC4: Opening file:',gridfilename)
     def get_pdim(self,comm,name): return parallel.PDim(self.nc.dimensions[name], comm)
@@ -310,6 +317,7 @@
     """ Helps class Unstructured_Mesh to read a MPAS mesh file."""
     start_index=1 # indexing starts at 1 as in Fortran
+    meshtype='unstructured'
     translate= {
-        'primal_num':'nCells', 'edge_num':'nEdges', 'dual_num':'nVertices',
+        'primal_cell':'nCells', 'edge':'nEdges', 'dual_cell':'nVertices',
         'primal_deg':'nEdgesOnCell', 'trisk_deg':'nEdgesOnEdge',
         'primal_edge':'edgesOnCell', 'dual_edge':'edgesOnVertex', 
@@ -452,5 +460,5 @@
         
 class Unstructured_PMesh(Abstract_Mesh): # Mesh data distributed across MPI processes
-    def __init__(self,comm, gridfile, meshtype='unstructured'):
+    def __init__(self,comm, gridfile):
         self.comm, self.gridfile, get_pvars = comm, gridfile, gridfile.get_pvars
         dim_primal, dim_edge, dim_dual = gridfile.get_pdims(comm,
@@ -463,4 +471,7 @@
         dual_deg, dual_vertex, dual_edge, Riv2, lon_v, lat_v, Av = get_pvars(
             dim_dual, 'dual_deg', 'dual_vertex', 'dual_edge', 'Riv2', 'lon_v', 'lat_v', 'Av')
+        if gridfile.meshtype == 'curvilinear':
+            self.primal_i, self.primal_j = get_pvars(dim_primal, 'primal_i','primal_j')
+
         # Let indices start at 0
         for x in (primal_vertex, primal_edge, dual_vertex,dual_edge,
@@ -478,5 +489,5 @@
         primal_owner = partition_from_stencil(edge_owner, primal_deg, primal_edge)
         primal_owner = parallel.LocPArray1D(dim_primal, primal_owner)
-        self.set_members(locals(), 'meshtype', 'dim_primal', 'dim_dual', 'dim_edge',
+        self.set_members(locals(), 'dim_primal', 'dim_dual', 'dim_edge',
                          'edge_owner', 'primal_owner', 'primal_deg', 'primal_vertex', 
                          'lon_v', 'lat_v', 'Av', 'lon_i', 'lat_i', 'Ai', 
@@ -530,6 +541,9 @@
         primal_num, dual_num, edge_num = map(len, (cells_C1, vertices_V1, edges_E2))
 
-        if pmesh.meshtype == 'curvilinear' :
-            self.primal_i, self.primal_j = pmesh.get(get_own_cells, 'primal_i', 'primal_j')
+        self.meshtype = pmesh.gridfile.meshtype
+        if self.meshtype == 'curvilinear' :
+            self.nx, self.ny = pmesh.gridfile.nx, pmesh.gridfile.ny
+            self.primal_i, self.primal_j = pmesh.get(get_all_cells, 'primal_i', 'primal_j')
+
 
         # construct local stencils from global stencils
@@ -566,4 +580,5 @@
         primal_own_all = comm.allgather(primal_own_num)
         displ = sum(primal_own_all[0:comm.Get_rank()]) # required by XIOS
+
         print 'local primal mesh :', primal_vertex.min(), primal_vertex.max(), primal_own_all, displ
         
@@ -742,107 +757,2 @@
     # example : edge_list = progressive_list(E0,E1,E2) with E0,E1,E2 increasing lists of cell edges
     return list(progressive_iter([], cell_lists))
-from dynamico.meshes import zeros
-import numpy as np
-import netCDF4 as cdf
-import argparse
-
-
-if False:
-    def __init__(self,nx,ny,llm,nqdyn,Lx,Ly):
-        dx,dy = Lx/nx, Ly/ny
-        self.dx, self.dy = dx,dy
-        self.nx, self.ny, self.llm, self.nqdyn = nx,ny,llm,nqdyn
-        self.field_z = self.field_mass
-        # 1D coordinate arrays
-        x=(np.arange(nx)-nx/2.)*Lx/nx
-        y=(np.arange(ny)-ny/2.)*Ly/ny
-        lev=np.arange(llm)
-        levp1=np.arange(llm+1)
-        self.x, self.y, self.lev, self.levp1 = x,y,lev,levp1
-        # 3D coordinate arrays
-        self.xx,self.yy,self.ll = np.meshgrid(x,y,lev, indexing='ij')
-        self.xxp1,self.yyp1,self.llp1 = np.meshgrid(x,y,levp1, indexing='ij')
-        # beware conventions for indexing
-        # Fortran order : llm,nx*ny,nqdyn  / indices start at 1
-        # Python order : nqdyn,ny,nx,llm   / indices start at 0
-        # indices below follow Fortran while x,y follow Python/C
-        index=lambda x,y : ((x+(nx*(y+2*ny)))%(nx*ny))+1
-        indexu=lambda x,y : 2*index(x,y)-1
-        indexv=lambda x,y : 2*index(x,y)
-        indices = lambda shape : np.zeros(shape,np.int32)
-
-        for x in range(nx):
-            for y in range(ny):
-                # NB : Fortran indices start at 1
-                # primal cells
-                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) ))
-                # dual cells
-                put(index(x,y)-1,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) ))
-                # edges :
-                # left and right are adjacent primal cells
-                # flux is positive when going from left to right
-                # up and down are adjacent dual cells (vertices)
-                # circulation is positive when going from down to up
-
-        Aiv=np.zeros(nx*ny)+dx*dy  
-        Ai=Av=np.zeros(nx*ny)+dx*dy
-
-        self.llm=llm
-        self.nqdyn=nqdyn
-        self.nx=nx
-        self.ny=ny
-        self.primal_deg=primal_deg
-        self.primal_edge=primal_edge
-        self.primal_ne=primal_ne
-        self.dual_deg=dual_deg
-        self.dual_edge=dual_edge
-        self.dual_ne=dual_ne
-        self.dual_vertex=dual_vertex
-        self.primal_vertex=primal_vertex
-        self.left=left
-        self.right=right
-        self.down=down
-        self.up=up
-        self.trisk_deg=trisk_deg
-        self.trisk=trisk
-        self.Aiv=Aiv
-        self.Ai=Ai
-        self.Av=Av
-        self.le=le
-        self.de=de
-        self.angle_e=angle_e
-        self.Riv2=Riv2
-        self.wee=wee
-        self.lon_i = lon_i
-        self.lon_v = lon_v
-        self.lon_e = lon_e
-        #self.lon_ev = indices(2*nx*ny)
-        self.lat_i = lat_i
-        self.lat_v = lat_v
-        self.lat_e = lat_e
-        #self.lat_ev = indices(2*nx*ny)
-      
-
-    def field_theta(self,n=1): return zeros((n,self.nqdyn,self.ny,self.nx,self.llm))
-    def field_mass(self,n=1): return zeros((n,self.ny,self.nx,self.llm))
-    def field_z(self,n=1): return zeros((n,self.ny,self.nx,self.llm))
-    def field_w(self,n=1): return zeros((n,self.ny,self.nx,self.llm+1))
-    def field_u(self,n=1):
-        if n==1:
-            return np.zeros((self.ny,2*self.nx,self.llm))
-        else:
-            return np.zeros((n,self.ny,2*self.nx,self.llm))
-    def field_ps(self,n=1): return zeros((n,self.ny,self.nx))
-    def ucomp(self,u): return u[:,range(0,2*self.nx,2),:]
-    def set_ucomp(self,uv,u): uv[:,range(0,2*self.nx,2),:]=u
-    def vcomp(self,u): return u[:,range(1,2*self.nx,2),:]
-
-
Index: codes/icosagcm/devel/Python/dynamico/xios.py
===================================================================
--- codes/icosagcm/devel/Python/dynamico/xios.py	(revision 758)
+++ codes/icosagcm/devel/Python/dynamico/xios.py	(revision 759)
@@ -4,4 +4,15 @@
 from dynamico.meshes import radian
 from cxios import finalize
+
+#-----------------------------------------------------------------------------
+
+def setup_curvilinear(cat,id,nx,ny,displ,index_own_i,index_own_j,lon,lat):
+    mesh = cxios.Handle(cat,id)
+    print 'setup_curvilinear : index_own_i.shape', index_own_i.shape
+    print 'setup_curvilinear : index_own_j.shape', index_own_j.shape
+    mesh.set_attr(type='curvilinear',
+                  ni_glo=nx, i_index=index_own_i,
+                  nj_glo=ny, j_index=index_own_j)
+    mesh.set_attr(lonvalue_1d=lon.flatten(), latvalue_1d=lat.flatten())
 
 def cf_boundaries(degree,points,lon,lat):
@@ -19,15 +30,20 @@
     return bnd_lon, bnd_lat
 
-def setup_domain(cat,id, degree,vertex,lon,lat,lon_v,lat_v, ncell_glo, index_glo, displ): # cat in 'domain','domaingroup'
+def setup_unstructured(cat,id, degree,vertex,lon,lat,lon_v,lat_v, ncell_glo, index_own_glo, displ): # cat in 'domain','domaingroup'
+    # ncell_glo (int) : global number of cells
+    # index_own_glo : global indices of cells own by this MPI process
+    # displ : min of index_own across all MPI processes
     bnd_lon, bnd_lat = cf_boundaries(degree, vertex, lon_v, lat_v)
     ncell, nvertex = bnd_lon.shape
     mesh = cxios.Handle(cat,id)
     mesh.set_attr(type='unstructured')
-    mesh.set_attr(ni_glo=ncell_glo, ni=ncell, ibegin=displ, i_index=index_glo)
+    mesh.set_attr(ni_glo=ncell_glo, ni=ncell, ibegin=displ, i_index=index_own_glo)
     mesh.set_attr(lonvalue_1d=lon, latvalue_1d=lat)
     mesh.set_attr(nvertex=nvertex,bounds_lon_1d=bnd_lon, bounds_lat_1d=bnd_lat)
 
+#-----------------------------------------------------------------------------
+
 class XIOS_Context:
-    def __init__(self, pmesh, mesh,nqtot, step):
+    def init_llm(self, mesh, nqtot):
         self.mesh=mesh
         ker.dynamico_setup_xios()
@@ -38,11 +54,5 @@
         levp1.set_attr(n_glo=llm+1, value=np.arange(llm+1))
         nq.set_attr(n_glo=nqtot, value=np.arange(nqtot))
-        # primal mesh
-        lon_i, lat_i, lon_v, lat_v = [x*radian for x in mesh.lon_i, mesh.lat_i, mesh.lon_v, mesh.lat_v]
-        mesh_i=setup_domain('domaingroup','i', mesh.primal_own_deg, mesh.primal_vertex,
-                            lon_i[mesh.primal_own_loc], lat_i[mesh.primal_own_loc],
-                            lon_v, lat_v,
-                            pmesh.dim_primal.n, mesh.primal_own_glo, mesh.displ)
-        # calendar
+    def init_calendar(self, step):
         calendar=cxios.Handle('calendar_wrapper')
         dtime = cxios.Duration(second=step)
@@ -53,4 +63,42 @@
         freq_offset = cxios.Duration(second=0)
         std.set_attr(freq_offset=freq_offset)
+    def update_calendar(self, i):
+        ker.dynamico_xios_update_calendar(i)
+    def finalize(self):
+        cxios.context_finalize()
+
+class XIOS_Context_Curvilinear(XIOS_Context):
+    def __init__(self, mesh, nqtot, step):
+        self.init_llm(mesh, nqtot)
+        # primal mesh
+#        lat_i,lon_i = np.meshgrid(mesh.y,mesh.x, indexing='ij')
+#        ny,nx=lat_i.shape
+        own = mesh.primal_own_loc
+        lon_i, lat_i = [ x[own]*radian for x in mesh.lon_i, mesh.lat_i ]
+        primal_i, primal_j = [ x[own] for x in mesh.primal_i, mesh.primal_j ]
+        print 'setup_curvilinear : ', lon_i.size, lat_i.size, primal_i.size, primal_j.size
+        setup_curvilinear('domaingroup','i', mesh.nx, mesh.ny, mesh.displ,
+                          primal_i, primal_j, lon_i, lat_i)
+        # calendar
+        self.init_calendar(step)
+        print 'cxios.context_close_definition()'
+        cxios.context_close_definition()
+    def send_field_primal(self,name, data):
+        if data.ndim==2:
+            data = data.transpose() # XIOS expects contiguous horizontal slices
+        data = np.ascontiguousarray(data)
+        cxios.send_field(name, data)
+
+class XIOS_Context_Unstructured(XIOS_Context):
+    def __init__(self, pmesh, mesh,nqtot, step):
+        self.init_llm(mesh, nqtot)
+        # primal mesh
+        lon_i, lat_i, lon_v, lat_v = [x*radian for x in mesh.lon_i, mesh.lat_i, mesh.lon_v, mesh.lat_v]
+        setup_unstructured('domaingroup','i', mesh.primal_own_deg, mesh.primal_vertex,
+                           lon_i[mesh.primal_own_loc], lat_i[mesh.primal_own_loc],
+                           lon_v, lat_v,
+                           pmesh.dim_primal.n, mesh.primal_own_glo, mesh.displ)
+        # calendar
+        self.init_calendar(step)
         print 'cxios.context_close_definition()'
         cxios.context_close_definition()
@@ -64,6 +112,2 @@
         data = np.ascontiguousarray(data)
         cxios.send_field(name, data)
-    def update_calendar(self, i):
-        ker.dynamico_xios_update_calendar(i)
-    def finalize(self):
-        cxios.context_finalize()
Index: codes/icosagcm/devel/Python/src/cxios.pyx
===================================================================
--- codes/icosagcm/devel/Python/src/cxios.pyx	(revision 758)
+++ codes/icosagcm/devel/Python/src/cxios.pyx	(revision 759)
@@ -67,5 +67,5 @@
             fun = cxios.vardict[prefix+key]
             extra=value
-            if type(value) in (int,c_int):
+            if type(value) in (int,c_int,np.int64,np.int32):
                 fun(handle,c_int(value))
             elif type(value) is np.ndarray:
@@ -102,6 +102,7 @@
                        ['update_calendar',c_int] ])
 
-    import_set_attr(['domain','domaingroup'], [c_int,'ni_glo','ibegin','ni','nvertex','data_dim'], [str,'type'],
-                        [np.ndarray, 'i_index','lonvalue_1d','latvalue_1d','bounds_lat_1d','bounds_lon_1d'])
+    import_set_attr(['domain','domaingroup'], [c_int,'ni_glo','nj_glo','ibegin','jbegin','ni','nvertex','data_dim'], [str,'type'],
+                        [np.ndarray, 'i_index','j_index',
+			'lonvalue_1d','latvalue_1d','bounds_lat_1d','bounds_lon_1d'])
     import_set_attr(['field','fieldgroup'], [Duration,'freq_op','freq_offset'])
     import_set_attr(['axis'], [c_int,'n_glo'], [np.ndarray,'value'])
Index: codes/icosagcm/devel/Python/test/py/RSW_2D_mesh.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/RSW_2D_mesh.py	(revision 758)
+++ codes/icosagcm/devel/Python/test/py/RSW_2D_mesh.py	(revision 759)
@@ -15,5 +15,5 @@
 dx,dy=Lx/nx,Ly/ny
 
-filename, llm, nqdyn, g, f, radius = 'cart_128.nc', 1, 1, 1., 1., None
+filename, llm, nqdyn, g, f, radius = 'cart_%03d_%03d.nc'%(nx,ny), 1, 1, 1., 1., None
 unst.setvar('g',g)
 
@@ -27,6 +27,7 @@
 caldyn = unst.Caldyn_RSW(mesh)
 
-xx = (mesh.lat_i-nx/2.)*dx
-yy = (mesh.lon_i-ny/2.)*dy
+#xx = (mesh.lat_i-nx/2.)*dx
+#yy = (mesh.lon_i-ny/2.)*dy
+xx,yy = mesh.lat_i, mesh.lon_i
 
 x1,x2,yy = xx-1., xx+1., yy
Index: codes/icosagcm/devel/Python/test/py/test_xios.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/test_xios.py	(revision 758)
+++ codes/icosagcm/devel/Python/test/py/test_xios.py	(revision 759)
@@ -28,5 +28,5 @@
 #--------------------------------- write some data ----------------------------------------
 
-context=xios.XIOS_Context(pmesh,mesh,nqtot, 3600)
+context=xios.XIOS_Context_Unstructured(pmesh,mesh,nqtot, 3600)
 
 lat_i = radian*mesh.lat_i
Index: codes/icosagcm/devel/Python/test/py/test_xios_cartesian.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/test_xios_cartesian.py	(revision 759)
+++ codes/icosagcm/devel/Python/test/py/test_xios_cartesian.py	(revision 759)
@@ -0,0 +1,57 @@
+from dynamico import cxios
+from dynamico import xios
+from dynamico import unstructured as unst
+from dynamico import meshes
+ker=unst.ker
+
+from mpi4py import MPI
+comm = MPI.COMM_WORLD
+mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size()
+print '%d/%d starting'%(mpi_rank,mpi_size)
+
+import numpy as np
+
+nqdyn,nqtot,llm = 1,1,4
+
+if False:
+    nx,ny,llm = 64,64
+    Lx,Ly,f = 1.,1.,1.
+    mesh = meshes.Cartesian_Mesh(nx,ny,llm,nqdyn,Lx,Ly,f)
+else:
+    filename, g, f, radius = 'cart_128_128.nc', 1., 1., None
+    filename = 'cart_004_004.nc'
+    unst.setvar('g',g)
+    print 'Reading Cartesian mesh ...'
+    def coriolis(lon,lat):
+        return f+0.*lon
+    meshfile = meshes.DYNAMICO_Format(filename)
+    nx,ny=meshfile.nx,meshfile.ny
+
+    pmesh = meshes.Unstructured_PMesh(comm,meshfile)
+    mesh = meshes.Local_Mesh(pmesh, llm, nqdyn, radius, coriolis)
+
+#--------------------------------- write some data ----------------------------------------
+
+context=xios.XIOS_Context_Curvilinear(mesh,nqtot, 3600)
+#lon_i,lat_i = mesh.xx[:,:,0].flatten(), mesh.yy[:,:,0].flatten()
+lon_i,lat_i = mesh.lon_i, mesh.lat_i
+lat_ik, junk  = np.meshgrid(lat_i, np.arange(llm), indexing='ij')
+
+no_error=True
+try:
+    for i in range(3):
+        context.update_calendar(i)
+        print 'send_field', i, lat_ik.shape
+        context.send_field_primal('ps', lat_i)
+        context.send_field_primal('theta', lat_ik)
+except:
+    print 'An error has occured, trying to close XIOS properly'
+    no_error=False
+
+print 'xios.context_finalize()'
+context.finalize()
+print 'xios.finalize()'
+xios.finalize()
+print 'Done'
+
+if not no_error : raise
Index: codes/icosagcm/devel/Python/test/python.sh
===================================================================
--- codes/icosagcm/devel/Python/test/python.sh	(revision 758)
+++ codes/icosagcm/devel/Python/test/python.sh	(revision 759)
@@ -8,4 +8,10 @@
 {
     LD_PRELOAD=$PYTHON_PRELOAD python -u $*
+}
+
+function cmd_gdb()
+{
+    set -x
+    LD_PRELOAD=$PYTHON_PRELOAD gdb --args python -u $*
 }
 
Index: codes/icosagcm/devel/Python/test/xml/filedef_simple.xml
===================================================================
--- codes/icosagcm/devel/Python/test/xml/filedef_simple.xml	(revision 758)
+++ codes/icosagcm/devel/Python/test/xml/filedef_simple.xml	(revision 759)
@@ -44,6 +44,5 @@
     <variable name="time_frequency" type="string" > 24h </variable>    
 
-  </file>
--->
+  </file> -->
 
 </file_definition>
Index: codes/icosagcm/devel/Python/test/xml/link.sh
===================================================================
--- codes/icosagcm/devel/Python/test/xml/link.sh	(revision 758)
+++ codes/icosagcm/devel/Python/test/xml/link.sh	(revision 759)
@@ -7,5 +7,5 @@
 	ln -s xml/${name}_$1.xml ${name}.xml
     done
-    ln -s iodef.xml
+    ln -s xml/iodef.xml
 }
 
Index: codes/icosagcm/devel/src/unstructured/data_unstructured.F90
===================================================================
--- codes/icosagcm/devel/src/unstructured/data_unstructured.F90	(revision 758)
+++ codes/icosagcm/devel/src/unstructured/data_unstructured.F90	(revision 759)
@@ -164,9 +164,9 @@
     dual_ne(:,:) = dual_ne_(:,:)
     dual_vertex(:,:) = dual_vertex_(:,:)
-    IF(MINVAL(dual_deg)<3) THEN
-       STOP 'At least one dual cell has less than 3 vertices'
-    END IF
-    IF(MINVAL(primal_deg)<3) THEN
-       STOP 'At least one primal cell has less than 3 vertices'
+    IF(MINVAL(dual_deg)<2) THEN
+       STOP 'At least one dual cell has less than 2 vertices'
+    END IF
+    IF(MINVAL(primal_deg)<2) THEN
+       STOP 'At least one primal cell has less than 2 vertices'
     END IF
     left(:)=left_(:)
