Index: codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py	(revision 794)
+++ codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py	(revision 795)
@@ -1,3 +1,5 @@
 from __future__ import print_function
+from __future__ import division
+
 from dynamico import getargs
 
@@ -13,5 +15,13 @@
 getargs.add("--ny", type=int, default=30)
 getargs.add("--betaplane", action='store_true')
+
+getargs.add("--kappa_divgrad",  type=float, default=3.0e15)
+getargs.add("--kappa_curlcurl", type=float, default=3.0e15)
+
 getargs.defaults(dt=360., llm=50)
+
+getargs.config_log(filename='out.log',filemode='w') # must be done before calling Logger()
+#    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 = getargs.Logger('main')
@@ -136,12 +146,21 @@
     return thermo, mesh, params, prec.asnum([mass_ik,Sik,u_ek,Phi_il,Wil]), gas
 
-def diagnose(Phi,S,m,W):
+def diagnose(Phi,S,m,W,u):
     s=S/m
+    curl_vk = curl(mesh, u)
+    abs_vort_vk = mesh.field_z() # absolute vorticity
+    un = mesh.field_u()
+    v = mesh.field_mass() # specific volume
+    w = mesh.field_mass()
+    z = mesh.field_mass()
     for l in range(llm):
         v[:,l]=(Phi[:,l+1]-Phi[:,l])/(g*m[:,l])
         w[:,l]=.5*params.g*(W[:,l+1]+W[:,l])/m[:,l]
         z[:,l]=.5*(Phi[:,l+1]+Phi[:,l])/params.g
+        un[:,l]=u[:,l]/mesh.de
+        abs_vort_vk[:,l]=curl_vk[:,l] + mesh.fv
     gas = thermo.set_vs(v,s)
-    return gas, w, z
+    ps = gas.p[:,0]+ .5*g*m[:,0]
+    return gas, w, z, ps, un, curl_vk, abs_vort_vk
 
 class myDavies(Davies):
@@ -153,10 +172,8 @@
     comm = client.comm
     mpi_rank, mpi_size = comm.Get_rank(), comm.Get_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
+    g, Lx, Ly = 9.80616, 4e7, 6e6
     nx, ny, llm = args.nx, args.ny, args.llm
     dx = Lx/nx
@@ -170,5 +187,4 @@
     #    print('WARNING: Davies_N1 not equal to Davies_N2; average value will be used')
 
-
     unst.setvar('g',g)
     
@@ -198,21 +214,20 @@
         scheme = time_step.ARK2(None, dt)
         if args.hydrostatic:
-            caldyn_step = unst.caldyn_step_HPE(mesh,scheme,1, caldyn_thermo,caldyn_eta,thermo,params,params.g)
+            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)
+            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() )
         logging.debug('mask min/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
-            #Relaxation zone inside the domain
-            kappa = 3.26e+15 #nonhyd=>#Y3.3#Y3.27#Y3.26#Y3.255#dblchkd;N3.25#N3.2#N3.0#Y3.4#dx**4/432000.
-            #Relaxation zone outside the domain
-            #kappa = 3.6e+15 #nonhyd=>#Y3.3#Y3.27#Y3.26#Y3.255#dblchkd;N3.25#N3.2#N3.0#Y3.4#dx**4/432000.
-            #print('Kappa=',kappa)
             for i in range(nt):
                 caldyn_step.next()
@@ -226,6 +241,6 @@
                 unst.ker.dynamico_curl_laplacian(u,lapu)
                 unst.ker.dynamico_curl_laplacian(lapu,bilapu)
-                caldyn_step.theta_rhodz[:] = S - dt*kappa*bilaps*m # Euler step
-                caldyn_step.u[:] = u - dt*kappa*bilapu # Euler step
+                caldyn_step.theta_rhodz[:] = S - dt*args.kappa_divgrad*bilaps*m # Euler step
+                caldyn_step.u[:] = u - dt*args.kappa_curlcurl*bilapu # Euler step
 
             return (caldyn_step.mass.copy(), caldyn_step.theta_rhodz.copy(), caldyn_step.u.copy(),
@@ -241,21 +256,16 @@
         title_format = 'Specific entropy at t=%g s (J/K/kg)'
 
-    w=mesh.field_mass()
-    z=mesh.field_mass()
-
     #T, nslice, dt = 3600., 1, 3600.
     Nslice = args.N
 
+    temp_v = mesh.field_z(),     
+
     with xios.Context_Curvilinear(mesh,1, 24*3600) as context:
         # now XIOS knows about the mesh and we can write to disk
-        v = mesh.field_mass() # specific volume (diagnosed)
-        
         for i in range(Nslice+1):
             context.update_calendar(i+1)
 
             # 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)
+            gas, w, z, ps, un, zeta_vk, zeta_abs_vk = diagnose(Phi,S,m,W,u)
             KE_ik = KE(mesh,u)
             du_fast, du_slow = caldyn_step.du_fast[0,:,:], caldyn_step.du_slow[0,:,:] 
@@ -277,7 +287,7 @@
 
             context.send_field_primal('div_slow', div_slow)
-#            context.send_field_primal('div_slow', div(mesh,hflux) )
-
-            context.send_field_dual('curl', curl_vk)
+
+            context.send_field_dual('curl', zeta_vk)
+            context.send_field_dual('curl_abs', zeta_abs_vk)
 
             print( 'ptop, model top (m) :', unst.getvar('ptop'), Phi.max()/unst.getvar('g'))
