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 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/dom_oce.F90	(revision 10806)
@@ -154,11 +154,10 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:)   ::  gde3w  
    
-   !                                                      !  reference depths of water column
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0  !: t-depth              [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0  !: u-depth              [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hv_0  !: v-depth              [m]
-   !                                                      !  time-dependent depths of water column
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:) ::  hu, hv, r1_hu, r1_hv  
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:)   ::  ht
+   !                                                      !  ref. ! before  !   now   !  after  !
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0            ,    ht_n             !: t-depth              [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
 
    !! TEMPORARY POINTERS FOR DEVELOPMENT ONLY
@@ -173,9 +172,4 @@
    REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::    gdepw_b , gdepw_n           !: w- depth              [m]
    REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::    gde3w_n                     !: w- depth (sum of e3w) [m]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:)   ::    ht_n                        !: t-depth              [m]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:)   ::    hu_b ,    hu_n ,    hu_a    !: u-depth              [m]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:)   ::    hv_b ,    hv_n ,    hv_a    !: v-depth              [m]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:)   :: r1_hu_b , r1_hu_n , r1_hu_a    !: inverse of u-depth [1/m]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:)   :: r1_hv_b , r1_hv_n , r1_hv_a    !: inverse of v-depth [1/m]
    !! TEMPORARY POINTERS FOR DEVELOPMENT ONLY
 
@@ -295,7 +289,8 @@
          &      e3f  (jpi,jpj,jpk),     STAT=ierr(5) )                       
          !
-      ALLOCATE( ht_0(jpi,jpj) ,    hu_0(jpi,jpj)     ,    hv_0(jpi,jpj) ,                                           &
-         &      ht  (jpi,jpj) ,    hu  (jpi,jpj,jpt) ,    hv  (jpi,jpj,jpt) ,                                       &
-         &                      r1_hu  (jpi,jpj,jpt) , r1_hv  (jpi,jpj,jpt) , STAT=ierr(6) )
+      ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) ,                                           &
+         &                      hu_b(jpi,jpj) , hv_b(jpi,jpj) , r1_hu_b(jpi,jpj) , r1_hv_b(jpi,jpj) ,     &
+         &      ht_n(jpi,jpj) , hu_n(jpi,jpj) , hv_n(jpi,jpj) , r1_hu_n(jpi,jpj) , r1_hv_n(jpi,jpj) ,     &
+         &                      hu_a(jpi,jpj) , hv_a(jpi,jpj) , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) , STAT=ierr(6)  )
          !
          !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf.F90	(revision 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf.F90	(revision 10806)
@@ -43,5 +43,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_ldf( kt )
+   SUBROUTINE dyn_ldf( kt, ktlev1, ktlev2, pu_rhs, pv_rhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_ldf  ***
@@ -49,5 +49,7 @@
       !! ** Purpose :   compute the lateral ocean dynamics physics.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   kt               ! ocean time-step index
+      INTEGER, INTENT(in) ::   ktlev1, ktlev2   ! time level index for source terms
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv
@@ -58,19 +60,19 @@
       IF( l_trddyn )   THEN                      ! temporary save of momentum trends
          ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
-         ztrdu(:,:,:) = ua(:,:,:) 
-         ztrdv(:,:,:) = va(:,:,:) 
+         ztrdu(:,:,:) = pu_rhs(:,:,:) 
+         ztrdv(:,:,:) = pv_rhs(:,:,:) 
       ENDIF
 
       SELECT CASE ( nldf_dyn )                   ! compute lateral mixing trend and add it to the general trend
       !
-      CASE ( np_lap   )    ;   CALL dyn_ldf_lap( kt, ub, vb, ua, va, 1 )      ! iso-level    laplacian
+      CASE ( np_lap   )    ;   CALL dyn_ldf_lap( kt, ktlev1, ktlev2, uu(:,:,:,ktlev1), vv(:,:,:,ktlev1), pu_rhs, pv_rhs, 1 )      ! iso-level    laplacian
       CASE ( np_lap_i )    ;   CALL dyn_ldf_iso( kt )                         ! rotated      laplacian
-      CASE ( np_blp   )    ;   CALL dyn_ldf_blp( kt, ub, vb, ua, va    )      ! iso-level bi-laplacian
+      CASE ( np_blp   )    ;   CALL dyn_ldf_blp( kt, ktlev1, ktlev2, uu(:,:,:,ktlev1), vv(:,:,:,ktlev1), pu_rhs, pv_rhs    )      ! iso-level bi-laplacian
       !
       END SELECT
 
       IF( l_trddyn ) THEN                        ! save the horizontal diffusive trends for further diagnostics
-         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
-         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
+         ztrdu(:,:,:) = pu_rhs(:,:,:) - ztrdu(:,:,:)
+         ztrdv(:,:,:) = pv_rhs(:,:,:) - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_ldf, kt )
          DEALLOCATE ( ztrdu , ztrdv )
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_lap_blp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_lap_blp.F90	(revision 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_lap_blp.F90	(revision 10806)
@@ -35,5 +35,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_ldf_lap( kt, pub, pvb, pua, pva, kpass )
+   SUBROUTINE dyn_ldf_lap( kt, ktlev1, ktlev2, pu, pv, pu_rhs, pva_rhs, kpass )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE dyn_ldf_lap  ***
@@ -45,10 +45,11 @@
       !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) 
       !!
-      !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb.
+      !! ** Action : - pu_rhs, pva_rhs increased by the harmonic operator applied on pu, pv.
       !!----------------------------------------------------------------------
-      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                         , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                         , INTENT(in   ) ::   ktlev1, ktlev2  ! time level index for scale factors
       INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity  [m/s]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! velocity trend   [m/s2]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv   ! before velocity  [m/s]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pva_rhs   ! velocity trend   [m/s2]
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -76,12 +77,12 @@
 !!gm open question here : e3f  at before or now ?    probably now...
 !!gm note that ahmf has already been multiplied by fmask
-               zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &
-                  &     * (  e2v(ji  ,jj-1) * pvb(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk)  &
-                  &        - e1u(ji-1,jj  ) * pub(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk)  )
+               zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &
+                  &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  &
+                  &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  )
                !                                      ! ahm * div        (computed from 2 to jpi/jpj)
 !!gm note that ahmt has already been multiplied by tmask
-               zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk)                                         &
-                  &     * (  e2u(ji,jj)*e3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk)  &
-                  &        + e1v(ji,jj)*e3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk)  )
+               zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev1)                                         &
+                  &     * (  e2u(ji,jj)*e3u(ji,jj,jk,ktlev1) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,ktlev1) * pu(ji-1,jj,jk)  &
+                  &        + e1v(ji,jj)*e3v(ji,jj,jk,ktlev1) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,ktlev1) * pv(ji,jj-1,jk)  )
             END DO  
          END DO  
@@ -89,10 +90,10 @@
          DO jj = 2, jpjm1                             ! - curl( curl) + grad( div )
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               pua(ji,jj,jk) = pua(ji,jj,jk) + zsign * (                                                 &
-                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk)   &
+               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * (                                                 &
+                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,ktlev2)   &
                   &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     )
                   !
-               pva(ji,jj,jk) = pva(ji,jj,jk) + zsign * (                                                 &
-                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk)   &
+               pva_rhs(ji,jj,jk) = pva_rhs(ji,jj,jk) + zsign * (                                                 &
+                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,ktlev2)   &
                   &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     )
             END DO
@@ -105,5 +106,5 @@
 
 
-   SUBROUTINE dyn_ldf_blp( kt, pub, pvb, pua, pva )
+   SUBROUTINE dyn_ldf_blp( kt, ktlev1, ktlev2, pu, pv, pu_rhs, pva_rhs )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE dyn_ldf_blp  ***
@@ -116,9 +117,10 @@
       !!      It is computed by two successive calls to dyn_ldf_lap routine
       !!
-      !! ** Action :   pta   updated with the before rotated bilaplacian diffusion
+      !! ** Action :   pt_rhs   updated with the before rotated bilaplacian diffusion
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! momentum trend
+      INTEGER                         , INTENT(in   ) ::   ktlev1, ktlev2  ! time level index for scale factors
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv   ! before velocity fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pva_rhs   ! momentum trend
       !
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point
@@ -134,9 +136,9 @@
       zvlap(:,:,:) = 0._wp
       !
-      CALL dyn_ldf_lap( kt, pub, pvb, zulap, zvlap, 1 )   ! rotated laplacian applied to ptb (output in zlap)
+      CALL dyn_ldf_lap( kt, ktlev1, ktlev2, pu, pv, zulap, zvlap, 1 )   ! rotated laplacian applied to pt (output in zlap)
       !
       CALL lbc_lnk_multi( 'dynldf_lap_blp', zulap, 'U', -1., zvlap, 'V', -1. )             ! Lateral boundary conditions
       !
