Index: branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 1360)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 1361)
@@ -4,5 +4,18 @@
    !! Ocean dynamics: time stepping
    !!======================================================================
-   
+   !!======================================================================
+   !! History :  OPA  !  1987-02  (P. Andrich, D. L Hostis)  Original code
+   !!                 !  1990-10  (C. Levy, G. Madec)
+   !!            7.0  !  1993-03  (M. Guyon)  symetrical conditions
+   !!            8.0  !  1997-02  (G. Madec & M. Imbard)  opa, release 8.0
+   !!            8.2  !  1997-04  (A. Weaver)  Euler forward step
+   !!             -   !  1997-06  (G. Madec)  lateral boudary cond., lbc routine
+   !!    NEMO    1.0  !  2002-08  (G. Madec)  F90: Free form and module
+   !!             -   !  2002-10  (C. Talandier, A-M. Treguier) Open boundary cond.
+   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization
+   !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines. 
+   !!            3.1  !  2009-02  (G. Madec)  re-introduce the vvl option
+   !!----------------------------------------------------------------------
+  
    !!----------------------------------------------------------------------
    !!   dyn_nxt      : update the horizontal velocity from the momentum trend
@@ -34,4 +47,8 @@
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
+   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! $Id$ 
+   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !!----------------------------------------------------------------------
 
 CONTAINS
@@ -59,34 +76,9 @@
       !!             - Update un,vn arrays, the now horizontal velocity
       !!
-      !! History :
-      !!        !  87-02  (P. Andrich, D. L Hostis)  Original code
-      !!        !  90-10  (C. Levy, G. Madec)
-      !!        !  93-03  (M. Guyon)  symetrical conditions
-      !!        !  97-02  (G. Madec & M. Imbard)  opa, release 8.0
-      !!        !  97-04  (A. Weaver)  Euler forward step
-      !!        !  97-06  (G. Madec)  lateral boudary cond., lbc routine
-      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
-      !!        !  02-10  (C. Talandier, A-M. Treguier) Open boundary cond.
-      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
-      !!    "   !  07-07  (D. Storkey) Calls to BDY routines. 
       !!----------------------------------------------------------------------
-      !! * Arguments
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-
-      !! * Local declarations
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      !!
+      INTEGER  ::   jk   ! dummy loop indices
       REAL(wp) ::   z2dt         ! temporary scalar
-      REAL(wp) ::   zsshun1, zsshvn1         ! temporary scalar
-      !! Variable volume
-      REAL(wp), DIMENSION(jpi,jpj)     ::  & ! 2D workspace
-         zsshub, zsshun, zsshua,           &
-         zsshvb, zsshvn, zsshva
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  &
-         zfse3ub, zfse3un, zfse3ua, &        ! 3D workspace
-         zfse3vb, zfse3vn, zfse3va
-      !!----------------------------------------------------------------------
-      !!  OPA 9.0 , LOCEAN-IPSL (2005) 
-      !! $Id$ 
-      !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
       !!----------------------------------------------------------------------
 
@@ -102,72 +94,4 @@
 
       !! Explicit physics with thickness weighted updates
