Index: codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90
===================================================================
--- codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 638)
+++ 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 638)
+++ 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
