Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/eosbn2.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/eosbn2.F90	(revision 14664)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/eosbn2.F90	(revision 14665)
@@ -577,4 +577,14 @@
 
    SUBROUTINE eos_insitu_pot_2d( pts, prhop )
+      !!
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
+      !                                                     ! 2 : salinity               [psu]
+      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   prhop  ! potential density (surface referenced)
+      !!
+      CALL eos_insitu_pot_2d_t( pts, is_tile(pts), prhop, is_tile(prhop) )
+   END SUBROUTINE eos_insitu_pot_2d
+
+
+   SUBROUTINE eos_insitu_pot_2d_t( pts, ktts, prhop, ktrhop )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE eos_insitu_pot  ***
@@ -589,7 +599,8 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
+      INTEGER                              , INTENT(in   ) ::   ktts, ktrhop
+      REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
       !                                                                ! 2 : salinity               [psu]
-      REAL(wp), DIMENSION(jpi,jpj     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
+      REAL(wp), DIMENSION(A2D_T(ktrhop)   ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
       !
       INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices
@@ -606,26 +617,26 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-            DO_2D( 1, 1, 1, 1 )
-               !
-               zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
-               zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
-               ztm = tmask(ji,jj,1)                                         ! tmask
-               !
-               zn0 = (((((EOS060*zt   &
-                  &   + EOS150*zs+EOS050)*zt   &
-                  &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
-                  &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
-                  &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
-                  &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
-                  &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
-                  !
-               !
-               prhop(ji,jj) = zn0 * ztm                           ! potential density referenced at the surface
-               !
-            END_2D
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+            !
+            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
+            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+            ztm = tmask(ji,jj,1)                                         ! tmask
+            !
+            zn0 = (((((EOS060*zt   &
+               &   + EOS150*zs+EOS050)*zt   &
+               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
+               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
+               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
+               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
+               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+               !
+            !
+            prhop(ji,jj) = zn0 * ztm                           ! potential density referenced at the surface
+            !
+         END_2D
 
       CASE( np_seos )                !==  simplified EOS  ==!
          !
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             zt  = pts  (ji,jj,jp_tem) - 10._wp
             zs  = pts  (ji,jj,jp_sal) - 35._wp
@@ -646,5 +657,5 @@
       IF( ln_timing )   CALL timing_stop('eos-pot')
       !
-   END SUBROUTINE eos_insitu_pot_2d
+   END SUBROUTINE eos_insitu_pot_2d_t
 
 
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfddm.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfddm.F90	(revision 14664)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfddm.F90	(revision 14665)
@@ -83,5 +83,5 @@
       REAL(dp) ::          zavfs    !   -      -
       REAL(wp) ::   zavdt, zavds    !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
+      REAL(wp), DIMENSION(A2D(nn_hls)) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
       !!----------------------------------------------------------------------
       !
@@ -96,5 +96,5 @@
          
          ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )           !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==!
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )           !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==!
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )           !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==!
             zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
 !!gm please, use e3w at Kmm below 
@@ -113,5 +113,5 @@
 
          ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )           !==  indicators  ==!
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )           !==  indicators  ==!
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )           !==  indicators  ==!
             ! stability indicator: msks=1 if rn2>0; 0 elsewhere
             IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
@@ -137,5 +137,5 @@
          END_2D
          ! mask zmsk in order to have avt and avs masked
-         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk)
+         zmsks(:,:) = zmsks(:,:) * wmask(A2D(nn_hls),jk)
 
 
@@ -144,5 +144,5 @@
          ! Constant eddy coefficient: reset to the background value
          ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             zinr = 1._wp / zrau(ji,jj)
             ! salt fingering
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfdrg.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfdrg.F90	(revision 14664)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfdrg.F90	(revision 14665)
@@ -118,5 +118,5 @@
       IF( l_log_not_linssh ) THEN     !==  "log layer"  ==!   compute Cd and -Cd*|U|
          ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             imk = k_mk(ji,jj)          ! ocean bottom level at t-points
             zut = uu(ji,jj,imk,Kmm) + uu(ji-1,jj,imk,Kmm)     ! 2 x velocity at t-point
@@ -131,5 +131,5 @@
       ELSE                                            !==  standard Cd  ==!
          ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             imk = k_mk(ji,jj)    ! ocean bottom level at t-points
             zut = uu(ji,jj,imk,Kmm) + uu(ji-1,jj,imk,Kmm)     ! 2 x velocity at t-point
@@ -179,5 +179,5 @@
 
       ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          ikbu = mbku(ji,jj)          ! deepest wet ocean u- & v-levels
          ikbv = mbkv(ji,jj)
@@ -193,5 +193,5 @@
       IF( ln_isfcav ) THEN        ! ocean cavities
          ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels
             ikbv = mikv(ji,jj)
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfevd.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfevd.F90	(revision 14664)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfevd.F90	(revision 14665)
@@ -62,16 +62,24 @@
       !
       INTEGER ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zavt_evd, zavm_evd
+      ! NOTE: [tiling] use a SAVE array to store diagnostics, then send after all tiles are finished. This is necessary because p_avt/p_avm are modified on adjacent tiles when using nn_hls > 1. zavt_evd/zavm_evd are then zero on some points when subsequently calculated for these tiles.
+      REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   zavt_evd, zavm_evd
       !!----------------------------------------------------------------------
       !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'zdf_evd : Enhanced Vertical Diffusion (evd)'
-         IF(lwp) WRITE(numout,*) '~~~~~~~ '
-         IF(lwp) WRITE(numout,*)
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+         IF( kt == nit000 ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'zdf_evd : Enhanced Vertical Diffusion (evd)'
+            IF(lwp) WRITE(numout,*) '~~~~~~~ '
+            IF(lwp) WRITE(numout,*)
+         ENDIF
+
+         ALLOCATE( zavt_evd(jpi,jpj,jpk) )
+         IF( nn_evdm == 1 ) ALLOCATE( zavm_evd(jpi,jpj,jpk) )
       ENDIF
       !
       !
-      zavt_evd(:,:,:) = p_avt(:,:,:)         ! set avt prior to evd application
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
+         zavt_evd(ji,jj,jk) = p_avt(ji,jj,jk)         ! set avt prior to evd application
+      END_3D
       !
       SELECT CASE ( nn_evdm )
@@ -79,5 +87,7 @@
       CASE ( 1 )           !==  enhance tracer & momentum Kz  ==!   (if rn2<-1.e-12)
          !
-         zavm_evd(:,:,:) = p_avm(:,:,:)      ! set avm prior to evd application
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
+            zavm_evd(ji,jj,jk) = p_avm(ji,jj,jk)      ! set avm prior to evd application
+         END_3D
          !
 !! change last digits results
@@ -88,5 +98,5 @@
          !
          ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
             IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN
                p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)
@@ -95,6 +105,11 @@
          END_3D
          !
-         zavm_evd(:,:,:) = p_avm(:,:,:) - zavm_evd(:,:,:)   ! change in avm due to evd
-         CALL iom_put( "avm_evd", zavm_evd )                ! output this change
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
+            zavm_evd(ji,jj,jk) = p_avm(ji,jj,jk) - zavm_evd(ji,jj,jk)   ! change in avm due to evd
+         END_3D
+         IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
+            CALL iom_put( "avm_evd", zavm_evd )                ! output this change
+            DEALLOCATE( zavm_evd )
+         ENDIF
          !
       CASE DEFAULT         !==  enhance tracer Kz  ==!   (if rn2<-1.e-12) 
@@ -105,5 +120,5 @@
 
          ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
             IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 )   &
                p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)
@@ -112,6 +127,11 @@
       END SELECT 
       !
-      zavt_evd(:,:,:) = p_avt(:,:,:) - zavt_evd(:,:,:)   ! change in avt due to evd
-      CALL iom_put( "avt_evd", zavt_evd )              ! output this change
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
+         zavt_evd(ji,jj,jk) = p_avt(ji,jj,jk) - zavt_evd(ji,jj,jk)   ! change in avt due to evd
+      END_3D
+      IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
+         CALL iom_put( "avt_evd", zavt_evd )              ! output this change
+         DEALLOCATE( zavt_evd )
+      ENDIF
       IF( l_trdtra ) CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_evd, zavt_evd )
       !
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfgls.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfgls.F90	(revision 14664)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfgls.F90	(revision 14665)
@@ -137,8 +137,8 @@
       USE zdf_oce , ONLY : en, avtb, avmb   ! ocean vertical physics
       !!
-      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)
+      INTEGER                             , INTENT(in   ) ::   kt             ! ocean time step
+      INTEGER                             , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   p_sh2          ! shear production term
+      REAL(wp), DIMENSION(:,:,:)          , INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points)
       !
       INTEGER  ::   ji, jj, jk    ! dummy loop arguments
@@ -151,34 +151,36 @@
       REAL(wp) ::   gh, gm, shr, dif, zsqen, zavt, zavm !   -      -
       REAL(wp) ::   zmsku, zmskv                        !   -      -
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zdep
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zkar
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zflxs       ! Turbulence fluxed induced by internal waves
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zhsro       ! Surface roughness (surface waves)
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra    ! Tapering of wave breaking under sea ice
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eb          ! tke at time before
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   hmxl_b      ! mixing length at time before
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eps         ! dissipation rate
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwall_psi   ! Wall function use in the wb case (ln_sigpsi)
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   psi         ! psi at time now
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zd_lw, zd_up, zdiag   ! lower, upper  and diagonal of the matrix
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zstt, zstm  ! stability function on tracer and momentum
+      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zdep
+      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zkar
+      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zflxs                 ! Turbulence fluxed induced by internal waves
+      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zhsro                 ! Surface roughness (surface waves)
+      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zice_fra              ! Tapering of wave breaking under sea ice
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   eb                    ! tke at time before
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   hmxl_b                ! mixing length at time before
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   eps                   ! dissipation rate
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zwall_psi             ! Wall function use in the wb case (ln_sigpsi)
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   psi                   ! psi at time now
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zd_lw, zd_up, zdiag   ! lower, upper  and diagonal of the matrix
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zstt, zstm            ! stability function on tracer and momentum
       !!--------------------------------------------------------------------
       !
       ! Preliminary computing
