Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIU/step_diu.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIU/step_diu.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIU/step_diu.F90	(revision 10922)
@@ -49,4 +49,5 @@
       INTEGER ::   indic    ! error indicator if < 0 
       REAL(wp), DIMENSION(jpi,jpj) :: z_fvel_bkginc, z_hflux_bkginc     
+      INTEGER :: Nbb, Nnn, Naa, Nrhs    ! local definitions as placeholders for now
       !! --------------------------------------------------------------------- 
       
@@ -60,5 +61,5 @@
          ENDIF
        
-            CALL sbc    ( kstp )                      ! Sea Boundary Conditions 
+            CALL sbc    ( kstp, Nbb, Nnn )            ! Sea Boundary Conditions 
       ENDIF
      
@@ -78,5 +79,6 @@
 
       IF( ln_diurnal_only ) THEN
-         IF( ln_diaobs )         CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
+         ! WILL HAVE TO INCREMENT Nbb and Nnn here in ln_diurnal_only case !
+         IF( ln_diaobs )         CALL dia_obs( kstp, Nnn )    ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
      
          !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
@@ -84,5 +86,5 @@
          !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
          IF( kstp == nit000   )   CALL iom_close( numror )     ! close input  ocean restart file 
-         IF( lrst_oce         )   CALL rst_write    ( kstp )   ! write output ocean restart file
+         IF( lrst_oce         )   CALL rst_write    ( kstp, Nbb, Nnn )   ! write output ocean restart file
      
          IF( ln_timing .AND.  kstp == nit000  )   CALL timing_reset 
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/istate.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/istate.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/istate.F90	(revision 10922)
@@ -51,5 +51,5 @@
 CONTAINS
 
-   SUBROUTINE istate_init
+   SUBROUTINE istate_init( Kbb, Kmm )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE istate_init  ***
@@ -57,4 +57,5 @@
       !! ** Purpose :   Initialization of the dynamics and tracer fields.
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
       INTEGER ::   ji, jj, jk   ! dummy loop indices
 !!gm see comment further down
@@ -85,5 +86,5 @@
       IF( ln_rstart ) THEN                    ! Restart from a file
          !                                    ! -------------------
-         CALL rst_read                        ! Read the restart file
+         CALL rst_read( Kbb, Kmm )            ! Read the restart file
          IF (ln_iscpl)       CALL iscpl_stp   ! extrapolate restart to wet and dry
          CALL day_init                        ! model calendar (using both namelist and restart infos)
@@ -124,5 +125,5 @@
          vn   (:,:,:)   = vb  (:,:,:)
          hdivn(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level
-         CALL div_hor( 0 )                    ! compute interior hdivn value  
+         CALL div_hor( 0, Kmm )               ! compute interior hdivn value  
 !!gm                                    hdivn(:,:,:) = 0._wp
 
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/divhor.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/divhor.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/divhor.F90	(revision 10922)
@@ -48,5 +48,5 @@
 CONTAINS
 
-   SUBROUTINE div_hor( kt )
+   SUBROUTINE div_hor( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE div_hor  ***
@@ -61,4 +61,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kmm  ! ocean time-level index
       !
       INTEGER  ::   ji, jj, jk    ! dummy loop indices
@@ -94,5 +95,5 @@
 #endif
       !
-      IF( ln_rnf )   CALL sbc_rnf_div( hdivn )              !==  runoffs    ==!   (update hdivn field)
+      IF( ln_rnf )   CALL sbc_rnf_div( hdivn, Kmm )                     !==  runoffs    ==!   (update hdivn field)
       !
 #if defined key_asminc 
@@ -100,7 +101,7 @@
       ! 
 #endif
-      IF( ln_isf )   CALL sbc_isf_div( hdivn )      !==  ice shelf  ==!   (update hdivn field)
+      IF( ln_isf )   CALL sbc_isf_div( hdivn, Kmm )                     !==  ice shelf  ==!   (update hdivn field)
       !
-      IF( ln_iscpl .AND. ln_hsb )   CALL iscpl_div( hdivn ) !==  ice sheet  ==!   (update hdivn field)
+      IF( ln_iscpl .AND. ln_hsb )   CALL iscpl_div( hdivn )             !==  ice sheet  ==!   (update hdivn field)
       !
       CALL lbc_lnk( 'divhor', hdivn, 'T', 1. )   !   (no sign change)
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/sshwzv.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/sshwzv.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/sshwzv.F90	(revision 10922)
@@ -54,5 +54,5 @@
 CONTAINS
 
-   SUBROUTINE ssh_nxt( kt )
+   SUBROUTINE ssh_nxt( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE ssh_nxt  ***
@@ -69,4 +69,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! time step
+      INTEGER, INTENT(in) ::   Kmm  ! time level index
       ! 
       INTEGER  ::   jk            ! dummy loop indice
@@ -94,5 +95,5 @@
       ENDIF
 
-      CALL div_hor( kt )                               ! Horizontal divergence
+      CALL div_hor( kt, Kmm )                          ! Horizontal divergence
       !
       zhdiv(:,:) = 0._wp
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/IOM/restart.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/IOM/restart.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/IOM/restart.F90	(revision 10922)
@@ -131,5 +131,5 @@
 
 
-   SUBROUTINE rst_write( kt )
+   SUBROUTINE rst_write( kt, Kbb, Kmm )
       !!---------------------------------------------------------------------
       !!                   ***  ROUTINE rstwrite  ***
@@ -140,5 +140,6 @@
       !!              file, save fields which are necessary for restart
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step
+      INTEGER, INTENT(in) ::   kt         ! ocean time-step
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
       !!----------------------------------------------------------------------
                      IF(lwxios) CALL iom_swap(      cwxios_context          )
@@ -147,14 +148,14 @@
 
       IF ( .NOT. ln_diurnal_only ) THEN
-                     CALL iom_rstput( kt, nitrst, numrow, 'ub'     , ub, ldxios = lwxios        )     ! before fields
-                     CALL iom_rstput( kt, nitrst, numrow, 'vb'     , vb, ldxios = lwxios        )
-                     CALL iom_rstput( kt, nitrst, numrow, 'tb'     , tsb(:,:,:,jp_tem), ldxios = lwxios )
-                     CALL iom_rstput( kt, nitrst, numrow, 'sb'     , tsb(:,:,:,jp_sal), ldxios = lwxios )
+                     CALL iom_rstput( kt, nitrst, numrow, 'ub'     , uu(:,:,:,Kbb), ldxios = lwxios        )     ! before fields
+                     CALL iom_rstput( kt, nitrst, numrow, 'vb'     , vv(:,:,:,Kbb), ldxios = lwxios        )
+                     CALL iom_rstput( kt, nitrst, numrow, 'tb'     , ts(:,:,:,jp_tem,Kbb), ldxios = lwxios )
+                     CALL iom_rstput( kt, nitrst, numrow, 'sb'     , ts(:,:,:,jp_sal,Kbb), ldxios = lwxios )
                      CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb, ldxios = lwxios      )
                      !
-                     CALL iom_rstput( kt, nitrst, numrow, 'un'     , un, ldxios = lwxios        )     ! now fields
-                     CALL iom_rstput( kt, nitrst, numrow, 'vn'     , vn, ldxios = lwxios        )
-                     CALL iom_rstput( kt, nitrst, numrow, 'tn'     , tsn(:,:,:,jp_tem), ldxios = lwxios )
-                     CALL iom_rstput( kt, nitrst, numrow, 'sn'     , tsn(:,:,:,jp_sal), ldxios = lwxios )
+                     CALL iom_rstput( kt, nitrst, numrow, 'un'     , uu(:,:,:,Kmm), ldxios = lwxios        )     ! now fields
+                     CALL iom_rstput( kt, nitrst, numrow, 'vn'     , vv(:,:,:,Kmm), ldxios = lwxios        )
+                     CALL iom_rstput( kt, nitrst, numrow, 'tn'     , ts(:,:,:,jp_tem,Kmm), ldxios = lwxios )
+                     CALL iom_rstput( kt, nitrst, numrow, 'sn'     , ts(:,:,:,jp_sal,Kmm), ldxios = lwxios )
                      CALL iom_rstput( kt, nitrst, numrow, 'sshn'   , sshn, ldxios = lwxios      )
                      CALL iom_rstput( kt, nitrst, numrow, 'rhop'   , rhop, ldxios = lwxios      )
@@ -165,8 +166,8 @@
                      CALL iom_rstput( kt, nitrst, numrow, 'vmask'  , vmask, ldxios = lwxios ) ! need to correct barotropic velocity
                      CALL iom_rstput( kt, nitrst, numrow, 'smask'  , ssmask, ldxios = lwxios) ! need to correct barotropic velocity
-                     CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios )   ! need to compute temperature correction
-                     CALL iom_rstput( kt, nitrst, numrow, 'e3u_n', e3u_n(:,:,:), ldxios = lwxios )   ! need to compute bt conservation
-                     CALL iom_rstput( kt, nitrst, numrow, 'e3v_n', e3v_n(:,:,:), ldxios = lwxios )   ! need to compute bt conservation
-                     CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', gdepw_n(:,:,:), ldxios = lwxios ) ! need to compute extrapolation if vvl
+                     CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios )   ! need to compute temperature correction
+                     CALL iom_rstput( kt, nitrst, numrow, 'e3u_n', e3u(:,:,:,Kmm), ldxios = lwxios )   ! need to compute bt conservation
+                     CALL iom_rstput( kt, nitrst, numrow, 'e3v_n', e3v(:,:,:,Kmm), ldxios = lwxios )   ! need to compute bt conservation
+                     CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', gdepw(:,:,:,Kmm), ldxios = lwxios ) ! need to compute extrapolation if vvl
                   END IF
       ENDIF
@@ -238,5 +239,5 @@
 
 
-   SUBROUTINE rst_read
+   SUBROUTINE rst_read( Kbb, Kmm )
       !!---------------------------------------------------------------------- 
       !!                   ***  ROUTINE rst_read  ***
@@ -246,4 +247,5 @@
       !! ** Method  :   Read in restart.nc file fields which are necessary for restart
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
       REAL(wp) ::   zrdt
       INTEGER  ::   jk
@@ -268,38 +270,38 @@
          rhop = rau0
          CALL iom_get( numror, jpdom_autoglo, 'tn'     , w3d, ldxios = lrxios ) 
-         tsn(:,:,1,jp_tem) = w3d(:,:,1)
+         ts(:,:,1,jp_tem,Kmm) = w3d(:,:,1)
          RETURN 
       ENDIF  
       
       IF( iom_varid( numror, 'ub', ldstop = .FALSE. ) > 0 ) THEN
-         CALL iom_get( numror, jpdom_autoglo, 'ub'     , ub, ldxios = lrxios                )   ! before fields
-         CALL iom_get( numror, jpdom_autoglo, 'vb'     , vb, ldxios = lrxios                )
-         CALL iom_get( numror, jpdom_autoglo, 'tb'     , tsb(:,:,:,jp_tem), ldxios = lrxios )
-         CALL iom_get( numror, jpdom_autoglo, 'sb'     , tsb(:,:,:,jp_sal), ldxios = lrxios )
-         CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb, ldxios = lrxios              )
+         CALL iom_get( numror, jpdom_autoglo, 'ub'     , uu(:,:,:,Kbb), ldxios = lrxios        )   ! before fields
+         CALL iom_get( numror, jpdom_autoglo, 'vb'     , vv(:,:,:,Kbb), ldxios = lrxios        )
+         CALL iom_get( numror, jpdom_autoglo, 'tb'     , ts(:,:,:,jp_tem,Kbb), ldxios = lrxios )
+         CALL iom_get( numror, jpdom_autoglo, 'sb'     , ts(:,:,:,jp_sal,Kbb), ldxios = lrxios )
+         CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb, ldxios = lrxios                 )
       ELSE
          neuler = 0
       ENDIF
       !
-      CALL iom_get( numror, jpdom_autoglo, 'un'     , un, ldxios = lrxios )   ! now    fields
-      CALL iom_get( numror, jpdom_autoglo, 'vn'     , vn, ldxios = lrxios )
-      CALL iom_get( numror, jpdom_autoglo, 'tn'     , tsn(:,:,:,jp_tem), ldxios = lrxios )
-      CALL iom_get( numror, jpdom_autoglo, 'sn'     , tsn(:,:,:,jp_sal), ldxios = lrxios )
+      CALL iom_get( numror, jpdom_autoglo, 'un'     , uu(:,:,:,Kmm), ldxios = lrxios        )   ! now    fields
+      CALL iom_get( numror, jpdom_autoglo, 'vn'     , vv(:,:,:,Kmm), ldxios = lrxios        )
+      CALL iom_get( numror, jpdom_autoglo, 'tn'     , ts(:,:,:,jp_tem,Kmm), ldxios = lrxios )
+      CALL iom_get( numror, jpdom_autoglo, 'sn'     , ts(:,:,:,jp_sal,Kmm), ldxios = lrxios )
       CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios )
       IF( iom_varid( numror, 'rhop', ldstop = .FALSE. ) > 0 ) THEN
          CALL iom_get( numror, jpdom_autoglo, 'rhop'   , rhop, ldxios = lrxios )   ! now    potential density
       ELSE
-         CALL eos( tsn, rhd, rhop, gdept_n(:,:,:) )   
+         CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, gdept(:,:,:,Kmm) )   
       ENDIF
       !
       IF( neuler == 0 ) THEN                                  ! Euler restart (neuler=0)
-         tsb  (:,:,:,:) = tsn  (:,:,:,:)                             ! all before fields set to now values
-         ub   (:,:,:)   = un   (:,:,:)
-         vb   (:,:,:)   = vn   (:,:,:)
-         sshb (:,:)     = sshn (:,:)
+         ts   (:,:,:,:,Kbb) = ts   (:,:,:,:,Kmm)              ! all before fields set to now values
+         uu   (:,:,:,Kbb)   = uu   (:,:,:,Kmm)
+         vv   (:,:,:,Kbb)   = vv   (:,:,:,Kmm)
+         sshb (:,:)         = sshn (:,:)
          !
          IF( .NOT.ln_linssh ) THEN
             DO jk = 1, jpk
-               e3t_b(:,:,jk) = e3t_n(:,:,jk)
+               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm)
             END DO
          ENDIF
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/LDF/ldfdyn.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/LDF/ldfdyn.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/LDF/ldfdyn.F90	(revision 10922)
@@ -339,5 +339,5 @@
 
 
-   SUBROUTINE ldf_dyn( kt )
+   SUBROUTINE ldf_dyn( kt, Kbb )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ldf_dyn  ***
@@ -357,4 +357,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! time step index
+      INTEGER, INTENT(in) ::   Kbb  ! ocean time level indices
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -373,7 +374,7 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
-                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
-                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
+                     zu2pv2_ij_p1 = uu(ji  ,jj+1,jk,Kbb) * uu(ji  ,jj+1,jk,Kbb) + vv(ji+1,jj  ,jk,Kbb) * vv(ji+1,jj  ,jk,Kbb)
+                     zu2pv2_ij    = uu(ji  ,jj  ,jk,Kbb) * uu(ji  ,jj  ,jk,Kbb) + vv(ji  ,jj  ,jk,Kbb) * vv(ji  ,jj  ,jk,Kbb)
+                     zu2pv2_ij_m1 = uu(ji-1,jj  ,jk,Kbb) * uu(ji-1,jj  ,jk,Kbb) + vv(ji  ,jj-1,jk,Kbb) * vv(ji  ,jj-1,jk,Kbb)
                      zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) )
                      zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) )
