Index: codes/icosagcm/devel/Python/dynamico/meshes.py
===================================================================
--- codes/icosagcm/devel/Python/dynamico/meshes.py	(revision 786)
+++ codes/icosagcm/devel/Python/dynamico/meshes.py	(revision 787)
@@ -478,4 +478,5 @@
         if gridfile.meshtype == 'curvilinear':
             self.primal_i, self.primal_j = get_pvars(dim_primal, 'primal_i','primal_j')
+            self.dual_i, self.dual_j = get_pvars(dim_dual, 'primal_i','primal_j')
             self.edge_i, self.edge_j = get_pvars(dim_edge, 'edge_i','edge_j')
 
@@ -494,12 +495,12 @@
                          'lon_v', 'lat_v', 'Av', 'lon_i', 'lat_i', 'Ai', 
                          'le', 'de', 'lon_e', 'lat_e', 'angle_e')
-    def partition(self, edge_owner):
+    def partition_metis(self):
+        print 'Partitioning with ParMetis...'
+        edge_owner = unst.partition_mesh(self.trisk_deg, self.trisk, self.comm.Get_size())
         self.edge_owner = parallel.LocPArray1D(self.dim_edge, edge_owner)
         primal_owner = partition_from_stencil(self.edge_owner, self.primal_deg, self.primal_edge)
         self.primal_owner = parallel.LocPArray1D(self.dim_primal, primal_owner)
-
-    def partition_metis(self):
-        print 'Partitioning with ParMetis...'
-        self.partition(unst.partition_mesh(self.trisk_deg, self.trisk, self.comm.Get_size()))
+        dual_owner = partition_from_stencil(self.edge_owner, self.dual_deg, self.dual_edge)
+        self.dual_owner = parallel.LocPArray1D(self.dim_dual, dual_owner)
 
     def partition_curvilinear(self, ni,nj):
@@ -513,4 +514,5 @@
         self.edge_owner   = owner(self.dim_edge,   self.edge_i,   self.edge_j)
         self.primal_owner = owner(self.dim_primal, self.primal_i, self.primal_j)
+        self.dual_owner = self.primal_owner
         print 'partition_curvilinear %d :'%(self.comm.Get_rank()), self.primal_owner.data
 
@@ -541,6 +543,6 @@
 class Local_Mesh(Abstract_Mesh):
     def __init__(self, pmesh, llm, nqdyn, radius, f):
-        comm, dim_primal, primal_owner, dim_edge, edge_owner, dim_dual = pmesh.members(
-            'comm', 'dim_primal', 'primal_owner', 'dim_edge', 'edge_owner', 'dim_dual')
+        comm, dim_primal, primal_owner, dim_edge, edge_owner, dim_dual, dual_owner = pmesh.members(
+            'comm', 'dim_primal', 'primal_owner', 'dim_edge', 'edge_owner', 'dim_dual', 'dual_owner' )
         print 'Constructing halos ...'
         edges_E0 = find_my_cells(edge_owner)
@@ -564,4 +566,5 @@
             self.nx, self.ny = pmesh.gridfile.nx, pmesh.gridfile.ny
             self.primal_i, self.primal_j = pmesh.get(get_all_cells, 'primal_i', 'primal_j')
+            self.dual_i, self.dual_j = pmesh.get(get_all_duals, 'dual_i', 'dual_j')
 
 
@@ -599,6 +602,10 @@
         primal_own_all = comm.allgather(primal_own_num)
         displ = sum(primal_own_all[0:comm.Get_rank()]) # required by XIOS
-
         print 'local primal mesh :', primal_vertex.min(), primal_vertex.max(), primal_own_all, displ
+
+        # construct own dual mesh for XIOS output
+        dual_own_glo = find_my_cells(dual_owner)
+        dual_own_glo = np.asarray(primal_own_glo, dtype=np.int32)
+        dual_own_loc = [dict_V1[i] for i in dual_own_glo]
         
         # compute_ne(...) and normalize(...) expect indices starting at 1
@@ -610,5 +617,5 @@
         self.set_members(locals(), 'nqdyn', 'llm', 'primal_num', 'dual_num', 'edge_num',
                          'primal_own_glo', 'primal_own_loc', 'primal_own_deg',
-                         'primal_deg', 'primal_vertex', 'displ',
+                         'primal_deg', 'primal_vertex', 'displ', 'dual_own_loc', 
                          'trisk_deg', 'trisk', 'wee', 
                          'dual_deg', 'dual_edge', 'dual_ne', 'dual_vertex', 'Riv2',
Index: codes/icosagcm/devel/Python/dynamico/xios.py
===================================================================
--- codes/icosagcm/devel/Python/dynamico/xios.py	(revision 786)
+++ codes/icosagcm/devel/Python/dynamico/xios.py	(revision 787)
@@ -4,44 +4,4 @@
 from dynamico.meshes import radian
 from mpi4py import MPI
-
-#-----------------------------------------------------------------------------
-
-def setup_curvilinear(cat,id,nx,ny,displ,index_own_i,index_own_j,lon,lat):
-    mesh = cxios.Handle(cat,id)
-    def log(name,data) : print 'setup_curvilinear : %s shape min max'%name, data.shape, data.min(), data.max()
-    log('index_own_i',index_own_i)
-    log('index_own_j',index_own_j)
-    log('lon',lon)
-    log('lat',lat)
-    mesh.set_attr(type='curvilinear',
-                  ni_glo=nx, i_index=index_own_i,
-                  nj_glo=ny, j_index=index_own_j)
-    mesh.set_attr(lonvalue_1d=lon.flatten(), latvalue_1d=lat.flatten())
-
-def cf_boundaries(degree,points,lon,lat):
-    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)
 
 #-----------------------------------------------------------------------------
