Index: codes/icosagcm/devel/Python/test/py/NH_3D_DCMIP31.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/NH_3D_DCMIP31.py	(revision 632)
+++ codes/icosagcm/devel/Python/test/py/NH_3D_DCMIP31.py	(revision 632)
@@ -0,0 +1,147 @@
+from dynamico import unstructured as unst
+from dynamico import dyn
+from dynamico import time_step
+from dynamico import DCMIP
+from dynamico import meshes
+from dynamico.meshes import MPAS_Mesh as Mesh
+import math as math
+import matplotlib.pyplot as plt
+import numpy as np
+import time
+
+#------------------------ initial condition -------------------------
+
+Cpd, Rd, g = 1004.5, 287., 9.81
+N, Teq, peq = 0.01, 300., 1e5   # background 
+lonc, d2 = 2.*math.pi/3., 25e6  # perturbation
+N2, g2 = N*N, g*g
+N2_g2, g2_N2 = N2/g2, g2/N2
+G, kappa = g2/(N2*Cpd), Rd/Cpd
+
+def psTs(lat) : 
+    cos2latm1 = np.cos(2*lat)-1.
+    Ts = G+(Teq-G)*np.exp(-.25*N2_g2*u0*u0*cos2latm1)
+    ps = peq*np.exp(u0/(4.*G*Rd)*u0*cos2latm1)*((Ts/Teq)**(1./kappa))
+    return ps, Ts
+    
+def Phi_top(lat) : 
+    ps, Ts = psTs(lat)
+    return -g2_N2*np.log( 1+(Ts/G)*( ((ptop/ps)**kappa) - 1 ) )
+    
+def T0_p0(lon,lat,Phi):
+    ps, Ts = psTs(lat)
+    pb = ps*( (G/Ts)*(np.exp(-N2_g2*Phi)-1.)+1. )**(1./kappa)
+    Tb = G + (Ts-G)*np.exp(N2_g2*Phi)
+    r  = radius*np.arccos(np.cos(lat)*np.cos(lon-lonc))
+    Theta = dT*np.sin(math.pi*Phi/(g*ztop))*d2/(d2+r*r)
+    T = Tb + Theta*(pb/peq)**kappa
+    return T, pb
+
+def DCMIP31_3D(grid):
+    def Phi(eta, lat): return eta*Phi_top(lat)
+    def ulon(lat): return u0*np.cos(lat)
+    def ulat(lat): return 0.*lat
+    def f(lon,lat): return 0.*np.sin(lat) # Coriolis parameter
+
+
+    nqdyn, preff, Treff = 1, 1e5, 300.
+    thermo = dyn.Ideal_perfect(Cpd, Rd, preff, Treff)
+    mesh = Mesh('grids/x1.%d.grid.nc'%grid, llm, nqdyn, radius, f)
+    lev,levp1 = np.arange(llm),np.arange(llm+1)
+    lat_ik,k_i = np.meshgrid(mesh.lat_i,lev, indexing='ij')
+    lon_ik,k_i = np.meshgrid(mesh.lon_i,lev, indexing='ij')
+    lat_il,l_i = np.meshgrid(mesh.lat_i,levp1, indexing='ij')            
+    lon_il,l_i = np.meshgrid(mesh.lon_i,levp1, indexing='ij')
+    lat_ek,k_e = np.meshgrid(mesh.lat_e,lev, indexing='ij')
+
+    Phi_ik, Phi_il = Phi(k_i/float(llm),lat_ik), Phi(l_i/float(llm),lat_il)
+    Tb_ik, pb_ik = T0_p0(lon_ik, lat_ik, Phi_ik)
+    Tb_il, pb_il = T0_p0(lon_il, lat_il, Phi_il)
+    gas = thermo.set_pT(pb_ik,Tb_ik)
+    mass_ik = mesh.field_mass()
+    for l in range(llm):
+        mass_ik[:,l]=(pb_il[:,l]-pb_il[:,l+1])/g
+    Sik  = gas.s*mass_ik    
+
+# start at hydrostatic geopotential
+    pb_ik = .5*(pb_il[:,1:]+pb_il[:,0:-1]) # pressure at full layers
+    gas = thermo.set_ps(pb_ik,gas.s)
+    for l in range(llm):
+        Phi_il[:,l+1] = Phi_il[:,l] + g*mass_ik[:,l]*gas.v[:,l]
+
+    ujk, Wil = mesh.ucov3D(ulon(lat_ek),ulat(lat_ek)), mesh.field_w()
+    
+    print 'ztop (m) = ', Phi_il[0,-1]/g, ztop
+    print 'ptop (Pa) = ', gas.p[0,-1], ptop
+    dx=mesh.de.min()
+    params=dyn.Struct()
+    params.g, params.ztop, params.ptop = g, ztop, ptop
+    params.dx, params.dx_g0 = dx, dx/g
+    params.pbot, params.rho_bot = 1e5+0.*mesh.lat_i, 1e6+0.*mesh.lat_i
+    return thermo, mesh, params, (mass_ik,Sik,ujk,Phi_il,Wil), gas
+
+#------------------------ main program -------------------------
+
+ztop, u0, dT, Xfactor =  1e4, 20., 1., 125.
+radius, grid, llm = 6.37122e6/Xfactor, 10242, 20
+T, Nslice, courant = 360., 10, 3.0
+
+junk, ptop = T0_p0(0.,0.,g*ztop)
+
+thermo, mesh, params, flow0, gas0 = DCMIP31_3D(grid)
+llm, dx = mesh.llm, params.dx
+dz = params.ztop/llm
+print 'llm, nb_hex, dx, dz =', llm, mesh.Ai.size, dx, dz
+
+caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_lag
+caldyn = unst.Caldyn_NH(caldyn_thermo,caldyn_eta, mesh,params,thermo,params,params.g)
+
+mesh.plot_e(mesh.le/mesh.de) ; plt.title('le/de')
+plt.savefig('fig_NH_3D_DCMIP31/le_de.png') ;plt.close()
+
+mesh.plot_i(mesh.Ai) ; plt.title('Ai')
+plt.savefig('fig_NH_3D_DCMIP31/Ai.png'); plt.close()
+
+#dt = courant*.5/np.sqrt(gas0.c2.max()*(dx**-2+dz**-2))
+dt = courant*.5*dx/np.sqrt(gas0.c2.max())
+
+nt = int(math.ceil(T/dt))
+dt = T/nt
+print 'Time step : %g s' % dt
+
+#scheme = time_step.RK4(caldyn.bwd_fast_slow, dt)
+scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, a32=0.7)
+
+flow=flow0
+w=mesh.field_mass()
+z=mesh.field_mass()
+
+def plots(it):
+    junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
+    m,S,u,Phi,W = flow
+    s=S/m ; s=.5*(s+abs(s))
+    for l in range(llm):
+        w[:,l]=.5*params.g*(W[:,l+1]+W[:,l])/m[:,l]
+        z[:,l]=.5*(Phi[:,l+1]+Phi[:,l])/params.g
+        
+    print( 'ptop, model top (m), max(w), max(W) (m/s) :', 
+          unst.getvar('ptop'), Phi.max()/unst.getvar('g'), w.max(), W.max() )
+
+    mesh.plot_i(w[:,llm/2])
+    plt.title('vertical velocity at t=%d'%(it*T))
+    plt.savefig('fig_NH_3D_DCMIP31/w%02d.png'%it)
+    plt.close()
+
+plots(0)
+
+for it in range(Nslice):
+    unst.setvar('debug_hevi_solver',True)
+    time1, elapsed1 =time.time(), unst.getvar('elapsed')
+    flow = scheme.advance(flow,nt)
+    time2, elapsed2 =time.time(), unst.getvar('elapsed')
+    factor = 1000./(3*nt)
+    print 'ms per call to caldyn_hevi : ', factor*(time2-time1), factor*(elapsed2-elapsed1)
+    factor = 1e9/(4*nt*w.size)
+    print 'nanosec per gridpoint per call to caldyn_hevi : ', factor*(time2-time1), factor*(elapsed2-elapsed1)
+    plots(it+1)
+    
Index: codes/icosagcm/devel/Python/test/py/NH_3D_bubble.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/NH_3D_bubble.py	(revision 632)
+++ codes/icosagcm/devel/Python/test/py/NH_3D_bubble.py	(revision 632)
@@ -0,0 +1,121 @@
+from dynamico import unstructured as unst
+from dynamico import dyn
+from dynamico import time_step
+from dynamico import DCMIP
+from dynamico import meshes
+from dynamico.meshes import Cartesian_mesh as Mesh
+import math as math
+import matplotlib.pyplot as plt
+import numpy as np
+import time
+
+def thermal_bubble_3D(Lx,nx,Ly,ny,llm,ztop=1000., zc=350.,
+                      rc=250, thetac=0.5, x0=0., y0=0.):
+    Cpd, Rd, g, p0,theta0, T0 = 1004.5, 287.,9.81, 1e5, 300., 300.
+    nqdyn = 1
+    
+    Phi = lambda eta : g*ztop*eta
+    p=lambda Phi : p0*np.exp(-Phi/(Rd*T0))
+    zz = lambda p: -(Rd*T0*np.log(p/p0))/g
+    rr = lambda x,y,p: np.sqrt((x-x0)**2 + (y-y0)**2 + (zz(p)-zc)**2)
+    sa = lambda x,y,p: rr(x,y,p) < rc
+    deform = lambda x,y,p: (0.5*thetac*(1+np.cos(np.pi*rr(x,y,p)/rc)))*sa(x,y,p)
+    temp =  lambda p: theta0*(p/p0)**(Rd/Cpd)
+    T = lambda x,y,p: deform(x,y,p) + temp(p)
+
+    mesh = Mesh(nx,ny,llm,nqdyn,Lx,Ly,0.)
+    thermo = dyn.Ideal_perfect(Cpd, Rd, p0, T0)
+
+    Phi_il = Phi(mesh.llp1/float(llm))
+    Phi_ik = Phi((mesh.ll+.5)/llm)
+    p_ik = p(Phi_ik)
+    T_ik = T(mesh.xx, mesh.yy, p_ik)
+
+    gas = thermo.set_pT(p_ik,T_ik)
+    mass_ik = mesh.field_mass()
+    for l in range(llm): 
+        mass_ik[:,:,l]=(Phi_il[:,:,l+1]-Phi_il[:,:,l])/(g*gas.v[:,:,l])
+    Sik, ujk, Wil  = gas.s*mass_ik, mesh.field_u(), mesh.field_w()
+    
+    print 'ztop (m) = ', Phi_il[0,0,-1]/g, ztop
+    ptop = p(g*ztop)
+    print 'ptop (Pa) = ', gas.p[0,0,-1], ptop
+    params=dyn.Struct()
+    params.ptop=ptop
+    params.dx=dx
+    params.dx_g0=dx/g
+    params.g = g
+    pbot = p(Phi_il[:,:,0])
+    gas_bot = thermo.set_pT(pbot, temp(pbot))
+    params.pbot = gas_bot.p
+    params.rho_bot = 1e6/gas_bot.v
+    
+    return thermo, mesh, params, (mass_ik,Sik,ujk,Phi_il,Wil), gas
+
+Lx, nx, llm, thetac, T, Nslice, courant = 2000., 100, 50, 30., 5., 10, 2.8
+#Lx, nx, llm, thetac, T, Nslice, courant = 2000., 50, 25, 1., 50., 10, 2.8
+#Lx, nx, llm, thetac, T, Nslice, courant = 3000., 75, 25, -30, 5., 10, 2.8
+
+nqdyn, dx = 1, Lx/nx
+Ly,ny,dy = Lx,nx,dx
+thermo, mesh, params, flow0, gas0 = thermal_bubble_3D(Lx,nx,Ly,ny,llm, thetac=thetac)
+
+caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass
+caldyn = unst.Caldyn_NH(caldyn_thermo,caldyn_eta, mesh,params,thermo,params,params.g)
+
+# compute hybrid coefs from initial distribution of mass
+mass_bl,mass_dak,mass_dbk = meshes.compute_hybrid_coefs(flow0[0])
+unst.ker.dynamico_init_hybrid(mass_bl,mass_dak,mass_dbk)
+print mass_dbk
+
+dz = flow0[3].max()/(params.g*llm)
+#dt = courant*.5/np.sqrt(gas0.c2.max()*(dx**-2+dy**-2+dz**-2))
+dt = courant*.5/np.sqrt(gas0.c2.max()*(dx**-2+dy**-2))
+nt = int(math.ceil(T/dt))
+dt = T/nt
+print 'Time step : %g s' % dt
+
+scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt)
+
+flow=flow0
+w=mesh.field_mass()
+z=mesh.field_mass()
+for it in range(Nslice):
+    junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
+    m,S,u,Phi,W = flow
+    s=S/m ; # s=.5*(s+abs(s))
+    for l in range(llm):
+        w[:,:,l]=.5*params.g*(W[:,:,l+1]+W[:,:,l])/m[:,:,l]
+        z[:,:,l]=.5*(Phi[:,:,l+1]+Phi[:,:,l])/params.g
+        
+    print 'ptop, model top (m) :', unst.getvar('ptop'), Phi.max()/unst.getvar('g')
+    jj=ny/2
+    xx,zz,ss,ww = mesh.xx[:,jj,:]/1000., z[:,jj,:]/1000., s[:,jj,:], w[:,jj,:]
+    
+    f, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(12,4))
+
+    c=ax1.contourf(xx,zz,ss,20)
+    ax1.set_xlim((-.5,.5)), ax1.set_xlabel('x (km)')
+    ax1.set_ylim((0.,1.)), ax1.set_ylabel('z (km)')        
+    plt.colorbar(c,ax=ax1)
+    ax1.set_title('Specific entropy at t=%g s (J/K/kg)' % (it*T,))
+#    plt.show()
+    
+#    plt.figure(figsize=(12,5))
+    c=ax2.contourf(xx,zz,ww,20)
+    ax2.set_xlim((-.5,.5)), ax2.set_xlabel('x (km)')
+    ax2.set_ylim((0.,1.)), ax2.set_ylabel('z (km)')        
+    plt.colorbar(c,ax=ax2)
+    ax2.set_title('Vertical velocity at t=%g s (m/s)' % (it*T,))
+#    plt.tight_layout()
+#    plt.show()
+    plt.savefig('fig_NH_3D_bubble/%02d.png'%it)
+    
+    time1, elapsed1 =time.time(), unst.getvar('elapsed')
+    flow = scheme.advance(flow,nt)
+    time2, elapsed2 =time.time(), unst.getvar('elapsed')
+    factor = 1000./(4*nt)
+    print 'ms per call to caldyn_hevi : ', factor*(time2-time1), factor*(elapsed2-elapsed1)
+    factor = 1e9/(4*nt*nx*ny*llm)
+    print 'nanosec per gridpoint per call to caldyn_hevi : ', factor*(time2-time1), factor*(elapsed2-elapsed1)
+        
Index: codes/icosagcm/devel/Python/test/py/RSW_2D.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/RSW_2D.py	(revision 631)
+++ codes/icosagcm/devel/Python/test/py/RSW_2D.py	(revision 632)
@@ -8,6 +8,4 @@
 import numpy as np
 import time
-
-unst.setvar('nb_threads', 1)
 
 nx,ny,llm,nqdyn=128,128,1,1
Index: codes/icosagcm/devel/Python/test/py/RSW_MPAS_W02.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/RSW_MPAS_W02.py	(revision 631)
+++ codes/icosagcm/devel/Python/test/py/RSW_MPAS_W02.py	(revision 632)
@@ -79,2 +79,4 @@
     print i, h.min(), h.max()
     flow = scheme.advance(flow, nt)
+
+print 'Time spent in DYNAMICO (s) : ', unst.getvar('elapsed')
