Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/LDF/ldfslp.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/LDF/ldfslp.F90	(revision 14575)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/LDF/ldfslp.F90	(revision 14576)
@@ -371,5 +371,6 @@
          !
          ip = jl   ;   jp = jl                ! guaranteed nonzero gradients ( absolute value larger than repsln)
-         DO_3D( 1, 0, 1, 0, 1, jpkm1 )        ! done each pair of triad ! NB: not masked ==>  a minimum value is set
+         ! [comm_cleanup] ! DO_3D( 1, 0, 1, 0, 1, jpkm1 )        ! done each pair of triad ! NB: not masked ==>  a minimum value is set
+         DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )        ! done each pair of triad ! NB: not masked ==>  a minimum value is set
             zdit = ( ts(ji+1,jj,jk,jp_tem,Kbb) - ts(ji,jj,jk,jp_tem,Kbb) )    ! i-gradient of T & S at u-point
             zdis = ( ts(ji+1,jj,jk,jp_sal,Kbb) - ts(ji,jj,jk,jp_sal,Kbb) )
@@ -383,5 +384,6 @@
          !
          IF( ln_zps .AND. l_grad_zps ) THEN     ! partial steps: correction of i- & j-grad on bottom
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                iku  = mbku(ji,jj)          ;   ikv  = mbkv(ji,jj)             ! last ocean level (u- & v-points)
                zdit = gtsu(ji,jj,jp_tem)   ;   zdjt = gtsv(ji,jj,jp_tem)      ! i- & j-gradient of Temperature
@@ -397,5 +399,6 @@
 
       DO kp = 0, 1                            !==  unmasked before density i- j-, k-gradients  ==!
-         DO_3D( 1, 1, 1, 1, 1, jpkm1 )        ! done each pair of triad ! NB: not masked ==>  a minimum value is set
+         ! [comm_cleanup] ! DO_3D( 1, 1, 1, 1, 1, jpkm1 )        ! done each pair of triad ! NB: not masked ==>  a minimum value is set
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )      ! done each pair of triad ! NB: not masked ==>  a minimum value is set
             IF( jk+kp > 1 ) THEN              ! k-gradient of T & S a jk+kp
                zdkt = ( ts(ji,jj,jk+kp-1,jp_tem,Kbb) - ts(ji,jj,jk+kp,jp_tem,Kbb) )
@@ -412,5 +415,6 @@
       END DO
       !
-      DO_2D( 1, 1, 1, 1 )                     !==  Reciprocal depth of the w-point below ML base  ==!
+      ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )                     !==  Reciprocal depth of the w-point below ML base  ==!
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )                   !== Reciprocal depth of the w-point below ML base  ==!
          jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1     ! MIN in case ML depth is the ocean depth
          z1_mlbw(ji,jj) = 1._wp / gdepw(ji,jj,jk,Kmm)
@@ -432,5 +436,6 @@
       DO jl = 0, 1                            ! calculate slope of the 4 triads immediately ONE level below mixed-layer base
          DO kp = 0, 1                         ! with only the slope-max limit   and   MASKED
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                ip = jl   ;   jp = jl
                !
@@ -469,5 +474,6 @@
                ! Must mask contribution to slope from dz/dx at constant s for triads jk=1,kp=0 that poke up though ocean surface
                znot_thru_surface = REAL( 1-1/(jk+kp), wp )  !jk+kp=1,=0.; otherwise=1.0
-               DO_2D( 1, 0, 1, 0 )
+               ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+               DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                   !
                   ! Calculate slope relative to geopotentials used for GM skew fluxes
@@ -552,4 +558,12 @@
       CALL lbc_lnk( 'ldfslp', wslp2, 'W', 1.0_wp )      ! lateral boundary confition on wslp2 only   ==>>> gm : necessary ? to be checked
       !
