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 10962)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/MY_TRC/trcsms_my_trc.F90	(revision 10963)
@@ -32,5 +32,5 @@
 CONTAINS
 
-   SUBROUTINE trc_sms_my_trc( kt )
+   SUBROUTINE trc_sms_my_trc( kt, Kmm, Krhs )
       !!----------------------------------------------------------------------
       !!                     ***  trc_sms_my_trc  ***
@@ -42,4 +42,5 @@
       !
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kmm, Krhs  ! time level indices
       INTEGER ::   jn   ! dummy loop index
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrmyt
@@ -54,5 +55,5 @@
       IF( l_trdtrc )  ALLOCATE( ztrmyt(jpi,jpj,jpk) )
 
-      CALL trc_bc ( kt )       ! tracers: surface and lateral Boundary Conditions
+      CALL trc_bc ( kt, Kmm, Krhs )       ! tracers: surface and lateral Boundary Conditions
 
       ! add here the call to BGC model
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 10962)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 10963)
@@ -46,5 +46,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_sms( kt )
+   SUBROUTINE p4z_sms( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_sms  ***
@@ -59,4 +59,5 @@
       !
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index      
+      INTEGER, INTENT( in ) ::   Kmm     ! time level index
       !!
       INTEGER ::   ji, jj, jk, jnt, jn, jl
@@ -155,5 +156,5 @@
       IF( ln_sediment ) THEN 
          !
-         CALL sed_model( kt )     !  Main program of Sediment model
+         CALL sed_model( kt, Kmm )     !  Main program of Sediment model
          !
          IF( ln_top_euler ) THEN
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/SED/sedmodel.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/SED/sedmodel.F90	(revision 10962)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/SED/sedmodel.F90	(revision 10963)
@@ -16,5 +16,5 @@
 CONTAINS
 
-   SUBROUTINE sed_model ( kt )
+   SUBROUTINE sed_model ( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE sed_model  ***
@@ -30,4 +30,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt       ! number of iteration
+      INTEGER, INTENT(in) ::   Kmm      ! time level indices
 
 
@@ -35,5 +36,5 @@
 
       IF( kt == nittrc000 ) CALL sed_initrc       ! Initialization of sediment model
-                            CALL sed_stp( kt )  ! Time stepping of Sediment model
+                            CALL sed_stp( kt, Kmm )  ! Time stepping of Sediment model
 
       IF( ln_timing )  CALL timing_stop('sed_model')
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/SED/sedstp.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/SED/sedstp.F90	(revision 10962)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/SED/sedstp.F90	(revision 10963)
@@ -29,5 +29,5 @@
 CONTAINS
 
-   SUBROUTINE sed_stp ( kt )
+   SUBROUTINE sed_stp ( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE sed_stp  ***
@@ -45,4 +45,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt       ! number of iteration
+      INTEGER, INTENT(in) ::   Kmm      ! time level indices
       INTEGER :: ji,jk,js,jn,jw
       !!----------------------------------------------------------------------
@@ -52,5 +53,5 @@
       IF( lrst_sed )            CALL sed_rst_cal  ( kt, 'WRITE' )   ! calenda
 
-      IF(ln_sediment_offline)   CALL trc_dmp_sed  ( kt )
+      IF(ln_sediment_offline)   CALL trc_dmp_sed  ( kt, Kmm )
 
       dtsed  = r2dttrc
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/SED/trcdmp_sed.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/SED/trcdmp_sed.F90	(revision 10962)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/SED/trcdmp_sed.F90	(revision 10963)
@@ -54,5 +54,5 @@
 
 
-   SUBROUTINE trc_dmp_sed( kt )
+   SUBROUTINE trc_dmp_sed( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE trc_dmp_sed  ***
@@ -74,4 +74,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kmm  ! time level index
       !
       INTEGER ::   ji, jj, jk, jn, jl, ikt   ! dummy loop indices
@@ -90,5 +91,5 @@
                !
                jl = n_trc_index(jn) 
-               CALL trc_dta( kt, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000
+               CALL trc_dta( kt, Kmm, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000
                !
                DO jj = 1, jpj
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 10962)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/trcsms_pisces.F90	(revision 10963)
@@ -25,5 +25,5 @@
 CONTAINS
 
-   SUBROUTINE trc_sms_pisces( kt )
+   SUBROUTINE trc_sms_pisces( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE trc_sms_pisces  ***
@@ -35,7 +35,8 @@
       !
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index      
+      INTEGER, INTENT( in ) ::   Kmm     ! time level index
       !!---------------------------------------------------------------------
       !
-      IF( ln_p4z .OR. ln_p5z ) THEN  ;   CALL p4z_sms( kt )   !  PISCES
+      IF( ln_p4z .OR. ln_p5z ) THEN  ;   CALL p4z_sms( kt, Kmm )   !  PISCES
       ELSE                           ;   CALL p2z_sms( kt )   !  LOBSTER
       ENDIF
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 10962)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcdmp.F90	(revision 10963)
@@ -106,5 +106,5 @@
                !
                jl = n_trc_index(jn) 
-               CALL trc_dta( kt, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000
+               CALL trc_dta( kt, Kmm, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000
                !
                SELECT CASE ( nn_zdmp_tr )
@@ -225,5 +225,5 @@
 
 
-   SUBROUTINE trc_dmp_clo( kt )
+   SUBROUTINE trc_dmp_clo( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trc_dmp_clo  ***
@@ -238,4 +238,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kmm     ! time level indices
       !
       INTEGER :: ji , jj, jk, jn, jl, jc                    ! dummy loop indicesa
@@ -355,5 +356,5 @@
             IF( ln_trc_ini(jn) ) THEN      ! update passive tracers arrays with input data read from file
                 jl = n_trc_index(jn)
-                CALL trc_dta( kt, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000
+                CALL trc_dta( kt, Kmm, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000
                 DO jc = 1, npncts
                    DO jk = 1, jpkm1
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 10962)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcnxt.F90	(revision 10963)
@@ -100,5 +100,5 @@
       CALL lbc_lnk( 'trcnxt', tra(:,:,:,:), 'T', 1. )   
 
-      IF( ln_bdy )  CALL trc_bdy( kt )
+      IF( ln_bdy )  CALL trc_bdy( kt, Kbb, Kmm, Krhs )
 
       IF( l_trdtrc )  THEN             ! trends: store now fields before the Asselin filter application
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 10962)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90	(revision 10963)
@@ -64,5 +64,5 @@
          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 )                 ! BDY 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 )  
@@ -79,6 +79,6 @@
                                 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 )                   ! internal damping trends on closed seas only
+         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
 
          !
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/oce_trc.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/oce_trc.F90	(revision 10962)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/oce_trc.F90	(revision 10963)
@@ -34,17 +34,22 @@
 
    !* ocean fields: here now and after fields *
-   USE oce , ONLY :   un      =>    un      !: i-horizontal velocity (m s-1) 
-   USE oce , ONLY :   vn      =>    vn      !: j-horizontal velocity (m s-1)
-   USE oce , ONLY :   wn      =>    wn      !: vertical velocity (m s-1)  
-   USE oce , ONLY :   tsn     =>    tsn     !: 4D array contaning ( tn, sn )
-   USE oce , ONLY :   tsb     =>    tsb     !: 4D array contaning ( tb, sb )
-   USE oce , ONLY :   tsa     =>    tsa     !: 4D array contaning ( ta, sa )
-   USE oce , ONLY :   rhop    =>    rhop    !: potential volumic mass (kg m-3) 
-   USE oce , ONLY :   rhd     =>    rhd     !: in situ density anomalie rhd=(rho-rau0)/rau0 (no units)
-   USE oce , ONLY :   hdivn   =>    hdivn   !: horizontal divergence (1/s)
-   USE oce , ONLY :   sshn    =>    sshn    !: sea surface height at t-point [m]   
-   USE oce , ONLY :   sshb    =>    sshb    !: sea surface height at t-point [m]   
-   USE oce , ONLY :   ssha    =>    ssha    !: sea surface height at t-point [m]   
-   USE oce , ONLY :   rab_n   =>    rab_n   !: local thermal/haline expension ratio at T-points
+   USE oce , ONLY :   tsn     =>    tsn     !: 4D array contaning ( tn, sn ) !TEMPORARY
+   USE oce , ONLY :   tsb     =>    tsb     !: 4D array contaning ( tb, sb ) !TEMPORARY
+   USE oce , ONLY :   tsa     =>    tsa     !: 4D array contaning ( ta, sa ) !TEMPORARY
+   USE oce , ONLY :   sshn    =>    sshn    !: sea surface height at t-point [m]    !TEMPORARY
+   USE oce , ONLY :   sshb    =>    sshb    !: sea surface height at t-point [m]    !TEMPORARY
+   USE oce , ONLY :   ssha    =>    ssha    !: sea surface height at t-point [m]    !TEMPORARY
+   USE oce , ONLY :    un     =>     un     !: 4D array  !TEMPORARY
+   USE oce , ONLY :    vn     =>     vn     !: 4D array  !TEMPORARY
+   USE oce , ONLY :    wn     =>     wn     !: 4D array  !TEMPORARY
+   USE oce , ONLY :   uu     =>    uu     !: i-horizontal velocity (m s-1) 
+   USE oce , ONLY :   vv     =>    vv     !: j-horizontal velocity (m s-1)
+   USE oce , ONLY :   ww     =>    ww     !: vertical velocity (m s-1)  
+   USE oce , ONLY :   ts     =>    ts     !: 4D array contaning ( tn, sn )
+   USE oce , ONLY :   rhop   =>    rhop   !: potential volumic mass (kg m-3) 
+   USE oce , ONLY :   rhd    =>    rhd    !: in situ density anomalie rhd=(rho-rau0)/rau0 (no units)
+   USE oce , ONLY :   hdiv   =>    hdiv   !: horizontal divergence (1/s)
+   USE oce , ONLY :   ssh    =>    ssh    !: sea surface height at t-point [m]   
+   USE oce , ONLY :   rab_n  =>    rab_n  !: local thermal/haline expension ratio at T-points
 
    !* surface fluxes *
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcbc.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcbc.F90	(revision 10962)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcbc.F90	(revision 10963)
@@ -55,5 +55,5 @@
 CONTAINS
 
-   SUBROUTINE trc_bc_ini( ntrc )
+   SUBROUTINE trc_bc_ini( ntrc, Kmm )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE trc_bc_ini  ***
@@ -64,5 +64,6 @@
       !!              - allocates passive tracer BC data structure 
       !!----------------------------------------------------------------------
-      INTEGER,INTENT(in) :: ntrc                           ! number of tracers
+      INTEGER, INTENT(in) :: ntrc                          ! number of tracers
+      INTEGER, INTENT(in) ::   Kmm                         ! time level index
       !
       INTEGER            :: jl, jn , ib, ibd, ii, ij, ik   ! dummy loop indices
@@ -267,5 +268,5 @@
                         ii = idx_bdy(ib)%nbi(ibd,igrd)
                         ij = idx_bdy(ib)%nbj(ibd,igrd)
-                        trcdta_bdy(jn,ib)%trc(ibd,ik) = trn(ii,ij,ik,jn) * tmask(ii,ij,ik)
+                        trcdta_bdy(jn,ib)%trc(ibd,ik) = tr(ii,ij,ik,jn,Kmm) * tmask(ii,ij,ik)
                      END DO
                   END DO
@@ -333,5 +334,5 @@
 
 
-   SUBROUTINE trc_bc(kt, jit)
+   SUBROUTINE trc_bc(kt, Kmm, Krhs, jit)
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE trc_bc  ***
@@ -344,6 +345,7 @@
       USE fldread
       !!      
-      INTEGER, INTENT(in)           ::   kt    ! ocean time-step index
-      INTEGER, INTENT(in), OPTIONAL ::   jit   ! subcycle time-step index (for timesplitting option)
+      INTEGER, INTENT(in)           ::   kt        ! ocean time-step index
+      INTEGER, INTENT(in)           ::   Kmm, Krhs ! time level indices
+      INTEGER, INTENT(in), OPTIONAL ::   jit       ! subcycle time-step index (for timesplitting option)
       !!
       INTEGER  :: ji, jj, jk, jn, jl             ! Loop index
@@ -412,5 +414,5 @@
                   DO jk = 1, nk_rnf(ji,jj)
                      zrnf = (rnf(ji,jj) + rnf_b(ji,jj)) * 0.5_wp * r1_rau0 / h_rnf(ji,jj)
-                     tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn)  + (trn(ji,jj,jk,jn) * zrnf)
+                     tr(ji,jj,jk,jn,Krhs) = tr(ji,jj,jk,jn,Krhs)  + (tr(ji,jj,jk,jn,Kmm) * zrnf)
                   END DO
                END DO
@@ -425,6 +427,6 @@
             DO jj = 2, jpj
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zfact = 1. / ( e3t_n(ji,jj,1) * rn_bc_time )
-                  tra(ji,jj,1,jn) = tra(ji,jj,1,jn) + rf_trsfac(jl) * sf_trcsbc(jl)%fnow(ji,jj,1) * zfact
+                  zfact = 1. / ( e3t(ji,jj,1,Kmm) * rn_bc_time )
+                  tr(ji,jj,1,jn,Krhs) = tr(ji,jj,1,jn,Krhs) + rf_trsfac(jl) * sf_trcsbc(jl)%fnow(ji,jj,1) * zfact
                END DO
             END DO
@@ -438,5 +440,5 @@
                   DO jk = 1, nk_rnf(ji,jj)
                      zfact = rn_rfact / ( e1e2t(ji,jj) * h_rnf(ji,jj) * rn_bc_time ) 
-                     tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + rf_trcfac(jl) * sf_trccbc(jl)%fnow(ji,jj,1) * zfact
+                     tr(ji,jj,jk,jn,Krhs) = tr(ji,jj,jk,jn,Krhs) + rf_trcfac(jl) * sf_trccbc(jl)%fnow(ji,jj,1) * zfact
                   END DO
                END DO
@@ -455,9 +457,11 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE trc_bc_ini( ntrc )        ! Empty routine
-      INTEGER,INTENT(IN) :: ntrc                           ! number of tracers
+   SUBROUTINE trc_bc_ini( ntrc, Kmm )        ! Empty routine
+      INTEGER, INTENT(IN) :: ntrc                           ! number of tracers
+      INTEGER, INTENT(in) :: Kmm                            ! time level index
       WRITE(*,*) 'trc_bc_ini: You should not have seen this print! error?', kt
    END SUBROUTINE trc_bc_ini
-   SUBROUTINE trc_bc( kt )        ! Empty routine
+   SUBROUTINE trc_bc( kt, Kmm, Krhs )        ! Empty routine
+      INTEGER, INTENT(in) :: Kmm, Krhs ! time level indices
       WRITE(*,*) 'trc_bc: You should not have seen this print! error?', kt
    END SUBROUTINE trc_bc
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcbdy.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcbdy.F90	(revision 10962)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcbdy.F90	(revision 10963)
@@ -37,5 +37,5 @@
 CONTAINS
 
-   SUBROUTINE trc_bdy( kt )
+   SUBROUTINE trc_bdy( kt, Kbb, Kmm, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE trc_bdy  ***
@@ -44,5 +44,6 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) :: kt     ! Main time step counter
+      INTEGER, INTENT( in ) :: kt              ! Main time step counter
+      INTEGER, INTENT( in ) :: Kbb, Kmm, Krhs  ! time level indices
       !!
       INTEGER                           :: ib_bdy ,jn ,igrd ! Loop indeces
@@ -63,13 +64,13 @@
             SELECT CASE( TRIM(trcdta_bdy(jn,ib_bdy)%cn_obc) )
             CASE('none'        )   ;   CYCLE
-            CASE('frs'         )   ;   CALL bdy_frs( idx_bdy(ib_bdy),                tra(:,:,:,jn), ztrc*zfac )
-            CASE('specified'   )   ;   CALL bdy_spe( idx_bdy(ib_bdy),                tra(:,:,:,jn), ztrc*zfac )
-            CASE('neumann'     )   ;   CALL bdy_nmn( idx_bdy(ib_bdy), igrd         , tra(:,:,:,jn) )
-            CASE('orlanski'    )   ;   CALL bdy_orl( idx_bdy(ib_bdy), trb(:,:,:,jn), tra(:,:,:,jn), ztrc*zfac, ll_npo=.false. )
-            CASE('orlanski_npo')   ;   CALL bdy_orl( idx_bdy(ib_bdy), trb(:,:,:,jn), tra(:,:,:,jn), ztrc*zfac, ll_npo=.true. )
+            CASE('frs'         )   ;   CALL bdy_frs( idx_bdy(ib_bdy), tr(:,:,:,jn,Krhs), ztrc*zfac )
+            CASE('specified'   )   ;   CALL bdy_spe( idx_bdy(ib_bdy), tr(:,:,:,jn,Krhs), ztrc*zfac )
+            CASE('neumann'     )   ;   CALL bdy_nmn( idx_bdy(ib_bdy), igrd            , tr(:,:,:,jn,Krhs) )
+            CASE('orlanski'    )   ;   CALL bdy_orl( idx_bdy(ib_bdy), tr(:,:,:,jn,Kbb), tr(:,:,:,jn,Krhs), ztrc*zfac, ll_npo=.false. )
+            CASE('orlanski_npo')   ;   CALL bdy_orl( idx_bdy(ib_bdy), tr(:,:,:,jn,Kbb), tr(:,:,:,jn,Krhs), ztrc*zfac, ll_npo=.true. )
             CASE DEFAULT           ;   CALL ctl_stop( 'trc_bdy : unrecognised option for open boundaries for passive tracers' )
             END SELECT
             ! Boundary points should be updated
-            CALL lbc_bdy_lnk( 'trcbdy', tra(:,:,:,jn), 'T', 1., ib_bdy )
+            CALL lbc_bdy_lnk( 'trcbdy', tr(:,:,:,jn,Krhs), 'T', 1., ib_bdy )
             !
          END DO
@@ -81,5 +82,5 @@
 
 
-   SUBROUTINE trc_bdy_dmp( kt )
+   SUBROUTINE trc_bdy_dmp( kt, Kbb, Krhs )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE trc_bdy_dmp  ***
@@ -90,4 +91,5 @@
       !!----------------------------------------------------------------------
       INTEGER,         INTENT(in) ::   kt
+      INTEGER,         INTENT(in) ::   Kbb, Krhs  ! time level indices
       !! 
       INTEGER  ::   jn             ! Tracer index
@@ -110,6 +112,6 @@
                   zwgt = idx_bdy(ib_bdy)%nbd(ib,igrd)
                   DO ik = 1, jpkm1
-                     zta = zwgt * ( trcdta_bdy(jn, ib_bdy)%trc(ib,ik) - trb(ii,ij,ik,jn) ) * tmask(ii,ij,ik)
-                     tra(ii,ij,ik,jn) = tra(ii,ij,ik,jn) + zta
+                     zta = zwgt * ( trcdta_bdy(jn, ib_bdy)%trc(ib,ik) - tr(ii,ij,ik,jn,Kbb) ) * tmask(ii,ij,ik)
+                     tr(ii,ij,ik,jn,Krhs) = tr(ii,ij,ik,jn,Krhs) + zta
                   END DO
                END DO
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcdta.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcdta.F90	(revision 10962)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcdta.F90	(revision 10963)
@@ -154,5 +154,5 @@
 
 
-   SUBROUTINE trc_dta( kt, sf_trcdta, ptrcfac, ptrcdta)
+   SUBROUTINE trc_dta( kt, Kmm, sf_trcdta, ptrcfac, ptrcdta)
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE trc_dta  ***
@@ -167,4 +167,5 @@
       !!----------------------------------------------------------------------
       INTEGER                          , INTENT(in   )   ::   kt         ! ocean time-step
+      INTEGER                          , INTENT(in   )   ::   Kmm        ! time level index
       TYPE(FLD), DIMENSION(1)          , INTENT(inout)   ::   sf_trcdta  ! array of information on the field to read
       REAL(wp)                         , INTENT(in   )   ::   ptrcfac    ! multiplication factor
@@ -194,5 +195,5 @@
                DO ji = 1, jpi
                   DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
-                     zl = gdept_n(ji,jj,jk)
+                     zl = gdept(ji,jj,jk,Kmm)
                      IF(     zl < gdept_1d(1  ) ) THEN         ! above the first level of data
                         ztp(jk) = ptrcdta(ji,jj,1)
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcini.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcini.F90	(revision 10962)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcini.F90	(revision 10963)
@@ -83,9 +83,9 @@
       ENDIF
       !
-      CALL trc_ini_state  !  passive tracers initialisation : from a restart or from clim
+      CALL trc_ini_state( Kbb, Kmm, Kaa )  !  passive tracers initialisation : from a restart or from clim
       IF( nn_dttrc /= 1 ) &
-      CALL trc_sub_ini    ! Initialize variables for substepping passive tracers
-      !
-      CALL trc_ini_inv   ! Inventories
+      CALL trc_sub_ini( Kmm )    ! Initialize variables for substepping passive tracers
+      !
+      CALL trc_ini_inv( Kmm )    ! Inventories
       !
       IF( ln_timing )   CALL timing_stop('trc_init')
@@ -110,10 +110,11 @@
 
 
-   SUBROUTINE trc_ini_inv
+   SUBROUTINE trc_ini_inv( Kmm )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE trc_ini_stat  ***
       !! ** Purpose :      passive tracers inventories at initialsation phase
       !!----------------------------------------------------------------------
-      INTEGER ::  jk, jn    ! dummy loop indices
+      INTEGER, INTENT(in) ::   Kmm    ! time level index
+      INTEGER             ::  jk, jn  ! dummy loop indices
       CHARACTER (len=25) :: charout
       !!----------------------------------------------------------------------
@@ -125,5 +126,5 @@
       !                          ! masked grid volume
       DO jk = 1, jpk
-         cvol(:,:,jk) = e1e2t(:,:) * e3t_n(:,:,jk) * tmask(:,:,jk)
+         cvol(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
       END DO
       !                          ! total volume of the ocean 
@@ -132,5 +133,5 @@
       trai(:) = 0._wp            ! initial content of all tracers
       DO jn = 1, jptra
-         trai(jn) = trai(jn) + glob_sum( 'trcini', trn(:,:,:,jn) * cvol(:,:,:)   )
+         trai(jn) = trai(jn) + glob_sum( 'trcini', tr(:,:,:,jn,Kmm) * cvol(:,:,:)   )
       END DO
 
@@ -151,5 +152,5 @@
          WRITE(charout, FMT="('ini ')")
          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
 9000  FORMAT('      tracer nb : ',i2,'      name :',a10,'      initial content :',e18.10)
@@ -229,5 +230,5 @@
 
 
-   SUBROUTINE trc_ini_state
+   SUBROUTINE trc_ini_state( Kbb, Kmm, Kaa )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE trc_ini_state ***
@@ -238,15 +239,16 @@
       USE trcdta          ! initialisation from files
       !
-      INTEGER :: jn, jl   ! dummy loop indices
-      !!----------------------------------------------------------------------
-      !
-      IF( ln_trcdta )   CALL trc_dta_ini( jptra )      ! set initial tracers values
-      !
-      IF( ln_my_trc )   CALL trc_bc_ini ( jptra )      ! set tracers Boundary Conditions
+      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa   ! time level index
+      INTEGER             :: jn, jl          ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_trcdta )   CALL trc_dta_ini( jptra )           ! set initial tracers values
+      !
+      IF( ln_my_trc )   CALL trc_bc_ini ( jptra, Kmm )      ! set tracers Boundary Conditions
       !
       !
       IF( ln_rsttr ) THEN              ! restart from a file
         !
-        CALL trc_rst_read
+        CALL trc_rst_read( Kbb, Kmm )
         !
       ELSE                             ! Initialisation of tracer from a file that may also be used for damping
@@ -257,5 +259,5 @@
                IF( ln_trc_ini(jn) ) THEN
                   jl = n_trc_index(jn) 
-                  CALL trc_dta( nit000, sf_trcdta(jl), rf_trfac(jl), trn(:,:,:,jn) )
+                  CALL trc_dta( nit000, Kmm, sf_trcdta(jl), rf_trfac(jl), tr(:,:,:,jn,Kmm) )
                   !
                   ! deallocate data structure if data are not used for damping
@@ -271,10 +273,10 @@
         ENDIF
         !
-        trb(:,:,:,:) = trn(:,:,:,:)
+        tr(:,:,:,:,Kbb) = tr(:,:,:,:,Kmm)
         ! 
       ENDIF
       !
-      tra(:,:,:,:) = 0._wp
-      !                                                         ! Partial top/bottom cell: GRADh(trn)
+      tr(:,:,:,:,Kaa) = 0._wp
+      !                                                         ! Partial top/bottom cell: GRADh(tr(Kmm))
    END SUBROUTINE trc_ini_state
 
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcrst.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcrst.F90	(revision 10962)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcrst.F90	(revision 10963)
@@ -94,5 +94,5 @@
    END SUBROUTINE trc_rst_opn
 
-   SUBROUTINE trc_rst_read
+   SUBROUTINE trc_rst_read( Kbb, Kmm )
       !!----------------------------------------------------------------------
       !!                    ***  trc_rst_opn  ***
@@ -100,4 +100,5 @@
       !! ** purpose  :   read passive tracer fields in restart files
       !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   Kbb, Kmm  ! time level indices
       INTEGER  ::  jn     
 
@@ -110,9 +111,9 @@
       ! READ prognostic variables and computes diagnostic variable
       DO jn = 1, jptra
-         CALL iom_get( numrtr, jpdom_autoglo, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) )
-      END DO
-
-      DO jn = 1, jptra
-         CALL iom_get( numrtr, jpdom_autoglo, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) )
+         CALL iom_get( numrtr, jpdom_autoglo, 'TRN'//ctrcnm(jn), tr(:,:,:,jn,Kmm) )
+      END DO
+
+      DO jn = 1, jptra
+         CALL iom_get( numrtr, jpdom_autoglo, 'TRB'//ctrcnm(jn), tr(:,:,:,jn,Kbb) )
       END DO
       !
@@ -121,5 +122,5 @@
    END SUBROUTINE trc_rst_read
 
-   SUBROUTINE trc_rst_wri( kt )
+   SUBROUTINE trc_rst_wri( kt, Kbb, Kmm, Krhs )
       !!----------------------------------------------------------------------
       !!                    ***  trc_rst_wri  ***
@@ -127,5 +128,6 @@
       !! ** purpose  :   write passive tracer fields in restart files
       !!----------------------------------------------------------------------
-      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  :: jn
@@ -136,9 +138,9 @@
       ! -------------------- 
       DO jn = 1, jptra
-         CALL iom_rstput( kt, nitrst, numrtw, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) )
-      END DO
-
-      DO jn = 1, jptra
-         CALL iom_rstput( kt, nitrst, numrtw, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) )
+         CALL iom_rstput( kt, nitrst, numrtw, 'TRN'//ctrcnm(jn), tr(:,:,:,jn,Kmm) )
+      END DO
+
+      DO jn = 1, jptra
+         CALL iom_rstput( kt, nitrst, numrtw, 'TRB'//ctrcnm(jn), tr(:,:,:,jn,Kbb) )
       END DO
       !
@@ -146,5 +148,5 @@
     
       IF( kt == nitrst ) THEN
-          CALL trc_rst_stat            ! statistics
+          CALL trc_rst_stat( Kmm, Krhs )             ! statistics
           CALL iom_close( numrtw )     ! close the restart file (only at last time step)
 #if ! defined key_trdmxl_trc
@@ -297,5 +299,5 @@
 
 
-   SUBROUTINE trc_rst_stat
+   SUBROUTINE trc_rst_stat( Kmm, Krhs )
       !!----------------------------------------------------------------------
       !!                    ***  trc_rst_stat  ***
@@ -303,4 +305,5 @@
       !! ** purpose  :   Compute tracers statistics
       !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   Kmm, Krhs  ! time level indices
       INTEGER  :: jk, jn
       REAL(wp) :: ztraf, zmin, zmax, zmean, zdrift
@@ -315,11 +318,11 @@
       !
       DO jk = 1, jpk
-         zvol(:,:,jk) = e1e2t(:,:) * e3t_a(:,:,jk) * tmask(:,:,jk)
-      END DO
-      !
-      DO jn = 1, jptra
-         ztraf = glob_sum( 'trcrst', trn(:,:,:,jn) * zvol(:,:,:) )
-         zmin  = MINVAL( trn(:,:,:,jn), mask= ((tmask*SPREAD(tmask_i,DIM=3,NCOPIES=jpk).NE.0.)) )
-         zmax  = MAXVAL( trn(:,:,:,jn), mask= ((tmask*SPREAD(tmask_i,DIM=3,NCOPIES=jpk).NE.0.)) )
+         zvol(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Krhs) * tmask(:,:,jk)
+      END DO
+      !
+      DO jn = 1, jptra
+         ztraf = glob_sum( 'trcrst', tr(:,:,:,jn,Kmm) * zvol(:,:,:) )
+         zmin  = MINVAL( tr(:,:,:,jn,Kmm), mask= ((tmask*SPREAD(tmask_i,DIM=3,NCOPIES=jpk).NE.0.)) )
+         zmax  = MAXVAL( tr(:,:,:,jn,Kmm), mask= ((tmask*SPREAD(tmask_i,DIM=3,NCOPIES=jpk).NE.0.)) )
          IF( lk_mpp ) THEN
             CALL mpp_min( 'trcrst', zmin )      ! min over the global domain
@@ -341,8 +344,10 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE trc_rst_read                      ! Empty routines
+   SUBROUTINE trc_rst_read( Kbb, Kmm)                      ! Empty routines
+      INTEGER, INTENT( in ) :: Kbb, Kmm  ! time level indices
    END SUBROUTINE trc_rst_read
-   SUBROUTINE trc_rst_wri( kt )
-      INTEGER, INTENT ( in ) :: kt
+   SUBROUTINE trc_rst_wri( kt, Kbb, Kmm, Krhs )
+      INTEGER, INTENT( in ) :: kt
+      INTEGER, INTENT( in ) :: Kbb, Kmm, Krhs  ! time level indices
       WRITE(*,*) 'trc_rst_wri: You should not have seen this print! error?', kt
    END SUBROUTINE trc_rst_wri   
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 10962)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcsms.F90	(revision 10963)
@@ -34,5 +34,5 @@
 CONTAINS
 
-   SUBROUTINE trc_sms( kt )
+   SUBROUTINE trc_sms( kt, Kmm , Krhs )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE trc_sms  ***
@@ -42,5 +42,6 @@
       !! ** Method  : -  call the main routine of of each defined tracer model
       !! -------------------------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index      
+      INTEGER, INTENT( in ) ::   kt        ! ocean time-step index      
+      INTEGER, INTENT( in ) ::   Kmm, Krhs ! time level indices
       !!
       CHARACTER (len=25) :: charout
@@ -49,14 +50,14 @@
       IF( ln_timing )   CALL timing_start('trc_sms')
       !
-      IF( ln_pisces  )   CALL trc_sms_pisces ( kt )    ! main program of PISCES 
+      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 )    ! MY_TRC  tracers
+      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)
          WRITE(charout, FMT="('sms ')")
          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
       !
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 10962)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcstp.F90	(revision 10963)
@@ -77,5 +77,5 @@
       IF( .NOT.ln_linssh ) THEN                                           ! update ocean volume due to ssh temporal evolution
          DO jk = 1, jpk
