Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfddm.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfddm.F90	(revision 10954)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfddm.F90	(revision 10955)
@@ -38,5 +38,5 @@
 CONTAINS
 
-   SUBROUTINE zdf_ddm( kt, p_avm, p_avt, p_avs )
+   SUBROUTINE zdf_ddm( kt, Kmm, p_avm, p_avt, p_avs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_ddm  ***
@@ -68,5 +68,6 @@
       !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step indexocean time step
+      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
+      INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm   !  Kz on momentum    (w-points)
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avt   !  Kz on temperature (w-points)
@@ -93,7 +94,7 @@
          DO jj = 1, jpj                !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==!
             DO ji = 1, jpi
-               zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
-!!gm please, use e3w_n below 
-                  &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) ) 
+               zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
+!!gm please, use e3w(:,:,:,Kmm) below 
+                  &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 
                !
                zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  )  &
@@ -102,6 +103,6 @@
                    &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
                !
-               zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )
-               zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) 
+               zdt = zaw * ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) )
+               zds = zbw * ( ts(ji,jj,jk-1,jp_sal,Kmm) - ts(ji,jj,jk,jp_sal,Kmm) ) 
                IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp
                zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfdrg.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfdrg.F90	(revision 10954)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfdrg.F90	(revision 10955)
@@ -81,5 +81,5 @@
 CONTAINS
 
-   SUBROUTINE zdf_drg( kt, k_mk, pCdmin, pCdmax, pz0, pke0, pCd0,   &   ! <<== in 
+   SUBROUTINE zdf_drg( kt, Kmm, k_mk, pCdmin, pCdmax, pz0, pke0, pCd0,   &   ! <<== in 
       &                                                     pCdU )      ! ==>> out : bottom drag [m/s]
       !!----------------------------------------------------------------------
@@ -99,4 +99,5 @@
       !!----------------------------------------------------------------------
       INTEGER                 , INTENT(in   ) ::   kt       ! ocean time-step index
+      INTEGER                 , INTENT(in   ) ::   Kmm      ! ocean time level index
       !                       !               !!         !==  top or bottom variables  ==!
       INTEGER , DIMENSION(:,:), INTENT(in   ) ::   k_mk     ! wet level (1st or last)
@@ -117,7 +118,7 @@
             DO ji = 2, jpim1
                imk = k_mk(ji,jj)          ! ocean bottom level at t-points
-               zut = un(ji,jj,imk) + un(ji-1,jj,imk)     ! 2 x velocity at t-point
-               zvt = vn(ji,jj,imk) + vn(ji,jj-1,imk)
-               zzz = 0.5_wp * e3t_n(ji,jj,imk)           ! altitude below/above (top/bottom) the boundary
+               zut = uu(ji,jj,imk,Kmm) + uu(ji-1,jj,imk,Kmm)     ! 2 x velocity at t-point
+               zvt = vv(ji,jj,imk,Kmm) + vv(ji,jj-1,imk,Kmm)
+               zzz = 0.5_wp * e3t(ji,jj,imk,Kmm)           ! altitude below/above (top/bottom) the boundary
                !
 !!JC: possible WAD implementation should modify line below if layers vanish
@@ -131,6 +132,6 @@
             DO ji = 2, jpim1
                imk = k_mk(ji,jj)    ! ocean bottom level at t-points
-               zut = un(ji,jj,imk) + un(ji-1,jj,imk)     ! 2 x velocity at t-point
-               zvt = vn(ji,jj,imk) + vn(ji,jj-1,imk)
+               zut = uu(ji,jj,imk,Kmm) + uu(ji-1,jj,imk,Kmm)     ! 2 x velocity at t-point
+               zvt = vv(ji,jj,imk,Kmm) + vv(ji,jj-1,imk,Kmm)
                !                                                           ! here pCd0 = mask*boost * drag
                pCdU(ji,jj) = - pCd0(ji,jj) * SQRT(  0.25 * ( zut*zut + zvt*zvt ) + pke0  )
@@ -183,6 +184,6 @@
             !
             ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)
