Index: /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/ASM/asminc.F90
===================================================================
--- /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/ASM/asminc.F90	(revision 13517)
+++ /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/ASM/asminc.F90	(revision 13518)
@@ -520,12 +520,15 @@
       INTEGER  :: it
       REAL(wp) :: zincwgt  ! IAU weight for current time step
-      REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values
+      REAL(wp), DIMENSION(ST_2D(nn_hls),jpk) :: fzptnz ! 3d freezing point values
       !!----------------------------------------------------------------------
       !
       ! freezing point calculation taken from oc_fz_pt (but calculated for all depths) 
       ! used to prevent the applied increments taking the temperature below the local freezing point 
-      DO jk = 1, jpkm1
-        CALL eos_fzp( pts(:,:,jk,jp_sal,Kmm), fzptnz(:,:,jk), gdept(:,:,jk,Kmm) )
-      END DO
+      ! TODO: NOT TESTED- logical is forced to False
+      IF( ln_temnofreeze ) THEN
+         DO jk = 1, jpkm1
+           CALL eos_fzp( pts(:,:,jk,jp_sal,Kmm), fzptnz(:,:,jk), gdept(:,:,jk,Kmm) )
+         END DO
+      ENDIF
          !
          !                             !--------------------------------------
@@ -538,30 +541,37 @@
             zincwgt = wgtiau(it) / rn_Dt   ! IAU weight for the current time step
             !
-            IF(lwp) THEN
-               WRITE(numout,*) 
-               WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
-               WRITE(numout,*) '~~~~~~~~~~~~'
+            IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+               IF(lwp) THEN
+                  WRITE(numout,*)
+                  WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
+                  WRITE(numout,*) '~~~~~~~~~~~~'
+               ENDIF
             ENDIF
             !
             ! Update the tracer tendencies
+            ! TODO: NOT TESTED- logical is forced to False
             DO jk = 1, jpkm1
                IF (ln_temnofreeze) THEN
                   ! Do not apply negative increments if the temperature will fall below freezing
-                  WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. &
-                     &   pts(:,:,jk,jp_tem,Kmm) + pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 
-                     pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt  
+                  WHERE(t_bkginc(ST_2D(0),jk) > 0.0_wp .OR. &
+                     &   pts(ST_2D(0),jk,jp_tem,Kmm) + pts(ST_2D(0),jk,jp_tem,Krhs) + t_bkginc(ST_2D(0),jk) * wgtiau(it) > fzptnz(:,:,jk) )
+                     pts(ST_2D(0),jk,jp_tem,Krhs) = pts(ST_2D(0),jk,jp_tem,Krhs) + t_bkginc(ST_2D(0),jk) * zincwgt
                   END WHERE
                ELSE
-                  pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt  
+                  DO_2D( 0, 0, 0, 0 )
+                     pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + t_bkginc(ji,jj,jk) * zincwgt
+                  END_2D
                ENDIF
                IF (ln_salfix) THEN
                   ! Do not apply negative increments if the salinity will fall below a specified
                   ! minimum value salfixmin
-                  WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. &
-                     &   pts(:,:,jk,jp_sal,Kmm) + pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 
-                     pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt
+                  WHERE(s_bkginc(ST_2D(0),jk) > 0.0_wp .OR. &
+                     &   pts(ST_2D(0),jk,jp_sal,Kmm) + pts(ST_2D(0),jk,jp_sal,Krhs) + s_bkginc(ST_2D(0),jk) * wgtiau(it) > salfixmin )
+                     pts(ST_2D(0),jk,jp_sal,Krhs) = pts(ST_2D(0),jk,jp_sal,Krhs) + s_bkginc(ST_2D(0),jk) * zincwgt
                   END WHERE
                ELSE
-                  pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt
+                  DO_2D( 0, 0, 0, 0 )
+                     pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs) + s_bkginc(ji,jj,jk) * zincwgt
+                  END_2D
                ENDIF
             END DO
@@ -569,7 +579,9 @@
          ENDIF
          !
-         IF ( kt == nitiaufin_r + 1  ) THEN   ! For bias crcn to work
-            DEALLOCATE( t_bkginc )
-            DEALLOCATE( s_bkginc )
+         IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile
+            IF ( kt == nitiaufin_r + 1  ) THEN   ! For bias crcn to work
+               DEALLOCATE( t_bkginc )
+               DEALLOCATE( s_bkginc )
+            ENDIF
          ENDIF
          !                             !--------------------------------------
@@ -582,23 +594,30 @@
             !
             ! Initialize the now fields with the background + increment
+            ! TODO: NOT TESTED- logical is forced to False
             IF (ln_temnofreeze) THEN
                ! Do not apply negative increments if the temperature will fall below freezing
-               WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_tem,Kmm) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 
-                  pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
+               WHERE( t_bkginc(ST_2D(0),:) > 0.0_wp .OR. pts(ST_2D(0),:,jp_tem,Kmm) + t_bkginc(ST_2D(0),:) > fzptnz(:,:,:) )
+                  pts(ST_2D(0),:,jp_tem,Kmm) = t_bkg(ST_2D(0),:) + t_bkginc(ST_2D(0),:)
                END WHERE
             ELSE
-               pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
+               DO_3D( 0, 0, 0, 0, 1, jpk )
+                  pts(ji,jj,jk,jp_tem,Kmm) = t_bkg(ji,jj,jk) + t_bkginc(ji,jj,jk)
+               END_3D
             ENDIF
             IF (ln_salfix) THEN
                ! Do not apply negative increments if the salinity will fall below a specified
                ! minimum value salfixmin
-               WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_sal,Kmm) + s_bkginc(:,:,:) > salfixmin ) 
-                  pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
+               WHERE( s_bkginc(ST_2D(0),:) > 0.0_wp .OR. pts(ST_2D(0),:,jp_sal,Kmm) + s_bkginc(ST_2D(0),:) > salfixmin )
+                  pts(ST_2D(0),:,jp_sal,Kmm) = s_bkg(ST_2D(0),:) + s_bkginc(ST_2D(0),:)
                END WHERE
             ELSE
