Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/dom_oce.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/dom_oce.F90	(revision 10879)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/dom_oce.F90	(revision 10880)
@@ -38,9 +38,4 @@
    LOGICAL , PUBLIC ::   ln_iscpl       !: coupling with ice sheet
    LOGICAL , PUBLIC ::   ln_crs         !: Apply grid coarsening to dynamical model output or online passive tracers
-
-   !!----------------------------------------------------------------------
-   !! time level indices
-   !!----------------------------------------------------------------------
-   INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs
 
    !! Free surface parameters
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv.F90	(revision 10879)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv.F90	(revision 10880)
@@ -75,5 +75,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv( kt )
+   SUBROUTINE tra_adv( kt, Kbb, Kmm, pts, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv  ***
@@ -81,10 +81,12 @@
       !! ** Purpose :   compute the ocean tracer advection trend.
       !!
-      !! ** Method  : - Update (ua,va) with the advection term following nadv
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !! ** Method  : - Update (uu(:,:,:,Krhs),vv(:,:,:,Krhs)) with the advection term following nadv
+      !!----------------------------------------------------------------------
+      INTEGER                                  , INTENT(in)    :: kt             ! ocean time-step index
+      INTEGER                                  , INTENT(in)    :: Kbb, Kmm, Krhs ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers and RHS of tracer equation
       !
       INTEGER ::   jk   ! dummy loop index
-      REAL(wp), DIMENSION(jpi,jpj,jpk)        :: zun, zvn, zwn   ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk)        :: zuu, zvv, zww   ! 3D workspace
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds
       !!----------------------------------------------------------------------
@@ -98,47 +100,47 @@
       !
       !                                         !==  effective transport  ==!
-      zun(:,:,jpk) = 0._wp
-      zvn(:,:,jpk) = 0._wp
-      zwn(:,:,jpk) = 0._wp
+      zuu(:,:,jpk) = 0._wp
+      zvv(:,:,jpk) = 0._wp
+      zww(:,:,jpk) = 0._wp
       IF( ln_wave .AND. ln_sdw )  THEN
          DO jk = 1, jpkm1                                                       ! eulerian transport + Stokes Drift
-            zun(:,:,jk) = e2u  (:,:) * e3u_n(:,:,jk) * ( un(:,:,jk) + usd(:,:,jk) )
-            zvn(:,:,jk) = e1v  (:,:) * e3v_n(:,:,jk) * ( vn(:,:,jk) + vsd(:,:,jk) )
-            zwn(:,:,jk) = e1e2t(:,:)                 * ( wn(:,:,jk) + wsd(:,:,jk) )
+            zuu(:,:,jk) = e2u  (:,:) * e3u(:,:,jk,Kmm) * ( uu(:,:,jk,Kmm) + usd(:,:,jk) )
+            zvv(:,:,jk) = e1v  (:,:) * e3v(:,:,jk,Kmm) * ( vv(:,:,jk,Kmm) + vsd(:,:,jk) )
+            zww(:,:,jk) = e1e2t(:,:)                 * ( ww(:,:,jk) + wsd(:,:,jk) )
          END DO
       ELSE
          DO jk = 1, jpkm1
-            zun(:,:,jk) = e2u  (:,:) * e3u_n(:,:,jk) * un(:,:,jk)               ! eulerian transport only
-            zvn(:,:,jk) = e1v  (:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
-            zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk)
+            zuu(:,:,jk) = e2u  (:,:) * e3u(:,:,jk,Kmm) * uu(:,:,jk,Kmm)               ! eulerian transport only
+            zvv(:,:,jk) = e1v  (:,:) * e3v(:,:,jk,Kmm) * vv(:,:,jk,Kmm)
+            zww(:,:,jk) = e1e2t(:,:)                 * ww(:,:,jk)
          END DO
       ENDIF
       !
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN                                ! add z-tilde and/or vvl corrections
-         zun(:,:,:) = zun(:,:,:) + un_td(:,:,:)
-         zvn(:,:,:) = zvn(:,:,:) + vn_td(:,:,:)
-      ENDIF
-      !
-      zun(:,:,jpk) = 0._wp                                                      ! no transport trough the bottom
-      zvn(:,:,jpk) = 0._wp
-      zwn(:,:,jpk) = 0._wp
+         zuu(:,:,:) = zuu(:,:,:) + un_td(:,:,:)
+         zvv(:,:,:) = zvv(:,:,:) + vn_td(:,:,:)
+      ENDIF
+      !
+      zuu(:,:,jpk) = 0._wp                                                      ! no transport trough the bottom
+      zvv(:,:,jpk) = 0._wp
+      zww(:,:,jpk) = 0._wp
       !
       IF( ln_ldfeiv .AND. .NOT. ln_traldf_triad )   &
-         &              CALL ldf_eiv_trp( kt, nit000, zun, zvn, zwn, 'TRA' )   ! add the eiv transport (if necessary)
-      !
-      IF( ln_mle    )   CALL tra_mle_trp( kt, nit000, zun, zvn, zwn, 'TRA' )   ! add the mle transport (if necessary)
-      !
-      CALL iom_put( "uocetr_eff", zun )                                        ! output effective transport      
-      CALL iom_put( "vocetr_eff", zvn )
-      CALL iom_put( "wocetr_eff", zwn )
+         &              CALL ldf_eiv_trp( kt, nit000, zuu, zvv, zww, 'TRA' )   ! add the eiv transport (if necessary)
+      !
+      IF( ln_mle    )   CALL tra_mle_trp( kt, nit000, zuu, zvv, zww, 'TRA' )   ! add the mle transport (if necessary)
+      !
+      CALL iom_put( "uocetr_eff", zuu )                                        ! output effective transport      
+      CALL iom_put( "vocetr_eff", zvv )
+      CALL iom_put( "wocetr_eff", zww )
       !
 !!gm ???
-      IF( ln_diaptr )   CALL dia_ptr( zvn )                                    ! diagnose the effective MSF 
+      IF( ln_diaptr )   CALL dia_ptr( zvv )                                    ! diagnose the effective MSF 
 !!gm ???
       !
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
          ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) )
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
+         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
+         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
       ENDIF
       !
@@ -146,13 +148,13 @@
       !
       CASE ( np_CEN )                                 ! Centered scheme : 2nd / 4th order
-         CALL tra_adv_cen    ( kt, nit000, 'TRA',         zun, zvn, zwn     , tsn, tsa, jpts, nn_cen_h, nn_cen_v )
+         CALL tra_adv_cen    ( kt, nit000, 'TRA',         zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v )
       CASE ( np_FCT )                                 ! FCT scheme      : 2nd / 4th order
-         CALL tra_adv_fct    ( kt, nit000, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts, nn_fct_h, nn_fct_v )
+         CALL tra_adv_fct    ( kt, nit000, 'TRA', r2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v )
       CASE ( np_MUS )                                 ! MUSCL
