Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynadv.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynadv.F90	(revision 14711)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynadv.F90	(revision 14712)
@@ -76,8 +76,6 @@
          CALL dyn_zad     ( kt,                 Kmm, puu, pvv, Krhs )    ! vector form : vertical advection
       CASE( np_FLX_c2  ) 
-         ! [comm_cleanup: dyn_adv_cen2 NOT TESTED] 
          CALL dyn_adv_cen2( kt,                 Kmm, puu, pvv, Krhs )    ! 2nd order centered scheme
       CASE( np_FLX_ubs )   
-         ! [comm_cleanup: dyn_adv_ubs NOT TESTED] 
          CALL dyn_adv_ubs ( kt,            Kbb, Kmm, puu, pvv, Krhs )    ! 3rd order UBS      scheme (UP3)
       END SELECT
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynadv_cen2.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynadv_cen2.F90	(revision 14711)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynadv_cen2.F90	(revision 14712)
@@ -72,6 +72,5 @@
          zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm)
          zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
-         ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )              ! horizontal momentum fluxes (at T- and F-point)
-         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )              ! horizontal momentum fluxes (at T- and F-point)
+         DO_2D( 1, 0, 1, 0 )              ! horizontal momentum fluxes (at T- and F-point)
             zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj  ,jk,Kmm) )
             zfv_f(ji  ,jj  ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) )
@@ -79,6 +78,5 @@
             zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji  ,jj+1,jk,Kmm) )
          END_2D
-         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )              ! divergence of horizontal momentum fluxes
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )              ! divergence of horizontal momentum fluxes
+         DO_2D( 0, 0, 0, 0 )              ! divergence of horizontal momentum fluxes
             puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
                &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj)   &
@@ -100,12 +98,10 @@
       !                             !==  Vertical advection  ==!
       !
-      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )                          ! surface/bottom advective fluxes set to zero
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                          ! surface/bottom advective fluxes set to zero
+      DO_2D( 0, 0, 0, 0 )                 ! surface/bottom advective fluxes set to zero
          zfu_uw(ji,jj,jpk) = 0._wp   ;   zfv_vw(ji,jj,jpk) = 0._wp
          zfu_uw(ji,jj, 1 ) = 0._wp   ;   zfv_vw(ji,jj, 1 ) = 0._wp
       END_2D
       IF( ln_linssh ) THEN                ! linear free surface: advection through the surface
-         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D( 0, 0, 0, 0 )
             zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
             zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji,jj+1) * ww(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
@@ -113,16 +109,13 @@
       ENDIF
       DO jk = 2, jpkm1                    ! interior advective fluxes
-         ! [comm_cleanup] ! DO_2D( 0, 1, 0, 1 )                  ! 1/4 * Vertical transport
-         DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )                  ! 1/4 * Vertical transport 
+         DO_2D( 0, 1, 0, 1 )                  ! 1/4 * Vertical transport
             zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
          END_2D
-         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D( 0, 0, 0, 0 )
             zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj  ,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) )
             zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji  ,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) )
          END_2D
       END DO
-      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )       ! divergence of vertical momentum flux divergence
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )       ! divergence of vertical momentum flux divergence
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )       ! divergence of vertical momentum flux divergence
          puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj)   &
             &                                      / e3u(ji,jj,jk,Kmm)
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynadv_ubs.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynadv_ubs.F90	(revision 14711)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynadv_ubs.F90	(revision 14712)
@@ -108,6 +108,7 @@
          zfv(:,:,jk) = e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
          !            
-         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )                       ! laplacian
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                       ! laplacian
+         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )                       ! laplacian 
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                       ! laplacia
+
             zlu_uu(ji,jj,jk,1) = ( puu (ji+1,jj  ,jk,Kbb) - 2.*puu (ji,jj,jk,Kbb) + puu (ji-1,jj  ,jk,Kbb) ) * umask(ji,jj,jk)
             zlv_vv(ji,jj,jk,1) = ( pvv (ji  ,jj+1,jk,Kbb) - 2.*pvv (ji,jj,jk,Kbb) + pvv (ji  ,jj-1,jk,Kbb) ) * vmask(ji,jj,jk)
@@ -137,6 +138,5 @@
          zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
          !
-         ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )              ! horizontal momentum fluxes at T- and F-point
-         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )              ! horizontal momentum fluxes at T- and F-point
+         DO_2D( 1, 0, 1, 0 )                       ! horizontal momentum fluxes at T- and F-point
             zui = ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj  ,jk,Kmm) )
             zvj = ( pvv(ji,jj,jk,Kmm) + pvv(ji  ,jj+1,jk,Kmm) )
@@ -170,6 +170,5 @@
                &                * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) - gamma1 * zl_v )
          END_2D
-         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )              ! divergence of horizontal momentum fluxes
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )              ! divergence of horizontal momentum fluxes
+         DO_2D( 0, 0, 0, 0 )                       ! divergence of horizontal momentum fluxes
             puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
                &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj)   &
@@ -190,6 +189,5 @@
       !                                      !  Vertical advection  !
       !                                      ! ==================== !
