Index: codes/icosagcm/devel/Python/dynamico/meshes.py
===================================================================
--- codes/icosagcm/devel/Python/dynamico/meshes.py	(revision 680)
+++ codes/icosagcm/devel/Python/dynamico/meshes.py	(revision 681)
@@ -10,5 +10,5 @@
 import unstructured as unst
 import parallel
-from util import list_stencil, Base_class
+from util import list_stencil, Base_class, inverse_list
 
 radian=180/math.pi
@@ -275,5 +275,5 @@
             if left[edge]==cell+1: ne[cell,iedge]+=1
             if right[edge]==cell+1: ne[cell,iedge]-=1
-            if ne[cell,iedge]==0 : print 'error at cell,iedge', cell, iedge
+#            if ne[cell,iedge]==0 : print 'compute_ne error at cell,iedge', cell, iedge
     return ne
 
@@ -299,5 +299,5 @@
         plt.colorbar()
         plt.ylim((-90,90))
-        plt.xlim((0,360))
+        plt.xlim((-180,180))
     def plot_i(self,data, **kwargs):
         self.plot(self.primal,data, **kwargs)
@@ -309,5 +309,4 @@
 class Unstructured_Mesh(Abstract_Mesh):
     def __init__(self, gridfile, llm, nqdyn, radius, f):
-        # open mesh file, get main dimensions
         getvars = gridfile.getvars
         # get positions, lengths, surfaces and weights
@@ -318,7 +317,5 @@
         primal_num, dual_num, edge_num = [ x.size for x in Ai, Av, le ]
         gridfile.primal_num, gridfile.dual_num = primal_num, dual_num
-
-        primal_deg, dual_deg, trisk_deg = getvars('primal_deg','dual_deg','trisk_deg')
-        
+        primal_deg, dual_deg, trisk_deg = getvars('primal_deg','dual_deg','trisk_deg')        
         # get indices for stencils 
         # primal <-> edges
@@ -349,7 +346,7 @@
                   left,right,down,up,trisk_deg,trisk,
                   self.Ai,self.Av,self.fv,le/de,Riv2,wee)
-        self.primal  = tri.Triangulation(lon_i*180./math.pi, lat_i*180./math.pi)
-        self.dual    = tri.Triangulation(lon_v*180./math.pi, lat_v*180./math.pi)
-        self.triedge = tri.Triangulation(lon_e*180./math.pi, lat_e*180./math.pi)        
+        self.primal  = tri.Triangulation(lon_i*radian, lat_i*radian)
+        self.dual    = tri.Triangulation(lon_v*radian, lat_v*radian)
+        self.triedge = tri.Triangulation(lon_e*radian, lat_e*radian)        
         
 class Unstructured_PMesh(Abstract_Mesh): # Mesh data distributed across MPI processes
@@ -358,10 +355,11 @@
         dim_primal, dim_edge, dim_dual = gridfile.get_pdims(comm,
             'primal_num','edge_num','dual_num')
