Index: trunk/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 1146)
+++ trunk/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 1438)
@@ -14,12 +14,16 @@
    !!            2.0  !  2006-02  (L. Debreu, C. Mazauric) Agrif implementation
    !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf
+   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
    !!   tra_nxt       : time stepping on temperature and salinity
+   !!   tra_nxt_fix   : time stepping on temperature and salinity : fixed    volume case
+   !!   tra_nxt_vvl   : time stepping on temperature and salinity : variable volume case
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers variables
    USE dom_oce         ! ocean space and time domain variables 
    USE zdf_oce         ! ???
+   USE domvvl          ! variable volume
    USE dynspg_oce      ! surface pressure gradient variables
    USE trdmod_oce      ! ocean variables trends
@@ -39,4 +43,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 +73,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 +84,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 +108,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 +120,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
 
@@ -165,8 +130,8 @@
 #endif      
 
-      ! trends diagnostics :  Asselin filter trend : (tb filtered - tb)/2dt
-      IF( l_trdtra ) THEN     
+      ! trends computation
+      IF( l_trdtra ) THEN              ! trend of the Asselin filter (tb filtered - tb)/dt     
          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
@@ -174,4 +139,5 @@
          CALL trd_mod( ztrdt, ztrds, jptra_trd_atf, 'TRA', kt )
       END IF
+
       !                        ! control print
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=tn, clinfo1=' nxt  - Tn: ', mask1=tmask,   &
@@ -180,4 +146,241 @@
    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        ! Euler time-stepping at first time-step
+            DO jk = 1, jpkm1                              ! (only swap)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tb <-- tn
+                     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
+                     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        ! Euler time-stepping at first time-step
+            DO jk = 1, jpkm1
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     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
+                     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)
+                  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
+      !!     
+      INTEGER  ::   ji, jj, jk             ! dummy loop indices
+      REAL(wp) ::   ztm , ztc_f , ztf , ztca, ztcn, ztcb   ! temporary scalar
+      REAL(wp) ::   zsm , zsc_f , zsf , zsca, zscn, zscb   !    -         -
+      REAL(wp) ::   ze3mr, ze3fr                           !    -         -
+      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f         !    -         -
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      ENDIF
+
+      !                                              ! ----------------------- !
+      IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case !
+         !                                           ! ----------------------- !
+         !
+         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
+            DO jk = 1, jpkm1                              ! (only swap)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     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
+            DO jk = 1, jpkm1
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ze3t_b = fse3t_b(ji,jj,jk)
+                     ze3t_n = fse3t_n(ji,jj,jk)
+                     ze3t_a = fse3t_a(ji,jj,jk)
+                     !                                         ! tracer content at Before, now and after
+                     ztcb = tb(ji,jj,jk) *  ze3t_b   ;   zscb = sb(ji,jj,jk) * ze3t_b
+                     ztcn = tn(ji,jj,jk) *  ze3t_n   ;   zscn = sn(ji,jj,jk) * ze3t_n
+                     ztca = ta(ji,jj,jk) *  ze3t_a   ;   zsca = sa(ji,jj,jk) * ze3t_a
+                     !
+                     !                                         ! Asselin filter on thickness and tracer content
+                     ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b )
+                     ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   ) 
+                     zsc_f  = atfp * ( zsca   - 2.* zscn   + zscb   ) 
+                     !
+                     !                                         ! filtered tracer including the correction 
+                     ze3fr = 1.e0  / ( ze3t_n + ze3t_f )
+                     ztf   = ze3fr * ( ztcn   + ztc_f  )
+                     zsf   = ze3fr * ( zscn   + zsc_f  )
+                     !                                         ! mean thickness and tracer
+                     ze3mr = 1.e0  / ( ze3t_a + 2.* ze3t_n + ze3t_b )
+                     ztm   = ze3mr * ( ztca   + 2.* ztcn   + ztcb   )
+                     zsm   = ze3mr * ( zsca   + 2.* zscn   + zscb   )
+!!gm mean e3t have to be saved and used in dynhpg  or it can be recomputed in dynhpg !!
+!!gm                 e3t_m(ji,jj,jk) = 0.25 / ze3mr
+                     !                                         ! swap of arrays
+                     tb(ji,jj,jk) = ztf                             ! tb <-- tn + filter
+                     sb(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                              ! No filter nor thickness weighting computation required
+               DO jj = 1, jpj                             ! ONLY swap
+                  DO ji = 1, jpi
+                     tn(ji,jj,jk) = ta(ji,jj,jk)                                 ! tn <-- ta
+                     sn(ji,jj,jk) = sa(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            !                                             ! general case (Leapfrog + Asselin filter)
+         ELSE                                             ! apply filter on thickness weighted tracer and swap
+            DO jk = 1, jpkm1
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ze3t_b = fse3t_b(ji,jj,jk)
+                     ze3t_n = fse3t_n(ji,jj,jk)
+                     ze3t_a = fse3t_a(ji,jj,jk)
+                     !                                         ! tracer content at Before, now and after
+                     ztcb = tb(ji,jj,jk) *  ze3t_b   ;   zscb = sb(ji,jj,jk) * ze3t_b
+                     ztcn = tn(ji,jj,jk) *  ze3t_n   ;   zscn = sn(ji,jj,jk) * ze3t_n
+                     ztca = ta(ji,jj,jk) *  ze3t_a   ;   zsca = sa(ji,jj,jk) * ze3t_a
+                     !
+                     !                                         ! Asselin filter on thickness and tracer content
+                     ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b )
+                     ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   ) 
+                     zsc_f  = atfp * ( zsca   - 2.* zscn   + zscb   ) 
+                     !
+                     !                                         ! filtered tracer including the correction 
+                     ze3fr = 1.e0 / ( ze3t_n + ze3t_f )
+                     ztf   =  ( ztcn  + ztc_f ) * ze3fr
+                     zsf   =  ( zscn  + zsc_f ) * ze3fr
+                     !                                         ! swap of arrays
+                     tb(ji,jj,jk) = ztf                             ! tb <-- tn filtered
+                     sb(ji,jj,jk) = zsf
+                     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_vvl
+
    !!======================================================================
 END MODULE tranxt
