Index: NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfddm.F90
===================================================================
--- NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfddm.F90	(revision 14598)
+++ NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfddm.F90	(revision 14601)
@@ -95,5 +95,6 @@
 !!gm                            and many acces in memory
          
-         DO_2D( 1, 1, 1, 1 )           !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==!
+         ! [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])  ==!
             zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
 !!gm please, use e3w at Kmm below 
@@ -111,5 +112,6 @@
          END_2D
 
-         DO_2D( 1, 1, 1, 1 )           !==  indicators  ==!
+         ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )           !==  indicators  ==!
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )           !==  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
@@ -141,5 +143,6 @@
          ! ------------------
          ! Constant eddy coefficient: reset to the background value
-         DO_2D( 1, 1, 1, 1 )
+         ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             zinr = 1._wp / zrau(ji,jj)
             ! salt fingering
Index: NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfdrg.F90
===================================================================
--- NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfdrg.F90	(revision 14598)
+++ NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfdrg.F90	(revision 14601)
@@ -117,5 +117,6 @@
       !
       IF( l_log_not_linssh ) THEN     !==  "log layer"  ==!   compute Cd and -Cd*|U|
-         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 )
             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
@@ -129,5 +130,6 @@
          END_2D
       ELSE                                            !==  standard Cd  ==!
-         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 )
             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
@@ -176,5 +178,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 )
          ikbu = mbku(ji,jj)          ! deepest wet ocean u- & v-levels
          ikbv = mbkv(ji,jj)
@@ -189,5 +192,6 @@
       !
       IF( ln_isfcav ) THEN        ! ocean cavities
-         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 )
             ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels
             ikbv = mikv(ji,jj)
@@ -432,5 +436,6 @@
             l_log_not_linssh = .FALSE.    !- don't update Cd at each time step
             !
-            DO_2D( 1, 1, 1, 1 )              ! pCd0 = mask (and boosted) logarithmic drag coef.
+            ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )              ! pCd0 = mask (and boosted) logarithmic drag coef.
+            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )              ! pCd0 = mask (and boosted) logarithmic drag coef.
                zzz =  0.5_wp * e3t_0(ji,jj,k_mk(ji,jj))
                zcd = (  vkarmn / LOG( zzz / rn_z0 )  )**2
Index: NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfevd.F90
===================================================================
--- NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfevd.F90	(revision 14598)
+++ NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfevd.F90	(revision 14601)
@@ -87,5 +87,6 @@
 !         END WHERE
          !
-         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 )
             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)
@@ -103,5 +104,6 @@
 !         END WHERE
 
-         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 )
             IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 )   &
                p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)
Index: NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfgls.F90
===================================================================
--- NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfgls.F90	(revision 14598)
+++ NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfgls.F90	(revision 14601)
@@ -179,5 +179,6 @@
 
       ! Compute surface, top and bottom friction at T-points
-      DO_2D( 0, 0, 0, 0 )          !==  surface ocean friction
+      ! [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 
          ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1)   ! surface friction
       END_2D
@@ -186,5 +187,6 @@
       !
       IF( .NOT.ln_drg_OFF ) THEN     !== top/bottom friction   (explicit before friction)
-         DO_2D( 0, 0, 0, 0 )         ! 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)
             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)
@@ -193,5 +195,6 @@
          END_2D
          IF( ln_isfcav ) THEN
-            DO_2D( 0, 0, 0, 0 )      ! top friction
+            ! [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
                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)
@@ -220,5 +223,6 @@
       zhsro(:,:) = ( (1._wp-zice_fra(:,:)) * zhsro(:,:) + zice_fra(:,:) * rn_hsri )*tmask(:,:,1)  + (1._wp - tmask(:,:,1))*rn_hsro
       !
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )  !==  Compute dissipation rate  ==!
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )  !==  Compute dissipation rate  ==!
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )  !==  Compute dissipation rate  ==!
          eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk)
       END_3D
@@ -229,5 +233,6 @@
 
       IF( nn_clos == 0 ) THEN    ! Mellor-Yamada
-         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 )
             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) )
@@ -250,5 +255,6 @@
       ! Warning : after this step, en : right hand side of the matrix
 
-      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 )
          !
          buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk)     ! stratif. destruction
@@ -327,5 +333,6 @@
       ! at k=2, set de/dz=Fw
       !cbr
-      DO_2D( 0, 0, 0, 0 )   ! zdiag zd_lw not defined/used on the halo
+      ! [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
@@ -348,5 +355,6 @@
          !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin
          !                      ! Balance between the production and the dissipation terms
-         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 )
 !!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 ???
@@ -366,5 +374,6 @@
          !
          IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
-            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 )
                itop   = mikt(ji,jj)       ! k   top w-point
                itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
@@ -384,5 +393,6 @@
       CASE ( 1 )             ! Neumman boundary condition
          !
-         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 )
             ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
             ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -399,5 +409,6 @@
          END_2D
          IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
-            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 )
                itop   = mikt(ji,jj)       ! k   top w-point
                itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
@@ -420,11 +431,14 @@
       ! ----------------------------------------------------------
       !
-      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)
       END_3D
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
          zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1)
       END_3D
-      DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )           ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+      ! [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
          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
@@ -441,20 +455,24 @@
       !
       CASE( 0 )               ! k-kl  (Mellor-Yamada)
-         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 )
             psi(ji,jj,jk)  = eb(ji,jj,jk) * hmxl_b(ji,jj,jk)
          END_3D
          !
       CASE( 1 )               ! k-eps
-         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 )
             psi(ji,jj,jk)  = eps(ji,jj,jk)
          END_3D
          !
       CASE( 2 )               ! k-w
-         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 )
             psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) )
          END_3D
          !
       CASE( 3 )               ! generic
-         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 )
             psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn
          END_3D
@@ -469,5 +487,6 @@
       ! Warning : after this step, en : right hand side of the matrix
 
-      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 )
          !
          ! psi / k
@@ -541,5 +560,6 @@
          !
          ! Neumann condition at k=2
-         DO_2D( 0, 0, 0, 0 )   ! zdiag zd_lw not defined/used on the halo
+         ! [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
@@ -569,5 +589,6 @@
          !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot
          !                      ! Balance between the production and the dissipation terms
-         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 )
             ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
             ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -588,5 +609,6 @@
       CASE ( 1 )             ! Neumman boundary condition
          !
