Index: codes/icosagcm/devel/Python/dev/__init__.py
===================================================================
--- codes/icosagcm/devel/Python/dev/__init__.py	(revision 977)
+++ codes/icosagcm/devel/Python/dev/__init__.py	(revision 977)
@@ -0,0 +1,6 @@
+# import openmp before NetCDF4 numpy because it uses numpy hence MKL and messes up OMP_NUM_PROCS
+import sys
+if 'numpy' in sys.modules:
+    raise ImportError('Module dynamico must be imported before NetCDF/numpy because numpy may mess up OpenMP')
+import libs
+libs.libicosa.dynamico_show_openmp()
Index: codes/icosagcm/devel/Python/dev/cxios.pyx
===================================================================
--- codes/icosagcm/devel/Python/dev/cxios.pyx	(revision 977)
+++ codes/icosagcm/devel/Python/dev/cxios.pyx	(revision 977)
@@ -0,0 +1,114 @@
+from dynamico.dev.libs import libxios
+from dynamico.dev import wrap
+import numpy as np
+
+from ctypes import c_int, c_double, c_bool, c_void_p, c_char_p, byref, POINTER, Structure
+c_void_pp=POINTER(c_void_p) # used in prototype
+c_void_p_byref=type(byref(c_void_p())) # used in py2c because byref creates an object of this type, not c_void_pp
+
+from dynamico import getargs
+log_master, log_world = getargs.getLogger(__name__)
+INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
+INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
+
+# ------------- Direct cython interfaces ----------------- #
+
+cdef extern from "functions.h":
+     cdef void cxios_context_close_definition()
+     cdef void cxios_context_finalize()
+     cdef void cxios_finalize()
+     cdef void cxios_write_data_k81(char*, int, double*, int)
+     cdef void cxios_write_data_k82(char*, int, double*, int,int)
+     cdef void cxios_write_data_k83(char*, int, double*, int,int,int)
+
+def context_close_definition(): cxios_context_close_definition()
+def context_finalize(): cxios_context_finalize()
+def finalize(): cxios_finalize()
+
+cdef send_field1(char* id, int idlen, double[:] data):     cxios_write_data_k81(id, idlen, &data[0],     data.shape[0])
+cdef send_field2(char* id, int idlen, double[:,:] data):   cxios_write_data_k82(id, idlen, &data[0,0],   data.shape[0],data.shape[1])
+cdef send_field3(char* id, int idlen, double[:,:,:] data): cxios_write_data_k83(id, idlen, &data[0,0,0], data.shape[0],data.shape[1],data.shape[2])
+
+def send_field(bytes id, data):
+    cdef char* idptr = id
+    cdef int idlen = len(id)
+#    print 'cxios.write_data', id, data.shape
+    ndim=data.ndim
+    if ndim==1: send_field1(idptr,idlen,data)
+    if ndim==2: send_field2(idptr,idlen,data)
+    if ndim==3: send_field3(idptr,idlen,data)
+
+# -------------------------------------------------------- #
+
+lib=wrap.Struct()
+cxios = wrap.SharedLib(vars(lib), libxios, prefix_so='cxios_')
+
+class Duration(Structure):
+    _fields_ = [(name, c_double) for name in ('year', 'month', 'day', 'hour', 'minute', 'second', 'timestep') ]
+class Date(Structure):
+    _fields_ = [(name, c_int) for name in ('year', 'month', 'day', 'hour', 'minute', 'second') ]
+
+for x in c_void_p_byref, Duration, Date : wrap.py2c[x]=wrap.noop # accept arguments of type (void*)&, Duration, ...
+
+def SHAPE(data): return np.asarray(list(reversed(data.shape)),c_int) # reverse shape to follow Fortran convention
+
+class Handle:
+    def __init__(self, cat, id=None):
+        self.handle = c_void_p()
+        if id is None:
+            id='default'
+            create_fun = cxios.vardict['get_current_'+cat]
+            create_fun(byref(self.handle))
+        else:
+            create_fun = cxios.vardict[cat+'_handle_create']
+            create_fun(byref(self.handle), id, c_int(len(id)))
+        self.cat, self.prefix_set, self.id = cat, 'set_'+cat+'_', id
+#        print 'Handle', cat, self.id, self.handle
+    def set_attr(self, **kwargs):
+        prefix, handle = self.prefix_set, self.handle
+        for key, value in kwargs.iteritems():
+            fun = cxios.vardict[prefix+key]
+            extra=value
+            if type(value) in (int,c_int,np.int64,np.int32):
+                fun(handle,c_int(value))
+            elif type(value) is np.ndarray:
+                extra=SHAPE(value)
+                fun(handle,value,extra)
+            elif type(value) is str:
+                fun(handle,value,c_int(len(value)))
+            elif type(value) in (Duration, Date):
+                fun(handle,value)
+            else:
+                raise(TypeError)
+            DEBUG('%s.%s %s %s' % (self.id, prefix+key, type(value), extra) )
+    def update_timestep(self):
+        lib.update_calendar_timestep(self.handle)
+
+def import_set_attr(cats, *args): #attr_int, attr_array, attr_str):
+    for arg in args:
+        argtype = arg.pop(0)
+        if argtype is str:
+            argtype = [c_void_p, c_char_p, c_int] #char *str, int strlen
+        elif argtype is np.ndarray:
+            argtype = [c_void_p, c_void_p, c_void_p] # void *handle, double *data, int* SHAPE(data)
+        else:
+            argtype = [c_void_p, argtype] #void *handle, argtype data
+        for cat in cats:
+            cxios.import_funs( [arg + argtype], prefix='set_'+cat+'_')
+
+def import_cxios_functions():
+    cats=['axis','domain','domaingroup','field','fieldgroup'] # XIOS categories
+    cxios.import_funs([
+                       [cat+'_handle_create' for cat in cats]+[c_void_pp,c_char_p,c_int],
+                       ['get_current_'+'calendar_wrapper',c_void_pp],
+                       ['update_calendar_timestep',c_void_p],
+                       ['update_calendar',c_int] ])
+
+    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'])
+    import_set_attr(['calendar_wrapper'], [Duration,'timestep'])
+
+import_cxios_functions()
Index: codes/icosagcm/devel/Python/dev/functions.h
===================================================================
--- codes/icosagcm/devel/Python/dev/functions.h	(revision 977)
+++ codes/icosagcm/devel/Python/dev/functions.h	(revision 977)
@@ -0,0 +1,37 @@
+/* KERNELS */
+
+#ifdef CPP_MIXED_PREC
+typedef float number;
+#else
+typedef double number;
+#endif
+
+enum {max_nb_stage=5};
+extern number tauj[max_nb_stage];
+extern number cslj[max_nb_stage][max_nb_stage], cflj[max_nb_stage][max_nb_stage];
+
+void dynamico_init_params(void);
+
+void dynamico_ARK_step(int nstep,
+		       number *mass_col, number *rhodz, number *theta_rhodz, 
+		       number *u, number *geopot, number *w,
+		       number *theta, number *ps, number *pk, number *hflux, number *qv,
+		       number *dmass_col, number *drhodz, number *dtheta_rhodz,
+		       number *du_fast, number *du_slow,
+		       number *dPhi_fast, number *dPhi_slow, 
+		       number *dW_fast, number *dW_slow);
+
+void dynamico_remap(number *rhodz, number *theta_rhodz, number *u);
+
+/* KERNELS -> XIOS */
+void dynamico_setup_xios(void);
+void dynamico_xios_set_timestep(double);
+void dynamico_xios_update_calendar(int);
+/* CXIOS */
+void cxios_context_close_definition(void);
+void cxios_context_finalize(void);
+void cxios_finalize(void);
+
+void cxios_write_data_k81(char*, int, double*, int);
+void cxios_write_data_k82(char*, int, double*, int, int);
+void cxios_write_data_k83(char*, int, double*, int, int, int);
Index: codes/icosagcm/devel/Python/dev/libs.py
===================================================================
--- codes/icosagcm/devel/Python/dev/libs.py	(revision 977)
+++ codes/icosagcm/devel/Python/dev/libs.py	(revision 977)
@@ -0,0 +1,2 @@
+from ctypes import cdll
+libxios ,libicosa = [ cdll.LoadLibrary(name) for name in 'libxios.so', 'libicosa.so' ]
Index: codes/icosagcm/devel/Python/dev/meshes.py
===================================================================
--- codes/icosagcm/devel/Python/dev/meshes.py	(revision 977)
+++ codes/icosagcm/devel/Python/dev/meshes.py	(revision 977)
@@ -0,0 +1,24 @@
+from __future__ import print_function
+
+from dynamico import meshes
+from dynamico.meshes import Unstructured_PMesh, compute_hybrid_coefs
+from dynamico.meshes import DYNAMICO_Format, MPAS_Format
+from dynamico.dev import unstructured as unst
+from dynamico.dev import numba, parallel
+from dynamico.dev.numba import int32, int64, float64
+
+class MeshData(numba.NumbaData):
+    def to_dynamico(self):
+        max_primal_deg, max_dual_deg, max_trisk_deg = [x.shape[1] for x in self.primal_edge, self.dual_edge, self.trisk]
+        unst.init_mesh(self.llm, self.nqdyn, self.edge_num, self.primal_num, self.dual_num,
+                  max_trisk_deg, max_primal_deg, max_dual_deg,
+                  self.primal_deg, self.primal_edge, self.primal_ne,
+                  self.dual_deg, self.dual_edge, self.dual_ne, self.dual_vertex,
+                  self.left, self.right, self.down, self.up, self.trisk_deg, self.trisk,
+                  self.Ai, self.Av, self.fv, self.le_de, self.Riv2, self.wee)
+
+# note inheritance order : inherit first from DevMesh to override empty to_dynamico() from AbstractMesh
+class Cartesian_Mesh(MeshData, meshes.Cartesian_Mesh): pass
+class Unstructured_Mesh(MeshData, meshes.Unstructured_Mesh): pass
+class Local_Mesh(MeshData, meshes.Local_Mesh):
+    Halo_Xchange = parallel.Halo_Xchange # dynamico.dev.parallel instead pf dynamico.parallel
Index: codes/icosagcm/devel/Python/dev/numba.py
===================================================================
--- codes/icosagcm/devel/Python/dev/numba.py	(revision 977)
+++ codes/icosagcm/devel/Python/dev/numba.py	(revision 977)
@@ -0,0 +1,47 @@
+from __future__ import absolute_import     
+from __future__ import print_function
+
+import numpy as np
+import numba
+from numba import int32, int64, float64
+
+class NumbaData(object):
+    """A base class to extract data from derived class instances and use it as argument to @jit functions. Derived classes must set the 'signature' attribute."""
+    def data(self):
+        """Returns a jitclass instance containing attributes copied from self, using self.signature which is of the form type,names,type,names ... where names is a string 'attr1 attr2 attr3' containing space-separated names of attributes of self. Those attributes are declared to numba with the type preceding them. The result of data() can be used as argument to a @jit function."""
+        cls = self.__class__.__name__
+        spec = []
+        for name in dir(self):
+            attr = getattr(self,name)
+            if not name.startswith('__') and not callable(attr):
+                tp = None
+                if isinstance(attr, int): tp=int64
+                if isinstance(attr, float): tp=float64
+                if isinstance(attr, np.ndarray):
+                    dtype=attr.dtype
+                    if dtype == np.int32 : dtype=int32
+                    elif dtype == np.float64 : dtype=float64
+                    else: dtype=None
+                    if dtype is None: print('Unknown dtype ', attr.dtype)
+                        
+                    if   len(attr.shape)==1 : tp = dtype[:]
+                    elif len(attr.shape)==2 : tp = dtype[:,:]
+                    elif len(attr.shape)==3 : tp = dtype[:,:,:]
+                    else: print('%s.%s is a numpy array with unsupported rank >3'%(cls,name))
+                        
+                if tp is None:
+                    print('Type of attribute %s.%s is not recognized'%(cls,name), type(attr))
+                else:
+                    spec.append( (name, tp) )
+        
+        @numba.jitclass(spec)
+        class JitClass(object):
+            def __init__(self): pass
+
+        data=JitClass()
+        for name,thetype in spec: 
+            print( 'Making %s.%s available @jit functions : '%(self.__class__.__name__, name), type(getattr(self,name)))
+            setattr(data, name, getattr(self,name))
+        return data
+
+jit=numba.jit(nopython=True, nogil=True, error_model='numpy', fastmath=True)
Index: codes/icosagcm/devel/Python/dev/parallel.py
===================================================================
--- codes/icosagcm/devel/Python/dev/parallel.py	(revision 977)
+++ codes/icosagcm/devel/Python/dev/parallel.py	(revision 977)
@@ -0,0 +1,9 @@
+from __future__ import absolute_import     
+from __future__ import print_function
+
+from dynamico import parallel
+from dynamico.dev.unstructured import ker
+
+class Halo_Xchange(parallel.Halo_Xchange):
+    def dynamico_init_transfer(self, *args):
+        ker.dynamico_init_transfer(*args)
Index: codes/icosagcm/devel/Python/dev/precision.py
===================================================================
--- codes/icosagcm/devel/Python/dev/precision.py	(revision 977)
+++ codes/icosagcm/devel/Python/dev/precision.py	(revision 977)
@@ -0,0 +1,10 @@
+"Helper function to manage mixed precision"
+import numpy as np
+from unstructured import np_num
+
+zero=np_num(0.)
+
+# cast array to float/double depending on CPP_MIXED_PREC
+def asnum(data): return [ x.astype(np_num) for x in data] if type(data) is list else data.astype(np_num)
+
+def zeros(dims): return np.zeros([n for n in dims if n>1], dtype=np_num)
Index: codes/icosagcm/devel/Python/dev/setup.py
===================================================================
--- codes/icosagcm/devel/Python/dev/setup.py	(revision 977)
+++ codes/icosagcm/devel/Python/dev/setup.py	(revision 977)
@@ -0,0 +1,19 @@
+from distutils.core import setup
+from distutils.extension import Extension
+from Cython.Build import cythonize
+
+from ctypes import cdll, c_bool
+libicosa = cdll.LoadLibrary('libicosa.so')
+mixed_prec=c_bool.in_dll(libicosa, 'mixed_precision').value
+
+if mixed_prec:
+    extra_flags=["-DCPP_MIXED_PREC"]
+else:
+    extra_flags=None
+
+setup(  name = "DYNAMICO kernels for unstructured mesh",
+        ext_modules = cythonize([
+            Extension("cxios",["cxios.pyx"],libraries=["xios"]),
+            Extension("unstructured",["unstructured.pyx"],libraries=["icosa","xios"], extra_compile_args=extra_flags)
+            ], compile_time_env={'CPP_MIXED_PREC': mixed_prec} )
+        )
Index: codes/icosagcm/devel/Python/dev/unstructured.pyx
===================================================================
--- codes/icosagcm/devel/Python/dev/unstructured.pyx	(revision 977)
+++ codes/icosagcm/devel/Python/dev/unstructured.pyx	(revision 977)
@@ -0,0 +1,339 @@
+import time
+import math
+from dynamico.dev import wrap
+from dynamico.dev.libs import libicosa
+from dynamico.util import list_stencil
+
+from ctypes import c_void_p, c_int, c_double, c_float, c_bool
+from numpy cimport ndarray
+cimport numpy as np
+import numpy as np
+
+from dynamico import getargs
+log_master, log_world = getargs.getLogger(__name__)
+INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
+INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
+
+#--------------   choose precision of kernel computations  ------------#
+
+# the compile-time constant CPP_MIXED_PREC is set by setup.py
+# based on the value of mixed_prec defined in module data_unstructured
+ 
+IF CPP_MIXED_PREC:
+   c_num=c_float
+   ctypedef float num
+   np_num=np.float32
+ELSE:
+   c_num=c_double
+   ctypedef double num
+   np_num=np.float64
+
+ctypedef num *num_ptr
+
+#------------- direct Cython interface to DYNAMICO routines -------------#
+
+
+cdef enum: max_nb_stage=5
+cdef extern :
+    cdef num tauj[max_nb_stage]
+    cdef num cslj[max_nb_stage][max_nb_stage]
+    cdef num cflj[max_nb_stage][max_nb_stage]
+    cdef int nb_stage[1]
+
+cdef extern from "functions.h":
+     cdef void dynamico_ARK_step(int nstep,
+                                 num *mass_col, num *rhodz, num *theta_rhodz, 
+                                 num *u, num *geopot, num *w,
+                                 num *theta, num *ps, num *pk, num *hflux, num *qv,
+                                 num *dmass_col, num *drhodz, num *dtheta_rhodz,
+                                 num *du_fast, num *du_slow,
+                                 num *dPhi_fast, num *dPhi_slow, 
+                                 num *dW_fast, num *dW_slow)
+     cdef void dynamico_remap(num *rhodz, num *theta_rhodz, num *u)
+     cdef void dynamico_init_params()
+     cpdef void dynamico_setup_xios()
+     cpdef void dynamico_xios_set_timestep(double)
+     cpdef void dynamico_xios_update_calendar(int)
+
+#------------- import and wrap DYNAMICO routines -------------#
+
+ker=wrap.Struct() # store imported fun X as funs.X
+
+check_args = False # use True instead of False for debugging, probably with some overhead
+
+try:    
+   kernels = wrap.SharedLib(vars(ker), libicosa, check_args=check_args) 
+   setvar, setvars, getvar, getvars = kernels.setvar, kernels.setvars, kernels.getvar, kernels.getvars
+except OSError:
+   ERROR("Unable to load shared library 'libicosa.so' !")
+   raise
+
+# providing a full prototype enables type-checking when calling
+# if a number n is present in the prototype, the previous type is repeated n times
+kernels.import_funs([
+                     ['dynamico_setup_xios',None],
+                     ['dynamico_print_trace',None],
+                     ['dynamico_xios_set_timestep',c_double],
+                     ['dynamico_xios_update_calendar',c_int],
+                     ['dynamico_init_mesh',c_void_p,13],
+                     ['dynamico_init_metric', c_void_p,6],
+                     ['dynamico_init_hybrid', c_void_p,3],
+                     ['dynamico_caldyn_unstructured', c_num, c_void_p,20],
+                     ['dynamico_scalar_laplacian', c_void_p,2],
+                     ['dynamico_curl_laplacian', c_void_p,2],
+                     ['dynamico_init_transfer', c_int, c_int,2,c_void_p,3, c_int,2,c_void_p,3],
+                     ['dynamico_update_halo', c_int,3,c_void_p] ])
+
+# set/get global variables
+eta_mass,eta_lag=(1,2)
+thermo_theta,thermo_entropy,thermo_moist,thermo_boussinesq=(1,2,3,4)
+
+kernels.addvars(
+    c_bool,'hydrostatic','debug_hevi_solver', 'is_mpi_master',
+    c_int,'llm','nqdyn','primal_num','max_primal_deg',
+    'dual_num','max_dual_deg','edge_num','max_trisk_deg',
+    'caldyn_thermo','caldyn_eta','nb_threads',
+    c_double,'elapsed',
+    c_num, 'g', 'ptop', 'cpp', 'cppv',
+    'Rd', 'Rv', 'preff', 'Treff', 'pbot', 'rho_bot', 'Phi_bot')
+
+elapsed=0.
+
+#------------------------ Extension type performing a full ARK time step ----------------------
+
+cdef num_ptr ptr1(num[:] data) except *: return &data[0]
+cdef num_ptr ptr2(num[:,:] data) except *: return &data[0,0]
+cdef num_ptr ptr3(num[:,:,:] data) except *: return &data[0,0,0]
+cdef num_ptr ptr4(num[:,:,:,:] data) except *: return &data[0,0,0,0]
+cdef num_ptr ptr(data) except * :
+    n=data.ndim
+    if n==1 : return ptr1(data)
+    if n==2 : return ptr2(data)
+    if n==3 : return ptr3(data)
+    if n==4 : return ptr4(data)
+    if n>4: raise IndexError
+        
+cdef alloc(num_ptr *p, allocator, n=1):
+    data=allocator(n)
+    p[0]=ptr(data)
+    return data
+
+cdef class Caldyn_step:
+    # number of time steps to do at each invocation of advance()
+    cdef int nstep
+    # pointer to allocated arrays
+    cdef num_ptr p_mass, p_theta_rhodz, p_u, p_geopot, p_W                   # prognostic
+    cdef num_ptr p_mass_col, p_dmass_col, p_ps, p_theta, p_pk, p_hflux, p_qv # diagnostic
+    cdef num_ptr p_drhodz, p_dtheta_rhodz, p_du_fast, p_du_slow                # tendencies
+    cdef num_ptr p_dPhi_fast, p_dPhi_slow, p_dW_fast, p_dW_slow                # tendencies
+    # allocated arrays, must remain referenced or segfault
+    cdef readonly ndarray mass, theta_rhodz, u, geopot, W
+    cdef readonly ndarray mass_col, dmass_col, ps, theta, pk, hflux, qv
+    cdef readonly ndarray drhodz, dtheta_rhodz, du_fast, du_slow
+    cdef readonly ndarray dPhi_fast, dPhi_slow, dW_fast, dW_slow
+
+    def __init__(self,mesh,time_scheme, nstep):
+        self.nstep=nstep
+        #        self.mesh=mesh
+        fps, ftheta, fmass = mesh.field_ps, mesh.field_theta, mesh.field_mass
+        fw, fu, fz = mesh.field_w, mesh.field_u, mesh.field_z
+        # collect coefficients of time scheme
+        cdef double[:]   tauj_ = time_scheme.tauj
+        cdef double[:,:] cslj_ = time_scheme.csjl
+        cdef double[:,:] cflj_ = time_scheme.cfjl
+        ns = time_scheme.nstage
+        nb_stage[0]=ns
+
+        cdef int i,j
+        for i in range(ns):
+            tauj[i]=tauj_[i]
+            for j in range(ns):
+                cslj[i][j]=cslj_[i,j]
+                cflj[i][j]=cflj_[i,j]
+        # allocate arrays, store pointers to avoid overhead when calling dynamico
+        #    prognostic/diagnostic
+        self.ps                       = alloc(&self.p_ps, fps) 
+        self.mass_col, self.dmass_col = alloc(&self.p_mass_col, fps), alloc(&self.p_dmass_col, fps,ns), 
+        self.mass, self.theta_rhodz   = alloc(&self.p_mass, fmass), alloc(&self.p_theta_rhodz, fmass),
+        self.theta, self.pk           = alloc(&self.p_theta, fmass), alloc(&self.p_pk, fmass), 
+        self.geopot, self.W           = alloc(&self.p_geopot, fw), alloc(&self.p_W, fw),
+        self.hflux, self.u            = alloc(&self.p_hflux, fu), alloc(&self.p_u, fu) 
+        self.qv                       = alloc(&self.p_qv,fz)
+        #    tendencies
+        self.drhodz, self.dtheta_rhodz  = alloc(&self.p_drhodz,fmass,ns), alloc(&self.p_dtheta_rhodz,fmass,ns) 
+        self.du_fast, self.du_slow      = alloc(&self.p_du_fast,fu,ns), alloc(&self.p_du_slow,fu,ns)
+        self.dPhi_fast, self.dPhi_slow  = alloc(&self.p_dPhi_fast,fw,ns), alloc(&self.p_dPhi_slow,fw,ns)
+        self.dW_fast, self.dW_slow      = alloc(&self.p_dW_fast,fw,ns), alloc(&self.p_dW_slow,fw,ns)
+    def next(self):
+        #        global elapsed
+        # time1=time.time()
+        dynamico_ARK_step(self.nstep,
+                          self.p_mass_col, self.p_mass, self.p_theta_rhodz,
+                          self.p_u, self.p_geopot, self.p_W,
+                          self.p_theta, self.p_ps, self.p_pk, self.p_hflux, self.p_qv,
+                          self.p_dmass_col, self.p_drhodz, self.p_dtheta_rhodz,
+                          self.p_du_fast, self.p_du_slow,
+                          self.p_dPhi_fast, self.p_dPhi_slow,
+                          self.p_dW_fast, self.p_dW_slow)
+        #time2=time.time()
+        #if time2>time1: elapsed=elapsed+time2-time1
+    def remap(self):
+        dynamico_remap(self.p_mass, self.p_theta_rhodz, self.p_u)
+
+def caldyn_step_TRSW(mesh,time_scheme,nstep):
+    setvars(('hydrostatic','caldyn_thermo','caldyn_eta'),
+            (True,thermo_boussinesq,eta_lag))
+    dynamico_init_params()
+    return Caldyn_step(mesh,time_scheme, nstep)
+def caldyn_step_HPE(mesh,time_scheme,nstep, caldyn_thermo,caldyn_eta, thermo,BC,g):
+    setvars(('hydrostatic','caldyn_thermo','caldyn_eta',
+             'g','ptop','Rd','cpp','preff','Treff'),
+             (True,caldyn_thermo,caldyn_eta,
+              g,BC.ptop,thermo.Rd,thermo.Cpd,thermo.p0,thermo.T0))
+    dynamico_init_params()
+    return Caldyn_step(mesh,time_scheme, nstep)
+def caldyn_step_NH(mesh,time_scheme,nstep, caldyn_thermo, caldyn_eta, thermo,BC,g):
+    setvars(('hydrostatic','caldyn_thermo','caldyn_eta',
+             'g','ptop','Rd','cpp','preff','Treff','pbot','rho_bot'),
+             (False,caldyn_thermo,caldyn_eta,
+              g,BC.ptop,thermo.Rd,thermo.Cpd,thermo.p0,thermo.T0,
+              BC.pbot.max(), BC.rho_bot.max()))
+    dynamico_init_params()
+    return Caldyn_step(mesh,time_scheme, nstep)
+    
+#----------------------------- Base class for dynamics ------------------------
+
+class Caldyn:
+    def __init__(self,mesh):
+        self.mesh=mesh
+        fps, ftheta, fmass = mesh.field_ps, mesh.field_theta, mesh.field_mass
+        fw, fu, fz = mesh.field_w, mesh.field_u, mesh.field_z
+        self.ps, self.ms, self.dms       = fps(), fps(), fps()
+        self.s, self.hs, self.dhs        = ftheta(), ftheta(), ftheta()
+        self.pk, self.berni, self.geopot, self.hflux = fmass(),fmass(),fw(),fu()
+        self.qu, self.qv                 = fu(),fz()
+        self.fmass, self.ftheta, self.fu, self.fw = fmass, ftheta, fu, fw
+    def bwd_fast_slow(self, flow, tau):
+        global elapsed
+        time1=time.time()
+        flow,fast,slow = self._bwd_fast_slow_(flow,tau)
+        time2=time.time()
+        elapsed=elapsed+time2-time1
+        return flow,fast,slow
+
+# when calling caldyn_unstructured, arrays for tendencies must be re-created each time
+# to avoid overwriting in the same memory space when time scheme is multi-stage
+
+#-------------------------- Shallow-water dynamics ---------------------
+
+class Caldyn_RSW(Caldyn):
+    def __init__(self,mesh):
+        Caldyn.__init__(self,mesh)
+        setvars(('hydrostatic','caldyn_thermo','caldyn_eta'),
+                (True,thermo_boussinesq,eta_lag))
+        self.dhs = self.fmass()
+        dynamico_init_params()
+    def _bwd_fast_slow_(self, flow, tau):
+        h,u = flow
+        # h*s = h => uniform buoyancy s=1 => shallow-water
+        dh, du_slow, du_fast, hs, buf = self.fmass(), self.fu(), self.fu(), h.copy(), self.geopot
+        assert type(tau) is np_num, 'tau must be of type unstructured.np_num.'
+        ker.dynamico_caldyn_unstructured(tau, self.ms, h, hs, u, self.geopot, buf,
+                  self.s, self.ps, self.pk, self.hflux, self.qv,
+                  self.dms, dh, self.dhs, du_fast, du_slow,
+                  buf, buf, buf, buf)
+        return (h,u), (0.,du_fast), (dh,du_slow)
+
+#----------------------------------- HPE ------------------------------------
+
+class Caldyn_HPE(Caldyn):
+    def __init__(self,caldyn_thermo,caldyn_eta, mesh,thermo,BC,g):
+        Caldyn.__init__(self,mesh)
+        setvars(('hydrostatic','caldyn_thermo','caldyn_eta',
+                 'g','ptop','Rd','cpp','preff','Treff'),
+                (True,caldyn_thermo,caldyn_eta,
+                 g,BC.ptop,thermo.Rd,thermo.Cpd,thermo.p0,thermo.T0))
+        dynamico_init_params()
+    def _bwd_fast_slow_(self, flow, tau):
+        dm, dS, du_slow, du_fast, buf = self.fmass(), self.ftheta(), self.fu(), self.fu(), self.geopot
+        m,S,u = flow
+        ker.dynamico_caldyn_unstructured(tau, self.ms, m, S, u, self.geopot, buf,
+                  self.s, self.ps, self.pk, self.hflux, self.qv,
+                  self.dms, dm, dS, du_fast, du_slow,
+                  buf, buf, buf, buf)
+        return (m,S,u), (0.,0.,du_fast), (dm,dS,du_slow)
+
+#----------------------------------- NH ------------------------------------
+
+class Caldyn_NH(Caldyn):
+    def __init__(self,caldyn_thermo,caldyn_eta, mesh,thermo,BC,g):
+        Caldyn.__init__(self,mesh)
+        setvars(('hydrostatic','caldyn_thermo','caldyn_eta',
+                 'g','ptop','Rd','cpp','preff','Treff',
+                 'pbot','rho_bot'),
+                (False,caldyn_thermo,caldyn_eta,
+                 g,BC.ptop,thermo.Rd,thermo.Cpd,thermo.p0,thermo.T0,
+                 BC.pbot.max(), BC.rho_bot.max()))
+        dynamico_init_params()
+    def bwd_fast_slow(self, flow, tau):
+        ftheta, fmass, fu, fw = self.ftheta, self.fmass, self.fu, self.fw 
+        dm, dS, du_slow, du_fast = fmass(), ftheta(), fu(), fu()
+        dPhi_slow, dPhi_fast, dW_slow, dW_fast = fw(), fw(), fw(), fw()
+        m,S,u,Phi,W = flow
+        ker.dynamico_caldyn_unstructured(tau, self.ms, m, S, u, Phi, W,
+                  self.s, self.ps, self.pk, self.hflux, self.qv,
+                  self.dms, dm, dS, du_fast, du_slow,
+                  dPhi_fast, dPhi_slow, dW_fast, dW_slow)
+        return ((m,S,u,Phi,W), (0.,0.,du_fast,dPhi_fast,dW_fast), 
+                (dm,dS,du_slow,dPhi_slow,dW_slow))
+
+#------------------------ Copy mesh info to Fortran side -------------------
+
+def init_mesh(llm, nqdyn, edge_num, primal_num, dual_num,
+              max_trisk_deg, max_primal_deg, max_dual_deg,
+              primal_nb, primal_edge, primal_ne,
+              dual_nb,dual_edge,dual_ne,dual_vertex, 
+              left,right,down,up,trisk_deg,trisk,
+              Ai, Av, fv, le_de, Riv2, wee):
+
+    DEBUG('Types of Ai, Av, fv, le_de, Riv2, wee : %s' % 
+          ((Ai.dtype,Av.dtype,fv.dtype,le_de.dtype,Riv2.dtype,wee.dtype),) )
+    for var,varname in zip((Ai,Av,fv,le_de,Riv2,wee), ('Ai','Av','fv','le_de','Riv2','wee')):
+        assert var.dtype == np.float64, '%s must be double precision'%varname
+
+    setvars( ('llm','nqdyn','edge_num','primal_num','dual_num',
+              'max_trisk_deg','max_primal_deg','max_dual_deg'),
+             (llm, nqdyn, edge_num, primal_num, dual_num, 
+              max_trisk_deg, max_primal_deg, max_dual_deg) )        
+    INFO('init_mesh ...')
+    ker.dynamico_init_mesh(primal_nb,primal_edge,primal_ne,
+              dual_nb,dual_edge,dual_ne,dual_vertex,
+              left,right,down,up,trisk_deg,trisk)
+    INFO('...done')
+    INFO('init_metric ...')
+    ker.dynamico_init_metric(Ai,Av,fv,le_de,Riv2,wee)
+    INFO('...done')
+
+#------------------------ Mesh partitioning ------------------------
+
+# Helper functions and interface to ParMETIS
+# loc_stencil returns the start/end indices (vtxdist) expected by ParMETIS
+# i.e. index[start:end] with start=vtxdist[cell], end=vtxdist[cell+1] lists the edges of cell 'cell'
+
+def loc_stencil(degree, stencil):
+    loc=0
+    for i in range(degree.size):
+        yield loc
+        loc=loc+degree[i]
+    yield loc
+
+def partition_mesh(degree, stencil, nparts): 
+    # arguments : PArray1D and PArray2D describing mesh, number of desired partitions
+    dim_cell, degree, stencil = degree.dim, degree.data, stencil.data
+    comm, vtxdist, idx_start, idx_end = dim_cell.comm, dim_cell.vtxdist, dim_cell.start, dim_cell.end
+    mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size()
+    adjncy_loc, xadj_loc = list_stencil(degree, stencil), loc_stencil(degree, stencil)
+    adjncy_loc, xadj_loc = [np.asarray(list(x), dtype=np.int32) for x in (adjncy_loc, xadj_loc)]
+    owner = np.zeros(idx_end-idx_start, dtype=np.int32);
+    ker.dynamico_partition_graph(mpi_rank, mpi_size, vtxdist, xadj_loc, adjncy_loc, nparts, owner)
+    return owner
Index: codes/icosagcm/devel/Python/dev/wrap.py
===================================================================
--- codes/icosagcm/devel/Python/dev/wrap.py	(revision 977)
+++ codes/icosagcm/devel/Python/dev/wrap.py	(revision 977)
@@ -0,0 +1,78 @@
+#------------------- wrap C/Fortran routines -------------------#                                                                                           
+from dynamico import getargs
+log_master, log_world = getargs.getLogger(__name__)
+INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
+INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
+
+
+import numpy as np 
+from ctypes import *
+
+class Struct: pass
+
+# converts numpy array to void * for calling C/Fortran routines                                                                                             
+def loc(array): return array.ctypes.data_as(c_void_p)
+
+# defines converter functions depending on type of Python argument
+def noop(x): return x
+py2c = { np.ndarray:loc, np.float64:c_double, np.float32:c_float, int:c_int, 
+         c_void_p:noop, str:noop, c_int:noop }
+
+class Wrap:
+    def __init__(self, fun, check_args, argtypes):
+        fun.restype=None
+        if check_args:
+            fun.argtypes=argtypes
+        self.fun=fun
+        self.convert=None
+    def __call__(self, *args):
+        if self.convert is None:
+            self.convert=[ py2c[type(arg)] for arg in args]
+        args = [ conv(arg) for conv,arg in zip(self.convert,args) ]
+        return self.fun(*args)
+
+def typelist(argtypes):
+    lst=[]
+    for x in argtypes:
+        if type(x) is int:
+            lst+=(x-1)*[prev] # x-1 because prev has already been appended once                                                             
+        else:
+            lst.append(x)
+        prev=x
+    return lst
+
+class SharedLib:
+    def __init__(self,vardict,lib,check_args=True, prefix_so='',prefix_py=''):
+        self.vardict, self.check_args = vardict, check_args
+        self.prefix_so, self.prefix_py = prefix_so, prefix_py
+        self.lib=lib
+        self.vars={}
+    def import_funs(self,prototypes,prefix=''):
+        for proto in prototypes:
+            names=[]
+            while type(proto[0]) is str:
+                names.append(proto.pop(0))
+            if proto[0] is None:
+                proto=None
+            else:
+                proto = typelist(proto)
+            for name in names:
+                name = prefix+name
+                soname = self.prefix_so+name
+                DEBUG('wrap.Sharedlib : importing %s' % soname)
+                self.vardict[self.prefix_py+name] = Wrap(self.lib[soname], self.check_args, proto)
+    def addvars(self,*types_and_names):
+        for x in types_and_names:
+            if type(x) is str:
+                self.vars[x]=ctype.in_dll(self.lib, x.lower())
+            else:
+                ctype = x
+    def getvar(self,name):
+        return self.vars[name].value
+    def getvars(self,*names):
+        return [self.vars[name].value for name in names]
+    def setvar(self,name,val):
+        self.vars[name].value=val
+    def setvars(self,names,vals):
+        for name,val in zip(names,vals):
+            self.vars[name].value=val
Index: codes/icosagcm/devel/Python/dev/xios.py
===================================================================
--- codes/icosagcm/devel/Python/dev/xios.py	(revision 977)
+++ codes/icosagcm/devel/Python/dev/xios.py	(revision 977)
@@ -0,0 +1,132 @@
+from __future__ import print_function
+import numpy as np
+import cxios
+from unstructured import ker
+from dynamico.meshes import radian
+from mpi4py import MPI
+
+from dynamico import getargs
+args=getargs.parse()
+log_master, log_world = getargs.getLogger(__name__)
+INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
+INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
+
+#-----------------------------------------------------------------------------
+
+class Client:
+    def __enter__(self):
+        ker.dynamico_setup_xios()
+        self.comm=MPI.COMM_WORLD # FIXME : we should use the communicator returned by XIOS
+        return self
+    def __exit__(self, type, value, traceback):
+        INFO('xios_finalize()')
+        DEBUG_ALL('xios_finalize')
+        cxios.finalize()
+    def min(self,data): return self.comm.allreduce(data, op=MPI.MIN)
+    def max(self,data): return self.comm.allreduce(data, op=MPI.MAX)
+    
+class Context:
+    def __enter__(self):
+        INFO('cxios.context_close_definition()')
+        cxios.context_close_definition()
+        return self
+    def __exit__(self, type, value, traceback):
+        INFO('xios_context_finalize()')
+        cxios.context_finalize()
+    def init_llm(self, mesh, nqtot):
+        self.mesh=mesh
+        # vertical axis, nq
+        llm = mesh.llm
+        lev, levp1, nq = [cxios.Handle('axis',name) for name in 'lev', 'levp1', 'nq']
+        lev.set_attr(n_glo=llm, value=np.arange(llm))
+        levp1.set_attr(n_glo=llm+1, value=np.arange(llm+1))
+        nq.set_attr(n_glo=nqtot, value=np.arange(nqtot))
+    def init_calendar(self, step):
+        calendar=cxios.Handle('calendar_wrapper')
+        dtime = cxios.Duration(second=step)
+        calendar.set_attr(timestep=dtime)
+        calendar.update_timestep()
+        std=cxios.Handle('fieldgroup','standard_output')
+        std.set_attr(freq_op=dtime)
+        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 send_field_primal(self,name, data):
+        self.send_field(self.mesh.primal_own_loc, name, data)
+    def send_field_dual(self,name, data):
+        self.send_field(self.mesh.dual_own_loc, name, data)
+    def send_field(self,own_loc, name, data):
+        if data.ndim==1:
+            data = data[own_loc]
+        if data.ndim==2:
+            data = data[own_loc,:]
+            data = data.transpose() # XIOS expects contiguous horizontal slices
+        data = np.ascontiguousarray(data)
+        cxios.send_field(name, data)
+
+#-----------------------------------------------------------------------------
+
+def setup_curvilinear(cat,id,nx,ny,own,i_index,j_index,lon,lat):
+    mesh = cxios.Handle(cat,id)
+    def log(name,data) : 
+        DEBUG('setup_curvilinear : %s shape min max %s %s %s'%(name, data.shape, data.min(), data.max()) )
+    i_index, j_index, lon, lat = [ x[own] for x in i_index, j_index, lon, lat ]
+    log('i_index',i_index)
+    log('j_index',j_index)
+    log('lon',lon)
+    log('lat',lat)
+    mesh.set_attr(type='curvilinear', ni_glo=nx, i_index=i_index, nj_glo=ny, j_index=j_index)
+    mesh.set_attr(lonvalue_1d=lon.flatten(), latvalue_1d=lat.flatten())
+
+class Context_Curvilinear(Context):
+    def __init__(self, mesh, nqtot, step):
+        self.init_llm(mesh, nqtot)
+        # primal mesh
+        setup_curvilinear('domaingroup','i', mesh.nx, mesh.ny, mesh.primal_own_loc, 
+                          mesh.primal_i, mesh.primal_j, mesh.lon_i, mesh.lat_i)
+        # dual mesh
+        setup_curvilinear('domaingroup','v', mesh.nx, mesh.ny, mesh.dual_own_loc, 
+                          mesh.dual_i, mesh.dual_j, mesh.lon_v, mesh.lat_v)
+        # calendar
+        self.init_calendar(step)
+
+#-----------------------------------------------------------------------------
+
+def cf_boundaries(degree,points,lon,lat):
+    n, nvertex = len(degree), degree.max()
+    bnd_lon, bnd_lat = np.zeros((n,nvertex)),np.zeros((n,nvertex))
+    for ij in range(n):
+        nb=degree[ij]
+        for k in range(nb):
+            vertex = points[ij,k]
+            bnd_lon[ij,k]=lon[vertex]
+            bnd_lat[ij,k]=lat[vertex]
+        for k in range(nb,nvertex):   # repeat last vertex if nb<nvertex
+            bnd_lon[ij,k]=bnd_lon[ij,nb-1]
+            bnd_lat[ij,k]=bnd_lat[ij,nb-1]
+    return bnd_lon, bnd_lat
+
+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_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 Context_Unstructured(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)
Index: codes/icosagcm/devel/Python/dynamico/meshes.py
===================================================================
--- codes/icosagcm/devel/Python/dynamico/meshes.py	(revision 976)
+++ codes/icosagcm/devel/Python/dynamico/meshes.py	(revision 977)
@@ -1,5 +1,3 @@
 from __future__ import absolute_import     