-        (primal_deg, primal_vertex, primal_edge, lon_i, 
-         lat_i) = get_pvars(dim_primal, 'primal_deg', 'primal_vertex', 'primal_edge', 'lon_i', 'lat_i' )
-        (edge_deg, trisk_deg, left_right, down_up, 
-         trisk) = get_pvars(dim_edge, 'edge_deg', 'trisk_deg', 'left_right', 'down_up', 'trisk')
-        dual_deg, dual_vertex, dual_edge, Riv2, lon_v, lat_v = get_pvars(dim_dual, 
-            'dual_deg', 'dual_vertex', 'dual_edge', 'Riv2', 'lon_v', 'lat_v')
+        primal_deg, primal_vertex, primal_edge, lon_i, lat_i, Ai = get_pvars(
+            dim_primal, 'primal_deg', 'primal_vertex', 'primal_edge', 'lon_i', 'lat_i', 'Ai' )
+        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', 'left_right', 'down_up', 
+             'trisk', 'le', 'de', 'wee', 'lon_e', 'lat_e', 'angle_e')
+        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')
         # Let indices start at 0
         for x in (primal_vertex, primal_edge, dual_vertex,dual_edge,
@@ -371,5 +369,6 @@
         self.edge2vertex = Stencil_glob(edge_deg, down_up)
         self.vertex2edge = Stencil_glob(dual_deg, dual_edge)
-        self.edge2edge   = Stencil_glob(trisk_deg, trisk)
+        self.vertex2cell = Stencil_glob(dual_deg, dual_vertex, Riv2)
+        self.edge2edge   = Stencil_glob(trisk_deg, trisk, wee)
         print 'Partitioning ...'
         edge_owner = unst.partition_mesh(trisk_deg, trisk, comm.Get_size())
@@ -377,8 +376,9 @@
         primal_owner = partition_from_stencil(edge_owner, primal_deg, primal_edge)
         primal_owner = parallel.LocPArray1D(dim_primal, primal_owner)
-        self.set_members(locals(), 'dim_primal', 'dim_dual', 'edge_owner', 
-                         'primal_owner', 'primal_deg', 'primal_vertex', 
-                         'lon_v', 'lat_v', 'lon_i', 'lat_i')
-        #        self.dim_primal, self.edge_owner, self.primal_owner = dim_primal, edge_owner, primal_owner
+        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', 
+                         'le', 'de', 'lon_e', 'lat_e', 'angle_e')
+    def get(self, getter, *names): return [ getter(self.__dict__[x]) for x in names ]
 
 def chain(start, links):
@@ -404,7 +404,8 @@
         yield polygon
 
-class Local_Mesh(Base_class):
-    def __init__(self, pmesh):
-        dim_primal, edge_owner, primal_owner = pmesh.dim_primal, pmesh.edge_owner, pmesh.primal_owner
+class Local_Mesh(Abstract_Mesh):
+    def __init__(self, pmesh, llm, nqdyn, radius, f):
+        dim_primal, primal_owner, dim_edge, edge_owner, dim_dual = pmesh.members(
+            'dim_primal', 'primal_owner', 'dim_edge', 'edge_owner', 'dim_dual')
         print 'Constructing halos ...'
         edges_E0 = find_my_cells(edge_owner)
@@ -414,8 +415,64 @@
         cells_C0, cells_C1 = progressive_list(cells_C0, cells_C1)
         print 'E2,E1,E0 ; C1,C0 : ', map(len, (edges_E2, edges_E1, edges_E0, cells_C1, cells_C0))
