Index: NEMO/trunk/src/OCE/DOM/dom_oce.F90
===================================================================
--- NEMO/trunk/src/OCE/DOM/dom_oce.F90	(revision 13970)
+++ NEMO/trunk/src/OCE/DOM/dom_oce.F90	(revision 13982)
@@ -74,4 +74,14 @@
    LOGICAL, PUBLIC ::   l_Iperio, l_Jperio   !   should we explicitely take care I/J periodicity 
 
+   ! Tiling namelist
+   LOGICAL, PUBLIC ::   ln_tile
+   INTEGER         ::   nn_ltile_i, nn_ltile_j
+
+   ! Domain tiling (all tiles)
+   INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ntsi_a       !: start of internal part of tile domain
+   INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ntsj_a       !
+   INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ntei_a       !: end of internal part of tile domain
+   INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ntej_a       !
+
    !                             !: domain MPP decomposition parameters
    INTEGER             , PUBLIC ::   nimpp, njmpp     !: i- & j-indexes for mpp-subdomain left bottom
@@ -87,4 +97,8 @@
    INTEGER, PUBLIC ::   noea, nowe        !: index of the local neighboring processors in
    INTEGER, PUBLIC ::   noso, nono        !: east, west, south and north directions
+   INTEGER, PUBLIC ::   nones, nonws        !: north-east, north-west directions for sending 
+   INTEGER, PUBLIC ::   noses, nosws        !: south-east, south-west directions for sending
+   INTEGER, PUBLIC ::   noner, nonwr        !: north-east, north-west directions for receiving
+   INTEGER, PUBLIC ::   noser, noswr        !: south-east, south-west directions for receiving
    INTEGER, PUBLIC ::   nidom             !: ???
 
@@ -296,9 +310,9 @@
       ALLOCATE( e3t(jpi,jpj,jpk,jpt) , e3u (jpi,jpj,jpk,jpt) , e3v (jpi,jpj,jpk,jpt) , e3f(jpi,jpj,jpk) ,      &
          &      e3w(jpi,jpj,jpk,jpt) , e3uw(jpi,jpj,jpk,jpt) , e3vw(jpi,jpj,jpk,jpt)                    ,  STAT=ierr(ii) )
-#endif  
+#endif
          !
       ii = ii+1
       ALLOCATE( r3t  (jpi,jpj,jpt)   , r3u  (jpi,jpj,jpt)    , r3v  (jpi,jpj,jpt)    , r3f  (jpi,jpj) ,  &
-         &      r3t_f(jpi,jpj)       , r3u_f(jpi,jpj)        , r3v_f(jpi,jpj)                         ,  STAT=ierr(ii) )       
+         &      r3t_f(jpi,jpj)       , r3u_f(jpi,jpj)        , r3v_f(jpi,jpj)                         ,  STAT=ierr(ii) )
          !
       ii = ii+1
@@ -317,5 +331,5 @@
          !
       ii = ii+1
-      ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(ii)  ) 
+      ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(ii)  )
          !
       ii = ii+1
@@ -323,5 +337,5 @@
          !
       ii = ii+1
-      ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                        & 
+      ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                        &
          &      ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) ,     &
          &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) ,                    STAT=ierr(ii) )
@@ -331,5 +345,5 @@
          !
       ii = ii+1
-      ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) ,     & 
+      ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) ,     &
          &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk) , STAT=ierr(ii) )
          !
Index: NEMO/trunk/src/OCE/DOM/domain.F90
===================================================================
--- NEMO/trunk/src/OCE/DOM/domain.F90	(revision 13970)
+++ NEMO/trunk/src/OCE/DOM/domain.F90	(revision 13982)
@@ -45,4 +45,5 @@
    USE closea , ONLY : dom_clo ! closed seas
    !
+   USE prtctl         ! Print control (prt_ctl_info routine)
    USE in_out_manager ! I/O manager
    USE iom            ! I/O library
@@ -55,4 +56,5 @@
    PUBLIC   dom_init     ! called by nemogcm.F90
    PUBLIC   domain_cfg   ! called by nemogcm.F90
+   PUBLIC   dom_tile     ! called by step.F90
 
    !!-------------------------------------------------------------------------
@@ -124,6 +126,8 @@
       !           !==  Reference coordinate system  ==!
       !