-      IF( lk_vvl .AND. .NOT. lk_dynspg_flt ) THEN
-
-         ! Sea surface elevation time stepping
-         ! -----------------------------------
-         !
-         DO jj = 1, jpjm1
-            DO ji = 1,jpim1
-
-               ! Sea Surface Height at u-point before
-               zsshub(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )     &
-                  &          * ( e1t(ji  ,jj  ) * e2t(ji  ,jj  ) * sshbb(ji  ,jj  )   &
-                  &            + e1t(ji+1,jj  ) * e2t(ji+1,jj  ) * sshbb(ji+1,jj  ) )
-
-               ! Sea Surface Height at v-point before
-               zsshvb(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )     &
-                  &          * ( e1t(ji  ,jj  ) * e2t(ji  ,jj  ) * sshbb(ji  ,jj  )   &
-                  &            + e1t(ji  ,jj+1) * e2t(ji  ,jj+1) * sshbb(ji  ,jj+1) )
-
-               ! Sea Surface Height at u-point after
-               zsshua(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )     &
-                  &          * ( e1t(ji  ,jj  ) * e2t(ji  ,jj  ) * ssha(ji  ,jj  )    &
-                  &            + e1t(ji+1,jj  ) * e2t(ji+1,jj  ) * ssha(ji+1,jj  ) )
-
-               ! Sea Surface Height at v-point after
-               zsshva(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )     &
-                  &          * ( e1t(ji  ,jj  ) * e2t(ji  ,jj  ) * ssha(ji  ,jj  )    &
-                  &            + e1t(ji  ,jj+1) * e2t(ji  ,jj+1) * ssha(ji  ,jj+1) )
-
-            END DO
-         END DO
-
-         ! Boundaries conditions
-         CALL lbc_lnk( zsshua, 'U', 1. )   ;    CALL lbc_lnk( zsshva, 'V', 1. )
-         CALL lbc_lnk( zsshub, 'U', 1. )   ;    CALL lbc_lnk( zsshvb, 'V', 1. )
-
-         ! Scale factors at before and after time step
-         ! -------------------------------------------
-         CALL dom_vvl_sf( zsshub, 'U', zfse3ub ) ;    CALL dom_vvl_sf( zsshua, 'U', zfse3ua ) 
-         CALL dom_vvl_sf( zsshvb, 'V', zfse3vb ) ;    CALL dom_vvl_sf( zsshva, 'V', zfse3va ) 
-
-         ! Asselin filtered scale factor at now time step
-         ! ----------------------------------------------
-         IF( (neuler == 0 .AND. kt == nit000) .OR. lk_dynspg_ts ) THEN
-            CALL dom_vvl_sf_ini( 'U', zfse3un ) ;   CALL dom_vvl_sf_ini( 'V', zfse3vn ) 
-         ELSE
-            zsshun(:,:) = atfp * ( zsshub(:,:) + zsshua(:,:) ) + atfp1 * sshu(:,:)
-            zsshvn(:,:) = atfp * ( zsshvb(:,:) + zsshva(:,:) ) + atfp1 * sshv(:,:)
-            CALL dom_vvl_sf( zsshun, 'U', zfse3un ) ;   CALL dom_vvl_sf( zsshvn, 'V', zfse3vn ) 
-         ENDIF
-
-         ! Thickness weighting
-         ! -------------------
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ua(ji,jj,jk) = ua(ji,jj,jk) * fse3u(ji,jj,jk)
-                  va(ji,jj,jk) = va(ji,jj,jk) * fse3v(ji,jj,jk)
-
-                  un(ji,jj,jk) = un(ji,jj,jk) * fse3u(ji,jj,jk)
-                  vn(ji,jj,jk) = vn(ji,jj,jk) * fse3v(ji,jj,jk)
-
-                  ub(ji,jj,jk) = ub(ji,jj,jk) * zfse3ub(ji,jj,jk)
-                  vb(ji,jj,jk) = vb(ji,jj,jk) * zfse3vb(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-
-      ENDIF
 
       ! Lateral boundary conditions on ( ua, va )
@@ -175,42 +99,28 @@
       CALL lbc_lnk( va, 'V', -1. )
 
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         ! Next velocity
-         ! -------------
+      ! Next velocity
+      ! -------------
 #if defined key_dynspg_flt
-         ! Leap-frog time stepping already done in dynspg.F routine
+      ! Leap-frog time stepping already done in dynspg_flt.F routine
 #else
-         DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
-            DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
-               ! Leap-frog time stepping
-               ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk)
-               va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk)
-            END DO
-         END DO
-
-         IF( lk_vvl ) THEN
-            ! Unweight velocities prior to updating open boundaries.
-
-            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
-               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
-                  ua(ji,jj,jk) = ua(ji,jj,jk) / fse3u(ji,jj,jk)
-                  va(ji,jj,jk) = va(ji,jj,jk) / fse3v(ji,jj,jk)
-
-                  un(ji,jj,jk) = un(ji,jj,jk) / fse3u(ji,jj,jk)
-                  vn(ji,jj,jk) = vn(ji,jj,jk) / fse3v(ji,jj,jk)
-
-                  ub(ji,jj,jk) = ub(ji,jj,jk) / zfse3ub(ji,jj,jk)
-                  vb(ji,jj,jk) = vb(ji,jj,jk) / zfse3vb(ji,jj,jk)
-               END DO
-            END DO
-
-         ENDIF
+      IF( lk_vvl ) THEN                                  ! Varying levels
+         !RB_vvl scale factors are wrong at this point
+         DO jk = 1, jpkm1
+            ua(ji,jj,jk) = (        ub(:,:,jk) * fse3u(:,:,jk)     &
+               &           + z2dt * ua(:,:,jk) * fse3u(:,:,jk) )   &
+               &         / fse3u(:,:,jk) * umask(:,:,jk)
+            va(ji,jj,jk) = (        vb(:,:,jk) * fse3v(:,:,jk)     &
+               &           + z2dt * va(:,:,jk) * fse3v(:,:,jk) )   &
+               &         / fse3v(:,:,jk) * vmask(:,:,jk)
+         END DO
+      ELSE
+         DO jk = 1, jpkm1
+                  ! Leap-frog time stepping
+            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
+            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
+         END DO
+      ENDIF
 
 # if defined key_obc
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
       ! Update (ua,va) along open boundaries (only in the rigid-lid case)
       CALL obc_dyn( kt )