@@ -387,7 +388,7 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
-                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
-                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
+                     zu2pv2_ij_p1 = uu(ji  ,jj+1,jk,Kbb) * uu(ji  ,jj+1,jk,Kbb) + vv(ji+1,jj  ,jk,Kbb) * vv(ji+1,jj  ,jk,Kbb)
+                     zu2pv2_ij    = uu(ji  ,jj  ,jk,Kbb) * uu(ji  ,jj  ,jk,Kbb) + vv(ji  ,jj  ,jk,Kbb) * vv(ji  ,jj  ,jk,Kbb)
+                     zu2pv2_ij_m1 = uu(ji-1,jj  ,jk,Kbb) * uu(ji-1,jj  ,jk,Kbb) + vv(ji  ,jj-1,jk,Kbb) * vv(ji  ,jj-1,jk,Kbb)
                      zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) )
                      zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) )
@@ -415,7 +416,7 @@
                DO jj = 2, jpj
                   DO ji = 2, jpi
-                     zdb = ( (  ub(ji,jj,jk) * r1_e2u(ji,jj) -  ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) )  &
+                     zdb = ( (  uu(ji,jj,jk,Kbb) * r1_e2u(ji,jj) -  uu(ji-1,jj,jk,Kbb) * r1_e2u(ji-1,jj) )  &
                           &                  * r1_e1t(ji,jj) * e2t(ji,jj)                           &
-                          & - ( vb(ji,jj,jk) * r1_e1v(ji,jj) -  vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) )  &
+                          & - ( vv(ji,jj,jk,Kbb) * r1_e1v(ji,jj) -  vv(ji,jj-1,jk,Kbb) * r1_e1v(ji,jj-1) )  &
                           &                  * r1_e2t(ji,jj) * e1t(ji,jj)    ) * tmask(ji,jj,jk)
                      dtensq(ji,jj) = zdb * zdb
@@ -425,7 +426,7 @@
                DO jj = 1, jpjm1
                   DO ji = 1, jpim1
-                     zdb = ( (  ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) -  ub(ji,jj,jk) * r1_e1u(ji,jj) )  &
+                     zdb = ( (  uu(ji,jj+1,jk,Kbb) * r1_e1u(ji,jj+1) -  uu(ji,jj,jk,Kbb) * r1_e1u(ji,jj) )  &
                           &                    * r1_e2f(ji,jj)   * e1f(ji,jj)                       &
-                          & + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) -  vb(ji,jj,jk) * r1_e2v(ji,jj) )  &
+                          & + ( vv(ji+1,jj,jk,Kbb) * r1_e2v(ji+1,jj) -  vv(ji,jj,jk,Kbb) * r1_e2v(ji,jj) )  &
                           &                    * r1_e1f(ji,jj)   * e2f(ji,jj)  ) * fmask(ji,jj,jk)
                      dshesq(ji,jj) = zdb * zdb
@@ -436,7 +437,7 @@
                   DO ji = fs_2, fs_jpim1
                      !
-                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
-                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
-                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
+                     zu2pv2_ij_p1 = uu(ji  ,jj+1,jk,Kbb) * uu(ji  ,jj+1,jk,Kbb) + vv(ji+1,jj  ,jk,Kbb) * vv(ji+1,jj  ,jk,Kbb)
+                     zu2pv2_ij    = uu(ji  ,jj  ,jk,Kbb) * uu(ji  ,jj  ,jk,Kbb) + vv(ji  ,jj  ,jk,Kbb) * vv(ji  ,jj  ,jk,Kbb)
+                     zu2pv2_ij_m1 = uu(ji-1,jj  ,jk,Kbb) * uu(ji-1,jj  ,jk,Kbb) + vv(ji  ,jj-1,jk,Kbb) * vv(ji  ,jj-1,jk,Kbb)
                                                      ! T-point value
                      zdelta         = zcmsmag * esqt(ji,jj)                                        ! L^2 * (C_smag/pi)^2
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/LDF/ldfslp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/LDF/ldfslp.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/LDF/ldfslp.F90	(revision 10922)
@@ -81,5 +81,5 @@
 CONTAINS
 
-   SUBROUTINE ldf_slp( kt, prd, pn2 )
+   SUBROUTINE ldf_slp( kt, prd, pn2, Kbb, Kmm )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE ldf_slp  ***
@@ -107,4 +107,5 @@
       !!----------------------------------------------------------------------
       INTEGER , INTENT(in)                   ::   kt    ! ocean time-step index
+      INTEGER , INTENT(in)                   ::   Kbb, Kmm   ! ocean time level indices
       REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   prd   ! in situ density
       REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   pn2   ! Brunt-Vaisala frequency (locally ref.)
@@ -171,5 +172,5 @@
       !
       !                          !==   Slopes just below the mixed layer   ==!
-      CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr )        ! output: uslpml, vslpml, wslpiml, wslpjml
+      CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr, Kmm )        ! output: uslpml, vslpml, wslpiml, wslpjml
 
 
@@ -205,14 +206,14 @@
                !                                      ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
                !                                      ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
-               zbu = MIN(  zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u_n(ji,jj,jk)* ABS( zau )  )
-               zbv = MIN(  zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v_n(ji,jj,jk)* ABS( zav )  )
+               zbu = MIN(  zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u(ji,jj,jk,Kmm)* ABS( zau )  )
+               zbv = MIN(  zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v(ji,jj,jk,Kmm)* ABS( zav )  )
                !                                      ! uslp and vslp output in zwz and zww, resp.
                zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )
                zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) )
                ! thickness of water column between surface and level k at u/v point
-               zdepu = 0.5_wp * ( ( gdept_n (ji,jj,jk) + gdept_n (ji+1,jj,jk) )                            &
-                                - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj) ) - e3u_n(ji,jj,miku(ji,jj))   )
-               zdepv = 0.5_wp * ( ( gdept_n (ji,jj,jk) + gdept_n (ji,jj+1,jk) )                            &
-                                - 2 * MAX( risfdep(ji,jj), risfdep(ji,jj+1) ) - e3v_n(ji,jj,mikv(ji,jj))   )
+               zdepu = 0.5_wp * ( ( gdept (ji,jj,jk,Kmm) + gdept (ji+1,jj,jk,Kmm) )                            &
+                                - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj) ) - e3u(ji,jj,miku(ji,jj),Kmm)   )
+               zdepv = 0.5_wp * ( ( gdept (ji,jj,jk,Kmm) + gdept (ji,jj+1,jk,Kmm) )                            &
+                                - 2 * MAX( risfdep(ji,jj), risfdep(ji,jj+1) ) - e3v(ji,jj,mikv(ji,jj),Kmm)   )
                !
                zwz(ji,jj,jk) = ( ( 1._wp - zfi) * zau / ( zbu - zeps )                                     &
@@ -224,6 +225,6 @@
 !               zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp )
 !               zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp )
-!               zci = 0.5 * ( gdept_n(ji+1,jj,jk)+gdept_n(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) )
-!               zcj = 0.5 * ( gdept_n(ji,jj+1,jk)+gdept_n(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) )
+!               zci = 0.5 * ( gdept(ji+1,jj,jk,Kmm)+gdept(ji,jj,jk,Kmm) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) )
+!               zcj = 0.5 * ( gdept(ji,jj+1,jk,Kmm)+gdept(ji,jj,jk,Kmm) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) )
 !               zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk)
 !               zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk)
@@ -296,9 +297,9 @@
                !                                        ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
                !                                        ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
-               zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/e3w_n(ji,jj,jk)* ABS( zai )  )
-               zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w_n(ji,jj,jk)* ABS( zaj )  )
+               zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/e3w(ji,jj,jk,Kmm)* ABS( zai )  )
+               zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w(ji,jj,jk,Kmm)* ABS( zaj )  )
                !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.)
                zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )   ! zfk=1 in the ML otherwise zfk=0
-               zck = ( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj) - gdepw_n(ji,jj,mikt(ji,jj)), 10._wp )
+               zck = ( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm) ) / MAX( hmlp(ji,jj) - gdepw(ji,jj,mikt(ji,jj),Kmm), 10._wp )
                zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk  ) * wmask(ji,jj,jk)
                zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk  ) * wmask(ji,jj,jk)
@@ -375,5 +376,5 @@
 
 
-   SUBROUTINE ldf_slp_triad ( kt )
+   SUBROUTINE ldf_slp_triad ( kt, Kbb, Kmm )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE ldf_slp_triad  ***
@@ -390,4 +391,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt             ! ocean time-step index
+      INTEGER , INTENT(in)  ::   Kbb, Kmm       ! ocean time level indices
       !!
       INTEGER  ::   ji, jj, jk, jl, ip, jp, kp  ! dummy loop indices
@@ -419,8 +421,8 @@
             DO jj = 1, jpjm1                  ! NB: not masked ==>  a minimum value is set
                DO ji = 1, fs_jpim1            ! vector opt.
-                  zdit = ( tsb(ji+1,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )    ! i-gradient of T & S at u-point
-                  zdis = ( tsb(ji+1,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
-                  zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )    ! j-gradient of T & S at v-point
-                  zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
+                  zdit = ( ts(ji+1,jj,jk,jp_tem,Kbb) - ts(ji,jj,jk,jp_tem,Kbb) )    ! i-gradient of T & S at u-point
+                  zdis = ( ts(ji+1,jj,jk,jp_sal,Kbb) - ts(ji,jj,jk,jp_sal,Kbb) )
+                  zdjt = ( ts(ji,jj+1,jk,jp_tem,Kbb) - ts(ji,jj,jk,jp_tem,Kbb) )    ! j-gradient of T & S at v-point
+                  zdjs = ( ts(ji,jj+1,jk,jp_sal,Kbb) - ts(ji,jj,jk,jp_sal,Kbb) )
                   zdxrho_raw = ( - rab_b(ji+ip,jj   ,jk,jp_tem) * zdit + rab_b(ji+ip,jj   ,jk,jp_sal) * zdis ) * r1_e1u(ji,jj)
                   zdyrho_raw = ( - rab_b(ji   ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji   ,jj+jp,jk,jp_sal) * zdjs ) * r1_e2v(ji,jj)
@@ -452,6 +454,6 @@
                DO ji = 1, jpi                 ! vector opt.
                   IF( jk+kp > 1 ) THEN        ! k-gradient of T & S a jk+kp
-                     zdkt = ( tsb(ji,jj,jk+kp-1,jp_tem) - tsb(ji,jj,jk+kp,jp_tem) )
-                     zdks = ( tsb(ji,jj,jk+kp-1,jp_sal) - tsb(ji,jj,jk+kp,jp_sal) )
+                     zdkt = ( ts(ji,jj,jk+kp-1,jp_tem,Kbb) - ts(ji,jj,jk+kp,jp_tem,Kbb) )
+                     zdks = ( ts(ji,jj,jk+kp-1,jp_sal,Kbb) - ts(ji,jj,jk+kp,jp_sal,Kbb) )
                   ELSE
                      zdkt = 0._wp                                             ! 1st level gradient set to zero
@@ -460,5 +462,5 @@
                   zdzrho_raw = ( - rab_b(ji,jj,jk+kp,jp_tem) * zdkt & 
                              &   + rab_b(ji,jj,jk+kp,jp_sal) * zdks &
-                             & ) / e3w_n(ji,jj,jk+kp)  
+                             & ) / e3w(ji,jj,jk+kp,Kmm)  
                   zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw )    ! force zdzrho >= repsln
                  END DO
@@ -470,5 +472,5 @@
          DO ji = 1, jpi
             jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1     ! MIN in case ML depth is the ocean depth
-            z1_mlbw(ji,jj) = 1._wp / gdepw_n(ji,jj,jk)
+            z1_mlbw(ji,jj) = 1._wp / gdepw(ji,jj,jk,Kmm)
          END DO
       END DO
@@ -499,6 +501,6 @@
                      ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth)
                      zti_g_raw = (  zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp)      &
-                        &          - ( gdept_n(ji+1,jj,jk-kp) - gdept_n(ji,jj,jk-kp) ) * r1_e1u(ji,jj)  ) * umask(ji,jj,jk)
-                     ze3_e1    =  e3w_n(ji+ip,jj,jk-kp) * r1_e1u(ji,jj) 
+                        &          - ( gdept(ji+1,jj,jk-kp,Kmm) - gdept(ji,jj,jk-kp,Kmm) ) * r1_e1u(ji,jj)  ) * umask(ji,jj,jk)
+                     ze3_e1    =  e3w(ji+ip,jj,jk-kp,Kmm) * r1_e1u(ji,jj) 
                      zti_mlb(ji+ip,jj   ,1-ip,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1  , ABS( zti_g_raw ) ), zti_g_raw )
                   ENDIF
@@ -509,6 +511,6 @@
                   ELSE
                      ztj_g_raw = (  zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp)      &
-                        &      - ( gdept_n(ji,jj+1,jk-kp) - gdept_n(ji,jj,jk-kp) ) / e2v(ji,jj)  ) * vmask(ji,jj,jk)
-                     ze3_e2    =  e3w_n(ji,jj+jp,jk-kp) / e2v(ji,jj)
+                        &      - ( gdept(ji,jj+1,jk-kp,Kmm) - gdept(ji,jj,jk-kp,Kmm) ) / e2v(ji,jj)  ) * vmask(ji,jj,jk)
+                     ze3_e2    =  e3w(ji,jj+jp,jk-kp,Kmm) / e2v(ji,jj)
                      ztj_mlb(ji   ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2  , ABS( ztj_g_raw ) ), ztj_g_raw )
                   ENDIF
@@ -542,11 +544,11 @@
                      !
                      ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface
-                     zti_coord = znot_thru_surface * ( gdept_n(ji+1,jj  ,jk) - gdept_n(ji,jj,jk) ) * r1_e1u(ji,jj)
-                     ztj_coord = znot_thru_surface * ( gdept_n(ji  ,jj+1,jk) - gdept_n(ji,jj,jk) ) * r1_e2v(ji,jj)     ! unmasked
+                     zti_coord = znot_thru_surface * ( gdept(ji+1,jj  ,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj)
+                     ztj_coord = znot_thru_surface * ( gdept(ji  ,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj)     ! unmasked
                      zti_g_raw = zti_raw - zti_coord      ! ref to geopot surfaces
                      ztj_g_raw = ztj_raw - ztj_coord
                      ! additional limit required in bilaplacian case
-                     ze3_e1    = e3w_n(ji+ip,jj   ,jk+kp) * r1_e1u(ji,jj)
-                     ze3_e2    = e3w_n(ji   ,jj+jp,jk+kp) * r1_e2v(ji,jj)
+                     ze3_e1    = e3w(ji+ip,jj   ,jk+kp,Kmm) * r1_e1u(ji,jj)
+                     ze3_e2    = e3w(ji   ,jj+jp,jk+kp,Kmm) * r1_e2v(ji,jj)
                      ! NB: hard coded factor 5 (can be a namelist parameter...)
                      zti_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS( zti_g_raw ) ), zti_g_raw )
@@ -561,8 +563,8 @@
                      zti_g_lim =          ( zfacti   * zti_g_lim                       &
                         &      + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp)   &
-                        &                           * gdepw_n(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp)
+                        &                           * gdepw(ji+ip,jj,jk+kp,Kmm) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp)
                      ztj_g_lim =          ( zfactj   * ztj_g_lim                       &
                         &      + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp)   &
-                        &                           * gdepw_n(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp)
+                        &                           * gdepw(ji,jj+jp,jk+kp,Kmm) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp)
                      !
                      triadi_g(ji+ip,jj   ,jk,1-ip,kp) = zti_g_lim
@@ -596,8 +598,8 @@
                      triadj(ji   ,jj+jp,jk,1-jp,kp) = ztj_lim * zjsw
                      !
