Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DOM/domqco.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DOM/domqco.F90	(revision 14679)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DOM/domqco.F90	(revision 14680)
@@ -154,5 +154,6 @@
 #if ! defined key_qcoTest_FluxForm
       !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
-         DO_2D( 0, 0, 0, 0 )
+         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
+         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
             pr3u(ji,jj) = 0.5_wp * (  e1e2t(ji  ,jj) * pssh(ji  ,jj)  &
                &                    + e1e2t(ji+1,jj) * pssh(ji+1,jj)  ) * r1_hu_0(ji,jj) * r1_e1e2u(ji,jj)
@@ -162,5 +163,6 @@
 !!st      ELSE                                         !- Flux Form   (simple averaging)
 #else
-         DO_2D( 0, 0, 0, 0 )
+         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
+         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
             pr3u(ji,jj) = 0.5_wp * (  pssh(ji,jj) + pssh(ji+1,jj  )  ) * r1_hu_0(ji,jj)
             pr3v(ji,jj) = 0.5_wp * (  pssh(ji,jj) + pssh(ji  ,jj+1)  ) * r1_hv_0(ji,jj)
@@ -170,5 +172,5 @@
       !
       IF( .NOT.PRESENT( pr3f ) ) THEN              !- lbc on ratio at u-, v-points only
-         CALL lbc_lnk( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp )
+         IF (nn_hls.eq.1) CALL lbc_lnk( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp )
          !
          !
@@ -179,20 +181,31 @@
          !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
 
-            DO_2D( 0, 0, 0, 0 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line
-               pr3f(ji,jj) = 0.25_wp * (  e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )  &
-                  &                     + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )  &
-                  &                     + e1e2t(ji  ,jj+1) * pssh(ji  ,jj+1)  &
-                  &                     + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj)
+            ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line
+               ! round brackets added to fix the order of floating point operations
+               ! needed to ensure halo 1 - halo 2 compatibility
+               pr3f(ji,jj) = 0.25_wp * ( ( e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )   &
+                  &                      + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )   &
+                  &                      )                                      & ! bracket for halo 1 - halo 2 compatibility
+                  &                     + ( e1e2t(ji  ,jj+1) * pssh(ji  ,jj+1)  &
+                  &                       + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  &
+                  &                       )                                     & ! bracket for halo 1 - halo 2 compatibility
+                  &                    ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj)
             END_2D
 !!st         ELSE                                      !- Flux Form   (simple averaging)
 #else
-            DO_2D( 0, 0, 0, 0 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line
-               pr3f(ji,jj) = 0.25_wp * (  pssh(ji,jj  ) + pssh(ji+1,jj  )  &
-                  &                     + pssh(ji,jj+1) + pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj)
+            ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
+               ! round brackets added to fix the order of floating point operations
+               ! needed to ensure halo 1 - halo 2 compatibility
+               pr3f(ji,jj) = 0.25_wp * ( ( pssh(ji,jj  ) + pssh(ji+1,jj  ) ) &
+                  &                     + ( pssh(ji,jj+1) + pssh(ji+1,jj+1)  & 
+                  &                       )                                  & ! bracket for halo 1 - halo 2 compatibility
+                  &                    ) * r1_hf_0(ji,jj)
             END_2D
 !!st         ENDIF
 #endif
          !                                                 ! lbc on ratio at u-,v-,f-points
-         CALL lbc_lnk( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp, pr3f, 'F', 1._wp )
+         IF (nn_hls.eq.1) CALL lbc_lnk( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp, pr3f, 'F', 1._wp )
          !
       ENDIF
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/divhor.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/divhor.F90	(revision 14679)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/divhor.F90	(revision 14680)
@@ -64,6 +64,4 @@
       !
       INTEGER  ::   ji, jj, jk    ! dummy loop indices
-      REAL(wp) ::   zraur, zdep   ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj) :: ztmp
       !!----------------------------------------------------------------------
       !
@@ -77,10 +75,15 @@
       ENDIF
       !
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )                                    !==  Horizontal divergence  ==!
-         hdiv(ji,jj,jk) = (   e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * uu(ji  ,jj,jk,Kmm)      &
-            &               - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm)      &
-            &               + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * vv(ji,jj  ,jk,Kmm)      &
-            &               - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm)  )   &
-            &            * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )                                    !==  Horizontal divergence  ==!
+      DO_3D( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 )                                    !==  Horizontal divergence  ==!
+         ! round brackets added to fix the order of floating point operations
+         ! needed to ensure halo 1 - halo 2 compatibility
+         hdiv(ji,jj,jk) = (  ( e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * uu(ji  ,jj,jk,Kmm)     &
+            &                - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm)     &
+            &                )                                                             & ! bracket for halo 1 - halo 2 compatibility
+            &              + ( e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * vv(ji,jj  ,jk,Kmm)     &
+            &                - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm)     & 
+            &                )                                                             & ! bracket for halo 1 - halo 2 compatibility
+            &             )  * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
       END_3D
       !
@@ -91,8 +94,8 @@
       ! 
 #endif
-      !
+      ! WED025 + isomip true 
       IF( ln_isf )                      CALL isf_hdiv( kt, Kmm, hdiv )           !==  ice shelf         ==!   (update hdiv field)
       !
-      CALL lbc_lnk( 'divhor', hdiv, 'T', 1.0_wp )   !   (no sign change)
+      IF (nn_hls.eq.1) CALL lbc_lnk( 'divhor', hdiv, 'T', 1.0_wp )   !   (no sign change)
       !
       IF( ln_timing )   CALL timing_stop('div_hor')
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 14679)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynadv.F90	(revision 14680)
@@ -76,6 +76,8 @@
          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 14679)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynadv_cen2.F90	(revision 14680)
@@ -72,5 +72,6 @@
          zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm)
          zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
-         DO_2D( 1, 0, 1, 0 )              ! horizontal momentum fluxes (at T- and F-point)
+         ! [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)
             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) )
@@ -78,5 +79,6 @@
             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
-         DO_2D( 0, 0, 0, 0 )              ! divergence of horizontal momentum fluxes
+         ! [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
             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)   &
@@ -98,10 +100,12 @@
       !                             !==  Vertical advection  ==!
       !
-      DO_2D( 0, 0, 0, 0 )                 ! surface/bottom advective fluxes set to zero
+      ! [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
          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
-         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 )
             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)
@@ -109,13 +113,16 @@
       ENDIF
       DO jk = 2, jpkm1                    ! interior advective fluxes
-         DO_2D( 0, 1, 0, 1 )                  ! 1/4 * Vertical transport
+         ! [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 
             zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
          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 )
             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
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )       ! divergence of vertical momentum flux divergence
+      ! [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
          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 14679)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynadv_ubs.F90	(revision 14680)
@@ -108,5 +108,6 @@
          zfv(:,:,jk) = e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
          !            
-         DO_2D( 0, 0, 0, 0 )                       ! laplacian
+         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )                       ! laplacian
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                       ! laplacian
             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)
@@ -124,8 +125,8 @@
          END_2D
       END DO
-      CALL lbc_lnk( 'dynadv_ubs', zlu_uu(:,:,:,1), 'U', 1.0_wp , zlu_uv(:,:,:,1), 'U', 1.0_wp,  &
-         &                        zlu_uu(:,:,:,2), 'U', 1.0_wp , zlu_uv(:,:,:,2), 'U', 1.0_wp,  & 
-         &                        zlv_vv(:,:,:,1), 'V', 1.0_wp , zlv_vu(:,:,:,1), 'V', 1.0_wp,  &
-         &                        zlv_vv(:,:,:,2), 'V', 1.0_wp , zlv_vu(:,:,:,2), 'V', 1.0_wp   )
+      IF (nn_hls.eq.1) CALL lbc_lnk( 'dynadv_ubs', zlu_uu(:,:,:,1), 'U', 1.0_wp , zlu_uv(:,:,:,1), 'U', 1.0_wp,  &
+                          &                        zlu_uu(:,:,:,2), 'U', 1.0_wp , zlu_uv(:,:,:,2), 'U', 1.0_wp,  & 
+                          &                        zlv_vv(:,:,:,1), 'V', 1.0_wp , zlv_vu(:,:,:,1), 'V', 1.0_wp,  &
+                          &                        zlv_vv(:,:,:,2), 'V', 1.0_wp , zlv_vu(:,:,:,2), 'V', 1.0_wp   )
       !
       !                                      ! ====================== !