-        mycells, halo_cells = cells_C0, cells_C1
-        get_mycells, get_halo_cells = dim_primal.getter(mycells), dim_primal.getter(halo_cells)
-        self.com_primal = parallel.Halo_Xchange(42, dim_primal, halo_cells, get_halo_cells(primal_owner))
-    def plot(self, ax, clim, degree, bounds, lon, lat, data):
+        get_own_cells, get_all_cells = dim_primal.getter(cells_C0), dim_primal.getter(cells_C1)
+        get_all_duals, get_all_edges = dim_dual.getter(vertices_V1), dim_edge.getter(edges_E2)
+        self.com_primal = parallel.Halo_Xchange(
+            42, dim_primal, cells_C1, get_all_cells(primal_owner))
+        self.com_edges = parallel.Halo_Xchange(
+            73, dim_edge, edges_E2, get_all_edges(edge_owner))
+        
+        self.lon_i, self.lat_i, self.Ai = pmesh.get(get_all_cells, 'lon_i', 'lat_i', 'Ai')
+        self.lon_v, self.lat_v, self.Av = pmesh.get(get_all_duals, 'lon_v', 'lat_v', 'Av')
+        self.fv = f(self.lon_v,self.lat_v) # Coriolis parameter
+        self.le, self.de, self.lon_e, self.lat_e, self.angle_e = pmesh.get(
+            get_all_edges, 'le', 'de', 'lon_e', 'lat_e', 'angle_e')
+        primal_num, dual_num, edge_num = map(len, (cells_C1, vertices_V1, edges_E2))
+        
+        dict_E2, dict_C1, dict_V1 = map(inverse_list, (edges_E2, cells_C1, vertices_V1))
+        down_up = pmesh.edge2vertex(edges_E2, dict_V1)
+        left_right = pmesh.edge2cell(edges_E2, dict_C1)
+        primal_edge = pmesh.cell2edge(cells_C1, dict_E2)
+        dual_edge = pmesh.vertex2edge(vertices_V1, dict_E2)
+        trisk = pmesh.edge2edge(edges_E2, dict_E2)
+        Riv2  = pmesh.vertex2cell(vertices_V1, dict_C1)
+        print 'Edge-vertex degree min/max %d %d'%(down_up.degree.min(), down_up.degree.max() )
+        print 'Edge-cell degree min/max %d %d'%(left_right.degree.min(), left_right.degree.max() )
+        print 'Primal degree min/max : %d %d'%(primal_edge.degree.min(), primal_edge.degree.max() )
+        print 'Dual degree min/max : %d %d'%(dual_edge.degree.min(), dual_edge.degree.max() )
+        print 'TRiSK degree min/max : %d %d'%(trisk.degree.min(), trisk.degree.max() )
+        print 'Riv2 degree min/max : %d %d'%(Riv2.degree.min(), Riv2.degree.max() )
+        primal_deg, primal_edge = primal_edge.degree, primal_edge.neigh_loc
+        dual_deg, dual_edge = dual_edge.degree, dual_edge.neigh_loc
+        trisk_deg, trisk, wee = trisk.degree, trisk.neigh_loc, trisk.weight
+        dual_deg, dual_vertex, Riv2 = Riv2.degree, Riv2.neigh_loc, Riv2.weight
+        edge_deg, left, right = left_right.degree, left_right.neigh_loc[:,0], left_right.neigh_loc[:,1]
+        edge_deg, down, up = down_up.degree, down_up.neigh_loc[:,0], down_up.neigh_loc[:,1]
+        for edge in range(edge_num): 
+            if edge_deg[edge]<2: up[edge]=down[edge]
+        max_primal_deg, max_dual_deg, max_trisk_deg = [x.shape[1] for x in primal_edge, dual_edge, trisk]
+        # compute_ne(...) and normalize(...) expect indices starting at 1
+        trisk, dual_vertex, primal_edge, dual_edge =  trisk+1, dual_vertex+1, primal_edge+1, dual_edge+1
+        left, right, down, up = left+1, right+1, down+1, up+1
+        primal_ne = compute_ne(primal_num,primal_deg,primal_edge,left,right)
+        dual_ne = compute_ne(dual_num,dual_deg,dual_edge,up,down)
+        self.set_members(locals(), 'nqdyn', 'llm', 'primal_num', 'dual_num', 'edge_num',
+                         'trisk_deg', 'trisk', 'wee', 'dual_deg', 'dual_vertex', 'Riv2')
+
+        pmesh.gridfile.normalize(self,radius)
+        unst.init_mesh(llm, nqdyn, self.edge_num, self.primal_num, self.dual_num,
+                  max_trisk_deg, max_primal_deg, max_dual_deg,
+                  primal_deg,primal_edge,primal_ne,
+                  dual_deg,dual_edge,dual_ne,dual_vertex,
+                  left,right,down,up,trisk_deg,trisk,
+                  self.Ai,self.Av,self.fv,self.le/self.de,Riv2,wee)
+        # llm must be set before we call set_dynamico_transfer
+        self.com_primal.set_dynamico_transfer('primal')
+        self.com_edges.set_dynamico_transfer('edge')
+
+        self.primal  = tri.Triangulation(self.lon_i*radian, self.lat_i*radian)
+        self.dual    = tri.Triangulation(self.lon_v*radian, self.lat_v*radian)
+        self.triedge = tri.Triangulation(self.lon_e*radian, self.lat_e*radian)        
+
+    def plot_patches(self, ax, clim, degree, bounds, lon, lat, data):
         nb_vertex = lon.size # global
         p = list(patches(degree, bounds, lon, lat))