-      CALL dyn_ldf_lap( kt, zulap, zvlap, pua, pva, 2 )   ! rotated laplacian applied to zlap (output in pta)
+      CALL dyn_ldf_lap( kt, ktlev1, ktlev2, zulap, zvlap, pu_rhs, pva_rhs, 2 )   ! rotated laplacian applied to zlap (output in pt_rhs)
       !
    END SUBROUTINE dyn_ldf_blp
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 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynvor.F90	(revision 10806)
@@ -119,6 +119,6 @@
          ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )
          !
-         ztrdu(:,:,:) = ua(:,:,:)            !* planetary vorticity trend (including Stokes-Coriolis force)
-         ztrdv(:,:,:) = va(:,:,:)
+         ztrdu(:,:,:) = pu_rhs(:,:,:)            !* planetary vorticity trend (including Stokes-Coriolis force)
+         ztrdv(:,:,:) = pv_rhs(:,:,:)
          SELECT CASE( nvor_scheme )
          CASE( np_ENS )           ;   CALL vor_ens( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! enstrophy conserving scheme
@@ -133,11 +133,11 @@
             IF( ln_stcor )            CALL vor_een( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! add the Stokes-Coriolis trend
          END SELECT
-         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
-         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
+         ztrdu(:,:,:) = pu_rhs(:,:,:) - ztrdu(:,:,:)
+         ztrdv(:,:,:) = pv_rhs(:,:,:) - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          !
          IF( n_dynadv /= np_LIN_dyn ) THEN   !* relative vorticity or metric trend (only in non-linear case)
-            ztrdu(:,:,:) = ua(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:)
+            ztrdu(:,:,:) = pu_rhs(:,:,:)
+            ztrdv(:,:,:) = pv_rhs(:,:,:)
             SELECT CASE( nvor_scheme )
             CASE( np_ENT )           ;   CALL vor_enT( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )  ! energy conserving scheme (T-pts)
@@ -147,6 +147,6 @@
             CASE( np_EEN )           ;   CALL vor_een( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )  ! energy & enstrophy scheme
             END SELECT
-            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
+            ztrdu(:,:,:) = pu_rhs(:,:,:) - ztrdu(:,:,:)
+            ztrdv(:,:,:) = pv_rhs(:,:,:) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
          ENDIF
@@ -410,5 +410,5 @@
 
          IF( ln_sco ) THEN
-            zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)
+            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
             zwx(:,:) = e2u(:,:) * e3u(:,:,jk,ktlev) * pu(:,:,jk)
             zwy(:,:) = e1v(:,:) * e3v(:,:,jk,ktlev) * pv(:,:,jk)
@@ -518,5 +518,5 @@
          !
          IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
-            zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)
+            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
             zwx(:,:) = e2u(:,:) * e3u(:,:,jk,ktlev) * pu(:,:,jk)
             zwy(:,:) = e1v(:,:) * e3v(:,:,jk,ktlev) * pv(:,:,jk)
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 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv.F90	(revision 10806)
@@ -75,5 +75,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv( kt, ktlev1, ktlev2, ktlev3, pts_rhs )
+   SUBROUTINE tra_adv( kt, ktlev1, ktlev2, ktlev3, kt2lev, pts_rhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv  ***
@@ -84,5 +84,6 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt                       ! ocean time-step index
-      INTEGER, INTENT(in) ::   ktlev1, ktlev2, ktlev3   ! time level indices for source terms
+      INTEGER, INTENT(in) ::   ktlev1, ktlev2, ktlev3   ! time level indices for 3-time-level source terms
+      INTEGER, INTENT(in) ::   kt2lev                   ! time level index for 2-time-level source terms
       REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends
       !
@@ -153,5 +154,5 @@
      &                         ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev2), pts_rhs, jpts, nn_fct_h, nn_fct_v )
       CASE ( np_MUS )                                 ! MUSCL
-         CALL tra_adv_mus    ( kt, nit000, ktlev2, 'TRA', r2dt, zun, zvn, zwn, ts(:,:,:,:,ktlev1),      pts_rhs, jpts        , ln_mus_ups ) 
+         CALL tra_adv_mus    ( kt, nit000, ktlev2, kt2lev, 'TRA', r2dt, zun, zvn, zwn, ts(:,:,:,:,ktlev1),      pts_rhs, jpts        , ln_mus_ups ) 
       CASE ( np_UBS )                                 ! UBS
          CALL tra_adv_ubs    ( kt, nit000, ktlev2, 'TRA', r2dt, zun, zvn, zwn, ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev2), pts_rhs, jpts        , nn_ubs_v   )
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 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_cen.F90	(revision 10806)
@@ -44,5 +44,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_cen( kt, kit000, ktlev, cdtype, pu, pv, pwn,     &
+   SUBROUTINE tra_adv_cen( kt, kit000, ktlev, cdtype, pu, pv, pw,     &
       &                                               pt, pt_rhs, kjpt, kn_cen_h, kn_cen_v ) 
       !!----------------------------------------------------------------------
@@ -70,5 +70,5 @@
       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, pv, pwn   ! 3 ocean velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pw   ! 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 
@@ -151,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) * ( pt(ji,jj,jk,jn) + pt(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)
+                     zwz(ji,jj,jk) = 0.5 * pw(ji,jj,jk) * ( pt(ji,jj,jk,jn) + pt(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)
                   END DO
                END DO
@@ -161,5 +161,5 @@
                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) = pw(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
                   END DO
                END DO
@@ -172,9 +172,9 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn) 
+                     zwz(ji,jj, mikt(ji,jj) ) = pw(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) * pt(:,:,1,jn)
+               zwz(:,:,1) = pw(:,:,1) * pt(:,:,1,jn)
             ENDIF
          ENDIF
@@ -194,5 +194,5 @@
             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) )
+            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pw, 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 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_fct.F90	(revision 10806)
@@ -52,5 +52,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_fct( kt, kit000, ktlev1, ktlev2, ktlev3, cdtype, p2dt, pu, pv, pwn,       &
+   SUBROUTINE tra_adv_fct( kt, kit000, ktlev1, ktlev2, ktlev3, cdtype, p2dt, pu, pv, pw,       &
       &                                                                pt_lev1, pt_lev2, pt_rhs, kjpt, kn_fct_h, kn_fct_v )
       !!----------------------------------------------------------------------
@@ -77,5 +77,5 @@
       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, pv, pwn   ! 3 ocean velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pw   ! 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 
@@ -139,6 +139,6 @@
             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) )
+                  zfp_wk = pw(ji,jj,jk) + ABS( pw(ji,jj,jk) )
+                  zfm_wk = pw(ji,jj,jk) - ABS( pw(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
@@ -149,9 +149,9 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * pt_lev1(ji,jj,mikt(ji,jj),jn)   ! linear free surface 
+                     zwz(ji,jj, mikt(ji,jj) ) = pw(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) * pt_lev1(:,:,1,jn)
+               zwz(:,:,1) = pw(:,:,1) * pt_lev1(:,:,1,jn)
             ENDIF
          ENDIF
@@ -258,5 +258,5 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     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) =  (  pw(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
@@ -269,5 +269,5 @@
                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) = ( pw(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
                   END DO
                END DO
@@ -306,5 +306,5 @@
                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) )
+               CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pw, pt_lev2(:,:,:,jn) )
             ENDIF
             !                             ! heat/salt transport
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 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_mus.F90	(revision 10806)
@@ -54,6 +54,6 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_mus( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pwn,             &
-      &                                                     pt, pt_rhs, kjpt, ld_msc_ups )
+   SUBROUTINE tra_adv_mus( kt, kit000, ktlev, kt2lev, cdtype, p2dt, pu, pv, pw,             &
+      &                                                             pt, pt_rhs, kjpt, ld_msc_ups )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE tra_adv_mus  ***
@@ -75,10 +75,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
+      INTEGER                              , INTENT(in   ) ::   ktlev           ! time level index for 3-time-level source terms
+      INTEGER                              , INTENT(in   ) ::   kt2lev          ! time level index for 2-time-level 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   ) ::   pu, pv, pwn   ! 3 ocean velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pw   ! 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 
@@ -240,10 +241,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, pw(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)
+                  zw  = z0w - 0.5 * pw(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,kt2lev)
                   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)
+                  zwx(ji,jj,jk+1) = pw(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)
                END DO 
             END DO
@@ -253,9 +254,9 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zwx(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn)
+                     zwx(ji,jj, mikt(ji,jj) ) = pw(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) * pt(:,:,1,jn)
+               zwx(:,:,1) = pw(:,:,1) * pt(:,:,1,jn)
             ENDIF
          ENDIF
@@ -269,5 +270,5 @@
          END DO
          !                                ! send trends for diagnostic
