Index: codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py	(revision 787)
+++ codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py	(revision 790)
@@ -1,2 +1,21 @@
+from __future__ import print_function
+from dynamico import getargs
+
+getargs.add("--T", type=float, default=3600.,
+                    help="Length of time slice in seconds")
+getargs.add("--perturb", type=float, default=1.,
+                    help="Amplitude of wind perturbation in m/s")
+getargs.add("--N", type=int, default=48,
+                    help="Number of time slices")
+getargs.add("--Davies_N1", type=int, default=3)
+getargs.add("--Davies_N2", type=int, default=3)
+getargs.add("--nx", type=int, default=200)
+getargs.add("--ny", type=int, default=30)
+getargs.add("--betaplane", action='store_true')
+getargs.defaults(dt=360., llm=50)
+
+logging = getargs.Logger('main')
+args = getargs.parse()
+
 from dynamico import unstructured as unst
 from dynamico import dyn
@@ -11,32 +30,7 @@
 import numpy as np
 import time
-import argparse
-import logging
 from numpy import pi, log, exp, sin, cos
 
 # Baroclinic instability test based on Ullrich et al. 2015, QJRMS
-
-parser = argparse.ArgumentParser()
-
-parser.add_argument("--mpi_ni", type=int, default=1,
-                    help="number of x processors")
-parser.add_argument("--mpi_nj", type=int, default=1,
-                    help="number of y processors")
-parser.add_argument("--T", type=float, default=3600.,
-                    help="Length of time slice in seconds")
-parser.add_argument("--N", type=int, default=48,
-                    help="Number of time slices")
-parser.add_argument("--log",type=str,default='DEBUG',choices=('DEBUG', 'INFO', 'WARNING', 'ERROR', 'CRITICAL'),
-                    help="Set verbosity level of logging")
-parser.add_argument("--Davies_N1", type=int, default=3)
-parser.add_argument("--Davies_N2", type=int, default=3)
-parser.add_argument("--llm", type=int, default=50)
-parser.add_argument("--nx", type=int, default=200)
-parser.add_argument("--ny", type=int, default=30)
-parser.add_argument("--dt", type=float, default=360., help='Time step in seconds')
-parser.add_argument("--hydrostatic", action='store_true')
-parser.add_argument("--betaplane", action='store_true')
-
-args = parser.parse_args()
 
 
@@ -52,11 +46,11 @@
     Cpd   = 1004.5
     p0    = 1e5
-    up    = 1. # amplitude (m/s)
+    up    = args.perturb # amplitude (m/s)
     xc,yc,lp = 0.,Ly/2.,600000.
 
     omega  = 7.292e-5  # Angular velocity of the Earth (s^-1)    
-    phi0   = 45.*np.pi/180.0 # Reference latitude North pi/4 (deg)
-    f0     = 2*omega*np.sin(phi0) 
-    beta0  = 2*omega*np.cos(phi0)/a if args.betaplane else 0.
+    phi0   = 45.*pi/180.0 # Reference latitude North pi/4 (deg)
+    f0     = 2*omega*sin(phi0) 
+    beta0  = 2*omega*cos(phi0)/a if args.betaplane else 0.
     fb     = f0 - beta0*y0
 
@@ -69,5 +63,5 @@
     def ulon(x,y,eta):
         u = -u0*(sin(pi*y/Ly)**2)*log(eta)*(eta**(-log(eta)/(b*b))) 
-#        u = u + up*exp(-(((x-xc)**2+(y-yc)**2)/lp**2))
+        u = u + up*exp(-(((x-xc)**2+(y-yc)**2)/lp**2))
         return  u
 
@@ -103,5 +97,5 @@
     eta_ik = eta(alpha_ik)
     eta_ek = eta(alpha_ek)
-    print('min max eta_il', np.min(eta_il),np.max(eta_il))
+#    print('min max eta_il', np.min(eta_il),np.max(eta_il))
 
     Phi_il = Phi_xyeta(y_il, eta_il)