-                     zbu  = e1e2u(ji   ,jj   ) * e3u_n(ji   ,jj   ,jk   )
-                     zbv  = e1e2v(ji   ,jj   ) * e3v_n(ji   ,jj   ,jk   )
-                     zbti = e1e2t(ji+ip,jj   ) * e3w_n(ji+ip,jj   ,jk+kp)
-                     zbtj = e1e2t(ji   ,jj+jp) * e3w_n(ji   ,jj+jp,jk+kp)
+                     zbu  = e1e2u(ji   ,jj   ) * e3u(ji   ,jj   ,jk   ,Kmm)
+                     zbv  = e1e2v(ji   ,jj   ) * e3v(ji   ,jj   ,jk   ,Kmm)
+                     zbti = e1e2t(ji+ip,jj   ) * e3w(ji+ip,jj   ,jk+kp,Kmm)
+                     zbtj = e1e2t(ji   ,jj+jp) * e3w(ji   ,jj+jp,jk+kp,Kmm)
                      !
                      wslp2(ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim*zti_g_lim      ! masked
@@ -618,5 +620,5 @@
 
 
-   SUBROUTINE ldf_slp_mxl( prd, pn2, p_gru, p_grv, p_dzr )
+   SUBROUTINE ldf_slp_mxl( prd, pn2, p_gru, p_grv, p_dzr, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ldf_slp_mxl  ***
@@ -638,4 +640,5 @@
       REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   p_gru, p_grv   ! i- & j-gradient of density (u- & v-pts)
       REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   p_dzr          ! z-gradient of density      (T-point)
+      INTEGER , INTENT(in)                   ::   Kmm            ! ocean time level indices
       !!
       INTEGER  ::   ji , jj , jk                   ! dummy loop indices
@@ -694,6 +697,6 @@
             !                        !- bound the slopes: abs(zw.)<= 1/100 and zb..<0
             !                           kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
-            zbu = MIN(  zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u_n(ji,jj,iku)* ABS( zau )  )
-            zbv = MIN(  zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v_n(ji,jj,ikv)* ABS( zav )  )
+            zbu = MIN(  zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u(ji,jj,iku,Kmm)* ABS( zau )  )
+            zbv = MIN(  zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v(ji,jj,ikv,Kmm)* ABS( zav )  )
             !                        !- Slope at u- & v-points (uslpml, vslpml)
             uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku)
@@ -717,6 +720,6 @@
             !                        !- bound the slopes: abs(zw.)<= 1/100 and zb..<0.
             !                           kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
-            zbi = MIN(  zbw , -100._wp* ABS( zai ) , -7.e+3_wp/e3w_n(ji,jj,ik)* ABS( zai )  )
-            zbj = MIN(  zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w_n(ji,jj,ik)* ABS( zaj )  )
+            zbi = MIN(  zbw , -100._wp* ABS( zai ) , -7.e+3_wp/e3w(ji,jj,ik,Kmm)* ABS( zai )  )
+            zbj = MIN(  zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w(ji,jj,ik,Kmm)* ABS( zaj )  )
             !                        !- i- & j-slope at w-points (wslpiml, wslpjml)
             wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik)
@@ -786,8 +789,8 @@
 !               DO jj = 2, jpjm1
 !                  DO ji = fs_2, fs_jpim1   ! vector opt.
-!                     uslp (ji,jj,jk) = - ( gdept_n(ji+1,jj,jk) - gdept_n(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
-!                     vslp (ji,jj,jk) = - ( gdept_n(ji,jj+1,jk) - gdept_n(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
-!                     wslpi(ji,jj,jk) = - ( gdepw_n(ji+1,jj,jk) - gdepw_n(ji-1,jj,jk) ) * r1_e1t(ji,jj) * wmask(ji,jj,jk) * 0.5
-!                     wslpj(ji,jj,jk) = - ( gdepw_n(ji,jj+1,jk) - gdepw_n(ji,jj-1,jk) ) * r1_e2t(ji,jj) * wmask(ji,jj,jk) * 0.5
+!                     uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kmm) - gdept(ji ,jj ,jk,Kmm) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
+!                     vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kmm) - gdept(ji ,jj ,jk,Kmm) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
+!                     wslpi(ji,jj,jk) = - ( gdepw(ji+1,jj,jk,Kmm) - gdepw(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj) * wmask(ji,jj,jk) * 0.5
+!                     wslpj(ji,jj,jk) = - ( gdepw(ji,jj+1,jk,Kmm) - gdepw(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj) * wmask(ji,jj,jk) * 0.5
 !                  END DO
 !               END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/LDF/ldftra.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/LDF/ldftra.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/LDF/ldftra.F90	(revision 10922)
@@ -382,5 +382,5 @@
 
 
-   SUBROUTINE ldf_tra( kt )
+   SUBROUTINE ldf_tra( kt, Kbb, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ldf_tra  ***
@@ -403,4 +403,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! time step
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -411,5 +412,5 @@
          !                                ! =F(growth rate of baroclinic instability)
          !                                ! max value aeiv_0 ; decreased to 0 within 20N-20S
-         CALL ldf_eiv( kt, aei0, aeiu, aeiv )
+         CALL ldf_eiv( kt, aei0, aeiu, aeiv, Kmm )
       ENDIF
       !
@@ -424,5 +425,5 @@
             ahtv(:,:,1) = aeiv(:,:,1)
          ELSE                                            ! compute aht. 
-            CALL ldf_eiv( kt, aht0, ahtu, ahtv )
+            CALL ldf_eiv( kt, aht0, ahtu, ahtv, Kmm )
          ENDIF
          !
@@ -448,11 +449,11 @@
          IF( ln_traldf_lap     ) THEN          !   laplacian operator |u| e /12
             DO jk = 1, jpkm1
-               ahtu(:,:,jk) = ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12   ! n.b. ub,vb are masked
-               ahtv(:,:,jk) = ABS( vb(:,:,jk) ) * e2v(:,:) * r1_12
+               ahtu(:,:,jk) = ABS( uu(:,:,jk,Kbb) ) * e1u(:,:) * r1_12   ! n.b. uu,vv are masked
+               ahtv(:,:,jk) = ABS( vv(:,:,jk,Kbb) ) * e2v(:,:) * r1_12
             END DO
          ELSEIF( ln_traldf_blp ) THEN      ! bilaplacian operator      sqrt( |u| e^3 /12 ) = sqrt( |u| e /12 ) * e
             DO jk = 1, jpkm1
-               ahtu(:,:,jk) = SQRT(  ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12  ) * e1u(:,:)
-               ahtv(:,:,jk) = SQRT(  ABS( vb(:,:,jk) ) * e2v(:,:) * r1_12  ) * e2v(:,:)
+               ahtu(:,:,jk) = SQRT(  ABS( uu(:,:,jk,Kbb) ) * e1u(:,:) * r1_12  ) * e1u(:,:)
+               ahtv(:,:,jk) = SQRT(  ABS( vv(:,:,jk,Kbb) ) * e2v(:,:) * r1_12  ) * e2v(:,:)
             END DO
          ENDIF
@@ -625,5 +626,5 @@
 
 
-   SUBROUTINE ldf_eiv( kt, paei0, paeiu, paeiv )
+   SUBROUTINE ldf_eiv( kt, paei0, paeiu, paeiv, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ldf_eiv  ***
@@ -637,4 +638,5 @@
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step index
+      INTEGER                         , INTENT(in   ) ::   Kmm            ! ocean time level indices
       REAL(wp)                        , INTENT(inout) ::   paei0          ! max value            [m2/s]
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   paeiu, paeiv   ! eiv coefficient      [m2/s]
@@ -658,9 +660,9 @@
                   ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 
                   zn2 = MAX( rn2b(ji,jj,jk), 0._wp )
-                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk)
+                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w(ji,jj,jk,Kmm)
                   ! Compute elements required for the inverse time scale of baroclinic
                   ! eddies using the isopycnal slopes calculated in ldfslp.F : 
                   ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
-                  ze3w = e3w_n(ji,jj,jk) * tmask(ji,jj,jk)
+                  ze3w = e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
                   zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w
                   zhw(ji,jj) = zhw(ji,jj) + ze3w
@@ -676,9 +678,9 @@
                   ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 
                   zn2 = MAX( rn2b(ji,jj,jk), 0._wp )
-                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk)
+                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w(ji,jj,jk,Kmm)
                   ! Compute elements required for the inverse time scale of baroclinic
                   ! eddies using the isopycnal slopes calculated in ldfslp.F : 
                   ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
-                  ze3w = e3w_n(ji,jj,jk) * tmask(ji,jj,jk)
+                  ze3w = e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
                   zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
                      &                            + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w
@@ -725,5 +727,5 @@
 
 
-   SUBROUTINE ldf_eiv_trp( kt, kit000, pun, pvn, pwn, cdtype )
+   SUBROUTINE ldf_eiv_trp( kt, kit000, pu, pv, pw, cdtype, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ldf_eiv_trp  ***
@@ -741,12 +743,13 @@
       !!                                    velocity and heat transport (call ldf_eiv_dia)
       !!
-      !! ** Action  : pun, pvn increased by the eiv transport
+      !! ** Action  : pu, pv increased by the eiv transport
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   kt       ! ocean time-step index
       INTEGER                         , INTENT(in   ) ::   kit000   ! first time step index
+      INTEGER                         , INTENT(in   ) ::   Kmm      ! ocean time level indices
       CHARACTER(len=3)                , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun      ! in : 3 ocean transport components   [m3/s]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pvn      ! out: 3 ocean transport components   [m3/s]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pwn      ! increased by the eiv                [m3/s]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu      ! in : 3 ocean transport components   [m3/s]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pv      ! out: 3 ocean transport components   [m3/s]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pw      ! increased by the eiv                [m3/s]
       !!
       INTEGER  ::   ji, jj, jk                 ! dummy loop indices
@@ -780,6 +783,6 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.               
-               pun(ji,jj,jk) = pun(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
-               pvn(ji,jj,jk) = pvn(ji,jj,jk) - ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
+               pu(ji,jj,jk) = pu(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
+               pv(ji,jj,jk) = pv(ji,jj,jk) - ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
             END DO
          END DO
@@ -788,5 +791,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               pwn(ji,jj,jk) = pwn(ji,jj,jk) + (  zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj  ,jk)   &
+               pw(ji,jj,jk) = pw(ji,jj,jk) + (  zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj  ,jk)   &
                   &                             + zpsi_vw(ji,jj,jk) - zpsi_vw(ji  ,jj-1,jk) )
             END DO
@@ -795,10 +798,10 @@
       !
       !                              ! diagnose the eddy induced velocity and associated heat transport
-      IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' )   CALL ldf_eiv_dia( zpsi_uw, zpsi_vw )
+      IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' )   CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm )
       !
     END SUBROUTINE ldf_eiv_trp
 
 
-   SUBROUTINE ldf_eiv_dia( psi_uw, psi_vw )
+   SUBROUTINE ldf_eiv_dia( psi_uw, psi_vw, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ldf_eiv_dia  ***
@@ -811,4 +814,5 @@
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   psi_uw, psi_vw   ! streamfunction   [m3/s]
+      INTEGER                         , INTENT(in   ) ::   Kmm   ! ocean time level indices
       !
       INTEGER  ::   ji, jj, jk    ! dummy loop indices
@@ -831,10 +835,10 @@
       !
       DO jk = 1, jpkm1                                         ! e2u e3u u_eiv = -dk[psi_uw]
-         zw3d(:,:,jk) = ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) / ( e2u(:,:) * e3u_n(:,:,jk) )
+         zw3d(:,:,jk) = ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) / ( e2u(:,:) * e3u(:,:,jk,Kmm) )
       END DO
       CALL iom_put( "uoce_eiv", zw3d )
       !
       DO jk = 1, jpkm1                                         ! e1v e3v v_eiv = -dk[psi_vw]
-         zw3d(:,:,jk) = ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) / ( e1v(:,:) * e3v_n(:,:,jk) )
+         zw3d(:,:,jk) = ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) / ( e1v(:,:) * e3v(:,:,jk,Kmm) )
       END DO
       CALL iom_put( "voce_eiv", zw3d )
@@ -859,6 +863,6 @@
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
-                 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   &
-                    &                            * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji+1,jj,jk,jp_tem) ) 
+                 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1)          - psi_uw(ji  ,jj,jk)            )   &
+                    &                            * ( ts    (ji,jj,jk,jp_tem,Kmm) + ts    (ji+1,jj,jk,jp_tem,Kmm) ) 
                  zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
               END DO
@@ -875,6 +879,6 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   &
-                  &                            * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji,jj+1,jk,jp_tem) ) 
+               zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)          - psi_vw(ji,jj  ,jk)            )   &
+                  &                            * ( ts    (ji,jj,jk,jp_tem,Kmm) + ts    (ji,jj+1,jk,jp_tem,Kmm) ) 
                zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
             END DO
@@ -894,6 +898,6 @@
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
-                 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   &
-                    &                            * ( tsn   (ji,jj,jk,jp_sal) + tsn   (ji+1,jj,jk,jp_sal) ) 
+                 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1)          - psi_uw(ji  ,jj,jk)            )   &
+                    &                            * ( ts    (ji,jj,jk,jp_sal,Kmm) + ts    (ji+1,jj,jk,jp_sal,Kmm) ) 
                  zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
               END DO
@@ -903,5 +907,5 @@
         CALL lbc_lnk( 'ldftra', zw3d, 'U', -1. )
         CALL iom_put( "ueiv_salttr", zztmp * zw2d )                  ! salt transport in i-direction
-        CALL iom_put( "ueiv_salttr3d", zztmp * zw3d )                  ! salt transport in i-direction
+        CALL iom_put( "ueiv_salttr3d", zztmp * zw3d )                ! salt transport in i-direction
       ENDIF
       zw2d(:,:) = 0._wp 
@@ -910,6 +914,6 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   &
-                  &                            * ( tsn   (ji,jj,jk,jp_sal) + tsn   (ji,jj+1,jk,jp_sal) ) 
+               zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)          - psi_vw(ji,jj  ,jk)            )   &
+                  &                            * ( ts    (ji,jj,jk,jp_sal,Kmm) + ts    (ji,jj+1,jk,jp_sal,Kmm) ) 
                zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
             END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/OBS/diaobs.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/OBS/diaobs.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/OBS/diaobs.F90	(revision 10922)
@@ -103,5 +103,5 @@
 CONTAINS
 
-   SUBROUTINE dia_obs_init
+   SUBROUTINE dia_obs_init( Kmm )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE dia_obs_init  ***
@@ -114,5 +114,6 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER, PARAMETER ::   jpmaxnfiles = 1000    ! Maximum number of files for each obs type
+      INTEGER, INTENT(in)                ::   Kmm                      ! ocean time level indices
+      INTEGER, PARAMETER                 ::   jpmaxnfiles = 1000       ! Maximum number of files for each obs type
       INTEGER, DIMENSION(:), ALLOCATABLE ::   ifilesprof, ifilessurf   ! Number of profile & surface files
       INTEGER :: ios             ! Local integer output status for namelist read
@@ -429,5 +430,5 @@
                &               jpi, jpj, jpk, &
                &               zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2,  &
-               &               ln_nea, ln_bound_reject, &
+               &               ln_nea, ln_bound_reject, Kmm, &
                &               kdailyavtypes = nn_profdavtypes )
          END DO
@@ -483,5 +484,5 @@
 
 
-   SUBROUTINE dia_obs( kstp )
+   SUBROUTINE dia_obs( kstp, Kmm )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE dia_obs  ***
@@ -496,7 +497,7 @@
       !! ** Action  :
       !!----------------------------------------------------------------------
-      USE dom_oce, ONLY : gdept_n, gdept_1d   ! Ocean space and time domain variables
+      USE dom_oce, ONLY : gdept, gdept_1d     ! Ocean space domain variables (Kmm time-level only)
       USE phycst , ONLY : rday                ! Physical constants
