Index: codes/icosagcm/devel/Python/dynamico/meshes.py
===================================================================
--- codes/icosagcm/devel/Python/dynamico/meshes.py	(revision 641)
+++ codes/icosagcm/devel/Python/dynamico/meshes.py	(revision 642)
@@ -148,10 +148,14 @@
                   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 field_theta(self,n=1): return squeeze((n,self.nqdyn,self.ny,self.nx,self.llm))
+    def field_mass(self,n=1): return squeeze((n,self.ny,self.nx,self.llm))
+    def field_z(self,n=1): return squeeze((n,self.ny,self.nx,self.llm))
+    def field_w(self,n=1): return squeeze((n,self.ny,self.nx,self.llm+1))
+    def field_u(self,n=1): 
+        if n==1:
+            return np.zeros((self.ny,2*self.nx,self.llm))
+        else:
+            return np.zeros((n,self.ny,2*self.nx,self.llm))
+    def field_ps(self,n=1): return squeeze((n,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
Index: codes/icosagcm/devel/Python/src/functions.h
===================================================================
--- codes/icosagcm/devel/Python/src/functions.h	(revision 641)
+++ codes/icosagcm/devel/Python/src/functions.h	(revision 642)
@@ -7,5 +7,6 @@
 void dynamico_init_params(void);
 
-void dynamico_ARK_step(double *mass_col, double *rhodz, double *theta_rhodz, 
+void dynamico_ARK_step(int nstep,
+		       double *mass_col, double *rhodz, double *theta_rhodz, 
 		       double *u, double *geopot, double *w,
 		       double *theta, double *ps, double *pk, double *hflux, double *qv,
Index: codes/icosagcm/devel/Python/src/unstructured.pyx
===================================================================
--- codes/icosagcm/devel/Python/src/unstructured.pyx	(revision 641)
+++ codes/icosagcm/devel/Python/src/unstructured.pyx	(revision 642)
@@ -18,5 +18,6 @@
 
 cdef extern from "functions.h":
-     cdef void dynamico_ARK_step(double *mass_col, double *rhodz, double *theta_rhodz, 
+     cdef void dynamico_ARK_step(int nstep,
+                                 double *mass_col, double *rhodz, double *theta_rhodz, 
                                  double *u, double *geopot, double *w,
                                  double *theta, double *ps, double *pk, double *hflux, double *qv,
@@ -78,4 +79,5 @@
 cdef dbl_ptr ptr2(double[:,:] data) : return &data[0,0]
 cdef dbl_ptr ptr3(double[:,:,:] data) : return &data[0,0,0]
+cdef dbl_ptr ptr4(double[:,:,:,:] data) : return &data[0,0,0,0]
 cdef dbl_ptr ptr(data):
     n=data.ndim
@@ -83,5 +85,7 @@
     if n==2 : return ptr2(data)
     if n==3 : return ptr3(data)
-
+    if n==4 : return ptr4(data)
+    if n>4: raise IndexError
+        
 cdef alloc(dbl_ptr *p, allocator, n=1):
     data=allocator(n)
@@ -93,16 +97,19 @@
         
 cdef class Caldyn_step:
+    # number of time steps to do at each invocation of advance()
+    cdef int nstep
     # pointer to allocated arrays
-    cdef dbl_ptr p_mass, p_theta_rhodz, p_u, p_geopot, p_w                   # prognostic
+    cdef dbl_ptr p_mass, p_theta_rhodz, p_u, p_geopot, p_W                   # prognostic
     cdef dbl_ptr p_mass_col, p_dmass_col, p_ps, p_theta, p_pk, p_hflux, p_qv # diagnostic
     cdef dbl_ptr p_drhodz, p_dtheta_rhodz, p_du_fast, p_du_slow                # tendencies
     cdef dbl_ptr p_dPhi_fast, p_dPhi_slow, p_dW_fast, p_dW_slow                # tendencies
     # allocated arrays, must remain referenced or segfault
-    cdef readonly np.ndarray mass, theta_rhodz, u, geopot, w
+    cdef readonly np.ndarray mass, theta_rhodz, u, geopot, W
     cdef readonly np.ndarray mass_col, dmass_col, ps, theta, pk, hflux, qv
     cdef readonly np.ndarray drhodz, dtheta_rhodz, du_fast, du_slow
     cdef readonly np.ndarray dPhi_fast, dPhi_slow, dW_fast, dW_slow
 
-    def __init__(self,mesh,time_scheme):
+    def __init__(self,mesh,time_scheme, nstep):
+        self.nstep=nstep
         #        self.mesh=mesh
         fps, ftheta, fmass = mesh.field_ps, mesh.field_theta, mesh.field_mass
@@ -113,5 +120,6 @@
         cdef double[:,:] cflj_ = time_scheme.cfjl
         ns = time_scheme.nstage
-        nb_stage[0]= ns
+        nb_stage[0]=ns
+
         cdef int i,j
         for i in range(ns):
@@ -123,8 +131,8 @@
         #    prognostic/diagnostic
         self.ps                       = alloc(&self.p_ps, fps) 
-        self.mass_col, self.dmass_col = alloc(&self.p_mass_col, fps), alloc(&self.p_dmass_col, fps), 
+        self.mass_col, self.dmass_col = alloc(&self.p_mass_col, fps), alloc(&self.p_dmass_col, fps,ns), 
         self.mass, self.theta_rhodz   = alloc(&self.p_mass, fmass), alloc(&self.p_theta_rhodz, fmass),
         self.theta, self.pk           = alloc(&self.p_theta, fmass), alloc(&self.p_pk, fmass), 
-        self.geopot, self.w           = alloc(&self.p_geopot, fw), alloc(&self.p_w, fw),
+        self.geopot, self.W           = alloc(&self.p_geopot, fw), alloc(&self.p_W, fw),
         self.hflux, self.u            = alloc(&self.p_hflux, fu), alloc(&self.p_u, fu) 
         self.qv                       = alloc(&self.p_qv,fz)
@@ -137,8 +145,7 @@
         #        global elapsed
         # time1=time.time()
-        check_ptr('mass', self.p_mass, self.mass)
-        check_ptr('u', self.p_u, self.u)
-        dynamico_ARK_step(self.p_mass_col, self.p_mass, self.p_theta_rhodz,
-                          self.p_u, self.p_geopot, self.p_w,
+        dynamico_ARK_step(self.nstep,
+                          self.p_mass_col, self.p_mass, self.p_theta_rhodz,
+                          self.p_u, self.p_geopot, self.p_W,
                           self.p_theta, self.p_ps, self.p_pk, self.p_hflux, self.p_qv,
                           self.p_dmass_col, self.p_drhodz, self.p_dtheta_rhodz,
@@ -148,23 +155,26 @@
         #time2=time.time()
         #if time2>time1: elapsed=elapsed+time2-time1
-    def setup_TRSW(self):
-        setvars(('hydrostatic','caldyn_thermo','caldyn_eta'),
-                (True,thermo_boussinesq,eta_lag))
-        dynamico_init_params()
-    def setup_HPE(self, caldyn_thermo, caldyn_eta, g, BC,thermo):
-        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))
-        dynamico_init_params()
-    def setup_NH(self, caldyn_thermo, caldyn_eta, g, BC,thermo):
-        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()))
-        dynamico_init_params()
-
+
+def caldyn_step_TRSW(mesh,time_scheme,nstep):
+    setvars(('hydrostatic','caldyn_thermo','caldyn_eta'),
+            (True,thermo_boussinesq,eta_lag))
+    dynamico_init_params()
+    return Caldyn_step(mesh,time_scheme, nstep)
+def caldyn_step_HPE(mesh,time_scheme,nstep, caldyn_thermo,caldyn_eta, thermo,BC,g):
+    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))
+    dynamico_init_params()
+    return Caldyn_step(mesh,time_scheme, nstep)
+def caldyn_step_NH(mesh,time_scheme,nstep, caldyn_thermo, caldyn_eta, thermo,BC,g):
+    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()))
+    dynamico_init_params()
+    return Caldyn_step(mesh,time_scheme, nstep)
+    
 #----------------------------- Base class for dynamics ------------------------
 