-
-from dynamico.dev import unstructured as unst
 
 import time
@@ -15,6 +13,5 @@
 from dynamico.util import list_stencil, BaseClass, inverse_list
 from dynamico import getargs
-# TODO from dynamico.parmetis import partition_mesh
-from dynamico.dev.unstructured import partition_mesh
+from dynamico.partition import partition_mesh
 
 log_master, log_world = getargs.getLogger(__name__)
@@ -419,5 +416,5 @@
                     vertex_lat[cell, ivertex] = vlat
             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.primal_bounds_lon, mesh.primal_bounds_lat = bounds(mesh.primal_own_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)
 
@@ -766,5 +763,5 @@
     nn, xx, yy, zz = x.size, toint(x), toint(y), toint(z)
     idx=np.zeros(nn, dtype=np.int64)
-    unst.ker.dynamico_morton_encode(nn, xx,yy,zz, idx)
+    partition.dynamico_morton_encode(nn, xx,yy,zz, idx)
     return idx
 
Index: codes/icosagcm/devel/Python/dynamico/partition.h
===================================================================
--- codes/icosagcm/devel/Python/dynamico/partition.h	(revision 977)
+++ codes/icosagcm/devel/Python/dynamico/partition.h	(revision 977)
@@ -0,0 +1,55 @@
+#include <mpi.h>
+#include <stdio.h>
+#include <parmetis.h>
+
+static inline void dynamico_partition_graph(int mpi_rank, int mpi_size, 
+                     idx_t vtxdist[], idx_t xadj_loc[], idx_t adjncy_loc[],
+                     idx_t nparts, idx_t part[])
+{
+  MPI_Comm comm = MPI_COMM_WORLD; /* avoid passing a communicator as argument due to Fortran/C issues */
+  idx_t edgecut;
+  idx_t options[]={0}, wgtflag=0, numflag=0, ncon=1;
+  real_t tpwgts[2*nparts], ubvec[]={1.05};
+
+  /* int ParMETIS V3 PartKway (
+     idx t *vtxdist, idx t *xadj, idx t *adjncy, idx t *vwgt, idx t *adjwgt, 
+     idx t *wgtflag, idx t *numflag, idx t *ncon, idx t *nparts, 
+     real t *tpwgts, real t *ubvec, idx t *options, idx t *edgecut, idx t *part, 
+     MPI Comm *comm )
+  */
+
+/*  printf("sizeof : %d %d \n", sizeof(int), sizeof(idx_t));
+  printf("partition_graph : %d %d %d\n", mpi_rank, mpi_size, nparts); 
+  printf("vtxdist : \n"); 
+  for(int i=0; i<mpi_size+1; i++) printf("%d %d %d\n", mpi_rank, i, vtxdist[i]); */
+      
+  for(int i=0; i<2*nparts; i++) tpwgts[i] = 1./nparts;
+  ParMETIS_V3_PartKway(vtxdist, xadj_loc, adjncy_loc, NULL, NULL, 
+		       &wgtflag, &numflag, &ncon, &nparts, 
+		       tpwgts, ubvec, options, &edgecut, part, &comm); 
+}
+
+
+// http://www.forceflow.be/2013/10/07/morton-encodingdecoding-through-bit-interleaving-implementations/ 
+// method to seperate bits from a given integer 3 positions apart
+static inline uint64_t splitBy3(unsigned int a){
+  uint64_t x = a & 0x1fffff; // we only look at the first 21 bits
+  x = (x | x << 32) & 0x1f00000000ffff;  // shift left 32 bits, OR with self, and 00011111000000000000000000000000000000001111111111111111
+  x = (x | x << 16) & 0x1f0000ff0000ff;  // shift left 32 bits, OR with self, and 00011111000000000000000011111111000000000000000011111111
+  x = (x | x << 8) & 0x100f00f00f00f00f; // shift left 32 bits, OR with self, and 0001000000001111000000001111000000001111000000001111000000000000
+  x = (x | x << 4) & 0x10c30c30c30c30c3; // shift left 32 bits, OR with self, and 0001000011000011000011000011000011000011000011000011000100000000
+  x = (x | x << 2) & 0x1249249249249249;
+  return x;
+}
+ 
+static inline uint64_t mortonEncode_magicbits(unsigned int x, unsigned int y, unsigned int z){
+  uint64_t answer = 0;
+  answer |= splitBy3(x) | splitBy3(y) << 1 | splitBy3(z) << 2;
+  return answer;
+}
+
+static inline void dynamico_morton_encode(int n, const int *x, const int *y, const int *z, uint64_t *m)
+{
+  for(int i=0; i<n; i++)
+    m[i]=mortonEncode_magicbits(x[i],y[i],z[i]);
+}
Index: codes/icosagcm/devel/Python/dynamico/partition.pyx
===================================================================
--- codes/icosagcm/devel/Python/dynamico/partition.pyx	(revision 977)
+++ codes/icosagcm/devel/Python/dynamico/partition.pyx	(revision 977)
@@ -0,0 +1,52 @@
+
+from ctypes import c_void_p, c_int, c_double, c_float, c_bool
+from numpy cimport ndarray
+cimport numpy as np
+import numpy as np
+
+ctypedef int idx_t
+
+cdef extern from "partition.h":
+    void dynamico_partition_graph(int mpi_rank, int mpi_size,
+                     idx_t vtxdist[], idx_t xadj_loc[], idx_t adjncy_loc[],
+                     idx_t nparts, idx_t part[])
+
+
+#------------------------ Mesh partitioning ------------------------
+
+# Helper functions and interface to ParMETIS
+# loc_stencil returns the start/end indices (vtxdist) expected by ParMETIS
+# i.e. index[start:end] with start=vtxdist[cell], end=vtxdist[cell+1] lists the edges of cell 'cell'
+# list_stencil converts an adjacency graph from array format index[num_cells, MAX_EDGES] to compressed format
+
+def list_stencil(degree, stencil, cond=lambda x:True):
+    for i in range(degree.size):
+        for j in range(degree[i]):
+            s=stencil[i,j]
+            if cond(s): yield stencil[i,j]
+
+def loc_stencil(degree, stencil):
+    loc=0
+    for i in range(degree.size):
+        yield loc
+        loc=loc+degree[i]
+    yield loc
+        
+def partition_mesh(degree, stencil, nparts): 
+    # arguments : PArray1D and PArray2D describing mesh, number of desired partitions
+    cdef int[::1] adjncy_loc, xadj_loc, vtxdist, owner
+    dim_cell, degree, stencil = degree.dim, degree.data, stencil.data
+    comm, vtxdist, idx_start, idx_end = dim_cell.comm, dim_cell.vtxdist, dim_cell.start, dim_cell.end
+    adjncy_loc, xadj_loc = [np.asarray(list(x), dtype=np.int32)
+                            for x in (list_stencil(degree, stencil), loc_stencil(degree, stencil)) ]
+    owner = np.zeros(idx_end-idx_start, dtype=np.int32)
+    cdef int mpi_rank = comm.Get_rank()
+    cdef int mpi_size = comm.Get_size()
+    dynamico_partition_graph(mpi_rank, mpi_size, &vtxdist[0], &xadj_loc[0], &adjncy_loc[0], nparts, &owner[0])
+    return owner
+
+def partition_graph(comm, int[::1] vtxdist, int[::1] xadj_loc, int[::1] adjncy_loc, int[::1] owner, int nparts):
+    cdef int mpi_rank = comm.Get_rank()
+    cdef int mpi_size = comm.Get_size()
+    dynamico_partition_graph(mpi_rank, mpi_size, &vtxdist[0], &xadj_loc[0], &adjncy_loc[0], nparts, &owner[0])
+    return owner
Index: codes/icosagcm/devel/Python/dynamico/setup.py
===================================================================
--- codes/icosagcm/devel/Python/dynamico/setup.py	(revision 977)
+++ codes/icosagcm/devel/Python/dynamico/setup.py	(revision 977)
@@ -0,0 +1,7 @@
+from distutils.core import setup
+from distutils.extension import Extension
+from Cython.Build import cythonize
+
+setup(  name = "Simple Cython interface to parmetis for DYNAMICO",
+        ext_modules = cythonize([ Extension("partition",["partition.pyx"],libraries=["parmetis"])] )
+    )
Index: codes/icosagcm/devel/Python/test/py/partition.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/partition.py	(revision 976)
+++ codes/icosagcm/devel/Python/test/py/partition.py	(revision 977)
@@ -3,7 +3,6 @@
 from dynamico import meshes
 from dynamico import parallel
