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 14574)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfosm.F90	(revision 14636)
@@ -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  )  &
