Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRD/trdtra.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRD/trdtra.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRD/trdtra.F90	(revision 10966)
@@ -172,5 +172,5 @@
          END SELECT
          !                            ! send trend to trd_trc
-         CALL trd_trc( ztrds, ktra, ktrd, kt ) 
+         CALL trd_trc( ztrds, ktra, ktrd, kt, Kmm ) 
          !
       ENDIF
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRD/trdtrc.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRD/trdtrc.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRD/trdtrc.F90	(revision 10966)
@@ -9,6 +9,7 @@
 CONTAINS
 
-   SUBROUTINE trd_trc( ptrtrd, kjn, ktrd, kt )
+   SUBROUTINE trd_trc( ptrtrd, kjn, ktrd, kt, Kmm )
       INTEGER ::   kt, kjn, ktrd   
+      INTEGER ::   Kmm            ! time level index
       REAL    ::   ptrtrd(:,:,:)  
       WRITE(*,*) 'trd_trc : You should not have seen this print! error?', ptrtrd(1,1,1)
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/AGE/trcsms_age.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/AGE/trcsms_age.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/AGE/trcsms_age.F90	(revision 10966)
@@ -37,5 +37,5 @@
 CONTAINS
 
-   SUBROUTINE trc_sms_age( kt )
+   SUBROUTINE trc_sms_age( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                     ***  trc_sms_age  ***
@@ -46,4 +46,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kmm  ! ocean time level
       INTEGER ::   jn, jk   ! dummy loop index
       !!----------------------------------------------------------------------
@@ -67,5 +68,5 @@
       END DO
       !
-      IF( l_trdtrc ) CALL trd_trc( tra(:,:,:,jp_age), jn, jptra_sms, kt )   ! save trends
+      IF( l_trdtrc ) CALL trd_trc( tra(:,:,:,jp_age), jn, jptra_sms, kt, Kmm )   ! save trends
       !
       IF( ln_timing )   CALL timing_stop('trc_sms_age')
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/C14/trcsms_c14.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/C14/trcsms_c14.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/C14/trcsms_c14.F90	(revision 10966)
@@ -33,5 +33,5 @@
 CONTAINS
 
-   SUBROUTINE trc_sms_c14( kt )
+   SUBROUTINE trc_sms_c14( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trc_sms_c14  ***
@@ -51,4 +51,5 @@
       !
       INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kmm   ! ocean time level
       !
       INTEGER  :: ji, jj, jk         ! dummy loop indices 
@@ -157,5 +158,5 @@
       ENDIF
 
-      IF( l_trdtrc )  CALL trd_trc( tra(:,:,:,jp_c14), 1, jptra_sms, kt )   ! save trends
+      IF( l_trdtrc )  CALL trd_trc( tra(:,:,:,jp_c14), 1, jptra_sms, kt, Kmm )   ! save trends
       !
       IF( ln_timing )   CALL timing_stop('trc_sms_c14')
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/CFC/trcsms_cfc.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/CFC/trcsms_cfc.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/CFC/trcsms_cfc.F90	(revision 10966)
@@ -54,5 +54,5 @@
 CONTAINS
 
-   SUBROUTINE trc_sms_cfc( kt )
+   SUBROUTINE trc_sms_cfc( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE trc_sms_cfc  ***
@@ -71,4 +71,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kmm   ! ocean time level
       !
       INTEGER  ::   ji, jj, jn, jl, jm
@@ -191,5 +192,5 @@
       IF( l_trdtrc ) THEN
           DO jn = jp_cfc0, jp_cfc1
-            CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt )   ! save trends
+            CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt, Kmm )   ! save trends
           END DO
       END IF
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/MY_TRC/trcsms_my_trc.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/MY_TRC/trcsms_my_trc.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/MY_TRC/trcsms_my_trc.F90	(revision 10966)
@@ -63,5 +63,5 @@
           DO jn = jp_myt0, jp_myt1
             ztrmyt(:,:,:) = tra(:,:,:,jn)
-            CALL trd_trc( ztrmyt, jn, jptra_sms, kt )   ! save trends
+            CALL trd_trc( ztrmyt, jn, jptra_sms, kt, Kmm )   ! save trends
           END DO
           DEALLOCATE( ztrmyt )
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P2Z/p2zsms.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P2Z/p2zsms.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P2Z/p2zsms.F90	(revision 10966)
@@ -35,5 +35,5 @@
 CONTAINS
 
-   SUBROUTINE p2z_sms( kt )
+   SUBROUTINE p2z_sms( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p2z_sms  ***
@@ -45,4 +45,5 @@
       !! --------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index      
+      INTEGER, INTENT( in ) ::   Kmm     ! ocean time level index      
       !
       INTEGER ::   jn   ! dummy loop index
@@ -58,5 +59,5 @@
       IF( l_trdtrc ) THEN
          DO jn = jp_pcs0, jp_pcs1
-           CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt )   ! save trends
+           CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt, Kmm )   ! save trends
          END DO
       END IF
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zbio.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zbio.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zbio.F90	(revision 10966)
@@ -45,5 +45,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_bio ( kt, knt )
+   SUBROUTINE p4z_bio ( kt, knt, Kbb, Kmm )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_bio  ***
@@ -56,4 +56,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) :: kt, knt
+      INTEGER, INTENT(in) :: Kbb, Kmm  ! time level indices
       !
       INTEGER             :: ji, jj, jk, jn
@@ -78,5 +79,5 @@
 
       CALL p4z_opt     ( kt, knt )     ! Optic: PAR in the water column
-      CALL p4z_sink    ( kt, knt )     ! vertical flux of particulate organic matter
+      CALL p4z_sink    ( kt, knt, Kbb, Kmm )     ! vertical flux of particulate organic matter
       CALL p4z_fechem  ( kt, knt )     ! Iron chemistry/scavenging
       !
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zsink.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zsink.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zsink.F90	(revision 10966)
@@ -49,5 +49,5 @@
    !!----------------------------------------------------------------------
 
-   SUBROUTINE p4z_sink ( kt, knt )
+   SUBROUTINE p4z_sink ( kt, knt, Kbb, Kmm )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_sink  ***
@@ -59,4 +59,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) :: kt, knt
+      INTEGER, INTENT(in) :: Kbb, Kmm  ! time level indices
       INTEGER  ::   ji, jj, jk
       CHARACTER (len=25) :: charout
@@ -104,10 +105,10 @@
       !   Compute the sedimentation term using p4zsink2 for all the sinking particles
       !   -----------------------------------------------------
-      CALL trc_sink( kt, wsbio3, sinking , jppoc, rfact2 )
-      CALL trc_sink( kt, wsbio3, sinkfer , jpsfe, rfact2 )
-      CALL trc_sink( kt, wsbio4, sinking2, jpgoc, rfact2 )
-      CALL trc_sink( kt, wsbio4, sinkfer2, jpbfe, rfact2 )
-      CALL trc_sink( kt, wsbio4, sinksil , jpgsi, rfact2 )
-      CALL trc_sink( kt, wsbio4, sinkcal , jpcal, rfact2 )
+      CALL trc_sink( kt, Kbb, Kmm, wsbio3, sinking , jppoc, rfact2 )
+      CALL trc_sink( kt, Kbb, Kmm, wsbio3, sinkfer , jpsfe, rfact2 )
+      CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinking2, jpgoc, rfact2 )
+      CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinkfer2, jpbfe, rfact2 )
+      CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinksil , jpgsi, rfact2 )
+      CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinkcal , jpcal, rfact2 )
 
       IF( ln_p5z ) THEN
@@ -119,8 +120,8 @@
          !   Compute the sedimentation term using p4zsink2 for all the sinking particles
          !   -----------------------------------------------------
-         CALL trc_sink( kt, wsbio3, sinkingn , jppon, rfact2 )
-         CALL trc_sink( kt, wsbio3, sinkingp , jppop, rfact2 )
-         CALL trc_sink( kt, wsbio4, sinking2n, jpgon, rfact2 )
-         CALL trc_sink( kt, wsbio4, sinking2p, jpgop, rfact2 )
+         CALL trc_sink( kt, Kbb, Kmm, wsbio3, sinkingn , jppon, rfact2 )
+         CALL trc_sink( kt, Kbb, Kmm, wsbio3, sinkingp , jppop, rfact2 )
+         CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinking2n, jpgon, rfact2 )
+         CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinking2p, jpgop, rfact2 )
       ENDIF
 
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zsms.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 10966)
@@ -46,5 +46,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_sms( kt, Kmm )
+   SUBROUTINE p4z_sms( kt, Kbb, Kmm )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_sms  ***
@@ -58,6 +58,6 @@
       !!---------------------------------------------------------------------
       !
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index      
-      INTEGER, INTENT( in ) ::   Kmm     ! time level index
+      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index      
+      INTEGER, INTENT( in ) ::   Kbb, Kmm     ! time level index
       !!
       INTEGER ::   ji, jj, jk, jnt, jn, jl
