Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfgls.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfgls.F90	(revision 10880)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfgls.F90	(revision 10883)
@@ -124,5 +124,5 @@
 
 
-   SUBROUTINE zdf_gls( kt, p_sh2, p_avm, p_avt )
+   SUBROUTINE zdf_gls( kt, Kbb, Kmm, p_sh2, p_avm, p_avt )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE zdf_gls  ***
@@ -134,4 +134,5 @@
       !!
       INTEGER                   , INTENT(in   ) ::   kt             ! ocean time step
+      INTEGER                   , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points)
@@ -176,6 +177,6 @@
           zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
           zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0)
-          ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2  &
-             &                                         + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  )
+          ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2  &
+             &                                         + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2  )
          END DO
       END DO
@@ -185,6 +186,6 @@
                zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
                zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )     ! (CAUTION: CdU<0)
-               ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2  &
-                  &                                         + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  )
+               ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )**2  &
+                  &                                         + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2  )
             END DO
          END DO
@@ -222,6 +223,6 @@
             DO jj = 2, jpjm1 
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zup   = hmxl_n(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1)
-                  zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) )
+                  zup   = hmxl_n(ji,jj,jk) * gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)
+                  zdown = vkarmn * gdepw(ji,jj,jk,Kmm) * ( -gdepw(ji,jj,jk,Kmm) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) )
                   zcoef = ( zup / MAX( zdown, rsmall ) )
                   zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk)
@@ -276,7 +277,7 @@
                zcof = rfact_tke * tmask(ji,jj,jk)
                !                                        ! lower diagonal, in fact not used for jk = 2 (see surface conditions)
-               zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )
+               zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) )
                !                                        ! upper diagonal, in fact not used for jk = ibotm1 (see bottom conditions)
-               zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) )
+               zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk,Kmm) )
                !                                        ! diagonal
                zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk)  + rdt * zdiss * wmask(ji,jj,jk) 
@@ -306,5 +307,5 @@
       ! 
       ! One level below
-      en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2))   &
+      en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw(:,:,2,Kmm))   &
          &                 / zhsro(:,:) )**(1.5_wp*ra_sf)  )**(2._wp/3._wp)                      , rn_emin   )
       zd_lw(:,:,2) = 0._wp 
@@ -325,9 +326,9 @@
       zdiag(:,:,2) = zdiag(:,:,2) +  zd_lw(:,:,2) ! Remove zd_lw from zdiag
       zd_lw(:,:,2) = 0._wp
-      zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) ))
+      zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:)) ))
       zflxs(:,:)   = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) &
-          &                    * (  ( zhsro(:,:)+gdept_n(:,:,1) ) / zhsro(:,:)  )**(1.5_wp*ra_sf)
-!!gm why not   :                        * ( 1._wp + gdept_n(:,:,1) / zhsro(:,:) )**(1.5_wp*ra_sf)
-      en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w_n(:,:,2)
+          &                    * (  ( zhsro(:,:)+gdept(:,:,1,Kmm) ) / zhsro(:,:)  )**(1.5_wp*ra_sf)
+!!gm why not   :                        * ( 1._wp + gdept(:,:,1,Kmm) / zhsro(:,:) )**(1.5_wp*ra_sf)
+      en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm)
       !
       !
@@ -526,7 +527,7 @@
                zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
                !                                               ! lower diagonal
-               zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )
+               zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) )
                !                                               ! upper diagonal
-               zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) )
+               zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk,Kmm) )
                !                                               ! diagonal
                zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk)
@@ -554,6 +555,6 @@
          !
          ! One level below
-         zkar    (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) )))
-         zdep    (:,:)   = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:)
+         zkar    (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(:,:,2,Kmm)/zhsro(:,:) )))
+         zdep    (:,:)   = (zhsro(:,:) + gdepw(:,:,2,Kmm)) * zkar(:,:)
          psi     (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
          zd_lw(:,:,2) = 0._wp
@@ -575,11 +576,11 @@
          !
          ! Set psi vertical flux at the surface:
-         zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope
-         zdep (:,:)   = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)
+         zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:) )) ! Lengh scale slope
+         zdep (:,:)   = ((zhsro(:,:) + gdept(:,:,1,Kmm)) / zhsro(:,:))**(rmm*ra_sf)
          zflxs(:,:)   = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp)
          zdep (:,:)   = rsbc_psi1 * (zwall_psi(:,:,1)*p_avm(:,:,1)+zwall_psi(:,:,2)*p_avm(:,:,2)) * &
-            &           ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.)
+            &           ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept(:,:,1,Kmm))**(rnn-1.)
          zflxs(:,:)   = zdep(:,:) * zflxs(:,:)
-         psi  (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2)
+         psi  (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm)
          !
       END SELECT
@@ -607,5 +608,5 @@
                !
                ! Just above last level, Dirichlet condition again (GOTM like)
-               zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t_n(ji,jj,ibotm1) )
+               zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t(ji,jj,ibotm1,Kmm) )
                psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn
                zd_lw(ji,jj,ibotm1) = 0._wp
@@ -635,8 +636,8 @@
                !
                ! Set psi vertical flux at the bottom:
-               zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t_n(ji,jj,ibotm1)
+               zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t(ji,jj,ibotm1,Kmm)
                zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) )   &
                   &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp)
-               psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w_n(ji,jj,ibotm1)
+               psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w(ji,jj,ibotm1,Kmm)
             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 10880)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfosm.F90	(revision 10883)
@@ -42,5 +42,5 @@
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and active tracers
-                      ! uses wn from previous time step (which is now wb) to calculate hbl
+                      ! uses ww from previous time step (which is now wb) to calculate hbl
    USE dom_oce        ! ocean space and time domain
    USE zdf_oce        ! ocean vertical physics