-         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, pt(:,:,:,jn) )
+         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pw, 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 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_qck.F90	(revision 10806)
@@ -47,5 +47,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_qck ( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pwn,      &
+   SUBROUTINE tra_adv_qck ( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pw,      &
       &                                                pt_lev1, pt_lev2, pt_rhs, kjpt )
       !!----------------------------------------------------------------------
@@ -90,5 +90,5 @@
       INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
       REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pwn   ! 3 ocean velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pw   ! 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 
@@ -113,5 +113,5 @@
 
       !        ! vertical fluxes are computed with the 2nd order centered scheme
-      CALL tra_adv_cen2_k( kt, ktlev, cdtype, pwn,         pt_lev2, pt_rhs, kjpt )
+      CALL tra_adv_cen2_k( kt, ktlev, cdtype, pw,         pt_lev2, pt_rhs, kjpt )
       !
    END SUBROUTINE tra_adv_qck
@@ -351,5 +351,5 @@
 
 
-   SUBROUTINE tra_adv_cen2_k( kt, ktlev, cdtype, pwn,           &
+   SUBROUTINE tra_adv_cen2_k( kt, ktlev, cdtype, pw,           &
      &                                    pt_lev2, pt_rhs, kjpt )
       !!----------------------------------------------------------------------
@@ -360,5 +360,5 @@
       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     ), INTENT(in   ) ::   pw      ! vertical velocity 
       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 
@@ -378,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) * ( pt_lev2(ji,jj,jk-1,jn) + pt_lev2(ji,jj,jk,jn) ) * wmask(ji,jj,jk)
+                  zwz(ji,jj,jk) = 0.5 * pw(ji,jj,jk) * ( pt_lev2(ji,jj,jk-1,jn) + pt_lev2(ji,jj,jk,jn) ) * wmask(ji,jj,jk)
                END DO
             END DO
@@ -386,9 +386,9 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * pt_lev2(ji,jj,mikt(ji,jj),jn)   ! linear free surface 
+                     zwz(ji,jj, mikt(ji,jj) ) = pw(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) * pt_lev2(:,:,1,jn)
+               zwz(:,:,1) = pw(:,:,1) * pt_lev2(:,:,1,jn)
             ENDIF
          ENDIF
@@ -403,5 +403,5 @@
          END DO
          !                                 ! Send trends for diagnostic
-         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, pt_lev2(:,:,:,jn) )
+         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pw, 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 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_ubs.F90	(revision 10806)
@@ -46,5 +46,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_ubs( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pwn,          &
+   SUBROUTINE tra_adv_ubs( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pw,          &
       &                                                pt_lev1, pt_lev2, pt_rhs, kjpt, kn_ubs_v )
       !!----------------------------------------------------------------------
@@ -91,5 +91,5 @@
       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, pv, pwn   ! 3 ocean transport components
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pw   ! 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 
@@ -200,6 +200,6 @@
                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) )
+                     zfp_wk = pw(ji,jj,jk) + ABS( pw(ji,jj,jk) )
+                     zfm_wk = pw(ji,jj,jk) - ABS( pw(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
@@ -210,9 +210,9 @@
                   DO jj = 1, jpj
                      DO ji = 1, jpi
-                        ztw(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * pt_lev1(ji,jj,mikt(ji,jj),jn)   ! linear free surface 
+                        ztw(ji,jj, mikt(ji,jj) ) = pw(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) * pt_lev1(:,:,1,jn)
+                  ztw(:,:,1) = pw(:,:,1) * pt_lev1(:,:,1,jn)
                ENDIF
             ENDIF
@@ -233,5 +233,5 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     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) = (   0.5_wp * pw(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
@@ -248,9 +248,9 @@
                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) * pt_lev2(:,:,1,jn)     !!gm ISF & 4th COMPACT doesn't work
+                     ztw(ji,jj,jk) = pw(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            IF( ln_linssh )   ztw(:,:, 1 ) = pw(:,:,1) * pt_lev2(:,:,1,jn)     !!gm ISF & 4th COMPACT doesn't work
             !
          END SELECT
@@ -269,5 +269,5 @@
                   DO ji = fs_2, fs_jpim1   ! vector opt.
                      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)  )   &
+                        &           + pt_lev2(ji,jj,jk,jn) * (  pw(ji,jj,jk) - pw(ji,jj,jk+1)  )   &
                         &                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev)
                   END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbc.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbc.F90	(revision 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbc.F90	(revision 10806)
@@ -51,5 +51,5 @@
 CONTAINS
 
-   SUBROUTINE tra_bbc( kt )
+   SUBROUTINE tra_bbc( kt, ktlev, pts_rhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_bbc  ***
@@ -74,4 +74,6 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   ktlev  ! time level index for source terms
+      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends
       !
       INTEGER  ::   ji, jj    ! dummy loop indices
@@ -83,17 +85,17 @@
       IF( l_trdtra )   THEN         ! Save the input temperature trend
          ALLOCATE( ztrdt(jpi,jpj,jpk) )
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
+         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem)
       ENDIF
       !                             !  Add the geothermal trend on temperature
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
-            tsa(ji,jj,mbkt(ji,jj),jp_tem) = tsa(ji,jj,mbkt(ji,jj),jp_tem) + qgh_trd0(ji,jj) / e3t_n(ji,jj,mbkt(ji,jj))
+            pts_rhs(ji,jj,mbkt(ji,jj),jp_tem) = pts_rhs(ji,jj,mbkt(ji,jj),jp_tem) + qgh_trd0(ji,jj) / e3t(ji,jj,mbkt(ji,jj),ktlev)
          END DO
       END DO
       !
-      CALL lbc_lnk( 'trabbc', tsa(:,:,:,jp_tem) , 'T', 1. )
+      CALL lbc_lnk( 'trabbc', pts_rhs(:,:,:,jp_tem) , 'T', 1. )
       !
       IF( l_trdtra ) THEN        ! Send the trend for diagnostics
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
+         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:)
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbc, ztrdt )
          DEALLOCATE( ztrdt )
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbl.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbl.F90	(revision 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbl.F90	(revision 10806)
@@ -89,5 +89,5 @@
 
 
-   SUBROUTINE tra_bbl( kt )
+   SUBROUTINE tra_bbl( kt, ktlev1, ktlev2, kt2lev, pts_rhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE bbl  ***
@@ -101,5 +101,8 @@
       !!              is added to the general tracer trend
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step
+      INTEGER, INTENT( in ) ::   kt              ! ocean time-step
+      INTEGER, INTENT( in ) ::   ktlev1, ktlev2  ! time level indices for 3-time-level source terms
+      INTEGER, INTENT( in ) ::   kt2lev          ! time level index for 2-time-level source terms
+      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
@@ -110,13 +113,13 @@
       IF( l_trdtra )   THEN                         !* Save the T-S input trends
          ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
-      ENDIF
-
-      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl)
+         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem)
+         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal)
+      ENDIF
+
+      IF( l_bbl )   CALL bbl( kt, nit000, ktlev1, ktlev2, kt2lev, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl)
 
       IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl
          !
-         CALL tra_bbl_dif( tsb, tsa, jpts )
+         CALL tra_bbl_dif( ts(:,:,:,:,ktlev1), e3t(:,:,:,ktlev2), pts_rhs, jpts )
          IF( ln_ctl )  &
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
@@ -131,5 +134,5 @@
       IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl
          !
-         CALL tra_bbl_adv( tsb, tsa, jpts )
+         CALL tra_bbl_adv( ts(:,:,:,:,ktlev1), e3t(:,:,:,ktlev2), pts_rhs, jpts )
          IF(ln_ctl)   &
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
@@ -143,6 +146,6 @@
 
       IF( l_trdtra )   THEN                      ! send the trends for further diagnostics
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
+         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:)
+         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) - ztrds(:,:,:)
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds )
@@ -155,5 +158,5 @@
 
 
-   SUBROUTINE tra_bbl_dif( ptb, pta, kjpt )
+   SUBROUTINE tra_bbl_dif( pt, pe3t, pt_rhs, kjpt )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_bbl_dif  ***
@@ -171,5 +174,5 @@
       !!      convection is satified)
       !!
-      !! ** Action  :   pta   increased by the bbl diffusive trend
+      !! ** Action  :   pt_rhs   increased by the bbl diffusive trend
       !!
       !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
@@ -177,6 +180,7 @@
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! 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     ! tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pe3t   ! thickness fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs ! tracer trend
       !
       INTEGER  ::   ji, jj, jn   ! dummy loop indices
@@ -191,5 +195,5 @@
             DO ji = 1, jpi
                ik = mbkt(ji,jj)                             ! bottom T-level index
-               zptb(ji,jj) = ptb(ji,jj,ik,jn)               ! bottom before T and S
+               zptb(ji,jj) = pt(ji,jj,ik,jn)               ! bottom before T and S
             END DO
          END DO
@@ -198,10 +202,10 @@
             DO ji = 2, jpim1
                ik = mbkt(ji,jj)                            ! bottom T-level index
