Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_ubs.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_ubs.F90	(revision 10799)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_ubs.F90	(revision 10802)
@@ -234,5 +234,5 @@
       ENDIF
       !                                         ! Control print
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' uu(:,:,:,ktlev1)s2 adv - Ua: ', mask1=umask,   &
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
          &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynvor.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynvor.F90	(revision 10799)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynvor.F90	(revision 10802)
@@ -878,5 +878,5 @@
       IF( ln_dynvor_msk )   CALL ctl_stop( 'dyn_vor_init:   masked vorticity is not currently not available')
 
-!!gm  this should be removed when choosing a uu(:,:,:,ktlev)ique strategy for fmask at the coast
+!!gm  this should be removed when choosing a unique strategy for fmask at the coast
       ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
       ! at angles with three ocean points and one land point
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 10799)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv.F90	(revision 10802)
@@ -75,5 +75,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv( kt )
+   SUBROUTINE tra_adv( kt, ktlev1, ktlev2, ktlev3, pts_rhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv  ***
@@ -81,7 +81,9 @@
       !! ** 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 (pu_rhs,pv_rhs) with the advection term following nadv
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt                       ! ocean time-step index
+      INTEGER, INTENT(in) ::   ktlev1, ktlev2, ktlev3   ! time level indices for source terms
+      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends
       !
       INTEGER ::   jk   ! dummy loop index
@@ -103,13 +105,13 @@
       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) )
+            zun(:,:,jk) = e2u  (:,:) * e3u(:,:,jk,ktlev2) * ( uu(:,:,jk,ktlev2) + usd(:,:,jk) )
+            zvn(:,:,jk) = e1v  (:,:) * e3v(:,:,jk,ktlev2) * ( vv(:,:,jk,ktlev2) + vsd(:,:,jk) )
+            zwn(:,:,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)
+            zun(:,:,jk) = e2u  (:,:) * e3u(:,:,jk,ktlev2) * uu(:,:,jk,ktlev2)               ! eulerian transport only
+            zvn(:,:,jk) = e1v  (:,:) * e3v(:,:,jk,ktlev2) * vv(:,:,jk,ktlev2)
+            zwn(:,:,jk) = e1e2t(:,:)                 * ww(:,:,jk)
          END DO
       ENDIF
@@ -139,6 +141,6 @@
       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_rhs(:,:,:,jp_tem)
+         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal)
       ENDIF
       !
@@ -146,13 +148,14 @@
       !
       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, ktlev2, 'TRA',         zun, zvn, zwn     , ts(:,:,:,:,ktlev2), pts_rhs, jpts, 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, ktlev1, ktlev2, ktlev3, 'TRA', r2dt, zun, zvn, zwn,                      &
+     &                         ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev2), pts_rhs, jpts, 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, ktlev2, 'TRA', r2dt, zun, zvn, zwn, ts(:,:,:,:,ktlev1),      pts_rhs, jpts        , 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, ktlev2, 'TRA', r2dt, zun, zvn, zwn, ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev2), pts_rhs, jpts        , 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, ktlev2, 'TRA', r2dt, zun, zvn, zwn, ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev2), pts_rhs, jpts                     )
       !
       END SELECT
@@ -160,6 +163,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_rhs(:,:,jk,jp_tem) - ztrdt(:,:,jk)
+            ztrds(:,:,jk) = pts_rhs(:,:,jk,jp_sal) - 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 10799)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_cen.F90	(revision 10802)
@@ -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, ktlev, cdtype, pu, pv, pwn,     &
+      &                                               pt, pt_rhs, kjpt, kn_cen_h, kn_cen_v ) 
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_cen  ***
@@ -59,5 +59,5 @@
       !!                = 4  ==>> 4th order COMPACT  scheme     -      -
       !!
-      !! ** Action : - update pta  with the now advective tracer trends
+      !! ** Action : - update pt_rhs  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)
@@ -65,11 +65,12 @@
       INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
       INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