@@ -486,5 +543,5 @@
     return owner1
             
-def find_my_cells(owner): # a PArray1D containing the data returned by partition_mesh
+def find_my_cells(owner): # owner : a PArray1D containing the data returned by partition_mesh
     dim, comm, owner = owner.dim, owner.dim.comm, owner.data
     mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size()
@@ -511,5 +568,6 @@
     for i in range(degree.size):
         for j in range(degree[i]):
-            bounds[i,j] = vertex_dict[bounds[i,j]]    
+            bounds[i,j] = vertex_dict[bounds[i,j]]
+    return bounds
 
 class Stencil_glob:
@@ -546,5 +604,6 @@
         myneigh_loc = reindex(neigh_dict, mydegree, myneigh)
         self.degree, self.neigh_glob, self.neigh_loc = mydegree, myneigh, myneigh_loc
-
+        if weight is not None: self.weight = myweight
+            
 def progressive_iter(mylist, cell_lists):
     for thelist in cell_lists: 
Index: codes/icosagcm/devel/Python/dynamico/parallel.py
===================================================================
--- codes/icosagcm/devel/Python/dynamico/parallel.py	(revision 680)
+++ codes/icosagcm/devel/Python/dynamico/parallel.py	(revision 681)
@@ -9,5 +9,7 @@
         order = [order[i] for i in range(len(lst))] # dict->list
         return lst,order  
-    
+
+from unstructured import ker
+
 #----------------------- Parallel access to distributed array ------------------------------#
 
@@ -143,6 +145,19 @@
         # asssociate local index to global index
         own_dict = { glob:loc for glob,loc in zip(own_global, own_local) }
-        self.send_list = [(rank, [own_dict[i] for i in send[rank]]) for rank in send_rank]                    
-        
+        self.send_list = [(rank, [own_dict[i] for i in send[rank]]) for rank in send_rank]
+    def set_dynamico_transfer(self,name):
+        def list_to_args(lst):
+            ranks, lst = zip(*lst) # list of tuples => tuple of lists
+            lens = map(len, lst)
+            return len(ranks), ranks, lens, sum(lens), sum(lst,[])
+        index = { 'primal':1, 'edge':2, 'dual':3 }[name]
+        self.index = index
+        send_num, send_rank, send_len, send_size, send_list = list_to_args(self.send_list) 
+        recv_num, recv_rank, recv_len, recv_size, recv_list = list_to_args(self.recv_list) 
+        send_rank, send_len, send_list, recv_rank, recv_len, recv_list = [ np.asarray(x, dtype=np.int32)
+            for x in send_rank, send_len, send_list, recv_rank, recv_len, recv_list ]
+        send_list, recv_list = send_list+1, recv_list+1 # Fortran expects that indices start at 1
+        ker.dynamico_init_transfer(index, send_num, send_size, send_rank, send_len, send_list,
+                                       recv_num, recv_size, recv_rank, recv_len, recv_list)
 class LocalArray: # a base class for arrays with halos
     def read_own(self, parray): self.put(self.halo.own_local, self.halo.get_own(parray))
Index: codes/icosagcm/devel/Python/dynamico/util.py
===================================================================
--- codes/icosagcm/devel/Python/dynamico/util.py	(revision 680)
+++ codes/icosagcm/devel/Python/dynamico/util.py	(revision 681)
@@ -13,3 +13,4 @@
         for name in names :
             self.__dict__[name]=loc[name]
-
+    def members(self, *names):
+        return [ self.__dict__[name] for name in names ]
Index: codes/icosagcm/devel/Python/src/unstructured.pyx
===================================================================
--- codes/icosagcm/devel/Python/src/unstructured.pyx	(revision 680)
+++ codes/icosagcm/devel/Python/src/unstructured.pyx	(revision 681)
@@ -60,4 +60,6 @@
                      ['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_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],
                      ['dynamico_morton_encode', c_int,c_void_p,4]
                      ])