-         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 )
             ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
             ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -616,11 +638,14 @@
       ! ----------------
       !
-      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)
       END_3D
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
          zd_lw(ji,jj,jk) = psi(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1)
       END_3D
-      DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+      ! [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
          psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
       END_3D
@@ -632,15 +657,18 @@
       !
       CASE( 0 )               ! k-kl  (Mellor-Yamada)
-         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 )
             eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / MAX( psi(ji,jj,jk), rn_epsmin)
          END_3D
          !
       CASE( 1 )               ! k-eps
-         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 )
             eps(ji,jj,jk) = psi(ji,jj,jk)
          END_3D
          !
       CASE( 2 )               ! k-w
-         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 )
             eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk)
          END_3D
@@ -650,5 +678,6 @@
          zex1  =      ( 1.5_wp + rmm/rnn )
          zex2  = -1._wp / rnn
-         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 )
             eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2
          END_3D
@@ -658,5 +687,6 @@
       ! Limit dissipation rate under stable stratification
       ! --------------------------------------------------
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! Note that this set boundary conditions on hmxl_n at the same time
+      ! [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
          ! limitation
          eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
@@ -674,5 +704,6 @@
       !
       CASE ( 0 , 1 )             ! Galperin or Kantha-Clayson stability functions
-         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 )
             ! zcof =  l²/q²
             zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
@@ -691,5 +722,6 @@
          !
       CASE ( 2, 3 )               ! Canuto stability functions
-         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 )
             ! zcof =  l²/q²
             zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )
@@ -723,5 +755,6 @@
       ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0)
       zstm(:,:,jpk) = 0.
-      DO_2D( 0, 0, 0, 0 )             ! update bottom with good values
+      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )             ! update bottom with good values
+      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )             ! update bottom with good values
          zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj))
       END_2D
@@ -738,5 +771,6 @@
       !     later overwritten by surface/bottom boundaries conditions, so we don't really care of p_avm(:,:1) and p_avm(:,:jpk)
       !     for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0)
-      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 )
          zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk)
          zavt  = zsqen * zstt(ji,jj,jk)
Index: NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfiwm.F90
===================================================================
--- NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfiwm.F90	(revision 14598)
+++ NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfiwm.F90	(revision 14601)
@@ -143,15 +143,18 @@
       ! Set to zero the 1st and last vertical levels of appropriate variables
       IF( iom_use("emix_iwm") ) 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 )
             zemx_iwm (ji,jj,1) = 0._wp   ;   zemx_iwm (ji,jj,jpk) = 0._wp
          END_2D
       ENDIF
       IF( iom_use("av_ratio") ) 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 )
             zav_ratio(ji,jj,1) = 0._wp   ;   zav_ratio(ji,jj,jpk) = 0._wp
          END_2D
       ENDIF
       IF( iom_use("av_wave") .OR. sn_cfctl%l_prtctl ) 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 )
             zav_wave (ji,jj,1) = 0._wp   ;   zav_wave (ji,jj,jpk) = 0._wp
          END_2D
@@ -164,5 +167,6 @@
       !                       !* Critical slope mixing: distribute energy over the time-varying ocean depth,
       !                                                 using an exponential decay from the seafloor.
-      DO_2D( 0, 0, 0, 0 )             ! part independent of the level
+      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )             ! part independent of the level
+      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )             ! part independent of the level
          zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
          zfact(ji,jj) = rho0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  )
@@ -170,5 +174,6 @@
       END_2D
 !!gm gde3w ==>>>  check for ssh taken into account.... seem OK gde3w_n=gdept(:,:,:,Kmm) - ssh(:,:,Kmm)
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )   ! complete with the level-dependent part
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )   ! complete with the level-dependent part
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )   ! complete with the level-dependent part
          IF ( zfact(ji,jj) == 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization
             zemx_iwm(ji,jj,jk) = 0._wp
@@ -190,16 +195,20 @@
       CASE ( 1 )               ! Dissipation scales as N (recommended)
          !
-         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 )
             zfact(ji,jj) = 0._wp
          END_2D
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! part independent of the level
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! part independent of the level
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! part independent of the level
             zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(ji,jj,jk) )  ) * wmask(ji,jj,jk)
          END_3D
          !
-         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 )
             IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) )
          END_2D
          !
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! complete with the level-dependent part
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! complete with the level-dependent part
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! complete with the level-dependent part
             zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * SQRT(  MAX( 0._wp, rn2(ji,jj,jk) )  ) * wmask(ji,jj,jk)
          END_3D
@@ -207,16 +216,20 @@
       CASE ( 2 )               ! Dissipation scales as N^2
          !
-         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 )
             zfact(ji,jj) = 0._wp
          END_2D
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! part independent of the level
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! part independent of the level
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! part independent of the level
             zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * MAX( 0._wp, rn2(ji,jj,jk) ) * wmask(ji,jj,jk)
          END_3D
          !
-         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 )
             IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) )
          END_2D
          !
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! complete with the level-dependent part
+         ! [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 )
             zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * MAX( 0._wp, rn2(ji,jj,jk) ) * wmask(ji,jj,jk)
          END_3D
@@ -227,23 +240,29 @@
       !                        !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot)
       !
-      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 )
          zwkb(ji,jj,1) = 0._wp
       END_2D
-      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 )
          zwkb(ji,jj,jk) = zwkb(ji,jj,jk-1) + e3w(ji,jj,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(ji,jj,jk) )  ) * wmask(ji,jj,jk)
       END_3D
-      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 )
          zfact(ji,jj) = zwkb(ji,jj,jpkm1)
       END_2D
       !
-      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 )
          IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   &
             &                                     * wmask(ji,jj,jk) / zfact(ji,jj)
       END_3D
-      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 )
          zwkb (ji,jj,1) = zhdep(ji,jj) * wmask(ji,jj,1)
       END_2D
       !
-      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 )
          IF ( rn2(ji,jj,jk) <= 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization: EXP coast a lot
             zweight(ji,jj,jk) = 0._wp
@@ -254,16 +273,20 @@
       END_3D
       !
-      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 )
          zfact(ji,jj) = 0._wp
       END_2D
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! part independent of the level
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! part independent of the level
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! part independent of the level
          zfact(ji,jj) = zfact(ji,jj) + zweight(ji,jj,jk)
       END_3D
       !
