Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ASM/asminc.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ASM/asminc.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ASM/asminc.F90	(revision 10954)
@@ -102,5 +102,5 @@
 CONTAINS
 
-   SUBROUTINE asm_inc_init
+   SUBROUTINE asm_inc_init( Kmm )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE asm_inc_init  ***
@@ -112,4 +112,5 @@
       !! ** Action  : 
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kmm  ! time level index
       INTEGER :: ji, jj, jk, jt  ! dummy loop indices
       INTEGER :: imid, inum      ! local integers
@@ -496,5 +497,5 @@
          IF( ln_bkgwri )   CALL asm_bkg_wri( nit000 - 1 )      ! Output background fields
          IF( ln_asmdin ) THEN                                  ! Direct initialization
-            IF( ln_trainc )   CALL tra_asm_inc( nit000 - 1 )      ! Tracers
+            IF( ln_trainc )   CALL tra_asm_inc( nit000 - 1, Kmm )      ! Tracers
             IF( ln_dyninc )   CALL dyn_asm_inc( nit000 - 1 )      ! Dynamics
             IF( ln_sshinc )   CALL ssh_asm_inc( nit000 - 1 )      ! SSH
@@ -505,5 +506,5 @@
    
    
-   SUBROUTINE tra_asm_inc( kt )
+   SUBROUTINE tra_asm_inc( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE tra_asm_inc  ***
@@ -516,4 +517,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(IN) ::   kt   ! Current time step
+      INTEGER, INTENT(IN) ::   Kmm  ! Current time level index
       !
       INTEGER  :: ji, jj, jk
@@ -607,10 +609,10 @@
 !!gm
 
-            IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)      &
-               &  CALL zps_hde    ( kt, jpts, tsb, gtsu, gtsv,        &  ! Partial steps: before horizontal gradient
-               &                              rhd, gru , grv          )  ! of t, s, rd at the last ocean level
-            IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)      &
-               &  CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi,    &    ! Partial steps for top cell (ISF)
-               &                                  rhd, gru , grv , grui, grvi     )    ! of t, s, rd at the last ocean level
+            IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)           &
+               &  CALL zps_hde    ( kt, Kmm, jpts, tsb, gtsu, gtsv,        &  ! Partial steps: before horizontal gradient
+               &                              rhd, gru , grv               )  ! of t, s, rd at the last ocean level
+            IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)                       &
+               &  CALL zps_hde_isf( nit000, Kmm, jpts, tsb, gtsu, gtsv, gtui, gtvi,    &  ! Partial steps for top cell (ISF)
+               &                                  rhd, gru , grv , grui, grvi          )  ! of t, s, rd at the last ocean level
 
             DEALLOCATE( t_bkginc )
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcisf.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcisf.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcisf.F90	(revision 10954)
@@ -698,5 +698,5 @@
                   zdep        = gdepw(ji,jj,ikt,Kmm)
                   !
-                  CALL eos_rab( zts, zdep, zab )
+                  CALL eos_rab( zts, zdep, zab, Kmm )
                   !
                   !! compute length scale 
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/eosbn2.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/eosbn2.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/eosbn2.F90	(revision 10954)
@@ -564,5 +564,5 @@
 
 
-   SUBROUTINE rab_3d( pts, pab )
+   SUBROUTINE rab_3d( pts, pab, Kmm )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE rab_3d  ***
@@ -574,4 +574,5 @@
       !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
       !!----------------------------------------------------------------------
+      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! pot. temperature & salinity
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab   ! thermal/haline expansion ratio
@@ -592,5 +593,5 @@
                DO ji = 1, jpi
                   !
-                  zh  = gdept_n(ji,jj,jk) * r1_Z0                                ! depth
+                  zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth
                   zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
                   zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
@@ -650,5 +651,5 @@
                   zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
                   zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
-                  zh  = gdept_n(ji,jj,jk)                ! depth in meters at t-point
+                  zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters at t-point
                   ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask
                   !
@@ -677,5 +678,5 @@
 
 
-   SUBROUTINE rab_2d( pts, pdep, pab )
+   SUBROUTINE rab_2d( pts, pdep, pab, Kmm )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE rab_2d  ***
@@ -685,4 +686,5 @@
       !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
       !!----------------------------------------------------------------------
+      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
       REAL(wp), DIMENSION(jpi,jpj,jpts)    , INTENT(in   ) ::   pts    ! pot. temperature & salinity
       REAL(wp), DIMENSION(jpi,jpj)         , INTENT(in   ) ::   pdep   ! depth                  [m]
@@ -791,5 +793,5 @@
 
 
-   SUBROUTINE rab_0d( pts, pdep, pab )
+   SUBROUTINE rab_0d( pts, pdep, pab, Kmm )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE rab_0d  ***
@@ -799,4 +801,5 @@
       !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
       !!----------------------------------------------------------------------
+      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
       REAL(wp), DIMENSION(jpts)    , INTENT(in   ) ::   pts    ! pot. temperature & salinity
       REAL(wp),                      INTENT(in   ) ::   pdep   ! depth                  [m]
@@ -887,5 +890,5 @@
 
 
-   SUBROUTINE bn2( pts, pab, pn2 )
+   SUBROUTINE bn2( pts, pab, pn2, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE bn2  ***
@@ -901,4 +904,5 @@
       !!
       !!----------------------------------------------------------------------
+      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celsius,psu]
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pab   ! thermal/haline expansion coef.  [Celsius-1,psu-1]
@@ -914,6 +918,6 @@
          DO jj = 1, jpj          ! surface and bottom value set to zero one for all in istate.F90
             DO ji = 1, jpi
-               zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
-                  &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) ) 
+               zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
+                  &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 
                   !
                zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 
@@ -922,5 +926,5 @@
                pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     &
                   &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  &
-                  &            / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)
+                  &            / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
             END DO
          END DO
@@ -1091,5 +1095,5 @@
 
 
-   SUBROUTINE eos_pen( pts, pab_pe, ppen )
+   SUBROUTINE eos_pen( pts, pab_pe, ppen, Kmm )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE eos_pen  ***
@@ -1111,4 +1115,5 @@
       !!                    pab_pe(:,:,:,jp_sal) is beta_pe
       !!----------------------------------------------------------------------
+      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts     ! pot. temperature & salinity
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab_pe  ! alpha_pe and beta_pe
@@ -1130,5 +1135,5 @@
                DO ji = 1, jpi
                   !
-                  zh  = gdept_n(ji,jj,jk) * r1_Z0                                ! depth
+                  zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth
                   zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
                   zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
@@ -1194,5 +1199,5 @@
                   zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0)
                   zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0)
-                  zh  = gdept_n(ji,jj,jk)              ! depth in meters  at t-point
+                  zh  = gdept(ji,jj,jk,Kmm)              ! depth in meters  at t-point
                   ztm = tmask(ji,jj,jk)                ! tmask
                   zn  = 0.5_wp * zh * r1_rau0 * ztm
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv.F90	(revision 10954)
@@ -129,5 +129,5 @@
          &              CALL ldf_eiv_trp( kt, nit000, zuu, zvv, zww, 'TRA', Kmm, Krhs )   ! add the eiv transport (if necessary)
       !
-      IF( ln_mle    )   CALL tra_mle_trp( kt, nit000, zuu, zvv, zww, 'TRA' )   ! add the mle transport (if necessary)
+      IF( ln_mle    )   CALL tra_mle_trp( kt, nit000, zuu, zvv, zww, 'TRA', Kmm       )   ! add the mle transport (if necessary)
       !
       CALL iom_put( "uocetr_eff", zuu )                                        ! output effective transport      
@@ -170,6 +170,6 @@
       ENDIF
       !                                              ! print mean trends (used for debugging)
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv  - Ta: ', mask1=tmask,               &
-         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' adv  - Ta: ', mask1=tmask,               &
+         &                       tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
       IF( ln_timing )   CALL timing_stop( 'tra_adv' )
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbc.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbc.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbc.F90	(revision 10954)
@@ -84,22 +84,22 @@
       IF( l_trdtra )   THEN         ! Save the input temperature trend
          ALLOCATE( ztrdt(jpi,jpj,jpk) )
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
+         ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs)
       ENDIF
       !                             !  Add the geothermal trend on temperature
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
-            tsa(ji,jj,mbkt(ji,jj),jp_tem) = tsa(ji,jj,mbkt(ji,jj),jp_tem) + qgh_trd0(ji,jj) / e3t_n(ji,jj,mbkt(ji,jj))
+            ts(ji,jj,mbkt(ji,jj),jp_tem,Krhs) = ts(ji,jj,mbkt(ji,jj),jp_tem,Krhs) + qgh_trd0(ji,jj) / e3t(ji,jj,mbkt(ji,jj),Kmm)
          END DO
       END DO
       !
-      CALL lbc_lnk( 'trabbc', tsa(:,:,:,jp_tem) , 'T', 1. )
+      CALL lbc_lnk( 'trabbc', ts(:,:,:,jp_tem,Krhs) , 'T', 1. )
       !
       IF( l_trdtra ) THEN        ! Send the trend for diagnostics
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
+         ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_bbc, ztrdt )
          DEALLOCATE( ztrdt )
       ENDIF
       !
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbc  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Krhs), clinfo1=' bbc  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
       !
       IF( ln_timing )   CALL timing_stop('tra_bbc')
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbl.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbl.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbl.F90	(revision 10954)
@@ -89,5 +89,5 @@
 
 
-   SUBROUTINE tra_bbl( kt, Kmm, Krhs )
+   SUBROUTINE tra_bbl( kt, Kbb, Kmm, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE bbl  ***
@@ -102,5 +102,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step
-      INTEGER, INTENT( in ) ::   Kmm, Krhs  ! time level indices
+      INTEGER, INTENT( in ) ::   Kbb, Kmm, Krhs  ! time level indices
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
@@ -111,16 +111,16 @@
       IF( l_trdtra )   THEN                         !* Save the T-S input trends
          ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
-      ENDIF
-
-      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl)
+         ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs)
+         ztrds(:,:,:) = ts(:,:,:,jp_sal,Krhs)
+      ENDIF
+
+      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA', Kbb, Kmm )   !* bbl coef. and transport (only if not already done in trcbbl)
 
       IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl
          !