-      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )                          ! surface/bottom advective fluxes set to zero
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                          ! surface/bottom advective fluxes set to zero
+      DO_2D( 0, 0, 0, 0 )                          ! surface/bottom advective fluxes set to zero
          zfu_uw(ji,jj,jpk) = 0._wp
          zfv_vw(ji,jj,jpk) = 0._wp
@@ -198,6 +196,5 @@
       END_2D
       IF( ln_linssh ) THEN                         ! constant volume : advection through the surface
-         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D( 0, 0, 0, 0 )
             zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
             zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji,jj+1) * ww(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
@@ -205,16 +202,13 @@
       ENDIF
       DO jk = 2, jpkm1                          ! interior fluxes
-         ! [comm_cleanup] ! DO_2D( 0, 1, 0, 1 )
-         DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls ) 
+         DO_2D( 0, 1, 0, 1 )
             zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
          END_2D
-         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D( 0, 0, 0, 0 )
             zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) )
             zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) )
          END_2D
       END DO
-      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )       ! divergence of vertical momentum flux divergence
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )       ! divergence of vertical momentum flux divergence
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )             ! divergence of vertical momentum flux divergence
          puu(ji,jj,jk,Krhs) =  puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj)   &
             &                                       / e3u(ji,jj,jk,Kmm)
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynatf.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynatf.F90	(revision 14711)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynatf.F90	(revision 14712)
@@ -169,5 +169,5 @@
 # endif
       !
-      CALL lbc_lnk( 'dynatf', puu(:,:,:,Kaa), 'U', -1.0_wp, pvv(:,:,:,Kaa), 'V', -1.0_wp )     !* local domain boundaries
+      IF (nn_hls.eq.1) CALL lbc_lnk( 'dynatf', puu(:,:,:,Kaa), 'U', -1.0_wp, pvv(:,:,:,Kaa), 'V', -1.0_wp )     !* local domain boundaries
       !
       !                                !* BDY open boundaries
@@ -201,5 +201,6 @@
          IF( ln_linssh ) THEN             ! Fixed volume !
             !                             ! =============!
-            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 )
                puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + rn_atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )
                pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + rn_atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )
@@ -237,5 +238,6 @@
                CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3u(:,:,:,Kmm), 'U' )
                CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3v(:,:,:,Kmm), 'V' )
-               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 )
                   puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + rn_atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )
                   pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + rn_atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )
@@ -248,5 +250,6 @@
                CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3u_f, 'U' )
                CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3v_f, 'V' )
-               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 )
                   zue3a = pe3u(ji,jj,jk,Kaa) * puu(ji,jj,jk,Kaa)
                   zve3a = pe3v(ji,jj,jk,Kaa) * pvv(ji,jj,jk,Kaa)
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynhpg.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynhpg.F90	(revision 14711)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynhpg.F90	(revision 14712)
@@ -462,5 +462,7 @@
           END IF
         END_2D
-        CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )
+        ! [ comm_cleanup ] ! I think lbc_lnk can be deleted for halo 1 case, too:
+        ! zcpx and zcpy are written and used only in the inner domain - can't test it 
+        ! CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )
       END IF
       !
@@ -689,5 +691,7 @@
           END IF
         END_2D
-        CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )
+        ! [ comm_cleanup ] ! I think lbc_lnk can be deleted for halo 1 case, too:
+        ! zcpx and zcpy are written and used only in the inner domain - can't test it 
+        ! CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )
       END IF
 
@@ -786,5 +790,6 @@
       !----------------------------------------------------------------------------------------
 
-      DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+      ! [ comm_cleanup ] ! DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
          zdrhox(ji,jj,jk) =   rhd    (ji+1,jj  ,jk) - rhd    (ji,jj,jk  )
          zdzx  (ji,jj,jk) = - gde3w(ji+1,jj  ,jk) + gde3w(ji,jj,jk  )
@@ -1043,5 +1048,7 @@
             ENDIF
          END_2D
-         CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )
+         ! [ comm_cleanup ] ! I think lbc_lnk can be deleted for halo 1 case, too:
+         ! zcpx and zcpy are written and used only in the inner domain - can't test it 
+         ! CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )
       ENDIF
 
@@ -1113,5 +1120,7 @@
       END_2D
 
-      CALL lbc_lnk ('dynhpg', zsshu_n, 'U', 1.0_wp, zsshv_n, 'V', 1.0_wp )
+      ! [ comm_cleanup ] ! I think lbc_lnk can be deleted for halo 1 case, too:
+      ! zcpx and zcpy are written and used only in the inner domain 
+      ! CALL lbc_lnk ('dynhpg', zsshu_n, 'U', 1.0_wp, zsshv_n, 'V', 1.0_wp )
 
       DO_2D( 0, 0, 0, 0 )
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynkeg.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynkeg.F90	(revision 14711)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynkeg.F90	(revision 14712)
@@ -101,6 +101,5 @@
       !
       CASE ( nkeg_C2 )                          !--  Standard scheme  --!