-from dynamico.dev import unstructured as unst
 from dynamico import maps
-from dynamico.parmetis import partition_graph
+from dynamico import partition
 
 from mpi4py import MPI
@@ -11,6 +10,4 @@
 mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size()
 print '%d/%d starting'%(mpi_rank,mpi_size)
-
-#from dynamico import partition
 
 import math as math
@@ -50,6 +47,5 @@
     part = 0*xadj_loc[0:-1];
 
-#    unst.ker.dynamico_partition_graph(mpi_rank, mpi_size, vtxdist, xadj_loc, adjncy_loc, 4, part)
-    partition_graph(comm, vtxdist, xadj_loc, adjncy_loc, part, nparts=4)
+    partition.partition_graph(comm, vtxdist, xadj_loc, adjncy_loc, part, nparts=4)
 
     for i in range(len(part)):
@@ -66,5 +62,5 @@
     mydegree, mybounds = [get_mycells(x) for x in primal_deg, primal_vertex]
     print '%d : len(mydegree)=%d'%(mpi_rank, len(mydegree))
-    vertex_list = sorted(set(unst.list_stencil(mydegree,mybounds))) 
+    vertex_list = sorted(set(partition.list_stencil(mydegree,mybounds))) 
     print '%d : len(vertex_list))=%d'%(mpi_rank, len(vertex_list))
     get_myvertices = parallel.Get_Indices(dim_vertex, vertex_list)