-         CALL tra_bbl_dif( tsb, tsa, jpts )
+         CALL tra_bbl_dif( ts(:,:,:,:,Kbb), ts(:,:,:,:,Krhs), jpts, Kmm )
          IF( ln_ctl )  &
-         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
-            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
+         CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
+            &          tab3d_2=ts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
          ! lateral boundary conditions ; just need for outputs
          CALL lbc_lnk_multi( 'trabbl', ahu_bbl, 'U', 1. , ahv_bbl, 'V', 1. )
@@ -132,8 +132,8 @@
       IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl
          !
-         CALL tra_bbl_adv( tsb, tsa, jpts )
+         CALL tra_bbl_adv( ts(:,:,:,:,Kbb), ts(:,:,:,:,Krhs), jpts, Kmm )
          IF(ln_ctl)   &
-         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
-            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
+         CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
+            &          tab3d_2=ts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
          ! lateral boundary conditions ; just need for outputs
          CALL lbc_lnk_multi( 'trabbl', utr_bbl, 'U', 1. , vtr_bbl, 'V', 1. )
@@ -144,6 +144,6 @@
 
       IF( l_trdtra )   THEN                      ! send the trends for further diagnostics
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
+         ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
+         ztrds(:,:,:) = ts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_bbl, ztrdt )
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_bbl, ztrds )
@@ -156,5 +156,5 @@
 
 
-   SUBROUTINE tra_bbl_dif( ptb, pta, kjpt )
+   SUBROUTINE tra_bbl_dif( ptb, pta, kjpt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_bbl_dif  ***
@@ -180,4 +180,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
+      INTEGER                              , INTENT(in   ) ::   Kmm    ! time level indices
       !
       INTEGER  ::   ji, jj, jn   ! dummy loop indices
@@ -204,5 +205,5 @@
                   &                + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )     &
                   &                - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )  )  &
-                  &             * r1_e1e2t(ji,jj) / e3t_n(ji,jj,ik)
+                  &             * r1_e1e2t(ji,jj) / e3t(ji,jj,ik,Kmm)
             END DO
          END DO
@@ -213,5 +214,5 @@
 
 
-   SUBROUTINE tra_bbl_adv( ptb, pta, kjpt )
+   SUBROUTINE tra_bbl_adv( ptb, pta, kjpt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trc_bbl  ***
@@ -231,4 +232,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
+      INTEGER                              , INTENT(in   ) ::   Kmm    ! time level indices
       !
       INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
@@ -251,15 +253,15 @@
                   !
                   !                                               ! up  -slope T-point (shelf bottom point)
-                  zbtr = r1_e1e2t(iis,jj) / e3t_n(iis,jj,ikus)
+                  zbtr = r1_e1e2t(iis,jj) / e3t(iis,jj,ikus,Kmm)
                   ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr
                   pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra
                   !
                   DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
-                     zbtr = r1_e1e2t(iid,jj) / e3t_n(iid,jj,jk)
+                     zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,jk,Kmm)
                      ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr
                      pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra
                   END DO
                   !
-                  zbtr = r1_e1e2t(iid,jj) / e3t_n(iid,jj,ikud)
+                  zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,ikud,Kmm)
                   ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr
                   pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra
@@ -273,15 +275,15 @@
                   !
                   ! up  -slope T-point (shelf bottom point)
-                  zbtr = r1_e1e2t(ji,ijs) / e3t_n(ji,ijs,ikvs)
+                  zbtr = r1_e1e2t(ji,ijs) / e3t(ji,ijs,ikvs,Kmm)
                   ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
                   pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra
                   !
                   DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
-                     zbtr = r1_e1e2t(ji,ijd) / e3t_n(ji,ijd,jk)
+                     zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,jk,Kmm)
                      ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
                      pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra
                   END DO
                   !                                               ! down-slope T-point (deep bottom point)
-                  zbtr = r1_e1e2t(ji,ijd) / e3t_n(ji,ijd,ikvd)
+                  zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,ikvd,Kmm)
                   ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
                   pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra
@@ -296,5 +298,5 @@
 
 
-   SUBROUTINE bbl( kt, kit000, cdtype )
+   SUBROUTINE bbl( kt, kit000, cdtype, Kbb, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE bbl  ***
@@ -325,4 +327,5 @@
       INTEGER         , INTENT(in   ) ::   kit000   ! first time step index
       CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
+      INTEGER         , INTENT(in   ) ::   Kbb, Kmm ! ocean time level index
       !
       INTEGER  ::   ji, jj                    ! dummy loop indices
@@ -345,14 +348,14 @@
          DO ji = 1, jpi
             ik = mbkt(ji,jj)                             ! bottom T-level index
-            zts (ji,jj,jp_tem) = tsb(ji,jj,ik,jp_tem)    ! bottom before T and S
-            zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal)
+            zts (ji,jj,jp_tem) = ts(ji,jj,ik,jp_tem,Kbb)    ! bottom before T and S
+            zts (ji,jj,jp_sal) = ts(ji,jj,ik,jp_sal,Kbb)
             !
-            zdep(ji,jj) = gdept_n(ji,jj,ik)              ! bottom T-level reference depth
-            zub (ji,jj) = un(ji,jj,mbku(ji,jj))          ! bottom velocity
-            zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj))
+            zdep(ji,jj) = gdept(ji,jj,ik,Kmm)              ! bottom T-level reference depth
+            zub (ji,jj) = uu(ji,jj,mbku(ji,jj),Kmm)          ! bottom velocity
+            zvb (ji,jj) = vv(ji,jj,mbkv(ji,jj),Kmm)
          END DO
       END DO
       !
-      CALL eos_rab( zts, zdep, zab )
+      CALL eos_rab( zts, zdep, zab, Kmm )
       !
       !                                   !-------------------!
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tradmp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tradmp.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tradmp.F90	(revision 10954)
@@ -72,5 +72,5 @@
 
 
-   SUBROUTINE tra_dmp( kt, Kmm, Krhs )
+   SUBROUTINE tra_dmp( kt, Kbb, Kmm, Krhs )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tra_dmp  ***
@@ -91,5 +91,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
       !
       INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
@@ -102,5 +102,5 @@
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
          ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) ) 
-         ztrdts(:,:,:,:) = tsa(:,:,:,:) 
+         ztrdts(:,:,:,:) = ts(:,:,:,:,Krhs) 
       ENDIF
       !                           !==  input T-S data at kt  ==!
@@ -114,5 +114,5 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - tsb(ji,jj,jk,jn) )
+                     ts(ji,jj,jk,jn,Krhs) = ts(ji,jj,jk,jn,Krhs) + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - ts(ji,jj,jk,jn,Kbb) )
                   END DO
                END DO
@@ -125,8 +125,8 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   IF( avt(ji,jj,jk) <= avt_c ) THEN
-                     tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
-                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )
-                     tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
-                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
+                     ts(ji,jj,jk,jp_tem,Krhs) = ts(ji,jj,jk,jp_tem,Krhs)   &
+                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - ts(ji,jj,jk,jp_tem,Kbb) )
+                     ts(ji,jj,jk,jp_sal,Krhs) = ts(ji,jj,jk,jp_sal,Krhs)   &
+                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - ts(ji,jj,jk,jp_sal,Kbb) )
                   ENDIF
                END DO
@@ -138,9 +138,9 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  IF( gdept_n(ji,jj,jk) >= hmlp (ji,jj) ) THEN
-                     tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
-                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )
-                     tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
-                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
+                  IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN
+                     ts(ji,jj,jk,jp_tem,Krhs) = ts(ji,jj,jk,jp_tem,Krhs)   &
+                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - ts(ji,jj,jk,jp_tem,Kbb) )
+                     ts(ji,jj,jk,jp_sal,Krhs) = ts(ji,jj,jk,jp_sal,Krhs)   &
+                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - ts(ji,jj,jk,jp_sal,Kbb) )
                   ENDIF
                END DO
@@ -151,5 +151,5 @@
       !
       IF( l_trdtra )   THEN       ! trend diagnostic
-         ztrdts(:,:,:,:) = tsa(:,:,:,:) - ztrdts(:,:,:,:)
+         ztrdts(:,:,:,:) = ts(:,:,:,:,Krhs) - ztrdts(:,:,:,:)
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) )
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) )
@@ -157,6 +157,6 @@
       ENDIF
       !                           ! Control print
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' dmp  - Ta: ', mask1=tmask,   &
-         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Krhs), clinfo1=' dmp  - Ta: ', mask1=tmask,   &
+         &                       tab3d_2=ts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
       IF( ln_timing )   CALL timing_stop('tra_dmp')
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf.F90	(revision 10954)
@@ -47,5 +47,5 @@
 CONTAINS
 
-   SUBROUTINE tra_ldf( kt, Kmm, Krhs )
+   SUBROUTINE tra_ldf( kt, Kbb, Kmm, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf  ***
@@ -54,5 +54,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
-      INTEGER, INTENT( in ) ::   Kmm, Krhs  ! ocean time level indices
+      INTEGER, INTENT( in ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
       !!
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
@@ -63,22 +63,22 @@
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
          ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
+         ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs) 
+         ztrds(:,:,:) = ts(:,:,:,jp_sal,Krhs)
       ENDIF
       !
       SELECT CASE ( nldf_tra )                 !* compute lateral mixing trend and add it to the general trend
       CASE ( np_lap   )                                  ! laplacian: iso-level operator
-         CALL tra_ldf_lap  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb,      tsa, jpts,  1   )
+         CALL tra_ldf_lap  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,Kbb), ts(:,:,:,:,Krhs),                   jpts,  1, Kmm )
       CASE ( np_lap_i )                                  ! laplacian: standard iso-neutral operator (Madec)
-         CALL tra_ldf_iso  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   )
+         CALL tra_ldf_iso  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,Kbb), ts(:,:,:,:,Kbb), ts(:,:,:,:,Krhs), jpts,  1, Kmm )
       CASE ( np_lap_it )                                 ! laplacian: triad iso-neutral operator (griffies)
-         CALL tra_ldf_triad( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1, Kmm   )
+         CALL tra_ldf_triad( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,Kbb), ts(:,:,:,:,Kbb), ts(:,:,:,:,Krhs), jpts,  1, Kmm )
       CASE ( np_blp , np_blp_i , np_blp_it )             ! bilaplacian: iso-level & iso-neutral operators