@@ -235,11 +145,5 @@
       ENDIF
 
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
 # elif defined key_bdy 
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
       ! Update (ua,va) along open boundaries (for exp or ts options).
       IF ( lk_dynspg_exp .or. lk_dynspg_ts ) THEN
@@ -279,82 +183,32 @@
       ENDIF
 
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
 # endif
 # if defined key_agrif
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-      ! Update (ua,va) along open boundaries (only in the rigid-lid case)
       CALL Agrif_dyn( kt )
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
 # endif
 #endif
 
-         ! Time filter and swap of dynamics arrays
-         ! ------------------------------------------
-         IF( neuler == 0 .AND. kt == nit000 ) THEN
-            IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN      ! Varying levels
-               ! caution: don't use (:,:) for this loop
-               ! it causes optimization problems on NEC in auto-tasking
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zsshun1 = umask(ji,jj,jk) / fse3u(ji,jj,jk)
-                     zsshvn1 = vmask(ji,jj,jk) / fse3v(ji,jj,jk)
-                     ub(ji,jj,jk) = un(ji,jj,jk) * zsshun1 * umask(ji,jj,jk)
-                     vb(ji,jj,jk) = vn(ji,jj,jk) * zsshvn1 * vmask(ji,jj,jk)
-                     zsshun1 = umask(ji,jj,jk) / zfse3ua(ji,jj,jk)
-                     zsshvn1 = vmask(ji,jj,jk) / zfse3va(ji,jj,jk)
-                     un(ji,jj,jk) = ua(ji,jj,jk) * zsshun1 * umask(ji,jj,jk)
-                     vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn1 * vmask(ji,jj,jk)
-                  END DO
-               END DO
-            ELSE                                               ! Fixed levels
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     ! Euler (forward) time stepping
-                     ub(ji,jj,jk) = un(ji,jj,jk)
-                     vb(ji,jj,jk) = vn(ji,jj,jk)
-                     un(ji,jj,jk) = ua(ji,jj,jk)
-                     vn(ji,jj,jk) = va(ji,jj,jk)
-                  END DO
-               END DO
-            ENDIF
-         ELSE
-            IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN      ! Varying levels
-               ! caution: don't use (:,:) for this loop
-               ! it causes optimization problems on NEC in auto-tasking
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zsshun1 = umask(ji,jj,jk) / zfse3un(ji,jj,jk)
-                     zsshvn1 = vmask(ji,jj,jk) / zfse3vn(ji,jj,jk)
-                     ub(ji,jj,jk) = ( atfp  * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) &
-                       &            + atfp1 * un(ji,jj,jk) ) * zsshun1
-                     vb(ji,jj,jk) = ( atfp  * ( vb(ji,jj,jk) + va(ji,jj,jk) ) &
-                       &            + atfp1 * vn(ji,jj,jk) ) * zsshvn1
-                     zsshun1 = umask(ji,jj,jk) / zfse3ua(ji,jj,jk)
-                     zsshvn1 = vmask(ji,jj,jk) / zfse3va(ji,jj,jk)
-                     un(ji,jj,jk) = ua(ji,jj,jk) * zsshun1
-                     vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn1
-                  END DO
-               END DO
-            ELSE                                               ! Fixed levels
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     ! Leap-frog time stepping
-                     ub(ji,jj,jk) = atfp  * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk)
-                     vb(ji,jj,jk) = atfp  * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk)
-                     un(ji,jj,jk) = ua(ji,jj,jk)
-                     vn(ji,jj,jk) = va(ji,jj,jk)
-                  END DO
-               END DO
-            ENDIF
+      ! Time filter and swap of dynamics arrays
+      ! ------------------------------------------
+      IF( neuler == 0 .AND. kt == nit000 ) THEN
+         DO jk = 1, jpkm1                                 
+            ub(:,:,jk) = un(:,:,jk) 
+            vb(:,:,jk) = vn(:,:,jk)
+            un(:,:,jk) = ua(:,:,jk)
+            vn(:,:,jk) = va(:,:,jk)
+         END DO
+      ELSE
+         IF( lk_vvl ) THEN                                ! Varying levels
+            ! Not done
+         ELSE                                             ! Fixed levels
+!RB_vvl : should be done as in tranxt ?
+            DO jk = 1, jpkm1                              ! filter applied on velocities
+               ub(:,:,jk) = atfp * ( ub(:,:,jk) + ua(:,:,jk) ) + atfp1 * un(:,:,jk)
+               vb(:,:,jk) = atfp * ( vb(:,:,jk) + va(:,:,jk) ) + atfp1 * vn(:,:,jk)
+               un(:,:,jk) = ua(:,:,jk)
+               vn(:,:,jk) = va(:,:,jk)
+            END DO
          ENDIF
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
+      ENDIF
 
       IF(ln_ctl)   THEN