-               pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
-            ENDIF
-
-            pts(:,:,:,:,Kbb) = pts(:,:,:,:,Kmm)                 ! Update before fields
+               DO_3D( 0, 0, 0, 0, 1, jpk )
+                  pts(ji,jj,jk,jp_sal,Kmm) = s_bkg(ji,jj,jk) + s_bkginc(ji,jj,jk)
+               END_3D
+            ENDIF
+
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               pts(ji,jj,jk,:,Kbb) = pts(ji,jj,jk,:,Kmm)             ! Update before fields
+            END_3D
 
             CALL eos( pts(:,:,:,:,Kbb), rhd, rhop, gdept_0(:,:,:) )  ! Before potential and in situ densities
@@ -607,18 +626,25 @@
 !!gm
 
-            IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)           &
-               &  CALL zps_hde    ( kt, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv,        &  ! Partial steps: before horizontal gradient
-               &                              rhd, gru , grv               )  ! of t, s, rd at the last ocean level
-            IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)                       &
-               &  CALL zps_hde_isf( nit000, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, gtui, gtvi,    &  ! Partial steps for top cell (ISF)
-               &                                  rhd, gru , grv , grui, grvi          )  ! of t, s, rd at the last ocean level
-
-            DEALLOCATE( t_bkginc )
-            DEALLOCATE( s_bkginc )
-            DEALLOCATE( t_bkg    )
-            DEALLOCATE( s_bkg    )
+            ! TEMP: This change not necessary after extra haloes development (lbc_lnk removed from zps_hde*)
+            IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
+               IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)           &
+                  &  CALL zps_hde    ( kt, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv,        &  ! Partial steps: before horizontal gradient
+                  &                              rhd, gru , grv               )  ! of t, s, rd at the last ocean level
+               IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)                       &
+                  &  CALL zps_hde_isf( nit000, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, gtui, gtvi,    &  ! Partial steps for top cell (ISF)
+                  &                                  rhd, gru , grv , grui, grvi          )  ! of t, s, rd at the last ocean level
+            ENDIF
+
+            IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile
+               DEALLOCATE( t_bkginc )
+               DEALLOCATE( s_bkginc )
+               DEALLOCATE( t_bkg    )
+               DEALLOCATE( s_bkg    )
+            ENDIF
+         !
          ENDIF
          !  
       ENDIF
+      ! TODO: NOT TESTED- logical is forced to False
       ! Perhaps the following call should be in step
       IF ( ln_seaiceinc  )   CALL seaice_asm_inc ( kt )   ! apply sea ice concentration increment
@@ -829,8 +855,9 @@
       INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation
       !
+      INTEGER  ::   ji, jj
       INTEGER  ::   it
       REAL(wp) ::   zincwgt   ! IAU weight for current time step
 #if defined key_si3
-      REAL(wp), DIMENSION(jpi,jpj) ::   zofrld, zohicif, zseaicendg, zhicifinc
+      REAL(wp), DIMENSION(ST_2D(nn_hls)) ::   zofrld, zohicif, zseaicendg, zhicifinc
       REAL(wp) ::   zhicifmin = 0.5_wp      ! ice minimum depth in metres
 #endif
@@ -847,8 +874,10 @@
             ! note this is not a tendency so should not be divided by rn_Dt (as with the tracer and other increments)
             !
-            IF(lwp) THEN
-               WRITE(numout,*) 
-               WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
-               WRITE(numout,*) '~~~~~~~~~~~~'
+            IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+               IF(lwp) THEN
+                  WRITE(numout,*)
+                  WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
+                  WRITE(numout,*) '~~~~~~~~~~~~'
+               ENDIF
             ENDIF
             !
@@ -856,16 +885,18 @@
             !
 #if defined key_si3