-      USE oce    , ONLY : tsn, un, vn, sshn   ! Ocean dynamics and tracers variables
+      USE oce    , ONLY : ts, uu, vv, sshn    ! Ocean dynamics and tracers variables (Kmm time-level only)
       USE phycst , ONLY : rday                ! Physical constants
 #if defined  key_si3
@@ -511,4 +512,5 @@
       !! * Arguments
       INTEGER, INTENT(IN) :: kstp  ! Current timestep
+      INTEGER, INTENT(in) :: Kmm   ! ocean time level indices
       !! * Local declarations
       INTEGER :: idaystp           ! Number of timesteps per day
@@ -551,6 +553,6 @@
             SELECT CASE ( TRIM(cobstypesprof(jtype)) )
             CASE('prof')
-               zprofvar1(:,:,:) = tsn(:,:,:,jp_tem)
-               zprofvar2(:,:,:) = tsn(:,:,:,jp_sal)
+               zprofvar1(:,:,:) = ts(:,:,:,jp_tem,Kmm)
+               zprofvar2(:,:,:) = ts(:,:,:,jp_sal,Kmm)
                zprofmask1(:,:,:) = tmask(:,:,:)
                zprofmask2(:,:,:) = tmask(:,:,:)
@@ -560,6 +562,6 @@
                zgphi2(:,:) = gphit(:,:)
             CASE('vel')
-               zprofvar1(:,:,:) = un(:,:,:)
-               zprofvar2(:,:,:) = vn(:,:,:)
+               zprofvar1(:,:,:) = uu(:,:,:,Kmm)
+               zprofvar2(:,:,:) = vv(:,:,:,Kmm)
                zprofmask1(:,:,:) = umask(:,:,:)
                zprofmask2(:,:,:) = vmask(:,:,:)
@@ -575,5 +577,5 @@
                &               nit000, idaystp,                         &
                &               zprofvar1, zprofvar2,                    &
-               &               gdept_n(:,:,:), gdepw_n(:,:,:),            & 
+               &               gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm),      & 
                &               zprofmask1, zprofmask2,                  &
                &               zglam1, zglam2, zgphi1, zgphi2,          &
@@ -594,9 +596,9 @@
             SELECT CASE ( TRIM(cobstypessurf(jtype)) )
             CASE('sst')
-               zsurfvar(:,:) = tsn(:,:,1,jp_tem)
+               zsurfvar(:,:) = ts(:,:,1,jp_tem,Kmm)
             CASE('sla')
                zsurfvar(:,:) = sshn(:,:)
             CASE('sss')
-               zsurfvar(:,:) = tsn(:,:,1,jp_sal)
+               zsurfvar(:,:) = ts(:,:,1,jp_sal,Kmm)
             CASE('sic')
                IF ( kstp == 0 ) THEN
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/OBS/obs_prep.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/OBS/obs_prep.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/OBS/obs_prep.F90	(revision 10922)
@@ -244,5 +244,5 @@
       &                     kpi, kpj, kpk, &
       &                     zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2,  &
-      &                     ld_nea, ld_bound_reject, kdailyavtypes,  kqc_cutoff )
+      &                     ld_nea, ld_bound_reject, Kmm, kdailyavtypes,  kqc_cutoff )
 
 !!----------------------------------------------------------------------
@@ -274,4 +274,5 @@
       LOGICAL, INTENT(IN) :: ld_bound_reject      ! Switch for rejecting observations near the boundary
       INTEGER, INTENT(IN) :: kpi, kpj, kpk        ! Local domain sizes
+      INTEGER, INTENT(IN) :: Kmm                  ! time-level index
       INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
          & kdailyavtypes                          ! Types for daily averages
@@ -420,5 +421,5 @@
          &                 inlav1obs,             ld_nea,               &
          &                 ibdyv1obs,             ld_bound_reject,      &
-         &                 iqc_cutoff       )
+         &                 iqc_cutoff,            Kmm                 )
 
       CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp )
@@ -442,5 +443,5 @@
          &                 inlav2obs,             ld_nea,               &
          &                 ibdyv2obs,             ld_bound_reject,      &
-         &                 iqc_cutoff       )
+         &                 iqc_cutoff,            Kmm                 )
 
       CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp )
@@ -1094,5 +1095,5 @@
       &                       klanobs, knlaobs, ld_nea,         &
       &                       kbdyobs, ld_bound_reject,         &
-      &                       kqc_cutoff                        )
+      &                       kqc_cutoff,       Kmm             )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_coo_spc_3d  ***
@@ -1116,9 +1117,9 @@
       !!----------------------------------------------------------------------
       !! * Modules used
-      USE dom_oce, ONLY : &       ! Geographical information
+      USE dom_oce, ONLY : &       ! Geographical information 
          & gdepw_1d,      &
          & gdepw_0,       &                       
-         & gdepw_n,       &
-         & gdept_n,       &
+         & gdepw,         &
+         & gdept,         &
          & ln_zco,        &
          & ln_zps             
@@ -1160,4 +1161,5 @@
       LOGICAL, INTENT(IN) :: ld_bound_reject  ! Flag observations near open boundary
       INTEGER, INTENT(IN) :: kqc_cutoff     ! Cutoff QC value
+      INTEGER, INTENT(IN) :: Kmm            ! time-level index
 
       !! * Local declarations
@@ -1230,5 +1232,5 @@
       CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam )
       CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
-      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), &
+      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw(:,:,:,Kmm), &
         &                     zgdepw )
 
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/fldread.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/fldread.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/fldread.F90	(revision 10922)
@@ -46,4 +46,5 @@
    PUBLIC   fld_clopn
 
+   INTEGER :: nfld_Nnn = 1
    TYPE, PUBLIC ::   FLD_N      !: Namelist field informations
       CHARACTER(len = 256) ::   clname      ! generic name of the NetCDF flux file
@@ -902,5 +903,5 @@
                      WRITE(ibstr,"(I10.10)") map%ptr(ib) 
                      CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
-                 !   IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1),  ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj
+                 !   IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t(zij,zjj,:,nfld_Nnn), mask=tmask(zij,zjj,:)==1),  ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj
                   ENDIF
                CASE(2)
@@ -908,5 +909,5 @@
                      WRITE(ibstr,"(I10.10)") map%ptr(ib) 
                      CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
-                     IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u_n(zij,zjj,:), mask=umask(zij,zjj,:)==1),  sum(umask(zij,zjj,:)), &
+                     IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u(zij,zjj,:,nfld_Nnn), mask=umask(zij,zjj,:)==1),  sum(umask(zij,zjj,:)), &
                        &                hu_n(zij,zjj), map%ptr(ib), ib, zij, zjj, narea-1  , &
                         &                dta_read(map%ptr(ib),1,:)
@@ -921,16 +922,16 @@
                SELECT CASE( igrd )                       
                   CASE(1)
-                     zl =  gdept_n(zij,zjj,ik)                                          ! if using in step could use fsdept instead of gdept_n?
+                     zl =  gdept(zij,zjj,ik,nfld_Nnn)                                          ! if using in step could use fsdept instead of gdept_n?
                   CASE(2)
                      IF(ln_sco) THEN
-                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
+                        zl =  ( gdept(zij,zjj,ik,nfld_Nnn) + gdept(zij+1,zjj,ik,nfld_Nnn) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
                      ELSE
-                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) ) 
+                        zl =  MIN( gdept(zij,zjj,ik,nfld_Nnn), gdept(zij+1,zjj,ik,nfld_Nnn) ) 
                      ENDIF
                   CASE(3)
                      IF(ln_sco) THEN
-                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
+                        zl =  ( gdept(zij,zjj,ik,nfld_Nnn) + gdept(zij,zjj+1,ik,nfld_Nnn) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
                      ELSE
-                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) )
+                        zl =  MIN( gdept(zij,zjj,ik,nfld_Nnn), gdept(zij,zjj+1,ik,nfld_Nnn) )
                      ENDIF
                END SELECT
@@ -940,5 +941,5 @@
                   dta(ib,1,ik) =  dta_read(map%ptr(ib),1,MAXLOC(dta_read_z(map%ptr(ib),1,:),1))
                ELSE                                                                                ! inbetween : vertical interpolation between ikk & ikk+1
-                  DO ikk = 1, jpkm1_bdy                                                            ! when  gdept_n(ikk) < zl < gdept_n(ikk+1)
+                  DO ikk = 1, jpkm1_bdy                                                            ! when  gdept(ikk,nfld_Nnn) < zl < gdept(ikk+1,nfld_Nnn)
                      IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp) &
                     &    .AND. (dta_read_z(map%ptr(ib),1,ikk+1) /= fv_alt)) THEN
@@ -964,5 +965,5 @@
               ENDDO
               DO ik = 1, ipk                                ! calculate transport on model grid
-                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik)
+                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3u(zij,zjj,ik,nfld_Nnn) * umask(zij,zjj,ik)
               ENDDO
               DO ik = 1, ipk                                ! make transport correction
@@ -989,5 +990,5 @@
               ENDDO
               DO ik = 1, ipk                                ! calculate transport on model grid
-                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik)
+                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3v(zij,zjj,ik,nfld_Nnn) * vmask(zij,zjj,ik)
               ENDDO
               DO ik = 1, ipk                                ! make transport correction
@@ -1027,5 +1028,5 @@
                      WRITE(ibstr,"(I10.10)") map%ptr(ib) 
                      CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
-                 !   IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1),  ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj
+                 !   IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t(zij,zjj,:,nfld_Nnn), mask=tmask(zij,zjj,:)==1),  ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj
                   ENDIF
                CASE(2)
@@ -1043,16 +1044,16 @@
                SELECT CASE( igrd )                          ! coded for sco - need zco and zps option using min
                   CASE(1)
-                     zl =  gdept_n(zij,zjj,ik)              ! if using in step could use fsdept instead of gdept_n?
+                     zl =  gdept(zij,zjj,ik,nfld_Nnn)              ! if using in step could use fsdept instead of gdept_n?
                   CASE(2)
                      IF(ln_sco) THEN
-                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
+                        zl =  ( gdept(zij,zjj,ik,nfld_Nnn) + gdept(zij+1,zjj,ik,nfld_Nnn) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
                      ELSE
-                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) )
+                        zl =  MIN( gdept(zij,zjj,ik,nfld_Nnn), gdept(zij+1,zjj,ik,nfld_Nnn) )
                      ENDIF
                   CASE(3)
                      IF(ln_sco) THEN
-                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
+                        zl =  ( gdept(zij,zjj,ik,nfld_Nnn) + gdept(zij,zjj+1,ik,nfld_Nnn) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
                      ELSE
-                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) )
+                        zl =  MIN( gdept(zij,zjj,ik,nfld_Nnn), gdept(zij,zjj+1,ik,nfld_Nnn) )
                      ENDIF
                END SELECT
@@ -1062,5 +1063,5 @@
                   dta(ib,1,ik) =  dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1))
                ELSE                                                                     ! inbetween : vertical interpolation between ikk & ikk+1
-                  DO ikk = 1, jpkm1_bdy                                                 ! when  gdept_n(ikk) < zl < gdept_n(ikk+1)
+                  DO ikk = 1, jpkm1_bdy                                                 ! when  gdept(ikk,nfld_Nnn) < zl < gdept(ikk+1,nfld_Nnn)
                      IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp) &
                     &    .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN
@@ -1088,5 +1089,5 @@
                ENDDO
                DO ik = 1, ipk                                ! calculate transport on model grid
-                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik)
+                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3u(zij,zjj,ik,nfld_Nnn) * umask(zij,zjj,ik)
                ENDDO
                DO ik = 1, ipk                                ! make transport correction
@@ -1113,5 +1114,5 @@
                ENDDO
                DO ik = 1, ipk                                ! calculate transport on model grid
-                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik)
+                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3v(zij,zjj,ik,nfld_Nnn) * vmask(zij,zjj,ik)
                ENDDO
                DO ik = 1, ipk                                ! make transport correction
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbccpl.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbccpl.F90	(revision 10922)
@@ -32,5 +32,5 @@
    USE cpl_oasis3     ! OASIS3 coupling
    USE geo2ocean      ! 
-   USE oce     , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev
+   USE oce     , ONLY : ts, uu, vv, sshn, sshb, fraqsr_1lev
    USE ocealb         ! 
    USE eosbn2         ! 
@@ -1049,5 +1049,5 @@
 
 
-   SUBROUTINE sbc_cpl_rcv( kt, k_fsbc, k_ice )     
+   SUBROUTINE sbc_cpl_rcv( kt, k_fsbc, k_ice, Kbb, Kmm )     
       !!----------------------------------------------------------------------
       !!             ***  ROUTINE sbc_cpl_rcv  ***
@@ -1099,4 +1099,5 @@
       INTEGER, INTENT(in) ::   k_fsbc      ! frequency of sbc (-> ice model) computation 
       INTEGER, INTENT(in) ::   k_ice       ! ice management in the sbc (=0/1/2/3)
+      INTEGER, INTENT(in) ::   Kbb, Kmm    ! ocean model time level indices
       !!
       LOGICAL  ::   llnewtx, llnewtau      ! update wind stress components and module??
@@ -1302,5 +1303,5 @@
          IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction &
                                       .OR. srcv(jpr_hsig)%laction   .OR. srcv(jpr_wfreq)%laction) THEN
-            CALL sbc_stokes()
+            CALL sbc_stokes( Kmm )
          ENDIF
       ENDIF
