Index: NEMO/branches/2020/dev_r13787_SI3-10_EAP/src/OCE/TRA/trabbl.F90
===================================================================
--- NEMO/branches/2020/dev_r13787_SI3-10_EAP/src/OCE/TRA/trabbl.F90	(revision 13796)
+++ NEMO/branches/2020/dev_r13787_SI3-10_EAP/src/OCE/TRA/trabbl.F90	(revision 13991)
@@ -106,4 +106,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts             ! active tracers and RHS of tracer equation
       !
+      INTEGER  ::   ji, jj, jk   ! Dummy loop indices
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
       !!----------------------------------------------------------------------
@@ -112,5 +113,5 @@
       !
       IF( l_trdtra )   THEN                         !* Save the T-S input trends
-         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
+         ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) )
          ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
          ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
@@ -125,8 +126,8 @@
          CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
             &          tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
-         ! lateral boundary conditions ; just need for outputs
-         CALL lbc_lnk_multi( 'trabbl', ahu_bbl, 'U', 1.0_wp , ahv_bbl, 'V', 1.0_wp )
-         CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef
-         CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
+         IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
+            CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef
+            CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
+         ENDIF
          !
       ENDIF
@@ -136,10 +137,13 @@
          CALL tra_bbl_adv( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm )
          IF(sn_cfctl%l_prtctl)   &
-         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
+         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_adv  - Ta: ', mask1=tmask, &
             &          tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
-         ! lateral boundary conditions ; just need for outputs
-         CALL lbc_lnk_multi( '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
+         IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
+            ! lateral boundary conditions ; just need for outputs
+            ! NOTE: [tiling-comms-merge] The diagnostic results change along the north fold if this is removed
+            CALL lbc_lnk_multi( '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
+         ENDIF
          !
       ENDIF
@@ -187,5 +191,5 @@
       INTEGER  ::   ik           ! local integers
       REAL(wp) ::   zbtr         ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::   zptb   ! workspace
+      REAL(wp), DIMENSION(A2D(nn_hls)) ::   zptb   ! workspace
       !!----------------------------------------------------------------------
       !
@@ -235,62 +239,63 @@
       INTEGER  ::   iis , iid , ijs , ijd    ! local integers
       INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
+      INTEGER  ::   isi, isj                 !   -       -
       REAL(wp) ::   zbtr, ztra               ! local scalars
       REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
       !!----------------------------------------------------------------------
       !
+      IF( ntsi == Nis0 ) THEN ; isi = 1 ; ELSE ; isi = 0 ; ENDIF    ! Avoid double-counting when using tiling
+      IF( ntsj == Njs0 ) THEN ; isj = 1 ; ELSE ; isj = 0 ; ENDIF
       !                                                          ! ===========
       DO jn = 1, kjpt                                            ! tracer loop
          !                                                       ! ===========
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west
-               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
-                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
-                  iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
-                  ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
-                  zu_bbl = ABS( utr_bbl(ji,jj) )
-                  !
-                  !                                               ! up  -slope T-point (shelf bottom point)
-                  zbtr = r1_e1e2t(iis,jj) / e3t(iis,jj,ikus,Kmm)
-                  ztra = zu_bbl * ( pt(iid,jj,ikus,jn) - pt(iis,jj,ikus,jn) ) * zbtr
-                  pt_rhs(iis,jj,ikus,jn) = pt_rhs(iis,jj,ikus,jn) + ztra
-                  !
-                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
-                     zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,jk,Kmm)
-                     ztra = zu_bbl * ( pt(iid,jj,jk+1,jn) - pt(iid,jj,jk,jn) ) * zbtr
-                     pt_rhs(iid,jj,jk,jn) = pt_rhs(iid,jj,jk,jn) + ztra
-                  END DO
-                  !
-                  zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,ikud,Kmm)
-                  ztra = zu_bbl * ( pt(iis,jj,ikus,jn) - pt(iid,jj,ikud,jn) ) * zbtr
-                  pt_rhs(iid,jj,ikud,jn) = pt_rhs(iid,jj,ikud,jn) + ztra
-               ENDIF
-               !
-               IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
-                  ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
-                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
-                  ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
-                  zv_bbl = ABS( vtr_bbl(ji,jj) )
-                  !
-                  ! up  -slope T-point (shelf bottom point)
-                  zbtr = r1_e1e2t(ji,ijs) / e3t(ji,ijs,ikvs,Kmm)
-                  ztra = zv_bbl * ( pt(ji,ijd,ikvs,jn) - pt(ji,ijs,ikvs,jn) ) * zbtr
-                  pt_rhs(ji,ijs,ikvs,jn) = pt_rhs(ji,ijs,ikvs,jn) + ztra
-                  !
-                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
-                     zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,jk,Kmm)
-                     ztra = zv_bbl * ( pt(ji,ijd,jk+1,jn) - pt(ji,ijd,jk,jn) ) * zbtr
-                     pt_rhs(ji,ijd,jk,jn) = pt_rhs(ji,ijd,jk,jn)  + ztra
-                  END DO
-                  !                                               ! down-slope T-point (deep bottom point)
-                  zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,ikvd,Kmm)
-                  ztra = zv_bbl * ( pt(ji,ijs,ikvs,jn) - pt(ji,ijd,ikvd,jn) ) * zbtr
-                  pt_rhs(ji,ijd,ikvd,jn) = pt_rhs(ji,ijd,ikvd,jn) + ztra
-               ENDIF
-            END DO
+         ! NOTE: [tiling-comms-merge] Bug fix- correct order of indices
+         DO_2D( isj, 0, isi, 0 )            ! CAUTION start from i=1 to update i=2 when cyclic east-west
+            IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
+               ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
+               iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
+               ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
+               zu_bbl = ABS( utr_bbl(ji,jj) )
+               !
+               !                                               ! up  -slope T-point (shelf bottom point)
+               zbtr = r1_e1e2t(iis,jj) / e3t(iis,jj,ikus,Kmm)
+               ztra = zu_bbl * ( pt(iid,jj,ikus,jn) - pt(iis,jj,ikus,jn) ) * zbtr
+               pt_rhs(iis,jj,ikus,jn) = pt_rhs(iis,jj,ikus,jn) + ztra
+               !
+               DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
+                  zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,jk,Kmm)
+                  ztra = zu_bbl * ( pt(iid,jj,jk+1,jn) - pt(iid,jj,jk,jn) ) * zbtr
+                  pt_rhs(iid,jj,jk,jn) = pt_rhs(iid,jj,jk,jn) + ztra
+               END DO
+               !
+               zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,ikud,Kmm)
+               ztra = zu_bbl * ( pt(iis,jj,ikus,jn) - pt(iid,jj,ikud,jn) ) * zbtr
+               pt_rhs(iid,jj,ikud,jn) = pt_rhs(iid,jj,ikud,jn) + ztra
+            ENDIF
             !