-            zofrld (:,:) = 1._wp - at_i(:,:)
-            zohicif(:,:) = hm_i(:,:)
-            !
-            at_i  (:,:) = 1. - MIN( MAX( 1.-at_i  (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
-            at_i_b(:,:) = 1. - MIN( MAX( 1.-at_i_b(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
-            fr_i(:,:) = at_i(:,:)        ! adjust ice fraction
-            !
-            zseaicendg(:,:) = zofrld(:,:) - (1. - at_i(:,:))   ! find out actual sea ice nudge applied
+            DO_2D( 0, 0, 0, 0 )
+               zofrld (ji,jj) = 1._wp - at_i(ji,jj)
+               zohicif(ji,jj) = hm_i(ji,jj)
+               !
+               at_i  (ji,jj) = 1. - MIN( MAX( 1.-at_i  (ji,jj) - seaice_bkginc(ji,jj) * zincwgt, 0.0_wp), 1.0_wp)
+               at_i_b(ji,jj) = 1. - MIN( MAX( 1.-at_i_b(ji,jj) - seaice_bkginc(ji,jj) * zincwgt, 0.0_wp), 1.0_wp)
+               fr_i(ji,jj) = at_i(ji,jj)        ! adjust ice fraction
+               !
+               zseaicendg(ji,jj) = zofrld(ji,jj) - (1. - at_i(ji,jj))   ! find out actual sea ice nudge applied
+            END_2D
             !
             ! Nudge sea ice depth to bring it up to a required minimum depth
-            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin ) 
-               zhicifinc(:,:) = (zhicifmin - hm_i(:,:)) * zincwgt    
+            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(ST_2D(0)) < zhicifmin )
+               zhicifinc(:,:) = (zhicifmin - hm_i(ST_2D(0))) * zincwgt
             ELSEWHERE
                zhicifinc(:,:) = 0.0_wp
@@ -873,5 +904,7 @@
             !
             ! nudge ice depth
-            hm_i (:,:) = hm_i (:,:) + zhicifinc(:,:)
+            DO_2D( 0, 0, 0, 0 )
+               hm_i (ji,jj) = hm_i (ji,jj) + zhicifinc(ji,jj)
+            END_2D
             !
             ! seaice salinity balancing (to add)
@@ -880,9 +913,13 @@
 #if defined key_cice && defined key_asminc
             ! Sea-ice : CICE case. Pass ice increment tendency into CICE
-            ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rn_Dt
-#endif
-            !
-            IF ( kt == nitiaufin_r ) THEN
-               DEALLOCATE( seaice_bkginc )
+            DO_2D( 0, 0, 0, 0 )
+               ndaice_da(ji,jj) = seaice_bkginc(ji,jj) * zincwgt / rn_Dt
+            END_2D
+#endif
+            !
+            IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile
+               IF ( kt == nitiaufin_r ) THEN
+                  DEALLOCATE( seaice_bkginc )
+               ENDIF
             ENDIF
             !
@@ -890,5 +927,7 @@
             !
 #if defined key_cice && defined key_asminc
-            ndaice_da(:,:) = 0._wp        ! Sea-ice : CICE case. Zero ice increment tendency into CICE
+            DO_2D( 0, 0, 0, 0 )
+               ndaice_da(ji,jj) = 0._wp        ! Sea-ice : CICE case. Zero ice increment tendency into CICE
+            END_2D
 #endif
             !
@@ -905,17 +944,19 @@
             !
 #if defined key_si3
-            zofrld (:,:) = 1._wp - at_i(:,:)
-            zohicif(:,:) = hm_i(:,:)
-            ! 
-            ! Initialize the now fields the background + increment
-            at_i(:,:) = 1. - MIN( MAX( 1.-at_i(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp)
-            at_i_b(:,:) = at_i(:,:) 
-            fr_i(:,:) = at_i(:,:)        ! adjust ice fraction
-            !
-            zseaicendg(:,:) = zofrld(:,:) - (1. - at_i(:,:))   ! find out actual sea ice nudge applied
+            DO_2D( 0, 0, 0, 0 )
+               zofrld (ji,jj) = 1._wp - at_i(ji,jj)
+               zohicif(ji,jj) = hm_i(ji,jj)
+               !
+               ! Initialize the now fields the background + increment
+               at_i(ji,jj) = 1. - MIN( MAX( 1.-at_i(ji,jj) - seaice_bkginc(ji,jj), 0.0_wp), 1.0_wp)
+               at_i_b(ji,jj) = at_i(ji,jj)
+               fr_i(ji,jj) = at_i(ji,jj)        ! adjust ice fraction
+               !
+               zseaicendg(ji,jj) = zofrld(ji,jj) - (1. - at_i(ji,jj))   ! find out actual sea ice nudge applied
+            END_2D
             !
             ! Nudge sea ice depth to bring it up to a required minimum depth
-            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin ) 
-               zhicifinc(:,:) = zhicifmin - hm_i(:,:)
+            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(ST_2D(0)) < zhicifmin )
+               zhicifinc(:,:) = zhicifmin - hm_i(ST_2D(0))
             ELSEWHERE
                zhicifinc(:,:) = 0.0_wp
@@ -923,5 +964,7 @@
             !
             ! nudge ice depth
-            hm_i (:,:) = hm_i (:,:) + zhicifinc(:,:)
+            DO_2D( 0, 0, 0, 0 )
+               hm_i(ji,jj) = hm_i (ji,jj) + zhicifinc(ji,jj)
+            END_2D
             !
             ! seaice salinity balancing (to add)
@@ -930,14 +973,20 @@
 #if defined key_cice && defined key_asminc
             ! Sea-ice : CICE case. Pass ice increment tendency into CICE
-           ndaice_da(:,:) = seaice_bkginc(:,:) / rn_Dt
-#endif
-            IF ( .NOT. PRESENT(kindic) ) THEN
-               DEALLOCATE( seaice_bkginc )
-            END IF
+            DO_2D( 0, 0, 0, 0 )
+               ndaice_da(ji,jj) = seaice_bkginc(ji,jj) / rn_Dt
+            END_2D
+#endif
+            IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile
+               IF ( .NOT. PRESENT(kindic) ) THEN
+                  DEALLOCATE( seaice_bkginc )
+               END IF
+            ENDIF
             !
          ELSE
             !
 #if defined key_cice && defined key_asminc
-            ndaice_da(:,:) = 0._wp     ! Sea-ice : CICE case. Zero ice increment tendency into CICE
+            DO_2D( 0, 0, 0, 0 )
+               ndaice_da(ji,jj) = 0._wp     ! Sea-ice : CICE case. Zero ice increment tendency into CICE
+            END_2D
 #endif
             !
Index: /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/BDY/bdytra.F90
===================================================================
--- /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/BDY/bdytra.F90	(revision 13517)
+++ /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/BDY/bdytra.F90	(revision 13518)
@@ -13,5 +13,5 @@
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers variables
-   USE dom_oce        ! ocean space and time domain variables 
+   USE dom_oce        ! ocean space and time domain variables
    USE bdy_oce        ! ocean open boundary conditions
    USE bdylib         ! for orlanski library routines
@@ -157,4 +157,8 @@
       INTEGER  ::   ib_bdy         ! Loop index
       !!----------------------------------------------------------------------
+      ! TODO: TO BE TILED
+      ! TODO: NOT TESTED- requires bdy
+      ! NOTE: Tiling these BDY loops is nontrivial; IF statements to check whether a point is in the current tile won't work (will be for every ib, every tile). The idx_bdy structure might require modifying to include a %nblen and list of ib indices for the current tile.
+      IF( ntile /= 0 .AND. ntile /= 1 ) RETURN                        ! Do only for the full domain
       !
       IF( ln_timing )   CALL timing_start('bdy_tra_dmp')
Index: /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/DOM/dtatsd.F90
===================================================================
--- /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/DOM/dtatsd.F90	(revision 13517)
+++ /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/DOM/dtatsd.F90	(revision 13518)
@@ -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,6 +136,6 @@
       !! ** 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(ST_2D(nn_hls),jpk,jpts), INTENT(  out) ::   ptsd   ! T & S data
       !
       INTEGER ::   ji, jj, jk, jl, jkk   ! dummy loop indicies
@@ -144,44 +145,53 @@
       !!----------------------------------------------------------------------
       !
-      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
+         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   ==!
+         ! TODO: NOT TESTED- requires orca2
+         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 = 1 )            ! Revert to tile domain
+      ENDIF
+      !
+      DO_3D( 0, 0, 0, 0, 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'
+         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
          !
@@ -215,7 +225,10 @@
       ELSE                                !==   z- or zps- coordinate   ==!
          !                             
-         ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) * tmask(:,:,:)    ! Mask
-         ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) * tmask(:,:,:)
-         !
+         DO_3D( 0, 0, 0, 0, 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
+         !
+         ! TODO: NOT TESTED- requires zps
          IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
             DO_2D( 1, 1, 1, 1 )
Index: /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/trabbc.F90
===================================================================
--- /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/trabbc.F90	(revision 13517)
+++ /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/trabbc.F90	(revision 13518)
@@ -17,4 +17,6 @@
    USE oce            ! ocean variables
    USE dom_oce        ! domain: ocean
+   ! TEMP: This change not necessary after trd_tra is tiled
+   USE domain, ONLY : dom_tile
    USE phycst         ! physical constants
    USE trd_oce        ! trends: ocean variables
@@ -80,13 +82,20 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts        ! active tracers and RHS of tracer equation
       !
-      INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt   ! 3D workspace
+      INTEGER  ::   ji, jj, jk    ! dummy loop indices
+      ! TEMP: This change not necessary after trd_tra is tiled
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ztrdt   ! 3D workspace
       !!----------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('tra_bbc')
       !
-      IF( l_trdtra )   THEN         ! Save the input temperature trend
-         ALLOCATE( ztrdt(jpi,jpj,jpk) )
-         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
+      IF( l_trdtra ) THEN           ! Save the input temperature trend
+         IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+            ! TEMP: This can be ST_2D(nn_hls) after trd_tra is tiled
+            ALLOCATE( ztrdt(jpi,jpj,jpk) )
+         ENDIF
+
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            ztrdt(ji,jj,jk) = pts(ji,jj,jk,jp_tem,Krhs)
+         END_3D
       ENDIF
       !                             !  Add the geothermal trend on temperature
@@ -96,14 +105,28 @@
       END_2D
       !
-      CALL lbc_lnk( 'trabbc', pts(:,:,:,jp_tem,Krhs) , 'T', 1.0_wp )
-      !
-      IF( l_trdtra ) THEN        ! Send the trend for diagnostics
-         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
-         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_bbc, ztrdt )
-         DEALLOCATE( ztrdt )
-      ENDIF
-      !
-      CALL iom_put ( "hfgeou" , rho0_rcp * qgh_trd0(:,:) )
-      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbc  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
+      ! TEMP: These changes not necessary after trd_tra is tiled, lbc_lnk not necessary if using XIOS (subdomain support, will not output haloes)
+      IF( l_trdtra ) THEN
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            ztrdt(ji,jj,jk) = pts(ji,jj,jk,jp_tem,Krhs) - ztrdt(ji,jj,jk)
+         END_3D
+      ENDIF
+
+      IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
+         !
+         IF( l_trdtra ) THEN        ! Send the trend for diagnostics
+            IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )         ! Use full domain
+
+            ! TODO: TO BE TILED- trd_tra
+            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_bbc, ztrdt )
+            DEALLOCATE( ztrdt )
+
+            IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = nijtile )   ! Revert to tile domain
+         ENDIF
+         !
+         CALL iom_put ( "hfgeou" , rho0_rcp * qgh_trd0(:,:) )
+      ENDIF
+
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbc  - Ta: ', mask1=tmask, &
+         &                                  clinfo3='tra-ta' )
       !
       IF( ln_timing )   CALL timing_stop('tra_bbc')
