Index: /branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 4735)
+++ /branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 4736)
@@ -684,4 +684,5 @@
    ln_traadv_qck    =  .false.  !  QUICKEST scheme
    ln_traadv_msc_ups=  .false.  !  use upstream scheme within muscl
+   ln_traadv_tvd_zts=  .false.  !  TVD scheme with sub-timestepping of vertical tracer advection
 /
 !-----------------------------------------------------------------------
@@ -758,4 +759,5 @@
    ln_dynadv_cen2= .false. !  flux form - 2nd order centered scheme
    ln_dynadv_ubs = .false. !  flux form - 3rd order UBS      scheme
+   ln_dynzad_zts = .false. !  Use (T) sub timestepping for vertical momentum advection
 /
 !-----------------------------------------------------------------------
Index: /branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
===================================================================
--- /branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90	(revision 4735)
+++ /branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90	(revision 4736)
@@ -30,4 +30,5 @@
    LOGICAL, PUBLIC ::   ln_dynadv_cen2  !: flux form - 2nd order centered scheme flag
    LOGICAL, PUBLIC ::   ln_dynadv_ubs   !: flux form - 3rd order UBS scheme flag
+   LOGICAL, PUBLIC ::   ln_dynzad_zts   !: vertical advection with sub-timestepping (requires vector form)
    
    INTEGER ::   nadv   ! choice of the formulation and scheme for the advection
@@ -64,7 +65,10 @@
                       CALL dyn_keg     ( kt )    ! vector form : horizontal gradient of kinetic energy
                       CALL dyn_zad     ( kt )    ! vector form : vertical advection
-      CASE ( 1 ) 
+      CASE ( 1 )     
+                      CALL dyn_keg     ( kt )    ! vector form : horizontal gradient of kinetic energy
+                      CALL dyn_zad_zts ( kt )    ! vector form : vertical advection with sub-timestepping
+      CASE ( 2 ) 
                       CALL dyn_adv_cen2( kt )    ! 2nd order centered scheme
-      CASE ( 2 )   
+      CASE ( 3 )   
                       CALL dyn_adv_ubs ( kt )    ! 3rd order UBS      scheme
       !
@@ -91,5 +95,5 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       !!
-      NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs
+      NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs, ln_dynzad_zts
       !!----------------------------------------------------------------------
 
@@ -111,4 +115,5 @@
          WRITE(numout,*) '          2nd order centred advection scheme ln_dynadv_cen2 = ', ln_dynadv_cen2
          WRITE(numout,*) '          3rd order UBS advection scheme     ln_dynadv_ubs  = ', ln_dynadv_ubs
+         WRITE(numout,*) '      Sub timestepping of vertical advection ln_dynzad_zts  = ', ln_dynzad_zts
       ENDIF
 
@@ -120,9 +125,12 @@
 
       IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namdyn_adv' )
+      IF( ln_dynzad_zts .AND. .NOT. ln_dynadv_vec )   &
+          CALL ctl_stop( 'Sub timestepping of vertical advection requires vector form; set ln_dynadv_vec = .TRUE.' )
 
       !                               ! Set nadv
       IF( ln_dynadv_vec  )   nadv =  0 
-      IF( ln_dynadv_cen2 )   nadv =  1
-      IF( ln_dynadv_ubs  )   nadv =  2
+      IF( ln_dynzad_zts  )   nadv =  1
+      IF( ln_dynadv_cen2 )   nadv =  2
+      IF( ln_dynadv_ubs  )   nadv =  3
       IF( lk_esopa       )   nadv = -1
 
@@ -130,6 +138,7 @@
          WRITE(numout,*)
          IF( nadv ==  0 )   WRITE(numout,*) '         vector form : keg + zad + vor is used'