-
-      ustar2_surf(:,:) = 0._wp   ;         psi(:,:,:) = 0._wp
-      ustar2_top (:,:) = 0._wp   ;   zwall_psi(:,:,:) = 0._wp
-      ustar2_bot (:,:) = 0._wp
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         ustar2_surf(ji,jj) = 0._wp   ;   ustar2_top(ji,jj) = 0._wp   ;   ustar2_bot(ji,jj) = 0._wp
+      END_2D
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
+         psi(ji,jj,jk) = 0._wp   ;   zwall_psi(ji,jj,jk) = 0._wp
+      END_3D
 
       SELECT CASE ( nn_z0_ice )
       CASE( 0 )   ;   zice_fra(:,:) = 0._wp
-      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(:,:) * 10._wp )
-      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(:,:)
-      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp )
+      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(A2D(nn_hls)) * 10._wp )
+      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(A2D(nn_hls))
+      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(A2D(nn_hls)) , 1._wp )
       END SELECT
 
       ! Compute surface, top and bottom friction at T-points
       ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )          !==  surface ocean friction
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )          !==  surface ocean friction 
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )          !==  surface ocean friction
          ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1)   ! surface friction
       END_2D
@@ -188,5 +190,5 @@
       IF( .NOT.ln_drg_OFF ) THEN     !== top/bottom friction   (explicit before friction)
          ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )         ! bottom friction (explicit before friction)
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )          ! bottom friction (explicit before friction)
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )          ! bottom friction (explicit before friction)
             zmsku = 0.5_wp * ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
             zmskv = 0.5_wp * ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0)
@@ -196,5 +198,5 @@
          IF( ln_isfcav ) THEN
             ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )      ! top friction
-            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )      ! top friction
+            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )      ! top friction
                zmsku = 0.5_wp * ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
                zmskv = 0.5_wp * ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )     ! (CAUTION: CdU<0)
@@ -209,17 +211,22 @@
          zhsro(:,:) = rn_hsro
       CASE ( 1 )             ! Standard Charnock formula
-         zhsro(:,:) = MAX( rsbc_zs1 * ustar2_surf(:,:) , rn_hsro )
+         zhsro(:,:) = MAX( rsbc_zs1 * ustar2_surf(A2D(nn_hls)) , rn_hsro )
       CASE ( 2 )             ! Roughness formulae according to Rascle et al., Ocean Modelling (2008)
 !!gm faster coding : the 2 comment lines should be used
 !!gm         zcof = 2._wp * 0.6_wp / 28._wp
 !!gm         zdep(:,:)  = 30._wp * TANH(  zcof/ SQRT( MAX(ustar2_surf(:,:),rsmall) )  )       ! Wave age (eq. 10)
-         zdep (:,:) = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(:,:),rsmall))) )         ! Wave age (eq. 10)
-         zhsro(:,:) = MAX(rsbc_zs2 * ustar2_surf(:,:) * zdep(:,:)**1.5, rn_hsro)          ! zhsro = rn_frac_hs * Hsw (eq. 11)
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            zcof = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(ji,jj),rsmall))) )          ! Wave age (eq. 10)
+            zhsro(ji,jj) = MAX(rsbc_zs2 * ustar2_surf(ji,jj) * zcof**1.5, rn_hsro)        ! zhsro = rn_frac_hs * Hsw (eq. 11)
+         END_2D
       CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file)
-         zhsro(:,:) = MAX(rn_frac_hs * hsw(:,:), rn_hsro)   ! (rn_frac_hs=1.6 see Eq. (5) of Rascle et al. 2008 )
+         zhsro(:,:) = MAX(rn_frac_hs * hsw(A2D(nn_hls)), rn_hsro)   ! (rn_frac_hs=1.6 see Eq. (5) of Rascle et al. 2008 )
       END SELECT
       !
       ! adapt roughness where there is sea ice
-      zhsro(:,:) = ( (1._wp-zice_fra(:,:)) * zhsro(:,:) + zice_fra(:,:) * rn_hsri )*tmask(:,:,1)  + (1._wp - tmask(:,:,1))*rn_hsro
+      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * rn_hsri )*tmask(ji,jj,1)  + &
+            &           (1._wp - tmask(ji,jj,1))*rn_hsro
+      END_2D
       !
       ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )  !==  Compute dissipation rate  ==!
@@ -229,10 +236,12 @@
 
       ! Save tke at before time step
-      eb    (:,:,:) = en    (:,:,:)
-      hmxl_b(:,:,:) = hmxl_n(:,:,:)
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
+         eb    (ji,jj,jk) = en    (ji,jj,jk)
+         hmxl_b(ji,jj,jk) = hmxl_n(ji,jj,jk)
+      END_3D
 
       IF( nn_clos == 0 ) THEN    ! Mellor-Yamada
          ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             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) )
@@ -256,5 +265,5 @@
 
       ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
          !
          buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk)     ! stratif. destruction
@@ -309,40 +318,43 @@
       !
       CASE ( 0 )             ! Dirichlet boundary condition (set e at k=1 & 2)
-      ! First level
-      en   (:,:,1) = MAX(  rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3  )
-      zd_lw(:,:,1) = en(:,:,1)
-      zd_up(:,:,1) = 0._wp
-      zdiag(:,:,1) = 1._wp
-      !
-      ! One level below
-      en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1 * ((zhsro(:,:)+gdepw(:,:,2,Kmm)) &
-         &                 / zhsro(:,:) )**(1.5_wp*ra_sf)  )**(2._wp/3._wp) , rn_emin   )
-      zd_lw(:,:,2) = 0._wp
-      zd_up(:,:,2) = 0._wp
-      zdiag(:,:,2) = 1._wp
-      !
-      !
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            ! First level
+            en   (ji,jj,1) = MAX(  rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1)**r2_3  )
+            zd_lw(ji,jj,1) = en(ji,jj,1)
+            zd_up(ji,jj,1) = 0._wp
+            zdiag(ji,jj,1) = 1._wp
+            !
+            ! One level below
+            en   (ji,jj,2) =  MAX( rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1          &
+               &                             * ((zhsro(ji,jj)+gdepw(ji,jj,2,Kmm)) / zhsro(ji,jj) )**(1.5_wp*ra_sf)  )**r2_3 )
+            zd_lw(ji,jj,2) = 0._wp
+            zd_up(ji,jj,2) = 0._wp
+            zdiag(ji,jj,2) = 1._wp
+         END_2D
+         !
+         !
       CASE ( 1 )             ! Neumann boundary condition (set d(e)/dz)
-      !
-      ! Dirichlet conditions at k=1
-      en   (:,:,1) = MAX(  rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3 , rn_emin  )
-      zd_lw(:,:,1) = en(:,:,1)
-      zd_up(:,:,1) = 0._wp
-      zdiag(:,:,1) = 1._wp
-      !
-      ! at k=2, set de/dz=Fw
-      !cbr
-      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )   ! zdiag zd_lw not defined/used on the halo
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )   ! zdiag zd_lw not defined/used on the halo
-         zdiag(ji,jj,2) = zdiag(ji,jj,2) +  zd_lw(ji,jj,2) ! Remove zd_lw from zdiag
-         zd_lw(ji,jj,2) = 0._wp
-      END_2D
-      zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:)) ))
-      zflxs(:,:)   = rsbc_tke2 * (1._wp-zice_fra(:,:)) * ustar2_surf(:,:)**1.5_wp * zkar(:,:) &
-          &                    * (  ( zhsro(:,:)+gdept(:,:,1,Kmm) ) / zhsro(:,:)  )**(1.5_wp*ra_sf)
+         !
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            ! Dirichlet conditions at k=1
+            en   (ji,jj,1) = MAX(  rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1)**r2_3  )
+            zd_lw(ji,jj,1) = en(ji,jj,1)
+            zd_up(ji,jj,1) = 0._wp
+            zdiag(ji,jj,1) = 1._wp
+            !
+            ! at k=2, set de/dz=Fw
+            !cbr
+            ! zdiag zd_lw not defined/used on the halo
+            zdiag(ji,jj,2) = zdiag(ji,jj,2) +  zd_lw(ji,jj,2) ! Remove zd_lw from zdiag
+            zd_lw(ji,jj,2) = 0._wp
+            !
+            zkar (ji,jj)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(ji,jj,1,Kmm)/zhsro(ji,jj)) ))
+            zflxs(ji,jj)   = rsbc_tke2 * (1._wp-zice_fra(ji,jj)) * ustar2_surf(ji,jj)**1.5_wp * zkar(ji,jj) &
+                &                    * (  ( zhsro(ji,jj)+gdept(ji,jj,1,Kmm) ) / zhsro(ji,jj)  )**(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)
-      !
-      !
+            en(ji,jj,2) = en(ji,jj,2) + zflxs(ji,jj) / e3w(ji,jj,2,Kmm)
+         END_2D
+         !
+         !
       END SELECT
 
@@ -356,5 +368,5 @@
          !                      ! Balance between the production and the dissipation terms
          ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
 !!gm This means that bottom and ocean w-level above have a specified "en" value.   Sure ????
 !!   With thick deep ocean level thickness, this may be quite large, no ???
@@ -373,7 +385,8 @@
          END_2D
          !
+         ! NOTE: [tiling] not tested- ctl_stop with ln_isfcav!
          IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
             ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                itop   = mikt(ji,jj)       ! k   top w-point
                itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
@@ -394,5 +407,5 @@
          !
          ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
             ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -408,7 +421,8 @@
             en   (ji,jj,ibot) = z_en
          END_2D
+         ! NOTE: [tiling] not tested- ctl_stop with ln_isfcav!
          IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
             ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                itop   = mikt(ji,jj)       ! k   top w-point
                itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
@@ -431,5 +445,5 @@
       ! ----------------------------------------------------------
       !
-      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
       DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
          zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
@@ -440,9 +454,11 @@
       END_3D
       ! [comm_cleanup] ! DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
-      DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+      DO_3DS_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
          en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
       END_3D
       !                                            ! set the minimum value of tke
-      en(:,:,:) = MAX( en(:,:,:), rn_emin )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
+         en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin )
+      END_3D
 
       !!----------------------------------------!!
