Changeset 14834 for NEMO/trunk/src/OCE/TRA

Timestamp:

2021-05-11T11:24:44+02:00 (3 years ago)

Author:

hadcv

Message:

#2600: Merge in dev_r14273_HPC-02_Daley_Tiling

Location:

NEMO/trunk/src/OCE/TRA

Files:

• eosbn2.F90 (modified) (4 diffs)
• traadv.F90 (modified) (15 diffs)
• traadv_cen.F90 (modified) (6 diffs)
• traadv_cen_lf.F90 (copied) (copied from NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traadv_cen_lf.F90)
• traadv_fct.F90 (modified) (11 diffs)
• traadv_fct_lf.F90 (deleted)
• traadv_mus.F90 (modified) (6 diffs)
• traadv_mus_lf.F90 (deleted)
• traadv_qck.F90 (modified) (13 diffs)
• traadv_qck_lf.F90 (copied) (copied from NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traadv_qck_lf.F90)
• traadv_ubs.F90 (modified) (5 diffs)
• traadv_ubs_lf.F90 (copied) (copied from NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traadv_ubs_lf.F90)
• trabbc.F90 (modified) (1 diff)
• trabbl.F90 (modified) (8 diffs)
• traisf.F90 (modified) (2 diffs)
• traldf.F90 (modified) (3 diffs)
• traldf_iso.F90 (modified) (21 diffs)
• traldf_lap_blp.F90 (modified) (7 diffs)
• traldf_triad.F90 (modified) (26 diffs)
• tramle.F90 (modified) (18 diffs)
• tranpc.F90 (modified) (4 diffs)
• traqsr.F90 (modified) (15 diffs)
• trasbc.F90 (modified) (8 diffs)
• trazdf.F90 (modified) (1 diff)
• zpshde.F90 (modified) (15 diffs)

NEMO/trunk/src/OCE/TRA/eosbn2.F90

@@ -577 +577 @@
 
    SUBROUTINE eos_insitu_pot_2d( pts, prhop )
+      !!
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
+      !                                                     ! 2 : salinity               [psu]
+      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   prhop  ! potential density (surface referenced)
+      !!
+      CALL eos_insitu_pot_2d_t( pts, is_tile(pts), prhop, is_tile(prhop) )
+   END SUBROUTINE eos_insitu_pot_2d
+
+
+   SUBROUTINE eos_insitu_pot_2d_t( pts, ktts, prhop, ktrhop )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE eos_insitu_pot  ***
@@ -589 +599 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
+      INTEGER                              , INTENT(in   ) ::   ktts, ktrhop
+      REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
       !                                                                ! 2 : salinity               [psu]
-      REAL(wp), DIMENSION(jpi,jpj     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
+      REAL(wp), DIMENSION(A2D_T(ktrhop)   ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
       !
       INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices
@@ -606 +617 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-            DO_2D( 1, 1, 1, 1 )
-               !
-               zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
-               zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
-               ztm = tmask(ji,jj,1)                                         ! tmask
-               !
-               zn0 = (((((EOS060*zt   &
-                  &   + EOS150*zs+EOS050)*zt   &
-                  &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
-                  &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
-                  &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
-                  &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
-                  &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
-                  !
-               !
-               prhop(ji,jj) = zn0 * ztm                           ! potential density referenced at the surface
-               !
-            END_2D
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+            !
+            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
+            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+            ztm = tmask(ji,jj,1)                                         ! tmask
+            !
+            zn0 = (((((EOS060*zt   &
+               &   + EOS150*zs+EOS050)*zt   &
+               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
+               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
+               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
+               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
+               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+               !
+            !
+            prhop(ji,jj) = zn0 * ztm                           ! potential density referenced at the surface
+            !
+         END_2D
 
       CASE( np_seos )                !==  simplified EOS  ==!
          !
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             zt  = pts  (ji,jj,jp_tem) - 10._wp
             zs  = pts  (ji,jj,jp_sal) - 35._wp
@@ -646 +657 @@
       IF( ln_timing )   CALL timing_stop('eos-pot')
       !
-   END SUBROUTINE eos_insitu_pot_2d
+   END SUBROUTINE eos_insitu_pot_2d_t
 
 

NEMO/trunk/src/OCE/TRA/traadv.F90

@@ -18 +18 @@
    USE oce            ! ocean dynamics and active tracers
    USE dom_oce        ! ocean space and time domain
-   ! TEMP: [tiling] This change not necessary after extended haloes development
+   ! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
    USE domtile
    USE domvvl         ! variable vertical scale factors
@@ -25 +25 @@
    USE traadv_cen     ! centered scheme            (tra_adv_cen  routine)
    USE traadv_fct     ! FCT      scheme            (tra_adv_fct  routine)
-   USE traadv_fct_lf  ! FCT      scheme            (tra_adv_fct  routine - loop fusion version)
    USE traadv_mus     ! MUSCL    scheme            (tra_adv_mus  routine)
-   USE traadv_mus_lf  ! MUSCL    scheme            (tra_adv_mus  routine - loop fusion version)
    USE traadv_ubs     ! UBS      scheme            (tra_adv_ubs  routine)
    USE traadv_qck     ! QUICKEST scheme            (tra_adv_qck  routine)
@@ -61 +59 @@
    LOGICAL ::   ln_traadv_qck    ! QUICKEST scheme flag
 
-   INTEGER ::   nadv             ! choice of the type of advection scheme
+   INTEGER, PUBLIC ::   nadv             ! choice of the type of advection scheme
    !                             ! associated indices:
-   INTEGER, PARAMETER ::   np_NO_adv  = 0   ! no T-S advection
-   INTEGER, PARAMETER ::   np_CEN     = 1   ! 2nd/4th order centered scheme
-   INTEGER, PARAMETER ::   np_FCT     = 2   ! 2nd/4th order Flux Corrected Transport scheme
-   INTEGER, PARAMETER ::   np_MUS     = 3   ! MUSCL scheme
-   INTEGER, PARAMETER ::   np_UBS     = 4   ! 3rd order Upstream Biased Scheme
-   INTEGER, PARAMETER ::   np_QCK     = 5   ! QUICK scheme
+   INTEGER, PARAMETER, PUBLIC ::   np_NO_adv  = 0   ! no T-S advection
+   INTEGER, PARAMETER, PUBLIC ::   np_CEN     = 1   ! 2nd/4th order centered scheme
+   INTEGER, PARAMETER, PUBLIC ::   np_FCT     = 2   ! 2nd/4th order Flux Corrected Transport scheme
+   INTEGER, PARAMETER, PUBLIC ::   np_MUS     = 3   ! MUSCL scheme
+   INTEGER, PARAMETER, PUBLIC ::   np_UBS     = 4   ! 3rd order Upstream Biased Scheme
+   INTEGER, PARAMETER, PUBLIC ::   np_QCK     = 5   ! QUICK scheme
 
    !! * Substitutions
@@ -93 +91 @@
       !
       INTEGER ::   ji, jj, jk   ! dummy loop index
-      ! TEMP: [tiling] This change not necessary and can be A2D(nn_hls) if using XIOS (subdomain support)
+      ! TEMP: [tiling] This change not necessary and can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zuu, zvv, zww   ! 3D workspace
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds
-      ! TEMP: [tiling] This change not necessary after extra haloes development
+      ! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
       LOGICAL :: lskip
       !!----------------------------------------------------------------------
@@ -104 +102 @@
       lskip = .FALSE.
 
-      ! TEMP: [tiling] These changes not necessary if using XIOS (subdomain support)
-      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+      ! TEMP: [tiling] These changes not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
          ALLOCATE( zuu(jpi,jpj,jpk), zvv(jpi,jpj,jpk), zww(jpi,jpj,jpk) )
       ENDIF
 
-      ! TEMP: [tiling] These changes not necessary after extra haloes development (lbc_lnk removed from tra_adv_*) and if XIOS has subdomain support (ldf_eiv_dia)
-      IF( nadv /= np_CEN .OR. (nadv == np_CEN .AND. nn_cen_h == 4) .OR. ln_ldfeiv_dia )  THEN
-         IF( ln_tile ) THEN
-            IF( ntile == 1 ) THEN
-               CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )
-            ELSE
-               lskip = .TRUE.
-            ENDIF
+      ! TEMP: [tiling] These changes not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
+      IF( ln_tile .AND. nadv == np_FCT )  THEN
+         IF( ntile == 1 ) THEN
+            CALL dom_tile_stop( ldhold=.TRUE. )
+         ELSE
+            lskip = .TRUE.
          ENDIF
       ENDIF
@@ -122 +118 @@
          !                                         !==  effective transport  ==!
          IF( ln_wave .AND. ln_sdw )  THEN
-            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+            DO_3D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )
                zuu(ji,jj,jk) = e2u  (ji,jj) * e3u(ji,jj,jk,Kmm) * ( uu(ji,jj,jk,Kmm) + usd(ji,jj,jk) )
                zvv(ji,jj,jk) = e1v  (ji,jj) * e3v(ji,jj,jk,Kmm) * ( vv(ji,jj,jk,Kmm) + vsd(ji,jj,jk) )
@@ -128 +124 @@
             END_3D
          ELSE
-            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+            DO_3D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )
                zuu(ji,jj,jk) = e2u  (ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm)               ! eulerian transport only
                zvv(ji,jj,jk) = e1v  (ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm)
@@ -136 +132 @@
          !
          IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN                                ! add z-tilde and/or vvl corrections
-            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+            DO_3D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )
                zuu(ji,jj,jk) = zuu(ji,jj,jk) + un_td(ji,jj,jk)
                zvv(ji,jj,jk) = zvv(ji,jj,jk) + vn_td(ji,jj,jk)
@@ -142 +138 @@
          ENDIF
          !
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
             zuu(ji,jj,jpk) = 0._wp                                                      ! no transport trough the bottom
             zvv(ji,jj,jpk) = 0._wp
@@ -148 +144 @@
          END_2D
          !
-         ! TEMP: [tiling] These changes not necessary if using XIOS (subdomain support)
          IF( ln_ldfeiv .AND. .NOT. ln_traldf_triad )   &
-            &              CALL ldf_eiv_trp( kt, nit000, zuu(A2D(nn_hls),:), zvv(A2D(nn_hls),:), zww(A2D(nn_hls),:), &
-            &                                'TRA', Kmm, Krhs )   ! add the eiv transport (if necessary)
-         !
-         IF( ln_mle    )   CALL tra_mle_trp( kt, nit000, zuu(A2D(nn_hls),:), zvv(A2D(nn_hls),:), zww(A2D(nn_hls),:), &
-            &                                'TRA', Kmm       )   ! add the mle transport (if necessary)
-         !
-         ! TEMP: [tiling] This change not necessary if using XIOS (subdomain support)
-         IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile
+            &              CALL ldf_eiv_trp( kt, nit000, zuu, zvv, zww, 'TRA', Kmm, Krhs )   ! add the eiv transport (if necessary)
+         !
+         IF( ln_mle    )   CALL tra_mle_trp( kt, nit000, zuu, zvv, zww, 'TRA', Kmm       )   ! add the mle transport (if necessary)
+         !
+         ! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
+         IF( .NOT. l_istiled .OR. ntile == nijtile )  THEN                ! Do only on the last tile
             CALL iom_put( "uocetr_eff", zuu )                                        ! output effective transport
             CALL iom_put( "vocetr_eff", zvv )
@@ -163 +156 @@
          ENDIF
          !
-   !!gm ???
-         ! TEMP: [tiling] This change not necessary if using XIOS (subdomain support)
+!!gm ???
+         ! TEMP: [tiling] This copy-in not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
          CALL dia_ptr( kt, Kmm, zvv(A2D(nn_hls),:) )                                    ! diagnose the effective MSF
-   !!gm ???
+!!gm ???
          !
 
@@ -178 +171 @@
          !
          CASE ( np_CEN )                                 ! Centered scheme : 2nd / 4th order
-            IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kmm), 'T', 1. )
             CALL tra_adv_cen    ( kt, nit000, 'TRA',         zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v )
          CASE ( np_FCT )                                 ! FCT scheme      : 2nd / 4th order
-            IF (nn_hls.EQ.2) THEN
-               CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1., pts(:,:,:,:,Kmm), 'T', 1.)
-               CALL lbc_lnk( 'traadv', zuu(:,:,:), 'U', -1., zvv(:,:,:), 'V', -1., zww(:,:,:), 'W', 1.)
-#if defined key_loop_fusion
-               CALL tra_adv_fct_lf ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v )
-#else
                CALL tra_adv_fct    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v )
-#endif
-            ELSE
-               CALL tra_adv_fct    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v )
-            END IF
          CASE ( np_MUS )                                 ! MUSCL
-            IF (nn_hls.EQ.2) THEN
-                CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)
-#if defined key_loop_fusion
-                CALL tra_adv_mus_lf ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )
-#else
                 CALL tra_adv_mus    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )
-#endif
-            ELSE
-                CALL tra_adv_mus    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )
-            END IF
          CASE ( np_UBS )                                 ! UBS
-            IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)
             CALL tra_adv_ubs    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v   )
          CASE ( np_QCK )                                 ! QUICKEST