-      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 )
          IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rho0 * zfact(ji,jj) )
       END_2D
       !
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! complete with the level-dependent part
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! complete with the level-dependent part
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! complete with the level-dependent part
          zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zweight(ji,jj,jk) * zfact(ji,jj) * wmask(ji,jj,jk)   &
             &                                                        / ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) )
@@ -273,10 +296,12 @@
 !!gm  this is to be replaced by just a constant value znu=1.e-6 m2/s
       ! Calculate molecular kinematic viscosity
-      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 )
          znu_t(ji,jj,jk) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * ts(ji,jj,jk,jp_tem,Kmm)   &
             &                                     + 0.00694_wp * ts(ji,jj,jk,jp_tem,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)  &
             &                                     + 0.02305_wp * ts(ji,jj,jk,jp_sal,Kmm)  ) * tmask(ji,jj,jk) * r1_rho0
       END_3D
-      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 )
          znu_w(ji,jj,jk) = 0.5_wp * ( znu_t(ji,jj,jk-1) + znu_t(ji,jj,jk) ) * wmask(ji,jj,jk)
       END_3D
@@ -284,15 +309,18 @@
       !
       ! Calculate turbulence intensity parameter Reb
-      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 )
          zReb(ji,jj,jk) = zemx_iwm(ji,jj,jk) / MAX( 1.e-20_wp, znu_w(ji,jj,jk) * rn2(ji,jj,jk) )
       END_3D
       !
       ! Define internal wave-induced diffusivity
-      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 )
          zav_wave(ji,jj,jk) = znu_w(ji,jj,jk) * zReb(ji,jj,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6
       END_3D
       !
       IF( ln_mevar ) THEN                ! Variable mixing efficiency case : modify zav_wave in the
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )   ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )   ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )   ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes
             IF( zReb(ji,jj,jk) > 480.00_wp ) THEN
                zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
@@ -303,5 +331,6 @@
       ENDIF
       !
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )      ! Bound diffusivity by molecular value and 100 cm2/s
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )      ! Bound diffusivity by molecular value and 100 cm2/s
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )      ! Bound diffusivity by molecular value and 100 cm2/s
          zav_wave(ji,jj,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(ji,jj,jk) ), 1.e-2_wp  ) * wmask(ji,jj,jk)
       END_3D
@@ -310,5 +339,6 @@
          zztmp = 0._wp
 !!gm used of glosum 3D....
-         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 )
             zztmp = zztmp + e3w(ji,jj,jk,Kmm) * e1e2t(ji,jj)   &
                &          * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
@@ -332,5 +362,6 @@
       IF( ln_tsdiff ) THEN                !* Option for differential mixing of salinity and temperature
          ztmp1 = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10( 1.e-20_wp ) - 0.60_wp ) )
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! Calculate S/T diffusivity ratio as a function of Reb
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! Calculate S/T diffusivity ratio as a function of Reb
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! Calculate S/T diffusivity ratio as a function of Reb
             ztmp2 = zReb(ji,jj,jk) * 5._wp * r1_6
             IF ( ztmp2 > 1.e-20_wp .AND. wmask(ji,jj,jk) == 1._wp ) THEN
@@ -341,5 +372,6 @@
          END_3D
          CALL iom_put( "av_ratio", zav_ratio )
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )    !* update momentum & tracer diffusivity with wave-driven mixing
+         ! [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
             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)
@@ -348,5 +380,6 @@
          !
       ELSE                                !* update momentum & tracer diffusivity with wave-driven mixing
-         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 )
             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)
Index: NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfmfc.F90
===================================================================
--- NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfmfc.F90	(revision 14598)
+++ NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfmfc.F90	(revision 14601)
@@ -218,5 +218,6 @@
          WHERE(zrautbm1 .NE. 0.) zfbuo(:,:)  =  grav * (zraupl(:,:) - zrautbm1(:,:)) / zrautbm1(:,:)
 
-         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 )
 
             ! Compute Environment of Plume. Interpolation T/S (before time step) on W-points
@@ -376,5 +377,6 @@
       !
       !
-      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.)
+      ! [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.)
       !
    END SUBROUTINE tra_mfc
Index: NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfmxl.F90
===================================================================
--- NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfmxl.F90	(revision 14598)
+++ NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfmxl.F90	(revision 14601)
@@ -99,5 +99,6 @@
       hmlp(:,:)  = 0._wp                  ! here hmlp used as a dummy variable, integrating vertically N^2
       zN2_c = grav * rho_c * r1_rho0      ! convert density criteria into N^2 criteria
-      DO_3D( 1, 1, 1, 1, nlb10, jpkm1 )   ! Mixed layer level: w-level
+      ! [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
          ikt = mbkt(ji,jj)
          hmlp(ji,jj) =   &
@@ -108,9 +109,11 @@
       ! w-level of the turbocline and mixing layer (iom_use)
       imld(:,:) = mbkt(:,:) + 1                ! Initialization to the number of w ocean point
-      DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )   ! from the bottom to nlb10
+      ! [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
-      DO_2D( 1, 1, 1, 1 )
+      ! [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)
Index: NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfosm.F90
===================================================================
--- NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfosm.F90	(revision 14598)
+++ NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfosm.F90	(revision 14601)
@@ -400,5 +400,6 @@
      zz0 =       rn_abs       ! surface equi-partition in 2-bands
      zz1 =  1. - rn_abs
-     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 )
         ! Surface downward irradiance (so always +ve)
         zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp
@@ -410,5 +411,6 @@
      END_2D
      ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL
-     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 )
         zthermal = rab_n(ji,jj,1,jp_tem)
         zbeta    = rab_n(ji,jj,1,jp_sal)
@@ -437,5 +439,6 @@
      ! Assume constant La#=0.3
      CASE(0)
-        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 )
            zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2
            zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2
@@ -446,5 +449,6 @@
      ! Assume Pierson-Moskovitz wind-wave spectrum
      CASE(1)
-        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 )
            ! Use wind speed wndm included in sbc_oce module
            zustke(ji,jj) =  MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )
@@ -455,5 +459,6 @@
         zfac =  2.0_wp * rpi / 16.0_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 )
            IF (hsw(ji,jj) > 1.e-4) THEN
               ! Use  wave fields
@@ -472,5 +477,6 @@
      IF (ln_zdfosm_ice_shelter) THEN
         ! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice
-        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 )
            zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - fr_i(ji,jj))
            dstokes(ji,jj) = dstokes(ji,jj) * (1.0_wp - fr_i(ji,jj))
@@ -494,5 +500,6 @@
         z_two_thirds = 2.0_wp / 3.0_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 )
            zthickness = rn_osm_hblfrac*hbl(ji,jj)
            z2k_times_thickness =  zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp )
@@ -509,5 +516,6 @@
         zsqrtpi = SQRT(rpi)
 
-        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 )
            zthickness = rn_osm_hblfrac*hbl(ji,jj)
            z2k_times_thickness =  zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp )
@@ -530,5 +538,6 @@
      ! Langmuir velocity scale (zwstrl), La # (zla)
      ! mixed scale (zvstr), convective velocity scale (zwstrc)
-     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 )
         ! Langmuir velocity scale (zwstrl), at T-point
         zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird
@@ -563,5 +572,6 @@
       hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,4,Kmm) )
       ibld(:,:) = 4
-      DO_3D( 1, 1, 1, 1, 5, jpkm1 )
+      ! [comm_cleanup] ! DO_3D( 1, 1, 1, 1, 5, jpkm1 )
+      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 5, jpkm1 )
          IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN
             ibld(ji,jj) = MIN(mbkt(ji,jj), jk)
@@ -570,5 +580,6 @@
      ! ##########################################################################
 
-      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 )
          zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm)
          imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t(ji, jj, ibld(ji,jj), Kmm )) , 1 ))
@@ -590,5 +601,6 @@
 ! Fox-Kemper Scheme
          mld_prof = 4
-         DO_3D( 0, 0, 0, 0, 5, jpkm1 )
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 5, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 5, jpkm1 )
          IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
          END_3D
@@ -596,5 +608,6 @@
         CALL zdf_osm_vertical_average(mld_prof, jp_ext_mle, zt_mle, zs_mle, zb_mle, zu_mle, zv_mle, zdt_mle, zds_mle, zdb_mle, zdu_mle, zdv_mle)
 
-         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 )
            zhmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm)
          END_2D
@@ -611,5 +624,6 @@
          lflux(:,:) = .FALSE.
          lmle(:,:) = .FALSE.
-         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 )
           IF ( lconv(ji,jj) .AND. zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE.
          END_2D
@@ -617,5 +631,6 @@
 
 ! Test if pycnocline well resolved
-      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 )
        IF (lconv(ji,jj) ) THEN
           ztmp = 0.2 * zhbl(ji,jj) / e3w(ji,jj,ibld(ji,jj),Kmm)
@@ -638,5 +653,6 @@
 ! Rate of change of hbl
       CALL zdf_osm_calculate_dhdt( zdhdt, zddhdt )
-      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 )
        zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need ww here, so subtract it
             ! adjustment to represent limiting by ocean bottom
@@ -650,5 +666,6 @@
       ibld(:,:) = 4
 
-      DO_3D( 0, 0, 0, 0, 4, jpkm1 )
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 4, jpkm1 )
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 4, jpkm1 )
          IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN
             ibld(ji,jj) = jk
@@ -669,5 +686,6 @@
       CALL zdf_osm_pycnocline_thickness( dh, zdh )
 
-      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 )
        IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) + jp_ext(ji,jj) >= mbkt(ji,jj) .or. ibld(ji,jj)-imld(ji,jj) == 1 ) lpyc(ji,jj) = .FALSE.
       END_2D
@@ -709,5 +727,6 @@
 
 
-       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 )
          IF ( lconv(ji,jj) ) THEN
            DO jk = 2, imld(ji,jj)
@@ -744,5 +763,6 @@
           IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv )
        END IF
-       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 )
           IF ( lconv(ji,jj) ) THEN
              DO jk = 2, imld(ji,jj)
@@ -776,5 +796,6 @@
        ENDWHERE
 
-       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 )
           IF (lconv(ji,jj) ) THEN
              DO jk = 2, imld(ji,jj)
@@ -838,5 +859,6 @@
        ENDWHERE
 
-       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 )
           IF ( lconv(ji,jj) ) THEN
              DO jk = 2 , imld(ji,jj)
@@ -856,5 +878,6 @@
        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 )
         IF ( lpyc(ji,jj) ) THEN
           IF ( j_ddh(ji,jj) == 0 ) THEN
@@ -891,5 +914,6 @@
 ! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ]
 
-       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 ( lconv(ji,jj) ) THEN
@@ -907,5 +931,6 @@
        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 )
          IF ( lconv(ji,jj) ) THEN
             DO jk = 2, imld(ji,jj)
@@ -954,5 +979,6 @@
        ENDWHERE
 
-       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 )
           IF ( lconv(ji,jj) ) THEN
             DO jk = 2, imld(ji,jj)
@@ -1001,5 +1027,6 @@
  ! Make surface forced velocity non-gradient terms go to zero at the base of the boundary layer.
 
-      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 )
          IF ( .not. lconv(ji,jj) ) THEN
             DO jk = 2, ibld(ji,jj)
@@ -1017,5 +1044,6 @@
 
       ! pynocline contributions
-       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 )
          IF ( .not. lconv(ji,jj) ) THEN
           IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
@@ -1035,5 +1063,6 @@
        END IF
 
-       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 )
           ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
           ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
@@ -1060,5 +1089,6 @@
        ! rotate non-gradient velocity terms back to model reference frame
 
-       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 )
           DO jk = 2, ibld(ji,jj)
              ztemp = ghamu(ji,jj,jk)
@@ -1076,5 +1106,7 @@
 ! KPP-style Ri# mixing
        IF( ln_kpprimix) THEN