-               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                  &
+               pt_rhs(ji,jj,ik,jn) = pt_rhs(ji,jj,ik,jn)                                                  &
                   &             + (  ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )     &
                   &                - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )     &
                   &                + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )     &
                   &                - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )  )  &
-                  &             * r1_e1e2t(ji,jj) / e3t_n(ji,jj,ik)
+                  &             * r1_e1e2t(ji,jj) / pe3t(ji,jj,ik)
             END DO
          END DO
@@ -212,5 +216,5 @@
 
 
-   SUBROUTINE tra_bbl_adv( ptb, pta, kjpt )
+   SUBROUTINE tra_bbl_adv( pt, pe3t, pt_rhs, kjpt )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trc_bbl  ***
@@ -228,6 +232,7 @@
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! 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    ! before and now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pe3t   ! thickness fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs    ! tracer trend
       !
       INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
@@ -250,17 +255,17 @@
                   !
                   !                                               ! up  -slope T-point (shelf bottom point)
-                  zbtr = r1_e1e2t(iis,jj) / e3t_n(iis,jj,ikus)
-                  ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr
-                  pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra
+                  zbtr = r1_e1e2t(iis,jj) / pe3t(iis,jj,ikus)
+                  ztra = zu_bbl * ( pt(iid,jj,ikus,jn) - pt(iis,jj,ikus,jn) ) * zbtr
+                  pt_rhs(iis,jj,ikus,jn) = pt_rhs(iis,jj,ikus,jn) + ztra
                   !
                   DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
-                     zbtr = r1_e1e2t(iid,jj) / e3t_n(iid,jj,jk)
-                     ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr
-                     pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra
+                     zbtr = r1_e1e2t(iid,jj) / pe3t(iid,jj,jk)
+                     ztra = zu_bbl * ( pt(iid,jj,jk+1,jn) - pt(iid,jj,jk,jn) ) * zbtr
+                     pt_rhs(iid,jj,jk,jn) = pt_rhs(iid,jj,jk,jn) + ztra
                   END DO
                   !
-                  zbtr = r1_e1e2t(iid,jj) / e3t_n(iid,jj,ikud)
-                  ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr
-                  pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra
+                  zbtr = r1_e1e2t(iid,jj) / pe3t(iid,jj,ikud)
+                  ztra = zu_bbl * ( pt(iis,jj,ikus,jn) - pt(iid,jj,ikud,jn) ) * zbtr
+                  pt_rhs(iid,jj,ikud,jn) = pt_rhs(iid,jj,ikud,jn) + ztra
                ENDIF
                !
@@ -272,17 +277,17 @@
                   !
                   ! up  -slope T-point (shelf bottom point)
-                  zbtr = r1_e1e2t(ji,ijs) / e3t_n(ji,ijs,ikvs)
-                  ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
-                  pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra
+                  zbtr = r1_e1e2t(ji,ijs) / pe3t(ji,ijs,ikvs)
+                  ztra = zv_bbl * ( pt(ji,ijd,ikvs,jn) - pt(ji,ijs,ikvs,jn) ) * zbtr
+                  pt_rhs(ji,ijs,ikvs,jn) = pt_rhs(ji,ijs,ikvs,jn) + ztra
                   !
                   DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
-                     zbtr = r1_e1e2t(ji,ijd) / e3t_n(ji,ijd,jk)
-                     ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
-                     pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra
+                     zbtr = r1_e1e2t(ji,ijd) / pe3t(ji,ijd,jk)
+                     ztra = zv_bbl * ( pt(ji,ijd,jk+1,jn) - pt(ji,ijd,jk,jn) ) * zbtr
+                     pt_rhs(ji,ijd,jk,jn) = pt_rhs(ji,ijd,jk,jn)  + ztra
                   END DO
                   !                                               ! down-slope T-point (deep bottom point)
-                  zbtr = r1_e1e2t(ji,ijd) / e3t_n(ji,ijd,ikvd)
-                  ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
-                  pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra
+                  zbtr = r1_e1e2t(ji,ijd) / pe3t(ji,ijd,ikvd)
+                  ztra = zv_bbl * ( pt(ji,ijs,ikvs,jn) - pt(ji,ijd,ikvd,jn) ) * zbtr
+                  pt_rhs(ji,ijd,ikvd,jn) = pt_rhs(ji,ijd,ikvd,jn) + ztra
                ENDIF
             END DO
@@ -295,5 +300,5 @@
 
 
-   SUBROUTINE bbl( kt, kit000, cdtype )
+   SUBROUTINE bbl( kt, kit000, ktlev1, ktlev2, kt2lev, cdtype )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE bbl  ***
@@ -323,4 +328,6 @@
       INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
       INTEGER         , INTENT(in   ) ::   kit000   ! first time step index
+      INTEGER         , INTENT(in   ) ::   ktlev1, ktlev2  ! time level indices for 3-time-levelsource terms
+      INTEGER         , INTENT(in   ) ::   kt2lev          ! time level index for 2-time-level source terms
       CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
       !
@@ -344,10 +351,10 @@
          DO ji = 1, jpi
             ik = mbkt(ji,jj)                             ! bottom T-level index
-            zts (ji,jj,jp_tem) = tsb(ji,jj,ik,jp_tem)    ! bottom before T and S
-            zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal)
+            zts (ji,jj,jp_tem) = ts(ji,jj,ik,jp_tem,ktlev1)    ! bottom before T and S
+            zts (ji,jj,jp_sal) = ts(ji,jj,ik,jp_sal,ktlev1)
             !
-            zdep(ji,jj) = gdept_n(ji,jj,ik)              ! bottom T-level reference depth
-            zub (ji,jj) = un(ji,jj,mbku(ji,jj))          ! bottom velocity
-            zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj))
+            zdep(ji,jj) = gdept(ji,jj,ik,kt2lev)              ! bottom T-level reference depth
+            zub (ji,jj) = uu(ji,jj,mbku(ji,jj),ktlev2)          ! bottom velocity
+            zvb (ji,jj) = vv(ji,jj,mbkv(ji,jj),ktlev2)
          END DO
       END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tradmp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tradmp.F90	(revision 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tradmp.F90	(revision 10806)
@@ -72,5 +72,5 @@
 
 
-   SUBROUTINE tra_dmp( kt )
+   SUBROUTINE tra_dmp( kt, ktlev, kt2lev, pts_rhs )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tra_dmp  ***
@@ -90,5 +90,8 @@
       !! ** Action  : - tsa: tracer trends updated with the damping trend
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
+      INTEGER, INTENT(in) ::   ktlev  ! time level index for 3-time-level source terms
+      INTEGER, INTENT(in) ::   kt2lev ! time level index for 2-time-level source terms
+      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends
       !
       INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
@@ -101,5 +104,5 @@
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
          ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) ) 
-         ztrdts(:,:,:,:) = tsa(:,:,:,:) 
+         ztrdts(:,:,:,:) = pts_rhs(:,:,:,:) 
       ENDIF
       !                           !==  input T-S data at kt  ==!
@@ -113,5 +116,5 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - tsb(ji,jj,jk,jn) )
+                     pts_rhs(ji,jj,jk,jn) = pts_rhs(ji,jj,jk,jn) + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - ts(ji,jj,jk,jn,ktlev) )
                   END DO
                END DO
@@ -124,8 +127,8 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   IF( avt(ji,jj,jk) <= avt_c ) THEN
-                     tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
-                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )
-                     tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
-                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
+                     pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem)   &
+                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - ts(ji,jj,jk,jp_tem,ktlev) )
+                     pts_rhs(ji,jj,jk,jp_sal) = pts_rhs(ji,jj,jk,jp_sal)   &
+                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - ts(ji,jj,jk,jp_sal,ktlev) )
                   ENDIF
                END DO
@@ -137,9 +140,9 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  IF( gdept_n(ji,jj,jk) >= hmlp (ji,jj) ) THEN
-                     tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
-                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )
-                     tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
-                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
+                  IF( gdept(ji,jj,jk,kt2lev) >= hmlp (ji,jj) ) THEN
+                     pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem)   &
+                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - ts(ji,jj,jk,jp_tem,ktlev) )
+                     pts_rhs(ji,jj,jk,jp_sal) = pts_rhs(ji,jj,jk,jp_sal)   &
+                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - ts(ji,jj,jk,jp_sal,ktlev) )
                   ENDIF
                END DO
@@ -150,5 +153,5 @@
       !
       IF( l_trdtra )   THEN       ! trend diagnostic