Index: /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/trabbl.F90
===================================================================
--- /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/trabbl.F90	(revision 13517)
+++ /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/trabbl.F90	(revision 13518)
@@ -26,4 +26,6 @@
    USE oce            ! ocean dynamics and active tracers
    USE dom_oce        ! ocean space and time domain
+   ! TEMP: This change not necessary after trd_tra is tiled
+   USE domain, ONLY : dom_tile
    USE phycst         ! physical constant
    USE eosbn2         ! equation of state
@@ -106,5 +108,7 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts             ! active tracers and RHS of tracer equation
       !
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
+      INTEGER  ::   ji, jj, jk   ! Dummy loop indices
+      ! TEMP: This change not necessary after trd_tra is tiled
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ztrdt, ztrds
       !!----------------------------------------------------------------------
       !
@@ -112,7 +116,13 @@
       !
       IF( l_trdtra )   THEN                         !* Save the T-S input trends
-         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
-         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
-         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
+         IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+            ! TEMP: This can be ST_2D(nn_hls) after trd_tra is tiled
+            ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) )
+         ENDIF
+
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            ztrdt(ji,jj,jk) = pts(ji,jj,jk,jp_tem,Krhs)
+            ztrds(ji,jj,jk) = pts(ji,jj,jk,jp_sal,Krhs)
+         END_3D
       ENDIF
 
@@ -125,8 +135,10 @@
          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
+            ! 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
+         ENDIF
          !
       ENDIF
@@ -136,19 +148,32 @@
          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
-         !
-      ENDIF
-
+         IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
+            ! 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
+         ENDIF
+         !
+      ENDIF
+
+      ! TEMP: These changes not necessary after trd_tra is tiled
       IF( l_trdtra )   THEN                      ! send the trends for further diagnostics