+      INTEGER                              , INTENT(in   ) ::   ktlev           ! time level index for source terms
       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 
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pwn   ! 3 ocean velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt             ! now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs             ! tracer trend 
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -106,6 +107,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(ji,jj,jk) * ( pt(ji,jj,jk,jn) + pt(ji+1,jj  ,jk,jn) )
+                     zwy(ji,jj,jk) = 0.5_wp * pv(ji,jj,jk) * ( pt(ji,jj,jk,jn) + pt(ji  ,jj+1,jk,jn) )
                   END DO
                END DO
@@ -118,6 +119,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) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
+                     ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
                   END DO
                END DO
@@ -128,12 +129,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) + pt(ji+1,jj  ,jk,jn)   ! C2 interpolation of T at u- & v-points (x2)
+                     zC2t_v = pt(ji,jj,jk,jn) + pt(ji  ,jj+1,jk,jn)
                      !                                                  ! 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(ji,jj,jk) * zC4t_u
+                     zwy(ji,jj,jk) =  0.5_wp * pv(ji,jj,jk) * zC4t_v
                   END DO
                END DO
@@ -150,5 +151,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 * pwn(ji,jj,jk) * ( pt(ji,jj,jk,jn) + pt(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)
                   END DO
                END DO
@@ -156,5 +157,5 @@
             !
          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) , ztw )      ! ztw = interpolated value of T at w-point
             DO jk = 2, jpkm1
                DO jj = 2, jpjm1
@@ -171,9 +172,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) ) = pwn(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn) 
                   END DO
                END DO   
             ELSE                                   ! no ice-shelf cavities (only ocean surface)
-               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
+               zwz(:,:,1) = pwn(:,:,1) * pt(:,:,1,jn)
             ENDIF
          ENDIF
@@ -182,8 +183,8 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn)    &
+                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    &
                      &             - (  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,ktlev)
                END DO
             END DO
@@ -191,7 +192,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, pt(:,:,:,jn) )
+            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pv, pt(:,:,:,jn) )
+            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, pt(:,:,:,jn) )
          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 10799)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_fct.F90	(revision 10802)
@@ -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, ktlev1, ktlev2, ktlev3, cdtype, p2dt, pu, pv, pwn,       &
+      &                                                                pt_lev1, pt_lev2, pt_rhs, kjpt, kn_fct_h, kn_fct_v )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_fct  ***
@@ -65,9 +65,10 @@
       !!               - corrected flux (monotonic correction) 
       !!
-      !! ** Action : - update pta  with the now advective tracer trends
+      !! ** Action : - update pt_rhs  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   ) ::   ktlev1, ktlev2, ktlev3   ! time level indices for source terms
       INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
@@ -76,7 +77,7 @@
       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 
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pwn   ! 3 ocean velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev1, pt_lev2        ! before and now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs             ! tracer trend 
       !
       INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices  
@@ -125,10 +126,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(ji,jj,jk) + ABS( pu(ji,jj,jk) )
+                  zfm_ui = pu(ji,jj,jk) - ABS( pu(ji,jj,jk) )
+                  zfp_vj = pv(ji,jj,jk) + ABS( pv(ji,jj,jk) )
+                  zfm_vj = pv(ji,jj,jk) - ABS( pv(ji,jj,jk) )
+                  zwx(ji,jj,jk) = 0.5 * ( zfp_ui * pt_lev1(ji,jj,jk,jn) + zfm_ui * pt_lev1(ji+1,jj  ,jk,jn) )
+                  zwy(ji,jj,jk) = 0.5 * ( zfp_vj * pt_lev1(ji,jj,jk,jn) + zfm_vj * pt_lev1(ji  ,jj+1,jk,jn) )
                END DO
             END DO
@@ -140,5 +141,5 @@
                   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)
+                  zwz(ji,jj,jk) = 0.5 * ( zfp_wk * pt_lev1(ji,jj,jk,jn) + zfm_wk * pt_lev1(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)
                END DO
             END DO
@@ -148,9 +149,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) ) = pwn(ji,jj,mikt(ji,jj)) * pt_lev1(ji,jj,mikt(ji,jj),jn)   ! linear free surface 
                   END DO
                END DO   
             ELSE                             ! no cavities: only at the ocean surface