@@ -136,5 +137,6 @@
          zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
          !
-         DO_2D( 1, 0, 1, 0 )                       ! horizontal momentum fluxes at T- and F-point
+         ! [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
             zui = ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj  ,jk,Kmm) )
             zvj = ( pvv(ji,jj,jk,Kmm) + pvv(ji  ,jj+1,jk,Kmm) )
@@ -168,5 +170,6 @@
                &                * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) - gamma1 * zl_v )
          END_2D
-         DO_2D( 0, 0, 0, 0 )                       ! divergence of horizontal momentum fluxes
+         ! [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
             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)   &
@@ -187,5 +190,6 @@
       !                                      !  Vertical advection  !
       !                                      ! ==================== !
-      DO_2D( 0, 0, 0, 0 )                          ! surface/bottom advective fluxes set to zero
+      ! [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
          zfu_uw(ji,jj,jpk) = 0._wp
          zfv_vw(ji,jj,jpk) = 0._wp
@@ -194,5 +198,6 @@
       END_2D
       IF( ln_linssh ) THEN                         ! constant volume : advection through the surface
-         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 )
             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)
@@ -200,13 +205,16 @@
       ENDIF
       DO jk = 2, jpkm1                          ! interior fluxes
-         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 ) 
             zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
          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 )
             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
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )             ! divergence of vertical momentum flux divergence
+      ! [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
          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/dynkeg.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynkeg.F90	(revision 14679)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynkeg.F90	(revision 14680)
@@ -101,5 +101,6 @@
       !
       CASE ( nkeg_C2 )                          !--  Standard scheme  --!
-         DO_3D( 0, 1, 0, 1, 1, jpkm1 )
+         ! [comm_cleanup] ! DO_3D( 0, 1, 0, 1, 1, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 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)
@@ -109,5 +110,6 @@
          END_3D
       CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         ! [comm_cleanup: Hollingsworth scheme NOT TESTED] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
             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) )  &
@@ -121,9 +123,10 @@
             zhke(ji,jj,jk) = r1_48 * ( zv + zu )
          END_3D
-         CALL lbc_lnk( 'dynkeg', zhke, 'T', 1.0_wp )
+         IF (nn_hls.eq.1) CALL lbc_lnk( 'dynkeg', zhke, 'T', 1.0_wp )
          !
       END SELECT 
       !
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )       !==  grad( KE ) added to the general momentum trends  ==!
+      ! [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  ==!
          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/dynzad.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynzad.F90	(revision 14679)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynzad.F90	(revision 14680)
@@ -79,5 +79,6 @@
 
       DO jk = 2, jpkm1                ! Vertical momentum advection at level w and u- and v- vertical
-         DO_2D( 0, 1, 0, 1 )              ! vertical fluxes
+         ! [comm_cleanup] ! DO_2D( 0, 1, 0, 1 )              ! vertical fluxes 
+         DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )              ! vertical fluxes
           IF( ln_vortex_force ) THEN
             zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wsd(ji,jj,jk) )
@@ -86,5 +87,6 @@
           ENDIF
          END_2D
-         DO_2D( 0, 0, 0, 0 )              ! vertical momentum advection at w-point
+         ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )              ! vertical momentum advection at w-point
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )              ! vertical momentum advection at w-point
             zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( puu(ji,jj,jk-1,Kmm) - puu(ji,jj,jk,Kmm) )
             zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( pvv(ji,jj,jk-1,Kmm) - pvv(ji,jj,jk,Kmm) )
@@ -93,5 +95,6 @@
       !
       ! Surface and bottom advective fluxes set to zero
-      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 )
          zwuw(ji,jj, 1 ) = 0._wp
          zwvw(ji,jj, 1 ) = 0._wp
@@ -100,5 +103,6 @@
       END_2D
       !
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! Vertical momentum advection at u- and v-points
+      ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! Vertical momentum advection at u- and v-points
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )   ! Vertical momentum advection at u- and v-points
          puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zwuw(ji,jj,jk) + zwuw(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/sshwzv.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/sshwzv.F90	(revision 14679)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/sshwzv.F90	(revision 14680)
@@ -78,5 +78,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! sea-surface height
       ! 
-      INTEGER  ::   jk      ! dummy loop index
+      INTEGER  ::   ji, jj, jk      ! dummy loop index
       REAL(wp) ::   zcoef   ! local scalar
       REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv   ! 2D workspace
@@ -103,12 +103,16 @@
       !
       zhdiv(:,:) = 0._wp
-      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
-        zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
-      END DO
+      ! [comm_cleanup] ! DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports 
+      DO_3D( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 )                                 ! Horizontal divergence of barotropic transports 
+        zhdiv(ji,jj) = zhdiv(ji,jj) + e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk)
+      END_3D
       !                                                ! Sea surface elevation time stepping
       ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
       ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
       ! 
-      pssh(:,:,Kaa) = (  pssh(:,:,Kbb) - rDt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
+      ! [comm_cleanup] 
+      DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
+         pssh(ji,jj,Kaa) = (  pssh(ji,jj,Kbb) - rDt * ( zcoef * ( emp_b(ji,jj) + emp(ji,jj) ) + zhdiv(ji,jj) )  ) * ssmask(ji,jj)
+      END_2D
       !
 #if defined key_agrif
@@ -119,5 +123,6 @@
       IF ( .NOT.ln_dynspg_ts ) THEN
          IF( ln_bdy ) THEN
-            CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp )    ! Not sure that's necessary
+            ! [comm_cleanup] 
+            IF (nn_hls.eq.1) CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp )    ! Not sure that's necessary
             CALL bdy_ssh( pssh(:,:,Kaa) )             ! Duplicate sea level across open boundaries
          ENDIF
@@ -178,9 +183,10 @@
             ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
             ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
-            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 )
                zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
             END_2D
          END DO
-         CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
+         IF (nn_hls.eq.1) CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
          !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
          !                             ! Same question holds for hdiv. Perhaps just for security
@@ -357,5 +363,6 @@
       zdt = 2._wp * rn_Dt                            ! 2*rn_Dt and not rDt (for restartability)
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
-         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 )
             z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
             Cu_adv(ji,jj,jk) =   zdt *                                                         &
@@ -374,5 +381,6 @@
          END_3D
       ELSE
-         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 )
             z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
             Cu_adv(ji,jj,jk) =   zdt *                                                      &
@@ -387,10 +395,11 @@
          END_3D
       ENDIF
-      CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1.0_wp )
+      IF (nn_hls.eq.1) CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1.0_wp )
       !
       CALL iom_put("Courant",Cu_adv)
       !
       IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
-         DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )             ! or scan Courant criterion and partition ! w where necessary
+         ! [comm_cleanup] ! DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )             ! or scan Courant criterion and partition ! w where necessary
+         DO_3DS( nn_hls, nn_hls, nn_hls, nn_hls, jpkm1, 2, -1 )             ! or scan Courant criterion and partition ! w where necessary
             !
             zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ISF/isfhdiv.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ISF/isfhdiv.F90	(revision 14679)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ISF/isfhdiv.F90	(revision 14680)
@@ -100,5 +100,6 @@
       !
       ! update divergence at each level affected by ice shelf top boundary layer