@@ -116,5 +112,5 @@
 plt.xlim(-190.,190.)
 plt.ylim(-90.,90.)
-plt.savefig('fig_partition/A%03d.pdf'%mpi_rank, dpi=1600)
+plt.savefig('fig_partition/A%03d.png'%mpi_rank, dpi=160)
 
 fig, ax = plt.subplots()
@@ -124,3 +120,3 @@
 plt.xlim(-190.,190.)
 plt.ylim(-90.,90.)
-plt.savefig('fig_partition/B%03d.pdf'%mpi_rank, dpi=1600)
+plt.savefig('fig_partition/B%03d.png'%mpi_rank, dpi=160)
Index: codes/icosagcm/devel/make_python
===================================================================
--- codes/icosagcm/devel/make_python	(revision 976)
+++ codes/icosagcm/devel/make_python	(revision 977)
@@ -13,6 +13,4 @@
     rm -rf *.o *.so build/
 }
-
-#------------------ Default action , invoked by : './make_python' --------------------
 
 function make_so()
@@ -33,10 +31,8 @@
 }
 
-
-function cmd_()
+function check_python_modules()
 {
-    echo "XIOS_ROOTDIR : $XIOS_ROOTDIR"
     MISSING=""
-    for module in netCDF4 jinja2 mpi4py; do
+    for module in $*; do
 	echo "import $module"
 	python -c "import $module" || MISSING="$MISSING $module"
@@ -49,24 +45,62 @@
 	exit
     fi
+}
 