@@ -1354,12 +1355,12 @@
       IF( srcv(jpr_ocx1)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
          ssu_m(:,:) = frcv(jpr_ocx1)%z3(:,:,1)
-         ub (:,:,1) = ssu_m(:,:)                             ! will be used in icestp in the call of ice_forcing_tau
-         un (:,:,1) = ssu_m(:,:)                             ! will be used in sbc_cpl_snd if atmosphere coupling
+         uu(:,:,1,Kbb) = ssu_m(:,:)                          ! will be used in icestp in the call of ice_forcing_tau
+         uu(:,:,1,Kmm) = ssu_m(:,:)                          ! will be used in sbc_cpl_snd if atmosphere coupling
          CALL iom_put( 'ssu_m', ssu_m )
       ENDIF
       IF( srcv(jpr_ocy1)%laction ) THEN
          ssv_m(:,:) = frcv(jpr_ocy1)%z3(:,:,1)
-         vb (:,:,1) = ssv_m(:,:)                             ! will be used in icestp in the call of ice_forcing_tau
-         vn (:,:,1) = ssv_m(:,:)                             ! will be used in sbc_cpl_snd if atmosphere coupling
+         vv(:,:,1,Kbb) = ssv_m(:,:)                          ! will be used in icestp in the call of ice_forcing_tau
+         vv(:,:,1,Kmm) = ssv_m(:,:)                          ! will be used in sbc_cpl_snd if atmosphere coupling
          CALL iom_put( 'ssv_m', ssv_m )
       ENDIF
@@ -2036,5 +2037,5 @@
    
    
-   SUBROUTINE sbc_cpl_snd( kt )
+   SUBROUTINE sbc_cpl_snd( kt, Kmm )
       !!----------------------------------------------------------------------
       !!             ***  ROUTINE sbc_cpl_snd  ***
@@ -2046,4 +2047,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt
+      INTEGER, INTENT(in) ::   Kmm    ! ocean model time level index
       !
       INTEGER ::   ji, jj, jl   ! dummy loop indices
@@ -2063,9 +2065,9 @@
          
          IF ( nn_components == jp_iam_opa ) THEN
-            ztmp1(:,:) = tsn(:,:,1,jp_tem)   ! send temperature as it is (potential or conservative) -> use of l_useCT on the received part
+            ztmp1(:,:) = ts(:,:,1,jp_tem,Kmm)   ! send temperature as it is (potential or conservative) -> use of l_useCT on the received part
          ELSE
             ! we must send the surface potential temperature 
-            IF( l_useCT )  THEN    ;   ztmp1(:,:) = eos_pt_from_ct( tsn(:,:,1,jp_tem), tsn(:,:,1,jp_sal) )
-            ELSE                   ;   ztmp1(:,:) = tsn(:,:,1,jp_tem)
+            IF( l_useCT )  THEN    ;   ztmp1(:,:) = eos_pt_from_ct( ts(:,:,1,jp_tem,Kmm), ts(:,:,1,jp_sal,Kmm) )
+            ELSE                   ;   ztmp1(:,:) = ts(:,:,1,jp_tem,Kmm)
             ENDIF
             !
@@ -2095,5 +2097,5 @@
                CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
                END SELECT
-            CASE( 'oce and weighted ice')    ;   ztmp1(:,:) =   tsn(:,:,1,jp_tem) + rt0  
+            CASE( 'oce and weighted ice')    ;   ztmp1(:,:) =   ts(:,:,1,jp_tem,Kmm) + rt0  
                SELECT CASE( sn_snd_temp%clcat ) 
                CASE( 'yes' )    
@@ -2316,6 +2318,6 @@
          !                                                               i      i+1 (for I)
          IF( nn_components == jp_iam_opa ) THEN
-            zotx1(:,:) = un(:,:,1)  
-            zoty1(:,:) = vn(:,:,1)  
+            zotx1(:,:) = uu(:,:,1,Kmm)  
+            zoty1(:,:) = vv(:,:,1,Kmm)  
          ELSE        
             SELECT CASE( TRIM( sn_snd_crt%cldes ) )
@@ -2323,6 +2325,6 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )
-                     zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji  ,jj-1,1) ) 
+                     zotx1(ji,jj) = 0.5 * ( uu(ji,jj,1,Kmm) + uu(ji-1,jj  ,1,Kmm) )
+                     zoty1(ji,jj) = 0.5 * ( vv(ji,jj,1,Kmm) + vv(ji  ,jj-1,1,Kmm) ) 
                   END DO
                END DO
@@ -2330,8 +2332,8 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)  
-                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)
-                     zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
-                     zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
+                     zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Kmm) + uu   (ji-1,jj  ,1,Kmm) ) * zfr_l(ji,jj)  
+                     zoty1(ji,jj) = 0.5 * ( vv   (ji,jj,1,Kmm) + vv   (ji  ,jj-1,1,Kmm) ) * zfr_l(ji,jj)
+                     zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  )     + u_ice(ji-1,jj    )     ) *  fr_i(ji,jj)
+                     zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  )     + v_ice(ji  ,jj-1  )     ) *  fr_i(ji,jj)
                   END DO
                END DO
@@ -2340,8 +2342,8 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   &
-                        &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
-                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   &
-                        &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
+                     zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Kmm) + uu   (ji-1,jj  ,1,Kmm) ) * zfr_l(ji,jj)   &
+                        &         + 0.5 * ( u_ice(ji,jj  )     + u_ice(ji-1,jj    )     ) *  fr_i(ji,jj)
+                     zoty1(ji,jj) = 0.5 * ( vv   (ji,jj,1,Kmm) + vv   (ji  ,jj-1,1,Kmm) ) * zfr_l(ji,jj)   &
+                        &         + 0.5 * ( v_ice(ji,jj  )     + v_ice(ji  ,jj-1  )     ) *  fr_i(ji,jj)
                   END DO
                END DO
@@ -2406,6 +2408,6 @@
              DO jj = 2, jpjm1 
                 DO ji = fs_2, fs_jpim1   ! vector opt. 
-                   zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) ) 
-                   zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji , jj-1,1) )  
+                   zotx1(ji,jj) = 0.5 * ( uu(ji,jj,1,Kmm) + uu(ji-1,jj  ,1,Kmm) ) 
+                   zoty1(ji,jj) = 0.5 * ( vv(ji,jj,1,Kmm) + vv(ji , jj-1,1,Kmm) )  
                 END DO 
              END DO 
@@ -2413,6 +2415,6 @@
              DO jj = 2, jpjm1 
                 DO ji = fs_2, fs_jpim1   ! vector opt. 
-                   zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   
-                   zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj) 
+                   zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Kmm) + uu   (ji-1,jj  ,1,Kmm) ) * zfr_l(ji,jj)   
+                   zoty1(ji,jj) = 0.5 * ( vv   (ji,jj,1,Kmm) + vv   (ji  ,jj-1,1,Kmm) ) * zfr_l(ji,jj) 
                    zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj) 
                    zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj) 
@@ -2423,7 +2425,7 @@
              DO jj = 2, jpjm1 
                 DO ji = fs_2, fs_jpim1   ! vector opt. 
-                   zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   & 
+                   zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Kmm) + uu   (ji-1,jj  ,1,Kmm) ) * zfr_l(ji,jj)   & 
                       &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj) 
-                   zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
+                   zoty1(ji,jj) = 0.5 * ( vv   (ji,jj,1,Kmm) + vv   (ji  ,jj-1,1,Kmm) ) * zfr_l(ji,jj)   & 
                       &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj) 
                 END DO
@@ -2497,9 +2499,9 @@
       !                                                        ! SSS
       IF( ssnd(jps_soce  )%laction )  THEN
-         CALL cpl_snd( jps_soce  , isec, RESHAPE ( tsn(:,:,1,jp_sal), (/jpi,jpj,1/) ), info )
+         CALL cpl_snd( jps_soce  , isec, RESHAPE ( ts(:,:,1,jp_sal,Kmm), (/jpi,jpj,1/) ), info )
       ENDIF
       !                                                        ! first T level thickness 
       IF( ssnd(jps_e3t1st )%laction )  THEN
-         CALL cpl_snd( jps_e3t1st, isec, RESHAPE ( e3t_n(:,:,1)   , (/jpi,jpj,1/) ), info )
+         CALL cpl_snd( jps_e3t1st, isec, RESHAPE ( e3t(:,:,1,Kmm)   , (/jpi,jpj,1/) ), info )
       ENDIF
       !                                                        ! Qsr fraction
@@ -2524,5 +2526,5 @@
       !                                                      ! ------------------------- !
       ! needed by Met Office
-      CALL eos_fzp(tsn(:,:,1,jp_sal), sstfrz)
+      CALL eos_fzp(ts(:,:,1,jp_sal,Kmm), sstfrz)
       ztmp1(:,:) = sstfrz(:,:) + rt0
       IF( ssnd(jps_sstfrz)%laction )  CALL cpl_snd( jps_sstfrz, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info)
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcice_cice.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcice_cice.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcice_cice.F90	(revision 10922)
@@ -169,5 +169,5 @@
       ! Values from a CICE restart file would overwrite this
       IF ( .NOT. ln_rstart ) THEN    
-         CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.) 
+         CALL nemo2cice( ts(:,:,1,jp_tem,Kmm) , sst , 'T' , 1.) 
       ENDIF  
 #endif
@@ -194,6 +194,6 @@
 ! Ensure ocean temperatures are nowhere below freezing if not a NEMO restart
       IF( .NOT. ln_rstart ) THEN
-         tsn(:,:,:,jp_tem) = MAX (tsn(:,:,:,jp_tem),Tocnfrz)
-         tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
+         ts(:,:,:,jp_tem,Kmm) = MAX (ts(:,:,:,jp_tem,Kmm),Tocnfrz)
+         ts(:,:,:,jp_tem,Kbb) = ts(:,:,:,jp_tem,Kmm)
       ENDIF
 
@@ -235,33 +235,33 @@
                !
                DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
-                  e3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
-                  e3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
+                  e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
+                  e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
                ENDDO
-               e3t_a(:,:,:) = e3t_b(:,:,:)
+               e3t(:,:,:,Krhs) = e3t(:,:,:,Kbb)
                ! Reconstruction of all vertical scale factors at now and before time-steps
                ! =============================================================================
                ! Horizontal scale factor interpolations
                ! --------------------------------------
-               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
-               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
-               CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
-               CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
-               CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
+               CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )
+               CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )
+               CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
+               CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
+               CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )
                ! Vertical scale factor interpolations
                ! ------------------------------------
-               CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
-               CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
-               CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
-               CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
-               CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
+               CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )
+               CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
+               CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
+               CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
+               CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
                ! t- and w- points depth
                ! ----------------------
-               gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
-               gdepw_n(:,:,1) = 0.0_wp
-               gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
+               gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
+               gdepw(:,:,1,Kmm) = 0.0_wp
+               gde3w(:,:,1)     = gdept(:,:,1,Kmm) - sshn(:,:)
                DO jk = 2, jpk
-                  gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk)
-                  gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
-                  gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn   (:,:)
+                  gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk,Kmm)
+                  gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm)
+                  gde3w(:,:,jk)     = gdept(:,:,jk  ,Kmm) - sshn   (:,:)
                END DO
             ENDIF
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcice_if.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcice_if.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcice_if.F90	(revision 10922)
@@ -42,5 +42,5 @@
 CONTAINS
 
-   SUBROUTINE sbc_ice_if( kt )
+   SUBROUTINE sbc_ice_if( kt, Kbb, Kmm )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE sbc_ice_if  ***
@@ -59,4 +59,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
       !
       INTEGER  ::   ji, jj     ! dummy loop indices
@@ -118,5 +119,5 @@
                ENDIF
 
-               tsn(ji,jj,1,jp_tem) = MAX( tsn(ji,jj,1,jp_tem), zt_fzp )     ! avoid over-freezing point temperature
+               ts(ji,jj,1,jp_tem,Kmm) = MAX( ts(ji,jj,1,jp_tem,Kmm), zt_fzp )     ! avoid over-freezing point temperature
 
                qsr(ji,jj) = ( 1. - zfr_obs ) * qsr(ji,jj)   ! solar heat flux : zero below observed ice cover
@@ -125,6 +126,6 @@
                !      # ztrp*(t-(tgel-1.))  if observed ice and no opa ice   (zfr_obs=1 fr_i=0)
                !      # ztrp*min(0,t-tgel)  if observed ice and opa ice      (zfr_obs=1 fr_i=1)
-               zqri = ztrp * ( tsb(ji,jj,1,jp_tem) - ( zt_fzp - 1.) )
-               zqrj = ztrp * MIN( 0., tsb(ji,jj,1,jp_tem) - zt_fzp )
+               zqri = ztrp * ( ts(ji,jj,1,jp_tem,Kbb) - ( zt_fzp - 1.) )
+               zqrj = ztrp * MIN( 0., ts(ji,jj,1,jp_tem,Kbb) - zt_fzp )
                zqrp = ( zfr_obs * ( (1. - fr_i(ji,jj) ) * zqri    &
                  &                 +      fr_i(ji,jj)   * zqrj ) ) * tmask(ji,jj,1)
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 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcisf.F90	(revision 10922)
@@ -75,5 +75,5 @@
 CONTAINS
  
-  SUBROUTINE sbc_isf( kt )
+  SUBROUTINE sbc_isf( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE sbc_isf  ***
@@ -89,4 +89,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
+      INTEGER, INTENT(in) ::   Kmm  ! ocean time level indices
       !
       INTEGER ::   ji, jj, jk   ! loop index
@@ -102,8 +103,8 @@
          CASE ( 1 )    ! realistic ice shelf formulation
             ! compute T/S/U/V for the top boundary layer
-            CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T')
-            CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T')
-            CALL sbc_isf_tbl(un(:,:,:)        ,utbl(:,:),'U')
-            CALL sbc_isf_tbl(vn(:,:,:)        ,vtbl(:,:),'V')
+            CALL sbc_isf_tbl(ts(:,:,:,jp_tem,Kmm),ttbl(:,:),'T',Kmm)
+            CALL sbc_isf_tbl(ts(:,:,:,jp_sal,Kmm),stbl(:,:),'T',Kmm)
+            CALL sbc_isf_tbl(uu(:,:,:,Kmm)       ,utbl(:,:),'U',Kmm)
+            CALL sbc_isf_tbl(vv(:,:,:,Kmm)       ,vtbl(:,:),'V',Kmm)
             ! iom print
             CALL iom_put('ttbl',ttbl(:,:))
@@ -113,5 +114,5 @@
             ! compute fwf and heat flux
             ! compute fwf and heat flux
-            IF( .NOT.l_isfcpl ) THEN    ;   CALL sbc_isf_cav (kt)
+            IF( .NOT.l_isfcpl ) THEN    ;   CALL sbc_isf_cav (kt, Kmm)
             ELSE                        ;   qisf(:,:)  = fwfisf(:,:) * rLfusisf  ! heat        flux
             ENDIF
@@ -119,5 +120,5 @@
          CASE ( 2 )    ! Beckmann and Goosse parametrisation 
             stbl(:,:)   = soce
-            CALL sbc_isf_bg03(kt)
+            CALL sbc_isf_bg03(kt, Kmm)
             !
          CASE ( 3 )    ! specified runoff in depth (Mathiot et al., XXXX in preparation)
@@ -179,14 +180,14 @@
                   ikb = misfkb(ji,jj)
                   DO jk = ikt, ikb - 1
-                     zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
-                     zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
-                     zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
+                     zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * e3t(ji,jj,jk,Kmm)
+                     zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * e3t(ji,jj,jk,Kmm)
+                     zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * e3t(ji,jj,jk,Kmm)
                   END DO
                   zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj)   & 
-                     &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk)
+                     &                                                                   * ralpha(ji,jj) * e3t(ji,jj,jk,Kmm)
                   zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj)   & 
-                     &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk)
+                     &                                                                   * ralpha(ji,jj) * e3t(ji,jj,jk,Kmm)
                   zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj)   &  
-                     &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk)
+                     &                                                                   * ralpha(ji,jj) * e3t(ji,jj,jk,Kmm)
                END DO
             END DO
@@ -251,5 +252,5 @@
 
 
-  SUBROUTINE sbc_isf_init
+  SUBROUTINE sbc_isf_init( Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE sbc_isf_init  ***
@@ -263,4 +264,5 @@
       !!                        4 : specified fwf and heat flux forcing beneath the ice shelf
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kmm  ! ocean time level indices
       INTEGER               :: ji, jj, jk           ! loop index
       INTEGER               :: ik                   ! current level index
@@ -355,5 +357,5 @@
                 ik = 2
 !!gm potential bug: use gdepw_0 not _n
-                DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_n(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO
+                DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw(ji,jj,ik,Kmm) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO
                 misfkt(ji,jj) = ik-1
             END DO
@@ -386,16 +388,16 @@
             ikb = misfkt(ji,jj)
             ! thickness of boundary layer at least the top level thickness
-            rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt))
+            rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t(ji,jj,ikt,Kmm))
 
             ! determine the deepest level influenced by the boundary layer
             DO jk = ikt+1, mbkt(ji,jj)
-               IF( (SUM(e3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) )   ikb = jk
-            END DO
-            rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
+               IF( (SUM(e3t(ji,jj,ikt:jk-1,Kmm)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) )   ikb = jk
+            END DO
+            rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t(ji,jj,ikt:ikb,Kmm))) ! limit the tbl to water thickness.
             misfkb(ji,jj) = ikb                                                   ! last wet level of the tbl
             r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
 