-               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
+               zwz(:,:,1) = pwn(:,:,1) * pt_lev1(:,:,1,jn)
             ENDIF
          ENDIF
@@ -164,6 +165,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_rhs(ji,jj,jk,jn) =                     pt_rhs(ji,jj,jk,jn) +        ztra   / e3t(ji,jj,jk,ktlev2) * tmask(ji,jj,jk)
+                  zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,ktlev1) * pt_lev1(ji,jj,jk,jn) + p2dt * ztra ) / e3t(ji,jj,jk,ktlev3) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -184,6 +185,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(ji,jj,jk) * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk)
+                     zwy(ji,jj,jk) = 0.5_wp * pv(ji,jj,jk) * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk)
                   END DO
                END DO
@@ -196,6 +197,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_lev2(ji+1,jj  ,jk,jn) - pt_lev2(ji,jj,jk,jn) ) * umask(ji,jj,jk)
+                     ztv(ji,jj,jk) = ( pt_lev2(ji  ,jj+1,jk,jn) - pt_lev2(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
                   END DO
                END DO
@@ -212,9 +213,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_lev2(ji,jj,jk,jn) + pt_lev2(ji+1,jj  ,jk,jn)   ! 2 x C2 interpolation of T at u- & v-points
+                     zC2t_v = pt_lev2(ji,jj,jk,jn) + pt_lev2(ji  ,jj+1,jk,jn)
                      !                                                  ! 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(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(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk)
                   END DO
                END DO
@@ -227,6 +228,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_lev2(ji+1,jj  ,jk,jn) - pt_lev2(ji,jj,jk,jn) ) * umask(ji,jj,jk)
+                     ztv(ji,jj,jk) = ( pt_lev2(ji  ,jj+1,jk,jn) - pt_lev2(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
                   END DO
                END DO
@@ -237,12 +238,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_lev2(ji,jj,jk,jn) + pt_lev2(ji+1,jj  ,jk,jn)   ! 2 x C2 interpolation of T at u- & v-points (x2)
+                     zC2t_v = pt_lev2(ji,jj,jk,jn) + pt_lev2(ji  ,jj+1,jk,jn)
                      !                                                  ! 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(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk)
+                     zwy(ji,jj,jk) =  0.5_wp * pv(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk)
                   END DO
                END DO
@@ -257,5 +258,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) =  (  pwn(ji,jj,jk) * 0.5_wp * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji,jj,jk-1,jn) )   &
                         &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk)
                   END DO
@@ -264,5 +265,5 @@
             !
          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_lev2(:,:,:,jn) , ztw )   ! zwt = COMPACT interpolation of T at w-point
             DO jk = 2, jpkm1
                DO jj = 2, jpjm1
@@ -282,5 +283,5 @@
          !        !==  monotonicity algorithm  ==!
          !
-         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
+         CALL nonosc( pt_lev1(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt, e3t(:,:,:,ktlev2) )
          !
          !        !==  final trend with corrected fluxes  ==!
@@ -289,8 +290,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_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) - (  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,ktlev2)
                END DO
             END DO
@@ -303,7 +304,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, pt_lev2(:,:,:,jn) )
+               CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pv, pt_lev2(:,:,:,jn) )
+               CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, pt_lev2(:,:,:,jn) )
             ENDIF
             !                             ! heat/salt transport
@@ -328,5 +329,5 @@
 
 
-   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
+   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt, pe3t )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE nonosc  ***
@@ -341,7 +342,7 @@
       !!       in-space based differencing for fluid
       !!----------------------------------------------------------------------
-      REAL(wp)                         , INTENT(in   ) ::   p2dt            ! tracer time-step
-      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field
-      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) ::   paa, pbb, pcc   ! monotonic fluxes in the 3 directions
+      REAL(wp)                         , INTENT(in   ) ::   p2dt             ! tracer time-step
+      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft, pe3t ! before & after field, now e3t field
+      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) ::   paa, pbb, pcc    ! monotonic fluxes in the 3 directions
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -392,5 +393,5 @@
 
                ! up & down beta terms