-            IF (nn_hls.EQ.2) THEN
-               CALL lbc_lnk( 'traadv', zuu(:,:,:), 'U', -1., zvv(:,:,:), 'V', -1.)
-               CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)
-            END IF
             CALL tra_adv_qck    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs )
          !
@@ -225 +193 @@
          ENDIF
 
-         ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed from tra_adv_*) and if XIOS has subdomain support (ldf_eiv_dia)
-         IF( ln_tile .AND. ntile == 0 ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 1 )
-
+         ! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
+         IF( ln_tile .AND. .NOT. l_istiled ) CALL dom_tile_start( ldhold=.TRUE. )
       ENDIF
       !                                              ! print mean trends (used for debugging)
@@ -233 +200 @@
          &                                  tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
 
-      ! TEMP: [tiling] This change not necessary if using XIOS (subdomain support)
-      IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
+      ! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
+      IF( .NOT. l_istiled .OR. ntile == nijtile )  THEN                ! Do only for the full domain
          DEALLOCATE( zuu, zvv, zww )
       ENDIF
@@ -306 +273 @@
         CALL ctl_stop( 'tra_adv_init: FCT scheme, choose 2nd or 4th order' )
       ENDIF
+      ! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
+      IF( ln_traadv_fct .AND. ln_tile ) THEN
+         CALL ctl_warn( 'tra_adv_init: FCT scheme does not yet work with tiling' )
+      ENDIF
       IF( ln_traadv_ubs .AND. ( nn_ubs_v /= 2 .AND. nn_ubs_v /= 4 )   ) THEN     ! UBS
         CALL ctl_stop( 'tra_adv_init: UBS scheme, choose 2nd or 4th order' )

NEMO/trunk/src/OCE/TRA/traadv_cen.F90

@@ -23 +23 @@
    USE trc_oce        ! share passive tracers/Ocean variables
    USE lib_mpp        ! MPP library
+#if defined key_loop_fusion
+   USE traadv_cen_lf  ! centered scheme            (tra_adv_cen  routine - loop fusion version)
+#endif
 
    IMPLICIT NONE
@@ -71 +74 @@
       INTEGER                                  , INTENT(in   ) ::   kn_cen_h        ! =2/4 (2nd or 4th order scheme)
       INTEGER                                  , INTENT(in   ) ::   kn_cen_v        ! =2/4 (2nd or 4th order scheme)
-      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
+      ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
       REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
@@ -82 +85 @@
       !!----------------------------------------------------------------------
       !
-      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+#if defined key_loop_fusion
+      CALL tra_adv_cen_lf    ( kt, nit000, cdtype, pU, pV, pW, Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v )
+#else
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
          IF( kt == kit000 )  THEN
             IF(lwp) WRITE(numout,*)
@@ -119 +125 @@
                ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
             END_3D
-            IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_cen', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp )   ! Lateral boundary cond.
+            IF (nn_hls==1) CALL lbc_lnk( 'traadv_cen', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp, ld4only= .TRUE. )   ! Lateral boundary cond.
             !
             DO_3D( nn_hls-1, 0, nn_hls-1, 0, 1, jpkm1 )           ! Horizontal advective fluxes
@@ -131 +137 @@
                zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * zC4t_v
             END_3D
-            IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_cen', zwx, 'U', -1. , zwy, 'V', -1. )
+            IF (nn_hls==1) CALL lbc_lnk( 'traadv_cen', zwx, 'U', -1. , zwy, 'V', -1. )
             !
          CASE DEFAULT
@@ -184 +190 @@
       END DO
       !
+#endif
    END SUBROUTINE tra_adv_cen
 

NEMO/trunk/src/OCE/TRA/traadv_fct.F90

@@ -34 +34 @@
    PUBLIC   tra_adv_fct        ! called by traadv.F90
    PUBLIC   interp_4th_cpt     ! called by traadv_cen.F90
-   PUBLIC   tridia_solver      ! called by traadv_fct_lf.F90
-   PUBLIC   nonosc             ! called by traadv_fct_lf.F90 - key_agrif
 
    LOGICAL  ::   l_trd   ! flag to compute trends
@@ -81 +79 @@
       INTEGER                                  , INTENT(in   ) ::   kn_fct_v        ! order of the FCT scheme (=2 or 4)
       REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
-      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
+      ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case
       REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
@@ -95 +93 @@
       !!----------------------------------------------------------------------
       !
-      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+#if defined key_loop_fusion
+      CALL tra_adv_fct_lf ( kt, nit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs, kn_fct_h, kn_fct_v )
+#else
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
          IF( kt == kit000 )  THEN
             IF(lwp) WRITE(numout,*)
@@ -136 +137 @@
       ! If adaptive vertical advection, check if it is needed on this PE at this time
       IF( ln_zad_Aimp ) THEN
-         IF( MAXVAL( ABS( wi(A2D(nn_hls),:) ) ) > 0._wp ) ll_zAimp = .TRUE.
+         IF( MAXVAL( ABS( wi(A2D(1),:) ) ) > 0._wp ) ll_zAimp = .TRUE.
       END IF
       ! If active adaptive vertical advection, build tridiagonal matrix
@@ -239 +240 @@
             END DO
             ! NOTE [ comm_cleanup ] : need to change sign to ensure halo 1 - halo 2 compatibility
-            CALL lbc_lnk( 'traadv_fct', zltu, 'T', -1.0_wp , zltv, 'T', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
-            !
-            DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+            CALL lbc_lnk( 'traadv_fct', zltu, 'T', -1.0_wp , zltv, 'T', -1.0_wp, ld4only= .TRUE. )   ! Lateral boundary cond. (unchanged sgn)
+            !
+            DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )
                zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points
                zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
@@ -254 +255 @@
                              &                          ) - zwy(ji,jj,jk)
             END_3D
-            IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv_fct', zwx, 'U', -1.0_wp, zwy, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
             !
          CASE(  41 )                   !- 4th order centered       ==>>   !!gm coding attempt   need to be tested
             ztu(:,:,jpk) = 0._wp             ! Bottom value : flux set to zero
             ztv(:,:,jpk) = 0._wp
-            DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )    ! 1st derivative (gradient)
+            DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )    ! 1st derivative (gradient)
                ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
                ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
             END_3D
-            IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
+            IF (nn_hls==1) CALL lbc_lnk( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp, ld4only= .TRUE. )   ! Lateral boundary cond. (unchanged sgn)
             !
             DO_3D( 0, 0, 0, 0, 1, jpkm1 )    ! Horizontal advective fluxes
@@ -275 +275 @@
                zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk)
             END_3D
-            IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv_fct', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
+            IF (nn_hls==2) CALL lbc_lnk( 'traadv_fct', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
             !
          END SELECT
@@ -298 +298 @@
          ENDIF
          !
-         IF (nn_hls.EQ.1) THEN
+         IF (nn_hls==1) THEN
             CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp, zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp, zwz, 'T', 1.0_wp )
          ELSE
@@ -381 +381 @@
       ENDIF
       !
+#endif
    END SUBROUTINE tra_adv_fct
 
@@ -456 +457 @@
          END_2D
       END DO
-      IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_fct', zbetup, 'T', 1.0_wp , zbetdo, 'T', 1.0_wp )   ! lateral boundary cond. (unchanged sign)
+      IF (nn_hls==1) CALL lbc_lnk( 'traadv_fct', zbetup, 'T', 1.0_wp , zbetdo, 'T', 1.0_wp, ld4only= .TRUE. )   ! lateral boundary cond. (unchanged sign)
 
       ! 3. monotonic flux in the i & j direction (paa & pbb)
@@ -677 +678 @@
    END SUBROUTINE tridia_solver
 