-         CALL tra_ldf_blp  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb      , tsa, jpts, nldf_tra, Kmm )
+         CALL tra_ldf_blp  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,Kbb), ts(:,:,:,:,Krhs),             jpts, nldf_tra, Kmm )
       END SELECT
       !
       IF( l_trdtra )   THEN                    !* save the horizontal diffusive trends for further diagnostics
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
+         ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
+         ztrds(:,:,:) = ts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_ldf, ztrdt )
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_ldf, ztrds )
@@ -86,6 +86,6 @@
       ENDIF
       !                                        !* print mean trends (used for debugging)
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf  - Ta: ', mask1=tmask,               &
-         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Krhs), clinfo1=' ldf  - Ta: ', mask1=tmask,               &
+         &                       tab3d_2=ts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
       IF( ln_timing )   CALL timing_stop('tra_ldf')
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_iso.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_iso.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_iso.F90	(revision 10954)
@@ -50,5 +50,5 @@
   SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
       &                                                   pgui, pgvi,   &
-      &                                       ptb , ptbb, pta , kjpt, kpass )
+      &                                       ptb , ptbb, pta , kjpt, kpass, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf_iso  ***
@@ -96,4 +96,5 @@
       INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
       INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
+      INTEGER                              , INTENT(in   ) ::   Kmm        ! ocean time level index
       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
       REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
@@ -182,5 +183,5 @@
                      DO ji = 1, fs_jpim1
                         akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   &
-                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  )
+                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  )
                      END DO
                   END DO
@@ -190,5 +191,5 @@
                   DO jj = 1, jpjm1
                      DO ji = 1, fs_jpim1
-                        ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk)
+                        ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
                         zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
                         akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
@@ -255,6 +256,6 @@
             DO jj = 1 , jpjm1            !==  Horizontal fluxes
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)
-                  zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)
+                  zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm)
+                  zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
                   !
                   zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   &
@@ -280,5 +281,5 @@
                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      &
                      &                                           + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   &
-                     &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                     &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                END DO
             END DO
@@ -325,5 +326,5 @@
                DO jj = 1, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)   &
+                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)   &
                         &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
                         &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
@@ -340,5 +341,5 @@
                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    &
                            &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   &
-                           &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)
+                           &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
                      END DO
                   END DO
@@ -348,5 +349,5 @@
                   DO jj = 1, jpjm1
                      DO ji = fs_2, fs_jpim1
-                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)                      &
+                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)                      &
                            &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   &
                            &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   )
@@ -361,5 +362,5 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   &
-                     &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                     &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                END DO
             END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_lap_blp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_lap_blp.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_lap_blp.F90	(revision 10954)
@@ -47,5 +47,5 @@
    SUBROUTINE tra_ldf_lap( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
       &                                                    pgui, pgvi,   &
-      &                                        ptb , pta , kjpt, kpass ) 
+      &                                        ptb , pta , kjpt, kpass, Kmm ) 
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf_lap  ***
@@ -70,4 +70,5 @@
       INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
       INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
+      INTEGER                              , INTENT(in   ) ::   Kmm        ! ocean time level index
       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
       REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
@@ -100,6 +101,6 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)   !!gm   * umask(ji,jj,jk) pah masked!
-               zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)   !!gm   * vmask(ji,jj,jk)
+               zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm)   !!gm   * umask(ji,jj,jk) pah masked!
+               zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm)   !!gm   * vmask(ji,jj,jk)
             END DO
          END DO
@@ -140,5 +141,5 @@
                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     &
                      &                                   + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )   &
-                     &                                / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) )
+                     &                                / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
                END DO
             END DO
@@ -207,7 +208,7 @@
       !
       CASE ( np_blp    )               ! iso-level bilaplacian
-         CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb,      zlap, kjpt, 1 )
+         CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb,      zlap, kjpt, 1, Kmm )
       CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec)
-         CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 )
+         CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1, Kmm )
       CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies)
          CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1, Kmm )
@@ -216,6 +217,6 @@
       CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1. )     ! Lateral boundary conditions (unchanged sign)
       !                                               ! Partial top/bottom cell: GRADh( zlap )  
-      IF( ln_isfcav .AND. ln_zps ) THEN   ;   CALL zps_hde_isf( kt, kjpt, zlap, zglu, zglv, zgui, zgvi )  ! both top & bottom
-      ELSEIF(             ln_zps ) THEN   ;   CALL zps_hde    ( kt, kjpt, zlap, zglu, zglv )              ! only bottom 
+      IF( ln_isfcav .AND. ln_zps ) THEN   ;   CALL zps_hde_isf( kt, Kmm, kjpt, zlap, zglu, zglv, zgui, zgvi )  ! both top & bottom
+      ELSEIF(             ln_zps ) THEN   ;   CALL zps_hde    ( kt, Kmm, kjpt, zlap, zglu, zglv )              ! only bottom 
       ENDIF
       !
@@ -223,7 +224,7 @@
       !
       CASE ( np_blp    )               ! iso-level bilaplacian
-         CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pta,      kjpt, 2 )
+         CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pta,      kjpt, 2, Kmm )
       CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec)
-         CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 )
+         CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2, Kmm )
       CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies)
          CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2, Kmm )
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_triad.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_triad.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_triad.F90	(revision 10954)
@@ -143,10 +143,10 @@
                   DO jj = 1, jpjm1
                      DO ji = 1, fs_jpim1
-                        ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp)
-                        zbu   = e1e2u(ji,jj) * e3u_n(ji,jj,jk)
+                        ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
+                        zbu   = e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
                         zah   = 0.25_wp * pahu(ji,jj,jk)
                         zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
                         ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth)
-                        zslope2 = zslope_skew + ( gdept_n(ji+1,jj,jk) - gdept_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
+                        zslope2 = zslope_skew + ( gdept(ji+1,jj,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
                         zslope2 = zslope2 *zslope2
                         ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji+ip,jj) * zslope2
@@ -167,11 +167,11 @@
                   DO jj = 1, jpjm1
                      DO ji = 1, fs_jpim1
-                        ze3wr = 1.0_wp / e3w_n(ji,jj+jp,jk+kp)
-                        zbv   = e1e2v(ji,jj) * e3v_n(ji,jj,jk)
+                        ze3wr = 1.0_wp / e3w(ji,jj+jp,jk+kp,Kmm)
+                        zbv   = e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
                         zah   = 0.25_wp * pahv(ji,jj,jk)
                         zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
                         ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces
                         !    (do this by *adding* gradient of depth)
-                        zslope2 = zslope_skew + ( gdept_n(ji,jj+1,jk) - gdept_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
+                        zslope2 = zslope_skew + ( gdept(ji,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
                         zslope2 = zslope2 * zslope2
                         ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj+jp) * zslope2
@@ -194,5 +194,5 @@
                      DO ji = 1, fs_jpim1
                         akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   &
-                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  )
+                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  )
                      END DO
                   END DO
@@ -202,5 +202,5 @@
                   DO jj = 1, jpjm1
                      DO ji = 1, fs_jpim1
-                        ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk)
+                        ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
                         zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
                         akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
@@ -274,10 +274,10 @@
                            ze1ur = r1_e1u(ji,jj)
                            zdxt  = zdit(ji,jj,jk) * ze1ur
-                           ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp)
+                           ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
                            zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr
                            zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
                            zslope_iso  = triadi  (ji+ip,jj,jk,1-ip,kp)
                            !
-                           zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk)
+                           zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
                            ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????   ahu is masked....
                            zah = pahu(ji,jj,jk)
@@ -297,9 +297,9 @@
                            ze2vr = r1_e2v(ji,jj)
                            zdyt  = zdjt(ji,jj,jk) * ze2vr
-                           ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp)
+                           ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,Kmm)
                            zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr
                            zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
                            zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp)
-                           zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk)
+                           zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
                            ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????  ahv is masked...
                            zah = pahv(ji,jj,jk)
@@ -321,10 +321,10 @@
                            ze1ur = r1_e1u(ji,jj)
                            zdxt  = zdit(ji,jj,jk) * ze1ur
-                           ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp)
+                           ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
                            zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr
                            zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
                            zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp)
                            !
-                           zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk)
+                           zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
                            ! ln_botmix_triad is .F. mask zah for bottom half cells
                            zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp)         ! pahu(ji+ip,jj,jk)   ===>>  ????
@@ -344,9 +344,9 @@
                            ze2vr = r1_e2v(ji,jj)
                            zdyt  = zdjt(ji,jj,jk) * ze2vr
-                           ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp)
+                           ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,Kmm)
                            zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr
                            zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
                            zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp)
-                           zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk)
+                           zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
                            ! ln_botmix_triad is .F. mask zah for bottom half cells
                            zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp)         ! pahv(ji,jj+jp,jk)  ????
@@ -365,5 +365,5 @@
                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       &
                      &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   &
-                     &                                        / (  e1e2t(ji,jj) * e3t_n(ji,jj,jk)  )
+                     &                                        / (  e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)  )
                END DO
             END DO
@@ -376,5 +376,5 @@
                DO jj = 1, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)   &
+                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)   &
                         &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
                         &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
@@ -388,5 +388,5 @@
                   DO jj = 1, jpjm1
                      DO ji = fs_2, fs_jpim1
-                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)             &
+                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)             &
                            &                            * ah_wslp2(ji,jj,jk) * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
                      END DO
@@ -397,5 +397,5 @@
                   DO jj = 1, jpjm1
                      DO ji = fs_2, fs_jpim1
-                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)                      &
+                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)                      &
                            &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   &
                            &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   )
@@ -410,5 +410,5 @@
                DO ji = fs_2, fs_jpim1  ! vector opt.
                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
-                     &                                        / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) )
+                     &                                        / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
                END DO
             END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tramle.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tramle.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tramle.F90	(revision 10954)
@@ -56,5 +56,5 @@
 CONTAINS
 
-   SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype )
+   SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_mle_trp  ***
@@ -71,5 +71,5 @@
       !!                p.n = p.n + z._mle
       !!
-      !! ** Action  : - (pun,pvn,pwn) increased by the mle transport
+      !! ** Action  : - (pu,pv,pw) increased by the mle transport
       !!                CAUTION, the transport is not updated at the last line/raw
       !!                         this may be a problem for some advection schemes
@@ -80,4 +80,5 @@
       INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
       INTEGER                         , INTENT(in   ) ::   kit000     ! first time step index