@@ -111,5 +111,5 @@
       DO jnt = 1, nrdttrc          ! Potential time splitting if requested
          !
-         CALL p4z_bio( kt, jnt )   ! Biology
+         CALL p4z_bio( kt, jnt, Kbb, Kmm )   ! Biology
          CALL p4z_lys( kt, jnt )   ! Compute CaCO3 saturation
          CALL p4z_sed( kt, jnt )   ! Surface and Bottom boundary conditions
@@ -149,5 +149,5 @@
       IF( l_trdtrc ) THEN
          DO jn = jp_pcs0, jp_pcs1
-           CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt )   ! save trends
+           CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt, Kmm )   ! save trends
          END DO
       END IF
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/trcsms_pisces.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/trcsms_pisces.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/trcsms_pisces.F90	(revision 10966)
@@ -25,5 +25,5 @@
 CONTAINS
 
-   SUBROUTINE trc_sms_pisces( kt, Kmm )
+   SUBROUTINE trc_sms_pisces( kt, Kbb, Kmm )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE trc_sms_pisces  ***
@@ -34,10 +34,10 @@
       !!---------------------------------------------------------------------
       !
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index      
-      INTEGER, INTENT( in ) ::   Kmm     ! time level index
+      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index      
+      INTEGER, INTENT( in ) ::   Kbb, Kmm     ! time level index
       !!---------------------------------------------------------------------
       !
-      IF( ln_p4z .OR. ln_p5z ) THEN  ;   CALL p4z_sms( kt, Kmm )   !  PISCES
-      ELSE                           ;   CALL p2z_sms( kt )   !  LOBSTER
+      IF( ln_p4z .OR. ln_p5z ) THEN  ;   CALL p4z_sms( kt, Kbb, Kmm )   !  PISCES
+      ELSE                           ;   CALL p2z_sms( kt,      Kmm )   !  LOBSTER
       ENDIF
 
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcadv.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcadv.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcadv.F90	(revision 10966)
@@ -74,5 +74,5 @@
       !! ** Purpose :   compute the ocean tracer advection trend.
       !!
-      !! ** Method  : - Update after tracers (tra) with the advection term following nadv
+      !! ** Method  : - Update after tracers (tr(Krhs)) with the advection term following nadv
       !!----------------------------------------------------------------------
       INTEGER                                   , INTENT(in)    :: kt   ! ocean time-step index
@@ -82,5 +82,5 @@
       INTEGER ::   jk   ! dummy loop index
       CHARACTER (len=22) ::   charout
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zun, zvn, zwn  ! effective velocity
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zuu, zvv, zww  ! effective velocity
       !!----------------------------------------------------------------------
       !
@@ -89,34 +89,34 @@
       !                                         !==  effective transport  ==!
       IF( l_offline ) THEN
-         zun(:,:,:) = un(:,:,:)                    ! already in (un,vn,wn)
-         zvn(:,:,:) = vn(:,:,:)
-         zwn(:,:,:) = wn(:,:,:)
+         zuu(:,:,:) = uu(:,:,:,Kmm)                ! already in (uu(Kmm),vv(Kmm),ww)
+         zvv(:,:,:) = vv(:,:,:,Kmm)
+         zww(:,:,:) = ww(:,:,:)
       ELSE                                         ! build the effective transport
-         zun(:,:,jpk) = 0._wp
-         zvn(:,:,jpk) = 0._wp
-         zwn(:,:,jpk) = 0._wp
+         zuu(:,:,jpk) = 0._wp
+         zvv(:,:,jpk) = 0._wp
+         zww(:,:,jpk) = 0._wp
          IF( ln_wave .AND. ln_sdw )  THEN
             DO jk = 1, jpkm1                                                       ! eulerian transport + Stokes Drift
-               zun(:,:,jk) = e2u  (:,:) * e3u_n(:,:,jk) * ( un(:,:,jk) + usd(:,:,jk) )
-               zvn(:,:,jk) = e1v  (:,:) * e3v_n(:,:,jk) * ( vn(:,:,jk) + vsd(:,:,jk) )
-               zwn(:,:,jk) = e1e2t(:,:)                 * ( wn(:,:,jk) + wsd(:,:,jk) )
+               zuu(:,:,jk) = e2u  (:,:) * e3u(:,:,jk,Kmm) * ( uu(:,:,jk,Kmm) + usd(:,:,jk) )
+               zvv(:,:,jk) = e1v  (:,:) * e3v(:,:,jk,Kmm) * ( vv(:,:,jk,Kmm) + vsd(:,:,jk) )
+               zww(:,:,jk) = e1e2t(:,:)                   * ( ww(:,:,jk) + wsd(:,:,jk) )
             END DO
          ELSE
             DO jk = 1, jpkm1
-               zun(:,:,jk) = e2u  (:,:) * e3u_n(:,:,jk) * un(:,:,jk)                   ! eulerian transport
-               zvn(:,:,jk) = e1v  (:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
-               zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk)
+               zuu(:,:,jk) = e2u  (:,:) * e3u(:,:,jk,Kmm) * uu(:,:,jk,Kmm)                   ! eulerian transport
+               zvv(:,:,jk) = e1v  (:,:) * e3v(:,:,jk,Kmm) * vv(:,:,jk,Kmm)
+               zww(:,:,jk) = e1e2t(:,:)                   * ww(:,:,jk)
             END DO
          ENDIF
          !
          IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN                                 ! add z-tilde and/or vvl corrections
-            zun(:,:,:) = zun(:,:,:) + un_td(:,:,:)
-            zvn(:,:,:) = zvn(:,:,:) + vn_td(:,:,:)
+            zuu(:,:,:) = zuu(:,:,:) + un_td(:,:,:)
+            zvv(:,:,:) = zvv(:,:,:) + vn_td(:,:,:)
          ENDIF
          !
          IF( ln_ldfeiv .AND. .NOT. ln_traldf_triad )   & 
-            &              CALL ldf_eiv_trp( kt, nittrc000, zun, zvn, zwn, 'TRC', Kmm, Krhs )  ! add the eiv transport
-         !
-         IF( ln_mle    )   CALL tra_mle_trp( kt, nittrc000, zun, zvn, zwn, 'TRC', Kmm       )  ! add the mle transport
+            &              CALL ldf_eiv_trp( kt, nittrc000, zuu, zvv, zww, 'TRC', Kmm, Krhs )  ! add the eiv transport
+         !
+         IF( ln_mle    )   CALL tra_mle_trp( kt, nittrc000, zuu, zvv, zww, 'TRC', Kmm       )  ! add the mle transport
          !
       ENDIF
@@ -125,13 +125,13 @@
       !
       CASE ( np_CEN )                                 ! Centered : 2nd / 4th order
-         CALL tra_adv_cen( kt, nittrc000,'TRC',          zun, zvn, zwn,      Kmm, ptr, jptra, Krhs, nn_cen_h, nn_cen_v )
+         CALL tra_adv_cen( kt, nittrc000,'TRC',          zuu, zvv, zww,      Kmm, ptr, jptra, Krhs, nn_cen_h, nn_cen_v )
       CASE ( np_FCT )                                 ! FCT      : 2nd / 4th order
-         CALL tra_adv_fct( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, Kbb, Kmm, ptr, jptra, Krhs, nn_fct_h, nn_fct_v )
+         CALL tra_adv_fct( kt, nittrc000,'TRC', r2dttrc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, nn_fct_h, nn_fct_v )
       CASE ( np_MUS )                                 ! MUSCL
-         CALL tra_adv_mus( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, Kbb, Kmm, ptr, jptra, Krhs, ln_mus_ups         ) 
+         CALL tra_adv_mus( kt, nittrc000,'TRC', r2dttrc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, ln_mus_ups         ) 
       CASE ( np_UBS )                                 ! UBS
-         CALL tra_adv_ubs( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, Kbb, Kmm, ptr, jptra, Krhs, nn_ubs_v           )
+         CALL tra_adv_ubs( kt, nittrc000,'TRC', r2dttrc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, nn_ubs_v           )
       CASE ( np_QCK )                                 ! QUICKEST
-         CALL tra_adv_qck( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, Kbb, Kmm, ptr, jptra, Krhs                     )
+         CALL tra_adv_qck( kt, nittrc000,'TRC', r2dttrc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs                     )
       !
       END SELECT