-      CALL dom_glo                     ! global domain versus local domain
-      CALL dom_nam                     ! read namelist ( namrun, namdom )
+      CALL dom_glo                            ! global domain versus local domain
+      CALL dom_nam                            ! read namelist ( namrun, namdom )
+      CALL dom_tile( ntsi, ntsj, ntei, ntej ) ! Tile domain
+
       !
       CALL dom_hgr                      ! Horizontal mesh
@@ -273,4 +277,87 @@
 
 
+   SUBROUTINE dom_tile( ktsi, ktsj, ktei, ktej, ktile )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE dom_tile  ***
+      !!
+      !! ** Purpose :   Set tile domain variables
+      !!
+      !! ** Action  : - ktsi, ktsj     : start of internal part of domain
+      !!              - ktei, ktej     : end of internal part of domain
+      !!              - ntile          : current tile number
+      !!              - nijtile        : total number of tiles
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(out) :: ktsi, ktsj, ktei, ktej      ! Tile domain indices
+      INTEGER, INTENT(in), OPTIONAL :: ktile              ! Tile number
+      INTEGER ::   jt                                     ! dummy loop argument
+      INTEGER ::   iitile, ijtile                         ! Local integers
+      CHARACTER (len=11) ::   charout
+      !!----------------------------------------------------------------------
+      IF( PRESENT(ktile) .AND. ln_tile ) THEN
+         ntile = ktile                 ! Set domain indices for tile
+         ktsi = ntsi_a(ktile)
+         ktsj = ntsj_a(ktile)
+         ktei = ntei_a(ktile)
+         ktej = ntej_a(ktile)
+
+         IF(sn_cfctl%l_prtctl) THEN
+            WRITE(charout, FMT="('ntile =', I4)") ktile
+            CALL prt_ctl_info( charout )
+         ENDIF
+      ELSE
+         ntile = 0                     ! Initialise to full domain
+         nijtile = 1
+         ktsi = Nis0
+         ktsj = Njs0
+         ktei = Nie0
+         ktej = Nje0
+
+         IF( ln_tile ) THEN            ! Calculate tile domain indices
+            iitile = Ni_0 / nn_ltile_i       ! Number of tiles
+            ijtile = Nj_0 / nn_ltile_j
+            IF( MOD( Ni_0, nn_ltile_i ) /= 0 ) iitile = iitile + 1
+            IF( MOD( Nj_0, nn_ltile_j ) /= 0 ) ijtile = ijtile + 1
+
+            nijtile = iitile * ijtile
+            ALLOCATE( ntsi_a(0:nijtile), ntsj_a(0:nijtile), ntei_a(0:nijtile), ntej_a(0:nijtile) )
+
+            ntsi_a(0) = ktsi                 ! Full domain
+            ntsj_a(0) = ktsj
+            ntei_a(0) = ktei
+            ntej_a(0) = ktej
+
+            DO jt = 1, nijtile               ! Tile domains
+               ntsi_a(jt) = Nis0 + nn_ltile_i * MOD(jt - 1, iitile)
+               ntsj_a(jt) = Njs0 + nn_ltile_j * ((jt - 1) / iitile)
+               ntei_a(jt) = MIN(ntsi_a(jt) + nn_ltile_i - 1, Nie0)
+               ntej_a(jt) = MIN(ntsj_a(jt) + nn_ltile_j - 1, Nje0)
+            ENDDO
+         ENDIF
+
+         IF(lwp) THEN                  ! control print
+            WRITE(numout,*)
+            WRITE(numout,*) 'dom_tile : Domain tiling decomposition'
+            WRITE(numout,*) '~~~~~~~~'
+            IF( ln_tile ) THEN
+               WRITE(numout,*) iitile, 'tiles in i'
+               WRITE(numout,*) '    Starting indices'
+               WRITE(numout,*) '        ', (ntsi_a(jt), jt=1, iitile)
+               WRITE(numout,*) '    Ending indices'
+               WRITE(numout,*) '        ', (ntei_a(jt), jt=1, iitile)
+               WRITE(numout,*) ijtile, 'tiles in j'
+               WRITE(numout,*) '    Starting indices'
+               WRITE(numout,*) '        ', (ntsj_a(jt), jt=1, nijtile, iitile)
+               WRITE(numout,*) '    Ending indices'
+               WRITE(numout,*) '        ', (ntej_a(jt), jt=1, nijtile, iitile)
+            ELSE
+               WRITE(numout,*) 'No domain tiling'
+               WRITE(numout,*) '    i indices =', ktsi, ':', ktei
+               WRITE(numout,*) '    j indices =', ktsj, ':', ktej
+            ENDIF
+         ENDIF
+      ENDIF
+   END SUBROUTINE dom_tile
+
+
    SUBROUTINE dom_nam
       !!----------------------------------------------------------------------