-            zCdu = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / e3u_n(ji,jj,ikbu)
-            zCdv = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / e3v_n(ji,jj,ikbv)
+            zCdu = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / e3u(ji,jj,ikbu,Kmm)
+            zCdv = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / e3v(ji,jj,ikbv,Kmm)
             !
             pua(ji,jj,ikbu) = pua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * pub(ji,jj,ikbu)
@@ -198,6 +199,6 @@
                !
                ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)
-               zCdu = 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / e3u_n(ji,jj,ikbu)    ! NB: Cdtop masked
-               zCdv = 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / e3v_n(ji,jj,ikbv)
+               zCdu = 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / e3u(ji,jj,ikbu,Kmm)    ! NB: Cdtop masked
+               zCdv = 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / e3v(ji,jj,ikbv,Kmm)
                !
                pua(ji,jj,ikbu) = pua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * pub(ji,jj,ikbu)
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfiwm.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfiwm.F90	(revision 10954)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfiwm.F90	(revision 10955)
@@ -69,5 +69,5 @@
 
 
-   SUBROUTINE zdf_iwm( kt, p_avm, p_avt, p_avs )
+   SUBROUTINE zdf_iwm( kt, Kmm, p_avm, p_avt, p_avs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_iwm  ***
@@ -118,4 +118,5 @@
       !!----------------------------------------------------------------------
       INTEGER                    , INTENT(in   ) ::   kt             ! ocean time step
+      INTEGER                    , INTENT(in   ) ::   Kmm            ! time level index
       REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm          ! momentum Kz (w-points)
       REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avt, p_avs   ! tracer   Kz (w-points)
@@ -155,5 +156,5 @@
          END DO
       END DO
-!!gm gde3w ==>>>  check for ssh taken into account.... seem OK gde3w_n=gdept_n - sshn
+!!gm gde3w ==>>>  check for ssh taken into account.... seem OK gde3w_n=gdept(:,:,:,Kmm) - ssh(:,:,Kmm)
       DO jk = 2, jpkm1              ! complete with the level-dependent part
          DO jj = 1, jpj             
@@ -162,11 +163,11 @@
                   zemx_iwm(ji,jj,jk) = 0._wp
                ELSE
-                  zemx_iwm(ji,jj,jk) = zfact(ji,jj) * (  EXP( ( gde3w_n(ji,jj,jk  ) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) )     &
-                       &                               - EXP( ( gde3w_n(ji,jj,jk-1) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) ) )   &
-                       &                            / ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) )
+                  zemx_iwm(ji,jj,jk) = zfact(ji,jj) * (  EXP( ( gde3w(ji,jj,jk  ) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) )     &
+                       &                               - EXP( ( gde3w(ji,jj,jk-1) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) ) )   &
+                       &                            / ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) )
                ENDIF
             END DO
          END DO
-!!gm delta(gde3w_n) = e3t_n  !!  Please verify the grid-point position w versus t-point
+!!gm delta(gde3w) = e3t(:,:,:,Kmm)  !!  Please verify the grid-point position w versus t-point
 !!gm it seems to me that only 1/hcri_iwm  is used ==>  compute it one for all
 
@@ -182,5 +183,5 @@
          zfact(:,:) = 0._wp
          DO jk = 2, jpkm1              ! part independent of the level
-            zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
+            zfact(:,:) = zfact(:,:) + e3w(:,:,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
          END DO
          !
@@ -199,5 +200,5 @@
          zfact(:,:) = 0._wp
          DO jk = 2, jpkm1              ! part independent of the level
-            zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk)
+            zfact(:,:) = zfact(:,:) + e3w(:,:,jk,Kmm) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk)
          END DO
          !
@@ -220,10 +221,10 @@
       zfact(:,:)   = 0._wp
       DO jk = 2, jpkm1