@@ -140,5 +140,5 @@
          WRITE(charout, FMT="('adv ')")
          CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         CALL prt_ctl_trc( tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
       END IF
       !
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcbbl.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcbbl.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcbbl.F90	(revision 10966)
@@ -61,5 +61,5 @@
       IF( l_trdtrc )  THEN
          ALLOCATE( ztrtrd(jpi,jpj,jpk,jptra) ) ! temporary save of trends
-         ztrtrd(:,:,:,:)  = tra(:,:,:,:)
+         ztrtrd(:,:,:,:)  = tr(:,:,:,:,Krhs)
       ENDIF
 
@@ -67,8 +67,8 @@
       IF( nn_bbl_ldf == 1 ) THEN
          !
-         CALL tra_bbl_dif( trb, tra, jptra, Kmm )  
+         CALL tra_bbl_dif( tr(:,:,:,:,Kbb), tr(:,:,:,:,Krhs), jptra, Kmm )  
          IF( ln_ctl )   THEN
             WRITE(charout, FMT="(' bbl_dif')")  ;  CALL prt_ctl_trc_info(charout)
-            CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+            CALL prt_ctl_trc( tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
          ENDIF
          !
@@ -78,8 +78,8 @@
       IF( nn_bbl_adv /= 0 ) THEN
          !
-         CALL tra_bbl_adv( trb, tra, jptra, Kmm )  
+         CALL tra_bbl_adv( tr(:,:,:,:,Kbb), tr(:,:,:,:,Krhs), jptra, Kmm )  
          IF( ln_ctl )   THEN
             WRITE(charout, FMT="(' bbl_adv')")  ;  CALL prt_ctl_trc_info(charout)
-            CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+            CALL prt_ctl_trc( tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
          ENDIF
          !
@@ -88,5 +88,5 @@
       IF( l_trdtrc )   THEN                      ! save the horizontal diffusive trends for further diagnostics
         DO jn = 1, jptra
-           ztrtrd(:,:,:,jn) = tra(:,:,:,jn) - ztrtrd(:,:,:,jn)
+           ztrtrd(:,:,:,jn) = tr(:,:,:,jn,Krhs) - ztrtrd(:,:,:,jn)
            CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_bbl, ztrtrd(:,:,:,jn) )
         END DO
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcdmp.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcdmp.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcdmp.F90	(revision 10966)
@@ -63,5 +63,5 @@
 
 
-   SUBROUTINE trc_dmp( kt, Kmm, Krhs )
+   SUBROUTINE trc_dmp( kt, Kbb, Kmm, Krhs )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE trc_dmp  ***
@@ -73,15 +73,15 @@
       !! ** Method  :   Newtonian damping towards trdta computed 
       !!      and add to the general tracer trends:
-      !!                     trn = tra + restotr * (trdta - trb)
+      !!                     tr(Kmm) = tr(Krhs) + restotr * (trdta - tr(Kbb))
       !!         The trend is computed either throughout the water column
       !!      (nlmdmptr=0) or in area of weak vertical mixing (nlmdmptr=1) or
       !!      below the well mixed layer (nlmdmptr=2)
       !!
-      !! ** Action  : - update the tracer trends tra with the newtonian 
+      !! ** Action  : - update the tracer trends tr(:,:,:,:,Krhs) with the newtonian 
       !!                damping trends.
       !!              - save the trends ('key_trdmxl_trc')
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      INTEGER, INTENT(in) ::   Kmm, Krhs  ! time level indices
+      INTEGER, INTENT(in) ::   kt              ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! time level indices
       !
       INTEGER ::   ji, jj, jk, jn, jl   ! dummy loop indices
@@ -101,5 +101,5 @@
          DO jn = 1, jptra                                           ! tracer loop
             !                                                       ! ===========
-            IF( l_trdtrc ) ztrtrd(:,:,:) = tra(:,:,:,jn)    ! save trends 
+            IF( l_trdtrc ) ztrtrd(:,:,:) = tr(:,:,:,jn,Krhs)    ! save trends 
             !
             IF( ln_trc_ini(jn) ) THEN      ! update passive tracers arrays with input data read from file
@@ -114,5 +114,5 @@
                      DO jj = 2, jpjm1
                         DO ji = fs_2, fs_jpim1   ! vector opt.
-                           tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + restotr(ji,jj,jk) * ( ztrcdta(ji,jj,jk) - trb(ji,jj,jk,jn) )
+                           tr(ji,jj,jk,jn,Krhs) = tr(ji,jj,jk,jn,Krhs) + restotr(ji,jj,jk) * ( ztrcdta(ji,jj,jk) - tr(ji,jj,jk,jn,Kbb) )
                         END DO
                      END DO
@@ -124,5 +124,5 @@
                         DO ji = fs_2, fs_jpim1   ! vector opt.
                            IF( avt(ji,jj,jk) <= avt_c )  THEN 
-                              tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + restotr(ji,jj,jk) * ( ztrcdta(ji,jj,jk) - trb(ji,jj,jk,jn) )
+                              tr(ji,jj,jk,jn,Krhs) = tr(ji,jj,jk,jn,Krhs) + restotr(ji,jj,jk) * ( ztrcdta(ji,jj,jk) - tr(ji,jj,jk,jn,Kbb) )
                            ENDIF
                         END DO
@@ -134,6 +134,6 @@
                      DO jj = 2, jpjm1
                         DO ji = fs_2, fs_jpim1   ! vector opt.
-                           IF( gdept_n(ji,jj,jk) >= hmlp (ji,jj) ) THEN
-                              tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + restotr(ji,jj,jk) * ( ztrcdta(ji,jj,jk) - trb(ji,jj,jk,jn) )
+                           IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN
+                              tr(ji,jj,jk,jn,Krhs) = tr(ji,jj,jk,jn,Krhs) + restotr(ji,jj,jk) * ( ztrcdta(ji,jj,jk) - tr(ji,jj,jk,jn,Kbb) )
                            END IF
                         END DO
@@ -146,5 +146,5 @@
             !
             IF( l_trdtrc ) THEN
-               ztrtrd(:,:,:) = tra(:,:,:,jn) -  ztrtrd(:,:,:)
+               ztrtrd(:,:,:) = tr(:,:,:,jn,Krhs) -  ztrtrd(:,:,:)
                CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_dmp, ztrtrd )
             END IF
@@ -160,5 +160,5 @@
          WRITE(charout, FMT="('dmp ')")
          CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         CALL prt_ctl_trc( tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
       ENDIF
       !
@@ -225,5 +225,5 @@
 
 
-   SUBROUTINE trc_dmp_clo( kt, Kmm )
+   SUBROUTINE trc_dmp_clo( kt, Kbb, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trc_dmp_clo  ***
@@ -237,6 +237,6 @@
       !!                nctsi2(), nctsj2() : north-east Closed sea limits (i,j)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-      INTEGER, INTENT( in ) ::   Kmm     ! time level indices
+      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kbb, Kmm     ! time level indices
       !
       INTEGER :: ji , jj, jk, jn, jl, jc                    ! dummy loop indicesa
@@ -361,6 +361,6 @@
                       DO jj = nctsj1(jc), nctsj2(jc)
                          DO ji = nctsi1(jc), nctsi2(jc)
-                            trn(ji,jj,jk,jn) = ztrcdta(ji,jj,jk)
-                            trb(ji,jj,jk,jn) = trn(ji,jj,jk,jn)
+                            tr(ji,jj,jk,jn,Kmm) = ztrcdta(ji,jj,jk)
+                            tr(ji,jj,jk,jn,Kbb) = tr(ji,jj,jk,jn,Kmm)
                          END DO
                       END DO
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcldf.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcldf.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcldf.F90	(revision 10966)
@@ -51,5 +51,5 @@
 CONTAINS
 
-   SUBROUTINE trc_ldf( kt, Kmm, Krhs )
+   SUBROUTINE trc_ldf( kt, Kbb, Kmm, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf  ***
@@ -58,6 +58,6 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
-      INTEGER, INTENT( in ) ::   Kmm, Krhs  ! ocean time-level index
+      INTEGER, INTENT( in ) ::   kt              ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kbb, Kmm, Krhs  ! ocean time-level index
       !
       INTEGER            :: ji, jj, jk, jn
@@ -74,5 +74,5 @@
       IF( l_trdtrc )  THEN
          ALLOCATE( ztrtrd(jpi,jpj,jpk,jptra) )
-         ztrtrd(:,:,:,:)  = tra(:,:,:,:)
+         ztrtrd(:,:,:,:)  = tr(:,:,:,:,Krhs)
       ENDIF
       !                                  !* set the lateral diffusivity coef. for passive tracer      
@@ -83,6 +83,6 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( gdept_n(ji,jj,jk) > 200. .AND. gphit(ji,jj) < 5. .AND. gphit(ji,jj) > -5. ) THEN
-                  zdep = MAX( gdept_n(ji,jj,jk) - 1000., 0. ) / 1000.
+               IF( gdept(ji,jj,jk,Kmm) > 200. .AND. gphit(ji,jj) < 5. .AND. gphit(ji,jj) > -5. ) THEN
+                  zdep = MAX( gdept(ji,jj,jk,Kmm) - 1000., 0. ) / 1000.
                   zahu(ji,jj,jk) = zahu(ji,jj,jk) * MAX( 1., rn_fact_lap * EXP( -zdep ) )
                ENDIF
@@ -93,17 +93,21 @@
       SELECT CASE ( nldf_trc )                 !* compute lateral mixing trend and add it to the general trend
       !
-      CASE ( np_lap   )                               ! iso-level laplacian
-         CALL tra_ldf_lap  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb,      tra, jptra,    1    , Kmm )
-      CASE ( np_lap_i )                               ! laplacian : standard iso-neutral operator (Madec)
-         CALL tra_ldf_iso  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb, trb, tra, jptra,    1    , Kmm )
-      CASE ( np_lap_it )                              ! laplacian : triad iso-neutral operator (griffies)
-         CALL tra_ldf_triad( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb, trb, tra, jptra,    1    , Kmm )
-      CASE ( np_blp , np_blp_i , np_blp_it )          ! bilaplacian: all operator (iso-level, -neutral)
-         CALL tra_ldf_blp  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb     , tra, jptra, nldf_trc, Kmm )
+      CASE ( np_lap   )                                                                                    ! iso-level laplacian
+         CALL tra_ldf_lap  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi,               &
+           &                     tr(:,:,:,:,Kbb), tr(:,:,:,:,Krhs),                 jptra, 1 , Kmm )
+      CASE ( np_lap_i )                                                                                    ! laplacian : standard iso-neutral operator (Madec)
+         CALL tra_ldf_iso  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi,               &
+           &                     tr(:,:,:,:,Kbb), tr(:,:,:,:,Kbb), tr(:,:,:,:,Krhs), jptra, 1 , Kmm )
+      CASE ( np_lap_it )                                                                                   ! laplacian : triad iso-neutral operator (griffies)
+         CALL tra_ldf_triad( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi,               &
+           &                     tr(:,:,:,:,Kbb), tr(:,:,:,:,Kbb), tr(:,:,:,:,Krhs), jptra, 1 , Kmm )
+      CASE ( np_blp , np_blp_i , np_blp_it )                                                               ! bilaplacian: all operator (iso-level, -neutral)
+         CALL tra_ldf_blp  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi,               &
+           &                     tr(:,:,:,:,Kbb) , tr(:,:,:,:,Krhs),                 jptra, nldf_trc, Kmm )
       END SELECT
       !
       IF( l_trdtrc )   THEN                    ! send the trends for further diagnostics
         DO jn = 1, jptra