Index: branches/dev_004_VVL/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 1360)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 1361)
@@ -39,4 +39,6 @@
    PUBLIC   tra_nxt    ! routine called by step.F90
 
+   REAL(wp), DIMENSION(jpk) ::   r2dt_t   ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler)
+
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -67,15 +69,6 @@
       !!             at the AGRIF zoom     boundaries (lk_agrif=T)
       !!
-      !!              - Apply the Asselin time filter on now fields,
-      !!             save in (ta,sa) an average over the three time levels 
-      !!             which will be used to compute rdn and thus the semi-implicit
-      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T), and
-      !!             swap tracer fields to prepare the next time_step.
-      !!                This can be summurized for tempearture as:
-      !!                    zt = (ta+2tn+tb)/4       ln_dynhpg_imp = T
-      !!                    zt = 0                   otherwise
-      !!                    tb = tn + atfp*[ tb - 2 tn + ta ]
-      !!                    tn = ta 
-      !!                    ta = zt        (NB: reset to 0 after eos_bn2 call)
+      !!              - Update lateral boundary conditions on AGRIF children
+      !!             domains (lk_agrif=T)
       !!
       !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
@@ -87,6 +80,6 @@
       INTEGER, INTENT(in) ::   kt    ! ocean time-step index
       !!
-      INTEGER  ::   ji, jj, jk    ! dummy loop indices
-      REAL(wp) ::   zt, zs, zfact ! temporary scalars
+      INTEGER  ::   jk    ! dummy loop indices
+      REAL(wp) ::   zfact ! temporary scalars
       !!----------------------------------------------------------------------
 