-         IF( nadv ==  1 )   WRITE(numout,*) '         flux form   : 2nd order scheme is used'
-         IF( nadv ==  2 )   WRITE(numout,*) '         flux form   : UBS       scheme is used'
+         IF( nadv ==  1 )   WRITE(numout,*) '         vector form : keg + zad_zts + vor is used'
+         IF( nadv ==  2 )   WRITE(numout,*) '         flux form   : 2nd order scheme is used'
+         IF( nadv ==  3 )   WRITE(numout,*) '         flux form   : UBS       scheme is used'
          IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection formulation'
       ENDIF
Index: /branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
===================================================================
--- /branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 4735)
+++ /branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 4736)
@@ -27,5 +27,6 @@
    PRIVATE
    
-   PUBLIC   dyn_zad   ! routine called by step.F90
+   PUBLIC   dyn_zad       ! routine called by dynadv.F90
+   PUBLIC   dyn_zad_zts   ! routine called by dynadv.F90
 
    !! * Substitutions
@@ -83,5 +84,5 @@
          DO jj = 2, jpj                   ! vertical fluxes 
             DO ji = fs_2, jpi             ! vector opt.
-               zww(ji,jj) = 0.25 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
+               zww(ji,jj) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
             END DO
          END DO
@@ -95,8 +96,8 @@
       DO jj = 2, jpjm1              ! Surface and bottom values set to zero
          DO ji = fs_2, fs_jpim1           ! vector opt.
-            zwuw(ji,jj, 1 ) = 0.e0
-            zwvw(ji,jj, 1 ) = 0.e0
-            zwuw(ji,jj,jpk) = 0.e0
-            zwvw(ji,jj,jpk) = 0.e0
+            zwuw(ji,jj, 1 ) = 0._wp
+            zwvw(ji,jj, 1 ) = 0._wp
+            zwuw(ji,jj,jpk) = 0._wp
+            zwvw(ji,jj,jpk) = 0._wp
          END DO  
       END DO
@@ -132,4 +133,146 @@
    END SUBROUTINE dyn_zad
 
