Index: /codes/icosagcm/devel/Python/codegen
===================================================================
--- /codes/icosagcm/devel/Python/codegen	(revision 615)
+++ /codes/icosagcm/devel/Python/codegen	(revision 615)
@@ -0,0 +1,65 @@
+#!/bin/bash
+
+function filter()
+{
+    PASS=$1 ; shift
+    cpp -P -DPASS_$PASS $* -include ../macros.jin
+}
+
+function log()
+{
+    echo $*
+    $*
+}
+
+function main()
+{
+    VARIANT=$SRCDIR/$1
+    cd $VARIANT
+    rm -f *~ tmp/*
+    mkdir -p tmp
+    cd tmp
+    cp $SRCDIR/macros.jin . # jinja macros common to all variants
+    cpp -P -DPASS_JINJA ../macros.jin >> macros.jin # variant-specific macros
+
+    # first expand user-defined CPP and jinja macros
+    cat $SRCDIR/kernels_*.jin | cpp -P > caldyn_all.jin
+    $PYTHON -u $ROOT/jinja caldyn_all.jin > caldyn_all.dysl
+    # then expand the code-generating CPP and jinja macros
+    cat caldyn_all.dysl | filter PRE1 | filter PRE2 >> macros.jin
+    $PYTHON -u $ROOT/jinja macros.jin > caldyn_all.k90
+    # extract each kernel and apply post-jinja CPP macros
+    KERNELS=$(grep KERNEL caldyn_all.dysl | cpp -P -D'KERNEL(x)=x')
+    for NAME in $KERNELS ; do
+	cat caldyn_all.k90 | filter POST1 -DKERNEL_${NAME} | filter POST2 | sed '/^$/d' > ../${NAME}.k90
+	head ../${NAME}.k90
+    done
+
+    cd ..
+    grep ERROR *.k90
+
+    cd $VARIANT
+    wc *.k90
+    cd $ROOT
+}
+
+if [ -z "$CODEGEN_ENV_SET" ] ; then
+    cat <<EOF
+Before running this script, you must source one of the *.env scripts found in the env/ directory.
+This will set up appropriate, machine-dependent values for the following environment variables :
+PYTHON F90 C99 CFLAGS FFLAGS LINK SOFLAGS LDFLAGS LIBS ...
+EOF
+else
+    set +x
+    
+    ROOT=$PWD
+    LIBDIR=$ROOT/lib
+    SRCDIR=$ROOT/src
+    cd $SRCDIR
+
+    for x in $* ; do
+	echo "Generating code for variant : $x"
+	grep FIXME *.jin
+	main $x
+    done
+fi
Index: /codes/icosagcm/devel/Python/dynamico/DCMIP.py
===================================================================
--- /codes/icosagcm/devel/Python/dynamico/DCMIP.py	(revision 615)
+++ /codes/icosagcm/devel/Python/dynamico/DCMIP.py	(revision 615)
@@ -0,0 +1,193 @@
+import matplotlib.pyplot as plt
+import time as time
+import math as math
+import numpy as np
+from numpy import cos, sin, exp, log
+from math import pi
+import dyn as dyn
+
+Linf = lambda x : abs(x).max()
+
+def run_with_scalars(model,dyn_scheme,transport,replace, plotter, 
+                     courant,flow0,scalars, T1, N1, N2):
+    metric, BC = model.metric, model.BC
+    mik, gas0, ujk, Phi_il, Wil = model.diagnose(*flow0)
+    ztop = Phi_il.max()/metric.g0
+    Phi_ik = metric.mk_int(Phi_il)
+    vol, jac = metric.volume(Phi_il[:,-1])
+    internal0 = (mik*(gas0.e+Phi_ik)).sum() + BC.ptop*vol.sum()
+
+    dtmax, dt_hydro = model.max_time_step(gas0,mik)
+    dt=courant*dtmax # RK4 stability limit is courant=2.sqrt(2)
+    
+    N0 = 2*int(math.ceil(.5*float(T1)/dt))
+    dt=T1/N0 # adjust time step to fraction of T1
+    step = dyn_scheme(dt)
+    print 'dt (s)', step.dt, 'Horizontal acoustic Courant number', dt/dt_hydro
+    
+    Sik = mik*gas0.s
+    flow0 = mik, Sik, ujk, Phi_il, Wil, 0.,0.
+    elapsed=time.time()
+    flow=step.advance(flow0,10)
+    elapsed = (time.time()-elapsed)/10
+    print "Elapsed time per time step (ms)", elapsed*1e3
+    print "Simulated time / elapsed time  (s/s)", step.dt/elapsed
+        
+    plotter(metric, (mik, Sik, ujk, Phi_il, Wil), scalars, 0.)
+
+    out=[flow0]
+    kinetic=[]
+    internal=[]
+    for i in range(N2):
+        for j in range(N1):
+            for k in range(N0/2):
+                mik_old = mik
+                mik, Sik, ujk, Phi_il, Wil, Fjk, Fil = step.next((mik, Sik, ujk, Phi_il, Wil, 0.,0.))
+                mik, Sik, ujk, Phi_il, Wil, Fjk, Fil = step.next((mik, Sik, ujk, Phi_il, Wil, Fjk, Fil))
+                if type(Fil) is type(Sik):# fix Lagrangian coordinate
+                    # transport scalars using time-integrated mass fluxes Fjk, Fil
+                    mik_new, scalars = transport.next(Fjk,Fil,mik_old,scalars)
+                    Sik = replace(Sik,scalars[0]) # replace (or not) entropy by FV-transported scalar
+                    # print Linf(mik-mik_new)/Linf(mik)
+            mik, gas_ik, ujk, Phi_il, Wil = model.diagnose(mik, Sik, ujk, Phi_il, Wil)
+            Phi_ik = metric.mk_int(Phi_il)
+            vol, jac = metric.volume(Phi_il[:,-1])
+            kinetic.append(model.kinetic(mik, ujk, Wil, Phi_il))
+            internal.append((mik*(gas_ik.e+Phi_ik)).sum() + BC.ptop*vol.sum())
+            print '\r', i*N1+j+1, '/', N1*N2,
+        plotter(metric, (mik, Sik, ujk, Phi_il, Wil), scalars, (i+1)*T1*N1)
+        
+    kinetic=np.array(kinetic)
+    internal=np.array(internal)
+    dE=kinetic+internal-internal0
+    plt.figure(figsize=(12,3))
+    plt.plot(T1*np.arange(N1*N2),dE/kinetic.max())
+    plt.xlabel('time')
+    plt.ylabel('(E(t)-E(0))/Kmax')
+    
+    plt.figure(figsize=(12,3))
+    plt.plot(T1*np.arange(N1*N2),kinetic/kinetic.max())
+    plt.xlabel('time')
+    plt.ylabel('K(t)/Kmax')
+    plt.show()
+
+    return model.diagnose(mik, Sik, ujk, Phi_il, Wil)
+
+def DCMIP21(Lx,Nx,llm, ztop=3e4, h0=250., 
+            dd=5000., xi=4000., Ueq=20.,
+            deform=lambda x,eta:0 ) :
+    """ DCMIP2.1 experiment at equator (phi=0)"""
+    Cpd, Rd, pref, Tref, = 1004.5, 287., 1e5, 273.15
+    g, p0, peq, Teq = 9.81, 1e5, 1e5, 300.
+
+    p=lambda x,Phi : peq*np.exp(-Phi/(Rd*Teq))
+    ptop = p(0,g*ztop)
+    print 'ptop (Pa) =', ptop
+
+    xmesh  = dyn.Periodic_FD(Nx)
+    zmesh  = dyn.Lorenz(llm)
+    mesh   = dyn.Slice(xmesh,zmesh)
+    metric = dyn.Shallow(mesh, Lx/Nx,g)
+    thermo = dyn.Ideal_perfect(Cpd, Rd, pref, Tref)
+
+    T=lambda x,Phi : 0*x+Teq
+    gauss=lambda x: np.exp(-x*x)
+    cos2=lambda x: cos(pi*x)**2
+    Phis = lambda x:g*h0*gauss(x/dd)*cos2(x/xi)
+    Phi = lambda x,eta : g*ztop*eta + (1.-eta)*Phis(x)
+    Phi_deform = lambda x,eta : Phi(x,eta+deform(x/Lx,eta))
+    mik, gas, Phi_il = metric.inicond_NH(thermo, Phi_deform,p,T)
+    ujk = 0*mik + Ueq*(Lx/Nx)
+    
+    print 'ztop (m) = ', Phi_il[0,-1]/g, ztop
+    print 'ptop (Pa) = ', gas.p[0,-1], ptop
+    Phi_bot = Phi_il[:,0:1]
+    gas_bot=thermo.set_pT(p(metric.xil[:,0:1], Phi_bot), Teq)
+    BC = dyn.BoundCond(ptop, Phi_bot, gas_bot.p, 1./gas_bot.v)
+    return metric, thermo, BC, mik, gas, Phi_il, ujk
+
+def DCMIP31(Lx,Nx,llm,deform=lambda x,eta:0, dTheta=1., d2=25e6) :
+    """ DCMIP3.1 initial condition with Phi(x,eta)=g.ztop*(eta + deform(eta,x/L)) """
+    # DCMIP 3.1 at equator (phi=0)
+    Cpd, Rd, pref, Tref, = 1004.5, 287., 1e5, 273.15
+    g, p0, ptop = 9.81, 1e5, 27391.9
+    peq, N2, Teq, ztop = 1e5, 1e-4, 300., 1e4
+
+    G=g**2/(Cpd*N2)
+    inv_G=1./G
+    inv_Teq=1./Teq
+    G_Teq = G/Teq
+
+    p=lambda x,Phi : peq*((1-G_Teq*(1-np.exp(-N2*Phi/g**2)))**(Cpd/Rd))
+    ptop = p(0,g*ztop)
+
+    xmesh  = dyn.Periodic_FD(Nx)
+    zmesh  = dyn.Lorenz(llm)
+    mesh   = dyn.Slice(xmesh,zmesh)
+    metric = dyn.Shallow(mesh, Lx/Nx,g)
+    thermo = dyn.Ideal_perfect(Cpd, Rd, pref, Tref)
+
+    ps = lambda xi: peq+0*xi
+    fac = lambda p : (p/peq)**(Rd/Cpd)
+    T = lambda x,p : fac(p)/(inv_G*(fac(p)-1)+inv_Teq)
+#    T = lambda x,p : 0.*p+Teq
+
+    Phi = lambda x,eta : g*ztop*(eta+deform(x/Lx,eta))
+    m0ik, gas, Phi0_il = metric.inicond_NH(thermo, Phi,p,T)
+    u0_jk = 0*m0ik
+#    u0_jk = 0*m0ik + 100.*Lx/Nx
+    
+    Tb  = gas.T # unperturbed temperature
+
+    print 'ztop (m) = ', Phi0_il[0,-1]/g, ztop
+    print 'ptop (Pa) = ', gas.p[0,-1], ptop
+    Phi_bot = Phi0_il[:,0:1]
+    pbot = p(metric.xil[:,0:1], Phi_bot)
+    gas_bot=thermo.set_pT(pbot, T(metric.xil[:,0:1],pbot))
+    BC = dyn.BoundCond(ptop, Phi_bot, gas_bot.p, 1./gas_bot.v)
+    print BC.rho_bot.shape
+    
+    # now add localized temperature perturbation
+    xik = metric.xik
+    pik = gas.p
+    zik = zmesh.mk_int(Phi0_il)/g
+    shape = d2/(d2+xik*xik)
+    shape = shape*np.sin(math.pi*zik/ztop)
+    shape = shape*((pik/p0)**(Rd/Cpd))
+    gas0 = thermo.set_pT(pik, Tb + dTheta*shape)
+
+    return metric, thermo, BC, m0ik, gas0, Phi0_il, u0_jk
+
+def thermal_bubble(Lx,Nx,llm, ztop=1000., zc=200., rc=250, thetac=0.5, x0=0.0):
+    """Hi"""
+    Cpd, Rd, g, p0,theta0, T0 = 1004.5, 287.,9.81, 1e5, 300., 300.
+
+    Phis = lambda x: 0.0
+    Phi = lambda x,eta : g*ztop*eta + (1.-eta)*Phis(x)
+    p=lambda x, Phi : p0*np.exp(-Phi/(Rd*T0))
+    ptop = p(0,g*ztop)
+    print 'ptop (Pa) =', ptop
+    xmesh  = dyn.Periodic_FD(Nx)
+    zmesh  = dyn.Lorenz(llm)
+    mesh   = dyn.Slice(xmesh,zmesh)
+    metric = dyn.Shallow(mesh, Lx/Nx,g)
+    thermo = dyn.Ideal_perfect(Cpd, Rd, p0, T0)
+    
+    zz = lambda p: -(Rd*T0*np.log(p/p0))/g
+    rr =   lambda x, p: np.sqrt((x-x0)**2 + (zz(p)-zc)**2)
+    sa = lambda x,p: rr(x,p) < rc
+    deform = lambda x,p: (0.5*thetac*(1+cos(pi*rr(x,p)/rc)))*sa(x,p)
+    temp =  lambda x, p: theta0*(p/p0)**(Rd/Cpd)
+    T = lambda x,p: deform(x,p) + temp(x,p)
+       
+    #T = lambda x, Phi: T0
+    mik, gas, Phi_il = metric.inicond_NH(thermo, Phi,p,T)
+    ujk = 0*mik
+        
+    print 'ztop (m) = ', Phi_il[0,-1]/g, ztop
+    print 'ptop (Pa) = ', gas.p[0,-1], ptop
+    Phi_bot = Phi_il[:,0:1]
+    gas_bot=thermo.set_pT(p(metric.xil[:,0:1], Phi_bot), T0) #### T ot T0 
+    BC = dyn.BoundCond(ptop, Phi_bot, gas_bot.p, 1./gas_bot.v)
+    return metric, thermo, BC, mik, gas, Phi_il, ujk
+
Index: /codes/icosagcm/devel/Python/dynamico/advect.py
===================================================================
--- /codes/icosagcm/devel/Python/dynamico/advect.py	(revision 615)
+++ /codes/icosagcm/devel/Python/dynamico/advect.py	(revision 615)
@@ -0,0 +1,168 @@
+import numpy as np
+from scipy import sparse as sparse
+
+#-------------------------------- Strang-split advection --------------------------------
+
+class SplitTransport:
+    def __init__(self, metric, Horiz, Vert):
+        # Horiz and Vert are classes for horizontal and vertical transport, respectively
+        self.horiz = Horiz(metric)
+        self.vert = Vert(metric)
+        self.di, self.dk = metric.di, metric.dk
+    def next(self, Fjk, Fil,  # time-integrated mass fluxes
+             mik,             # mass field before advection
+             scalars):        # tuple of 2D mass-weighted scalars to be advected
+        scalars =  [Qik for Qik in self.vert.advect(.5*Fil,mik,scalars)]
+        mik = mik - .5*self.dk(Fil)
+        scalars =  [Qik for Qik in self.horiz.advect(Fjk,mik,scalars)]
+        mik = mik - self.di(Fjk)
+        scalars =  [Qik for Qik in self.vert.advect(.5*Fil,mik,scalars)]
+        mik = mik - .5*self.dk(Fil)
+        return mik, scalars
+    
+#---------------------------- Horizontal advection schemes ----------------------------
+
+class HorizGodunov:
+    def __init__(self,metric):
+        self.di, self.dj = metric.di, metric.dj
+        self.mi, self.mj = metric.mi, metric.mj
+    # member function reconstruct() (to be overloaded by derived classes)
+    # yields for each point j the mean and difference of right and left reconstructions 
+    def reconstruct(self,Fjk,maik,scalars):
+        for Qik in scalars:
+            qik = Qik/maik
+            yield Qik, self.mj(qik), self.dj(qik) # piecewise-constant reconstruction
+    def advect(self,Fjk,maik,scalars):
+        for Qik, mean_qr_ql, qr_minus_ql in self.reconstruct(Fjk,maik,scalars):
+            qjk = mean_qr_ql - 0.5*np.sign(Fjk)*qr_minus_ql  # Upwind sub-grid reconstruction 
+            yield Qik - self.di(qjk*Fjk)
+            
+class HorizVanLeer(HorizGodunov):
+    def reconstruct(self,Fjk,mik,scalars):
+        fjk = .5*Fjk/self.mj(mik) # displacement, same for all scalars
+        for Qik in scalars:
+            qik = Qik/mik
+            delqik = self.di(self.mj(qik))
+            mean_qr_ql =   self.mj(qik) - .25*self.dj(delqik) - fjk*self.mj(delqik)
+            qr_minus_ql =  self.dj(qik) - self.mj(delqik)     - fjk*self.dj(delqik)
+            yield Qik, mean_qr_ql, qr_minus_ql
+            
+#----------------------------- Vertical advection schemes ---------------------------------
+
+class VertGodunov:
+    def __init__(self,metric):
+        self.dk, self.dl, self.ml = metric.dk, metric.dl, metric.ml_int
+    # member function reconstruct() (to be overloaded by derived classes)
+    # yields for each point j the mean and difference of right and left reconstructions 
+    def advect(self,Fil,mik,scalars):
+        for Qik, mean_qr_ql, qr_minus_ql in self.reconstruct(Fil,mik,scalars):
+            qil = mean_qr_ql - 0.5*np.sign(Fil)*qr_minus_ql  # Upwind sub-grid reconstruction 
+            yield Qik - self.dk(qil*Fil)
+    def reconstruct(self,Fil,mik,scalars):
+        for Qik in scalars:
+            qik = Qik/mik
+            yield Qik, self.ml(qik), self.dl(0.,qik,0.)
+
+class VertVanLeer(VertGodunov):
+    def reconstruct(self,Fil,mik,scalars):
+        fil = .5*np.array(Fil/self.ml(mik)) # displacement, same for all scalars
+        for Qik in scalars:
+            qik = Qik/mik
+            qil = self.ml(qik)
+            delqik = self.dk(qil) # subgrid slope, FIXME at top and bottom
+            mean_qr_ql  = self.ml(qik)     - .25*self.dl(0,delqik,0) - fil*self.ml(delqik)
+            qr_minus_ql = self.dl(0,qik,0) - self.ml(delqik)         - fil*self.dl(0,delqik,0)
+            yield Qik, mean_qr_ql, qr_minus_ql
+    
+#--------------------- Sarvesh original code -------------------------
+
+class adv:
+    def __init__(self,metric,maik,qik,Uajk,Fail):
+        self.metric = metric
+        self.maik=maik
+        self.qik=qik 
+        self.Uajk=Uajk
+        self.Fail=Fail 
+        self.di, self.dj = metric.di, metric.dj
+        self.dk, self.dl = metric.dk, metric.dl
+        self.mi, self.mj = metric.mi, metric.mj
+        Nx, Ny = np.shape(maik)
+        ones=np.ones(Nx);
+        left = sparse.spdiags([ones],[-1],Nx,Nx, format='csr')
+        left[0,-1]=1
+        right = sparse.spdiags([ones],[1],Nx,Nx, format='csr')
+        right[-1,0]=1
+        self.left =lambda f:left.dot(f)
+        self.right =lambda f:right.dot(f)
+        self.mk_ext, self.ml_ext = metric.mk_ext, metric.ml_ext
+        self.mk_int, self.ml_int = metric.mk_int, metric.ml_int
+
+    def adv_horizon(self,maik,qik,Uajk):
+        maik_old = maik  # Mass at the start of the advection (before N-adv dynamical steps)  
+        qik_r = self.right(qik)
+        maik_r = self.right(maik)
+        qik_l = self.left(qik)
+        delqik = 0.5*(qik_r-qik_l)
+        qr = qik + 0.5*(1 - Uajk/maik)*delqik
+        ql = qik_r - 0.5*(1 + Uajk/maik_r)*self.right(delqik)
+        qUPjk= 0.5*((qr+ql) + np.sign(Uajk)*(qr-ql))  ##Upwind sub-grid reconstruction 
+        qxflx = qUPjk*Uajk
+        maik = maik - self.di(Uajk) ## Updated Mass 
+        qik = (qik*maik_old - self.di(qxflx))/maik
+        return maik,qik
+    
+    def adv_horizon1(self,maik,qik,Uajk):
+        maik_old = maik  # Mass at the start of the advection (before N-adv dynamical steps)  
+        qjk = self.mj(qik)
+        delqik = self.di(qjk)
+        qik_r = qik + 0.5*(1 - Uajk/maik)*delqik
+        qik_l = np.zeros(np.shape(qik))
+        qik_l[:-1,:] = qik[1:,:]-0.5*(1.0 - Uajk[:-1,:]/maik[1:,:])*delqik[1:,:]
+        qik_l[-1,:] = qik[0,:] - 0.5*(1.0 - Uajk[-1,:]/maik[0,:])*delqik[0,:]
+        qUPjk= 0.5*((qik_r+qik_l) + np.sign(Uajk[:,:])*(qik_r-qik_l))
+        qxflx = qUPjk*Uajk
+        maik = maik - self.di(Uajk) ## Updated Mass 
+        qik = (qik*maik_old - self.di(qxflx))/maik
+        return maik,qik
+    
+    def adv_godnov(self,maik,qik,Uajk):
+        maik_old = maik  # Mass at the start of the advection (before N-adv dynamical steps)  
+        qjk = self.mj(qik)
+        qik_r = qik 
+        qik_l = np.zeros(np.shape(qik))
+        qik_l[:-1,:] = qik[1:,:]
+        qik_l[-1,:] = qik[0,:] 
+        qUPjk= 0.5*((qik_r+qik_l) + np.sign(Uajk[:,:])*(qik_r-qik_l))
+        qxflx = qUPjk*Uajk
+        maik = maik - self.di(Uajk) ## Updated Mass 
+        qik = (qik*maik_old - self.di(qxflx))/maik
+        return maik,qik
+    
+    def adv_vert(self,maik,qik,Fil):
+        maik_old = maik
+        qil = self.ml_int(qik)
+        delqik = self.dk(qil) # subgrid slope 
+        qrik = qik[:,:-1] + 0.5*(1 - Fil[:,1:-1]/maik[:,:-1])*delqik[:,:-1]
+        qlik = qik[:,1:] - 0.5*(1 - Fil[:,1:-1]/maik[:,1:])*delqik[:,1:]
+      
+        qUPil = np.zeros(np.shape(Fil))
+        qUPil[:,1:-1]= 0.5*((qrik+qlik) + np.sign(Fil[:,1:-1])*(qrik-qlik))
+        qUPil[:,0]= 0.0 #qik[:,0] - 0.5*(1 - Fil[:,0]/maik[:,0])*delqik[:,0]
+        qUPil[:,-1]= 0.0 #qik[:,-1] + 0.5*(1 - Fil[:,-1]/maik[:,-1])*delqik[:,-1]
+        qzflx = qUPil*Fil
+        maik = maik - self.dk(Fil)
+        qik = (qik*maik_old - self.dk(qzflx))/maik
+        return maik,qik
+    
+    def adv_vgod(self,maik,qik,Fil):
+        maik_old = maik
+        qrik = qik[:,:-1]
+        qlik = qik[:,1:]       
+        qUPil = np.zeros(np.shape(Fil))
+        qUPil[:,1:-1]= 0.5*((qrik+qlik) + np.sign(Fil[:,1:-1])*(qrik-qlik))
+        qUPil[:,0]= 0.0
+        qUPil[:,-1]= 0.0
+        qzflx = qUPil*Fil
+        maik = maik - self.dk(Fil)
+        qik = (qik - self.dk(qzflx)) #/maik
+        return maik,qik
Index: /codes/icosagcm/devel/Python/dynamico/dyn.py
===================================================================
--- /codes/icosagcm/devel/Python/dynamico/dyn.py	(revision 615)
+++ /codes/icosagcm/devel/Python/dynamico/dyn.py	(revision 615)
@@ -0,0 +1,497 @@
+import numpy as np
+from scipy import sparse as sparse
+import math as math
+import time_step as time_step
+
+class Struct: pass
+
+def Linf(data):
+    return abs(data).max()
+def compare(msg, data1, data2):
+    print msg, Linf(data1-data2)/Linf(data1)    
+
+def col_vec(data):
+    vec=np.zeros((data.size,1))
+    vec[:,0]=data
+    return vec
+def row_vec(data):
+    vec=np.zeros((1,data.size))
+    vec[0,:]=data
+    return vec
+
+zero=np.zeros((1))
+
+####################### Horizontal mesh #########################
+
+class Periodic_FD:
+    # staggering
+    # m0  m1  m2  m3
+    #   u0  u1  u2  u3
+    def __init__(self, Nx):
+        self.Nx=Nx
+        ones=np.ones(Nx);
+
+        right = sparse.spdiags([ones],[0],Nx,Nx, format='csr')
+        left = sparse.spdiags([ones],[-1],Nx,Nx, format='csr')
+        left[0,-1]=1
+        di = right-left
+        mi = .5*(right+left)
+        self.di=lambda f:di.dot(f)
+        self.mi=lambda f:mi.dot(f)
+        self.xi = col_vec(np.linspace(0,Nx-1,Nx))-Nx/2.
+        self.onesi = np.ones((Nx,1))
+        self.zeroi = np.zeros((Nx,1))
+        
+        right = sparse.spdiags([ones],[1],Nx,Nx, format='csr')
+        right[-1,0]=1
+        left = sparse.spdiags([ones],[0],Nx,Nx, format='csr')
+        dj = right-left
+        mj = .5*(right+left)
+        self.dj=lambda f:dj.dot(f)
+        self.mj=lambda f:mj.dot(f)
+        self.xj = col_vec(np.linspace(0.5,Nx-.5,Nx))-Nx/2.
+        self.onesj = np.ones((Nx,1))
+
+############################## Vertical mesh ############################
+
+class Lorenz:
+    def __init__(self, llm):
+        self.llm = llm
+        self.onesk = np.ones((1,llm))
+        self.onesl = np.ones((1,llm+1))
+        self.etak = row_vec(np.linspace(.5/llm,1.-.5/llm,llm))
+        self.etal = row_vec(np.linspace(0.,1.,llm+1))
+    def fieldk(self,Nx):
+        return np.zeros((Nx,self.llm))
+    def fieldl(self,Nx):
+        return np.zeros((Nx,self.llm+1))
+    def dk(self, fl):
+        return fl[:,1:]-fl[:,:-1]
+    def dl(self, fbot, fk, ftop): 
+    # vertical FD from k to l with extra info at top/bottom 
+        llm=self.llm
+        gbot = fk[:,0:1]-fbot 
+        gl = fk[:,1:]-fk[:,:-1]
+        gtop = ftop-fk[:,llm-1:llm]
+        return np.array(np.hstack((gbot,gl,gtop)))
+    # Vertical averaging, fac=1 or 1/2 is weight of boundary contributions
+    # mk and ml with same fac are adjoint
+    # ml_ext, mk_int use fac=1/2
+    # mk_ext, ml_int use fac=1
+    def ml(self, fk, fac): # fk has size llm, last index llm-1
+        llm=self.llm
+        fbot, ftop = fk[:,0:1], fk[:,llm-1:llm]
+        fl = .5*(fk[:,1:]+fk[:,:-1]) # size llm-1
+        return np.array(np.hstack((fac*fbot,fl,fac*ftop))) # size llm+1
+    def mk(self, fl, fac): # fl has size llm+1, last index llm
+        llm=self.llm
+        fbot = fac*fl[:,0:1] + .5*fl[:,1:2]
+        ftop = fac*fl[:,llm:llm+1] + .5*fl[:,llm-1:llm]
+        fk = .5*(fl[:,2:-1]+fl[:,1:-2]) # size llm-2
+        return np.array(np.hstack((fbot,fk,ftop)))  # size llm
+    # Averaging of extensive fields
+    def ml_ext(self, fk): return self.ml(fk,.5)
+    def mk_ext(self, fl): return self.mk(fl,1.)
+    # Averaging of intensive fields
+    def ml_int(self, fk): return self.ml(fk,1.)
+    def mk_int(self, fl): return self.mk(fl,.5)
+    
+############################## Slice #############################
+
+class Slice:
+    def __init__(self, xmesh, zmesh):
+        self.xmesh, self.zmesh = xmesh, zmesh
+        self.xik = xmesh.xi*zmesh.onesk
+        self.xjk = xmesh.xj*zmesh.onesk
+        self.xil = xmesh.xi*zmesh.onesl
+        self.xjl = xmesh.xj*zmesh.onesl
+        self.etaik = xmesh.onesi*zmesh.etak
+        self.etail = xmesh.onesi*zmesh.etal
+        self.etajk = xmesh.onesj*zmesh.etak
+        self.etajl = xmesh.onesj*zmesh.etal
+        self.ones_ik=xmesh.onesi*zmesh.onesk
+    def field_ik(self):
+        return self.zmesh.fieldk(self.xmesh.Nx)
+    def field_il(self):
+        return self.zmesh.fieldl(self.xmesh.Nx)
+    def field_jk(self):
+        return self.zmesh.fieldk(self.xmesh.Nx)
+    def field_jl(self):
+        return self.zmesh.fieldl(self.xmesh.Nx)
+    def cumsum_il(self, f0, dkf):
+        f_il = np.hstack((f0, dkf))
+        return np.cumsum(f_il,axis=1)
+
+########################## Thermodynamics #########################
+
+class Gas:
+    def dp(self, dvol_ik):
+        """ Variations of p,T,g at constant s """
+        dpik = -self.c2*dvol_ik/(self.v*self.v) # dp/dv = -(1/v2)dp/drho
+        dTik = -self.dpds*dvol_ik # dp/ds = -d2e/dvds = -dT/dv
+        dgik = self.v*dpik-self.s*dTik # dG = vdp - sdT
+        return dpik, dTik, dgik
+    
+class Ideal_perfect:
+    def __init__(self,Cpd,Rd,p0,T0):
+        # make sure everything is floating-point
+        self.Cpd = 1.*Cpd
+        self.Cvd = 1.*Cpd-Rd
+        self.Rd  = 1.*Rd
+        self.T0  = 1.*T0
+        self.p0  = 1.*p0
+        self.v0  = Rd*T0/(1.*p0)
+    def set_ps(self,p,s):
+        Cpd = self.Cpd ; Rd = self.Rd
+        T = self.T0*np.exp( (s + Rd*np.log(p/self.p0))/Cpd )
+        return self.energies(p,Rd*T/p,T,s)
+    def set_vs(self,v,s):
+        Cvd = self.Cvd ; Rd = self.Rd
+        T = self.T0*np.exp( (s - Rd*np.log(v/self.v0))/Cvd )
+        return self.energies(Rd*T/v,v,T,s)
+    def set_pT(self,p,T):
+        Cpd = self.Cpd ; Rd = self.Rd
+        s = Cpd*np.log(T/self.T0) - Rd*np.log(p/self.p0)
+        return self.energies(p,Rd*T/p,T,s)
+    def set_vT(self,v,T):
+        Cvd = self.Cvd ; Rd = self.Rd
+        s = Cvd*np.log(T/self.T0) + Rd*np.log(v/self.v0)
+        return self.energies(Rd*T/v,v,T,s)
+    def energies(self,p,v,T,s):
+        gas = Gas();
+        gas.p, gas.v, gas.T, gas.s = p, v, T, s
+        gas.dpds = p/self.Cvd
+        gas.c2=((self.Cpd/self.Cvd)*self.Rd)*T
+        gas.e = self.Cvd*T
+        gas.h = self.Cpd*T
+        gas.g = gas.h-T*s
+        return gas
+
+################################# Dynamics ##############################
+
+class BoundCond:
+    def __init__(self, ptop, Phi_bot, pbot, rho_bot):
+        self.ptop, self.Phi_bot = ptop, Phi_bot # standard BC
+        self.pbot, self.rho_bot = pbot, rho_bot # soft ground
+        self.rigid_bottom, self.rigid_top = True, True # by default
+        self.sponge = lambda Phi,ujk,dujk,tau : (ujk+tau*dujk,dujk)
+        
+class Dynamics:
+    def __init__(self, thermo, metric, BC, fac=(1,1,1)):
+        self.thermo, self.metric = thermo, metric
+        self.BC, self.fac = BC, fac
+        self.di, self.dj = metric.di, metric.dj
+        self.dk, self.dl = metric.dk, metric.dl
+        self.mi, self.mj = metric.mi, metric.mj
+        self.mk_ext, self.ml_ext = metric.mk_ext, metric.ml_ext
+        self.mk_int, self.ml_int = metric.mk_int, metric.ml_int
+    def kinetic(self,mik,ujk,Wil,Phi_il):
+        djPhil = self.dj(Phi_il)
+#        ujk, Wil, djPhil = self.mult_by_fac((ujk, Wil, djPhil))
+        return self.metric.kinetic(mik, ujk, Wil, djPhil)
+    def time_steps(self,gas,mik) :
+        dx = self.metric.dx
+        cmax, dz = np.sqrt(gas.c2.max()), gas.v*mik/dx
+        print 'Min/max layer thickness (m)', dz.min(), dz.max()
+        dz=dz.min()
+        dt_hydro=.5*dx/cmax
+        dt_NH=.5/cmax/math.sqrt(dx**-2+dz**-2)
+        print 'Horizontal resolution (m)', dx
+        print 'Max temperature (K), sound speed (m/s)', gas.T.max(), cmax
+        print 'Courant number 1 times steps for HPE and NH (s) : ', dt_hydro, dt_NH
+        return dt_NH, dt_hydro
+
+############################## Hydrostatic dynamics ###########################
+
+class Shallow_hydro: # hydrostatic routines of the Shallow class defined at the end
+    def hydrostatic_pressure(self, ptop, mik):
+        M_il=np.cumsum(mik,1)
+        Mtop_il=col_vec(M_il[:,-1])*self.zmesh.onesk
+        M_ik = Mtop_il - M_il + .5*mik
+        return ptop + M_ik*self.g0_dx
+    def geopotential(self, Phi0_i, mik, gas):
+        dkphi = mik*gas.v*self.g0_dx
+        phi_il = self.mesh.cumsum_il(Phi0_i, dkphi)
+        return phi_il
+    def hydrostatic_adjustment(self, thermo, BC, mik, sik):
+        pik=self.hydrostatic_pressure(BC.ptop, mik)
+        gas=thermo.set_ps(pik, sik)
+        Phi_il=self.geopotential(BC.Phi_bot, mik, gas)
+        return Phi_il, gas
+    def dK_hydro(self, mik, ujk):
+        """Partial derivatives of HPE kinetic energy"""
+        Ujk=self.mj(mik)*ujk/self.dx2
+        Kik=(.5/self.dx2)*self.mi(ujk*ujk)
+        return Kik, Ujk
+    def inicond(self, hybrid, thermo, ps, T, u):
+        """Hydrostatic initialization given
+        surface pressure, temperature, horiz vel""" 
+        Mi=(ps(self.xi)-self.ptop)*self.dx_g0
+        mik=hybrid.mik(Mi)
+        pik=self.hydrostatic_pressure(mik)
+        Tik=T(self.xik, pik)
+        gas=thermo.set_pT(pik, Tik)
+        Phi_il=self.geopotential(mik, gas)
+        pjk=self.mj(pik)
+        ujk=self.dx*u(self.xjk,pjk)
+        return Mi, mik, gas, Phi_il, ujk
+
+class Hydrostatic(Dynamics):
+    def __init__(self,  thermo, metric, BC):
+        Dynamics.__init__(self, thermo, metric, BC, (1,0,0))
+    def max_time_step(self,gas,mik):
+         dt_NH, dt_hydro = self.time_steps(gas,mik)
+         return dt_hydro, dt_hydro
+    def diagnose(self,mik,Sik,ujk,Phi_il,Wil):
+        Phi_il, gas_ik = self.metric.hydrostatic_adjustment(
+                                    self.thermo, self.BC, mik, Sik/mik)
+        return mik,gas_ik,ujk,Phi_il,0*Phi_il
+    def dH(self, mik, sik, ujk):
+        Phi_il, gas_ik = self.metric.hydrostatic_adjustment(
+                                    self.thermo, self.BC, mik, sik)
+        Kik, Ujk = self.metric.dK_hydro(mik, ujk)
+        Bik = Kik + self.mk_int(Phi_il) + gas_ik.g
+        return Bik, gas_ik, Ujk, Phi_il
+    def bwd_fast_slow(self, mik, Sik, ujk, tau):
+        sik = Sik/mik
+        Phi_il, gas_ik = \
+            self.metric.hydrostatic_adjustment(self.thermo, self.BC, mik, sik)
+        sjk = self.mj(sik)
+        Bik = self.mk_int(Phi_il) + gas_ik.g
+        dujk_fast = -self.dj(Bik)-sjk*self.dj(gas_ik.T)
+        ujk, dujk_fast = self.BC.sponge(Phi_il, ujk, dujk_fast, tau) # updates ujk, dujk_dt
+        Kik, Fjk = self.metric.dK_hydro(mik, ujk)
+        dSik = -self.di(sjk*Fjk)
+        dujk = -self.dj(Kik)
+#        print 'Linf(du/dt) (fast,slow) : ', Linf(dujk_fast),Linf(dujk)
+        return ujk, dujk_fast, Fjk, dSik, dujk 
+
+############################## NH dynamics ###########################
+
+class Shallow_NH: # NH routines of the Shallow class defined at the end
+    def wil(self, mik, Wil):
+        mil=self.ml_ext(mik)
+        return mil, Wil/mil
+    def lambda_il(self,djPhil):
+        return self.mi(djPhil*djPhil)/self.dx2
+    def kinetic(self, mik, ujk, Wil, djPhil):
+        """NH kinetic energy - coincides with HPE if Wil=0"""
+        dx2,mi,mj,ml_int=self.dx2, self.mi, self.mj, self.ml_int
+        mil,wil = self.wil(mik,Wil)
+        lambda_il = (self.lambda_il(djPhil) + self.g2)*wil
+        Kjk = (.5/dx2)*mj(mik)*ujk**2
+        Kjl = (1./dx2)*djPhil*ml_int(ujk)*mj(Wil)
+        Kil = .5*lambda_il*Wil
+        return Kjk.sum()-Kjl.sum()+Kil.sum()
+    def dK_NH_horiz(self, mik, ujk, Wil, djPhil):
+        """Partial derivatives of horizontal kinetic energy"""
+        dx2, mi, mj = self.dx2, self.mi, self.mj
+        mk_ext, mk_int, ml_int = self.mk_int, self.mk_int, self.ml_int
+        mil, wil = self.wil(mik,Wil)
+        lambda_il = wil*self.lambda_il(djPhil)
+        djPhi_W = mk_ext(djPhil*mj(Wil))
+        dK_dujk = (mj(mik)*ujk-djPhi_W)/dx2
+        dK_dmik = (.5/dx2)*mi(ujk*ujk) - .5*mk_int(lambda_il*wil)
+        ujl = ml_int(ujk)
+        dK_dWil = lambda_il - mi(djPhil*ujl)/dx2
+        dK_djPhil = (djPhil*mj(Wil*wil)-mj(Wil)*ujl)/dx2
+        return dK_dmik, dK_dujk, dK_dWil, dK_djPhil
+    def dK_NH_vert(self, mik, Wil):
+        """Partial derivatives of vertical kinetic energy,  
+        also returns mil,wil to be used by HEVI solver"""
+        mil, wil = self.wil(mik,Wil)
+        dK_dWil = self.g2*wil
+        dK_dmik = -.5*self.mk_int(dK_dWil*wil)
+        return dK_dmik, dK_dWil, mil, wil
+    def inicond_NH(self, thermo, Phi, p, T):
+        """ Initialization given Phi(x,eta), p(x,Phi), T(x,p) """
+        Phi_il = Phi(self.xil, self.mesh.etail)
+        Phi_ik = Phi(self.xik, self.mesh.etaik)
+        pik = p(self.xik, Phi_ik)
+        gas = thermo.set_pT(pik, T(self.xik,pik))
+        vol, Jac = self.volume(Phi_il)
+        mik = self.dk(vol) / gas.v
+        return mik, gas, Phi_il
+    def solver_NH(self, BC, dt, mik,Jac_il,gas_ik,mil):
+        dx_g0, g2 = self.dx_g0, self.g2
+        Ak = (dt*dx_g0/gas_ik.v)**2
+        Ak = gas_ik.c2/mik*Ak
+        ml_g2 = mil/g2
+        Bl = 2*self.ml_ext(Ak) + ml_g2
+        Bl[:,0:1] = Bl[:,0:1] + (dx_g0*dt*dt)*BC.rho_bot
+        return time_step.SymTriDiag(Bl,Ak), ml_g2
+
+class Nonhydro(Dynamics):
+    def max_time_step(self,gas,mik):
+         dt_NH, dt_hydro = self.time_steps(gas,mik)
+         return dt_hydro, dt_hydro
+    def diagnose(self,mik,Sik,ujk,Phi_il,Wil):
+        vol_il,dvol_il = self.metric.volume(Phi_il)
+        gas_ik = self.thermo.set_vs(self.dk(vol_il)/mik, Sik/mik)
+        return mik,gas_ik,ujk,Phi_il,Wil
+    def pressure(self, mik,sik,Phi_il):
+        BC, metric = self.BC, self.metric
+        vol_il, Jac_il = metric.volume(Phi_il)
+        gas_ik = self.thermo.set_vs(self.dk(vol_il)/mik, sik)
+        pbot = BC.pbot + BC.rho_bot*(BC.Phi_bot-Phi_il[:,0:1])
+        ptop = BC.ptop*metric.mesh.ones_ik[:,0:1]
+        return gas_ik, Jac_il, self.dl(pbot, gas_ik.p, ptop)
+    def bwd_fast_slow(self, mik,Sik,ujk,Phi_il,Wil, tau, verbose=False):
+        di, dj, dk, dl = self.di, self.dj, self.dk, self.dl
+        BC, metric = self.BC, self.metric
+        sik = Sik/mik
+        mil, wil = metric.wil(mik,Wil)
+        # solve Phi_il - (tau**2)*X.dl(pk) = Phi_star_il
+        if tau>0.:
+            Phi_star_il = Phi_il + tau*metric.g2*(wil - tau)
+            # Newton iteration : Phi_il contains current guess value
+            for ii in range(2):
+                gas_ik, Jac_il, dl_pi = self.pressure(mik,sik,Phi_il)
+                tridiag, ml_g2 = metric.solver_NH(BC, tau, mik,Jac_il,gas_ik,mil)
+                residual = (ml_g2*(Phi_il-Phi_star_il)
+                            +tau*tau*Jac_il*dl_pi)
+                dPhi_il = tridiag.solve(residual)
+                if verbose :
+                    print 'NH solver : Linf(dPhi_il)=', abs(dPhi_il).max(axis=0)
+                    print 'NH solver : Linf(residual)=', abs(residual).max(axis=0)
+                Phi_il = Phi_il - dPhi_il
+        # update Wil
+        gas_ik, Jac_il, dl_pi = self.pressure(mik,sik,Phi_il)
+        dH_Phiil = mil + Jac_il*dl_pi
+        Wil = Wil - tau*dH_Phiil
+        # update ujk
+        dK_mik, dPhiil_fast, mil, wil = metric.dK_NH_vert(mik, Wil)
+        dH_mik = dK_mik + self.mk_int(Phi_il) + gas_ik.g
+        sjk = self.mj(sik)
+        dujk_fast = -dj(dH_mik)-sjk*dj(gas_ik.T)
+        ujk, dujk_fast = self.BC.sponge(Phi_il, ujk, dujk_fast, tau)
+        # slow trends
+        djPhil = dj(Phi_il)
+        Bik, Fjk, dPhiil_slow, dK_djPhil = self.metric.dK_NH_horiz(mik, ujk, Wil, djPhil)
+        dSik_slow = -di(sjk*Fjk)
+        dujk_slow = -dj(Bik)
+        dWil_slow = di(dK_djPhil)
+        return  ( ujk,Phi_il,Wil,                               # updated flow state
+                dujk_fast,dPhiil_fast,-dH_Phiil,                # fast tendencies
+                Fjk,dSik_slow,dujk_slow,dPhiil_slow,dWil_slow ) # slow tendencies
+
+################## Fully equipped shallow-atmosphere metric  ##################
+
+class Shallow(Shallow_hydro, Shallow_NH):
+    def __init__(self, mesh, dx, g0):
+        self.dx, self.g0  = dx, g0
+        self.dx2, self.g2 = dx*dx, g0*g0
+        self.dx_g0, self.g0_dx = dx/g0, g0/dx
+        self.mesh, xmesh, zmesh = mesh, mesh.xmesh, mesh.zmesh
+        self.xmesh, self.zmesh = xmesh, zmesh
+        self.di, self.dj = xmesh.di, xmesh.dj
+        self.mi, self.mj = xmesh.mi, xmesh.mj
+        self.dk, self.dl = zmesh.dk, zmesh.dl
+        self.mk_int, self.ml_int = zmesh.mk_int, zmesh.ml_int
+        self.mk_ext, self.ml_ext = zmesh.mk_ext, zmesh.ml_ext
+        self.xik = mesh.xik*dx
+        self.xjk = mesh.xjk*dx
+        self.xjl = mesh.xjl*dx
+        self.xil = mesh.xil*dx
+        self.xi  = xmesh.xi*dx
+    def volume(self, phi):  # V(Phi), J=dV/dPhi
+        return phi*self.dx_g0, self.dx_g0
+
+######################## Lagrangian vertical coordinate #####################
+
+class Coord:
+    def __init__(self, dyn):
+        self.dyn, self.metric = dyn, dyn.metric
+        
+class LagrangianHydro(Coord):
+    def bwd_fast_slow(self, flow, tau):
+        dyn=self.dyn
+        mik, Sik, ujk, Phiil, Wil, Mjk, Mil =  flow # Phiil, Wil ignored
+        ujk, dujk_fast, Fjk, dSik, dujk = dyn.bwd_fast_slow(mik, Sik, ujk, tau)
+        flow = (mik, Sik, ujk, 0., 0., Mjk, Mil)
+        fast = (0., 0., dujk_fast, 0., 0., 0., 0.)
+        dmik = -dyn.di(Fjk)
+        slow = (dmik, dSik, dujk, 0., 0., Fjk, 0.)
+        return flow, fast, slow
+
+class LagrangianNH(Coord):
+    def bwd_fast_slow(self, flow, tau):
+        dyn=self.dyn
+        mik, Sik, ujk, Phiil, Wil, Mjk, Mil =  flow
+        ujk, Phiil, Wil, \
+            dujk_fast, dPhiil_fast, dWil_fast, \
+            Fjk, dSik, dujk, dPhiil, dWil = dyn.bwd_fast_slow(mik, Sik, ujk, Phiil, Wil, tau)
+        flow = (mik, Sik, ujk, Phiil, Wil, Mjk, Mil)
+        fast = (0., 0., dujk_fast, dPhiil_fast, dWil_fast, 0., 0.)
+        dmik = -dyn.di(Fjk)
+        slow = (dmik, dSik, dujk, dPhiil, dWil, Fjk, 0.)
+        return flow, fast, slow
+
+######################## Mass-based vertical coordinate #####################
+
+class GeneralCoord(Coord):
+    # this class computes the terms due to vertical transport (eta_dot)
+    # given the vertical mass flux Fil
+    def trend_vert_hydro(self, mik,Sik,ujk, Fil):
+        dyn=self.dyn
+        sik=Sik/mik
+        sil=dyn.ml_int(sik) # centered
+        dSik_dt = - dyn.dk(Fil*sil)
+        dluj = dyn.dl(0.,ujk,0.)
+        dujk_dt = - dyn.mk_ext(dluj*dyn.mj(Fil))/dyn.mj(mik)
+        return dSik_dt, dujk_dt
+    def trend_vert_NH(self, mik,Sik,ujk,Phi_il,Wil, Fil):
+        dyn=self.dyn
+        dSik, dujk = self.trend_vert_hydro(mik,Sik,ujk, Fil)
+        Fik = dyn.mk_int(Fil)
+        wik = dyn.mk_ext(Wil)/mik 
+        fil = Fil/dyn.ml_ext(mik) # fil = eta_dot
+        dujk = dujk - dyn.mj(wik*dyn.dk(Phi_il))*dyn.dj(Fik/mik)
+        dWil = - dyn.dl(0.,Fik*wik,0.) # centered vertical flux of vertical momentum
+        dPhiil = - fil*dyn.dl(0., dyn.mk_int(Phi_il), 0.)
+        return dSik, dujk, dPhiil, dWil
+
+class MassCoord(GeneralCoord):
+    # computes the vertical mass flux for a mass-based coordinate
+    def massflux(self,mik,Fjk):
+        def sum(fik): return col_vec(fik.sum(1))*onesk
+        mesh=self.metric.mesh
+        onesk=mesh.zmesh.onesk
+        # mik = Aik*Mi
+        Mik=sum(mik)
+        Aik = mik/Mik
+        # dmik_mass = Aik*dMi = Aik*sumk(dmik_lag)
+        dmik_lag = -self.dyn.di(Fjk)
+        dMik = sum(dmik_lag)
+        dmik_mass = Aik*dMik 
+        #    dmik_lag  + di(Fjk)            = 0
+        #    dmik_mass + di(Fjk) + dk(Fil)  = 0
+        # => dmik_mass - dmik_lag + dk(Fil) = 0
+        dkFil = dmik_lag - dmik_mass
+        Fil = mesh.cumsum_il(mesh.xmesh.zeroi, dkFil)
+        return Fil, dmik_mass
+
+class MassBasedHydro(MassCoord):
+    def bwd_fast_slow(self, flow, tau):
+        mik, Sik, ujk, Phiil, Wil, Mjk, Mil =  flow # Phiil, Wil ignored
+        ujk, dujk_fast, Fjk, dSik, dujk = self.dyn.bwd_fast_slow(mik, Sik, ujk, tau)
+        flow = (mik, Sik, ujk, 0., 0., Mjk, Mil)
+        fast = (0., 0., dujk_fast, 0., 0., 0., 0.)
+        Fil, dmik = self.massflux(mik,Fjk)
+        dSik_vert, dujk_vert = self.trend_vert_hydro(mik,Sik,ujk, Fil)
+        slow = (dmik, dSik+dSik_vert, dujk+dujk_vert, 0., 0., Fjk, Fil)
+        return flow, fast, slow
+
+class MassBasedNH(MassCoord):
+    def bwd_fast_slow(self, flow, tau, verbose=False):
+        dyn=self.dyn
+        mik, Sik, ujk, Phiil, Wil, Mjk, Mil =  flow
+        ujk,Phiil,Wil, \
+            dujk_fast,dPhiil_fast,dWil_fast, \
+            Fjk, dSik, dujk, dPhiil, dWil = dyn.bwd_fast_slow(mik,Sik,ujk,Phiil,Wil, tau, verbose)
+        flow = (mik, Sik, ujk, Phiil, Wil, Mjk, Mil)
+        fast = (0.,0., dujk_fast,dPhiil_fast,dWil_fast, 0.,0.)
+        Fil, dmik = self.massflux(mik,Fjk)
+        dSik_vert, dujk_vert, dPhiil_vert, dWil_vert = self.trend_vert_NH(mik,Sik,ujk,Phiil,Wil,Fil)
+        slow = (dmik, dSik+dSik_vert, dujk+dujk_vert, dPhiil+dPhiil_vert, dWil+dWil_vert, Fjk, Fil)
+        return flow, fast, slow
Index: /codes/icosagcm/devel/Python/dynamico/kessler.py
===================================================================
--- /codes/icosagcm/devel/Python/dynamico/kessler.py	(revision 615)
+++ /codes/icosagcm/devel/Python/dynamico/kessler.py	(revision 615)
@@ -0,0 +1,27 @@
+#------------- import and wrap Kessler routines and variables -------------#
+
+from ctypes import *
+lib = cdll.LoadLibrary("Kessler/libkessler.so")
+
+libvars={}
+def addvars(ctype,names):
+    for name in names:
+        libvars[name]=ctype.in_dll(lib, name)
+def getvar(name):
+    return libvars[name].value
+def setvar(name,val):
+    libvars[name].value=val
+def setvars(names,vals):
+    for name,val in zip(names,vals):
+        libvars[name].value=val
+
+addvars(c_int,('nx','nz','nz1'))
+addvars(c_double,('g', 'p0', 'dt', 'cp', 'r', 'dt', 'time', 'tsfc', 'qc'))
+
+ptr_arg = POINTER(c_int)
+def loc(array): return array.ctypes.data_as(ptr_arg)
+def locs(arrays): return [loc(v) for v in arrays]
+def wrap(routine): return lambda *args : routine(*locs(args))
+
+create_moist_sounding = wrap(lib.create_moist_sounding)
+kessler = wrap(lib.kessler)
Index: /codes/icosagcm/devel/Python/dynamico/parallel.py
===================================================================
--- /codes/icosagcm/devel/Python/dynamico/parallel.py	(revision 615)
+++ /codes/icosagcm/devel/Python/dynamico/parallel.py	(revision 615)
@@ -0,0 +1,167 @@
+from itertools import groupby
+import numpy as np
+
+def inverse_list(lst): return {j:i for i,j in enumerate(lst)}
+def ordered_list(lst): 
+        lst = sorted(zip( lst, range(len(lst)) ))
+        lst, order = zip(*lst)
+        order = inverse_list(order)
+        order = [order[i] for i in range(len(lst))] # dict->list
+        return lst,order  
+    
+#----------------------- Parallel access to distributed array ------------------------------#
+
+# The classes PDim, Get_Indices and PArrayX are used to read NetCDF arrays in parallel
+# Each MPI process reads a contiguous chunk of the array once.
+# Classes CstPArrayX emulate constant arrays that are not stored in a NetCDF file.
+# With class LocArray1D, local data is copied from a user-provided array, rather than from a NetCDF file.
+# Using the class Get_Indices, each MPI process can obtain data for an arbitrary list of indices
+# Data is not re-read from file but scattered via MPI all2all
+
+class PDim:
+    def __init__(self, ncdim, comm):
+        mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size()
+        self.ncdim, n = ncdim, len(ncdim)
+        vtxdist = [(n*i)/mpi_size for i in range(mpi_size+1)]
+        self.n, self.vtxdist = n, np.asarray(vtxdist, dtype=np.int32)
+        self.comm, self.start, self.end = comm, vtxdist[mpi_rank], vtxdist[mpi_rank+1]
+    def getter(self, index_list): return Get_Indices(self, index_list)
+    def get(self, index_list, data):
+        getter=self.getter(index_list)
+        return getter(data)
+class Get_Indices:
+    def __init__(self, dim, index_list): # list MUST be a list of unique indices
+        index_list, self.order = ordered_list(index_list)
+        comm, vtxdist = dim.comm, dim.vtxdist
+        list_local = [ [ j-vtxdist[i] for j in index_list if j>=vtxdist[i] and j<vtxdist[i+1]] 
+                      for i in range(len(vtxdist)-1) ]
+        list_send = comm.alltoall(list_local)
+        self.comm, self.n, self.list_n = comm, len(index_list), [len(i) for i in list_local]
+        self.list_local, self.list_send, self.dict = list_local, list_send, inverse_list(index_list)
+        self.vtxdist=vtxdist
+    def get(self, get_data, put_data):
+        data_send = [get_data(indices) for indices in self.list_send]
+        data_recv = self.comm.alltoall(data_send)
+        start=0
+        for data,n in zip(data_recv,self.list_n):
+            if n>0:
+                end = start + n
+                put_data(start,end,data)
+                start=end
+    def __call__(indices, self): # self is a PArrayND
+        shape = list(self.data.shape)
+        shape[0] = indices.n
+        self.data_out = np.zeros(shape, dtype=self.data.dtype)
+        indices.get(self.get_data, self.put_data )
+#        return self.data_out
+        return self.reorder(indices.order)
+
+class PArrayND:
+    def init_data(self, dim, data):
+        self.dim, self.data = dim, np.asarray(data) # local chunk of data
+    def max(self):
+        max_loc = self.data.max()
+        self.dim.comm.allreduce
+
+class PArray1D(PArrayND):
+    def __init__(self, dim, data): self.init_data(dim, data[dim.start:dim.end])
+    def get_data(self, indices): return np.array(self.data[indices])
+    def put_data(self, start, end, data): self.data_out[start:end]=data
+    def reorder(self, order): return self.data_out[order]
+class LocPArray1D(PArray1D):
+    def __init__(self, dim, data): self.init_data(dim,data)
+class CstPArray1D(PArray1D):
+    def __init__(self, dim, dtype, val): 
+        self.init_data(dim,np.full( (dim.end-dim.start,), val, dtype=dtype))
+
+class PArray2D(PArrayND):
+    def __init__(self, dim, data): self.init_data(dim, data[dim.start:dim.end,:])
+    def get_data(self, indices): return self.data[indices,:]
+    def put_data(self, start, end, data): self.data_out[start:end]=data
+    def reorder(self, order): return self.data_out[order,:]
+
+def PArray(dim, data): return {1:PArray1D, 2:PArray2D}[len(data.shape)](dim,data)
+
+#-------------------------------------- Halo management -----------------------------------#
+
+# Classes LocalDim, LocalArrayX and Halo_Xchange are used to 
+# store local arrays with halos and update halos, respectively
+# A LocalDim instance is created using a list of (global indices of) cells. The instance contains
+# a lookup table associating a local index to a global index. 
+# A Halo_Xchange instance is created first based on a LocalDim and a partitioning table giving MPI process owning each cell.
+# This instance can the be used to create instances of Halo_Xchange. 
+# LocalArray data can be initialized by reading from a PArray or a NumPy array containing local values. 
+# The update() methods updates the halo using point-to-point MPI communications (send/recv) 
+
+class LocalDim:
+    def __init__(self, dim, cells):
+        # dim : a PDim instance ; cells : a list of (global indices of) cells, in any order
+        self.dim, self.loc2glob, self.glob2loc = dim, cells, inverse_list(list)
+        
+def dict2list(size, thedict, default):
+    lst=[default]*size
+    for rank,val in thedict.items() : lst[rank]=val
+    return lst
+
+class Halo_Xchange:
+    def __init__(self, tag, dim, cells, part, reorder_cells=False): 
+        # cells = global indices, part = rank owning each cell
+        # if reorder_cells is True, cells will be reordered so that each halo is contiguous
+        # reordering is not recommended if cells is a telescopic sum of increasing sets C0<C1<C2 ...
+        comm = dim.comm
+        mpi_size, mpi_rank = comm.Get_size(), comm.Get_rank()
+        self.comm, self.tag, self.reorder_cell = comm, tag, reorder_cells
+        # sort cells by rank, keep track of their global and local index
+        recv = sorted(zip(part,cells,range(len(cells))))
+        # group by MPI rank ; lst is a list of (rank,global,local) tuples so zip(*it) is a tuple of lists (rank,global,local)
+        recv = { rank : zip(*lst) for rank,lst in groupby(recv,lambda x:x[0]) }
+        # remove mpi_rank from dict and get list of own cells
+        junk, own_global, own_local = recv.pop(mpi_rank)
+        own_len = len(own_local)
+        if reorder_cells :  own_local = range(own_len)
+        self.own_len, self.own_local, self.get_own = own_len, list(own_local), Get_Indices(dim, own_global) 
+        # figure out data we want to receive from other CPUs
+        recv_rank = sorted(recv.keys())
+        recv_glob = { rank : recv[rank][1] for rank in recv_rank} # global indices of cells to receive
+        recv_len = { rank : len(recv_glob[rank]) for rank in recv_rank}
+        if reorder_cells :
+            recv_loc, start = {}, own_len
+            for rank in recv_rank:
+                end = start+recv_len[rank]
+                recv_loc[rank] = range(start,end)
+                start=end
+            cells = sum([recv_glob[rank] for rank in recv_rank], own_global)
+        else:
+            recv_loc = { rank : recv[rank][2] for rank in recv_rank}
+        self.recv_list = [(rank, list(recv_loc[rank])) for rank in recv_rank]
+        self.cells, self.get_all = cells, Get_Indices(dim, cells)
+        # now figure out the data we must send
+        send=comm.alltoall(dict2list(mpi_size, recv_glob, []))
+        send_len = [len(lst) for lst in send]
+        send_rank = [rank for rank in range(mpi_size) if (send_len[rank]>0)] # CPUs asking for data
+        print 'send_rank %d'%mpi_rank, send_rank
+        # asssociate local index to global index
+        own_dict = { glob:loc for glob,loc in zip(own_global, own_local) }
+        self.send_list = [(rank, [own_dict[i] for i in send[rank]]) for rank in send_rank]                    
+        
+class LocalArray: # a base class for arrays with halos
+    def read_own(self, parray): self.put(self.halo.own_local, self.halo.get_own(parray))
+    def read_all(self, parray): self.data = self.halo.get_all(parray)
+    def update(self):
+        halo=self.halo
+        comm, tag = halo.comm, halo.tag
+        for rank,loc in halo.send_list:
+            comm.send(self.get(loc), dest=rank, tag=tag)
+        for rank,loc in halo.recv_list:
+            self.put(loc, comm.recv(source=rank, tag=tag) )
+        
+class LocalArray1(LocalArray): # a 1D array with halo
+    def __init__(self, halo):
+        self.halo, self.data = halo, np.zeros((len(halo.cells),))
+    def get(self, cells): return np.asarray([self.data[i] for i in cells])
+    def put(self,cells,data): self.data[cells]=data
+class LocalArray2(LocalArray): # a 2D array with halo
+    def __init__(self, halo, llm):
+        self.halo, self.llm, self.data = halo, llm, np.array((halo.local_len,llm))
+    def get(self, cells): return self.data[cells,:]
+    def put(self,cells,data): self.data[cells,:]=data
Index: /codes/icosagcm/devel/Python/dynamico/time_step.py
===================================================================
--- /codes/icosagcm/devel/Python/dynamico/time_step.py	(revision 615)
+++ /codes/icosagcm/devel/Python/dynamico/time_step.py	(revision 615)
@@ -0,0 +1,174 @@
+import numpy as np
+from scipy import sparse as sparse
+import math as math
+
+class Struct: pass
+
+def Linf(data):
+    return abs(data).max()
+def compare(msg, data1, data2):
+    print msg, Linf(data1-data2)/Linf(data1)    
+
+############################# Tridiagonal solver #############################
+
+class SymTriDiag:
+    """ Solve -c(i-1)x(i-1) + b(i)x(i) - c(i)x(i+1) = d(i)
+    Forward sweep :
+    C(i) =        -c(i)       / (b(i)+c(i-1)C(i-1)), C(0)=0
+    D(i) = (d(i)+c(i-1)D(i-1)) / (b(i)+c(i-1)C(i-1)), D(0)=0
+    Back substituion :
+    x(i) = D(i)-C(i)x(i+1), x(N+1)=0
+    """
+    def __init__(self, b,c):
+        M, N = b.shape
+        self.M, self.N, self.b, self.c = M,N,b,c
+        B,C = np.zeros((M,N)), np.zeros((M,N-1))
+        Bi = b[:,0]
+        B[:,0]=Bi
+        for i in xrange(1,N):
+            Ci = -c[:,i-1]/Bi
+            C[:,i-1]=Ci
+            Bi = b[:,i]+c[:,i-1]*Ci
+            B[:,i]=Bi
+        self.B, self.C = B,C
+    def dot(self,x):
+        M,N,b,c = self.M, self.N, self.b, self.c
+        y = np.zeros((M,N))
+        y[:,0] = b[:,0]*x[:,0] - c[:,0]*x[:,1]
+        for i in xrange(1,N-1):
+            y[:,i] = -c[:,i-1]*x[:,i-1] + b[:,i]*x[:,i] - c[:,i]*x[:,i+1]
+        y[:,N-1] = -c[:,N-2]*x[:,N-2] + b[:,N-1]*x[:,N-1]
+        return y
+    def solve(self,d):
+        M,N,c,B,C = self.M, self.N, self.c, self.B, self.C
+        D,x = np.zeros((M,N)), np.zeros((M,N))
+        Di = d[:,0]/B[:,0]
+        D[:,0]=Di
+        for i in xrange(1,N):
+            Di = (d[:,i]+c[:,i-1]*Di)/B[:,i]
+            D[:,i]=Di
+        xi=D[:,N-1]
+        x[:,N-1]=xi
+        for i in xrange(N-2,-1,-1):
+            xi = D[:,i]-C[:,i]*xi
+            x[:,i]=xi
+        return x
+
+########################### ARK time schemes ##########################
+
+class SIRK:
+    """ 
+    Generic SIRK scheme with fast/slow Butcher tableaux af,as
+    Following Weller et al. (2013)
+    sj=slow(yj), fj=fast(yj)
+    zj = y(n) + dt*( sum(l<j)asjl*sl + sum(l<j)afjl*fl )
+    yj = zj + dt*afjj*fast(yj)
+    y(n+1) = y(n)+dt*sum(wsj*sj+wsf*fj)
+    Reformulation : remarking that 
+    y(j)   = y(n) + dt*sum_(l<=j)[as(j,l)*s(l)   + af(j,l)*f(l) ]
+    z(j+1) = y(n) + dt*sum_(l<=j)[as(j+1,l)*s(l) + af(j+1,l)*f(l) ]
+           = y(j) + dt*sum_(l<=j)[cs(j,l)*s(l)   + cf(j,l)*f(l) ]
+    where b(j,l)=a(j,l) ; b(nstage,l)=w(l) ; c(j) = b(j+1)-b(j) => c0=b1
+    => implementation
+    z(0) = y(n)
+    DO j=0, nstage-1
+        [y(j), fj, sj] = bwd_fast_slow(z(j),0.)
+        z(j+1) = y(j) + dt*sum_(l<=j)[cs(j,l)*s(l) + cf(j,l)*f(l) ]
+    END DO
+    y(n+1) = z(nstage)
+
+    where z(j) and y(j) are actually stored in a single variable zj
+    and [y,f,s] = bwd_fast_slow(z,tau) 
+        * solves y - tau*fast(y) = z
+        * computes f=fast(y) and s=slow(y)
+    """
+    def cjl(self,b): # converts Butcher tableau b(j,l) to update coefficients c(j,l)
+        return np.array([ b[j+1,:]-b[j,:] for j in range(self.nstage)])
+    def __init__(self,bwd_fast_slow, dt,nstage, bsjl,bfjl):
+        self.bwd_fast_slow = bwd_fast_slow
+        self.dt, self.nstage = dt, nstage
+        self.csjl = dt*self.cjl(bsjl)
+        self.cfjl = dt*self.cjl(bfjl)
+        self.tauj = np.array([dt*bfjl[j,j] for j in range(self.nstage)])
+    def next(self,zj):
+        csjl, cfjl, tauj = self.csjl, self.cfjl, self.tauj
+        fastj, slowj = [], []
+        for j in range(self.nstage): # here zj is z(j)
+            zj, fj, sj = self.bwd_fast_slow(zj, tauj[j])
+            fastj.append(fj)
+            slowj.append(sj)
+            for l in range(0,j+1): # here zj is initially y(j) and finally z(j+1)
+                zj = [z + csjl[j,l]*s + cfjl[j,l]*f for z,s,f in zip(zj, slowj[l], fastj[l]) ]
+        return zj        
+    def advance(self, flow, N):
+        for k in range(N):
+            flow=self.next(flow)
+        return flow
+    
+def SIRK1(bwd_fast_slow, dt): # untested
+    """ Forward-Backward""" 
+    wj  = np.array([1])
+    bsjl = np.array([[0],wj])
+    bfjl = np.array([[1],wj])
+    return SIRK(bwd_fast_slow,dt,1,bsjl,bfjl)
+        
+def ARK2(bwd_fast_slow, dt, a32 = (3.+2.*np.sqrt(2.))/6.):
+    """ ARK2 scheme by Giraldo et al. (2013), noted ARK2(2,3,2) in Weller et al. (2013) """ 
+    gamma = 1.-1./np.sqrt(2.)
+    delta = 1./(2.*np.sqrt(2.))
+    wj  = np.array([delta,delta,gamma])
+    bsjl = np.array([[0,0,0],[2.*gamma,0,0],[1.-a32,a32,0],wj])
+    bfjl = np.array([[0,0,0],[gamma,gamma,0],[delta,delta,gamma],wj])
+    return SIRK(bwd_fast_slow,dt,3,bsjl,bfjl)
+        
+def UJ3(bwd_fast_slow, dt):
+    """ Strang-carryover scheme by Ullrich et al. (2013), noted UJ3(1,3,2) in Weller et al. (2013) """
+    zz   = [0,0,0,0,0,0]
+    wsj  = [0,1/6.,1/6.,2/3.,0,0]
+    h    = [.5,0,0,0,0,0]
+    wfj  = [.5,0,0,0,0,.5]
+    bsjl = np.array([zz,zz,[0,1,0,0,0,0],
+                     [0,.25,.25,0,0,0],wsj,wsj,wsj])
+    bfjl = np.array([zz,h,h,h,h,wfj,wfj])
+    return SIRK(bwd_fast_slow,dt,6,bsjl,bfjl)
+
+######### Explicit RK schemes written as ARK with two identical Butcher tableaus ########
+
+def RK4(bwd_fast_slow, dt):
+    """ 4-stage 4th-order RK """ 
+    bjl = np.array([[0,0,0,0],
+                    [.5,0,0,0],
+                    [0,.5,0,0],
+                    [0.,0.,1.,0],
+                    [1/6.,1/3.,1/3.,1/6.]])
+    return SIRK(bwd_fast_slow,dt,4,bjl,bjl)
+
+######### Simplified explicit RK schemes (2nd-order) ########
+
+def RK_simple(bwd_fast_slow, dt, tau):
+    # take coefficients of simplified RK scheme
+    # y(n) := x + tau*f(y(n-1))
+    # and computes the equivalent Butcher tableau
+    tau=np.array(tau)
+    n = tau.size # number of stages
+    bjl = np.zeros((n+1,n)) # a[stage, coef]
+    for i in range(n):
+        bjl[i+1,i]=tau[i] 
+    return SIRK(bwd_fast_slow,dt,n,bjl,bjl)
+
+def RKn_simple(n, bwd_fast_slow, dt):
+    tau = [1./(n-i) for i in range(n)]
+    return RK_simple(bwd_fast_slow, dt, tau)
+
+def RK_Kinnmark(K, bwd_fast_slow, dt):
+    """ RK schemes with max stability along imaginary axis
+        Kinnmark & Gray, 1984 """
+    ee=np.zeros((K+1,K+1))
+    ee[2,0], ee[2,1], ee[2,2] = 1,1,1
+    ee[3,0], ee[3,1], ee[3,2], ee[3,3] = 1,2,2,2
+    for k in range(4,K+1):
+        ee[k,0]=1
+        for n in range(1,k+1):
+            ee[k,n]=ee[k-2,n]+2*ee[k-1,n-1]
+    tau = [ee[K,K-n]/ee[K,K-n-1]/(K-1) for n in range(K)]
+    return RK_simple(bwd_fast_slow, dt, tau)
Index: /codes/icosagcm/devel/Python/dynamico/wrap.py
===================================================================
--- /codes/icosagcm/devel/Python/dynamico/wrap.py	(revision 615)
+++ /codes/icosagcm/devel/Python/dynamico/wrap.py	(revision 615)
@@ -0,0 +1,73 @@
+#------------------- wrap C/Fortran routines -------------------#                                                                                           
+
+import numpy as np 
+from ctypes import *
+
+class Struct: pass
+
+# converts numpy array to void * for calling C/Fortran routines                                                                                             
+def loc(array): return array.ctypes.data_as(c_void_p)
+
+# defines converter functions depending on type of Python argument
+def noop(x): return x
+py2c = { np.ndarray:loc, np.float64:c_double, float:c_double, int:c_int, 
+         c_void_p:noop, str:noop, c_int:noop }
+
+class Wrap:
+    def __init__(self, fun, check_args, argtypes):
+        fun.restype=None
+        if check_args:
+            fun.argtypes=argtypes
+        self.fun=fun
+        self.convert=None
+    def __call__(self, *args):
+        if self.convert is None:
+            self.convert=[ py2c[type(arg)] for arg in args]
+        args = [ conv(arg) for conv,arg in zip(self.convert,args) ]
+        return self.fun(*args)
+
+def typelist(argtypes):
+    lst=[]
+    for x in argtypes:
+        if type(x) is int:
+            lst+=(x-1)*[prev] # x-1 because prev has already been appended once                                                             
+        else:
+            lst.append(x)
+        prev=x
+    return lst
+
+class SharedLib:
+    def __init__(self,vardict,name,check_args=True, prefix_so='',prefix_py=''):
+        self.vardict, self.check_args = vardict, check_args
+        self.prefix_so, self.prefix_py = prefix_so, prefix_py
+        self.lib=cdll.LoadLibrary(name)
+        self.vars={}
+    def import_funs(self,prototypes,prefix=''):
+        for proto in prototypes:
+            names=[]
+            while type(proto[0]) is str:
+                names.append(proto.pop(0))
+            if proto[0] is None:
+                proto=None
+            else:
+                proto = typelist(proto)
+            for name in names:
+                name = prefix+name
+                soname = self.prefix_so+name
+                print 'wrap.Sharedlib : importing', soname
+                self.vardict[self.prefix_py+name] = Wrap(self.lib[soname], self.check_args, proto)
+    def addvars(self,*types_and_names):
+        for x in types_and_names:
+            if type(x) is str:
+                self.vars[x]=ctype.in_dll(self.lib, x.lower())
+            else:
+                ctype = x
+    def getvar(self,name):
+        return self.vars[name].value
+    def getvars(self,*names):
+        return [self.vars[name].value for name in names]
+    def setvar(self,name,val):
+        self.vars[name].value=val
+    def setvars(self,names,vals):
+        for name,val in zip(names,vals):
+            self.vars[name].value=val
Index: /codes/icosagcm/devel/Python/env/CICLAD.env
===================================================================
--- /codes/icosagcm/devel/Python/env/CICLAD.env	(revision 615)
+++ /codes/icosagcm/devel/Python/env/CICLAD.env	(revision 615)
@@ -0,0 +1,46 @@
+module purge
+module load intel/15.0.6.233
+module load openmpi/1.6.5-ifort
+module load hdf5/1.8.18-parallel-ifort
+module load netcdf4/4.4.1.1-parallel-ifort
+module load python/2.7-anaconda
+
+# If not done yet, certain python packages must be installed in a dynamico-only conda environment :
+#   conda create -n dynamico python
+#   source activate dynamico
+#   pip install scipy mpi4py jinja2 matplotlib cython cmake
+
+# netCDF4 must be installed manually in order to use the netCDF library from the modules (requires Cython, see above)
+# for this download and unpack 
+#   https://pypi.python.org/packages/6a/37/49fd7cc21af0fb173e40de5e15f7fdd48d521429922a90347219c77b0c36/netCDF4-1.2.9.tar.gz#md5=e320491d52c42e937e6df47b56a2579c
+# then
+#   python setup.py build
+#   python setup.py install
+
+# CICLAD does not have a parmetis module
+# compiling parmetis uses cmake, installed above using pip because the system-wide cmake might not work with the modules
+# then download parmetis from 
+#   http://glaros.dtc.umn.edu/gkhome/fetch/sw/parmetis/parmetis-4.0.3.tar.gz
+# unpack, follow workaround given in 
+#   https://stackoverflow.com/questions/35058709/building-shared-parmetis-4-0-3
+# then
+#   make config shared=1 prefix=$HOME/.conda/envs/dynamico
+#   make install
+# manually install libmetis.so and metis.h ...
+#   cp $(find . -name 'libmetis.so') $HOME/.conda/envs/dynamico/lib/
+#   cp $(find . -name 'metis.h') $HOME/.conda/envs/dynamico/include/
+
+source activate dynamico
+PARMETIS_CFLAGS="-I$HOME/.conda/envs/dynamico/include"
+PARMETIS_LDFLAGS="-L$HOME/.conda/envs/dynamico/lib"
+
+export NETCDF_LDFLAGS="$(nf-config --flibs) -lhdf5 -lz -ldl -lm"
+
+# check that packages work
+for x in netCDF4 mpi4py numpy jinja2 ; do python -c "import $x" ; done
+python -c "import ctypes ; parmetis=ctypes.cdll.LoadLibrary('libparmetis.so')"
+
+export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$HOME/.conda/envs/dynamico/lib
+
+BASEDIR=$(dirname ${BASH_SOURCE[0]})
+. $BASEDIR/linux.env
Index: /codes/icosagcm/devel/Python/env/OSX.env
===================================================================
--- /codes/icosagcm/devel/Python/env/OSX.env	(revision 615)
+++ /codes/icosagcm/devel/Python/env/OSX.env	(revision 615)
@@ -0,0 +1,38 @@
+BASEDIR=$(dirname ${BASH_SOURCE[0]})
+export BASEDIR=$(cd $BASEDIR/../.. ; pwd)
+export PYTHONPATH="$BASEDIR"
+export XIOS_ROOT="$BASEDIR/MODELS/XIOS"
+
+export PATH=/opt/local/bin:/opt/local/sbin:/usr/bin:/bin:/usr/sbin:/sbin:/usr/local/bin:/opt/X11/bin:/usr/texbin
+export LD_LIBRARY_PATH="$BASEDIR/kernels/lib"
+
+export PYTHON=/opt/local/bin/python
+export PYTHON_ROOT=/opt/local/Library/Frameworks/Python.framework/Versions/2.7/
+
+export C99=/usr/bin/clang
+export C99FLAGS="-std=c99 -I /opt/local/include/openmpi-gcc48 -I /opt/local/include $CLFAGS"
+
+export F90=gfortran
+export FFLAGS_BASE="-I $XIOS_ROOT/inc -fPIC -g -openmp"
+export FFLAGS_PROD="-O3"
+export FFLAGS_DBG="-O3 -fcheck=bounds"
+export FFLAGS="$FFLAGS_BASE $FFLAGS_PROD"
+
+export LINK=mpicxx
+export LDFLAGS="-fopenmp"
+export SOFLAGS="-dynamiclib -undefined dynamic_lookup -Wl,-headerpad_max_install_names"
+export LIBS="-L /opt/local/lib/openmpi-gcc48/ -L /opt/local/lib -L/opt/local/lib/gcc48 /opt/local/lib/gcc48/libgfortran.a -lgcc_ext.10.5 -lmpi -lparmetis"
+
+export MPIRUN=mpirun-openmpi-gcc48
+
+export WHOLE_ARCHIVE='-Wl,-all_load'
+export NOWHOLE_ARCHIVE=''
+
+export NETCDF_LDFLAGS="-L/opt/local/lib -lnetcdff -lnetcdf -lstdc++"
+export LDD="otool -L"
+
+# for the record
+#export CC=/usr/bin/clang
+#export CFLAGS="-fno-strict-aliasing -fno-common -dynamic -pipe -Os -fwrapv -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -I/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/core/include -I/opt/local/Library/Frameworks/Python.framework/Versions/2.7/include/python2.7"
+
+export CODEGEN_ENV_SET=yes
Index: /codes/icosagcm/devel/Python/env/X64_CURIE.env
===================================================================
--- /codes/icosagcm/devel/Python/env/X64_CURIE.env	(revision 615)
+++ /codes/icosagcm/devel/Python/env/X64_CURIE.env	(revision 615)
@@ -0,0 +1,5 @@
+BASEDIR=$(dirname ${BASH_SOURCE[0]})
+. $BASEDIR/linux.env
+
+# Needed for Intel OpenMP
+KERNELS_LDFLAGS="-liomp5"
Index: /codes/icosagcm/devel/Python/env/linux.env
===================================================================
--- /codes/icosagcm/devel/Python/env/linux.env	(revision 615)
+++ /codes/icosagcm/devel/Python/env/linux.env	(revision 615)
@@ -0,0 +1,42 @@
+#----------- Useful directories ---------
+BASEDIR=$(dirname ${BASH_SOURCE[0]})
+export BASEDIR=$(cd $BASEDIR/../.. ; pwd)
+
+#--------- To create shared libraries -------
+export F90=mpif90
+export LINK=$F90
+export SOFLAGS="-shared"
+export WHOLE_ARCHIVE='-Wl,--whole-archive'
+export NO_WHOLE_ARCHIVE='-Wl,--no-whole-archive'
+export LDD="ldd"
+
+#----------- To compile kernels ------------
+
+export C99=mpicc
+export CFLAGS="-fPIC $PARMETIS_CFLAGS"
+export C99FLAGS="-std=c99 $CFLAGS"
+export FFLAGS_BASE="-I $XIOS_INCDIR -fPIC -g -openmp"
+export FFLAGS_PROD="-O3"
+export FFLAGS_DBG="-O3 -fcheck=bounds"
+export FFLAGS="$FFLAGS_BASE $FFLAGS_PROD"
+
+#----------- To compile with cython ---------
+
+#---------------- To execute ----------------
+
+export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$DYNAMICO_ROOT/lib
+#export PYTHONPATH="$BASEDIR"
+export PYTHON=$(which python)
+#export PYTHON_ROOT=$(basename $(basename $PYTHON))
+
+
+#--------------- Link --------------
+#export LDFLAGS="$PARMETIS_LDFLAGS"
+#export LIBS="-lparmetis"
+
+
+#-------------- Execute --------------
+# disable Intel OpenMP in MKL (loaded by numpy) to avoid conflict with our openMP
+export MKL_THREADING_LAYER=SEQUENTIAL
+
+export CODEGEN_ENV_SET=yes
Index: /codes/icosagcm/devel/Python/jinja
===================================================================
--- /codes/icosagcm/devel/Python/jinja	(revision 615)
+++ /codes/icosagcm/devel/Python/jinja	(revision 615)
@@ -0,0 +1,42 @@
+from jinja2 import Environment, FileSystemLoader
+from sys import argv
+from textwrap import dedent
+from re import search # regular expressions
+
+def pretty(line,indent):
+    line = dedent(line)
+    if search('^#', line): # check if line is a preprocessor directive
+       print line
+    else:
+       print (3*indent)*' ' + line
+    	
+env=Environment(loader=FileSystemLoader('.'))
+tpl=env.get_template(argv[1])
+source = tpl.render()
+source = [line for line in source.split('\n') if line.strip() != '']
+
+# regular expressions matching Fortran keywords
+# \b detects word boundaries
+# $ and ( must be escaped : \$, \(
+# THIS|THAT means THIS or THAT
+# raw Python strings r'blabla' avoid doubling '\'
+noop      = r'!\$OMP'
+increase  = r'\bTHEN\b|\bDO\b|\bSELECT CASE\b'
+inter     = r'\bELSE\b|\bCASE\(|\bCASE DEFAULT\b'
+decrease  = r'\bEND IF\b|\bENDIF\b|\bEND DO\b|\bENDDO\b|\bEND SELECT\b' 
+
+indent=1
+
+for line in source :
+    if search(noop, line) :
+        pretty(line,indent)
+    elif search(decrease, line) :
+        indent = indent-1
+        pretty(line,indent)
+    elif search(increase,line):
+        pretty(line,indent)
+        indent = indent+1
+    elif search(inter,line):
+        pretty(line,indent-1)
+    else:
+        pretty(line,indent)
Index: /codes/icosagcm/devel/Python/src/functions.h
===================================================================
--- /codes/icosagcm/devel/Python/src/functions.h	(revision 615)
+++ /codes/icosagcm/devel/Python/src/functions.h	(revision 615)
@@ -0,0 +1,12 @@
+/* KERNELS */
+void init_params(void);
+
+/* KERNELS -> XIOS */
+void setup_xios(void);
+void call_xios_set_timestep(double);
+void call_xios_update_calendar(int);
+
+/* CXIOS */
+void cxios_context_close_definition(void);
+void cxios_context_finalize(void);
+void cxios_finalize(void);
Index: /codes/icosagcm/devel/Python/src/hexagonal/macros.jin
===================================================================
--- /codes/icosagcm/devel/Python/src/hexagonal/macros.jin	(revision 615)
+++ /codes/icosagcm/devel/Python/src/hexagonal/macros.jin	(revision 615)
@@ -0,0 +1,274 @@
+#ifdef PASS_PRE1
+
+#define _AND_ &&
+#define _OR_ ||
+#define _NOT_ !
+#define IS_INNER_LAYER (_NOT_ (IS_BOTTOM_LEVEL _OR_ IS_TOP_LAYER))
+#define IS_INNER_INTERFACE (_NOT_ (IS_BOTTOM_LEVEL _OR_ IS_TOP_INTERFACE))
+
+#define CELL ij,LL
+#define VERTEX1 ij+{{vertex1}},LL
+#define VERTEX2 ij+{{vertex2}},LL
+#define EDGE ij+{{edge}},LL
+#define EDGE_TRISK ij+{{edge_trisk}},LL
+#define DUAL_CELL ij+{{triangle}},LL
+#define CELL1 ij,LL
+#define CELL2 ij+{{neighbour}},LL
+#define VERTEX {{vertex}},LL
+#define SIGN {{sign}}
+#define LE_DE le_de(ij+{{edge}})
+#define WEE wee(ij+{{edge}},{{itrisk}})
+#define DE de(ij+{{edge}})
+#define RIV2 Riv2({{vertex}},{{vertex_dir}})
+#define AI Ai(ij)
+#define FV fv(ij+{{triangle}})
+#define AV Av(ij+{{triangle}})
+
+#endif
+
+#ifdef PASS_PRE2
+
+#define BARRIER !$OMP BARRIER
+
+#define KERNEL(name) {% call define_kernel(#name) %}
+#define FORALL_CELLS(...) {% call forall_cells(__VA_ARGS__) %}
+#define FORALL_CELLS_EXT(...) {% call forall_cells_ext(__VA_ARGS__) %}
+#define ON_PRIMAL {% call(forall_edges) on_hexagons() %}
+#define ON_DUAL {% call(triangle,forall_edges,forall_vertices) on_triangles() %}
+#define ON_EDGES {% call(edge,neighbour,sign,vertex1,vertex2,forall_trisk) on_edges() %}
+#define FORALL_EDGES {% call(edge,sign,vertex1,vertex2) forall_edges() %}
+#define FORALL_VERTICES {% call(vertex,vertex_dir) forall_vertices() %}
+#define FORALL_TRISK {% call(itrisk,edge_trisk) forall_trisk() %}
+#define END_BLOCK {% endcall %}
+#define SEQUENCE {% call(ext) sequence() %}
+#define SEQUENCE_EXT {% call(ext) sequence_ext() %}
+#define PROLOGUE(level) {% call at_level(level,ext) %}
+#define BODY(range){% call body(range,ext) %}
+#define EPILOGUE(level){% call at_level(level,ext) %}
+
+#define CST_IF(condition, action) {{ cst_if(#condition, #action) }}
+#define CST_IFTHEN(condition) {{flat}}if condition
+#define CST_ELSEIF(condition) {{flat}}elif condition
+#define CST_ELSE {{flat}}else
+#define CST_ENDIF {{flat}}endif
+
+#endif
+
+#if PASS_JINJA
+
+{% macro sequence() -%}
+{{ caller('') }}
+{%- endmacro %}
+
+{% macro sequence_ext() -%}
+{{ caller('_ext') }}
+{%- endmacro %}
+
+{% macro at_level(lev,ext) -%}
+{{ define('LL', lev) }}
+DO ij=ij_omp_begin{{ext}},ij_omp_end{{ext}}
+{{ caller() }}
+END DO
+{%- endmacro %}
+
+{% macro body(range,ext) -%}
+{{ define('LL', 'l') }}
+DO l = {{ range }}
+  DO ij=ij_omp_begin{{ext}},ij_omp_end{{ext}}
+    {{ caller() }}
+  END DO
+END DO
+{%- endmacro %}
+
+{#---------------------- big macro generating nested DO loops -------------------------#}
+{# specialized code is generated for l=1 and l=llm,llm+1 if necessary #}
+{# start/end values for l are replaced by corresponding OpenMP values #}
+
+{% macro define_levels(ll,kdown,kup,thecode) %}
+  {{ define('LL', ll) }}
+  {{ define_cppmacro('KDOWN','(ij,ll)', kdown) if 'KDOWN' in thecode }}
+  {{ define_cppmacro('KUP','(ij,ll)', kup) if 'KUP'   in thecode }}
+{%- endmacro %}
+
+{% macro forall(start,end,ij_range,thecode) %}
+{% set start_omp={'1':'ll_begin' , '2':'ll_beginp1'}[start] %}  
+{% set end_omp={'llm-1':'ll_endm1' , 'llm':'ll_end', 'llm+1':'ll_endp1'}[end] %}
+{% set is_top_layer, is_top_inter = False, False %}
+
+{{ define('IS_TOP_LAYER', '_FALSE_') }}
+{{ define('IS_TOP_INTERFACE', '_FALSE_') }}
+
+{% if 'IS_BOTTOM_LEVEL' in thecode or 'KDOWN' in thecode%}
+{{ 'ERROR : using IS_BOTTOM_LEVEL in a loop starting at l=2' if start=='2' }}
+{# the code in the loop checks whether l==1, for the sake of performance 
+we shall write special code for l=1 and start the loop at l=2 #}
+
+IF (ll_begin==1) THEN
+  {{ define('IS_BOTTOM_LEVEL', '_TRUE_') }}
+  {{ define_levels('1', 'ij,1', 'ij,1', thecode) }}
+   !DIR$ SIMD
+    DO ij={{ij_range}}
+      {{ thecode }}
+    END DO
+END IF
+{% set start_omp='ll_beginp1' %}
+{% endif %}
+
+{% if 'IS_TOP_LAYER' in thecode %}
+{{ 'ERROR : using _IS_TOP_LAYER_ in a loop ending at l=llm+1' if end=='llm+1' }}
+{# the code checks whether l==llm, write special code for l=llm and end the loop at l=llm-1 #}
+{% set end_omp, is_top_layer = 'll_endm1', True %}
+{% endif %}
+
+{% if 'IS_TOP_INTERFACE' in thecode or 'KUP' in thecode %}
+{# the code checks whether l==llm+1, write special code for l=llm+1 and end the loop at l=llm #}
+{{ 'ERROR : using _IS_TOP_INTERFACE_ in a loop ending at l=llm' if end=='llm' }}
+{% set end_omp, is_top_inter ='ll_end', True %}
+{% endif %}
+
+{{ define('IS_BOTTOM_LEVEL', '_FALSE_') }}
+{{ define_levels('l', 'ij,l-1', 'ij,l', thecode) }}
+DO l = {{start_omp}}, {{end_omp}}
+!DIR$ SIMD
+  DO ij={{ij_range}}
+  {{ thecode }}
+  END DO
+END DO
+
+{% if is_top_layer %}
+IF(ll_end==llm) THEN
+  {{ define('IS_TOP_LAYER', '_TRUE_') }}
+  {{ define_levels('llm', 'ij,llm-1', 'ij,llm', thecode) }}
+  !DIR$ SIMD
+  DO ij={{ij_range}}
+    {{ thecode }}
+  END DO
+END IF
+{% endif %}
+
+{% if is_top_inter %}
+IF(ll_endp1==llm+1) THEN
+  {{ define('IS_TOP_INTERFACE', '_TRUE_') }}
+  {{ define_levels('llm+1', 'ij,llm', 'ij,llm', thecode) }}
+  !DIR$ SIMD
+  DO ij={{ij_range}}
+    {{ thecode }}
+  END DO
+END IF
+{% endif %}
+
+{% endmacro %}
+
+{#-------------------------------------------------------------------------------------#}
+
+{% macro forall_cells(start='1', end='llm') -%}
+{{ forall(start,end,'ij_begin, ij_end', caller()) }}
+{% endmacro %}
+
+{% macro forall_cells_ext(start='1', end='llm') -%}
+{{ forall(start,end,'ij_begin_ext, ij_end_ext', caller()) }}
+{% endmacro %}
+
+{% macro forall_edges_hex() -%}
+{{ caller('u_rup','ne_rup', 'z_rup','z_up') }}
+{{ caller('u_lup','ne_lup', 'z_lup','z_up') }}
+{{ caller('u_left','ne_left','z_lup','z_ldown') }}
+{{ caller('u_ldown','ne_ldown', 'z_ldown','z_down') }}
+{{ caller('u_rdown','ne_rdown', 'z_rdown','z_down') }}
+{{ caller('u_right','ne_right', 'z_rup','z_rdown') }}
+{%- endmacro %}
+
+{% macro on_hexagons() -%}
+{{ caller(forall_edges_hex) }}
+{%- endmacro %}
+
+{% macro forall_edges_zup() -%}
+{{ caller('u_rup','ne_rup', 'TODO','TODO') }}
+{{ caller('t_rup+u_left','ne_left', 'TODO','TODO') }}
+{{ caller('u_lup','(-ne_lup)', 'TODO','TODO') }}
+{%- endmacro %}
+
+{% macro forall_edges_zdown() -%}
+{{ caller('u_rdown','(-ne_rdown)', 'TODO', 'TODO') }}
+{{ caller('t_ldown+u_right','ne_right', 'TODO', 'TODO') }}
+{{ caller('u_ldown','ne_ldown', 'TODO', 'TODO') }}
+{%- endmacro %}
+
+{% macro forall_vertices_zup() -%}
+{{ caller('ij','vup') }}
+{{ caller('ij+t_rup','vldown') }}
+{{ caller('ij+t_lup','vrdown') }}
+{%- endmacro %}
+
+{% macro forall_vertices_zdown() -%}
+{{ caller('ij','vdown') }}
+{{ caller('ij+t_ldown','vrup') }}
+{{ caller('ij+t_rdown','vlup') }}
+{%- endmacro %}
+
+{% macro on_triangles() -%}
+{{ caller('z_up',forall_edges_zup,forall_vertices_zup) }}
+{{ caller('z_down',forall_edges_zdown,forall_vertices_zdown) }}
+{%- endmacro %}
+
+{% macro trisk_right() -%}
+{{ caller('1,1', 'u_rup') }}
+{{ caller('2,1', 'u_lup') }}
+{{ caller('3,1', 'u_left') }}
+{{ caller('4,1', 'u_ldown') }}
+{{ caller('5,1', 'u_rdown') }}
+{{ caller('1,2', 't_right+u_ldown') }}
+{{ caller('2,2', 't_right+u_rdown') }}
+{{ caller('3,2', 't_right+u_right') }}
+{{ caller('4,2', 't_right+u_rup') }}
+{{ caller('5,2', 't_right+u_lup') }}
+{%- endmacro %}
+
+{% macro trisk_lup() -%}
+{{ caller('1,1', 'u_left') }}
+{{ caller('2,1', 'u_ldown') }}
+{{ caller('3,1', 'u_rdown') }}
+{{ caller('4,1', 'u_right') }}
+{{ caller('5,1', 'u_rup') }}
+{{ caller('1,2', 't_lup+u_right') }}
+{{ caller('2,2', 't_lup+u_rup') }}
+{{ caller('3,2', 't_lup+u_lup') }}
+{{ caller('4,2', 't_lup+u_left') }}
+{{ caller('5,2', 't_lup+u_ldown') }}
+{%- endmacro %}
+
+{% macro trisk_ldown() -%}
+{{ caller('1,1', 'u_rdown') }}
+{{ caller('2,1', 'u_right') }}
+{{ caller('3,1', 'u_rup') }}
+{{ caller('4,1', 'u_lup') }}
+{{ caller('5,1', 'u_left') }}
+{{ caller('1,2', 't_ldown+u_lup') }}
+{{ caller('2,2', 't_ldown+u_left') }}
+{{ caller('3,2', 't_ldown+u_ldown') }}
+{{ caller('4,2', 't_ldown+u_rdown') }}
+{{ caller('5,2', 't_ldown+u_right') }}
+{%- endmacro %}
+
+{% macro on_edges() -%}
+{{ caller('u_right','t_right','ne_right','z_rdown','z_rup',trisk_right) }}
+{{ caller('u_lup','t_lup','ne_lup','z_up','z_lup', trisk_lup) }}
+{{ caller('u_ldown','t_ldown','ne_ldown','z_ldown','z_down', trisk_ldown) }}
+{%- endmacro %}
+
+{% set llm = 'llm' %}
+
+#endif
+
+#ifdef PASS_POST1
+#define _TRUE_ (0==0)
+#define _FALSE_ (0==1)
+#define CELL ij,LL
+#endif
+
+#ifdef PASS_POST2
+#define UP(ij,l) ij,l+1
+#define DOWN(ij,l) ij,l-1
+#define HIDX(ij,l) ij
+#define VIDX(ij,l) l
+#endif
Index: /codes/icosagcm/devel/Python/src/kernels_caldyn_NH.jin
===================================================================
--- /codes/icosagcm/devel/Python/src/kernels_caldyn_NH.jin	(revision 615)
+++ /codes/icosagcm/devel/Python/src/kernels_caldyn_NH.jin	(revision 615)
@@ -0,0 +1,312 @@
+{% macro compute_p_and_Aik() %}
+{% set p_and_c2_from_theta=caller %}
+    SEQUENCE_EXT
+      BODY('1,llm')
+          rho_ij = (g*m_ik(CELL))/(Phi_il(UP(CELL))-Phi_il(CELL))
+          {{ p_and_c2_from_theta() }}
+          A_ik(CELL) = c2_mik*(tau/g*rho_ij)**2
+      END_BLOCK
+    END_BLOCK
+{%- endmacro %}
+
+KERNEL(compute_NH_geopot)
+  tau2_g=tau*tau/g
+  g2=g*g
+  gm2 = 1./g2
+  gamma = 1./(1.-kappa)
+    
+  BARRIER
+
+  ! compute Phi_star
+  SEQUENCE_EXT
+    BODY('1,llm+1')
+      Phi_star_il(CELL) = Phi_il(CELL) + tau*g2*(W_il(CELL)/m_il(CELL)-tau)
+    END_BLOCK
+  END_BLOCK
+
+  ! Newton-Raphson iteration : Phi_il contains current guess value
+  DO iter=1,2 ! 2 iterations should be enough
+    ! Compute pressure, A_ik
+    SELECT CASE(caldyn_thermo)
+    CASE(thermo_theta)
+      {% call() compute_p_and_Aik() %}
+        X_ij = (cpp/preff)*kappa*theta(CELL)*rho_ij
+        p_ik(CELL) = preff*(X_ij**gamma)
+        c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho
+      {% endcall %}
+    CASE(thermo_entropy)
+      {% call() compute_p_and_Aik() %}
+        X_ij = Treff*exp(theta(CELL)/cpp) ! theta = Tref.exp(s/Cp)
+        X_ij = (cpp/preff)*kappa*X_ij*rho_ij
+        p_ik(CELL) = preff*(X_ij**gamma)
+        c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho
+      {% endcall %}
+    CASE DEFAULT
+        PRINT *, 'caldyn_thermo not supported by compute_NH_geopot', caldyn_thermo
+        STOP
+    END SELECT
+    
+    ! NB : A(1), A(llm), R(1), R(llm+1) = 0 => x(l)=0 at l=1,llm+1 => flat, rigid top and bottom
+    ! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm
+
+    SEQUENCE_EXT
+      ! Compute residual R_il and B_il
+      PROLOGUE(1)
+        ! bottom interface l=1
+        ml_g2 = gm2*m_il(CELL) 
+        B_il(CELL) = A_ik(CELL) + ml_g2 + tau2_g*rho_bot
+        R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  &
+                  + tau2_g*( p_ik(CELL)-pbot+rho_bot*(Phi_il(CELL)-PHI_BOT(HIDX(CELL))) )
+      END_BLOCK
+      BODY('2,llm')
+        ! inner interfaces
+        ml_g2 = gm2*m_il(CELL) 
+        B_il(CELL) = A_ik(CELL)+A_ik(DOWN(CELL)) + ml_g2
+        R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  &
+                  + tau2_g*(p_ik(CELL)-p_ik(DOWN(CELL)))
+        ! consistency check : if Wil=0 and initial state is in hydrostatic balance
+        ! then Phi_star_il(CELL) = Phi_il(CELL) - tau^2*g^2
+        ! and residual = tau^2*(ml+(1/g)dl_pi)=0
+      END_BLOCK
+      EPILOGUE('llm+1')
+        ! top interface l=llm+1
+        ml_g2 = gm2*m_il(CELL) 
+        B_il(CELL) = A_ik(DOWN(CELL)) + ml_g2
+        R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  &
+                  + tau2_g*( ptop-p_ik(DOWN(CELL)) )
+      END_BLOCK
+      !
+      ! Forward sweep :
+      ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)), 
+      ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1))
+      PROLOGUE(1)
+        X_ij = 1./B_il(CELL)
+        C_ik(CELL) = -A_ik(CELL) * X_ij 
+        D_il(CELL) =  R_il(CELL) * X_ij
+      END_BLOCK
+      BODY('2,llm')
+        X_ij = 1./( B_il(CELL) + A_ik(DOWN(CELL))*C_ik(DOWN(CELL)) )
+        C_ik(CELL) = -A_ik(CELL) * X_ij 
+        D_il(CELL) =  (R_il(CELL)+A_ik(DOWN(CELL))*D_il(DOWN(CELL))) * X_ij
+      END_BLOCK
+      EPILOGUE('llm+1')
+        X_ij = 1./( B_il(CELL) + A_ik(DOWN(CELL))*C_ik(DOWN(CELL)) )
+        D_il(CELL) =  (R_il(CELL)+A_ik(DOWN(CELL))*D_il(DOWN(CELL))) * X_ij
+        ! Back substitution :
+        ! x(i) = D(i)-C(i)x(i+1), x(llm+1)=0
+        ! + Newton-Raphson update
+        ! top interface l=llm+1
+        x_il(CELL) = D_il(CELL)
+        Phi_il(CELL) = Phi_il(CELL) - x_il(CELL)
+      END_BLOCK
+      BODY('llm,1,-1')
+        ! Back substitution at lower interfaces
+        x_il(CELL) = D_il(CELL) - C_ik(CELL)*x_il(UP(CELL))
+        Phi_il(CELL) = Phi_il(CELL) - x_il(CELL)
+      END_BLOCK
+    END_BLOCK
+
+    IF(debug_hevi_solver) THEN
+       PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il))
+       PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il))
+       DO l=1,llm+1
+          WRITE(*,'(A,I2.1,I3.2,E9.2)'), '[hevi_solver] x_il', iter,l, MAXVAL(ABS(x_il(l,:)))
+       END DO
+       DO l=1,llm+1
+          WRITE(*,'(A,I2.1,I3.2,E9.2)'), '[hevi_solver] R_il', iter,l, MAXVAL(ABS(R_il(l,:))) 
+       END DO
+    END IF
+
+  END DO ! Newton-Raphson
+
+  BARRIER
+
+  debug_hevi_solver=.FALSE.
+END_BLOCK
+
+KERNEL(caldyn_solver)
+  FORALL_CELLS_EXT('1', 'llm+1')
+    ON_PRIMAL
+       CST_IF(IS_BOTTOM_LEVEL,    m_il(CELL) = .5*rhodz(CELL)                     )
+       CST_IF(IS_INNER_INTERFACE, m_il(CELL) = .5*(rhodz(CELL)+rhodz(DOWN(CELL))) )
+       CST_IF(IS_TOP_INTERFACE,   m_il(CELL) = .5*rhodz(DOWN(CELL))               )
+    END_BLOCK
+  END_BLOCK
+  !
+  IF(tau>0) THEN ! solve implicit problem for geopotential
+    CALL compute_NH_geopot(tau,PHI_BOT_VAR, rhodz, m_il, theta, W, geopot)
+  END IF
+  ! 
+  ! Compute pressure (pres) and Exner function (pk)
+  ! kappa = R/Cp
+  ! 1-kappa = Cv/Cp
+  ! Cp/Cv = 1/(1-kappa)
+  gamma = 1./(1.-kappa)
+  vreff    = Rd*Treff/preff ! reference specific volume
+  Cvd   = cpp-Rd
+  FORALL_CELLS_EXT()
+    SELECT CASE(caldyn_thermo)
+    CASE(thermo_theta)
+      ON_PRIMAL
+        rho_ij = g*rhodz(CELL)/(geopot(UP(CELL))-geopot(CELL))
+        X_ij = (cpp/preff)*kappa*theta(CELL,1)*rho_ij
+        ! kappa.theta.rho = p/exner
+        ! => X = (p/p0)/(exner/Cp)
+        !      = (p/p0)^(1-kappa)
+        pres(CELL) = preff*(X_ij**gamma)  ! pressure
+        ! Compute Exner function (needed by compute_caldyn_fast)
+        ! other formulae possible if exponentiation is slow
+        pk(CELL)   = cpp*((pres(CELL)/preff)**kappa) ! Exner
+      END_BLOCK
+    CASE(thermo_entropy)
+      ON_PRIMAL
+        rho_ij = g*rhodz(CELL)/(geopot(UP(CELL))-geopot(CELL))
+        T_ij = Treff*exp( (theta(CELL,1)+Rd*log(vreff*rho_ij))/Cvd )
+        pres(CELL) = rho_ij*Rd*T_ij
+        pk(CELL)   = T_ij
+      END_BLOCK
+    CASE DEFAULT
+      STOP
+    END SELECT
+  END_BLOCK
+
+  ! We need a barrier here because we compute pres above and do a vertical difference below
+  BARRIER
+
+  FORALL_CELLS_EXT('1', 'llm+1')
+    ON_PRIMAL
+      CST_IFTHEN(IS_BOTTOM_LEVEL)
+        ! Lower BC
+        dW(CELL)   = (1./g)*(pbot-rho_bot*(geopot(CELL)-PHI_BOT(HIDX(CELL)))-pres(CELL)) - m_il(CELL)
+      CST_ELSEIF(IS_TOP_INTERFACE)
+        ! Top BC
+        dW(CELL)   = (1./g)*(pres(DOWN(CELL))-ptop) - m_il(CELL)
+      CST_ELSE
+        dW(CELL)   = (1./g)*(pres(DOWN(CELL))-pres(CELL)) - m_il(CELL)
+      CST_ENDIF
+      W(CELL)    = W(CELL)+tau*dW(CELL) ! update W
+      dPhi(CELL) = g*g*W(CELL)/m_il(CELL)
+    END_BLOCK
+  END_BLOCK
+
+  ! We need a barrier here because we update W above and do a vertical average below
+  BARRIER
+
+  FORALL_CELLS_EXT()
+    ON_PRIMAL
+      ! compute du = -0.5*g^2.grad(w^2)
+      berni(CELL) = (-.25*g*g)*((W(CELL)/m_il(CELL))**2 + (W(UP(CELL))/m_il(UP(CELL)))**2 )
+    END_BLOCK
+  END_BLOCK
+  FORALL_CELLS_EXT()
+    ON_EDGES
+      du(EDGE) = SIGN*(berni(CELL1)-berni(CELL2))
+    END_BLOCK
+  END_BLOCK
+END_BLOCK
+
+KERNEL(caldyn_vert_NH)
+  FORALL_CELLS_EXT()
+    ON_PRIMAL
+      CST_IF(IS_BOTTOM_LEVEL,     w_ij = ( W(CELL)+.5*W(UP(CELL)) )/mass(CELL) )
+      CST_IF(IS_INNER_LAYER , w_ij = .5*( W(CELL)+W(UP(CELL)) )/mass(CELL) )
+      CST_IF(IS_TOP_LAYER,        w_ij = ( .5*W(CELL)+W(UP(CELL)) )/mass(CELL) )
+      wflux_ij = .5*(wflux(CELL)+wflux(UP(CELL)))
+      W_etadot(CELL) = wflux_ij*w_ij
+      eta_dot(CELL) = wflux_ij / mass(CELL)
+      wcov(CELL) = w_ij*(geopot(UP(CELL))-geopot(CELL))
+    END_BLOCK
+  END_BLOCK
+  ! add NH term to du
+  FORALL_CELLS()
+    ON_EDGES
+      du(EDGE) = du(EDGE) - .5*(wcov(CELL1)+wcov(CELL2))*SIGN*(eta_dot(CELL2)-eta_dot(CELL1))
+    END_BLOCK
+  END_BLOCK
+  ! add NH terms to dW, dPhi
+  ! FIXME : TODO top and bottom
+  FORALL_CELLS('2','llm')
+    ON_PRIMAL
+      dPhi(CELL)=dPhi(CELL)-wflux(CELL)*(geopot(UP(CELL))-geopot(DOWN(CELL)))/(mass(DOWN(CELL))+mass(CELL))
+    END_BLOCK
+  END_BLOCK
+
+  ! We need a barrier here because we compute W_etadot above and do a vertical difference below
+  BARRIER
+
+  FORALL_CELLS('1','llm+1')
+    ON_PRIMAL
+      CST_IFTHEN(IS_BOTTOM_LEVEL)
+         dW(CELL) = dW(CELL) - W_etadot(CELL)
+      CST_ELSEIF(IS_TOP_INTERFACE)
+         dW(CELL) = dW(CELL) + W_etadot(DOWN(CELL))
+      CST_ELSE
+         dW(CELL) = dW(CELL) + W_etadot(DOWN(CELL)) - W_etadot(CELL)
+      CST_ENDIF
+    END_BLOCK
+  END_BLOCK
+END_BLOCK
+
+KERNEL(caldyn_slow_NH)
+  FORALL_CELLS_EXT('1', 'llm+1')
+     ON_PRIMAL
+        CST_IF(IS_INNER_INTERFACE, w_il(CELL) = 2.*W(CELL)/(rhodz(KDOWN(CELL))+rhodz(KUP(CELL))) )
+        CST_IF(IS_BOTTOM_LEVEL,    w_il(CELL) = 2.*W(CELL)/rhodz(KUP(CELL)) )
+        CST_IF(IS_TOP_INTERFACE,   w_il(CELL) = 2.*W(CELL)/rhodz(KDOWN(CELL)) )
+     END_BLOCK
+  END_BLOCK
+  FORALL_CELLS_EXT('1', 'llm+1')
+     ON_EDGES
+       ! compute DePhi, v_el, G_el, F_el
+       ! v_el, W2_el and therefore G_el incorporate metric factor le_de
+       ! while DePhil, W_el and F_el do not
+       W_el         = .5*( W(CELL2)+W(CELL1) )
+       DePhil(EDGE) = SIGN*(Phi(CELL2)-Phi(CELL1))
+       F_el(EDGE)   = DePhil(EDGE)*W_el
+       W2_el        = .5*LE_DE * ( W(CELL1)*w_il(CELL1) + W(CELL2)*w_il(CELL2) )
+       v_el(EDGE)   = .5*LE_DE*(u(KUP(EDGE))+u(KDOWN(EDGE))) ! checked
+       G_el(EDGE)   = v_el(EDGE)*W_el - DePhil(EDGE)*W2_el
+     END_BLOCK
+  END_BLOCK
+
+  FORALL_CELLS_EXT('1', 'llm+1')
+    ! compute GradPhi2, dPhi, dW
+    ON_PRIMAL
+       gPhi2=0.
+       dP=0.
+       divG=0
+       FORALL_EDGES       
+         gPhi2 = gPhi2 + LE_DE*DePhil(EDGE)**2
+         dP = dP + LE_DE*DePhil(EDGE)*v_el(EDGE)
+         divG = divG + SIGN*G_el(EDGE) ! -div(G_el), G_el already has le_de
+       END_BLOCK
+       gradPhi2(CELL) = 1./(2.*AI) * gPhi2
+       dPhi(CELL) = gradPhi2(CELL)*w_il(CELL) - 1./(2.*AI)*dP
+       dW(CELL) = (-1./AI)*divG
+    END_BLOCK
+  END_BLOCK
+
+  ! We need a barrier here because we compute gradPhi2, F_el and w_il above and do a vertical average below
+  BARRIER
+
+  FORALL_CELLS_EXT()
+    ! Compute berni at scalar points
+    ON_PRIMAL
+       u2=0.
+       FORALL_EDGES       
+         u2 = u2 + LE_DE*u(EDGE)**2
+       END_BLOCK
+       berni(CELL) = 1./(4.*AI) * u2 - .25*( gradPhi2(CELL)*w_il(CELL)**2 + gradPhi2(UP(CELL))*w_il(UP(CELL))**2 )
+    END_BLOCK
+  END_BLOCK
+
+  FORALL_CELLS_EXT()
+    ON_EDGES
+    ! Compute mass flux and grad(berni)
+      uu = .5*(rhodz(CELL1)+rhodz(CELL2))*u(EDGE) - .5*( F_el(EDGE)+F_el(UP(EDGE)) )
+      hflux(EDGE) = LE_DE*uu
+      du(EDGE) = SIGN*(berni(CELL1)-berni(CELL2))
+    END_BLOCK
+  END_BLOCK
+
+END_BLOCK
Index: /codes/icosagcm/devel/Python/src/kernels_caldyn_hevi.jin
===================================================================
--- /codes/icosagcm/devel/Python/src/kernels_caldyn_hevi.jin	(revision 615)
+++ /codes/icosagcm/devel/Python/src/kernels_caldyn_hevi.jin	(revision 615)
@@ -0,0 +1,202 @@
+KERNEL(caldyn_vert)
+  SEQUENCE
+    BODY('llm-1,1,-1')
+      ! cumulate mass flux convergence from top to bottom
+      convm(CELL) = convm(CELL) + convm(UP(CELL))
+    END_BLOCK
+    EPILOGUE(1)
+      dmass_col(HIDX(CELL)) = convm(CELL)
+    END_BLOCK
+    BODY('2,llm')
+      ! Compute vertical mass flux (l=1,llm+1 set to zero at init)
+      wflux(CELL) = mass_bl(CELL) * dmass_col(HIDX(CELL)) - convm(CELL) 
+    END_BLOCK
+  END_BLOCK
+
+  ! make sure wflux is up to date
+  BARRIER
+  
+  FORALL_CELLS('2','llm')
+    ON_EDGES
+      wwuu(EDGE) = .5*(wflux(CELL1)+wflux(CELL2))*(u(EDGE)-u(DOWN(EDGE)))
+    END_BLOCK
+  END_BLOCK
+
+  ! make sure wwuu is up to date
+  BARRIER
+
+  FORALL_CELLS()
+    ON_PRIMAL
+      convm(CELL) = mass_dbk(CELL) * dmass_col(HIDX(CELL)) ! FIXME : we shoud prognose mass_col
+    END_BLOCK
+    ON_EDGES
+      du(EDGE) = du(EDGE) - (wwuu(EDGE)+wwuu(UP(EDGE))) / (rhodz(CELL1)+rhodz(CELL2))
+    END_BLOCK
+  END_BLOCK
+
+  DO iq=1,nqdyn
+    FORALL_CELLS('2', 'llm')
+      ON_PRIMAL
+        dtheta_rhodz(CELL,iq) = dtheta_rhodz(CELL,iq) + 0.5*(theta(CELL,iq)+theta(DOWN(CELL),iq))*wflux(CELL)  
+      END_BLOCK
+    END_BLOCK
+    FORALL_CELLS('1', 'llm-1')
+      ON_PRIMAL
+        dtheta_rhodz(CELL,iq) = dtheta_rhodz(CELL,iq) - 0.5*(theta(CELL,iq)+theta(UP(CELL),iq))*wflux(UP(CELL))  
+      END_BLOCK
+    END_BLOCK
+  END DO
+END_BLOCK
+
+KERNEL(gradient)
+  FORALL_CELLS_EXT()
+    ON_EDGES
+      grad(EDGE) = SIGN*(b(CELL2)-b(CELL1))
+    END_BLOCK
+  END_BLOCK
+END_BLOCK
+
+KERNEL(div)
+  FORALL_CELLS_EXT()
+    ON_PRIMAL
+      div_ij=0.
+      FORALL_EDGES
+        div_ij = div_ij + SIGN*LE_DE*u(EDGE)
+      END_BLOCK
+      divu(CELL) = div_ij / AI
+    END_BLOCK
+  END_BLOCK
+END_BLOCK
+
+KERNEL(theta)
+  IF(caldyn_eta==eta_mass) THEN ! Compute mass
+    ! FIXME : here mass_col is computed from rhodz
+    ! so that the DOFs are the same whatever caldyn_eta
+    ! in DYNAMICO mass_col is prognosed rather than rhodz
+    SEQUENCE_EXT
+      PROLOGUE(0)
+        mass_col(HIDX(CELL))=0.
+      END_BLOCK
+      BODY('1,llm')
+          mass_col(HIDX(CELL)) = mass_col(HIDX(CELL)) + rhodz(CELL)
+      END_BLOCK
+    END_BLOCK
+    FORALL_CELLS_EXT()
+      ON_PRIMAL
+        ! FIXME : formula below (used in DYNAMICO) is for dak, dbk based on pressure rather than mass
+        !        m = mass_dak(CELL)+(mass_col(HIDX(CELL))*g+ptop)*mass_dbk(CELL)
+      	!        rhodz(CELL) = m/g
+        rhodz(CELL) = mass_dak(CELL) + mass_col(HIDX(CELL))*mass_dbk(CELL)
+      END_BLOCK
+    END_BLOCK
+  END IF
+  DO iq=1,nqdyn
+    FORALL_CELLS_EXT()
+      ON_PRIMAL
+        theta(CELL,iq) = theta_rhodz(CELL,iq)/rhodz(CELL)
+      END_BLOCK
+    END_BLOCK
+  END DO
+END_BLOCK
+
+KERNEL(caldyn_fast)
+!
+  SELECT CASE(caldyn_thermo)
+  CASE(thermo_boussinesq)
+    FORALL_CELLS()
+      ON_PRIMAL
+        berni(CELL) = pk(CELL)
+        ! from now on pk contains the vertically-averaged geopotential
+        pk(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
+      END_BLOCK
+    END_BLOCK
+  CASE(thermo_theta)
+    FORALL_CELLS()
+      ON_PRIMAL
+        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
+      END_BLOCK
+    END_BLOCK
+  CASE(thermo_entropy)
+    FORALL_CELLS()
+      ON_PRIMAL
+        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
+        berni(CELL) = berni(CELL) + pk(CELL)*(cpp-theta(CELL,1)) ! Gibbs = Cp.T-Ts = T(Cp-s)
+      END_BLOCK
+    END_BLOCK
+  CASE DEFAULT
+    PRINT *, 'Unsupported value of caldyn_thermo : ',caldyn_thermo  ! FIXME
+    STOP
+  END SELECT
+!
+  FORALL_CELLS()
+    ON_EDGES
+      due = .5*(theta(CELL1,1)+theta(CELL2,1))*(pk(CELL2)-pk(CELL1)) + berni(CELL2)-berni(CELL1)
+      du(EDGE) = du(EDGE) - SIGN*due
+      u(EDGE) = u(EDGE) + tau*du(EDGE)
+    END_BLOCK
+  END_BLOCK  
+!
+END_BLOCK
+
+KERNEL(pvort_only)
+  FORALL_CELLS_EXT()
+    ON_DUAL
+      etav = 0.d0
+      FORALL_EDGES
+         etav = etav + SIGN*u(EDGE)
+      END_BLOCK
+      hv=0.
+      FORALL_VERTICES
+        hv = hv + RIV2*rhodz(VERTEX)
+      END_BLOCK
+      qv(DUAL_CELL) = (etav + FV*AV )/(hv*AV)
+    END_BLOCK
+  END_BLOCK
+
+  FORALL_CELLS()
+    ON_EDGES
+      qu(EDGE)=0.5d0*(qv(VERTEX1)+qv(VERTEX2))
+    END_BLOCK
+  END_BLOCK
+END_BLOCK
+
+KERNEL(coriolis)
+!
+  DO iq=1,nqdyn
+    FORALL_CELLS_EXT()
+      ON_EDGES
+        Ftheta(EDGE) = .5*(theta(CELL1,iq)+theta(CELL2,iq))*hflux(EDGE)
+      END_BLOCK
+    END_BLOCK
+    FORALL_CELLS()
+      ON_PRIMAL
+        divF=0.
+	FORALL_EDGES
+	  divF = divF + Ftheta(EDGE)*SIGN
+	END_BLOCK
+        dtheta_rhodz(CELL,iq) = -divF / AI
+      END_BLOCK
+    END_BLOCK
+  END DO ! iq
+!
+  FORALL_CELLS()
+    ON_PRIMAL
+      divF=0.
+      FORALL_EDGES
+        divF = divF + hflux(EDGE)*SIGN
+      END_BLOCK
+      convm(CELL) = -divF / AI
+    END_BLOCK
+  END_BLOCK
+!
+  FORALL_CELLS() 
+    ON_EDGES
+      du_trisk=0.
+      FORALL_TRISK
+        du_trisk = du_trisk + WEE*hflux(EDGE_TRISK)*(qu(EDGE)+qu(EDGE_TRISK))
+      END_BLOCK
+      du(EDGE) = du(EDGE) + .5*du_trisk
+    END_BLOCK
+  END_BLOCK
+
+END_BLOCK
Index: /codes/icosagcm/devel/Python/src/kernels_caldyn_hydro.jin
===================================================================
--- /codes/icosagcm/devel/Python/src/kernels_caldyn_hydro.jin	(revision 615)
+++ /codes/icosagcm/devel/Python/src/kernels_caldyn_hydro.jin	(revision 615)
@@ -0,0 +1,112 @@
+#define THECELL {{ thecell }}
+
+{# ---------------- macro to generate code computing pressure top-down ---------------
+   formula = formula to compute 'gravitational' mass
+   = rhodz (dry) rhodz*theta (boussinesq) rhodz*(1+qv) (moist)         #}
+
+#define BALANCE(formula) {% call(thecell) balance() %} formula {% endcall %}
+{% macro balance() %} {% set formula=caller %}
+SEQUENCE_EXT
+  PROLOGUE(llm)
+    pk(CELL) = ptop + .5*g*{{ formula('CELL') }}
+  END_BLOCK
+  BODY('llm-1,1,-1') 
+    pk(CELL) = pk(UP(CELL)) + (.5*g)*({{ formula('CELL') }}+{{ formula('UP(CELL)') }})
+  END_BLOCK
+  IF(caldyn_eta == eta_lag) THEN
+    EPILOGUE(1)
+      ps(HIDX(CELL)) = pk(CELL) + .5*g*{{ formula('CELL') }}
+    END_BLOCK
+  END IF
+END_BLOCK
+{%- endmacro %}
+
+{# ------------ macro to generate code computing geopotential bottom-up --------------
+   var = variable to be stored in pk(CELL)
+   caller() computes gv = g*v where v = specific volume
+   details depend on caldyn_thermo #}
+
+#define GEOPOT(var) {% call geopot(var) %}
+{% macro geopot(var) %} {% set formula=caller %}
+  SEQUENCE_EXT
+    BODY('1,llm')
+      p_ik = pk(CELL)
+      {{ formula() }}
+      pk(CELL) = {{ var }}
+      geopot(UP(CELL)) = geopot(CELL) + gv*rhodz(CELL)
+    END_BLOCK
+  END_BLOCK
+{%- endmacro %}
+
+#define END_GEOPOT {% endcall %}
+
+KERNEL(compute_geopot)
+  SELECT CASE(caldyn_thermo)
+  CASE(thermo_boussinesq)
+    ! use hydrostatic balance with theta*rhodz to find pk (=Lagrange multiplier=pressure)
+    BALANCE( theta(THECELL,1)*rhodz(THECELL) )
+    ! now pk contains the Lagrange multiplier (pressure)
+    ! specific volume 1 = dphi/g/rhodz
+    SEQUENCE_EXT
+      BODY('1,llm')
+        geopot(UP(CELL)) = geopot(CELL) + g*rhodz(CELL)
+      END_BLOCK
+    END_BLOCK
+  CASE(thermo_theta)
+    BALANCE( rhodz(THECELL) )
+    GEOPOT('exner_ik')
+      exner_ik = cpp * (p_ik/preff) ** kappa
+      gv = (g*kappa)*theta(CELL,1)*exner_ik/p_ik
+    END_GEOPOT
+  CASE(thermo_entropy)
+    BALANCE( rhodz(THECELL) )
+    GEOPOT('temp_ik')
+      temp_ik = Treff*exp((theta(CELL,1) + Rd*log(p_ik/preff))/cpp)
+      gv = (g*Rd)*temp_ik/p_ik   ! specific volume v = Rd*T/p
+    END_GEOPOT
+  CASE(thermo_moist)
+    BALANCE( rhodz(THECELL)*(1.+theta(THECELL,2)) )
+    GEOPOT('temp_ik')
+      qv = theta(CELL,2) ! water vaper mixing ratio = mv/md
+      Rmix = Rd+qv*Rv
+      chi = ( theta(CELL,1) + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff)
+      temp_ik = Treff*exp(chi)
+      ! specific volume v = R*T/p
+      ! R = (Rd + qv.Rv)/(1+qv)
+      gv = g*Rmix*temp_ik/(p_ik*(1+qv))
+    END_GEOPOT
+  END SELECT
+END_BLOCK
+
+KERNEL(caldyn_slow_hydro)
+  FORALL_CELLS_EXT()
+    ON_EDGES
+      uu = .5*(rhodz(CELL1)+rhodz(CELL2))*u(EDGE)
+      hflux(EDGE) = uu*LE_DE
+    END_BLOCK
+  END_BLOCK
+
+  FORALL_CELLS()
+    ON_PRIMAL
+      ke=0.d0
+      FORALL_EDGES
+        ke = ke + LE_DE*u(EDGE)**2
+      END_BLOCK
+      BERNI(CELL)=ke*(.25/AI)
+    END_BLOCK
+  END_BLOCK
+  IF(zero) THEN
+    FORALL_CELLS()
+      ON_EDGES
+        du(EDGE) = SIGN*(berni(CELL1)-berni(CELL2)) ! minus gradient
+      END_BLOCK
+    END_BLOCK
+  ELSE
+    FORALL_CELLS()
+      ON_EDGES
+        du(EDGE) = du(EDGE) + SIGN*(berni(CELL1)-berni(CELL2)) ! minus gradient
+      END_BLOCK
+    END_BLOCK
+  END IF
+
+END_BLOCK
Index: /codes/icosagcm/devel/Python/src/kernels_transport.jin
===================================================================
--- /codes/icosagcm/devel/Python/src/kernels_transport.jin	(revision 615)
+++ /codes/icosagcm/devel/Python/src/kernels_transport.jin	(revision 615)
@@ -0,0 +1,226 @@
+KERNEL(advect_vert)
+
+  FORALL_CELLS_EXT('2','llm')
+    ON_PRIMAL
+      dzqw(CELL) = q(CELL)-q(DOWN(CELL))
+    END_BLOCK
+  END_BLOCK
+
+  ! We need a barrier here because we compute dzqw above and do a vertical average below
+  BARRIER  
+
+  FORALL_CELLS_EXT()
+    ON_PRIMAL
+      CST_IFTHEN(IS_BOTTOM_LEVEL || IS_TOP_LAYER)
+        dzq(CELL)=0.
+      CST_ELSE
+        dzq_down=dzqw(CELL)
+        dzq_up=dzqw(UP(CELL)
+        IF(dzq_up*dzq_down>0.) THEN
+           dzqm = .5*(dzqup+dzq_down)
+           dzqmax = max_slope * min( ABS(dzqw_down), ABS(dzq_up) )
+           dzq(CELL) = sign( min(abs(dzqm),dzqmax) , dzqm )  ! sign(a,b)=a*sign(b) 
+        ELSE
+          dzq(CELL)=0.
+        END IF
+      CST_ENDIF
+    END_BLOCK
+  END_BLOCK
+   
+  ! We need a barrier here because we compute dzq above and do a vertical average below
+  BARRIER
+
+  FORALL_CELLS_EXT('1','llm+1')
+    ON_PRIMAL
+      CST_IFTHEN(IS_BOTTOM_LEVEL || IS_TOP_INTERFACE)
+        wq(CELL)=0.
+      CST_ELSE
+        w = fac*wfluxt(CELL)
+	IF(w>0.) THEN  ! upward transport, upwind side is at level l
+	  sigw = w/mass(DOWN(CELL))
+	  qq   = q(DOWN(CELL) + 0.5*(1.-sigw)*dzq(DOWN(CELL)) ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0
+        ELSE           ! downward transport, upwind side is at level l+1
+          sigw = w/mass(CELL)
+          qq   = q(CELL)      -0.5*(1.+sigw)*dzq(ij,l) ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0               
+	ENDIF
+	wq(ij,l) = w*qq																 
+      CST_ENDIF
+    END_BLOCK
+  END_BLOCK
+
+  ! We need a barrier here because we compute wq above and do a vertical difference below
+  BARRIER
+
+  FORALL_CELLS_EXT()
+    ON_PRIMAL
+      newmass = mass(CELL) + fac*(wflux(CELL)-wflux(UP(CELL)))
+      q(CELL) = ( q(CELL)*mass(CELL) + wq(CELL)-wq(UP(CELL)) ) / newmass
+      IF(update_mass) mass(CELL) = newmass
+    END_BLOCK
+  END_BLOCK
+
+END_BLOCK
+
+KERNEL(advect_horiz)
+  ! evaluate tracer field at point cc using piecewise linear reconstruction
+  ! q(cc)= q0 + gradq.(cc-xyz_i)  with xi centroid of hexagon
+  ! sign*hfluxt>0 iff outgoing
+  FORALL_CELLS_EXT()
+    ON_EDGES
+      IF(SIGN*hfluxt(EDGE)>0.) THEN
+        qe = qi(CELL1)
+        qe = qe + (cc(EDGE,1)-xyz_i(HIDX(CELL1),1))*gradq3d(CELL1,1) 
+        qe = qe + (cc(EDGE,2)-xyz_i(HIDX(CELL1),2))*gradq3d(CELL1,2) 
+        qe = qe + (cc(EDGE,3)-xyz_i(HIDX(CELL1),3))*gradq3d(CELL1,3) 
+      ELSE
+        qe = qi(CELL2)
+        qe = qe + (cc(EDGE,1)-xyz_i(HIDX(CELL2),1))*gradq3d(CELL2,1) 
+        qe = qe + (cc(EDGE,2)-xyz_i(HIDX(CELL2),2))*gradq3d(CELL2,2) 
+        qe = qe + (cc(EDGE,3)-xyz_i(HIDX(CELL2),3))*gradq3d(CELL2,3) 
+      END IF
+      qflux(EDGE) = hfluxt(EDGE)*qe
+      IF(diagflux_on) qfluxt(EDGE) = qfluxt(EDGE)+qflux(EDGE)
+    END_BLOCK
+  END_BLOCK
+
+  ! update q and, if update_mass, update mass
+  FORALL_CELLS()
+    ON_PRIMAL
+      dmass=0.
+      dq=0.
+      FORALL_EDGES
+        dmass = dmass + SIGN*hfluxt(EDGE)
+        dq = dq + SIGN*qflux(EDGE)
+      END_BLOCK
+      newmass = mass(CELL) - dmass/AI
+      qi(CELL) = (qi(CELL)*mass(CELL)-dq/AI) / newmass
+      IF(update_mass) mass(CELL)=newmass
+    END_BLOCK
+  END_BLOCK
+
+END_BLOCK
+
+{% macro perot_reconstruction() %}
+  ! Perot reconstruction based on Gauss theorem
+  ! u = sum( u.edge_normal * edge_length * (edge_midpoint-cell_centroid) ) /cell_area
+  FORALL_CELLS()
+    ON_PRIMAL
+      cx=centroid(HIDX(CELL),1)
+      cy=centroid(HIDX(CELL),2)
+      cz=centroid(HIDX(CELL),3)
+      ux=0. ; uy=0. ; uz=0.
+      FORALL_EDGES
+        {{ caller() }}
+	ux = ux + ue_le*(.5*(xyz_v(HIDX(VERTEX1),1)+xyz_v(HIDX(VERTEX2),1))-cx)
+	uy = uy + ue_le*(.5*(xyz_v(HIDX(VERTEX1),2)+xyz_v(HIDX(VERTEX2),2))-cy)
+	uz = uz + ue_le*(.5*(xyz_v(HIDX(VERTEX1),3)+xyz_v(HIDX(VERTEX2),3))-cz)
+      END_BLOCK
+      fac = scale*(1./AI)
+      ucenter(CELL,1)=ux*fac
+      ucenter(CELL,2)=uy*fac
+      ucenter(CELL,3)=uz*fac
+    END_BLOCK
+  END_BLOCK
+{% endmacro %}
+
+KERNEL(wind_centered)
+{% call() perot_reconstruction() %}
+ue_le = SIGN*ue(EDGE)*le(HIDX(EDGE))
+{% endcall %}
+END_BLOCK
+
+KERNEL(flux_centered)
+! NB : here the input data is a flux and has already the factor l_e in it
+! Input data is rescaled by factor scale
+{% call() perot_reconstruction() %}
+ue_le = SIGN*ue(EDGE)
+{% endcall %}
+END_BLOCK
+
+/*------------------------------------- Energetics --------------------------------*/
+
+{% macro energy_flux(energy) %}
+  FORALL_CELLS_EXT()
+    ON_PRIMAL
+      {{ caller() }}
+      {{ energy }}(CELL) = {{ energy }}(CELL) + frac*rhodz(CELL)*energy
+      pk(CELL) = energy
+    END_BLOCK
+  END_BLOCK
+  FORALL_CELLS_EXT()
+    ON_EDGES
+      {{ energy }}_flux(EDGE) = {{ energy }}_flux(EDGE) + .5*massflux(EDGE)*(pk(CELL1)+pk(CELL2))
+    END_BLOCK
+  END_BLOCK
+{% endmacro %}
+
+KERNEL(energy_fluxes)
+  ! First diagnose geopotential and temperature, column-wise
+  BARRIER
+  SEQUENCE_EXT
+    PROLOGUE(llm)
+      pk(CELL) = ptop + .5*g*rhodz(CELL)
+    END_BLOCK
+    BODY('llm-1,1,-1')
+      pk(CELL) = pk(UP(CELL)) + (.5*g)*( rhodz(CELL)+rhodz(UP(CELL)) )
+    END_BLOCK
+  END_BLOCK
+  ! NB : at this point pressure is stored in array pk
+  !      pk then serves as buffer to store temperature
+  SELECT CASE(caldyn_thermo)
+  CASE(thermo_theta)
+    GEOPOT('temp_ik')
+      theta_ik = theta_rhodz(CELL,1)/rhodz(CELL)
+      theta(CELL) = theta_ik
+      temp_ik = theta_ik*(p_ik/preff)**kappa
+      gv = (g*Rd)*temp_ik/p_ik
+    END_GEOPOT
+  CASE(thermo_entropy)
+    GEOPOT('temp_ik')
+      theta_ik = theta_rhodz(CELL,1)/rhodz(CELL)
+      temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp)
+      theta(CELL) = Treff*exp(theta_ik/cpp)
+      gv = (g*Rd)*temp_ik/p_ik   ! specific volume v = Rd*T/p
+    END_GEOPOT
+  END SELECT
+  BARRIER
+  ! Now accumulate energies and energy fluxes
+  ! NB : at this point temperature is stored in array pk
+  !      pk then serves as buffer to store energy
+  ! enthalpy
+  {% call energy_flux('enthalpy') %}
+    energy = cpp*pk(CELL)
+  {% endcall %}
+  ! potential energy
+  {% call energy_flux('epot') %}
+    energy = .5*(geopot(UP(CELL))+geopot(CELL))
+  {% endcall %}
+  ! theta
+  {% call energy_flux('thetat') %}
+    energy = theta(CELL)
+  {% endcall %}
+  ! kinetic energy
+  {% call energy_flux('ekin') %}
+    energy=0.d0
+    FORALL_EDGES
+      energy = energy + le(HIDX(EDGE))*de(HIDX(EDGE))*ue(EDGE)**2
+    END_BLOCK
+    energy = energy * (.25/AI)
+  {% endcall %}
+  ! ulon
+  {% call energy_flux('ulon') %}
+    cx=centroid(HIDX(CELL),1)
+    cy=centroid(HIDX(CELL),2)
+    cz=centroid(HIDX(CELL),3)
+    ux=0. ; uy=0. ; uz=0.
+    FORALL_EDGES
+      ue_le = SIGN*ue(EDGE)*le(HIDX(EDGE))
+      ux = ux + ue_le*(.5*(xyz_v(HIDX(VERTEX1),1)+xyz_v(HIDX(VERTEX2),1))-cx)
+      uy = uy + ue_le*(.5*(xyz_v(HIDX(VERTEX1),2)+xyz_v(HIDX(VERTEX2),2))-cy)
+      uz = uz + ue_le*(.5*(xyz_v(HIDX(VERTEX1),3)+xyz_v(HIDX(VERTEX2),3))-cz)
+    END_BLOCK
+    ulon_i = ux*elon_i(HIDX(CELL),1) + uy*elon_i(HIDX(CELL),2) + uz*elon_i(HIDX(CELL),3)
+    energy = ulon_i*(1./AI)
+  {% endcall %}
+
+END_BLOCK
Index: /codes/icosagcm/devel/Python/src/macros.jin
===================================================================
--- /codes/icosagcm/devel/Python/src/macros.jin	(revision 615)
+++ /codes/icosagcm/devel/Python/src/macros.jin	(revision 615)
@@ -0,0 +1,39 @@
+{#  ---------------------------------- Macros common to all variants ---------------------------------------- #}
+
+{% set flat='#' %}
+
+{% macro cst_if(condition, action) %}
+{{ flat }}if {{ condition }}
+{{ action }}
+{{ flat }}endif
+{%- endmacro %}
+ 
+{% macro define(key,val) %} 
+{{ flat }}undef {{ key }} 
+{{ flat }}define {{ key }} {{ val }} 
+{%- endmacro %}
+
+{% macro define_cppmacro(key,args,val) %} 
+{{ flat }}undef {{ key }} 
+{{ flat }}define {{ key+args }} {{ val }} 
+{%- endmacro %}
+
+{% macro undef(key) %} 
+{{ flat }}undef {{ key }}
+{{ flat }}define {{ key }} {{ key }}_UNDEFINED_HERE
+{%- endmacro %}
+
+{% macro dummy() -%}
+  {{ caller() }}
+{%- endmacro %}
+
+{% macro define_kernel(name) %} 
+{{ flat }}ifdef KERNEL_{{name}}
+!--------------------------------------------------------------------------
+!---------------------------- {{ name }} ----------------------------------
+  {{ caller() }}
+!---------------------------- {{ name }} ----------------------------------
+!--------------------------------------------------------------------------
+
+{{ flat }}endif
+{%- endmacro %}
Index: /codes/icosagcm/devel/Python/src/setup.py
===================================================================
--- /codes/icosagcm/devel/Python/src/setup.py	(revision 615)
+++ /codes/icosagcm/devel/Python/src/setup.py	(revision 615)
@@ -0,0 +1,4 @@
+from distutils.core import setup
+from Cython.Build import cythonize
+setup(  name = "DYNAMICO kernels for unstructured mesh",
+        ext_modules = cythonize(["*.pyx"]) )
Index: /codes/icosagcm/devel/Python/src/unstructured.pyx
===================================================================
--- /codes/icosagcm/devel/Python/src/unstructured.pyx	(revision 615)
+++ /codes/icosagcm/devel/Python/src/unstructured.pyx	(revision 615)
@@ -0,0 +1,544 @@
+import time
+import math
+import numpy as np
+import netCDF4 as cdf
+import matplotlib.tri as tri
+import matplotlib.pyplot as plt
+
+import dynamico.wrap as wrap
+from ctypes import c_void_p, c_int, c_double, c_bool
+radian=180/math.pi
+
+#------------- direct Cython interface to DYNAMICO routines -------------#
+
+cdef extern from "functions.h":
+     cpdef void init_params()
+     cpdef void setup_xios()
+     cpdef void call_xios_set_timestep(double)
+     cpdef void call_xios_update_calendar(int)
+
+#------------- import and wrap DYNAMICO routines -------------#
+
+ker=wrap.Struct() # store imported fun X as funs.X
+
+check_args = False # use True instead of False for debugging, probably with some overhead
+
+try:    
+    kernels = wrap.SharedLib(vars(ker), 'libkernels.so', check_args=check_args) 
+    setvar, setvars, getvar, getvars = kernels.setvar, kernels.setvars, kernels.getvar, kernels.getvars
+except OSError:
+    print """
+Unable to load shared library 'libkernels.so' !
+    """
+    raise
+
+# providing a full prototype enables type-checking when calling
+# if a number n is present in the prototype, the previous type is repeated n times
+kernels.import_funs([
+#                     ['setup_xios',None],
+#                     ['call_xios_set_timestep',c_double],
+#                     ['call_xios_update_calendar',c_int]
+#                     ['init_params',c_double],
+                     ['init_mesh',c_void_p,13],
+                     ['init_metric', c_void_p,6],
+                     ['init_hybrid', c_void_p,3],
+                     ['caldyn_hevi', c_double, c_void_p,20],
+                     ['partition_graph', c_int,2, c_void_p,3, c_int, c_void_p],
+                     ])
+
+# set/get global variables
+eta_mass,eta_lag=(1,2)
+thermo_theta,thermo_entropy,thermo_moist,thermo_boussinesq=(1,2,3,4)
+
+kernels.addvars(
+    c_bool,'hydrostatic','debug_hevi_solver','rigid',
+    c_int,'llm','nqdyn','primal_num','max_primal_deg',
+    'dual_num','max_dual_deg','edge_num','max_trisk_deg',
+    'caldyn_thermo','caldyn_eta','nb_threads',
+    c_double,'elapsed','g', 'ptop', 'cpp', 'cppv',
+    'Rd', 'Rv', 'preff', 'Treff', 'pbot', 'rho_bot', 'Phi_bot')
+
+elapsed=0.
+
+#----------------------------- Base class for dynamics ------------------------
+
+class Caldyn:
+    def __init__(self,mesh):
+        self.mesh=mesh
+        fps, ftheta, fmass = mesh.field_ps, mesh.field_theta, mesh.field_mass
+        fw, fu, fz = mesh.field_w, mesh.field_u, mesh.field_z
+        self.ps, self.ms, self.dms       = fps(), fps(), fps()
+        self.s, self.hs, self.dhs        = ftheta(), ftheta(), ftheta()
+        self.pk, self.berni, self.geopot, self.hflux = fmass(),fmass(),fw(),fu()
+        self.qu, self.qv                 = fu(),fz()
+        self.fmass, self.ftheta, self.fu, self.fw = fmass, ftheta, fu, fw
+    def bwd_fast_slow(self, flow, tau):
+        global elapsed
+        time1=time.time()
+        flow,fast,slow = self._bwd_fast_slow_(flow,tau)
+        time2=time.time()
+        elapsed=elapsed+time2-time1
+        return flow,fast,slow
+
+# when calling caldyn_hevi, arrays for tendencies must be re-created each time
+# to avoid overwriting in the same memory space when time scheme is multi-stage
+
+#-------------------------- Shallow-water dynamics ---------------------
+
+class Caldyn_RSW(Caldyn):
+    def __init__(self,mesh):
+        Caldyn.__init__(self,mesh)
+        setvars(('hydrostatic','caldyn_thermo','caldyn_eta'),
+                (True,thermo_boussinesq,eta_lag))
+        self.dhs = self.fmass()
+        ker.init_params()
+    def _bwd_fast_slow_(self, flow, tau):
+        h,u = flow
+        # h*s = h => uniform buoyancy s=1 => shallow-water
+        dh, du_slow, du_fast, hs, buf = self.fmass(), self.fu(), self.fu(), h.copy(), self.geopot
+        ker.caldyn_hevi(tau, self.ms, h, hs, u, self.geopot, buf,
+                  self.s, self.ps, self.pk, self.hflux, self.qv,
+                  self.dms, dh, self.dhs, du_fast, du_slow,
+                  buf, buf, buf, buf)
+        return (h,u), (0.,du_fast), (dh,du_slow)
+
+#----------------------------------- HPE ------------------------------------
+
+class Caldyn_HPE(Caldyn):
+    def __init__(self,caldyn_thermo,caldyn_eta, mesh,metric,thermo,BC,g):
+        Caldyn.__init__(self,mesh)
+        setvars(('hydrostatic','caldyn_thermo','caldyn_eta',
+                 'g','ptop','Rd','cpp','preff','Treff'),
+                (True,caldyn_thermo,caldyn_eta,
+                 g,BC.ptop,thermo.Rd,thermo.Cpd,thermo.p0,thermo.T0))
+        ker.init_params()
+    def _bwd_fast_slow_(self, flow, tau):
+        dm, dS, du_slow, du_fast, buf = self.fmass(), self.ftheta(), self.fu(), self.fu(), self.geopot
+        m,S,u = flow
+        ker.caldyn_hevi(tau, self.ms, m, S, u, self.geopot, buf,
+                  self.s, self.ps, self.pk, self.hflux, self.qv,
+                  self.dms, dm, dS, du_fast, du_slow,
+                  buf, buf, buf, buf)
+        return (m,S,u), (0.,0.,du_fast), (dm,dS,du_slow)
+
+#----------------------------------- NH ------------------------------------
+
+class Caldyn_NH(Caldyn):
+    def __init__(self,caldyn_thermo,caldyn_eta, mesh,metric,thermo,BC,g):
+        Caldyn.__init__(self,mesh)
+        setvars(('hydrostatic','caldyn_thermo','caldyn_eta',
+                 'g','ptop','Rd','cpp','preff','Treff',
+                 'pbot','rho_bot'),
+                (False,caldyn_thermo,caldyn_eta,
+                 g,BC.ptop,thermo.Rd,thermo.Cpd,thermo.p0,thermo.T0,
+                 BC.pbot.max(), BC.rho_bot.max()))
+        ker.init_params()
+    def bwd_fast_slow(self, flow, tau):
+        ftheta, fmass, fu, fw = self.ftheta, self.fmass, self.fu, self.fw 
+        dm, dS, du_slow, du_fast = fmass(), ftheta(), fu(), fu()
+        dPhi_slow, dPhi_fast, dW_slow, dW_fast = fw(), fw(), fw(), fw()
+        m,S,u,Phi,W = flow
+        ker.caldyn_hevi(tau, self.ms, m, S, u, Phi, W,
+                  self.s, self.ps, self.pk, self.hflux, self.qv,
+                  self.dms, dm, dS, du_fast, du_slow,
+                  dPhi_fast, dPhi_slow, dW_fast, dW_slow)
+        return ((m,S,u,Phi,W), (0.,0.,du_fast,dPhi_fast,dW_fast), 
+                (dm,dS,du_slow,dPhi_slow,dW_slow))
+
+#-------------------------------- Hybrid mass-based coordinate -------------
+
+# compute hybrid coefs from distribution of mass                                                                                                                
+def compute_hybrid_coefs(mass):
+    nx,llm=mass.shape
+    mass_dak = np.zeros((nx,llm))
+    mass_dbk = np.zeros((nx,llm))
+    mass_bl = np.zeros((nx,llm+1))+1.
+    for i in range(nx):
+        m_i = mass[i,:]
+        dbk_i = m_i/sum(m_i)
+        mass_dbk[i,:] = dbk_i
+        mass_bl[i,1:]= 1.-dbk_i.cumsum()
+    return mass_bl, mass_dak, mass_dbk
+
+#----------------------- Cartesian mesh -----------------------
+
+def squeeze(dims):
+#    return np.zeros(dims)
+    return np.zeros([n for n in dims if n>1])
+
+# arrays is a list of arrays
+# vals is a list of tuples
+# each tuple is stored in each array
+def put(ij, deg, arrays, vals):
+    k=0
+    for vv in vals: # vv is a tuple of values to be stored in arrays
+        for array,v in zip(arrays,vv):
+            array[ij,k]=v
+        k=k+1        
+    deg[ij]=k
+
+class Cartesian_mesh:
+    def __init__(self,nx,ny,llm,nqdyn,Lx,Ly,f):
+        dx,dy = Lx/nx, Ly/ny
+        self.dx, self.dy, self.f = dx,dy,f
+        self.nx, self.ny, self.llm, self.nqdyn = nx,ny,llm,nqdyn
+        self.field_z = self.field_mass
+        # 1D coordinate arrays
+        x=(np.arange(nx)-nx/2.)*Lx/nx
+        y=(np.arange(ny)-ny/2.)*Ly/ny
+        lev=np.arange(llm)
+        levp1=np.arange(llm+1)
+        self.x, self.y, self.lev, self.levp1 = x,y,lev,levp1
+        # 3D coordinate arrays
+        self.xx,self.yy,self.ll = np.meshgrid(x,y,lev, indexing='ij')       
+        self.xxp1,self.yyp1,self.llp1 = np.meshgrid(x,y,levp1, indexing='ij')
+        # beware conventions for indexing
+        # Fortran order : llm,nx*ny,nqdyn  / indices start at 1
+        # Python order : nqdyn,ny,nx,llm   / indices start at 0
+        # indices below follow Fortran while x,y follow Python/C
+        index=lambda x,y : ((x+(nx*(y+2*ny)))%(nx*ny))+1
+        indexu=lambda x,y : 2*index(x,y)-1
+        indexv=lambda x,y : 2*index(x,y)
+        indices = lambda shape : np.zeros(shape,np.int32)
+
+        primal_nb = indices(nx*ny)
+        primal_edge = indices((nx*ny,4))
+        primal_ne = indices((nx*ny,4))
+    
+        dual_nb = indices(nx*ny)
+        dual_edge = indices((nx*ny,4))
+        dual_ne = indices((nx*ny,4))
+        dual_vertex = indices((nx*ny,4))
+        Riv2 = np.zeros((nx*ny,4))
+    
+        left = indices(2*nx*ny)
+        right = indices(2*nx*ny)
+        up = indices(2*nx*ny)
+        down = indices(2*nx*ny)
+        le_de = np.zeros(2*nx*ny)
+    
+        trisk_deg = indices(2*nx*ny)
+        trisk = indices((2*nx*ny,4)) # 4 TRiSK coefs per edge
+        wee = np.zeros((2*nx*ny,4))
+    
+        for x in range(nx):
+            for y in range(ny):
+                # NB : Fortran indices start at 1
+                # primal cells
+                put(index(x,y)-1,primal_nb,(primal_edge,primal_ne),
+                    ((indexu(x,y),1),
+                     (indexv(x,y),1),
+                     (indexu(x-1,y),-1),
+                     (indexv(x,y-1),-1) ))
+                # dual cells
+                put(index(x,y)-1,dual_nb,(dual_edge,dual_vertex,dual_ne,Riv2),
+                   ((indexv(x+1,y),index(x,y),1,.25),
+                    (indexu(x,y+1),index(x+1,y),-1,.25),
+                    (indexv(x,y),index(x+1,y+1),-1,.25),
+                    (indexu(x,y),index(x,y+1),1,.25) ))
+                # edges :
+                # left and right are adjacent primal cells
+                # flux is positive when going from left to right
+                # up and down are adjacent dual cells (vertices)
+                # circulation is positive when going from down to up
+                # u-points
+                ij =indexu(x,y)-1 
+                left[ij]=index(x,y)
+                right[ij]=index(x+1,y)
+                down[ij]=index(x,y-1)
+                up[ij]=index(x,y)
+                le_de[ij]=dy/dx
+                put(ij,trisk_deg,(trisk,wee),(
+                    (indexv(x,y),-.25),
+                    (indexv(x+1,y),-.25),
+                    (indexv(x,y-1),-.25),
+                    (indexv(x+1,y-1),-.25)))
+                # v-points
+                ij = indexv(x,y)-1
+                left[ij]=index(x,y)
+                right[ij]=index(x,y+1)
+                down[ij]=index(x,y)
+                up[ij]=index(x-1,y)
+                le_de[ij]=dx/dy
+                put(ij,trisk_deg,(trisk,wee),(
+                    (indexu(x,y),.25),
+                    (indexu(x-1,y),.25),
+                    (indexu(x,y+1),.25),
+                    (indexu(x-1,y+1),.25)))
+        setvars(('llm','nqdyn','edge_num','primal_num','dual_num',
+                 'max_trisk_deg','max_primal_deg','max_dual_deg'),
+                (llm,nqdyn,2*nx*ny,nx*ny,nx*ny,4,4,4) )        
+        ker.init_mesh(primal_nb,primal_edge,primal_ne,
+                  dual_nb,dual_edge,dual_ne,dual_vertex,
+                  left,right,down,up,trisk_deg,trisk)
+
+        Aiv=np.zeros(nx*ny)+dx*dy # Ai=Av=dx*dy
+        ker.init_metric(Aiv,Aiv,f+0.*Aiv,le_de,Riv2,-wee)
+    def field_theta(self): return squeeze((self.nqdyn,self.ny,self.nx,self.llm))
+    def field_mass(self): return squeeze((self.ny,self.nx,self.llm))
+    def field_z(self): return squeeze((self.ny,self.nx,self.llm))
+    def field_w(self): return squeeze((self.ny,self.nx,self.llm+1))
+    def field_u(self): return np.zeros((self.ny,2*self.nx,self.llm))
+    def field_ps(self): return squeeze((self.ny,self.nx))
+    def ucomp(self,u): return u[:,range(0,2*self.nx,2),:]
+    def set_ucomp(self,uv,u): uv[:,range(0,2*self.nx,2),:]=u
+    def vcomp(self,u): return u[:,range(1,2*self.nx,2),:]
+
+#---------------------- MPAS fully unstructured mesh ------------------------
+
+def compute_ne(num,deg,edges,left,right):
+    ne = np.zeros_like(edges)
+    for cell in range(num):
+        for iedge in range(deg[cell]):
+            edge = edges[cell,iedge]-1
+            if left[edge]==cell+1: ne[cell,iedge]+=1
+            if right[edge]==cell+1: ne[cell,iedge]-=1
+            if ne[cell,iedge]==0 : print 'error at cell,iedge', cell, iedge
+    return ne
+
+def plot(tri,data):
+    plt.figure(figsize=(12,4))
+    plt.gca().set_aspect('equal')
+    plt.tricontourf(tri, data, 20)
+    plt.colorbar()
+    plt.ylim((-90,90))
+    plt.xlim((0,360))
+
+class MPAS_Mesh:
+    def __init__(self, gridfile, llm, nqdyn, radius, f):
+        self.gridfile, self.llm, self.nqdyn = gridfile,llm,nqdyn
+        self.radius, self.f = radius, f
+        # open mesh file, get main dimensions
+        try:
+            nc = cdf.Dataset(gridfile, "r")
+        except RuntimeError:
+            print """
+Unable to open grid file %s, maybe you forgot to download it ?
+To do so, go to the 'Python/' dir and execute './get_MPAS_grids.sh'.
+            """ % gridfile
+            raise
+
+        def getdims(*names): return [len(nc.dimensions[name]) for name in names]
+        def getvars(*names): 
+            for name in names : print "getvar %s ..."%name
+            time1=time.time()
+            ret=[nc.variables[name][:] for name in names]
+            print "... Done in %f seconds"%(time.time()-time1)
+            return ret
+        primal_num, edge_num, dual_num = getdims('nCells','nEdges','nVertices')                                   
+        print 'Number of primal cells, dual cells and edges : %d, %d, %d '%(primal_num,dual_num,edge_num)
+        primal_deg, trisk_deg = getvars('nEdgesOnCell','nEdgesOnEdge')
+        dual_deg = [3 for i in range(dual_num) ]
+        dual_deg = np.ascontiguousarray(dual_deg,dtype=np.int32) # NB : Fortran code expects 32-bit ints
+
+        # get indices for stencils 
+        # primal -> vertices (unused)
+        primal_vertex, dual_vertex = getvars('verticesOnCell','cellsOnVertex')
+        # primal <-> edges
+        primal_edge, left_right = getvars('edgesOnCell','cellsOnEdge')
+        left,right = left_right[:,0], left_right[:,1]
+        primal_ne = compute_ne(primal_num,primal_deg,primal_edge,left,right)
+        # dual <-> edges
+        dual_edge, down_up = getvars('edgesOnVertex','verticesOnEdge')
+        down,up = down_up[:,0], down_up[:,1]
+        dual_ne = compute_ne(dual_num,dual_deg,dual_edge,up,down)
+        # primal <-> dual, edges <-> edges (TRiSK)
+        dual_vertex, trisk = getvars('cellsOnVertex','edgesOnEdge')
+        # get positions, lengths, surfaces and weights
+        le,de,Ai,Av = getvars('dvEdge','dcEdge','areaCell','areaTriangle')
+        lat_i,lon_i = getvars('latCell','lonCell')
+        lat_v,lon_v = getvars('latVertex','lonVertex')
+        lat_e,lon_e,angle_e = getvars('latEdge','lonEdge','angleEdge')
+        wee,Riv2 = getvars('weightsOnEdge','kiteAreasOnVertex')
+
+        # fix normalization of wee and Riv2 weights
+        for edge1 in range(edge_num):
+            for i in range(trisk_deg[edge1]):
+                edge2=trisk[edge1,i]-1  # NB Fortran vs C/Python indexing
+                wee[edge1,i] = de[edge1]*wee[edge1,i]/le[edge2]
+        for ivertex in range(dual_num):
+            for j in range(3):
+                Riv2[ivertex,j]=Riv2[ivertex,j]/Av[ivertex]
+            r=Riv2[ivertex,:]
+            r=sum(r)
+            if abs(r-1.)>1e-6 : print 'error with Riv2 at vertex ', ivertex
+            
+        max_primal_deg, max_dual_deg, max_trisk_deg = getdims('maxEdges','vertexDegree','maxEdges2')
+
+        # CRITICAL : force arrays left, etc. to be contiguous in memory:
+        left,right,up,down = [np.ascontiguousarray(x) for x in (left,right,up,down)]
+        trisk,wee,primal_edge = [np.ascontiguousarray(x) for x in (trisk,wee,primal_edge)]
+        
+        print ('Max stencil sizes (div,curl,trisk) : %d, %d, %d' 
+               % (max_primal_deg, max_dual_deg, max_trisk_deg) )
+
+        r2 = radius**2
+        Av = r2*Av
+        fv = f(lon_v,lat_v) # Coriolis parameter
+        self.Ai, self.Av, self.fv = r2*Ai,Av,fv
+        self.le, self.de, self.le_de = radius*le, radius*de, le/de
+        self.trisk_deg, self.trisk, self.wee = trisk_deg, trisk, wee
+        setvars(('llm','nqdyn','edge_num','primal_num','dual_num',
+                      'max_trisk_deg','max_primal_deg','max_dual_deg'),
+                     (llm, nqdyn, edge_num, primal_num,dual_num,
+                      max_trisk_deg, max_primal_deg, max_dual_deg) )
+
+        ker.init_mesh(primal_deg,primal_edge,primal_ne,
+                       dual_deg,dual_edge,dual_ne,dual_vertex,
+                       left,right,down,up,trisk_deg,trisk)
+        ker.init_metric(self.Ai,self.Av,self.fv,le/de,Riv2,wee)
+        
+        for edge in range(edge_num):
+            iedge = trisk[edge,0:trisk_deg[edge]]
+            if iedge.min()<1 :
+                print 'error'
+                          
+        self.primal_num, self.edge_num, self.dual_num = primal_num, edge_num, dual_num
+        def period(x) : return (x+2*math.pi)%(2*math.pi)
+        lon_i, lon_v, lon_e = map(period, (lon_i,lon_v,lon_e))
+        self.lon_i, self.lat_i = lon_i, lat_i
+        self.lon_v, self.lat_v = lon_v, lat_v
+        self.lon_e, self.lat_e, self.angle_e = lon_e, lat_e, angle_e
+        self.primal_deg, self.primal_vertex = primal_deg, primal_vertex
+        self.primal = tri.Triangulation(lon_i*180./math.pi, lat_i*180./math.pi)
+        self.dual_deg, self.dual_vertex = dual_deg, dual_vertex
+        self.dual = tri.Triangulation(lon_v*180./math.pi, lat_v*180./math.pi)
+        self.triedge = tri.Triangulation(lon_e*180./math.pi, lat_e*180./math.pi)
+        self.dx = de.min()
+        self.lon3D_i, self.ll3D = np.meshgrid(lon_i, range(llm))
+        self.lat3D_i, self.ll3D = np.meshgrid(lat_i, range(llm))
+    def field_theta(self): return squeeze((self.nqdyn,self.primal_num,self.llm))
+    def field_mass(self): return squeeze((self.primal_num,self.llm))
+    def field_z(self): return squeeze((self.dual_num,self.llm))
+    def field_w(self): return squeeze((self.primal_num,self.llm+1))
+    def field_u(self): return squeeze((self.edge_num,self.llm))
+    def field_ps(self): return squeeze((self.primal_num,))
+    def ucov2D(self, ulon, ulat): 
+        return self.de*(ulon*np.cos(self.angle_e)+ulat*np.sin(self.angle_e))
+    def ucov3D(self, ulon, ulat):
+        ucov = np.zeros((self.edge_num,ulon.shape[1]))
+        for edge in range(self.edge_num):
+            angle=self.angle_e[edge]
+            ucov[edge,:] = self.de[edge]*(ulon[edge,:]*math.cos(angle)+ulat[edge,:]*math.sin(angle))
+        return ucov
+    def plot_i(self,data):
+        plot(self.primal,data)
+    def plot_v(self,data):
+        plot(self.dual,data)
+    def plot_e(self,data):
+        plot(self.triedge,data)
+
+#-------------------------------------- Mesh partitioning ------------------------------------------#
+
+# Helper functions and interface to ParMETIS
+# list_stencil converts an adjacency graph from array format index[num_cells, MAX_EDGES] to compressed format
+# loc_stencil returns the start/end indices (vtxdist) expected by ParMETIS
+# i.e. index[start:end] with start=vtxdist[cell], end=vtxdist[cell+1] lists the edges of cell 'cell'
+
+def list_stencil(degree, stencil, cond=lambda x:True):
+    for i in range(degree.size):
+        for j in range(degree[i]):
+            s=stencil[i,j]
+            if cond(s): yield stencil[i,j]
+                
+def loc_stencil(degree, stencil):
+    loc=0
+    for i in range(degree.size):
+        yield loc
+        loc=loc+degree[i]
+    yield loc
+
+def partition_mesh(degree, stencil, nparts): 
+    # arguments : PArray1D and PArray2D describing mesh, number of desired partitions
+    dim_cell, degree, stencil = degree.dim, degree.data, stencil.data
+    comm, vtxdist, idx_start, idx_end = dim_cell.comm, dim_cell.vtxdist, dim_cell.start, dim_cell.end
+    mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size()
+    adjncy_loc, xadj_loc = list_stencil(degree, stencil), loc_stencil(degree, stencil)
+    adjncy_loc, xadj_loc = [np.asarray(list(x), dtype=np.int32) for x in (adjncy_loc, xadj_loc)]
+    owner = np.zeros(idx_end-idx_start, dtype=np.int32);
+    ker.partition_graph(mpi_rank, mpi_size, vtxdist, xadj_loc, adjncy_loc, nparts, owner)
+    return owner
+
+def partition_from_stencil(owner2, degree, stencil):
+    # given a stencil dim1->dim2 and owner2 on dim2, define owner[i] on dim1 as min(stencil[i,:] if i is even, max if odd
+    dim1, dim2= degree.dim, owner2.dim
+    degree, stencil, n = degree.data, stencil.data, dim1.end-dim1.start
+    cells2 = list_stencil(degree, stencil)
+    cells2 = sorted(list(set(list(cells2)))) # list of cells for which we need to know owner2
+    get2 = dim2.getter(cells2)
+    owner1, owner2, glob2loc = np.zeros(n, dtype=np.int32), get2(owner2), get2.dict
+    for i in range(n):
+        owners = [ owner2[glob2loc[stencil[i,j]]] for j in range(degree[i]) ]
+        if i%2 == 0 : owner1[i] = min(owners)
+        else : owner1[i] = max(owners)
+    return owner1
+            
+def find_my_cells(owner): # a PArray1D containing the data returned by partition_mesh
+    dim, comm, owner = owner.dim, owner.dim.comm, owner.data
+    mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size()
+    cells=[set() for i in range(mpi_size)]
+    for i in range(len(owner)):
+        cells[owner[i]].add(dim.start+i)
+    cells = [sorted(list(cells[i])) for i in range(mpi_size)]
+    mycells = comm.alltoall(cells)
+    mycells = sorted(sum(mycells, [])) # concatenate into a single list
+    return mycells
+
+#---------------------------------- Stencil management ---------------------------------------#
+
+# Class Stencil represents an adjacency relationship (e.g. cell=>edges)
+# using adjacency information read from PArrays
+# It computes a list of "edges" adjacent to a given list of "cells"
+# This is used to form the sets E0 -> C0 -> E1 -> V1 -> E2 -> C1
+# which are then used to form lists of global indices for V1,C1,E2 such that
+# C0, E0, E1 form contiguous subsets of C1, E2 starting from 0
+
+def reindex(vertex_dict, degree, bounds):
+    for i in range(degree.size):
+        for j in range(degree[i]):
+            bounds[i,j] = vertex_dict[bounds[i,j]]    
+
+class Stencil_glob:
+    def __init__(self, degree, neigh, weight=None):
+        self.degree, self.neigh, self.weight = degree, neigh, weight
+    def __call__(self, cells, neigh_dict=None):
+        return Stencil(cells, self.degree, self.neigh, neigh_dict, self.weight)
+            
+class Stencil:
+    def __init__(self, cells, degree, neigh, neigh_dict, weight=None):
+        get = degree.dim.getter(cells)
+        mydegree, myneigh = [get(x) for x in (degree, neigh)]
+        if not weight is None : myweight = get(weight)
+        if neigh_dict is None : 
+            keep = lambda n : True
+        else : # if neigh_dict is present, only neighbours in dict are retained
+            keep = lambda n : n in neigh_dict
+        neigh_set = list_stencil(mydegree, myneigh, keep)
+        self.neigh_set = list(set(list(neigh_set)     ))                         
+        rej=0
+        for i in range(len(mydegree)): # keep only elements in neigh_dict, in-place
+            k=0
+            for j in range(mydegree[i]):
+                n=myneigh[i,j]
+                if keep(n):
+                    myneigh[i,k]=n
+                    if not weight is None : myweight[i,k]=myweight[i,j]
+                    k=k+1
+                else:
+                    rej=rej+1
+            mydegree[i]=k
+        if neigh_dict is None:
+            neigh_dict = {j:i for i,j in enumerate(self.neigh_set)}
+        myneigh_loc = reindex(neigh_dict, mydegree, myneigh)
+        self.degree, self.neigh_glob, self.neigh_loc = mydegree, myneigh, myneigh_loc
+
+def progressive_iter(mylist, cell_lists):
+    for thelist in cell_lists: 
+        mylist = mylist + list(set(thelist)-set(mylist))
+        yield mylist
+def progressive_list(*cell_lists):
+    # cell_lists : a tuple of lists of global indices, with each list a subset of the next
+    # returns : a list 'mylist' such that for each list 'thelist' in cell_lists, thelist = mylist[0:len(thelist)]
+    # example : edge_list = progressive_list(E0,E1,E2) with E0,E1,E2 increasing lists of cell edges
+    return list(progressive_iter([], cell_lists))
+
Index: /codes/icosagcm/devel/Python/src/unstructured/macros.jin
===================================================================
--- /codes/icosagcm/devel/Python/src/unstructured/macros.jin	(revision 615)
+++ /codes/icosagcm/devel/Python/src/unstructured/macros.jin	(revision 615)
@@ -0,0 +1,250 @@
+#ifdef PASS_PRE1
+#define _AND_ &&
+#define _OR_ ||
+#define _NOT_ !
+#define IS_INNER_LAYER (_NOT_ (IS_BOTTOM_LEVEL _OR_ IS_TOP_LAYER))
+#define IS_INNER_INTERFACE (_NOT_ (IS_BOTTOM_LEVEL _OR_ IS_TOP_INTERFACE))
+#endif
+
+#ifdef PASS_PRE2
+
+#define BARRIER !$OMP BARRIER
+#define IS_INNER_INTERFACE (.NOT. (IS_BOTTOM_LEVEL .OR. IS_TOP_INTERFACE))
+#define KERNEL(name) {% call() define_kernel(#name) %}
+#define FORALL_CELLS(...) {% call(lrange,flags) forall_cells(__VA_ARGS__) %}
+#define FORALL_CELLS_EXT(...) {% call(lrange,flags) forall_cells(__VA_ARGS__) %}
+#define ON_PRIMAL {% call(mesh,lrange) on_mesh('primal',lrange,flags) %}
+#define ON_DUAL {% call(mesh,lrange) on_mesh('dual',lrange,flags) %}
+#define ON_EDGES {% call(mesh,lrange) on_edges('edge',lrange,flags) %}
+#define FORALL_EDGES {% call forall_edges(mesh) %}
+#define FORALL_VERTICES {% call forall_vertices(mesh) %}
+#define FORALL_TRISK {% call forall_trisk() %}
+#define END_BLOCK {% endcall %}
+#define SEQUENCE {% call sequence() %} 
+#define SEQUENCE_EXT {% call sequence() %} 
+#define PROLOGUE(level) {% call at_level(level) %}
+#define EPILOGUE(level) {% call at_level(level) %}
+#define BODY(range) {% call(cell_up,cell_down) body(range) %}
+
+#define CST_IF(condition, action) {{ cst_if(#condition, #action) }}
+#define CST_IFTHEN(condition) {{flat}}if condition
+#define CST_ELSEIF(condition) {{flat}}elif condition
+#define CST_ELSE {{flat}}else
+#define CST_ENDIF {{flat}}endif
+
+#endif
+
+#ifdef PASS_JINJA
+
+{#  -------------------------- Design notes ------------------------------
+
+Using its optional argument 'options', forall_cells() sets variable 'flags' 
+and passes it on to on_mesh() / on_edges()
+
+Similarly variable 'lrange' is passed from forall_cells to on_mesh() / on_egdes()  
+
+-------------------------------------------------------------------------- #}
+
+{# ----------------------- SEQUENCE, FORALL ----------------------- #}
+
+{% macro sequence() -%}
+  
+!$OMP DO SCHEDULE(STATIC)                                                                            
+  DO ij=1,primal_num
+    {{ caller() }}
+  END DO
+!$OMP END DO
+
+{%- endmacro %}
+
+{% macro body(range) -%}
+  {{ define('CELL','l,ij') }}
+    DO l = {{ range }}
+      {{ caller() }}
+    END DO
+  {{ undef('CELL') }}
+{%- endmacro %}
+
+{% macro at_level(lev) -%}
+  {{ define('CELL','l,ij') }}
+  l={{ lev }}
+  {{ caller() }}
+  {{ undef('CELL') }}
+{%- endmacro %}
+
+{% macro forall_cells(start='1', end='llm', options=[]) -%}
+  {% set flags=[options] if options is string else options %}
+  {{ caller((start,end),flags) }}
+{%- endmacro %}
+
+{# ------------------------------ MESHES --------------------------- #}
+
+{% macro vloop(mesh, lrange, flags, code) %}
+{% set start,end = lrange %}
+{% set thecode, is_top_layer, is_top_inter = code(mesh), False, False %}
+
+{{ define('IS_TOP_LAYER', '_FALSE_') }}
+{{ define('IS_TOP_INTERFACE', '_FALSE_') }}
+
+{% if 'IS_BOTTOM_LEVEL' in thecode or 'KDOWN' in thecode%}
+{# the code in the loop checks whether l==1, for the sake of performance 
+we shall write special code for l=1 and start the loop at l=2 #}
+{{ define('IS_BOTTOM_LEVEL', '_TRUE_') }}
+{{ 'kdown = 1' if 'KDOWN' in thecode }}
+{{ 'kup   = 1' if 'KUP'   in thecode }}
+l=1
+{{ thecode }}
+{% set start='2' %}
+{% endif %}
+
+{% if 'IS_TOP_LAYER' in thecode %}
+{{ 'ERROR : using IS_TOP_LAYER in a loop ending at l=llm+1' if end=='llm+1' }}
+{# the code checks whether l==llm, write special code for l=llm and end the loop at l=llm-1 #}
+{% set end, is_top_layer = 'llm-1', True %}
+{% endif %}
+
+{% if 'IS_TOP_INTERFACE' in thecode or 'KUP' in thecode %}
+{# the code checks whether l==llm+1, write special code for l=llm+1 and end the loop at l=llm #}
+{{ 'ERROR : using IS_TOP_INTERFACE in a loop ending at l=llm' if end=='llm' }}
+{% set end, is_top_inter ='llm', True %}
+{% endif %}
+
+{{ define('IS_BOTTOM_LEVEL', '_FALSE_') }}
+  DO l = {{start}}, {{end}}
+    {{ 'kdown = l-1' if 'KDOWN' in thecode }}
+    {{ 'kup   = l' if 'KUP'   in thecode }}
+    {{ thecode }}
+  END DO
+
+{% if is_top_layer %}
+{{ define('IS_TOP_LAYER', '_TRUE_') }}
+{{ 'kdown = llm-1' if 'KDOWN' in thecode }}
+{{ 'kup   = llm' if 'KUP'   in thecode }}
+l=llm
+{{ thecode }}
+{% endif %}
+
+{% if is_top_inter %}
+{{ define('IS_TOP_INTERFACE', '_TRUE_') }}
+{{ 'kdown = llm' if 'KDOWN' in thecode }}
+{{ 'kup   = llm' if 'KUP'   in thecode }}
+l=llm+1
+{{ thecode }}
+{% endif %}
+
+{% endmacro %}
+
+{% macro on_mesh(mesh,lrange,flags) -%}
+{{ define('CELL','l,ij') if mesh=='primal' }}
+{{ define('DUAL_CELL', 'l,ij') if mesh=='dual'}}
+!$OMP DO SCHEDULE(STATIC)                                                                            
+DO ij = 1, {{ mesh }}_num
+{{ vloop(mesh, lrange, flags, caller) }}
+END DO
+!$OMP END DO
+{{ undef('CELL') }}
+{{ undef('DUAL_CELL') }}
+{%- endmacro %}
+
+{# ------------------------------ STENCILS --------------------------- #}
+
+{% macro on_edges(mesh,lrange,flags) -%}
+{% set thecode = caller(mesh) %}
+{{ define('EDGE', 'l,edge') }}
+{{ define('LE_DE', 'le_de(edge)') }}
+{{ define('SIGN', '1.') }}
+{{ define('CELL1', 'l,ij_left') }}
+{{ define('CELL2',  'l,ij_right') }}
+{{ define('VERTEX1', 'l,ij_down') }}
+{{ define('VERTEX2', 'l,ij_up') }}
+!$OMP DO SCHEDULE(STATIC)
+DO edge = 1, edge_num
+{{ 'ij_left = left(edge)'   if 'CELL1'    in thecode }}
+{{ 'ij_right = right(edge)' if 'CELL2'   in thecode }}
+{{ 'ij_up = up(edge)'       if 'VERTEX1' in thecode }}
+{{ 'ij_down = down(edge)'   if 'VERTEX2' in thecode }}
+{{ vloop(mesh, lrange, flags, caller) }}
+END DO
+!$OMP END DO
+{{ undef('EDGE') }}
+{{ undef('LE_DE') }}
+{{ undef('SIGN') }}
+{{ undef('CELL1') }}
+{{ undef('CELL2') }}
+{{ undef('VERTEX1') }}
+{{ undef('VERTEX2') }}
+{%- endmacro %}
+
+{% macro forall_edges(mesh) -%}
+{% set thecode = caller() %}
+{{ define('EDGE', 'l,edge') }}
+{{ define('VERTEX1', 'l,ij_down') }}
+{{ define('VERTEX2', 'l,ij_up') }}
+{{ define('SIGN', mesh + '_ne(iedge,ij)') }}
+{{ define('LE_DE', 'le_de(edge)') }}
+DO iedge = 1, {{ mesh }}_deg(ij)
+  edge = {{ mesh }}_edge(iedge,ij)
+{{ 'ij_up = up(edge)'       if 'VERTEX1' in thecode }}
+{{ 'ij_down = down(edge)'   if 'VERTEX2' in thecode }}
+  {{ thecode }}
+END DO
+{{ undef('EDGE') }}
+{{ undef('SIGN') }}
+{{ undef('LE_DE') }}
+{%- endmacro %}
+
+{% macro forall_trisk() -%}
+{{ define('EDGE_TRISK', 'l,edge_trisk') }}
+DO itrisk = 1, trisk_deg(edge)
+  edge_trisk = trisk(itrisk,edge)
+  {{ caller() }}
+END DO
+{{ undef('EDGE_TRISK') }}
+{%- endmacro %}
+
+{% macro forall_vertices(mesh) -%}
+{{ define('VERTEX', 'l,vertex') }}
+DO ivertex = 1, {{ mesh }}_deg(ij)
+  vertex = {{ mesh }}_vertex(ivertex,ij)
+  {{ caller() }}
+END DO
+{{ undef('VERTEX') }}
+{%- endmacro %}
+
+{# --------------------------------------------------------- #}
+
+{% set llm='llm' %}
+
+{{ undef('SIGN') }}
+{{ undef('CELL') }}
+{{ undef('CELL1') }}
+{{ undef('CELL2') }}
+{{ undef('VERTEX1') }}
+{{ undef('VERTEX2') }}
+{{ undef('EDGE_TRISK') }}
+
+{# --------------------- END JINJA ------------------------- #}
+
+#endif
+
+#ifdef PASS_POST1
+
+#define _TRUE_ (0==0)
+#define _FALSE_ (0==1)
+
+#define RIV2 Riv2(ivertex,ij)
+#define AI Ai(ij)
+#define AV Av(ij)
+#define FV fv(ij)
+#define WEE wee(itrisk,edge)
+#define edge_ne(iedge,ij) 1.
+#endif
+
+#ifdef PASS_POST2
+#define KUP(l,ij) kup,ij
+#define KDOWN(l,ij) kdown,ij
+#define DOWN(l,ij) l-1,ij
+#define UP(l,ij) l+1,ij
+#define HIDX(l,ij) ij
+#define VIDX(l,ij) l
+#endif
Index: /codes/icosagcm/devel/Python/src/xios.pyx
===================================================================
--- /codes/icosagcm/devel/Python/src/xios.pyx	(revision 615)
+++ /codes/icosagcm/devel/Python/src/xios.pyx	(revision 615)
@@ -0,0 +1,109 @@
+import numpy as np
+import dynamico.wrap as wrap
+from ctypes import c_int, c_double, c_bool, c_void_p, c_char_p, byref, POINTER, Structure
+c_void_pp=POINTER(c_void_p) # used in prototype
+c_void_p_byref=type(byref(c_void_p())) # used in py2c because byref creates an object of this type, not c_void_pp
+
+# ------------- Direct cython interfaces ----------------- #
+
+cdef extern from "functions.h":
+     cpdef void context_close_definition "cxios_context_close_definition"()
+     cpdef void context_finalize "cxios_context_finalize"()
+     cpdef void finalize "cxios_finalize"()
+
+# -------------------------------------------------------- #
+
+lib=wrap.Struct()
+cxios = wrap.SharedLib(vars(lib), "libxios.so", prefix_so='cxios_')
+
+class Duration(Structure):
+    _fields_ = [(name, c_double) for name in ('year', 'month', 'day', 'hour', 'minute', 'second', 'timestep') ]
+class Date(Structure):
+    _fields_ = [(name, c_int) for name in ('year', 'month', 'day', 'hour', 'minute', 'second') ]
+
+print Duration.from_param
+
+for x in c_void_p_byref, Duration, Date : wrap.py2c[x]=wrap.noop # accept arguments of type (void*)&, Duration, ...
+
+def SHAPE(data): return np.asarray(list(reversed(data.shape)),c_int) # reverse shape to follow Fortran convention
+
+class Handle:
+    def __init__(self, cat, id=None):
+        self.handle = c_void_p()
+        if id is None:
+            id='default'
+            create_fun = cxios.vardict['get_current_'+cat]
+            create_fun(byref(self.handle))
+        else:
+            create_fun = cxios.vardict[cat+'_handle_create']
+            create_fun(byref(self.handle), id, c_int(len(id)))
+        self.cat, self.prefix_set, self.id = cat, 'set_'+cat+'_', id
+        print 'Handle', cat, self.id, self.handle
+    def set_attr(self, **kwargs):
+        prefix, handle = self.prefix_set, self.handle
+        for key, value in kwargs.iteritems():
+            fun = cxios.vardict[prefix+key]
+            extra=value
+            if type(value) in (int,c_int):
+                fun(handle,c_int(value))
+            elif type(value) is np.ndarray:
+                extra=SHAPE(value)
+                fun(handle,value,extra)
+            elif type(value) is str:
+                fun(handle,value,c_int(len(value)))
+            elif type(value) in (Duration, Date):
+                fun(handle,value)
+            else:
+                raise(TypeError)
+            print self.id+'.'+prefix+key, type(value), extra
+    def update_timestep(self):
+        lib.update_calendar_timestep(self.handle)
+
+def import_set_attr(cats, *args): #attr_int, attr_array, attr_str):
+    for arg in args:
+        argtype = arg.pop(0)
+        if argtype is str:
+            argtype = [c_void_p, c_char_p, c_int] #char *str, int strlen
+        elif argtype is np.ndarray:
+            argtype = [c_void_p, c_void_p, c_void_p] # void *handle, double *data, int* SHAPE(data)
+        else:
+            argtype = [c_void_p, argtype] #void *handle, argtype data
+        for cat in cats:
+            cxios.import_funs( [arg + argtype], prefix='set_'+cat+'_')
+
+def import_cxios_functions():
+    cats=['axis','domain','domaingroup','field','fieldgroup'] # XIOS categories
+    cxios.import_funs([
+#                       ['finalize','context_finalize','context_close_definition', None],
+                       [cat+'_handle_create' for cat in cats]+[c_void_pp,c_char_p,c_int],
+                       ['get_current_'+'calendar_wrapper',c_void_pp],
+                       ['update_calendar_timestep',c_void_p],
+                       ['update_calendar',c_int] ])
+
+    import_set_attr(['domain','domaingroup'], [c_int,'ni_glo','ibegin','ni','nvertex','data_dim'], [str,'type'],
+                        [np.ndarray, 'i_index','lonvalue_1d','latvalue_1d','bounds_lat_1d','bounds_lon_1d'])
+    import_set_attr(['field','fieldgroup'], [Duration,'freq_op','freq_offset'])
+    import_set_attr(['axis'], [c_int,'n_glo'], [np.ndarray,'value'])
+    import_set_attr(['calendar_wrapper'], [Duration,'timestep'])
+
+    cxios.import_funs([['write_data_k80', c_char_p, c_int, c_void_p],
+                       ['write_data_k81', c_char_p, c_int, c_void_p, c_int],
+                       ['write_data_k82', c_char_p, c_int, c_void_p, c_int, c_int],
+                       ['write_data_k83', c_char_p, c_int, c_void_p, c_int, c_int, c_int]])
+
+def send_field(id, data):
+    if type(data) is np.ndarray:
+        shp=data.shape
+        print 'cxios.write_data', id, shp
+        ndim=len(shp)
+        if ndim==1:
+            lib.write_data_k81(id, c_int(len(id)), data, shp[0])
+        if ndim==2:
+            lib.write_data_k82(id, c_int(len(id)), data, shp[0], shp[1]) # reverse order of dimensions ??
+        if ndim==3:
+            lib.write_data_k83(id, c_int(len(id)), data, shp[0], shp[1], shp[2])
+    else:
+        print id, type(data)
+        raise(TypeError)
+
+import_cxios_functions()
Index: /codes/icosagcm/devel/arch/arch-X64_CURIE.env
===================================================================
--- /codes/icosagcm/devel/arch/arch-X64_CURIE.env	(revision 614)
+++ /codes/icosagcm/devel/arch/arch-X64_CURIE.env	(revision 615)
@@ -3,4 +3,8 @@
 module load mpi/bullxmpi/1.2.8.4
 module load netcdf/4.2_hdf5_parallel
-module load vampir
 module load vampirtrace
+
+# for make_python
+module load parmetis/4.0.3
+module load python/2.7.12
+export CURRENT_ARCH=X64_CURIE
Index: /codes/icosagcm/devel/arch/arch-X64_CURIE.fcm
===================================================================
--- /codes/icosagcm/devel/arch/arch-X64_CURIE.fcm	(revision 614)
+++ /codes/icosagcm/devel/arch/arch-X64_CURIE.fcm	(revision 615)
@@ -6,5 +6,5 @@
 %FPP_FLAGS           -P
 %FPP_DEF             KEY_NONE CPP_NETCDF4
-%BASE_FFLAGS         -i4 -r8 -auto -align all -I${MKLROOT}/include -fPIC
+%BASE_FFLAGS         -i4 -r8 -auto -align all -I${MKLROOT}/include
 %PROD_FFLAGS         -g -traceback -O3 -vec-report2 -xAVX
 %DEV_FFLAGS          -g -O1 -traceback
Index: /codes/icosagcm/devel/arch/arch-X64_CURIE.path
===================================================================
--- /codes/icosagcm/devel/arch/arch-X64_CURIE.path	(revision 614)
+++ /codes/icosagcm/devel/arch/arch-X64_CURIE.path	(revision 615)
@@ -11,4 +11,5 @@
 HDF5_LIB="-lhdf5_hl -lhdf5 -lhdf5 -lz -lcurl"
 
+XIOS_ROOTDIR="$PWD/../XIOS"
 XIOS_INCDIR="-I $PWD/../XIOS/inc"
 XIOS_LIBDIR="-L $PWD/../XIOS/lib"
Index: /codes/icosagcm/devel/bld.cfg
===================================================================
--- /codes/icosagcm/devel/bld.cfg	(revision 614)
+++ /codes/icosagcm/devel/bld.cfg	(revision 615)
@@ -15,21 +15,23 @@
 dir::root            $PWD
 
-src::ioipsl	%IOIPSL_PATH
-src::src	src
-src::vertical	src/vertical
-src::base	src/base
-src::sphere	src/sphere
-src::parallel	src/parallel
-src::output	src/output
-src::initial	src/initial
-src::dcmip	src/dcmip
-src::physics	src/physics
-src::dissip	src/dissip
-src::dynamics	src/dynamics
-src::transport	src/transport
+src::ioipsl		%IOIPSL_PATH
+src::src		src
+src::vertical		src/vertical
+src::base		src/base
+src::sphere		src/sphere
+src::parallel		src/parallel
+src::output		src/output
+src::initial		src/initial
+src::dcmip		src/dcmip
+src::physics		src/physics
+src::dissip		src/dissip
+src::dynamics		src/dynamics
+src::transport		src/transport
 src::diagnostics	src/diagnostics
-src::time	src/time
+src::time		src/time
+src::unstructured	src/unstructured
 src::kernels_hex	src/kernels_hex
-search_src           false
+src::kernels_unst	src/kernels_unst
+search_src		false
 #bld::lib::ioipsl ioipsl
 #bld::target libioipsl.a
Index: /codes/icosagcm/devel/make_python
===================================================================
--- /codes/icosagcm/devel/make_python	(revision 615)
+++ /codes/icosagcm/devel/make_python	(revision 615)
@@ -0,0 +1,86 @@
+#!/bin/bash
+
+function log()
+{
+    echo $*
+    $*
+}
+
+function full_dir()
+{
+    (cd -P $1 ; pwd)
+}
+
+#--------------- Generate kernels, invoked by : './make_python kernels' ---------------
+
+function update()
+{
+ diff $2 $1 || log cp -pr $1 $2
+}
+
+function cmd_kernels() # this function is invoked by : './make_python kernels'
+{
+    cd $KERNELS
+    ./codegen hexagonal unstructured
+
+    cd $KERNELS/src/hexagonal
+    for x in *.k90 ; do
+	update "$x" "$DYNAMICO_ROOT/src/kernels_hex/$x"
+    done
+
+    cd $KERNELS/src/unstructured
+    for x in *.k90 ; do
+	update "$x" "$DYNAMICO_ROOT/src/kernels_unst/$x"
+    done
+
+# make sure to recompile dynamics
+    cd $DYNAMICO_ROOT
+    touch src/dynamics/*.F90 src/diagnostics/*.F90 src/transport/*.F90 src/unstructured/*.F90}
+}
+
+#------------------ Default action , invoked by : './make_python' --------------------
+
+function make_so()
+{
+    ARLIB=$1
+    SOLIB=$2
+    ls -l $ARLIB
+    ls -l $SOLIB
+    if [[ $ARLIB -nt $SOLIB ]]; then
+        echo $ARLIB is newer than $SOLIB
+        log $LINK $SOFLAGS $WHOLE_ARCHIVE $ARLIB $NO_WHOLE_ARCHIVE $3 -o $SOLIB
+        $LDD $SOLIB
+    fi
+}
+
+
+function cmd_()
+{
+    make_so $XIOS_ROOTDIR/lib/libxios.a $LIBDIR/libxios.so "$NETCDF_LDFLAGS"
+    make_so $LIBDIR/libicosa.a $LIBDIR/libicosa.so "-L$LIBDIR -lxios"
+
+    # Cythonize Python modules
+    cd $DYNAMICO_ROOT/Python/src
+    CYTHONFLAGS="$SOFLAGS $LDFLAGS -L$LIBDIR -licosa"
+    LDSHARED="$LINK" CFLAGS="-fno-strict-aliasing -O2" LDFLAGS="$CYTHONFLAGS" log python setup.py build_ext --inplace
+    mv *.so $DYNAMICO_ROOT/Python/dynamico
+    cd ..
+    python -c 'import dynamico ; help(dynamico)'
+}
+
+#------------------------- script starts here ------------------------
+
+DYNAMICO_ROOT=$PWD
+LIBDIR=$DYNAMICO_ROOT/lib
+KERNELS=$DYNAMICO_ROOT/Python
+
+. arch.env
+. arch.path
+
+echo "Selected architecture : $CURRENT_ARCH"
+
+if [[ "$CURRENT_ARCH" == "" ]] ; then echo "Current architecture undefined or unsupported" ; exit ; fi
+
+. Python/env/${CURRENT_ARCH}.env
+
+cmd_$1
Index: /codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90
===================================================================
--- /codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 615)
+++ /codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 615)
@@ -0,0 +1,445 @@
+MODULE caldyn_unstructured_mod
+  USE ISO_C_BINDING
+  USE OMP_LIB
+  USE xios
+  IMPLICIT NONE
+  PRIVATE
+  SAVE
+
+#define BINDC(thename) BIND(C, name='dynamico_' // #thename)
+!
+#define DBL REAL(C_DOUBLE)
+#define DOUBLE1(m) DBL, DIMENSION(m)
+#define DOUBLE2(m,n) DBL, DIMENSION(m,n)
+#define DOUBLE3(m,n,p) DBL, DIMENSION(m,n,p)
+#define INDEX INTEGER(C_INT)
+#define PHI_BOT(ij) Phi_bot
+#define PHI_BOT_VAR 0.d0
+!
+  INTEGER, PARAMETER :: eta_mass=1, eta_lag=2, thermo_theta=1, thermo_entropy=2, thermo_moist=3, thermo_boussinesq=4
+  INTEGER(C_INT),  BIND(C) :: caldyn_thermo=thermo_theta, caldyn_eta=eta_lag, nb_threads=1
+  LOGICAL(C_BOOL), BIND(C) :: hydrostatic=.TRUE., debug_hevi_solver=.TRUE., rigid=.TRUE.
+!
+  INDEX, BIND(C) :: llm, nqdyn, edge_num, primal_num, dual_num, &
+       max_primal_deg, max_dual_deg, max_trisk_deg
+  INDEX, ALLOCATABLE, PRIVATE :: & ! deg(ij) = nb of vertices = nb of edges of primal/dual cell ij
+       primal_deg(:), primal_edge(:,:), primal_vertex(:,:), primal_ne(:,:), & 
+       dual_deg(:), dual_edge(:,:), dual_vertex(:,:), dual_ne(:,:), &
+       trisk_deg(:), trisk(:,:), &
+       left(:), right(:), up(:), down(:)
+  ! left and right are adjacent primal cells
+  ! flux is positive when going from left to right
+  ! up and down are adjacent dual cells
+  ! circulation is positive when going from down to up
+! 
+  DBL, BIND(C) :: elapsed, g, ptop, cpp, cppv, Rd, Rv, preff, Treff, &
+      pbot, Phi_bot, rho_bot
+  DBL :: kappa
+  DOUBLE1(:), ALLOCATABLE, PRIVATE   :: le_de, fv, Av, Ai
+  DOUBLE2(:,:), ALLOCATABLE, PRIVATE :: Riv2, wee, ap,bp, mass_bl, mass_dak, mass_dbk
+!
+  INTEGER(C_INT), BIND(C) :: comm_icosa
+  TYPE(xios_context) :: ctx_hdl
+
+CONTAINS
+!
+!-------------------------------------- KERNELS --------------------------------------
+!
+!
+#define DECLARE_INDICES INTEGER ij,l,iq,iedge,edge,ivertex,vertex,ij_left,ij_right,ij_up,ij_down,itrisk,edge_trisk,kup,kdown
+
+#define FIELD_PS      DOUBLE1(primal_num)
+#define FIELD_MASS    DOUBLE2(llm, primal_num)
+#define FIELD_Z       DOUBLE2(llm, dual_num)
+#define FIELD_U       DOUBLE2(llm, edge_num)
+#define FIELD_UL      DOUBLE2(llm+1, edge_num)
+#define FIELD_THETA   DOUBLE3(llm, primal_num, nqdyn)
+#define FIELD_GEOPOT  DOUBLE2(llm+1, primal_num)
+!
+#define HASNAN(field) (ANY(.NOT.ABS(field)<1e20))
+!
+SUBROUTINE caldyn_unstructured(tau, mass_col,rhodz,theta_rhodz,u,geopot,w, & ! IN : flow state
+                               theta,ps,pk,hflux,qv, &    ! OUT : diags (except surface geopot : IN)
+                               dmass_col,drhodz,dtheta_rhodz,du_fast,du_slow, &
+                               dPhi_fast, dPhi_slow, dW_fast, dW_slow) BINDC(caldyn_unstructured) ! OUT : tendencies
+  DBL, VALUE :: tau
+  FIELD_MASS   :: rhodz, drhodz, pk, berni         ! IN, OUT, DIAG
+  FIELD_THETA  :: theta_rhodz, dtheta_rhodz, theta ! IN, OUT, DIAG
+  FIELD_GEOPOT :: wflux, w, geopot, &              ! DIAG, INOUT
+       dPhi_fast, dPhi_slow, dW_fast, dW_slow      ! OUT
+  FIELD_U      :: u,du_fast,du_slow,hflux,qu       ! INOUT,OUT,OUT,DIAG
+  FIELD_Z      :: qv                               ! DIAG
+  FIELD_PS     :: ps,dmass_col,mass_col            ! OUT,OUT,IN (if eta_mass) or OUT,UNUSED,UNUSED (if eta_lag)
+  DOUBLE2(llm+1, edge_num) :: wwuu
+  DBL          :: time1,time2
+  INTEGER :: ij
+
+  !CALL CPU_TIME(time1)
+  time1=OMP_GET_WTIME()
+
+  IF(hydrostatic) THEN
+
+    !$OMP PARALLEL NUM_THREADS(nb_threads)
+    !$OMP DO SCHEDULE(STATIC)
+    DO ij=1,edge_num
+      du_fast(:,ij)=0.
+      du_slow(:,ij)=0.
+    END DO
+    !$OMP END DO
+    CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
+    CALL compute_geopot(rhodz,theta, ps,pk,geopot)
+
+    CALL compute_caldyn_fast(tau, pk,berni,theta,geopot, du_fast,u)
+
+    CALL compute_pvort_only(rhodz,u,qv,qu)
+    CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow)
+    CALL compute_coriolis(hflux,theta,qu, drhodz,dtheta_rhodz,du_slow)
+    IF(caldyn_eta == eta_mass) THEN
+       CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,wwuu)
+    END IF
+    !$OMP END PARALLEL
+
+  ELSE ! NH
+    DO ij=1,edge_num
+      du_fast(:,ij)=0.
+      du_slow(:,ij)=0.
+    END DO
+    DO ij=1,primal_num
+      wflux(1,ij)=0.
+      wflux(llm+1,ij)=0.
+    END DO
+    CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
+    CALL compute_caldyn_solver(tau,rhodz,theta,pk,geopot,W,dPhi_fast,dW_fast,du_fast)
+    CALL compute_caldyn_fast(tau, pk,berni,theta,geopot, du_fast,u)
+    CALL compute_pvort_only(rhodz,u,qv,qu)
+    CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, hflux,du_slow,dPhi_slow,dW_slow) 
+    CALL compute_coriolis(hflux,theta,qu, drhodz,dtheta_rhodz,du_slow)
+    IF(caldyn_eta == eta_mass) THEN
+       CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,wwuu)
+       CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, du_slow,dPhi_slow,dW_slow)
+    END IF
+  END IF
+ 
+  time2=OMP_GET_WTIME()
+!  CALL CPU_TIME(time2)
+  elapsed = elapsed + time2-time1
+END SUBROUTINE caldyn_unstructured
+!
+!----------------------------- Non-Hydrostatic -----------------------------
+!
+SUBROUTINE compute_NH_geopot(tau,dummy, m_ik, m_il, theta, W_il, Phi_il)
+  FIELD_MASS   :: m_ik, theta, p_ik, A_ik, C_ik   ! IN*2,LOCAL*3
+  FIELD_GEOPOT :: m_il, W_il, Phi_il, Phi_star_il, R_il, x_il, B_il, D_il  ! IN,INOUT*2, LOCAL*5 
+  DBL :: tau, dummy, gamma, rho_ij, X_ij, Y_ij, wil, tau2_g, g2, gm2, ml_g2, c2_mik
+  DECLARE_INDICES
+  INTEGER :: iter
+#include "../kernels_unst/compute_NH_geopot.k90"
+END SUBROUTINE compute_NH_geopot
+!
+SUBROUTINE compute_caldyn_slow_NH(u,rhodz,Phi,W, hflux,du,dPhi,dW)
+  FIELD_U      :: u, hflux, du   ! IN, OUT, OUT
+  FIELD_MASS   :: rhodz, berni   ! IN, LOCAL
+  FIELD_GEOPOT :: Phi,W,dPhi,dW, w_il, gradPhi2  ! IN,IN, OUT,OUT, LOCAL 
+  FIELD_UL     :: DePhil, v_el, G_el, F_el ! LOCAL
+  DECLARE_INDICES
+  DBL :: W_el, W2_el, gPhi2, dP, divG, u2, uu
+#include "../kernels_unst/caldyn_slow_NH.k90"
+END SUBROUTINE compute_caldyn_slow_NH
+!
+SUBROUTINE compute_caldyn_solver(tau,rhodz,theta,pk,geopot,W,dPhi,dW,du)
+  DBL, INTENT(IN) :: tau
+  FIELD_MASS   :: rhodz,pk,berni,pres    ! IN, OUT, LOCAL
+  FIELD_THETA  :: theta                  ! IN
+  FIELD_GEOPOT :: geopot,W,dPhi,dW, m_il ! INOUT,INOUT, OUT,OUT, LOCAL 
+  FIELD_U      :: du                     ! OUT
+  DECLARE_INDICES
+  DBL :: X_ij, rho_ij, T_ij, gamma, Cvd, vreff
+#include "../kernels_unst/caldyn_solver.k90"
+END SUBROUTINE compute_caldyn_solver
+!
+SUBROUTINE compute_caldyn_vert_NH(mass,geopot,W,wflux, du,dPhi,dW)
+  FIELD_MASS   :: mass, eta_dot, wcov, W_etadot
+  FIELD_GEOPOT :: geopot,W,wflux,dPhi,dW
+  FIELD_U      :: du
+  DECLARE_INDICES
+  DBL :: w_ij, wflux_ij
+#include "../kernels_unst/caldyn_vert_NH.k90"
+END SUBROUTINE compute_caldyn_vert_NH
+!
+!----------------------------- Hydrostatic -------------------------------
+!
+SUBROUTINE compute_geopot(rhodz,theta,ps,pk,geopot)
+  FIELD_MASS  :: rhodz,pk   ! IN, OUT
+  FIELD_THETA :: theta      ! IN
+  FIELD_GEOPOT :: geopot    ! IN(l=1)/OUT(l>1)
+  FIELD_PS     :: ps        ! OUT
+  DECLARE_INDICES
+  DBL :: gdz, ke, uu, chi, gv, exner_ik, temp_ik, p_ik, qv, Rmix
+#include "../kernels_unst/compute_geopot.k90"
+END SUBROUTINE compute_geopot
+!
+SUBROUTINE compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du)
+  FIELD_MASS  :: rhodz,berni! IN
+  FIELD_THETA :: theta      ! IN
+  FIELD_U     :: u,hflux,du ! IN, OUT, OUT
+  DECLARE_INDICES
+  LOGICAL, PARAMETER :: zero=.TRUE.
+  DBL :: ke, uu
+#include "../kernels_unst/caldyn_slow_hydro.k90"
+END SUBROUTINE compute_caldyn_slow_hydro
+!
+!---------------------------------- Generic ------------------------------
+!
+SUBROUTINE caldyn_vert(convm,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du, wwuu)
+  FIELD_PS    :: dmass_col
+  FIELD_MASS  :: convm, rhodz
+  FIELD_U     :: u,du
+  FIELD_THETA :: dtheta_rhodz, theta
+  FIELD_GEOPOT :: wflux
+  DOUBLE2(llm+1, edge_num) :: wwuu
+  DECLARE_INDICES
+  wwuu=0.
+#include "../kernels_unst/caldyn_vert.k90"
+END SUBROUTINE caldyn_vert
+!
+SUBROUTINE compute_coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
+  FIELD_U     :: hflux, Ftheta, qu, du
+  FIELD_MASS  :: convm
+  FIELD_THETA :: theta, dtheta_rhodz
+  DECLARE_INDICES
+  DBL :: divF, du_trisk
+#include "../kernels_unst/coriolis.k90"
+END SUBROUTINE
+!
+SUBROUTINE compute_theta(mass_col,rhodz,theta_rhodz, theta)
+  FIELD_PS :: mass_col
+  FIELD_MASS :: rhodz
+  FIELD_THETA :: theta, theta_rhodz
+  DECLARE_INDICES
+  DBL :: m
+#include "../kernels_unst/theta.k90"
+END SUBROUTINE
+!
+SUBROUTINE compute_pvort_only(rhodz,u,qv,qu)
+  FIELD_MASS :: rhodz
+  FIELD_U    :: u,qu
+  FIELD_Z    :: qv
+  DECLARE_INDICES
+  DBL :: etav, hv
+#include "../kernels_unst/pvort_only.k90"
+END SUBROUTINE compute_pvort_only
+!
+SUBROUTINE compute_caldyn_fast(tau, pk,berni,theta,geopot, du,u)
+  DBL, INTENT(IN) :: tau
+  FIELD_MASS   :: pk,berni  ! INOUT, OUT
+  FIELD_THETA  :: theta     ! IN
+  FIELD_GEOPOT :: geopot    ! IN
+  FIELD_U      :: u,du      ! INOUT,INOUT
+  DECLARE_INDICES
+  DBL          :: due
+
+#include "../kernels_unst/caldyn_fast.k90"
+
+END SUBROUTINE compute_caldyn_fast
+!
+!----------------------------- Unused -----------------------------
+!
+SUBROUTINE gradient(b,grad) BINDC(gradient)
+  DOUBLE2(llm,primal_num) :: b
+  DOUBLE2(llm,edge_num)  :: grad
+  DECLARE_INDICES
+#include "../kernels_unst/gradient.k90"
+END SUBROUTINE
+!
+SUBROUTINE div(u,divu) BINDC(div)
+  DOUBLE2(llm,primal_num) :: divu
+  DOUBLE2(llm,edge_num)  :: u
+  DECLARE_INDICES
+  DBL :: div_ij
+  !$OMP PARALLEL NUM_THREADS(nb_threads)
+#include "../kernels_unst/div.k90"
+  !$OMP END PARALLEL
+END SUBROUTINE
+!
+!---------------------------- CONTEXT INITIALIZATION --------------------------
+!
+#define ALLOC1(v,n1) IF(ALLOCATED(v)) DEALLOCATE(v) ; ALLOCATE(v(n1))
+#define ALLOC2(v,n1,n2) IF(ALLOCATED(v)) DEALLOCATE(v) ; ALLOCATE(v(n1,n2))
+
+SUBROUTINE init_mesh( & 
+     primal_deg_, primal_edge_, primal_ne_, &
+     dual_deg_, dual_edge_, dual_ne_, dual_vertex_, &
+     left_, right_, up_, down_ ,&
+     trisk_deg_, trisk_) BINDC(init_mesh)
+  INDEX :: primal_deg_(primal_num), primal_edge_(max_primal_deg,primal_num), &
+       primal_ne_(max_primal_deg,primal_num), &
+       dual_deg_(dual_num), dual_edge_(max_dual_deg,dual_num), &
+       dual_ne_(max_dual_deg,dual_num), &
+       dual_vertex_(max_dual_deg,dual_num), &
+       trisk_deg_(edge_num), trisk_(max_trisk_deg, edge_num)
+  INDEX, DIMENSION(edge_num) :: left_, right_, down_, up_
+
+  PRINT *, 'init_mesh ...'
+  PRINT *, 'Primal mesh : ', primal_num, max_primal_deg
+  PRINT *, 'Dual mesh   : ', dual_num, max_dual_deg
+  PRINT *, 'Edge mesh   : ', edge_num, max_trisk_deg
+  PRINT *, 'Vertical levels :', llm
+  ALLOC1(primal_deg, primal_num)
+  ALLOC2(primal_edge, max_primal_deg,primal_num)
+  ALLOC2(primal_ne, max_primal_deg,primal_num)
+  ALLOC1(dual_deg,dual_num)
+  ALLOC2(dual_edge, max_dual_deg,dual_num)
+  ALLOC2(dual_ne, max_dual_deg,dual_num)
+  ALLOC2(dual_vertex, max_dual_deg,dual_num)
+  ALLOC1(trisk_deg, edge_num)
+  ALLOC2(trisk, max_trisk_deg, edge_num)
+  PRINT *, SHAPE(trisk), edge_num
+  ALLOC1(left, edge_num)
+  ALLOC1(right, edge_num)
+  ALLOC1(up, edge_num)
+  ALLOC1(down, edge_num)
+  primal_deg(:) = primal_deg_(:)
+  primal_edge(:,:) = primal_edge_(:,:)
+  primal_ne(:,:) = primal_ne_(:,:)
+  dual_deg(:) = dual_deg_(:)
+  dual_edge(:,:) = dual_edge_(:,:)
+  dual_ne(:,:) = dual_ne_(:,:)
+  dual_vertex(:,:) = dual_vertex_(:,:)
+  IF(MINVAL(dual_deg)<3) THEN
+     STOP 'At least one dual cell has less than 3 vertices'
+  END IF
+  IF(MINVAL(primal_deg)<3) THEN
+     STOP 'At least one primal cell has less than 3 vertices'
+  END IF
+  left(:)=left_(:)
+  right(:)=right_(:)
+  down(:)=down_(:)
+  up=up_(:)
+  trisk_deg(:)=trisk_deg_(:)
+  trisk(:,:)=trisk_(:,:)
+  PRINT *, MAXVAL(primal_edge), edge_num
+  PRINT *, MAXVAL(dual_edge), edge_num
+  PRINT *, MAXVAL(dual_vertex), dual_num
+  PRINT *, MAXVAL(trisk), edge_num
+  PRINT *, MAX(MAXVAL(left),MAXVAL(right)), primal_num
+  PRINT *, MAX(MAXVAL(up),MAXVAL(down)), dual_num
+  PRINT *, SHAPE(trisk), edge_num
+  PRINT *,' ... Done.'
+END SUBROUTINE init_mesh
+
+SUBROUTINE init_metric(Ai_, Av_, fv_, le_de_, Riv2_, wee_) BINDC(init_metric)
+  DOUBLE1(primal_num) :: Ai_
+  DOUBLE1(dual_num)   :: Av_, fv_
+  DOUBLE1(edge_num)   :: le_de_
+  DOUBLE2(max_dual_deg,dual_num) :: Riv2_
+  DOUBLE2(max_trisk_deg,edge_num) :: wee_
+  PRINT *, 'init_metric ...'
+  ALLOC1(Ai,primal_num)
+  ALLOC1(Av,dual_num)
+  ALLOC1(fv,dual_num)
+  ALLOC1(le_de,edge_num)
+  ALLOC2(Riv2, max_dual_deg, dual_num)
+  ALLOC2(wee, max_trisk_deg, edge_num)
+  Ai(:) = Ai_(:)
+  Av(:) = Av_(:)
+  fv(:) = fv_(:)
+  le_de(:) = le_de_(:)
+  Riv2(:,:)=Riv2_(:,:)
+  wee(:,:) = wee_(:,:)
+  PRINT *, MAXVAL(ABS(Ai))
+  PRINT *, MAXVAL(ABS(Av))
+  PRINT *, MAXVAL(ABS(fv))
+  PRINT *, MAXVAL(ABS(le_de))
+  PRINT *, MAXVAL(ABS(Riv2))
+  PRINT *, MAXVAL(ABS(wee))
+  PRINT *, MINVAL(right), MAXVAL(right)
+  PRINT *, MINVAL(right), MAXVAL(left)
+  PRINT *,' ... Done.'
+  PRINT *,'OpenMP : max_threads, num_procs, nb_threads', OMP_GET_MAX_THREADS(), omp_get_num_procs(), nb_threads
+
+END SUBROUTINE init_metric
+!
+SUBROUTINE init_params() BINDC(init_params)
+  PRINT *, 'Setting physical parameters ...'
+  IF(hydrostatic) THEN
+     PRINT *, 'Hydrostatic dynamics (HPE)'
+  ELSE
+     PRINT *, 'Non-hydrostatic dynamics (Euler)'
+  END IF
+  kappa = Rd/cpp
+  PRINT *, 'g = ',g
+  PRINT *, 'preff = ',preff
+  PRINT *, 'Treff = ',Treff
+  PRINT *, 'Rd = ',Rd
+  PRINT *, 'cpp = ',cpp
+  PRINT *, 'kappa = ',kappa
+  PRINT *, '... Done'
+END SUBROUTINE init_params
+!
+SUBROUTINE init_hybrid(bl,dak,dbk) BINDC(init_hybrid)
+  DOUBLE2(llm+1, primal_num) :: bl
+  DOUBLE2(llm, primal_num) :: dak,dbk
+  PRINT *, 'Setting hybrid coefficients ...'
+  ALLOC2(mass_bl, llm+1, primal_num)
+  ALLOC2(mass_dak, llm, primal_num)
+  ALLOC2(mass_dbk, llm, primal_num)
+  mass_bl(:,:)  = bl(:,:)
+  mass_dak(:,:) = dak(:,:)
+  mass_dbk(:,:) = dbk(:,:)
+  PRINT *, '... Done, llm = ', llm
+END SUBROUTINE Init_hybrid
+
+!---------------------------------------------- XIOS -----------------------------------------
+
+  SUBROUTINE setup_xios() BINDC(setup_xios)
+    ! MPI_INIT / MPI_finalize are assumed to be called BEFORE/AFTER this routine
+    ! This is the case when calling from Python after importing mpi4py
+    INTEGER :: ierr, mpi_threading_mode
+    
+    PRINT *, 'Initialize XIOS and obtain our communicator'
+    CALL xios_initialize("icosagcm",return_comm=comm_icosa)
+    
+    PRINT *, 'Initialize our XIOS context'
+    
+    CALL xios_context_initialize("icosagcm",comm_icosa)
+    CALL xios_get_handle("icosagcm",ctx_hdl)
+    CALL xios_set_current_context(ctx_hdl)
+    !   CALL xios_set_axis_attr("lev",n_glo=llm ,value=lev_value) ;
+    !   CALL xios_set_axis_attr("levp1",n_glo=llm+1 ,value=lev_valuep1) ;
+    !   CALL xios_set_axis_attr("nq",n_glo=nqtot, value=nq_value) ;
+    !   CALL xios_set_domaingroup_attr("i",ni_glo=ncell_tot, ibegin=displ, ni=ncell)
+    !   CALL xios_set_domaingroup_attr("i", data_dim=1, type='unstructured' , nvertex=6, i_index=ind_glo)
+    !   CALL xios_set_domaingroup_attr("i",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat)
+    !   CALL xios_set_domain_attr("u",ni_glo=ncell_tot, ibegin=displ, ni=ncell)
+    !   CALL xios_set_domain_attr("u", data_dim=1, type='unstructured' , nvertex=2, i_index=ind_glo)
+    !  CALL xios_set_domain_attr("u",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat)
+    !   CALL xios_set_domain_attr("v",ni_glo=ncell_tot, ibegin=displ, ni=ncell)
+    !   CALL xios_set_domain_attr("v", data_dim=1, type='unstructured' , nvertex=3)
+    !   CALL xios_set_domain_attr("v",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat)
+    !   CALL xios_set_timestep(dtime)  
+    !   CALL xios_set_fieldgroup_attr("standard_output", freq_op=itau_out*xios_timestep, freq_offset=(itau_out-1)*xios_timestep)    
+    ! CALL xios_close_context_definition()
+    
+    ! PRINT *, 'Read data'
+    ! CALL xios_recv_field(name,field)
+    
+    ! PRINT *,'Write data'
+    ! CALL xios_send_field(name,field_tmp)
+    
+   !   PRINT *, 'Flush to disk, clean up and die'
+   !   CALL xios_context_finalize
+  END SUBROUTINE setup_xios
+!
+  SUBROUTINE call_xios_set_timestep(dt) BINDC(xios_set_timestep)
+    DBL, VALUE :: dt
+    TYPE(xios_duration)      :: dtime
+    dtime%second=dt
+    CALL xios_set_timestep(dtime)
+  END SUBROUTINE call_xios_set_timestep
+
+  SUBROUTINE call_xios_update_calendar(step) BINDC(xios_update_calendar)
+    INDEX, value :: step
+    CALL xios_update_calendar(step)
+  END SUBROUTINE call_xios_update_calendar
+
+END MODULE caldyn_unstructured_mod