@@ -122,5 +122,5 @@
 
 
-   SUBROUTINE zdf_osm( kt, p_avm, p_avt )
+   SUBROUTINE zdf_osm( kt, Kbb, Kmm, p_avm, p_avt )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE zdf_osm  ***
@@ -157,5 +157,6 @@
       !!         the equation number. (LMD94, here after)
       !!----------------------------------------------------------------------
-      INTEGER                   , INTENT(in   ) ::   kt            ! ocean time step
+      INTEGER                   , INTENT(in   ) ::  kt             ! ocean time step
+      INTEGER                   , INTENT(in   ) ::  Kbb, Kmm       ! ocean time level indices
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::  p_avm, p_avt   ! momentum and tracer Kz (w-points)
       !!
@@ -314,5 +315,5 @@
            zwth0(ji,jj) = - qns(ji,jj) * r1_rau0_rcp * tmask(ji,jj,1)
            ! Upwards surface salinity flux for non-local term
-           zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal)  + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1)
+           zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm)  + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1)
            ! Non radiative upwards surface buoyancy flux
            zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) -  grav * zbeta * zws0(ji,jj)
@@ -407,10 +408,10 @@
      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      ! BL must be always 2 levels deep.
-      hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,3) )
+      hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,3,Kmm) )
       ibld(:,:) = 3
       DO jk = 4, jpkm1
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               IF ( hbl(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN
+               IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN
                   ibld(ji,jj) = MIN(mbkt(ji,jj), jk)
                ENDIF
@@ -430,12 +431,12 @@
                zthick=0._wp
                DO jm = 2, ibld(ji,jj)
-                  zthick=zthick+e3t_n(ji,jj,jm)
-                  zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem)
-                  zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal)
-                  zu   = zu  + e3t_n(ji,jj,jm) &
-                     &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) &
+                  zthick=zthick+e3t(ji,jj,jm,Kmm)
+                  zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm)
+                  zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm)
+                  zu   = zu  + e3t(ji,jj,jm,Kmm) &
+                     &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) &
                      &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
-                  zv   = zv  + e3t_n(ji,jj,jm) &
-                     &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) &
+                  zv   = zv  + e3t(ji,jj,jm,Kmm) &
+                     &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) &
                      &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
                END DO
@@ -444,9 +445,9 @@
                zu_bl(ji,jj) = zu / zthick
                zv_bl(ji,jj) = zv / zthick
-               zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem)
-               zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal)
-               zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) &
+               zdt_bl(ji,jj) = zt_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm)
+               zds_bl(ji,jj) = zs_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm)
+               zdu_bl(ji,jj) = zu_bl(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) &
                      &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
-               zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) &
+               zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) &
                      &   / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) )
                zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj)
@@ -487,12 +488,12 @@
       ibld(:,:) = 3
 
-      zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it
+      zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - ww(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it
       zhbl_t(:,:) = MIN(zhbl_t(:,:), ht_n(:,:))
-      zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + wn(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom
+      zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + ww(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom
 
       DO jk = 4, jpkm1
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN
+               IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN
                   ibld(ji,jj) =  MIN(mbkt(ji,jj), jk)
                ENDIF
@@ -520,11 +521,11 @@
 
                   DO jk = imld(ji,jj), ibld(ji,jj)
-                     zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) &
-                          & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) + zvel_max
-
-                     zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w_n(ji,jj,jk) )
+                     zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) &
+                          & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) + zvel_max
+
+                     zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w(ji,jj,jk,Kmm) )
                      zhbl_s = MIN(zhbl_s, ht_n(ji,jj))
 
-                     IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1
+                     IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1
                   END DO
                   hbl(ji,jj) = zhbl_s
@@ -534,6 +535,6 @@
 ! stable
                   DO jk = imld(ji,jj), ibld(ji,jj)
-                     zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )          &
-                          &               - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) &
+                     zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) )          &
+                          &               - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) &
                           & + 2.0 * zwstrl(ji,jj)**2 / zhbl_s
 
@@ -543,10 +544,10 @@
                           &               + ( ( 0.32 / 3.0 )           * EXP( -  2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) )   &
                           &               -   ( 0.32 / 3.0  - 0.0485 ) * EXP( - 12.5 * ( hbli(ji,jj) / zhbl_s      ) ) ) &
-                          &          * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w_n(ji,jj,jk) / zdhdt(ji,jj)  ! ALMG to investigate whether need to include wn here
+                          &          * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w(ji,jj,jk,Kmm) / zdhdt(ji,jj)  ! ALMG to investigate whether need to include ww here
 
                      zhbl_s = MIN(zhbl_s, ht_n(ji,jj))
-                     IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1
+                     IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1
                   END DO
-                  hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,3) )
+                  hbl(ji,jj) = MAX(zhbl_s, gdepw(ji,jj,3,Kmm) )
                   ibld(ji,jj) = MAX(jm, 3 )
                   IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj)
@@ -558,9 +559,9 @@
                   hbli(ji,jj) = hbl(ji,jj)
                ELSE
-                  hbl(ji,jj) = MAX(hbl(ji,jj), gdepw_n(ji,jj,3) )
+                  hbl(ji,jj) = MAX(hbl(ji,jj), gdepw(ji,jj,3,Kmm) )
                   IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj)
                ENDIF
             ENDIF
-            zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
+            zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm)
          END DO
       END DO
@@ -581,12 +582,12 @@
                zthick=0._wp
                DO jm = 2, ibld(ji,jj)
-                  zthick=zthick+e3t_n(ji,jj,jm)
-                  zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem)
-                  zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal)
-                  zu   = zu  + e3t_n(ji,jj,jm) &
-                     &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) &
+                  zthick=zthick+e3t(ji,jj,jm,Kmm)
+                  zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm)
+                  zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm)
+                  zu   = zu  + e3t(ji,jj,jm,Kmm) &
+                     &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) &
                      &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