@@ -170,4 +164,17 @@
     return div_ik
 
+def KE(mesh, u_ek):
+    llm=u_ek.shape[1]
+    primal_num, primal_deg, primal_edge = mesh.primal_num, mesh.primal_deg, mesh.primal_edge
+    primal_ne, le_de, Ai = mesh.primal_ne, mesh.le_de, mesh.Ai
+    div_ik = prec.zeros((primal_num, llm))
+    for i in range(primal_num):
+        div_i = 0.
+        for iedge in range(primal_deg[i]):
+            e = primal_edge[i,iedge]-1 # Fortran => Python indexing
+            div_i = div_i + le_de[e]*(u_ek[e,:]**2)
+        div_ik[i,:] = div_i * (.25/ Ai[i])
+    return div_ik
+
 def curl(mesh, u_ek):
     llm=u_ek.shape[1]
@@ -194,5 +201,5 @@
         for i in range(c.size):
             ci=c[i]
-            m[i] = (1.+np.cos((ci-c2)*np.pi/(N2*delta)))/2.0
+            m[i] = (1.+cos((ci-c2)*pi/(N2*delta)))/2.0
             if ci < c2 : m[i]=1.
             if ci > c1 : m[i]=0.
@@ -211,20 +218,14 @@
 class myDavies(Davies):
     def mask(self,x,y):
-        logging.debug('davies dy = %f' % dy)
+        logging.debug('davies dy = %f'% dy)
         return self.mask0(y,.5*Ly,dy)*self.mask0(-y,.5*Ly,dy)
-
-
         
 with xios.Client() as client: # setup XIOS which creates the DYNAMICO communicator
     comm = client.comm
     mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size()
-
-    #-------------Logging verbosity configuration---------------------------
-    myloglevel = args.log
-    myloglevel = getattr(logging, myloglevel.upper())
-    logging.basicConfig(filename='out.log',filemode='w',level=myloglevel, 
-        format='%(processName)s:%(name)s:%(filename)s:%(module)s:%(funcName)s:%(lineno)d:%(levelname)s:%(message)s' )
-
-    logging.debug('%d/%d starting'%(mpi_rank,mpi_size))
+#    getargs.config_log(filename='out.log',filemode='w',
+#                    format='%(processName)s:%(name)s:%(filename)s:%(module)s:%(funcName)s:%(lineno)d:%(levelname)s:%(message)s' )
+
+    logging.info('%d/%d starting'%(mpi_rank,mpi_size))
 
     g, Lx, Ly = 9.81, 4e7, 6e6
@@ -237,5 +238,4 @@
 
     mass_bl,mass_dak,mass_dbk = meshes.compute_hybrid_coefs(flow0[0])
-    print( 'Type of mass_bl, mass_dak, mass_dbk : ', [x.dtype for x in mass_bl, mass_dak, mass_dbk])
     unst.ker.dynamico_init_hybrid(mass_bl,mass_dak,mass_dbk)
 
@@ -247,6 +247,4 @@
     logging.info( 'Time step : %d x %g = %g s' % (nt,dt,nt*dt))
     
-    dt=0. # FIXME
-
     caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass
 
@@ -264,13 +262,13 @@
             caldyn_step = unst.caldyn_step_NH(mesh,scheme,1, caldyn_thermo,caldyn_eta,thermo,params,params.g)
         davies = myDavies(args.Davies_N1, args.Davies_N2, mesh.lon_i, mesh.lat_i, mesh.lon_e,mesh.lat_e)
-        print('mask min/max :', davies.beta_i.min(), davies.beta_i.max() )
+#        print('mask min/max :', davies.beta_i.min(), davies.beta_i.max() )
         logging.debug('mask min/max :')
