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