+#if defined key_loop_fusion
+#define tracer_flux_i(out,zfp,zfm,ji,jj,jk) \
+        zfp = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) ) ; \
+        zfm = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) ) ; \
+        out = 0.5 * ( zfp * pt(ji,jj,jk,jn,Kbb) + zfm * pt(ji+1,jj,jk,jn,Kbb) )
+
+#define tracer_flux_j(out,zfp,zfm,ji,jj,jk) \
+        zfp = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) ) ; \
+        zfm = pV(ji,jj,jk) - ABS( pV(ji,jj,jk) ) ; \
+        out = 0.5 * ( zfp * pt(ji,jj,jk,jn,Kbb) + zfm * pt(ji,jj+1,jk,jn,Kbb) )
+
+   SUBROUTINE tra_adv_fct_lf( kt, kit000, cdtype, p2dt, pU, pV, pW,       &
+      &                    Kbb, Kmm, pt, kjpt, Krhs, kn_fct_h, kn_fct_v )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_adv_fct  ***
+      !!
+      !! **  Purpose :   Compute the now trend due to total advection of tracers
+      !!               and add it to the general trend of tracer equations
+      !!
+      !! **  Method  : - 2nd or 4th FCT scheme on the horizontal direction
+      !!               (choice through the value of kn_fct)
+      !!               - on the vertical the 4th order is a compact scheme
+      !!               - corrected flux (monotonic correction)
+      !!
+      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
+      !!             - send trends to trdtra module for further diagnostics (l_trdtra=T)
+      !!             - poleward advective heat and salt transport (ln_diaptr=T)
+      !!----------------------------------------------------------------------
+      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
+      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   ) ::   kn_fct_h        ! order of the FCT scheme (=2 or 4)
+      INTEGER                                  , INTENT(in   ) ::   kn_fct_v        ! order of the FCT scheme (=2 or 4)
+      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
+      !
+      INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices
+      REAL(wp) ::   ztra                                     ! local scalar
+      REAL(wp) ::   zwx_im1, zfp_ui, zfp_ui_m1, zfp_vj, zfp_vj_m1, zfp_wk, zC2t_u, zC4t_u   !   -      -
+      REAL(wp) ::   zwy_jm1, zfm_ui, zfm_ui_m1, zfm_vj, zfm_vj_m1, zfm_wk, zC2t_v, zC4t_v   !   -      -
+      REAL(wp) ::   ztu, ztv, ztu_im1, ztu_ip1, ztv_jm1, ztv_jp1
+      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwx_3d, zwy_3d, zwz, ztw, zltu_3d, zltv_3d
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdx, ztrdy, ztrdz, zptry
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zwinf, zwdia, zwsup
+      LOGICAL  ::   ll_zAimp                                 ! flag to apply adaptive implicit vertical advection
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == kit000 )  THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_adv_fct_lf : FCT advection scheme on ', cdtype
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      ENDIF
+      !! -- init to 0
+      zwx_3d(:,:,:) = 0._wp
+      zwy_3d(:,:,:) = 0._wp
+      zwz(:,:,:) = 0._wp
+      zwi(:,:,:) = 0._wp
+      !
+      l_trd = .FALSE.            ! set local switches
+      l_hst = .FALSE.
+      l_ptr = .FALSE.
+      ll_zAimp = .FALSE.
+      IF( ( cdtype == 'TRA' .AND. l_trdtra  ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
+      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) )    l_ptr = .TRUE.
+      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR.  &
+         &                         iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
+      !
+      IF( l_trd .OR. l_hst )  THEN
+         ALLOCATE( ztrdx(jpi,jpj,jpk), ztrdy(jpi,jpj,jpk), ztrdz(jpi,jpj,jpk) )
+         ztrdx(:,:,:) = 0._wp   ;    ztrdy(:,:,:) = 0._wp   ;   ztrdz(:,:,:) = 0._wp
+      ENDIF
+      !
+      IF( l_ptr ) THEN
+         ALLOCATE( zptry(jpi,jpj,jpk) )
+         zptry(:,:,:) = 0._wp
+      ENDIF
+      !
+      ! If adaptive vertical advection, check if it is needed on this PE at this time
+      IF( ln_zad_Aimp ) THEN
+         IF( MAXVAL( ABS( wi(:,:,:) ) ) > 0._wp ) ll_zAimp = .TRUE.
+      END IF
+      ! If active adaptive vertical advection, build tridiagonal matrix
+      IF( ll_zAimp ) THEN
+         ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk))
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+            zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )   &
+            &                               / e3t(ji,jj,jk,Krhs)
+            zwinf(ji,jj,jk) =  p2dt * MIN( wi(ji,jj,jk  ) , 0._wp ) / e3t(ji,jj,jk,Krhs)
+            zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Krhs)
+         END_3D
+      END IF
+      !
+      DO jn = 1, kjpt            !==  loop over the tracers  ==!
+         !
+         !        !==  upstream advection with initial mass fluxes & intermediate update  ==!
+         !                               !* upstream tracer flux in the k direction *!
+         DO_3D( 1, 1, 1, 1, 2, jpkm1 )      ! Interior value ( multiplied by wmask)
+            zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) )
+            zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) )
+            zwz(ji,jj,jk) = 0.5 * ( zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb) ) * wmask(ji,jj,jk)
+         END_3D
+         IF( ln_linssh ) THEN               ! top ocean value (only in linear free surface as zwz has been w-masked)
+            IF( ln_isfcav ) THEN                        ! top of the ice-shelf cavities and at the ocean surface
+               DO_2D( 1, 1, 1, 1 )
+                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface
+               END_2D
+            ELSE                                        ! no cavities: only at the ocean surface
+               DO_2D( 1, 1, 1, 1 )
+                  zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb)
+               END_2D
+            ENDIF
+         ENDIF
+         !
+         !                    !* upstream tracer flux in the i and j direction
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj-1
+               tracer_flux_i(zwx_3d(1,jj,jk),zfp_ui,zfm_ui,1,jj,jk)
+               tracer_flux_j(zwy_3d(1,jj,jk),zfp_vj,zfm_vj,1,jj,jk)
+            END DO
+            DO ji = 1, jpi-1
+               tracer_flux_i(zwx_3d(ji,1,jk),zfp_ui,zfm_ui,ji,1,jk)
+               tracer_flux_j(zwy_3d(ji,1,jk),zfp_vj,zfm_vj,ji,1,jk)
+            END DO
+            DO_2D( 1, 1, 1, 1 )
+               tracer_flux_i(zwx_3d(ji,jj,jk),zfp_ui,zfm_ui,ji,jj,jk)
+               tracer_flux_i(zwx_im1,zfp_ui_m1,zfm_ui_m1,ji-1,jj,jk)
+               tracer_flux_j(zwy_3d(ji,jj,jk),zfp_vj,zfm_vj,ji,jj,jk)
+               tracer_flux_j(zwy_jm1,zfp_vj_m1,zfm_vj_m1,ji,jj-1,jk)
+               ztra = - ( zwx_3d(ji,jj,jk) - zwx_im1 + zwy_3d(ji,jj,jk) - zwy_jm1 + zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) * r1_e1e2t(ji,jj)
+               !                               ! update and guess with monotonic sheme
+               pt(ji,jj,jk,jn,Krhs) =                   pt(ji,jj,jk,jn,Krhs) +       ztra   &
+                  &                                  / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk)
+               zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) &
+                  &                                  / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
+            END_2D
+         END DO
+
+         IF ( ll_zAimp ) THEN
+            CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 )
+            !
+            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ;
+            DO_3D( 1, 1, 1, 1, 2, jpkm1 )       ! Interior value ( multiplied by wmask)
+               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
+               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
+               ztw(ji,jj,jk) =  0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk)
+               zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes
+            END_3D
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) &
+                  &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+            END_3D
+            !
+         END IF
+         !
+         IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes)
+            ztrdx(:,:,:) = zwx_3d(:,:,:)   ;   ztrdy(:,:,:) = zwy_3d(:,:,:)   ;   ztrdz(:,:,:) = zwz(:,:,:)
+         END IF
+         !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes)
+         IF( l_ptr )   zptry(:,:,:) = zwy_3d(:,:,:)
+         !
+         !        !==  anti-diffusive flux : high order minus low order  ==!
+         !
+         SELECT CASE( kn_fct_h )    !* horizontal anti-diffusive fluxes
+         !
+         CASE(  2  )                   !- 2nd order centered
+            DO_3D( 2, 1, 2, 1, 1, jpkm1 )
+               zwx_3d(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj,jk,jn,Kmm) ) - zwx_3d(ji,jj,jk)
+               zwy_3d(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj+1,jk,jn,Kmm) ) - zwy_3d(ji,jj,jk)
+            END_3D
+            !
+         CASE(  4  )                   !- 4th order centered
+            zltu_3d(:,:,jpk) = 0._wp            ! Bottom value : flux set to zero
+            zltv_3d(:,:,jpk) = 0._wp
+            !                          ! Laplacian
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )                 ! 2nd derivative * 1/ 6
+                  !             ! 1st derivative (gradient)
+                  ztu = ( pt(ji+1,jj,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
+                  ztu_im1 = ( pt(ji,jj,jk,jn,Kmm) - pt(ji-1,jj,jk,jn,Kmm) ) * umask(ji-1,jj,jk)
+                  ztv = ( pt(ji,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
+                  ztv_jm1 = ( pt(ji,jj,jk,jn,Kmm) - pt(ji,jj-1,jk,jn,Kmm) ) * vmask(ji,jj-1,jk)
+                  !             ! 2nd derivative * 1/ 6
+                  zltu_3d(ji,jj,jk) = (  ztu + ztu_im1  ) * r1_6
+                  zltv_3d(ji,jj,jk) = (  ztv + ztv_jm1  ) * r1_6
+               END_2D
+            END DO
+            ! NOTE [ comm_cleanup ] : need to change sign to ensure halo 1 - halo 2 compatibility
+            CALL lbc_lnk( 'traadv_fct', zltu_3d, 'T', -1.0_wp , zltv_3d, 'T', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
+            !
+            DO_3D( 2, 1, 2, 1, 1, jpkm1 )
+               zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points
+               zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
+               !                                                        ! C4 minus upstream advective fluxes
+               ! round brackets added to fix the order of floating point operations
+               ! needed to ensure halo 1 - halo 2 compatibility
+               zwx_3d(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * ( zC2t_u + ( zltu_3d(ji,jj,jk) - zltu_3d(ji+1,jj,jk)   &
+                             &                                        )                                           & ! bracket for halo 1 - halo 2 compatibility
+                             &                             ) - zwx_3d(ji,jj,jk)
+               zwy_3d(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * ( zC2t_v + ( zltv_3d(ji,jj,jk) - zltv_3d(ji,jj+1,jk)   &
+                             &                                        )                                           & ! bracket for halo 1 - halo 2 compatibility
+                             &                             ) - zwy_3d(ji,jj,jk)
+            END_3D
+            !
+         CASE(  41 )                   !- 4th order centered       ==>>   !!gm coding attempt   need to be tested
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )    ! Horizontal advective fluxes
+               ztu_im1 = ( pt(ji  ,jj  ,jk,jn,Kmm) - pt(ji-1,jj,jk,jn,Kmm) ) * umask(ji-1,jj,jk)
+               ztu_ip1 = ( pt(ji+2,jj  ,jk,jn,Kmm) - pt(ji+1,jj,jk,jn,Kmm) ) * umask(ji+1,jj,jk)
+
+               ztv_jm1 = ( pt(ji,jj  ,jk,jn,Kmm) - pt(ji,jj-1,jk,jn,Kmm) ) * vmask(ji,jj-1,jk)
+               ztv_jp1 = ( pt(ji,jj+2,jk,jn,Kmm) - pt(ji,jj+1,jk,jn,Kmm) ) * vmask(ji,jj+1,jk)
+               zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points (x2)
+               zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
+               !                                                  ! C4 interpolation of T at u- & v-points (x2)
+               zC4t_u =  zC2t_u + r1_6 * ( ztu_im1 - ztu_ip1 )
+               zC4t_v =  zC2t_v + r1_6 * ( ztv_jm1 - ztv_jp1 )
+               !                                                  ! C4 minus upstream advective fluxes
+               zwx_3d(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx_3d(ji,jj,jk)
+               zwy_3d(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy_3d(ji,jj,jk)
+            END_3D
+            CALL lbc_lnk( 'traadv_fct', zwx_3d, 'U', -1.0_wp , zwy_3d, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
+            !
+         END SELECT
+         !
+         SELECT CASE( kn_fct_v )    !* vertical anti-diffusive fluxes (w-masked interior values)
+         !
+         CASE(  2  )                   !- 2nd order centered
+            DO_3D( 1, 1, 1, 1, 2, jpkm1 )
+               zwz(ji,jj,jk) =  (  pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   &
+                  &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk)
+            END_3D
+            !
+         CASE(  4  )                   !- 4th order COMPACT
+            CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )   ! zwt = COMPACT interpolation of T at w-point
+            DO_3D( 1, 1, 1, 1, 2, jpkm1 )
+               zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
+            END_3D
+            !
+         END SELECT
+         IF( ln_linssh ) THEN    ! top ocean value: high order = upstream  ==>>  zwz=0
+            zwz(:,:,1) = 0._wp   ! only ocean surface as interior zwz values have been w-masked
+         ENDIF
+         !
+         CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp)
+         !
+         IF ( ll_zAimp ) THEN
+            DO_3D( 1, 1, 1, 1, 1, jpkm1 )    !* trend and after field with monotonic scheme
+               !                                                ! total intermediate advective trends
+               ztra = - (  zwx_3d(ji,jj,jk) - zwx_3d(ji-1,jj  ,jk  )   &
+                  &      + zwy_3d(ji,jj,jk) - zwy_3d(ji  ,jj-1,jk  )   &
+                  &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
+               ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
+            END_3D
+            !
+            CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 )
+            !
+            DO_3D( 1, 1, 1, 1, 2, jpkm1 )       ! Interior value ( multiplied by wmask)
+               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
+               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
+               zwz(ji,jj,jk) = zwz(ji,jj,jk) + 0.5 * e1e2t(ji,jj) * ( zfp_wk * ztw(ji,jj,jk) + zfm_wk * ztw(ji,jj,jk-1) ) * wmask(ji,jj,jk)
+            END_3D
+         END IF
+         !
+         !        !==  monotonicity algorithm  ==!
+         !
+         CALL nonosc( Kmm, pt(:,:,:,jn,Kbb), zwx_3d, zwy_3d, zwz, zwi, p2dt )
+         !
+         !        !==  final trend with corrected fluxes  ==!
+         !
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            ztra = - (  zwx_3d(ji,jj,jk) - zwx_3d(ji-1,jj  ,jk  )   &
+               &      + zwy_3d(ji,jj,jk) - zwy_3d(ji  ,jj-1,jk  )   &
+               &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
+            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm)
+            zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
+         END_3D
+         !
+         IF ( ll_zAimp ) THEN
+            !
+            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )      ! Interior value ( multiplied by wmask)
+               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
+               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
+               ztw(ji,jj,jk) = - 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk)
+               zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic
+            END_3D
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) &
+                  &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+            END_3D
+         END IF
+         ! NOT TESTED - NEED l_trd OR l_hst TRUE
+         IF( l_trd .OR. l_hst ) THEN   ! trend diagnostics // heat/salt transport
+            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx_3d(:,:,:)  ! <<< add anti-diffusive fluxes
+            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy_3d(:,:,:)  !     to upstream fluxes
+            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  !
+            !
+            IF( l_trd ) THEN              ! trend diagnostics
+               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztrdx, pU, pt(:,:,:,jn,Kmm) )
+               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztrdy, pV, pt(:,:,:,jn,Kmm) )
+               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, ztrdz, pW, pt(:,:,:,jn,Kmm) )
+            ENDIF
+            !                             ! heat/salt transport
+            IF( l_hst )   CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) )
+            !
+         ENDIF
+         ! NOT TESTED - NEED l_ptr TRUE
+         IF( l_ptr ) THEN              ! "Poleward" transports
+            zptry(:,:,:) = zptry(:,:,:) + zwy_3d(:,:,:)  ! <<< add anti-diffusive fluxes
+            CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
+         ENDIF
+         !
+      END DO                     ! end of tracer loop
+      !
+      IF ( ll_zAimp ) THEN
+         DEALLOCATE( zwdia, zwinf, zwsup )
+      ENDIF
+      IF( l_trd .OR. l_hst ) THEN
+         DEALLOCATE( ztrdx, ztrdy, ztrdz )
+      ENDIF
+      IF( l_ptr ) THEN
+         DEALLOCATE( zptry )
+      ENDIF
+      !
+   END SUBROUTINE tra_adv_fct_lf
+#endif
    !!======================================================================
 END MODULE traadv_fct

NEMO/trunk/src/OCE/TRA/traadv_mus.F90

@@ -81 +81 @@
       LOGICAL                                  , INTENT(in   ) ::   ld_msc_ups      ! use upstream scheme within muscl
       REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
-      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
+      ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
       REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
@@ -93 +93 @@
       !!----------------------------------------------------------------------
       !
-      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
          IF( kt == kit000 )  THEN
             IF(lwp) WRITE(numout,*)