@@ -281,4 +368,5 @@
       !! ** input   : - namrun namelist
       !!              - namdom namelist
+      !!              - namtile namelist
       !!              - namnc4 namelist   ! "key_netcdf4" only
       !!----------------------------------------------------------------------
@@ -293,4 +381,5 @@
          &             ln_cfmeta, ln_xios_read, nn_wxios
       NAMELIST/namdom/ ln_linssh, rn_Dt, rn_atfp, ln_crs, ln_meshmask
+      NAMELIST/namtile/ ln_tile, nn_ltile_i, nn_ltile_j
 #if defined key_netcdf4
       NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
@@ -429,4 +518,24 @@
       r1_Dt = 1._wp / rDt
 
+      READ  ( numnam_ref, namtile, IOSTAT = ios, ERR = 905 )
+905   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtile in reference namelist' )
+      READ  ( numnam_cfg, namtile, IOSTAT = ios, ERR = 906 )
+906   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtile in configuration namelist' )
+      IF(lwm) WRITE( numond, namtile )
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*)    '   Namelist : namtile   ---   Domain tiling decomposition'
+         WRITE(numout,*)    '      Tiling (T) or not (F)                ln_tile    = ', ln_tile
+         WRITE(numout,*)    '      Length of tile in i                  nn_ltile_i = ', nn_ltile_i
+         WRITE(numout,*)    '      Length of tile in j                  nn_ltile_j = ', nn_ltile_j
+         WRITE(numout,*)
+         IF( ln_tile ) THEN
+            WRITE(numout,*) '      The domain will be decomposed into tiles of size', nn_ltile_i, 'x', nn_ltile_j
+         ELSE
+            WRITE(numout,*) '      Domain tiling will NOT be used'
+         ENDIF
+      ENDIF
+
       IF( TRIM(Agrif_CFixed()) == '0' ) THEN
          lrxios = ln_xios_read.AND.ln_rstart
Index: NEMO/trunk/src/OCE/DOM/domutl.F90
===================================================================
--- NEMO/trunk/src/OCE/DOM/domutl.F90	(revision 13970)
+++ NEMO/trunk/src/OCE/DOM/domutl.F90	(revision 13982)
@@ -21,6 +21,11 @@
    PRIVATE
 
+   INTERFACE is_tile
+      MODULE PROCEDURE is_tile_2d, is_tile_3d, is_tile_4d
+   END INTERFACE is_tile
+
    PUBLIC dom_ngb    ! routine called in iom.F90 module
    PUBLIC dom_uniq   ! Called by dommsk and domwri
+   PUBLIC is_tile
 
    !!----------------------------------------------------------------------
@@ -109,5 +114,44 @@
       !
    END SUBROUTINE dom_uniq
-   
+
+
+   FUNCTION is_tile_2d( pt )
+      !!
+      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pt
+      INTEGER :: is_tile_2d
+      !!
+      IF( ln_tile .AND. (SIZE(pt, 1) < jpi .OR. SIZE(pt, 2) < jpj) ) THEN
+         is_tile_2d = 1
+      ELSE
+         is_tile_2d = 0
+      ENDIF
+   END FUNCTION is_tile_2d
+
+
+   FUNCTION is_tile_3d( pt )
+      !!
+      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   pt
+      INTEGER :: is_tile_3d
+      !!
+      IF( ln_tile .AND. (SIZE(pt, 1) < jpi .OR. SIZE(pt, 2) < jpj) ) THEN
+         is_tile_3d = 1
+      ELSE
+         is_tile_3d = 0
+      ENDIF
+   END FUNCTION is_tile_3d
+
+
+   FUNCTION is_tile_4d( pt )
+      !!
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(in) ::   pt
+      INTEGER :: is_tile_4d
+      !!
+      IF( ln_tile .AND. (SIZE(pt, 1) < jpi .OR. SIZE(pt, 2) < jpj) ) THEN
+         is_tile_4d = 1
+      ELSE
+         is_tile_4d = 0
+      ENDIF
+   END FUNCTION is_tile_4d
+
    !!======================================================================
 END MODULE domutl