@@ -535,44 +551,45 @@
       CASE ( 0 )             ! Dirichlet boundary conditions
          !
-         ! Surface value
-         zdep    (:,:)   = zhsro(:,:) * rl_sf ! Cosmetic
-         psi     (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-         zd_lw(:,:,1) = psi(:,:,1)
-         zd_up(:,:,1) = 0._wp
-         zdiag(:,:,1) = 1._wp
-         !
-         ! One level below
-         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
-         zd_up(:,:,2) = 0._wp
-         zdiag(:,:,2) = 1._wp
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            ! Surface value
+            zdep    (ji,jj)   = zhsro(ji,jj) * rl_sf ! Cosmetic
+            psi     (ji,jj,1) = rc0**rpp * en(ji,jj,1)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1)
+            zd_lw(ji,jj,1) = psi(ji,jj,1)
+            zd_up(ji,jj,1) = 0._wp
+            zdiag(ji,jj,1) = 1._wp
+            !
+            ! One level below
+            zkar    (ji,jj)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(ji,jj,2,Kmm)/zhsro(ji,jj) )))
+            zdep    (ji,jj)   = (zhsro(ji,jj) + gdepw(ji,jj,2,Kmm)) * zkar(ji,jj)
+            psi     (ji,jj,2) = rc0**rpp * en(ji,jj,2)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1)
+            zd_lw(ji,jj,2) = 0._wp
+            zd_up(ji,jj,2) = 0._wp
+            zdiag(ji,jj,2) = 1._wp
+         END_2D
          !
       CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
          !
-         ! Surface value: Dirichlet
-         zdep    (:,:)   = zhsro(:,:) * rl_sf
-         psi     (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-         zd_lw(:,:,1) = psi(:,:,1)
-         zd_up(:,:,1) = 0._wp
-         zdiag(:,:,1) = 1._wp
-         !
-         ! Neumann condition at k=2
-         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )   ! zdiag zd_lw not defined/used on the halo
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )   ! zdiag zd_lw not defined/used on the halo
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            ! Surface value: Dirichlet
+            zdep    (ji,jj)   = zhsro(ji,jj) * rl_sf
+            psi     (ji,jj,1) = rc0**rpp * en(ji,jj,1)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1)
+            zd_lw(ji,jj,1) = psi(ji,jj,1)
+            zd_up(ji,jj,1) = 0._wp
+            zdiag(ji,jj,1) = 1._wp
+            !
+            ! Neumann condition at k=2, zdiag zd_lw not defined/used on the halo
             zdiag(ji,jj,2) = zdiag(ji,jj,2) +  zd_lw(ji,jj,2) ! Remove zd_lw from zdiag
             zd_lw(ji,jj,2) = 0._wp
+            !
+            ! Set psi vertical flux at the surface:
+            zkar (ji,jj)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(ji,jj,1,Kmm)/zhsro(ji,jj) )) ! Lengh scale slope
+            zdep (ji,jj)   = ((zhsro(ji,jj) + gdept(ji,jj,1,Kmm)) / zhsro(ji,jj))**(rmm*ra_sf)
+            zflxs(ji,jj)   = (rnn + (1._wp-zice_fra(ji,jj))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(ji,jj)) &
+               &           *(1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1*zdep(ji,jj))**(2._wp*rmm/3._wp-1_wp)
+            zdep (ji,jj)   = rsbc_psi1 * (zwall_psi(ji,jj,1)*p_avm(ji,jj,1)+zwall_psi(ji,jj,2)*p_avm(ji,jj,2)) * &
+               &           ustar2_surf(ji,jj)**rmm * zkar(ji,jj)**rnn * (zhsro(ji,jj) + gdept(ji,jj,1,Kmm))**(rnn-1.)
+            zflxs(ji,jj)   = zdep(ji,jj) * zflxs(ji,jj)
+            psi  (ji,jj,2) = psi(ji,jj,2) + zflxs(ji,jj) / e3w(ji,jj,2,Kmm)
          END_2D
-         !
-         ! Set psi vertical flux at the surface:
-         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 + (1._wp-zice_fra(:,:))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:)) &
-            &           *(1._wp + (1._wp-zice_fra(:,:))*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(:,:,1,Kmm))**(rnn-1.)
-         zflxs(:,:)   = zdep(:,:) * zflxs(:,:)
-         psi  (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm)
          !
       END SELECT
@@ -638,5 +655,5 @@
       ! ----------------
       !
-      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
       DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
          zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
@@ -688,12 +705,15 @@
       ! --------------------------------------------------
       ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! Note that this set boundary conditions on hmxl_n at the same time
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )   ! Note that this set boundary conditions on hmxl_n at the same time
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )   ! Note that this set boundary conditions on hmxl_n at the same time
          ! limitation
          eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
          hmxl_n(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
-         ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)
-         zrn2 = MAX( rn2(ji,jj,jk), rsmall )
-         IF( ln_length_lim )   hmxl_n(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) )
-      END_3D
+      END_3D
+      IF( ln_length_lim ) THEN        ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
+            zrn2 = MAX( rn2(ji,jj,jk), rsmall )
+            hmxl_n(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) )
+         END_3D
+      ENDIF
 
       !
@@ -760,6 +780,6 @@
       END_2D
 
-      zstt(:,:,  1) = wmask(:,:,  1)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0)
-      zstt(:,:,jpk) = wmask(:,:,jpk)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0)
+      zstt(:,:,  1) = wmask(A2D(nn_hls),  1)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0)
+      zstt(:,:,jpk) = wmask(A2D(nn_hls),jpk)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0)
 
 !!gm should be done for ISF (top boundary cond.)
@@ -772,5 +792,5 @@
       !     for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0)
       ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpk )
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
          zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk)
          zavt  = zsqen * zstt(ji,jj,jk)
@@ -779,5 +799,5 @@
          p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom
       END_3D
-      p_avt(:,:,1) = 0._wp
+      p_avt(A2D(nn_hls),1) = 0._wp
       !
       IF(sn_cfctl%l_prtctl) THEN
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfiwm.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfiwm.F90	(revision 14664)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfiwm.F90	(revision 14665)
@@ -125,15 +125,17 @@
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zztmp, ztmp1, ztmp2        ! scalar workspace
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zfact       ! Used for vertical structure
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zhdep       ! Ocean depth
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwkb        ! WKB-stretched height above bottom
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zweight     ! Weight for high mode vertical distribution
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   znu_t       ! Molecular kinematic viscosity (T grid)
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   znu_w       ! Molecular kinematic viscosity (W grid)
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zReb        ! Turbulence intensity parameter
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zemx_iwm    ! local energy density available for mixing (W/kg)
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zav_ratio   ! S/T diffusivity ratio (only for ln_tsdiff=T)
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zav_wave    ! Internal wave-induced diffusivity
+      ! NOTE: [tiling] declare zztmp with SAVE for mpp_sum
+      REAL(wp), SAVE :: zztmp
+      REAL(wp)       :: ztmp1, ztmp2        ! scalar workspace
+      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zfact       ! Used for vertical structure
+      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zhdep       ! Ocean depth
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zwkb        ! WKB-stretched height above bottom
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zweight     ! Weight for high mode vertical distribution
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   znu_t       ! Molecular kinematic viscosity (T grid)
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   znu_w       ! Molecular kinematic viscosity (W grid)
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zReb        ! Turbulence intensity parameter
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zemx_iwm    ! local energy density available for mixing (W/kg)
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zav_ratio   ! S/T diffusivity ratio (only for ln_tsdiff=T)
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zav_wave    ! Internal wave-induced diffusivity
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z3d  ! 3D workspace used for iom_put 
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   z2d  ! 2D     -      -    -     -
@@ -337,5 +339,5 @@
       !
       IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave
-         zztmp = 0._wp
+         IF( .NOT. l_istiled .OR. ntile == 1 ) zztmp = 0._wp                    ! Do only on the first tile
 !!gm used of glosum 3D....
          ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
@@ -344,13 +346,16 @@
                &          * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
          END_3D
-         CALL mpp_sum( 'zdfiwm', zztmp )
-         zztmp = rho0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing 
-         !
-         IF(lwp) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) 'zdf_iwm : Internal wave-driven mixing (iwm)'
-            WRITE(numout,*) '~~~~~~~ '
-            WRITE(numout,*)
-            WRITE(numout,*) '      Total power consumption by av_wave =  ', zztmp * 1.e-12_wp, 'TW'
+
+         IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
+            CALL mpp_sum( 'zdfiwm', zztmp )
+            zztmp = rho0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing
+            !
+            IF(lwp) THEN
+               WRITE(numout,*)
+               WRITE(numout,*) 'zdf_iwm : Internal wave-driven mixing (iwm)'
+               WRITE(numout,*) '~~~~~~~ '
+               WRITE(numout,*)
+               WRITE(numout,*) '      Total power consumption by av_wave =  ', zztmp * 1.e-12_wp, 'TW'
+            ENDIF
          ENDIF
       ENDIF
@@ -373,5 +378,5 @@
          CALL iom_put( "av_ratio", zav_ratio )
          ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )    !* update momentum & tracer diffusivity with wave-driven mixing
-         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )    !* update momentum & tracer diffusivity with wave-driven mixing
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )    !* update momentum & tracer diffusivity with wave-driven mixing
             p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) * zav_ratio(ji,jj,jk)
             p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk)
@@ -381,5 +386,5 @@
       ELSE                                !* update momentum & tracer diffusivity with wave-driven mixing
          ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk)
             p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk)
@@ -394,13 +399,9 @@
                                           !  vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm)
       IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN
-         ALLOCATE( z2d(jpi,jpj) , z3d(jpi,jpj,jpk) )
+         ALLOCATE( z2d(A2D(nn_hls)) , z3d(A2D(nn_hls),jpk) )
          ! Initialisation for iom_put
          DO_2D( 0, 0, 0, 0 )
             z3d(ji,jj,1) = 0._wp   ;   z3d(ji,jj,jpk) = 0._wp
          END_2D