-            zhk           = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1
-            ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
+            zhk           = SUM( e3t(ji, jj, ikt:ikb - 1,Kmm)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1
+            ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t(ji,jj,ikb,Kmm)  ! proportion of bottom cell influenced by boundary layer
          END DO
       END DO
@@ -411,5 +413,5 @@
 
 
-  SUBROUTINE sbc_isf_bg03(kt)
+  SUBROUTINE sbc_isf_bg03( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE sbc_isf_bg03  ***
@@ -426,4 +428,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT ( in ) :: kt
+      INTEGER, INTENT ( in ) :: Kmm  ! ocean time level indices
       !
       INTEGER  :: ji, jj, jk ! dummy loop index
@@ -444,7 +447,7 @@
                DO jk = misfkt(ji,jj),misfkb(ji,jj)
                   ! Calculate freezing temperature
-                  zpress = grav*rau0*gdept_n(ji,jj,ik)*1.e-04
+                  zpress = grav*rau0*gdept(ji,jj,ik,Kmm)*1.e-04
                   CALL eos_fzp(stbl(ji,jj), zt_frz, zpress) 
-                  zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk)  ! sum temp
+                  zt_sum = zt_sum + (ts(ji,jj,jk,jp_tem,Kmm)-zt_frz) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)  ! sum temp
                END DO
                zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
@@ -466,5 +469,5 @@
 
 
-  SUBROUTINE sbc_isf_cav( kt )
+  SUBROUTINE sbc_isf_cav( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE sbc_isf_cav  ***
@@ -480,4 +483,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
+      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
       !
       INTEGER  ::   ji, jj     ! dummy loop indices
@@ -520,5 +524,5 @@
 
             ! compute gammat every where (2d)
-            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
+            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, Kmm)
             
             ! compute upward heat flux zhtflx and upward water flux zwflx
@@ -536,5 +540,5 @@
          CASE ( 2 )  ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)
             ! compute gammat every where (2d)
-            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
+            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, Kmm)
 
             ! compute upward heat flux zhtflx and upward water flux zwflx
@@ -600,5 +604,5 @@
 
 
-   SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf )
+   SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf, Kmm )
       !!----------------------------------------------------------------------
       !! ** Purpose    : compute the coefficient echange for heat flux
@@ -611,4 +615,5 @@
       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pgt   , pgs      ! 
       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pqhisf, pqwisf   ! 
+      INTEGER                 , INTENT(in   ) ::   Kmm  ! ocean time level indices
       !
       INTEGER  :: ji, jj                     ! loop index
@@ -679,10 +684,10 @@
 !!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation
 !!gm moreover, use Max(rn2,0) to take care of static instabilities....
-                  zcoef = 0.5_wp / e3w_n(ji,jj,ikt+1)
+                  zcoef = 0.5_wp / e3w(ji,jj,ikt+1,Kmm)
                   !                                            ! shear of horizontal velocity
-                  zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  &
-                     &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
-                  zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )  &
-                     &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
+                  zdku = zcoef * (  uu(ji-1,jj  ,ikt  ,Kmm) + uu(ji,jj,ikt  ,Kmm)  &
+                     &             -uu(ji-1,jj  ,ikt+1,Kmm) - uu(ji,jj,ikt+1,Kmm)  )
+                  zdkv = zcoef * (  vv(ji  ,jj-1,ikt  ,Kmm) + vv(ji,jj,ikt  ,Kmm)  &
+                     &             -vv(ji  ,jj-1,ikt+1,Kmm) - vv(ji,jj,ikt+1,Kmm)  )
                   !                                            ! richardson number (minimum value set to zero)
                   zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps )
@@ -691,5 +696,5 @@
                   zts(jp_tem) = ttbl(ji,jj)
                   zts(jp_sal) = stbl(ji,jj)
-                  zdep        = gdepw_n(ji,jj,ikt)
+                  zdep        = gdepw(ji,jj,ikt,Kmm)
                   !
                   CALL eos_rab( zts, zdep, zab )
@@ -700,5 +705,5 @@
                   !! compute Monin Obukov Length
                   ! Maximum boundary layer depth
-                  zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp
+                  zhmax = gdept(ji,jj,mbkt(ji,jj),Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm) -0.001_wp
                   ! Compute Monin obukhov length scale at the surface and Ekman depth:
                   zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
@@ -727,5 +732,5 @@
 
 
-   SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin )
+   SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE sbc_isf_tbl  ***
@@ -737,4 +742,5 @@
       REAL(wp), DIMENSION(:,:)  , INTENT(  out) :: pvarout
       CHARACTER(len=1),           INTENT(in   ) :: cd_ptin ! point of variable in/out
+      INTEGER                   , INTENT(in   ) :: Kmm  ! ocean time level indices
       !
       INTEGER ::   ji, jj, jk                ! loop index
@@ -753,20 +759,20 @@
                ikt = miku(ji,jj) ; ikb = miku(ji,jj)
                ! thickness of boundary layer at least the top level thickness
-               zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj) , e3u_n(ji,jj,ikt) )
+               zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj) , e3u(ji,jj,ikt,Kmm) )
 
                ! determine the deepest level influenced by the boundary layer
                DO jk = ikt+1, mbku(ji,jj)
-                  IF ( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
-               END DO
-               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
+                  IF ( (SUM(e3u(ji,jj,ikt:jk-1,Kmm)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
+               END DO
+               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u(ji,jj,ikt:ikb,Kmm)))  ! limit the tbl to water thickness.
 
                ! level fully include in the ice shelf boundary layer
                DO jk = ikt, ikb - 1
-                  ze3 = e3u_n(ji,jj,jk)
+                  ze3 = e3u(ji,jj,jk,Kmm)
                   pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
                END DO
 
                ! level partially include in ice shelf boundary layer 
-               zhk = SUM( e3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
+               zhk = SUM( e3u(ji, jj, ikt:ikb - 1,Kmm)) / zhisf_tbl(ji,jj)
                pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
             END DO
@@ -785,20 +791,20 @@
                ikt = mikv(ji,jj) ; ikb = mikv(ji,jj)
                ! thickness of boundary layer at least the top level thickness
-               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v_n(ji,jj,ikt))
+               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v(ji,jj,ikt,Kmm))
 
                ! determine the deepest level influenced by the boundary layer
                DO jk = ikt+1, mbkv(ji,jj)
-                  IF ( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
-               END DO
-               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
+                  IF ( (SUM(e3v(ji,jj,ikt:jk-1,Kmm)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
+               END DO
+               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v(ji,jj,ikt:ikb,Kmm)))  ! limit the tbl to water thickness.
 
                ! level fully include in the ice shelf boundary layer
                DO jk = ikt, ikb - 1
-                  ze3 = e3v_n(ji,jj,jk)
+                  ze3 = e3v(ji,jj,jk,Kmm)
                   pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
                END DO
 
                ! level partially include in ice shelf boundary layer 
-               zhk = SUM( e3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
+               zhk = SUM( e3v(ji, jj, ikt:ikb - 1,Kmm)) / zhisf_tbl(ji,jj)
                pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
             END DO
@@ -820,10 +826,10 @@
                ! level fully include in the ice shelf boundary layer
                DO jk = ikt, ikb - 1
-                  ze3 = e3t_n(ji,jj,jk)
+                  ze3 = e3t(ji,jj,jk,Kmm)
                   pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
                END DO
 
                ! level partially include in ice shelf boundary layer 
-               zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
+               zhk = SUM( e3t(ji, jj, ikt:ikb - 1,Kmm)) * r1_hisf_tbl(ji,jj)
                pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
             END DO
@@ -837,5 +843,5 @@
       
 
-   SUBROUTINE sbc_isf_div( phdivn )
+   SUBROUTINE sbc_isf_div( phdivn, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE sbc_isf_div  ***
@@ -850,4 +856,5 @@
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdivn   ! horizontal divergence
+      INTEGER                   , INTENT( in    ) ::   Kmm      ! ocean time level indices
       ! 
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -865,16 +872,16 @@
                ikb = misfkt(ji,jj)
                ! thickness of boundary layer at least the top level thickness
-               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt))
+               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t(ji,jj,ikt,Kmm))
 
                ! determine the deepest level influenced by the boundary layer
                DO jk = ikt, mbkt(ji,jj)
-                  IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
-               END DO
-               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
+                  IF ( (SUM(e3t(ji,jj,ikt:jk-1,Kmm)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
+               END DO
+               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t(ji,jj,ikt:ikb,Kmm)))  ! limit the tbl to water thickness.
                misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
                r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
 
-               zhk           = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
-               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
+               zhk           = SUM( e3t(ji, jj, ikt:ikb - 1,Kmm)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
+               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t(ji,jj,ikb,Kmm)  ! proportion of bottom cell influenced by boundary layer
             END DO
          END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcmod.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcmod.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcmod.F90	(revision 10922)
@@ -76,5 +76,5 @@
 CONTAINS
 
-   SUBROUTINE sbc_init
+   SUBROUTINE sbc_init( Kbb, Kmm )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE sbc_init ***
@@ -88,4 +88,5 @@
       !!              - nsbc: type of sbc
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kbb, Kmm              ! ocean time level indices
       INTEGER ::   ios, icpt                         ! local integer
       LOGICAL ::   ll_purecpl, ll_opa, ll_not_nemo   ! local logical
@@ -323,15 +324,15 @@
       !                       !**  associated modules : initialization
       !
-                          CALL sbc_ssm_init            ! Sea-surface mean fields initialization
-      !
-      IF( ln_blk      )   CALL sbc_blk_init            ! bulk formulae initialization
-
-      IF( ln_ssr      )   CALL sbc_ssr_init            ! Sea-Surface Restoring initialization
-      !
-      IF( ln_isf      )   CALL sbc_isf_init            ! Compute iceshelves
-      !
-                          CALL sbc_rnf_init            ! Runof initialization
-      !
-      IF( ln_apr_dyn )    CALL sbc_apr_init            ! Atmo Pressure Forcing initialization
+                          CALL sbc_ssm_init ( Kbb, Kmm ) ! Sea-surface mean fields initialization
+      !
+      IF( ln_blk      )   CALL sbc_blk_init              ! bulk formulae initialization
+
+      IF( ln_ssr      )   CALL sbc_ssr_init              ! Sea-Surface Restoring initialization
+      !
+      IF( ln_isf      )   CALL sbc_isf_init( Kmm )       ! Compute iceshelves
+      !
+                          CALL sbc_rnf_init( Kmm )       ! Runof initialization
+      !
+      IF( ln_apr_dyn )    CALL sbc_apr_init              ! Atmo Pressure Forcing initialization
       !
 #if defined key_si3
@@ -359,5 +360,5 @@
 
 
-   SUBROUTINE sbc( kt )
+   SUBROUTINE sbc( kt, Kbb, Kmm )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE sbc  ***
@@ -376,4 +377,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
       !
       LOGICAL ::   ll_sas, ll_opa   ! local logical
@@ -411,6 +413,6 @@
       ll_opa = nn_components == jp_iam_opa
       !
-      IF( .NOT.ll_sas )   CALL sbc_ssm ( kt )            ! mean ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
-      IF( ln_wave     )   CALL sbc_wave( kt )            ! surface waves
+      IF( .NOT.ll_sas )   CALL sbc_ssm ( kt, Kbb, Kmm )  ! mean ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
+      IF( ln_wave     )   CALL sbc_wave( kt, Kmm )       ! surface waves
 
       !
@@ -419,23 +421,23 @@
       SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition
       !                                                  ! (i.e. utau,vtau, qns, qsr, emp, sfx)
-      CASE( jp_usr   )     ;   CALL usrdef_sbc_oce( kt )                    ! user defined formulation 
-      CASE( jp_flx     )   ;   CALL sbc_flx       ( kt )                    ! flux formulation
+      CASE( jp_usr   )     ;   CALL usrdef_sbc_oce( kt )                             ! user defined formulation 
+      CASE( jp_flx     )   ;   CALL sbc_flx       ( kt )                             ! flux formulation
       CASE( jp_blk     )
-         IF( ll_sas    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: SAS receiving fields from OPA
-                               CALL sbc_blk       ( kt )                    ! bulk formulation for the ocean
+         IF( ll_sas    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )  ! OPA-SAS coupling: SAS receiving fields from OPA
+                               CALL sbc_blk       ( kt )                             ! bulk formulation for the ocean
                                !
-      CASE( jp_purecpl )   ;   CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! pure coupled formulation
+      CASE( jp_purecpl )   ;   CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )  ! pure coupled formulation
       CASE( jp_none    )
-         IF( ll_opa    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: OPA receiving fields from SAS
+         IF( ll_opa    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )  ! OPA-SAS coupling: OPA receiving fields from SAS
       END SELECT
       !
-      IF( ln_mixcpl )          CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing
-      !
-      IF ( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( )      ! Wind stress provided by waves 
+      IF( ln_mixcpl )          CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )  ! forced-coupled mixed formulation after forcing
+      !
+      IF ( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( )              ! Wind stress provided by waves 
       !
       !                                            !==  Misc. Options  ==!
       !
       SELECT CASE( nn_ice )                                       ! Update heat and freshwater fluxes over sea-ice areas
-      CASE(  1 )   ;         CALL sbc_ice_if   ( kt )             ! Ice-cover climatology ("Ice-if" model)
+      CASE(  1 )   ;         CALL sbc_ice_if   ( kt, Kbb, Kmm )   ! Ice-cover climatology ("Ice-if" model)
 #if defined key_si3
       CASE(  2 )   ;         CALL ice_stp  ( kt, nsbc )           ! SI3 ice model
@@ -451,5 +453,5 @@
       ENDIF
 
-      IF( ln_isf         )   CALL sbc_isf( kt )                   ! compute iceshelves
+      IF( ln_isf         )   CALL sbc_isf( kt, Kmm )              ! compute iceshelves
 
       IF( ln_rnf         )   CALL sbc_rnf( kt )                   ! add runoffs to fresh water fluxes
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcrnf.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcrnf.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcrnf.F90	(revision 10922)
@@ -173,5 +173,5 @@
 
 
-   SUBROUTINE sbc_rnf_div( phdivn )
+   SUBROUTINE sbc_rnf_div( phdivn, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE sbc_rnf  ***
@@ -185,4 +185,5 @@
       !! ** Action  :   phdivn   decreased by the runoff inflow
       !!----------------------------------------------------------------------
+      INTEGER                   , INTENT(in   ) ::   Kmm      ! ocean time level index
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
       !!
@@ -207,5 +208,5 @@
                   h_rnf(ji,jj) = 0._wp
                   DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres
-                     h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk)   ! to the bottom of the relevant grid box
+                     h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)   ! to the bottom of the relevant grid box
                   END DO
                   !                          ! apply the runoff input flow
@@ -217,6 +218,6 @@
          ENDIF
       ELSE                       !==   runoff put only at the surface   ==!
-         h_rnf (:,:)   = e3t_n (:,:,1)        ! update h_rnf to be depth of top box
-         phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:) + rnf_b(:,:) ) * zfact * r1_rau0 / e3t_n(:,:,1)
+         h_rnf (:,:)   = e3t (:,:,1,Kmm)        ! update h_rnf to be depth of top box
+         phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:) + rnf_b(:,:) ) * zfact * r1_rau0 / e3t(:,:,1,Kmm)
       ENDIF
       !
@@ -224,5 +225,5 @@
 
 
-   SUBROUTINE sbc_rnf_init
+   SUBROUTINE sbc_rnf_init( Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE sbc_rnf_init  ***
@@ -234,4 +235,5 @@
       !! ** Action  : - read parameters
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kmm           ! ocean time level index
       CHARACTER(len=32) ::   rn_dep_file   ! runoff file name
       INTEGER           ::   ji, jj, jk, jm    ! dummy loop indices
@@ -356,5 +358,5 @@
                h_rnf(ji,jj) = 0._wp
                DO jk = 1, nk_rnf(ji,jj)
-                  h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk)
+                  h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)
                END DO
             END DO
@@ -411,5 +413,5 @@
                h_rnf(ji,jj) = 0._wp
                DO jk = 1, nk_rnf(ji,jj)