@@ -212,5 +222,5 @@
 
 class Caldyn_HPE(Caldyn):
-    def __init__(self,caldyn_thermo,caldyn_eta, mesh,metric,thermo,BC,g):
+    def __init__(self,caldyn_thermo,caldyn_eta, mesh,thermo,BC,g):
         Caldyn.__init__(self,mesh)
         setvars(('hydrostatic','caldyn_thermo','caldyn_eta',
@@ -231,5 +241,5 @@
 
 class Caldyn_NH(Caldyn):
-    def __init__(self,caldyn_thermo,caldyn_eta, mesh,metric,thermo,BC,g):
+    def __init__(self,caldyn_thermo,caldyn_eta, mesh,thermo,BC,g):
         Caldyn.__init__(self,mesh)
         setvars(('hydrostatic','caldyn_thermo','caldyn_eta',
Index: codes/icosagcm/devel/Python/test/py/NH_3D_bubble.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/NH_3D_bubble.py	(revision 641)
+++ codes/icosagcm/devel/Python/test/py/NH_3D_bubble.py	(revision 642)
@@ -53,6 +53,6 @@
     return thermo, mesh, params, (mass_ik,Sik,ujk,Phi_il,Wil), gas
 
-Lx, nx, llm, thetac, T, Nslice, courant = 2000., 100, 50, 30., 5., 10, 2.8
-#Lx, nx, llm, thetac, T, Nslice, courant = 2000., 50, 25, 1., 50., 10, 2.8
+#Lx, nx, llm, thetac, T, Nslice, courant = 2000., 100, 50, 30., 5., 10, 2.8
+Lx, nx, llm, thetac, T, Nslice, courant = 2000., 50, 25, 30, 5., 10, 2.8
 #Lx, nx, llm, thetac, T, Nslice, courant = 3000., 75, 25, -30, 5., 10, 2.8
 
@@ -61,11 +61,7 @@
 thermo, mesh, params, flow0, gas0 = thermal_bubble_3D(Lx,nx,Ly,ny,llm, thetac=thetac)
 
-caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass
-caldyn = unst.Caldyn_NH(caldyn_thermo,caldyn_eta, mesh,params,thermo,params,params.g)
-
 # compute hybrid coefs from initial distribution of mass
 mass_bl,mass_dak,mass_dbk = meshes.compute_hybrid_coefs(flow0[0])
 unst.ker.dynamico_init_hybrid(mass_bl,mass_dak,mass_dbk)
-print mass_dbk
 
 dz = flow0[3].max()/(params.g*llm)
@@ -74,15 +70,34 @@
 nt = int(math.ceil(T/dt))
 dt = T/nt
-print 'Time step : %g s' % dt
+print 'Time step : %d x %g s' % (nt,dt)
 
-scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt)
+caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass
+#caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_lag
 
-flow=flow0
+if False: # time stepping in Python
+    caldyn = unst.Caldyn_NH(caldyn_thermo,caldyn_eta, mesh,thermo,params,params.g)
+    scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt)
+    def next_flow(m,S,u,Phi,W):
+        # junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
+        return scheme.advance((m,S,u,Phi,W),nt)
+
+else: # time stepping in Fortran
+    scheme = time_step.ARK2(None, dt)
+    caldyn_step = unst.caldyn_step_NH(mesh,scheme,nt, caldyn_thermo,caldyn_eta,
+                                      thermo,params,params.g)
+    def next_flow(m,S,u,Phi,W):
+        # junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
+        caldyn_step.mass[:,:,:], caldyn_step.theta_rhodz[:,:,:], caldyn_step.u[:,:,:] = m,S,u
+        caldyn_step.geopot[:,:,:], caldyn_step.W[:,:,:] = Phi,W
+        caldyn_step.next()
+        return (caldyn_step.mass.copy(), caldyn_step.theta_rhodz.copy(), caldyn_step.u.copy(),
+                caldyn_step.geopot.copy(), caldyn_step.W.copy())
+    
+m,S,u,Phi,W=flow0
 w=mesh.field_mass()
 z=mesh.field_mass()
+
 for it in range(Nslice):
-    junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
-    m,S,u,Phi,W = flow
-    s=S/m ; # s=.5*(s+abs(s))
+    s=S/m ; s=.5*(s+abs(s))
     for l in range(llm):
         w[:,:,l]=.5*params.g*(W[:,:,l+1]+W[:,:,l])/m[:,:,l]
@@ -92,7 +107,7 @@
     jj=ny/2
     xx,zz,ss,ww = mesh.xx[:,jj,:]/1000., z[:,jj,:]/1000., s[:,jj,:], w[:,jj,:]
+
+    f, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(12,4))
     
-    f, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(12,4))
-
     c=ax1.contourf(xx,zz,ss,20)
     ax1.set_xlim((-.5,.5)), ax1.set_xlabel('x (km)')