Index: codes/icosagcm/devel/Python/test/py/RSW2_MPAS_W02.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/RSW2_MPAS_W02.py	(revision 680)
+++ codes/icosagcm/devel/Python/test/py/RSW2_MPAS_W02.py	(revision 681)
@@ -1,5 +1,13 @@
+print 'Starting'
+
+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)
+prefix='fig_RSW2_MPAS_W02/%02d'%mpi_rank
+
 print 'Loading DYNAMICO modules ...'
 from dynamico import unstructured as unst
-from dynamico.meshes import MPAS_Format, Unstructured_Mesh as Mesh
+from dynamico.meshes import MPAS_Format, Unstructured_PMesh as PMesh, Local_Mesh as Mesh
 from dynamico import time_step
 print '...Done'
@@ -11,7 +19,7 @@
 print '...Done'
 
-grid, llm, nqdyn = 10242, 1,1 # 2562, 10242, 40962
+grid, llm, nqdyn = 2562, 1,1 # 2562, 10242, 40962
 Omega, radius, g, gh0 = 2.*np.pi/86400., 6.4e6, 1., 2.94e4
-N, T, courant = 40, 10800., 1.2 # simulation length = N*T
+N, T, courant = 40, 10400., 1.2 # simulation length = N*T
 
 print 'Omega, planetary PV', Omega, 2*Omega/gh0
@@ -20,5 +28,6 @@
 print 'Reading MPAS mesh ...'
 meshfile = MPAS_Format('grids/x1.%d.grid.nc'%grid)
-mesh = Mesh(meshfile, llm, nqdyn, radius, f)
+pmesh = PMesh(comm,meshfile)
+mesh = Mesh(pmesh, llm, nqdyn, radius, f)
 print '...Done'
 lon, lat = mesh.lon_i, mesh.lat_i
@@ -30,4 +39,5 @@
 dx = mesh.de.min()
 dt = courant*dx/c0
+print dx, dt
 nt = int(math.ceil(T/dt))
 dt = T/nt
@@ -61,4 +71,6 @@
 flow=gh,u
 
+unst.ker.dynamico_update_halo(mesh.com_edges.index, 1, u.shape[0], u) 
+                                  
 for i in range(20):
     if True:
@@ -67,5 +79,5 @@
         print i, gh.shape, gh.min(), gh.max(), u.min(), u.max()
         mesh.plot_i(gh-Phi0) ; plt.title('err(gh)');
-        plt.savefig('fig_RSW2_MPAS_W02/err_gh_%02d.png'%i)
+        plt.savefig('%s_err_gh_%02d.png'%(prefix,i))
         plt.close()
     else:
@@ -100,5 +112,5 @@
         ax8.scatter(u_nowref,u_now-u_nowref) ; ax8.set_title('u_now');
         f.subplots_adjust(hspace=1.)
-        plt.savefig('fig_RSW2_MPAS_W02/scatter_%02d.png'%i)
+        plt.savefig('%s_scatter_%02d.png'%(prefix,i))
         plt.close()
 
Index: codes/icosagcm/devel/src/unstructured/data_unstructured.F90
===================================================================
--- codes/icosagcm/devel/src/unstructured/data_unstructured.F90	(revision 680)
+++ codes/icosagcm/devel/src/unstructured/data_unstructured.F90	(revision 681)
@@ -18,5 +18,5 @@
        thermo_theta=1, thermo_entropy=2, thermo_moist=3, thermo_boussinesq=4, &
        caldyn_vert_cons=1, max_nb_stage=5
-  INTEGER(C_INT),  BIND(C) :: caldyn_thermo=thermo_theta, caldyn_eta=eta_lag, &
+  INDEX,  BIND(C) :: caldyn_thermo=thermo_theta, caldyn_eta=eta_lag, &
        caldyn_vert_variant=caldyn_vert_cons, nb_threads=0, nb_stage=0
   LOGICAL(C_BOOL), BIND(C) :: hydrostatic=.TRUE., debug_hevi_solver=.TRUE.