-                  h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk)
+                  h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)
                END DO
             END DO
@@ -424,5 +426,5 @@
       ELSE                                       ! runoffs applied at the surface
          nk_rnf(:,:) = 1
-         h_rnf (:,:) = e3t_n(:,:,1)
+         h_rnf (:,:) = e3t(:,:,1,Kmm)
       ENDIF
       !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcssm.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcssm.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcssm.F90	(revision 10922)
@@ -39,5 +39,5 @@
 CONTAINS
 
-   SUBROUTINE sbc_ssm( kt )
+   SUBROUTINE sbc_ssm( kt, Kbb, Kmm )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE sbc_oce  ***
@@ -53,4 +53,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
       !
       INTEGER  ::   ji, jj               ! loop index
@@ -62,6 +63,6 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            zts(ji,jj,jp_tem) = tsn(ji,jj,mikt(ji,jj),jp_tem)
-            zts(ji,jj,jp_sal) = tsn(ji,jj,mikt(ji,jj),jp_sal)
+            zts(ji,jj,jp_tem) = ts(ji,jj,mikt(ji,jj),jp_tem,Kmm)
+            zts(ji,jj,jp_sal) = ts(ji,jj,mikt(ji,jj),jp_sal,Kmm)
          END DO
       END DO
@@ -69,6 +70,6 @@
       IF( nn_fsbc == 1 ) THEN                             !   Instantaneous surface fields        !
          !                                                ! ---------------------------------------- !
-         ssu_m(:,:) = ub(:,:,1)
-         ssv_m(:,:) = vb(:,:,1)
+         ssu_m(:,:) = uu(:,:,1,Kbb)
+         ssv_m(:,:) = vv(:,:,1,Kbb)
          IF( l_useCT )  THEN    ;   sst_m(:,:) = eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
          ELSE                    ;   sst_m(:,:) = zts(:,:,jp_tem)
@@ -80,5 +81,5 @@
          ENDIF
          !
-         e3t_m(:,:) = e3t_n(:,:,1)
+         e3t_m(:,:) = e3t(:,:,1,Kmm)
          !
          frq_m(:,:) = fraqsr_1lev(:,:)
@@ -92,6 +93,6 @@
             IF(lwp) WRITE(numout,*) '~~~~~~~   '
             zcoef = REAL( nn_fsbc - 1, wp )
-            ssu_m(:,:) = zcoef * ub(:,:,1)
-            ssv_m(:,:) = zcoef * vb(:,:,1)
+            ssu_m(:,:) = zcoef * uu(:,:,1,Kbb)
+            ssv_m(:,:) = zcoef * vv(:,:,1,Kbb)
             IF( l_useCT )  THEN    ;   sst_m(:,:) = zcoef * eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
             ELSE                    ;   sst_m(:,:) = zcoef * zts(:,:,jp_tem)
@@ -103,5 +104,5 @@
             ENDIF
             !
-            e3t_m(:,:) = zcoef * e3t_n(:,:,1)
+            e3t_m(:,:) = zcoef * e3t(:,:,1,Kmm)
             !
             frq_m(:,:) = zcoef * fraqsr_1lev(:,:)
@@ -120,6 +121,6 @@
          !                                                !        Cumulate at each time step        !
          !                                                ! ---------------------------------------- !
-         ssu_m(:,:) = ssu_m(:,:) + ub(:,:,1)
-         ssv_m(:,:) = ssv_m(:,:) + vb(:,:,1)
+         ssu_m(:,:) = ssu_m(:,:) + uu(:,:,1,Kbb)
+         ssv_m(:,:) = ssv_m(:,:) + vv(:,:,1,Kbb)
          IF( l_useCT )  THEN     ;   sst_m(:,:) = sst_m(:,:) + eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
          ELSE                    ;   sst_m(:,:) = sst_m(:,:) + zts(:,:,jp_tem)
@@ -131,5 +132,5 @@
          ENDIF
          !
-         e3t_m(:,:) = e3t_m(:,:) + e3t_n(:,:,1)
+         e3t_m(:,:) = e3t_m(:,:) + e3t(:,:,1,Kmm)
          !
          frq_m(:,:) = frq_m(:,:) + fraqsr_1lev(:,:)
@@ -184,5 +185,5 @@
 
 
-   SUBROUTINE sbc_ssm_init
+   SUBROUTINE sbc_ssm_init( Kbb, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE sbc_ssm_init  ***
@@ -192,4 +193,5 @@
       !! ** Action  : - read parameters
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
       REAL(wp) ::   zcoef, zf_sbc   ! local scalar
       !!----------------------------------------------------------------------
@@ -242,12 +244,12 @@
          !
          IF(lwp) WRITE(numout,*) '   default initialisation of ss._m arrays'
-         ssu_m(:,:) = ub(:,:,1)
-         ssv_m(:,:) = vb(:,:,1)
-         IF( l_useCT )  THEN    ;   sst_m(:,:) = eos_pt_from_ct( tsn(:,:,1,jp_tem), tsn(:,:,1,jp_sal) )
-         ELSE                   ;   sst_m(:,:) = tsn(:,:,1,jp_tem)
-         ENDIF
-         sss_m(:,:) = tsn  (:,:,1,jp_sal)
-         ssh_m(:,:) = sshn (:,:)
-         e3t_m(:,:) = e3t_n(:,:,1)
+         ssu_m(:,:) = uu(:,:,1,Kbb)
+         ssv_m(:,:) = vv(:,:,1,Kbb)
+         IF( l_useCT )  THEN    ;   sst_m(:,:) = eos_pt_from_ct( ts(:,:,1,jp_tem,Kmm), ts(:,:,1,jp_sal,Kmm) )
+         ELSE                   ;   sst_m(:,:) = ts(:,:,1,jp_tem,Kmm)
+         ENDIF
+         sss_m(:,:) = ts  (:,:,1,jp_sal,Kmm)
+         ssh_m(:,:) = sshn(:,:)
+         e3t_m(:,:) = e3t (:,:,1,Kmm)
          frq_m(:,:) = 1._wp
          !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcwave.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcwave.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcwave.F90	(revision 10922)
@@ -80,5 +80,5 @@
 CONTAINS
 
-   SUBROUTINE sbc_stokes( )
+   SUBROUTINE sbc_stokes( Kmm )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE sbc_stokes  ***
@@ -92,4 +92,5 @@
       !! ** action  
       !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kmm ! ocean time level index
       INTEGER  ::   jj, ji, jk   ! dummy loop argument
       INTEGER  ::   ik           ! local integer 
@@ -152,6 +153,6 @@
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
-                  zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) )
-                  zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) )
+                  zdep_u = 0.5_wp * ( gdept(ji,jj,jk,Kmm) + gdept(ji+1,jj,jk,Kmm) )
+                  zdep_v = 0.5_wp * ( gdept(ji,jj,jk,Kmm) + gdept(ji,jj+1,jk,Kmm) )
                   !                          
                   zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth
@@ -179,11 +180,11 @@
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
-                  zbot_u = ( gdepw_n(ji,jj,jk+1) + gdepw_n(ji+1,jj,jk+1) )  ! 2 * bottom depth
-                  zbot_v = ( gdepw_n(ji,jj,jk+1) + gdepw_n(ji,jj+1,jk+1) )  ! 2 * bottom depth
+                  zbot_u = ( gdepw(ji,jj,jk+1,Kmm) + gdepw(ji+1,jj,jk+1,Kmm) )  ! 2 * bottom depth
+                  zbot_v = ( gdepw(ji,jj,jk+1,Kmm) + gdepw(ji,jj+1,jk+1,Kmm) )  ! 2 * bottom depth
                   zkb_u  = zk_u(ji,jj) * zbot_u                             ! 2 * k * bottom depth
                   zkb_v  = zk_v(ji,jj) * zbot_v                             ! 2 * k * bottom depth
                   !
-                  zke3_u = MAX(1.e-8_wp, 2.0_wp * zk_u(ji,jj) * e3u_n(ji,jj,jk))     ! 2k * thickness
-                  zke3_v = MAX(1.e-8_wp, 2.0_wp * zk_v(ji,jj) * e3v_n(ji,jj,jk))     ! 2k * thickness
+                  zke3_u = MAX(1.e-8_wp, 2.0_wp * zk_u(ji,jj) * e3u(ji,jj,jk,Kmm))     ! 2k * thickness
+                  zke3_v = MAX(1.e-8_wp, 2.0_wp * zk_v(ji,jj) * e3v(ji,jj,jk,Kmm))     ! 2k * thickness
 
                   ! Depth attenuation .... do u component first..
@@ -223,8 +224,8 @@
          DO jj = 2, jpj
             DO ji = fs_2, jpi
-               ze3divh(ji,jj,jk) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * usd(ji  ,jj,jk)    &
-                  &                 - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * usd(ji-1,jj,jk)    &
-                  &                 + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vsd(ji,jj  ,jk)    &
-                  &                 - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vsd(ji,jj-1,jk)  ) * r1_e1e2t(ji,jj)
+               ze3divh(ji,jj,jk) = (  e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * usd(ji  ,jj,jk)    &
+                  &                 - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * usd(ji-1,jj,jk)    &
+                  &                 + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * vsd(ji,jj  ,jk)    &
+                  &                 - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vsd(ji,jj-1,jk)  ) * r1_e1e2t(ji,jj)
             END DO
          END DO
@@ -307,5 +308,5 @@
 
 
-   SUBROUTINE sbc_wave( kt )
+   SUBROUTINE sbc_wave( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE sbc_wave  ***
@@ -322,4 +323,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in   ) ::   kt   ! ocean time step
+      INTEGER, INTENT(in   ) ::   Kmm  ! ocean time index
       !!---------------------------------------------------------------------
       !
@@ -361,5 +363,5 @@
          !
          IF( ( ll_st_bv_li   .AND. jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0 ) .OR. &
-           & ( ll_st_peakfr  .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0                ) ) CALL sbc_stokes()
+           & ( ll_st_peakfr  .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0                ) ) CALL sbc_stokes( Kmm )
          !
       ENDIF
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 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv.F90	(revision 10922)
@@ -127,5 +127,5 @@
       !
       IF( ln_ldfeiv .AND. .NOT. ln_traldf_triad )   &
-         &              CALL ldf_eiv_trp( kt, nit000, zuu, zvv, zww, 'TRA' )   ! add the eiv transport (if necessary)
+         &              CALL ldf_eiv_trp( kt, nit000, zuu, zvv, zww, 'TRA', Kmm )   ! 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)
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 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf.F90	(revision 10922)
@@ -47,5 +47,5 @@
 CONTAINS
 
-   SUBROUTINE tra_ldf( kt )
+   SUBROUTINE tra_ldf( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf  ***
@@ -54,4 +54,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kmm  ! ocean time level indices
       !!
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
@@ -72,7 +73,7 @@
          CALL tra_ldf_iso  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   )
       CASE ( np_lap_it )                                 ! laplacian: triad iso-neutral operator (griffies)
-         CALL tra_ldf_triad( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   )
+         CALL tra_ldf_triad( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, 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 )
+         CALL tra_ldf_blp  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb      , tsa, jpts, nldf_tra, Kmm )
       END SELECT
       !
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 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_lap_blp.F90	(revision 10922)
@@ -161,5 +161,5 @@
    SUBROUTINE tra_ldf_blp( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
       &                                                    pgui, pgvi,   &
-      &                                                    ptb , pta , kjpt, kldf )
+      &                                                    ptb , pta , kjpt, kldf, Kmm )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE tra_ldf_blp  ***
@@ -179,4 +179,5 @@
       INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
       INTEGER                              , INTENT(in   ) ::   kldf       ! type of operator used
+      INTEGER                              , INTENT(in   ) ::   Kmm        ! ocean time level indices
       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
@@ -210,5 +211,5 @@
          CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 )
       CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies)
-         CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 )
+         CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1, Kmm )
       END SELECT
       !
@@ -226,5 +227,5 @@
          CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 )
       CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies)
-         CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 )
+         CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2, Kmm )
       END SELECT
       !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_triad.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_triad.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_triad.F90	(revision 10922)
@@ -50,5 +50,5 @@
   SUBROUTINE tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
       &                                                     pgui, pgvi,   &
-      &                                         ptb , ptbb, pta , kjpt, kpass )
+      &                                         ptb , ptbb, pta , kjpt, kpass, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf_triad  ***
@@ -75,4 +75,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 indices
       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
@@ -213,5 +214,5 @@
          ENDIF
          !
-         IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' )   CALL ldf_eiv_dia( zpsi_uw, zpsi_vw )
+         IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' )   CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm )
          !
       ENDIF                                  !==  end 1st pass only  ==!
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 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90	(revision 10922)
@@ -428,10 +428,10 @@
       !                            
       IF( ln_diurnal_only ) THEN                   ! diurnal only: a subset of the initialisation routines
-         CALL  istate_init                            ! ocean initial state (Dynamics and tracers)
-         CALL     sbc_init                            ! Forcings : surface module
+         CALL  istate_init( Nbb, Nnn )                ! ocean initial state (Dynamics and tracers)
+         CALL     sbc_init( Nbb, Nnn )                ! Forcings : surface module
          CALL tra_qsr_init                            ! penetrative solar radiation qsr
          IF( ln_diaobs ) THEN                         ! Observation & model comparison
-            CALL dia_obs_init                            ! Initialize observational data
-            CALL dia_obs( nit000 - 1 )                   ! Observation operator for restart
+            CALL dia_obs_init( Nnn )                     ! Initialize observational data
+            CALL dia_obs( nit000 - 1, Nnn )              ! Observation operator for restart
          ENDIF     
          IF( lk_asminc )   CALL asm_inc_init          ! Assimilation increments
@@ -440,10 +440,10 @@
       ENDIF
       
-                           CALL  istate_init    ! ocean initial state (Dynamics and tracers)
+                           CALL  istate_init( Nbb, Nnn ) ! ocean initial state (Dynamics and tracers)
 
       !                                      ! external forcing 
-                           CALL    tide_init    ! tidal harmonics
-                           CALL     sbc_init    ! surface boundary conditions (including sea-ice)
-                           CALL     bdy_init    ! Open boundaries initialisation
+                           CALL    tide_init             ! tidal harmonics
+                           CALL     sbc_init( Nbb, Nnn ) ! surface boundary conditions (including sea-ice)
+                           CALL     bdy_init             ! Open boundaries initialisation
 
       !                                      ! Ocean physics
@@ -491,8 +491,8 @@
                            CALL dia_hsb_init    ! heat content, salt content and volume budgets
                            CALL     trd_init    ! Mixed-layer/Vorticity/Integral constraints trends
-                           CALL dia_obs_init    ! Initialize observational data
+                           CALL dia_obs_init( Nnn )        ! Initialize observational data
                            CALL dia_tmb_init    ! TMB outputs
                            CALL dia_25h_init    ! 25h mean  outputs