-         CALL tra_adv_mus    ( kt, nit000, 'TRA', r2dt, zun, zvn, zwn, tsb,      tsa, jpts        , ln_mus_ups ) 
+         CALL tra_adv_mus    ( kt, nit000, 'TRA', r2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 
       CASE ( np_UBS )                                 ! UBS
-         CALL tra_adv_ubs    ( kt, nit000, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts        , nn_ubs_v   )
+         CALL tra_adv_ubs    ( kt, nit000, 'TRA', r2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v   )
       CASE ( np_QCK )                                 ! QUICKEST
-         CALL tra_adv_qck    ( kt, nit000, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts                     )
+         CALL tra_adv_qck    ( kt, nit000, 'TRA', r2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs )
       !
       END SELECT
@@ -160,6 +162,6 @@
       IF( l_trdtra )   THEN                      ! save the advective trends for further diagnostics
          DO jk = 1, jpkm1
-            ztrdt(:,:,jk) = tsa(:,:,jk,jp_tem) - ztrdt(:,:,jk)
-            ztrds(:,:,jk) = tsa(:,:,jk,jp_sal) - ztrds(:,:,jk)
+            ztrdt(:,:,jk) = pts(:,:,jk,jp_tem,Krhs) - ztrdt(:,:,jk)
+            ztrds(:,:,jk) = pts(:,:,jk,jp_sal,Krhs) - ztrds(:,:,jk)
          END DO
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_totad, ztrdt )
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_cen.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_cen.F90	(revision 10879)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_cen.F90	(revision 10880)
@@ -44,6 +44,6 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pun, pvn, pwn,     &
-      &                                        ptn, pta, kjpt, kn_cen_h, kn_cen_v ) 
+   SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pu_mm, pv_mm, pww,     &
+      &                    Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v ) 
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_cen  ***
@@ -59,17 +59,17 @@
       !!                = 4  ==>> 4th order COMPACT  scheme     -      -
       !!
-      !! ** Action : - update pta  with the now advective tracer trends
+      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
       !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
       !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
       !!----------------------------------------------------------------------
-      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
-      INTEGER                              , INTENT(in   ) ::   kn_cen_h        ! =2/4 (2nd or 4th order scheme)
-      INTEGER                              , INTENT(in   ) ::   kn_cen_v        ! =2/4 (2nd or 4th order scheme)
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn             ! now tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend 
+      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) ::   Kmm, Krhs       ! ocean time level indices
+      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
+      INTEGER                                  , INTENT(in   ) ::   kn_cen_h        ! =2/4 (2nd or 4th order scheme)
+      INTEGER                                  , INTENT(in   ) ::   kn_cen_v        ! =2/4 (2nd or 4th order scheme)
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pu_mm, pv_mm, pww   ! 3 ocean velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -106,6 +106,6 @@
                DO jj = 1, jpjm1
                   DO ji = 1, fs_jpim1   ! vector opt.
-                     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) )
+                     zwx(ji,jj,jk) = 0.5_wp * pu_mm(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm) )
+                     zwy(ji,jj,jk) = 0.5_wp * pv_mm(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) )
                   END DO
                END DO
@@ -118,6 +118,6 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
-                     ztv(ji,jj,jk) = ( ptn(ji  ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
+                     ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
+                     ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
                   END DO
                END DO
@@ -128,12 +128,12 @@
                DO jj = 2, jpjm1
                   DO ji = 1, fs_jpim1   ! vector opt.
-                     zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn)   ! C2 interpolation of T at u- & v-points (x2)
-                     zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn)
+                     zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! C2 interpolation of T at u- & v-points (x2)
+                     zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
                      !                                                  ! C4 interpolation of T at u- & v-points (x2)
                      zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj,jk) - ztu(ji+1,jj,jk) )
                      zC4t_v =  zC2t_v + r1_6 * ( ztv(ji,jj-1,jk) - ztv(ji,jj+1,jk) )
                      !                                                  ! C4 fluxes
-                     zwx(ji,jj,jk) =  0.5_wp * pun(ji,jj,jk) * zC4t_u
-                     zwy(ji,jj,jk) =  0.5_wp * pvn(ji,jj,jk) * zC4t_v
+                     zwx(ji,jj,jk) =  0.5_wp * pu_mm(ji,jj,jk) * zC4t_u
+                     zwy(ji,jj,jk) =  0.5_wp * pv_mm(ji,jj,jk) * zC4t_v
                   END DO
                END DO
@@ -150,5 +150,5 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)
+                     zwz(ji,jj,jk) = 0.5 * pww(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) ) * wmask(ji,jj,jk)
                   END DO
                END DO
@@ -156,9 +156,9 @@
             !
          CASE(  4  )                         !* 4th order compact
-            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )      ! ztw = interpolated value of T at w-point
+            CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )      ! ztw = interpolated value of T at w-point
             DO jk = 2, jpkm1
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     zwz(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
+                     zwz(ji,jj,jk) = pww(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
                   END DO
                END DO
@@ -171,9 +171,9 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn) 
+                     zwz(ji,jj, mikt(ji,jj) ) = pww(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm) 
                   END DO
                END DO   
             ELSE                                   ! no ice-shelf cavities (only ocean surface)
-               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
+               zwz(:,:,1) = pww(:,:,1) * pt(:,:,1,jn,Kmm)
             ENDIF
          ENDIF
@@ -182,8 +182,8 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn)    &
+                  pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs)    &
                      &             - (  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)  ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                     &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                END DO
             END DO
@@ -191,7 +191,7 @@
          !                             ! trend diagnostics
          IF( l_trd ) THEN
-            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
-            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
-            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
+            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pu_mm, pt(:,:,:,jn,Kmm) )
+            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pv_mm, pt(:,:,:,jn,Kmm) )
+            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pww, pt(:,:,:,jn,Kmm) )
          END IF
          !                                 ! "Poleward" heat and salt transports 
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_fct.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_fct.F90	(revision 10879)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_fct.F90	(revision 10880)
@@ -52,6 +52,6 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, pun, pvn, pwn,       &
-      &                                              ptb, ptn, pta, kjpt, kn_fct_h, kn_fct_v )
+   SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, pu_mm, pv_mm, pww,       &
+      &                    Kbb, Kmm, pt, kjpt, Krhs, kn_fct_h, kn_fct_v )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_fct  ***
@@ -65,18 +65,18 @@
       !!               - corrected flux (monotonic correction) 
       !!
-      !! ** Action : - update pta  with the now advective tracer trends
+      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
       !!             - send trends to trdtra module for further diagnostics (l_trdtra=T)
       !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
       !!----------------------------------------------------------------------