-                  zv   = zv  + e3t_n(ji,jj,jm) &
-                     &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) &
+                  zv   = zv  + e3t(ji,jj,jm,Kmm) &
+                     &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) &
                      &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
                END DO
@@ -595,12 +596,12 @@
                zu_bl(ji,jj) = zu / zthick
                zv_bl(ji,jj) = zv / zthick
-               zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem)
-               zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal)
-               zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) &
+               zdt_bl(ji,jj) = zt_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm)
+               zds_bl(ji,jj) = zs_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm)
+               zdu_bl(ji,jj) = zu_bl(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) &
                       &   / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
-               zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) &
+               zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) &
                       &  / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) )
                zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj)
-               zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
+               zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm)
                IF ( lconv(ji,jj) ) THEN
                   IF ( zdb_bl(ji,jj) > 0._wp )THEN
@@ -616,11 +617,11 @@
                         zwb_ent(ji,jj) = 0._wp
                      ENDIF
-                     inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 )
+                     inhml = MAX( INT( zari * zhbl(ji,jj) / e3t(ji,jj,ibld(ji,jj),Kmm) ) , 1 )
                      imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1)
-                     zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
+                     zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)
                      zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
                   ELSE  ! IF (zdb_bl)
                      imld(ji,jj) = ibld(ji,jj) - 1
-                     zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
+                     zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)
                      zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
                   ENDIF
@@ -632,11 +633,11 @@
                         zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
                           & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 )
-                        inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 )
+                        inhml = MAX( INT( zari * zhbl(ji,jj) / e3t(ji,jj,ibld(ji,jj),Kmm) ) , 1 )
                         imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1)
-                        zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
+                        zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)
                         zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
                      ELSE
                         imld(ji,jj) = ibld(ji,jj) - 1
-                        zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
+                        zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)
                         zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
                      ENDIF ! IF (zdb_bl > 0.0)
@@ -665,12 +666,12 @@
                zthick=0._wp
                DO jm = 2, imld(ji,jj)
-                  zthick=zthick+e3t_n(ji,jj,jm)
-                  zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem)
-                  zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal)
-                  zu   = zu  + e3t_n(ji,jj,jm) &
-                     &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) &
+                  zthick=zthick+e3t(ji,jj,jm,Kmm)
+                  zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm)
+                  zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm)
+                  zu   = zu  + e3t(ji,jj,jm,Kmm) &
+                     &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) &
                      &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
-                  zv   = zv  + e3t_n(ji,jj,jm) &
-                     &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) &
+                  zv   = zv  + e3t(ji,jj,jm,Kmm) &
+                     &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) &
                      &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
                END DO
@@ -679,9 +680,9 @@
                zu_ml(ji,jj) = zu / zthick
                zv_ml(ji,jj) = zv / zthick
-               zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem)
-               zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal)
-               zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) &
+               zdt_ml(ji,jj) = zt_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm)
+               zds_ml(ji,jj) = zs_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm)
+               zdu_ml(ji,jj) = zu_ml(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) &
                      &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
-               zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) &
+               zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) &
                      &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) )
                zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj)
@@ -696,12 +697,12 @@
                   zthick=0._wp
                   DO jm = 2, imld(ji,jj)
-                     zthick=zthick+e3t_n(ji,jj,jm)
-                     zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem)
-                     zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal)
-                     zu   = zu  + e3t_n(ji,jj,jm) &
-                        &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) &
+                     zthick=zthick+e3t(ji,jj,jm,Kmm)
+                     zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm)
+                     zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm)
+                     zu   = zu  + e3t(ji,jj,jm,Kmm) &
+                        &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) &
                         &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
-                     zv   = zv  + e3t_n(ji,jj,jm) &
-                        &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) &
+                     zv   = zv  + e3t(ji,jj,jm,Kmm) &
+                        &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) &
                         &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
                   END DO
@@ -710,9 +711,9 @@
                   zu_ml(ji,jj) = zu / zthick
                   zv_ml(ji,jj) = zv / zthick
-                  zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem)
-                  zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal)
-                  zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) &
+                  zdt_ml(ji,jj) = zt_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm)
+                  zds_ml(ji,jj) = zs_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm)
+                  zdu_ml(ji,jj) = zu_ml(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) &
                         &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
-                  zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) &
+                  zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) &
                         &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) )
                   zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj)
@@ -775,5 +776,5 @@
                    zbgrad = ( zdb_ml(ji,jj) / zdh(ji,jj) )
                    DO jk = 2 , ibld(ji,jj)
-                      znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
+                      znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj)
                       zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
                       zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
@@ -791,5 +792,5 @@
                          zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj)
                          DO jk = 2, ibld(ji,jj)
-                            znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
+                            znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
                             zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
                             zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
@@ -801,5 +802,5 @@
                          zbgrad = zdb_bl(ji,jj) / zdh(ji,jj)
                          DO jk = 2, ibld(ji,jj)
-                            znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
+                            znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj)
                             zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
                             zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
@@ -824,5 +825,5 @@
               & ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
                 DO jk = 2 , ibld(ji,jj)-1
-                   znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
+                   znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj)
                    zdudz_pyc(ji,jj,jk) =  zugrad * EXP( -1.75 * ( znd + 0.75 )**2 )
                    zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
@@ -833,5 +834,5 @@
                 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj)
                 DO jk = 2, ibld(ji,jj)
-                   znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
+                   znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
                    IF ( znd < 1.0 ) THEN
                       zdudz_pyc(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 )
@@ -880,5 +881,5 @@
              IF ( lconv(ji,jj) ) THEN
                 DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity
-                    zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
+                    zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
                     !
                     zdiffut(ji,jj,jk) = 0.8   * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml    )**1.5