@@ -139 +139 @@
             zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( pt(ji,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
          END_3D
-         ! lateral boundary conditions   (changed sign)
-         IF ( nn_hls.EQ.1 ) CALL lbc_lnk( 'traadv_mus', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp )
          !                                !-- Slopes of tracer
          zslpx(:,:,jpk) = 0._wp                 ! bottom values
          zslpy(:,:,jpk) = 0._wp
-         DO_3D( nn_hls-1, 1, nn_hls-1, 1, 1, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
             zslpx(ji,jj,jk) =                       ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   &
                &            * ( 0.25 + SIGN( 0.25_wp, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) )
@@ -151 +149 @@
          END_3D
          !
-         DO_3D( nn_hls-1, 1, nn_hls-1, 1, 1, jpkm1 )    !-- Slopes limitation
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )    !-- Slopes limitation
             zslpx(ji,jj,jk) = SIGN( 1.0_wp, zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   &
                &                                                     2.*ABS( zwx  (ji-1,jj,jk) ),   &
@@ -159 +157 @@
                &                                                     2.*ABS( zwy  (ji,jj  ,jk) ) )
          END_3D
-         !
-         DO_3D( nn_hls-1, 0, nn_hls-1, 0, 1, jpkm1 )    !-- MUSCL horizontal advective fluxes
+         ! NOTE [ comm_cleanup ] : need to change sign to ensure halo 1 - halo 2 compatibility
+         IF ( nn_hls==1 ) CALL lbc_lnk( 'traadv_mus', zslpx, 'T', -1.0_wp , zslpy, 'T', -1.0_wp )   ! lateral boundary conditions   (changed sign)
+         !
+         DO_3D( 1, 0, 1, 0, 1, jpkm1 )    !-- MUSCL horizontal advective fluxes
             ! MUSCL fluxes
             z0u = SIGN( 0.5_wp, pU(ji,jj,jk) )
@@ -176 +176 @@
             zwy(ji,jj,jk) = pV(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
          END_3D
-         IF ( nn_hls.EQ.1 ) CALL lbc_lnk( 'traadv_mus', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp )   ! lateral boundary conditions   (changed sign)
          !
          DO_3D( 0, 0, 0, 0, 1, jpkm1 )    !-- Tracer advective trend

NEMO/trunk/src/OCE/TRA/traadv_qck.F90

@@ -27 +27 @@
    USE lbclnk          ! ocean lateral boundary condition (or mpp link)
    USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
+#if defined key_loop_fusion
+   USE traadv_qck_lf   ! QCK    scheme            (tra_adv_qck  routine - loop fusion version)
+#endif
 
    IMPLICIT NONE
@@ -91 +94 @@
       INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers
       REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
-      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
+      ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
       REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume transport components
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
       !!----------------------------------------------------------------------
       !
-      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+#if defined key_loop_fusion
+      CALL tra_adv_qck_lf ( kt, kit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs )
+#else
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
          IF( kt == kit000 )  THEN
             IF(lwp) WRITE(numout,*)
@@ -117 +123 @@
       CALL tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs )
       !
+#endif
    END SUBROUTINE tra_adv_qck
 
@@ -129 +136 @@
       INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers
       REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step
-      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
+      ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
       REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU        ! i-velocity components
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
@@ -149 +156 @@
             zfd(ji,jj,jk) = pt(ji+1,jj,jk,jn,Kbb)        ! Downstream in the x-direction for the tracer
          END_3D
-         IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp )   ! Lateral boundary conditions
+         IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, ld4only= .TRUE. )   ! Lateral boundary conditions
 
          !
@@ -167 +174 @@
          END_3D
          !--- Lateral boundary conditions
-         IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp,  zwx(:,:,:), 'T', 1.0_wp )
+         IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp,  zwx(:,:,:), 'T', 1.0_wp )
 
          !--- QUICKEST scheme
@@ -176 +183 @@
             zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
          END_3D
-         IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp )      ! Lateral boundary conditions
+         IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp, ld4only= .TRUE. )      ! Lateral boundary conditions
 
          !
@@ -214 +221 @@
       INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers
       REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step
-      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
+      ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
       REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pV        ! j-velocity components
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
@@ -229 +236 @@
          zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0
          !
-         DO jk = 1, jpkm1
-            !
-            !--- Computation of the ustream and downstream value of the tracer and the mask
-            DO_2D( 0, 0, nn_hls-1, nn_hls-1 )
-               ! Upstream in the x-direction for the tracer
-               zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb)
-               ! Downstream in the x-direction for the tracer
-               zfd(ji,jj,jk) = pt(ji,jj+1,jk,jn,Kbb)
-            END_2D
-         END DO
-         IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp )   ! Lateral boundary conditions
+         !--- Computation of the ustream and downstream value of the tracer and the mask
+         DO_3D( 0, 0, nn_hls-1, nn_hls-1, 1, jpkm1 )
+            ! Upstream in the x-direction for the tracer
+            zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb)
+            ! Downstream in the x-direction for the tracer
+            zfd(ji,jj,jk) = pt(ji,jj+1,jk,jn,Kbb)
+         END_3D
+
+         IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, ld4only= .TRUE. )   ! Lateral boundary conditions
 
          !
@@ -259 +264 @@
 
          !--- Lateral boundary conditions
-         IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwy(:,:,:), 'T', 1.0_wp )
+         IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwy(:,:,:), 'T', 1.0_wp )
 
          !--- QUICKEST scheme
@@ -268 +273 @@
             zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
          END_3D
-         IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp )    !--- Lateral boundary conditions
+         IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp, ld4only= .TRUE. )    !--- Lateral boundary conditions
          !
          ! Tracer flux on the x-direction
@@ -306 +311 @@
       CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
       INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers
-      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
+      ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
       REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pW      ! vertical velocity
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
@@ -365 +370 @@
       !----------------------------------------------------------------------
       !
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+      DO_3D( 1, 0, 1, 0, 1, jpkm1 )
          zc     = puc(ji,jj,jk)                         ! Courant number
          zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)

NEMO/trunk/src/OCE/TRA/traadv_ubs.F90

@@ -26 +26 @@
    USE lbclnk         ! ocean lateral boundary condition (or mpp link)
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
+#if defined key_loop_fusion
+   USE traadv_ubs_lf  ! UBS      scheme            (tra_adv_ubs  routine - loop fusion version)
+#endif
 
    IMPLICIT NONE
@@ -92 +95 @@
       INTEGER                                  , INTENT(in   ) ::   kn_ubs_v        ! number of tracers
       REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
-      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
+      ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
       REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume transport components
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
@@ -103 +106 @@
       !!----------------------------------------------------------------------
       !
-      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+#if defined key_loop_fusion
+      CALL tra_adv_ubs_lf    ( kt, kit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs, kn_ubs_v )
+#else
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
          IF( kt == kit000 )  THEN
             IF(lwp) WRITE(numout,*)
@@ -140 +146 @@
             !
          END DO
-         IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1.0_wp, zltv, 'T', 1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
+         IF (nn_hls==1) CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1.0_wp, zltv, 'T', 1.0_wp, ld4only= .TRUE. )   ! Lateral boundary cond. (unchanged sgn)
          !
          DO_3D( 1, 0, 1, 0, 1, jpkm1 )   !==  Horizontal advective fluxes  ==!     (UBS)
@@ -260 +266 @@
       END DO
       !
+#endif
    END SUBROUTINE tra_adv_ubs
 

NEMO/trunk/src/OCE/TRA/trabbc.F90

@@ -102 +102 @@
       ENDIF
       !
-      IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
-         CALL iom_put ( "hfgeou" , rho0_rcp * qgh_trd0(:,:) )
-      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' )
       !

NEMO/trunk/src/OCE/TRA/trabbl.F90

@@ -126 +126 @@
          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' )
-         IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
-            CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef
-            CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
-         ENDIF
+         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
@@ -139 +137 @@
          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' )
-         IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
-            ! lateral boundary conditions ; just need for outputs
-            CALL lbc_lnk( 'trabbl', utr_bbl, 'U', -1.0_wp , vtr_bbl, 'V', -1.0_wp )
-            CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport
-            CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
-         ENDIF
+         CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport
+         CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
          !
       ENDIF
@@ -215 +209 @@
 
 
+   ! NOTE: [tiling] tiling changes the results, but only the order of floating point operations is different
    SUBROUTINE tra_bbl_adv( pt, pt_rhs, kjpt, Kmm )
       !!----------------------------------------------------------------------
@@ -238 +233 @@
       INTEGER  ::   iis , iid , ijs , ijd    ! local integers
       INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
-      INTEGER  ::   isi, isj                 !   -       -
       REAL(wp) ::   zbtr, ztra               ! local scalars
       REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
       !!----------------------------------------------------------------------
-      !
-      IF( ntsi == Nis0 ) THEN ; isi = 1 ; ELSE ; isi = 0 ; ENDIF    ! Avoid double-counting when using tiling
-      IF( ntsj == Njs0 ) THEN ; isj = 1 ; ELSE ; isj = 0 ; ENDIF
       !                                                          ! ===========
       DO jn = 1, kjpt                                            ! tracer loop
          !                                                       ! ===========
-         DO_2D( isi, 0, isj, 0 )            ! CAUTION start from i=1 to update i=2 when cyclic east-west
+         DO_2D_OVR( 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)
@@ -340 +331 @@
       !!----------------------------------------------------------------------
       !
-      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
          IF( kt == kit000 )  THEN
             IF(lwp)  WRITE(numout,*)
@@ -363 +354 @@
       IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
          !                                !-------------------!
-         DO_2D( 1, 0, 1, 0 )                   ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
+         DO_2D_OVR( 1, 0, 1, 0 )                   ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
             !                                                   ! i-direction
             za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
@@ -393 +384 @@
          !
          CASE( 1 )                                   != use of upper velocity
-            DO_2D( 1, 0, 1, 0 )                              ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
+            DO_2D_OVR( 1, 0, 1, 0 )                              ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
                !                                                  ! i-direction
                za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
@@ -422 +413 @@
          CASE( 2 )                                 != bbl velocity = F( delta rho )
             zgbbl = grav * rn_gambbl
-            DO_2D( 1, 0, 1, 0 )                         ! criteria: rho_up > rho_down
+            DO_2D_OVR( 1, 0, 1, 0 )                         ! criteria: rho_up > rho_down
                !                                                  ! i-direction
                ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)

NEMO/trunk/src/OCE/TRA/traisf.F90

@@ -47 +47 @@
       IF( ln_timing )   CALL timing_start('tra_isf')
       !
-      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
          IF( kt == nit000 ) THEN
             IF(lwp) WRITE(numout,*)
@@ -79 +79 @@
       !
       IF ( ln_isfdebug ) THEN
-         IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only for the full domain
+         IF( .NOT. l_istiled .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))

NEMO/trunk/src/OCE/TRA/traldf.F90

@@ -17 +17 @@
    USE oce            ! ocean dynamics and tracers
    USE dom_oce        ! ocean space and time domain
-   ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed from tra_ldf_blp, zps_hde*)
-   USE domtile
    USE phycst         ! physical constants
    USE ldftra         ! lateral diffusion: eddy diffusivity & EIV coeff.
@@ -58 +56 @@
       !!
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
-      ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed from tra_ldf_blp, zps_hde*)
-      LOGICAL :: lskip
       !!----------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('tra_ldf')
       !
-      lskip = .FALSE.
-
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
          ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
@@ -71 +65 @@
          ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
       ENDIF
-
-      ! TEMP: [tiling] These changes not necessary after extra haloes development (lbc_lnk removed from tra_ldf_blp, zps_hde*)
-      IF( nldf_tra == np_blp .OR. nldf_tra == np_blp_i .OR. nldf_tra == np_blp_it )  THEN
-         IF( ln_tile ) THEN
-            IF( ntile == 1 ) THEN
-               CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )
-            ELSE
-               lskip = .TRUE.
-            ENDIF
-         ENDIF
-      ENDIF
-      IF( .NOT. lskip ) THEN
-         !
-         SELECT CASE ( nldf_tra )                 !* compute lateral mixing trend and add it to the general trend
-         CASE ( np_lap   )                                  ! laplacian: iso-level operator
-            CALL tra_ldf_lap  ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs),                   jpts,  1 )
-         CASE ( np_lap_i )                                  ! laplacian: standard iso-neutral operator (Madec)
-            CALL tra_ldf_iso  ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts,  1 )
-         CASE ( np_lap_it )                                 ! laplacian: triad iso-neutral operator (griffies)
-            CALL tra_ldf_triad( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts,  1 )
-         CASE ( np_blp , np_blp_i , np_blp_it )             ! bilaplacian: iso-level & iso-neutral operators
-            IF(nn_hls.EQ.2) CALL lbc_lnk( 'tra_ldf', pts(:,:,:,:,Kbb), 'T',1.)
-            CALL tra_ldf_blp  ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs),             jpts, nldf_tra )
-         END SELECT
-         !
-         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_ldf, ztrdt )
-            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_ldf, ztrds )
-            DEALLOCATE( ztrdt, ztrds )
-         ENDIF
-
-         ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed from tra_ldf_blp, zps_hde*)
-         IF( ln_tile .AND. ntile == 0 ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 1 )
+      !
+      SELECT CASE ( nldf_tra )                 !* compute lateral mixing trend and add it to the general trend
+      CASE ( np_lap   )                                  ! laplacian: iso-level operator
+         CALL tra_ldf_lap  ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs),                   jpts,  1 )
+      CASE ( np_lap_i )                                  ! laplacian: standard iso-neutral operator (Madec)
+         CALL tra_ldf_iso  ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts,  1 )
+      CASE ( np_lap_it )                                 ! laplacian: triad iso-neutral operator (griffies)
+         CALL tra_ldf_triad( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts,  1 )
+      CASE ( np_blp , np_blp_i , np_blp_it )             ! bilaplacian: iso-level & iso-neutral operators
+         CALL tra_ldf_blp  ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs),             jpts, nldf_tra )
+      END SELECT
+      !
+      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_ldf, ztrdt )
+         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_ldf, ztrds )
+         DEALLOCATE( ztrdt, ztrds )
       ENDIF
       !                                        !* print mean trends (used for debugging)

NEMO/trunk/src/OCE/TRA/traldf_iso.F90

@@ -132 +132 @@
       INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
       INTEGER  ::  ikt
-      INTEGER  ::  ierr             ! local integer
+      INTEGER  ::  ierr, iij        ! local integer
       REAL(wp) ::  zmsku, zahu_w, zabe1, zcof1, zcoef3   ! local scalars
       REAL(wp) ::  zmskv, zahv_w, zabe2, zcof2, zcoef4   !   -      -