@@ -54,4 +54,14 @@
        'pvort_only', 'slow_hydro', 'fast      ', 'coriolis  ', 'theta     ', 'geopot    ', 'vert      ', &
        'solver    ', 'slow_NH   ', 'NH_geopot ', 'vert_NH   ',  'update    ' /)
+
+  INTEGER, PARAMETER ::transfer_primal=1, transfer_edge=2, transfer_dual=3, transfer_max=3
+  TYPE Halo_transfer
+     INTEGER :: ranks ! size of arrays rank, len
+     INTEGER, ALLOCATABLE :: rank(:), & ! MPI ranks to communicate with
+          num(:), & ! number of cells to send to / receive from other MPI ranks
+          cells(:) ! local indices of cells to send/receive
+     DBL, ALLOCATABLE :: buf1(:), buf2(:,:), buf3(:,:,:)
+  END TYPE Halo_transfer
+  TYPE(Halo_transfer), TARGET :: send_info(transfer_max), recv_info(transfer_max)
 
 CONTAINS
Index: codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90
===================================================================
--- codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90	(revision 680)
+++ codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90	(revision 681)
@@ -6,4 +6,5 @@
 #endif
   USE caldyn_unstructured_mod
+  USE transfer_unstructured_mod
   IMPLICIT NONE
   PRIVATE
@@ -151,4 +152,7 @@
        DO stage=1, nb_stage
           
+          CALL update_halo(transfer_primal, rhodz)
+          CALL update_halo(transfer_primal, theta_rhodz(:,:,1)) ! FIXME
+
           IF(hydrostatic) THEN
              
@@ -163,6 +167,8 @@
              
              CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage), u)
-             
+             CALL update_halo(transfer_edge, u)
              CALL compute_pvort_only(rhodz,u,qv,qu)
+             CALL update_halo(transfer_edge, qu)
+
              CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow(:,:,stage))
              CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz(:,:,stage), dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage))