@@ -890,5 +891,5 @@
                 IF ( zdh(ji,jj) > 0._wp ) THEN
                    DO jk = imld(ji,jj)+1 , ibld(ji,jj)
-                       zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
+                       zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj)
                        !
                        zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 + zznd_pyc )
@@ -897,6 +898,6 @@
                    END DO
                 ENDIF
-                ! Temporay fix to ensure zdiffut is +ve; won't be necessary with wn taken out
-                zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t_n(ji,jj,ibld(ji,jj))
+                ! Temporay fix to ensure zdiffut is +ve; won't be necessary with ww taken out
+                zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t(ji,jj,ibld(ji,jj),Kmm)
                 ! could be taken out, take account of entrainment represents as a diffusivity
                 ! should remove w from here, represents entrainment
@@ -904,5 +905,5 @@
              ! stable conditions
                 DO jk = 2, ibld(ji,jj)
-                   zznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
+                   zznd_ml = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
                    zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5
                    zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 )
@@ -932,5 +933,5 @@
             IF ( lconv(ji,jj) ) THEN
               DO jk = 2, imld(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.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
                  !
@@ -967,5 +968,5 @@
              IF ( lconv(ji,jj) ) THEN
                 DO jk = 2, imld(ji,jj)
-                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
+                   zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
                    ghamu(ji,jj,jk) = ghamu(ji,jj,jk) +      ( -0.05 * EXP ( -0.4 * zznd_d )   * zsc_uw_1(ji,jj)   &
                         &          +                        0.00125 * EXP (      - zznd_d )   * zsc_uw_2(ji,jj) ) &
@@ -978,5 +979,5 @@
 ! Stable conditions
                 DO jk = 2, ibld(ji,jj) ! corrected to ibld
-                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
+                   zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
                    ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 *   1.3 * EXP ( -0.5 * zznd_d )                       &
                         &                                   * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj)
@@ -1001,5 +1002,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)
                    ! calculate turbulent length scale
                    zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           &
@@ -1035,5 +1036,5 @@
              IF ( lconv(ji,jj) ) THEN
                 DO jk = 2 , imld(ji,jj)
-                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
+                   zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
                    ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) +   0.125 * EXP( -0.5 * zznd_d )     &
                         &                                                            * (   1.0 - EXP( -0.5 * zznd_d ) )   &
@@ -1078,6 +1079,6 @@
             ELSE
                DO jk = 2, ibld(ji,jj)
-                  zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
-                  znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
+                  zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
+                  znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
                &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj)
@@ -1104,6 +1105,6 @@
              IF ( lconv(ji,jj) ) THEN
                DO jk = 2, imld(ji,jj)
-                  zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
-                  zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
+                  zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
+                  zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
                        & + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml**4 ) - EXP ( -8.0 * zznd_ml ) ) * zsc_uw_1(ji,jj)
@@ -1114,6 +1115,6 @@
              ELSE
                DO jk = 2, ibld(ji,jj)
-                  znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
-                  zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
+                  znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
+                  zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
                   IF ( zznd_d <= 2.0 ) THEN
                      ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 &
@@ -1141,5 +1142,5 @@
             IF ( lconv(ji,jj) ) THEN
                DO jk = 2, ibld(ji,jj)
-                  znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about
+                  znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about
                   IF ( znd >= 0.0 ) THEN
                      ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) )
@@ -1152,5 +1153,5 @@
             ELSE
                DO jk = 2, ibld(ji,jj)
-                  znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about
+                  znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about
                   IF ( znd >= 0.0 ) THEN
                      ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
@@ -1171,5 +1172,5 @@
           DO ji = 2, jpim1
              DO jk= 2, ibld(ji,jj)
-                znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
+                znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk)
                 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk)
@@ -1194,5 +1195,5 @@
                END DO
                DO jk = imld(ji,jj), ibld(ji,jj)
-                  znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
+                  znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj)
                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd )
                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd )
@@ -1239,10 +1240,10 @@
              DO jj = 1, jpjm1
                 DO ji = 1, jpim1   ! vector opt.
-                   z3du(ji,jj,jk) = 0.5 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )   &
-                        &                 * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) * wumask(ji,jj,jk) &
-                        &                 / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
-                   z3dv(ji,jj,jk) = 0.5 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )   &
-                        &                 * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) * wvmask(ji,jj,jk) &
-                        &                 / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
+                   z3du(ji,jj,jk) = 0.5 * (  uu(ji,jj,jk-1,Kmm) -  uu(ji  ,jj,jk,Kmm) )   &
+                        &                 * (  uu(ji,jj,jk-1,Kbb) -  uu(ji  ,jj,jk,Kbb) ) * wumask(ji,jj,jk) &
+                        &                 / (  e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) )
+                   z3dv(ji,jj,jk) = 0.5 * (  vv(ji,jj,jk-1,Kmm) -  vv(ji,jj  ,jk,Kmm) )   &
+                        &                 * (  vv(ji,jj,jk-1,Kbb) -  vv(ji,jj  ,jk,Kbb) ) * wvmask(ji,jj,jk) &
+                        &                 / (  e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) )
                 END DO
              END DO
@@ -1364,5 +1365,5 @@
 
 
-   SUBROUTINE zdf_osm_init
+   SUBROUTINE zdf_osm_init( Kmm ) 
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE zdf_osm_init  ***
@@ -1376,4 +1377,6 @@
      !! ** input   :   Namlist namosm
      !!----------------------------------------------------------------------
+     INTEGER, INTENT(in)    :: Kmm ! time level index (middle)
+     !
      INTEGER  ::   ios            ! local integer
      INTEGER  ::   ji, jj, jk     ! dummy loop indices