-         ztrdts(:,:,:,:) = tsa(:,:,:,:) - ztrdts(:,:,:,:)
+         ztrdts(:,:,:,:) = pts_rhs(:,:,:,:) - ztrdts(:,:,:,:)
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) )
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf.F90	(revision 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf.F90	(revision 10806)
@@ -47,5 +47,5 @@
 CONTAINS
 
-   SUBROUTINE tra_ldf( kt )
+   SUBROUTINE tra_ldf( kt, ktlev1, ktlev2, kt2lev, pts_rhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf  ***
@@ -53,5 +53,8 @@
       !! ** Purpose :   compute the lateral ocean tracer physics.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      INTEGER, INTENT( in ) ::   kt              ! ocean time-step index
+      INTEGER, INTENT( in ) ::   ktlev1, ktlev2  ! time level indices for 3-time-level source terms
+      INTEGER, INTENT( in ) ::   kt2lev          ! time level index for 2-time-level source terms
+      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends
       !!
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
@@ -62,22 +65,22 @@
       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
       !
       SELECT CASE ( nldf_tra )                 !* compute lateral mixing trend and add it to the general trend
       CASE ( np_lap   )                                  ! laplacian: iso-level operator
-         CALL tra_ldf_lap  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb,      tsa, jpts,  1   )
+         CALL tra_ldf_lap  ( kt, nit000, ktlev2, kt2lev, 'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,ktlev1),      pts_rhs, jpts,  1   )
       CASE ( np_lap_i )                                  ! laplacian: standard iso-neutral operator (Madec)
-         CALL tra_ldf_iso  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   )
+         CALL tra_ldf_iso  ( kt, nit000, ktlev2, kt2lev, 'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev1), pts_rhs, jpts,  1   )
       CASE ( np_lap_it )                                 ! laplacian: triad iso-neutral operator (griffies)
-         CALL tra_ldf_triad( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   )
+         CALL tra_ldf_triad( kt, nit000, ktlev2, kt2lev, 'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev1), pts_rhs, jpts,  1   )
       CASE ( np_blp , np_blp_i , np_blp_it )             ! bilaplacian: iso-level & iso-neutral operators
-         CALL tra_ldf_blp  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb      , tsa, jpts, nldf_tra )
+         CALL tra_ldf_blp  ( kt, nit000, ktlev2, kt2lev, 'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,ktlev1)      , pts_rhs, jpts, nldf_tra )
       END SELECT
       !
       IF( l_trdtra )   THEN                    !* save the horizontal diffusive trends for further diagnostics
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
+         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:)
+         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) - ztrds(:,:,:)
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrdt )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrds )
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_iso.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_iso.F90	(revision 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_iso.F90	(revision 10806)
@@ -48,7 +48,7 @@
 CONTAINS
 
-  SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
+  SUBROUTINE tra_ldf_iso( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu , pgv ,   &
       &                                                   pgui, pgvi,   &
-      &                                       ptb , ptbb, pta , kjpt, kpass )
+      &                                       pt , pt_lev0, pt_rhs , kjpt, kpass )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf_iso  ***
@@ -87,10 +87,12 @@
       !!         difft = 1/(e1e2t*e3t) dk[ zftw ]
       !!      Add this trend to the general trend (ta,sa):
-      !!         pta = pta + difft
-      !!
-      !! ** Action :   Update pta arrays with the before rotated diffusion
+      !!         pt_rhs = pt_rhs + difft
+      !!
+      !! ** Action :   Update pt_rhs arrays with the before rotated diffusion
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
       INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index
+      INTEGER                              , INTENT(in   ) ::   ktlev      ! time level index for e3t
+      INTEGER                              , INTENT(in   ) ::   kt2lev     ! time level index for 2-time-level thicknesses
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
@@ -99,7 +101,7 @@
       REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! tracer (kpass=1) or laplacian of tracer (kpass=2)
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptbb       ! tracer (only used in kpass=2)
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt        ! tracer (kpass=1) or laplacian of tracer (kpass=2)
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev0       ! tracer (only used in kpass=2)
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs        ! tracer trend
       !
       INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
@@ -182,5 +184,5 @@
                      DO ji = 1, fs_jpim1
                         akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   &
-                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  )
+                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,kt2lev) * e3w(ji,jj,jk,kt2lev) )  )
                      END DO
                   END DO
@@ -190,5 +192,5 @@
                   DO jj = 1, jpjm1
                      DO ji = 1, fs_jpim1
-                        ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk)
+                        ze3w_2 = e3w(ji,jj,jk,kt2lev) * e3w(ji,jj,jk,kt2lev)
                         zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
                         akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
@@ -219,6 +221,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
-                  zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
+                  zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
+                  zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
                END DO
             END DO
@@ -248,13 +250,13 @@
             !
             !                             !== Vertical tracer gradient
-            zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * wmask(:,:,jk+1)     ! level jk+1
+            zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * wmask(:,:,jk+1)     ! level jk+1
             !
             IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:)                          ! surface: zdkt(jk=1)=zdkt(jk=2)
-            ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk)
+            ELSE                 ;   zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * wmask(:,:,jk)
             ENDIF
             DO jj = 1 , jpjm1            !==  Horizontal fluxes
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)
-                  zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)
+                  zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,ktlev)
+                  zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,ktlev)
                   !
                   zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   &
@@ -276,9 +278,9 @@
             END DO
             !
-            DO jj = 2 , jpjm1          !== horizontal divergence and add to pta
+            DO jj = 2 , jpjm1          !== horizontal divergence and add to pt_rhs
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      &
+                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      &
                      &                                           + zftv(ji,jj,jk) - zftv(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
@@ -325,7 +327,7 @@
                DO jj = 1, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)   &
+                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,kt2lev) * wmask(ji,jj,jk)   &
                         &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
-                        &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
+                        &                            * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
                   END DO
                END DO
@@ -340,15 +342,15 @@
                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    &
                            &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   &
-                           &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)
+                           &           * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) / e3w(ji,jj,jk,kt2lev) * wmask(ji,jj,jk)
                      END DO
                   END DO
                END DO 
-            CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp.
+            CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt_lev0 gradients, resp.
                DO jk = 2, jpkm1 
                   DO jj = 1, jpjm1
                      DO ji = fs_2, fs_jpim1
-                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)                      &
-                           &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   &
-                           &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   )
+                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,kt2lev) * wmask(ji,jj,jk)                      &
+                           &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   &
+                           &                               + akz     (ji,jj,jk) * ( pt_lev0(ji,jj,jk-1,jn) - pt_lev0(ji,jj,jk,jn) )   )
                      END DO
                   END DO
@@ -357,9 +359,9 @@
          ENDIF
          !         
-         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==!
+         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pt_rhs  ==!
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(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) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   &
+                     &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev)
                END DO
             END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_lap_blp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_lap_blp.F90	(revision 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_lap_blp.F90	(revision 10806)
@@ -45,7 +45,7 @@
 CONTAINS
 
-   SUBROUTINE tra_ldf_lap( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
+   SUBROUTINE tra_ldf_lap( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu , pgv ,   &
       &                                                    pgui, pgvi,   &
-      &                                        ptb , pta , kjpt, kpass ) 
+      &                                        pt , pt_rhs , kjpt, kpass ) 
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf_lap  ***
@@ -59,12 +59,14 @@
       !!          difft = 1/(e1e2t*e3t) {  di-1[ pahu e2u*e3u/e1u di(tb) ]
       !!                                 + dj-1[ pahv e1v*e3v/e2v dj(tb) ] }
-      !!      Add this trend to the general tracer trend pta :
-      !!          pta = pta + difft
-      !!
-      !! ** Action  : - Update pta arrays with the before iso-level 
+      !!      Add this trend to the general tracer trend pt_rhs :
+      !!          pt_rhs = pt_rhs + difft
+      !!
+      !! ** Action  : - Update pt_rhs arrays with the before iso-level 
       !!                harmonic mixing trend.
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
       INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index
+      INTEGER                              , INTENT(in   ) ::   ktlev      ! time level index for e3t
+      INTEGER                              , INTENT(in   ) ::   kt2lev     ! time level index for 2-time-level thicknesses
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
@@ -73,6 +75,6 @@
       REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! 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        ! 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
@@ -100,6 +102,6 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)   !!gm   * umask(ji,jj,jk) pah masked!
-               zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)   !!gm   * vmask(ji,jj,jk)
+               zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,ktlev)   !!gm   * umask(ji,jj,jk) pah masked!
+               zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,ktlev)   !!gm   * vmask(ji,jj,jk)
             END DO
          END DO
@@ -113,6 +115,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
-                  ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
+                  ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) )
+                  ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) )
                END DO
             END DO
@@ -138,7 +140,7 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     &
+                  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) )   &
-                     &                                / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) )
+                     &                                / ( e1e2t(ji,jj) * e3t(ji,jj,jk,ktlev) )
                END DO
             END DO
@@ -159,7 +161,7 @@
    
 
-   SUBROUTINE tra_ldf_blp( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
+   SUBROUTINE tra_ldf_blp( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu , pgv ,   &
       &                                                    pgui, pgvi,   &
-      &                                                    ptb , pta , kjpt, kldf )
+      &                                                    pt , pt_rhs , kjpt, kldf )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE tra_ldf_blp  ***
@@ -172,8 +174,10 @@
       !!      It is computed by two successive calls to laplacian routine
       !!