-               zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt
+               zbt = e1e2t(ji,jj) * pe3t(ji,jj,jk) / 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
@@ -634,5 +635,5 @@
       !!        The tri-diagonals matrix is given as input 3D arrays:   pD, pU, pL 
       !!        (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal).
-      !!        The solution is pta.
+      !!        The solution is pt_rhs.
       !!        The 3d array zwt is used as a work space array.
       !!----------------------------------------------------------------------
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 10799)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_mus.F90	(revision 10802)
@@ -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, ktlev, cdtype, p2dt, pu, pv, pwn,             &
+      &                                                     pt, pt_rhs, kjpt, 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_rhs  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)
@@ -75,11 +75,12 @@
       INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
       INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
+      INTEGER                              , INTENT(in   ) ::   ktlev           ! time level index for source terms
       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 
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pwn   ! 3 ocean velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt             ! before tracer field
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs             ! tracer trend 
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -134,6 +135,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) - pt(ji,jj,jk,jn) )
+                  zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( pt(ji,jj+1,jk,jn) - pt(ji,jj,jk,jn) )
                END DO
            END DO
@@ -172,17 +173,17 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! MUSCL fluxes
-                  z0u = SIGN( 0.5, pun(ji,jj,jk) )
+                  z0u = SIGN( 0.5, pu(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(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,ktlev)
+                  zzwx = pt(ji+1,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk)
+                  zzwy = pt(ji  ,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk)
+                  zwx(ji,jj,jk) = pu(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
                   !
-                  z0v = SIGN( 0.5, pvn(ji,jj,jk) )
+                  z0v = SIGN( 0.5, pv(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(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,ktlev)
+                  zzwx = pt(ji,jj+1,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk)
+                  zzwy = pt(ji,jj  ,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk)
+                  zwy(ji,jj,jk) = pv(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
                END DO
             END DO
@@ -193,7 +194,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_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) - ( 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,ktlev)
                END DO
            END DO
@@ -201,6 +202,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, pt(:,:,:,jn) )
+            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pv, pt(:,:,:,jn) )
          END IF
          !                                 ! "Poleward" heat and salt transports 
@@ -215,5 +216,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) - pt(:,:,jk,jn) )
          END DO
          !                                !-- Slopes of tracer
@@ -242,6 +243,6 @@
                   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  )
+                  zzwx = pt(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1)
+                  zzwy = pt(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)
                END DO 
@@ -252,9 +253,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) ) = pwn(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn)
                   END DO
                END DO   
             ELSE                                      ! no cavities: only at the ocean surface
-               zwx(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
+               zwx(:,:,1) = pwn(:,:,1) * pt(:,:,1,jn)
             ENDIF
          ENDIF
@@ -263,10 +264,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_rhs(ji,jj,jk,jn) =  pt_rhs(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev)
                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, pwn, pt(:,:,:,jn) )
          !
       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 10799)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_qck.F90	(revision 10802)
@@ -47,6 +47,6 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
-      &                                       ptb, ptn, pta, kjpt )
+   SUBROUTINE tra_adv_qck ( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pwn,      &
+      &                                                pt_lev1, pt_lev2, pt_rhs, kjpt )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_qck  ***
@@ -72,5 +72,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_lev2
       !!
       !!               The fluxes are bounded by the ULTIMATE limiter to
@@ -78,5 +78,5 @@
       !!            prevent the appearance of spurious numerical oscillations
       !!
-      !! ** Action : - update pta  with the now advective tracer trends
+      !! ** Action : - update pt_rhs  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)
@@ -86,10 +86,11 @@
       INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
       INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
+      INTEGER                              , INTENT(in   ) ::   ktlev           ! time level index for source terms
       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 
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pwn   ! 3 ocean velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev1, pt_lev2        ! before and now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs             ! tracer trend 
       !!----------------------------------------------------------------------
       !