@@ -1423,5 +1426,5 @@
      IF( zdf_osm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' )
 
-     call osm_rst( nit000, 'READ' ) !* read or initialize hbl
+     call osm_rst( nit000, Kmm, 'READ' ) !* read or initialize hbl
 
      IF( ln_zdfddm) THEN
@@ -1517,5 +1520,5 @@
 
 
-   SUBROUTINE osm_rst( kt, cdrw )
+   SUBROUTINE osm_rst( kt, Kmm, cdrw )
      !!---------------------------------------------------------------------
      !!                   ***  ROUTINE osm_rst  ***
@@ -1527,5 +1530,6 @@
      !!----------------------------------------------------------------------
 
-     INTEGER, INTENT(in) :: kt
+     INTEGER         , INTENT(in) ::   kt     ! ocean time step index
+     INTEGER         , INTENT(in) ::   Kmm    ! ocean time level index (middle)
      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
 
@@ -1545,9 +1549,9 @@
         id1 = iom_varid( numror, 'wn'   , ldstop = .FALSE. )
         IF( id1 > 0 ) THEN                       ! 'wn' exists; read
-           CALL iom_get( numror, jpdom_autoglo, 'wn', wn, ldxios = lrxios )
-           WRITE(numout,*) ' ===>>>> :  wn read from restart file'
+           CALL iom_get( numror, jpdom_autoglo, 'wn', ww, ldxios = lrxios )
+           WRITE(numout,*) ' ===>>>> :  ww read from restart file'
         ELSE
-           wn(:,:,:) = 0._wp
-           WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially'
+           ww(:,:,:) = 0._wp
+           WRITE(numout,*) ' ===>>>> :  ww not in restart file, set to zero initially'
         END IF
         id1 = iom_varid( numror, 'hbl'   , ldstop = .FALSE. )
@@ -1568,5 +1572,5 @@
      IF( TRIM(cdrw) == 'WRITE') THEN     !* Write hbli into the restart file, then return
         IF(lwp) WRITE(numout,*) '---- osm-rst ----'
-         CALL iom_rstput( kt, nitrst, numrow, 'wn'     , wn  , ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'wn'     , ww  , ldxios = lwxios )
          CALL iom_rstput( kt, nitrst, numrow, 'hbl'    , hbl , ldxios = lwxios )
          CALL iom_rstput( kt, nitrst, numrow, 'hbli'   , hbli, ldxios = lwxios )
@@ -1580,6 +1584,6 @@
      ALLOCATE( imld_rst(jpi,jpj) )
      ! w-level of the mixing and mixed layers
-     CALL eos_rab( tsn, rab_n )
-     CALL bn2(tsn, rab_n, rn2)
+     CALL eos_rab( ts(:,:,:,:,Kmm), rab_n )
+     CALL bn2(ts(:,:,:,:,Kmm), rab_n, rn2)
      imld_rst(:,:)  = nlb10         ! Initialization to the number of w ocean point
      hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
@@ -1591,5 +1595,5 @@
            DO ji = 1, jpi
               ikt = mbkt(ji,jj)
-              hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk)
+              hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm)
               IF( hbl(ji,jj) < zN2_c )   imld_rst(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
            END DO
@@ -1600,5 +1604,5 @@
         DO ji = 1, jpi
            iiki = imld_rst(ji,jj)
-           hbl (ji,jj) = gdepw_n(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth
+           hbl (ji,jj) = gdepw(ji,jj,iiki  ,Kmm) * ssmask(ji,jj)    ! Turbocline depth
         END DO
      END DO
@@ -1610,5 +1614,5 @@
 
 
-   SUBROUTINE tra_osm( kt )
+   SUBROUTINE tra_osm( kt, Kmm, pts, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_osm  ***
@@ -1620,5 +1624,8 @@
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) :: kt
+      INTEGER                                  , INTENT(in)    :: kt        ! time step index
+      INTEGER                                  , INTENT(in)    :: Kmm, Krhs ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts       ! active tracers and RHS of tracer equation
+      !
       INTEGER :: ji, jj, jk
       !
@@ -1630,6 +1637,6 @@
 
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
-         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
-         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsa(:,:,:,jp_sal)
+         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
+         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
       ENDIF
 
@@ -1638,10 +1645,10 @@
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               tsa(ji,jj,jk,jp_tem) =  tsa(ji,jj,jk,jp_tem)                      &
+               pts(ji,jj,jk,jp_tem,Krhs) =  pts(ji,jj,jk,jp_tem,Krhs)                      &
                   &                 - (  ghamt(ji,jj,jk  )  &
-                  &                    - ghamt(ji,jj,jk+1) ) /e3t_n(ji,jj,jk)
-               tsa(ji,jj,jk,jp_sal) =  tsa(ji,jj,jk,jp_sal)                      &
+                  &                    - ghamt(ji,jj,jk+1) ) /e3t(ji,jj,jk,Kmm)
+               pts(ji,jj,jk,jp_sal,Krhs) =  pts(ji,jj,jk,jp_sal,Krhs)                      &
                   &                 - (  ghams(ji,jj,jk  )  &
-                  &                    - ghams(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
+                  &                    - ghams(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
             END DO
          END DO
@@ -1651,6 +1658,6 @@
       ! save the non-local tracer flux trends for diagnostic
       IF( l_trdtra )   THEN
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
+         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
+         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
 !!bug gm jpttdzdf ==> jpttosm
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )
@@ -1660,6 +1667,6 @@
 
       IF(ln_ctl) THEN
-         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' osm  - Ta: ', mask1=tmask,   &
-         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' osm  - Ta: ', mask1=tmask,   &
+         &             tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       ENDIF
       !
@@ -1684,5 +1691,5 @@
 
 
-   SUBROUTINE dyn_osm( kt )
+   SUBROUTINE dyn_osm( kt, Kmm, puu, pvv, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_osm  ***
@@ -1693,5 +1700,7 @@
       !! ** Method  :   ???
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   !
+      INTEGER                             , INTENT( in )  ::  kt          ! ocean time step index
+      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation
       !
       INTEGER :: ji, jj, jk   ! dummy loop indices
@@ -1708,10 +1717,10 @@
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               ua(ji,jj,jk) =  ua(ji,jj,jk)                      &
+               puu(ji,jj,jk,Krhs) =  puu(ji,jj,jk,Krhs)                      &
                   &                 - (  ghamu(ji,jj,jk  )  &
-                  &                    - ghamu(ji,jj,jk+1) ) / e3u_n(ji,jj,jk)
-               va(ji,jj,jk) =  va(ji,jj,jk)                      &
+                  &                    - ghamu(ji,jj,jk+1) ) / e3u(ji,jj,jk,Kmm)
+               pvv(ji,jj,jk,Krhs) =  pvv(ji,jj,jk,Krhs)                      &
                   &                 - (  ghamv(ji,jj,jk  )  &
-                  &                    - ghamv(ji,jj,jk+1) ) / e3v_n(ji,jj,jk)
+                  &                    - ghamv(ji,jj,jk+1) ) / e3v(ji,jj,jk,Kmm)
             END DO
          END DO
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 10880)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfphy.F90	(revision 10883)
@@ -61,5 +61,5 @@
 CONTAINS
 
-   SUBROUTINE zdf_phy_init
+   SUBROUTINE zdf_phy_init( Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_phy_init  ***
@@ -70,4 +70,6 @@
       !!                set horizontal shape and vertical profile of background mixing coef.
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in)    :: Kmm ! time level index (middle)
+      !
       INTEGER ::   jk            ! dummy loop indices
       INTEGER ::   ioptio, ios   ! local integers
@@ -193,5 +195,5 @@
       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_zdfosm ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_OSM   ;   CALL zdf_osm_init   ;   ENDIF
+      IF( ln_zdfosm ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_OSM   ;   CALL zdf_osm_init( Kmm )   ;   ENDIF
       !
       IF( ioptio /= 1 )    CALL ctl_stop( 'zdf_phy_init: one and only one vertical diffusion option has to be defined ' )
@@ -218,5 +220,5 @@
 
 
-   SUBROUTINE zdf_phy( kt )
+   SUBROUTINE zdf_phy( kt, Kbb, Kmm )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE zdf_phy  ***
@@ -230,5 +232,6 @@
       !!                bottom stress.....                               <<<<====verifier !
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   kt         ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
       !
       INTEGER ::   ji, jj, jk   ! dummy loop indice
@@ -254,12 +257,12 @@
       !
       IF( l_zdfsh2 )   &         !* shear production at w-points (energy conserving form)
-         CALL zdf_sh2( ub, vb, un, vn, avm_k,   &     ! <<== in
-            &                           zsh2    )     ! ==>> out : shear production
+         CALL zdf_sh2( Kbb, Kmm, avm_k,   &     ! <<== in
+            &                     zsh2    )     ! ==>> out : shear production
       !
       SELECT CASE ( nzdf_phy )                  !* Vertical eddy viscosity and diffusivity coefficients at w-points
-      CASE( np_RIC )   ;   CALL zdf_ric( kt, gdept_n, zsh2, avm_k, avt_k )    ! Richardson number dependent Kz
-      CASE( np_TKE )   ;   CALL zdf_tke( kt         , zsh2, avm_k, avt_k )    ! TKE closure scheme for Kz
-      CASE( np_GLS )   ;   CALL zdf_gls( kt         , zsh2, avm_k, avt_k )    ! GLS closure scheme for Kz
-      CASE( np_OSM )   ;   CALL zdf_osm( kt               , avm_k, avt_k )    ! OSMOSIS closure scheme for Kz
+      CASE( np_RIC )   ;   CALL zdf_ric( kt,      Kmm, zsh2, avm_k, avt_k )    ! Richardson number dependent Kz
+      CASE( np_TKE )   ;   CALL zdf_tke( kt, Kbb, Kmm, zsh2, avm_k, avt_k )    ! TKE closure scheme for Kz
+      CASE( np_GLS )   ;   CALL zdf_gls( kt, Kbb, Kmm, zsh2, avm_k, avt_k )    ! GLS closure scheme for Kz
+      CASE( np_OSM )   ;   CALL zdf_osm( kt, Kbb, Kmm      , avm_k, avt_k )    ! OSMOSIS closure scheme for Kz
 !     CASE( np_CST )                                  ! Constant Kz (reset avt, avm to the background value)
 !         ! avt_k and avm_k set one for all at initialisation phase
@@ -318,5 +321,5 @@
          IF( ln_zdfgls )   CALL gls_rst( kt, 'WRITE' )
          IF( ln_zdfric )   CALL ric_rst( kt, 'WRITE' ) 
-         ! NB. OSMOSIS restart (osm_rst) will be called in step.F90 after wn has been updated
+         ! NB. OSMOSIS restart (osm_rst) will be called in step.F90 after ww has been updated
       ENDIF
       !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfric.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfric.F90	(revision 10880)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfric.F90	(revision 10883)
@@ -112,5 +112,5 @@
 
 
-   SUBROUTINE zdf_ric( kt, pdept, p_sh2, p_avm, p_avt )
+   SUBROUTINE zdf_ric( kt, Kmm, p_sh2, p_avm, p_avt )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE zdfric  ***
@@ -125,5 +125,5 @@
       !!                    avt = avm0 / (1 + rn_alp*ri)
       !!                with ri  = N^2 / dz(u)**2
-      !!                         = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ]
+      !!                         = e3w**2 * rn2/[ mi( dk(uu(:,:,:,Kbb)) )+mj( dk(vv(:,:,:,Kbb)) ) ]
       !!                    avm0= rn_avmri / (1 + rn_alp*Ri)**nn_ric
       !!                where ri is the before local Richardson number,
@@ -152,5 +152,5 @@
       !!----------------------------------------------------------------------
       INTEGER                   , INTENT(in   ) ::   kt             ! ocean time-step
-      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdept          ! depth of t-point  [m]
+      INTEGER                   , INTENT(in   ) ::   Kmm            ! ocean time level index
       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   ! momentum and tracer Kz (w-points)
@@ -189,5 +189,5 @@
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
-                  IF( pdept(ji,jj,jk) < zh_ekm(ji,jj) ) THEN
+                  IF( gdept(ji,jj,jk,Kmm) < zh_ekm(ji,jj) ) THEN
                      p_avm(ji,jj,jk) = MAX(  p_avm(ji,jj,jk), rn_wvmix  ) * wmask(ji,jj,jk)
                      p_avt(ji,jj,jk) = MAX(  p_avt(ji,jj,jk), rn_wtmix  ) * wmask(ji,jj,jk)
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfsh2.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfsh2.F90	(revision 10880)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfsh2.F90	(revision 10883)
@@ -11,4 +11,5 @@
    !!   zdf_sh2       : compute mixing the shear production term of TKE
    !!----------------------------------------------------------------------
+   USE oce
    USE dom_oce        ! domain: ocean
    !
@@ -28,5 +29,5 @@
 CONTAINS
 
-   SUBROUTINE zdf_sh2( pub, pvb, pun, pvn, p_avm, p_sh2  ) 
+   SUBROUTINE zdf_sh2( Kbb, Kmm, p_avm, p_sh2  ) 
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE zdf_sh2  ***
@@ -47,5 +48,5 @@
       !! References :   Bruchard, OM 2002
       !! ---------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   pub, pvb, pun, pvn   ! before, now horizontal velocities
+      INTEGER                    , INTENT(in   ) ::   Kbb, Kmm             ! ocean time level indices
       REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm                ! vertical eddy viscosity (w-points)
       REAL(wp), DIMENSION(:,:,:) , INTENT(  out) ::   p_sh2                ! shear production of TKE (w-points)
@@ -59,9 +60,9 @@
             DO ji = 1, jpim1
                zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) &