-           ztrtrd(:,:,:,jn) = tra(:,:,:,jn) - ztrtrd(:,:,:,jn)
+           ztrtrd(:,:,:,jn) = tr(:,:,:,jn,Krhs) - ztrtrd(:,:,:,jn)
            CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_ldf, ztrtrd(:,:,:,jn) )
         END DO
@@ -114,5 +118,5 @@
          WRITE(charout, FMT="('ldf ')")
          CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         CALL prt_ctl_trc( tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
       ENDIF
       !
@@ -168,5 +172,5 @@
       IF( ln_trcldf_OFF  ) THEN   ;   nldf_trc = np_no_ldf   ;   ioptio = ioptio + 1   ;   ENDIF
       IF( ln_trcldf_tra  ) THEN   ;   nldf_trc = nldf_tra    ;   ioptio = ioptio + 1   ;   ENDIF
-      IF( ioptio /=  1   )   CALL ctl_stop( 'trc_ldf_ini: use ONE of the 2 operator options (NONE/tra)' )
+      IF( ioptio /=  1   )   CALL ctl_stop( 'trc_ldf_ini: use ONE of the 2 operator options (OFF/tra)' )
       
       !                                ! multiplier : passive/active tracers ration
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcnxt.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcnxt.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcnxt.F90	(revision 10966)
@@ -61,20 +61,20 @@
       !!      next time-step from their temporal trends and swap the fields.
       !! 
-      !! ** Method  :   Apply lateral boundary conditions on (ua,va) through 
+      !! ** Method  :   Apply lateral boundary conditions on (uu(Krhs),vv(Krhs)) through 
       !!      call to lbc_lnk routine
       !!   default:
       !!      arrays swap
-      !!         (trn) = (tra) ; (tra) = (0,0)
-      !!         (trb) = (trn) 
+      !!         (tr(Kmm)) = (tr(Krhs)) ; (tr(Krhs)) = (0,0)
+      !!         (tr(Kbb)) = (tr(Kmm)) 
       !!
       !!   For Arakawa or TVD Scheme : 
-      !!      A Asselin time filter applied on now tracers (trn) to avoid
+      !!      A Asselin time filter applied on now tracers (tr(:,:,:,:,Kmm)) to avoid
       !!      the divergence of two consecutive time-steps and tr arrays
       !!      to prepare the next time_step:
-      !!         (trb) = (trn) + atfp [ (trb) + (tra) - 2 (trn) ]
-      !!         (trn) = (tra) ; (tra) = (0,0)
-      !!
-      !!
-      !! ** Action  : - update trb, trn
+      !!         (tr(Kbb)) = (tr(Kmm)) + atfp [ (tr(Kbb)) + (tr(Krhs)) - 2 (tr(Kmm)) ]
+      !!         (tr(Kmm)) = (tr(Krhs)) ; (tr(Krhs)) = (0,0)
+      !!
+      !!
+      !! ** Action  : - update tr(Kbb), tr(Kmm)
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
@@ -98,5 +98,5 @@
 #endif
       ! Update after tracer on domain lateral boundaries
-      CALL lbc_lnk( 'trcnxt', tra(:,:,:,:), 'T', 1. )   
+      CALL lbc_lnk( 'trcnxt', tr(:,:,:,:,Krhs), 'T', 1. )   
 
       IF( ln_bdy )  CALL trc_bdy( kt, Kbb, Kmm, Krhs )
@@ -113,9 +113,9 @@
          ! total trend for the non-time-filtered variables. 
          zfact = 1.0 / rdttrc
-         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from tsn terms
+         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from ts(Kmm) terms
          IF( ln_linssh ) THEN       ! linear sea surface height only
             DO jn = 1, jptra
                DO jk = 1, jpkm1
-                  ztrdt(:,:,jk,jn) = ( tra(:,:,jk,jn)*e3t_a(:,:,jk) / e3t_n(:,:,jk) - trn(:,:,jk,jn)) * zfact
+                  ztrdt(:,:,jk,jn) = ( tr(:,:,jk,jn,Krhs)*e3t(:,:,jk,Krhs) / e3t(:,:,jk,Kmm) - tr(:,:,jk,jn,Kmm)) * zfact
                END DO
             END DO
@@ -123,5 +123,5 @@
             DO jn = 1, jptra
                DO jk = 1, jpkm1
-                  ztrdt(:,:,jk,jn) = ( tra(:,:,jk,jn) - trn(:,:,jk,jn) ) * zfact
+                  ztrdt(:,:,jk,jn) = ( tr(:,:,jk,jn,Krhs) - tr(:,:,jk,jn,Kmm) ) * zfact
                END DO
             END DO
@@ -135,5 +135,5 @@
             ! Store now fields before applying the Asselin filter 
             ! in order to calculate Asselin filter trend later.
-            ztrdt(:,:,:,:) = trn(:,:,:,:) 
+            ztrdt(:,:,:,:) = tr(:,:,:,:,Kmm) 
          ENDIF
 
@@ -143,6 +143,6 @@
          DO jn = 1, jptra
             DO jk = 1, jpkm1
-               trn(:,:,jk,jn) = tra(:,:,jk,jn)
-               trb(:,:,jk,jn) = trn(:,:,jk,jn)  
+               tr(:,:,jk,jn,Kmm) = tr(:,:,jk,jn,Krhs)
+               tr(:,:,jk,jn,Kbb) = tr(:,:,jk,jn,Kmm)  
             END DO
          END DO
@@ -157,13 +157,15 @@
       ELSE     
          IF( .NOT. l_offline ) THEN ! Leap-Frog + Asselin filter time stepping
-            IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt,      Kmm,       nittrc000,         'TRC', trb, trn, tra, jptra )  !     linear ssh
-            ELSE                   ;   CALL tra_nxt_vvl( kt, Kbb, Kmm, Krhs, nittrc000, rdttrc, 'TRC', trb, trn, tra,      &
-              &                                                                                   sbc_trc, sbc_trc_b, jptra )  ! non-linear ssh
+            IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt,      Kmm,       nittrc000,         'TRC',                &
+              &                                              tr(:,:,:,:,Kbb), tr(:,:,:,:,Kmm), tr(:,:,:,:,Krhs), jptra )  !     linear ssh
+            ELSE                   ;   CALL tra_nxt_vvl( kt, Kbb, Kmm, Krhs, nittrc000, rdttrc, 'TRC',                &
+              &                                              tr(:,:,:,:,Kbb), tr(:,:,:,:,Kmm), tr(:,:,:,:,Krhs),      &
+              &                                              sbc_trc, sbc_trc_b, jptra )                                  ! non-linear ssh
             ENDIF
          ELSE
