Index: codes/icosagcm/devel/Python/dynamico/meshes.py
===================================================================
--- codes/icosagcm/devel/Python/dynamico/meshes.py	(revision 673)
+++ codes/icosagcm/devel/Python/dynamico/meshes.py	(revision 674)
@@ -5,5 +5,5 @@
 import matplotlib.tri as tri
 import matplotlib.pyplot as plt
-from unstructured import init_mesh, list_stencil
+from unstructured import init_mesh, list_stencil, ker
 
 radian=180/math.pi
@@ -174,8 +174,8 @@
     return ne
 
-def plot(tri,data):
+def plot(tri,data, **kwargs):
     plt.figure(figsize=(12,4))
     plt.gca().set_aspect('equal')
-    plt.tricontourf(tri, data, 20)
+    plt.tricontourf(tri, data, 20, **kwargs)
     plt.colorbar()
     plt.ylim((-90,90))
@@ -210,6 +210,4 @@
 
         # get indices for stencils 
-        # primal -> vertices (unused)
-        primal_vertex, dual_vertex = getvars('verticesOnCell','cellsOnVertex')
         # primal <-> edges
         primal_edge, left_right = getvars('edgesOnCell','cellsOnEdge')
@@ -221,5 +219,5 @@
         dual_ne = compute_ne(dual_num,dual_deg,dual_edge,up,down)
         # primal <-> dual, edges <-> edges (TRiSK)
-        dual_vertex, trisk = getvars('cellsOnVertex','edgesOnEdge')
+        primal_vertex, dual_vertex, trisk = getvars('verticesOnCell','cellsOnVertex','edgesOnEdge')
         # get positions, lengths, surfaces and weights
         le,de,Ai,Av = getvars('dvEdge','dcEdge','areaCell','areaTriangle')
@@ -296,12 +294,57 @@
             ucov[edge,:] = self.de[edge]*(ulon[edge,:]*math.cos(angle)+ulat[edge,:]*math.sin(angle))
         return ucov
-    def plot_i(self,data):
-        plot(self.primal,data)
-    def plot_v(self,data):
-        plot(self.dual,data)
-    def plot_e(self,data):
-        plot(self.triedge,data)
+    def plot_i(self,data, **kwargs):
+        plot(self.primal,data, **kwargs)
+    def plot_v(self,data, **kwargs):
+        plot(self.dual,data, **kwargs)
+    def plot_e(self,data, **kwargs):
+        plot(self.triedge,data, **kwargs)
+
+#--------------------------------------  Mesh reordering  ------------------------------------------#
+
+# NB : Indices start at 0
+# permutations map an old index to a new index : perm[old_index]=new_index
+# enumerate(perm) is the list [ (0,perm[0]) , (1,perm[1]), ...] = [(old, new) ... ]
+
+def reorder_values(perm, *datas): 
+    # data in datas : ndarray containing indices ; perm : permutation to apply to values
+    for data in datas:
+        for x in np.nditer(data, op_flags=['readwrite']):
+            x[...]=perm[x]
+
+def reorder_indices(perm, *datas): 
+    # data in datas : ndarray containing values ; perm : permutation to apply to indices
+    for data in datas:
+        data_out, ndim = data.copy(), data.ndim
+        if ndim == 1:
+            for old,new in enumerate(perm):
+                data_out[new]=data[old]
+        if ndim == 2:
+            for old,new in enumerate(perm):
+                data_out[new,:]=data[old,:]
+        data[:]=data_out[:]
+
+def key_to_perm(keys): 
+    # keys : maps an old index to a key
+    # returns : a list mapping old indices to new indices obtained by sorting keys in ascending order
+    ordered = [old for key,old in sorted(zip(keys,range(len(keys))))] 
+    # ordered maps new indices to old indices
+    # but we want to map old indices to new indices
+    perm = list(range(len(keys)))
+    for new,old in enumerate(ordered):
+        perm[old]=new
+    return perm  # maps old indices to new indices
+
+def toint(x): return np.int32((x+1.)*65536)
+    
+def morton_index(x,y,z):
+    nn, xx, yy, zz = x.size, toint(x), toint(y), toint(z)
+    idx=np.zeros(nn, dtype=np.int64)
+    ker.dynamico_morton_encode(nn, xx,yy,zz, idx)
+    return idx
 
 #-------------------------------------- Mesh partitioning ------------------------------------------#
+
+# NB : Indices start at 0
 
 def partition_from_stencil(owner2, degree, stencil):
@@ -331,4 +374,6 @@
 
 #---------------------------------- Stencil management ---------------------------------------#
+
+# NB : Indices start at 0
 
 # Class Stencil represents an adjacency relationship (e.g. cell=>edges)
Index: codes/icosagcm/devel/Python/src/partition.c
===================================================================
--- codes/icosagcm/devel/Python/src/partition.c	(revision 673)
+++ codes/icosagcm/devel/Python/src/partition.c	(revision 674)
@@ -29,2 +29,27 @@
 		       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;
+}
+
+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/src/unstructured.pyx
===================================================================
--- codes/icosagcm/devel/Python/src/unstructured.pyx	(revision 673)
+++ codes/icosagcm/devel/Python/src/unstructured.pyx	(revision 674)
@@ -44,5 +44,5 @@
 except OSError:
     print """
-Unable to load shared library 'libkernels.so' !
+Unable to load shared library 'libicosa.so' !
     """
     raise