@@ -55,5 +15,7 @@
         print '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):
@@ -84,5 +46,8 @@
         ker.dynamico_xios_update_calendar(i)
     def send_field_primal(self,name, data):
-        own_loc = self.mesh.primal_own_loc
+        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]
@@ -93,15 +58,56 @@
         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) : print 'setup_curvilinear : %s shape min max'%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
-        own = mesh.primal_own_loc
-        lon_i, lat_i = [ x[own] for x in mesh.lon_i, mesh.lat_i ]
-        primal_i, primal_j = [ x[own] for x in mesh.primal_i, mesh.primal_j ]
-        setup_curvilinear('domaingroup','i', mesh.nx, mesh.ny, mesh.displ,
-                          primal_i, primal_j, lon_i, lat_i)
+        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):
Index: codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py	(revision 786)
+++ codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py	(revision 787)
@@ -31,5 +31,10 @@
 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()
@@ -44,5 +49,5 @@
     a     = 6.371229e6 # Radius of the Earth (m)
     ptop  = 2000.
-    y0    = Ly*0.5
+    y0    = .5*Ly
     Cpd   = 1004.5
     p0    = 1e5
@@ -51,8 +56,8 @@
 
     omega  = 7.292e-5  # Angular velocity of the Earth (s^-1)    
-    phi0   = 90.*np.pi/180.0 # Reference latitude North pi/4 (deg)
+    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
-    fb     = 2*omega*np.sin(phi0) - y0*2*omega*np.cos(phi0)/a
+    beta0  = 2*omega*np.cos(phi0)/a if args.betaplane else 0.
+    fb     = f0 - beta0*y0
 
     def Phi_xy(y):
@@ -64,5 +69,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
 
@@ -71,6 +76,6 @@
     def p(eta): return p0*eta     # eta  = p/p_s
 
-    def eta(alpha) : return (1-(lap*ztop*alpha/(T0)))**(g/(Rd*lap)) # roughly equispaced levels
-    def coriolis(x,y): return f0+beta0*y
+    def eta(alpha) : return (1.-(lap*ztop*alpha/T0))**(g/(Rd*lap)) # roughly equispaced levels
+    def coriolis(x,y): return f0+beta0*y # here y is 0 at domain center
 
     filename = 'cart_%03d_%03d.nc'%(nx,ny)
@@ -82,13 +87,12 @@
     mesh  = meshes.Local_Mesh(pmesh, llm, nqdyn, radius, coriolis)
 
-
     alpha_k = (np.arange(llm)  +.5)/llm
     alpha_l = (np.arange(llm+1)+ 0.)/llm
     x_ik, alpha_ik = np.meshgrid(mesh.lon_i,       alpha_k, indexing='ij')
-    y_ik, alpha_ik = np.meshgrid(mesh.lat_i+.5*Ly, alpha_k, indexing='ij')
+    y_ik, alpha_ik = np.meshgrid(mesh.lat_i+y0, alpha_k, indexing='ij')
     x_il, alpha_il = np.meshgrid(mesh.lon_i,       alpha_l, indexing='ij')
-    y_il, alpha_il = np.meshgrid(mesh.lat_i+.5*Ly, alpha_l, indexing='ij')
+    y_il, alpha_il = np.meshgrid(mesh.lat_i+y0, alpha_l, indexing='ij')
     x_ek, alpha_ek = np.meshgrid(mesh.lon_e,       alpha_k, indexing='ij')
-    y_ek, alpha_ek = np.meshgrid(mesh.lat_e+.5*Ly, alpha_k, indexing='ij')
+    y_ek, alpha_ek = np.meshgrid(mesh.lat_e+y0, alpha_k, indexing='ij')
 
     print('ztop(ptop) according to Eq. 7:', T0/lap*(1.-(ptop/p0)**(Rd*lap/g))) 
@@ -103,5 +107,5 @@
     Phi_il = Phi_xyeta(y_il, eta_il)
     Phi_ik = Phi_xyeta(y_ik, eta_ik)