-         zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
+         zfact(:,:) = zfact(:,:) + e3w(:,:,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
          zwkb(:,:,jk) = zfact(:,:)
       END DO
 !!gm even better:
 !      DO jk = 2, jpkm1
-!         zwkb(:,:) = zwkb(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  )
+!         zwkb(:,:) = zwkb(:,:) + e3w(:,:,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  )
 !      END DO
 !      zfact(:,:) = zwkb(:,:,jpkm1)
@@ -267,12 +268,12 @@
       DO jk = 2, jpkm1              ! complete with the level-dependent part
          zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   &
-            &                                / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) )
-!!gm  use of e3t_n just above?
+            &                                / ( gde3w(:,:,jk) - gde3w(:,:,jk-1) )
+!!gm  use of e3t(:,:,:,Kmm) just above?
       END DO
       !
 !!gm  this is to be replaced by just a constant value znu=1.e-6 m2/s
       ! Calculate molecular kinematic viscosity
-      znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(:,:,:,jp_tem) + 0.00694_wp * tsn(:,:,:,jp_tem) * tsn(:,:,:,jp_tem)  &
-         &                                  + 0.02305_wp * tsn(:,:,:,jp_sal)  ) * tmask(:,:,:) * r1_rau0
+      znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * ts(:,:,:,jp_tem,Kmm) + 0.00694_wp * ts(:,:,:,jp_tem,Kmm) * ts(:,:,:,jp_tem,Kmm)  &
+         &                                  + 0.02305_wp * ts(:,:,:,jp_sal,Kmm)  ) * tmask(:,:,:) * r1_rau0
       DO jk = 2, jpkm1
          znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk)
@@ -314,5 +315,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zztmp = zztmp + e3w_n(ji,jj,jk) * e1e2t(ji,jj)   &
+                  zztmp = zztmp + e3w(ji,jj,jk,Kmm) * e1e2t(ji,jj)   &
                      &          * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
                END DO
@@ -374,5 +375,5 @@
          z2d(:,:) = 0._wp
          DO jk = 2, jpkm1