-                                       CALL trc_nxt_off( kt )       ! offline 
-         ENDIF
-         !
-         CALL lbc_lnk_multi( 'trcnxt', trb(:,:,:,:), 'T', 1._wp, trn(:,:,:,:), 'T', 1._wp, tra(:,:,:,:), 'T', 1._wp )
+                                       CALL trc_nxt_off( kt, Kbb, Kmm, Krhs )       ! offline 
+         ENDIF
+         !
+         CALL lbc_lnk_multi( 'trcnxt', tr(:,:,:,:,Kbb), 'T', 1._wp, tr(:,:,:,:,Kmm), 'T', 1._wp, tr(:,:,:,:,Krhs), 'T', 1._wp )
       ENDIF
       !
@@ -172,5 +174,5 @@
             DO jk = 1, jpkm1
                zfact = 1._wp / r2dttrc  
-               ztrdt(:,:,jk,jn) = ( trb(:,:,jk,jn) - ztrdt(:,:,jk,jn) ) * zfact 
+               ztrdt(:,:,jk,jn) = ( tr(:,:,jk,jn,Kbb) - ztrdt(:,:,jk,jn) ) * zfact 
             END DO
             CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
@@ -182,5 +184,5 @@
          WRITE(charout, FMT="('nxt')")
          CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=trn, mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Kmm), mask=tmask, clinfo=ctrcnm)
       ENDIF
       !
@@ -190,5 +192,5 @@
 
 
-   SUBROUTINE trc_nxt_off( kt )
+   SUBROUTINE trc_nxt_off( kt, Kbb, Kmm, Krhs )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tra_nxt_vvl  ***
@@ -204,8 +206,8 @@
       !!                This can be summurized for tempearture as:
       !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
-      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
+      !!                  /( e3t(:,:,:,Kmm)    + rbcp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Krhs)    ] )   
       !!             ztm = 0                                                       otherwise
       !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
-      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
+      !!                  /( e3t(:,:,:,Kmm)    + atfp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Krhs)    ] )
       !!             tn  = ta 
       !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
@@ -214,5 +216,6 @@
       !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
       !!----------------------------------------------------------------------
-      INTEGER , INTENT(in   )   ::  kt       ! ocean time-step index
+      INTEGER, INTENT(in   )   ::  kt              ! ocean time-step index
+      INTEGER, INTENT(in   )   ::  Kbb, Kmm, Krhs  ! time level indices
       !!     
       INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
@@ -236,11 +239,11 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ze3t_b = e3t_b(ji,jj,jk)
-                  ze3t_n = e3t_n(ji,jj,jk)
-                  ze3t_a = e3t_a(ji,jj,jk)
+                  ze3t_b = e3t(ji,jj,jk,Kbb)
+                  ze3t_n = e3t(ji,jj,jk,Kmm)
+                  ze3t_a = e3t(ji,jj,jk,Krhs)
                   !                                         ! tracer content at Before, now and after
-                  ztc_b  = trb(ji,jj,jk,jn) * ze3t_b
-                  ztc_n  = trn(ji,jj,jk,jn) * ze3t_n
-                  ztc_a  = tra(ji,jj,jk,jn) * ze3t_a
+                  ztc_b  = tr(ji,jj,jk,jn,Kbb)  * ze3t_b
+                  ztc_n  = tr(ji,jj,jk,jn,Kmm)  * ze3t_n
+                  ztc_a  = tr(ji,jj,jk,jn,Krhs) * ze3t_a
                   !
                   ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
@@ -256,6 +259,6 @@
 
                   ze3t_f = 1.e0 / ze3t_f
-                  trb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
-                  trn(ji,jj,jk,jn) = tra(ji,jj,jk,jn)     ! ptn <-- pta
+                  tr(ji,jj,jk,jn,Kbb) = ztc_f * ze3t_f           ! pt(:,:,:,:,Kbb) <-- pt(:,:,:,:,Kmm) filtered
+                  tr(ji,jj,jk,jn,Kmm) = tr(ji,jj,jk,jn,Krhs)     ! pt(:,:,:,:,Kmm) <-- pt(:,:,:,:,Krhs)
                   !
                END DO
@@ -272,6 +275,7 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE trc_nxt( kt )  
+   SUBROUTINE trc_nxt( kt, Kbb, Kmm, Krhs )  
       INTEGER, INTENT(in) :: kt
+      INTEGER, INTENT(in) :: Kbb, Kmm, Krhs  ! time level indices
       WRITE(*,*) 'trc_nxt: You should not have seen this print! error?', kt
    END SUBROUTINE trc_nxt
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcrad.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcrad.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcrad.F90	(revision 10966)
@@ -37,5 +37,5 @@
 CONTAINS
 