-         z3d(           1:nn_hls,:,:) = 0._wp   ;   z3d(:,           1:nn_hls,:) = 0._wp
-         z3d(jpi-nn_hls+1:jpi   ,:,:) = 0._wp   ;   z3d(:,jpj-nn_hls+1:   jpj,:) = 0._wp
-         z2d(           1:nn_hls,:  ) = 0._wp   ;   z2d(:,           1:nn_hls  ) = 0._wp
-         z2d(jpi-nn_hls+1:jpi   ,:  ) = 0._wp   ;   z2d(:,jpj-nn_hls+1:   jpj  ) = 0._wp
 
          DO_3D( 0, 0, 0, 0, 2, jpkm1 )
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfmfc.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfmfc.F90	(revision 14664)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfmfc.F90	(revision 14665)
@@ -96,24 +96,24 @@
       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
-      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   ztsp         ! T/S of the plume
-      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   ztse         ! T/S at W point
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrwp          !
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrwp2         !
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zapp          !
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zedmf         !
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zepsT, zepsW  !
-      !
-      REAL(wp), DIMENSION(jpi,jpj) :: zustar, zustar2   !
-      REAL(wp), DIMENSION(jpi,jpj) :: zuws, zvws, zsws, zfnet          !
-      REAL(wp), DIMENSION(jpi,jpj) :: zfbuo, zrautbm1, zrautb, zraupl
-      REAL(wp), DIMENSION(jpi,jpj) :: zwpsurf            !
-      REAL(wp), DIMENSION(jpi,jpj) :: zop0 , zsp0 !
-      REAL(wp), DIMENSION(jpi,jpj) :: zrwp_0, zrwp2_0  !
-      REAL(wp), DIMENSION(jpi,jpj) :: zapp0           !
-      REAL(wp), DIMENSION(jpi,jpj) :: zphp, zph, zphpm1, zphm1, zNHydro
-      REAL(wp), DIMENSION(jpi,jpj) :: zhcmo          !
-      !
-      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zn2    ! N^2
-      REAL(wp), DIMENSION(jpi,jpj,2  ) ::   zab, zabm1, zabp ! alpha and beta
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk,2) ::   ztsp         ! T/S of the plume
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk,2) ::   ztse         ! T/S at W point
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zrwp          !
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zrwp2         !
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zapp          !
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zedmf         !
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zepsT, zepsW  !
+      !
+      REAL(wp), DIMENSION(A2D(nn_hls)) :: zustar, zustar2   !
+      REAL(wp), DIMENSION(A2D(nn_hls)) :: zuws, zvws, zsws, zfnet          !
+      REAL(wp), DIMENSION(A2D(nn_hls)) :: zfbuo, zrautbm1, zrautb, zraupl
+      REAL(wp), DIMENSION(A2D(nn_hls)) :: zwpsurf            !
+      REAL(wp), DIMENSION(A2D(nn_hls)) :: zop0 , zsp0 !
+      REAL(wp), DIMENSION(A2D(nn_hls)) :: zrwp_0, zrwp2_0  !
+      REAL(wp), DIMENSION(A2D(nn_hls)) :: zapp0           !
+      REAL(wp), DIMENSION(A2D(nn_hls)) :: zphp, zph, zphpm1, zphm1, zNHydro
+      REAL(wp), DIMENSION(A2D(nn_hls)) :: zhcmo          !
+      !
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk)   ::   zn2    ! N^2
+      REAL(wp), DIMENSION(A2D(nn_hls),2  ) ::   zab, zabm1, zabp ! alpha and beta
      
       REAL(wp), PARAMETER :: zepsilon = 1.e-30                 ! local small value
@@ -136,35 +136,37 @@
       zcd          = 1._wp
 
-      !------------------------------------------------------------------
-      ! Surface boundary condition
-      !------------------------------------------------------------------
-      ! surface Stress
-      !--------------------
-      zuws(:,:) = utau(:,:) * r1_rho0 
-      zvws(:,:) = vtau(:,:) * r1_rho0 
-      zustar2(:,:) = SQRT(zuws(:,:)*zuws(:,:)+zvws(:,:)*zvws(:,:))
-      zustar(:,:)  = SQRT(zustar2(:,:))
-
-      ! Heat Flux
-      !--------------------
-      zfnet(:,:) = qns(:,:) + qsr(:,:)
-      zfnet(:,:) = zfnet(:,:) / (rho0 * rcp)
-
-      ! Water Flux
-      !---------------------
-      zsws(:,:) = emp(:,:)
-
-      !-------------------------------------------
-      ! Initialisation of prognostic variables
-      !-------------------------------------------
-      zrwp (:,:,:) =  0._wp ; zrwp2(:,:,:) =  0._wp ; zedmf(:,:,:) =  0._wp
-      zph  (:,:)   =  0._wp ; zphm1(:,:)   =  0._wp ; zphpm1(:,:)  =  0._wp
-      ztsp(:,:,:,:)=  0._wp
-
-      ! Tracers inside plume (ztsp) and environment (ztse)
-      ztsp(:,:,1,jp_tem) = pts(:,:,1,jp_tem,Kmm) * tmask(:,:,1)
-      ztsp(:,:,1,jp_sal) = pts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
-      ztse(:,:,1,jp_tem) = pts(:,:,1,jp_tem,Kmm) * tmask(:,:,1)
-      ztse(:,:,1,jp_sal) = pts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         !------------------------------------------------------------------
+         ! Surface boundary condition
+         !------------------------------------------------------------------
+         ! surface Stress
+         !--------------------
+         zuws(ji,jj) = utau(ji,jj) * r1_rho0
+         zvws(ji,jj) = vtau(ji,jj) * r1_rho0
+         zustar2(ji,jj) = SQRT(zuws(ji,jj)*zuws(ji,jj)+zvws(ji,jj)*zvws(ji,jj))
+         zustar(ji,jj)  = SQRT(zustar2(ji,jj))
+
+         ! Heat Flux
+         !--------------------
+         zfnet(ji,jj) = qns(ji,jj) + qsr(ji,jj)
+         zfnet(ji,jj) = zfnet(ji,jj) / (rho0 * rcp)
+
+         ! Water Flux
+         !---------------------
+         zsws(ji,jj) = emp(ji,jj)
+
+         !-------------------------------------------
+         ! Initialisation of prognostic variables
+         !-------------------------------------------
+         zrwp (ji,jj,:) =  0._wp ; zrwp2(ji,jj,:) =  0._wp ; zedmf(ji,jj,:) =  0._wp
+         zph  (ji,jj)   =  0._wp ; zphm1(ji,jj)   =  0._wp ; zphpm1(ji,jj)  =  0._wp
+         ztsp(ji,jj,:,:)=  0._wp
+
+         ! Tracers inside plume (ztsp) and environment (ztse)
+         ztsp(ji,jj,1,jp_tem) = pts(ji,jj,1,jp_tem,Kmm) * tmask(ji,jj,1)
+         ztsp(ji,jj,1,jp_sal) = pts(ji,jj,1,jp_sal,Kmm) * tmask(ji,jj,1)
+         ztse(ji,jj,1,jp_tem) = pts(ji,jj,1,jp_tem,Kmm) * tmask(ji,jj,1)
+         ztse(ji,jj,1,jp_sal) = pts(ji,jj,1,jp_sal,Kmm) * tmask(ji,jj,1)
+      END_2D
 
       CALL eos( ztse(:,:,1,:) ,  zrautb(:,:) )
@@ -174,8 +176,10 @@
       ! Boundary Condition of Mass Flux (plume velo.; convective area, entrain/detrain)
       !-------------------------------------------
-      zhcmo(:,:) = e3t(:,:,1,Kmm)
+      zhcmo(:,:) = e3t(A1Di(nn_hls),A1Dj(nn_hls),1,Kmm)
       zfbuo(:,:)   = 0._wp
       WHERE ( ABS(zrautb(:,:)) > 1.e-20 ) zfbuo(:,:)   =   &
-         &      grav * ( 2.e-4_wp *zfnet(:,:) - 7.6E-4_wp*pts(:,:,1,jp_sal,Kmm)*zsws(:,:)/zrautb(:,:)) * zhcmo(:,:)
+         &      grav * ( 2.e-4_wp *zfnet(:,:)              &
+         &      - 7.6E-4_wp*pts(A2D(nn_hls),1,jp_sal,Kmm)  &
+         &      * zsws(:,:)/zrautb(:,:)) * zhcmo(:,:)
 
       zedmf(:,:,1) = -0.065_wp*(ABS(zfbuo(:,:)))**(1._wp/3._wp)*SIGN(1.,zfbuo(:,:))
@@ -211,13 +215,14 @@
          CALL eos( ztsp(:,:,jk-1,:    ) ,  zraupl(:,:)   )
 
-         zphm1(:,:)  = zphm1(:,:)  + grav * zrautbm1(:,:) * e3t(:,:,jk-1, Kmm)
-         zphpm1(:,:) = zphpm1(:,:) + grav * zraupl(:,:)   * e3t(:,:,jk-1, Kmm)
-         zph(:,:)    = zphm1(:,:)  + grav * zrautb(:,:)   * e3t(:,:,jk  , Kmm)
-         zph(:,:)    = MAX( zph(:,:), zepsilon)
+         DO_2D( 0, 0, 0, 0 )
+            zphm1(ji,jj)  = zphm1(ji,jj)  + grav * zrautbm1(ji,jj) * e3t(ji,jj,jk-1, Kmm)
+            zphpm1(ji,jj) = zphpm1(ji,jj) + grav * zraupl(ji,jj)   * e3t(ji,jj,jk-1, Kmm)
+            zph(ji,jj)    = zphm1(ji,jj)  + grav * zrautb(ji,jj)   * e3t(ji,jj,jk  , Kmm)
+            zph(ji,jj)    = MAX( zph(ji,jj), zepsilon)
+         END_2D
 
          WHERE(zrautbm1 .NE. 0.) zfbuo(:,:)  =  grav * (zraupl(:,:) - zrautbm1(:,:)) / zrautbm1(:,:)
 
-         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D( 0, 0, 0, 0 )
 
             ! Compute Environment of Plume. Interpolation T/S (before time step) on W-points
@@ -323,8 +328,10 @@
 
       ! Compute Mass Flux on T-point