+   SUBROUTINE dyn_zad_zts ( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dynzad_zts  ***
+      !! 
+      !! ** Purpose :   Compute the now vertical momentum advection trend and 
+      !!      add it to the general trend of momentum equation. This version
+      !!      uses sub-timesteps for improved numerical stability with small
+      !!      vertical grid sizes. This is especially relevant when using 
+      !!      embedded ice with thin surface boxes.
+      !!
+      !! ** Method  :   The now vertical advection of momentum is given by:
+      !!         w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ]
+      !!         w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ]
+      !!      Add this trend to the general trend (ua,va):
+      !!         (ua,va) = (ua,va) + w dz(u,v)
+      !!
+      !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends
+      !!              - Save the trends in (ztrdu,ztrdv) ('key_trddyn')
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
+      !
+      INTEGER  ::   ji, jj, jk, jl  ! dummy loop indices
+      INTEGER  ::   jnzts = 5       ! number of sub-timesteps for vertical advection
+      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps
+      REAL(wp) ::   zua, zva        ! temporary scalars
+      REAL(wp) ::   zr_rdt          ! temporary scalar
+      REAL(wp) ::   z2dtzts         ! length of Euler forward sub-timestep for vertical advection
+      REAL(wp) ::   zts             ! length of sub-timestep for vertical advection
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zwuw , zwvw, zww
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztrdu, ztrdv
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zus , zvs
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('dyn_zad_zts')
+      !
+      CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 
+      CALL wrk_alloc( jpi,jpj,jpk,3, zus, zvs ) 
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp)WRITE(numout,*)
+         IF(lwp)WRITE(numout,*) 'dyn_zad_zts : arakawa advection scheme with sub-timesteps'
+      ENDIF
+
+      IF( l_trddyn )   THEN         ! Save ua and va trends
+         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
+         ztrdu(:,:,:) = ua(:,:,:) 
+         ztrdv(:,:,:) = va(:,:,:) 
+      ENDIF
+      
+      IF( neuler == 0 .AND. kt == nit000 ) THEN
+          z2dtzts =         rdt / REAL( jnzts, wp )   ! = rdt (restart with Euler time stepping)
+      ELSE
+          z2dtzts = 2._wp * rdt / REAL( jnzts, wp )   ! = 2 rdt (leapfrog)
+      ENDIF
+      
+      DO jk = 2, jpkm1                    ! Calculate and store vertical fluxes
+         DO jj = 2, jpj                   
+            DO ji = fs_2, jpi             ! vector opt.
+               zww(ji,jj,jk) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      DO jj = 2, jpjm1                    ! Surface and bottom advective fluxes set to zero
+         DO ji = fs_2, fs_jpim1           ! vector opt.
+            zwuw(ji,jj, 1 ) = 0._wp
+            zwvw(ji,jj, 1 ) = 0._wp
+            zwuw(ji,jj,jpk) = 0._wp
+            zwvw(ji,jj,jpk) = 0._wp
+         END DO  
+      END DO
+
+! Start with before values and use sub timestepping to reach after values
+
+      zus(:,:,:,1) = ub(:,:,:)
+      zvs(:,:,:,1) = vb(:,:,:)
+
+      DO jl = 1, jnzts                   ! Start of sub timestepping loop
+
+         IF( jl == 1 ) THEN              ! Euler forward to kick things off
+           jtb = 1   ;   jtn = 1   ;   jta = 2
+           zts = z2dtzts
+         ELSEIF( jl == 2 ) THEN          ! First leapfrog step
+           jtb = 1   ;   jtn = 2   ;   jta = 3
+           zts = 2._wp * z2dtzts
+         ELSE                            ! Shuffle pointers for subsequent leapfrog steps
+           jtb = MOD(jtb,3) + 1
+           jtn = MOD(jtn,3) + 1
+           jta = MOD(jta,3) + 1
+         ENDIF
+
+         DO jk = 2, jpkm1           ! Vertical momentum advection at level w and u- and v- vertical
+            DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
+               DO ji = fs_2, fs_jpim1        ! vector opt.
+                  zwuw(ji,jj,jk) = ( zww(ji+1,jj  ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) )
+                  zwvw(ji,jj,jk) = ( zww(ji  ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) )
+               END DO  
+            END DO   
+         END DO
+         DO jk = 1, jpkm1           ! Vertical momentum advection at u- and v-points
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1       ! vector opt.
+                  !                         ! vertical momentum advective trends
+                  zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
+                  zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
+                  zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts
+                  zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts
+               END DO  
+            END DO  
+         END DO
+
+      END DO      ! End of sub timestepping loop
+
+      zr_rdt = 1._wp / ( REAL( jnzts, wp ) * z2dtzts )
+      DO jk = 1, jpkm1              ! Recover trends over the outer timestep
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1       ! vector opt.
+               !                         ! vertical momentum advective trends
+               !                         ! add the trends to the general momentum trends
+               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zus(ji,jj,jk,jta) - ub(ji,jj,jk)) * zr_rdt
+               va(ji,jj,jk) = va(ji,jj,jk) + ( zvs(ji,jj,jk,jta) - vb(ji,jj,jk)) * zr_rdt
+            END DO  
+         END DO  
+      END DO
+
+      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
+         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
+         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
+         CALL trd_mod(ztrdu, ztrdv, jpdyn_trd_zad, 'DYN', kt)
+         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
+      ENDIF
+      !                             ! Control print
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
+         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      !
+      CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 
+      CALL wrk_dealloc( jpi,jpj,jpk,3, zus, zvs ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('dyn_zad_zts')
+      !
+   END SUBROUTINE dyn_zad_zts
+
    !!======================================================================
 END MODULE dynzad
Index: /branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- /branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 4735)
+++ /branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 4736)
@@ -43,4 +43,5 @@
    LOGICAL ::   ln_traadv_cen2     ! 2nd order centered scheme flag
    LOGICAL ::   ln_traadv_tvd      ! TVD scheme flag
+   LOGICAL ::   ln_traadv_tvd_zts  ! TVD scheme flag with vertical sub time-stepping
    LOGICAL ::   ln_traadv_muscl    ! MUSCL scheme flag
    LOGICAL ::   ln_traadv_muscl2   ! MUSCL2 scheme flag