Index: codes/icosagcm/devel/src/unstructured/transfer_unstructured.F90
===================================================================
--- codes/icosagcm/devel/src/unstructured/transfer_unstructured.F90	(revision 681)
+++ codes/icosagcm/devel/src/unstructured/transfer_unstructured.F90	(revision 681)
@@ -0,0 +1,114 @@
+MODULE transfer_unstructured_mod
+  USE data_unstructured_mod
+  IMPLICIT NONE
+  SAVE
+  PRIVATE
+
+  INCLUDE 'mpif.h'
+  
+  PUBLIC :: update_halo
+  
+CONTAINS
+
+  SUBROUTINE init_halo_transfer(halo, rank, num, cells)
+    TYPE(Halo_transfer) :: halo
+    INTEGER :: ranks, total, rank(:), num(:), cells(:)
+    ranks=SIZE(rank)
+    total=SIZE(cells)
+    IF(total /= SUM(num) ) STOP 'init_halo_transfer : SIZE(cells)<>SUM(num)'
+    halo%ranks = ranks
+    ALLOCATE(halo%rank(ranks))
+    ALLOCATE(halo%num(ranks))
+    ALLOCATE(halo%cells(total))
+    ALLOCATE(halo%buf2(llm+1,total))
+    halo%rank=rank
+    halo%num=num
+    halo%cells=cells
+  END SUBROUTINE init_halo_transfer
+
+  SUBROUTINE init_transfer(index, &
+       send_num, send_size, send_rank, send_len, send_list, &
+       recv_num, recv_size, recv_rank, recv_len, recv_list) &
+       BIND(C, name='dynamico_init_transfer')
+    INTEGER(C_INT), VALUE :: index, send_num, send_size, recv_num, recv_size
+    INTEGER(C_INT) :: send_rank(send_num), send_len(send_num), send_list(send_size), &
+         recv_rank(recv_num), recv_len(recv_num), recv_list(recv_size)
+    PRINT *, 'init_transfer', index, &
+       send_num, send_size, '/', send_rank, '/', send_len, '/', &
+       recv_num, recv_size, '/', recv_rank, '/', recv_len, '/' 
+    CALL init_halo_transfer(send_info(index), send_rank, send_len, send_list)
+    CALL init_halo_transfer(recv_info(index), recv_rank, recv_len, recv_list)
+  END SUBROUTINE init_transfer
+
+  !-------------------------------------------------------------------------------------------------!
+  
+  SUBROUTINE copy_in(lldata, ndata, llbuf, nbuf, cells, data, buf)
+    INTEGER :: ndata, lldata, nbuf, llbuf, cells(nbuf), i
+    REAL(C_DOUBLE) :: data(lldata,ndata), buf(llbuf,nbuf)
+    !$OMP DO
+    DO i=1,nbuf
+       buf(1:lldata,i) = data(1:lldata, cells(i))
+    END DO
+    !$OMP END DO
+  END SUBROUTINE copy_in
+
+  SUBROUTINE copy_out(lldata, ndata, llbuf, nbuf, cells, data, buf)
+    INTEGER :: ndata, lldata, nbuf, llbuf, cells(nbuf), i
+    REAL(C_DOUBLE) :: data(lldata,ndata), buf(llbuf,nbuf)
+    !$OMP DO SCHEDULE(STATIC)
+    DO i=1,nbuf
+       data(1:lldata, cells(i)) = buf(1:lldata,i) 
+    END DO
+    !$OMP END DO
+  END SUBROUTINE copy_out
+
+  SUBROUTINE update_halo(index, data)
+    INTEGER :: index ! transfer_XXX, cf data_unstructured.F90
+    REAL(C_DOUBLE) :: data(:,:)    
+    TYPE(Halo_transfer), POINTER :: halo
+    REAL(C_DOUBLE), POINTER :: buf(:,:)
+    INTEGER :: req(send_info(index)%ranks+recv_info(index)%ranks), &
+         status(MPI_STATUS_SIZE, send_info(index)%ranks+recv_info(index)%ranks)
+    INTEGER :: llbuf, num, i, istart, ireq, ierr
+    halo => send_info(index)
+    buf => halo%buf2
+    llbuf = SIZE(buf,1)
+    CALL copy_in(SIZE(data,1), SIZE(data,2), llbuf, SIZE(buf,2), halo%cells, data, buf) 
+
+    !$OMP MASTER
+    ireq=0
+    istart=1
+    DO i = 1,halo%ranks
+       num = halo%num(i)
+       ireq = ireq+1
+       call MPI_ISEND(buf(1,istart), llbuf*num, MPI_REAL8, halo%rank(i), &
+            42, MPI_COMM_WORLD, req(ireq), ierr)
+       istart = istart+num
+    END DO
+    !$OMP END MASTER
+    halo => recv_info(index)
+    buf => halo%buf2
+    llbuf = SIZE(buf,1)
+    !$OMP MASTER
+    istart=1
+    DO i = 1,halo%ranks
+       num = halo%num(i)
+       ireq = ireq+1
+       call MPI_IRECV(buf(1,istart), llbuf*num, MPI_REAL8, halo%rank(i), &
+            42, MPI_COMM_WORLD, req(ireq), ierr)
+       istart = istart+num
+    END DO
+    CALL MPI_WAITALL(ireq, req, status, ierr)
+    !$OMP END MASTER
+    !$OMP BARRIER
+
+    CALL copy_out(SIZE(data,1), SIZE(data,2), SIZE(buf,1), SIZE(buf,2), halo%cells, data, buf)
+  END SUBROUTINE update_halo
+
+  SUBROUTINE update_halo_c(index, lldata, ndata, data) BIND(C, name='dynamico_update_halo')
+    INTEGER, VALUE :: index, lldata, ndata
+    REAL(C_DOUBLE) :: data(lldata, ndata)
+    CALL update_halo(index, data)
+  END SUBROUTINE update_halo_c
+  
+END MODULE transfer_unstructured_mod