+      INTEGER                         , INTENT(in   ) ::   Kmm        ! ocean time level index
       CHARACTER(len=3)                , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu         ! in : 3 ocean transport components
@@ -115,5 +116,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points
+               zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points
                zmld(ji,jj) = zmld(ji,jj) + zc
                zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0
@@ -147,5 +148,5 @@
       END SELECT
       !                                                ! convert density into buoyancy
-      zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) )
+      zbm(:,:) = + grav * zbm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) )
       !
       !
@@ -205,6 +206,6 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zcuw = 1._wp - ( gdepw_n(ji+1,jj,jk) + gdepw_n(ji,jj,jk) ) * zhu(ji,jj)
-               zcvw = 1._wp - ( gdepw_n(ji,jj+1,jk) + gdepw_n(ji,jj,jk) ) * zhv(ji,jj)
+               zcuw = 1._wp - ( gdepw(ji+1,jj,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhu(ji,jj)
+               zcvw = 1._wp - ( gdepw(ji,jj+1,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhv(ji,jj)
                zcuw = zcuw * zcuw
                zcvw = zcvw * zcvw
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tranpc.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tranpc.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tranpc.F90	(revision 10954)
@@ -85,6 +85,6 @@
          IF( l_trdtra )   THEN                    !* Save initial after fields
             ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
-            ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
-            ztrds(:,:,:) = tsa(:,:,:,jp_sal)
+            ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs) 
+            ztrds(:,:,:) = ts(:,:,:,jp_sal,Krhs)
          ENDIF
          !
@@ -96,6 +96,6 @@
          ENDIF
          !
-         CALL eos_rab( tsa, zab )         ! after alpha and beta (given on T-points)
-         CALL bn2    ( tsa, zab, zn2 )    ! after Brunt-Vaisala  (given on W-points)
+         CALL eos_rab( ts(:,:,:,:,Krhs), zab, Kmm )         ! after alpha and beta (given on T-points)
+         CALL bn2    ( ts(:,:,:,:,Krhs), zab, zn2, Kmm )    ! after Brunt-Vaisala  (given on W-points)
          !
          inpcc = 0
@@ -106,6 +106,6 @@
                IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points
                   !                                     ! consider one ocean column 
-                  zvts(:,jp_tem) = tsa(ji,jj,:,jp_tem)      ! temperature
-                  zvts(:,jp_sal) = tsa(ji,jj,:,jp_sal)      ! salinity
+                  zvts(:,jp_tem) = ts(ji,jj,:,jp_tem,Krhs)      ! temperature
+                  zvts(:,jp_sal) = ts(ji,jj,:,jp_sal,Krhs)      ! salinity
                   !
                   zvab(:,jp_tem)  = zab(ji,jj,:,jp_tem)     ! Alpha 
@@ -187,5 +187,5 @@
                            DO jk = ikup, ikbot      ! Inside the instable (and overlying neutral) portion of the column
                               !
-                              zdz       = e3t_n(ji,jj,jk)
+                              zdz       = e3t(ji,jj,jk,Kmm)
                               zsum_temp = zsum_temp + zvts(jk,jp_tem)*zdz
                               zsum_sali = zsum_sali + zvts(jk,jp_sal)*zdz
@@ -236,6 +236,6 @@
 
                               !! Interpolating alfa and beta at W point:
-                              zrw =  (gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk)) &
-                                 & / (gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk))
+                              zrw =  (gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm)) &
+                                 & / (gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm))
                               zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw
                               zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw
@@ -244,5 +244,5 @@
                               zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     &
                                  &            - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )  &
-                                 &       / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)
+                                 &       / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
 
                               !! OR, faster  => just considering the vertical gradient of density
@@ -288,6 +288,6 @@
 
                   !! Updating tsa:
-                  tsa(ji,jj,:,jp_tem) = zvts(:,jp_tem)
-                  tsa(ji,jj,:,jp_sal) = zvts(:,jp_sal)
+                  ts(ji,jj,:,jp_tem,Krhs) = zvts(:,jp_tem)
+                  ts(ji,jj,:,jp_sal,Krhs) = zvts(:,jp_sal)
 
                   !! LB:  Potentially some other global variable beside theta and S can be treated here
@@ -303,6 +303,6 @@
          IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic
             z1_r2dt = 1._wp / (2._wp * rdt)
-            ztrdt(:,:,:) = ( tsa(:,:,:,jp_tem) - ztrdt(:,:,:) ) * z1_r2dt
-            ztrds(:,:,:) = ( tsa(:,:,:,jp_sal) - ztrds(:,:,:) ) * z1_r2dt
+            ztrdt(:,:,:) = ( ts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) ) * z1_r2dt
+            ztrds(:,:,:) = ( ts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) ) * z1_r2dt
             CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_npc, ztrdt )
             CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_npc, ztrds )
@@ -310,5 +310,5 @@
          ENDIF
          !
-         CALL lbc_lnk_multi( 'tranpc', tsa(:,:,:,jp_tem), 'T', 1., tsa(:,:,:,jp_sal), 'T', 1. )
+         CALL lbc_lnk_multi( 'tranpc', ts(:,:,:,jp_tem,Krhs), 'T', 1., ts(:,:,:,jp_sal,Krhs), 'T', 1. )
          !
          IF( lwp .AND. l_LB_debug ) THEN
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tranxt.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tranxt.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tranxt.F90	(revision 10954)
@@ -64,5 +64,5 @@
 CONTAINS
 
-   SUBROUTINE tra_nxt( kt, Kmm, Krhs )
+   SUBROUTINE tra_nxt( kt, Kbb, Kmm, Krhs )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tranxt  ***
@@ -84,8 +84,8 @@
       !!             domains (lk_agrif=T)
       !!
-      !! ** Action  : - tsb & tsn ready for the next time step
+      !! ** Action  : - ts(Kbb) & ts(Kmm) ready for the next time step
       !!----------------------------------------------------------------------
       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
       !!
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -108,5 +108,5 @@
 #endif
       !                                              ! local domain boundaries  (T-point, unchanged sign)
-      CALL lbc_lnk_multi( 'tranxt', tsa(:,:,:,jp_tem), 'T', 1., tsa(:,:,:,jp_sal), 'T', 1. )
+      CALL lbc_lnk_multi( 'tranxt', ts(:,:,:,jp_tem,Krhs), 'T', 1., ts(:,:,:,jp_sal,Krhs), 'T', 1. )
       !
       IF( ln_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries
@@ -128,8 +128,8 @@
          ! total trend for the non-time-filtered variables. 
          zfact = 1.0 / rdt
-         ! 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
          DO jk = 1, jpkm1
-            ztrdt(:,:,jk) = ( tsa(:,:,jk,jp_tem)*e3t_a(:,:,jk) / e3t_n(:,:,jk) - tsn(:,:,jk,jp_tem)) * zfact
-            ztrds(:,:,jk) = ( tsa(:,:,jk,jp_sal)*e3t_a(:,:,jk) / e3t_n(:,:,jk) - tsn(:,:,jk,jp_sal)) * zfact
+            ztrdt(:,:,jk) = ( ts(:,:,jk,jp_tem,Krhs)*e3t(:,:,jk,Krhs) / e3t(:,:,jk,Kmm) - ts(:,:,jk,jp_tem,Kmm)) * zfact
+            ztrds(:,:,jk) = ( ts(:,:,jk,jp_sal,Krhs)*e3t(:,:,jk,Krhs) / e3t(:,:,jk,Kmm) - ts(:,:,jk,jp_sal,Kmm)) * zfact
          END DO
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_tot, ztrdt )
@@ -138,6 +138,6 @@
             ! Store now fields before applying the Asselin filter 
             ! in order to calculate Asselin filter trend later.
-            ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
-            ztrds(:,:,:) = tsn(:,:,:,jp_sal)
+            ztrdt(:,:,:) = ts(:,:,:,jp_tem,Kmm) 
+            ztrds(:,:,:) = ts(:,:,:,jp_sal,Kmm)
          ENDIF
       ENDIF
@@ -146,5 +146,5 @@
          DO jn = 1, jpts
             DO jk = 1, jpkm1
-               tsn(:,:,jk,jn) = tsa(:,:,jk,jn)    
+               ts(:,:,jk,jn,Kmm) = ts(:,:,jk,jn,Krhs)    
             END DO
          END DO
@@ -159,12 +159,14 @@
       ELSE                                            ! Leap-Frog + Asselin filter time stepping
          !
-         IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt, nit000,      'TRA', tsb, tsn, tsa, jpts )  ! linear free surface 
-         ELSE                   ;   CALL tra_nxt_vvl( kt, Kmm, Krhs, nit000, rdt, 'TRA', tsb, tsn, tsa,   &
-           &                                                                sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface
-         ENDIF
-         !
-         CALL lbc_lnk_multi( 'tranxt', tsb(:,:,:,jp_tem), 'T', 1., tsb(:,:,:,jp_sal), 'T', 1., &
-                  &          tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., &
-                  &          tsa(:,:,:,jp_tem), 'T', 1., tsa(:,:,:,jp_sal), 'T', 1.  )
+         IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt,      Kmm,       nit000,      'TRA',                           &
+           &                                                   ts(:,:,:,:,Kbb), ts(:,:,:,:,Kmm), ts(:,:,:,:,Krhs), jpts )  ! linear free surface 
+         ELSE                   ;   CALL tra_nxt_vvl( kt, Kbb, Kmm, Krhs, nit000, rdt, 'TRA',                           &
+           &                                                   ts(:,:,:,:,Kbb), ts(:,:,:,:,Kmm), ts(:,:,:,:,Krhs),      &
+           &                                                   sbc_tsc        , sbc_tsc_b                        , jpts )  ! non-linear free surface
+         ENDIF
+         !
+         CALL lbc_lnk_multi( 'tranxt', ts(:,:,:,jp_tem,Kbb) , 'T', 1., ts(:,:,:,jp_sal,Kbb) , 'T', 1., &
+                  &                    ts(:,:,:,jp_tem,Kmm) , 'T', 1., ts(:,:,:,jp_sal,Kmm) , 'T', 1., &
+                  &                    ts(:,:,:,jp_tem,Krhs), 'T', 1., ts(:,:,:,jp_sal,Krhs), 'T', 1.  )
          !
       ENDIF     