@@ -111,4 +104,9 @@
       CALL Agrif_tra                   ! AGRIF zoom boundaries
 #endif
+ 
+      ! set time step size (Euler/Leapfrog)
+      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt_t(:) =     rdttra(:)      ! at nit000             (Euler)
+      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt_t(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
+      ENDIF
 
       ! trends computation initialisation
@@ -118,44 +116,7 @@
       ENDIF
 
-      ! Asselin time filter and swap of arrays
-      !                                             
-      IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler 1st time step : swap only
-         DO jk = 1, jpkm1
-            tb(:,:,jk) = tn(:,:,jk)                       ! ta, sa remain at their values which
-            sb(:,:,jk) = sn(:,:,jk)                       ! correspond to tn, sn after the sawp
-            tn(:,:,jk) = ta(:,:,jk)                     
-            sn(:,:,jk) = sa(:,:,jk)
-         END DO
-         !                                           
-      ELSE                                           ! Leap-Frog : filter + swap
-         !                                            
-         IF( ln_dynhpg_imp ) THEN                         ! semi-implicite hpg case
-            DO jk = 1, jpkm1                              ! (save the averaged of the 3 time steps 
-               DO jj = 1, jpj                             !                   in the after fields)
-                  DO ji = 1, jpi
-                     zt = ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) ) * 0.25
-                     zs = ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ) * 0.25
-                     tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk)
-                     sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk)
-                     tn(ji,jj,jk) = ta(ji,jj,jk)
-                     sn(ji,jj,jk) = sa(ji,jj,jk)
-                     ta(ji,jj,jk) = zt
-                     sa(ji,jj,jk) = zs
-                  END DO
-               END DO
-            END DO
-         ELSE                                            ! explicit hpg case
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk)
-                     sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk)
-                     tn(ji,jj,jk) = ta(ji,jj,jk)
-                     sn(ji,jj,jk) = sa(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-         ENDIF
-         !
+      ! Leap-Frog + Asselin filter time stepping
+      IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl( kt )      ! variable volume level (vvl)
+      ELSE                  ;   CALL tra_nxt_fix( kt )      ! fixed    volume level
       ENDIF
 
@@ -168,5 +129,5 @@
       IF( l_trdtra ) THEN     
          DO jk = 1, jpkm1
-            zfact = 1.e0 / ( 2.*rdttra(jk) )             ! NB: euler case, (tb filtered - tb)=0 so 2dt always OK
+            zfact = 1.e0 / r2dt_t(jk)             
             ztrdt(:,:,jk) = ( tb(:,:,jk) - ztrdt(:,:,jk) ) * zfact
             ztrds(:,:,jk) = ( sb(:,:,jk) - ztrds(:,:,jk) ) * zfact
@@ -180,4 +141,147 @@
    END SUBROUTINE tra_nxt
 
+
+   SUBROUTINE tra_nxt_fix( kt )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE tra_nxt_fix  ***
+      !!
+      !! ** Purpose :   fixed volume: apply the Asselin time filter and 
+      !!                swap the tracer fields.
+      !! 
+      !! ** Method  : - Apply a Asselin time filter on now fields.
+      !!              - save in (ta,sa) an average over the three time levels 
+      !!             which will be used to compute rdn and thus the semi-implicit
+      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
+      !!              - swap tracer fields to prepare the next time_step.
+      !!                This can be summurized for tempearture as:
+      !!             ztm = (ta+2tn+tb)/4       ln_dynhpg_imp = T
+      !!             ztm = 0                   otherwise
+      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
+      !!             tn  = ta 
+      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
+      !!
+      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
+      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   ztm, ztf     ! temporary scalars
+      REAL(wp) ::   zsm, zsf     !    -         -
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      ENDIF
+      !
+      !                                              ! ----------------------- !
+      IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case !
+         !                                           ! ----------------------- !
+         !
+         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! case of Euler time-stepping at first time-step
+            DO jk = 1, jpkm1
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztm = 0.25 * ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) )      ! mean t
+                     zsm = 0.25 * ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) )
+                     tb(ji,jj,jk) = tn(ji,jj,jk)                                           ! tb <-- tn
+                     sb(ji,jj,jk) = sn(ji,jj,jk)
+                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tb <-- tn
+                     sn(ji,jj,jk) = sa(ji,jj,jk)
+                     ta(ji,jj,jk) = ztm                                                    ! ta <-- mean t
+                     sa(ji,jj,jk) = zsm
+                  END DO
+               END DO
+            END DO
+         ELSE                                             ! general case (Leapfrog + Asselin filter
+            DO jk = 1, jpkm1
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztm = 0.25 * ( ta(ji,jj,jk) + 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )       ! mean t
+                     zsm = 0.25 * ( sa(ji,jj,jk) + 2.* sn(ji,jj,jk) + sb(ji,jj,jk) )
+                     ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )       ! Asselin filter on t 
+                     zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) )
+                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                     ! tb <-- filtered tn 
+                     sb(ji,jj,jk) = sn(ji,jj,jk) + zsf
+                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta
+                     sn(ji,jj,jk) = sa(ji,jj,jk)
+                     ta(ji,jj,jk) = ztm                                                    ! ta <-- mean t
+                     sa(ji,jj,jk) = zsm
+                  END DO
+               END DO
+            END DO
+         ENDIF
+         !                                           ! ----------------------- !
+      ELSE                                           !    explicit hpg case    !
+         !                                           ! ----------------------- !
+         !
+         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! case of Euler time-stepping at first time-step
+            DO jk = 1, jpkm1
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     tb(ji,jj,jk) = tn(ji,jj,jk)                                           ! tb <-- tn
+                     sb(ji,jj,jk) = sn(ji,jj,jk)
+                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta
+                     sn(ji,jj,jk) = sa(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         ELSE                                             ! general case (Leapfrog + Asselin filter
+            DO jk = 1, jpkm1
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+!RBvvl for reproducibility
+!                     ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )       ! Asselin filter on t 
+!                     zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) )
+!                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                     ! tb <-- filtered tn 
+!                     sb(ji,jj,jk) = sn(ji,jj,jk) + zsf
+                     tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk)
+                     sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk)
+                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta
+                     sn(ji,jj,jk) = sa(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         ENDIF
+      ENDIF
+      !
+   END SUBROUTINE tra_nxt_fix
+
+
+   SUBROUTINE tra_nxt_vvl( kt )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE tra_nxt_vvl  ***
+      !!
+      !! ** Purpose :   Time varying volume: apply the Asselin time filter  
+      !!                and swap the tracer fields.
+      !! 
+      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
+      !!              - save in (ta,sa) a thickness weighted average over the three 
+      !!             time levels which will be used to compute rdn and thus the semi-
+      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
+      !!              - swap tracer fields to prepare the next time_step.
+      !!                This can be summurized for tempearture as:
+      !!             ztm = (e3t_a*ta+2*e3t_n*tn+e3t_b*tb)   ln_dynhpg_imp = T
+      !!                  /(e3t_a   +2*e3t_n   +e3t_b   )    
+      !!             ztm = 0                                otherwise
+      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
+      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
+      !!             tn  = ta 
+      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
+      !!
+      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
+      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      !!     
+      
+      ! Not yet ready
+      WRITE(*,*) 'tra_next_vvl : you should not be there'
+      STOP
+      !
+   END SUBROUTINE tra_nxt_vvl
+
    !!======================================================================
 END MODULE tranxt