-         ! [comm_cleanup] ! DO_3D( 0, 1, 0, 1, 1, jpkm1 )
-         DO_3D( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 )
+         DO_3D( 0, 1, 0, 1, 1, jpkm1 )
             zu =    puu(ji-1,jj  ,jk,Kmm) * puu(ji-1,jj  ,jk,Kmm)   &
                &  + puu(ji  ,jj  ,jk,Kmm) * puu(ji  ,jj  ,jk,Kmm)
@@ -110,15 +109,19 @@
          END_3D
       CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
-         ! [comm_cleanup: Hollingsworth scheme NOT TESTED] ! 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 )
+            ! round brackets added to fix the order of floating point operations
+            ! needed to ensure halo 1 - halo 2 compatibility
             zu = 8._wp * ( puu(ji-1,jj  ,jk,Kmm) * puu(ji-1,jj  ,jk,Kmm)    &
                &         + puu(ji  ,jj  ,jk,Kmm) * puu(ji  ,jj  ,jk,Kmm) )  &
-               &   +     ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) ) * ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) )   &
-               &   +     ( puu(ji  ,jj-1,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) ) * ( puu(ji  ,jj-1,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) )
+               &   +     ( ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) ) * ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) )   &
+               &   +       ( puu(ji  ,jj-1,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) ) * ( puu(ji  ,jj-1,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) )   &
+               &         )                                                               ! bracket for halo 1 - halo 2 compatibility
                !
             zv = 8._wp * ( pvv(ji  ,jj-1,jk,Kmm) * pvv(ji  ,jj-1,jk,Kmm)    &
                &         + pvv(ji  ,jj  ,jk,Kmm) * pvv(ji  ,jj  ,jk,Kmm) )  &
-               &  +      ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) ) * ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) )   &
-               &  +      ( pvv(ji-1,jj  ,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) ) * ( pvv(ji-1,jj  ,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) )
+               &  +      ( ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) ) * ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) )  &
+               &  +        ( pvv(ji-1,jj  ,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) ) * ( pvv(ji-1,jj  ,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) )  & 
+               &         )                                                               ! bracket for halo 1 - halo 2 compatibility
             zhke(ji,jj,jk) = r1_48 * ( zv + zu )
          END_3D
@@ -127,6 +130,5 @@
       END SELECT 
       !
-      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )       !==  grad( KE ) added to the general momentum trends  ==!
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )       !==  grad( KE ) added to the general momentum trends  ==!
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )       !==  grad( KE ) added to the general momentum trends  ==!
          puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
          pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynldf_iso.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynldf_iso.F90	(revision 14711)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynldf_iso.F90	(revision 14712)
@@ -128,5 +128,6 @@
       IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
          !
-         DO_3D( 0, 0, 0, 0, 1, jpk )      ! set the slopes of iso-level
+         ! [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 )      ! set the slopes of iso-level 
             uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
             vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
@@ -135,5 +136,5 @@
          END_3D
          ! Lateral boundary conditions on the slopes
-         CALL lbc_lnk( 'dynldf_iso', uslp , 'U', -1.0_wp, vslp , 'V', -1.0_wp, wslpi, 'W', -1.0_wp, wslpj, 'W', -1.0_wp )
+         IF (nn_hls.eq.1) CALL lbc_lnk( 'dynldf_iso', uslp , 'U', -1.0_wp, vslp , 'V', -1.0_wp, wslpi, 'W', -1.0_wp, wslpj, 'W', -1.0_wp )
          !
        ENDIF
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynldf_lap_blp.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynldf_lap_blp.F90	(revision 14711)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynldf_lap_blp.F90	(revision 14712)
@@ -83,5 +83,6 @@
          DO jk = 1, jpkm1                                 ! Horizontal slab
             !
-            DO_2D( 0, 1, 0, 1 )
+            ! [comm_cleanup] ! DO_2D( 0, 1, 0, 1 ) 
+            DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
                !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1)
                zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &   ! ahmf already * by fmask
@@ -94,5 +95,6 @@
             END_2D
             !