@@ -108,25 +109,26 @@
       !
       !        ! 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, ktlev, cdtype, p2dt, pu, pt_lev1, pt_lev2, pt_rhs, kjpt ) 
+      CALL tra_adv_qck_j( kt, ktlev, cdtype, p2dt, pv, pt_lev1, pt_lev2, pt_rhs, kjpt ) 
 
       !        ! 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, ktlev, cdtype, pwn,         pt_lev2, pt_rhs, kjpt )
       !
    END SUBROUTINE tra_adv_qck
 
 
-   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun,                  &
-      &                                        ptb, ptn, pta, kjpt   )
+   SUBROUTINE tra_adv_qck_i( kt, ktlev, cdtype, p2dt, pu,                  &
+      &                                        pt_lev1, pt_lev2, pt_rhs, kjpt   )
       !!----------------------------------------------------------------------
       !!
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   ktlev           ! time level index for source terms
       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 
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu        ! i-velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev1, pt_lev2   ! before and now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs        ! tracer trend 
       !!
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -145,6 +147,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_lev1(ji-1,jj,jk,jn)        ! Upstream   in the x-direction for the tracer
+                  zfd(ji,jj,jk) = pt_lev1(ji+1,jj,jk,jn)        ! Downstream in the x-direction for the tracer
                END DO
             END DO
@@ -158,5 +160,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(ji,jj,jk) )   ! if pu > 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 +169,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(ji,jj,jk) )   ! if pu > 0 : zdir = 1 otherwise zdir = 0 
+                  zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,ktlev)
+                  zwx(ji,jj,jk)  = ABS( pu(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
+                  zfc(ji,jj,jk)  = zdir * pt_lev1(ji  ,jj,jk,jn) + ( 1. - zdir ) * pt_lev1(ji+1,jj,jk,jn)  ! FC in the x-direction for T
+                  zfd(ji,jj,jk)  = zdir * pt_lev1(ji+1,jj,jk,jn) + ( 1. - zdir ) * pt_lev1(ji  ,jj,jk,jn)  ! FD in the x-direction for T
                END DO
             END DO
@@ -197,10 +199,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(ji,jj,jk) )   ! if pu > 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(ji,jj,jk)
                END DO
             END DO
@@ -213,14 +215,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,ktlev)
                   ! 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_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + 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, pt_lev2(:,:,:,jn) )
          !
       END DO
@@ -229,16 +231,17 @@
 
 
-   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                &
-      &                                        ptb, ptn, pta, kjpt )
+   SUBROUTINE tra_adv_qck_j( kt, ktlev, cdtype, p2dt, pv,                &
+      &                                        pt_lev1, pt_lev2, pt_rhs, kjpt )
       !!----------------------------------------------------------------------
       !!
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   ktlev           ! time level index for source terms
       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 
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pv        ! j-velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev1, pt_lev2   ! before and now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs        ! tracer trend 
       !!
       INTEGER  :: ji, jj, jk, jn                ! dummy loop indices
@@ -259,7 +262,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_lev1(ji,jj-1,jk,jn)
                   ! Downstream in the x-direction for the tracer
-                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
+                  zfd(ji,jj,jk) = pt_lev1(ji,jj+1,jk,jn)
                END DO
             END DO
@@ -275,5 +278,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(ji,jj,jk) )   ! if pu > 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 +287,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(ji,jj,jk) )   ! if pu > 0 : zdir = 1 otherwise zdir = 0 
+                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,ktlev)
+                  zwy(ji,jj,jk)  = ABS( pv(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
+                  zfc(ji,jj,jk)  = zdir * pt_lev1(ji,jj  ,jk,jn) + ( 1. - zdir ) * pt_lev1(ji,jj+1,jk,jn)  ! FC in the x-direction for T
+                  zfd(ji,jj,jk)  = zdir * pt_lev1(ji,jj+1,jk,jn) + ( 1. - zdir ) * pt_lev1(ji,jj  ,jk,jn)  ! FD in the x-direction for T
                END DO
             END DO
@@ -314,10 +317,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(ji,jj,jk) )   ! if pu > 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(ji,jj,jk)
                END DO
             END DO
@@ -330,14 +333,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,ktlev)
                   ! 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_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + 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, pt_lev2(:,:,:,jn) )
          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
          IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