-          DO_3D( 1, 0, 1, 0, 2, jpkm1 )      !* Shear production at uw- and vw-points (energy conserving form)
+          
+          ! [comm_cleanup] ! DO_3D( 1, 0, 1, 0, 2, jpkm1 )      !* Shear production at uw- and vw-points (energy conserving form)
+          DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 2, jpkm1 )
              z3du(ji,jj,jk) = 0.5 * (  uu(ji,jj,jk-1,Kmm) -  uu(ji  ,jj,jk,Kmm) )   &
                   &                 * (  uu(ji,jj,jk-1,Kbb) -  uu(ji  ,jj,jk,Kbb) ) * wumask(ji,jj,jk) &
@@ -1085,5 +1117,6 @@
           END_3D
       !
-         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 )
             !                                          ! shear prod. at w-point weightened by mask
             zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
@@ -1096,5 +1129,6 @@
          END_3D
 
-          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 )
              DO jk = ibld(ji,jj) + 1, jpkm1
                 zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
@@ -1107,5 +1141,6 @@
 ! KPP-style set diffusivity large if unstable below BL
        IF( ln_convmix) 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 )
              DO jk = ibld(ji,jj) + 1, jpkm1
                IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv
@@ -1117,5 +1152,6 @@
 
        IF ( ln_osm_mle ) THEN  ! set up diffusivity and non-gradient mixing
-          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 )
               IF ( lflux(ji,jj) ) THEN ! MLE mixing extends below boundary layer
              ! Calculate MLE flux contribution from surface fluxes
@@ -1158,11 +1194,12 @@
        ! GN 25/8: need to change tmask --> wmask
 
-     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 )
           p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
           p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk)
      END_3D
       ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid  (sign unchanged), needed to caclulate gham[uv] on u and v grids
-     CALL lbc_lnk( 'zdfosm', p_avt, 'W', 1.0_wp , p_avm, 'W', 1.0_wp,   &
-        &                    ghamu, 'W', 1.0_wp , ghamv, 'W', 1.0_wp )
+     IF (nn_hls.eq.1) CALL lbc_lnk( 'zdfosm', p_avt, 'W', 1.0_wp , p_avm, 'W', 1.0_wp,   &
+                         &                    ghamu, 'W', 1.0_wp , ghamv, 'W', 1.0_wp )
        DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) &
@@ -1176,9 +1213,11 @@
        END_3D
         ! Lateral boundary conditions on final outputs for hbl,  on T-grid (sign unchanged)
-        CALL lbc_lnk( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. )
+        ! [comm_cleanup] ! no need lbc_lnk for output
+        ! 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)
-        CALL lbc_lnk( 'zdfosm', ghamt, 'W',  1.0_wp , ghams, 'W',  1.0_wp,   &
-           &                    ghamu, 'U', -1.0_wp , ghamv, 'V', -1.0_wp )
+        ! [comm_cleanup] ! no need lbc_lnk for output
+        ! CALL lbc_lnk( 'zdfosm', ghamt, 'W',  1.0_wp , ghams, 'W',  1.0_wp,   &
+        !   &                    ghamu, 'U', -1.0_wp , ghamv, 'V', -1.0_wp )
 
       IF(ln_dia_osm) THEN
@@ -1280,5 +1319,6 @@
       REAL(wp), PARAMETER :: rn_vispyc_shr = 0.15
 
-      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 )
           IF ( lconv(ji,jj) ) THEN
 
@@ -1323,5 +1363,6 @@
       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 )
           IF ( lconv(ji,jj) ) THEN
              DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity
@@ -1422,5 +1463,6 @@
 
 ! Determins stability and set flag lconv
-     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 )
       IF ( zhol(ji,jj) < 0._wp ) THEN
          lconv(ji,jj) = .TRUE.
@@ -1439,5 +1481,6 @@
      j_ddh(:,:) = 1
 
-     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 )
       IF ( lconv(ji,jj) ) THEN
          IF ( zdb_bl(ji,jj) > 0._wp ) THEN
@@ -1476,5 +1519,6 @@
 ! Calculate entrainment buoyancy flux due to surface fluxes.
 
-     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 )
       IF ( lconv(ji,jj) ) THEN
         zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln
@@ -1501,5 +1545,6 @@
      zwb_min(:,:) = 0._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 )
       IF ( lshear(ji,jj) ) THEN
         IF ( lconv(ji,jj) ) THEN
@@ -1562,5 +1607,6 @@
         zu   = 0._wp
         zv   = 0._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 )
          zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
          zbeta    = rab_n(ji,jj,1,jp_sal)
@@ -1619,5 +1665,6 @@
         REAL(wp) :: ztemp
 
-        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 )
            ztemp = zu(ji,jj)
            zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj)
@@ -1653,9 +1700,11 @@
       znd_param(:,:) = 0._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 )
           ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
           zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * zdbds_mle(ji,jj) * zdbds_mle(ji,jj)
         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 )
                  !
          IF ( lconv(ji,jj) ) THEN
@@ -1681,5 +1730,6 @@
 
 ! Diagnosis
-        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 )
           IF ( lconv(ji,jj) ) THEN
             zwb_ent = - 2.0 * 0.2 * zwbav(ji,jj) &
@@ -1751,5 +1801,6 @@
 
 
-     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 )
         IF ( jbase(ji,jj)+1 < mbkt(ji,jj) ) THEN
            zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
@@ -1781,5 +1832,6 @@
      REAL(wp), PARAMETER :: zgamma_b = 2.25, zzeta_sh = 0.15
 
-     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 )
         IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
            IF ( lconv(ji,jj) ) THEN  ! convective conditions
@@ -1866,5 +1918,6 @@
       REAL(wp) :: zzeta_v = 0.45
       !
-      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 )
          !
          IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
@@ -1929,5 +1982,6 @@
     REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth
 
-  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 )
 
     IF ( lshear(ji,jj) ) THEN
@@ -2072,5 +2126,6 @@
     REAL(wp) :: zthermal, zbeta
 
-     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 )
         IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN
 !
@@ -2176,5 +2231,6 @@
       REAL, PARAMETER :: a_ddh_2 = 3.5 ! also in pycnocline_depth
 
-    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 )
 
       IF ( lshear(ji,jj) ) THEN