-      DO_2D( 1, 1, 1, 1 )
+      ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
          ikt = ktop(ji,jj)
          ikb = kbot(ji,jj)
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/LDF/ldfc1d_c2d.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/LDF/ldfc1d_c2d.F90	(revision 14679)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/LDF/ldfc1d_c2d.F90	(revision 14680)
@@ -135,5 +135,5 @@
       !
       CASE( 'DYN' )                       ! T- and F-points
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             pah1(ji,jj,1) = pUfac * MAX( e1t(ji,jj) , e2t(ji,jj) )**knn
             pah2(ji,jj,1) = pUfac * MAX( e1f(ji,jj) , e2f(ji,jj) )**knn
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/LDF/ldftra.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/LDF/ldftra.F90	(revision 14679)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/LDF/ldftra.F90	(revision 14680)
@@ -633,9 +633,9 @@
       INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step index
       INTEGER                         , INTENT(in   ) ::   Kmm            ! ocean time level indices
-      REAL(wp)                        , INTENT(inout) ::   paei0          ! max value            [m2/s]
+      REAL(wp)                        , INTENT(in   ) ::   paei0          ! max value            [m2/s]
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   paeiu, paeiv   ! eiv coefficient      [m2/s]
       !
       INTEGER  ::   ji, jj, jk    ! dummy loop indices
-      REAL(wp) ::   zfw, ze3w, zn2, z1_f20, zaht, zaht_min, zzaei    ! local scalars
+      REAL(wp) ::   zfw, ze3w, zn2, z1_f20, zzaei    ! local scalars
       REAL(wp), DIMENSION(jpi,jpj) ::   zn, zah, zhw, zRo, zaeiw   ! 2D workspace
       !!----------------------------------------------------------------------
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/SBC/sbcrnf.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/SBC/sbcrnf.F90	(revision 14679)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/SBC/sbcrnf.F90	(revision 14680)
@@ -206,5 +206,6 @@
       IF( ln_rnf_depth .OR. ln_rnf_depth_ini ) THEN      !==   runoff distributed over several levels   ==!
          IF( ln_linssh ) THEN    !* constant volume case : just apply the runoff input flow
-            DO_2D( 1, 1, 1, 1 )
+            ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
+            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
                DO jk = 1, nk_rnf(ji,jj)
                   phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rho0 / h_rnf(ji,jj)
@@ -212,5 +213,6 @@
             END_2D
          ELSE                    !* variable volume case
-            DO_2D( 1, 1, 1, 1 )              ! update the depth over which runoffs are distributed
+            ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )              ! update the depth over which runoffs are distributed
+            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )              ! update the depth over which runoffs are distributed
                h_rnf(ji,jj) = 0._wp
                DO jk = 1, nk_rnf(ji,jj)                             ! recalculates h_rnf to be the depth in metres
@@ -358,5 +360,6 @@
          !
          nk_rnf(:,:) = 0                               ! set the number of level over which river runoffs are applied
-         DO_2D( 1, 1, 1, 1 )
+         ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             IF( h_rnf(ji,jj) > 0._wp ) THEN
                jk = 2
@@ -371,5 +374,6 @@
             ENDIF
          END_2D
-         DO_2D( 1, 1, 1, 1 )                           ! set the associated depth
+         ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )                           ! set the associated depth
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )                           ! set the associated depth
             h_rnf(ji,jj) = 0._wp
             DO jk = 1, nk_rnf(ji,jj)
@@ -401,5 +405,6 @@
          WHERE( zrnfcl(:,:,1) > 0._wp )  h_rnf(:,:) = zacoef * zrnfcl(:,:,1)   ! compute depth for all runoffs
          !
-         DO_2D( 1, 1, 1, 1 )                ! take in account min depth of ocean rn_hmin
+         ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )                ! take in account min depth of ocean rn_hmin
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )                ! take in account min depth of ocean rn_hmin
             IF( zrnfcl(ji,jj,1) > 0._wp ) THEN
                jk = mbkt(ji,jj)
@@ -409,5 +414,6 @@
          !
          nk_rnf(:,:) = 0                       ! number of levels on which runoffs are distributed
-         DO_2D( 1, 1, 1, 1 )
+         ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             IF( zrnfcl(ji,jj,1) > 0._wp ) THEN
                jk = 2
@@ -420,5 +426,6 @@
          END_2D
          !
-         DO_2D( 1, 1, 1, 1 )                          ! set the associated depth
+         ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )                          ! set the associated depth
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )                          ! set the associated depth
             h_rnf(ji,jj) = 0._wp
             DO jk = 1, nk_rnf(ji,jj)
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf.F90	(revision 14679)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf.F90	(revision 14680)
@@ -92,6 +92,4 @@
             CALL tra_ldf_triad( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts,  1 )
          CASE ( np_blp , np_blp_i , np_blp_it )             ! bilaplacian: iso-level & iso-neutral operators
-            ! [comm_cleanup]
-            ! IF (nn_hls.EQ.2) CALL lbc_lnk( 'tra_ldf', pts(:,:,:,:,Kbb), 'T',1.)
             CALL tra_ldf_blp  ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs),             jpts, nldf_tra )
          END SELECT
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_iso.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_iso.F90	(revision 14679)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_iso.F90	(revision 14680)
@@ -185,11 +185,14 @@
             zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          &
                &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  )
-               !
-            ! NOTE: [halo1-halo2]
-            ! Extra brackets required to ensure consistent floating point arithmetic for different nn_hls for bilaplacian
-            zahu_w = (   (pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk))    &
-               &       + (pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk))  ) * zmsku
-            zahv_w = (   (pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk))    &
-               &       + (pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk))  ) * zmskv
+            ! round brackets added to fix the order of floating point operations
+            ! needed to ensure halo 1 - halo 2 compatibility
+            zahu_w = ( (  pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)                    &
+               &       )                                                           & ! bracket for halo 1 - halo 2 compatibility
+               &       + ( pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)                   &
+               &         ) ) * zmsku                                                 ! bracket for halo 1 - halo 2 compatibility
+            zahv_w = ( (  pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)                    &
+               &       )                                                           & ! bracket for halo 1 - halo 2 compatibility
+               &       + ( pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)                   &
+               &         ) ) * zmskv                                                 ! bracket for halo 1 - halo 2 compatibility
                !
             ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
@@ -200,11 +203,13 @@
             ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
-               ! NOTE: [halo1-halo2]
-               ! Extra brackets required to ensure consistent floating point arithmetic for different nn_hls for bilaplacian
-               akz(ji,jj,jk) = 0.25_wp * (                                                                       &
-                  &              (( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )    &
-                  &            + (  pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) ))   &
-                  &            + (( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) )    &
-                  &            + (  pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ))   )
+               ! round brackets added to fix the order of floating point operations
+               ! needed to ensure halo 1 - halo 2 compatibility
+               akz(ji,jj,jk) = 0.25_wp * (                                                                     &
+                  &            ( ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   &
+                  &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   &
+                  &            )                                                                               & ! bracket for halo 1 - halo 2 compatibility
+                  &            + ( ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) ) &
+                  &              + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) &
+                  &              ) )                                                                             ! bracket for halo 1 - halo 2 compatibility
             END_3D
             !
@@ -296,21 +301,30 @@
                zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
                !