-                  &         * (   pun(ji,jj,jk-1) -   pun(ji,jj,jk) ) &
-                  &         * (   pub(ji,jj,jk-1) -   pub(ji,jj,jk) ) / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) * wumask(ji,jj,jk)
+                  &         * (   uu(ji,jj,jk-1,Kmm) -   uu(ji,jj,jk,Kmm) ) &
+                  &         * (   uu(ji,jj,jk-1,Kbb) -   uu(ji,jj,jk,Kbb) ) / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) * wumask(ji,jj,jk)
                zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) &
-                  &         * (   pvn(ji,jj,jk-1) -   pvn(ji,jj,jk) ) &
-                  &         * (   pvb(ji,jj,jk-1) -   pvb(ji,jj,jk) ) / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) * wvmask(ji,jj,jk)
+                  &         * (   vv(ji,jj,jk-1,Kmm) -   vv(ji,jj,jk,Kmm) ) &
+                  &         * (   vv(ji,jj,jk-1,Kbb) -   vv(ji,jj,jk,Kbb) ) / ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) * wvmask(ji,jj,jk)
             END DO
          END DO
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 10880)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdftke.F90	(revision 10883)
@@ -109,5 +109,5 @@
 
 
-   SUBROUTINE zdf_tke( kt, p_sh2, p_avm, p_avt )
+   SUBROUTINE zdf_tke( kt, Kbb, Kmm, p_sh2, p_avm, p_avt )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE zdf_tke  ***
@@ -155,16 +155,17 @@
       !!----------------------------------------------------------------------
       INTEGER                   , INTENT(in   ) ::   kt             ! ocean time step