@@ -2322,10 +2378,12 @@
       zmld(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2
       zN2_c = grav * rn_osm_mle_rho_c * r1_rho0   ! convert density criteria into N^2 criteria
-      DO_3D( 1, 1, 1, 1, nlb10, jpkm1 )
+      ! [comm_cleanup] ! DO_3D( 1, 1, 1, 1, nlb10, jpkm1 )
+      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, nlb10, jpkm1 )
          ikt = mbkt(ji,jj)
          zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm)
          IF( zmld(ji,jj) < zN2_c )   mld_prof(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
       END_3D
-      DO_2D( 1, 1, 1, 1 )
+      ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj))
          zmld(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm)
@@ -2337,5 +2395,6 @@
       ztm(:,:) = 0._wp
       zsm(:,:) = 0._wp
-      DO_3D( 1, 1, 1, 1, 1, ikmax )
+      ! [comm_cleanup] ! DO_3D( 1, 1, 1, 1, 1, ikmax )
+      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, ikmax )
          zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points
          ztm(ji,jj) = ztm(ji,jj) + zc * ts(ji,jj,jk,jp_tem,Kmm)
@@ -2347,5 +2406,6 @@
       ! calculate horizontal gradients at u & v points
 
-      DO_2D( 1, 0, 0, 0 )
+      ! [comm_cleanup] ! DO_2D( 1, 0, 0, 0 )
+      DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 )
          zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj)
          zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj)
@@ -2355,5 +2415,6 @@
       END_2D
 
-      DO_2D( 0, 0, 1, 0 )
+      ! [comm_cleanup] ! DO_2D( 0, 0, 1, 0 )
+      DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 )
          zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
          zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
@@ -2366,12 +2427,15 @@
       CALL eos_rab(ztsm_midv, zmld_midv, zabv, Kmm)
 
-      DO_2D( 1, 0, 0, 0 )
+      ! [comm_cleanup] ! DO_2D( 1, 0, 0, 0 )
+      DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 )
          dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal))
       END_2D
-      DO_2D( 0, 0, 1, 0 )
+      ! [comm_cleanup] ! DO_2D( 0, 0, 1, 0 )
+      DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 )
          dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal))
       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 )
         ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
         zdbds_mle(ji,jj) = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) &
@@ -2400,5 +2464,6 @@
    ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE.
 
-      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 )
        IF ( lconv(ji,jj) ) THEN
           ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
@@ -2409,5 +2474,6 @@
       END_2D
    ! Timestep mixed layer eddy depth.
-      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 )
         IF ( lmle(ji,jj) ) THEN  ! MLE layer growing.
 ! Buoyancy gradient at base of MLE layer.
@@ -2433,8 +2499,11 @@
 
       mld_prof = 4
-      DO_3D( 0, 0, 0, 0, 5, jpkm1 )
+
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 5, jpkm1 )
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 5, jpkm1 )
       IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
       END_3D
-      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 )
          zhmle(ji,jj) = gdepw(ji,jj, mld_prof(ji,jj),Kmm)
       END_2D
@@ -2585,5 +2654,6 @@
          !                                ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case)
          z1_t2 = 2.e-5
-         DO_2D( 1, 1, 1, 1 )
+         ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2)
          END_2D
@@ -2630,5 +2700,6 @@
         etmean(:,:,:) = 0.e0
 
-        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 )
            etmean(ji,jj,jk) = tmask(ji,jj,jk)                     &
                 &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
@@ -2644,5 +2715,6 @@
         etmean(:,:,:) = 0.e0
 
-        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 )
            etmean(ji,jj,jk) = tmask(ji, jj,jk)                           &
                 & / MAX( 1., 2.* tmask(ji,jj,jk)                           &
@@ -2759,5 +2831,6 @@
      !
      hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
-     DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+     ! [comm_cleanup] ! DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+     DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
         ikt = mbkt(ji,jj)
         hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm)
@@ -2765,5 +2838,6 @@
      END_3D
      !
-     DO_2D( 1, 1, 1, 1 )
+     ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
+     DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
         iiki = MAX(4,imld_rst(ji,jj))
         hbl (ji,jj) = gdepw(ji,jj,iiki,Kmm  )    ! Turbocline depth
@@ -2812,5 +2886,6 @@
       ENDIF
 
-      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 )
          pts(ji,jj,jk,jp_tem,Krhs) =  pts(ji,jj,jk,jp_tem,Krhs)                      &
             &                 - (  ghamt(ji,jj,jk  )  &
@@ -2879,5 +2954,6 @@
       !code saving tracer trends removed, replace with trdmxl_oce
 
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )       ! add non-local u and v fluxes
+      ! [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
          puu(ji,jj,jk,Krhs) =  puu(ji,jj,jk,Krhs)                      &
             &                 - (  ghamu(ji,jj,jk  )  &
Index: NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfphy.F90
===================================================================
--- NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfphy.F90	(revision 14598)
+++ NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfphy.F90	(revision 14601)
@@ -54,5 +54,5 @@
    INTEGER, PARAMETER ::   np_OSM = 5   ! OSMOSIS-OBL closure scheme for Kz
 
-   LOGICAL ::   l_zdfsh2   ! shear production term flag (=F for CST, =T otherwise (i.e. TKE, GLS, RIC))
+   LOGICAL, PUBLIC ::   l_zdfsh2   ! shear production term flag (=F for CST, =T otherwise (i.e. TKE, GLS, RIC))
 
    !!----------------------------------------------------------------------
@@ -284,4 +284,5 @@
       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)
@@ -322,14 +323,17 @@
 
       !                                         !* Lateral boundary conditions (sign unchanged)
-      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
+      ! [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
+            ENDIF
          ENDIF
       ENDIF
Index: NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfric.F90
===================================================================
--- NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfric.F90	(revision 14598)
+++ NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfric.F90	(revision 14601)
@@ -156,5 +156,6 @@
       !
       !                       !==  avm and avt = F(Richardson number)  ==!
-      DO_3D( 1, 0, 1, 0, 2, jpkm1 )       ! coefficient = F(richardson number) (avm-weighted Ri)
+      ! [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)
          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
@@ -169,10 +170,12 @@
       IF( ln_mldw ) THEN      !==  set a minimum value in the Ekman layer  ==!
          !
-         DO_2D( 0, 0, 0, 0 )             !* Ekman depth
+         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )             !* Ekman depth
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
             zustar = SQRT( taum(ji,jj) * r1_rho0 )
             zhek   = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )   ! Ekman depth
             zh_ekm(ji,jj) = MAX(  rn_mldmin , MIN( zhek , rn_mldmax )  )   ! set allowed range
          END_2D
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )   !* minimum mixing coeff. within the Ekman layer
+         ! [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
             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_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfsh2.F90
===================================================================
--- NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfsh2.F90	(revision 14598)
+++ NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfsh2.F90	(revision 14601)
@@ -65,5 +65,6 @@
       DO jk = 2, jpkm1                 !* Shear production at uw- and vw-points (energy conserving form)
          IF ( cpl_sdrftx .AND. ln_stshear )  THEN       ! Surface Stokes Drift available  ===>>>  shear + stokes drift contibution
-            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 )
                zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) )        &
                   &         * ( uu (ji,jj,jk-1,Kmm) -   uu (ji,jj,jk,Kmm)    &
@@ -78,5 +79,6 @@
             END_2D
          ELSE
-            DO_2D( 1, 0, 1, 0 )     !* 2 x shear production at uw- and vw-points (energy conserving form)
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )     !* 2 x shear production at uw- and vw-points (energy conserving form)
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )     !* 2 x shear production at uw- and vw-points (energy conserving form)
                zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) &
                   &         * (   uu(ji,jj,jk-1,Kmm) -   uu(ji,jj,jk,Kmm) ) &