@@ -100,7 +115,7 @@
     plt.colorbar(c,ax=ax1)
     ax1.set_title('Specific entropy at t=%g s (J/K/kg)' % (it*T,))
-#    plt.show()
+    #    plt.show()
     
-#    plt.figure(figsize=(12,5))
+    #    plt.figure(figsize=(12,5))
     c=ax2.contourf(xx,zz,ww,20)
     ax2.set_xlim((-.5,.5)), ax2.set_xlabel('x (km)')
@@ -108,10 +123,10 @@
     plt.colorbar(c,ax=ax2)
     ax2.set_title('Vertical velocity at t=%g s (m/s)' % (it*T,))
-#    plt.tight_layout()
-#    plt.show()
+    #    plt.tight_layout()
+    #    plt.show()
     plt.savefig('fig_NH_3D_bubble/%02d.png'%it)
     
     time1, elapsed1 =time.time(), unst.getvar('elapsed')
-    flow = scheme.advance(flow,nt)
+    m,S,u,Phi,W = next_flow(m,S,u,Phi,W)
     time2, elapsed2 =time.time(), unst.getvar('elapsed')
     factor = 1000./(4*nt)
Index: codes/icosagcm/devel/Python/test/py/RSW2_MPAS_W02.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/RSW2_MPAS_W02.py	(revision 641)
+++ codes/icosagcm/devel/Python/test/py/RSW2_MPAS_W02.py	(revision 642)
@@ -38,6 +38,5 @@
 #scheme = time_step.RKn_simple(1,caldyn.bwd_fast_slow, dt) # forward Euler scheme
 print dt, scheme.csjl, scheme.cfjl