-    log mpicc -c -std=c99 -fPIC $PARMETIS_CFLAGS Python/src/partition.c -o obj/partition.o
-    make_so $XIOS_ROOTDIR/lib/libxios.a $LIBDIR/libxios.so "$NETCDF_LDFLAGS -lstdc++"
-    make_so $LIBDIR/libicosa.a $LIBDIR/libicosa.so "-L$LIBDIR -lxios $KERNELS_LDFLAGS" "physics_external_mod.o physics_mod.o timeloop_gcm_mod.o icosa_init_mod.o" obj/partition.o
-
-    # Cythonize Python modules
-    cd $DYNAMICO_ROOT/Python/src
-    NUMPY=$(python -c 'import numpy as np ; print np.get_include()')
-    CYTHONFLAGS="$SOFLAGS $LDFLAGS -L$LIBDIR -licosa $KERNELS_LDFLAGS"
-    LDSHARED="$LINK" CFLAGS="-fno-strict-aliasing -O2 -I$NUMPY" LDFLAGS="$CYTHONFLAGS" log python setup.py build_ext --inplace || exit
-
-    cp -fpr $DYNAMICO_ROOT/Python/dynamico $DYNAMICO_ROOT/lib
-    cp -fp *.so $DYNAMICO_ROOT/lib/dynamico/dev/
+function check_dynamico_modules()
+{
     cd $DYNAMICO_ROOT
-
-    for module in dev.cxios dev.unstructured meshes dyn time_step DCMIP; do
+    for module in $*; do
 	echo "import dynamico.$module"
 	LD_PRELOAD=$PYTHON_PRELOAD python -c "import dynamico.$module"
     done
 }