@@ -91,5 +93,6 @@
             END_2D
          ENDIF
-         DO_2D( 0, 0, 0, 0 )     !* shear production at w-point ! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked)
+         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )     !* shear production at w-point ! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked)
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )     !* shear production at w-point ! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked)
             p_sh2(ji,jj,jk) = 0.25 * (   ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2. - umask(ji-1,jj,jk) * umask(ji,jj,jk) )   &
                &                       + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2. - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )   )
Index: NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfswm.F90
===================================================================
--- NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfswm.F90	(revision 14598)
+++ NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfswm.F90	(revision 14601)
@@ -63,5 +63,6 @@
       !
       zcoef = 1._wp * 0.353553_wp
-      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 )
          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_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdftke.F90
===================================================================
--- NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdftke.F90	(revision 14598)
+++ NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdftke.F90	(revision 14601)
@@ -241,5 +241,6 @@
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !
-      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 )
          en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) )
          zdiag(ji,jj,1) = 1._wp/en(ji,jj,1)
@@ -258,5 +259,6 @@
       IF( .NOT.ln_drg_OFF ) THEN    !== friction used as top/bottom boundary condition on TKE
          !
-         DO_2D( 0, 0, 0, 0 )        ! bottom friction
+         ! [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
             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)) )
@@ -267,5 +269,6 @@
          END_2D
          IF( ln_isfcav ) THEN
-            DO_2D( 0, 0, 0, 0 )     ! top friction
+            ! [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 
                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)) )
@@ -294,5 +297,6 @@
 !!gm  ! PS: currently we don't have neither the 2 stress components at t-point !nor the angle between u* and u_s
 !!gm  ! so we will overestimate the LC velocity....   !!gm I will do the work if !LC have an effect !
-            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 )
 !!XC                  zWlc2(ji,jj) = 0.5_wp * SQRT( taum(ji,jj) * r1_rho0 * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 )  )
                   zWlc2(ji,jj) = 0.5_wp *  ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 )
@@ -301,5 +305,6 @@
 !  Projection of Stokes drift in the wind stress direction
 !
-            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 )
                   ztaui   = 0.5_wp * ( utau(ji,jj) + utau(ji-1,jj) )
                   ztauj   = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj-1) )
@@ -307,5 +312,6 @@
                   zWlc2(ji,jj) = 0.5_wp * z1_norm * ( MAX( ut0sd(ji,jj)*ztaui + vt0sd(ji,jj)*ztauj, 0._wp ) )**2
             END_2D
-         CALL lbc_lnk      ( 'zdftke', zWlc2, 'T', 1. )
+            ! [comm_cleanup]
+            IF (nn_hls.eq.1) CALL lbc_lnk      ( 'zdftke', zWlc2, 'T', 1. )
 !
          ELSE                          ! Surface Stokes drift deduced from surface stress
@@ -315,5 +321,6 @@
             !                                ! 1/2 Wlc^2 = 0.5 * 0.016 * 0.016 |tau| /( rho_air Cdrag )
             zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )      ! to convert stress in 10m wind using a constant drag
-            DO_2D( 1, 1, 1, 1 )
+            ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
+            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
                zWlc2(ji,jj) = zcof * taum(ji,jj)
             END_2D
@@ -331,18 +338,22 @@
          !                             !- compare LHS to RHS of Eq.47
          imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
-         DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )
+         ! [comm_cleanup] ! DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )
+         DO_3DS( nn_hls, nn_hls, nn_hls, nn_hls, jpkm1, 2, -1 )
             IF( zpelc(ji,jj,jk) > zWlc2(ji,jj) )   imlc(ji,jj) = jk
          END_3D
          !                               ! finite LC depth
-         DO_2D( 1, 1, 1, 1 )
+         ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm)
          END_2D
          !
          zcof = 0.016 / SQRT( zrhoa * zcdrag )
-         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 )
             zus = SQRT( 2. * zWlc2(ji,jj) )             ! Stokes drift
             zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok
          END_2D
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )                  !* TKE Langmuir circulation source term added to en
+         ! [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
             IF ( zus3(ji,jj) /= 0._wp ) THEN
                IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN
@@ -365,5 +376,6 @@
       !
       IF( nn_pdl == 1 ) THEN          !* Prandtl number = F( Ri )
-         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 )
             !                             ! local Richardson number
             IF (rn2b(ji,jj,jk) <= 0.0_wp) then
@@ -377,5 +389,6 @@
       ENDIF
       !
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )   !* Matrix and right hand side in en
+      ! [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
          zcof   = zfact1 * tmask(ji,jj,jk)
          !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical
@@ -406,5 +419,6 @@
 
          CASE ( 0 ) ! Dirichlet BC
-            DO_2D( 0, 0, 0, 0 )    ! en(1)   = rn_ebb taum / rho0  (min value rn_emin0)
+            ! [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)
                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)