@@ -141 +141 @@
       !
       IF( kpass == 1 .AND. kt == kit000 )  THEN
-         IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+         IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype
@@ -147 +147 @@
          ENDIF
          !
-         DO_3D( 0, 0, 0, 0, 1, jpk )
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
             akz     (ji,jj,jk) = 0._wp
             ah_wslp2(ji,jj,jk) = 0._wp
@@ -153 +153 @@
       ENDIF
       !
-      IF( ntile == 0 .OR. ntile == 1 )  THEN                           ! Do only on the first tile
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                           ! Do only on the first tile
          l_hst = .FALSE.
          l_ptr = .FALSE.
@@ -161 +161 @@
       ENDIF
       !
+      ! Define pt_rhs halo points for multi-point haloes in bilaplacian case
+      IF( nldf_tra == np_blp_i .AND. kpass == 1 ) THEN ; iij = nn_hls
+      ELSE                                             ; iij = 1
+      ENDIF
+
       !
       IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
@@ -172 +177 @@
       IF( kpass == 1 ) THEN                  !==  first pass only  ==!
          !
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             !
             zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
@@ -185 +190 @@
                &       + ( pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)                   &
                &         )                                                         & ! bracket for halo 1 - halo 2 compatibility
-               &     ) * zmsku                                       
+               &     ) * zmsku
             zahv_w = ( (  pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)                    &
                &       )                                                           & ! bracket for halo 1 - halo 2 compatibility
-               &       + ( pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)                   & 
+               &       + ( pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)                   &
                &         )                                                         & ! bracket for halo 1 - halo 2 compatibility
-               &     ) * zmskv                                               
+               &     ) * zmskv
                !
             ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
@@ -197 +202 @@
          !
          IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient
-            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+            DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
                ! round brackets added to fix the order of floating point operations
                ! needed to ensure halo 1 - halo 2 compatibility
@@ -205 +210 @@
                   &            )                                                                               & ! bracket for halo 1 - halo 2 compatibility
                   &            + ( ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) ) &
-                  &              + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) &  
+                  &              + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) &
                   &              )                                                                             & ! bracket for halo 1 - halo 2 compatibility
-                  &                      )                                                                         
+                  &                      )
             END_3D
             !
             IF( ln_traldf_blp ) THEN                ! bilaplacian operator
-               DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+               DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
                   akz(ji,jj,jk) = 16._wp   &
                      &   * ah_wslp2   (ji,jj,jk)   &
@@ -219 +224 @@
                END_3D
             ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
-               DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+               DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 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  )
@@ -227 +232 @@
            !
          ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
-            DO_3D( 0, 0, 0, 0, 1, jpk )
+            DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
                akz(ji,jj,jk) = ah_wslp2(ji,jj,jk)
             END_3D
@@ -240 +245 @@
          !!   I - masked horizontal derivative
          !!----------------------------------------------------------------------
-!!gm : bug.... why (x,:,:)?   (1,jpj,:) and (jpi,1,:) should be sufficient....
-         zdit (ntsi-nn_hls,:,:) = 0._wp     ;     zdit (ntei+nn_hls,:,:) = 0._wp
-         zdjt (ntsi-nn_hls,:,:) = 0._wp     ;     zdjt (ntei+nn_hls,:,:) = 0._wp
-         !!end
+         zdit(:,:,:) = 0._wp
+         zdjt(:,:,:) = 0._wp
 
          ! Horizontal tracer gradient
-         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+         DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 )
             zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
             zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
          END_3D
          IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient
-            DO_2D( 1, 0, 1, 0 )           ! bottom correction (partial bottom cell)
+            DO_2D( iij, iij-1, iij, iij-1 )            ! bottom correction (partial bottom cell)
                zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)
                zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
             END_2D
             IF( ln_isfcav ) THEN      ! first wet level beneath a cavity
-               DO_2D( 1, 0, 1, 0 )
+               DO_2D( iij, iij-1, iij, iij-1 )
                   IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)
                   IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)
@@ -269 +272 @@
          DO jk = 1, jpkm1                                 ! Horizontal slab
             !
-            DO_2D( 1, 1, 1, 1 )
+            DO_2D( iij, iij, iij, iij )
                !                             !== Vertical tracer gradient
                zdk1t(ji,jj) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1)     ! level jk+1
@@ -278 +281 @@
             END_2D
             !
-            DO_2D( 1, 0, 1, 0 )           !==  Horizontal fluxes
+            DO_2D( iij, iij-1, iij, iij-1 )           !==  Horizontal fluxes
                zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm)
                zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
@@ -296 +299 @@
                   &               + zcof1 * ( ( zdkt (ji+1,jj) + zdk1t(ji,jj)    &
                   &                           )                                  & ! bracket for halo 1 - halo 2 compatibility
-                  &                         + ( zdk1t(ji+1,jj) + zdkt (ji,jj)    & 
+                  &                         + ( zdk1t(ji+1,jj) + zdkt (ji,jj)    &
                   &                           )                                  & ! bracket for halo 1 - halo 2 compatibility
                   &                         ) ) * umask(ji,jj,jk)
@@ -304 +307 @@
                   &                         + ( zdk1t(ji,jj+1) + zdkt (ji,jj)    &
                   &                           )                                  & ! bracket for halo 1 - halo 2 compatibility
-                  &                         ) ) * vmask(ji,jj,jk)      
+                  &                         ) ) * vmask(ji,jj,jk)
             END_2D
             !
-            DO_2D( 0, 0, 0, 0 )           !== horizontal divergence and add to pta
+            DO_2D( iij-1, iij-1, iij-1, iij-1 )           !== horizontal divergence and add to pta
                ! round brackets added to fix the order of floating point operations
                ! needed to ensure halo 1 - halo 2 compatibility
@@ -328 +331 @@
          ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp
 
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )    ! interior (2=<jk=<jpk-1)
+         DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 )    ! interior (2=<jk=<jpk-1)
             !
             zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
@@ -348 +351 @@
                   &                     )                                              & ! bracket for halo 1 - halo 2 compatibility
                   &                   + ( zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)    &
-                  &                     )                                              & ! bracket for halo 1 - halo 2 compatibility 
+                  &                     )                                              & ! bracket for halo 1 - halo 2 compatibility
                   &                   )                                                &
                   &        + zcoef4 * ( ( zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)    &
@@ -358 +361 @@
          !                                !==  add the vertical 33 flux  ==!
          IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
-            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+            DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 )
                ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)   &
                   &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )               &
@@ -367 +370 @@
             SELECT CASE( kpass )
             CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
-               DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+               DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 )
                   ztfw(ji,jj,jk) =   &
                      &  ztfw(ji,jj,jk) + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   &
@@ -381 +384 @@
          ENDIF
          !
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )    !==  Divergence of vertical fluxes added to pta  ==!
+         DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 )    !==  Divergence of vertical fluxes added to pta  ==!
             pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * r1_e1e2t(ji,jj)   &
                &                                             / e3t(ji,jj,jk,Kmm)

NEMO/trunk/src/OCE/TRA/traldf_lap_blp.F90

@@ -103 +103 @@
       !
       INTEGER  ::   ji, jj, jk, jn      ! dummy loop indices
-      INTEGER  ::   isi, iei, isj, iej  ! local integers
+      INTEGER  ::   iij
       REAL(wp) ::   zsign               ! local scalars
       REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   ztu, ztv, zaheeu, zaheev
       !!----------------------------------------------------------------------
       !
-      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
          IF( kt == nit000 .AND. lwp )  THEN
             WRITE(numout,*)
@@ -122 +122 @@
       ENDIF
       !
+      ! Define pt_rhs halo points for multi-point haloes in bilaplacian case
+      IF( nldf_tra == np_blp .AND. kpass == 1 ) THEN ; iij = nn_hls
+      ELSE                                           ; iij = 1
+      ENDIF
+
       !                                !==  Initialization of metric arrays used for all tracers  ==!
       IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
@@ -127 +132 @@
       ENDIF
 
-      IF( ntsi == Nis0 ) THEN ; isi = nn_hls - 1 ; ELSE ; isi = 0 ; ENDIF    ! Avoid double-counting when using tiling
-      IF( ntsj == Njs0 ) THEN ; isj = nn_hls - 1 ; ELSE ; isj = 0 ; ENDIF
-      IF( ntei == Nie0 ) THEN ; iei = nn_hls - 1 ; ELSE ; iei = 0 ; ENDIF
-      IF( ntej == Nje0 ) THEN ; iej = nn_hls - 1 ; ELSE ; iej = 0 ; ENDIF
-
-      DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )            !== First derivative (gradient)  ==!
+      DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 )            !== First derivative (gradient)  ==!
          zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm)   !!gm   * umask(ji,jj,jk) pah masked!
          zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm)   !!gm   * vmask(ji,jj,jk)
@@ -141 +141 @@
          !                          ! =========== !
          !
-         DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )            !== First derivative (gradient)  ==!
+         DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 )            !== First derivative (gradient)  ==!
             ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) )
             ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) )
          END_3D
          IF( ln_zps ) THEN                             ! set gradient at bottom/top ocean level
-            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )                              ! bottom
+            DO_2D( iij, iij-1, iij, iij-1 )                              ! bottom
                ztu(ji,jj,mbku(ji,jj)) = zaheeu(ji,jj,mbku(ji,jj)) * pgu(ji,jj,jn)
                ztv(ji,jj,mbkv(ji,jj)) = zaheev(ji,jj,mbkv(ji,jj)) * pgv(ji,jj,jn)
             END_2D
             IF( ln_isfcav ) THEN                             ! top in ocean cavities only
-               DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
+               DO_2D( iij, iij-1, iij, iij-1 )
                   IF( miku(ji,jj) > 1 )   ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn)
                   IF( mikv(ji,jj) > 1 )   ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn)
@@ -158 +158 @@
          ENDIF
          !
-         DO_3D( isi, iei, isj, iej, 1, jpkm1 )            !== Second derivative (divergence) added to the general tracer trends  ==!
+         DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 )            !== Second derivative (divergence) added to the general tracer trends  ==!
             ! round brackets added to fix the order of floating point operations
             ! needed to ensure halo 1 - halo 2 compatibility
@@ -215 +215 @@
       !!---------------------------------------------------------------------
       !
-      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
          IF( kt == kit000 .AND. lwp )  THEN
             WRITE(numout,*)
@@ -239 +239 @@
       END SELECT
       !
-      CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1.0_wp )     ! Lateral boundary conditions (unchanged sign)
+      IF (nn_hls==1) CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1.0_wp )     ! Lateral boundary conditions (unchanged sign)
       !                                               ! Partial top/bottom cell: GRADh( zlap )
       IF( ln_isfcav .AND. ln_zps ) THEN   ;   CALL zps_hde_isf( kt, Kmm, kjpt, zlap, zglu, zglv, zgui, zgvi )  ! both top & bottom

NEMO/trunk/src/OCE/TRA/traldf_triad.F90

@@ -13 +13 @@
    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 domtile
    USE domutl, ONLY : is_tile
    USE phycst         ! physical constants
@@ -109 +107 @@
       REAL(wp), DIMENSION(A2D_T(ktt_rhs),JPK,KJPT), INTENT(inout) ::   pt_rhs     ! tracer trend
       !
-      INTEGER  ::  ji, jj, jk, jn, kp   ! dummy loop indices
+      INTEGER  ::  ji, jj, jk, jn, kp, iij   ! dummy loop indices
       REAL(wp) ::  zcoef0, ze3w_2, zsign          !   -      -
       !
@@ -115 +113 @@
       REAL(wp) ::   ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt, zdyt_jp1, ze3wr_jp1, zdzt_jp1, zah_slp1, zah_slp_jp1, zaei_slp_jp1
       REAL(wp) ::   zah_slp, zaei_slp, zdxt_ip1, ze3wr_ip1, zdzt_ip1, zah_slp_ip1, zaei_slp_ip1, zaei_slp1
-      REAL(wp), DIMENSION(A2D(nn_hls),0:1)     ::   zdkt3d                         ! vertical tracer gradient at 2 levels
-      REAL(wp), DIMENSION(A2D(nn_hls)        ) ::   z2d                            ! 2D workspace
-      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
+      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, zpsi_uw, zpsi_vw   ! 3D     -
       !!----------------------------------------------------------------------
       !
-      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
          IF( kpass == 1 .AND. kt == kit000 )  THEN
             IF(lwp) WRITE(numout,*)
@@ -138 +134 @@
       ENDIF
       !
+      ! Define pt_rhs halo points for multi-point haloes in bilaplacian case
+      IF( nldf_tra == np_blp_it .AND. kpass == 1 ) THEN ; iij = nn_hls
+      ELSE                                              ; iij = 1
+      ENDIF
+
+      !
       IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
       ELSE                    ;   zsign = -1._wp
@@ -148 +150 @@
       IF( kpass == 1 ) THEN         !==  first pass only  and whatever the tracer is  ==!
          !
-         DO_3D( 0, 0, 0, 0, 1, jpk )
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
             akz     (ji,jj,jk) = 0._wp
             ah_wslp2(ji,jj,jk) = 0._wp
@@ -154 +156 @@
          !
          DO kp = 0, 1                            ! i-k triads