-      IF( ln_diaobs    )   CALL dia_obs( nit000-1 )   ! Observation operator for restart
+      IF( ln_diaobs    )   CALL dia_obs( nit000-1, Nnn )   ! Observation operator for restart
 
       !                                      ! Assimilation increments
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 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10922)
@@ -119,5 +119,5 @@
       IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                   ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib) 
       IF( ln_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
-                         CALL sbc     ( kstp )                   ! Sea Boundary Condition (including sea-ice)
+                         CALL sbc     ( kstp, Nbb, Nnn )         ! Sea Boundary Condition (including sea-ice)
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -152,12 +152,12 @@
             &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
          IF( ln_traldf_triad ) THEN 
-                         CALL ldf_slp_triad( kstp )                       ! before slope for triad operator
+                         CALL ldf_slp_triad( kstp, Nbb, Nnn )             ! before slope for triad operator
          ELSE     
-                         CALL ldf_slp     ( kstp, rhd, rn2b )             ! before slope for standard operator
+                         CALL ldf_slp     ( kstp, rhd, rn2b, Nbb, Nnn )   ! before slope for standard operator
          ENDIF
       ENDIF
-      !                                                                   ! eddy diffusivity coeff.
-      IF( l_ldftra_time .OR. l_ldfeiv_time )   CALL ldf_tra( kstp )       !       and/or eiv coeff.
-      IF( l_ldfdyn_time                    )   CALL ldf_dyn( kstp )       ! eddy viscosity coeff. 
+      !                                                                        ! eddy diffusivity coeff.
+      IF( l_ldftra_time .OR. l_ldfeiv_time )   CALL ldf_tra( kstp, Nbb, Nnn )  !       and/or eiv coeff.
+      IF( l_ldfdyn_time                    )   CALL ldf_dyn( kstp, Nbb )       ! eddy viscosity coeff. 
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -165,5 +165,5 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
-                            CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_hor)
+                            CALL ssh_nxt       ( kstp, Nnn )  ! after ssh (includes call to div_hor)
       IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
                             CALL wzv           ( kstp )  ! now cross-level velocity 
@@ -202,5 +202,5 @@
                                                       ! With split-explicit free surface, since now transports have been updated and ssha as well
       IF( ln_dynspg_ts ) THEN                         ! vertical scale factors and vertical velocity need to be updated
-                            CALL div_hor    ( kstp )              ! Horizontal divergence  (2nd call in time-split case)
+                            CALL div_hor    ( kstp, Nnn )         ! Horizontal divergence  (2nd call in time-split case)
          IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component)
                             CALL wzv        ( kstp )              ! now cross-level velocity 
@@ -256,5 +256,5 @@
       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 )  ! lateral mixing
+                         CALL tra_ldf       ( kstp, Nnn )  ! lateral mixing
 
 !!gm : why CALL to dia_ptr has been moved here??? (use trends info?)
@@ -293,5 +293,5 @@
 !!jc: That would be better, but see comment above
 !!
-      IF( lrst_oce   )   CALL rst_write    ( kstp )   ! write output ocean restart file
+      IF( lrst_oce   )   CALL rst_write    ( kstp, Nbb, Nnn )   ! write output ocean restart file
       IF( ln_sto_eos )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters
 
@@ -304,5 +304,5 @@
                          IF( Agrif_NbStepint() == 0 ) CALL Agrif_update_all( ) ! Update all components
 #endif
-      IF( ln_diaobs  )   CALL dia_obs      ( kstp )      ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
+      IF( ln_diaobs  )   CALL dia_obs      ( kstp, Nnn )      ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -321,5 +321,5 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 !!gm why lk_oasis and not lk_cpl ????
-      IF( lk_oasis   )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
+      IF( lk_oasis   )   CALL sbc_cpl_snd( kstp, Nnn )     ! coupled mode : field exchanges
       !
 #if defined key_iomput
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 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/dtadyn.F90	(revision 10922)
@@ -122,5 +122,5 @@
       CALL fld_read( kt, 1, sf_dyn )      !=  read data at kt time step   ==!
       !
-      IF( l_ldfslp .AND. .NOT.lk_c1d )   CALL  dta_dyn_slp( kt )    ! Computation of slopes
+      IF( l_ldfslp .AND. .NOT.lk_c1d )   CALL  dta_dyn_slp( kt, Kbb, Kmm )    ! Computation of slopes
       !
       ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature
@@ -679,5 +679,5 @@
 
 
-   SUBROUTINE dta_dyn_slp( kt )
+   SUBROUTINE dta_dyn_slp( kt, Kbb, Kmm )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE dta_dyn_slp  ***
@@ -687,4 +687,5 @@
       !!---------------------------------------------------------------------
       INTEGER,  INTENT(in) :: kt       ! time step
+      INTEGER,  INTENT(in) :: Kbb, Kmm ! ocean time level indices
       !
       INTEGER  ::   ji, jj     ! dummy loop indices
@@ -702,5 +703,5 @@
             zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity 
             avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:)   ! vertical diffusive coef.
-            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
+            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm )
             uslpdta (:,:,:,1) = zuslp (:,:,:) 
             vslpdta (:,:,:,1) = zvslp (:,:,:) 
@@ -711,5 +712,5 @@
             zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity 
             avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef.
-            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
+            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm )
             uslpdta (:,:,:,2) = zuslp (:,:,:) 
             vslpdta (:,:,:,2) = zvslp (:,:,:) 
@@ -730,5 +731,5 @@
               zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity 
               avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef.
-              CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
+              CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm )
               !
               uslpdta (:,:,:,2) = zuslp (:,:,:) 
@@ -754,5 +755,5 @@
          zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:)   ! salinity 
          avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)   ! vertical diffusive coef.
-         CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
+         CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm )
          !
          IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)
@@ -767,5 +768,5 @@
 
 
-   SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj )
+   SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj, Kbb, Kmm )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE dta_dyn_slp  ***
@@ -779,4 +780,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpi   ! zonal diapycnal slopes
       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpj   ! meridional diapycnal slopes
+      INTEGER ,                              INTENT(in ) :: Kbb, Kmm ! ocean time level indices
       !!---------------------------------------------------------------------
       !
@@ -796,5 +798,5 @@
          rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl
          CALL zdf_mxl( kt )            ! mixed layer depth
-         CALL ldf_slp( kt, rhd, rn2 )  ! slopes
+         CALL ldf_slp( kt, rhd, rn2, Kbb, Kmm )  ! slopes
          puslp (:,:,:) = uslp (:,:,:)
          pvslp (:,:,:) = vslp (:,:,:)
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90	(revision 10922)
@@ -317,5 +317,5 @@
                            CALL  istate_init    ! ocean initial state (Dynamics and tracers)
 
-                           CALL     sbc_init    ! Forcings : surface module
+                           CALL     sbc_init( Nbb, Nnn )    ! Forcings : surface module
 
       !                                      ! Tracer physics
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/nemogcm.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/nemogcm.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/nemogcm.F90	(revision 10922)
@@ -335,4 +335,9 @@
       CALL nemo_alloc()
 
+      ! Initialise time level indices
+      Nbb = 1; Nnn = 2; Naa = 3; Nrhs = Naa
+
+      ! Initialisation of temporary pointers (to be deleted after development finished)
+      CALL update_pointers( Nbb, Nnn, Naa )
       !                             !-------------------------------!
       !                             !  NEMO general initialization  !
@@ -354,5 +359,5 @@
 
       !                                      ! external forcing 
-                           CALL sbc_init        ! Forcings : surface module 
+                           CALL sbc_init( Nbb, Nnn )  ! Forcings : surface module 
 
       ! ==> clem: open boundaries init. is mandatory for sea-ice because ice BDY is not decoupled from  
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/sbcssm.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/sbcssm.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/sbcssm.F90	(revision 10922)
@@ -62,5 +62,5 @@
 CONTAINS
 
-   SUBROUTINE sbc_ssm( kt )
+   SUBROUTINE sbc_ssm( kt, Kbb, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE sbc_ssm  ***
@@ -73,4 +73,6 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
+                          ! (not needed for SAS but needed to keep a consistent interface in sbcmod.F90)
       !
       INTEGER  ::   ji, jj     ! dummy loop indices
@@ -156,5 +158,5 @@
 
 
-   SUBROUTINE sbc_ssm_init
+   SUBROUTINE sbc_ssm_init( Kbb, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE sbc_ssm_init  ***
@@ -162,4 +164,6 @@
       !! ** Purpose :   Initialisation of sea surface mean data     
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices 
+                          ! (not needed for SAS but needed to keep a consistent interface in sbcmod.F90)
       INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code
       INTEGER  :: ifpr                               ! dummy loop indice
@@ -311,5 +315,5 @@
       ENDIF
       !
-      CALL sbc_ssm( nit000 )   ! need to define ss?_m arrays used in iceistate
+      CALL sbc_ssm( nit000, Kbb, Kmm )   ! need to define ss?_m arrays used in iceistate
       l_initdone = .TRUE.
       !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/step.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/step.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/step.F90	(revision 10922)
@@ -47,5 +47,10 @@
 
    PUBLIC   stp   ! called by nemogcm.F90
+   PUBLIC   update_pointers ! called by nemo_init
 
+   !!----------------------------------------------------------------------
+   !! time level indices
+   !!----------------------------------------------------------------------
+   INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs          !! used by nemo_init
    !!----------------------------------------------------------------------
    !! NEMO/SAS 4.0 , NEMO Consortium (2018)
@@ -98,5 +103,5 @@
       IF( ln_bdy     )       CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
       ! ==>
-                             CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
+                             CALL sbc    ( kstp, Nbb, Nnn )         ! Sea Boundary Condition (including sea-ice)
 
                              CALL dia_wri( kstp )         ! ocean model: outputs
@@ -128,5 +133,5 @@
       ! Coupled mode
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( lk_oasis    )  CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges if OASIS-coupled ice
+      IF( lk_oasis    )  CALL sbc_cpl_snd( kstp, Nnn )     ! coupled mode : field exchanges if OASIS-coupled ice
 
 #if defined key_iomput
@@ -148,4 +153,42 @@
    END SUBROUTINE stp
 
+   SUBROUTINE update_pointers( Kbb, Kmm, Kaa )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE update_pointers  ***
+      !!
+      !! ** Purpose :   Associate temporary pointer arrays.
+      !!                For IMMERSE development phase only - to be deleted
+      !!
+      !! ** Method  :
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices
+
+      ub => uu(:,:,:,Kbb); un => uu(:,:,:,Kmm); ua => uu(:,:,:,Kaa)
+      vb => vv(:,:,:,Kbb); vn => vv(:,:,:,Kmm); va => vv(:,:,:,Kaa)
+      wn => ww(:,:,:)
+      hdivn => hdiv(:,:,:)
+
+      sshb =>  ssh(:,:,Kbb); sshn =>  ssh(:,:,Kmm); ssha =>  ssh(:,:,Kaa)
+      ub_b => uu_b(:,:,Kbb); un_b => uu_b(:,:,Kmm); ua_b => uu_b(:,:,Kaa)
+      vb_b => vv_b(:,:,Kbb); vn_b => vv_b(:,:,Kmm); va_b => vv_b(:,:,Kaa)
+
+      tsb => ts(:,:,:,:,Kbb); tsn => ts(:,:,:,:,Kmm); tsa => ts(:,:,:,:,Kaa)
+
+      e3t_b => e3t(:,:,:,Kbb); e3t_n => e3t(:,:,:,Kmm); e3t_a => e3t(:,:,:,Kaa)
+      e3u_b => e3u(:,:,:,Kbb); e3u_n => e3u(:,:,:,Kmm); e3u_a => e3u(:,:,:,Kaa)
+      e3v_b => e3v(:,:,:,Kbb); e3v_n => e3v(:,:,:,Kmm); e3v_a => e3v(:,:,:,Kaa)
+
+      e3f_n => e3f(:,:,:)
+
+      e3w_b  => e3w (:,:,:,Kbb); e3w_n  => e3w (:,:,:,Kmm)
+      e3uw_b => e3uw(:,:,:,Kbb); e3uw_n => e3uw(:,:,:,Kmm)
+      e3vw_b => e3vw(:,:,:,Kbb); e3vw_n => e3vw(:,:,:,Kmm)
+
+      gdept_b => gdept(:,:,:,Kbb); gdept_n => gdept(:,:,:,Kmm) 
+      gdepw_b => gdepw(:,:,:,Kbb); gdepw_n => gdepw(:,:,:,Kmm) 
+      gde3w_n => gde3w(:,:,:)
+
+   END SUBROUTINE update_pointers
+
    !!======================================================================
 END MODULE step
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 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcadv.F90	(revision 10922)
@@ -116,5 +116,5 @@
          !
          IF( ln_ldfeiv .AND. .NOT. ln_traldf_triad )   & 
-            &              CALL ldf_eiv_trp( kt, nittrc000, zun, zvn, zwn, 'TRC' )  ! add the eiv transport
+            &              CALL ldf_eiv_trp( kt, nittrc000, zun, zvn, zwn, 'TRC', Kmm )  ! add the eiv transport
          !
          IF( ln_mle    )   CALL tra_mle_trp( kt, nittrc000, zun, zvn, zwn, 'TRC' )  ! add the mle transport
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 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcldf.F90	(revision 10922)
@@ -51,5 +51,5 @@
 CONTAINS
 
-   SUBROUTINE trc_ldf( kt )
+   SUBROUTINE trc_ldf( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf  ***
@@ -59,4 +59,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kmm  ! ocean time-level index
       !
       INTEGER            :: ji, jj, jk, jn
@@ -97,7 +98,7 @@
          CALL tra_ldf_iso  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb, trb, tra, jptra,    1     )
       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     )
+         CALL tra_ldf_triad( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb, trb, tra, jptra,    1    , Kmm )
       CASE ( np_blp , np_blp_i , np_blp_it )          ! bilaplacian: all operator (iso-level, -neutral)
-         CALL tra_ldf_blp  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb     , tra, jptra, nldf_trc )
+         CALL tra_ldf_blp  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb     , tra, jptra, nldf_trc, Kmm )
       END SELECT
       !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90	(revision 10922)
@@ -73,5 +73,5 @@
          ENDIF
          !                                                      
-                                CALL trc_ldf    ( kt )      ! lateral mixing
+                                CALL trc_ldf    ( kt, Kmm ) ! lateral mixing
 #if defined key_agrif
          IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_trc       ! tracers sponge
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcstp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcstp.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcstp.F90	(revision 10922)
@@ -87,5 +87,5 @@
       IF( l_trcdm2dc )   CALL trc_mean_qsr( kt )
       !    
-      IF( nn_dttrc /= 1 )   CALL trc_sub_stp( kt )  ! averaging physical variables for sub-stepping
+      IF( nn_dttrc /= 1 )   CALL trc_sub_stp( kt, Kmm )  ! averaging physical variables for sub-stepping
       !    
       IF( MOD( kt , nn_dttrc ) == 0 ) THEN      ! only every nn_dttrc time step
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcsub.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcsub.F90	(revision 10921)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcsub.F90	(revision 10922)
@@ -82,5 +82,5 @@
 CONTAINS
 
-   SUBROUTINE trc_sub_stp( kt )
+   SUBROUTINE trc_sub_stp( kt, Kmm )
       !!-------------------------------------------------------------------
       !!                     ***  ROUTINE trc_stp  ***
@@ -92,4 +92,5 @@
       !!-------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kmm  ! ocean time-level index
       !
       INTEGER ::   ji, jj, jk   ! dummy loop indices
@@ -282,5 +283,5 @@
          ENDIF
          !
-         CALL trc_sub_ssh( kt )         ! after ssh & vertical velocity
+         CALL trc_sub_ssh( kt, Kmm )         ! after ssh & vertical velocity
          !
       ENDIF
@@ -445,5 +446,5 @@
 
 
-   SUBROUTINE trc_sub_ssh( kt ) 
+   SUBROUTINE trc_sub_ssh( kt, Kmm ) 
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE trc_sub_ssh  ***
@@ -464,4 +465,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! time step
+      INTEGER, INTENT(in) ::   Kmm  ! ocean time-level index
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -484,5 +486,5 @@
       !
 !!gm BUG here !   hdivn will include the runoff divergence at the wrong timestep !!!!
-      CALL div_hor( kt )                              ! Horizontal divergence & Relative vorticity
+      CALL div_hor( kt, Kmm )                         ! Horizontal divergence & Relative vorticity
       !
       z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