-      !! ** Action :   pta   updated with the before rotated bilaplacian diffusion
+      !! ** Action :   pt_rhs   updated with the before rotated bilaplacian diffusion
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
       INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index
+      INTEGER                              , INTENT(in   ) ::   ktlev      ! time level index for e3t
+      INTEGER                              , INTENT(in   ) ::   kt2lev     ! time level index for 2-time-level thicknesses
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
@@ -182,6 +186,6 @@
       REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top levels
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! 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        ! 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
@@ -203,12 +207,12 @@
       zlap(:,:,:,:) = 0._wp
       !
-      SELECT CASE ( kldf )       !==  1st laplacian applied to ptb (output in zlap)  ==!
+      SELECT CASE ( kldf )       !==  1st laplacian applied to pt (output in zlap)  ==!
       !
       CASE ( np_blp    )               ! iso-level bilaplacian
-         CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb,      zlap, kjpt, 1 )
+         CALL tra_ldf_lap  ( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt,      zlap, kjpt, 1 )
       CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec)
-         CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 )
+         CALL tra_ldf_iso  ( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt, pt, zlap, kjpt, 1 )
       CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies)
-         CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 )
+         CALL tra_ldf_triad( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt, pt, zlap, kjpt, 1 )
       END SELECT
       !
@@ -219,12 +223,12 @@
       ENDIF
       !
-      SELECT CASE ( kldf )       !==  2nd laplacian applied to zlap (output in pta)  ==!
+      SELECT CASE ( kldf )       !==  2nd laplacian applied to zlap (output in pt_rhs)  ==!
       !
       CASE ( np_blp    )               ! iso-level bilaplacian
-         CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pta,      kjpt, 2 )
+         CALL tra_ldf_lap  ( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pt_rhs,      kjpt, 2 )
       CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec)
-         CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 )
+         CALL tra_ldf_iso  ( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pt, pt_rhs, kjpt, 2 )
       CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies)
-         CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 )
+         CALL tra_ldf_triad( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pt, pt_rhs, kjpt, 2 )
       END SELECT
       !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_triad.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_triad.F90	(revision 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_triad.F90	(revision 10806)
@@ -48,7 +48,7 @@
 CONTAINS
 
-  SUBROUTINE tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
+  SUBROUTINE tra_ldf_triad( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu , pgv ,   &
       &                                                     pgui, pgvi,   &
-      &                                         ptb , ptbb, pta , kjpt, kpass )
+      &                                         pt , pt_lev0, pt_rhs , kjpt, kpass )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf_triad  ***
@@ -66,5 +66,5 @@
       !!      see documentation for the desciption
       !!
-      !! ** Action :   pta   updated with the before rotated diffusion
+      !! ** Action :   pt_rhs   updated with the before rotated diffusion
       !!               ah_wslp2 ....
       !!               akz   stabilizing vertical diffusivity coefficient (used in trazdf_imp)
@@ -72,4 +72,6 @@
       INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
       INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index
+      INTEGER                              , INTENT(in   ) ::   ktlev      ! time level index for e3t
+      INTEGER                              , INTENT(in   ) ::   kt2lev     ! time level index for 2-time-level thicknesses
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
@@ -78,7 +80,7 @@
       REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu , pgv  ! tracer gradient at pstep levels
       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! tracer (kpass=1) or laplacian of tracer (kpass=2)
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptbb       ! tracer (only used in kpass=2)
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt        ! tracer (kpass=1) or laplacian of tracer (kpass=2)
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev0       ! tracer (only used in kpass=2)
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs        ! tracer trend
       !
       INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
@@ -142,10 +144,10 @@
                   DO jj = 1, jpjm1
                      DO ji = 1, fs_jpim1
-                        ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp)
-                        zbu   = e1e2u(ji,jj) * e3u_n(ji,jj,jk)
+                        ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,kt2lev)
+                        zbu   = e1e2u(ji,jj) * e3u(ji,jj,jk,ktlev)
                         zah   = 0.25_wp * pahu(ji,jj,jk)
                         zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
                         ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth)
-                        zslope2 = zslope_skew + ( gdept_n(ji+1,jj,jk) - gdept_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
+                        zslope2 = zslope_skew + ( gdept(ji+1,jj,jk,kt2lev) - gdept(ji,jj,jk,kt2lev) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
                         zslope2 = zslope2 *zslope2
                         ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji+ip,jj) * zslope2
@@ -166,11 +168,11 @@
                   DO jj = 1, jpjm1
                      DO ji = 1, fs_jpim1
-                        ze3wr = 1.0_wp / e3w_n(ji,jj+jp,jk+kp)
-                        zbv   = e1e2v(ji,jj) * e3v_n(ji,jj,jk)
+                        ze3wr = 1.0_wp / e3w(ji,jj+jp,jk+kp,kt2lev)
+                        zbv   = e1e2v(ji,jj) * e3v(ji,jj,jk,ktlev)
                         zah   = 0.25_wp * pahv(ji,jj,jk)
                         zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
                         ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces
                         !    (do this by *adding* gradient of depth)
-                        zslope2 = zslope_skew + ( gdept_n(ji,jj+1,jk) - gdept_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
+                        zslope2 = zslope_skew + ( gdept(ji,jj+1,jk,kt2lev) - gdept(ji,jj,jk,kt2lev) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
                         zslope2 = zslope2 * zslope2
                         ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj+jp) * zslope2
@@ -193,5 +195,5 @@
                      DO ji = 1, fs_jpim1
                         akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   &
-                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  )
+                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,kt2lev) * e3w(ji,jj,jk,kt2lev) )  )
                      END DO
                   END DO
@@ -201,5 +203,5 @@
                   DO jj = 1, jpjm1
                      DO ji = 1, fs_jpim1
-                        ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk)
+                        ze3w_2 = e3w(ji,jj,jk,kt2lev) * e3w(ji,jj,jk,kt2lev)
                         zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
                         akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
@@ -229,6 +231,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
-                  zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
+                  zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
+                  zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
                END DO
             END DO
@@ -257,9 +259,9 @@
          DO jk = 1, jpkm1
             !                    !==  Vertical tracer gradient at level jk and jk+1
-            zdkt3d(:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
+            zdkt3d(:,:,1) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
             !
             !                    ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1)
             IF( jk == 1 ) THEN   ;   zdkt3d(:,:,0) = zdkt3d(:,:,1)
-            ELSE                 ;   zdkt3d(:,:,0) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk)
+            ELSE                 ;   zdkt3d(:,:,0) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)
             ENDIF
             !
@@ -273,10 +275,10 @@
                            ze1ur = r1_e1u(ji,jj)
                            zdxt  = zdit(ji,jj,jk) * ze1ur
-                           ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp)
+                           ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,kt2lev)
                            zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr
                            zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
                            zslope_iso  = triadi  (ji+ip,jj,jk,1-ip,kp)
                            !
-                           zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk)
+                           zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,ktlev)
                            ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????   ahu is masked....
                            zah = pahu(ji,jj,jk)
@@ -296,9 +298,9 @@
                            ze2vr = r1_e2v(ji,jj)
                            zdyt  = zdjt(ji,jj,jk) * ze2vr
-                           ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp)
+                           ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,kt2lev)
                            zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr
                            zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
                            zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp)
-                           zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk)
+                           zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,ktlev)
                            ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????  ahv is masked...
                            zah = pahv(ji,jj,jk)
@@ -320,10 +322,10 @@
                            ze1ur = r1_e1u(ji,jj)
                            zdxt  = zdit(ji,jj,jk) * ze1ur
-                           ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp)
+                           ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,kt2lev)
                            zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr
                            zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
                            zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp)
                            !
-                           zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk)
+                           zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,ktlev)
                            ! ln_botmix_triad is .F. mask zah for bottom half cells
                            zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp)         ! pahu(ji+ip,jj,jk)   ===>>  ????
@@ -343,9 +345,9 @@
                            ze2vr = r1_e2v(ji,jj)
                            zdyt  = zdjt(ji,jj,jk) * ze2vr
-                           ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp)
+                           ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,kt2lev)
                            zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr
                            zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
                            zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp)
-                           zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk)
+                           zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,ktlev)
                            ! ln_botmix_triad is .F. mask zah for bottom half cells
                            zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp)         ! pahv(ji,jj+jp,jk)  ????
@@ -362,7 +364,7 @@
             DO jj = 2 , jpjm1
                DO ji = fs_2, fs_jpim1  ! vector opt.
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       &
+                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       &
                      &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   &
-                     &                                        / (  e1e2t(ji,jj) * e3t_n(ji,jj,jk)  )
+                     &                                        / (  e1e2t(ji,jj) * e3t(ji,jj,jk,ktlev)  )
                END DO
             END DO
@@ -375,7 +377,7 @@
                DO jj = 1, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)   &
+                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,kt2lev) * tmask(ji,jj,jk)   &
                         &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