-      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
-      INTEGER                              , INTENT(in   ) ::   kn_fct_h        ! order of the FCT scheme (=2 or 4)
-      INTEGER                              , INTENT(in   ) ::   kn_fct_v        ! order of the FCT scheme (=2 or 4)
-      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! 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 
+      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
+      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
+      INTEGER                                  , INTENT(in   ) ::   kn_fct_h        ! order of the FCT scheme (=2 or 4)
+      INTEGER                                  , INTENT(in   ) ::   kn_fct_v        ! order of the FCT scheme (=2 or 4)
+      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pu_mm, pv_mm, pww   ! 3 ocean velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
       !
       INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices  
@@ -125,10 +125,10 @@
                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 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn) )
-                  zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn) )
+                  zfp_ui = pu_mm(ji,jj,jk) + ABS( pu_mm(ji,jj,jk) )
+                  zfm_ui = pu_mm(ji,jj,jk) - ABS( pu_mm(ji,jj,jk) )
+                  zfp_vj = pv_mm(ji,jj,jk) + ABS( pv_mm(ji,jj,jk) )
+                  zfm_vj = pv_mm(ji,jj,jk) - ABS( pv_mm(ji,jj,jk) )
+                  zwx(ji,jj,jk) = 0.5 * ( zfp_ui * pt(ji,jj,jk,jn,Kbb) + zfm_ui * pt(ji+1,jj  ,jk,jn,Kbb) )
+                  zwy(ji,jj,jk) = 0.5 * ( zfp_vj * pt(ji,jj,jk,jn,Kbb) + zfm_vj * pt(ji  ,jj+1,jk,jn,Kbb) )
                END DO
             END DO
@@ -138,7 +138,7 @@
             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 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)
+                  zfp_wk = pww(ji,jj,jk) + ABS( pww(ji,jj,jk) )
+                  zfm_wk = pww(ji,jj,jk) - ABS( pww(ji,jj,jk) )
+                  zwz(ji,jj,jk) = 0.5 * ( zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb) ) * wmask(ji,jj,jk)
                END DO
             END DO
@@ -148,9 +148,9 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface 
+                     zwz(ji,jj, mikt(ji,jj) ) = pww(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface 
                   END DO
                END DO   
             ELSE                             ! no cavities: only at the ocean surface
-               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
+               zwz(:,:,1) = pww(:,:,1) * pt(:,:,1,jn,Kbb)
             ENDIF
          ENDIF
@@ -164,6 +164,6 @@
                      &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
                   !                             ! update and guess with monotonic sheme
-                  pta(ji,jj,jk,jn) =                     pta(ji,jj,jk,jn) +        ztra   / e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
-                  zwi(ji,jj,jk)    = ( e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * ztra ) / e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
+                  pt(ji,jj,jk,jn,Krhs) =                     pt(ji,jj,jk,jn,Krhs) +        ztra   / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
+                  zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -184,6 +184,6 @@
                DO jj = 1, jpjm1
                   DO ji = 1, fs_jpim1   ! vector opt.
-                     zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk)
-                     zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk)
+                     zwx(ji,jj,jk) = 0.5_wp * pu_mm(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj,jk,jn,Kmm) ) - zwx(ji,jj,jk)
+                     zwy(ji,jj,jk) = 0.5_wp * pv_mm(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj+1,jk,jn,Kmm) ) - zwy(ji,jj,jk)
                   END DO
                END DO
@@ -196,6 +196,6 @@
                DO jj = 1, jpjm1                    ! 1st derivative (gradient)
                   DO ji = 1, fs_jpim1   ! vector opt.
-                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
-                     ztv(ji,jj,jk) = ( ptn(ji  ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
+                     ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
+                     ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
                   END DO
                END DO
@@ -212,9 +212,9 @@
                DO jj = 1, jpjm1
                   DO ji = 1, fs_jpim1   ! vector opt.
-                     zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn)   ! 2 x C2 interpolation of T at u- & v-points
-                     zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn)
+                     zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points
+                     zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
                      !                                                  ! C4 minus upstream advective fluxes 
-                     zwx(ji,jj,jk) =  0.5_wp * pun(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk)
-                     zwy(ji,jj,jk) =  0.5_wp * pvn(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk)
+                     zwx(ji,jj,jk) =  0.5_wp * pu_mm(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk)
+                     zwy(ji,jj,jk) =  0.5_wp * pv_mm(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk)
                   END DO
                END DO
@@ -227,6 +227,6 @@
                DO jj = 1, jpjm1
                   DO ji = 1, fs_jpim1   ! vector opt.
-                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
-                     ztv(ji,jj,jk) = ( ptn(ji  ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
+                     ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
+                     ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
                   END DO
                END DO
@@ -237,12 +237,12 @@
                DO jj = 2, jpjm1
                   DO ji = 2, fs_jpim1   ! vector opt.
-                     zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn)   ! 2 x C2 interpolation of T at u- & v-points (x2)
-                     zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn)
+                     zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points (x2)
+                     zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
                      !                                                  ! C4 interpolation of T at u- & v-points (x2)
                      zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj  ,jk) - ztu(ji+1,jj  ,jk) )
                      zC4t_v =  zC2t_v + r1_6 * ( ztv(ji  ,jj-1,jk) - ztv(ji  ,jj+1,jk) )
                      !                                                  ! C4 minus upstream advective fluxes 
-                     zwx(ji,jj,jk) =  0.5_wp * pun(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk)
-                     zwy(ji,jj,jk) =  0.5_wp * pvn(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk)
+                     zwx(ji,jj,jk) =  0.5_wp * pu_mm(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk)
+                     zwy(ji,jj,jk) =  0.5_wp * pv_mm(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk)
                   END DO
                END DO
@@ -257,5 +257,5 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     zwz(ji,jj,jk) =  (  pwn(ji,jj,jk) * 0.5_wp * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )   &
+                     zwz(ji,jj,jk) =  (  pww(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   &
                         &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk)
                   END DO
@@ -264,9 +264,9 @@
             !
          CASE(  4  )                   !- 4th order COMPACT
-            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )   ! zwt = COMPACT interpolation of T at w-point
+            CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )   ! zwt = COMPACT interpolation of T at w-point
             DO jk = 2, jpkm1
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
+                     zwz(ji,jj,jk) = ( pww(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
                   END DO
                END DO
@@ -282,5 +282,5 @@
          !        !==  monotonicity algorithm  ==!
          !
-         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
+         CALL nonosc( Kmm, pt(:,:,:,jn,Kbb), zwx, zwy, zwz, zwi, p2dt )
          !
          !        !==  final trend with corrected fluxes  ==!
@@ -289,8 +289,8 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.  
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                  pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - (  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) ) &
-                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                     &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                END DO
             END DO
@@ -303,7 +303,7 @@
             !
             IF( l_trd ) THEN              ! trend diagnostics
-               CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
-               CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )
-               CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )
+               CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pu_mm, pt(:,:,:,jn,Kmm) )
+               CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pv_mm, pt(:,:,:,jn,Kmm) )
+               CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pww, pt(:,:,:,jn,Kmm) )
             ENDIF
             !                             ! heat/salt transport