-   SUBROUTINE trc_rad( kt, Kmm, Krhs )
+   SUBROUTINE trc_rad( kt, Kbb, Kmm, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trc_rad  ***
@@ -52,6 +52,6 @@
       !!                (the total CFC content is not strictly preserved)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      INTEGER, INTENT(in) ::   Kmm, Krhs  ! time level indices
+      INTEGER, INTENT(in) ::   kt              ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! time level indices
       !
       CHARACTER (len=22) :: charout
@@ -60,14 +60,14 @@
       IF( ln_timing )   CALL timing_start('trc_rad')
       !
-      IF( ln_age     )   CALL trc_rad_sms( kt, Kmm, Krhs, trb, trn, jp_age , jp_age                )  !  AGE
-      IF( ll_cfc     )   CALL trc_rad_sms( kt, Kmm, Krhs, trb, trn, jp_cfc0, jp_cfc1               )  !  CFC model
-      IF( ln_c14     )   CALL trc_rad_sms( kt, Kmm, Krhs, trb, trn, jp_c14 , jp_c14                )  !  C14
-      IF( ln_pisces  )   CALL trc_rad_sms( kt, Kmm, Krhs, trb, trn, jp_pcs0, jp_pcs1, cpreserv='Y' )  !  PISCES model
-      IF( ln_my_trc  )   CALL trc_rad_sms( kt, Kmm, Krhs, trb, trn, jp_myt0, jp_myt1               )  !  MY_TRC model
+      IF( ln_age     )   CALL trc_rad_sms( kt, Kmm, Krhs, tr(:,:,:,:,Kbb), tr(:,:,:,:,Kmm), jp_age , jp_age                )  !  AGE
+      IF( ll_cfc     )   CALL trc_rad_sms( kt, Kmm, Krhs, tr(:,:,:,:,Kbb), tr(:,:,:,:,Kmm), jp_cfc0, jp_cfc1               )  !  CFC model
+      IF( ln_c14     )   CALL trc_rad_sms( kt, Kmm, Krhs, tr(:,:,:,:,Kbb), tr(:,:,:,:,Kmm), jp_c14 , jp_c14                )  !  C14
+      IF( ln_pisces  )   CALL trc_rad_sms( kt, Kmm, Krhs, tr(:,:,:,:,Kbb), tr(:,:,:,:,Kmm), jp_pcs0, jp_pcs1, cpreserv='Y' )  !  PISCES model
+      IF( ln_my_trc  )   CALL trc_rad_sms( kt, Kmm, Krhs, tr(:,:,:,:,Kbb), tr(:,:,:,:,Kmm), jp_myt0, jp_myt1               )  !  MY_TRC model
       !
       IF(ln_ctl) THEN      ! print mean trends (used for debugging)
          WRITE(charout, FMT="('rad')")
          CALL prt_ctl_trc_info( charout )
-         CALL prt_ctl_trc( tab4d=trn, mask=tmask, clinfo=ctrcnm )
+         CALL prt_ctl_trc( tab4d=tr(:,:,:,:,Kmm), mask=tmask, clinfo=ctrcnm )
       ENDIF
       !
@@ -160,5 +160,5 @@
          DO jn = jp_sms0, jp_sms1
             !
-            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrb(:,:,:,jn)                            ! save input trb for trend computation           
+            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrb(:,:,:,jn)                            ! save input tr(:,:,:,:,Kbb) for trend computation           
             !
             DO jk = 1, jpkm1
@@ -210,5 +210,5 @@
          DO jn = jp_sms0, jp_sms1
             !
-            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrn(:,:,:,jn)                            ! save input trb for trend computation
+            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrn(:,:,:,jn)                            ! save input tr for trend computation
             !
             DO jk = 1, jpkm1
@@ -257,5 +257,5 @@
          DO jn = jp_sms0, jp_sms1  
             !
-            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrb(:,:,:,jn)                        ! save input trb for trend computation
+            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrb(:,:,:,jn)                        ! save input tr for trend computation
             !
             WHERE( ptrb(:,:,:,jn) < 0. )   ptrb(:,:,:,jn) = 0.
@@ -266,5 +266,5 @@
             ENDIF
             !
-            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrn(:,:,:,jn)                        ! save input trn for trend computation
+            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrn(:,:,:,jn)                        ! save input tr for trend computation
             !
             WHERE( ptrn(:,:,:,jn) < 0. )   ptrn(:,:,:,jn) = 0.
@@ -288,6 +288,7 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE trc_rad( kt )              ! Empty routine
+   SUBROUTINE trc_rad( kt, Kbb, Kmm, Krhs )              ! Empty routine
       INTEGER, INTENT(in) ::   kt
+      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! time level indices
       WRITE(*,*) 'trc_rad: You should not have seen this print! error?', kt
    END SUBROUTINE trc_rad
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcsbc.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcsbc.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcsbc.F90	(revision 10966)
@@ -49,10 +49,10 @@
       !!            The surface freshwater flux modify the ocean volume
       !!         and thus the concentration of a tracer as :
-      !!            tra = tra + emp * trn / e3t   for k=1
+      !!            tr(Krhs) = tr(Krhs) + emp * tr(Kmm) / e3t   for k=1
       !!         where emp, the surface freshwater budget (evaporation minus
       !!         precipitation ) given in kg/m2/s is divided
       !!         by 1035 kg/m3 (density of ocean water) to obtain m/s.
       !!
-      !! ** Action  : - Update the 1st level of tra with the trend associated
+      !! ** Action  : - Update the 1st level of tr(:,:,:,:,Krhs) with the trend associated
       !!                with the tracer surface boundary condition 
       !!
@@ -103,5 +103,5 @@
       ENDIF
 
-      ! Coupling online : river runoff is added to the horizontal divergence (hdivn) in the subroutine sbc_rnf_div 
+      ! Coupling online : river runoff is added to the horizontal divergence (hdiv) in the subroutine sbc_rnf_div 
       ! one only consider the concentration/dilution effect due to evaporation minus precipitation + freezing/melting of sea-ice
       ! Coupling offline : runoff are in emp which contains E-P-R
@@ -119,5 +119,5 @@
             DO jj = 2, jpj
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  sbc_trc(ji,jj,jn) = zsfx(ji,jj) * r1_rau0 * trn(ji,jj,1,jn)
+                  sbc_trc(ji,jj,jn) = zsfx(ji,jj) * r1_rau0 * tr(ji,jj,1,jn,Kmm)
                END DO
             END DO
@@ -129,5 +129,5 @@
             DO jj = 2, jpj
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zse3t = 1. / e3t_n(ji,jj,1)
+                  zse3t = 1. / e3t(ji,jj,1,Kmm)
                   ! tracer flux at the ice/ocean interface (tracer/m2/s)
                   zftra = - trc_i(ji,jj,jn) * fmmflx(ji,jj) ! uptake of tracer in the sea ice
@@ -138,7 +138,7 @@
                   ztfx  = zftra                             ! net tracer flux
                   !
-                  zdtra = r1_rau0 * ( ztfx + zsfx(ji,jj) * trn(ji,jj,1,jn) ) 
+                  zdtra = r1_rau0 * ( ztfx + zsfx(ji,jj) * tr(ji,jj,1,jn,Kmm) ) 
                   IF ( zdtra < 0. ) THEN
-                     zratio = -zdtra * zse3t * r2dttrc / ( trn(ji,jj,1,jn) + zrtrn )
+                     zratio = -zdtra * zse3t * r2dttrc / ( tr(ji,jj,1,jn,Kmm) + zrtrn )
                      zdtra = MIN(1.0, zratio) * zdtra ! avoid negative concentrations to arise
                   ENDIF
@@ -153,15 +153,15 @@
       DO jn = 1, jptra
          !
-         IF( l_trdtrc )   ztrtrd(:,:,:) = tra(:,:,:,jn)  ! save trends
+         IF( l_trdtrc )   ztrtrd(:,:,:) = tr(:,:,:,jn,Krhs)  ! save trends
          !
          DO jj = 2, jpj
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zse3t = zfact / e3t_n(ji,jj,1)
-               tra(ji,jj,1,jn) = tra(ji,jj,1,jn) + ( sbc_trc_b(ji,jj,jn) + sbc_trc(ji,jj,jn) ) * zse3t
+               zse3t = zfact / e3t(ji,jj,1,Kmm)
+               tr(ji,jj,1,jn,Krhs) = tr(ji,jj,1,jn,Krhs) + ( sbc_trc_b(ji,jj,jn) + sbc_trc(ji,jj,jn) ) * zse3t
             END DO
          END DO
          !
          IF( l_trdtrc ) THEN
-            ztrtrd(:,:,:) = tra(:,:,:,jn) - ztrtrd(:,:,:)
+            ztrtrd(:,:,:) = tr(:,:,:,jn,Krhs) - ztrtrd(:,:,:)
             CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_nsr, ztrtrd )
          END IF
@@ -184,5 +184,5 @@
       IF( ln_ctl )   THEN
          WRITE(charout, FMT="('sbc ')") ;  CALL prt_ctl_trc_info(charout)