-      DO jk=1,jpk-1
-         edmfm(:,:,jk) = (zedmf(:,:,jk+1)  + zedmf(:,:,jk) )*0.5_wp
-      END DO
-      edmfm(:,:,jpk) = zedmf(:,:,jpk) 
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         edmfm(ji,jj,jk) = (zedmf(ji,jj,jk+1)  + zedmf(ji,jj,jk) )*0.5_wp
+      END_3D
+      DO_2D( 0, 0, 0, 0 )
+         edmfm(ji,jj,jpk) = zedmf(ji,jj,jpk)
+      END_2D
 
       ! Save variable (on T point)
@@ -339,44 +346,48 @@
       !  Computation of a tridiagonal matrix and right hand side terms of the linear system
       !=================================================================================
-      edmfa(:,:,:)     = 0._wp
-      edmfb(:,:,:)     = 0._wp
-      edmfc(:,:,:)     = 0._wp
-      edmftra(:,:,:,:) = 0._wp
+      DO_3D( 0, 0, 0, 0, 1, jpk )
+         edmfa(ji,jj,jk)     = 0._wp
+         edmfb(ji,jj,jk)     = 0._wp
+         edmfc(ji,jj,jk)     = 0._wp
+         edmftra(ji,jj,jk,:) = 0._wp
+      END_3D
 
       !---------------------------------------------------------------
       ! Diagonal terms 
       !---------------------------------------------------------------
-      DO jk=1,jpk-1
-         edmfa(:,:,jk) =  0._wp
-         edmfb(:,:,jk) = -edmfm(:,:,jk  ) / e3w(:,:,jk+1,Kmm)
-         edmfc(:,:,jk) =  edmfm(:,:,jk+1) / e3w(:,:,jk+1,Kmm)
-      END DO
-      edmfa(:,:,jpk)   = -edmfm(:,:,jpk-1) / e3w(:,:,jpk,Kmm)
-      edmfb(:,:,jpk)   =  edmfm(:,:,jpk  ) / e3w(:,:,jpk,Kmm)
-      edmfc(:,:,jpk)   =  0._wp
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         edmfa(ji,jj,jk) =  0._wp
+         edmfb(ji,jj,jk) = -edmfm(ji,jj,jk  ) / e3w(ji,jj,jk+1,Kmm)
+         edmfc(ji,jj,jk) =  edmfm(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
+      END_3D
+      DO_2D( 0, 0, 0, 0 )
+         edmfa(ji,jj,jpk)   = -edmfm(ji,jj,jpk-1) / e3w(ji,jj,jpk,Kmm)
+         edmfb(ji,jj,jpk)   =  edmfm(ji,jj,jpk  ) / e3w(ji,jj,jpk,Kmm)
+         edmfc(ji,jj,jpk)   =  0._wp
+      END_2D
 
       !---------------------------------------------------------------
       ! right hand side term for Temperature
       !---------------------------------------------------------------
-      DO jk=1,jpk-1
-        edmftra(:,:,jk,1) = - edmfm(:,:,jk  ) * ztsp(:,:,jk  ,jp_tem) / e3w(:,:,jk+1,Kmm) &
-                          & + edmfm(:,:,jk+1) * ztsp(:,:,jk+1,jp_tem) / e3w(:,:,jk+1,Kmm)
-      END DO
-      edmftra(:,:,jpk,1) = - edmfm(:,:,jpk-1) * ztsp(:,:,jpk-1,jp_tem) / e3w(:,:,jpk,Kmm) &
-                         & + edmfm(:,:,jpk  ) * ztsp(:,:,jpk  ,jp_tem) / e3w(:,:,jpk,Kmm)
-                        
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+        edmftra(ji,jj,jk,1) = - edmfm(ji,jj,jk  ) * ztsp(ji,jj,jk  ,jp_tem) / e3w(ji,jj,jk+1,Kmm) &
+                            & + edmfm(ji,jj,jk+1) * ztsp(ji,jj,jk+1,jp_tem) / e3w(ji,jj,jk+1,Kmm)
+      END_3D
+      DO_2D( 0, 0, 0, 0 )
+         edmftra(ji,jj,jpk,1) = - edmfm(ji,jj,jpk-1) * ztsp(ji,jj,jpk-1,jp_tem) / e3w(ji,jj,jpk,Kmm) &
+                              & + edmfm(ji,jj,jpk  ) * ztsp(ji,jj,jpk  ,jp_tem) / e3w(ji,jj,jpk,Kmm)
+      END_2D
+
       !---------------------------------------------------------------
       ! Right hand side term for Salinity
       !---------------------------------------------------------------
-      DO jk=1,jpk-1
-         edmftra(:,:,jk,2) =  - edmfm(:,:,jk  ) * ztsp(:,:,jk  ,jp_sal) / e3w(:,:,jk+1,Kmm) &
-                           &  + edmfm(:,:,jk+1) * ztsp(:,:,jk+1,jp_sal) / e3w(:,:,jk+1,Kmm)
-      END DO
-      edmftra(:,:,jpk,2) = - edmfm(:,:,jpk-1) * ztsp(:,:,jpk-1,jp_sal) / e3w(:,:,jpk,Kmm) &
-                         & + edmfm(:,:,jpk  ) * ztsp(:,:,jpk  ,jp_sal) / e3w(:,:,jpk,Kmm)
-      !
-      !
-      ! [comm_cleanup] 
-      IF (nn_hls.eq.1) CALL lbc_lnk( 'zdfmfc', edmfm,'T',1., edmfa,'T',1., edmfb,'T',1., edmfc,'T',1., edmftra(:,:,:,jp_tem),'T',1., edmftra(:,:,:,jp_sal),'T',1.)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         edmftra(ji,jj,jk,2) =  - edmfm(ji,jj,jk  ) * ztsp(ji,jj,jk  ,jp_sal) / e3w(ji,jj,jk+1,Kmm) &
+                             &  + edmfm(ji,jj,jk+1) * ztsp(ji,jj,jk+1,jp_sal) / e3w(ji,jj,jk+1,Kmm)
+      END_3D
+      DO_2D( 0, 0, 0, 0 )
+         edmftra(ji,jj,jpk,2) = - edmfm(ji,jj,jpk-1) * ztsp(ji,jj,jpk-1,jp_sal) / e3w(ji,jj,jpk,Kmm) &
+                              & + edmfm(ji,jj,jpk  ) * ztsp(ji,jj,jpk  ,jp_sal) / e3w(ji,jj,jpk,Kmm)
+      END_2D
       !
    END SUBROUTINE tra_mfc
@@ -385,7 +396,7 @@
    SUBROUTINE diag_mfc( zdiagi, zdiagd, zdiags, p2dt, Kaa )
 
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  zdiagi, zdiagd, zdiags  ! inout: tridaig. terms 
-      REAL(wp)                        , INTENT(in   ) ::   p2dt                   ! tracer time-step
-      INTEGER                         , INTENT(in   ) ::   Kaa                    ! ocean time level indices
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) ::  zdiagi, zdiagd, zdiags  ! inout: tridaig. terms
+      REAL(wp)                            , INTENT(in   ) ::   p2dt                   ! tracer time-step
+      INTEGER                             , INTENT(in   ) ::   Kaa                    ! ocean time level indices
 
       INTEGER  ::   ji, jj, jk  ! dummy  loop arguments   
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfmxl.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfmxl.F90	(revision 14664)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfmxl.F90	(revision 14665)
@@ -26,5 +26,5 @@
    PRIVATE
 
-   PUBLIC   zdf_mxl   ! called by zdfphy.F90
+   PUBLIC   zdf_mxl, zdf_mxl_turb   ! called by zdfphy.F90
 
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmln    !: number of level in the mixed layer (used by LDF, ZDF, TRD, TOP)
@@ -65,16 +65,11 @@
       !!                  ***  ROUTINE zdfmxl  ***
       !!                   
-      !! ** Purpose :   Compute the turbocline depth and the mixed layer depth
-      !!              with density criteria.
+      !! ** Purpose :   Compute the mixed layer depth with density criteria.
       !!
       !! ** Method  :   The mixed layer depth is the shallowest W depth with 
       !!      the density of the corresponding T point (just bellow) bellow a
       !!      given value defined locally as rho(10m) + rho_c
-      !!               The turbocline depth is the depth at which the vertical
-      !!      eddy diffusivity coefficient (resulting from the vertical physics
-      !!      alone, not the isopycnal part, see trazdf.F) fall below a given
-      !!      value defined locally (avt_c here taken equal to 5 cm/s2 by default)
       !!
-      !! ** Action  :   nmln, hmld, hmlp, hmlpt
+      !! ** Action  :   nmln, hmlp, hmlpt
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
@@ -82,23 +77,26 @@
       !
       INTEGER  ::   ji, jj, jk      ! dummy loop indices
-      INTEGER  ::   iikn, iiki, ikt ! local integer
+      INTEGER  ::   iik, ikt        ! local integer
       REAL(wp) ::   zN2_c           ! local scalar
-      INTEGER, DIMENSION(jpi,jpj) ::   imld   ! 2D workspace
       !!----------------------------------------------------------------------
       !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth'
-         IF(lwp) WRITE(numout,*) '~~~~~~~ '
-         !                             ! allocate zdfmxl arrays
-         IF( zdf_mxl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl : unable to allocate arrays' )
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+         IF( kt == nit000 ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth'
+            IF(lwp) WRITE(numout,*) '~~~~~~~ '
+            !                             ! allocate zdfmxl arrays
+            IF( zdf_mxl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl : unable to allocate arrays' )
+         ENDIF
       ENDIF
       !
       ! w-level of the mixing and mixed layers
-      nmln(:,:)  = nlb10                  ! Initialization to the number of w ocean point
-      hmlp(:,:)  = 0._wp                  ! here hmlp used as a dummy variable, integrating vertically N^2
+      DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
+         nmln(ji,jj)  = nlb10                  ! Initialization to the number of w ocean point
+         hmlp(ji,jj)  = 0._wp                  ! here hmlp used as a dummy variable, integrating vertically N^2
+      END_2D
       zN2_c = grav * rho_c * r1_rho0      ! convert density criteria into N^2 criteria
       ! [comm_cleanup] ! DO_3D( 1, 1, 1, 1, nlb10, jpkm1 )   ! Mixed layer level: w-level
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, nlb10, jpkm1 )   ! Mixed layer level: w-level
+      DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, nlb10, jpkm1 )   ! Mixed layer level: w-level
          ikt = mbkt(ji,jj)
          hmlp(ji,jj) =   &
@@ -106,32 +104,16 @@
          IF( hmlp(ji,jj) < zN2_c )   nmln(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
       END_3D
-      !
-      ! w-level of the turbocline and mixing layer (iom_use)
-      imld(:,:) = mbkt(:,:) + 1                ! Initialization to the number of w ocean point
-      ! [comm_cleanup] ! DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )   ! from the bottom to nlb10
-      DO_3DS( nn_hls, nn_hls, nn_hls, nn_hls, jpkm1, nlb10, -1 )   ! from the bottom to nlb10
-         IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) )   imld(ji,jj) = jk      ! Turbocline 
-      END_3D
-      ! depth of the mixing and mixed layers
+      ! depth of the mixed layer
       ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-         iiki = imld(ji,jj)