+      ! [comm_cleanup] - this comm have to be deleted !!!!!
+      IF (nn_hls.eq.2) THEN
+         CALL lbc_lnk( 'ldfslp', triadi(:,:,:,:,0), 'U', 1.0_wp, triadi(:,:,:,:,1), 'U', 1.0_wp )
+         CALL lbc_lnk( 'ldfslp', triadj(:,:,:,:,0), 'V', 1.0_wp, triadj(:,:,:,:,1), 'V', 1.0_wp )
+         CALL lbc_lnk( 'ldfslp', triadi_g(:,:,:,:,0), 'U', 1.0_wp, triadi_g(:,:,:,:,1), 'U', 1.0_wp )
+         CALL lbc_lnk( 'ldfslp', triadj_g(:,:,:,:,0), 'V', 1.0_wp, triadj_g(:,:,:,:,1), 'V', 1.0_wp )
+      END IF
+
       IF( ln_timing )   CALL timing_stop('ldf_slp_triad')
       !
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/LDF/ldftra.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/LDF/ldftra.F90	(revision 14575)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/LDF/ldftra.F90	(revision 14576)
@@ -647,5 +647,6 @@
       !                       ! Compute lateral diffusive coefficient at T-point
       IF( ln_traldf_triad ) THEN
-         DO_3D( 0, 0, 0, 0, 1, jpk )
+         ! [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 )
             ! Take the max of N^2 and zero then take the vertical sum
             ! of the square root of the resulting N^2 ( required to compute
@@ -661,5 +662,6 @@
          END_3D
       ELSE
-         DO_3D( 0, 0, 0, 0, 1, jpk )
+         ! [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 )
             ! Take the max of N^2 and zero then take the vertical sum
             ! of the square root of the resulting N^2 ( required to compute
@@ -677,5 +679,6 @@
       ENDIF
 
-      DO_2D( 0, 0, 0, 0 )
+      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
+      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 )
          ! Rossby radius at w-point taken betwenn 2 km and  40km
@@ -687,5 +690,6 @@
       !                                         !==  Bound on eiv coeff.  ==!
       z1_f20 = 1._wp / (  2._wp * omega * sin( rad * 20._wp )  )
-      DO_2D( 0, 0, 0, 0 )
+      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
+      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj)     ! tropical decrease
          zaeiw(ji,jj) = MIN( zzaei , paei0 )                                  ! Max value = paei0
@@ -693,5 +697,6 @@
       CALL lbc_lnk( 'ldftra', zaeiw(:,:), 'W', 1.0_wp )       ! lateral boundary condition
       !
-      DO_2D( 0, 0, 0, 0 )                       !== aei at u- and v-points  ==!
+      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
+      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          paeiu(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji+1,jj  ) ) * umask(ji,jj,1)
          paeiv(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji  ,jj+1) ) * vmask(ji,jj,1)
@@ -751,5 +756,6 @@
       zpsi_uw(:,:,jpk) = 0._wp   ;   zpsi_vw(:,:,jpk) = 0._wp
       !
-      DO_3D( 1, 0, 1, 0, 2, jpkm1 )
+      ! [comm_cleanup] ! DO_3D( 1, 0, 1, 0, 2, jpkm1 )
+      DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 2, jpkm1 )
          zpsi_uw(ji,jj,jk) = - r1_4 * e2u(ji,jj) * ( wslpi(ji,jj,jk  ) + wslpi(ji+1,jj,jk) )   &
             &                                    * ( aeiu (ji,jj,jk-1) + aeiu (ji  ,jj,jk) ) * wumask(ji,jj,jk)
@@ -758,9 +764,11 @@
       END_3D
       !
-      DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+      ! [comm_cleanup] ! DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+      DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )
          pu(ji,jj,jk) = pu(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
          pv(ji,jj,jk) = pv(ji,jj,jk) - ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
       END_3D
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+      ! [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 )
          pw(ji,jj,jk) = pw(ji,jj,jk) + (  zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj  ,jk)   &
             &                             + zpsi_vw(ji,jj,jk) - zpsi_vw(ji  ,jj-1,jk) )
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traadv.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traadv.F90	(revision 14575)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traadv.F90	(revision 14576)
@@ -61,12 +61,12 @@
    LOGICAL ::   ln_traadv_qck    ! QUICKEST scheme flag
 
-   INTEGER ::   nadv             ! choice of the type of advection scheme
+   INTEGER, PUBLIC ::   nadv             ! choice of the type of advection scheme
    !                             ! associated indices:
-   INTEGER, PARAMETER ::   np_NO_adv  = 0   ! no T-S advection
-   INTEGER, PARAMETER ::   np_CEN     = 1   ! 2nd/4th order centered scheme
-   INTEGER, PARAMETER ::   np_FCT     = 2   ! 2nd/4th order Flux Corrected Transport scheme
-   INTEGER, PARAMETER ::   np_MUS     = 3   ! MUSCL scheme
-   INTEGER, PARAMETER ::   np_UBS     = 4   ! 3rd order Upstream Biased Scheme
-   INTEGER, PARAMETER ::   np_QCK     = 5   ! QUICK scheme
+   INTEGER, PARAMETER, PUBLIC ::   np_NO_adv  = 0   ! no T-S advection
+   INTEGER, PARAMETER, PUBLIC ::   np_CEN     = 1   ! 2nd/4th order centered scheme
+   INTEGER, PARAMETER, PUBLIC ::   np_FCT     = 2   ! 2nd/4th order Flux Corrected Transport scheme
+   INTEGER, PARAMETER, PUBLIC ::   np_MUS     = 3   ! MUSCL scheme
+   INTEGER, PARAMETER, PUBLIC ::   np_UBS     = 4   ! 3rd order Upstream Biased Scheme
+   INTEGER, PARAMETER, PUBLIC ::   np_QCK     = 5   ! QUICK scheme
 
    !! * Substitutions
@@ -178,10 +178,12 @@
          !
          CASE ( np_CEN )                                 ! Centered scheme : 2nd / 4th order
-            IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kmm), 'T', 1. )
+            ! [comm_cleanup]
+            ! IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kmm), 'T', 1. )
             CALL tra_adv_cen    ( kt, nit000, 'TRA',         zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v )
          CASE ( np_FCT )                                 ! FCT scheme      : 2nd / 4th order
             IF (nn_hls.EQ.2) THEN
-               CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1., pts(:,:,:,:,Kmm), 'T', 1.)
-               CALL lbc_lnk( 'traadv', zuu(:,:,:), 'U', -1., zvv(:,:,:), 'V', -1., zww(:,:,:), 'W', 1.)
+               ! [comm_cleanup] - lbc_lnk shifted into step
+               ! CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1., pts(:,:,:,:,Kmm), 'T', 1.)
+               ! CALL lbc_lnk( 'traadv', zuu(:,:,:), 'U', -1., zvv(:,:,:), 'V', -1., zww(:,:,:), 'W', 1.)
 #if defined key_loop_fusion
                CALL tra_adv_fct_lf ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v )
@@ -194,5 +196,6 @@
          CASE ( np_MUS )                                 ! MUSCL
             IF (nn_hls.EQ.2) THEN
-                CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)
+               ! [comm_cleanup] - lbc_lnk shifted into step
+               ! CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)
 #if defined key_loop_fusion
                 CALL tra_adv_mus_lf ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )
@@ -204,11 +207,13 @@
             END IF
          CASE ( np_UBS )                                 ! UBS
-            IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)
+            ! [comm_cleanup] - lbc_lnk shifted into step
+            ! IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)
             CALL tra_adv_ubs    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v   )
          CASE ( np_QCK )                                 ! QUICKEST
-            IF (nn_hls.EQ.2) THEN
-               CALL lbc_lnk( 'traadv', zuu(:,:,:), 'U', -1., zvv(:,:,:), 'V', -1.)
-               CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)
-            END IF
+            ! [comm_cleanup] - lbc_lnk shifted into step
+            ! IF (nn_hls.EQ.2) THEN
+            !   CALL lbc_lnk( 'traadv', zuu(:,:,:), 'U', -1., zvv(:,:,:), 'V', -1.)
+            !   CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)
+            ! END IF
             CALL tra_adv_qck    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs )
          !
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traatf.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traatf.F90	(revision 14575)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traatf.F90	(revision 14576)
@@ -110,5 +110,6 @@
 #endif
       !                                              ! local domain boundaries  (T-point, unchanged sign)
-      CALL lbc_lnk( 'traatf', pts(:,:,:,jp_tem,Kaa), 'T', 1.0_wp, pts(:,:,:,jp_sal,Kaa), 'T', 1.0_wp )
+      ! [comm_cleanup] ! lbc_lnk moved into stp
+      ! CALL lbc_lnk( 'traatf', pts(:,:,:,jp_tem,Kaa), 'T', 1.0_wp, pts(:,:,:,jp_sal,Kaa), 'T', 1.0_wp )
       !
       IF( ln_bdy )   CALL bdy_tra( kt, Kbb, pts, Kaa )  ! BDY open boundaries
@@ -156,5 +157,6 @@
          ENDIF
          !