-    p_ik = p(eta_ik)
+    p_ik, p_il = p(eta_ik), p(eta_il)
     T_ik = T(y_ik, eta_ik)  #ik full level(40), il(41)
 
@@ -110,4 +114,5 @@
     for l in range(llm): 
         mass_ik[:,l]=(Phi_il[:,l+1]-Phi_il[:,l])/(g*gas.v[:,l])
+#        mass_ik[:,l]=(p_il[:,l]-p_il[:,l+1])/g
     Sik, Wil  = gas.s*mass_ik, mesh.field_w()
     
@@ -224,5 +229,5 @@
 
     g, Lx, Ly = 9.81, 4e7, 6e6
-    nx, ny, llm = 200, 30, 30
+    nx, ny, llm = args.nx, args.ny, args.llm
     dx,dy = Lx/nx,Ly/ny
 
@@ -236,5 +241,5 @@
 
     T  = args.T
-    dt = 100.#180.#360.
+    dt = args.dt
     dz = flow0[3].max()/(params.g*llm)
     nt = int(math.ceil(T/dt))
@@ -242,4 +247,6 @@
     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
 
@@ -252,5 +259,8 @@
     else: # time stepping in Fortran
         scheme = time_step.ARK2(None, dt)
-        caldyn_step = unst.caldyn_step_NH(mesh,scheme,1, caldyn_thermo,caldyn_eta,thermo,params,params.g)
+        if args.hydrostatic:
+            caldyn_step = unst.caldyn_step_HPE(mesh,scheme,1, caldyn_thermo,caldyn_eta,thermo,params,params.g)
+        else:
+            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() )
@@ -292,5 +302,4 @@
     z=mesh.field_mass()
 
-
     #T, nslice, dt = 3600., 1, 3600.
     Nslice = args.N
@@ -305,7 +314,10 @@
             # Diagnose quantities of interest from prognostic variables m,S,u,Phi,W
             gas, w, z = diagnose(Phi,S,m,W)
+            ps = caldyn_step.ps
             curl_vk = curl(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)
             # write to disk
-            context.send_field_primal('ps', davies.beta_i)
+            context.send_field_primal('ps', ps)
             context.send_field_primal('temp', gas.T)
             context.send_field_primal('p', gas.p)
@@ -313,5 +325,8 @@
             context.send_field_primal('uz', w)
             context.send_field_primal('z', z)
-            context.send_field_primal('curl', curl_vk)
+            context.send_field_primal('div_fast', div_fast)
+            context.send_field_primal('div_slow', div_slow)
+            print curl_vk.size
+            context.send_field_dual('curl', curl_vk)
 
             print( 'ptop, model top (m) :', unst.getvar('ptop'), Phi.max()/unst.getvar('g'))
Index: codes/icosagcm/devel/Python/test/xml/context_icosa_simple.xml
===================================================================
--- codes/icosagcm/devel/Python/test/xml/context_icosa_simple.xml	(revision 786)
+++ codes/icosagcm/devel/Python/test/xml/context_icosa_simple.xml	(revision 787)
@@ -20,5 +20,9 @@
 
      <domain_group id="i">
-       <domain id="i" name="mesh"/>
+       <domain id="i" name="primal"/>
+     </domain_group>
+
+     <domain_group id="v">
+       <domain id="v" name="dual"/>
      </domain_group>
 
Index: codes/icosagcm/devel/Python/test/xml/field_def_simple.xml
===================================================================
--- codes/icosagcm/devel/Python/test/xml/field_def_simple.xml	(revision 786)
+++ codes/icosagcm/devel/Python/test/xml/field_def_simple.xml	(revision 787)
@@ -9,4 +9,10 @@
     <field id="mid_ap" axis_ref="lev"   long_name="hybrid A coefficient at midpoints" />
     <field id="mid_bp" axis_ref="lev"  long_name="hybrid B coefficient at midpoints" />
+
+    <field_group domain_ref="v">
+      <field_group axis_ref="lev"> 
+        <field id="curl" />
+      </field_group>
+    </field_group>
     
     <field_group domain_ref="i">
@@ -23,5 +29,6 @@
         <field id="z" />
         <field id="theta" />
-        <field id="curl" />
+        <field id="div_fast" />
+        <field id="div_slow" />
       </field_group>
      
Index: codes/icosagcm/devel/Python/test/xml/filedef_simple.xml
===================================================================
--- codes/icosagcm/devel/Python/test/xml/filedef_simple.xml	(revision 786)
+++ codes/icosagcm/devel/Python/test/xml/filedef_simple.xml	(revision 787)
@@ -21,4 +21,6 @@
       <field field_ref="theta" />
       <field field_ref="curl" />
+      <field field_ref="div_fast" />
+      <field field_ref="div_slow" />
       <field field_ref="phis" operation="once" freq_offset="0ts" name="PHIS"   standard_name="surface_geopotential" long_name="Surface geopotential" unit="m2/m2"/>      
     </field_group>