-         END DO
-         !                                                  ! ===========
-      END DO                                                ! end tracer
-      !                                                     ! ===========
+            IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
+               ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
+               ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
+               ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
+               zv_bbl = ABS( vtr_bbl(ji,jj) )
+               !
+               ! up  -slope T-point (shelf bottom point)
+               zbtr = r1_e1e2t(ji,ijs) / e3t(ji,ijs,ikvs,Kmm)
+               ztra = zv_bbl * ( pt(ji,ijd,ikvs,jn) - pt(ji,ijs,ikvs,jn) ) * zbtr
+               pt_rhs(ji,ijs,ikvs,jn) = pt_rhs(ji,ijs,ikvs,jn) + ztra
+               !
+               DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
+                  zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,jk,Kmm)
+                  ztra = zv_bbl * ( pt(ji,ijd,jk+1,jn) - pt(ji,ijd,jk,jn) ) * zbtr
+                  pt_rhs(ji,ijd,jk,jn) = pt_rhs(ji,ijd,jk,jn)  + ztra
+               END DO
+               !                                               ! down-slope T-point (deep bottom point)
+               zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,ikvd,Kmm)
+               ztra = zv_bbl * ( pt(ji,ijs,ikvs,jn) - pt(ji,ijd,ikvd,jn) ) * zbtr
+               pt_rhs(ji,ijd,ikvd,jn) = pt_rhs(ji,ijd,ikvd,jn) + ztra
+            ENDIF
+         END_2D
+         !                                                       ! ===========
+      END DO                                                     ! end tracer
+      !                                                          ! ===========
    END SUBROUTINE tra_bbl_adv
 
@@ -333,12 +338,14 @@
       REAL(wp) ::   za, zb, zgdrho            ! local scalars
       REAL(wp) ::   zsign, zsigna, zgbbl      !   -      -
-      REAL(wp), DIMENSION(jpi,jpj,jpts)   :: zts, zab         ! 3D workspace
-      REAL(wp), DIMENSION(jpi,jpj)        :: zub, zvb, zdep   ! 2D workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == kit000 )  THEN
-         IF(lwp)  WRITE(numout,*)
-         IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
-         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
+      REAL(wp), DIMENSION(A2D(nn_hls),jpts)   :: zts, zab         ! 3D workspace
+      REAL(wp), DIMENSION(A2D(nn_hls))        :: zub, zvb, zdep   ! 2D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+         IF( kt == kit000 )  THEN
+            IF(lwp)  WRITE(numout,*)
+            IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
+            IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
+         ENDIF
       ENDIF
       !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