-            DO_2D( 0, 0, 0, 0 )                       ! - curl( curl) + grad( div )
+            ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )                       ! - curl( curl) + grad( div ) 
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                       ! - curl( curl) + grad( div )
                pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * umask(ji,jj,jk) * (    &    ! * by umask is mandatory for dyn_ldf_blp use
                   &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   &
@@ -114,5 +116,6 @@
          DO jk = 1, jpkm1                                 ! Horizontal slab
             !
-            DO_2D( 0, 1, 0, 1 )
+            ! [comm_cleanup] ! DO_2D( 0, 1, 0, 1 ) 
+            DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
                !                                      ! shearing stress component (F-point)   NB : ahmf has already been multiplied by fmask
                zshe(ji-1,jj-1) = ahmf(ji-1,jj-1,jk)                                                              &
@@ -129,5 +132,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 )
                pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                               &
                   &    * (   (   zten(ji+1,jj  ) * e2t(ji+1,jj  )*e2t(ji+1,jj  ) * e3t(ji+1,jj  ,jk,Kmm)                       &
@@ -185,5 +189,5 @@
       CALL dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, zulap, zvlap, 1 )   ! rotated laplacian applied to pt (output in zlap,Kbb)
       !
-      CALL lbc_lnk( 'dynldf_lap_blp', zulap, 'U', -1.0_wp, zvlap, 'V', -1.0_wp )             ! Lateral boundary conditions
+      IF (nn_hls.eq.1) CALL lbc_lnk( 'dynldf_lap_blp', zulap, 'U', -1.0_wp, zvlap, 'V', -1.0_wp )             ! Lateral boundary conditions
       !
       CALL dyn_ldf_lap( kt, Kbb, Kmm, zulap, zvlap, pu_rhs, pv_rhs, 2 )   ! rotated laplacian applied to zlap (output in pt(:,:,:,:,Krhs))
Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynvor.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynvor.F90	(revision 14711)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynvor.F90	(revision 14712)
@@ -256,15 +256,17 @@
          ALLOCATE( zwz(jpi,jpj,jpk) )
          DO jk = 1, jpkm1                                ! Horizontal slab
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
                   &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
             END_2D
             IF( ln_dynvor_msk ) THEN                     ! mask relative vorticity
-               DO_2D( 1, 0, 1, 0 )
+               ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
+               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                   zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
                END_2D
             ENDIF
          END DO
-         CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
+         IF (nn_hls.eq.1) CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
          !
       END SELECT
@@ -625,5 +627,6 @@
          !
 #if defined key_qco   ||   defined key_linssh
-         DO_2D( 1, 0, 1, 0 )                 ! == reciprocal of e3 at F-point (key_qco)
+         ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )                 ! == reciprocal of e3 at F-point (key_qco) 
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                 ! == reciprocal of e3 at F-point (key_qco)
             z1_e3f(ji,jj) = 1._wp / e3f_vor(ji,jj,jk)
          END_2D
@@ -631,5 +634,6 @@
          SELECT CASE( nn_e3f_typ )           ! == reciprocal of e3 at F-point
          CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
                   &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
@@ -641,5 +645,6 @@
             END_2D
          CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
                   &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
@@ -658,35 +663,46 @@
          !
          CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj)
             END_2D
          CASE ( np_RVO )                           !* relative vorticity
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                zwz(ji,jj,jk) = ( e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
                   &            - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
             END_2D
             IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
-               DO_2D( 1, 0, 1, 0 )
+               ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
+               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                   zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
                END_2D
             ENDIF
          CASE ( np_MET )                           !* metric term
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                zwz(ji,jj,jk) = (   ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
                   &              - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
             END_2D
          CASE ( np_CRV )                           !* Coriolis + relative vorticity
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
-                  &                              - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  )   &
-                  &                           * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            ! round brackets added to fix the order of floating point operations
+            ! needed to ensure halo 1 - halo 2 compatibility
+               zwz(ji,jj,jk) = (  ff_f(ji,jj) + ( ( e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
+                  &                               )                                                                  & ! bracket for halo 1 - halo 2 compatibility
+                  &                             - ( e1u(ji  ,jj+1) * pu(ji,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk)      &
+                  &                               )                                                                  & ! bracket for halo 1 - halo 2 compatibility
+                  &                             ) * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
             END_2D
             IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
-               DO_2D( 1, 0, 1, 0 )
+               ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
+               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                   zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
                END_2D
             ENDIF
          CASE ( np_CME )                           !* Coriolis + metric
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                zwz(ji,jj,jk) = (   ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
                   &                            - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
@@ -699,5 +715,5 @@
       !                                                ! ===============
       !
-      CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
+      IF (nn_hls.eq.1) CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
       !
       !                                                ! ===============
@@ -776,9 +792,11 @@
          SELECT CASE( kvor )                 !==  vorticity considered  ==!
          CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                zwz(ji,jj,jk) = ff_f(ji,jj)
             END_2D
          CASE ( np_RVO )                           !* relative vorticity
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                zwz(ji,jj,jk) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
                   &             - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
@@ -786,15 +804,18 @@
             END_2D
             IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
-               DO_2D( 1, 0, 1, 0 )
+               ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
+               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                   zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
                END_2D
             ENDIF
          CASE ( np_MET )                           !* metric term
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                zwz(ji,jj,jk) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
                   &          - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
             END_2D
          CASE ( np_CRV )                           !* Coriolis + relative vorticity
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
                   &                              - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
@@ -802,10 +823,12 @@
             END_2D
             IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
-               DO_2D( 1, 0, 1, 0 )
+               ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
+               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                   zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
                END_2D
             ENDIF
          CASE ( np_CME )                           !* Coriolis + metric
-            DO_2D( 1, 0, 1, 0 )
+            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
                   &                        - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
@@ -819,5 +842,5 @@
       !                                                ! ===============
       !
-      CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
+      IF (nn_hls.eq.1) CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
       !
       !                                                ! ===============
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 14711)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traatf.F90	(revision 14712)
@@ -111,5 +111,5 @@
       !                                              ! local domain boundaries  (T-point, unchanged sign)
       ! [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 (nn_hls.eq.1) 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
@@ -158,5 +158,5 @@
          !
          ! [comm_cleanup] 
-         ! CALL lbc_lnk( 'traatf',  pts(:,:,:,jp_tem,Kmm) , 'T', 1.0_wp, pts(:,:,:,jp_sal,Kmm) , 'T', 1.0_wp )
+         IF (nn_hls.eq.1) 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/ZDF/zdfosm.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfosm.F90	(revision 14711)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfosm.F90	(revision 14712)
@@ -400,6 +400,5 @@
      zz0 =       rn_abs       ! surface equi-partition in 2-bands
      zz1 =  1. - rn_abs
-     ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-     DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+     DO_2D( 0, 0, 0, 0 )
         ! Surface downward irradiance (so always +ve)
         zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp
@@ -411,6 +410,5 @@
      END_2D
      ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL
-     ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-     DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+     DO_2D( 0, 0, 0, 0 )
         zthermal = rab_n(ji,jj,1,jp_tem)
         zbeta    = rab_n(ji,jj,1,jp_sal)
@@ -439,6 +437,5 @@
      ! Assume constant La#=0.3
      CASE(0)
-        ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-        DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+        DO_2D( 0, 0, 0, 0 )
            zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2
            zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2
@@ -449,6 +446,5 @@
      ! Assume Pierson-Moskovitz wind-wave spectrum
      CASE(1)
-        ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-        DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+        DO_2D( 0, 0, 0, 0 )
            ! Use wind speed wndm included in sbc_oce module
            zustke(ji,jj) =  MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )
@@ -459,6 +455,5 @@
         zfac =  2.0_wp * rpi / 16.0_wp
 
-        ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-        DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+        DO_2D( 0, 0, 0, 0 )
            IF (hsw(ji,jj) > 1.e-4) THEN
               ! Use  wave fields
@@ -477,6 +472,5 @@
      IF (ln_zdfosm_ice_shelter) THEN
         ! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice
-        ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-        DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+        DO_2D( 0, 0, 0, 0 )
            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))