+
+function cythonize()
+{
+    EXTRA=$1
+    NUMPY=$(python -c 'import numpy as np ; print np.get_include()')
+    CYTHONFLAGS="$SOFLAGS $LDFLAGS -L$LIBDIR $KERNELS_LDFLAGS $EXTRA"
+    CC=mpicc LDSHARED="$LINK" CFLAGS="-fno-strict-aliasing -O2 -I$NUMPY $PARMETIS_CFLAGS" LDFLAGS="$CYTHONFLAGS" log python setup.py build_ext --inplace || exit
+    rm -rf build *.c # temp files
+}
+
+#------------------ Default action , invoked by : './make_python' --------------------
+
+function make_user()
+{
+    check_python_modules netCDF4 mpi4py
+    # Cythonize Python modules
+    cd $DYNAMICO_ROOT/Python/dynamico
+    cythonize
+    rm -rf build *.c # temp files
+    cp -fpr $DYNAMICO_ROOT/Python/dynamico $DYNAMICO_ROOT/lib
+    check_dynamico_modules partition meshes dyn time_step DCMIP
+}
+
+function make_dev()
+{
+    check_python_modules jinja2
+    echo "XIOS_ROOTDIR : $XIOS_ROOTDIR"
+
+    make_so $XIOS_ROOTDIR/lib/libxios.a $LIBDIR/libxios.so "$NETCDF_LDFLAGS -lstdc++"
+    make_so $LIBDIR/libicosa.a $LIBDIR/libicosa.so "-L$LIBDIR -lxios $KERNELS_LDFLAGS" "physics_external_mod.o physics_mod.o timeloop_gcm_mod.o icosa_init_mod.o"
+
+    cd $DYNAMICO_ROOT/Python/dev
+    cythonize "-licosa"
+    cp -fpr $DYNAMICO_ROOT/Python/dev $DYNAMICO_ROOT/lib/dynamico/
+    check_dynamico_modules dev.cxios dev.unstructured
+}
+
+function cmd_()
+{
+    make_user
+}
+
+function cmd_dev()
+{
+    make_user
+    make_dev
+}
+
 
 #------------------------- script starts here ------------------------