@@ -173,6 +175,6 @@
          zfact = 1._wp / r2dt             
          DO jk = 1, jpkm1
-            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
-            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
+            ztrdt(:,:,jk) = ( ts(:,:,jk,jp_tem,Kbb) - ztrdt(:,:,jk) ) * zfact
+            ztrds(:,:,jk) = ( ts(:,:,jk,jp_sal,Kbb) - ztrds(:,:,jk) ) * zfact
          END DO
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_atf, ztrdt )
@@ -182,6 +184,6 @@
       !
       !                        ! control print
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
-         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
+         &                       tab3d_2=ts(:,:,:,jp_sal,Kmm), clinfo2=       ' Sn: ', mask2=tmask )
       !
       IF( ln_timing )   CALL timing_stop('tra_nxt')
@@ -190,5 +192,5 @@
 
 
-   SUBROUTINE tra_nxt_fix( kt, kit000, cdtype, ptb, ptn, pta, kjpt )
+   SUBROUTINE tra_nxt_fix( kt, Kmm, kit000, cdtype, ptb, ptn, pta, kjpt )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tra_nxt_fix  ***
@@ -200,7 +202,8 @@
       !!              - swap tracer fields to prepare the next time_step.
       !!
-      !! ** Action  : - tsb & tsn ready for the next time step
+      !! ** Action  : - ptb & ptn ready for the next time step
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   ) ::  kt        ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::  Kmm       ! time level index
       INTEGER                              , INTENT(in   ) ::  kit000    ! first time step index
       CHARACTER(len=3)                     , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator)
@@ -239,5 +242,5 @@
 
 
-   SUBROUTINE tra_nxt_vvl( kt, Kmm, Krhs, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )
+   SUBROUTINE tra_nxt_vvl( kt, Kbb, Kmm, Krhs, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tra_nxt_vvl  ***
@@ -248,12 +251,12 @@
       !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
       !!              - swap tracer fields to prepare the next time_step.
-      !!             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    ] )
+      !!             tb  = ( e3t(Kmm)*tn + atfp*[ e3t(Kbb)*tb - 2 e3t(Kmm)*tn + e3t_a*ta ] )
+      !!                  /( e3t(Kmm)    + atfp*[ e3t(Kbb)    - 2 e3t(Kmm)    + e3t(Krhs)    ] )
       !!             tn  = ta 
       !!
-      !! ** Action  : - tsb & tsn ready for the next time step
+      !! ** Action  : - ptb & ptn ready for the next time step
       !!----------------------------------------------------------------------
       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
       INTEGER                              , INTENT(in   ) ::  kit000    ! first time step index
       REAL(wp)                             , INTENT(in   ) ::  p2dt      ! time-step
@@ -300,7 +303,7 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-                  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  = ptb(ji,jj,jk,jn) * ze3t_b
@@ -323,5 +326,5 @@
                      IF( mikt(ji,jj) <=jk .and. jk <= nk_rnf(ji,jj)  ) THEN
                         ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj) - rnf(ji,jj)   )  ) &
-                    &                            * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) 
+                    &                            * ( e3t(ji,jj,jk,Kmm) / h_rnf(ji,jj) ) 
                      ENDIF
                   ELSE
@@ -339,5 +342,5 @@
                   IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
                      &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
-                     &                              * e3t_n(ji,jj,jk) / h_rnf(ji,jj)
+                     &                              * e3t(ji,jj,jk,Kmm) / h_rnf(ji,jj)
                      !
                   ! ice shelf
@@ -346,9 +349,9 @@
                      IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          &
                         ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
-                               &                 * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj)
+                               &                 * e3t(ji,jj,jk,Kmm) * r1_hisf_tbl (ji,jj)
                      ! level partially include in Losch_2008 ice shelf boundary layer 
                      IF ( jk == misfkb(ji,jj) )                                                   &
                         ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
-                               &                 * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)
+                               &                 * e3t(ji,jj,jk,Kmm) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)
                   END IF
                   !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traqsr.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traqsr.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traqsr.F90	(revision 10954)
@@ -127,5 +127,5 @@
       IF( l_trdtra ) THEN      ! trends diagnostic: save the input temperature trend
          ALLOCATE( ztrdt(jpi,jpj,jpk) ) 
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
+         ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs)
       ENDIF
       !
@@ -173,5 +173,5 @@
                      zze     = 568.2 * zCtot**(-0.746)
                      IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293)
-                     zpsi    = gdepw_n(ji,jj,jk) / zze
+                     zpsi    = gdepw(ji,jj,jk,Kmm) / zze
                      !
                      zlogc   = LOG( zchl )
@@ -219,8 +219,8 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-                  zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * xsi0r       )
-                  zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekb(ji,jj) )
-                  zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekg(ji,jj) )
-                  zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekr(ji,jj) )
+                  zc0 = ze0(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * xsi0r       )
+                  zc1 = ze1(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekb(ji,jj) )
+                  zc2 = ze2(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekg(ji,jj) )
+                  zc3 = ze3(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekr(ji,jj) )
                   ze0(ji,jj,jk) = zc0
                   ze1(ji,jj,jk) = zc1
@@ -249,6 +249,6 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-                  zc0 = zz0 * EXP( -gdepw_n(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk  )*xsi1r )
-                  zc1 = zz0 * EXP( -gdepw_n(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk+1)*xsi1r )
+                  zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r )
+                  zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r )
                   qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 
                END DO
@@ -262,6 +262,6 @@
          DO jj = 2, jpjm1        !-----------------------------!
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
-                  &                 + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t_n(ji,jj,jk)
+               ts(ji,jj,jk,jp_tem,Krhs) = ts(ji,jj,jk,jp_tem,Krhs)   &
+                  &                 + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t(ji,jj,jk,Kmm)
             END DO
          END DO
@@ -296,10 +296,10 @@
       !
       IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
+         ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt )
          DEALLOCATE( ztrdt ) 
       ENDIF
       !                       ! print mean trends (used for debugging)
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Krhs), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
       !
       IF( ln_timing )   CALL timing_stop('tra_qsr')
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trasbc.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trasbc.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trasbc.F90	(revision 10954)
@@ -63,5 +63,5 @@
       !!      (2) Fwe , tracer carried with the water that is exchanged with air+ice. 
       !!               The input forcing fields (emp, rnf, sfx, isf) contain Fext+Fwe,
-      !!             they are simply added to the tracer trend (tsa).
+      !!             they are simply added to the tracer trend (ts(Krhs)).
       !!               In linear free surface case (ln_linssh=T), the volume of the
       !!             ocean does not change with the water exchanges at the (air+ice)-sea
@@ -69,5 +69,5 @@
       !!             concentration/dilution effect associated with water exchanges.
       !!
-      !! ** Action  : - Update tsa with the surface boundary condition trend 
+      !! ** Action  : - Update ts(Krhs) with the surface boundary condition trend 
       !!              - send trends to trdtra module for further diagnostics(l_trdtra=T)
       !!----------------------------------------------------------------------
@@ -91,6 +91,6 @@
       IF( l_trdtra ) THEN                    !* Save ta and sa trends
          ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
+         ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs)
+         ztrds(:,:,:) = ts(:,:,:,jp_sal,Krhs)
       ENDIF
       !
@@ -132,10 +132,10 @@
          DO jj = 2, jpj                         !==>> add concentration/dilution effect due to constant volume cell
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_tem)
-               sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_sal)
+               sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * ts(ji,jj,1,jp_tem,Kmm)
+               sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * ts(ji,jj,1,jp_sal,Kmm)
             END DO
          END DO                                 !==>> output c./d. term
-         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) )
-         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) )
+         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * ts(:,:,1,jp_tem,Kmm) )
+         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * ts(:,:,1,jp_sal,Kmm) )
       ENDIF
       !
@@ -143,5 +143,5 @@
          DO jj = 2, jpj
             DO ji = fs_2, fs_jpim1   ! vector opt.  
-               tsa(ji,jj,1,jn) = tsa(ji,jj,1,jn) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t_n(ji,jj,1)
+               ts(ji,jj,1,jn,Krhs) = ts(ji,jj,1,jn,Krhs) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t(ji,jj,1,Kmm)
             END DO
          END DO
@@ -174,5 +174,5 @@
                DO jk = ikt, ikb - 1
                ! compute trend
-                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                                &
+                  ts(ji,jj,jk,jp_tem,Krhs) = ts(ji,jj,jk,jp_tem,Krhs)                                                &
                      &           + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             &
                      &           * r1_hisf_tbl(ji,jj)
@@ -181,5 +181,5 @@
                ! level partially include in ice shelf boundary layer 
                ! compute trend
-               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                                 &
+               ts(ji,jj,ikb,jp_tem,Krhs) = ts(ji,jj,ikb,jp_tem,Krhs)                                                 &
                   &              + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             &
                   &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)
@@ -200,7 +200,7 @@
                   zdep = zfact / h_rnf(ji,jj)
                   DO jk = 1, nk_rnf(ji,jj)
-                                        tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                 &
+                                        ts(ji,jj,jk,jp_tem,Krhs) = ts(ji,jj,jk,jp_tem,Krhs)                                 &
                                            &                 +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep
-                     IF( ln_rnf_sal )   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                 &
+                     IF( ln_rnf_sal )   ts(ji,jj,jk,jp_sal,Krhs) = ts(ji,jj,jk,jp_sal,Krhs)                                 &
                                            &                 +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 
                   END DO
@@ -210,6 +210,6 @@
       ENDIF
 
-      IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*tsn(:,:,1,jp_tem) )   ! runoff term on sst
-      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*tsn(:,:,1,jp_sal) )   ! runoff term on sss
+      IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*ts(:,:,1,jp_tem,Kmm) )   ! runoff term on sst
+      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*ts(:,:,1,jp_sal,Kmm) )   ! runoff term on sss
 
 #if defined key_asminc
@@ -224,7 +224,7 @@
             DO jj = 2, jpj 
                DO ji = fs_2, fs_jpim1
-                  ztim = ssh_iau(ji,jj) / e3t_n(ji,jj,1)
-                  tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + tsn(ji,jj,1,jp_tem) * ztim
-                  tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + tsn(ji,jj,1,jp_sal) * ztim
+                  ztim = ssh_iau(ji,jj) / e3t(ji,jj,1,Kmm)
+                  ts(ji,jj,1,jp_tem,Krhs) = ts(ji,jj,1,jp_tem,Krhs) + ts(ji,jj,1,jp_tem,Kmm) * ztim
+                  ts(ji,jj,1,jp_sal,Krhs) = ts(ji,jj,1,jp_sal,Krhs) + ts(ji,jj,1,jp_sal,Kmm) * ztim
                END DO
             END DO