-step = unst.Caldyn_step(mesh,scheme)
-step.setup_TRSW()
+step = unst.caldyn_step_TRSW(mesh,scheme,nt)
 
 u0 = Omega*radius/12. # cf Williamson (1991), p.13
@@ -63,6 +62,5 @@
 for i in range(20):
     if True:
-        for j in range(nt):
-            step.next()
+        step.next()
         gh,u = step.mass, step.u
         print i, gh.shape, gh.min(), gh.max(), u.min(), u.max()
Index: codes/icosagcm/devel/Python/test/py/slice_GW_hydro.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/slice_GW_hydro.py	(revision 641)
+++ codes/icosagcm/devel/Python/test/py/slice_GW_hydro.py	(revision 642)
@@ -15,5 +15,4 @@
     caldyn_thermo = unst.thermo_entropy
     g = mesh.dx/metric.dx_g0
-    caldyn = unst.Caldyn_HPE(caldyn_thermo,caldyn_eta, mesh,metric,thermo,BC,g)
 
     Sik = m0ik*gas0.s
@@ -25,5 +24,5 @@
         S0ik = m0ik*gas0.s
 
-    u0=mesh.field_u() 
+    u0=mesh.field_u()
     u0[:,range(0,2*nx,2),:] = 100.*mesh.dx # Doppler shift 100 m/s
     flow0=(m0ik,S0ik,u0)