-               ! NOTE: [halo1-halo2]
-               ! Extra brackets required to ensure consistent floating point arithmetic for different nn_hls for bilaplacian
-               zftu(ji,jj,jk) = (  zabe1 * zdit(ji,jj,jk)                          &
-                  &              + zcof1 * (  (zdkt (ji+1,jj) + zdk1t(ji,jj))      &
-                  &                         + (zdk1t(ji+1,jj) + zdkt (ji,jj))  )  ) * umask(ji,jj,jk)
-               zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)                          &
-                  &              + zcof2 * (  (zdkt (ji,jj+1) + zdk1t(ji,jj))      &
-                  &                         + (zdk1t(ji,jj+1) + zdkt (ji,jj))  )  ) * vmask(ji,jj,jk)
+               ! round brackets added to fix the order of floating point operations
+               ! needed to ensure halo 1 - halo 2 compatibility
+               zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)                       &
+                  &               + zcof1 * ( ( zdkt (ji+1,jj) + zdk1t(ji,jj)    &
+                  &                           )                                  & ! bracket for halo 1 - halo 2 compatibility
+                  &                         + ( zdk1t(ji+1,jj) + zdkt (ji,jj)    &
+                  &                           )                                  & ! bracket for halo 1 - halo 2 compatibility
+                  &                         ) ) * umask(ji,jj,jk)
+               zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)                        &
+                  &              + zcof2 * ( ( zdkt (ji,jj+1) + zdk1t(ji,jj)     &
+                  &                           )                                  & ! bracket for halo 1 - halo 2 compatibility
+                  &                         + ( zdk1t(ji,jj+1) + zdkt (ji,jj)    &
+                  &                           )                                  & ! bracket for halo 1 - halo 2 compatibility
+                  &                         ) ) * vmask(ji,jj,jk)
             END_2D
             !
             ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )           !== horizontal divergence and add to pta
             DO_2D( iij-1, iij-1, iij-1, iij-1 )           !== horizontal divergence and add to pta
-               ! NOTE: [halo1-halo2]
-               ! Extra brackets required to ensure consistent floating point arithmetic for different nn_hls for bilaplacian
-               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    &
-                  &       + zsign * (  (zftu(ji,jj,jk) - zftu(ji-1,jj,jk)) + &
-                  &                    (zftv(ji,jj,jk) - zftv(ji,jj-1,jk))  ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+               ! round brackets added to fix the order of floating point operations
+               ! needed to ensure halo 1 - halo 2 compatibility
+               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)                         &
+                  &       + zsign * ( ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk)        &
+                  &                   )                                          & ! bracket for halo 1 - halo 2 compatibility
+                  &                 + ( zftv(ji,jj,jk) - zftv(ji,jj-1,jk)        &
+                  &                   )                                          & ! bracket for halo 1 - halo 2 compatibility
+                  &                 ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
             END_2D
          END DO                                        !   End of slab
@@ -341,10 +355,16 @@
             zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk)
             !
-            ! NOTE: [halo1-halo2]
-            ! Extra brackets required to ensure consistent floating point arithmetic for different nn_hls for bilaplacian
-            ztfw(ji,jj,jk) = zcoef3 * (   (zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk))      &
-               &                        + (zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk))  )   &
-               &           + zcoef4 * (   (zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk))      &
-               &                        + (zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk))  )
+            ! round brackets added to fix the order of floating point operations
+            ! needed to ensure halo 1 - halo 2 compatibility
+            ztfw(ji,jj,jk) = zcoef3 * ( ( zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)    &
+                  &                     )                                              & ! bracket for halo 1 - halo 2 compatibility
+                  &                   + ( zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)    &
+                  &                     )                                              & ! bracket for halo 1 - halo 2 compatibility
+                  &                   )                                                &
+                  &        + zcoef4 * ( ( zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)    &
+                  &                     )                                              & ! bracket for halo 1 - halo 2 compatibility
+                  &                   + ( zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)    &
+                  &                     )                                              & ! bracket for halo 1 - halo 2 compatibility
+                  &                   )
          END_3D
          !                                !==  add the vertical 33 flux  ==!
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_lap_blp.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_lap_blp.F90	(revision 14679)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_lap_blp.F90	(revision 14680)
@@ -158,10 +158,12 @@
          ENDIF
          !
-         DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 )        !== Second derivative (divergence) added to the general tracer trends  ==!
-            ! NOTE: [halo1-halo2]
-            ! Extra brackets required to ensure consistent floating point arithmetic for different nn_hls for bilaplacian
-            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + (  (ztu(ji,jj,jk) - ztu(ji-1,jj,jk))     &
-               &                                      +    (ztv(ji,jj,jk) - ztv(ji,jj-1,jk)) )   &
-               &                                      / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
+         DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 )            !== Second derivative (divergence) added to the general tracer trends  ==!
+            ! round brackets added to fix the order of floating point operations
+            ! needed to ensure halo 1 - halo 2 compatibility
+            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + ( ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk)    &
+               &                                          )                                    & ! bracket for halo 1 - halo 2 compatibility
+               &                                      +   ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk)    &
+               &                                          )                                    & ! bracket for halo 1 - halo 2 compatibility
+               &                                        ) / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
          END_3D
          !
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_triad.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_triad.F90	(revision 14679)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_triad.F90	(revision 14680)
@@ -107,20 +107,13 @@
       REAL(wp), DIMENSION(A2D_T(ktt_rhs),JPK,KJPT), INTENT(inout) ::   pt_rhs     ! tracer trend
       !
-      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
-      INTEGER  ::  ip,jp,kp         ! dummy loop indices
-      INTEGER  ::  ierr, iij        ! local integer
-      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3    ! local scalars
-      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4    !   -      -
+      INTEGER  ::  ji, jj, jk, jn, kp, iij   ! dummy loop indices
       REAL(wp) ::  zcoef0, ze3w_2, zsign          !   -      -
       !
-      REAL(wp) ::   zslope_skew, zslope_iso, zslope2, zbu, zbv
-      REAL(wp) ::   ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt
-      REAL(wp) ::   zah, zah_slp, zaei_slp
-      REAL(wp), DIMENSION(A2D(nn_hls),0:1) ::   zdkt3d                                           ! vertical tracer gradient at 2 levels
-      REAL(wp), DIMENSION(A2D(nn_hls)    ) ::   z2d                                              ! 2D workspace
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zdit, zdjt, zpsi_uw, zpsi_vw   ! 3D     -
-      ! NOTE: [halo1-halo2]
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk,2) :: zftu, zftv
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk,2,2) :: ztfw
+      REAL(wp) ::   zslope2, zbu, zbv, zbu1, zbv1, zslope21, zah, zah1, zah_ip1, zah_jp1, zbu_ip1, zbv_jp1
+      REAL(wp) ::   ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt, zdyt_jp1, ze3wr_jp1, zdzt_jp1, zah_slp1, zah_slp_jp1, zaei_slp_jp1
+      REAL(wp) ::   zah_slp, zaei_slp, zdxt_ip1, ze3wr_ip1, zdzt_ip1, zah_slp_ip1, zaei_slp_ip1, zaei_slp1
+      REAL(wp), DIMENSION(A2D(nn_hls),0:1)     ::   zdkt3d                                           ! vertical tracer gradient at 2 levels
+      REAL(wp), DIMENSION(A2D(nn_hls)        ) ::   z2d                                              ! 2D workspace
+      REAL(wp), DIMENSION(A2D(nn_hls)    ,jpk) ::   zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw   ! 3D     -
       !!----------------------------------------------------------------------
       !
@@ -163,68 +156,53 @@
          END_3D
          !