@@ -348,15 +351,16 @@
 
 
-   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           &
-     &                                    ptn, pta, kjpt )
+   SUBROUTINE tra_adv_cen2_k( kt, ktlev, cdtype, pwn,           &
+     &                                    pt_lev2, pt_rhs, kjpt )
       !!----------------------------------------------------------------------
       !!
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   ktlev    ! time level index for source terms
       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 
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev2      ! before and now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs      ! tracer trend 
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -374,5 +378,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 * pwn(ji,jj,jk) * ( pt_lev2(ji,jj,jk-1,jn) + pt_lev2(ji,jj,jk,jn) ) * wmask(ji,jj,jk)
                END DO
             END DO
@@ -382,9 +386,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) ) = pwn(ji,jj,mikt(ji,jj)) * pt_lev2(ji,jj,mikt(ji,jj),jn)   ! linear free surface 
                   END DO
                END DO   
             ELSE                                   ! no ocean cavities (only ocean surface)
-               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
+               zwz(:,:,1) = pwn(:,:,1) * pt_lev2(:,:,1,jn)
             ENDIF
          ENDIF
@@ -393,11 +397,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_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
+                     &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev)
                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, pwn, pt_lev2(:,:,:,jn) )
          !
       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 10799)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_ubs.F90	(revision 10802)
@@ -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, ktlev, cdtype, p2dt, pu, pv, pwn,          &
+      &                                                pt_lev1, pt_lev2, pt_rhs, kjpt, 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  ) ,ktlev)   if uu(i,ktlev) >= 0
       !!          ztu = !  or 
-      !!                !  e2u e3u un ( mi(Tn) - zltu(i+1) )   if un(i) < 0
+      !!                !  e2u e3u uu ( mi(Tn) - zltu(i+1) ,ktlev)   if uu(i,ktlev) < 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_rhs  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)
@@ -86,11 +86,12 @@
       INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
       INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
+      INTEGER                              , INTENT(in   ) ::   ktlev           ! time level index for source terms
       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 
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pwn   ! 3 ocean transport components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev1, pt_lev2        ! before and now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs             ! tracer trend 
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -126,13 +127,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,ktlev) * umask(ji,jj,jk)
+                  zeev = e1_e2v(ji,jj) * e3v(ji,jj,jk,ktlev) * vmask(ji,jj,jk)
+                  ztu(ji,jj,jk) = zeeu * ( pt_lev1(ji+1,jj  ,jk,jn) - pt_lev1(ji,jj,jk,jn) )
+                  ztv(ji,jj,jk) = zeev * ( pt_lev1(ji  ,jj+1,jk,jn) - pt_lev1(ji,jj,jk,jn) )
                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,ktlev) )
                   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 +147,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(ji,jj,jk) + ABS( pu(ji,jj,jk) )      ! upstream transport (x2)
+                  zfm_ui = pu(ji,jj,jk) - ABS( pu(ji,jj,jk) )
+                  zfp_vj = pv(ji,jj,jk) + ABS( pv(ji,jj,jk) )
+                  zfm_vj = pv(ji,jj,jk) - ABS( pv(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(ji,jj,jk) * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji+1,jj  ,jk,jn) )
+                  zcenvt = pv(ji,jj,jk) * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji  ,jj+1,jk,jn) )
                   !                                                  ! 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 +161,12 @@
          END DO         
          !
-         zltu(:,:,:) = pta(:,:,:,jn)      ! store the initial trends before its update
+         zltu(:,:,:) = pt_rhs(:,:,:,jn)      ! 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_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)                        &
                      &             - (  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,ktlev)
                END DO
             END DO
@@ -173,10 +174,10 @@
          END DO
          !