@@ -500,6 +494,5 @@
         z_two_thirds = 2.0_wp / 3.0_wp
 
-        ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-        DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+        DO_2D( 0, 0, 0, 0 )
            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 )
@@ -516,6 +509,5 @@
         zsqrtpi = SQRT(rpi)
 
-        ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-        DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+        DO_2D( 0, 0, 0, 0 )
            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 )
@@ -538,6 +530,5 @@
      ! Langmuir velocity scale (zwstrl), La # (zla)
      ! mixed scale (zvstr), convective velocity scale (zwstrc)
-     ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-     DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+     DO_2D( 0, 0, 0, 0 )
         ! Langmuir velocity scale (zwstrl), at T-point
         zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird
@@ -572,6 +563,5 @@
       hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,4,Kmm) )
       ibld(:,:) = 4
-      ! [comm_cleanup] ! DO_3D( 1, 1, 1, 1, 5, jpkm1 )
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 5, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 5, jpkm1 )
          IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN
             ibld(ji,jj) = MIN(mbkt(ji,jj), jk)
@@ -580,6 +570,5 @@
      ! ##########################################################################
 
-      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D( 0, 0, 0, 0 )
          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 ))
@@ -601,6 +590,5 @@
 ! Fox-Kemper Scheme
          mld_prof = 4
-         ! [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 )
+         DO_3D( 0, 0, 0, 0, 5, jpkm1 )
          IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
          END_3D
@@ -608,6 +596,5 @@
         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)
 
-         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D( 0, 0, 0, 0 )
            zhmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm)
          END_2D
@@ -624,6 +611,5 @@
          lflux(:,:) = .FALSE.
          lmle(:,:) = .FALSE.
-         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D( 0, 0, 0, 0 )
           IF ( lconv(ji,jj) .AND. zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE.
          END_2D
@@ -631,6 +617,5 @@
 
 ! Test if pycnocline well resolved
-      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D( 0, 0, 0, 0 )
        IF (lconv(ji,jj) ) THEN
           ztmp = 0.2 * zhbl(ji,jj) / e3w(ji,jj,ibld(ji,jj),Kmm)
@@ -653,6 +638,5 @@
 ! Rate of change of hbl
       CALL zdf_osm_calculate_dhdt( zdhdt, zddhdt )
-      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D( 0, 0, 0, 0 )
        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
@@ -666,6 +650,5 @@
       ibld(:,:) = 4
 
-      ! [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 )
+      DO_3D( 0, 0, 0, 0, 4, jpkm1 )
          IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN
             ibld(ji,jj) = jk
@@ -686,6 +669,5 @@
       CALL zdf_osm_pycnocline_thickness( dh, zdh )
 
-      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D( 0, 0, 0, 0 )
        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
@@ -727,6 +709,5 @@
 
 
-       ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+       DO_2D( 0, 0, 0, 0 )
          IF ( lconv(ji,jj) ) THEN
            DO jk = 2, imld(ji,jj)
@@ -763,6 +744,5 @@
           IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv )
        END IF
-       ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+       DO_2D( 0, 0, 0, 0 )
           IF ( lconv(ji,jj) ) THEN
              DO jk = 2, imld(ji,jj)
@@ -796,6 +776,5 @@
        ENDWHERE
 