@@ -35,17 +34,30 @@
     dt=T/nt
     print 'nt,dt,dx', nt,dt,dx
-    scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, a32=0.7)
-    #scheme = time_step.ARK3(caldyn.bwd_fast_slow, dt)
 
-    flow=flow0
+    m,S,u=flow0
+
+    if False: # time stepping in Python
+        caldyn = unst.Caldyn_HPE(caldyn_thermo,caldyn_eta, mesh,thermo,BC,g)
+        scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, a32=0.7)
+        #scheme = time_step.ARK3(caldyn.bwd_fast_slow, dt)
+        def next_flow(m,S,u):
+            m,S,u = scheme.advance((m,S,u),nt)
+            return m,S,u,caldyn.geopot
+    else: # time stepping in Fortran
+        scheme = time_step.ARK2(None, dt, a32=0.7)
+        caldyn_step = unst.caldyn_step_HPE(mesh,scheme,nt, caldyn_thermo,caldyn_eta,thermo,BC,g)
+        def next_flow(m,S,u):
+            caldyn_step.mass[:,:], caldyn_step.theta_rhodz[:,:], caldyn_step.u[:,:] = m,S,u
+            caldyn_step.next()
+            return caldyn_step.mass, caldyn_step.theta_rhodz, caldyn_step.u, caldyn_step.geopot
+        
     for i in range(10):
         time1=time.time()
-        flow = scheme.advance(flow,nt)
+        m,S,u,geopot = next_flow(m,S,u)
         time2=time.time()
-        m,S,u = flow
-        print 'ms per time step : ', 1000*(time2-time1)/nt
-        print 'ptop, model top (m) :', unst.getvar('ptop'), caldyn.geopot.max()/unst.getvar('g')
-        junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
-        zz=caldyn.geopot[:,0:llm]/metric.g0/1000
+        print 'ms per time step : ', 1000*(time2-time1)/nt, 1000*unst.getvar('elapsed')/(nt*(i+1))
+        print 'ptop, model top (m) :', unst.getvar('ptop'), geopot.max()/unst.getvar('g')
+        #        junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
+        zz=geopot[:,0:llm]/metric.g0/1000
         plt.figure(figsize=(12,3))
         ucomp = mesh.ucomp(u)
Index: codes/icosagcm/devel/src/icosa_init.f90
===================================================================
--- codes/icosagcm/devel/src/icosa_init.f90	(revision 641)
+++ codes/icosagcm/devel/src/icosa_init.f90	(revision 642)
@@ -107,3 +107,7 @@
   END SUBROUTINE icosa_init
 
+  SUBROUTINE force_recompile_unstructured
+    USE timestep_unstructured_mod
+  END SUBROUTINE force_recompile_unstructured
+
 END MODULE icosa_init_mod
Index: codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90
===================================================================
--- codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 641)
+++ codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 642)
@@ -103,5 +103,5 @@
   FIELD_THETA :: dtheta_rhodz, theta
   FIELD_GEOPOT :: wflux
-  DOUBLE2(llm+1, edge_num) :: wwuu
+  FIELD_UL     :: wwuu
   DECLARE_INDICES
   DBL :: dF_deta, dFu_deta
Index: codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90
===================================================================
--- codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90	(revision 641)
+++ codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90	(revision 642)
@@ -36,79 +36,79 @@
 #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)