-         zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:)    ! Horizontal advective trend used in vertical 2nd order FCT case
+         zltu(:,:,:) = pt_rhs(:,:,:,jn) - 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, pt_lev2(:,:,:,jn) )
+             CALL trd_tra( kt, cdtype, jn, jptra_yad, ztv, pv, pt_lev2(:,:,:,jn) )
          END IF
          !     
@@ -193,5 +194,5 @@
          CASE(  2  )                   ! 2nd order FCT 
             !         
-            IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag.
+            IF( l_trd )   zltv(:,:,:) = pt_rhs(:,:,:,jn)          ! store pt_rhs if trend diag.
             !
             !                          !*  upstream advection with initial mass fluxes & intermediate update  ==!
@@ -201,5 +202,5 @@
                      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)
+                     ztw(ji,jj,jk) = 0.5_wp * (  zfp_wk * pt_lev1(ji,jj,jk,jn) + zfm_wk * pt_lev1(ji,jj,jk-1,jn)  ) * wmask(ji,jj,jk)
                   END DO
                END DO
@@ -209,9 +210,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) ) = pwn(ji,jj,mikt(ji,jj)) * pt_lev1(ji,jj,mikt(ji,jj),jn)   ! linear free surface 
                      END DO
                   END DO   
                ELSE                                ! no cavities: only at the ocean surface
-                  ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
+                  ztw(:,:,1) = pwn(:,:,1) * pt_lev1(:,:,1,jn)
                ENDIF
             ENDIF
@@ -220,7 +221,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,ktlev)
+                     pt_rhs(ji,jj,jk,jn) =   pt_rhs(ji,jj,jk,jn) +  ztak 
+                     zti(ji,jj,jk)    = ( pt_lev1(ji,jj,jk,jn) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
                   END DO
                END DO
@@ -232,5 +233,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 * pwn(ji,jj,jk) * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji,jj,jk-1,jn) )   &
                         &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk)
                   END DO
@@ -240,8 +241,8 @@
             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( pt_lev1(:,:,:,jn), ztw, zti, p2dt, e3t(:,:,:,ktlev) )      !  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_lev2(:,:,:,jn) , ztw )         ! 4th order compact interpolation of T at w-point
             DO jk = 2, jpkm1
                DO jj = 2, jpjm1
@@ -251,5 +252,5 @@
                END DO
             END DO
-            IF( ln_linssh )   ztw(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)     !!gm ISF & 4th COMPACT doesn't work
+            IF( ln_linssh )   ztw(:,:, 1 ) = pwn(:,:,1) * pt_lev2(:,:,1,jn)     !!gm ISF & 4th COMPACT doesn't work
             !
          END SELECT
@@ -258,5 +259,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_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev)
                END DO
             END DO
@@ -264,10 +265,10 @@
          !
          IF( l_trd )  THEN       ! vertical advective trend diagnostics
-            DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
+            DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + pt_lev2.dk[w])
                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_rhs(ji,jj,jk,jn) - zltv(ji,jj,jk)                          &
+                        &           + pt_lev2(ji,jj,jk,jn) * (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  )   &
+                        &                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev)
                   END DO
                END DO
@@ -281,5 +282,5 @@
 
 
-   SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt )
+   SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt, pe3t )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE nonosc_z  ***
@@ -296,4 +297,5 @@
       REAL(wp), INTENT(in   )                          ::   p2dt   ! tracer time-step
       REAL(wp),                DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field
+      REAL(wp), INTENT(in   ), DIMENSION (jpi,jpj,jpk) ::   pe3t   ! now cell thickness field
       REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paft   ! after field
       REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pcc    ! monotonic flux in the k direction