-         iikn = nmln(ji,jj)
-         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
+      DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
+         iik = nmln(ji,jj)
+         hmlp (ji,jj) = gdepw(ji,jj,iik  ,Kmm) * ssmask(ji,jj)    ! Mixed layer depth
+         hmlpt(ji,jj) = gdept(ji,jj,iik-1,Kmm) * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer
       END_2D
       !
-      IF( .NOT.l_offline ) THEN
-         IF( iom_use("mldr10_1") ) THEN
-            IF( ln_isfcav ) THEN  ;  CALL iom_put( "mldr10_1", hmlp - risfdep)   ! mixed layer thickness
-            ELSE                  ;  CALL iom_put( "mldr10_1", hmlp )            ! mixed layer depth
-            END IF
+      IF( .NOT.l_offline .AND. iom_use("mldr10_1") ) THEN
+         IF( ln_isfcav ) THEN  ;  CALL iom_put( "mldr10_1", hmlp - risfdep)   ! mixed layer thickness
+         ELSE                  ;  CALL iom_put( "mldr10_1", hmlp )            ! mixed layer depth
          END IF
-         IF( iom_use("mldkz5") ) THEN
-            IF( ln_isfcav ) THEN  ;  CALL iom_put( "mldkz5"  , hmld - risfdep )   ! turbocline thickness
-            ELSE                  ;  CALL iom_put( "mldkz5"  , hmld )             ! turbocline depth
-            END IF
-         ENDIF
       ENDIF
       !
@@ -140,4 +122,44 @@
    END SUBROUTINE zdf_mxl
 
+
+   SUBROUTINE zdf_mxl_turb( kt, Kmm )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_mxl_turb  ***
+      !!
+      !! ** Purpose :   Compute the turbocline depth.
+      !!
+      !! ** Method  :   The turbocline depth is the depth at which the vertical
+      !!      eddy diffusivity coefficient (resulting from the vertical physics
+      !!      alone, not the isopycnal part, see trazdf.F) fall below a given
+      !!      value defined locally (avt_c here taken equal to 5 cm/s2 by default)
+      !!
+      !! ** Action  :   hmld
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
+      !
+      INTEGER  ::   ji, jj, jk      ! dummy loop indices
+      INTEGER  ::   iik             ! local integer
+      INTEGER, DIMENSION(A2D(nn_hls)) ::   imld   ! 2D workspace
+      !!----------------------------------------------------------------------
+      !
+      ! w-level of the turbocline and mixing layer (iom_use)
+      imld(:,:) = mbkt(A2D(nn_hls)) + 1                ! Initialization to the number of w ocean point
+      DO_3DS( 0, 0, 0, 0, jpkm1, nlb10, -1 )   ! from the bottom to nlb10
+         IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) )   imld(ji,jj) = jk      ! Turbocline
+      END_3D
+      ! depth of the mixing layer
+      DO_2D( 0, 0, 0, 0 )
+         iik = imld(ji,jj)
+         hmld (ji,jj) = gdepw(ji,jj,iik  ,Kmm) * ssmask(ji,jj)    ! Turbocline depth
+      END_2D
+      !
+      IF( .NOT.l_offline .AND. iom_use("mldkz5") ) THEN
+         IF( ln_isfcav ) THEN  ;  CALL iom_put( "mldkz5"  , hmld - risfdep )   ! turbocline thickness
+         ELSE                  ;  CALL iom_put( "mldkz5"  , hmld )             ! turbocline depth
+         END IF
+      ENDIF
+      !
+   END SUBROUTINE zdf_mxl_turb
    !!======================================================================
 END MODULE zdfmxl
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfosm.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfosm.F90	(revision 14664)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfosm.F90	(revision 14665)
@@ -1214,5 +1214,6 @@
         ! Lateral boundary conditions on final outputs for hbl,  on T-grid (sign unchanged)
         ! [comm_cleanup] ! no need lbc_lnk for output
-        ! CALL lbc_lnk( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. )
+        ! NOTE: [tiling] still needed (at least for hmle in tra_mle_trp)
+         CALL lbc_lnk( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. )
         ! Lateral boundary conditions on final outputs for gham[ts],  on W-grid  (sign unchanged)
         ! Lateral boundary conditions on final outputs for gham[uv],  on [UV]-grid  (sign changed)
@@ -2874,5 +2875,5 @@
       !
       IF( kt == nit000 ) THEN
-         IF( ntile == 0 .OR. ntile == 1 ) THEN                    ! Do only on the first tile
+         IF( .NOT. l_istiled .OR. ntile == 1 ) THEN                    ! Do only on the first tile
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes'
@@ -2886,6 +2887,5 @@
       ENDIF
 
-      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          pts(ji,jj,jk,jp_tem,Krhs) =  pts(ji,jj,jk,jp_tem,Krhs)                      &
             &                 - (  ghamt(ji,jj,jk  )  &
@@ -2948,12 +2948,13 @@
       !
       IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity'
-         IF(lwp) WRITE(numout,*) '~~~~~~~   '
+         IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity'
+            IF(lwp) WRITE(numout,*) '~~~~~~~   '
+         ENDIF
       ENDIF
       !code saving tracer trends removed, replace with trdmxl_oce
 
-      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )       ! add non-local u and v fluxes
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )       ! add non-local u and v fluxes
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )       ! add non-local u and v fluxes
          puu(ji,jj,jk,Krhs) =  puu(ji,jj,jk,Krhs)                      &
             &                 - (  ghamu(ji,jj,jk  )  &
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfphy.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfphy.F90	(revision 14664)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfphy.F90	(revision 14665)
@@ -12,4 +12,6 @@
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers variables
+   ! TEMP: [tiling] This change not necessary after finalisation of zdf_osm
+   USE domtile
    USE zdf_oce        ! vertical physics: shared variables
    USE zdfdrg         ! vertical physics: top/bottom drag coef.
@@ -56,4 +58,8 @@
    LOGICAL, PUBLIC ::   l_zdfsh2   ! shear production term flag (=F for CST, =T otherwise (i.e. TKE, GLS, RIC))
 
+   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avm_k_n !: "Now" avm_k used for calculation of zsh2 with tiling
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -210,7 +216,8 @@
       ENDIF
       !                                ! shear production term flag
-      IF( ln_zdfcst ) THEN   ;   l_zdfsh2 = .FALSE.
-      ELSE                   ;   l_zdfsh2 = .TRUE.
-      ENDIF
+      IF( ln_zdfcst .OR. ln_zdfosm ) THEN   ;   l_zdfsh2 = .FALSE.
+      ELSE                                  ;   l_zdfsh2 = .TRUE.
+      ENDIF
+      IF( ln_tile .AND. l_zdfsh2 ) ALLOCATE( avm_k_n(jpi,jpj,jpk) )
       !                          !== Mass Flux Convectiive algorithm  ==!
       IF( ln_zdfmfc )   CALL zdf_mfc_init       ! Convection computed with eddy diffusivity mass flux
@@ -246,8 +253,21 @@
       !
       INTEGER ::   ji, jj, jk   ! dummy loop indice
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsh2   ! shear production
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zsh2   ! shear production
+      ! TEMP: [tiling] This change not necessary after finalisation of zdf_osm
+      LOGICAL :: lskip
       !! ---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('zdf_phy')
+
+      ! TEMP: [tiling] These changes not necessary after finalisation of zdf_osm
+      lskip = .FALSE.
+
+      IF( ln_tile .AND. nzdf_phy == np_OSM )  THEN
+         IF( ntile == 1 ) THEN
+            CALL dom_tile_stop( ldhold=.TRUE., cstr='zdfphy' )
+         ELSE
+            lskip = .TRUE.
+         ENDIF
+      ENDIF
       !
       IF( l_zdfdrg ) THEN     !==  update top/bottom drag  ==!   (non-linear cases)
@@ -264,42 +284,64 @@
       ENDIF
       !
+! TODO: [tiling] NOT TESTED- requires SI3 and isfcav
 #if defined key_si3
       IF ( ln_drgice_imp) THEN
          IF ( ln_isfcav ) THEN
-            rCdU_top(:,:) = rCdU_top(:,:) + ssmask(:,:) * tmask(:,:,1) * rCdU_ice(:,:)
+            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+               rCdU_top(ji,jj) = rCdU_top(ji,jj) + ssmask(ji,jj) * tmask(ji,jj,1) * rCdU_ice(ji,jj)
+            END_2D
          ELSE
-            rCdU_top(:,:) = rCdU_ice(:,:)
+            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+               rCdU_top(ji,jj) = rCdU_ice(ji,jj)
+            END_2D
          ENDIF
       ENDIF
 #endif
       !
-      !                       !==  Kz from chosen turbulent closure  ==!   (avm_k, avt_k)
-      !
-      IF( l_zdfsh2 )   &         !* shear production at w-points (energy conserving form)
-         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,      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
-      ! [comm_cleanup] ! modified but not tested - no ref config uses this scheme
-      CASE( np_OSM )   ;   CALL zdf_osm( kt, Kbb, Kmm, Krhs, 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
+      ! NOTE: [tiling] zdf_mxl has been split into zdf_mxl (nmln, hmlp, hmlpt) and zdf_mxl_turb (hmld), which depend on rn2b and avt respectively. zdf_mxl has been moved here from below (which is now zdf_mxl_turb) because zdf_tke with nn_etau = 2 depends on nmln. The previous solution was to call zdf_mxl in zdf_tke_init, but rn2b is not yet calculated here (0 everywhere), so this specific option was not restartable and tiling changed the results. Furthermore, the MLD diagnostics sent from zdf_mxl were incorrect for the first timestep. This bug fix has therefore changed the results when using zdf_tke with nn_etau = 2.
+      CALL zdf_mxl( kt, Kmm )                        !* mixed layer depth, and level
+
+      ! TEMP: [tiling] These changes not necessary after finalisation of zdf_osm
+      IF( .NOT. lskip ) THEN
+         !                       !==  Kz from chosen turbulent closure  ==!   (avm_k, avt_k)
+         !
+         ! NOTE: [tiling] the closure schemes (zdf_tke etc) will update avm_k. With tiling, the calculation of zsh2 on adjacent tiles then uses both updated (next timestep) and non-updated (current timestep) values of avm_k. To preserve results, we save a read-only copy of the "now" avm_k to use in the calculation of zsh2.
+         IF( l_zdfsh2 ) THEN        !* shear production at w-points (energy conserving form)
+            IF( ln_tile ) THEN
+               IF( ntile == 1 ) avm_k_n(:,:,:) = avm_k(:,:,:)     ! Preserve "now" avm_k for calculation of zsh2
+               CALL zdf_sh2( Kbb, Kmm, avm_k_n, &     ! <<== in
+                  &                     zsh2    )     ! ==>> out : shear production
+            ELSE
+               CALL zdf_sh2( Kbb, Kmm, avm_k,   &     ! <<== in
+                  &                     zsh2    )     ! ==>> out : shear production
+            ENDIF
+         ENDIF
+         !
+         SELECT CASE ( nzdf_phy )                  !* Vertical eddy viscosity and diffusivity coefficients at w-points
+         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, Krhs, 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
 !!gm         avt(2:jpim1,2:jpjm1,1:jpkm1) = rn_avt0 * wmask(2:jpim1,2:jpjm1,1:jpkm1)
 !!gm         avm(2:jpim1,2:jpjm1,1:jpkm1) = rn_avm0 * wmask(2:jpim1,2:jpjm1,1:jpkm1)
-      END SELECT
+         END SELECT
+
+         IF( ln_tile .AND. .NOT. l_istiled ) CALL dom_tile_start( ldhold=.TRUE., cstr='zdfphy' )
+      ENDIF
       !
       !                          !==  ocean Kz  ==!   (avt, avs, avm)
       !
       !                                         !* start from turbulent closure values
-      avt(:,:,2:jpkm1) = avt_k(:,:,2:jpkm1)
-      avm(:,:,2:jpkm1) = avm_k(:,:,2:jpkm1)
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
+         avt(ji,jj,jk) = avt_k(ji,jj,jk)
+         avm(ji,jj,jk) = avm_k(ji,jj,jk)
+      END_3D
       !
       IF( ln_rnf_mouth ) THEN                   !* increase diffusivity at rivers mouths
-         DO jk = 2, nkrnf
-            avt(:,:,jk) = avt(:,:,jk) + 2._wp * rn_avt_rnf * rnfmsk(:,:) * wmask(:,:,jk)
-         END DO
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, nkrnf )
+            avt(ji,jj,jk) = avt(ji,jj,jk) + 2._wp * rn_avt_rnf * rnfmsk(ji,jj) * wmask(ji,jj,jk)
+         END_3D
       ENDIF
       !
@@ -310,8 +352,11 @@
                         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)
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
+            avs(ji,jj,jk) = avt(ji,jj,jk)
+         END_3D
       ENDIF
       !
       !                                         !* wave-induced mixing