@@ -328,5 +328,5 @@
 
 
-   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
+   SUBROUTINE nonosc( Kmm, pbef, paa, pbb, pcc, paft, p2dt )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE nonosc  ***
@@ -341,4 +341,5 @@
       !!       in-space based differencing for fluid
       !!----------------------------------------------------------------------
+      INTEGER                          , INTENT(in   ) ::   Kmm             ! time level index 
       REAL(wp)                         , INTENT(in   ) ::   p2dt            ! tracer time-step
       REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field
@@ -392,5 +393,5 @@
 
                ! up & down beta terms
-               zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt
+               zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt
                zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
                zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_mus.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_mus.F90	(revision 10879)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_mus.F90	(revision 10880)
@@ -54,6 +54,6 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pun, pvn, pwn,             &
-      &                                              ptb, pta, kjpt, ld_msc_ups )
+   SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pu_mm, pv_mm, pww,             &
+      &                    Kbb, Kmm, pt, kjpt, Krhs, ld_msc_ups )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE tra_adv_mus  ***
@@ -66,5 +66,5 @@
       !!              ld_msc_ups=T : 
       !!
-      !! ** Action : - update pta  with the now advective tracer trends
+      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
       !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
       !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
@@ -73,13 +73,13 @@
       !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
       !!----------------------------------------------------------------------
-      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
-      LOGICAL                              , INTENT(in   ) ::   ld_msc_ups      ! use upstream scheme within muscl
-      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! 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             ! before tracer field
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend 
+      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
+      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
+      LOGICAL                                  , INTENT(in   ) ::   ld_msc_ups      ! use upstream scheme within muscl
+      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pu_mm, pv_mm, pww   ! 3 ocean velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -134,6 +134,6 @@
             DO jj = 1, jpjm1      
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwx(ji,jj,jk) = umask(ji,jj,jk) * ( ptb(ji+1,jj,jk,jn) - ptb(ji,jj,jk,jn) )
-                  zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( ptb(ji,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
+                  zwx(ji,jj,jk) = umask(ji,jj,jk) * ( pt(ji+1,jj,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
+                  zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( pt(ji,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
                END DO
            END DO
@@ -172,17 +172,17 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! MUSCL fluxes
-                  z0u = SIGN( 0.5, pun(ji,jj,jk) )
+                  z0u = SIGN( 0.5, pu_mm(ji,jj,jk) )
                   zalpha = 0.5 - z0u
-                  zu  = z0u - 0.5 * pun(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
-                  zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk)
-                  zzwy = ptb(ji  ,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk)
-                  zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
+                  zu  = z0u - 0.5 * pu_mm(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
+                  zzwx = pt(ji+1,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk)
+                  zzwy = pt(ji  ,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk)
+                  zwx(ji,jj,jk) = pu_mm(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
                   !
-                  z0v = SIGN( 0.5, pvn(ji,jj,jk) )
+                  z0v = SIGN( 0.5, pv_mm(ji,jj,jk) )
                   zalpha = 0.5 - z0v
-                  zv  = z0v - 0.5 * pvn(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
-                  zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk)
-                  zzwy = ptb(ji,jj  ,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk)
-                  zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
+                  zv  = z0v - 0.5 * pv_mm(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
+                  zzwx = pt(ji,jj+1,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk)
+                  zzwy = pt(ji,jj  ,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk)
+                  zwy(ji,jj,jk) = pv_mm(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
                END DO
             END DO
@@ -193,7 +193,7 @@
             DO jj = 2, jpjm1      
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       &
+                  pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       &
                   &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     &
-                  &                                   * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                  &                                   * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                END DO
            END DO
@@ -201,6 +201,6 @@
          !                                ! trend diagnostics
          IF( l_trd )  THEN
-            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) )
-            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) )
+            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pu_mm, pt(:,:,:,jn,Kbb) )
+            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pv_mm, pt(:,:,:,jn,Kbb) )
          END IF
          !                                 ! "Poleward" heat and salt transports 
@@ -215,5 +215,5 @@
          zwx(:,:,jpk) = 0._wp
          DO jk = 2, jpkm1                       ! interior values
-            zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) )
+            zwx(:,:,jk) = tmask(:,:,jk) * ( pt(:,:,jk-1,jn,Kbb) - pt(:,:,jk,jn,Kbb) )
          END DO
          !                                !-- Slopes of tracer
@@ -239,10 +239,10 @@
             DO jj = 2, jpjm1      
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )
+                  z0w = SIGN( 0.5, pww(ji,jj,jk+1) )
                   zalpha = 0.5 + z0w
-                  zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w_n(ji,jj,jk+1)
-                  zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1)
-                  zzwy = ptb(ji,jj,jk  ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  )
-                  zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)
+                  zw  = z0w - 0.5 * pww(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,Kmm)
+                  zzwx = pt(ji,jj,jk+1,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1)
+                  zzwy = pt(ji,jj,jk  ,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  )
+                  zwx(ji,jj,jk+1) = pww(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)
                END DO 
             END DO
@@ -252,9 +252,9 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zwx(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)
+                     zwx(ji,jj, mikt(ji,jj) ) = pww(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)
                   END DO
                END DO   
             ELSE                                      ! no cavities: only at the ocean surface
-               zwx(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
+               zwx(:,:,1) = pww(:,:,1) * pt(:,:,1,jn,Kbb)
             ENDIF
          ENDIF
@@ -263,10 +263,10 @@
             DO jj = 2, jpjm1      
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  pta(ji,jj,jk,jn) =  pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                  pt(ji,jj,jk,jn,Krhs) =  pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                END DO
             END DO
          END DO
          !                                ! send trends for diagnostic
-         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) )
+         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pww, pt(:,:,:,jn,Kbb) )
          !
       END DO                     ! end of tracer loop
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_qck.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_qck.F90	(revision 10879)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_qck.F90	(revision 10880)
@@ -47,6 +47,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
-      &                                       ptb, ptn, pta, kjpt )
+   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pu_mm, pv_mm, pww, Kbb, Kmm, pt, kjpt, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_qck  ***
@@ -72,5 +71,5 @@
       !!         dt = 2*rdtra and the scalar values are tb and sb
       !!
-      !!       On the vertical, the simple centered scheme used ptn
+      !!       On the vertical, the simple centered scheme used pt(:,:,:,:,Kmm)
       !!
       !!               The fluxes are bounded by the ULTIMATE limiter to
@@ -78,5 +77,5 @@
       !!            prevent the appearance of spurious numerical oscillations
       !!
-      !! ** Action : - update pta  with the now advective tracer trends
+      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
       !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
       !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