@@ -352,5 +354,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) * pe3t(ji,jj,jk) / 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/oce.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/oce.F90	(revision 10799)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/oce.F90	(revision 10802)
@@ -23,5 +23,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           wi             !: vertical vel. (adaptive-implicit) [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           hdivn          !: horizontal divergence        [s-1]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   tsb  ,  tsn   , tsa    !: 4D T-S fields                  [Celsius,psu] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:), TARGET :: ts             !: 4D T-S fields                  [Celsius,psu] 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   rab_b,  rab_n          !: thermal/haline expansion coef. [Celsius-1,psu-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   rn2b ,  rn2            !: brunt-vaisala frequency**2     [s-2]
@@ -72,4 +72,6 @@
    REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::   vb   ,  vn    , va     !: j-horizontal velocity        [m/s]
    REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::           wn             !: k-vertical   velocity        [m/s]
+   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:,:) ::   tsb  ,  tsn   , tsa    !: 4D T-S fields                [Celsius,psu]         
+   !! TEMPORARY POINTERS FOR DEVELOPMENT ONLY
 
    !!----------------------------------------------------------------------
@@ -90,5 +92,5 @@
       ALLOCATE( uu   (jpi,jpj,jpk, jpt) , vv   (jpi,jpj,jpk,jpt)  ,                             &
          &      ww   (jpi,jpj,jpk)      , hdivn(jpi,jpj,jpk)      ,                             &
-         &      tsb  (jpi,jpj,jpk,jpts) , tsn  (jpi,jpj,jpk,jpts) , tsa(jpi,jpj,jpk,jpts) ,     &
+         &      ts  (jpi,jpj,jpk,jpts,jpt) ,                                                    &
          &      rab_b(jpi,jpj,jpk,jpts) , rab_n(jpi,jpj,jpk,jpts) ,                             &
          &      rn2b (jpi,jpj,jpk)      , rn2  (jpi,jpj,jpk)      ,                             &
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 10799)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10802)
@@ -247,5 +247,5 @@
                &         CALL Agrif_Sponge_tra        ! tracers sponge
 #endif
-                         CALL tra_adv       ( kstp )  ! horizontal & vertical advection
+                         CALL tra_adv       ( kstp, Nm1, Nnn, Np1, ts(:,:,:,:,Nrhs) )  ! horizontal & vertical advection
       IF( ln_zdfosm  )   CALL tra_osm       ( kstp )  ! OSMOSIS non-local tracer fluxes
       IF( lrst_oce .AND. ln_zdfosm ) &
@@ -344,4 +344,6 @@
       wn => ww(:,:,:)
 
+      tsb => ts(:,:,:,:,Nm1); tsn => ts(:,:,:,:,Nnn); tsa => ts(:,:,:,:,Np1)
+
       e3t_b => e3t(:,:,:,Nm1); e3t_n => e3t(:,:,:,Nnn); e3t_a => e3t(:,:,:,Np1)
       e3u_b => e3u(:,:,:,Nm1); e3u_n => e3u(:,:,:,Nnn); e3u_a => e3u(:,:,:,Np1)
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 10799)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcadv.F90	(revision 10802)
@@ -68,5 +68,5 @@
 CONTAINS
 
-   SUBROUTINE trc_adv( kt )
+   SUBROUTINE trc_adv( kt, ktlev1, ktlev2, ktlev3 )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trc_adv  ***
@@ -77,4 +77,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   ktlev1, ktlev2, ktlev3   ! time level indices for source terms
       !
       INTEGER ::   jk   ! dummy loop index
@@ -123,13 +124,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, ktlev2, 'TRC',          zun, zvn, zwn     , trn, tra, jptra, 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, ktlev1, ktlev2, ktlev3, 'TRC', r2dttrc, zun, zvn, zwn, trb, trn, tra, jptra, 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, ktlev2, 'TRC', r2dttrc, zun, zvn, zwn, trb,      tra, jptra        , 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, ktlev2, 'TRC', r2dttrc, zun, zvn, zwn, trb, trn, tra, jptra          , 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, ktlev2, 'TRC', r2dttrc, zun, zvn, zwn, trb, trn, tra, jptra                     )
       !
       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 10799)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90	(revision 10802)
@@ -64,5 +64,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, Nm1, Nnn, Np1 )      ! horizontal & vertical advection 
          !                                                         ! Partial top/bottom cell: GRADh( trb )  
          IF( ln_zps ) THEN