+      ! TODO: [tiling] NOT TESTED- requires ln_wave and ln_sdw
       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)
@@ -319,31 +364,39 @@
 #if defined key_agrif
       ! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2)
-      IF( l_zdfsh2 )   CALL Agrif_avm
+      ! TODO: [tiling] NOT TESTED- requires AGRIF
+      IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
+         IF( l_zdfsh2 )   CALL Agrif_avm
+      ENDIF
 #endif
 
       !                                         !* Lateral boundary conditions (sign unchanged)
       ! [comm_cleanup] ! lbc_lnk shifted in stp
-      IF(nn_hls.eq.1) THEN 
-         IF( l_zdfsh2 ) THEN
-            CALL lbc_lnk( 'zdfphy', avm_k, 'W', 1.0_wp , avt_k, 'W', 1.0_wp,   &
-               &                    avm  , 'W', 1.0_wp , avt  , 'W', 1.0_wp , avs , 'W', 1.0_wp )
-         ELSE
-            CALL lbc_lnk( 'zdfphy', avm  , 'W', 1.0_wp , avt  , 'W', 1.0_wp , avs , 'W', 1.0_wp )
-         ENDIF
-      !
-         IF( l_zdfdrg ) THEN     ! drag  have been updated (non-linear cases)
-            IF( ln_isfcav ) THEN   ;  CALL lbc_lnk( 'zdfphy', rCdU_top, 'T', 1.0_wp , rCdU_bot, 'T', 1.0_wp )   ! top & bot drag
-            ELSE                   ;  CALL lbc_lnk( 'zdfphy', rCdU_bot, 'T', 1.0_wp )                       ! bottom drag only
+      IF(nn_hls.eq.1) THEN
+         ! NOTE: [tiling] this solution for lbc_lnk will not work if outer halo values are accessed before the end of the tiling loop
+         IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
+            IF( l_zdfsh2 ) THEN
+               CALL lbc_lnk( 'zdfphy', avm_k, 'W', 1.0_wp , avt_k, 'W', 1.0_wp,   &
+                     &                 avm  , 'W', 1.0_wp , avt  , 'W', 1.0_wp , avs , 'W', 1.0_wp )
+            ELSE
+               CALL lbc_lnk( 'zdfphy', avm  , 'W', 1.0_wp , avt  , 'W', 1.0_wp , avs , 'W', 1.0_wp )
             ENDIF
-         ENDIF
-      ENDIF
-      !
-      CALL zdf_mxl( kt, Kmm )                        !* mixed layer depth, and level
-      !
-      IF( lrst_oce ) THEN                       !* write TKE, GLS or RIC fields in the restart file
-         IF( ln_zdftke )   CALL tke_rst( kt, 'WRITE' )
-         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 ww has been updated
+            !
+            IF( l_zdfdrg ) THEN     ! drag  have been updated (non-linear cases)
+               IF( ln_isfcav ) THEN   ;  CALL lbc_lnk( 'zdfphy', rCdU_top, 'T', 1.0_wp , rCdU_bot, 'T', 1.0_wp )   ! top & bot drag
+               ELSE                   ;  CALL lbc_lnk( 'zdfphy', rCdU_bot, 'T', 1.0_wp )                           ! bottom drag only
+               ENDIF
+            ENDIF
+         ENDIF
+      ENDIF
+      !
+      CALL zdf_mxl_turb( kt, Kmm )                   !* turbocline depth
+      !
+      IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
+         IF( lrst_oce ) THEN                       !* write TKE, GLS or RIC fields in the restart file
+            IF( ln_zdftke )   CALL tke_rst( kt, 'WRITE' )
+            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 ww has been updated
+         ENDIF
       ENDIF
       !
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfric.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfric.F90	(revision 14664)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfric.F90	(revision 14665)
@@ -51,4 +51,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -145,17 +146,17 @@
       !!              PFJ Lermusiaux 2001.
       !!----------------------------------------------------------------------
-      INTEGER                   , INTENT(in   ) ::   kt             ! ocean time-step
-      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)
+      INTEGER                             , INTENT(in   ) ::   kt             ! ocean time-step
+      INTEGER                             , INTENT(in   ) ::   Kmm            ! ocean time level index
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   p_sh2          ! shear production term
+      REAL(wp), DIMENSION(:,:,:)          , INTENT(inout) ::   p_avm, p_avt   ! momentum and tracer Kz (w-points)
       !!
       INTEGER  ::   ji, jj, jk                  ! dummy loop indices
       REAL(wp) ::   zcfRi, zav, zustar, zhek    ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::   zh_ekm  ! 2D workspace
+      REAL(wp), DIMENSION(A2D(nn_hls)) ::   zh_ekm  ! 2D workspace
       !!----------------------------------------------------------------------
       !
       !                       !==  avm and avt = F(Richardson number)  ==!
       ! [comm_cleanup] ! DO_3D( 1, 0, 1, 0, 2, jpkm1 )       ! coefficient = F(richardson number) (avm-weighted Ri)
-      DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 2, jpkm1 )       ! coefficient = F(richardson number) (avm-weighted Ri)
+      DO_3D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 2, jpkm1 )       ! coefficient = F(richardson number) (avm-weighted Ri)
          zcfRi = 1._wp / (  1._wp + rn_alp * MAX(  0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( p_sh2(ji,jj,jk) + 1.e-20 ) )  )
          zav   = rn_avmri * zcfRi**nn_ric
@@ -177,5 +178,5 @@
          END_2D
          ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )   !* minimum mixing coeff. within the Ekman layer
-         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )   !* minimum mixing coeff. within the Ekman layer
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )   !* minimum mixing coeff. within the Ekman layer
             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)
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfsh2.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfsh2.F90	(revision 14664)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfsh2.F90	(revision 14665)
@@ -55,13 +55,14 @@
       !! References :   Bruchard, OM 2002
       !! ---------------------------------------------------------------------
-      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)
+      INTEGER                              , INTENT(in   ) ::   Kbb, Kmm             ! ocean time level indices
+      REAL(wp), DIMENSION(:,:,:)           , INTENT(in   ) ::   p_avm                ! vertical eddy viscosity (w-points)
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) , INTENT(  out) ::   p_sh2                ! shear production of TKE (w-points)
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop arguments
-      REAL(wp), DIMENSION(jpi,jpj) ::   zsh2u, zsh2v   ! 2D workspace
+      REAL(wp), DIMENSION(A2D(nn_hls)) ::   zsh2u, zsh2v   ! 2D workspace
       !!--------------------------------------------------------------------
       !
       DO jk = 2, jpkm1                 !* Shear production at uw- and vw-points (energy conserving form)