-         ! NOTE: [halo1-halo2]
-         zftu(:,:,:,:) = 0._wp
-         zftv(:,:,:,:) = 0._wp
-         !
-         DO ip = 0, 1                            ! i-k triads
-            DO kp = 0, 1
-               ! [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 ) 
-                  ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm)
-                  zbu   = e1e2u(ji-ip,jj) * e3u(ji-ip,jj,jk,Kmm)
-                  zah   = 0.25_wp * pahu(ji-ip,jj,jk)
-                  zslope_skew = triadi_g(ji,jj,jk,1-ip,kp)
-                  ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth)
-                  zslope2 = zslope_skew + ( gdept(ji-ip+1,jj,jk,Kmm) - gdept(ji-ip,jj,jk,Kmm) ) * r1_e1u(ji-ip,jj) * umask(ji-ip,jj,jk+kp)
-                  zslope2 = zslope2 *zslope2
-                  ! NOTE: [halo1-halo2]
-                  zftu(ji,jj,jk+kp,ip+1) = zftu(ji,jj,jk+kp,ip+1) + &
-                     &                     zah * zbu * ze3wr * r1_e1e2t(ji,jj) * zslope2
-                  zftv(ji,jj,jk+kp,ip+1) = zftv(ji,jj,jk+kp,ip+1) + &
-                     &                     zah * r1_e1u(ji-ip,jj) * r1_e1u(ji-ip,jj) * umask(ji-ip,jj,jk+kp)
-                     !
-               END_3D
-            END DO
+         DO kp = 0, 1                            ! i-k triads
+            ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
+               ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm)
+               zbu   = e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
+               zbu1  = e1e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)
+               zah   = 0.25_wp * pahu(ji,jj,jk)
+               zah1  = 0.25_wp * pahu(ji-1,jj,jk)
+               ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth)
+               zslope2 = triadi_g(ji,jj,jk,1,kp) + ( gdept(ji+1,jj,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
+               zslope2 = zslope2 *zslope2
+               zslope21 = triadi_g(ji,jj,jk,0,kp) + ( gdept(ji,jj,jk,Kmm) - gdept(ji-1,jj,jk,Kmm) ) * r1_e1u(ji-1,jj) * umask(ji-1,jj,jk+kp)
+               zslope21 = zslope21 *zslope21
+               ! round brackets added to fix the order of floating point operations
+               ! needed to ensure halo 1 - halo 2 compatibility
+               ah_wslp2(ji,jj,jk+kp) =  ah_wslp2(ji,jj,jk+kp) + ( zah * zbu * ze3wr * r1_e1e2t(ji,jj) * zslope2                    &
+                        &                                       + zah1 * zbu1 * ze3wr * r1_e1e2t(ji,jj) * zslope21                 &
+                        &                                       )                                                                  ! bracket for halo 1 - halo 2 compatibility
+               akz     (ji,jj,jk+kp) =  akz     (ji,jj,jk+kp) + ( zah * r1_e1u(ji,jj) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)         &
+                                                                + zah1 * r1_e1u(ji-1,jj) * r1_e1u(ji-1,jj) * umask(ji-1,jj,jk+kp)  &
+                        &                                       )                                                                  ! bracket for halo 1 - halo 2 compatibility
+            END_3D
          END DO
          !
-         ! NOTE: [halo1-halo2] Use DO_3D instead of DO_3D_OVR
-         ! ip loop contributions added here to ensure consistent floating point arithmetic for different nn_hls
-         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
-            ah_wslp2(ji,jj,jk) = ah_wslp2(ji,jj,jk) + (zftu(ji,jj,jk,1) + zftu(ji,jj,jk,2))
-            akz     (ji,jj,jk) = akz     (ji,jj,jk) + (zftv(ji,jj,jk,1) + zftv(ji,jj,jk,2))
-         END_3D
-         !
-         ! NOTE: [halo1-halo2]
-         zftu(:,:,:,:) = 0._wp
-         zftv(:,:,:,:) = 0._wp
-         !
-         DO jp = 0, 1                            ! j-k triads
-            DO kp = 0, 1
-               ! [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 ) 
-                  ze3wr = 1.0_wp / e3w(ji,jj,jk+kp,Kmm)
-                  zbv   = e1e2v(ji,jj-jp) * e3v(ji,jj-jp,jk,Kmm)
-                  zah   = 0.25_wp * pahv(ji,jj-jp,jk)
-                  zslope_skew = triadj_g(ji,jj,jk,1-jp,kp)
-                  ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces
-                  !    (do this by *adding* gradient of depth)
-                  zslope2 = zslope_skew + ( gdept(ji,jj-jp+1,jk,Kmm) - gdept(ji,jj-jp,jk,Kmm) ) * r1_e2v(ji,jj-jp) * vmask(ji,jj-jp,jk+kp)
-                  zslope2 = zslope2 * zslope2
-                  ! NOTE: [halo1-halo2]
-                  zftu(ji,jj,jk+kp,jp+1) = zftu(ji,jj,jk+kp,jp+1) + &
-                     &                     zah * zbv * ze3wr * r1_e1e2t(ji,jj) * zslope2
-                  zftv(ji,jj,jk+kp,jp+1) = zftv(ji,jj,jk+kp,jp+1) + &
-                     &                     zah * r1_e2v(ji,jj-jp) * r1_e2v(ji,jj-jp) * vmask(ji,jj-jp,jk+kp)
-                  !
-               END_3D
-            END DO
+         DO kp = 0, 1                            ! j-k triads
+            ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
+               ze3wr = 1.0_wp / e3w(ji,jj,jk+kp,Kmm)
+               zbv   = e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
+               zbv1   = e1e2v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)
+               zah   = 0.25_wp * pahv(ji,jj,jk)
+               zah1   = 0.25_wp * pahv(ji,jj-1,jk)
+               ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces
+               !    (do this by *adding* gradient of depth)
+               zslope2 = triadj_g(ji,jj,jk,1,kp) + ( gdept(ji,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
+               zslope2 = zslope2 * zslope2
+               zslope21 = triadj_g(ji,jj,jk,0,kp) + ( gdept(ji,jj,jk,Kmm) - gdept(ji,jj-1,jk,Kmm) ) * r1_e2v(ji,jj-1) * vmask(ji,jj-1,jk+kp)
+               zslope21 = zslope21 * zslope21
+               ! round brackets added to fix the order of floating point operations
+               ! needed to ensure halo 1 - halo 2 compatibility
+               ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + ( zah * zbv * ze3wr * r1_e1e2t(ji,jj) * zslope2                     &
+                        &                                      + zah1 * zbv1 * ze3wr * r1_e1e2t(ji,jj) * zslope21                  &
+                        &                                      )                                                                   ! bracket for halo 1 - halo 2 compatibility
+               akz     (ji,jj,jk+kp) = akz     (ji,jj,jk+kp) + ( zah * r1_e2v(ji,jj) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)          &
+                        &                                      + zah1 * r1_e2v(ji,jj-1) * r1_e2v(ji,jj-1) * vmask(ji,jj-1,jk+kp)   &
+                        &                                      )                                                                   ! bracket for halo 1 - halo 2 compatibility
+               !
+            END_3D
          END DO
-         !
-         ! NOTE: [halo1-halo2] Use DO_3D instead of DO_3D_OVR
-         ! jp loop contributions added here to ensure consistent floating point arithmetic for different nn_hls
-         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
-            ah_wslp2(ji,jj,jk) = ah_wslp2(ji,jj,jk) + (zftu(ji,jj,jk,1) + zftu(ji,jj,jk,2))
-            akz     (ji,jj,jk) = akz     (ji,jj,jk) + (zftv(ji,jj,jk,1) + zftv(ji,jj,jk,2))
-         END_3D
          !
          IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient
@@ -259,13 +237,17 @@
             zpsi_vw(:,:,:) = 0._wp
 
-            DO jp = 0, 1
-               DO kp = 0, 1
-                  DO_3D( 1, 0, 1, 0, 1, jpkm1 )
-                     zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) &
-                        & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+jp,jj,jk,1-jp,kp)
-                     zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp) &
-                        & + 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj+jp,jk,1-jp,kp)
-                  END_3D
-               END DO
+            DO kp = 0, 1
+               DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+                  ! round brackets added to fix the order of floating point operations
+                  ! needed to ensure halo 1 - halo 2 compatibility
+                  zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp)                                     &
+                     & + ( 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji,jj,jk,1,kp)        &
+                     &   + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+1,jj,jk,0,kp)      &
+                     &   )                                                                        ! bracket for halo 1 - halo 2 compatibility
+                  zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp)                                     &
+                     & + ( 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj,jk,1,kp)        &
+                     &   + 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj+1,jk,0,kp)      &
+                     &   )                                                                        ! bracket for halo 1 - halo 2 compatibility
+               END_3D
             END DO
             CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm )
@@ -279,8 +261,7 @@
          ! Zero fluxes for each tracer
 !!gm  this should probably be done outside the jn loop
