Index: NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/TRA/traldf_triad.F90
===================================================================
--- NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/TRA/traldf_triad.F90	(revision 13724)
+++ NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/TRA/traldf_triad.F90	(revision 14046)
@@ -13,4 +13,7 @@
    USE oce            ! ocean dynamics and active tracers
    USE dom_oce        ! ocean space and time domain
+   ! TEMP: [tiling] This change not necessary if XIOS has subdomain support
+   USE domain, ONLY : dom_tile
+   USE domutl, ONLY : is_tile
    USE phycst         ! physical constants
    USE trc_oce        ! share passive tracers/Ocean variables
@@ -33,6 +36,4 @@
    PUBLIC   tra_ldf_triad   ! routine called by traldf.F90
 
-   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE ::   zdkt3d   !: vertical tracer gradient at 2 levels
-
    LOGICAL  ::   l_ptr   ! flag to compute poleward transport
    LOGICAL  ::   l_hst   ! flag to compute heat transport
@@ -49,7 +50,30 @@
 CONTAINS
 
-  SUBROUTINE tra_ldf_triad( kt, Kmm, kit000, cdtype, pahu, pahv,               &
-      &                                              pgu , pgv  , pgui, pgvi , &
-      &                                         pt , pt2, pt_rhs, kjpt, kpass )
+   SUBROUTINE tra_ldf_triad( kt, Kmm, kit000, cdtype, pahu, pahv,             &
+      &                                               pgu , pgv , pgui, pgvi, &
+      &                                               pt, pt2, pt_rhs, kjpt, kpass )
+      !!
+      INTEGER                     , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                     , INTENT(in   ) ::   kit000     ! first time step index
+      CHARACTER(len=3)            , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
+      INTEGER                     , INTENT(in   ) ::   kjpt       ! number of tracers
+      INTEGER                     , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
+      INTEGER                     , INTENT(in   ) ::   Kmm        ! ocean time level indices
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pgu , pgv  ! tracer gradient at pstep levels
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pt         ! tracer (kpass=1) or laplacian of tracer (kpass=2)
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pt2        ! tracer (only used in kpass=2)
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pt_rhs     ! tracer trend
+      !!
+      CALL tra_ldf_triad_t( kt, Kmm, kit000, cdtype, pahu, pahv, is_tile(pahu),                            &
+      &                                              pgu , pgv , is_tile(pgu) , pgui, pgvi, is_tile(pgui), &
+      &                                              pt, is_tile(pt), pt2, is_tile(pt2), pt_rhs, is_tile(pt_rhs), kjpt, kpass )
+   END SUBROUTINE tra_ldf_triad
+
+
+  SUBROUTINE tra_ldf_triad_t( kt, Kmm, kit000, cdtype, pahu, pahv, ktah,                   &
+      &                                                pgu , pgv , ktg , pgui, pgvi, ktgi, &
+      &                                                pt, ktt, pt2, ktt2, pt_rhs, ktt_rhs, kjpt, kpass )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf_triad  ***
@@ -77,10 +101,11 @@
       INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
       INTEGER                              , INTENT(in)    ::   Kmm        ! ocean time level indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
-      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu , pgv  ! tracer gradient at pstep levels
-      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt         ! tracer (kpass=1) or laplacian of tracer (kpass=2)
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt2        ! tracer (only used in kpass=2)
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs     ! tracer trend
+      INTEGER                              , INTENT(in   ) ::   ktah, ktg, ktgi, ktt, ktt2, ktt_rhs
+      REAL(wp), DIMENSION(A2D_T(ktah),   JPK)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
+      REAL(wp), DIMENSION(A2D_T(ktg),        KJPT), INTENT(in   ) ::   pgu , pgv  ! tracer gradient at pstep levels
+      REAL(wp), DIMENSION(A2D_T(ktgi),       KJPT), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
+      REAL(wp), DIMENSION(A2D_T(ktt),    JPK,KJPT), INTENT(in   ) ::   pt         ! tracer (kpass=1) or laplacian of tracer (kpass=2)
+      REAL(wp), DIMENSION(A2D_T(ktt2),   JPK,KJPT), INTENT(in   ) ::   pt2        ! tracer (only used in kpass=2)
+      REAL(wp), DIMENSION(A2D_T(ktt_rhs),JPK,KJPT), INTENT(inout) ::   pt_rhs     ! tracer trend
       !
       INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
@@ -94,26 +119,25 @@
       REAL(wp) ::   ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt
       REAL(wp) ::   zah, zah_slp, zaei_slp