-         CALL lbc_lnk( 'traatf',  pts(:,:,:,jp_tem,Kmm) , 'T', 1.0_wp, pts(:,:,:,jp_sal,Kmm) , 'T', 1.0_wp )
+         ! [comm_cleanup] 
+         ! CALL lbc_lnk( 'traatf',  pts(:,:,:,jp_tem,Kmm) , 'T', 1.0_wp, pts(:,:,:,jp_sal,Kmm) , 'T', 1.0_wp )
 
       ENDIF
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf.F90	(revision 14575)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf.F90	(revision 14576)
@@ -92,5 +92,6 @@
             CALL tra_ldf_triad( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts,  1 )
          CASE ( np_blp , np_blp_i , np_blp_it )             ! bilaplacian: iso-level & iso-neutral operators
-            IF(nn_hls.EQ.2) CALL lbc_lnk( 'tra_ldf', pts(:,:,:,:,Kbb), 'T',1.)
+            ! [comm_cleanup]
+            ! IF (nn_hls.EQ.2) CALL lbc_lnk( 'tra_ldf', pts(:,:,:,:,Kbb), 'T',1.)
             CALL tra_ldf_blp  ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs),             jpts, nldf_tra )
          END SELECT
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_iso.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_iso.F90	(revision 14575)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_iso.F90	(revision 14576)
@@ -147,5 +147,6 @@
          ENDIF
          !
-         DO_3D( 0, 0, 0, 0, 1, jpk )
+         ! [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 ) 
             akz     (ji,jj,jk) = 0._wp
             ah_wslp2(ji,jj,jk) = 0._wp
@@ -172,5 +173,6 @@
       IF( kpass == 1 ) THEN                  !==  first pass only  ==!
          !
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         ! [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 ) 
             !
             zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
@@ -189,5 +191,6 @@
          !
          IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient
-            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+            ! [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 ) 
                akz(ji,jj,jk) = 0.25_wp * (                                                                     &
                   &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   &
@@ -198,5 +201,6 @@
             !
             IF( ln_traldf_blp ) THEN                ! bilaplacian operator
-               DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+               ! [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 ) 
                   akz(ji,jj,jk) = 16._wp   &
                      &   * ah_wslp2   (ji,jj,jk)   &
@@ -206,5 +210,6 @@
                END_3D
             ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
-               DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+               ! [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 ) 
                   ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
                   zcoef0 = rDt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
@@ -214,5 +219,6 @@
            !
          ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
-            DO_3D( 0, 0, 0, 0, 1, jpk )
+            ! [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 ) 
                akz(ji,jj,jk) = ah_wslp2(ji,jj,jk)
             END_3D
@@ -233,15 +239,18 @@
 
          ! Horizontal tracer gradient
-         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+         ! [comm_cleanup] ! DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+         DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 
             zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
             zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
          END_3D
          IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient
-            DO_2D( 1, 0, 1, 0 )           ! bottom correction (partial bottom cell)
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )           ! bottom correction (partial bottom cell)
+               DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )            ! bottom correction (partial bottom cell)
                zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)
                zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
             END_2D
             IF( ln_isfcav ) THEN      ! first wet level beneath a cavity
-               DO_2D( 1, 0, 1, 0 )
+               ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+               DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
                   IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)
                   IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)
@@ -256,5 +265,6 @@
          DO jk = 1, jpkm1                                 ! Horizontal slab
             !
-            DO_2D( 1, 1, 1, 1 )
+            ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
+            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
                !                             !== Vertical tracer gradient
                zdk1t(ji,jj) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1)     ! level jk+1
@@ -265,5 +275,6 @@
             END_2D
             !
-            DO_2D( 1, 0, 1, 0 )           !==  Horizontal fluxes
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )           !==  Horizontal fluxes
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )           !==  Horizontal fluxes
                zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm)
                zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
@@ -286,5 +297,6 @@
             END_2D
             !
-            DO_2D( 0, 0, 0, 0 )           !== horizontal divergence and add to pta
+            ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )           !== horizontal divergence and add to pta
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )           !== horizontal divergence and add to pta
                pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    &
                   &       + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   &
@@ -302,5 +314,6 @@
          ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp
 
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )    ! interior (2=<jk=<jpk-1)
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )    ! interior (2=<jk=<jpk-1)
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )    ! interior (2=<jk=<jpk-1)
             !
             zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
@@ -324,5 +337,6 @@
          !                                !==  add the vertical 33 flux  ==!
          IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