@@ -84,12 +83,12 @@
       !! ** Reference : Leonard (1979, 1991)
       !!----------------------------------------------------------------------
-      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)                             , INTENT(in   ) ::   p2dt            ! 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 
+      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
+      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)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pu_mm, pv_mm, pww   ! 3 ocean velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
       !!----------------------------------------------------------------------
       !
@@ -108,25 +107,24 @@
       !
       !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
-      CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt ) 
-      CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 
+      CALL tra_adv_qck_i( kt, cdtype, p2dt, pu_mm, Kbb, Kmm, pt, kjpt, Krhs ) 
+      CALL tra_adv_qck_j( kt, cdtype, p2dt, pv_mm, Kbb, Kmm, pt, kjpt, Krhs ) 
 
       !        ! vertical fluxes are computed with the 2nd order centered scheme
-      CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt )
+      CALL tra_adv_cen2_k( kt, cdtype, pww, Kmm, pt, kjpt, Krhs )
       !
    END SUBROUTINE tra_adv_qck
 
 
-   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun,                  &
-      &                                        ptb, ptn, pta, kjpt   )
-      !!----------------------------------------------------------------------
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
-      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
-      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
-      REAL(wp)                             , INTENT(in   ) ::   p2dt       ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun        ! i-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 
+   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pu_mm, Kbb, Kmm, pt, kjpt, Krhs )
+      !!----------------------------------------------------------------------
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
+      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
+      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers
+      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pu_mm        ! i-velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
       !!
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -145,6 +143,6 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)        ! Upstream   in the x-direction for the tracer
-                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)        ! Downstream in the x-direction for the tracer
+                  zfc(ji,jj,jk) = pt(ji-1,jj,jk,jn,Kbb)        ! Upstream   in the x-direction for the tracer
+                  zfd(ji,jj,jk) = pt(ji+1,jj,jk,jn,Kbb)        ! Downstream in the x-direction for the tracer
                END DO
             END DO
@@ -158,5 +156,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.         
-                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0 
+                  zdir = 0.5 + SIGN( 0.5, pu_mm(ji,jj,jk) )   ! if pu_mm > 0 : zdir = 1 otherwise zdir = 0 
                   zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T 
                END DO
@@ -167,9 +165,9 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.   
-                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0 
-                  zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u_n(ji,jj,jk)
-                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
-                  zfc(ji,jj,jk)  = zdir * ptb(ji  ,jj,jk,jn) + ( 1. - zdir ) * ptb(ji+1,jj,jk,jn)  ! FC in the x-direction for T
-                  zfd(ji,jj,jk)  = zdir * ptb(ji+1,jj,jk,jn) + ( 1. - zdir ) * ptb(ji  ,jj,jk,jn)  ! FD in the x-direction for T
+                  zdir = 0.5 + SIGN( 0.5, pu_mm(ji,jj,jk) )   ! if pu_mm > 0 : zdir = 1 otherwise zdir = 0 
+                  zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
+                  zwx(ji,jj,jk)  = ABS( pu_mm(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
+                  zfc(ji,jj,jk)  = zdir * pt(ji  ,jj,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji+1,jj,jk,jn,Kbb)  ! FC in the x-direction for T
+                  zfd(ji,jj,jk)  = zdir * pt(ji+1,jj,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji  ,jj,jk,jn,Kbb)  ! FD in the x-direction for T
                END DO
             END DO
@@ -197,10 +195,10 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.               
-                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0 
+                  zdir = 0.5 + SIGN( 0.5, pu_mm(ji,jj,jk) )   ! if pu_mm > 0 : zdir = 1 otherwise zdir = 0 
                   !--- If the second ustream point is a land point
                   !--- the flux is computed by the 1st order UPWIND scheme
                   zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
                   zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
-                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
+                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pu_mm(ji,jj,jk)
                END DO
             END DO
@@ -213,14 +211,14 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.  
-                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                  zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                   ! horizontal advective trends
                   ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
                   !--- add it to the general tracer trends
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
+                  pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra
                END DO
             END DO
          END DO
          !                                 ! trend diagnostics
-         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
+         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pu_mm, pt(:,:,:,jn,Kmm) )
          !
       END DO
@@ -229,16 +227,15 @@
 
 
-   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                &
-      &                                        ptb, ptn, pta, kjpt )
-      !!----------------------------------------------------------------------
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
-      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
-      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
-      REAL(wp)                             , INTENT(in   ) ::   p2dt       ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn        ! j-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 
+   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pv_mm, Kbb, Kmm, pt, kjpt, Krhs )
+      !!----------------------------------------------------------------------
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
+      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
+      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers
+      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pv_mm        ! j-velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
       !!
       INTEGER  :: ji, jj, jk, jn                ! dummy loop indices
@@ -259,7 +256,7 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! Upstream in the x-direction for the tracer
-                  zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
+                  zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb)
                   ! Downstream in the x-direction for the tracer
-                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
+                  zfd(ji,jj,jk) = pt(ji,jj+1,jk,jn,Kbb)
                END DO
             END DO
@@ -275,5 +272,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.         
-                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0 
+                  zdir = 0.5 + SIGN( 0.5, pv_mm(ji,jj,jk) )   ! if pu_mm > 0 : zdir = 1 otherwise zdir = 0 
                   zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T 
                END DO
@@ -284,9 +281,9 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.   
-                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0 
-                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v_n(ji,jj,jk)
-                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
-                  zfc(ji,jj,jk)  = zdir * ptb(ji,jj  ,jk,jn) + ( 1. - zdir ) * ptb(ji,jj+1,jk,jn)  ! FC in the x-direction for T
-                  zfd(ji,jj,jk)  = zdir * ptb(ji,jj+1,jk,jn) + ( 1. - zdir ) * ptb(ji,jj  ,jk,jn)  ! FD in the x-direction for T
+                  zdir = 0.5 + SIGN( 0.5, pv_mm(ji,jj,jk) )   ! if pu_mm > 0 : zdir = 1 otherwise zdir = 0 
+                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
+                  zwy(ji,jj,jk)  = ABS( pv_mm(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
+                  zfc(ji,jj,jk)  = zdir * pt(ji,jj  ,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji,jj+1,jk,jn,Kbb)  ! FC in the x-direction for T
+                  zfd(ji,jj,jk)  = zdir * pt(ji,jj+1,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji,jj  ,jk,jn,Kbb)  ! FD in the x-direction for T
                END DO
             END DO
@@ -314,10 +311,10 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.               
-                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0 
+                  zdir = 0.5 + SIGN( 0.5, pv_mm(ji,jj,jk) )   ! if pu_mm > 0 : zdir = 1 otherwise zdir = 0 
                   !--- If the second ustream point is a land point
                   !--- the flux is computed by the 1st order UPWIND scheme
                   zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
                   zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
-                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
+                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pv_mm(ji,jj,jk)
                END DO
             END DO
@@ -330,14 +327,14 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.  
-                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                  zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                   ! horizontal advective trends
                   ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
                   !--- add it to the general tracer trends
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
+                  pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra
                END DO
             END DO
          END DO
          !                                 ! trend diagnostics
-         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
+         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pv_mm, pt(:,:,:,jn,Kmm) )
          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
          IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
@@ -348,15 +345,14 @@
 
 
-   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           &
-     &                                    ptn, pta, kjpt )
-      !!----------------------------------------------------------------------
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
-      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
-      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity 
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn      ! before and now tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend 
+   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pww, Kmm, pt, kjpt, Krhs )
+      !!----------------------------------------------------------------------
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER                                  , INTENT(in   ) ::   kt       ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) ::   Kmm, Krhs  ! ocean time level indices
+      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
+      INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pww      ! vertical velocity 
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -374,5 +370,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) * wmask(ji,jj,jk)
+                  zwz(ji,jj,jk) = 0.5 * pww(ji,jj,jk) * ( pt(ji,jj,jk-1,jn,Kmm) + pt(ji,jj,jk,jn,Kmm) ) * wmask(ji,jj,jk)
                END DO
             END DO