@@ -120,4 +121,5 @@
       CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  UBS 
       CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  QUICKEST 
+      CASE ( 7 )   ;   CALL tra_adv_tvd_zts( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD ZTS
       !
       CASE (-1 )                                      !==  esopa: test all possibility with control print  ==!
@@ -166,5 +168,5 @@
          &                 ln_traadv_muscl, ln_traadv_muscl2,  &
          &                 ln_traadv_ubs  , ln_traadv_qck,     &
-         &                 ln_traadv_msc_ups
+         &                 ln_traadv_msc_ups, ln_traadv_tvd_zts
       !!----------------------------------------------------------------------
 
@@ -190,4 +192,5 @@
          WRITE(numout,*) '      QUICKEST advection scheme      ln_traadv_qck     = ', ln_traadv_qck
          WRITE(numout,*) '      upstream scheme within muscl   ln_traadv_msc_ups = ', ln_traadv_msc_ups
+         WRITE(numout,*) '      TVD advection scheme with zts  ln_traadv_tvd_zts = ', ln_traadv_tvd_zts
       ENDIF
 
@@ -199,4 +202,5 @@
       IF( ln_traadv_ubs    )   ioptio = ioptio + 1
       IF( ln_traadv_qck    )   ioptio = ioptio + 1
+      IF( ln_traadv_tvd_zts)   ioptio = ioptio + 1
       IF( lk_esopa         )   ioptio =          1
 
@@ -210,4 +214,5 @@
       IF( ln_traadv_ubs    )   nadv =  5
       IF( ln_traadv_qck    )   nadv =  6
+      IF( ln_traadv_tvd_zts)   nadv =  7
       IF( lk_esopa         )   nadv = -1
 
@@ -220,4 +225,5 @@
          IF( nadv ==  5 )   WRITE(numout,*) '         UBS       scheme is used'
          IF( nadv ==  6 )   WRITE(numout,*) '         QUICKEST  scheme is used'
+         IF( nadv ==  7 )   WRITE(numout,*) '         TVD ZTS   scheme is used'
          IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection scheme'
       ENDIF
Index: /branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
===================================================================
--- /branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 4735)
+++ /branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 4736)
@@ -37,5 +37,6 @@
    PRIVATE
 
-   PUBLIC   tra_adv_tvd    ! routine called by step.F90
+   PUBLIC   tra_adv_tvd        ! routine called by traadv.F90
+   PUBLIC   tra_adv_tvd_zts    ! routine called by traadv.F90
 
    LOGICAL ::   l_trd   ! flag to compute trends
@@ -247,4 +248,279 @@
    END SUBROUTINE tra_adv_tvd
 
