Index: /codes/icosagcm/devel/Python/test/clean.sh
===================================================================
--- /codes/icosagcm/devel/Python/test/clean.sh	(revision 617)
+++ /codes/icosagcm/devel/Python/test/clean.sh	(revision 617)
@@ -0,0 +1,2 @@
+#!/bin/bash
+rm -f fig_*/*
Index: /codes/icosagcm/devel/Python/test/get_MPAS_grids.sh
===================================================================
--- /codes/icosagcm/devel/Python/test/get_MPAS_grids.sh	(revision 617)
+++ /codes/icosagcm/devel/Python/test/get_MPAS_grids.sh	(revision 617)
@@ -0,0 +1,29 @@
+#!/bin/bash
+
+function get_grid()
+{
+  GRID=$1
+  if [ -f $GRID ];
+  then
+      echo "Grid file $GRID is already there."
+  else
+      FILE=$2
+      rm -f $FILE
+#      wget "http://www2.mmm.ucar.edu/projects/mpas/atmosphere_meshes/$FILE"
+      scp ${user}@forge.ipsl.fr:/home/dubos/grids/$FILE . && tar xzf $FILE $GRID
+  fi
+}
+
+echo "Usage : get_MPAS_grids.sh username SIZE1 SIZE2 ..."
+echo "where username is your username on forge.ipsl.fr ;"
+echo "optional arguments SIZEx are among : x1.163842 x1.655362 x1.2621442 x1.5898242 x4.163842"
+echo "By default, download small meshes for sizes 2562 10242 40962"
+
+user=$1 ; shift
+SIZES="x1.2562 x1.10242 x1.40962 x4.163842 $*"
+
+cd grids
+
+for SIZE in $SIZES ; do 
+    get_grid "$SIZE.grid.nc" "$SIZE.tar.gz" 
+done
Index: /codes/icosagcm/devel/Python/test/py/RSW_2D.py
===================================================================
--- /codes/icosagcm/devel/Python/test/py/RSW_2D.py	(revision 617)
+++ /codes/icosagcm/devel/Python/test/py/RSW_2D.py	(revision 617)
@@ -0,0 +1,63 @@
+from dynamico import unstructured as unst
+from dynamico import dyn
+from dynamico import time_step
+from dynamico import DCMIP
+import math as math
+import matplotlib.pyplot as plt
+import numpy as np
+import time
+
+unst.setvar('nb_threads', 1)
+
+nx,ny,llm,nqdyn=128,128,1,1
+Lx,Ly,g,f = 8.,8.,1.,1.
+dx,dy=Lx/nx,Ly/ny
+
+unst.setvar('g',g)
+mesh = unst.Cartesian_mesh(nx,ny,llm,nqdyn,Lx,Ly,f)
+caldyn = unst.Caldyn_RSW(mesh)
+
+x1,x2,yy = mesh.xx[:,:,0]-1., mesh.xx[:,:,0]+1., mesh.yy[:,:,0]
+h0 = 1+0.1*(np.exp(-2.*(x1*x1+yy*yy))+
+            np.exp(-2.*(x2*x2+yy*yy)))
+u0 = mesh.field_u()  # flow initally at rest
+flow0=(h0,u0)
+
+cfl = .8
+dt = cfl/math.sqrt((nx/Lx)**2+(ny/Ly)**2)
+
+T=10.
+N=int(T/dt)+1
+dt=T/N
+print N,dt,Lx/nx
+scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt)
+
+flow=flow0
+for i in range(10):
+    h,u=flow
+    caldyn.bwd_fast_slow(flow, 0.)
+    plt.figure(); plt.pcolor(mesh.x,mesh.y,caldyn.qv)
+    plt.colorbar(); plt.title('potential vorticity')
+    plt.savefig('fig_RSW_2D/%02d.png'%i)
+#    plt.figure(); plt.pcolor(mesh.x,mesh.y,h)
+#    plt.colorbar(); plt.title('h')
+#    plt.show()
+#    plt.pcolor(x,y,vcomp(u)/dx)
+#    plt.colorbar(); plt.title('v')
+#    plt.show()
+    for j in range(5):
+        unst.elapsed=0.
+        flow = scheme.advance(flow,N)
+        # flops
+        # mass flux : 1add+1mul per edge          =>  4
+        # div U     : 4 add per cell              =>  4
+        # KE        : 4*(2add+1mul) per cell      => 12
+        # grad KE   : 1 add per edge              =>  2
+        # grad h    : 1 add per edge              =>  2
+        # qv    : 4+4+1 add +4mul + 1div per cell => 10
+        # qu    : 1add+1mul per edge              =>  4
+        # TrisK : 4add+4mul+4add+1add per edge    => 26
+        # Total :                                    64 FLOPS/cell
+        print i,j, unst.elapsed, 100.*(1.-unst.getvar('elapsed')/unst.elapsed)
+        print 'GFlops', 64*(N*nx*ny)/unst.elapsed/1e9
+        unst.setvar('elapsed',0.)
Index: /codes/icosagcm/devel/Python/test/py/RSW_MPAS_W02.py
===================================================================
--- /codes/icosagcm/devel/Python/test/py/RSW_MPAS_W02.py	(revision 617)
+++ /codes/icosagcm/devel/Python/test/py/RSW_MPAS_W02.py	(revision 617)
@@ -0,0 +1,89 @@
+import sys
+print 'Loading modules ...'
+sys.stdout.flush()
+
+import math as math
+# select non-interactive backend, cf http://stackoverflow.com/questions/4931376/generating-matplotlib-graphs-without-a-running-x-server
+import matplotlib
+matplotlib.use('Agg') 
+import matplotlib.pyplot as plt
+import numpy as np
+
+print 'Loading DYNAMICO modules ...'
+sys.stdout.flush()
+from dynamico import unstructured as unst
+from dynamico import time_step
+print '...Done'
+sys.stdout.flush()
+
+grid, llm, nqdyn = 10242, 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
+
+print 'Omega, planetary PV', Omega, 2*Omega/gh0
+sys.stdout.flush()
+
+def f(lon,lat): return 2*Omega*np.sin(lat) # Coriolis parameter
+print 'Reading MPAS mesh ...'
+sys.stdout.flush()
+mesh = unst.MPAS_Mesh('grids/x1.%d.grid.nc'%grid, llm, nqdyn, radius, f)
+print '...Done'
+sys.stdout.flush()
+lon, lat = mesh.lon_i, mesh.lat_i
+x,y,z = np.cos(lat)*np.cos(lon), np.cos(lat)*np.sin(lon), np.sin(lat)
+
+unst.setvar('g',g)
+caldyn = unst.Caldyn_RSW(mesh)
+
+c0 = math.sqrt(gh0) # phase speed of barotropic mode
+dx = mesh.de.min()
+dt = courant*dx/c0
+nt = int(math.ceil(T/dt))
+dt = T/nt
+#scheme = time_step.RKn_simple(1,caldyn.bwd_fast_slow, dt)
+scheme = time_step.RK4(caldyn.bwd_fast_slow, dt)
+    
+print dx, dt, dt*c0/dx
+print T, nt
+
+#mesh.plot_e(mesh.le_de) ; plt.title('le/de') ; plt.show()
+#mesh.plot_i(mesh.Ai) ; plt.title('Ai') ; plt.show()
+
+#Phi0 = gh0*(1+0.1*np.exp(-2000.*(y+1.)**2))
+#flow = (Phi0,mesh.field_u())
+#for i in range(N):
+#    h,u=flow
+#    plot_i(h) ; plt.title('h'); plt.show()
+#    junk, fast, slow = caldyn.bwd_fast_slow(flow,0.)
+#    print i
+#    flow = scheme.advance(flow, nt)
+
+u0 = Omega*radius/12. # cf Williamson (1991), p.13
+gh1 = radius*Omega*u0+.5*u0**2
+print 'Williamson (1991) test 2, u0=', u0
+ulon = u0*np.cos(mesh.lat_e)
+Phi0 = gh0 - gh1*(np.sin(mesh.lat_i)**2)
+zeta0 = (2*u0/radius+2*Omega)*np.sin(mesh.lat_v)
+Phi0v = gh0 - (radius*Omega*u0+.5*u0**2)*(np.sin(mesh.lat_v)**2)
+q0 = zeta0/Phi0v
+fu_perp = mesh.ucov2D(0.*ulon,2*Omega*np.sin(mesh.lat_e)*ulon)
+
+flow = (Phi0,mesh.ucov2D(ulon,0.*ulon))
+for i in range(N):
+    h,u=flow
+    mesh.plot_i(h-Phi0) ; plt.title('err(gh)'); 
+    plt.savefig('fig_RSW_MPAS_W02/err_gh_%02d.png'%i)
+    plt.close()
+#    junk, fast, slow = caldyn.bwd_fast_slow(flow,0.)
+#    plot_i(slow[0]) ; plt.title('dh/dt'); plt.show()
+#    plot_e(slow[1]+fast[1]) ; plt.title('du/dt'); plt.show()
+#    plt.figure(); plt.plot(fu_perp,slow[1],'.');
+#    plt.xlabel('fu_perp'); plt.ylabel('u_slow'); plt.show()
+#    plt.figure(); plt.plot(fu_perp,fast[1],'.');
+#    plt.xlabel('fu_perp'); plt.ylabel('u_fast'); plt.show()
+#    plot_e(fast[1]) ; plt.title('u_fast'); plt.show()
+#    plot_e(slow[1]) ; plt.title('u_slow'); plt.show()
+#    plt.figure(); plt.plot(q0,caldyn.qv,'.');
+#    plt.xlabel('q0'); plt.ylabel('qv'); plt.show()
+    print i, h.min(), h.max()
+    flow = scheme.advance(flow, nt)
Index: /codes/icosagcm/devel/Python/test/py/bubble.py
===================================================================
--- /codes/icosagcm/devel/Python/test/py/bubble.py	(revision 617)
+++ /codes/icosagcm/devel/Python/test/py/bubble.py	(revision 617)
@@ -0,0 +1,193 @@
+import math as math
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.animation as manimation
+
+import dynamico.dyn as dyn
+import dynamico.time_step as time_step
+import dynamico.DCMIP as DCMIP
+import dynamico.advect as advect
+from dynamico import unstructured as unst
+
+def reldiff(a,b):
+    return (a-b)/(np.absolute(a).max()+np.absolute(b).max()+1e-20)
+
+def go(model,courant,scheme,start):
+    def plotter(metric, flow, scalars, t):
+        mik, gas_ik, ujk, Phi_il, Wil = model.diagnose(*flow)
+        Q1, Q2 = scalars
+        
+        mil = metric.ml_ext(mik)
+        wil = Wil/mil*metric.g0
+        wik = metric.mk_ext(Wil)/mik*metric.g0
+        xik = metric.xik/1000
+        xjl = metric.xjl/1000
+        xil = metric.xil/1000
+        zil = Phi_il/metric.g0/1000
+        zik = metric.mk_int(zil)
+        zjl = metric.mj(zil)
+        print 's ',np.max(gas_ik.s),np.min(gas_ik.s)
+        print 'ujk ',np.max(ujk),np.min(ujk)
+        print 'wil ',np.max(wil),np.min(wil)
+        sik = gas_ik.s
+        sik = .5*(abs(sik)+sik)
+
+        plt.clf()
+        plt.contourf(xik, zik, sik, 20 )
+#        plt.contourf(xik, zik, wik, 20 )
+        plt.xlim((-1,1)), plt.xlabel('x (km)')
+        plt.ylim((0,ztop/1000.)), plt.ylabel('z (km)')
+        plt.colorbar()
+        plt.title("t = %f s" % (t))
+        plt.savefig('fig_bubble/bubble.png',bbox_inches='tight')
+        writer.grab_frame()
+        
+    transport = advect.SplitTransport(model.metric, advect.HorizVanLeer, advect.VertVanLeer)
+    replace = lambda fd,fv : fv
+    return DCMIP.run_with_scalars(model,scheme,transport,replace, plotter, courant, start,scalars, T1, N1, N2)
+
+Lx, ztop, nx, llm = 2000.,1000., 100, 50
+nqdyn, ny, dx = 1, 1, Lx/nx
+T1, N1, N2 = 5., 1, 200
+#T1, N1, N2 = 5., 2, 5
+
+metric, thermo, BC, m0ik, gas0, Phi0_il, u0_jk = DCMIP.thermal_bubble(Lx,nx,llm, ztop=ztop, zc=350.)
+Sik = m0ik*gas0.s
+start = (m0ik, Sik, u0_jk, Phi0_il, 0*Phi0_il)
+scalars = (Sik, m0ik*metric.mk_int(Phi0_il))
+
+model = dyn.Nonhydro(thermo,metric,BC)
+solver, caldyn_eta = dyn.MassBasedNH(model), unst.eta_mass
+#solver, caldyn_eta = dyn.LagrangianNH(model), unst.eta_lag
+
+rho_bot=BC.rho_bot
+if False: # explicit time-stepping with rigid lid
+    def rigid(flow,tau): # impose rigid BCs at top and bottom
+        flow, fast, slow = solver.bwd_fast_slow(flow,0.)
+        fast[4][:,0]=0.
+        fast[4][:,-1]=0.
+        return flow,fast,slow
+    BC.rho_bot = 100*rho_bot # moderate stiffness
+    courant, time_scheme = 2.0, lambda dt : time_step.RK4(rigid, dt)
+else: # HEVI time-stepping with pressure BC
+    # Note : parameters BC.rigid_top and BC.rigid_bottom are True by default, but are ignored by the Python code...
+    BC.rigid_top=False
+    BC.rigid_bottom=False
+    BC.rho_bot = 1e6*rho_bot # large stiffness
+#    BC.rho_bot = 100*rho_bot # moderate stiffness
+    courant, time_scheme = 3.0, lambda dt : time_step.ARK2(solver.bwd_fast_slow, dt, a32=0.7)
+
+FFMpegWriter = manimation.writers['ffmpeg']
+metadata = dict(title='Aquaplanet', artist='DYNAMICO_LMDZ5',
+                comment='Movie support!')
+writer = FFMpegWriter(fps=15, metadata=metadata, bitrate=4000, codec="h263p")
+fig = plt.figure(figsize=(12,5))
+with writer.saving(fig, "fig_bubble/bubble.avi", 100):
+    m0ik, gas0_ik, u0jk, Phi0_il, W0il = go(model, courant, time_scheme, start)
+
+#-------------------------------------------------------------------------------------------------------
+# Now use final state to check that Fortran and Python implementations produce the exact same tendencies
+#-------------------------------------------------------------------------------------------------------
+
+#W0il[:,0]=0.
+#W0il[:,-1]=0.
+
+unst.setvar('nb_threads', 1)
+mesh = unst.Cartesian_mesh(nx,ny,llm,nqdyn,Lx,ny*dx,0.)
+xx_ik, xx_il, ll = mesh.xx[:,0,:]/1000, mesh.xxp1[:,0,:]/1000, mesh.ll[:,0,:]
+
+mass_bl,mass_dak,mass_dbk=unst.compute_hybrid_coefs(m0ik)
+unst.init_hybrid(mass_bl,mass_dak,mass_dbk)
+
+caldyn_thermo = unst.thermo_entropy
+#caldyn_thermo = unst.thermo_theta
+g = mesh.dx/metric.dx_g0
+caldyn = unst.Caldyn_NH(caldyn_thermo,caldyn_eta, mesh,metric,thermo,BC,g)
+unst.setvar('rigid',False)
+
+s0ik = gas0_ik.s
+if caldyn_thermo == unst.thermo_theta:
+    s0ik = thermo.T0*np.exp(s0ik/thermo.Cpd)
+
+# 2 momentum components in DYNAMICO vs 1 component in slice
+u0jk = u0jk + 100*dx # shift horizontal wind by 100 m/s
+u0_ker = mesh.field_u()
+mesh.set_ucomp(u0_ker,u0jk)
+
+tau=0.1
+# impose hybdrid distribution of mass (Fortran)
+flow = (m0ik/dx,s0ik*m0ik/dx,u0_ker,Phi0_il,W0il/dx) # NB : different definitions of mass field in DYNAMICO / Python
+flow,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
+# shift Phi from target value Phi_bot
+Phi0_il = Phi0_il+g*100.*np.exp(-10.*xx_il*xx_il)
+m,S,u,Phi,W = flow
+flow = m,S,u,Phi0_il,W
+# compute tendencies (0=Python)
+flow0 = (m*dx,S*dx,u0jk,Phi0_il,W0il, 0., 0.) # NB : different definitions of mass field in DYNAMICO / Python
+flow0, fast0, slow0 = solver.bwd_fast_slow(flow0,tau,True) 
+# compute tendencies (Fortran)
+flow,fast,slow = caldyn.bwd_fast_slow(flow, tau) 
+
+
+zz_il = Phi0_il/metric.g0/1000.
+zz_ik=.5*(zz_il[:,0:-1]+zz_il[:,1:])
+print
+print 'ptop, model top (km) :', unst.getvar('ptop'), zz_il.max()
+        
+for data0,data,dname in ((flow0,flow,'flow'),(slow0,slow,'slow'),(fast0,fast,'fast')):
+    def plot(xx,zz,var,vname):
+        title = "%s_%s" % (vname,dname)
+        filename = "bubble/bubble_%s.png" % title
+        var = 1.*var
+        if var.size>1:
+            plt.figure(figsize=(12,5))
+            plt.contourf(xx,zz,var,20)
+#                    plt.pcolor(xx,zz,var)
+            plt.xlim((-1.,1.)), plt.xlabel('x (km)')
+            plt.ylim((0.,1.)), plt.ylabel('z (km)')        
+            plt.colorbar()
+            plt.title(title)
+            plt.savefig(filename,bbox_inches='tight')
+            print title, var.min(), var.max()
+            if var.max()>1e100:
+                raise OverflowError
+                    
+    dm,dS,du,dPhi,dW = data    
+    du=mesh.ucomp(du)/dx
+    du=du[0,:,:]
+
+    dm0, dS0, du0, dPhi0, dW0, junk, junk = data0
+    dm0, dS0, du0, dW0 = dm0/dx, dS0/dx, du0/dx, dW0/dx
+
+    if dname == 'flow':
+        hflux=mesh.ucomp(caldyn.hflux)
+        plot(xx_ik, zz_ik, dm,'m')
+        plot(xx_ik, zz_ik, du,'u')
+        plot(xx_ik, zz_ik, du-du0,'uu')
+        plot(xx_ik, zz_ik, dS/m,'s')
+        plot(xx_il, zz_il, dW,'w')
+        plot(xx_il, zz_il, dW-dW0,'ww')
+    else:
+        dS=reldiff(dS,dS0)
+        dW=dW-dW0
+        dPhi=dPhi-dPhi0
+        du=reldiff(du,du0)
+        plot(xx_ik, zz_ik, dS,'dS')
+        plot(xx_il, zz_il, dPhi,'dPhi')
+        plot(xx_il, zz_il, dW,'dW')
+        plot(xx_ik, zz_ik, du,'du')
+
+if unst.getvar('rigid'):
+    print 'Fortran BCs : rigid'
+else:
+    print 'Fortran BCs : pressure/soft', unst.getvars('ptop','rho_bot','Phi_bot','pbot')
+
+if BC.rigid_top:
+    print 'Top BC : rigid'
+else:
+    print 'Top BC : pressure', BC.ptop
+
+if BC.rigid_bottom:
+    print 'Bottom BC : rigid'
+else:
+    print 'Bottom BC : soft', BC.rho_bot[0], BC.Phi_bot[0], BC.pbot[0]
Index: /codes/icosagcm/devel/Python/test/py/slice_GW_NH.py
===================================================================
--- /codes/icosagcm/devel/Python/test/py/slice_GW_NH.py	(revision 617)
+++ /codes/icosagcm/devel/Python/test/py/slice_GW_NH.py	(revision 617)
@@ -0,0 +1,105 @@
+from dynamico import unstructured as unst
+from dynamico import dyn
+from dynamico import time_step
+from dynamico import DCMIP
+
+import math as math
+import matplotlib.pyplot as plt
+import numpy as np
+import time
+
+unst.setvar('nb_threads', 1)
+
+def run(mesh,metric,thermo,BC,caldyn_eta,u, m0ik,gas0,Phi0_il, T=300., courant=2.9):
+    caldyn_thermo = unst.thermo_entropy
+    g = mesh.dx/metric.dx_g0
+    caldyn = unst.Caldyn_NH(caldyn_thermo,caldyn_eta, mesh,metric,thermo,BC,g)
+
+    Sik = m0ik*gas0.s
+    m0ik = m0ik/dx # NB : different definitions of mass field in DYNAMICO / Slice
+    if caldyn_thermo == unst.thermo_theta:
+        theta0 = thermo.T0*np.exp(gas0.s/thermo.Cpd)
+        S0ik = m0ik*theta0
+    else:
+        S0ik = m0ik*gas0.s
+
+    u0=mesh.field_u() 
+    u0[:,range(0,2*nx,2),:] = u*mesh.dx # Doppler shift by u
+    flow0=(m0ik,S0ik,u0,Phi0_il,0*Phi0_il)
+
+    dt = courant*.5*dx/np.sqrt(gas0.c2.max())
+    nt=int(T/dt)+1
+    dt=T/nt
+    print 'nt,dt,dx', nt,dt,dx
+    scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, a32=0.7)
+    #scheme = time_step.ARK3(caldyn.bwd_fast_slow, dt)
+
+#    nt=1 # FIXME
+    flow=flow0
+    for i in range(12):
+        time1=time.time()
+        flow = scheme.advance(flow,nt)
+        time2=time.time()
+        print 'ms per time step : ', 1000*(time2-time1)/nt
+
+        m,S,u,Phi,W = flow
+        w=.5*(W[:,0:llm]+W[:,1:llm+1])/m
+        junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
+#        zz=caldyn.geopot[:,0:llm]/metric.g0/1000
+        zz=Phi[:,0:llm]/metric.g0/1000
+        print 'ptop, model top (m) :', unst.getvar('ptop'), Phi.max()/unst.getvar('g')
+
+        def plt_axes():
+            plt.plot(xx,zz[:,llm/4])
+            plt.plot(xx,zz[:,llm/2])
+            plt.plot(xx,zz[:,3*llm/4])
+            plt.ylim(0,ztop/1000.)
+            plt.colorbar()
+
+        plt.figure(figsize=(12,3))
+        ucomp = mesh.ucomp(u)
+        plt.pcolor(xx,zz,caldyn.pk)
+        plt_axes()
+        plt.title("T(t = %f s)" % ((i+1)*dt*nt))
+        plt.savefig('fig_slice_GW_NH/T%02d.png'%i)
+        plt.close()
+
+        plt.figure(figsize=(12,3))
+        ucomp = mesh.ucomp(u)
+        plt.pcolor(xx,zz,ucomp[0,:,:]/dx)
+        plt_axes()
+        plt.title("u(t = %f s)" % ((i+1)*dt*nt))
+        plt.savefig('fig_slice_GW_NH/U%02d.png'%i)
+        plt.close()
+
+        plt.figure(figsize=(12,3))
+        plt.pcolor(xx,zz,w)
+        plt_axes()
+        plt.title("w(t = %f s)" % ((i+1)*dt*nt))
+        plt.savefig('fig_slice_GW_NH/W%02d.png'%i)
+        plt.close()
+
+def deform(x,eta): # x and eta go from 0 to 1
+    h0, dd, xi = 2000., 5000., 4000.
+    cos2=lambda x: np.cos(np.pi*x)**2
+    gauss=lambda x: np.exp(-x*x)
+    xx = x*Lx
+    phis = h0*gauss(xx/dd)*cos2(xx/xi)
+    return (phis/ztop)*np.sin(np.pi*eta)
+
+#Lx, nx, llm = 3e5, 100, 1e4, 30
+Lx, nx, ztop, llm = 2e5, 400, 3e4, 60
+nqdyn, ny, dx = 1, 1, Lx/nx
+mesh = unst.Cartesian_mesh(nx,ny,llm,nqdyn,nx*dx,ny*dx,0.)
+xx,ll = mesh.xx[:,0,:]/1000, mesh.ll[:,0,:]
+
+#metric, thermo, BC, m0ik, gas0, Phi0_il, u0_jk = DCMIP.DCMIP31(Lx, nx, llm, deform=deform, dTheta=1.)
+metric, thermo, BC, m0ik, gas0, Phi0_il, u0_jk = DCMIP.DCMIP21(Lx, nx, llm, deform=deform, h0=0.)
+
+BC.rho_bot = 1e6*BC.rho_bot # large stiffness
+
+mass_bl,mass_dak,mass_dbk = unst.compute_hybrid_coefs(m0ik)
+unst.init_hybrid(mass_bl,mass_dak,mass_dbk)
+
+run(mesh,metric,thermo,BC,unst.eta_mass,20.,  m0ik,gas0,Phi0_il, courant=2.) # lower courant number required if h0=2000m
+#run(mesh,metric,thermo,BC,unst.eta_lag,20.,  m0ik,gas0,Phi0_il)
Index: /codes/icosagcm/devel/Python/test/py/slice_GW_hydro.py
===================================================================
--- /codes/icosagcm/devel/Python/test/py/slice_GW_hydro.py	(revision 617)
+++ /codes/icosagcm/devel/Python/test/py/slice_GW_hydro.py	(revision 617)
@@ -0,0 +1,64 @@
+from dynamico import unstructured as unst
+from dynamico import dyn
+from dynamico import time_step
+from dynamico import DCMIP
+
+import math as math
+import matplotlib.pyplot as plt
+import numpy as np
+import time
+
+unst.setvar('nb_threads', 1)
+
+def run(mesh,metric,thermo,BC,caldyn_eta, m0ik,gas0):
+    caldyn_thermo = unst.thermo_entropy
+    g = mesh.dx/metric.dx_g0
+    caldyn = unst.Caldyn_HPE(caldyn_thermo,caldyn_eta, mesh,metric,thermo,BC,g)
+
+    Sik = m0ik*gas0.s
+    m0ik = m0ik/dx # NB : different definitions of mass field in DYNAMICO / Slice
+    if caldyn_thermo == unst.thermo_theta:
+        theta0 = thermo.T0*np.exp(gas0.s/thermo.Cpd)
+        S0ik = m0ik*theta0
+    else:
+        S0ik = m0ik*gas0.s
+
+    u0=mesh.field_u() 
+    u0[:,range(0,2*nx,2),:] = 100.*mesh.dx # Doppler shift 100 m/s
+    flow0=(m0ik,S0ik,u0)
+
+    #T, courant = 60., 1.
+    T, courant = 60., 2.9
+    dt = courant*.5*dx/np.sqrt(gas0.c2.max())
+    nt=int(T/dt)+1
+    dt=T/nt
+    print 'nt,dt,dx', nt,dt,dx
+    scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, a32=0.7)
+    #scheme = time_step.ARK3(caldyn.bwd_fast_slow, dt)
+
+    flow=flow0
+    for i in range(10):
+        time1=time.time()
+        flow = scheme.advance(flow,nt)
+        time2=time.time()
+        m,S,u = flow
+        print 'ms per time step : ', 1000*(time2-time1)/nt
+        print 'ptop, model top (m) :', unst.getvar('ptop'), caldyn.geopot.max()/unst.getvar('g')
+        junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
+        zz=caldyn.geopot[:,0:llm]/metric.g0/1000
+        plt.figure(figsize=(12,3))
+        ucomp = mesh.ucomp(u)
+        plt.pcolor(xx,zz,ucomp[0,:,:]/dx)
+        plt.ylim(0,10.)
+        plt.colorbar()
+        plt.title("u(t = %f s)" % ((i+1)*dt*nt))
+        plt.savefig('fig_slice_GW_hydro/%02d.png'%i)
+        plt.close()
+
+Lx, nx, llm = 3e5, 300, 20
+nqdyn, ny, dx = 1, 1, Lx/nx
+mesh = unst.Cartesian_mesh(nx,ny,llm,nqdyn,nx*dx,ny*dx,0.)
+xx,ll = mesh.xx[:,0,:]/1000, mesh.ll[:,0,:]
+metric, thermo, BC, m0ik, gas0, Phi0_il, u0_jk = DCMIP.DCMIP31(Lx, nx, llm)
+
+run(mesh,metric,thermo,BC,unst.eta_lag,  m0ik,gas0)
Index: /codes/icosagcm/devel/Python/test/py/slice_orographic.py
===================================================================
--- /codes/icosagcm/devel/Python/test/py/slice_orographic.py	(revision 617)
+++ /codes/icosagcm/devel/Python/test/py/slice_orographic.py	(revision 617)
@@ -0,0 +1,81 @@
+import math as math
+import matplotlib.pyplot as plt
+import numpy as np
+import dynamico.dyn as dyn
+import dynamico.time_step as time_step
+import dynamico.DCMIP as DCMIP
+import dynamico.advect as advect
+
+dx, Nx, llm, tau_relax = 500., 400, 60, 25.
+#dx, Nx, llm, tau_relax = 250., 800, 120, 25.
+T1, N1, N2 = 30., 6, 20 # check energy every 60s, plot every 10min, run 10h
+#dx, Nx, llm, tau_relax = 20000., 200, 200, 1e3
+#T1, N1, N2 = 1800., 4, 6 # check energy every 60s, plot every 10min, run 2h
+
+Lx=Nx*dx
+metric, thermo, BC, m0ik, gas0, Phi0_il, u0_jk = DCMIP.DCMIP21(
+    Lx, Nx, llm, h0=250., ztop=3e4, dd=4000, xi=3200)
+#    Lx, Nx, llm, h0=250., ztop=3e4, dd=10*dx, xi=8*dx)
+Sik = m0ik*gas0.s
+start =(m0ik, Sik, u0_jk, Phi0_il, 0*Phi0_il)
+scalars = (Sik, m0ik*metric.mk_int(Phi0_il))
+
+relax = lambda Phi : (1./tau_relax)*(1+np.tanh((Phi/9.81-23e3)/3e3))
+
+def sponge(Phi_il, vjk, dujk_dt, tau):
+    tau_rel = 1./relax(Phi_il[:,0:-1])
+    dujk = (vjk-u0_jk + tau*dujk_dt)/(tau+tau_rel)
+    ujk = u0_jk + tau_rel*dujk
+    dujk_dt = dujk_dt - dujk
+    return ujk, dujk_dt
+
+BC.sponge=sponge # will be called by fast_bwd
+BC.rigid_top = False
+BC.rigid_bottom = False
+rho_bot = BC.rho_bot
+BC.rho_bot = 1e6*rho_bot
+
+def go(model,courant,scheme,start):
+    def plotter(metric, flow, scalars, t):
+        mik, gas_ik, ujk, Phi_il, Wil = model.diagnose(*flow)
+        mil = metric.ml_ext(mik)
+        wil = metric.g0*Wil/mil
+        xil = metric.xil/1000
+        xjl = metric.xjl/1000
+        zil = Phi_il/metric.g0/1000
+        zjl = metric.mj(zil)
+        plt.figure(figsize=(12,5))
+#        plt.figure(figsize=(6,5))
+#        plt.pcolormesh(xjl, zjl, gas_ik.T-gas0.T)
+#        plt.contourf(xil, zil, wil, np.arange(-.2,.2,.02))
+        plt.contourf(xil, zil, wil, np.arange(-1.,1.,.1))
+        plt.plot(xil,zil[:,llm/10])
+        plt.plot(xil,zil[:,llm/5])
+        plt.plot(xil,zil[:,llm/4])
+        plt.plot(xil,zil[:,llm/2])
+        plt.xlabel('x (km)') 
+        plt.ylabel('z (km)')
+#        plt.ylim((-1,30.)), plt.xlim((-dx/20.,dx/10.))
+        plt.ylim((0.,10.)), plt.xlim((-10.,10.))
+        plt.colorbar()
+        plt.title("t = %f s" % (t))
+        plt.savefig("fig_slice_orographic/%04d"%t)
+        plt.close()
+    transport = advect.SplitTransport(model.metric, advect.HorizVanLeer, advect.VertVanLeer)
+    replace = lambda fd,fv : fv
+    DCMIP.run_with_scalars(model,scheme,transport,replace, plotter, courant, start,scalars, T1, N1, N2)
+
+model = dyn.Nonhydro(thermo,metric,BC)
+solver = dyn.MassBasedNH(model)
+#solver = dyn.LagrangianNH(model)
+
+def bwd_fast_slow(flow,tau):
+    flow,fast,slow=solver.bwd_fast_slow(flow,tau)
+    m,S,u,Phi,W,X,Y = flow
+    dm,dS,du,dPhi,dW,dX,dY = slow
+    dW = dW - relax(Phi)*W # damp vertical velocity
+    slow = dm,dS,du,dPhi,dW,dX,dY
+    return flow,fast,slow
+
+time_scheme = lambda dt : time_step.ARK2(bwd_fast_slow, dt, a32=0.7)
+go(model, 3.0, time_scheme, start )
Index: /codes/icosagcm/devel/Python/test/python.sh
===================================================================
--- /codes/icosagcm/devel/Python/test/python.sh	(revision 617)
+++ /codes/icosagcm/devel/Python/test/python.sh	(revision 617)
@@ -0,0 +1,7 @@
+#!/bin/bash
+# pleae do not move this script as relies on ../.. being DYNAMICO_ROOT
+DYNAMICO_ROOT=$(cd -P $(dirname $0)/../.. ; pwd)
+. $DYNAMICO_ROOT/arch.env
+. $DYNAMICO_ROOT/Python//env/$CURRENT_ARCH.env
+
+python -u $*