+      INTEGER                   , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points)
       !!----------------------------------------------------------------------
       !
-      CALL tke_tke( gdepw_n, e3t_n, e3w_n, p_sh2, p_avm, p_avt )   ! now tke (en)
-      !
-      CALL tke_avn( gdepw_n, e3t_n, e3w_n,        p_avm, p_avt )   ! now avt, avm, dissl
+      CALL tke_tke( Kbb, Kmm, p_sh2, p_avm, p_avt )   ! now tke (en)
+      !
+      CALL tke_avn( Kbb, Kmm,        p_avm, p_avt )   ! now avt, avm, dissl
       !
   END SUBROUTINE zdf_tke
 
 
-   SUBROUTINE tke_tke( pdepw, p_e3t, p_e3w, p_sh2, p_avm, p_avt )
+   SUBROUTINE tke_tke( Kbb, Kmm, p_sh2, p_avm, p_avt )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tke_tke  ***
@@ -186,6 +187,5 @@
       USE zdf_oce , ONLY : en   ! ocean vertical physics
       !!
-      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   pdepw          ! depth of w-points
-      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_e3t, p_e3w   ! level thickness (t- & w-points)
+      INTEGER                    , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
       REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_sh2          ! shear production term
       REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points)
@@ -243,6 +243,6 @@
                zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
                !                       ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0)
-               zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2  &
-                  &                                           + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  )
+               zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2  &
+                  &                                           + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2  )
                en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj)
             END DO
@@ -254,6 +254,6 @@
                   zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )
                   !                             ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0)
-                  zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2  &
-                     &                                           + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  )
+                  zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )**2  &
+                     &                                           + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2  )
                   en(ji,jj,mikt(ji,jj)) = MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1))   ! masked at ocean surface
                END DO
@@ -268,7 +268,7 @@
          !
          !                        !* total energy produce by LC : cumulative sum over jk
-         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * pdepw(:,:,1) * p_e3w(:,:,1)
+         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm)
          DO jk = 2, jpk