-            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+            ! [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 ) 
                ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)   &
                   &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )               &
@@ -333,5 +347,6 @@
             SELECT CASE( kpass )
             CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
-               DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+               ! [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 ) 
                   ztfw(ji,jj,jk) =   &
                      &  ztfw(ji,jj,jk) + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   &
@@ -339,5 +354,6 @@
                END_3D
             CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp.
-               DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+               ! [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 ) 
                   ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)                  &
                      &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   &
@@ -347,5 +363,6 @@
          ENDIF
          !
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )    !==  Divergence of vertical fluxes added to pta  ==!
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )    !==  Divergence of vertical fluxes added to pta  ==!
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )    !==  Divergence of vertical fluxes added to pta  ==! 
             pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * r1_e1e2t(ji,jj)   &
                &                                             / e3t(ji,jj,jk,Kmm)
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_lap_blp.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_lap_blp.F90	(revision 14575)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_lap_blp.F90	(revision 14576)
@@ -229,5 +229,6 @@
       END SELECT
       !
-      CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1.0_wp )     ! Lateral boundary conditions (unchanged sign)
+      ! [comm_cleanup] 
+      IF (nn_hls.EQ.1) CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1.0_wp )     ! Lateral boundary conditions (unchanged sign)
       !                                               ! Partial top/bottom cell: GRADh( zlap )
       IF( ln_isfcav .AND. ln_zps ) THEN   ;   CALL zps_hde_isf( kt, Kmm, kjpt, zlap, zglu, zglv, zgui, zgvi )  ! both top & bottom
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_triad.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_triad.F90	(revision 14575)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_triad.F90	(revision 14576)
@@ -152,5 +152,6 @@
       IF( kpass == 1 ) THEN         !==  first pass only  and whatever the tracer is  ==!
          !
-         DO_3D( 0, 0, 0, 0, 1, jpk )
+         ! [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 ) 
             akz     (ji,jj,jk) = 0._wp
             ah_wslp2(ji,jj,jk) = 0._wp
@@ -159,5 +160,6 @@
          DO ip = 0, 1                            ! i-k triads
             DO kp = 0, 1
-               DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+               ! [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 ) 
                   ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm)
                   zbu   = e1e2u(ji-ip,jj) * e3u(ji-ip,jj,jk,Kmm)
@@ -177,5 +179,6 @@
          DO jp = 0, 1                            ! j-k triads
             DO kp = 0, 1
-               DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+               ! [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 ) 
                   ze3wr = 1.0_wp / e3w(ji,jj,jk+kp,Kmm)
                   zbv   = e1e2v(ji,jj-jp) * e3v(ji,jj-jp,jk,Kmm)
@@ -197,5 +200,6 @@
             !
             IF( ln_traldf_blp ) THEN                ! bilaplacian operator
-               DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+               ! [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 ) 
                   akz(ji,jj,jk) = 16._wp           &
                      &   * ah_wslp2   (ji,jj,jk)   &
@@ -205,5 +209,6 @@
                END_3D
             ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
-               DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+               ! [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 ) 
                   ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
                   zcoef0 = rDt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
@@ -213,5 +218,6 @@
            !
          ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
-            DO_3D( 0, 0, 0, 0, 1, jpk )
+            ! [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 ) 
                akz(ji,jj,jk) = ah_wslp2(ji,jj,jk)
             END_3D
@@ -228,5 +234,6 @@
                DO jp = 0, 1
                   DO kp = 0, 1
-                     DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+                     ! [comm_cleanup] ! DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+                     DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 
                         zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) &
                            & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+jp,jj,jk,1-jp,kp)
@@ -253,15 +260,18 @@
          zftv(:,:,:) = 0._wp
          !
-         DO_3D( 1, 0, 1, 0, 1, jpkm1 )    !==  before lateral T & S gradients at T-level jk  ==!
+         ! [comm_cleanup] ! DO_3D( 1, 0, 1, 0, 1, jpkm1 )    !==  before lateral T & S gradients at T-level jk  ==!
+         DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 
             zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
             zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
          END_3D
          IF( ln_zps .AND. l_grad_zps ) THEN    ! partial steps: correction at top/bottom ocean level
-            DO_2D( 1, 0, 1, 0 )                    ! bottom level
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )                    ! bottom level
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
                zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)
                zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
             END_2D
             IF( ln_isfcav ) THEN                   ! top level (ocean cavities only)