-         ! NOTE: [halo1-halo2]
-         ztfw(:,:,:,:,:) = 0._wp
-         zftu(:,:,:,:) = 0._wp
-         zftv(:,:,:,:) = 0._wp
+         ztfw(:,:,:) = 0._wp
+         zftu(:,:,:) = 0._wp
+         zftv(:,:,:) = 0._wp
          !
          ! [comm_cleanup] ! DO_3D( 1, 0, 1, 0, 1, jpkm1 )    !==  before lateral T & S gradients at T-level jk  ==!
@@ -327,123 +308,167 @@
             !
             IF( ln_botmix_triad ) THEN
-               DO ip = 0, 1              !==  Horizontal & vertical fluxes
-                  DO kp = 0, 1
-                     ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
-                     DO_2D( iij, iij-1, iij, iij-1 )
-                        ze1ur = r1_e1u(ji,jj)
-                        zdxt  = zdit(ji,jj,jk) * ze1ur
-                        ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
-                        zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr
-                        zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
-                        zslope_iso  = triadi  (ji+ip,jj,jk,1-ip,kp)
-                        !
-                        zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
-                        ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????   ahu is masked....
-                        zah = pahu(ji,jj,jk)
-                        zah_slp  = zah * zslope_iso
-                        IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew
-                        ! NOTE: [halo1-halo2]
-                        zftu(ji   ,jj,jk,ip+1  ) = zftu(ji   ,jj,jk,ip+1  ) - &
-                           &                       ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
-                        ztfw(ji,jj,jk+kp,ip+1,1) = ztfw(ji,jj,jk+kp,ip+1,1) - &
-                           &                                      (zah_slp + zaei_slp) * zdxt   * zbu * ze3wr
-                     END_2D
-                  END DO
+               DO kp = 0, 1              !==  Horizontal & vertical fluxes
+                  ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+                  DO_2D( iij, iij-1, iij, iij-1 )
+                     ze1ur = r1_e1u(ji,jj)
+                     zdxt  = zdit(ji,jj,jk) * ze1ur
+                     zdxt_ip1  = zdit(ji+1,jj,jk) * r1_e1u(ji+1,jj)
+                     ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm)
+                     ze3wr_ip1 = 1._wp / e3w(ji+1,jj,jk+kp,Kmm)
+                     zdzt  = zdkt3d(ji,jj,kp) * ze3wr
+                     zdzt_ip1  = zdkt3d(ji+1,jj,kp) * ze3wr_ip1
+                     !
+                     zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
+                     zbu_ip1 = 0.25_wp * e1e2u(ji+1,jj) * e3u(ji+1,jj,jk,Kmm)
+                     ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????   ahu is masked....
+                     zah = pahu(ji,jj,jk)
+                     zah_ip1 = pahu(ji+1,jj,jk)
+                     zah_slp  = zah * triadi(ji,jj,jk,1,kp)
+                     zah_slp_ip1  = zah_ip1 * triadi(ji+1,jj,jk,1,kp)
+                     zah_slp1  = zah * triadi(ji+1,jj,jk,0,kp)
+                     IF( ln_ldfeiv )   THEN
+                        zaei_slp = aeiu(ji,jj,jk) * triadi_g(ji,jj,jk,1,kp)
+                        zaei_slp_ip1 = aeiu(ji+1,jj,jk) * triadi_g(ji+1,jj,jk,1,kp)
+                        zaei_slp1 = aeiu(ji,jj,jk) * triadi_g(ji+1,jj,jk,0,kp)
+                     ENDIF
+                     ! round brackets added to fix the order of floating point operations
+                     ! needed to ensure halo 1 - halo 2 compatibility
+                     zftu(ji   ,jj,jk  ) =  zftu(ji   ,jj,jk )                                                               &
+                                         &    - ( ( zah * zdxt + ( zah_slp - zaei_slp ) * zdzt ) * zbu * ze1ur               &
+                                         &      + ( zah * zdxt + zah_slp1 * zdzt_ip1 - zaei_slp1 * zdzt_ip1 ) * zbu * ze1ur  &
+                                         &      )                                                                            ! bracket for halo 1 - halo 2 compatibility
+                     ztfw(ji+1,jj,jk+kp) =  ztfw(ji+1,jj,jk+kp)                                                              &
+                                         &    - ( (zah_slp_ip1 + zaei_slp_ip1) * zdxt_ip1 * zbu_ip1 * ze3wr_ip1              &
+                                         &      + ( zah_slp1 + zaei_slp1) * zdxt * zbu * ze3wr_ip1                           &
+                                         &      )                                                                            ! bracket for halo 1 - halo 2 compatibility
+                  END_2D
                END DO
                !
-               DO jp = 0, 1
-                  DO kp = 0, 1
-                     ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
-                     DO_2D( iij, iij-1, iij, iij-1 )
-                        ze2vr = r1_e2v(ji,jj)
-                        zdyt  = zdjt(ji,jj,jk) * ze2vr
-                        ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,Kmm)
-                        zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr
-                        zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
-                        zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp)
-                        zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
-                        ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????  ahv is masked...
-                        zah = pahv(ji,jj,jk)
-                        zah_slp = zah * zslope_iso
-                        IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew
-                        ! NOTE: [halo1-halo2]
-                        zftv(ji,jj   ,jk,jp+1  ) = zftv(ji,jj   ,jk,jp+1  ) - &
-                           &                       ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
-                        ztfw(ji,jj,jk+kp,jp+1,2) = ztfw(ji,jj,jk+kp,jp+1,2) - &
-                           &                                      (zah_slp + zaei_slp) * zdyt   * zbv * ze3wr
-                     END_2D
-                  END DO
+               DO kp = 0, 1
+                  ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+                  DO_2D( iij, iij-1, iij, iij-1 )
+                     ze2vr = r1_e2v(ji,jj)
+                     zdyt  = zdjt(ji,jj,jk) * ze2vr
+                     zdyt_jp1  = zdjt(ji,jj+1,jk) * r1_e2v(ji,jj+1)
+                     ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm)
+                     ze3wr_jp1 = 1._wp / e3w(ji,jj+1,jk+kp,Kmm)
+                     zdzt  = zdkt3d(ji,jj,kp) * ze3wr
+                     zdzt_jp1  = zdkt3d(ji,jj+1,kp) * ze3wr_jp1
+                     zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
+                     zbv_jp1 = 0.25_wp * e1e2v(ji,jj+1) * e3v(ji,jj+1,jk,Kmm)
+                     ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????   ahu is masked....
+                     zah = pahv(ji,jj,jk)          ! pahv(ji,jj+jp,jk)  ????
+                     zah_jp1 = pahv(ji,jj+1,jk)
+                     zah_slp = zah * triadj(ji,jj,jk,1,kp)
+                     zah_slp1 = zah * triadj(ji,jj+1,jk,0,kp)
+                     zah_slp_jp1 = zah_jp1 * triadj(ji,jj+1,jk,1,kp)
+                     IF( ln_ldfeiv )   THEN
+                        zaei_slp = aeiv(ji,jj,jk) * triadj_g(ji,jj,jk,1,kp)
+                        zaei_slp_jp1 = aeiv(ji,jj+1,jk) * triadj_g(ji,jj+1,jk,1,kp)
+                        zaei_slp1 = aeiv(ji,jj,jk) * triadj_g(ji,jj+1,jk,0,kp)
+                     ENDIF
+                     ! round brackets added to fix the order of floating point operations
+                     ! needed to ensure halo 1 - halo 2 compatibility
+                     zftv(ji,jj  ,jk   ) =  zftv(ji,jj  ,jk   )                                                              &
+                                         &    - ( ( zah * zdyt + ( zah_slp - zaei_slp ) * zdzt ) * zbv * ze2vr               &
+                                         &      + ( zah * zdyt + zah_slp1 * zdzt_jp1 - zaei_slp1 * zdzt_jp1 ) * zbv * ze2vr  &
+                                         &      )                                                                            ! bracket for halo 1 - halo 2 compatibility
+                     ztfw(ji,jj+1,jk+kp) =  ztfw(ji,jj+1,jk+kp)                                                              &
+                                         &    - ( ( zah_slp_jp1 + zaei_slp_jp1) * zdyt_jp1 * zbv_jp1 * ze3wr_jp1             &
+                                         &      + ( zah_slp1 + zaei_slp1) * zdyt * zbv * ze3wr_jp1                           &
+                                         &      )                                                                            ! bracket for halo 1 - halo 2 compatibility
+                  END_2D
                END DO
                !
             ELSE
                !