-      REAL(wp), DIMENSION(jpi,jpj    ) ::   z2d                                              ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw   ! 3D     -
+      REAL(wp), DIMENSION(A2D(nn_hls),0:1)     ::   zdkt3d                         ! vertical tracer gradient at 2 levels
+      REAL(wp), DIMENSION(A2D(nn_hls)        ) ::   z2d                            ! 2D workspace
+      REAL(wp), DIMENSION(A2D(nn_hls)    ,jpk) ::   zdit, zdjt, zftu, zftv, ztfw   ! 3D     -
+      ! TEMP: [tiling] This can be A2D(nn_hls) if XIOS has subdomain support
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpsi_uw, zpsi_vw
       !!----------------------------------------------------------------------
       !
-      IF( .NOT.ALLOCATED(zdkt3d) )  THEN
-         ALLOCATE( zdkt3d(jpi,jpj,0:1) , STAT=ierr )
-         CALL mpp_sum ( 'traldf_triad', ierr )
-         IF( ierr > 0 )   CALL ctl_stop('STOP', 'tra_ldf_triad: unable to allocate arrays')
-      ENDIF
-     !
-      IF( kpass == 1 .AND. kt == kit000 )  THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_ldf_triad : rotated laplacian diffusion operator on ', cdtype
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
-      ENDIF
-      !   
-      l_hst = .FALSE.
-      l_ptr = .FALSE.
-      IF( cdtype == 'TRA' ) THEN
-         IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf') )      l_ptr = .TRUE. 
-         IF( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR.                   &
-         &   iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  )   l_hst = .TRUE.
+      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+         IF( kpass == 1 .AND. kt == kit000 )  THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'tra_ldf_triad : rotated laplacian diffusion operator on ', cdtype
+            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
+         ENDIF
+         !
+         l_hst = .FALSE.
+         l_ptr = .FALSE.
+         IF( cdtype == 'TRA' ) THEN
+            IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf') )      l_ptr = .TRUE.
+            IF( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR.                   &
+            &   iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  )   l_hst = .TRUE.
+         ENDIF
       ENDIF
       !
@@ -128,27 +152,23 @@
       IF( kpass == 1 ) THEN         !==  first pass only  and whatever the tracer is  ==!
          !
-         akz     (:,:,:) = 0._wp      
-         ah_wslp2(:,:,:) = 0._wp
-         IF( ln_ldfeiv_dia ) THEN
-            zpsi_uw(:,:,:) = 0._wp
-            zpsi_vw(:,:,:) = 0._wp
-         ENDIF
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            akz     (ji,jj,jk) = 0._wp
+            ah_wslp2(ji,jj,jk) = 0._wp
+         END_3D
          !
          DO ip = 0, 1                            ! i-k triads
             DO kp = 0, 1
-               DO_3D( 1, 0, 1, 0, 1, jpkm1 )
-                  ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
-                  zbu   = e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
-                  zah   = 0.25_wp * pahu(ji,jj,jk)
-                  zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
+               DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+                  ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm)
+                  zbu   = e1e2u(ji-ip,jj) * e3u(ji-ip,jj,jk,Kmm)
+                  zah   = 0.25_wp * pahu(ji-ip,jj,jk)
+                  zslope_skew = triadi_g(ji,jj,jk,1-ip,kp)
                   ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth)
-                  zslope2 = zslope_skew + ( gdept(ji+1,jj,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
+                  zslope2 = zslope_skew + ( gdept(ji-ip+1,jj,jk,Kmm) - gdept(ji-ip,jj,jk,Kmm) ) * r1_e1u(ji-ip,jj) * umask(ji-ip,jj,jk+kp)
                   zslope2 = zslope2 *zslope2
-                  ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji+ip,jj) * zslope2
-                  akz     (ji+ip,jj,jk+kp) = akz     (ji+ip,jj,jk+kp) + zah * r1_e1u(ji,jj)       &
-                     &                                                      * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
+                  ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji,jj) * zslope2
+                  akz     (ji,jj,jk+kp) = akz     (ji,jj,jk+kp) + zah * r1_e1u(ji-ip,jj)       &
+                     &                                                      * r1_e1u(ji-ip,jj) * umask(ji-ip,jj,jk+kp)
                      !
-                 IF( ln_ldfeiv_dia )   zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp)   &
-                     &                                       + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * zslope_skew
                END_3D
             END DO