-            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
                ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm)
                zbu   = e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
@@ -177 +179 @@
          !
          DO kp = 0, 1                            ! j-k triads
-            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
                ze3wr = 1.0_wp / e3w(ji,jj,jk+kp,Kmm)
                zbv   = e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
@@ -203 +205 @@
             !
             IF( ln_traldf_blp ) THEN                ! bilaplacian operator
-               DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+               DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
                   akz(ji,jj,jk) = 16._wp           &
                      &   * ah_wslp2   (ji,jj,jk)   &
@@ -211 +213 @@
                END_3D
             ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
-               DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+               DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 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  )
@@ -219 +221 @@
            !
          ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
-            DO_3D( 0, 0, 0, 0, 1, jpk )
+            DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 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 kp = 0, 1
-                     DO_3D( 1, 0, 1, 0, 1, jpkm1 )
-                        ! round brackets added to fix the order of floating point operations
-                        ! needed to ensure halo 1 - halo 2 compatibility [ NOT TESTED ]
-                        zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp)                                     &
-                           & + ( 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji,jj,jk,1,kp)        & 
-                           &   + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+1,jj,jk,0,kp)      &
-                           &   )                                                                        ! bracket for halo 1 - halo 2 compatibility
-                        zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp)                                     &
-                           & + ( 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj,jk,1,kp)        &
-                           &   + 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj+1,jk,0,kp)      &
-                           &   )                                                                        ! bracket for halo 1 - halo 2 compatibility
-                     END_3D
-                  END DO
-               CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm )
-
-               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = nijtile )
-            ENDIF
+         IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) THEN
+            zpsi_uw(:,:,:) = 0._wp
+            zpsi_vw(:,:,:) = 0._wp
+
+            DO kp = 0, 1
+               DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+                  ! round brackets added to fix the order of floating point operations
+                  ! needed to ensure halo 1 - halo 2 compatibility
+                  zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp)                                     &
+                     & + ( 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji,jj,jk,1,kp)        &
+                     &   + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+1,jj,jk,0,kp)      &
+                     &   )                                                                        ! bracket for halo 1 - halo 2 compatibility
+                  zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp)                                     &
+                     & + ( 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj,jk,1,kp)        &
+                     &   + 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj+1,jk,0,kp)      &
+                     &   )                                                                        ! bracket for halo 1 - halo 2 compatibility
+               END_3D
+            END DO
+            CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm )
          ENDIF
          !
@@ -262 +257 @@
          zftu(:,:,:) = 0._wp
          zftv(:,:,:) = 0._wp
-         !
-         DO_3D( 1, 0, 1, 0, 1, jpkm1 )    !==  before lateral T & S gradients at T-level jk  ==!
+         zdit(:,:,:) = 0._wp
+         zdjt(:,:,:) = 0._wp
+         !
+         DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 )    !==  before lateral T & S gradients at T-level jk  ==!
             zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
             zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
          END_3D
          IF( ln_zps .AND. l_grad_zps ) THEN    ! partial steps: correction at top/bottom ocean level
-            DO_2D( 1, 0, 1, 0 )                    ! bottom level
+            DO_2D( iij, iij-1, iij, iij-1 )                    ! bottom level
                zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)
                zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
             END_2D
             IF( ln_isfcav ) THEN                   ! top level (ocean cavities only)
-               DO_2D( 1, 0, 1, 0 )
+               DO_2D( iij, iij-1, iij, iij-1 )
                   IF( miku(ji,jj)  > 1 )   zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn)
                   IF( mikv(ji,jj)  > 1 )   zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn)
@@ -286 +283 @@
          DO jk = 1, jpkm1
             !                    !==  Vertical tracer gradient at level jk and jk+1
-            DO_2D( 1, 1, 1, 1 )
+            DO_2D( iij, iij, iij, iij )
                zdkt3d(ji,jj,1) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1)
             END_2D
@@ -293 +290 @@
             IF( jk == 1 ) THEN   ;   zdkt3d(:,:,0) = zdkt3d(:,:,1)
             ELSE
-               DO_2D( 1, 1, 1, 1 )
+               DO_2D( iij, iij, iij, iij )
                   zdkt3d(ji,jj,0) = ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
                END_2D
@@ -305 +302 @@
             IF( ln_botmix_triad ) THEN
                DO kp = 0, 1              !==  Horizontal & vertical fluxes
-                  DO_2D( 1, 0, 1, 0 )
+                  DO_2D( iij, iij-1, iij, iij-1 )
                      ze1ur = r1_e1u(ji,jj)
                      zdxt  = zdit(ji,jj,jk) * ze1ur
@@ -317 +314 @@
                      zbu_ip1 = 0.25_wp * e1e2u(ji+1,jj) * e3u(ji+1,jj,jk,Kmm)
                      ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????   ahu is masked....
-                     zah = pahu(ji,jj,jk) 
-                     zah_ip1 = pahu(ji+1,jj,jk) 
+                     zah = pahu(ji,jj,jk)
+                     zah_ip1 = pahu(ji+1,jj,jk)
                      zah_slp  = zah * triadi(ji,jj,jk,1,kp)
                      zah_slp_ip1  = zah_ip1 * triadi(ji+1,jj,jk,1,kp)
@@ -333 +330 @@
                                          &      + ( zah * zdxt + zah_slp1 * zdzt_ip1 - zaei_slp1 * zdzt_ip1 ) * zbu * ze1ur  &
                                          &      )                                                                            ! bracket for halo 1 - halo 2 compatibility
-                     ztfw(ji+1,jj,jk+kp) =  ztfw(ji+1,jj,jk+kp)                                                              & 
+                     ztfw(ji+1,jj,jk+kp) =  ztfw(ji+1,jj,jk+kp)                                                              &
                                          &    - ( (zah_slp_ip1 + zaei_slp_ip1) * zdxt_ip1 * zbu_ip1 * ze3wr_ip1              &
-                                         &      + ( zah_slp1 + zaei_slp1) * zdxt * zbu * ze3wr_ip1                           & 
-                                         &      )                                                                            ! bracket for halo 1 - halo 2 compatibility
-                     END_2D
-                  END DO
+                                         &      + ( zah_slp1 + zaei_slp1) * zdxt * zbu * ze3wr_ip1                           &
+                                         &      )                                                                            ! bracket for halo 1 - halo 2 compatibility
+                  END_2D
+               END DO
                !
                DO kp = 0, 1
-                  DO_2D( 1, 0, 1, 0 )
+                  DO_2D( iij, iij-1, iij, iij-1 )
                      ze2vr = r1_e2v(ji,jj)
                      zdyt  = zdjt(ji,jj,jk) * ze2vr
@@ -353 +350 @@
                      ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????   ahu is masked....
                      zah = pahv(ji,jj,jk)          ! pahv(ji,jj+jp,jk)  ????
-                     zah_jp1 = pahv(ji,jj+1,jk) 
+                     zah_jp1 = pahv(ji,jj+1,jk)
                      zah_slp = zah * triadj(ji,jj,jk,1,kp)
                      zah_slp1 = zah * triadj(ji,jj+1,jk,0,kp)
                      zah_slp_jp1 = zah_jp1 * triadj(ji,jj+1,jk,1,kp)
                      IF( ln_ldfeiv )   THEN
-                        zaei_slp = aeiv(ji,jj,jk) * triadj_g(ji,jj,jk,1,kp)  
-                        zaei_slp_jp1 = aeiv(ji,jj+1,jk) * triadj_g(ji,jj+1,jk,1,kp)  
+                        zaei_slp = aeiv(ji,jj,jk) * triadj_g(ji,jj,jk,1,kp)
+                        zaei_slp_jp1 = aeiv(ji,jj+1,jk) * triadj_g(ji,jj+1,jk,1,kp)
                         zaei_slp1 = aeiv(ji,jj,jk) * triadj_g(ji,jj+1,jk,0,kp)
                      ENDIF
@@ -377 +374 @@
             ELSE
                !
-               DO kp = 0, 1
-                  DO_2D( 1, 0, 1, 0 )
+               DO kp = 0, 1               !==  Horizontal & vertical fluxes
+                  DO_2D( iij, iij-1, iij, iij-1 )
                      ze1ur = r1_e1u(ji,jj)
                      zdxt  = zdit(ji,jj,jk) * ze1ur
@@ -391 +388 @@
                      ! ln_botmix_triad is .F. mask zah for bottom half cells
                      zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp)         ! pahu(ji+ip,jj,jk)   ===>>  ????
-                     zah_ip1 = pahu(ji+1,jj,jk) * umask(ji+1,jj,jk+kp) 
+                     zah_ip1 = pahu(ji+1,jj,jk) * umask(ji+1,jj,jk+kp)
                      zah_slp  = zah * triadi(ji,jj,jk,1,kp)
                      zah_slp_ip1  = zah_ip1 * triadi(ji+1,jj,jk,1,kp)
                      zah_slp1  = zah * triadi(ji+1,jj,jk,0,kp)
                      IF( ln_ldfeiv )   THEN
-                        zaei_slp = aeiu(ji,jj,jk) * triadi_g(ji,jj,jk,1,kp) 
+                        zaei_slp = aeiu(ji,jj,jk) * triadi_g(ji,jj,jk,1,kp)
                         zaei_slp_ip1 = aeiu(ji+1,jj,jk) * triadi_g(ji+1,jj,jk,1,kp)
                         zaei_slp1 = aeiu(ji,jj,jk) * triadi_g(ji+1,jj,jk,0,kp)
@@ -406 +403 @@
                                          &      + ( zah * zdxt + zah_slp1 * zdzt_ip1 - zaei_slp1 * zdzt_ip1 ) * zbu * ze1ur  &
                                          &      )                                                                            ! bracket for halo 1 - halo 2 compatibility
-                     ztfw(ji+1,jj,jk+kp) =  ztfw(ji+1,jj,jk+kp)                                                              & 
+                     ztfw(ji+1,jj,jk+kp) =  ztfw(ji+1,jj,jk+kp)                                                              &
                                          &    - ( (zah_slp_ip1 + zaei_slp_ip1) * zdxt_ip1 * zbu_ip1 * ze3wr_ip1              &
                                          &      + ( zah_slp1 + zaei_slp1) * zdxt * zbu * ze3wr_ip1                           &
@@ -414 +411 @@
                !
                DO kp = 0, 1
-                  DO_2D( 1, 0, 1, 0 )
+                  DO_2D( iij, iij-1, iij, iij-1 )
                      ze2vr = r1_e2v(ji,jj)
                      zdyt  = zdjt(ji,jj,jk) * ze2vr
@@ -431 +428 @@
                      zah_slp_jp1 = zah_jp1 * triadj(ji,jj+1,jk,1,kp)
                      IF( ln_ldfeiv )   THEN
-                        zaei_slp = aeiv(ji,jj,jk) * triadj_g(ji,jj,jk,1,kp) 
+                        zaei_slp = aeiv(ji,jj,jk) * triadj_g(ji,jj,jk,1,kp)
                         zaei_slp_jp1 = aeiv(ji,jj+1,jk) * triadj_g(ji,jj+1,jk,1,kp)
                         zaei_slp1 = aeiv(ji,jj,jk) * triadj_g(ji,jj+1,jk,0,kp)
@@ -449 +446 @@
             ENDIF
             !                             !==  horizontal divergence and add to the general trend  ==!
-            DO_2D( 0, 0, 0, 0 )
+            DO_2D( iij-1, iij-1, iij-1, iij-1 )
                ! round brackets added to fix the order of floating point operations
                ! needed to ensure halo 1 - halo 2 compatibility
@@ -464 +461 @@
          !                                !==  add the vertical 33 flux  ==!
          IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
-            DO_3D( 0, 0, 1, 0, 2, jpkm1 )
+            DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 )
                ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)   &
                   &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
@@ -472 +469 @@
             SELECT CASE( kpass )
             CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
-               DO_3D( 0, 0, 1, 0, 2, jpkm1 )
+               DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 )
                   ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)             &
                      &                            * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
                END_3D
             CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp.
-               DO_3D( 0, 0, 1, 0, 2, jpkm1 )
+               DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                   ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)                      &
                      &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   &
@@ -485 +482 @@
          ENDIF
          !
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )      !==  Divergence of vertical fluxes added to pta  ==!
+         DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 )      !==  Divergence of vertical fluxes added to pta  ==!
             pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    &
             &                                  + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &

NEMO/trunk/src/OCE/TRA/tramle.F90

@@ -87 +87 @@
       INTEGER                     , INTENT(in   ) ::   Kmm        ! ocean time level index
       CHARACTER(len=3)            , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) ::   pu         ! in : 3 ocean transport components
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) ::   pv         ! out: same 3  transport components
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) ::   pw         !   increased by the MLE induced transport
+      ! TEMP: [tiling] Can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu         ! in : 3 ocean transport components
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pv         ! out: same 3  transport components
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pw         !   increased by the MLE induced transport
       !
       INTEGER  ::   ji, jj, jk          ! dummy loop indices
@@ -96 +97 @@
       REAL(wp) ::   zcvw, zmvw          !   -      -
       INTEGER , DIMENSION(A2D(nn_hls))     :: inml_mle