-                                           CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+                                           CALL prt_ctl_trc( tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
       ENDIF
       IF( l_trdtrc )  DEALLOCATE( ztrtrd )
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcsink.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcsink.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcsink.F90	(revision 10966)
@@ -35,5 +35,5 @@
    !!----------------------------------------------------------------------
 
-   SUBROUTINE trc_sink ( kt, pwsink, psinkflx, jp_tra, rsfact )
+   SUBROUTINE trc_sink ( kt, Kbb, Kmm, pwsink, psinkflx, jp_tra, rsfact )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE trc_sink  ***
@@ -45,4 +45,5 @@
       !!---------------------------------------------------------------------
       INTEGER , INTENT(in)  :: kt
+      INTEGER , INTENT(in)  :: Kbb, Kmm
       INTEGER , INTENT(in)  :: jp_tra    ! tracer index index      
       REAL(wp), INTENT(in)  :: rsfact    ! time step duration
@@ -75,5 +76,5 @@
                DO jk = 1, jpkm1
                   IF( tmask(ji,jj,jk) == 1.0 ) THEN
-                      zwsmax =  0.5 * e3t_n(ji,jj,jk) * rday / rsfact
+                      zwsmax =  0.5 * e3t(ji,jj,jk,Kmm) * rday / rsfact
                       iiter(ji,jj) =  MAX( iiter(ji,jj), INT( pwsink(ji,jj,jk) / zwsmax ) )
                   ENDIF
@@ -88,5 +89,5 @@
             DO ji = 1, jpi
                IF( tmask(ji,jj,jk) == 1 ) THEN
-                 zwsmax = 0.5 * e3t_n(ji,jj,jk) * rday / rsfact
+                 zwsmax = 0.5 * e3t(ji,jj,jk,Kmm) * rday / rsfact
                  zwsink(ji,jj,jk) = MIN( pwsink(ji,jj,jk), zwsmax * REAL( iiter(ji,jj), wp ) )
                ENDIF
@@ -101,5 +102,5 @@
       !   Compute the sedimentation term using trc_sink2 for the considered sinking particle
       !   -----------------------------------------------------
-      CALL trc_sink2( zwsink, psinkflx, jp_tra, iiter, rsfact )
+      CALL trc_sink2( Kbb, Kmm, zwsink, psinkflx, jp_tra, iiter, rsfact )
       !
       IF( ln_timing )   CALL timing_stop('trc_sink')
@@ -107,5 +108,5 @@
    END SUBROUTINE trc_sink
 
-   SUBROUTINE trc_sink2( pwsink, psinkflx, jp_tra, kiter, rsfact )
+   SUBROUTINE trc_sink2( Kbb, Kmm, pwsink, psinkflx, jp_tra, kiter, rsfact )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE trc_sink2  ***
@@ -118,4 +119,5 @@
       !!      transport term, i.e.  div(u*tra).
       !!---------------------------------------------------------------------
+      INTEGER,  INTENT(in   )                         ::   Kbb, Kmm  ! time level indices
       INTEGER,  INTENT(in   )                         ::   jp_tra    ! tracer index index      
       REAL(wp), INTENT(in   )                         ::   rsfact    ! duration of time step
@@ -133,5 +135,5 @@
       ztraz(:,:,:) = 0.e0
       zakz (:,:,:) = 0.e0
-      ztrb (:,:,:) = trb(:,:,:,jp_tra)
+      ztrb (:,:,:) = tr(:,:,:,jp_tra,Kbb)
 
       DO jk = 1, jpkm1
@@ -150,5 +152,5 @@
                !              
                DO jk = 2, jpkm1
-                  ztraz(ji,jj,jk) = ( trb(ji,jj,jk-1,jp_tra) - trb(ji,jj,jk,jp_tra) ) * tmask(ji,jj,jk)
+                  ztraz(ji,jj,jk) = ( tr(ji,jj,jk-1,jp_tra,Kbb) - tr(ji,jj,jk,jp_tra,Kbb) ) * tmask(ji,jj,jk)
                END DO
                ztraz(ji,jj,1  ) = 0.0
@@ -169,7 +171,7 @@
                ! vertical advective flux
                DO jk = 1, jpkm1
-                  zigma = zwsink2(ji,jj,jk+1) * zstep / e3w_n(ji,jj,jk+1)
+                  zigma = zwsink2(ji,jj,jk+1) * zstep / e3w(ji,jj,jk+1,Kmm)
                   zew   = zwsink2(ji,jj,jk+1)
-                  psinkflx(ji,jj,jk+1) = -zew * ( trb(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
+                  psinkflx(ji,jj,jk+1) = -zew * ( tr(ji,jj,jk,jp_tra,Kbb) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
                END DO
                !
@@ -179,6 +181,6 @@
          
                DO jk=1,jpkm1
-                  zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
-                  trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + zflx
+                  zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
+                  tr(ji,jj,jk,jp_tra,Kbb) = tr(ji,jj,jk,jp_tra,Kbb) + zflx
                END DO
             END DO
@@ -189,5 +191,5 @@
          DO jj = 1,jpj
             DO ji = 1, jpi
-               zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
+               zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
                ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx
             END DO
@@ -195,5 +197,5 @@
       END DO
 
-      trb(:,:,:,jp_tra) = ztrb(:,:,:)
+      tr(:,:,:,jp_tra,Kbb) = ztrb(:,:,:)
       psinkflx(:,:,:)   = 2. * psinkflx(:,:,:)
       !
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90	(revision 10966)
@@ -61,32 +61,32 @@
       IF( .NOT. lk_c1d ) THEN
          !
-                                CALL trc_sbc    ( kt, Kmm, Krhs )      ! surface boundary condition
-         IF( ln_trabbl )        CALL trc_bbl    ( kt, Kbb, Kmm, Krhs ) ! advective (and/or diffusive) bottom boundary layer scheme
-         IF( ln_trcdmp )        CALL trc_dmp    ( kt, Kmm, Krhs )      ! internal damping trends
-         IF( ln_bdy )           CALL trc_bdy_dmp( kt, Kbb,      Krhs ) ! BDY damping trends
-                                CALL trc_adv    ( kt, Kbb, Kmm, tr, Krhs )      ! horizontal & vertical advection 
-         !                                                             ! Partial top/bottom cell: GRADh( trb )  
+                                CALL trc_sbc    ( kt,      Kmm, Krhs )      ! surface boundary condition
+         IF( ln_trabbl )        CALL trc_bbl    ( kt, Kbb, Kmm, Krhs )      ! advective (and/or diffusive) bottom boundary layer scheme
+         IF( ln_trcdmp )        CALL trc_dmp    ( kt, Kbb, Kmm, Krhs )      ! internal damping trends
+         IF( ln_bdy )           CALL trc_bdy_dmp( kt, Kbb,      Krhs )      ! BDY damping trends
+                                CALL trc_adv    ( kt, Kbb, Kmm, tr, Krhs )  ! horizontal & vertical advection 
+         !                                                             ! Partial top/bottom cell: GRADh( tr(Kbb) )  
          IF( ln_zps ) THEN
-           IF( ln_isfcav ) THEN ; CALL zps_hde_isf( kt, Kmm, jptra, trb, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )  ! both top & bottom
-           ELSE                 ; CALL zps_hde    ( kt, Kmm, jptra, trb, gtru, gtrv )                                      !  only bottom
+           IF( ln_isfcav ) THEN ; CALL zps_hde_isf( kt, Kmm, jptra, tr(:,:,:,:,Kbb), pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )  ! both top & bottom
+           ELSE                 ; CALL zps_hde    ( kt, Kmm, jptra, tr(:,:,:,:,Kbb), gtru, gtrv )                                      !  only bottom
            ENDIF
          ENDIF
          !                                                      
-                                CALL trc_ldf    ( kt, Kmm, Krhs ) ! lateral mixing
+                                CALL trc_ldf    ( kt, Kbb, Kmm, Krhs )     ! lateral mixing
 #if defined key_agrif
          IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_trc       ! tracers sponge
 #endif
-                                CALL trc_zdf( kt, Kbb, Kmm, Krhs, tr, Kaa  )  ! vert. mixing & after tracer	==> after
-                                CALL trc_nxt    ( kt, Kbb, Kmm, Krhs )        ! tracer fields at next time step     
-         IF( ln_trcrad )        CALL trc_rad    ( kt, Kmm, Krhs )             ! Correct artificial negative concentrations
-         IF( ln_trcdmp_clo )    CALL trc_dmp_clo( kt, Kmm )                   ! internal damping trends on closed seas only
+                                CALL trc_zdf    ( kt, Kbb, Kmm, Krhs, tr, Kaa  )  ! vert. mixing & after tracer	==> after
+                                CALL trc_nxt    ( kt, Kbb, Kmm, Krhs )            ! tracer fields at next time step     
+         IF( ln_trcrad )        CALL trc_rad    ( kt, Kbb, Kmm, Krhs )            ! Correct artificial negative concentrations
+         IF( ln_trcdmp_clo )    CALL trc_dmp_clo( kt, Kbb, Kmm )                  ! internal damping trends on closed seas only
 
          !
       ELSE                                               ! 1D vertical configuration
-                                CALL trc_sbc( kt, Kmm, Krhs )            ! surface boundary condition
-         IF( ln_trcdmp )        CALL trc_dmp( kt, Kmm, Krhs )            ! internal damping trends
+                                CALL trc_sbc( kt,      Kmm, Krhs )            ! surface boundary condition
+         IF( ln_trcdmp )        CALL trc_dmp( kt, Kbb, Kmm, Krhs )            ! internal damping trends
                                 CALL trc_zdf( kt, Kbb, Kmm, Krhs, tr, Kaa  )  ! vert. mixing & after tracer	==> after
-                                CALL trc_nxt( kt, Kbb, Kmm, Krhs )       ! tracer fields at next time step     
-          IF( ln_trcrad )       CALL trc_rad( kt, Kmm, Krhs )            ! Correct artificial negative concentrations
+                                CALL trc_nxt( kt, Kbb, Kmm, Krhs )            ! tracer fields at next time step     
+          IF( ln_trcrad )       CALL trc_rad( kt, Kbb, Kmm, Krhs )            ! Correct artificial negative concentrations
          !
       END IF
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trczdf.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trczdf.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trczdf.F90	(revision 10966)
@@ -70,5 +70,5 @@
          WRITE(charout, FMT="('zdf ')")
          CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         CALL prt_ctl_trc( tab4d=tr(:,:,:,:,Kaa), mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
       END IF
       !
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trdmxl_trc.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trdmxl_trc.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trdmxl_trc.F90	(revision 10966)
@@ -16,5 +16,5 @@
    !!   trd_mxl_trc_init : initialization step
    !!----------------------------------------------------------------------
-   USE trc               ! tracer definitions (trn, trb, tra, etc.)
+   USE trc               ! tracer definitions (tr etc.)
    USE trc_oce, ONLY :   nn_dttrc  ! frequency of step on passive tracers
    USE dom_oce           ! domain definition
@@ -70,5 +70,5 @@
 
 
-   SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn )
+   SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trd_mxl_trc_zint  ***
@@ -92,4 +92,5 @@
       !!
       INTEGER, INTENT( in ) ::   ktrd, kjn                        ! ocean trend index and passive tracer rank
+      INTEGER, INTENT( in ) ::   Kmm                              ! time level index
       CHARACTER(len=2), INTENT( in ) ::  ctype                    ! surface/bottom (2D) or interior (3D) physics
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::  ptrc_trdmxl ! passive tracer trend
@@ -150,5 +151,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  IF( jk - nmld_trc(ji,jj) < 0 )   wkx_trc(ji,jj,jk) = e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
+                  IF( jk - nmld_trc(ji,jj) < 0 )   wkx_trc(ji,jj,jk) = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -183,5 +184,5 @@
 
 
-   SUBROUTINE trd_mxl_trc( kt )
+   SUBROUTINE trd_mxl_trc( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trd_mxl_trc  ***
@@ -232,4 +233,5 @@
       !
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kmm                              ! time level index
       !
       INTEGER ::   ji, jj, jk, jl, ik, it, itmod, jn
@@ -267,6 +269,6 @@
                   ik = nmld_trc(ji,jj)
                   IF( ln_trdtrc(jn) )    &
-                  tmltrd_trc(ji,jj,jpmxl_trc_zdf,jn) = - avs(ji,jj,ik) / e3w_n(ji,jj,ik) * tmask(ji,jj,ik)  &
-                       &                    * ( trn(ji,jj,ik-1,jn) - trn(ji,jj,ik,jn) )            &
+                  tmltrd_trc(ji,jj,jpmxl_trc_zdf,jn) = - avs(ji,jj,ik) / e3w(ji,jj,ik,Kmm) * tmask(ji,jj,ik)  &
+                       &                    * ( tr(ji,jj,ik-1,jn,Kmm) - tr(ji,jj,ik,jn,Kmm) )            &
                        &                    / MAX( 1., rmld_trc(ji,jj) ) * tmask(ji,jj,1)
                END DO
@@ -322,5 +324,5 @@
          DO jn = 1, jptra
             IF( ln_trdtrc(jn) ) &
-               tml_trc(:,:,jn) = tml_trc(:,:,jn) + wkx_trc(:,:,jk) * trn(:,:,jk,jn)
+               tml_trc(:,:,jn) = tml_trc(:,:,jn) + wkx_trc(:,:,jk) * tr(:,:,jk,jn,Kmm)
          END DO
       END DO
@@ -970,10 +972,12 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE trd_mxl_trc( kt )                                   ! Empty routine
+   SUBROUTINE trd_mxl_trc( kt, Kmm )                                   ! Empty routine
       INTEGER, INTENT( in) ::   kt
+      INTEGER, INTENT( in) ::   Kmm            ! time level index
       WRITE(*,*) 'trd_mxl_trc: You should not have seen this print! error?', kt
    END SUBROUTINE trd_mxl_trc
-   SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn )
+   SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn, Kmm )
       INTEGER               , INTENT( in ) ::  ktrd, kjn              ! ocean trend index and passive tracer rank
+      INTEGER               , INTENT( in ) ::  Kmm                    ! time level index
       CHARACTER(len=2)      , INTENT( in ) ::  ctype                  ! surface/bottom (2D) or interior (3D) physics
       REAL, DIMENSION(:,:,:), INTENT( in ) ::  ptrc_trdmxl            ! passive trc trend
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trdtrc.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trdtrc.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trdtrc.F90	(revision 10966)
@@ -13,5 +13,5 @@
    !!   trdtrc      : passive tracer trends 
    !!----------------------------------------------------------------------
-   USE trc               ! tracer definitions (trn, trb, tra, etc.)
+   USE trc               ! tracer definitions (tr(:,:,:,:,Kmm), tr(:,:,:,:,Kbb), tr(:,:,:,:,Krhs), etc.)
    USE trd_oce
    USE trdtrc_oce       ! definition of main arrays used for trends computations
@@ -32,9 +32,10 @@
 CONTAINS
 
-   SUBROUTINE trd_trc( ptrtrd, kjn, ktrd, kt )
+   SUBROUTINE trd_trc( ptrtrd, kjn, ktrd, kt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trd_trc  ***
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in )  ::   kt                                  ! time step
+      INTEGER, INTENT( in )  ::   Kmm                                 ! time level index
       INTEGER, INTENT( in )  ::   kjn                                 ! tracer index
       INTEGER, INTENT( in )  ::   ktrd                                ! tracer trend index
@@ -56,21 +57,21 @@
          !
          SELECT CASE ( ktrd )
-         CASE ( jptra_xad     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_xad, '3D', kjn )
-         CASE ( jptra_yad     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_yad, '3D', kjn )
-         CASE ( jptra_zad     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_zad, '3D', kjn )
-         CASE ( jptra_ldf     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_ldf, '3D', kjn )
-         CASE ( jptra_bbl     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_bbl, '3D', kjn )
+         CASE ( jptra_xad     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_xad, '3D', kjn, Kmm )
+         CASE ( jptra_yad     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_yad, '3D', kjn, Kmm )
+         CASE ( jptra_zad     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_zad, '3D', kjn, Kmm )
+         CASE ( jptra_ldf     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_ldf, '3D', kjn, Kmm )
+         CASE ( jptra_bbl     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_bbl, '3D', kjn, Kmm )
          CASE ( jptra_zdf     )
             IF( ln_trcldf_iso ) THEN
-               CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_ldf, '3D', kjn )
+               CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_ldf, '3D', kjn, Kmm )
             ELSE
-               CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_zdf, '3D', kjn )
+               CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_zdf, '3D', kjn, Kmm )
             ENDIF