@@ -233,6 +233,6 @@
                DO ji = fs_2, fs_jpim1
                   ztim = ssh_iau(ji,jj) / ( ht_n(ji,jj) + 1. - ssmask(ji, jj) )
-                  tsa(ji,jj,:,jp_tem) = tsa(ji,jj,:,jp_tem) + tsn(ji,jj,:,jp_tem) * ztim
-                  tsa(ji,jj,:,jp_sal) = tsa(ji,jj,:,jp_sal) + tsn(ji,jj,:,jp_sal) * ztim
+                  ts(ji,jj,:,jp_tem,Krhs) = ts(ji,jj,:,jp_tem,Krhs) + ts(ji,jj,:,jp_tem,Kmm) * ztim
+                  ts(ji,jj,:,jp_sal,Krhs) = ts(ji,jj,:,jp_sal,Krhs) + ts(ji,jj,:,jp_sal,Kmm) * ztim
                END DO  
             END DO  
@@ -251,7 +251,7 @@
             DO jj = 2, jpj 
                DO ji = fs_2, fs_jpim1
-                  zdep = 1._wp / e3t_n(ji,jj,jk) 
-                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem) * zdep
-                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal) * zdep  
+                  zdep = 1._wp / e3t(ji,jj,jk,Kmm) 
+                  ts(ji,jj,jk,jp_tem,Krhs) = ts(ji,jj,jk,jp_tem,Krhs) - htsc_iscpl(ji,jj,jk,jp_tem) * zdep
+                  ts(ji,jj,jk,jp_sal,Krhs) = ts(ji,jj,jk,jp_sal,Krhs) - htsc_iscpl(ji,jj,jk,jp_sal) * zdep  
                END DO  
             END DO  
@@ -260,6 +260,6 @@
 
       IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
+         ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
+         ztrds(:,:,:) = ts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_nsr, ztrdt )
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_nsr, ztrds )
@@ -267,6 +267,6 @@
       ENDIF
       !
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' sbc  - Ta: ', mask1=tmask,   &
-         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Krhs), clinfo1=' sbc  - Ta: ', mask1=tmask,   &
+         &                       tab3d_2=ts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
       IF( ln_timing )   CALL timing_stop('tra_sbc')
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trazdf.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trazdf.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trazdf.F90	(revision 10954)
@@ -101,6 +101,6 @@
       ENDIF
       !                                          ! print mean trends (used for debugging)
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf  - Ta: ', mask1=tmask,               &
-         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kaa), clinfo1=' zdf  - Ta: ', mask1=tmask,               &
+         &                       tab3d_2=pts(:,:,:,jp_sal,Kaa), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
       IF( ln_timing )   CALL timing_stop('tra_zdf')
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/zpshde.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/zpshde.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/zpshde.F90	(revision 10954)
@@ -39,5 +39,5 @@
 CONTAINS
 
-   SUBROUTINE zps_hde( kt, kjpt, pta, pgtu, pgtv,   &
+   SUBROUTINE zps_hde( kt, Kmm, kjpt, pta, pgtu, pgtv,   &
       &                          prd, pgru, pgrv    )
       !!----------------------------------------------------------------------
@@ -85,4 +85,5 @@
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index
+      INTEGER                              , INTENT(in   )           ::  Kmm         ! ocean time level index
       INTEGER                              , INTENT(in   )           ::  kjpt        ! number of tracers
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta         ! 4D tracers fields
@@ -109,11 +110,11 @@
                iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
                ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
-!!gm BUG ? when applied to before fields, e3w_b should be used....
-               ze3wu = e3w_n(ji+1,jj  ,iku) - e3w_n(ji,jj,iku)
-               ze3wv = e3w_n(ji  ,jj+1,ikv) - e3w_n(ji,jj,ikv)
+!!gm BUG ? when applied to before fields, e3w(:,:,:,Kbb) should be used....
+               ze3wu = e3w(ji+1,jj  ,iku,Kmm) - e3w(ji,jj,iku,Kmm)
+               ze3wv = e3w(ji  ,jj+1,ikv,Kmm) - e3w(ji,jj,ikv,Kmm)
                !
                ! i- direction
                IF( ze3wu >= 0._wp ) THEN      ! case 1
-                  zmaxu =  ze3wu / e3w_n(ji+1,jj,iku)
+                  zmaxu =  ze3wu / e3w(ji+1,jj,iku,Kmm)
                   ! interpolated values of tracers
                   zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
@@ -121,5 +122,5 @@
                   pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
                ELSE                           ! case 2
-                  zmaxu = -ze3wu / e3w_n(ji,jj,iku)
+                  zmaxu = -ze3wu / e3w(ji,jj,iku,Kmm)
                   ! interpolated values of tracers
                   zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
@@ -130,5 +131,5 @@
                ! j- direction
                IF( ze3wv >= 0._wp ) THEN      ! case 1
-                  zmaxv =  ze3wv / e3w_n(ji,jj+1,ikv)
+                  zmaxv =  ze3wv / e3w(ji,jj+1,ikv,Kmm)
                   ! interpolated values of tracers
                   ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
@@ -136,5 +137,5 @@
                   pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
                ELSE                           ! case 2
-                  zmaxv =  -ze3wv / e3w_n(ji,jj,ikv)
+                  zmaxv =  -ze3wv / e3w(ji,jj,ikv,Kmm)
                   ! interpolated values of tracers
                   ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
@@ -155,11 +156,11 @@
                iku = mbku(ji,jj)
                ikv = mbkv(ji,jj)
-               ze3wu  = e3w_n(ji+1,jj  ,iku) - e3w_n(ji,jj,iku)
-               ze3wv  = e3w_n(ji  ,jj+1,ikv) - e3w_n(ji,jj,ikv)
-               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_n(ji  ,jj,iku)     ! i-direction: case 1
-               ELSE                        ;   zhi(ji,jj) = gdept_n(ji+1,jj,iku)     ! -     -      case 2
-               ENDIF
-               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_n(ji,jj  ,ikv)     ! j-direction: case 1
-               ELSE                        ;   zhj(ji,jj) = gdept_n(ji,jj+1,ikv)     ! -     -      case 2
+               ze3wu  = e3w(ji+1,jj  ,iku,Kmm) - e3w(ji,jj,iku,Kmm)
+               ze3wv  = e3w(ji  ,jj+1,ikv,Kmm) - e3w(ji,jj,ikv,Kmm)
+               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept(ji  ,jj,iku,Kmm)     ! i-direction: case 1
+               ELSE                        ;   zhi(ji,jj) = gdept(ji+1,jj,iku,Kmm)     ! -     -      case 2
+               ENDIF
+               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept(ji,jj  ,ikv,Kmm)     ! j-direction: case 1
+               ELSE                        ;   zhj(ji,jj) = gdept(ji,jj+1,ikv,Kmm)     ! -     -      case 2
                ENDIF
             END DO
@@ -173,6 +174,6 @@
                iku = mbku(ji,jj)
                ikv = mbkv(ji,jj)
-               ze3wu  = e3w_n(ji+1,jj  ,iku) - e3w_n(ji,jj,iku)
-               ze3wv  = e3w_n(ji  ,jj+1,ikv) - e3w_n(ji,jj,ikv)
+               ze3wu  = e3w(ji+1,jj  ,iku,Kmm) - e3w(ji,jj,iku,Kmm)
+               ze3wv  = e3w(ji  ,jj+1,ikv,Kmm) - e3w(ji,jj,ikv,Kmm)
                IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1
                ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2
@@ -192,5 +193,5 @@
 
 
-   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv, pgtui, pgtvi,  &
+   SUBROUTINE zps_hde_isf( kt, Kmm, kjpt, pta, pgtu, pgtv, pgtui, pgtvi,  &
       &                          prd, pgru, pgrv, pgrui, pgrvi )
       !!----------------------------------------------------------------------
@@ -241,4 +242,5 @@
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   )           ::  kt           ! ocean time-step index
+      INTEGER                              , INTENT(in   )           ::  Kmm          ! ocean time level index
       INTEGER                              , INTENT(in   )           ::  kjpt         ! number of tracers
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta          ! 4D tracers fields
@@ -270,10 +272,10 @@
                iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
                ikv = mbkv(ji,jj); ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
-               ze3wu = gdept_n(ji+1,jj,iku) - gdept_n(ji,jj,iku)
-               ze3wv = gdept_n(ji,jj+1,ikv) - gdept_n(ji,jj,ikv)
+               ze3wu = gdept(ji+1,jj,iku,Kmm) - gdept(ji,jj,iku,Kmm)
+               ze3wv = gdept(ji,jj+1,ikv,Kmm) - gdept(ji,jj,ikv,Kmm)
                !
                ! i- direction
                IF( ze3wu >= 0._wp ) THEN      ! case 1
-                  zmaxu =  ze3wu / e3w_n(ji+1,jj,iku)
+                  zmaxu =  ze3wu / e3w(ji+1,jj,iku,Kmm)
                   ! interpolated values of tracers
                   zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
