Index: NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/TRA/traldf_triad.F90
===================================================================
--- NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/TRA/traldf_triad.F90	(revision 14384)
+++ NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/TRA/traldf_triad.F90	(revision 15127)
@@ -13,6 +13,4 @@
    USE oce            ! ocean dynamics and active tracers
    USE dom_oce        ! ocean space and time domain
-   ! TEMP: [tiling] This change not necessary if XIOS has subdomain support
-   USE domtile
    USE domutl, ONLY : is_tile
    USE phycst         ! physical constants
@@ -109,22 +107,16 @@
       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            ! 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, zftu, zftv, ztfw   ! 3D     -
-      ! TEMP: [tiling] This can be A2D(nn_hls) if XIOS has subdomain support
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpsi_uw, zpsi_vw
+      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     -
       !!----------------------------------------------------------------------
       !
-      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
          IF( kpass == 1 .AND. kt == kit000 )  THEN
             IF(lwp) WRITE(numout,*)
@@ -142,4 +134,10 @@
       ENDIF
       !
+      ! Define pt_rhs halo points for multi-point haloes in bilaplacian case
+      IF( nldf_tra == np_blp_it .AND. kpass == 1 ) THEN ; iij = nn_hls
+      ELSE                                              ; iij = 1
+      ENDIF
+
+      !
       IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
       ELSE                    ;   zsign = -1._wp
@@ -152,44 +150,54 @@
       IF( kpass == 1 ) THEN         !==  first pass only  and whatever the tracer is  ==!
          !
-         DO_3D( 0, 0, 0, 0, 1, jpk )
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
             akz     (ji,jj,jk) = 0._wp
             ah_wslp2(ji,jj,jk) = 0._wp
          END_3D
          !
-         DO ip = 0, 1                            ! i-k triads
-            DO kp = 0, 1
-               DO_3D( 0, 0, 0, 0, 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
-                  ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji,jj) * zslope2
-                  akz     (ji,jj,jk+kp) = akz     (ji,jj,jk+kp) + 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
+            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
          !
-         DO jp = 0, 1                            ! j-k triads
-            DO kp = 0, 1
-               DO_3D( 0, 0, 0, 0, 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
-                  ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj) * zslope2
-                  akz     (ji,jj,jk+kp) = akz     (ji,jj,jk+kp) + 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
+            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
          !
@@ -197,13 +205,13 @@
             !
             IF( ln_traldf_blp ) THEN                ! bilaplacian operator
-               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 )
                   akz(ji,jj,jk) = 16._wp           &
                      &   * ah_wslp2   (ji,jj,jk)   &
                      &   * (  akz     (ji,jj,jk)   &
                      &      + ah_wslp2(ji,jj,jk)   &
-                     &        / ( e3w (ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  )
+                     &        / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  )
                END_3D
             ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
-               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 )
                   ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
                   zcoef0 = rDt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
@@ -213,31 +221,28 @@
            !
          ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
-            DO_3D( 0, 0, 0, 0, 1, jpk )
+            DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
                akz(ji,jj,jk) = ah_wslp2(ji,jj,jk)
             END_3D
          ENDIF
          !
-         ! TEMP: [tiling] These changes not necessary if XIOS has subdomain support
-         IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
-            IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) THEN
-               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )
-
-               zpsi_uw(:,:,:) = 0._wp
-               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
-               END DO
-               CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm )
-
-               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = nijtile )
-            ENDIF
+         IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) THEN
+            zpsi_uw(:,:,:) = 0._wp
+            zpsi_vw(:,:,:) = 0._wp
+
+            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 )
          ENDIF
          !
@@ -252,16 +257,18 @@
          zftu(:,:,:) = 0._wp
          zftv(:,:,:) = 0._wp
-         !
-         DO_3D( 1, 0, 1, 0, 1, jpkm1 )    !==  before lateral T & S gradients at T-level jk  ==!
+         zdit(:,:,:) = 0._wp
+         zdjt(:,:,:) = 0._wp
+         !
+         DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 )    !==  before lateral T & S gradients at T-level jk  ==!
             zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
             zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
          END_3D
          IF( ln_zps .AND. l_grad_zps ) THEN    ! partial steps: correction at top/bottom ocean level