@@ -157,19 +177,17 @@
          DO jp = 0, 1                            ! j-k triads 
             DO kp = 0, 1
-               DO_3D( 1, 0, 1, 0, 1, jpkm1 )
-                  ze3wr = 1.0_wp / e3w(ji,jj+jp,jk+kp,Kmm)
-                  zbv   = e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
-                  zah   = 0.25_wp * pahv(ji,jj,jk)
-                  zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
+               DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+                  ze3wr = 1.0_wp / e3w(ji,jj,jk+kp,Kmm)
+                  zbv   = e1e2v(ji,jj-jp) * e3v(ji,jj-jp,jk,Kmm)
+                  zah   = 0.25_wp * pahv(ji,jj-jp,jk)
+                  zslope_skew = triadj_g(ji,jj,jk,1-jp,kp)
                   ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces
                   !    (do this by *adding* gradient of depth)
-                  zslope2 = zslope_skew + ( gdept(ji,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
+                  zslope2 = zslope_skew + ( gdept(ji,jj-jp+1,jk,Kmm) - gdept(ji,jj-jp,jk,Kmm) ) * r1_e2v(ji,jj-jp) * vmask(ji,jj-jp,jk+kp)
                   zslope2 = zslope2 * zslope2
-                  ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj+jp) * zslope2
-                  akz     (ji,jj+jp,jk+kp) = akz     (ji,jj+jp,jk+kp) + zah * r1_e2v(ji,jj)     &
-                     &                                                      * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
+                  ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj) * zslope2
+                  akz     (ji,jj,jk+kp) = akz     (ji,jj,jk+kp) + zah * r1_e2v(ji,jj-jp)     &
+                     &                                                      * r1_e2v(ji,jj-jp) * vmask(ji,jj-jp,jk+kp)
                   !
-                  IF( ln_ldfeiv_dia )   zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp)   &
-                     &                                       + 0.25 * aeiv(ji,jj,jk) * e1v(ji,jj) * zslope_skew
                END_3D
             END DO
@@ -179,5 +197,5 @@
             !
             IF( ln_traldf_blp ) THEN                ! bilaplacian operator
-               DO_3D( 1, 0, 1, 0, 2, jpkm1 )
+               DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                   akz(ji,jj,jk) = 16._wp           &
                      &   * ah_wslp2   (ji,jj,jk)   &
@@ -187,5 +205,5 @@
                END_3D
             ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
-               DO_3D( 1, 0, 1, 0, 2, jpkm1 )
+               DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                   ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
                   zcoef0 = rDt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
@@ -195,8 +213,32 @@
            !
          ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
-            akz(:,:,:) = ah_wslp2(:,:,:)      
-         ENDIF
-         !
-         IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' )   CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm )
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               akz(ji,jj,jk) = ah_wslp2(ji,jj,jk)
+            END_3D
+         ENDIF
+         !
+         ! TEMP: [tiling] These changes not necessary if XIOS has subdomain support
+         IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
+            IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) THEN
+               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )
+
+               zpsi_uw(:,:,:) = 0._wp
+               zpsi_vw(:,:,:) = 0._wp
+
+               DO jp = 0, 1
+                  DO kp = 0, 1
+                     DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+                        zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) &
+                           & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+jp,jj,jk,1-jp,kp)
+                        zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp) &
+                           & + 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj+jp,jk,1-jp,kp)
+                     END_3D
+                  END DO
+               END DO
+               CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm )
+
+               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = nijtile )
+            ENDIF
+         ENDIF
          !
       ENDIF                                  !==  end 1st pass only  ==!
@@ -234,9 +276,14 @@
          DO jk = 1, jpkm1
             !                    !==  Vertical tracer gradient at level jk and jk+1
-            zdkt3d(:,:,1) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
+            DO_2D( 1, 1, 1, 1 )
+               zdkt3d(ji,jj,1) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1)
+            END_2D
             !
             !                    ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1)
             IF( jk == 1 ) THEN   ;   zdkt3d(:,:,0) = zdkt3d(:,:,1)
-            ELSE                 ;   zdkt3d(:,:,0) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)
+            ELSE
+               DO_2D( 1, 1, 1, 1 )
+                  zdkt3d(ji,jj,0) = ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
+               END_2D
             ENDIF
             !
@@ -380,5 +427,5 @@
       END DO                                                      ! end tracer loop
       !                                                           ! ===============
-   END SUBROUTINE tra_ldf_triad
+   END SUBROUTINE tra_ldf_triad_t
 
    !!==============================================================================