-                        &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
+                        &                            * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
                   END DO
                END DO
@@ -387,16 +389,16 @@
                   DO jj = 1, jpjm1
                      DO ji = fs_2, fs_jpim1
-                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)             &
-                           &                            * ah_wslp2(ji,jj,jk) * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
+                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,kt2lev) * tmask(ji,jj,jk)             &
+                           &                            * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
                      END DO
                   END DO
                END DO 
-            CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp.
+            CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt_lev0 gradients, resp.
                DO jk = 2, jpkm1 
                   DO jj = 1, jpjm1
                      DO ji = fs_2, fs_jpim1
-                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)                      &
-                           &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   &
-                           &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   )
+                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,kt2lev) * tmask(ji,jj,jk)                      &
+                           &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   &
+                           &                               + akz     (ji,jj,jk) * ( pt_lev0(ji,jj,jk-1,jn) - pt_lev0(ji,jj,jk,jn) )   )
                      END DO
                   END DO
@@ -405,9 +407,9 @@
          ENDIF
          !
-         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==!
+         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pt_rhs  ==!
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1  ! vector opt.
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
-                     &                                        / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) )
+                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
+                     &                                        / ( e1e2t(ji,jj) * e3t(ji,jj,jk,ktlev) )
                END DO
             END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traqsr.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traqsr.F90	(revision 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traqsr.F90	(revision 10806)
@@ -75,5 +75,5 @@
 CONTAINS
 
-   SUBROUTINE tra_qsr( kt )
+   SUBROUTINE tra_qsr( kt, ktlev, kt2lev, pts_rhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_qsr  ***
@@ -102,4 +102,7 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! ocean time-step
+      INTEGER, INTENT(in) ::   ktlev  ! time level index for 3-time-level source terms
+      INTEGER, INTENT(in) ::   kt2lev ! time level index for 2-time-level source terms
+      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends
       !
       INTEGER  ::   ji, jj, jk               ! dummy loop indices
@@ -126,5 +129,5 @@
       IF( l_trdtra ) THEN      ! trends diagnostic: save the input temperature trend
          ALLOCATE( ztrdt(jpi,jpj,jpk) ) 
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
+         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem)
       ENDIF
       !
@@ -172,5 +175,5 @@
                      zze     = 568.2 * zCtot**(-0.746)
                      IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293)
-                     zpsi    = gdepw_n(ji,jj,jk) / zze
+                     zpsi    = gdepw(ji,jj,jk,kt2lev) / zze
                      !
                      zlogc   = LOG( zchl )
@@ -218,8 +221,8 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-                  zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * xsi0r       )
-                  zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekb(ji,jj) )
-                  zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekg(ji,jj) )
-                  zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekr(ji,jj) )
+                  zc0 = ze0(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,ktlev) * xsi0r       )
+                  zc1 = ze1(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,ktlev) * zekb(ji,jj) )
+                  zc2 = ze2(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,ktlev) * zekg(ji,jj) )
+                  zc3 = ze3(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,ktlev) * zekr(ji,jj) )
                   ze0(ji,jj,jk) = zc0
                   ze1(ji,jj,jk) = zc1
@@ -248,6 +251,6 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-                  zc0 = zz0 * EXP( -gdepw_n(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk  )*xsi1r )
-                  zc1 = zz0 * EXP( -gdepw_n(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk+1)*xsi1r )
+                  zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,kt2lev)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,kt2lev)*xsi1r )
+                  zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,kt2lev)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,kt2lev)*xsi1r )
                   qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 
                END DO
@@ -261,6 +264,6 @@
          DO jj = 2, jpjm1        !-----------------------------!
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
-                  &                 + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t_n(ji,jj,jk)
+               pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem)   &
+                  &                 + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t(ji,jj,jk,ktlev)
             END DO
          END DO
@@ -295,5 +298,5 @@
       !
       IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
+         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:)
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt )
          DEALLOCATE( ztrdt ) 
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trasbc.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trasbc.F90	(revision 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trasbc.F90	(revision 10806)
@@ -51,5 +51,5 @@
 CONTAINS
 
-   SUBROUTINE tra_sbc ( kt )
+   SUBROUTINE tra_sbc ( kt, ktlev, pts_rhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_sbc  ***
@@ -63,5 +63,5 @@
       !!      (2) Fwe , tracer carried with the water that is exchanged with air+ice. 
       !!               The input forcing fields (emp, rnf, sfx, isf) contain Fext+Fwe,
-      !!             they are simply added to the tracer trend (tsa).
+      !!             they are simply added to the tracer trend (pts_rhs).
       !!               In linear free surface case (ln_linssh=T), the volume of the
       !!             ocean does not change with the water exchanges at the (air+ice)-sea
@@ -69,8 +69,10 @@
       !!             concentration/dilution effect associated with water exchanges.
       !!
-      !! ** Action  : - Update tsa with the surface boundary condition trend 
+      !! ** Action  : - Update pts_rhs with the surface boundary condition trend 
       !!              - send trends to trdtra module for further diagnostics(l_trdtra=T)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
+      INTEGER, INTENT(in) ::   ktlev  ! time level index for source terms
+      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends
       !
       INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices  
@@ -90,6 +92,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
       !
@@ -131,10 +133,10 @@
          DO jj = 2, jpj                         !==>> add concentration/dilution effect due to constant volume cell
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_tem)
-               sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_sal)
+               sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * ts(ji,jj,1,jp_tem,ktlev)
+               sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * ts(ji,jj,1,jp_sal,ktlev)
             END DO
          END DO                                 !==>> output c./d. term
-         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) )
-         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) )
+         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * ts(:,:,1,jp_tem,ktlev) )
+         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * ts(:,:,1,jp_sal,ktlev) )
       ENDIF
       !
@@ -142,5 +144,5 @@
          DO jj = 2, jpj
             DO ji = fs_2, fs_jpim1   ! vector opt.  
-               tsa(ji,jj,1,jn) = tsa(ji,jj,1,jn) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t_n(ji,jj,1)
+               pts_rhs(ji,jj,1,jn) = pts_rhs(ji,jj,1,jn) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t(ji,jj,1,ktlev)
             END DO
          END DO
@@ -173,5 +175,5 @@
                DO jk = ikt, ikb - 1
                ! compute trend
-                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                                &
+                  pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem)                                                &
                      &           + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             &
                      &           * r1_hisf_tbl(ji,jj)
@@ -180,5 +182,5 @@
                ! level partially include in ice shelf boundary layer 
                ! compute trend
-               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                                 &
+               pts_rhs(ji,jj,ikb,jp_tem) = pts_rhs(ji,jj,ikb,jp_tem)                                                 &
                   &              + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             &
                   &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)
@@ -199,7 +201,7 @@
                   zdep = zfact / h_rnf(ji,jj)
                   DO jk = 1, nk_rnf(ji,jj)
-                                        tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                 &
+                                        pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem)                                 &
                                            &                 +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep
-                     IF( ln_rnf_sal )   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                 &
+                     IF( ln_rnf_sal )   pts_rhs(ji,jj,jk,jp_sal) = pts_rhs(ji,jj,jk,jp_sal)                                 &
                                            &                 +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 
                   END DO
@@ -209,6 +211,6 @@
       ENDIF
 
-      IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*tsn(:,:,1,jp_tem) )   ! runoff term on sst
-      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*tsn(:,:,1,jp_sal) )   ! runoff term on sss
+      IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*ts(:,:,1,jp_tem,ktlev) )   ! runoff term on sst
+      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*ts(:,:,1,jp_sal,ktlev) )   ! runoff term on sss
 
 #if defined key_asminc
@@ -223,7 +225,7 @@
             DO jj = 2, jpj 
                DO ji = fs_2, fs_jpim1
-                  ztim = ssh_iau(ji,jj) / e3t_n(ji,jj,1)
-                  tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + tsn(ji,jj,1,jp_tem) * ztim
-                  tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + tsn(ji,jj,1,jp_sal) * ztim
+                  ztim = ssh_iau(ji,jj) / e3t(ji,jj,1,ktlev)
+                  pts_rhs(ji,jj,1,jp_tem) = pts_rhs(ji,jj,1,jp_tem) + ts(ji,jj,1,jp_tem,ktlev) * ztim
+                  pts_rhs(ji,jj,1,jp_sal) = pts_rhs(ji,jj,1,jp_sal) + ts(ji,jj,1,jp_sal,ktlev) * ztim
                END DO
             END DO
@@ -232,6 +234,6 @@
                DO ji = fs_2, fs_jpim1
                   ztim = ssh_iau(ji,jj) / ( ht_n(ji,jj) + 1. - ssmask(ji, jj) )