-        logging.debug(davies.beta_i.min())
-        logging.debug(davies.beta_i.max())
+        logging.debug('%f'% davies.beta_i.min())
+        logging.debug('%f'% davies.beta_i.max())
         def next_flow(m,S,u,Phi,W):
             # junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
             caldyn_step.mass[:,:], caldyn_step.theta_rhodz[:,:], caldyn_step.u[:,:] = m,S,u
             caldyn_step.geopot[:,:], caldyn_step.W[:,:] = Phi,W
-            kappa = dx**4/43200.
+            kappa = dx**4/432000.
             for i in range(nt):
                 caldyn_step.next()
@@ -282,8 +280,4 @@
                 unst.ker.dynamico_scalar_laplacian(s,laps)
                 unst.ker.dynamico_scalar_laplacian(laps,bilaps)
-#                grads = grad(mesh,s)
-#                laps = div(mesh,grads)
-#                gradlaps = grad(mesh,laps) # bilaplacian
-#                bilaps = div(mesh,gradlaps)                
                 caldyn_step.theta_rhodz[:] = S - dt*kappa*bilaps*m # Euler step
             return (caldyn_step.mass.copy(), caldyn_step.theta_rhodz.copy(), caldyn_step.u.copy(),
@@ -293,5 +287,5 @@
     if caldyn_thermo == unst.thermo_theta:
         s=S/m
-        theta = thermo.T0*np.exp(s/thermo.Cpd)
+        theta = thermo.T0*exp(s/thermo.Cpd)
         S=m*theta
         title_format = 'Potential temperature at t=%g s (K)'
@@ -316,16 +310,25 @@
             ps = caldyn_step.ps
             curl_vk = curl(mesh, u)
+            KE_ik = KE(mesh,u)
             du_fast, du_slow = caldyn_step.du_fast[0,:,:], caldyn_step.du_slow[0,:,:] 
             div_fast, div_slow = div(mesh,du_fast), div(mesh,du_slow)
+            drhodz, hflux = caldyn_step.drhodz[0,:,:], caldyn_step.hflux[:,:]
+
             # write to disk
             context.send_field_primal('ps', ps)
             context.send_field_primal('temp', gas.T)
+
             context.send_field_primal('p', gas.p)
+#            context.send_field_primal('p', KE_ik)
+#            context.send_field_primal('p', drhodz)
+
             context.send_field_primal('theta', gas.s)
             context.send_field_primal('uz', w)
             context.send_field_primal('z', z)
             context.send_field_primal('div_fast', div_fast)
+
             context.send_field_primal('div_slow', div_slow)
-            print curl_vk.size
+#            context.send_field_primal('div_slow', div(mesh,hflux) )
+
             context.send_field_dual('curl', curl_vk)
 
Index: codes/icosagcm/devel/Python/test/py/test_halo.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/test_halo.py	(revision 790)
+++ codes/icosagcm/devel/Python/test/py/test_halo.py	(revision 790)
@@ -0,0 +1,230 @@
+from dynamico import meshes
+from dynamico import unstructured as unst
+from dynamico import time_step
+from dynamico import precision as prec
+from dynamico import xios
+
+import math as math
+import numpy as np
+import argparse
+
+#------------------ Command-line arguments ------------------
+
+parser = argparse.ArgumentParser()
+
+parser.add_argument("--mpi_ni", type=int, default=1,
+                    help="number of x processors")
+parser.add_argument("--mpi_nj", type=int, default=1,
+                    help="number of y processors")
+parser.add_argument("-T", type=float, default=10.,
+                    help="Length of time slice")
+args = parser.parse_args()
+
+#------------------------- Functions ----------------------
+
+def globmin(x):
+   xloc = x[mesh.primal_own_loc]
+   locmin= xloc.min()
+   return client.min(locmin)
+
+def globmax(x):
+   xloc = x[mesh.primal_own_loc]
+   locmin= xloc.max()
+   return client.max(locmin)
+
+def minmax(name, x):
+   gmin,gmax  = globmin(x),globmax(x)
+   if(mpi_rank==0) : print('Min/max %s :'%name, gmin, gmax ) 
+
+def laplace(mesh, p):
+   left, right, ne = mesh.left, mesh.right, mesh.primal_ne
+   gradp = mesh.field_u().flatten()
+   for ij in range(mesh.left.size):
+      p_left = p[left[ij]-1] # left,right start at 1 (CHECK)
+      p_right = p[right[ij]-1] # left,right start at 1 (CHECK)
+      gradp[ij]= mesh.le_de[ij]*(p_right-p_left) # convert to contravariant
+   lap = mesh.field_mass().flatten()
+   for ij in range(lap.size):
+      lap_ij=0.
+      for k in range(mesh.primal_deg[ij]):
+         edge=mesh.primal_edge[ij,k]-1 # primal_edge starts at 1 (CHECK)
+         lap_ij = lap_ij + mesh.primal_ne[ij,k]*gradp[edge]
+      lap[ij]=lap_ij
+   return lap/mesh.Ai
+
+def curlgrad(mesh,p):
+   left, right, ne = mesh.left, mesh.right, mesh.primal_ne
+   gradp = mesh.field_u().flatten()
+   for ij in range(mesh.left.size):
+      p_left = p[left[ij]-1] # left,right start at 1 (CHECK)
+      p_right = p[right[ij]-1] # left,right start at 1 (CHECK)
+      gradp[ij]= p_right-p_left # covariant
+   curl = mesh.field_z().flatten()
+   for ij in range(curl.size):
+      curl_ij=0.
+      for k in range(mesh.dual_deg[ij]):
+         edge=mesh.dual_edge[ij,k]-1 # dual_edge starts at 1 (CHECK)
+         curl_ij = curl_ij + mesh.dual_ne[ij,k]*gradp[edge]
+      curl[ij]=curl_ij
+   return curl/mesh.Av
+
+def curl(mesh, ucov):
+   curl = mesh.field_z().flatten()
+   for ij in range(curl.size):
+      curl_ij=0.
+      for k in range(mesh.dual_deg[ij]):
+         edge=mesh.dual_edge[ij,k]-1 # dual_edge starts at 1 (CHECK)
+         curl_ij = curl_ij + mesh.dual_ne[ij,k]*ucov[edge]
+      curl[ij]=curl_ij
+   return curl/mesh.Av
+
+def test_laplace_seq():
+   p=p.flatten()
+   print p.shape
+   curl = curlgrad(mesh,p)
+   print 'curl(grad(p)) min max :', curl.min(), curl.max()       
+   lap = laplace(mesh,p)
+   p  = p+lap/vp
+   print 'p+lap(p)/vp min max :', p.min(), p.max()
+
+def test_laplace_parallel():
+   transfer_primal, transfer_edge = mesh.com_primal.index, mesh.com_edges.index
+   with xios.Context_Curvilinear(mesh,1, 24*3600) as context:
+      # now XIOS knows about the mesh and we can write to disk
+      context.update_calendar(0)
+      for it in range(4):
+         context.update_calendar(it+1)
+         unst.ker.dynamico_update_halo(mesh.com_primal.index, 1, mesh.primal_num, p)
+         lap = laplace(mesh,p)
+          # compute eigenvalue and save residual (should be 0.)
+         context.send_field_primal('p', p)
+         p  = p+lap/vp
+         context.send_field_primal('ps', p)
+          # next iteration will start from lap(p)
+         p=-lap/vp
+
+def periodize(z,nz): return (z+2*nz)%nz
+def index(x,y):      return periodize(x,nx)+nx*periodize(y,ny)
+def indexu(x,y):     return 2*index(x,y)
+def indexv(x,y):     return 2*index(x,y)+1
+
+def test_RSW():
+   x1,x2= xx-1., xx+1.
+   if False:
+      u0 = mesh.field_u()  # flow initally at rest
+      h0 = 1+0.1*(np.exp(-2.*(x1*x1+yy*yy))+
+               np.exp(-2.*(x2*x2+yy*yy)))
+   else: # uniform flow
+      ulon = 1.+0.*mesh.lat_e
+      ulat = 0.*ulon
+      u0 = mesh.ucov2D(ulon,ulat)
+      h0 = mesh.field_mass()
+      h0[:]=1.
+
+   dt=0.1
+   scheme = time_step.RK1(None, dt)
+   step = unst.caldyn_step_TRSW(mesh,scheme,1)
+
+   # check halo exchange
+   step.u[:] = mesh.edges_E2 # global index of edges
+   unst.ker.dynamico_update_halo(mesh.com_edges.index, 1, mesh.edge_num, step.u)
+   for i in range(step.u.size):
+      if step.u[i] != mesh.edges_E2[i]:
+         print '[%d] minmax halo exchanges mismatch : ', i, step.u[i], mesh.edges_E2[i]
+
+   step.mass[:]=h0
+   step.theta_rhodz[:]=h0
+   step.u[:]=u0
+
+   def minmax(name, data): print '[%d] minmax %s :'%(mpi_rank,name), data.min(), data.max()
+
+   print 'minmax dual_deg : ', mesh.dual_deg.min(), mesh.dual_deg.max()
+   with xios.Context_Curvilinear(mesh,1, 24*3600) as context:
+      for it in range(3):
+         step.next()
+         context.update_calendar(it+1)
+         u, dh, fast, slow = step.u[:], step.drhodz[:], step.du_fast[:], step.du_slow[:]
+         qv = step.qv[:]
+         # better exchange halos before evaluating min/max
+         unst.ker.dynamico_update_halo(mesh.com_primal.index, 1, mesh.primal_num, dh)
+         unst.ker.dynamico_update_halo(mesh.com_edges.index, 1, mesh.edge_num, fast)
+         unst.ker.dynamico_update_halo(mesh.com_edges.index, 1, mesh.edge_num, slow)
+#         unst.ker.dynamico_update_halo(mesh.com_edges.index, 1, mesh.edge_num, qu)
+#         unst.ker.dynamico_update_halo(mesh.com_dual.index, 1, mesh.dual_num, qv)
+         minmax('dh',dh)
+         minmax('du_fast',fast)
+         minmax('du_slow',slow)
+#         minmax('qu',qu)
+         minmax('qv',qv)
+         context.send_field_primal('p', dh )
+
+         unst.ker.dynamico_update_halo(mesh.com_edges.index, 1, mesh.edge_num, u)
+         zeta = curl(mesh, u)
+         minmax('zeta', zeta)
+         for i in range(zeta.size):
+            if zeta[i]<0. :
+               if mpi_rank==0:
+                  ij = mesh.vertices_V1[i] # global index of dual cell (starting at 0)
+                  x,y = ij%nx, ij/nx
+                  edges1 = [ indexv(x+1,y), indexu(x,y+1), indexv(x,y), indexu(x,y) ]
+                  edges2 = [ mesh.edges_E2[k-1] for k in mesh.dual_edge[i,:] ]
+                  if edges1!=edges2: print '[%d] minmax zeta '%mpi_rank, ij, x, y
+                  print '[%d] minmax zeta '%mpi_rank, [ u[edge-1] for edge in mesh.dual_edge[i,:] ]
+                  print '[%d] minmax zeta '%mpi_rank, mesh.dual_ne[i,:]
+
+#--------------------------------- Main program -----------------------------
+
+with xios.Client() as client: # setup XIOS which creates the DYNAMICO communicator
+    comm = client.comm
+    mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size()
+    print '%d/%d starting'%(mpi_rank,mpi_size)
+
+    #    nx,ny,Lx,Ly=200,30,4e7,6e6
+    nx,ny,Lx,Ly      = 32,32,8.,8.
+    llm, nqdyn, g, f = 1,1,1.,0.
+    dx,dy=Lx/nx,Ly/ny
+
+    cfl = .8
+    dt = cfl/math.sqrt((nx/Lx)**2+(ny/Ly)**2)
+    T=args.T
+    N=int(T/dt)+1
+    dt=T/N
+    print N,dt,Lx/nx
+
+    filename, radius = 'cart_%03d_%03d.nc'%(nx,ny), None
+    unst.setvar('g',g)
+
+    def coriolis(lon,lat):
+       return f+0.*lon
+
+    print 'Reading Cartesian mesh ...'
+    meshfile = meshes.DYNAMICO_Format(filename)
+
+    parallel=True
+    laplace=False
+
+    if parallel:
+       pmesh = meshes.Unstructured_PMesh(comm,meshfile)
+       pmesh.partition_curvilinear(args.mpi_ni,args.mpi_nj)
+    #    pmesh.partition_metis()
+       mesh  = meshes.Local_Mesh(pmesh, llm, nqdyn, radius, coriolis)
+    else:
+       mesh = meshes.Cartesian_mesh(nx,ny,llm,nqdyn,Lx,Ly,f)
+
+    xx,yy = mesh.lon_i, mesh.lat_i
+
+    if laplace:
+       p = np.cos(2*np.pi*xx/Lx)*np.cos(2.*np.pi*yy/Ly)
+
+    # analytic eigenvalue : vp = -((2*np/pi/Lx)**2+(2*np.pi/Ly)**2)
+       kx = 2.*np.sin(dx*np.pi/Lx)/dx
+       ky = 2.*np.sin(dy*np.pi/Ly)/dy
+       vp = kx**2+ky**2
+
+       if parallel:
+          test_laplace_parallel()
+       else:
+          test_laplace_seq()
+
+    else:
+       test_RSW()
Index: codes/icosagcm/devel/Python/test/py/write_Cartesian_mesh.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/write_Cartesian_mesh.py	(revision 787)
+++ codes/icosagcm/devel/Python/test/py/write_Cartesian_mesh.py	(revision 790)
@@ -1,2 +1,18 @@
+from dynamico import getargs
+getargs.add("--nx", type=int,
+                    default=64, choices=None,
+                    help="number of x points")
+getargs.add("--ny", type=int,
+                    default=64, choices=None,
+                    help="number of y points")
+getargs.add("--Lx", type=float,
+                    default=8., choices=None,
+                    help="Lx")
+getargs.add("--Ly", type=float,
+                    default=8., choices=None,
+                    help="Ly")
+
+args = getargs.parse()
+
 from dynamico.meshes import zeros
 from dynamico import meshes
@@ -4,20 +20,4 @@
 import netCDF4 as cdf
 import argparse
-
-parser = argparse.ArgumentParser()
-
-parser.add_argument("-nx", type=int,
-                    default=64, choices=None,
-                    help="number of x points")
-parser.add_argument("-ny", type=int,
-                    default=64, choices=None,
-                    help="number of y points")
-parser.add_argument("-Lx", type=float,
-                    default=8., choices=None,
-                    help="Lx")
-parser.add_argument("-Ly", type=float,
-                    default=8., choices=None,
-                    help="Ly")
-args = parser.parse_args()
 
 nx, ny, Lx, Ly, llm, nqdyn = args.nx, args.ny,args.Lx, args.Ly, 1, 1
Index: codes/icosagcm/devel/Python/test/python.sh
===================================================================
--- codes/icosagcm/devel/Python/test/python.sh	(revision 787)
+++ codes/icosagcm/devel/Python/test/python.sh	(revision 790)
@@ -20,5 +20,6 @@
     NB_MPI=$1
     shift
-    mpirun -np $NB_MPI python -u $*
+    rm -f xios_client*.*
+    mpirun -np $NB_MPI python -u $* | tee dynamico.log
 }
 
