Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfiwm.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfiwm.F90	(revision 14315)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfiwm.F90	(revision 14636)
@@ -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)