-            cvol(:,:,jk) = e1e2t(:,:) * e3t_n(:,:,jk) * tmask(:,:,jk)
+            cvol(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
          END DO
          IF ( ln_ctl .OR. kt == nitrst .OR. ( ln_check_mass .AND. kt == nitend )              &
@@ -87,5 +87,5 @@
       IF( l_trcdm2dc )   CALL trc_mean_qsr( kt )
       !    
-      IF( nn_dttrc /= 1 )   CALL trc_sub_stp( kt, Kmm )  ! averaging physical variables for sub-stepping
+      IF( nn_dttrc /= 1 )   CALL trc_sub_stp( kt, Kbb, Kmm, Krhs )  ! averaging physical variables for sub-stepping
       !    
       IF( MOD( kt , nn_dttrc ) == 0 ) THEN      ! only every nn_dttrc time step
@@ -96,19 +96,19 @@
          ENDIF
          !
-         tra(:,:,:,:) = 0.e0
+         tr(:,:,:,:,Krhs) = 0.e0
          !
                                    CALL trc_rst_opn  ( kt )       ! Open tracer restart file 
          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 )       ! tracers: sinks and sources
-                                   CALL trc_trp      ( kt, Kbb, Kmm, Krhs, Kaa )       ! transport of passive tracers
+                                   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_trp      ( kt, Kbb, Kmm, Krhs, Kaa )  ! transport of passive tracers
          IF( kt == nittrc000 ) THEN
             CALL iom_close( numrtr )       ! close input tracer restart file
             IF(lwm) CALL FLUSH( numont )   ! flush namelist output
          ENDIF
-         IF( lrst_trc )            CALL trc_rst_wri  ( kt )       ! write tracer restart file
-         IF( lk_trdmxl_trc  )      CALL trd_mxl_trc  ( kt )       ! trends: Mixed-layer
-         !
-         IF( nn_dttrc /= 1   )     CALL trc_sub_reset( kt )       ! resetting physical variables when sub-stepping
+         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( nn_dttrc /= 1   )     CALL trc_sub_reset( kt, Kbb, Kmm, Krhs )       ! resetting physical variables when sub-stepping
          !
       ENDIF
@@ -117,5 +117,5 @@
          ztrai = 0._wp                                                   !  content of all tracers
          DO jn = 1, jptra
-            ztrai = ztrai + glob_sum( 'trcstp', trn(:,:,:,jn) * cvol(:,:,:)   )
+            ztrai = ztrai + glob_sum( 'trcstp', tr(:,:,:,jn,Kmm) * cvol(:,:,:)   )
          END DO
          IF( lwm ) WRITE(numstr,9300) kt,  ztrai / areatot
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcsub.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcsub.F90	(revision 10962)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcsub.F90	(revision 10963)
@@ -43,5 +43,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::    un_tm ,   un_temp   !: i-horizontal velocity average     [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::    vn_tm ,   vn_temp   !: j-horizontal velocity average     [m/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::              wn_temp   !: hold current values of avt, un, vn, wn
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::              wn_temp   !: hold current values of avt, uu(Kmm), vv(Kmm), ww
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   tsn_tm ,  tsn_temp   !: t/s average     [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   avs_tm ,  avs_temp   !: vertical diffusivity coeff. at  w-point   [m2/s]
@@ -59,5 +59,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   wndm_tm ,   wndm_temp   !: 10m wind average [m]
    !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  sshb_hold   !:hold sshb from the beginning of each sub-stepping[m]  
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  sshb_hold   !: hold ssh(Kbb) from the beginning of each sub-stepping[m]  
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   sshb_temp, ssha_temp
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  hdivn_temp, rotn_temp
@@ -82,5 +82,5 @@
 CONTAINS
 
-   SUBROUTINE trc_sub_stp( kt, Kmm )
+   SUBROUTINE trc_sub_stp( kt, Kbb, Kmm, Krhs )
       !!-------------------------------------------------------------------
       !!                     ***  ROUTINE trc_stp  ***
@@ -91,6 +91,6 @@
       !!              on TOP steps, calculate averages.
       !!-------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
-      INTEGER, INTENT( in ) ::   Kmm  ! 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   ! dummy loop indices
@@ -105,5 +105,5 @@
            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
            !
-           sshb_hold  (:,:) = sshn  (:,:)
+           sshb_hold  (:,:) = ssh   (:,:,Kmm)
            emp_b_hold (:,:) = emp_b (:,:)
            !
@@ -114,10 +114,10 @@
       IF( MOD( kt , nn_dttrc ) /= 0 ) THEN
          !
-         un_tm   (:,:,:)        = un_tm   (:,:,:)        + un   (:,:,:)        * e3u_n(:,:,:) 
-         vn_tm   (:,:,:)        = vn_tm   (:,:,:)        + vn   (:,:,:)        * e3v_n(:,:,:) 
-         tsn_tm  (:,:,:,jp_tem) = tsn_tm  (:,:,:,jp_tem) + tsn  (:,:,:,jp_tem) * e3t_n(:,:,:)  
-         tsn_tm  (:,:,:,jp_sal) = tsn_tm  (:,:,:,jp_sal) + tsn  (:,:,:,jp_sal) * e3t_n(:,:,:)  
-         rhop_tm (:,:,:)        = rhop_tm (:,:,:)        + rhop (:,:,:)        * e3t_n(:,:,:)  
-         avs_tm  (:,:,:)        = avs_tm  (:,:,:)        + avs  (:,:,:)        * e3w_n(:,:,:)  
+         un_tm   (:,:,:)        = un_tm   (:,:,:)        + uu   (:,:,:,Kmm)        * e3u(:,:,:,Kmm) 
+         vn_tm   (:,:,:)        = vn_tm   (:,:,:)        + vv   (:,:,:,Kmm)        * e3v(:,:,:,Kmm) 
+         tsn_tm  (:,:,:,jp_tem) = tsn_tm  (:,:,:,jp_tem) + ts   (:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm)  
+         tsn_tm  (:,:,:,jp_sal) = tsn_tm  (:,:,:,jp_sal) + ts   (:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm)  
+         rhop_tm (:,:,:)        = rhop_tm (:,:,:)        + rhop (:,:,:)            * e3t(:,:,:,Kmm)  
+         avs_tm  (:,:,:)        = avs_tm  (:,:,:)        + avs  (:,:,:)            * e3w(:,:,:,Kmm)  
          IF( l_ldfslp ) THEN
             uslp_tm (:,:,:)      = uslp_tm (:,:,:)        + uslp (:,:,:)
@@ -137,5 +137,5 @@
          ENDIF 
          !
-         sshn_tm  (:,:)         = sshn_tm  (:,:)         + sshn  (:,:) 
+         sshn_tm  (:,:)         = sshn_tm  (:,:)         + ssh   (:,:,Kmm) 
          rnf_tm   (:,:)         = rnf_tm   (:,:)         + rnf   (:,:) 
          h_rnf_tm (:,:)         = h_rnf_tm (:,:)         + h_rnf (:,:) 
@@ -149,8 +149,8 @@
       ELSE                           !  It is time to substep 
          !   1. set temporary arrays to hold physics/dynamical variables
-         un_temp    (:,:,:)      = un    (:,:,:)
-         vn_temp    (:,:,:)      = vn    (:,:,:)
-         wn_temp    (:,:,:)      = wn    (:,:,:)
-         tsn_temp   (:,:,:,:)    = tsn   (:,:,:,:)
+         un_temp    (:,:,:)      = uu    (:,:,:,Kmm)
+         vn_temp    (:,:,:)      = vv    (:,:,:,Kmm)
+         wn_temp    (:,:,:)      = ww    (:,:,:)
+         tsn_temp   (:,:,:,:)    = ts    (:,:,:,:,Kmm)
          rhop_temp  (:,:,:)      = rhop  (:,:,:)    
          avs_temp   (:,:,:)      = avs   (:,:,:)
@@ -169,7 +169,7 @@
             ENDIF
          ENDIF 
-         sshn_temp  (:,:)        = sshn  (:,:)
-         sshb_temp  (:,:)        = sshb  (:,:)
-         ssha_temp  (:,:)        = ssha  (:,:)
+         sshn_temp  (:,:)        = ssh   (:,:,Kmm)
+         sshb_temp  (:,:)        = ssh   (:,:,Kbb)
+         ssha_temp  (:,:)        = ssh   (:,:,Krhs)
          rnf_temp   (:,:)        = rnf   (:,:)
          h_rnf_temp (:,:)        = h_rnf (:,:)
@@ -182,13 +182,13 @@
          wndm_temp  (:,:)        = wndm  (:,:)
          !                                    !  Variables reset in trc_sub_ssh
-         hdivn_temp (:,:,:)      = hdivn (:,:,:)
+         hdivn_temp (:,:,:)      = hdiv (:,:,:)
          !
          ! 2. Create averages and reassign variables
-         un_tm    (:,:,:)        = un_tm   (:,:,:)        + un   (:,:,:)        * e3u_n(:,:,:) 
-         vn_tm    (:,:,:)        = vn_tm   (:,:,:)        + vn   (:,:,:)        * e3v_n(:,:,:) 
-         tsn_tm   (:,:,:,jp_tem) = tsn_tm  (:,:,:,jp_tem) + tsn  (:,:,:,jp_tem) * e3t_n(:,:,:)  
-         tsn_tm   (:,:,:,jp_sal) = tsn_tm  (:,:,:,jp_sal) + tsn  (:,:,:,jp_sal) * e3t_n(:,:,:)  
-         rhop_tm (:,:,:)         = rhop_tm (:,:,:)        + rhop (:,:,:)        * e3t_n(:,:,:)  
-         avs_tm   (:,:,:)        = avs_tm  (:,:,:)        + avs  (:,:,:)        * e3w_n(:,:,:)  
+         un_tm    (:,:,:)        = un_tm   (:,:,:)        + uu   (:,:,:,Kmm)        * e3u(:,:,:,Kmm) 
+         vn_tm    (:,:,:)        = vn_tm   (:,:,:)        + vv   (:,:,:,Kmm)        * e3v(:,:,:,Kmm) 
+         tsn_tm   (:,:,:,jp_tem) = tsn_tm  (:,:,:,jp_tem) + ts   (:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm)  
+         tsn_tm   (:,:,:,jp_sal) = tsn_tm  (:,:,:,jp_sal) + ts   (:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm)  
+         rhop_tm  (:,:,:)        = rhop_tm (:,:,:)        + rhop (:,:,:)            * e3t(:,:,:,Kmm)  
+         avs_tm   (:,:,:)        = avs_tm  (:,:,:)        + avs  (:,:,:)            * e3w(:,:,:,Kmm)  
          IF( l_ldfslp ) THEN
             uslp_tm  (:,:,:)     = uslp_tm (:,:,:)        + uslp (:,:,:) 
@@ -207,5 +207,5 @@
             ENDIF
          ENDIF 
-         sshn_tm  (:,:)          = sshn_tm    (:,:)       + sshn  (:,:) 
+         sshn_tm  (:,:)          = sshn_tm    (:,:)       + ssh   (:,:,Kmm) 
          rnf_tm   (:,:)          = rnf_tm     (:,:)       + rnf   (:,:) 
          h_rnf_tm (:,:)          = h_rnf_tm   (:,:)       + h_rnf (:,:) 
@@ -217,11 +217,11 @@
          wndm_tm  (:,:)          = wndm_tm    (:,:)       + wndm  (:,:)
          !
-         sshn     (:,:)          = sshn_tm    (:,:) * r1_ndttrcp1 
-         sshb     (:,:)          = sshb_hold  (:,:)
+         ssh      (:,:,Kmm)      = sshn_tm    (:,:) * r1_ndttrcp1 
+         ssh      (:,:,Kbb)      = sshb_hold  (:,:)
          rnf      (:,:)          = rnf_tm     (:,:) * r1_ndttrcp1 
          h_rnf    (:,:)          = h_rnf_tm   (:,:) * r1_ndttrcp1 
          hmld     (:,:)          = hmld_tm    (:,:) * r1_ndttrcp1 
          !  variables that are initialized after averages
-         emp_b    (:,:) = emp_b_hold (:,:)
+         emp_b    (:,:)          = emp_b_hold (:,:)
          IF( kt == nittrc000 ) THEN
             wndm  (:,:)          = wndm_tm    (:,:) * r1_ndttrc 
@@ -261,16 +261,16 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  z1_ne3t = r1_ndttrcp1  / e3t_n(ji,jj,jk)
-                  z1_ne3u = r1_ndttrcp1  / e3u_n(ji,jj,jk)
-                  z1_ne3v = r1_ndttrcp1  / e3v_n(ji,jj,jk)
-                  z1_ne3w = r1_ndttrcp1  / e3w_n(ji,jj,jk)
+                  z1_ne3t = r1_ndttrcp1  / e3t(ji,jj,jk,Kmm)
+                  z1_ne3u = r1_ndttrcp1  / e3u(ji,jj,jk,Kmm)
+                  z1_ne3v = r1_ndttrcp1  / e3v(ji,jj,jk,Kmm)
+                  z1_ne3w = r1_ndttrcp1  / e3w(ji,jj,jk,Kmm)
                   !
-                  un   (ji,jj,jk)        = un_tm   (ji,jj,jk)        * z1_ne3u
-                  vn   (ji,jj,jk)        = vn_tm   (ji,jj,jk)        * z1_ne3v
-                  tsn  (ji,jj,jk,jp_tem) = tsn_tm  (ji,jj,jk,jp_tem) * z1_ne3t
-                  tsn  (ji,jj,jk,jp_sal) = tsn_tm  (ji,jj,jk,jp_sal) * z1_ne3t
-                  rhop (ji,jj,jk)        = rhop_tm (ji,jj,jk)        * z1_ne3t
+                  uu   (ji,jj,jk,Kmm)        = un_tm   (ji,jj,jk)        * z1_ne3u
+                  vv   (ji,jj,jk,Kmm)        = vn_tm   (ji,jj,jk)        * z1_ne3v
+                  ts   (ji,jj,jk,jp_tem,Kmm) = tsn_tm  (ji,jj,jk,jp_tem) * z1_ne3t
+                  ts   (ji,jj,jk,jp_sal,Kmm) = tsn_tm  (ji,jj,jk,jp_sal) * z1_ne3t
+                  rhop (ji,jj,jk)            = rhop_tm (ji,jj,jk)        * z1_ne3t
 !!gm : BUG ==>> for avs I don't understand the division by e3w
-                  avs  (ji,jj,jk)        = avs_tm  (ji,jj,jk)        * z1_ne3w
+                  avs  (ji,jj,jk)            = avs_tm  (ji,jj,jk)        * z1_ne3w
                END DO
             END DO
@@ -283,5 +283,5 @@
          ENDIF
          !
-         CALL trc_sub_ssh( kt, Kmm )         ! after ssh & vertical velocity
+         CALL trc_sub_ssh( kt, Kbb, Kmm, Krhs )         ! after ssh & vertical velocity
          !
       ENDIF
@@ -292,5 +292,5 @@
 
 
-   SUBROUTINE trc_sub_ini
+   SUBROUTINE trc_sub_ini( Kmm )
       !!-------------------------------------------------------------------
       !!                     ***  ROUTINE trc_sub_ini  ***
@@ -302,4 +302,5 @@
       !!-------------------------------------------------------------------
       INTEGER ::   ierr
+      INTEGER, INTENT( in ) ::   Kmm    ! time level index
       !!-------------------------------------------------------------------
       !
@@ -312,11 +313,11 @@
       IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'top_sub_alloc : unable to allocate standard ocean arrays' )
 
-      un_tm   (:,:,:)        = un   (:,:,:)        * e3u_n(:,:,:) 
-      vn_tm   (:,:,:)        = vn   (:,:,:)        * e3v_n(:,:,:) 
-      tsn_tm  (:,:,:,jp_tem) = tsn  (:,:,:,jp_tem) * e3t_n(:,:,:)  
-      tsn_tm  (:,:,:,jp_sal) = tsn  (:,:,:,jp_sal) * e3t_n(:,:,:)  
-      rhop_tm (:,:,:)        = rhop (:,:,:)        * e3t_n(:,:,:)  
+      un_tm   (:,:,:)        = uu   (:,:,:,Kmm)        * e3u(:,:,:,Kmm) 
+      vn_tm   (:,:,:)        = vv   (:,:,:,Kmm)        * e3v(:,:,:,Kmm) 
+      tsn_tm  (:,:,:,jp_tem) = ts   (:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm)  
+      tsn_tm  (:,:,:,jp_sal) = ts   (:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm)  
+      rhop_tm (:,:,:)        = rhop (:,:,:)            * e3t(:,:,:,Kmm)  
 !!gm : BUG? ==>> for avt & avs I don't understand the division by e3w
-      avs_tm  (:,:,:)        = avs  (:,:,:)        * e3w_n(:,:,:)  
+      avs_tm  (:,:,:)        = avs  (:,:,:)            * e3w(:,:,:,Kmm)  
       IF( l_ldfslp ) THEN
          wslpi_tm(:,:,:)     = wslpi(:,:,:)
@@ -325,5 +326,5 @@
          vslp_tm (:,:,:)     = vslp (:,:,:)
       ENDIF
-      sshn_tm  (:,:) = sshn  (:,:) 
+      sshn_tm  (:,:) = ssh   (:,:,Kmm) 
       rnf_tm   (:,:) = rnf   (:,:) 
       h_rnf_tm (:,:) = h_rnf (:,:) 
@@ -333,5 +334,5 @@
       fr_i_tm  (:,:) = 0._wp
       emp_tm   (:,:) = 0._wp
-      fmmflx_tm(:,:)  = 0._wp
+      fmmflx_tm(:,:) = 0._wp
       qsr_tm   (:,:) = 0._wp
       wndm_tm  (:,:) = 0._wp
@@ -350,5 +351,5 @@
 
 
-   SUBROUTINE trc_sub_reset( kt )
+   SUBROUTINE trc_sub_reset( kt, Kbb, Kmm, Krhs )
       !!-------------------------------------------------------------------
       !!                     ***  ROUTINE trc_sub_reset  ***
@@ -359,6 +360,7 @@
       !!              Compute the averages for sub-stepping
       !!-------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::  kt  ! ocean time-step index
-      INTEGER :: jk                 ! dummy loop indices
+      INTEGER, INTENT( in ) ::  kt             ! ocean time-step index
+      INTEGER, INTENT( in ) ::  Kbb, Kmm, Krhs ! time level index
+      INTEGER :: jk                            ! dummy loop indices
       !!-------------------------------------------------------------------
       !
@@ -366,21 +368,21 @@
       !
       !   restore physics variables
-      un    (:,:,:)   =  un_temp    (:,:,:)
-      vn    (:,:,:)   =  vn_temp    (:,:,:)
-      wn    (:,:,:)   =  wn_temp    (:,:,:)
-      tsn   (:,:,:,:) =  tsn_temp   (:,:,:,:)
-      rhop  (:,:,:)   =  rhop_temp  (:,:,:)
-      avs   (:,:,:)   =  avs_temp   (:,:,:)
+      uu   (:,:,:,Kmm)   =  un_temp  (:,:,:)
+      vv   (:,:,:,Kmm)   =  vn_temp  (:,:,:)
+      ww   (:,:,:)       =  wn_temp  (:,:,:)
+      ts   (:,:,:,:,Kmm) =  tsn_temp (:,:,:,:)
+      rhop (:,:,:)       =  rhop_temp(:,:,:)
+      avs  (:,:,:)       =  avs_temp (:,:,:)
       IF( l_ldfslp ) THEN
-         wslpi (:,:,:)=  wslpi_temp (:,:,:)
-         wslpj (:,:,:)=  wslpj_temp (:,:,:)
-         uslp  (:,:,:)=  uslp_temp  (:,:,:)
-         vslp  (:,:,:)=  vslp_temp  (:,:,:)
-      ENDIF
-      sshn  (:,:)     =  sshn_temp  (:,:)
-      sshb  (:,:)     =  sshb_temp  (:,:)
-      ssha  (:,:)     =  ssha_temp  (:,:)
-      rnf   (:,:)     =  rnf_temp   (:,:)
-      h_rnf (:,:)     =  h_rnf_temp (:,:)
+         wslpi (:,:,:)   =  wslpi_temp (:,:,:)
+         wslpj (:,:,:)   =  wslpj_temp (:,:,:)
+         uslp  (:,:,:)   =  uslp_temp  (:,:,:)
+         vslp  (:,:,:)   =  vslp_temp  (:,:,:)
+      ENDIF
+      ssh   (:,:,Kmm)     =  sshn_temp  (:,:)
+      ssh   (:,:,Kbb)     =  sshb_temp  (:,:)
+      ssh   (:,:,Krhs)    =  ssha_temp  (:,:)
+      rnf   (:,:)         =  rnf_temp   (:,:)
+      h_rnf (:,:)         =  h_rnf_temp (:,:)
       !
       hmld  (:,:)     =  hmld_temp  (:,:)
@@ -402,13 +404,13 @@
       ENDIF
       !
-      hdivn (:,:,:)   =  hdivn_temp (:,:,:)
+      hdiv (:,:,:)   =  hdivn_temp (:,:,:)
       !                                      
       ! Start new averages
-         un_tm   (:,:,:)        = un   (:,:,:)        * e3u_n(:,:,:) 
-         vn_tm   (:,:,:)        = vn   (:,:,:)        * e3v_n(:,:,:) 
-         tsn_tm  (:,:,:,jp_tem) = tsn  (:,:,:,jp_tem) * e3t_n(:,:,:)  
-         tsn_tm  (:,:,:,jp_sal) = tsn  (:,:,:,jp_sal) * e3t_n(:,:,:)  
-         rhop_tm (:,:,:)        = rhop (:,:,:)        * e3t_n(:,:,:)  
-         avs_tm  (:,:,:)        = avs  (:,:,:)        * e3w_n(:,:,:)  
+         un_tm   (:,:,:)        = uu   (:,:,:,Kmm)        * e3u(:,:,:,Kmm) 
+         vn_tm   (:,:,:)        = vv   (:,:,:,Kmm)        * e3v(:,:,:,Kmm) 
+         tsn_tm  (:,:,:,jp_tem) = ts   (:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm)  
+         tsn_tm  (:,:,:,jp_sal) = ts   (:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm)  
+         rhop_tm (:,:,:)        = rhop (:,:,:)            * e3t(:,:,:,Kmm)  
+         avs_tm  (:,:,:)        = avs  (:,:,:)            * e3w(:,:,:,Kmm)  
       IF( l_ldfslp ) THEN
          uslp_tm (:,:,:)        = uslp (:,:,:)
@@ -418,7 +420,7 @@
       ENDIF
       !
-      sshb_hold  (:,:) = sshn  (:,:)
+      sshb_hold  (:,:) = ssh   (:,:,Kmm )
       emp_b_hold (:,:) = emp   (:,:)
-      sshn_tm    (:,:) = sshn  (:,:) 
+      sshn_tm    (:,:) = ssh   (:,:,Kmm) 
       rnf_tm     (:,:) = rnf   (:,:) 
       h_rnf_tm   (:,:) = h_rnf (:,:) 
@@ -446,9 +448,9 @@
 
 
-   SUBROUTINE trc_sub_ssh( kt, Kmm ) 
+   SUBROUTINE trc_sub_ssh( kt, Kbb, Kmm, Krhs ) 
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE trc_sub_ssh  ***
       !!                   
-      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity
+      !! ** Purpose :   compute the after ssh (ssh(:,:,Krhs)), the now vertical velocity
       !!              and update the now vertical coordinate (ln_linssh=F).
       !!
@@ -458,12 +460,12 @@
       !!        The boundary conditions are w=0 at the bottom (no flux) and.
       !!
-      !! ** action  :   ssha    : after sea surface height
-      !!                wn      : now vertical velocity
+      !! ** action  :   ssh(:,:,Krhs)    : after sea surface height
+      !!                ww      : now vertical velocity
       !!                sshu_a, sshv_a, sshf_a  : after sea surface height (ln_linssh=F)
       !!
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! time step
-      INTEGER, INTENT(in) ::   Kmm  ! ocean time-level index
+      INTEGER, INTENT(in) ::   kt              ! time step
+      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! ocean time-level index
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -481,9 +483,9 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
          !
-         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
-         !
-      ENDIF
-      !
-!!gm BUG here !   hdivn will include the runoff divergence at the wrong timestep !!!!
+         ww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
+         !
+      ENDIF
+      !
+!!gm BUG here !   hdiv will include the runoff divergence at the wrong timestep !!!!
       CALL div_hor( kt, Kmm )                         ! Horizontal divergence & Relative vorticity
       !
@@ -496,5 +498,5 @@
       zhdiv(:,:) = 0._wp
       DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
-        zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
+        zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
       END DO
       !                                                ! Sea surface elevation time stepping
@@ -502,5 +504,5 @@
       ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
       z1_rau0 = 0.5 / rau0
-      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
+      ssh(:,:,Krhs) = (  ssh(:,:,Kbb) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
 
       IF( .NOT.ln_dynspg_ts ) THEN
@@ -511,6 +513,6 @@
 #endif
          IF( ln_bdy ) THEN
-            ssha(:,:) = ssha(:,:) * bdytmask(:,:)
-            CALL lbc_lnk( 'trcsub', ssha, 'T', 1. ) 
+            ssh(:,:,Krhs) = ssh(:,:,Krhs) * bdytmask(:,:)
+            CALL lbc_lnk( 'trcsub', ssh(:,:,Krhs), 'T', 1. ) 
          ENDIF
       ENDIF
@@ -522,8 +524,8 @@
       DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
          ! - ML - need 3 lines here because replacement of e3t by its expression yields too long lines otherwise
-         wn(:,:,jk) = wn(:,:,jk+1) -   e3t_n(:,:,jk) * hdivn(:,:,jk)        &
-            &                      - ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )    &
+         ww(:,:,jk) = ww(:,:,jk+1) -   e3t(:,:,jk,Kmm ) * hdiv(:,:,jk)          &
+            &                      - ( e3t(:,:,jk,Krhs) - e3t (:,:,jk,Kbb) )    &
             &                         * tmask(:,:,jk) * z1_2dt
-         IF( ln_bdy ) wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
+         IF( ln_bdy ) ww(:,:,jk) = ww(:,:,jk) * bdytmask(:,:)
       END DO
       !
@@ -588,8 +590,10 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE trc_sub_stp( kt )        ! Empty routine
+   SUBROUTINE trc_sub_stp( kt, Kbb, Kmm, Krhs )        ! Empty routine
+      INTEGER, INTENT( in ) ::   Kbb, Kmm, Krhs  ! time level indices
       WRITE(*,*) 'trc_sub_stp: You should not have seen this print! error?', kt
    END SUBROUTINE trc_sub_stp
-   SUBROUTINE trc_sub_ini        ! Empty routine
+   SUBROUTINE trc_sub_ini( Kmm )        ! Empty routine
+      INTEGER, INTENT( in ) ::   Kmm    ! time level index
       WRITE(*,*) 'trc_sub_ini: You should not have seen this print! error?', kt
    END SUBROUTINE trc_sub_ini