-            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * pdepw(:,:,jk) * p_e3w(:,:,jk)
+            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm)
          END DO
          !                        !* finite Langmuir Circulation depth
@@ -286,5 +286,5 @@
          DO jj = 1, jpj 
             DO ji = 1, jpi
-               zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj))
+               zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm)
             END DO
          END DO
@@ -302,7 +302,7 @@
                   IF ( zfr_i(ji,jj) /= 0. ) THEN               
                      ! vertical velocity due to LC   
-                     IF ( pdepw(ji,jj,jk) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN
+                     IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN
                         !                                           ! vertical velocity due to LC
-                        zwlc = rn_lc * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) )   ! warning: optimization: zus^3 is in zfr_i
+                        zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) )   ! warning: optimization: zus^3 is in zfr_i
                         !                                           ! TKE Langmuir circulation source term
                         en(ji,jj,jk) = en(ji,jj,jk) + rdt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)
@@ -342,7 +342,7 @@
                !                                   ! eddy coefficient (ensure numerical stability)
                zzd_up = zcof * MAX(  p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal
-                  &          /    (  p_e3t(ji,jj,jk  ) * p_e3w(ji,jj,jk  )  )
+                  &          /    (  e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk  ,Kmm)  )
                zzd_lw = zcof * MAX(  p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal
-                  &          /    (  p_e3t(ji,jj,jk-1) * p_e3w(ji,jj,jk  )  )
+                  &          /    (  e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk  ,Kmm)  )
                !
                zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
@@ -402,5 +402,5 @@
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 !!gm BUG : in the exp  remove the depth of ssh !!!
-!!gm       i.e. use gde3w in argument (pdepw)
+!!gm       i.e. use gde3w in argument (gdepw(:,:,:,Kmm))
       
       
@@ -409,5 +409,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
+                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
                      &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
@@ -418,5 +418,5 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                jk = nmln(ji,jj)
-               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
+               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
                   &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
             END DO
@@ -431,5 +431,5 @@
                   zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean 
                   zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
-                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
+                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
                      &                        * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
@@ -441,5 +441,5 @@
 
 
-   SUBROUTINE tke_avn( pdepw, p_e3t, p_e3w, p_avm, p_avt )
+   SUBROUTINE tke_avn( Kbb, Kmm, p_avm, p_avt )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tke_avn  ***
@@ -477,6 +477,5 @@
       USE zdf_oce , ONLY : en, avtb, avmb, avtb_2d   ! ocean vertical physics
       !!
-      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdepw          ! depth (w-points)
-      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_e3t, p_e3w   ! level thickness (t- & w-points)
+      INTEGER                   , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points)
       !
@@ -526,14 +525,14 @@
       !
  !!gm Not sure of that coding for ISF....
-      ! where wmask = 0 set zmxlm == p_e3w
+      ! where wmask = 0 set zmxlm == e3w(:,:,:,Kmm)
       CASE ( 0 )           ! bounded by the distance to surface and bottom
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   &
-                  &            pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) )
+                  zemxl = MIN( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm), zmxlm(ji,jj,jk),   &
+                  &            gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - gdepw(ji,jj,jk,Kmm) )
                   ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
-                  zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk))
-                  zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk))
+                  zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk))
+                  zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk))
                END DO
             END DO
@@ -544,5 +543,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) )
+                  zemxl = MIN( e3w(ji,jj,jk,Kmm), zmxlm(ji,jj,jk) )
                   zmxlm(ji,jj,jk) = zemxl
                   zmxld(ji,jj,jk) = zemxl
@@ -555,5 +554,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
+                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
                END DO
             END DO
@@ -562,5 +561,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
+                  zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
                   zmxlm(ji,jj,jk) = zemxl
                   zmxld(ji,jj,jk) = zemxl
@@ -573,5 +572,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
+                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
                END DO
             END DO
@@ -580,5 +579,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
+                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
                END DO
             END DO
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 10880)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90	(revision 10883)
@@ -448,5 +448,5 @@
 
       !                                      ! Ocean physics
-                           CALL zdf_phy_init    ! Vertical physics
+                           CALL zdf_phy_init( Nnn )    ! Vertical physics
                                      
       !                                         ! Lateral physics
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 10880)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10883)
@@ -137,5 +137,5 @@
 
       !  VERTICAL PHYSICS
-                         CALL zdf_phy( kstp )         ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
+                         CALL zdf_phy( kstp, Nbb, Nnn )   ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
 
       !  LATERAL  PHYSICS
@@ -196,5 +196,5 @@
                          CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity           	==> RHS
                          CALL dyn_ldf       ( kstp )  ! lateral mixing
-      IF( ln_zdfosm  )   CALL dyn_osm       ( kstp )  ! OSMOSIS non-local velocity fluxes
+      IF( ln_zdfosm  )   CALL dyn_osm( kstp,      Nnn      , uu, vv, Nrhs )  ! OSMOSIS non-local velocity fluxes ==> RHS
                          CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure
                          CALL dyn_spg       ( kstp )  ! surface pressure gradient
@@ -253,7 +253,7 @@
 #endif
                          CALL tra_adv( kstp, Nbb, Nnn      , ts, Nrhs )  ! hor. + vert. advection	==> RHS
-      IF( ln_zdfosm  )   CALL tra_osm       ( kstp )  ! OSMOSIS non-local tracer fluxes
+      IF( ln_zdfosm  )   CALL tra_osm( kstp,      Nnn      , ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS
       IF( lrst_oce .AND. ln_zdfosm ) &
-           &             CALL osm_rst( kstp, 'WRITE' )! write OSMOSIS outputs + wn (so must do here) to restarts
+           &             CALL osm_rst( kstp, Nnn, 'WRITE' )! write OSMOSIS outputs + wn (so must do here) to restarts
                          CALL tra_ldf       ( kstp )  ! lateral mixing
 