+         ! TODO: [tiling] NOT TESTED- requires key_oasis3
          IF ( cpl_sdrftx .AND. ln_stshear )  THEN       ! Surface Stokes Drift available  ===>>>  shear + stokes drift contibution
             ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfswm.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfswm.F90	(revision 14664)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfswm.F90	(revision 14665)
@@ -64,5 +64,5 @@
       zcoef = 1._wp * 0.353553_wp
       ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
          zqb = zcoef * hsw(ji,jj) * tsd2d(ji,jj) * EXP( -3. * wnum(ji,jj) * gdepw(ji,jj,jk,Kmm) ) * wmask(ji,jj,jk)
          !
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdftke.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdftke.F90	(revision 14664)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdftke.F90	(revision 14665)
@@ -168,8 +168,8 @@
       !!              Bruchard OM 2002
       !!----------------------------------------------------------------------
-      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)
+      INTEGER                             , INTENT(in   ) ::   kt             ! ocean time step
+      INTEGER                             , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   p_sh2          ! shear production term
+      REAL(wp), DIMENSION(:,:,:)          , INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points)
       !!----------------------------------------------------------------------
       !
@@ -201,7 +201,7 @@
       USE zdf_oce , ONLY : en   ! ocean vertical physics
       !!
-      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)
+      INTEGER                              , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) , INTENT(in   ) ::   p_sh2          ! shear production term
+      REAL(wp), DIMENSION(:,:,:)           , INTENT(in   ) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points)
       !
       INTEGER ::   ji, jj, jk                  ! dummy loop arguments
@@ -216,7 +216,7 @@
       REAL(wp) ::   zzd_up, zzd_lw             !   -      -
       REAL(wp) ::   ztaui, ztauj, z1_norm
-      INTEGER , DIMENSION(jpi,jpj)     ::   imlc
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra, zhlc, zus3, zWlc2
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc, zdiag, zd_up, zd_lw
+      INTEGER , DIMENSION(A2D(nn_hls))     ::   imlc
+      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zice_fra, zhlc, zus3, zWlc2
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zpelc, zdiag, zd_up, zd_lw
       !!--------------------------------------------------------------------
       !
@@ -232,7 +232,7 @@
       SELECT CASE ( nn_eice )
       CASE( 0 )   ;   zice_fra(:,:) = 0._wp
-      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(:,:) * 10._wp )
-      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(:,:)
-      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp )
+      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(A2D(nn_hls)) * 10._wp )
+      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(A2D(nn_hls))
+      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(A2D(nn_hls)) , 1._wp )
       END SELECT
       !
@@ -242,5 +242,5 @@
       !
       ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) )
          zdiag(ji,jj,1) = 1._wp/en(ji,jj,1)
@@ -260,5 +260,5 @@
          !
          ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )        ! bottom friction
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )        ! bottom friction
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )        ! bottom friction
             zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
             zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
@@ -270,5 +270,5 @@
          IF( ln_isfcav ) THEN
             ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )     ! top friction
-            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )     ! top friction 
+            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )     ! top friction
                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)) )
@@ -290,4 +290,5 @@
          !                       !* Langmuir velocity scale
          !
+         ! TODO: [tiling] NOT TESTED- requires key_oasis3
          IF ( cpl_sdrftx )  THEN       ! Surface Stokes Drift available
             !                                ! Craik-Leibovich velocity scale Wlc = ( u* u_s )^1/2    with u* = (taum/rho0)^1/2
@@ -312,7 +313,4 @@
                   zWlc2(ji,jj) = 0.5_wp * z1_norm * ( MAX( ut0sd(ji,jj)*ztaui + vt0sd(ji,jj)*ztauj, 0._wp ) )**2
             END_2D
-            ! [comm_cleanup]
-            IF (nn_hls.eq.1) CALL lbc_lnk      ( 'zdftke', zWlc2, 'T', 1. )
-!
          ELSE                          ! Surface Stokes drift deduced from surface stress
             !                                ! Wlc = u_s   with u_s = 0.016*U_10m, the surface stokes drift  (Axell 2002, Eq.44)
@@ -322,5 +320,5 @@
             zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )      ! to convert stress in 10m wind using a constant drag
             ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
-            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                zWlc2(ji,jj) = zcof * taum(ji,jj)
             END_2D
@@ -330,19 +328,21 @@
          !                       !* Depth of the LC circulation  (Axell 2002, Eq.47)
          !                             !- LHS of Eq.47
-         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 ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm)
-         END DO
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            zpelc(ji,jj,1) =  MAX( rn2b(ji,jj,1), 0._wp ) * gdepw(ji,jj,1,Kmm) * e3w(ji,jj,1,Kmm)
+         END_2D
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpk )
+            zpelc(ji,jj,jk)  = zpelc(ji,jj,jk-1) +   &
+               &          MAX( rn2b(ji,jj,jk), 0._wp ) * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
+         END_3D
          !
          !                             !- compare LHS to RHS of Eq.47
-         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
+         imlc(:,:) = mbkt(A2D(nn_hls)) + 1       ! Initialization to the number of w ocean point (=2 over land)
          ! [comm_cleanup] ! DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )
-         DO_3DS( nn_hls, nn_hls, nn_hls, nn_hls, jpkm1, 2, -1 )
+         DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 )
             IF( zpelc(ji,jj,jk) > zWlc2(ji,jj) )   imlc(ji,jj) = jk
          END_3D
          !                               ! finite LC depth
          ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm)
          END_2D
@@ -355,5 +355,5 @@
          END_2D
          ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )                  !* TKE Langmuir circulation source term added to en
-         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                  !* TKE Langmuir circulation source term added to en
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                  !* TKE Langmuir circulation source term added to en
             IF ( zus3(ji,jj) /= 0._wp ) THEN
                IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN
@@ -377,5 +377,5 @@
       IF( nn_pdl == 1 ) THEN          !* Prandtl number = F( Ri )
          ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             !                             ! local Richardson number
             IF (rn2b(ji,jj,jk) <= 0.0_wp) then
@@ -390,5 +390,5 @@
       !
       ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )   !* Matrix and right hand side in en
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )   !* Matrix and right hand side in en
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )   !* Matrix and right hand side in en
          zcof   = zfact1 * tmask(ji,jj,jk)
          !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical
@@ -415,4 +415,5 @@
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !
+      ! TODO: [tiling] NOT TESTED- requires key_oasis3 and waves
       IF ( cpl_phioc .and. ln_phioc )  THEN
          SELECT CASE (nn_bc_surf) ! Boundary Condition using surface TKE flux from waves
@@ -420,5 +421,5 @@
          CASE ( 0 ) ! Dirichlet BC
             ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )    ! en(1)   = rn_ebb taum / rho0  (min value rn_emin0)
-            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )    ! en(1)   = rn_ebb taum / rho0  (min value rn_emin0)
+            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )    ! en(1)   = rn_ebb taum / rho0  (min value rn_emin0)
                IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp
                en(ji,jj,1) = MAX( rn_emin0, .5 * ( 15.8 * phioc(ji,jj) / rho0 )**(2./3.) )  * tmask(ji,jj,1)
@@ -428,5 +429,5 @@
          CASE ( 1 ) ! Neumann BC
             ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp
                en(ji,jj,2)    = en(ji,jj,2) + ( rn_Dt * phioc(ji,jj) / rho0 ) /e3w(ji,jj,2,Kmm)
@@ -455,13 +456,13 @@
       END_3D
       ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )                          ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                          ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                          ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
          en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
       END_2D
       ! [comm_cleanup] ! DO_3DS( 0, 0, 0, 0, jpk-2, 2, -1 )
-      DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, 2, -1 )
+      DO_3DS_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, 2, -1 )
          en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
       END_3D
       ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! set the minimum value of tke
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! set the minimum value of tke
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! set the minimum value of tke
          en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
       END_3D
@@ -477,5 +478,5 @@
       IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
          ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             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._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
@@ -483,12 +484,13 @@
       ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
          ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             jk = nmln(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._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
          END_2D
+      ! TODO: [tiling] NOT TESTED- requires ln_cpl
       ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
          ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
             zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
@@ -547,5 +549,5 @@
       REAL(wp) ::   zdku,   zdkv, zsqen       !   -      -
       REAL(wp) ::   zemxl, zemlm, zemlp, zmaxice       !   -      -
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmxlm, zmxld   ! 3D workspace
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zmxlm, zmxld   ! 3D workspace
       !!--------------------------------------------------------------------
       !
@@ -560,4 +562,5 @@
       zmxld(:,:,:)  = rmxl_min
       !
+      ! TODO: [tiling] NOT TESTED- requires waves
       IF(ln_sdw .AND. ln_mxhsw) THEN
          zmxlm(:,:,1)= vkarmn * MAX ( 1.6 * hsw(:,:) , 0.02 )        ! surface mixing length = F(wave height)
@@ -576,4 +579,5 @@
             END_2D
 #else
+            ! TODO: [tiling] NOT TESTED- requires sea ice
             SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice
             !
@@ -698,5 +702,5 @@
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !* vertical eddy viscosity & diffivity at w-points
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )   !* vertical eddy viscosity & diffivity at w-points
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )   !* vertical eddy viscosity & diffivity at w-points
          zsqen = SQRT( en(ji,jj,jk) )
          zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
@@ -709,5 +713,5 @@
       IF( nn_pdl == 1 ) THEN          !* Prandtl number case: update avt
          ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             p_avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
          END_3D
@@ -825,7 +829,7 @@
       !
       !                               !* Check of some namelist values
-      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
-      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
-      IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
+      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1, 2 or 3' )
+      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1' )
+      IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0 or 1' )
       IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
       !
@@ -835,7 +839,4 @@
          rn_mxl0 = rmxl_min
       ENDIF
-
-      IF( nn_etau == 2  )   CALL zdf_mxl( nit000, Kmm )      ! Initialization of nmln
-
       !                               !* depth of penetration of surface tke
       IF( nn_etau /= 0 ) THEN