-               DO ip = 0, 1               !==  Horizontal & vertical fluxes
-                  DO kp = 0, 1
-                     ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
-                     DO_2D( iij, iij-1, iij, iij-1 )
-                        ze1ur = r1_e1u(ji,jj)
-                        zdxt  = zdit(ji,jj,jk) * ze1ur
-                        ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
-                        zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr
-                        zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
-                        zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp)
-                        !
-                        zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
-                        ! ln_botmix_triad is .F. mask zah for bottom half cells
-                        zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp)         ! pahu(ji+ip,jj,jk)   ===>>  ????
-                        zah_slp  = zah * zslope_iso
-                        IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew        ! aeit(ji+ip,jj,jk)*zslope_skew
-                        ! NOTE: [halo1-halo2]
-                        zftu(ji   ,jj,jk,ip+1  ) = zftu(ji   ,jj,jk,ip+1  ) - &
-                           &                       ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
-                        ztfw(ji,jj,jk+kp,ip+1,1) = ztfw(ji,jj,jk+kp,ip+1,1) - &
-                           &                                      (zah_slp + zaei_slp) * zdxt   * zbu * ze3wr
-                     END_2D
-                  END DO
+               DO kp = 0, 1               !==  Horizontal & vertical fluxes
+                  ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+                  DO_2D( iij, iij-1, iij, iij-1 )
+                     ze1ur = r1_e1u(ji,jj)
+                     zdxt  = zdit(ji,jj,jk) * ze1ur
+                     zdxt_ip1  = zdit(ji+1,jj,jk) * r1_e1u(ji+1,jj)
+                     ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm)
+                     ze3wr_ip1 = 1._wp / e3w(ji+1,jj,jk+kp,Kmm)
+                     zdzt  = zdkt3d(ji,jj,kp) * ze3wr
+                     zdzt_ip1  = zdkt3d(ji+1,jj,kp) * ze3wr_ip1
+                     !
+                     zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
+                     zbu_ip1 = 0.25_wp * e1e2u(ji+1,jj) * e3u(ji+1,jj,jk,Kmm)
+                     ! ln_botmix_triad is .F. mask zah for bottom half cells
+                     zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp)         ! pahu(ji+ip,jj,jk)   ===>>  ????
+                     zah_ip1 = pahu(ji+1,jj,jk) * umask(ji+1,jj,jk+kp)
+                     zah_slp  = zah * triadi(ji,jj,jk,1,kp)
+                     zah_slp_ip1  = zah_ip1 * triadi(ji+1,jj,jk,1,kp)
+                     zah_slp1  = zah * triadi(ji+1,jj,jk,0,kp)
+                     IF( ln_ldfeiv )   THEN
+                        zaei_slp = aeiu(ji,jj,jk) * triadi_g(ji,jj,jk,1,kp)
+                        zaei_slp_ip1 = aeiu(ji+1,jj,jk) * triadi_g(ji+1,jj,jk,1,kp)
+                        zaei_slp1 = aeiu(ji,jj,jk) * triadi_g(ji+1,jj,jk,0,kp)
+                     ENDIF
+!                     zftu(ji   ,jj,jk  ) = zftu(ji   ,jj,jk ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur - ( zah * zdxt + (zah_slp1 - zaei_slp1) * zdzt_ip1 ) * zbu * ze1ur
+!                     ztfw(ji+1,jj,jk+kp) = ztfw(ji+1,jj,jk+kp) - (zah_slp_ip1 + zaei_slp_ip1) * zdxt_ip1 * zbu_ip1 * ze3wr_ip1 - (zah_slp1 + zaei_slp1) * zdxt * zbu * ze3wr_ip1
+                     ! round brackets added to fix the order of floating point operations
+                     ! needed to ensure halo 1 - halo 2 compatibility
+                     zftu(ji   ,jj,jk  ) =  zftu(ji   ,jj,jk )                                                               &
+                                         &    - ( ( zah * zdxt + ( zah_slp - zaei_slp ) * zdzt ) * zbu * ze1ur               &
+                                         &      + ( zah * zdxt + zah_slp1 * zdzt_ip1 - zaei_slp1 * zdzt_ip1 ) * zbu * ze1ur  &
+                                         &      )                                                                            ! bracket for halo 1 - halo 2 compatibility
+                     ztfw(ji+1,jj,jk+kp) =  ztfw(ji+1,jj,jk+kp)                                                              &
+                                         &    - ( (zah_slp_ip1 + zaei_slp_ip1) * zdxt_ip1 * zbu_ip1 * ze3wr_ip1              &
+                                         &      + ( zah_slp1 + zaei_slp1) * zdxt * zbu * ze3wr_ip1                           &
+                                         &      )                                                                            ! bracket for halo 1 - halo 2 compatibility
+                  END_2D
                END DO
                !
-               DO jp = 0, 1
-                  DO kp = 0, 1
-                     ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
-                     DO_2D( iij, iij-1, iij, iij-1 )
-                        ze2vr = r1_e2v(ji,jj)
-                        zdyt  = zdjt(ji,jj,jk) * ze2vr
-                        ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,Kmm)
-                        zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr
-                        zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
-                        zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp)
-                        zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
-                        ! ln_botmix_triad is .F. mask zah for bottom half cells
-                        zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp)         ! pahv(ji,jj+jp,jk)  ????
-                        zah_slp = zah * zslope_iso
-                        IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew        ! aeit(ji,jj+jp,jk)*zslope_skew
-                        ! NOTE: [halo1-halo2]
-                        zftv(ji,jj,jk,   jp+1  ) = zftv(ji,jj,jk,   jp+1  ) - &
-                           &                       ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
-                        ztfw(ji,jj,jk+kp,jp+1,2) = ztfw(ji,jj,jk+kp,jp+1,2) - &
-                           &                                      (zah_slp + zaei_slp) * zdyt   * zbv * ze3wr
-                     END_2D
-                  END DO
+               DO kp = 0, 1
+                  ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
+                  DO_2D( iij, iij-1, iij, iij-1 )
+                     ze2vr = r1_e2v(ji,jj)
+                     zdyt  = zdjt(ji,jj,jk) * ze2vr
+                     zdyt_jp1  = zdjt(ji,jj+1,jk) * r1_e2v(ji,jj+1)
+                     ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm)
+                     ze3wr_jp1 = 1._wp / e3w(ji,jj+1,jk+kp,Kmm)
+                     zdzt  = zdkt3d(ji,jj,kp) * ze3wr
+                     zdzt_jp1  = zdkt3d(ji,jj+1,kp) * ze3wr_jp1
+                     zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
+                     zbv_jp1 = 0.25_wp * e1e2v(ji,jj+1) * e3v(ji,jj+1,jk,Kmm)
+                     ! ln_botmix_triad is .F. mask zah for bottom half cells
+                     zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp)         ! pahv(ji,jj+jp,jk)  ????
+                     zah_jp1 = pahv(ji,jj+1,jk) * vmask(ji,jj+1,jk+kp)
+                     zah_slp = zah * triadj(ji,jj,jk,1,kp)
+                     zah_slp1 = zah * triadj(ji,jj+1,jk,0,kp)
+                     zah_slp_jp1 = zah_jp1 * triadj(ji,jj+1,jk,1,kp)
+                     IF( ln_ldfeiv )   THEN
+                        zaei_slp = aeiv(ji,jj,jk) * triadj_g(ji,jj,jk,1,kp)
+                        zaei_slp_jp1 = aeiv(ji,jj+1,jk) * triadj_g(ji,jj+1,jk,1,kp)
+                        zaei_slp1 = aeiv(ji,jj,jk) * triadj_g(ji,jj+1,jk,0,kp)
+                     ENDIF
+!                     zftv(ji,jj  ,jk   ) = zftv(ji,jj  ,jk   ) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr - ( zah * zdyt + (zah_slp1 - zaei_slp1) * zdzt_jp1 ) * zbv * ze2vr
+!                     ztfw(ji,jj+1,jk+kp) = ztfw(ji,jj+1,jk+kp) - ( zah_slp_jp1 + zaei_slp_jp1) * zdyt_jp1 * zbv_jp1 * ze3wr_jp1 - (zah_slp1 + zaei_slp1) * zdyt * zbv * ze3wr_jp1
+                     ! round brackets added to fix the order of floating point operations
+                     ! needed to ensure halo 1 - halo 2 compatibility
+                     zftv(ji,jj  ,jk   ) =  zftv(ji,jj  ,jk   )                                                              &
+                                         &    - ( ( zah * zdyt + ( zah_slp - zaei_slp ) * zdzt ) * zbv * ze2vr               &
+                                         &      + ( zah * zdyt + zah_slp1 * zdzt_jp1 - zaei_slp1 * zdzt_jp1 ) * zbv * ze2vr  &
+                                         &      )                                                                            ! bracket for halo 1 - halo 2 compatibility
+                     ztfw(ji,jj+1,jk+kp) =  ztfw(ji,jj+1,jk+kp)                                                              &
+                                         &    - ( ( zah_slp_jp1 + zaei_slp_jp1) * zdyt_jp1 * zbv_jp1 * ze3wr_jp1             &
+                                         &      + ( zah_slp1 + zaei_slp1) * zdyt * zbv * ze3wr_jp1                           &
+                                         &      )                                                                            ! bracket for halo 1 - halo 2 compatibility
+                  END_2D
                END DO
             ENDIF
