Index: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfgls.F90
===================================================================
--- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfgls.F90	(revision 14315)
+++ NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfgls.F90	(revision 14636)
@@ -179,5 +179,6 @@
 
       ! Compute surface, top and bottom friction at T-points
-      DO_2D( 0, 0, 0, 0 )          !==  surface ocean friction
+      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )          !==  surface ocean friction
+      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )          !==  surface ocean friction 
          ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1)   ! surface friction
       END_2D
@@ -186,5 +187,6 @@
       !
       IF( .NOT.ln_drg_OFF ) THEN     !== top/bottom friction   (explicit before friction)
-         DO_2D( 0, 0, 0, 0 )         ! bottom friction (explicit before friction)
+         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )         ! bottom friction (explicit before friction)
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )          ! bottom friction (explicit before friction)
             zmsku = 0.5_wp * ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
             zmskv = 0.5_wp * ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0)
@@ -193,5 +195,6 @@
          END_2D
          IF( ln_isfcav ) THEN
-            DO_2D( 0, 0, 0, 0 )      ! top friction
+            ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )      ! top friction
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )      ! top friction
                zmsku = 0.5_wp * ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
                zmskv = 0.5_wp * ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )     ! (CAUTION: CdU<0)
@@ -220,5 +223,6 @@
       zhsro(:,:) = ( (1._wp-zice_fra(:,:)) * zhsro(:,:) + zice_fra(:,:) * rn_hsri )*tmask(:,:,1)  + (1._wp - tmask(:,:,1))*rn_hsro
       !
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )  !==  Compute dissipation rate  ==!
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )  !==  Compute dissipation rate  ==!
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )  !==  Compute dissipation rate  ==!
          eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk)
       END_3D
@@ -229,5 +233,6 @@
 
       IF( nn_clos == 0 ) THEN    ! Mellor-Yamada
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             zup   = hmxl_n(ji,jj,jk) * gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)
             zdown = vkarmn * gdepw(ji,jj,jk,Kmm) * ( -gdepw(ji,jj,jk,Kmm) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) )
@@ -250,5 +255,6 @@
       ! Warning : after this step, en : right hand side of the matrix
 
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
          !
          buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk)     ! stratif. destruction
@@ -327,5 +333,6 @@
       ! at k=2, set de/dz=Fw
       !cbr
-      DO_2D( 0, 0, 0, 0 )   ! zdiag zd_lw not defined/used on the halo
+      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )   ! zdiag zd_lw not defined/used on the halo
+      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )   ! zdiag zd_lw not defined/used on the halo
          zdiag(ji,jj,2) = zdiag(ji,jj,2) +  zd_lw(ji,jj,2) ! Remove zd_lw from zdiag
          zd_lw(ji,jj,2) = 0._wp
@@ -348,5 +355,6 @@
          !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin
          !                      ! Balance between the production and the dissipation terms
-         DO_2D( 0, 0, 0, 0 )
+         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
 !!gm This means that bottom and ocean w-level above have a specified "en" value.   Sure ????
 !!   With thick deep ocean level thickness, this may be quite large, no ???
@@ -366,5 +374,6 @@
          !
          IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
-            DO_2D( 0, 0, 0, 0 )
+            ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                itop   = mikt(ji,jj)       ! k   top w-point
                itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
@@ -384,5 +393,6 @@
       CASE ( 1 )             ! Neumman boundary condition
          !
-         DO_2D( 0, 0, 0, 0 )
+         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
             ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -399,5 +409,6 @@
          END_2D
          IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
-            DO_2D( 0, 0, 0, 0 )
+            ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                itop   = mikt(ji,jj)       ! k   top w-point
                itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
@@ -420,11 +431,14 @@
       ! ----------------------------------------------------------
       !
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
          zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
       END_3D
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
          zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1)
       END_3D
-      DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )           ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+      ! [comm_cleanup] ! DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+      DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
          en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
       END_3D
@@ -441,20 +455,24 @@
       !
       CASE( 0 )               ! k-kl  (Mellor-Yamada)
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             psi(ji,jj,jk)  = eb(ji,jj,jk) * hmxl_b(ji,jj,jk)
          END_3D
          !
       CASE( 1 )               ! k-eps
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             psi(ji,jj,jk)  = eps(ji,jj,jk)
          END_3D
          !
       CASE( 2 )               ! k-w
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) )
          END_3D
          !
       CASE( 3 )               ! generic
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn
          END_3D
@@ -469,5 +487,6 @@
       ! Warning : after this step, en : right hand side of the matrix
 
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
          !
          ! psi / k
@@ -541,5 +560,6 @@
          !
          ! Neumann condition at k=2
-         DO_2D( 0, 0, 0, 0 )   ! zdiag zd_lw not defined/used on the halo
+         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )   ! zdiag zd_lw not defined/used on the halo
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )   ! zdiag zd_lw not defined/used on the halo
             zdiag(ji,jj,2) = zdiag(ji,jj,2) +  zd_lw(ji,jj,2) ! Remove zd_lw from zdiag
             zd_lw(ji,jj,2) = 0._wp
@@ -569,5 +589,6 @@
          !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot
          !                      ! Balance between the production and the dissipation terms
-         DO_2D( 0, 0, 0, 0 )
+         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
             ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -588,5 +609,6 @@
       CASE ( 1 )             ! Neumman boundary condition
          !
-         DO_2D( 0, 0, 0, 0 )
+         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
             ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -616,11 +638,14 @@
       ! ----------------
       !
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
          zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
       END_3D
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
          zd_lw(ji,jj,jk) = psi(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1)
       END_3D
-      DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+      ! [comm_cleanup] ! DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+      DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
          psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
       END_3D
@@ -632,15 +657,18 @@
       !
       CASE( 0 )               ! k-kl  (Mellor-Yamada)
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
             eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / MAX( psi(ji,jj,jk), rn_epsmin)
          END_3D
          !
       CASE( 1 )               ! k-eps
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
             eps(ji,jj,jk) = psi(ji,jj,jk)
          END_3D
          !
       CASE( 2 )               ! k-w
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
             eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk)
          END_3D
@@ -650,5 +678,6 @@
          zex1  =      ( 1.5_wp + rmm/rnn )
          zex2  = -1._wp / rnn
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
             eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2
          END_3D
@@ -658,5 +687,6 @@
       ! Limit dissipation rate under stable stratification
       ! --------------------------------------------------
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! Note that this set boundary conditions on hmxl_n at the same time
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! Note that this set boundary conditions on hmxl_n at the same time
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )   ! Note that this set boundary conditions on hmxl_n at the same time
          ! limitation
          eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
@@ -674,5 +704,6 @@
       !
       CASE ( 0 , 1 )             ! Galperin or Kantha-Clayson stability functions
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             ! zcof =  l²/q²
             zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
@@ -691,5 +722,6 @@
          !
       CASE ( 2, 3 )               ! Canuto stability functions
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             ! zcof =  l²/q²
             zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )
@@ -723,5 +755,6 @@
       ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0)
       zstm(:,:,jpk) = 0.
-      DO_2D( 0, 0, 0, 0 )             ! update bottom with good values
+      ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )             ! update bottom with good values
+      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )             ! update bottom with good values
          zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj))
       END_2D
@@ -738,5 +771,6 @@
       !     later overwritten by surface/bottom boundaries conditions, so we don't really care of p_avm(:,:1) and p_avm(:,:jpk)
       !     for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0)
-      DO_3D( 0, 0, 0, 0, 1, jpk )
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpk )
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
          zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk)
          zavt  = zsqen * zstt(ji,jj,jk)