-         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
-         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
-         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_bbl, ztrdt )
-         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_bbl, ztrds )
-         DEALLOCATE( ztrdt, ztrds )
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            ztrdt(ji,jj,jk) = pts(ji,jj,jk,jp_tem,Krhs) - ztrdt(ji,jj,jk)
+            ztrds(ji,jj,jk) = pts(ji,jj,jk,jp_sal,Krhs) - ztrds(ji,jj,jk)
+         END_3D
+
+         IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
+            IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )         ! Use full domain
+
+            ! TODO: TO BE TILED- trd_tra
+            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_bbl, ztrdt )
+            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_bbl, ztrds )
+            DEALLOCATE( ztrdt, ztrds )
+
+            IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = nijtile )   ! Revert to tile domain
+         ENDIF
       ENDIF
       !
@@ -187,5 +212,5 @@
       INTEGER  ::   ik           ! local integers
       REAL(wp) ::   zbtr         ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::   zptb   ! workspace
+      REAL(wp), DIMENSION(ST_2D(nn_hls)) ::   zptb   ! workspace
       !!----------------------------------------------------------------------
       !
@@ -242,55 +267,52 @@
       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
+         DO_2D( 1, 0, 1, 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 +355,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(ST_2D(nn_hls),jpts)   :: zts, zab         ! 3D workspace
+      REAL(wp), DIMENSION(ST_2D(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)
Index: /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/tradmp.F90
===================================================================
--- /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/tradmp.F90	(revision 13517)
+++ /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/tradmp.F90	(revision 13518)
@@ -24,4 +24,6 @@
    USE oce            ! ocean: variables
    USE dom_oce        ! ocean: domain variables
+   ! TEMP: This change not necessary after trd_tra is tiled
+   USE domain, ONLY : dom_tile
    USE c1d            ! 1D vertical configuration
    USE trd_oce        ! trends: ocean variables
@@ -95,6 +97,7 @@
       !
       INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts)     ::  zts_dta
-      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  ztrdts
+      REAL(wp), DIMENSION(ST_2D(nn_hls),jpk,jpts)     ::  zts_dta
+      ! TEMP: This change not necessary after trd_tra is tiled
+      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, SAVE ::  ztrdts
       !!----------------------------------------------------------------------
       !
@@ -102,6 +105,12 @@
       !
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
-         ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) ) 
-         ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs) 
+         IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+            ! TEMP: This can be ST_2D(nn_hls) after trd_tra is tiled
+            ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) )
+         ENDIF
+
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            ztrdts(ji,jj,jk,:) = pts(ji,jj,jk,:,Krhs)
+         END_3D
       ENDIF
       !                           !==  input T-S data at kt  ==!
@@ -140,9 +149,20 @@
       END SELECT
       !
+      ! TEMP: These changes not necessary after trd_tra is tiled
       IF( l_trdtra )   THEN       ! trend diagnostic
-         ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs) - ztrdts(:,:,:,:)
-         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) )
-         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) )
-         DEALLOCATE( ztrdts ) 
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            ztrdts(ji,jj,jk,:) = pts(ji,jj,jk,:,Krhs) - ztrdts(ji,jj,jk,:)
+         END_3D
+
+         IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
+            IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )         ! Use full domain
+
+            ! TODO: TO BE TILED- trd_tra
+            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) )
+            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) )
+            DEALLOCATE( ztrdts )
+
+            IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = nijtile )   ! Revert to tile domain
+         ENDIF
       ENDIF
       !                           ! Control print
Index: /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/traisf.F90
===================================================================
--- /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/traisf.F90	(revision 13517)
+++ /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/traisf.F90	(revision 13518)
@@ -11,4 +11,5 @@
    !!----------------------------------------------------------------------
    USE isf_oce                                     ! Ice shelf variables
+   USE par_oce , ONLY : nijtile, ntile, ntsi, ntei, ntsj, ntej
    USE dom_oce                                     ! ocean space domain variables
    USE isfutils, ONLY : debug                      ! debug option
@@ -31,4 +32,5 @@
 CONTAINS
 
+   ! TODO: NOT TESTED- requires isf
    SUBROUTINE tra_isf ( kt, Kmm, pts, Krhs )
       !!----------------------------------------------------------------------
@@ -46,8 +48,10 @@
       IF( ln_timing )   CALL timing_start('tra_isf')
       !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_isf : Ice shelf heat fluxes'
-         IF(lwp) WRITE(numout,*) '~~~~~~~ '
+      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+         IF( kt == nit000 ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'tra_isf : Ice shelf heat fluxes'
+            IF(lwp) WRITE(numout,*) '~~~~~~~ '
+         ENDIF
       ENDIF
       !
@@ -76,6 +80,8 @@
       !
       IF ( ln_isfdebug ) THEN
-         CALL debug('tra_isf: pts(:,:,:,:,Krhs) T', pts(:,:,:,1,Krhs))
-         CALL debug('tra_isf: pts(:,:,:,:,Krhs) S', pts(:,:,:,2,Krhs))
+         IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only for the full domain
+            CALL debug('tra_isf: pts(:,:,:,:,Krhs) T', pts(:,:,:,1,Krhs))
+            CALL debug('tra_isf: pts(:,:,:,:,Krhs) S', pts(:,:,:,2,Krhs))
+         ENDIF
       END IF
       !
@@ -84,4 +90,5 @@
    END SUBROUTINE tra_isf
    !
+   ! TODO: NOT TESTED- requires isf
    SUBROUTINE tra_isf_mlt(ktop, kbot, phtbl, pfrac, ptsc, ptsc_b, pts)
       !!----------------------------------------------------------------------
@@ -101,9 +108,11 @@
       INTEGER                      :: ji,jj,jk  ! loop index   
       INTEGER                      :: ikt, ikb  ! top and bottom level of the tbl
-      REAL(wp), DIMENSION(jpi,jpj) :: ztc       ! total ice shelf tracer trend
+      REAL(wp), DIMENSION(ST_2D(nn_hls))     :: ztc       ! total ice shelf tracer trend
       !!----------------------------------------------------------------------
       !
       ! compute 2d total trend due to isf
-      ztc(:,:) = 0.5_wp * ( ptsc(:,:,jp_tem) + ptsc_b(:,:,jp_tem) ) / phtbl(:,:)
+      DO_2D( 1, 1, 1, 1 )
+         ztc(ji,jj) = 0.5_wp * ( ptsc(ji,jj,jp_tem) + ptsc_b(ji,jj,jp_tem) ) / phtbl(ji,jj)
+      END_2D
       !
       ! update pts(:,:,:,:,Krhs)
@@ -125,4 +134,5 @@
    END SUBROUTINE tra_isf_mlt
    !
+   ! TODO: NOT TESTED- requires isf
    SUBROUTINE tra_isf_cpl( Kmm, ptsc, ptsa )
       !!----------------------------------------------------------------------
@@ -137,13 +147,11 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) :: ptsc
       !!----------------------------------------------------------------------