@@ -382,9 +378,9 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn)   ! linear free surface 
+                     zwz(ji,jj, mikt(ji,jj) ) = pww(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)   ! linear free surface 
                   END DO
                END DO   
             ELSE                                   ! no ocean cavities (only ocean surface)
-               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
+               zwz(:,:,1) = pww(:,:,1) * pt(:,:,1,jn,Kmm)
             ENDIF
          ENDIF
@@ -393,11 +389,11 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
-                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                  pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
+                     &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                END DO
             END DO
          END DO
          !                                 ! Send trends for diagnostic
-         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
+         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pww, pt(:,:,:,jn,Kmm) )
          !
       END DO
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_ubs.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_ubs.F90	(revision 10879)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_ubs.F90	(revision 10880)
@@ -46,6 +46,6 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pun, pvn, pwn,          &
-      &                                                ptb, ptn, pta, kjpt, kn_ubs_v )
+   SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pu_mm, pv_mm, pww,          &
+      &                    Kbb, Kmm, pt, kjpt, Krhs, kn_ubs_v )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_ubs  ***
@@ -58,7 +58,7 @@
       !!      It is only used in the horizontal direction.
       !!      For example the i-component of the advective fluxes are given by :
-      !!                !  e2u e3u un ( mi(Tn) - zltu(i  ) )   if un(i) >= 0
+      !!                !  e2u e3u uu ( mi(Tn) - zltu(i  ) ,Kmm)   if uu(i,Kmm) >= 0
       !!          ztu = !  or 
-      !!                !  e2u e3u un ( mi(Tn) - zltu(i+1) )   if un(i) < 0
+      !!                !  e2u e3u uu ( mi(Tn) - zltu(i+1) ,Kmm)   if uu(i,Kmm) < 0
       !!      where zltu is the second derivative of the before temperature field:
       !!          zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ]
@@ -77,5 +77,5 @@
       !!      scheme (kn_ubs_v=4).
       !!
-      !! ** Action : - update pta  with the now advective tracer trends
+      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
       !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
       !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
@@ -84,13 +84,13 @@
       !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741. 
       !!----------------------------------------------------------------------
-      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
-      INTEGER                              , INTENT(in   ) ::   kn_ubs_v        ! number of tracers
-      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean transport 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 
+      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
+      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
+      INTEGER                                  , INTENT(in   ) ::   kn_ubs_v        ! number of tracers
+      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pu_mm, pv_mm, pww   ! 3 ocean transport components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -126,13 +126,13 @@
             DO jj = 1, jpjm1              ! First derivative (masked gradient)
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zeeu = e2_e1u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk)
-                  zeev = e1_e2v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk)
-                  ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
-                  ztv(ji,jj,jk) = zeev * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
+                  zeeu = e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
+                  zeev = e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
+                  ztu(ji,jj,jk) = zeeu * ( pt(ji+1,jj  ,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
+                  ztv(ji,jj,jk) = zeev * ( pt(ji  ,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
                END DO
             END DO
             DO jj = 2, jpjm1              ! Second derivative (divergence)
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zcoef = 1._wp / ( 6._wp * e3t_n(ji,jj,jk) )
+                  zcoef = 1._wp / ( 6._wp * e3t(ji,jj,jk,Kmm) )
                   zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
                   zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef
@@ -146,11 +146,11 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )      ! upstream transport (x2)
-                  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) )
+                  zfp_ui = pu_mm(ji,jj,jk) + ABS( pu_mm(ji,jj,jk) )      ! upstream transport (x2)
+                  zfm_ui = pu_mm(ji,jj,jk) - ABS( pu_mm(ji,jj,jk) )
+                  zfp_vj = pv_mm(ji,jj,jk) + ABS( pv_mm(ji,jj,jk) )
+                  zfm_vj = pv_mm(ji,jj,jk) - ABS( pv_mm(ji,jj,jk) )
                   !                                                  ! 2nd order centered advective fluxes (x2)
-                  zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
-                  zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
+                  zcenut = pu_mm(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm) )
+                  zcenvt = pv_mm(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) )
                   !                                                  ! UBS advective fluxes
                   ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) )
@@ -160,12 +160,12 @@
          END DO         
          !
-         zltu(:,:,:) = pta(:,:,:,jn)      ! store the initial trends before its update
+         zltu(:,:,:) = pt(:,:,:,jn,Krhs)      ! store the initial trends before its update
          !
          DO jk = 1, jpkm1        !==  add the horizontal advective trend  ==!
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                        &
+                  pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs)                        &
                      &             - (  ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk)    &
-                     &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                     &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                END DO
             END DO
@@ -173,10 +173,10 @@
          END DO
          !
-         zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:)    ! Horizontal advective trend used in vertical 2nd order FCT case
+         zltu(:,:,:) = pt(:,:,:,jn,Krhs) - zltu(:,:,:)    ! Horizontal advective trend used in vertical 2nd order FCT case
          !                                            ! and/or in trend diagnostic (l_trd=T) 
          !                
          IF( l_trd ) THEN                  ! trend diagnostics
-             CALL trd_tra( kt, cdtype, jn, jptra_xad, ztu, pun, ptn(:,:,:,jn) )
-             CALL trd_tra( kt, cdtype, jn, jptra_yad, ztv, pvn, ptn(:,:,:,jn) )
+             CALL trd_tra( kt, cdtype, jn, jptra_xad, ztu, pu_mm, pt(:,:,:,jn,Kmm) )
+             CALL trd_tra( kt, cdtype, jn, jptra_yad, ztv, pv_mm, pt(:,:,:,jn,Kmm) )
          END IF
          !     