-                  tsa(ji,jj,:,jp_tem) = tsa(ji,jj,:,jp_tem) + tsn(ji,jj,:,jp_tem) * ztim
-                  tsa(ji,jj,:,jp_sal) = tsa(ji,jj,:,jp_sal) + tsn(ji,jj,:,jp_sal) * ztim
+                  pts_rhs(ji,jj,:,jp_tem) = pts_rhs(ji,jj,:,jp_tem) + ts(ji,jj,:,jp_tem,ktlev) * ztim
+                  pts_rhs(ji,jj,:,jp_sal) = pts_rhs(ji,jj,:,jp_sal) + ts(ji,jj,:,jp_sal,ktlev) * ztim
                END DO  
             END DO  
@@ -250,7 +252,7 @@
             DO jj = 2, jpj 
                DO ji = fs_2, fs_jpim1
-                  zdep = 1._wp / e3t_n(ji,jj,jk) 
-                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem) * zdep
-                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal) * zdep  
+                  zdep = 1._wp / e3t(ji,jj,jk,ktlev) 
+                  pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem) * zdep
+                  pts_rhs(ji,jj,jk,jp_sal) = pts_rhs(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal) * zdep  
                END DO  
             END DO  
@@ -259,6 +261,6 @@
 
       IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
+         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:)
+         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) - ztrds(:,:,:)
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds )
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 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/oce.F90	(revision 10806)
@@ -22,8 +22,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:),   TARGET ::   ww             !: vertical velocity            [m/s]
    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(:,:,:), TARGET   ::   hdiv           !: horizontal divergence        [s-1]
    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]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:), TARGET :: r_ab           !: thermal/haline expansion coef. [Celsius-1,psu-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:),TARGET    :: r_n2           !: brunt-vaisala frequency**2     [s-2]
    !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rhd    !: in situ density anomalie rhd=(rho-rau0)/rau0  [no units]
@@ -72,4 +72,7 @@
    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(:,:,:)   ::           hdivn          !: horizontal divergence        [s-1]
+   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:,:) ::   rab_b,  rab_n          !: thermal/haline expansion coef. [Celsius-1,psu-1]
+   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::   rn2b ,  rn2            !: brunt-vaisala frequency**2     [s-2]
    REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:,:) ::   tsb  ,  tsn   , tsa    !: 4D T-S fields                [Celsius,psu]         
    !! TEMPORARY POINTERS FOR DEVELOPMENT ONLY
@@ -91,9 +94,8 @@
       ierr(:) = 0 
       ALLOCATE( uu   (jpi,jpj,jpk, jpt) , vv   (jpi,jpj,jpk,jpt)  ,                             &
-         &      ww   (jpi,jpj,jpk)      , hdivn(jpi,jpj,jpk)      ,                             &
+         &      ww   (jpi,jpj,jpk)      , hdiv(jpi,jpj,jpk)       ,                             &
          &      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)      ,                             &
-         &      rhd  (jpi,jpj,jpk)      , rhop (jpi,jpj,jpk)                              , STAT=ierr(1) )
+         &      r_ab (jpi,jpj,jpk,jpts,jpt), r_n2 (jpi,jpj,jpk,jpt),                            &
+         &      rhd  (jpi,jpj,jpk)         , rhop (jpi,jpj,jpk)                              , STAT=ierr(1) )
          !
       ALLOCATE( sshb(jpi,jpj)     , sshn(jpi,jpj)   , ssha(jpi,jpj)   ,     &
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 10802)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10806)
@@ -178,6 +178,6 @@
 !!jc: fs simplification
                             
-                         ua(:,:,:) = 0._wp            ! set dynamics trends to zero
-                         va(:,:,:) = 0._wp
+                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
+                         vv(:,:,:,Nrhs) = 0._wp
 
       IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   &
@@ -190,5 +190,5 @@
                          CALL dyn_adv       ( kstp, Nm1, Nnn, uu(:,:,:,Nrhs), vv(:,:,:,Nrhs) )  ! advection (vector or flux form)
                          CALL dyn_vor       ( kstp,      Nnn, uu(:,:,:,Nrhs), vv(:,:,:,Nrhs) )  ! vorticity term including Coriolis
-                         CALL dyn_ldf       ( kstp )  ! lateral mixing
+                         CALL dyn_ldf       ( kstp, Nm1, Nnn, uu(:,:,:,Nrhs), vv(:,:,:,Nrhs) )  ! lateral mixing
       IF( ln_zdfosm  )   CALL dyn_osm       ( kstp )  ! OSMOSIS non-local velocity fluxes
                          CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure
@@ -233,13 +233,13 @@
       ! Active tracers                              
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         tsa(:,:,:,:) = 0._wp         ! set tracer trends to zero
+                         ts(:,:,:,:,Nrhs) = 0._wp         ! set tracer trends to zero
 
       IF(  lk_asminc .AND. ln_asmiau .AND. &
          & ln_trainc )   CALL tra_asm_inc   ( kstp )  ! apply tracer assimilation increment
-                         CALL tra_sbc       ( kstp )  ! surface boundary condition
-      IF( ln_traqsr  )   CALL tra_qsr       ( kstp )  ! penetrative solar radiation qsr
-      IF( ln_trabbc  )   CALL tra_bbc       ( kstp )  ! bottom heat flux
-      IF( ln_trabbl  )   CALL tra_bbl       ( kstp )  ! advective (and/or diffusive) bottom boundary layer scheme
-      IF( ln_tradmp  )   CALL tra_dmp       ( kstp )  ! internal damping trends
+                         CALL tra_sbc       ( kstp, Nnn, ts(:,:,:,:,Nrhs) )  ! surface boundary condition
+      IF( ln_traqsr  )   CALL tra_qsr       ( kstp, Nnn, Nnn_2lev, ts(:,:,:,:,Nrhs) )  ! penetrative solar radiation qsr
+      IF( ln_trabbc  )   CALL tra_bbc       ( kstp, Nnn, ts(:,:,:,:,Nrhs) )  ! bottom heat flux
+      IF( ln_trabbl  )   CALL tra_bbl       ( kstp, Nm1, Nnn, Nnn_2lev, ts(:,:,:,:,Nrhs) )  ! advective (and/or diffusive) bottom boundary layer scheme
+      IF( ln_tradmp  )   CALL tra_dmp       ( kstp, Nm1, Nnn_2lev, ts(:,:,:,:,Nrhs) )  ! internal damping trends
       IF( ln_bdy     )   CALL bdy_tra_dmp   ( kstp )  ! bdy damping trends
 #if defined key_agrif
@@ -247,9 +247,9 @@
                &         CALL Agrif_Sponge_tra        ! tracers sponge
 #endif
-                         CALL tra_adv       ( kstp, Nm1, Nnn, Np1, ts(:,:,:,:,Nrhs) )  ! horizontal & vertical advection
+                         CALL tra_adv       ( kstp, Nm1, Nnn, Np1, Nnn_2lev, ts(:,:,:,:,Nrhs) )  ! horizontal & vertical advection
       IF( ln_zdfosm  )   CALL tra_osm       ( kstp )  ! OSMOSIS non-local tracer fluxes
       IF( lrst_oce .AND. ln_zdfosm ) &
            &             CALL osm_rst( kstp, 'WRITE' )! write OSMOSIS outputs + wn (so must do here) to restarts
-                         CALL tra_ldf       ( kstp )  ! lateral mixing
+                         CALL tra_ldf       ( kstp, Nm1, Nnn, Nnn_2lev, ts(:,:,:,:,Nrhs)  )  ! lateral mixing
 
 !!gm : why CALL to dia_ptr has been moved here??? (use trends info?)
@@ -343,4 +343,7 @@
       vb => vv(:,:,:,Nm1); vn => vv(:,:,:,Nnn); va => vv(:,:,:,Np1)
       wn => ww(:,:,:)
+      hdivn => hdiv(:,:,:)
+      rab_b  => r_ab  (:,:,:,:,Nm1_2lev); rab_n  => r_ab (:,:,:,:,Nnn_2lev)
+      rn2b   => r_n2 (:,:,:,Nm1_2lev)   ; rn2    => r_n2 (:,:,:,Nnn_2lev)
 
       tsb => ts(:,:,:,:,Nm1); tsn => ts(:,:,:,:,Nnn); tsa => ts(:,:,:,:,Np1)
@@ -360,10 +363,4 @@
       gde3w_n => gde3w(:,:,:)
 
-      ht_n => ht(:,:)
-      hu_b => hu(:,:,Nm1); hu_n => hu(:,:,Nnn); hu_a => hu(:,:,Np1)
-      hv_b => hv(:,:,Nm1); hv_n => hv(:,:,Nnn); hv_a => hv(:,:,Np1)
-      r1_hu_b => r1_hu(:,:,Nm1); r1_hu_n => r1_hu(:,:,Nnn); r1_hu_a => r1_hu(:,:,Np1)
-      r1_hv_b => r1_hv(:,:,Nm1); r1_hv_n => r1_hv(:,:,Nnn); r1_hv_a => r1_hv(:,:,Np1)
-
    END SUBROUTINE update_pointers
 