@@ -413,5 +427,6 @@
 
          CASE ( 1 ) ! Neumann BC
-            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 )
                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)
@@ -427,5 +442,6 @@
       !
       !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
-      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)
       END_3D
@@ -434,14 +450,18 @@
 !         zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
 !      END_2D
-      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 )
          zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1)
       END_3D
-      DO_2D( 0, 0, 0, 0 )                          ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+      ! [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
          en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
       END_2D
-      DO_3DS( 0, 0, 0, 0, jpk-2, 2, -1 )
+      ! [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 )
          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
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! set the minimum value of tke
+      ! [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
          en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
       END_3D
@@ -456,10 +476,12 @@
       !
       IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
-         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 )
             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_3D
       ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
-         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 )
             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) )   &
@@ -467,5 +489,6 @@
          END_2D
       ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
-         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 )
             ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
             zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
@@ -548,5 +571,6 @@
             zraug = vkarmn * 2.e5_wp / ( rho0 * grav )
 #if ! defined key_si3 && ! defined key_cice
-            DO_2D( 0, 0, 0, 0 )                  ! No sea-ice
+            ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )                  ! No sea-ice
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                  ! No sea-ice
                zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1)
             END_2D
@@ -555,10 +579,12 @@
             !
             CASE( 0 )                      ! No scaling under sea-ice
-               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 )
                   zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)
                END_2D
                !
             CASE( 1 )                      ! scaling with constant sea-ice thickness
-               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 )
                   zmxlm(ji,jj,1) =  ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
                      &                          fr_i(ji,jj)   * rn_mxlice           ) * tmask(ji,jj,1)
@@ -566,5 +592,6 @@
                !
             CASE( 2 )                      ! scaling with mean sea-ice thickness
-               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 )
 #if defined key_si3
                   zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
@@ -578,5 +605,6 @@
                !
             CASE( 3 )                      ! scaling with max sea-ice thickness
-               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 )
                   zmaxice = MAXVAL( h_i(ji,jj,:) )
                   zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
@@ -587,5 +615,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 )
                zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) )
             END_2D
@@ -596,5 +625,6 @@
       ENDIF
       !
-      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 )
          zrn2 = MAX( rn2(ji,jj,jk), rsmall )
          zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  )
@@ -611,5 +641,6 @@
       ! where wmask = 0 set zmxlm == e3w(:,:,:,Kmm)
       CASE ( 0 )           ! bounded by the distance to surface and bottom
-         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 )
             zemxl = MIN( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm), zmxlm(ji,jj,jk),   &
             &            gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - gdepw(ji,jj,jk,Kmm) )
@@ -622,5 +653,6 @@
          !
       CASE ( 1 )           ! bounded by the vertical scale factor
-         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 )
             zemxl = MIN( e3w(ji,jj,jk,Kmm), zmxlm(ji,jj,jk) )
             zmxlm(ji,jj,jk) = zemxl
@@ -629,9 +661,11 @@
          !
       CASE ( 2 )           ! |dk[xml]| bounded by e3t :
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )        ! from the surface to the bottom :
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! from the surface to the bottom :
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! from the surface to the bottom :
             zmxlm(ji,jj,jk) =   &
                &    MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
          END_3D
-         DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )   ! from the bottom to the surface :
+         ! [comm_cleanup] ! DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )   ! from the bottom to the surface :
+         DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 )   ! from the bottom to the surface :
             zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
             zmxlm(ji,jj,jk) = zemxl
@@ -640,13 +674,16 @@
          !
       CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )        ! from the surface to the bottom : lup
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )        ! from the surface to the bottom : lup
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )        ! from the surface to the bottom : lup
             zmxld(ji,jj,jk) =    &
                &    MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
          END_3D
-         DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )   ! from the bottom to the surface : ldown
+         ! [comm_cleanup] ! DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )   ! from the bottom to the surface : ldown
+         DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 )   ! from the bottom to the surface : ldown
             zmxlm(ji,jj,jk) =   &
                &    MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
          END_3D
-         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 )
             zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
             zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
@@ -660,5 +697,6 @@
       !                     !  Vertical eddy viscosity and diffusivity  (avm and avt)
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !* vertical eddy viscosity & diffivity at w-points
+      ! [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
          zsqen = SQRT( en(ji,jj,jk) )
          zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
@@ -670,5 +708,6 @@
       !
       IF( nn_pdl == 1 ) THEN          !* Prandtl number case: update avt
-         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 )
             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
Index: NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/step.F90
===================================================================
--- NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/step.F90	(revision 14598)
+++ NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/step.F90	(revision 14601)
@@ -168,4 +168,13 @@
                          CALL bn2    ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn  ) ! now    Brunt-Vaisala frequency
 
+      ! [comm_cleanup] 
+      IF (nn_hls.eq.2) THEN
+         IF( l_zdfsh2 ) THEN
+            CALL lbc_lnk( 'stp', 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( 'stp', avm  , 'W', 1.0_wp , avt  , 'W', 1.0_wp , avs , 'W', 1.0_wp )
+         ENDIF
+      ENDIF
       !  VERTICAL PHYSICS
                          CALL zdf_phy( kstp, Nbb, Nnn, Nrhs )   ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
Index: NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/stpmlf.F90
===================================================================
--- NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/stpmlf.F90	(revision 14598)
+++ NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/stpmlf.F90	(revision 14601)
@@ -177,4 +177,13 @@
 
       !  VERTICAL PHYSICS
+      ! [comm_cleanup] 
+      IF (nn_hls.eq.2) THEN
+         IF( l_zdfsh2 ) THEN
+            CALL lbc_lnk( 'stp_MLF', 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( 'stp_MLF', avm  , 'W', 1.0_wp , avt  , 'W', 1.0_wp , avs , 'W', 1.0_wp )
+         ENDIF
+      ENDIF
                          CALL zdf_phy( kstp, Nbb, Nnn, Nrhs )   ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
 
@@ -309,6 +318,6 @@
 #endif
 
-     ! [comm_cleanup]
-     IF (nn_hls.EQ.2) THEN
+      ! [comm_cleanup]
+      IF (nn_hls.EQ.2) THEN
          SELECT CASE ( nadv )
          CASE ( np_FCT )                                 ! FCT scheme : 2nd / 4th order