+  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)
+    FIELD_UL     :: 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
-!$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)
-  IF(time2>time1) elapsed = elapsed + time2-time1
-END SUBROUTINE caldyn_unstructured
-!
-!----------------------------- Time stepping -------------------------------
-
-!
-  SUBROUTINE ARK_step(mass_col,rhodz,theta_rhodz,u,geopot,w, & ! INOUT : flow state
+    
+    time2=OMP_GET_WTIME()
+    !  CALL CPU_TIME(time2)
+    IF(time2>time1) elapsed = elapsed + time2-time1
+  END SUBROUTINE caldyn_unstructured
+  !
+  !----------------------------- Time stepping -------------------------------
+  
+  !
+  SUBROUTINE ARK_step(nstep, mass_col,rhodz,theta_rhodz,u,geopot,w, & ! INOUT : 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(ARK_step) ! OUT : tendencies
+    INTEGER, VALUE :: nstep ! advance by nstep time steps
     FIELD_PS     :: mass_col, ps                     ! OUT,IN (if eta_mass) or OUT,UNUSED (if eta_lag)
     FIELD_MASS   :: rhodz, pk, berni                 ! IN, DIAG
@@ -123,146 +123,144 @@
     DOUBLE3(llm+1, primal_num, max_nb_stage) :: &
          dPhi_fast, dPhi_slow, dW_fast, dW_slow      ! OUT
-    FIELD_UL     :: DePhil, v_el, G_el               ! DIAG*3
-    DOUBLE2(llm+1, edge_num) :: wwuu
+    FIELD_UL     :: DePhil, v_el, G_el, wwuu         ! DIAG*4
     DBL       :: time1,time2
-    INTEGER :: stage, ij
-
+    INTEGER :: step, stage, ij
+    
     !CALL CPU_TIME(time1)
     time1=OMP_GET_WTIME()
-
+    
     !$OMP PARALLEL NUM_THREADS(nb_threads)
     
-    DO stage=1, nb_stage
-       
-       IF(hydrostatic) THEN
-
-          !$OMP DO SCHEDULE(STATIC)
-          DO ij=1,edge_num
-             du_fast(:,ij,stage)=0.
-             du_slow(:,ij,stage)=0.
-          END DO
-
-          CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
-          CALL compute_geopot(rhodz,theta, ps,pk,geopot)
-          
-          CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage), u)
-          
-          CALL compute_pvort_only(rhodz,u,qv,qu)
-          CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow(:,:,stage))
-          CALL compute_coriolis(hflux,theta,qu, drhodz(:,:,stage), dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage))
-          IF(caldyn_eta == eta_mass) THEN
-             STOP ! FIXME
-             CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, &
-                  dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage),wwuu)
+    DO step=1, nstep
+       DO stage=1, nb_stage
+          
+          IF(hydrostatic) THEN
+             
+             !$OMP DO SCHEDULE(STATIC)
+             DO ij=1,edge_num
+                du_fast(:,ij,stage)=0.
+                du_slow(:,ij,stage)=0.
+             END DO
+             
+             CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
+             CALL compute_geopot(rhodz,theta, ps,pk,geopot)
+             
+             CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage), u)
+             
+             CALL compute_pvort_only(rhodz,u,qv,qu)
+             CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow(:,:,stage))
+             CALL compute_coriolis(hflux,theta,qu, drhodz(:,:,stage), dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage))
+             IF(caldyn_eta == eta_mass) THEN
+                CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, &
+                     dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage),wwuu)
+             END IF
+             
+          ELSE ! NH
+             !$OMP DO SCHEDULE(STATIC)
+             DO ij=1,primal_num
+                wflux(1,ij)=0.
+                wflux(llm+1,ij)=0.
+             END DO
+             !$OMP END DO
+             CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
+             CALL compute_caldyn_solver(tauj(stage),rhodz,theta,pk,geopot,W, &
+                  dPhi_fast(:,:,stage), dW_fast(:,:,stage), du_fast(:,:,stage))
+             CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage),u)
+             CALL compute_pvort_only(rhodz,u,qv,qu)
+             CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, hflux, du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage))
+             CALL compute_coriolis(hflux,theta,qu, drhodz(:,:,stage),dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage))
+             IF(caldyn_eta == eta_mass) THEN
+                CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, &
+                     dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage), wwuu)
+                CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage) )
+             END IF
+          END IF ! NH
+          
+          ! FIXME : 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
+          
+#define UPDATE(clj,field,dfield) update(stage,SIZE(field),clj,field,dfield)
+          
+          CALL UPDATE(cslj, rhodz, drhodz)
+          CALL UPDATE(cslj, theta_rhodz, dtheta_rhodz)
+          CALL UPDATE(cslj, u, du_slow)
+          CALL UPDATE(cflj, u, du_fast)
+          
+          IF(.NOT.hydrostatic) THEN
+             CALL UPDATE(cslj, geopot, dPhi_slow)
+             CALL UPDATE(cflj, geopot, dPhi_fast)
+             CALL UPDATE(cslj, W, dW_slow)
+             CALL UPDATE(cflj, W, dW_fast)
           END IF
           