@@ -193,5 +193,5 @@
          CASE(  2  )                   ! 2nd order FCT 
             !         
-            IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag.
+            IF( l_trd )   zltv(:,:,:) = pt(:,:,:,jn,Krhs)          ! store pt(:,:,:,:,Krhs) if trend diag.
             !
             !                          !*  upstream advection with initial mass fluxes & intermediate update  ==!
@@ -199,7 +199,7 @@
                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) )
-                     ztw(ji,jj,jk) = 0.5_wp * (  zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn)  ) * wmask(ji,jj,jk)
+                     zfp_wk = pww(ji,jj,jk) + ABS( pww(ji,jj,jk) )
+                     zfm_wk = pww(ji,jj,jk) - ABS( pww(ji,jj,jk) )
+                     ztw(ji,jj,jk) = 0.5_wp * (  zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb)  ) * wmask(ji,jj,jk)
                   END DO
                END DO
@@ -209,9 +209,9 @@
                   DO jj = 1, jpj
                      DO ji = 1, jpi
-                        ztw(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface 
+                        ztw(ji,jj, mikt(ji,jj) ) = pww(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface 
                      END DO
                   END DO   
                ELSE                                ! no cavities: only at the ocean surface
-                  ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
+                  ztw(:,:,1) = pww(:,:,1) * pt(:,:,1,jn,Kbb)
                ENDIF
             ENDIF
@@ -220,7 +220,7 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
-                     pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak 
-                     zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
+                     ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+                     pt(ji,jj,jk,jn,Krhs) =   pt(ji,jj,jk,jn,Krhs) +  ztak 
+                     zti(ji,jj,jk)    = ( pt(ji,jj,jk,jn,Kbb) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
                   END DO
                END DO
@@ -232,5 +232,5 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     ztw(ji,jj,jk) = (   0.5_wp * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )   &
+                     ztw(ji,jj,jk) = (   0.5_wp * pww(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   &
                         &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk)
                   END DO
@@ -240,16 +240,16 @@
             IF( ln_linssh )   ztw(:,:, 1 ) = 0._wp       ! only ocean surface as interior zwz values have been w-masked
             !
-            CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm
+            CALL nonosc_z( Kmm, pt(:,:,:,jn,Kbb), ztw, zti, p2dt )      !  monotonicity algorithm
             !
          CASE(  4  )                               ! 4th order COMPACT
-            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )         ! 4th order compact interpolation of T at w-point
+            CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )         ! 4th order compact interpolation of T at w-point
             DO jk = 2, jpkm1
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     ztw(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-            IF( ln_linssh )   ztw(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)     !!gm ISF & 4th COMPACT doesn't work
+                     ztw(ji,jj,jk) = pww(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            IF( ln_linssh )   ztw(:,:, 1 ) = pww(:,:,1) * pt(:,:,1,jn,Kmm)     !!gm ISF & 4th COMPACT doesn't work
             !
          END SELECT
@@ -258,5 +258,5 @@
             DO jj = 2, jpjm1 
                DO ji = fs_2, fs_jpim1   ! vector opt.   
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                  pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                END DO
             END DO
@@ -267,7 +267,7 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk)                          &
-                        &           + ptn(ji,jj,jk,jn) * (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  )   &
-                        &                              * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                     zltv(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) - zltv(ji,jj,jk)                          &
+                        &           + pt(ji,jj,jk,jn,Kmm) * (  pww(ji,jj,jk) - pww(ji,jj,jk+1)  )   &
+                        &                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                   END DO
                END DO
@@ -281,5 +281,5 @@
 
 
-   SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt )
+   SUBROUTINE nonosc_z( Kmm, pbef, pcc, paft, p2dt )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE nonosc_z  ***
@@ -294,4 +294,5 @@
       !!       in-space based differencing for fluid
       !!----------------------------------------------------------------------
+      INTEGER , INTENT(in   )                          ::   Kmm    ! time level index
       REAL(wp), INTENT(in   )                          ::   p2dt   ! tracer time-step
       REAL(wp),                DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field
@@ -352,5 +353,5 @@
                zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
                ! up & down beta terms
-               zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt
+               zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt
                zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
                zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10879)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10880)
@@ -47,4 +47,9 @@
 
    !!----------------------------------------------------------------------
+   !! time level indices
+   !!----------------------------------------------------------------------
+   INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs          !! used by nemo_init
+
+   !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
    !! $Id$
@@ -247,5 +252,5 @@
                &         CALL Agrif_Sponge_tra        ! tracers sponge
 #endif
-                         CALL tra_adv       ( kstp )  ! horizontal & vertical advection
+                         CALL tra_adv( kstp, Nbb, Nnn      , ts, Nrhs )  ! hor. + vert. advection	==> RHS
       IF( ln_zdfosm  )   CALL tra_osm       ( kstp )  ! OSMOSIS non-local tracer fluxes
       IF( lrst_oce .AND. ln_zdfosm ) &
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step_oce.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step_oce.F90	(revision 10879)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step_oce.F90	(revision 10880)
@@ -107,5 +107,5 @@
 #endif
 #if defined key_top
-   USE trcstp           ! passive tracer time-stepping      (trc_stp routine)
+   USE trcstp, ONLY : trc_stp    ! passive tracer time-stepping      (trc_stp routine)
 #endif
    !!----------------------------------------------------------------------
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcadv.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcadv.F90	(revision 10879)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcadv.F90	(revision 10880)
@@ -68,5 +68,5 @@
 CONTAINS
 
-   SUBROUTINE trc_adv( kt )
+   SUBROUTINE trc_adv( kt, Kbb, Kmm, ptr, Krhs  )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trc_adv  ***
@@ -76,5 +76,7 @@
       !! ** Method  : - Update after tracers (tra) with the advection term following nadv
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER                                   , INTENT(in)    :: kt   ! ocean time-step index
+      INTEGER                                   , INTENT(in)    :: Kbb, Kmm, Krhs ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr            ! passive tracers and RHS of tracer equation
       !
       INTEGER ::   jk   ! dummy loop index
@@ -123,13 +125,13 @@
       !
       CASE ( np_CEN )                                 ! Centered : 2nd / 4th order
-         CALL tra_adv_cen( kt, nittrc000,'TRC',          zun, zvn, zwn     , trn, tra, jptra, nn_cen_h, nn_cen_v )
+         CALL tra_adv_cen( kt, nittrc000,'TRC',          zun, zvn, zwn,      Kmm, ptr, jptra, Krhs, nn_cen_h, nn_cen_v )
       CASE ( np_FCT )                                 ! FCT      : 2nd / 4th order
-         CALL tra_adv_fct( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, trb, trn, tra, jptra, nn_fct_h, nn_fct_v )
+         CALL tra_adv_fct( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, Kbb, Kmm, ptr, jptra, Krhs, nn_fct_h, nn_fct_v )
       CASE ( np_MUS )                                 ! MUSCL