-            DO_2D( 1, 0, 1, 0 )                    ! bottom level
+            DO_2D( iij, iij-1, iij, iij-1 )                    ! bottom level
                zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)
                zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
             END_2D
             IF( ln_isfcav ) THEN                   ! top level (ocean cavities only)
-               DO_2D( 1, 0, 1, 0 )
+               DO_2D( iij, iij-1, iij, iij-1 )
                   IF( miku(ji,jj)  > 1 )   zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn)
                   IF( mikv(ji,jj)  > 1 )   zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn)
@@ -276,5 +283,5 @@
          DO jk = 1, jpkm1
             !                    !==  Vertical tracer gradient at level jk and jk+1
-            DO_2D( 1, 1, 1, 1 )
+            DO_2D( iij, iij, iij, iij )
                zdkt3d(ji,jj,1) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1)
             END_2D
@@ -283,5 +290,5 @@
             IF( jk == 1 ) THEN   ;   zdkt3d(:,:,0) = zdkt3d(:,:,1)
             ELSE
-               DO_2D( 1, 1, 1, 1 )
+               DO_2D( iij, iij, iij, iij )
                   zdkt3d(ji,jj,0) = ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
                END_2D
@@ -289,96 +296,163 @@
             !
             zaei_slp = 0._wp
+            zaei_slp_ip1 = 0._wp
+            zaei_slp_jp1 = 0._wp
+            zaei_slp1 = 0._wp
             !
             IF( ln_botmix_triad ) THEN
-               DO ip = 0, 1              !==  Horizontal & vertical fluxes
-                  DO kp = 0, 1
-                     DO_2D( 1, 0, 1, 0 )
-                        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
-                        zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
-                        ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - ( zah_slp + zaei_slp) * zdxt                 * zbu * ze3wr
-                     END_2D
-                  END DO
+               DO kp = 0, 1              !==  Horizontal & vertical fluxes
+                  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
-                     DO_2D( 1, 0, 1, 0 )
-                        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
-                        zftv(ji,jj   ,jk   ) = zftv(ji,jj   ,jk   ) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
-                        ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - ( zah_slp + zaei_slp ) * zdyt                * zbv * ze3wr
-                     END_2D
-                  END DO
+               DO kp = 0, 1
+                  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
-                     DO_2D( 1, 0, 1, 0 )
-                        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
-                        zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
-                        ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr
-                     END_2D
-                  END DO
+               DO kp = 0, 1               !==  Horizontal & vertical fluxes
+                  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
+                     ! 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
-                     DO_2D( 1, 0, 1, 0 )
-                        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
-                        zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
-                        ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr
-                     END_2D
-                  END DO
+               DO kp = 0, 1
+                  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
+                     ! 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
             !                             !==  horizontal divergence and add to the general trend  ==!
-            DO_2D( 0, 0, 0, 0 )
-               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    &
-                  &                       + zsign * (  zftu(ji-1,jj  ,jk) - zftu(ji,jj,jk)       &
-                  &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   &
-                  &                                        / (  e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)  )
+            DO_2D( iij-1, iij-1, iij-1, iij-1 )
+               ! 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
             !
@@ -387,5 +461,5 @@
          !                                !==  add the vertical 33 flux  ==!
          IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
-            DO_3D( 0, 0, 1, 0, 2, jpkm1 )
+            DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 )
                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) )             &
@@ -395,10 +469,10 @@
             SELECT CASE( kpass )
             CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
-               DO_3D( 0, 0, 1, 0, 2, jpkm1 )
+               DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 )
                   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.
-               DO_3D( 0, 0, 1, 0, 2, jpkm1 )
+               DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                   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) )   &
@@ -408,5 +482,5 @@
          ENDIF
          !
-         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 )      !==  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)  )   &