@@ -59,4 +59,5 @@
                      ['dynamico_caldyn_unstructured', c_double, c_void_p,20],
                      ['dynamico_partition_graph', c_int,2, c_void_p,3, c_int, c_void_p],
+                     ['dynamico_morton_encode', c_int,c_void_p,4]
                      ])
 
Index: codes/icosagcm/devel/Python/test/py/reorder.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/reorder.py	(revision 674)
+++ codes/icosagcm/devel/Python/test/py/reorder.py	(revision 674)
@@ -0,0 +1,86 @@
+print 'Loading DYNAMICO modules ...'
+from dynamico.meshes import morton_index, key_to_perm, reorder_indices, reorder_values
+print '...Done'
+import numpy as np
+import netCDF4 as cdf
+import matplotlib.pyplot as plt
+import time
+
+grid = 10242 # 2562, 10242, 40962
+
+#-------------- read MPAS mesh ----------------#
+
+print 'Reading MPAS mesh ...'
+mesh=cdf.Dataset('grids_MPAS/x1.%d.grid.nc'%grid, "r")
+
+def getdims(*names): return [len(mesh.dimensions[name]) for name in names]
+def getvars(*names): 
+    for name in names : 
+        print "getvar %s ..."%name
+        time1=time.time()
+        ret=[mesh.variables[name][:] for name in names]
+        print "... Done in %f seconds"%(time.time()-time1)
+    return ret
+
+primal_num, edge_num, dual_num = getdims('nCells','nEdges','nVertices')                                   
+print 'Number of primal cells, dual cells and edges : %d, %d, %d '%(primal_num,dual_num,edge_num)
+# get degree of primal, dual cells
+primal_deg, trisk_deg = getvars('nEdgesOnCell','nEdgesOnEdge')
+# get indices for stencils 
+primal_edge, left_right = getvars('edgesOnCell','cellsOnEdge')
+dual_edge, down_up = getvars('edgesOnVertex','verticesOnEdge')
+primal_vertex, dual_vertex, trisk = getvars('verticesOnCell','cellsOnVertex','edgesOnEdge')
+# get positions, lengths, surfaces and weights
+le,de,Ai,Av = getvars('dvEdge','dcEdge','areaCell','areaTriangle')
+lat_i,lon_i = getvars('latCell','lonCell')
+lat_v,lon_v = getvars('latVertex','lonVertex')
+lat_e,lon_e,angle_e = getvars('latEdge','lonEdge','angleEdge')
+wee,Riv2 = getvars('weightsOnEdge','kiteAreasOnVertex')
+mesh.close()
+
+#------------ reorder mesh according to Morton index --------------#
+
+def morton_order(lon,lat):
+    x,y,z = np.cos(lat)*np.cos(lon), np.cos(lat)*np.sin(lon), np.sin(lat)
+    keys = morton_index(x,y,z)
+    return key_to_perm(keys)
+
+perm_i = morton_order(lon_i,lat_i)
+perm_v = morton_order(lon_v,lat_v)
+perm_e = morton_order(lon_e,lat_e)
+
+# MPAS => Python
+for idx in primal_edge, left_right, dual_edge, down_up, primal_vertex, dual_vertex, trisk : idx[:]=idx[:]-1
+
+reorder_indices(perm_i, primal_deg, primal_vertex, primal_edge, lon_i, lat_i, Ai)
+reorder_indices(perm_v, dual_vertex, dual_edge, lon_v, lat_v, Av, Riv2)
+reorder_indices(perm_e, trisk_deg, left_right, down_up, trisk, lon_e, lat_e, angle_e, le, de, wee)
+reorder_values(perm_i, dual_vertex, left_right)
+reorder_values(perm_v, primal_vertex, down_up)
+reorder_values(perm_e, primal_edge, dual_edge, trisk)
+
+# Python => MPAS
+for idx in primal_edge, left_right, dual_edge, down_up, primal_vertex, dual_vertex, trisk : idx[:]=idx[:]+1
+
+#--- Open copy of MPAS mesh file ; keep file structure but replace contents of variables ---#
+
+mesh=cdf.Dataset('grids/x1.%d.grid.nc'%grid, "r+")
+
+def putvars(datas,names): 
+    for name,data in zip(names,datas) :
+        print "putvar %s ..."%name
+        time1=time.time()
+        mesh.variables[name][:] = data[:]
+        print "... Done in %f seconds"%(time.time()-time1)
+    
+putvars( (primal_deg, trisk_deg), ('nEdgesOnCell','nEdgesOnEdge') )
+putvars( (primal_edge, left_right), ('edgesOnCell','cellsOnEdge') )
+putvars( (dual_edge, down_up), ('edgesOnVertex','verticesOnEdge') )
+putvars( (primal_vertex, dual_vertex, trisk), ('verticesOnCell','cellsOnVertex','edgesOnEdge') )
+putvars( (le,de,Ai,Av), ('dvEdge','dcEdge','areaCell','areaTriangle') )
+putvars( (lat_i,lon_i), ('latCell','lonCell') )
+putvars( (lat_v,lon_v), ('latVertex','lonVertex') )
+putvars( (lat_e,lon_e,angle_e), ('latEdge','lonEdge','angleEdge') )
+putvars( (wee,Riv2), ('weightsOnEdge','kiteAreasOnVertex') )
+
+mesh.close()