-               DO_2D( 1, 0, 1, 0 )
+               ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+               DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
                   IF( miku(ji,jj)  > 1 )   zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn)
                   IF( mikv(ji,jj)  > 1 )   zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn)
@@ -276,5 +286,6 @@
          DO jk = 1, jpkm1
             !                    !==  Vertical tracer gradient at level jk and jk+1
-            DO_2D( 1, 1, 1, 1 )
+            ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
+            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
                zdkt3d(ji,jj,1) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1)
             END_2D
@@ -283,5 +294,6 @@
             IF( jk == 1 ) THEN   ;   zdkt3d(:,:,0) = zdkt3d(:,:,1)
             ELSE
-               DO_2D( 1, 1, 1, 1 )
+               ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
+               DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
                   zdkt3d(ji,jj,0) = ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
                END_2D
@@ -293,5 +305,6 @@
                DO ip = 0, 1              !==  Horizontal & vertical fluxes
                   DO kp = 0, 1
-                     DO_2D( 1, 0, 1, 0 )
+                     ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+                     DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
                         ze1ur = r1_e1u(ji,jj)
                         zdxt  = zdit(ji,jj,jk) * ze1ur
@@ -314,5 +327,6 @@
                DO jp = 0, 1
                   DO kp = 0, 1
-                     DO_2D( 1, 0, 1, 0 )
+                     ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+                     DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
                         ze2vr = r1_e2v(ji,jj)
                         zdyt  = zdjt(ji,jj,jk) * ze2vr
@@ -336,5 +350,6 @@
                DO ip = 0, 1               !==  Horizontal & vertical fluxes
                   DO kp = 0, 1
-                     DO_2D( 1, 0, 1, 0 )
+                     ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+                     DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
                         ze1ur = r1_e1u(ji,jj)
                         zdxt  = zdit(ji,jj,jk) * ze1ur
@@ -357,5 +372,6 @@
                DO jp = 0, 1
                   DO kp = 0, 1
-                     DO_2D( 1, 0, 1, 0 )
+                     ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+                     DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
                         ze2vr = r1_e2v(ji,jj)
                         zdyt  = zdjt(ji,jj,jk) * ze2vr
@@ -376,5 +392,6 @@
             ENDIF
             !                             !==  horizontal divergence and add to the general trend  ==!
-            DO_2D( 0, 0, 0, 0 )
+            ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
                pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    &
                   &                       + zsign * (  zftu(ji-1,jj  ,jk) - zftu(ji,jj,jk)       &
@@ -387,5 +404,6 @@
          !                                !==  add the vertical 33 flux  ==!
          IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
-            DO_3D( 0, 0, 1, 0, 2, jpkm1 )
+            ! [comm_cleanup] ! DO_3D( 0, 0, 1, 0, 2, jpkm1 )
+            DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
                ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)   &
                   &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
@@ -395,10 +413,12 @@
             SELECT CASE( kpass )
             CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
-               DO_3D( 0, 0, 1, 0, 2, jpkm1 )
+               ! [comm_cleanup] ! DO_3D( 0, 0, 1, 0, 2, jpkm1 )
+               DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
                   ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)             &
                      &                            * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
                END_3D
             CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp.
-               DO_3D( 0, 0, 1, 0, 2, jpkm1 )
+               ! [comm_cleanup] ! DO_3D( 0, 0, 1, 0, 2, jpkm1 )
+               DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
                   ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)                      &
                      &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   &
@@ -408,5 +428,6 @@
          ENDIF
          !
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )      !==  Divergence of vertical fluxes added to pta  ==!
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )      !==  Divergence of vertical fluxes added to pta  ==!
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 
             pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    &
             &                                  + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/tramle.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/tramle.F90	(revision 14575)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/tramle.F90	(revision 14576)
@@ -110,15 +110,18 @@
          SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
          CASE ( 0 )                                               != min of the 2 neighbour MLDs
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zhu(ji,jj) = MIN( hmle(ji+1,jj), hmle(ji,jj) )
                zhv(ji,jj) = MIN( hmle(ji,jj+1), hmle(ji,jj) )
             END_2D
          CASE ( 1 )                                               != average of the 2 neighbour MLDs
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) )
                zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) )
             END_2D
          CASE ( 2 )                                               != max of the 2 neighbour MLDs
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) )
                zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) )