-       ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+       DO_2D( 0, 0, 0, 0 )
           IF (lconv(ji,jj) ) THEN
              DO jk = 2, imld(ji,jj)
@@ -859,6 +838,5 @@
        ENDWHERE
 
-       ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+       DO_2D( 0, 0, 0, 0 )
           IF ( lconv(ji,jj) ) THEN
              DO jk = 2 , imld(ji,jj)
@@ -878,6 +856,5 @@
        END_2D
 
-       ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+       DO_2D( 0, 0, 0, 0 )
         IF ( lpyc(ji,jj) ) THEN
           IF ( j_ddh(ji,jj) == 0 ) THEN
@@ -914,6 +891,5 @@
 ! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ]
 
-       ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
-       DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
+       DO_2D( 1, 0, 1, 0 )
 
          IF ( lconv(ji,jj) ) THEN
@@ -931,6 +907,5 @@
        END_2D
 
-       ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+       DO_2D( 0, 0, 0, 0 )
          IF ( lconv(ji,jj) ) THEN
             DO jk = 2, imld(ji,jj)
@@ -979,6 +954,5 @@
        ENDWHERE
 
-       ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+       DO_2D( 0, 0, 0, 0 )
           IF ( lconv(ji,jj) ) THEN
             DO jk = 2, imld(ji,jj)
@@ -1027,6 +1001,5 @@
  ! Make surface forced velocity non-gradient terms go to zero at the base of the boundary layer.
 
-      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D( 0, 0, 0, 0 )
          IF ( .not. lconv(ji,jj) ) THEN
             DO jk = 2, ibld(ji,jj)
@@ -1044,6 +1017,5 @@
 
       ! pynocline contributions
-       ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+       DO_2D( 0, 0, 0, 0 )
          IF ( .not. lconv(ji,jj) ) THEN
           IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
@@ -1063,6 +1035,5 @@
        END IF
 
-       ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+       DO_2D( 0, 0, 0, 0 )
           ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
           ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
@@ -1089,6 +1060,5 @@
        ! rotate non-gradient velocity terms back to model reference frame
 
-       ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+       DO_2D( 0, 0, 0, 0 )
           DO jk = 2, ibld(ji,jj)
              ztemp = ghamu(ji,jj,jk)
@@ -1106,7 +1076,5 @@
 ! KPP-style Ri# mixing
        IF( ln_kpprimix) THEN
-          
-          ! [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 )
+          DO_3D( 1, 0, 1, 0, 2, jpkm1 )      !* Shear production at uw- and vw-points (energy conserving form)
              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) &
@@ -1117,6 +1085,5 @@
           END_3D
       !
-         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 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) )   &
@@ -1129,6 +1096,5 @@
          END_3D
 
-          ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-          DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+          DO_2D( 0, 0, 0, 0 )
              DO jk = ibld(ji,jj) + 1, jpkm1
                 zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
@@ -1141,6 +1107,5 @@
 ! KPP-style set diffusivity large if unstable below BL
        IF( ln_convmix) THEN
-          ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-          DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+          DO_2D( 0, 0, 0, 0 )
              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
@@ -1152,6 +1117,5 @@
 
        IF ( ln_osm_mle ) THEN  ! set up diffusivity and non-gradient mixing
-          ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-          DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+          DO_2D( 0, 0, 0, 0 )
               IF ( lflux(ji,jj) ) THEN ! MLE mixing extends below boundary layer
              ! Calculate MLE flux contribution from surface fluxes
@@ -1194,12 +1158,11 @@
        ! GN 25/8: need to change tmask --> wmask
 
-     ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-     DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
+     DO_3D( 0, 0, 0, 0, 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
-     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 )
+     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) ) &
@@ -1213,12 +1176,9 @@
        END_3D
         ! Lateral boundary conditions on final outputs for hbl,  on T-grid (sign unchanged)
-        ! [comm_cleanup] ! no need lbc_lnk for output
-        ! NOTE: [tiling] still needed (at least for hmle in tra_mle_trp)
-         CALL lbc_lnk( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. )
+        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)
-        ! [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 )
+        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
@@ -1320,6 +1280,5 @@
       REAL(wp), PARAMETER :: rn_vispyc_shr = 0.15
 
-      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D( 0, 0, 0, 0 )
           IF ( lconv(ji,jj) ) THEN
 
@@ -1364,6 +1323,5 @@
       END_2D
 !
-       ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+       DO_2D( 0, 0, 0, 0 )
           IF ( lconv(ji,jj) ) THEN
              DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity
@@ -1464,6 +1422,5 @@
 
 ! Determins stability and set flag lconv
-     ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-     DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+     DO_2D( 0, 0, 0, 0 )
       IF ( zhol(ji,jj) < 0._wp ) THEN
          lconv(ji,jj) = .TRUE.
@@ -1482,6 +1439,5 @@
      j_ddh(:,:) = 1
 
-     ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-     DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+     DO_2D( 0, 0, 0, 0 )
       IF ( lconv(ji,jj) ) THEN
          IF ( zdb_bl(ji,jj) > 0._wp ) THEN
@@ -1520,6 +1476,5 @@
 ! Calculate entrainment buoyancy flux due to surface fluxes.
 