-            z2d(:,:) = z2d(:,:) + e3w_n(:,:,jk) * z3d(:,:,jk) * wmask(:,:,jk)
+            z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)
          END DO
          z2d(:,:) = rau0 * z2d(:,:)
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfmxl.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfmxl.F90	(revision 10954)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfmxl.F90	(revision 10955)
@@ -57,5 +57,5 @@
 
 
-   SUBROUTINE zdf_mxl( kt )
+   SUBROUTINE zdf_mxl( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdfmxl  ***
@@ -75,4 +75,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
       !
       INTEGER  ::   ji, jj, jk      ! dummy loop indices
@@ -98,5 +99,5 @@
             DO ji = 1, jpi
                ikt = mbkt(ji,jj)
-               hmlp(ji,jj) = hmlp(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk)
+               hmlp(ji,jj) = hmlp(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm)
                IF( hmlp(ji,jj) < zN2_c )   nmln(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
             END DO
@@ -118,7 +119,7 @@
             iiki = imld(ji,jj)
             iikn = nmln(ji,jj)
-            hmld (ji,jj) = gdepw_n(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth 
-            hmlp (ji,jj) = gdepw_n(ji,jj,iikn  ) * ssmask(ji,jj)    ! Mixed layer depth
-            hmlpt(ji,jj) = gdept_n(ji,jj,iikn-1) * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer
+            hmld (ji,jj) = gdepw(ji,jj,iiki  ,Kmm) * ssmask(ji,jj)    ! Turbocline depth 
+            hmlp (ji,jj) = gdepw(ji,jj,iikn  ,Kmm) * ssmask(ji,jj)    ! Mixed layer depth
+            hmlpt(ji,jj) = gdept(ji,jj,iikn-1,Kmm) * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer
          END DO
       END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfosm.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfosm.F90	(revision 10954)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfosm.F90	(revision 10955)
@@ -941,5 +941,5 @@
  ! Stable conditions
                DO jk = 2, ibld(ji,jj)
-                  zznd_d=gdepw_n(ji,jj,jk) / dstokes(ji,jj)
+                  zznd_d=gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) &
                        &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
@@ -1066,5 +1066,5 @@
             IF ( lconv(ji,jj) ) THEN
                DO jk = 2, imld(ji,jj)
-                  zznd_ml=gdepw_n(ji,jj,jk) / zhml(ji,jj)
+                  zznd_ml=gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj)                &
                        &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      &
@@ -1188,5 +1188,5 @@
              IF ( lconv(ji,jj) ) THEN
                DO jk = 1, imld(ji,jj) - 1
-                  znd=gdepw_n(ji,jj,jk) / zhml(ji,jj)
+                  znd=gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd
                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfphy.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfphy.F90	(revision 10954)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfphy.F90	(revision 10955)
@@ -192,7 +192,7 @@
       ioptio = 0 
       IF( ln_zdfcst ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_CST   ;   ENDIF
-      IF( ln_zdfric ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_RIC   ;   CALL zdf_ric_init   ;   ENDIF
-      IF( ln_zdftke ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_TKE   ;   CALL zdf_tke_init   ;   ENDIF
-      IF( ln_zdfgls ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_GLS   ;   CALL zdf_gls_init   ;   ENDIF
+      IF( ln_zdfric ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_RIC   ;   CALL zdf_ric_init          ;   ENDIF
+      IF( ln_zdftke ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_TKE   ;   CALL zdf_tke_init( Kmm )   ;   ENDIF
+      IF( ln_zdfgls ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_GLS   ;   CALL zdf_gls_init          ;   ENDIF
       IF( ln_zdfosm ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_OSM   ;   CALL zdf_osm_init( Kmm )   ;   ENDIF
       !
@@ -244,9 +244,9 @@
          !
          !                       !* bottom drag
-         CALL zdf_drg( kt, mbkt    , r_Cdmin_bot, r_Cdmax_bot,   &   ! <<== in 
+         CALL zdf_drg( kt, Kmm, mbkt , r_Cdmin_bot, r_Cdmax_bot,   &   ! <<== in 
             &              r_z0_bot,   r_ke0_bot,    rCd0_bot,   &
             &                                        rCdU_bot  )     ! ==>> out : bottom drag [m/s]
          IF( ln_isfcav ) THEN    !* top drag   (ocean cavities)
-            CALL zdf_drg( kt, mikt    , r_Cdmin_top, r_Cdmax_top,   &   ! <<== in 
+            CALL zdf_drg( kt, Kmm, mikt , r_Cdmin_top, r_Cdmax_top,   &   ! <<== in 
                &              r_z0_top,   r_ke0_top,    rCd0_top,   &
                &                                        rCdU_top  )     ! ==>> out : bottom drag [m/s]
@@ -287,5 +287,5 @@
       !                                         !* double diffusive mixing
       IF( ln_zdfddm ) THEN                            ! update avt and compute avs
-                        CALL zdf_ddm( kt, avm, avt, avs )
+                        CALL zdf_ddm( kt, Kmm,  avm, avt, avs )
       ELSE                                            ! same mixing on all tracers
          avs(2:jpim1,2:jpjm1,1:jpkm1) = avt(2:jpim1,2:jpjm1,1:jpkm1)
@@ -293,6 +293,6 @@
       !
       !                                         !* wave-induced mixing 
-      IF( ln_zdfswm )   CALL zdf_swm( kt, avm, avt, avs )   ! surface  wave (Qiao et al. 2004) 
-      IF( ln_zdfiwm )   CALL zdf_iwm( kt, avm, avt, avs )   ! internal wave (de Lavergne et al 2017)
+      IF( ln_zdfswm )   CALL zdf_swm( kt, Kmm, avm, avt, avs )   ! surface  wave (Qiao et al. 2004) 
+      IF( ln_zdfiwm )   CALL zdf_iwm( kt, Kmm, avm, avt, avs )   ! internal wave (de Lavergne et al 2017)
 
 #if defined key_agrif 
@@ -315,5 +315,5 @@
       ENDIF
       !
-      CALL zdf_mxl( kt )                        !* mixed layer depth, and level
+      CALL zdf_mxl( kt, Kmm )                        !* mixed layer depth, and level
       !
       IF( lrst_oce ) THEN                       !* write TKE, GLS or RIC fields in the restart file
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfswm.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfswm.F90	(revision 10954)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfswm.F90	(revision 10955)
@@ -34,5 +34,5 @@
 CONTAINS
 
-   SUBROUTINE zdf_swm( kt, p_avm, p_avt, p_avs )
+   SUBROUTINE zdf_swm( kt, Kmm, p_avm, p_avt, p_avs )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE zdf_swm ***
@@ -52,4 +52,5 @@
       !!---------------------------------------------------------------------
       INTEGER                    , INTENT(in   ) ::   kt             ! ocean time step
+      INTEGER                    , INTENT(in   ) ::   Kmm            ! time level index
       REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm          ! momentum Kz (w-points)
       REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avt, p_avs   ! tracer   Kz (w-points)
@@ -63,5 +64,5 @@
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               zqb = zcoef * hsw(ji,jj) * tsd2d(ji,jj) * EXP( -3. * wnum(ji,jj) * gdepw_n(ji,jj,jk) ) * wmask(ji,jj,jk)
+               zqb = zcoef * hsw(ji,jj) * tsd2d(ji,jj) * EXP( -3. * wnum(ji,jj) * gdepw(ji,jj,jk,Kmm) ) * wmask(ji,jj,jk)
                !
                p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zqb
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdftke.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdftke.F90	(revision 10954)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdftke.F90	(revision 10955)
@@ -630,5 +630,5 @@
 
 
-   SUBROUTINE zdf_tke_init
+   SUBROUTINE zdf_tke_init( Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_tke_init  ***
@@ -646,6 +646,7 @@
       USE zdf_oce , ONLY : ln_zdfiwm   ! Internal Wave Mixing flag
       !!
-      INTEGER ::   ji, jj, jk   ! dummy loop indices
-      INTEGER ::   ios
+      INTEGER, INTENT(in) ::   Kmm          ! time level index
+      INTEGER             ::   ji, jj, jk   ! dummy loop indices
+      INTEGER             ::   ios
       !!
       NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,          &
@@ -724,5 +725,5 @@
       ENDIF
       
-      IF( nn_etau == 2  )   CALL zdf_mxl( nit000 )      ! Initialization of nmln 
+      IF( nn_etau == 2  )   CALL zdf_mxl( nit000, Kmm )      ! Initialization of nmln 
 
       !                               !* depth of penetration of surface tke
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 10954)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/dtadyn.F90	(revision 10955)
@@ -162,6 +162,6 @@
       CALL bn2    ( ts(:,:,:,:,Kmm), rab_n, rn2, Kmm )  ! before Brunt-Vaisala frequency need for zdfmxl
 
-      rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl
-      CALL zdf_mxl( kt )                                                   ! In any case, we need mxl
+      rn2b(:,:,:) = rn2(:,:,:)         ! needed for zdfmxl
+      CALL zdf_mxl( kt, Kmm )          ! In any case, we need mxl
       !
       hmld(:,:)       = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1)    ! mixed layer depht
@@ -796,6 +796,6 @@
          &                                        rhd, gru , grv , grui, grvi )  ! of t, s, rd at the first ocean level
 
-         rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl
-         CALL zdf_mxl( kt )            ! mixed layer depth
+         rn2b(:,:,:) = rn2(:,:,:)                ! needed for zdfmxl
+         CALL zdf_mxl( kt, Kmm )                 ! mixed layer depth
          CALL ldf_slp( kt, rhd, rn2, Kbb, Kmm )  ! slopes
          puslp (:,:,:) = uslp (:,:,:)