-      INTEGER :: jk
+      INTEGER :: ji, jj, jk
       !!----------------------------------------------------------------------
       !
-      DO jk = 1,jpk
-         ptsa(:,:,jk,jp_tem) =   &
-            &  ptsa(:,:,jk,jp_tem) + ptsc(:,:,jk,jp_tem) * r1_e1e2t(:,:) / e3t(:,:,jk,Kmm)
-         ptsa(:,:,jk,jp_sal) =   &
-            &  ptsa(:,:,jk,jp_sal) + ptsc(:,:,jk,jp_sal) * r1_e1e2t(:,:) / e3t(:,:,jk,Kmm)
-      END DO
+      DO_3D( 0, 0, 0, 0, 1, jpk )
+         ptsa(ji,jj,jk,jp_tem) = ptsa(ji,jj,jk,jp_tem) + ptsc(ji,jj,jk,jp_tem) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+         ptsa(ji,jj,jk,jp_sal) = ptsa(ji,jj,jk,jp_sal) + ptsc(ji,jj,jk,jp_sal) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+      END_3D
       !
    END SUBROUTINE tra_isf_cpl
Index: /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/traqsr.F90
===================================================================
--- /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/traqsr.F90	(revision 13517)
+++ /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/traqsr.F90	(revision 13518)
@@ -22,4 +22,5 @@
    USE phycst         ! physical constants
    USE dom_oce        ! ocean space and time domain
+   USE domain, ONLY : dom_tile
    USE sbc_oce        ! surface boundary condition: ocean
    USE trc_oce        ! share SMS/Ocean variables
@@ -114,4 +115,6 @@
       REAL(wp) ::   zCb, zCmax, zpsi, zpsimax, zrdpsi, zCze
       REAL(wp) ::   zlogc, zlogze, zlogCtot, zlogCze
+      ! TEMP: These changes not necessary after trd_tra is tiled
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ztrdt
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: ze0, ze1, ze2, ze3
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, zetot, ztmp3d
@@ -120,13 +123,21 @@
       IF( ln_timing )   CALL timing_start('tra_qsr')
       !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation'
-         IF(lwp) WRITE(numout,*) '~~~~~~~'
+      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+         IF( kt == nit000 ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation'
+            IF(lwp) WRITE(numout,*) '~~~~~~~'
+         ENDIF
       ENDIF
       !
       IF( l_trdtra ) THEN      ! trends diagnostic: save the input temperature trend
-         ALLOCATE( ztrdt(jpi,jpj,jpk) ) 
-         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
+         IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+            ! TEMP: This can be ST_2D(nn_hls) after trd_tra is tiled
+            ALLOCATE( ztrdt(jpi,jpj,jpk) )
+         ENDIF
+
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            ztrdt(ji,jj,jk) = pts(ji,jj,jk,jp_tem,Krhs)
+         END_3D
       ENDIF
       !
@@ -136,14 +147,20 @@
       IF( kt == nit000 ) THEN          !==  1st time step  ==!
          IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0  .AND. .NOT.l_1st_euler ) THEN    ! read in restart
-            IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field read in the restart file'
             z1_2 = 0.5_wp
-            CALL iom_get( numror, jpdom_auto, 'qsr_hc_b', qsr_hc_b, ldxios = lrxios )   ! before heat content trend due to Qsr flux
+            IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+               IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field read in the restart file'
+               CALL iom_get( numror, jpdom_auto, 'qsr_hc_b', qsr_hc_b, ldxios = lrxios )   ! before heat content trend due to Qsr flux
+            ENDIF
          ELSE                                           ! No restart or restart not found: Euler forward time stepping
             z1_2 = 1._wp
-            qsr_hc_b(:,:,:) = 0._wp
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               qsr_hc_b(ji,jj,jk) = 0._wp
+            END_3D
          ENDIF
       ELSE                             !==  Swap of qsr heat content  ==!
          z1_2 = 0.5_wp
-         qsr_hc_b(:,:,:) = qsr_hc(:,:,:)
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            qsr_hc_b(ji,jj,jk) = qsr_hc(ji,jj,jk)
+         END_3D
       ENDIF
       !
@@ -154,16 +171,20 @@
       CASE( np_BIO )                   !==  bio-model fluxes  ==!
          !
-         DO jk = 1, nksr
-            qsr_hc(:,:,jk) = r1_rho0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
-         END DO
+         DO_3D( 0, 0, 0, 0, 1, nksr )
+            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) )
+         END_3D
          !
       CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==!
          !
-         ALLOCATE( ze0 (jpi,jpj)           , ze1 (jpi,jpj) ,   &
-            &      ze2 (jpi,jpj)           , ze3 (jpi,jpj) ,   &
-            &      ztmp3d(jpi,jpj,nksr + 1)                     )
+         ALLOCATE( ze0 (ST_2D(nn_hls))           , ze1 (ST_2D(nn_hls)) ,   &
+            &      ze2 (ST_2D(nn_hls))           , ze3 (ST_2D(nn_hls)) ,   &
+            &      ztmp3d(ST_2D(nn_hls),nksr + 1)                     )
          !
          IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll
-            CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step
+            IF( ntile == 0 .OR. ntile == 1 )  THEN                                         ! Do only for the full domain
+               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )            ! Use full domain
+               CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step
+               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 1 )            ! Revert to tile domain
+            ENDIF
             !
             ! Separation in R-G-B depending on the surface Chl
@@ -215,7 +236,7 @@
             ! Convert chlorophyll value to attenuation coefficient look-up table index
             zlui = 41 + 20.*LOG10(zchl) + 1.e-15
-            DO jk = 1, nksr + 1
-               ztmp3d(:,:,jk) = zlui 
-            END DO
+            DO_3D( 1, 1, 1, 1, 1, nksr + 1 )
+               ztmp3d(ji,jj,jk) = zlui
+            END_3D
          ENDIF
          !
@@ -277,27 +298,46 @@
          ENDIF
       END_2D
-      CALL lbc_lnk( 'traqsr', fraqsr_1lev(:,:), 'T', 1._wp )
-      !
-      IF( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution
-         ALLOCATE( zetot(jpi,jpj,jpk) )
-         zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero
-         DO jk = nksr, 1, -1
-            zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp
-         END DO         
-         CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation
-         DEALLOCATE( zetot ) 
-      ENDIF
-      !
-      IF( lrst_oce ) THEN     ! write in the ocean restart file
-         IF( lwxios ) CALL iom_swap(      cwxios_context          )
-         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc     , ldxios = lwxios )
-         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev, ldxios = lwxios ) 
-         IF( lwxios ) CALL iom_swap(      cxios_context          )
-      ENDIF
-      !
+      ! TEMP: This change not necessary after extra haloes development (lbc_lnk removed)
+      IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
+         CALL lbc_lnk( 'traqsr', fraqsr_1lev(:,:), 'T', 1._wp )
+      ENDIF
+      !
+      ! TEMP: This change not necessary and working array can use ST_2D(nn_hls) if using XIOS (subdomain support)
+      IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
+         IF( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution
+            ALLOCATE( zetot(jpi,jpj,jpk) )
+            zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero
+            DO jk = nksr, 1, -1
+               zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp
+            END DO
+            CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation
+            DEALLOCATE( zetot )
+         ENDIF
+      ENDIF
+      !
+      IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile
+         IF( lrst_oce ) THEN     ! write in the ocean restart file
+            IF( lwxios ) CALL iom_swap(      cwxios_context          )
+            CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc     , ldxios = lwxios )
+            CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev, ldxios = lwxios )
+            IF( lwxios ) CALL iom_swap(      cxios_context          )
+         ENDIF
+      ENDIF
+      !
+      ! TEMP: These changes not necessary after trd_tra is tiled
       IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
-         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
-         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt )
-         DEALLOCATE( ztrdt ) 
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            ztrdt(ji,jj,jk) = pts(ji,jj,jk,jp_tem,Krhs) - ztrdt(ji,jj,jk)
+         END_3D
+
+         IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
+            IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )         ! Use full domain
+
+            ! TODO: TO BE TILED- trd_tra
+            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt )
+            DEALLOCATE( ztrdt )
+
+            IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = nijtile )   ! Revert to tile domain
+         ENDIF
       ENDIF
       !                       ! print mean trends (used for debugging)
Index: /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/trasbc.F90
===================================================================
--- /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/trasbc.F90	(revision 13517)
+++ /NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/trasbc.F90	(revision 13518)
@@ -19,4 +19,6 @@
    USE sbc_oce        ! surface boundary condition: ocean
    USE dom_oce        ! ocean space domain variables
+   ! TEMP: This change not necessary after trd_tra is tiled
+   USE domain, ONLY : dom_tile
    USE phycst         ! physical constant
    USE eosbn2         ! Equation Of State
@@ -79,25 +81,36 @@
       INTEGER  ::   ikt, ikb                    ! local integers
       REAL(wp) ::   zfact, z1_e3t, zdep, ztim   ! local scalar
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdt, ztrds
+      ! TEMP: This change not necessary after trd_tra is tiled
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  ztrdt, ztrds
       !!----------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('tra_sbc')
       !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_sbc : TRAcer Surface Boundary Condition'
-         IF(lwp) WRITE(numout,*) '~~~~~~~ '
+      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+         IF( kt == nit000 ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'tra_sbc : TRAcer Surface Boundary Condition'
+            IF(lwp) WRITE(numout,*) '~~~~~~~ '
+         ENDIF
       ENDIF
       !
       IF( l_trdtra ) THEN                    !* Save ta and sa trends
-         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
-         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
-         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
+         IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+            ! TEMP: This can be ST_2D(nn_hls) after trd_tra is tiled
+            ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) )
+         ENDIF
+
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            ztrdt(ji,jj,jk) = pts(ji,jj,jk,jp_tem,Krhs)
+            ztrds(ji,jj,jk) = pts(ji,jj,jk,jp_sal,Krhs)
+         END_3D
       ENDIF
       !
 !!gm  This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist)
       IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
-         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
-         qsr(:,:) = 0._wp                     ! qsr set to zero
+         DO_2D( 0, 0, 0, 0 )
+            qns(ji,jj) = qns(ji,jj) + qsr(ji,jj)      ! total heat flux in qns
+            qsr(ji,jj) = 0._wp                        ! qsr set to zero
+         END_2D
       ENDIF
 
@@ -109,44 +122,53 @@
          IF( ln_rstart .AND.    &               ! Restart: read in restart file
               & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN
-            IF(lwp) WRITE(numout,*) '          nit000-1 sbc tracer content field read in the restart file'
             zfact = 0.5_wp
-            sbc_tsc(:,:,:) = 0._wp
-            CALL iom_get( numror, jpdom_auto, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem), ldxios = lrxios )   ! before heat content sbc trend
-            CALL iom_get( numror, jpdom_auto, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal), ldxios = lrxios )   ! before salt content sbc trend
+            IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+               IF(lwp) WRITE(numout,*) '          nit000-1 sbc tracer content field read in the restart file'
+               sbc_tsc(:,:,:) = 0._wp
+               CALL iom_get( numror, jpdom_auto, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem), ldxios = lrxios )   ! before heat content sbc trend
+               CALL iom_get( numror, jpdom_auto, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal), ldxios = lrxios )   ! before salt content sbc trend
+            ENDIF
          ELSE                                   ! No restart or restart not found: Euler forward time stepping
             zfact = 1._wp