-      REAL(wp), DIMENSION(A2D(nn_hls))     :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_MH
+      REAL(wp), DIMENSION(A2D(nn_hls))     :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH
       REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zpsi_uw, zpsi_vw
-      ! TEMP: [tiling] These changes not necessary if using XIOS (subdomain support)
-      REAL(wp), DIMENSION(:,:),   ALLOCATABLE, SAVE :: zLf_NH
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zpsiu_mle, zpsiv_mle
       !!----------------------------------------------------------------------
       !
@@ -110 +108 @@
          SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
          CASE ( 0 )                                               != min of the 2 neighbour MLDs
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zhu(ji,jj) = MIN( hmle(ji+1,jj), hmle(ji,jj) )
                zhv(ji,jj) = MIN( hmle(ji,jj+1), hmle(ji,jj) )
             END_2D
          CASE ( 1 )                                               != average of the 2 neighbour MLDs
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) )
                zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) )
             END_2D
          CASE ( 2 )                                               != max of the 2 neighbour MLDs
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) )
                zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) )
@@ -126 +124 @@
          END SELECT
          IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2u(ji,jj)                                            &
                     &           * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) )   &
@@ -137 +135 @@
             !
          ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2u(ji,jj)               &
                     &                  * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) )
@@ -149 +147 @@
          !                                      !==  MLD used for MLE  ==!
          !                                                ! compute from the 10m density to deal with the diurnal cycle
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             inml_mle(ji,jj) = mbkt(ji,jj) + 1                    ! init. to number of ocean w-level (T-level + 1)
          END_2D
          IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m
-           DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )        ! from the bottom to nlb10 (10m)
+           DO_3DS( nn_hls, nn_hls, nn_hls, nn_hls, jpkm1, nlb10, -1 )        ! from the bottom to nlb10 (10m)
               IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer
            END_3D
@@ -163 +161 @@
          zbm (:,:) = 0._wp
          zn2 (:,:) = 0._wp
-         DO_3D( 1, 1, 1, 1, 1, ikmax )                    ! MLD and mean buoyancy and N2 over the mixed layer
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, ikmax )                    ! MLD and mean buoyancy and N2 over the mixed layer
             zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points
             zmld(ji,jj) = zmld(ji,jj) + zc
@@ -172 +170 @@
          SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
          CASE ( 0 )                                               != min of the 2 neighbour MLDs
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) )
                zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) )
             END_2D
          CASE ( 1 )                                               != average of the 2 neighbour MLDs
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp
                zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp
             END_2D
          CASE ( 2 )                                               != max of the 2 neighbour MLDs
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) )
                zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) )
@@ -188 +186 @@
          END SELECT
          !                                                ! convert density into buoyancy
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             zbm(ji,jj) = + grav * zbm(ji,jj) / MAX( e3t(ji,jj,1,Kmm), zmld(ji,jj) )
          END_2D
@@ -201 +199 @@
          !
          IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2_e1u(ji,jj)                                            &
                     &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )   &
@@ -212 +210 @@
             !
          ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2_e1u(ji,jj)               &
                     &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )
@@ -222 +220 @@
          !
          IF( nn_conv == 1 ) THEN              ! No MLE in case of convection
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp
                IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp
@@ -230 +228 @@
       ENDIF  ! end of ln_osm_mle conditional
     !                                      !==  structure function value at uw- and vw-points  ==!
-    DO_2D( 1, 0, 1, 0 )
+    DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
        zhu(ji,jj) = 1._wp / MAX(zhu(ji,jj), rsmall)                   ! hu --> 1/hu
        zhv(ji,jj) = 1._wp / MAX(zhv(ji,jj), rsmall) 
@@ -238 +236 @@
     zpsi_vw(:,:,:) = 0._wp
     !
-      DO_3D( 1, 0, 1, 0, 2, ikmax )                ! start from 2 : surface value = 0
+      DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 2, ikmax )                ! start from 2 : surface value = 0
+      
          zcuw = 1._wp - ( gdepw(ji+1,jj,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhu(ji,jj)
          zcvw = 1._wp - ( gdepw(ji,jj+1,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhv(ji,jj)
@@ -252 +251 @@
       !                                      !==  transport increased by the MLE induced transport ==!
       DO jk = 1, ikmax
-         DO_2D( 1, 0, 1, 0 )                      ! CAUTION pu,pv must be defined at row/column i=1 / j=1
+         DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
             pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
             pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
          END_2D
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk)   &
                &                          + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) ) * wmask(ji,jj,1)
@@ -262 +261 @@
       END DO
 
-      ! TEMP: [tiling] These changes not necessary if using XIOS (subdomain support)
       IF( cdtype == 'TRA') THEN              !==  outputs  ==!
-         IF( ntile == 0 .OR. ntile == 1 ) THEN                             ! Do only on the first tile
-            ALLOCATE( zLf_NH(jpi,jpj), zpsiu_mle(jpi,jpj,jpk), zpsiv_mle(jpi,jpj,jpk) )
-            zpsiu_mle(:,:,:) = 0._wp ; zpsiv_mle(:,:,:) = 0._wp
-         ENDIF
          !
          IF (ln_osm_mle.and.ln_zdfosm) THEN
@@ -279 +273 @@
          ENDIF
          !
+         CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f
+         !
          ! divide by cross distance to give streamfunction with dimensions m^2/s
          DO_3D( 0, 0, 0, 0, 1, ikmax+1 )
-            zpsiu_mle(ji,jj,jk) = zpsi_uw(ji,jj,jk) * r1_e2u(ji,jj)
-            zpsiv_mle(ji,jj,jk) = zpsi_vw(ji,jj,jk) * r1_e1v(ji,jj)
+            zpsi_uw(ji,jj,jk) = zpsi_uw(ji,jj,jk) * r1_e2u(ji,jj)
+            zpsi_vw(ji,jj,jk) = zpsi_vw(ji,jj,jk) * r1_e1v(ji,jj)
          END_3D
-
-         IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
-            CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f
-            CALL iom_put( "psiu_mle", zpsiu_mle )    ! i-mle streamfunction
-            CALL iom_put( "psiv_mle", zpsiv_mle )    ! j-mle streamfunction
-            DEALLOCATE( zLf_NH, zpsiu_mle, zpsiv_mle )
-         ENDIF
+         CALL iom_put( "psiu_mle", zpsi_uw )    ! i-mle streamfunction
+         CALL iom_put( "psiv_mle", zpsi_vw )    ! j-mle streamfunction
       ENDIF
       !
@@ -375 +366 @@
          r1_ft(:,:) = 1._wp / SQRT(  ff_t(:,:) * ff_t(:,:) + z1_t2  )
          !
+         ! Specifically, dbdx_mle, dbdy_mle and mld_prof need to be defined for nn_hls = 2
+         IF( nn_hls == 2 .AND. ln_osm_mle .AND. ln_zdfosm ) THEN
+            CALL ctl_stop('nn_hls = 2 cannot be used with ln_mle = ln_osm_mle = ln_zdfosm = T (zdfosm not updated for nn_hls = 2)')
+         ENDIF
       ENDIF
       !

NEMO/trunk/src/OCE/TRA/tranpc.F90

@@ -17 +17 @@
    USE oce            ! ocean dynamics and active tracers
    USE dom_oce        ! ocean space and time domain
-   ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed)
-   USE domtile
    USE phycst         ! physical constants
    USE zdf_oce        ! ocean vertical physics
@@ -81 +79 @@
       LOGICAL, PARAMETER :: l_LB_debug = .FALSE. ! set to true if you want to follow what is
       INTEGER :: ilc1, jlc1, klc1, nncpu         ! actually happening in a water column at point "ilc1, jlc1"
-      INTEGER :: isi, isj, iei, iej
       LOGICAL :: lp_monitor_point = .FALSE.      ! in CPU domain "nncpu"
       !!----------------------------------------------------------------------
@@ -105 +102 @@
          CALL bn2    ( pts(:,:,:,:,Kaa), zab, zn2, Kmm )    ! after Brunt-Vaisala  (given on W-points)
          !
-         IF( ntile == 0 .OR. ntile == 1 ) nnpcc = 0         ! Do only on the first tile
-         !
-         IF( ntsi == Nis0 ) THEN ; isi = nn_hls ; ELSE ; isi = 0 ; ENDIF    ! Avoid double-counting when using tiling
-         IF( ntsj == Njs0 ) THEN ; isj = nn_hls ; ELSE ; isj = 0 ; ENDIF
-         IF( ntei == Nie0 ) THEN ; iei = nn_hls ; ELSE ; iei = 0 ; ENDIF
-         IF( ntej == Nje0 ) THEN ; iej = nn_hls ; ELSE ; iej = 0 ; ENDIF
-         !
-         DO_2D( isi, iei, isj, iej )                        ! interior column only
+         IF( .NOT. l_istiled .OR. ntile == 1 ) nnpcc = 0         ! Do only on the first tile
+         !
+         DO_2D_OVR( 0, 0, 0, 0 )                        ! interior column only
             !
             IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points
@@ -319 +311 @@
          ENDIF
          !
-         IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
+         IF( .NOT. l_istiled .OR. ntile == nijtile )  THEN                ! Do only for the full domain
             IF( lwp .AND. l_LB_debug ) THEN
                WRITE(numout,*) 'Exiting tra_npc , kt = ',kt,', => numb. of statically instable water-columns: ', nnpcc

NEMO/trunk/src/OCE/TRA/traqsr.F90

@@ -108 +108 @@
       !
       INTEGER  ::   ji, jj, jk               ! dummy loop indices
-      INTEGER  ::   irgb, isi, iei, isj, iej ! local integers
+      INTEGER  ::   irgb                     ! local integers
       REAL(wp) ::   zchl, zcoef, z1_2        ! local scalars
       REAL(wp) ::   zc0 , zc1 , zc2 , zc3    !    -         -
@@ -121 +121 @@
       IF( ln_timing )   CALL timing_start('tra_qsr')
       !
-      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
          IF( kt == nit000 ) THEN
             IF(lwp) WRITE(numout,*)
@@ -137 +137 @@
       !                         !  before qsr induced heat content  !
       !                         !-----------------------------------!
-      IF( ntsi == Nis0 ) THEN ; isi = nn_hls ; ELSE ; isi = 0 ; ENDIF    ! Avoid double-counting when using tiling
-      IF( ntsj == Njs0 ) THEN ; isj = nn_hls ; ELSE ; isj = 0 ; ENDIF
-      IF( ntei == Nie0 ) THEN ; iei = nn_hls ; ELSE ; iei = 0 ; ENDIF
-      IF( ntej == Nje0 ) THEN ; iej = nn_hls ; ELSE ; iej = 0 ; ENDIF
-
       IF( kt == nit000 ) THEN          !==  1st time step  ==!
          IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN    ! read in restart
             z1_2 = 0.5_wp
-            IF( ntile == 0 .OR. ntile == 1 )  THEN                        ! Do only on the first tile
+            IF( .NOT. l_istiled .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 )   ! before heat content trend due to Qsr flux
@@ -151 +146 @@
          ELSE                                           ! No restart or Euler forward at 1st time step
             z1_2 = 1._wp
-            DO_3D( isi, iei, isj, iej, 1, jpk )
+            DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
                qsr_hc_b(ji,jj,jk) = 0._wp
             END_3D
@@ -157 +152 @@
       ELSE                             !==  Swap of qsr heat content  ==!
          z1_2 = 0.5_wp
-         DO_3D( isi, iei, isj, iej, 1, jpk )
+         DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
             qsr_hc_b(ji,jj,jk) = qsr_hc(ji,jj,jk)
          END_3D
@@ -168 +163 @@
       CASE( np_BIO )                   !==  bio-model fluxes  ==!
          !
-         DO_3D( isi, iei, isj, iej, 1, nksr )
+         DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr )
             qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) )
          END_3D
@@ -179 +174 @@
          !
          IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll
-            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
+            IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                                         ! Do only for the full domain
+               IF( ln_tile ) CALL dom_tile_stop( ldhold=.TRUE. )             ! 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
+               IF( ln_tile ) CALL dom_tile_start( ldhold=.TRUE. )            ! Revert to tile domain
             ENDIF
             !
@@ -190 +185 @@
             ! most expensive calculations)
             !
-            DO_2D( isi, iei, isj, iej )
+            DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
                        ! zlogc = log(zchl)
                zlogc = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) )
@@ -209 +204 @@
 
 !
-            DO_3D( isi, iei, isj, iej, 1, nksr + 1 )
+            DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr + 1 )
                ! zchl    = ALOG( ze0(ji,jj) )
                zlogc = ze0(ji,jj)
@@ -239 +234 @@
          !
          zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B
-         DO_2D( isi, iei, isj, iej )
+         DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
             ze0(ji,jj) = rn_abs * qsr(ji,jj)
             ze1(ji,jj) = zcoef  * qsr(ji,jj)
@@ -250 +245 @@
          !
          !                                    !* interior equi-partition in R-G-B depending on vertical profile of Chl
-         DO_3D( isi, iei, isj, iej, 2, nksr + 1 )
+         DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 2, nksr + 1 )
             ze3t = e3t(ji,jj,jk-1,Kmm)
             irgb = NINT( ztmp3d(ji,jj,jk) )