+   SUBROUTINE tra_adv_tvd_zts ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
+      &                                       ptb, ptn, pta, kjpt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_adv_tvd_zts  ***
+      !! 
+      !! **  Purpose :   Compute the now trend due to total advection of 
+      !!       tracers and add it to the general trend of tracer equations
+      !!
+      !! **  Method  :   TVD ZTS scheme, i.e. 2nd order centered scheme with
+      !!       corrected flux (monotonic correction). This version use sub-
+      !!       timestepping for the vertical advection which increases stability
+      !!       when vertical metrics are small.
+      !!       note: - this advection scheme needs a leap-frog time scheme
+      !!
+      !! ** Action : - update (pta) with the now advective tracer trends
+      !!             - save the trends 
+      !!----------------------------------------------------------------------
+      USE oce     , ONLY:   zwx => ua        , zwy => va          ! (ua,va) used as workspace
+      !
+      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
+      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
+      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
+      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend 
+      !
+      REAL(wp), DIMENSION( jpk )                           ::   zts             ! length of sub-timestep for vertical advection
+      REAL(wp), DIMENSION( jpk )                           ::   zr_p2dt         ! reciprocal of tracer timestep
+      INTEGER  ::   ji, jj, jk, jl, jn       ! dummy loop indices  
+      INTEGER  ::   jnzts = 5       ! number of sub-timesteps for vertical advection
+      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps
+      INTEGER  ::   jtaken          ! toggle for collecting appropriate fluxes from sub timesteps
+      REAL(wp) ::   z_rzts          ! Fractional length of Euler forward sub-timestep for vertical advection
+      REAL(wp) ::   z2dtt, zbtr, ztra        ! local scalar
+      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !   -      -
+      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk   !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zwx_sav , zwy_sav
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwz, zhdiv, zwz_sav, zwzts
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrs
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('tra_adv_tvd_zts')
+      !
+      CALL wrk_alloc( jpi, jpj, zwx_sav, zwy_sav )
+      CALL wrk_alloc( jpi, jpj, jpk, zwi, zwz , zhdiv, zwz_sav, zwzts )
+      CALL wrk_alloc( jpi, jpj, jpk, 3, ztrs )
+      !
+      IF( kt == kit000 )  THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_adv_tvd_zts : TVD ZTS advection scheme on ', cdtype
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      ENDIF
+      !
+      l_trd = .FALSE.
+      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
+      !
+      IF( l_trd )  THEN
+         CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
+         ztrdx(:,:,:) = 0._wp  ;    ztrdy(:,:,:) = 0._wp  ;   ztrdz(:,:,:) = 0._wp
+      ENDIF
+      !
+      zwi(:,:,:) = 0._wp
+      z_rzts = 1._wp / REAL( jnzts, wp )
+      zr_p2dt(:) = 1._wp / p2dt(:)
+      !
+      !                                                          ! ===========
+      DO jn = 1, kjpt                                            ! tracer loop
+         !                                                       ! ===========
+         ! 1. Bottom value : flux set to zero
+         ! ----------------------------------
+         zwx(:,:,jpk) = 0._wp   ;    zwz(:,:,jpk) = 0._wp
+         zwy(:,:,jpk) = 0._wp   ;    zwi(:,:,jpk) = 0._wp
+
+         ! 2. upstream advection with initial mass fluxes & intermediate update
+         ! --------------------------------------------------------------------
+         ! upstream tracer flux in the i and j direction
+         DO jk = 1, jpkm1
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  ! upstream scheme
+                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
+                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
+                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
+                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
+                  zwx(ji,jj,jk) = 0.5_wp * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn) )
+                  zwy(ji,jj,jk) = 0.5_wp * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn) )
+               END DO
+            END DO
+         END DO
+
+         ! upstream tracer flux in the k direction
+         ! Surface value
+         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0._wp                        ! volume variable
+         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface 
+         ENDIF
+         ! Interior value
+         DO jk = 2, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
+                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
+                  zwz(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) )
+               END DO
+            END DO
+         END DO
+
+         ! total advective trend
+         DO jk = 1, jpkm1
+            z2dtt = p2dt(jk)
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                  ! total intermediate advective trends
+                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
+                     &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) )
+                  ! update and guess with monotonic sheme
+                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra
+                  zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !                             ! Lateral boundary conditions on zwi  (unchanged sign)
+         CALL lbc_lnk( zwi, 'T', 1. )  
+
+         !                                 ! trend diagnostics (contribution of upstream fluxes)
+         IF( l_trd )  THEN 
+            ! store intermediate advective trends
+            ztrdx(:,:,:) = zwx(:,:,:)   ;    ztrdy(:,:,:) = zwy(:,:,:)  ;   ztrdz(:,:,:) = zwz(:,:,:)
+         END IF
+         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
+         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN  
+           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
+           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
+         ENDIF
+
+         ! 3. antidiffusive flux : high order minus low order
+         ! --------------------------------------------------
+         ! antidiffusive flux on i and j
+
+
+         DO jk = 1, jpkm1
+
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwx_sav(ji,jj) = zwx(ji,jj,jk)
+                  zwy_sav(ji,jj) = zwy(ji,jj,jk)
+
+                  zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) )
+                  zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) )
+               END DO
+            END DO
+
+            DO jj = 2, jpjm1         ! partial horizontal divergence
+               DO ji = fs_2, fs_jpim1
+                  zhdiv(ji,jj,jk) = (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
+                     &               + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
+               END DO
+            END DO
+
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwx(ji,jj,jk) = zwx(ji,jj,jk)  - zwx_sav(ji,jj)
+                  zwy(ji,jj,jk) = zwy(ji,jj,jk)  - zwy_sav(ji,jj)
+               END DO
+            END DO
+         END DO
+      
+         ! antidiffusive flux on k
+         zwz(:,:,1) = 0._wp        ! Surface value
+         zwz_sav(:,:,:) = zwz(:,:,:)
+         !
+         ztrs(:,:,:,1) = ptb(:,:,:,jn)
+         zwzts(:,:,:) = 0._wp
+
+         DO jl = 1, jnzts                   ! Start of sub timestepping loop
+
+            IF( jl == 1 ) THEN              ! Euler forward to kick things off
+              jtb = 1   ;   jtn = 1   ;   jta = 2
+              zts(:) = p2dt(:) * z_rzts
+              jtaken = MOD( jnzts + 1 , 2)  ! Toggle to collect every second flux
+                                            ! starting at jl =1 if jnzts is odd; 
+                                            ! starting at jl =2 otherwise
+            ELSEIF( jl == 2 ) THEN          ! First leapfrog step
+              jtb = 1   ;   jtn = 2   ;   jta = 3
+              zts(:) = 2._wp * p2dt(:) * z_rzts
+            ELSE                            ! Shuffle pointers for subsequent leapfrog steps
+              jtb = MOD(jtb,3) + 1
+              jtn = MOD(jtn,3) + 1
+              jta = MOD(jta,3) + 1
+            ENDIF
+            DO jk = 2, jpkm1          ! Interior value
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) )
+                     IF( jtaken == 0 ) zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk)*zts(jk)           ! Accumulate time-weighted vertcal flux
+                  END DO
+               END DO
+            END DO
+
+            jtaken = MOD( jtaken + 1 , 2 )
+
+            DO jk = 2, jpkm1          ! Interior value
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                     ! total advective trends
+                     ztra = - zbtr * (  zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) )
+                     ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb) + zts(jk) * ztra
+                  END DO
+               END DO
+            END DO
+
+         END DO
+
+         DO jk = 2, jpkm1          ! Anti-diffusive vertical flux using average flux from the sub-timestepping
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  zwz(ji,jj,jk) = zwzts(ji,jj,jk) * zr_p2dt(jk) - zwz_sav(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions
+         CALL lbc_lnk( zwz, 'W',  1. )
+
+         ! 4. monotonicity algorithm
+         ! -------------------------
+         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
+
+
+         ! 5. final trend with corrected fluxes
+         ! ------------------------------------
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.  
+                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                  ! total advective trends
+                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
+                     &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) )
+                  ! add them to the general tracer trends
+                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
+               END DO
+            END DO
+         END DO
+
+         !                                 ! trend diagnostics (contribution of upstream fluxes)
+         IF( l_trd )  THEN 
+            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed
+            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
+            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed
+            
+            CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, ztrdx, pun, ptn(:,:,:,jn) )   
+            CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, ztrdy, pvn, ptn(:,:,:,jn) )  
+            CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 
+         END IF
+         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
+         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN  
+           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) ) + htr_adv(:)
+           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) ) + str_adv(:)
+         ENDIF
+         !
+      END DO
+      !
+                   CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz, zhdiv, zwz_sav, zwzts )
+                   CALL wrk_dealloc( jpi, jpj, jpk, 3, ztrs )
+                   CALL wrk_dealloc( jpi, jpj, zwx_sav, zwy_sav )
+      IF( l_trd )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_tvd_zts')
+      !
+   END SUBROUTINE tra_adv_tvd_zts
 
    SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
