Index: NEMO/trunk/src/OCE/TRA/traadv_fct.F90
===================================================================
--- NEMO/trunk/src/OCE/TRA/traadv_fct.F90	(revision 14718)
+++ NEMO/trunk/src/OCE/TRA/traadv_fct.F90	(revision 14820)
@@ -238,5 +238,6 @@
                END_2D
             END DO
-            CALL lbc_lnk( 'traadv_fct', zltu, 'T', 1.0_wp , zltv, 'T', 1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
+            ! NOTE [ comm_cleanup ] : need to change sign to ensure halo 1 - halo 2 compatibility
+            CALL lbc_lnk( 'traadv_fct', zltu, 'T', -1.0_wp , zltv, 'T', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
             !
             DO_3D( 1, 0, 1, 0, 1, jpkm1 )
@@ -244,6 +245,12 @@
                zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
                !                                                        ! C4 minus upstream advective fluxes
-               zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk)
-               zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk)
+               ! round brackets added to fix the order of floating point operations
+               ! needed to ensure halo 1 - halo 2 compatibility
+               zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * ( zC2t_u + ( zltu(ji,jj,jk) - zltu(ji+1,jj,jk)   &
+                             &                                     )                                     & ! bracket for halo 1 - halo 2 compatibility
+                             &                          ) - zwx(ji,jj,jk)
+               zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * ( zC2t_v + ( zltv(ji,jj,jk) - zltv(ji,jj+1,jk)   &
+                             &                                     )                                     & ! bracket for halo 1 - halo 2 compatibility
+                             &                          ) - zwy(ji,jj,jk)
             END_3D
             IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv_fct', zwx, 'U', -1.0_wp, zwy, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
Index: NEMO/trunk/src/OCE/TRA/trabbl.F90
===================================================================
--- NEMO/trunk/src/OCE/TRA/trabbl.F90	(revision 14718)
+++ NEMO/trunk/src/OCE/TRA/trabbl.F90	(revision 14820)
@@ -141,5 +141,5 @@
          IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
             ! lateral boundary conditions ; just need for outputs
-            CALL lbc_lnk( 'trabbl', utr_bbl, 'U', 1.0_wp , vtr_bbl, 'V', 1.0_wp )
+            CALL lbc_lnk( 'trabbl', utr_bbl, 'U', -1.0_wp , vtr_bbl, 'V', -1.0_wp )
             CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport
             CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
Index: NEMO/trunk/src/OCE/TRA/traldf_iso.F90
===================================================================
--- NEMO/trunk/src/OCE/TRA/traldf_iso.F90	(revision 14718)
+++ NEMO/trunk/src/OCE/TRA/traldf_iso.F90	(revision 14820)
@@ -179,8 +179,16 @@
                &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  )
                !
-            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)                   &
+               &         )                                                         & ! bracket for halo 1 - halo 2 compatibility
+               &     ) * zmsku                                       
+            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)                   & 
+               &         )                                                         & ! bracket for halo 1 - halo 2 compatibility
+               &     ) * zmskv                                               
                !
             ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
@@ -190,9 +198,14 @@
          IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient
             DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+               ! 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  ,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) )   )
+                  &            )                                                                               & ! 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
             !
@@ -278,16 +291,29 @@
                zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
                !
-               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
             !
             DO_2D( 0, 0, 0, 0 )           !== horizontal divergence and add to pta
-               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
@@ -317,8 +343,16 @@
             zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk)
             !
-            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/trunk/src/OCE/TRA/traldf_lap_blp.F90
===================================================================
--- NEMO/trunk/src/OCE/TRA/traldf_lap_blp.F90	(revision 14718)
+++ NEMO/trunk/src/OCE/TRA/traldf_lap_blp.F90	(revision 14820)
@@ -159,7 +159,11 @@
          !
          DO_3D( isi, iei, isj, iej, 1, jpkm1 )            !== Second derivative (divergence) added to the general tracer trends  ==!
-            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) )
+            ! 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/trunk/src/OCE/TRA/traldf_triad.F90
===================================================================
--- NEMO/trunk/src/OCE/TRA/traldf_triad.F90	(revision 14718)
+++ NEMO/trunk/src/OCE/TRA/traldf_triad.F90	(revision 14820)
@@ -109,14 +109,10 @@
       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   ! 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) ::   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
@@ -157,39 +153,49 @@
          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( 0, 0, 0, 0, 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( 0, 0, 0, 0, 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
          !
@@ -223,17 +229,21 @@
                IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )
 
-               zpsi_uw(:,:,:) = 0._wp
-               zpsi_vw(:,:,:) = 0._wp
-
-               DO jp = 0, 1
+                  zpsi_uw(:,:,:) = 0._wp
+                  zpsi_vw(:,:,:) = 0._wp
+
                   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)
+                        ! round brackets added to fix the order of floating point operations
+                        ! needed to ensure halo 1 - halo 2 compatibility [ NOT TESTED ]
+                        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
-               END DO
                CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm )
 
@@ -289,96 +299,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
+               DO kp = 0, 1              !==  Horizontal & vertical fluxes
+                  DO_2D( 1, 0, 1, 0 )
+                     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
-               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( 1, 0, 1, 0 )
+                     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
+                  DO_2D( 1, 0, 1, 0 )
+                     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( 1, 0, 1, 0 )
+                     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)  )
+               ! 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
             !