Index: NEMO/trunk/src/OCE/DOM/domvvl.F90
===================================================================
--- NEMO/trunk/src/OCE/DOM/domvvl.F90	(revision 13970)
+++ NEMO/trunk/src/OCE/DOM/domvvl.F90	(revision 13982)
@@ -422,5 +422,4 @@
          !                             (stored for tracer advction and continuity equation)
          CALL lbc_lnk_multi( 'domvvl', un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
-
          ! 4 - Time stepping of baroclinic scale factors
          ! ---------------------------------------------
Index: NEMO/trunk/src/OCE/DOM/dtatsd.F90
===================================================================
--- NEMO/trunk/src/OCE/DOM/dtatsd.F90	(revision 13970)
+++ NEMO/trunk/src/OCE/DOM/dtatsd.F90	(revision 13982)
@@ -18,4 +18,5 @@
    USE phycst          ! physical constants
    USE dom_oce         ! ocean space and time domain
+   USE domain, ONLY : dom_tile
    USE fldread         ! read input fields
    !
@@ -135,56 +136,67 @@
       !! ** Action  :   ptsd   T-S data on medl mesh and interpolated at time-step kt
       !!----------------------------------------------------------------------
-      INTEGER                              , INTENT(in   ) ::   kt     ! ocean time-step
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   ptsd   ! T & S data
+      INTEGER                          , INTENT(in   ) ::   kt     ! ocean time-step
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpts), INTENT(  out) ::   ptsd   ! T & S data
       !
       INTEGER ::   ji, jj, jk, jl, jkk   ! dummy loop indicies
       INTEGER ::   ik, il0, il1, ii0, ii1, ij0, ij1   ! local integers
+      INTEGER ::   itile
       REAL(wp)::   zl, zi                             ! local scalars
       REAL(wp), DIMENSION(jpk) ::  ztp, zsp   ! 1D workspace
       !!----------------------------------------------------------------------
       !
-      CALL fld_read( kt, 1, sf_tsd )      !==   read T & S data at kt time step   ==!
+      IF( ntile == 0 .OR. ntile == 1 )  THEN                                         ! Do only for the full domain
+         itile = ntile
+         IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )            ! Use full domain
+            CALL fld_read( kt, 1, sf_tsd )   !==   read T & S data at kt time step   ==!
       !
       !
 !!gm  This should be removed from the code   ===>>>>  T & S files has to be changed
-      !
-      !                                   !==   ORCA_R2 configuration and T & S damping   ==! 
-      IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN
-         IF( nn_cfg == 2 .AND. ln_tsd_dmp ) THEN    ! some hand made alterations
-            !
-            ij0 = 101 + nn_hls       ;   ij1 = 109 + nn_hls                       ! Reduced T & S in the Alboran Sea
-            ii0 = 141 + nn_hls - 1   ;   ii1 = 155 + nn_hls - 1
-            DO jj = mj0(ij0), mj1(ij1)
-               DO ji = mi0(ii0), mi1(ii1)
-                  sf_tsd(jp_tem)%fnow(ji,jj,13:13) = sf_tsd(jp_tem)%fnow(ji,jj,13:13) - 0.20_wp
-                  sf_tsd(jp_tem)%fnow(ji,jj,14:15) = sf_tsd(jp_tem)%fnow(ji,jj,14:15) - 0.35_wp
-                  sf_tsd(jp_tem)%fnow(ji,jj,16:25) = sf_tsd(jp_tem)%fnow(ji,jj,16:25) - 0.40_wp
-                  !
-                  sf_tsd(jp_sal)%fnow(ji,jj,13:13) = sf_tsd(jp_sal)%fnow(ji,jj,13:13) - 0.15_wp
-                  sf_tsd(jp_sal)%fnow(ji,jj,14:15) = sf_tsd(jp_sal)%fnow(ji,jj,14:15) - 0.25_wp
-                  sf_tsd(jp_sal)%fnow(ji,jj,16:17) = sf_tsd(jp_sal)%fnow(ji,jj,16:17) - 0.30_wp
-                  sf_tsd(jp_sal)%fnow(ji,jj,18:25) = sf_tsd(jp_sal)%fnow(ji,jj,18:25) - 0.35_wp
+         !
+         !                                   !==   ORCA_R2 configuration and T & S damping   ==!
+         IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN
+            IF( nn_cfg == 2 .AND. ln_tsd_dmp ) THEN    ! some hand made alterations
+               !
+               ij0 = 101 + nn_hls       ;   ij1 = 109 + nn_hls                       ! Reduced T & S in the Alboran Sea
+               ii0 = 141 + nn_hls - 1   ;   ii1 = 155 + nn_hls - 1
+               DO jj = mj0(ij0), mj1(ij1)
+                  DO ji = mi0(ii0), mi1(ii1)
+                     sf_tsd(jp_tem)%fnow(ji,jj,13:13) = sf_tsd(jp_tem)%fnow(ji,jj,13:13) - 0.20_wp
+                     sf_tsd(jp_tem)%fnow(ji,jj,14:15) = sf_tsd(jp_tem)%fnow(ji,jj,14:15) - 0.35_wp
+                     sf_tsd(jp_tem)%fnow(ji,jj,16:25) = sf_tsd(jp_tem)%fnow(ji,jj,16:25) - 0.40_wp
+                     !
+                     sf_tsd(jp_sal)%fnow(ji,jj,13:13) = sf_tsd(jp_sal)%fnow(ji,jj,13:13) - 0.15_wp
+                     sf_tsd(jp_sal)%fnow(ji,jj,14:15) = sf_tsd(jp_sal)%fnow(ji,jj,14:15) - 0.25_wp
+                     sf_tsd(jp_sal)%fnow(ji,jj,16:17) = sf_tsd(jp_sal)%fnow(ji,jj,16:17) - 0.30_wp
+                     sf_tsd(jp_sal)%fnow(ji,jj,18:25) = sf_tsd(jp_sal)%fnow(ji,jj,18:25) - 0.35_wp
+                  END DO
                END DO