@@ -281,5 +283,5 @@
                   pgtu(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
                ELSE                           ! case 2
-                  zmaxu = -ze3wu / e3w_n(ji,jj,iku)
+                  zmaxu = -ze3wu / e3w(ji,jj,iku,Kmm)
                   ! interpolated values of tracers
                   zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
@@ -290,5 +292,5 @@
                ! j- direction
                IF( ze3wv >= 0._wp ) THEN      ! case 1
-                  zmaxv =  ze3wv / e3w_n(ji,jj+1,ikv)
+                  zmaxv =  ze3wv / e3w(ji,jj+1,ikv,Kmm)
                   ! interpolated values of tracers
                   ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
@@ -296,5 +298,5 @@
                   pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
                ELSE                           ! case 2
-                  zmaxv =  -ze3wv / e3w_n(ji,jj,ikv)
+                  zmaxv =  -ze3wv / e3w(ji,jj,ikv,Kmm)
                   ! interpolated values of tracers
                   ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
@@ -318,12 +320,12 @@
                iku = mbku(ji,jj)
                ikv = mbkv(ji,jj)
-               ze3wu = gdept_n(ji+1,jj,iku) - gdept_n(ji,jj,iku)
-               ze3wv = gdept_n(ji,jj+1,ikv) - gdept_n(ji,jj,ikv)
-               !
-               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_n(ji  ,jj,iku)    ! i-direction: case 1
-               ELSE                        ;   zhi(ji,jj) = gdept_n(ji+1,jj,iku)    ! -     -      case 2
-               ENDIF
-               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_n(ji,jj  ,ikv)    ! j-direction: case 1
-               ELSE                        ;   zhj(ji,jj) = gdept_n(ji,jj+1,ikv)    ! -     -      case 2
+               ze3wu = gdept(ji+1,jj,iku,Kmm) - gdept(ji,jj,iku,Kmm)
+               ze3wv = gdept(ji,jj+1,ikv,Kmm) - gdept(ji,jj,ikv,Kmm)
+               !
+               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept(ji  ,jj,iku,Kmm)    ! i-direction: case 1
+               ELSE                        ;   zhi(ji,jj) = gdept(ji+1,jj,iku,Kmm)    ! -     -      case 2
+               ENDIF
+               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept(ji,jj  ,ikv,Kmm)    ! j-direction: case 1
+               ELSE                        ;   zhj(ji,jj) = gdept(ji,jj+1,ikv,Kmm)    ! -     -      case 2
                ENDIF
 
@@ -340,6 +342,6 @@
                iku = mbku(ji,jj)
                ikv = mbkv(ji,jj)
-               ze3wu = gdept_n(ji+1,jj,iku) - gdept_n(ji,jj,iku)
-               ze3wv = gdept_n(ji,jj+1,ikv) - gdept_n(ji,jj,ikv)
+               ze3wu = gdept(ji+1,jj,iku,Kmm) - gdept(ji,jj,iku,Kmm)
+               ze3wv = gdept(ji,jj+1,ikv,Kmm) - gdept(ji,jj,ikv,Kmm)
 
                IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1
@@ -369,10 +371,10 @@
                ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj
                ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0
-               ze3wu  =  gdept_n(ji,jj,iku) - gdept_n(ji+1,jj,iku)
-               ze3wv  =  gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv) 
+               ze3wu  =  gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm)
+               ze3wv  =  gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm) 
 
                ! i- direction
                IF( ze3wu >= 0._wp ) THEN      ! case 1
-                  zmaxu = ze3wu / e3w_n(ji+1,jj,ikup1)
+                  zmaxu = ze3wu / e3w(ji+1,jj,ikup1,Kmm)
                   ! interpolated values of tracers
                   zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikup1,jn) - pta(ji+1,jj,iku,jn) )
@@ -380,5 +382,5 @@
                   pgtui(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
                ELSE                           ! case 2
-                  zmaxu = - ze3wu / e3w_n(ji,jj,ikup1)
+                  zmaxu = - ze3wu / e3w(ji,jj,ikup1,Kmm)
                   ! interpolated values of tracers
                   zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikup1,jn) - pta(ji,jj,iku,jn) )
@@ -389,5 +391,5 @@
                ! j- direction
                IF( ze3wv >= 0._wp ) THEN      ! case 1
-                  zmaxv =  ze3wv / e3w_n(ji,jj+1,ikvp1)
+                  zmaxv =  ze3wv / e3w(ji,jj+1,ikvp1,Kmm)
                   ! interpolated values of tracers
                   ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvp1,jn) - pta(ji,jj+1,ikv,jn) )
@@ -395,5 +397,5 @@
                   pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
                ELSE                           ! case 2
-                  zmaxv =  - ze3wv / e3w_n(ji,jj,ikvp1)
+                  zmaxv =  - ze3wv / e3w(ji,jj,ikvp1,Kmm)
                   ! interpolated values of tracers
                   ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvp1,jn) - pta(ji,jj,ikv,jn) )
@@ -416,13 +418,13 @@
                iku = miku(ji,jj)
                ikv = mikv(ji,jj)
-               ze3wu  =  gdept_n(ji,jj,iku) - gdept_n(ji+1,jj,iku)
-               ze3wv  =  gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv) 
-               !
-               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_n(ji  ,jj,iku)    ! i-direction: case 1
-               ELSE                        ;   zhi(ji,jj) = gdept_n(ji+1,jj,iku)    ! -     -      case 2
-               ENDIF
-
-               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_n(ji,jj  ,ikv)    ! j-direction: case 1
-               ELSE                        ;   zhj(ji,jj) = gdept_n(ji,jj+1,ikv)    ! -     -      case 2
+               ze3wu  =  gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm)
+               ze3wv  =  gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm) 
+               !
+               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept(ji  ,jj,iku,Kmm)    ! i-direction: case 1
+               ELSE                        ;   zhi(ji,jj) = gdept(ji+1,jj,iku,Kmm)    ! -     -      case 2
+               ENDIF
+
+               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept(ji,jj  ,ikv,Kmm)    ! j-direction: case 1
+               ELSE                        ;   zhj(ji,jj) = gdept(ji,jj+1,ikv,Kmm)    ! -     -      case 2
                ENDIF
 
@@ -437,6 +439,6 @@
                iku = miku(ji,jj) 
                ikv = mikv(ji,jj) 
-               ze3wu  =  gdept_n(ji,jj,iku) - gdept_n(ji+1,jj,iku)
-               ze3wv  =  gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv) 
+               ze3wu  =  gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm)
+               ze3wv  =  gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm) 
 
                IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj      ) - prd(ji,jj,iku) ) ! i: 1
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRD/trdpen.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRD/trdpen.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRD/trdpen.F90	(revision 10954)
@@ -78,5 +78,5 @@
       IF( kt /= nkstp ) THEN     ! full eos: set partial derivatives at the 1st call of kt time step
          nkstp = kt
-         CALL eos_pen( ts(:,:,:,:,Kmm), rab_PE, zpe )
+         CALL eos_pen( ts(:,:,:,:,Kmm), rab_PE, zpe, Kmm )
          CALL iom_put( "alphaPE", rab_pe(:,:,:,jp_tem) )
          CALL iom_put( "betaPE" , rab_pe(:,:,:,jp_sal) )
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfosm.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfosm.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfosm.F90	(revision 10954)
@@ -1584,6 +1584,6 @@
      ALLOCATE( imld_rst(jpi,jpj) )
      ! w-level of the mixing and mixed layers
-     CALL eos_rab( ts(:,:,:,:,Kmm), rab_n )
-     CALL bn2(ts(:,:,:,:,Kmm), rab_n, rn2)
+     CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm )
+     CALL bn2(ts(:,:,:,:,Kmm), rab_n, rn2, Kmm)
      imld_rst(:,:)  = nlb10         ! Initialization to the number of w ocean point
      hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90	(revision 10954)
@@ -435,5 +435,5 @@
             CALL dia_obs( nit000 - 1, Nnn )              ! Observation operator for restart
          ENDIF     
-         IF( lk_asminc )   CALL asm_inc_init          ! Assimilation increments
+         IF( lk_asminc )   CALL asm_inc_init( Nnn )   ! Assimilation increments
          !
          RETURN                                       ! end of initialization
@@ -497,5 +497,5 @@
 
       !                                      ! Assimilation increments
-      IF( lk_asminc    )   CALL asm_inc_init    ! Initialize assimilation increments
+      IF( lk_asminc    )   CALL asm_inc_init( Nnn )        ! Initialize assimilation increments
       !
       IF(lwp) WRITE(numout,cform_aaa)           ! Flag AAAAAAA
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10954)
@@ -131,8 +131,8 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !  THERMODYNAMICS
-                         CALL eos_rab( tsb, rab_b )       ! before local thermal/haline expension ratio at T-points
-                         CALL eos_rab( tsn, rab_n )       ! now    local thermal/haline expension ratio at T-points
-                         CALL bn2    ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency
-                         CALL bn2    ( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency
+                         CALL eos_rab( tsb, rab_b, Nnn )       ! before local thermal/haline expension ratio at T-points
+                         CALL eos_rab( tsn, rab_n, Nnn )       ! now    local thermal/haline expension ratio at T-points
+                         CALL bn2    ( tsb, rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
+                         CALL bn2    ( tsn, rab_n, rn2, Nnn  ) ! now    Brunt-Vaisala frequency
 
       !  VERTICAL PHYSICS
@@ -144,11 +144,11 @@
                          CALL eos( tsb, rhd, gdept_0(:,:,:) )               ! before in situ density
 
-         IF( ln_zps .AND. .NOT. ln_isfcav)                               &
-            &            CALL zps_hde    ( kstp, jpts, tsb, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
-            &                                          rhd, gru , grv    )  ! of t, s, rd at the last ocean level
-
-         IF( ln_zps .AND.       ln_isfcav)                               &
-            &            CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
-            &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
+         IF( ln_zps .AND. .NOT. ln_isfcav)                                    &
+            &            CALL zps_hde    ( kstp, Nnn, jpts, tsb, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
+            &                                          rhd, gru , grv    )       ! of t, s, rd at the last ocean level
+
+         IF( ln_zps .AND.       ln_isfcav)                                                &
+            &            CALL zps_hde_isf( kstp, Nnn, jpts, tsb, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+            &                                          rhd, gru , grv , grui, grvi   )       ! of t, s, rd at the first ocean level
          IF( ln_traldf_triad ) THEN 
                          CALL ldf_slp_triad( kstp, Nbb, Nnn )             ! before slope for triad operator
@@ -175,10 +175,10 @@
 !!                                         but ensures reproductible results
 !!                                         with previous versions using split-explicit free surface          
-            IF( ln_zps .AND. .NOT. ln_isfcav )                               &
-               &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,   &  ! Partial steps: before horizontal gradient
-               &                                          rhd, gru , grv     )  ! of t, s, rd at the last ocean level
-            IF( ln_zps .AND.       ln_isfcav )                                          &
-               &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
-               &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
+            IF( ln_zps .AND. .NOT. ln_isfcav )                                    &
+               &            CALL zps_hde    ( kstp, Nnn, jpts, tsn, gtsu, gtsv,   &  ! Partial steps: before horizontal gradient
+               &                                          rhd, gru , grv     )       ! of t, s, rd at the last ocean level
+            IF( ln_zps .AND.       ln_isfcav )                                               &
+               &            CALL zps_hde_isf( kstp, Nnn, jpts, tsn, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+               &                                          rhd, gru , grv , grui, grvi   )       ! of t, s, rd at the first ocean level
 !!jc: fs simplification
                             
@@ -241,10 +241,10 @@
 
       IF(  lk_asminc .AND. ln_asmiau .AND. &
-         & ln_trainc )   CALL tra_asm_inc   ( kstp )  ! apply tracer assimilation increment
+         & ln_trainc )   CALL tra_asm_inc   ( kstp, Nnn       )  ! apply tracer assimilation increment
                          CALL tra_sbc       ( kstp, Nnn, Nrhs )  ! surface boundary condition
       IF( ln_traqsr  )   CALL tra_qsr       ( kstp, Nnn, Nrhs )  ! penetrative solar radiation qsr
       IF( ln_trabbc  )   CALL tra_bbc       ( kstp, Nnn, Nrhs )  ! bottom heat flux