-            ! NOTE: [halo1-halo2]
-            ! ip and jp loop contributions added here to ensure consistent floating point arithmetic for different nn_hls
-            DO_2D( iij, iij-1, iij, iij-1 )
-               zftu(ji,jj,jk,1) = zftu(ji,jj,jk,1) + zftu(ji,jj,jk,2)
-               zftv(ji,jj,jk,1) = zftv(ji,jj,jk,1) + zftv(ji,jj,jk,2)
-            END_2D
-            DO_2D( iij-1, iij-1, iij-1, iij-1 )
-               ztfw(ji,jj,jk,1,1) = (ztfw(ji,jj,jk,1,1) + ztfw(ji-1,jj,jk,2,1)) + &
-                  &                 (ztfw(ji,jj,jk,1,2) + ztfw(ji,jj-1,jk,2,2))
-            END_2D
             !                             !==  horizontal divergence and add to the general trend  ==!
             ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )
             DO_2D( iij-1, iij-1, iij-1, iij-1 )
-               ! NOTE: [halo1-halo2]
-               ! Extra brackets required to ensure consistent floating point arithmetic for different nn_hls for bilaplacian
-               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    &
-                  &                       + zsign * (  (zftu(ji-1,jj,jk,1) - zftu(ji,jj,jk,1))       &
-                  &                                  + (zftv(ji,jj-1,jk,1) - zftv(ji,jj,jk,1))   )   &
-                  &                               / (  e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)  )
+               ! round brackets added to fix the order of floating point operations
+               ! needed to ensure halo 1 - halo 2 compatibility
+               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)                                                &
+                  &                       + zsign * ( ( zftu(ji-1,jj  ,jk) - zftu(ji,jj,jk)             &
+                  &                                   )                                                 & ! bracket for halo 1 - halo 2 compatibility
+                  &                                 + ( zftv(ji,jj-1,jk) - zftv(ji,jj,jk)               &
+                  &                                   )                                                 & ! bracket for halo 1 - halo 2 compatibility
+                  &                                 ) / (  e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)  )
             END_2D
             !
@@ -454,5 +479,5 @@
             ! [comm_cleanup] ! DO_3D( 0, 0, 1, 0, 2, jpkm1 )
             DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 )
-               ztfw(ji,jj,jk,1,1) = ztfw(ji,jj,jk,1,1) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)   &
+               ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)   &
                   &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
                   &                            * (  pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
@@ -463,10 +488,11 @@
                ! [comm_cleanup] ! DO_3D( 0, 0, 1, 0, 2, jpkm1 )
                DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 )
-                  ztfw(ji,jj,jk,1,1) = ztfw(ji,jj,jk,1,1) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)             &
+                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)             &
                      &                            * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
                END_3D
             CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp.
+               ! [comm_cleanup] ! DO_3D( 0, 0, 1, 0, 2, jpkm1 )
                DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-                  ztfw(ji,jj,jk,1,1) = ztfw(ji,jj,jk,1,1) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)                      &
+                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)                      &
                      &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   &
                      &                               + akz     (ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) )   )
@@ -476,6 +502,7 @@
          !
          ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )      !==  Divergence of vertical fluxes added to pta  ==!
-         DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 )
-            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1,1,1) - ztfw(ji,jj,jk,1,1)  )   &
+         DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 )      !==  Divergence of vertical fluxes added to pta  ==!
+            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    &
+            &                                  + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
                &                                              / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
          END_3D
@@ -485,9 +512,7 @@
             !
             !                          ! "Poleward" diffusive heat or salt transports (T-S case only)
-            ! NOTE: [halo1-halo2]
-            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', zftv(:,:,:,1)  )
+            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', zftv(:,:,:)  )
             !                          ! Diffusive heat transports
-            ! NOTE: [halo1-halo2]
-            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', zftu(:,:,:,1), zftv(:,:,:,1) )
+            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', zftu(:,:,:), zftv(:,:,:) )
             !
          ENDIF                                                    !== end pass selection  ==!
Index: /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/zpshde.F90
===================================================================
--- /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/zpshde.F90	(revision 14679)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/zpshde.F90	(revision 14680)
@@ -175,5 +175,5 @@
          pgru(:,:) = 0._wp
          pgrv(:,:) = 0._wp                ! depth of the partial step level
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
             iku = mbku(ji,jj)
             ikv = mbkv(ji,jj)
@@ -191,5 +191,5 @@
          CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj
          !
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )              ! Gradient of density at the last level
+         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )              ! Gradient of density at the last level
             iku = mbku(ji,jj)
             ikv = mbkv(ji,jj)
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 14679)
+++ /NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/stpmlf.F90	(revision 14680)
@@ -187,6 +187,4 @@
             CALL lbc_lnk( 'stp_MLF', avm  , 'W', 1.0_wp , avt  , 'W', 1.0_wp , avs , 'W', 1.0_wp )
          ENDIF
-         ! CALL lbc_lnk( 'stp_MLF', uu(:,:,:,Nnn), 'U', -1.0_wp, vv(:,:,:,Nnn), 'V', -1.0_wp )
-         ! IF(.NOT.lk_linssh) CALL lbc_lnk( 'stp_MLF', r3u(:,:,Nnn), 'U', 1.0_wp, r3v(:,:,Nnn), 'V', 1.0_wp, r3t(:,:,Nnn), 'T', 1.0_wp )
       ENDIF
                          CALL zdf_phy( kstp, Nbb, Nnn, Nrhs )   ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
@@ -286,4 +284,6 @@
                          CALL ssh_atf    ( kstp, Nbb, Nnn, Naa, ssh )            ! time filtering of "now" sea surface height
       IF(.NOT.lk_linssh) CALL dom_qco_r3c( ssh(:,:,Nnn), r3t_f, r3u_f, r3v_f )   ! "now" ssh/h_0 ratio from filtrered ssh
+         ! [comm_cleanup] this should not be needed
+         IF(nn_hls.eq.2.AND..NOT.lk_linssh) CALL lbc_lnk( 'stp_MLF', r3u_f, 'U', 1.0_wp, r3v_f, 'V', 1.0_wp )
 #if defined key_top
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