-            END DO
-            ij0 =  87 + nn_hls       ;   ij1 =  96 + nn_hls                       ! Reduced temperature in Red Sea
-            ii0 = 148 + nn_hls - 1   ;   ii1 = 160 + nn_hls - 1
-            sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ,  4:10 ) = 7.0_wp
-            sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp
-            sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp
-         ENDIF
-      ENDIF
+               ij0 =  87 + nn_hls       ;   ij1 =  96 + nn_hls                       ! Reduced temperature in Red Sea
+               ii0 = 148 + nn_hls - 1   ;   ii1 = 160 + nn_hls - 1
+               sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ,  4:10 ) = 7.0_wp
+               sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp
+               sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp
+            ENDIF
+         ENDIF
 !!gm end
-      !
-      ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:)    ! NO mask
-      ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:) 
+         IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = itile )            ! Revert to tile domain
+      ENDIF
+      !
+      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+         ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jk)    ! NO mask
+         ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jk)
+      END_3D
       !
       IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==!
          !
-         IF( kt == nit000 .AND. lwp )THEN
-            WRITE(numout,*)
-            WRITE(numout,*) 'dta_tsd: interpolates T & S data onto the s- or mixed s-z-coordinate mesh'
-         ENDIF
-         !
-         DO_2D( 1, 1, 1, 1 )                  ! vertical interpolation of T & S
+         IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+            IF( kt == nit000 .AND. lwp )THEN
+               WRITE(numout,*)
+               WRITE(numout,*) 'dta_tsd: interpolates T & S data onto the s- or mixed s-z-coordinate mesh'
+            ENDIF
+         ENDIF
+         !
+         ! NOTE: [tiling-comms-merge] This fix was necessary to take out tra_adv lbc_lnk statements in the zps case, but did not work in the zco case
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )                  ! vertical interpolation of T & S
             DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
                zl = gdept_0(ji,jj,jk)
@@ -215,9 +227,12 @@
       ELSE                                !==   z- or zps- coordinate   ==!
          !                             
-         ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) * tmask(:,:,:)    ! Mask
-         ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) * tmask(:,:,:)
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+            ptsd(ji,jj,jk,jp_tem) = ptsd(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)    ! Mask
+            ptsd(ji,jj,jk,jp_sal) = ptsd(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
+         END_3D
          !
          IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
-            DO_2D( 1, 1, 1, 1 )
+            ! NOTE: [tiling-comms-merge] This fix was necessary to take out tra_adv lbc_lnk statements in the zps case
+            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
                ik = mbkt(ji,jj) 
                IF( ik > 1 ) THEN