-         CASE ( jptra_dmp     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_dmp , '3D', kjn )
-         CASE ( jptra_nsr     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_sbc , '2D', kjn )
-         CASE ( jptra_sms     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_sms , '3D', kjn )
-         CASE ( jptra_radb    )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_radb, '3D', kjn )
-         CASE ( jptra_radn    )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_radn, '3D', kjn )
-         CASE ( jptra_atf     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_atf , '3D', kjn )
+         CASE ( jptra_dmp     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_dmp , '3D', kjn, Kmm )
+         CASE ( jptra_nsr     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_sbc , '2D', kjn, Kmm )
+         CASE ( jptra_sms     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_sms , '3D', kjn, Kmm )
+         CASE ( jptra_radb    )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_radb, '3D', kjn, Kmm )
+         CASE ( jptra_radn    )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_radn, '3D', kjn, Kmm )
+         CASE ( jptra_atf     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_atf , '3D', kjn, Kmm )
          END SELECT
          !
@@ -110,6 +111,7 @@
 CONTAINS
 
-   SUBROUTINE trd_trc( ptrtrd, kjn, ktrd, kt )
+   SUBROUTINE trd_trc( ptrtrd, kjn, ktrd, kt, Kmm )
       INTEGER               , INTENT( in )     ::   kt      ! time step
+      INTEGER               , INTENT( in )     ::   Kmm     ! time level index
       INTEGER               , INTENT( in )     ::   kjn     ! tracer index
       INTEGER               , INTENT( in )     ::   ktrd    ! tracer trend index
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcsms.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcsms.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcsms.F90	(revision 10966)
@@ -34,5 +34,5 @@
 CONTAINS
 
-   SUBROUTINE trc_sms( kt, Kmm , Krhs )
+   SUBROUTINE trc_sms( kt, Kbb, Kmm , Krhs )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE trc_sms  ***
@@ -43,5 +43,5 @@
       !! -------------------------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt        ! ocean time-step index      
-      INTEGER, INTENT( in ) ::   Kmm, Krhs ! time level indices
+      INTEGER, INTENT( in ) ::   Kbb, Kmm, Krhs ! time level indices
       !!
       CHARACTER (len=25) :: charout
@@ -50,9 +50,9 @@
       IF( ln_timing )   CALL timing_start('trc_sms')
       !
-      IF( ln_pisces  )   CALL trc_sms_pisces ( kt, Kmm )    ! main program of PISCES 
-      IF( ll_cfc     )   CALL trc_sms_cfc    ( kt )    ! surface fluxes of CFC
-      IF( ln_c14     )   CALL trc_sms_c14    ( kt )    ! surface fluxes of C14
-      IF( ln_age     )   CALL trc_sms_age    ( kt )    ! Age tracer
-      IF( ln_my_trc  )   CALL trc_sms_my_trc ( kt, Kmm, Krhs )    ! MY_TRC  tracers
+      IF( ln_pisces  )   CALL trc_sms_pisces ( kt, Kbb, Kmm )          ! main program of PISCES 
+      IF( ll_cfc     )   CALL trc_sms_cfc    ( kt,      Kmm )          ! surface fluxes of CFC
+      IF( ln_c14     )   CALL trc_sms_c14    ( kt,      Kmm )          ! surface fluxes of C14
+      IF( ln_age     )   CALL trc_sms_age    ( kt,      Kmm )          ! Age tracer
+      IF( ln_my_trc  )   CALL trc_sms_my_trc ( kt,      Kmm, Krhs )    ! MY_TRC  tracers
 
       IF(ln_ctl) THEN      ! print mean trends (used for debugging)
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcstp.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcstp.F90	(revision 10965)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcstp.F90	(revision 10966)
@@ -101,5 +101,5 @@
          IF( lrst_trc )            CALL trc_rst_cal  ( kt, 'WRITE' )   ! calendar
                                    CALL trc_wri      ( kt )                       ! output of passive tracers with iom I/O manager
-                                   CALL trc_sms      ( kt,      Kmm, Krhs )       ! tracers: sinks and sources
+                                   CALL trc_sms      ( kt, Kbb, Kmm, Krhs )       ! tracers: sinks and sources
                                    CALL trc_trp      ( kt, Kbb, Kmm, Krhs, Kaa )  ! transport of passive tracers
          IF( kt == nittrc000 ) THEN
@@ -108,5 +108,5 @@
          ENDIF
          IF( lrst_trc )            CALL trc_rst_wri  ( kt, Kbb, Kmm, Krhs )       ! write tracer restart file
-         IF( lk_trdmxl_trc  )      CALL trd_mxl_trc  ( kt )                       ! trends: Mixed-layer
+         IF( lk_trdmxl_trc  )      CALL trd_mxl_trc  ( kt,      Kmm       )       ! trends: Mixed-layer
          !
          IF( nn_dttrc /= 1   )     CALL trc_sub_reset( kt, Kbb, Kmm, Krhs )       ! resetting physical variables when sub-stepping