-     ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-     DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+     DO_2D( 0, 0, 0, 0 )
       IF ( lconv(ji,jj) ) THEN
         zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln
@@ -1546,6 +1501,5 @@
      zwb_min(:,:) = 0._wp
 
-     ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-     DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+     DO_2D( 0, 0, 0, 0 )
       IF ( lshear(ji,jj) ) THEN
         IF ( lconv(ji,jj) ) THEN
@@ -1608,6 +1562,5 @@
         zu   = 0._wp
         zv   = 0._wp
-        ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-        DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+        DO_2D( 0, 0, 0, 0 )
          zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
          zbeta    = rab_n(ji,jj,1,jp_sal)
@@ -1666,6 +1619,5 @@
         REAL(wp) :: ztemp
 
-        ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-        DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+        DO_2D( 0, 0, 0, 0 )
            ztemp = zu(ji,jj)
            zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj)
@@ -1701,11 +1653,9 @@
       znd_param(:,:) = 0._wp
 
-        ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-        DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+        DO_2D( 0, 0, 0, 0 )
           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
-        ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-        DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+        DO_2D( 0, 0, 0, 0 )
                  !
          IF ( lconv(ji,jj) ) THEN
@@ -1731,6 +1681,5 @@
 
 ! Diagnosis
-        ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-        DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+        DO_2D( 0, 0, 0, 0 )
           IF ( lconv(ji,jj) ) THEN
             zwb_ent = - 2.0 * 0.2 * zwbav(ji,jj) &
@@ -1802,6 +1751,5 @@
 
 
-     ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-     DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+     DO_2D( 0, 0, 0, 0 )
         IF ( jbase(ji,jj)+1 < mbkt(ji,jj) ) THEN
            zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
@@ -1833,6 +1781,5 @@
      REAL(wp), PARAMETER :: zgamma_b = 2.25, zzeta_sh = 0.15
 
-     ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-     DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+     DO_2D( 0, 0, 0, 0 )
         IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
            IF ( lconv(ji,jj) ) THEN  ! convective conditions
@@ -1919,6 +1866,5 @@
       REAL(wp) :: zzeta_v = 0.45
       !
-      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D( 0, 0, 0, 0 )
          !
          IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
@@ -1983,6 +1929,5 @@
     REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth
 
-  ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-  DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+  DO_2D( 0, 0, 0, 0 )
 
     IF ( lshear(ji,jj) ) THEN
@@ -2127,6 +2072,5 @@
     REAL(wp) :: zthermal, zbeta
 
-     ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-     DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+     DO_2D( 0, 0, 0, 0 )
         IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN
 !
@@ -2232,6 +2176,5 @@
       REAL, PARAMETER :: a_ddh_2 = 3.5 ! also in pycnocline_depth
 
-    ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-    DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+    DO_2D( 0, 0, 0, 0 )
 
       IF ( lshear(ji,jj) ) THEN
@@ -2379,12 +2322,10 @@
       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
-      ! [comm_cleanup] ! DO_3D( 1, 1, 1, 1, nlb10, jpkm1 )
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, nlb10, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 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
-      ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 1, 1, 1, 1 )
          mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj))
          zmld(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm)
@@ -2396,6 +2337,5 @@
       ztm(:,:) = 0._wp
       zsm(:,:) = 0._wp
-      ! [comm_cleanup] ! DO_3D( 1, 1, 1, 1, 1, ikmax )
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, ikmax )
+      DO_3D( 1, 1, 1, 1, 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)
@@ -2407,6 +2347,5 @@
       ! calculate horizontal gradients at u & v points
 
-      ! [comm_cleanup] ! DO_2D( 1, 0, 0, 0 )
-      DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D( 1, 0, 0, 0 )
          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)
@@ -2416,6 +2355,5 @@
       END_2D
 
-      ! [comm_cleanup] ! DO_2D( 0, 0, 1, 0 )
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 )
+      DO_2D( 0, 0, 1, 0 )
          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)
@@ -2428,15 +2366,12 @@
       CALL eos_rab(ztsm_midv, zmld_midv, zabv, Kmm)
 
-      ! [comm_cleanup] ! DO_2D( 1, 0, 0, 0 )
-      DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D( 1, 0, 0, 0 )
          dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal))
       END_2D
-      ! [comm_cleanup] ! DO_2D( 0, 0, 1, 0 )
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 )
+      DO_2D( 0, 0, 1, 0 )
          dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal))
       END_2D
 
-      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D( 0, 0, 0, 0 )
         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) &
@@ -2465,6 +2400,5 @@
    ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE.
 
-      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D( 0, 0, 0, 0 )
        IF ( lconv(ji,jj) ) THEN
           ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
@@ -2475,6 +2409,5 @@
       END_2D
    ! Timestep mixed layer eddy depth.
-      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D( 0, 0, 0, 0 )
         IF ( lmle(ji,jj) ) THEN  ! MLE layer growing.
 ! Buoyancy gradient at base of MLE layer.
@@ -2500,11 +2433,8 @@
 
       mld_prof = 4
-
-      ! [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 )
+      DO_3D( 0, 0, 0, 0, 5, jpkm1 )
       IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
       END_3D