-         CALL tra_adv_mus( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, trb,      tra, jptra        , ln_mus_ups ) 
+         CALL tra_adv_mus( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, Kbb, Kmm, ptr, jptra, Krhs, ln_mus_ups         ) 
       CASE ( np_UBS )                                 ! UBS
-         CALL tra_adv_ubs( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, trb, trn, tra, jptra          , nn_ubs_v )
+         CALL tra_adv_ubs( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, Kbb, Kmm, ptr, jptra, Krhs, nn_ubs_v           )
       CASE ( np_QCK )                                 ! QUICKEST
-         CALL tra_adv_qck( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, trb, trn, tra, jptra                     )
+         CALL tra_adv_qck( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, Kbb, Kmm, ptr, jptra, Krhs                     )
       !
       END SELECT
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90	(revision 10879)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90	(revision 10880)
@@ -44,5 +44,5 @@
 CONTAINS
 
-   SUBROUTINE trc_trp( kt )
+   SUBROUTINE trc_trp( kt, Kbb, Kmm, Krhs, Kaa )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE trc_trp  ***
@@ -53,5 +53,6 @@
       !!              - Update the passive tracers
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::  kt  ! ocean time-step index
+      INTEGER, INTENT( in ) :: kt                  ! ocean time-step index
+      INTEGER, INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! time level indices
       !! ---------------------------------------------------------------------
       !
@@ -64,5 +65,5 @@
          IF( ln_trcdmp )        CALL trc_dmp    ( kt )      ! internal damping trends
          IF( ln_bdy )           CALL trc_bdy_dmp( kt )      ! BDY damping trends
-                                CALL trc_adv    ( kt )      ! horizontal & vertical advection 
+                                CALL trc_adv    ( kt, Kbb, Kmm, tr, Krhs )      ! horizontal & vertical advection 
          !                                                         ! Partial top/bottom cell: GRADh( trb )  
          IF( ln_zps ) THEN
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/oce_trc.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/oce_trc.F90	(revision 10879)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/oce_trc.F90	(revision 10880)
@@ -8,4 +8,5 @@
    !!----------------------------------------------------------------------
    !                                            !* Domain size *
+   USE par_oce , ONLY :   jpt      =>   jpt        !: time dimension
    USE par_oce , ONLY :   jpi      =>   jpi        !: first  dimension of grid --> i 
    USE par_oce , ONLY :   jpj      =>   jpj        !: second dimension of grid --> j  
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trc.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trc.F90	(revision 10879)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trc.F90	(revision 10880)
@@ -33,7 +33,5 @@
    REAL(wp), PUBLIC                                        ::  areatot        !: total volume 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:  ) ::  cvol           !: volume correction -degrad option- 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  trn            !: tracer concentration for now time step
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  tra            !: tracer concentration for next time step
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  trb            !: tracer concentration for before time step
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:), TARGET ::  tr             !: tracer concentration 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:  ) ::  sbc_trc_b      !: Before sbc fluxes for tracers
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:  ) ::  sbc_trc        !: Now sbc fluxes for tracers
@@ -42,4 +40,10 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:  ) ::  trc_o          !: prescribed tracer concentration in ocean for SBC
    INTEGER             , PUBLIC                            ::  nn_ice_tr      !: handling of sea ice tracers
+
+   !! TEMPORARY POINTERS - TO BE DELETED AFTER IMMERSE DEVELOPMENT COMPLETE
+   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:,:)   ::  trn            !: tracer concentration for now time step
+   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:,:)   ::  tra            !: tracer concentration for next time step
+   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:,:)   ::  trb            !: tracer concentration for before time step
+   !! TEMPORARY POINTERS - TO BE DELETED AFTER IMMERSE DEVELOPMENT COMPLETE
 
    !! interpolated gradient
@@ -147,5 +151,5 @@
       ierr(:) = 0
       !
-      ALLOCATE( trn(jpi,jpj,jpk,jptra), trb(jpi,jpj,jpk,jptra), tra(jpi,jpj,jpk,jptra),       &  
+      ALLOCATE( tr(jpi,jpj,jpk,jptra,jpt)                                             ,       &  
          &      trc_i(jpi,jpj,jptra)  , trc_o(jpi,jpj,jptra)                          ,       &
          &      gtru (jpi,jpj,jptra)  , gtrv (jpi,jpj,jptra)                          ,       &
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcini.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcini.F90	(revision 10879)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcini.F90	(revision 10880)
@@ -26,4 +26,5 @@
    USE trcice          ! tracers in sea ice
    USE trcbc,   only : trc_bc_ini ! generalized Boundary Conditions
+   USE trcstp          ! for time level indices (to be initialised)
  
    IMPLICIT NONE
@@ -61,4 +62,10 @@
       CALL trc_nam       ! read passive tracers namelists
       CALL top_alloc()   ! allocate TOP arrays
+
+      ! Initialise time level indices
+      Nbb = 1; Nnn = 2; Naa = 3; Nrhs = Naa
+
+      ! Initialisation of temporary pointers (to be deleted after development finished)
+      CALL update_pointers_trc()
       !
       IF(.NOT.ln_trcdta )   ln_trc_ini(:) = .FALSE.
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcstp.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcstp.F90	(revision 10879)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcstp.F90	(revision 10880)
@@ -30,4 +30,10 @@
 
    PUBLIC   trc_stp    ! called by step
+   PUBLIC   update_pointers_trc ! called in initialisation
+
+   !!----------------------------------------------------------------------
+   !! time level indices
+   !!----------------------------------------------------------------------
+   INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs               !! used by trc_init
 
    LOGICAL  ::   llnew                   ! ???
@@ -100,5 +106,5 @@
                                    CALL trc_wri      ( kt )       ! output of passive tracers with iom I/O manager
                                    CALL trc_sms      ( kt )       ! tracers: sinks and sources
-                                   CALL trc_trp      ( kt )       ! transport of passive tracers
+                                   CALL trc_trp      ( kt, Nbb, Nnn, Nrhs, Naa )       ! transport of passive tracers
          IF( kt == nittrc000 ) THEN
             CALL iom_close( numrtr )       ! close input tracer restart file
@@ -125,4 +131,17 @@
    END SUBROUTINE trc_stp
 
+   SUBROUTINE update_pointers_trc
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE update_pointers_trc  ***
+      !!
+      !! ** Purpose :   Associate temporary pointer arrays.
+      !!                For IMMERSE development phase only - to be deleted
+      !!
+      !! ** Method  :
+      !!----------------------------------------------------------------------
+
+      trb => tr(:,:,:,:,Nbb); trn => tr(:,:,:,:,Nnn); tra => tr(:,:,:,:,Naa)
+
+   END SUBROUTINE update_pointers_trc
 
    SUBROUTINE trc_mean_qsr( kt )