@@ -264 +259 @@
          END_3D
          !
-         DO_3D( isi, iei, isj, iej, 1, nksr )          !* now qsr induced heat content
+         DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr )          !* now qsr induced heat content
             qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) )
          END_3D
@@ -274 +269 @@
          zz0 =        rn_abs   * r1_rho0_rcp      ! surface equi-partition in 2-bands
          zz1 = ( 1. - rn_abs ) * r1_rho0_rcp
-         DO_3D( isi, iei, isj, iej, 1, nksr )          !* now qsr induced heat content
+         DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr )          !* now qsr induced heat content
             zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r )
             zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r )
@@ -292 +287 @@
       !
       ! sea-ice: store the 1st ocean level attenuation coefficient
-      DO_2D( isi, iei, isj, iej )
+      DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
          IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) )
          ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
@@ -298 +293 @@
       END_2D
       !
-      ! TEMP: [tiling] This change not necessary and working array can use A2D(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( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution
+         ALLOCATE( zetot(A2D(nn_hls),jpk) )
+         zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero
+         DO_3DS(0, 0, 0, 0, nksr, 1, -1)
+            zetot(ji,jj,jk) = zetot(ji,jj,jk+1) + qsr_hc(ji,jj,jk) * rho0_rcp
+         END_3D
+         CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation
+         DEALLOCATE( zetot )
+      ENDIF
+      !
+      IF( .NOT. l_istiled .OR. ntile == nijtile )  THEN                ! Do only on the last tile
          IF( lrst_oce ) THEN     ! write in the ocean restart file
             CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc      )

NEMO/trunk/src/OCE/TRA/trasbc.F90

@@ -77 +77 @@
       !
       INTEGER  ::   ji, jj, jk, jn               ! dummy loop indices
-      INTEGER  ::   ikt, ikb, isi, iei, isj, iej ! local integers
+      INTEGER  ::   ikt, ikb                     ! local integers
       REAL(wp) ::   zfact, z1_e3t, zdep, ztim    ! local scalar
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdt, ztrds
@@ -84 +84 @@
       IF( ln_timing )   CALL timing_start('tra_sbc')
       !
-      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
          IF( kt == nit000 ) THEN
             IF(lwp) WRITE(numout,*)
@@ -98 +98 @@
       ENDIF
       !
-      IF( ntsi == Nis0 ) THEN ; isi = nn_hls ; ELSE ; isi = 0 ; ENDIF    ! Avoid double-counting when using tiling
-      IF( ntsj == Njs0 ) THEN ; isj = nn_hls ; ELSE ; isj = 0 ; ENDIF
-      IF( ntei == Nie0 ) THEN ; iei = nn_hls ; ELSE ; iei = 0 ; ENDIF
-      IF( ntej == Nje0 ) THEN ; iej = nn_hls ; ELSE ; iej = 0 ; 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
-         DO_2D( isi, iei, isj, iej )
+         DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
             qns(ji,jj) = qns(ji,jj) + qsr(ji,jj)      ! total heat flux in qns
             qsr(ji,jj) = 0._wp                        ! qsr set to zero
@@ -118 +113 @@
          IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN      ! Restart: read in restart file
             zfact = 0.5_wp
-            IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+            IF( .NOT. l_istiled .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
@@ -126 +121 @@
          ELSE                                             ! No restart or restart not found: Euler forward time stepping
             zfact = 1._wp
-            DO_2D( isi, iei, isj, iej )
+            DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
                sbc_tsc(ji,jj,:) = 0._wp
                sbc_tsc_b(ji,jj,:) = 0._wp
@@ -133 +128 @@
       ELSE                                !* other time-steps: swap of forcing fields
          zfact = 0.5_wp
-         DO_2D( isi, iei, isj, iej )
+         DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
             sbc_tsc_b(ji,jj,:) = sbc_tsc(ji,jj,:)
          END_2D
       ENDIF
       !                             !==  Now sbc tracer content fields  ==!
-      DO_2D( isi, iei, isj, iej )
+      DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
          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( isi, iei, isj, iej )                    !==>> add concentration/dilution effect due to constant volume cell
+         DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )                    !==>> add concentration/dilution effect due to constant volume cell
             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                                 !==>> output c./d. term
-         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
+         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
       !
@@ -160 +153 @@
       END DO
       !
-      IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile
+      IF( .NOT. l_istiled .OR. ntile == nijtile )  THEN                ! Do only on the last tile
          IF( lrst_oce ) THEN           !==  write sbc_tsc in the ocean restart file  ==!
             CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem) )
@@ -186 +179 @@
       ENDIF
 
-      IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last 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( 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 defined key_asminc

NEMO/trunk/src/OCE/TRA/trazdf.F90

@@ -64 +64 @@
       !
       IF( kt == nit000 )  THEN
-         IF( ntile == 0 .OR. ntile == 1 )  THEN                   ! Do only on the first tile
+         IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                   ! Do only on the first tile
             IF(lwp)WRITE(numout,*)
             IF(lwp)WRITE(numout,*) 'tra_zdf : implicit vertical mixing on T & S'

NEMO/trunk/src/OCE/TRA/zpshde.F90

@@ -47 +47 @@
       INTEGER                     , INTENT(in   )           ::  Kmm         ! ocean time level index
       INTEGER                     , INTENT(in   )           ::  kjpt        ! number of tracers
-      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout)           ::  pta         ! 4D tracers fields
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   )           ::  pta         ! 4D tracers fields
       REAL(wp), DIMENSION(:,:,:)  , INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts
-      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout), OPTIONAL ::  prd         ! 3D density anomaly fields
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields
       REAL(wp), DIMENSION(:,:)    , INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad of prd at u- & v-pts (bottom)
       !
@@ -111 +111 @@
       INTEGER                                , INTENT(in   )           ::  kjpt        ! number of tracers
       INTEGER                                , INTENT(in   )           ::  ktta, ktgt, ktrd, ktgr
-      REAL(wp), DIMENSION(A2D_T(ktta),JPK,KJPT), INTENT(inout)           ::  pta         ! 4D tracers fields
+      REAL(wp), DIMENSION(A2D_T(ktta),JPK,KJPT), INTENT(in   )           ::  pta         ! 4D tracers fields
       REAL(wp), DIMENSION(A2D_T(ktgt)    ,KJPT), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts
-      REAL(wp), DIMENSION(A2D_T(ktrd),JPK     ), INTENT(inout), OPTIONAL ::  prd         ! 3D density anomaly fields
+      REAL(wp), DIMENSION(A2D_T(ktrd),JPK     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields
       REAL(wp), DIMENSION(A2D_T(ktgr)         ), INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad of prd at u- & v-pts (bottom)
       !
@@ -124 +124 @@
       !
       IF( ln_timing )   CALL timing_start( 'zps_hde')
-      IF (nn_hls.EQ.2) THEN
-         CALL lbc_lnk( 'zpshde', pta, 'T', 1.0_wp)
-         IF(PRESENT(prd)) CALL lbc_lnk( 'zpshde', prd, 'T', 1.0_wp)
-      END IF
       !
       pgtu(:,:,:) = 0._wp   ;   zti (:,:,:) = 0._wp   ;   zhi (:,:) = 0._wp
@@ -134 +130 @@
       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
          !
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )              ! Gradient of density at the last level
+         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )              ! Gradient of density at the last level
             iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
             ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
@@ -173 +169 @@
       END DO
       !
-      IF (nn_hls.EQ.1) CALL lbc_lnk( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp )   ! Lateral boundary cond.
+      IF (nn_hls==1) CALL lbc_lnk( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp )   ! Lateral boundary cond.
       !
       IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
@@ -206 +202 @@
             ENDIF
          END_2D
-         IF (nn_hls.EQ.1) CALL lbc_lnk( 'zpshde', pgru , 'U', -1.0_wp , pgrv , 'V', -1.0_wp )   ! Lateral boundary conditions
+         IF (nn_hls==1) CALL lbc_lnk( 'zpshde', pgru , 'U', -1.0_wp , pgrv , 'V', -1.0_wp )   ! Lateral boundary conditions
          !
       END IF
@@ -221 +217 @@
       INTEGER                     , INTENT(in   )           ::  Kmm          ! ocean time level index
       INTEGER                     , INTENT(in   )           ::  kjpt         ! number of tracers
-      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout)           ::  pta          ! 4D tracers fields
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   )           ::  pta          ! 4D tracers fields
       REAL(wp), DIMENSION(:,:,:)  , INTENT(  out)           ::  pgtu, pgtv   ! hor. grad. of ptra at u- & v-pts
       REAL(wp), DIMENSION(:,:,:)  , INTENT(  out)           ::  pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF)
-      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout), OPTIONAL ::  prd          ! 3D density anomaly fields
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ), OPTIONAL ::  prd          ! 3D density anomaly fields
       REAL(wp), DIMENSION(:,:)    , INTENT(  out), OPTIONAL ::  pgru, pgrv   ! hor. grad of prd at u- & v-pts (bottom)
       REAL(wp), DIMENSION(:,:)    , INTENT(  out), OPTIONAL ::  pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top)
@@ -291 +287 @@
       INTEGER                                , INTENT(in   )           ::  kjpt         ! number of tracers
       INTEGER                                , INTENT(in   )           ::  ktta, ktgt, ktgti, ktrd, ktgr, ktgri
-      REAL(wp), DIMENSION(A2D_T(ktta),JPK,KJPT), INTENT(inout)           ::  pta          ! 4D tracers fields
+      REAL(wp), DIMENSION(A2D_T(ktta),JPK,KJPT), INTENT(in   )           ::  pta          ! 4D tracers fields
       REAL(wp), DIMENSION(A2D_T(ktgt)    ,KJPT), INTENT(  out)           ::  pgtu, pgtv   ! hor. grad. of ptra at u- & v-pts
       REAL(wp), DIMENSION(A2D_T(ktgti)   ,KJPT), INTENT(  out)           ::  pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF)
-      REAL(wp), DIMENSION(A2D_T(ktrd),JPK     ), INTENT(inout), OPTIONAL ::  prd          ! 3D density anomaly fields
+      REAL(wp), DIMENSION(A2D_T(ktrd),JPK     ), INTENT(in   ), OPTIONAL ::  prd          ! 3D density anomaly fields
       REAL(wp), DIMENSION(A2D_T(ktgr)         ), INTENT(  out), OPTIONAL ::  pgru, pgrv   ! hor. grad of prd at u- & v-pts (bottom)
       REAL(wp), DIMENSION(A2D_T(ktgri)        ), INTENT(  out), OPTIONAL ::  pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top)
@@ -307 +303 @@
       IF( ln_timing )   CALL timing_start( 'zps_hde_isf')
       !
-      IF (nn_hls.EQ.2) THEN
-         CALL lbc_lnk( 'zpshde', pta, 'T', 1.0_wp)
-         IF (PRESENT(prd)) CALL lbc_lnk( 'zpshde', prd, 'T', 1.0_wp)
-      END IF
-
       pgtu (:,:,:) = 0._wp   ;   pgtv (:,:,:) =0._wp
       pgtui(:,:,:) = 0._wp   ;   pgtvi(:,:,:) =0._wp
@@ -319 +310 @@
       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
          !
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
 
             iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
@@ -359 +350 @@
       END DO
       !
-      IF (nn_hls.EQ.1) CALL lbc_lnk( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp )   ! Lateral boundary cond.
+      IF (nn_hls==1) CALL lbc_lnk( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp )   ! Lateral boundary cond.
 
       ! horizontal derivative of density anomalies (rd)
@@ -401 +392 @@
          END_2D
 
-         IF (nn_hls.EQ.1) CALL lbc_lnk( 'zpshde', pgru , 'U', -1.0_wp , pgrv , 'V', -1.0_wp )   ! Lateral boundary conditions
+         IF (nn_hls==1) CALL lbc_lnk( 'zpshde', pgru , 'U', -1.0_wp , pgrv , 'V', -1.0_wp )   ! Lateral boundary conditions
          !
       END IF
@@ -408 +399 @@
       !
       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            !
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
             iku = miku(ji,jj); ikup1 = miku(ji,jj) + 1
             ikv = mikv(ji,jj); ikvp1 = mikv(ji,jj) + 1
@@ -452 +443 @@
          !
       END DO
-      IF (nn_hls.EQ.1) CALL lbc_lnk( 'zpshde', pgtui(:,:,:), 'U', -1.0_wp , pgtvi(:,:,:), 'V', -1.0_wp )   ! Lateral boundary cond.
+      IF (nn_hls==1) CALL lbc_lnk( 'zpshde', pgtui(:,:,:), 'U', -1.0_wp , pgtvi(:,:,:), 'V', -1.0_wp )   ! Lateral boundary cond.
 
       IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
@@ -491 +482 @@
 
          END_2D
-         IF (nn_hls.EQ.1) CALL lbc_lnk( 'zpshde', pgrui, 'U', -1.0_wp , pgrvi, 'V', -1.0_wp )   ! Lateral boundary conditions
+         IF (nn_hls==1) CALL lbc_lnk( 'zpshde', pgrui, 'U', -1.0_wp , pgrvi, 'V', -1.0_wp )   ! Lateral boundary conditions
          !
       END IF