-            sbc_tsc(:,:,:) = 0._wp
-            sbc_tsc_b(:,:,:) = 0._wp
+            DO_2D( 0, 0, 0, 0 )
+               sbc_tsc(ji,jj,:) = 0._wp
+               sbc_tsc_b(ji,jj,:) = 0._wp
+            END_2D
          ENDIF
       ELSE                                !* other time-steps: swap of forcing fields
          zfact = 0.5_wp
-         sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:)
+         DO_2D( 0, 0, 0, 0 )
+            sbc_tsc_b(ji,jj,:) = sbc_tsc(ji,jj,:)
+         END_2D
       ENDIF
       !                             !==  Now sbc tracer content fields  ==!
-      DO_2D( 0, 1, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )
          sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj)   ! non solar heat flux
          sbc_tsc(ji,jj,jp_sal) = r1_rho0     * sfx(ji,jj)   ! salt flux due to freezing/melting
       END_2D
       IF( ln_linssh ) THEN                !* linear free surface  
-         DO_2D( 0, 1, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )
             sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm)
             sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_sal,Kmm)
          END_2D
-         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) )
-         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * pts(:,:,1,jp_sal,Kmm) )
+         IF( ntile == 0 .OR. ntile == nijtile )  THEN             ! Do only on the last tile
+            IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) )
+            IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * pts(:,:,1,jp_sal,Kmm) )
+         ENDIF
       ENDIF
       !
       DO jn = 1, jpts               !==  update tracer trend  ==!
-         DO_2D( 0, 1, 0, 0 )
-            pts(ji,jj,1,jn,Krhs) = pts(ji,jj,1,jn,Krhs) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) )    &
-               &                                                / e3t(ji,jj,1,Kmm)
+         DO_2D( 0, 0, 0, 0 )
+            pts(ji,jj,1,jn,Krhs) = pts(ji,jj,1,jn,Krhs) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t(ji,jj,1,Kmm)
          END_2D
       END DO
       !                  
-      IF( lrst_oce ) THEN           !==  write sbc_tsc in the ocean restart file  ==!
-         IF( lwxios ) CALL iom_swap(      cwxios_context          )
-         CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem), ldxios = lwxios )
-         CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_tsc(:,:,jp_sal), ldxios = lwxios )
-         IF( lwxios ) CALL iom_swap(      cxios_context          )
+      IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile
+         IF( lrst_oce ) THEN           !==  write sbc_tsc in the ocean restart file  ==!
+            IF( lwxios ) CALL iom_swap(      cwxios_context          )
+            CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem), ldxios = lwxios )
+            CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_tsc(:,:,jp_sal), ldxios = lwxios )
+            IF( lwxios ) CALL iom_swap(      cxios_context          )
+         ENDIF
       ENDIF
       !
@@ -157,5 +179,5 @@
       IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff 
          zfact = 0.5_wp
-         DO_2D( 0, 1, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )
             IF( rnf(ji,jj) /= 0._wp ) THEN
                zdep = zfact / h_rnf(ji,jj)
@@ -170,6 +192,8 @@
       ENDIF
 
-      IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*pts(:,:,1,jp_tem,Kmm) )   ! runoff term on sst
-      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*pts(:,:,1,jp_sal,Kmm) )   ! runoff term on sss
+      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+         IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*pts(:,:,1,jp_tem,Kmm) )   ! runoff term on sst
+         IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*pts(:,:,1,jp_sal,Kmm) )   ! runoff term on sss
+      ENDIF
 
 #if defined key_asminc
@@ -182,5 +206,5 @@
       	 !
          IF( ln_linssh ) THEN 
-            DO_2D( 0, 1, 0, 0 )
+            DO_2D( 0, 0, 0, 0 )
                ztim = ssh_iau(ji,jj) / e3t(ji,jj,1,Kmm)
                pts(ji,jj,1,jp_tem,Krhs) = pts(ji,jj,1,jp_tem,Krhs) + pts(ji,jj,1,jp_tem,Kmm) * ztim
@@ -188,5 +212,5 @@
             END_2D
          ELSE
-            DO_2D( 0, 1, 0, 0 )
+            DO_2D( 0, 0, 0, 0 )
                ztim = ssh_iau(ji,jj) / ( ht(ji,jj) + 1. - ssmask(ji, jj) )
                pts(ji,jj,:,jp_tem,Krhs) = pts(ji,jj,:,jp_tem,Krhs) + pts(ji,jj,:,jp_tem,Kmm) * ztim
@@ -199,10 +223,21 @@
 #endif
       !
+      ! TEMP: These changes not necessary after trd_tra is tiled
       IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
-         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
-         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
-         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_nsr, ztrdt )
-         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_nsr, ztrds )
-         DEALLOCATE( ztrdt , ztrds ) 
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            ztrdt(ji,jj,jk) = pts(ji,jj,jk,jp_tem,Krhs) - ztrdt(ji,jj,jk)
+            ztrds(ji,jj,jk) = pts(ji,jj,jk,jp_sal,Krhs) - ztrds(ji,jj,jk)
+         END_3D
+
+         IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
+            IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )         ! Use full domain
+
+            ! TODO: TO BE TILED- trd_tra
+            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_nsr, ztrdt )
+            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_nsr, ztrds )
+            DEALLOCATE( ztrdt , ztrds )
+
+            IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = nijtile )   ! Revert to tile domain
+         ENDIF
       ENDIF
       !