@@ -126,5 +129,6 @@
          END SELECT
          IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2u(ji,jj)                                            &
                     &           * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) )   &
@@ -137,5 +141,6 @@
             !
          ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2u(ji,jj)               &
                     &                  * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) )
@@ -149,9 +154,11 @@
          !                                      !==  MLD used for MLE  ==!
          !                                                ! compute from the 10m density to deal with the diurnal cycle
-         DO_2D( 1, 1, 1, 1 )
+         ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             inml_mle(ji,jj) = mbkt(ji,jj) + 1                    ! init. to number of ocean w-level (T-level + 1)
          END_2D
          IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m
-           DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )        ! from the bottom to nlb10 (10m)
+           ! [comm_cleanup] ! DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )        ! from the bottom to nlb10 (10m)
+           DO_3DS( nn_hls, nn_hls, nn_hls, nn_hls, jpkm1, nlb10, -1 )        ! from the bottom to nlb10 (10m)
               IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer
            END_3D
@@ -163,5 +170,6 @@
          zbm (:,:) = 0._wp
          zn2 (:,:) = 0._wp
-         DO_3D( 1, 1, 1, 1, 1, ikmax )                    ! MLD and mean buoyancy and N2 over the mixed layer
+         ! [comm_cleanup] ! DO_3D( 1, 1, 1, 1, 1, ikmax )                    ! MLD and mean buoyancy and N2 over the mixed layer
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, ikmax )                    ! MLD and mean buoyancy and N2 over the mixed layer
             zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points
             zmld(ji,jj) = zmld(ji,jj) + zc
@@ -172,15 +180,18 @@
          SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
          CASE ( 0 )                                               != min of the 2 neighbour MLDs
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) )
                zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) )
             END_2D
          CASE ( 1 )                                               != average of the 2 neighbour MLDs
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp
                zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp
             END_2D
          CASE ( 2 )                                               != max of the 2 neighbour MLDs
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) )
                zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) )
@@ -188,5 +199,6 @@
          END SELECT
          !                                                ! convert density into buoyancy
-         DO_2D( 1, 1, 1, 1 )
+         ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             zbm(ji,jj) = + grav * zbm(ji,jj) / MAX( e3t(ji,jj,1,Kmm), zmld(ji,jj) )
          END_2D
@@ -201,5 +213,6 @@
          !
          IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2_e1u(ji,jj)                                            &
                     &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )   &
@@ -212,5 +225,6 @@
             !
          ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2_e1u(ji,jj)               &
                     &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )
@@ -222,5 +236,6 @@
          !
          IF( nn_conv == 1 ) THEN              ! No MLE in case of convection
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp
                IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp
@@ -230,5 +245,6 @@
       ENDIF  ! end of ln_osm_mle conditional
     !                                      !==  structure function value at uw- and vw-points  ==!
-    DO_2D( 1, 0, 1, 0 )
+    ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+    DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
        zhu(ji,jj) = 1._wp / MAX(zhu(ji,jj), rsmall)                   ! hu --> 1/hu
        zhv(ji,jj) = 1._wp / MAX(zhv(ji,jj), rsmall) 
@@ -238,5 +254,7 @@
     zpsi_vw(:,:,:) = 0._wp
     !
-      DO_3D( 1, 0, 1, 0, 2, ikmax )                ! start from 2 : surface value = 0
+      ! [comm_cleanup] ! DO_3D( 1, 0, 1, 0, 2, ikmax )                ! start from 2 : surface value = 0
+      DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 2, ikmax )                ! start from 2 : surface value = 0
+      
          zcuw = 1._wp - ( gdepw(ji+1,jj,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhu(ji,jj)
          zcvw = 1._wp - ( gdepw(ji,jj+1,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhv(ji,jj)
@@ -252,9 +270,11 @@
       !                                      !==  transport increased by the MLE induced transport ==!
       DO jk = 1, ikmax
-         DO_2D( 1, 0, 1, 0 )                      ! CAUTION pu,pv must be defined at row/column i=1 / j=1
+         ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )                      ! CAUTION pu,pv must be defined at row/column i=1 / j=1
+         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
             pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
             pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
          END_2D
-         DO_2D( 0, 0, 0, 0 )
+         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk)   &
                &                          + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) ) * wmask(ji,jj,1)
@@ -270,9 +290,11 @@
          !
          IF (ln_osm_mle.and.ln_zdfosm) THEN