-       ELSE ! NH
-          STOP ! FIXME
-          !$OMP DO SCHEDULE(STATIC)
-          DO ij=1,primal_num
-             wflux(1,ij)=0.
-             wflux(llm+1,ij)=0.
-          END DO
-          !$OMP END DO
-          CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
-          CALL compute_caldyn_solver(tauj(stage),rhodz,theta,pk,geopot,W, &
-               dPhi_fast(:,:,stage), dW_fast(:,:,stage), du_fast(:,:,stage))
-          CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage),u)
-          CALL compute_pvort_only(rhodz,u,qv,qu)
-          CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, hflux, du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage))
-          CALL compute_coriolis(hflux,theta,qu, drhodz(:,:,stage),dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage))
-          IF(caldyn_eta == eta_mass) THEN
-             CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, &
-                  dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage), wwuu)
-             CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage) )
-          END IF
-       END IF ! NH
-       
-      ! FIXME : 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
-
-#define UPDATE(clj,field,dfield) update(stage,SIZE(field),clj,field,dfield)
-
-       CALL UPDATE(cslj, rhodz, drhodz)
-       CALL UPDATE(cslj, theta_rhodz, dtheta_rhodz)
-       CALL UPDATE(cslj, u, du_slow)
-       CALL UPDATE(cflj, u, du_fast)
-
-       IF(.NOT.hydrostatic) THEN
-          STOP ! FIXME
-          CALL UPDATE(cslj, geopot, dPhi_slow)
-          CALL UPDATE(cflj, geopot, dPhi_fast)
-          CALL UPDATE(cslj, W, dW_slow)
-          CALL UPDATE(cflj, W, dW_fast)
-       END IF
-
+       END DO
     END DO
-!$OMP END PARALLEL
-    
-  time2=OMP_GET_WTIME()
-!  CALL CPU_TIME(time2)
-  IF(time2>time1) elapsed = elapsed + time2-time1
-
+    !$OMP END PARALLEL
+    
+    time2=OMP_GET_WTIME()
+    !  CALL CPU_TIME(time2)
+    IF(time2>time1) elapsed = elapsed + time2-time1
+    
   END SUBROUTINE ARK_step