-      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D( 0, 0, 0, 0 )
          zhmle(ji,jj) = gdepw(ji,jj, mld_prof(ji,jj),Kmm)
       END_2D
@@ -2655,6 +2585,5 @@
          !                                ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case)
          z1_t2 = 2.e-5
-         ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         DO_2D( 1, 1, 1, 1 )
             r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2)
          END_2D
@@ -2701,6 +2630,5 @@
         etmean(:,:,:) = 0.e0
 
-        ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-        DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
+        DO_3D( 0, 0, 0, 0, 1, jpkm1 )
            etmean(ji,jj,jk) = tmask(ji,jj,jk)                     &
                 &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
@@ -2716,6 +2644,5 @@
         etmean(:,:,:) = 0.e0
 
-        ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-        DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
+        DO_3D( 0, 0, 0, 0, 1, jpkm1 )
            etmean(ji,jj,jk) = tmask(ji, jj,jk)                           &
                 & / MAX( 1., 2.* tmask(ji,jj,jk)                           &
@@ -2832,6 +2759,5 @@
      !
      hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
-     ! [comm_cleanup] ! DO_3D( 1, 1, 1, 1, 1, jpkm1 )
-     DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+     DO_3D( 1, 1, 1, 1, 1, jpkm1 )
         ikt = mbkt(ji,jj)
         hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm)
@@ -2839,6 +2765,5 @@
      END_3D
      !
-     ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
-     DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+     DO_2D( 1, 1, 1, 1 )
         iiki = MAX(4,imld_rst(ji,jj))
         hbl (ji,jj) = gdepw(ji,jj,iiki,Kmm  )    ! Turbocline depth
@@ -2875,5 +2800,5 @@
       !
       IF( kt == nit000 ) THEN
-         IF( .NOT. l_istiled .OR. ntile == 1 ) THEN                    ! Do only on the first tile
+         IF( ntile == 0 .OR. ntile == 1 ) THEN                    ! Do only on the first tile
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes'
@@ -2948,9 +2873,7 @@
       !
       IF( kt == nit000 ) THEN
-         IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity'
-            IF(lwp) WRITE(numout,*) '~~~~~~~   '
-         ENDIF
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity'
+         IF(lwp) WRITE(numout,*) '~~~~~~~   '
       ENDIF
       !code saving tracer trends removed, replace with trdmxl_oce
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 14711)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/step.F90	(revision 14712)
@@ -224,4 +224,11 @@
                &         CALL Agrif_Sponge_dyn        ! momentum sponge
 #endif
+      IF (nn_hls.eq.2) THEN
+         ! [comm_cleanup] ! needed from DYN 
+         CALL lbc_lnk( 'stp', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1., ww(:,:,:), 'W', 1.)
+         IF(.NOT.ln_linssh) CALL lbc_lnk( 'stp', e3f, 'F', 1. )
+         ! [comm_cleanup] ! needed from DYN dyn_ldf_blp 
+         CALL lbc_lnk( 'stp', uu(:,:,:,Nbb), 'U', -1., vv(:,:,:,Nbb), 'V', -1.)
+      ENDIF
                          CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)	==> RHS
                          CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity           	==> RHS
@@ -349,5 +356,7 @@
 !!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 )
+                         IF (nn_hls.eq.2) CALL lbc_lnk( 'stp', ts(:,:,:,jp_tem,Naa), 'T', 1.0_wp, ts(:,:,:,jp_sal,Naa), 'T', 1.0_wp, &
+                                                  &            uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.)
+
                          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 14711)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/stpmlf.F90	(revision 14712)
@@ -62,5 +62,4 @@
 #  include "do_loop_substitute.h90"
 #  include "domzgr_substitute.h90"
-#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -241,8 +240,17 @@
                &         CALL Agrif_Sponge_dyn        ! momentum sponge
 #endif
+      ! [comm_cleanup] ! lbc_lnk from DYN
+      IF (nn_hls.eq.2) THEN
+         CALL lbc_lnk( 'stp_MLF', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1., ww(:,:,:), 'W', 1.,       & 
+                 &                uu(:,:,:,Nbb), 'U', -1., vv(:,:,:,Nbb), 'V', -1.)
+         IF(.NOT.lk_linssh) CALL lbc_lnk( 'stp_MLF', r3u(:,:,Nnn), 'U', 1.0_wp, r3v(:,:,Nnn), 'V', 1.0_wp,    &
+                 &                                   r3u(:,:,Nbb), 'U', 1.0_wp, r3v(:,:,Nbb), 'V', 1.0_wp,    &
+                 &                                   r3t(:,:,Nbb), 'T', 1.0_wp )
+      ENDIF
                          CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)	==> RHS
                          CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity           	==> RHS
                          CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing
       IF( ln_zdfosm  )   CALL dyn_osm( kstp,      Nnn      , uu, vv, Nrhs )  ! OSMOSIS non-local velocity fluxes ==> RHS
+
                          CALL dyn_hpg( kstp,      Nnn      , uu, vv, Nrhs )  ! horizontal gradient of Hydrostatic pressure
                          CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa )  ! surface pressure gradient