-      IF( ln_trabbl  )   CALL tra_bbl       ( kstp, Nnn, Nrhs )  ! advective (and/or diffusive) bottom boundary layer scheme
-      IF( ln_tradmp  )   CALL tra_dmp       ( kstp, Nnn, Nrhs )  ! internal damping trends
+      IF( ln_trabbl  )   CALL tra_bbl       ( kstp, Nbb, Nnn, Nrhs )  ! advective (and/or diffusive) bottom boundary layer scheme
+      IF( ln_tradmp  )   CALL tra_dmp       ( kstp, Nbb, Nnn, Nrhs )  ! internal damping trends
       IF( ln_bdy     )   CALL bdy_tra_dmp   ( kstp )  ! bdy damping trends
 #if defined key_agrif
@@ -255,6 +255,6 @@
       IF( ln_zdfosm  )   CALL tra_osm( kstp,      Nnn      , ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS
       IF( lrst_oce .AND. ln_zdfosm ) &
-           &             CALL osm_rst( kstp, Nnn, 'WRITE' )! write OSMOSIS outputs + wn (so must do here) to restarts
-                         CALL tra_ldf( kstp, Nnn, Nrhs )   ! lateral mixing
+           &             CALL osm_rst( kstp,      Nnn, 'WRITE' )         ! write OSMOSIS outputs + wn (so must do here) to restarts
+                         CALL tra_ldf( kstp, Nbb, Nnn, Nrhs )            ! lateral mixing
 
 !!gm : why CALL to dia_ptr has been moved here??? (use trends info?)
@@ -281,6 +281,6 @@
 !! 
 !!jc2: dynnxt must be the latest call. e3t_b are indeed updated in that routine
-                         CALL tra_nxt       ( kstp, Nnn, Nrhs )  ! finalize (bcs) tracer fields at next time step and swap
-                         CALL dyn_nxt       ( kstp, Nnn )        ! finalize (bcs) velocities at next time step and swap (always called after tra_nxt)
+                         CALL tra_nxt       ( kstp, Nbb, Nnn, Nrhs )  ! finalize (bcs) tracer fields at next time step and swap
+                         CALL dyn_nxt       ( kstp,      Nnn )        ! finalize (bcs) velocities at next time step and swap (always called after tra_nxt)
                          CALL ssh_swp       ( kstp )  ! swap of sea surface height
       IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/dtadyn.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/dtadyn.F90	(revision 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/dtadyn.F90	(revision 10954)
@@ -159,6 +159,6 @@
       !
       CALL eos    ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
-      CALL eos_rab( ts(:,:,:,:,Kmm), rab_n )       ! now    local thermal/haline expension ratio at T-points
-      CALL bn2    ( ts(:,:,:,:,Kmm), rab_n, rn2 ) ! before Brunt-Vaisala frequency need for zdfmxl
+      CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm )       ! now    local thermal/haline expension ratio at T-points
+      CALL bn2    ( ts(:,:,:,:,Kmm), rab_n, rn2, Kmm )  ! before Brunt-Vaisala frequency need for zdfmxl
 
       rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl
@@ -785,13 +785,13 @@
       IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)
          CALL eos    ( pts, rhd, rhop, gdept_0(:,:,:) )
-         CALL eos_rab( pts, rab_n )       ! now local thermal/haline expension ratio at T-points
-         CALL bn2    ( pts, rab_n, rn2  ) ! now    Brunt-Vaisala
+         CALL eos_rab( pts, rab_n, Kmm )       ! now local thermal/haline expension ratio at T-points
+         CALL bn2    ( pts, rab_n, rn2, Kmm  ) ! now    Brunt-Vaisala
 
       ! Partial steps: before Horizontal DErivative
       IF( ln_zps  .AND. .NOT. ln_isfcav)                            &
-         &            CALL zps_hde    ( kt, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
+         &            CALL zps_hde    ( kt, Kmm, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
          &                                        rhd, gru , grv    )  ! of t, s, rd at the last ocean level
       IF( ln_zps .AND.        ln_isfcav)                            &
-         &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, gtui, gtvi, &  ! Partial steps for top cell (ISF)
+         &            CALL zps_hde_isf( kt, Kmm, jpts, pts, gtsu, gtsv, gtui, gtvi, &  ! Partial steps for top cell (ISF)
          &                                        rhd, gru , grv , grui, grvi )  ! of t, s, rd at the first ocean level
 
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 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcadv.F90	(revision 10954)
@@ -118,5 +118,5 @@
             &              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' )  ! add the mle transport
+         IF( ln_mle    )   CALL tra_mle_trp( kt, nittrc000, zun, zvn, zwn, 'TRC', Kmm       )  ! add the mle transport
          !
       ENDIF
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 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcbbl.F90	(revision 10954)
@@ -36,5 +36,5 @@
 CONTAINS
 
-   SUBROUTINE trc_bbl( kt, Kmm, Krhs )
+   SUBROUTINE trc_bbl( kt, Kbb, Kmm, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE bbl  ***
@@ -46,5 +46,5 @@
       !!----------------------------------------------------------------------  
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step 
-      INTEGER, INTENT( in ) ::   Kmm, Krhs  ! time level indices
+      INTEGER, INTENT( in ) ::   Kbb, Kmm, Krhs  ! time level indices
       INTEGER :: jn                   ! loop index
       CHARACTER (len=22) :: charout
@@ -55,6 +55,6 @@
       !
       IF( .NOT. l_offline .AND. nn_dttrc == 1 ) THEN
-         CALL bbl( kt, nittrc000, 'TRC' )      ! Online coupling with dynamics  : Computation of bbl coef and bbl transport
-         l_bbl = .FALSE.                       ! Offline coupling with dynamics : Read bbl coef and bbl transport from input files
+         CALL bbl( kt, nittrc000, 'TRC', Kbb, Kmm )  ! Online coupling with dynamics  : Computation of bbl coef and bbl transport
+         l_bbl = .FALSE.                             ! Offline coupling with dynamics : Read bbl coef and bbl transport from input files
       ENDIF
 
@@ -67,5 +67,5 @@
       IF( nn_bbl_ldf == 1 ) THEN
          !
-         CALL tra_bbl_dif( trb, tra, jptra )  
+         CALL tra_bbl_dif( trb, tra, jptra, Kmm )  
          IF( ln_ctl )   THEN
             WRITE(charout, FMT="(' bbl_dif')")  ;  CALL prt_ctl_trc_info(charout)
@@ -78,5 +78,5 @@
       IF( nn_bbl_adv /= 0 ) THEN
          !
-         CALL tra_bbl_adv( trb, tra, jptra )  
+         CALL tra_bbl_adv( trb, tra, jptra, Kmm )  
          IF( ln_ctl )   THEN
             WRITE(charout, FMT="(' bbl_adv')")  ;  CALL prt_ctl_trc_info(charout)
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 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcldf.F90	(revision 10954)
@@ -94,7 +94,7 @@
       !
       CASE ( np_lap   )                               ! iso-level laplacian
-         CALL tra_ldf_lap  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb,      tra, jptra,    1     )
+         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     )
+         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 )
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 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcnxt.F90	(revision 10954)
@@ -54,5 +54,5 @@
 CONTAINS
 
-   SUBROUTINE trc_nxt( kt, Kmm, Krhs )
+   SUBROUTINE trc_nxt( kt, Kbb, Kmm, Krhs )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE trcnxt  ***
@@ -79,5 +79,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
       !
       INTEGER  ::   jk, jn   ! dummy loop indices
@@ -157,7 +157,7 @@
       ELSE     
          IF( .NOT. l_offline ) THEN ! Leap-Frog + Asselin filter time stepping
-            IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt, nittrc000,         'TRC', trb, trn, tra, jptra )  !     linear ssh
-            ELSE                   ;   CALL tra_nxt_vvl( kt, 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', 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
             ENDIF
          ELSE
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 10946)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90	(revision 10954)
@@ -62,5 +62,5 @@
          !
                                 CALL trc_sbc    ( kt, Kmm, Krhs )      ! surface boundary condition
-         IF( ln_trabbl )        CALL trc_bbl    ( kt, Kmm, Krhs )      ! advective (and/or diffusive) bottom boundary layer scheme
+         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
@@ -68,6 +68,6 @@
          !                                                             ! Partial top/bottom cell: GRADh( trb )  
          IF( ln_zps ) THEN
-           IF( ln_isfcav ) THEN ; CALL zps_hde_isf( kt, jptra, trb, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )  ! both top & bottom
-           ELSE                 ; CALL zps_hde    ( kt, jptra, trb, gtru, gtrv )                                      !  only bottom
+           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
            ENDIF
          ENDIF
@@ -78,5 +78,5 @@
 #endif
                                 CALL trc_zdf( kt, Kbb, Kmm, Krhs, tr, Kaa  )  ! vert. mixing & after tracer	==> after
-                                CALL trc_nxt    ( kt, Kmm, Krhs )        ! tracer fields at next time step     
+                                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
@@ -87,5 +87,5 @@
          IF( ln_trcdmp )        CALL trc_dmp( kt, Kmm, Krhs )            ! internal damping trends
                                 CALL trc_zdf( kt, Kbb, Kmm, Krhs, tr, Kaa  )  ! vert. mixing & after tracer	==> after
-                                CALL trc_nxt( kt, Kmm, Krhs )            ! tracer fields at next time step     
+                                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
          !