-  !
-!
-SUBROUTINE update(j,sz,clj,field,dfield)
-  INTEGER :: j, sz ! stage in ARK scheme, field size
-  DOUBLE2(max_nb_stage,max_nb_stage) :: clj ! modified Butcher tableau
-  DOUBLE1(sz) :: field
-  DOUBLE2(sz, max_nb_stage) :: dfield
-  !
-  INTEGER :: l, ij
-!  PRINT *, clj(1:j,j)
-  SELECT CASE(j)
-  CASE(1)
-     !$OMP DO SCHEDULE(static)
-     DO ij=1,sz
-        field(ij) = field(ij) &
-             + clj(1,j)*dfield(ij,1)
-     END DO
-  CASE(2)
-     !$OMP DO SCHEDULE(static)
-     DO ij=1,sz
-        field(ij) = field(ij) &
-             + clj(1,j)*dfield(ij,1) &
-             + clj(2,j)*dfield(ij,2)
-     END DO
-  CASE(3)
-     !$OMP DO SCHEDULE(static)
-     DO ij=1,sz
-        field(ij) = field(ij) &
-             + clj(1,j)*dfield(ij,1) &
-             + clj(2,j)*dfield(ij,2) &
-             + clj(3,j)*dfield(ij,3)
-
-     END DO
-  CASE(4)
-     !$OMP DO SCHEDULE(static)
-     DO ij=1,sz
-        field(ij) = field(ij) &
-             + clj(1,j)*dfield(ij,1) &
-             + clj(2,j)*dfield(ij,2) &
-             + clj(3,j)*dfield(ij,3) &
-             + clj(4,j)*dfield(ij,4)
-     END DO
-  CASE(5)
-     !$OMP DO SCHEDULE(static)
-     DO ij=1,sz
-        field(ij) = field(ij) &
-             + clj(1,j)*dfield(ij,1) &
-             + clj(2,j)*dfield(ij,2) &
-             + clj(3,j)*dfield(ij,3) &
-             + clj(4,j)*dfield(ij,4) &
-             + clj(5,j)*dfield(ij,5)
-     END DO
-  END SELECT
-END SUBROUTINE update
-
-!---------------------------------------------- XIOS ----------------------------------------
+  
+  
+  SUBROUTINE update(j,sz,clj,field,dfield)
+    INTEGER :: j, sz ! stage in ARK scheme, field size
+    DOUBLE2(max_nb_stage,max_nb_stage) :: clj ! modified Butcher tableau
+    DOUBLE1(sz) :: field
+    DOUBLE2(sz, max_nb_stage) :: dfield
+    !
+    INTEGER :: l, ij
+    !  PRINT *, clj(1:j,j)
+    SELECT CASE(j)
+    CASE(1)
+       !$OMP DO SCHEDULE(static)
+       DO ij=1,sz
+          field(ij) = field(ij) &
+               + clj(1,j)*dfield(ij,1)
+       END DO
+    CASE(2)
+       !$OMP DO SCHEDULE(static)
+       DO ij=1,sz
+          field(ij) = field(ij) &
+               + clj(1,j)*dfield(ij,1) &
+               + clj(2,j)*dfield(ij,2)
+       END DO
+    CASE(3)
+       !$OMP DO SCHEDULE(static)
+       DO ij=1,sz
+          field(ij) = field(ij) &
+               + clj(1,j)*dfield(ij,1) &
+               + clj(2,j)*dfield(ij,2) &
+               + clj(3,j)*dfield(ij,3)
+          
+       END DO
+    CASE(4)
+       !$OMP DO SCHEDULE(static)
+       DO ij=1,sz
+          field(ij) = field(ij) &
+               + clj(1,j)*dfield(ij,1) &
+               + clj(2,j)*dfield(ij,2) &
+               + clj(3,j)*dfield(ij,3) &
+               + clj(4,j)*dfield(ij,4)
+       END DO
+    CASE(5)
+       !$OMP DO SCHEDULE(static)
+       DO ij=1,sz
+          field(ij) = field(ij) &
+               + clj(1,j)*dfield(ij,1) &
+               + clj(2,j)*dfield(ij,2) &
+               + clj(3,j)*dfield(ij,3) &
+               + clj(4,j)*dfield(ij,4) &
+               + clj(5,j)*dfield(ij,5)
+       END DO
+    END SELECT
+  END SUBROUTINE update
+  
+  !---------------------------------------------- XIOS ----------------------------------------
 
 #ifdef CPP_USING_XIOS
-
+  
   SUBROUTINE setup_xios() BINDC(setup_xios)
     ! MPI_INIT / MPI_finalize are assumed to be called BEFORE/AFTER this routine