-            DO_2D( 0, 0, 0, 0 )
+            ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                zLf_NH(ji,jj) = SQRT( rb_c * hmle(ji,jj) ) * r1_ft(ji,jj)      ! Lf = N H / f
             END_2D
          ELSE
-            DO_2D( 0, 0, 0, 0 )
+            ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                zLf_NH(ji,jj) = SQRT( rb_c * zmld(ji,jj) ) * r1_ft(ji,jj)      ! Lf = N H / f
             END_2D
@@ -280,5 +302,6 @@
          !
          ! divide by cross distance to give streamfunction with dimensions m^2/s
-         DO_3D( 0, 0, 0, 0, 1, ikmax+1 )
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, ikmax+1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, ikmax+1 )
             zpsiu_mle(ji,jj,jk) = zpsi_uw(ji,jj,jk) * r1_e2u(ji,jj)
             zpsiv_mle(ji,jj,jk) = zpsi_vw(ji,jj,jk) * r1_e1v(ji,jj)
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/step.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/step.F90	(revision 14575)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/step.F90	(revision 14576)
@@ -289,4 +289,20 @@
 #endif
 
+     ! [comm_cleanup]
+     IF (nn_hls.EQ.2) THEN
+         SELECT CASE ( nadv )
+         CASE ( np_FCT )                                 ! FCT scheme      : 2nd / 4th order
+               CALL lbc_lnk( 'stp', ts(:,:,:,:,Nbb), 'T', 1., ts(:,:,:,:,Nnn), 'T', 1.)
+               CALL lbc_lnk( 'stp', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1., ww(:,:,:), 'W', 1.)
+         CASE ( np_MUS )                                 ! MUSCL
+                CALL lbc_lnk( 'stp', ts(:,:,:,:,Nbb), 'T', 1.)
+         CASE ( np_UBS )                                 ! UBS
+                CALL lbc_lnk( 'stp', ts(:,:,:,:,Nbb), 'T', 1.)
+         CASE ( np_QCK )                                 ! QUICKEST
+               CALL lbc_lnk( 'stp', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.)
+               CALL lbc_lnk( 'stp', ts(:,:,:,:,Nbb), 'T', 1.)
+         END SELECT
+      ENDIF
+
       ! TEMP: [tiling] Separate loop over tile domains (due to tra_adv workarounds for tiling)
       DO jtile = 1, nijtile
@@ -323,4 +339,6 @@
 !!
 !!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine
+                         ! [comm_cleanup] 
+                         CALL lbc_lnk( 'stp', ts(:,:,:,jp_tem,Naa), 'T', 1.0_wp, ts(:,:,:,jp_sal,Naa), 'T', 1.0_wp )
                          CALL tra_atf       ( kstp, Nbb, Nnn, Naa, ts )                      ! time filtering of "now" tracer arrays
                          CALL dyn_atf       ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v  )  ! time filtering of "now" velocities and scale factors
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/stpmlf.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/stpmlf.F90	(revision 14575)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/stpmlf.F90	(revision 14576)
@@ -309,4 +309,20 @@
 #endif
 
+     ! [comm_cleanup]
+     IF (nn_hls.EQ.2) THEN
+         SELECT CASE ( nadv )
+         CASE ( np_FCT )                                 ! FCT scheme : 2nd / 4th order
+               CALL lbc_lnk( 'stp_MLF', ts(:,:,:,:,Nbb), 'T', 1., ts(:,:,:,:,Nnn), 'T', 1.)
+               CALL lbc_lnk( 'stp_MLF', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1., ww(:,:,:), 'W', 1.)
+         CASE ( np_MUS )                                 ! MUSCL
+                CALL lbc_lnk( 'stp_MLF', ts(:,:,:,:,Nbb), 'T', 1.)
+         CASE ( np_UBS )                                 ! UBS
+                CALL lbc_lnk( 'stp_MLF', ts(:,:,:,:,Nbb), 'T', 1.)
+         CASE ( np_QCK )                                 ! QUICKEST
+               CALL lbc_lnk( 'stp_MLF', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.)
+               CALL lbc_lnk( 'stp_MLF', ts(:,:,:,:,Nbb), 'T', 1.)
+         END SELECT
+      ENDIF
+
       ! TEMP: [tiling] Separate loop over tile domains (due to tra_adv workarounds for tiling)
       DO jtile = 1, nijtile
