Index: NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/ICE/icedyn_adv_pra.F90
===================================================================
--- NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/ICE/icedyn_adv_pra.F90	(revision 13925)
+++ NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/ICE/icedyn_adv_pra.F90	(revision 13946)
@@ -273,4 +273,20 @@
          ! --- Lateral boundary conditions --- !
          !     caution: for gradients (sx and sy) the sign changes
+#if defined key_mpi3
+         CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', z0ice , 'T', 1._wp, sxice , 'T', -1._wp, syice , 'T', -1._wp  & ! ice volume
+            &                                , sxxice, 'T', 1._wp, syyice, 'T',  1._wp, sxyice, 'T',  1._wp  &
+            &                                , z0snw , 'T', 1._wp, sxsn  , 'T', -1._wp, sysn  , 'T', -1._wp  & ! snw volume
+            &                                , sxxsn , 'T', 1._wp, syysn , 'T',  1._wp, sxysn , 'T',  1._wp  )
+         CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', z0smi , 'T', 1._wp, sxsal , 'T', -1._wp, sysal , 'T', -1._wp  & ! ice salinity
+            &                                , sxxsal, 'T', 1._wp, syysal, 'T',  1._wp, sxysal, 'T',  1._wp  &
+            &                                , z0ai  , 'T', 1._wp, sxa   , 'T', -1._wp, sya   , 'T', -1._wp  & ! ice concentration
+            &                                , sxxa  , 'T', 1._wp, syya  , 'T',  1._wp, sxya  , 'T',  1._wp  )
+         CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', z0oi  , 'T', 1._wp, sxage , 'T', -1._wp, syage , 'T', -1._wp  & ! ice age
+            &                                , sxxage, 'T', 1._wp, syyage, 'T',  1._wp, sxyage, 'T',  1._wp  )
+         CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', z0es  , 'T', 1._wp, sxc0  , 'T', -1._wp, syc0  , 'T', -1._wp  & ! snw enthalpy
+            &                                , sxxc0 , 'T', 1._wp, syyc0 , 'T',  1._wp, sxyc0 , 'T',  1._wp  ) 
+         CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', z0ei  , 'T', 1._wp, sxe   , 'T', -1._wp, sye   , 'T', -1._wp  & ! ice enthalpy
+            &                                , sxxe  , 'T', 1._wp, syye  , 'T',  1._wp, sxye  , 'T',  1._wp  )
+#else
          CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ice , 'T', 1._wp, sxice , 'T', -1._wp, syice , 'T', -1._wp  & ! ice volume
             &                                , sxxice, 'T', 1._wp, syyice, 'T',  1._wp, sxyice, 'T',  1._wp  &
@@ -287,12 +303,25 @@
          CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ei  , 'T', 1._wp, sxe   , 'T', -1._wp, sye   , 'T', -1._wp  & ! ice enthalpy
             &                                , sxxe  , 'T', 1._wp, syye  , 'T',  1._wp, sxye  , 'T',  1._wp  )
+#endif
          IF ( ln_pnd_LEV ) THEN
+#if defined key_mpi3
+            CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp  & ! melt pond fraction
+               &                                , sxxap, 'T', 1._wp, syyap, 'T',  1._wp, sxyap, 'T',  1._wp  &
+               &                                , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp  & ! melt pond volume
+               &                                , sxxvp, 'T', 1._wp, syyvp, 'T',  1._wp, sxyvp, 'T',  1._wp  ) 
+#else
             CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp  & ! melt pond fraction
                &                                , sxxap, 'T', 1._wp, syyap, 'T',  1._wp, sxyap, 'T',  1._wp  &
                &                                , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp  & ! melt pond volume
                &                                , sxxvp, 'T', 1._wp, syyvp, 'T',  1._wp, sxyvp, 'T',  1._wp  ) 
+#endif
             IF ( ln_pnd_lids ) THEN
+#if defined key_mpi3
+               CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', z0vl ,'T', 1._wp, sxvl ,'T', -1._wp, syvl ,'T', -1._wp  & ! melt pond lid volume
+                  &                                , sxxvl,'T', 1._wp, syyvl,'T',  1._wp, sxyvl,'T',  1._wp  ) 
+#else
                CALL lbc_lnk_multi( 'icedyn_adv_pra', z0vl ,'T', 1._wp, sxvl ,'T', -1._wp, syvl ,'T', -1._wp  & ! melt pond lid volume
                   &                                , sxxvl,'T', 1._wp, syyvl,'T',  1._wp, sxyvl,'T',  1._wp  ) 
+#endif
             ENDIF
          ENDIF
Index: NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/ICE/icedyn_adv_umx.F90
===================================================================
--- NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/ICE/icedyn_adv_umx.F90	(revision 13925)
+++ NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/ICE/icedyn_adv_umx.F90	(revision 13946)
@@ -368,14 +368,33 @@
          ! --- Lateral boundary conditions --- !
          IF    ( ln_pnd_LEV .AND. ln_pnd_lids ) THEN
+#if defined key_mpi3
+            CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp &
+               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp )
+#else
             CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp &
                &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp )
+#endif
          ELSEIF( ln_pnd_LEV .AND. .NOT.ln_pnd_lids ) THEN
+#if defined key_mpi3
+            CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp &
+               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp )
+#else
             CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp &
                &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp )
+#endif
          ELSE
+#if defined key_mpi3
+            CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp )
+#else
             CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp )
+#endif
          ENDIF
+#if defined key_mpi3
+         CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', pe_i, 'T', 1._wp )
+         CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', pe_s, 'T', 1._wp )
+#else
          CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp )
          CALL lbc_lnk( 'icedyn_adv_umx', pe_s, 'T', 1._wp )
+#endif
          !
          !== Open water area ==!
Index: NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/LBC/lbc_lnk_nc_generic.h90
===================================================================
--- NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/LBC/lbc_lnk_nc_generic.h90	(revision 13925)
+++ NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/LBC/lbc_lnk_nc_generic.h90	(revision 13946)
@@ -58,7 +58,7 @@
       !!
       INTEGER                          ::   kfld        ! number of elements that will be attributed
-      PTR_TYPE         , DIMENSION(11) ::   ptab_ptr    ! pointer array
-      CHARACTER(len=1) , DIMENSION(11) ::   cdna_ptr    ! nature of ptab_ptr grid-points
-      REAL(wp)         , DIMENSION(11) ::   psgn_ptr    ! sign used across the north fold boundary
+      PTR_TYPE         , DIMENSION(16) ::   ptab_ptr    ! pointer array
+      CHARACTER(len=1) , DIMENSION(16) ::   cdna_ptr    ! nature of ptab_ptr grid-points
+      REAL(wp)         , DIMENSION(16) ::   psgn_ptr    ! sign used across the north fold boundary
       !!---------------------------------------------------------------------
       !
@@ -79,9 +79,9 @@
       IF( PRESENT(psgn10) )   CALL ROUTINE_NC_LOAD( pt10, cdna10, psgn10, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
       IF( PRESENT(psgn11) )   CALL ROUTINE_NC_LOAD( pt11, cdna11, psgn11, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
-      IF( PRESENT(psgn11) )   CALL ROUTINE_NC_LOAD( pt12, cdna12, psgn12, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
-      IF( PRESENT(psgn11) )   CALL ROUTINE_NC_LOAD( pt13, cdna13, psgn13, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
-      IF( PRESENT(psgn11) )   CALL ROUTINE_NC_LOAD( pt14, cdna14, psgn14, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
-      IF( PRESENT(psgn11) )   CALL ROUTINE_NC_LOAD( pt15, cdna15, psgn15, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
-      IF( PRESENT(psgn11) )   CALL ROUTINE_NC_LOAD( pt16, cdna16, psgn16, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
+      IF( PRESENT(psgn12) )   CALL ROUTINE_NC_LOAD( pt12, cdna12, psgn12, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
+      IF( PRESENT(psgn13) )   CALL ROUTINE_NC_LOAD( pt13, cdna13, psgn13, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
+      IF( PRESENT(psgn14) )   CALL ROUTINE_NC_LOAD( pt14, cdna14, psgn14, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
+      IF( PRESENT(psgn15) )   CALL ROUTINE_NC_LOAD( pt15, cdna15, psgn15, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
+      IF( PRESENT(psgn16) )   CALL ROUTINE_NC_LOAD( pt16, cdna16, psgn16, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
       !
       CALL lbc_lnk_nc    ( cdname,               ptab_ptr, cdna_ptr, psgn_ptr, kfld, kfillmode, pfillval, lsend, lrecv, ncsten )
Index: NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/TRA/traadv.F90
===================================================================
--- NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/TRA/traadv.F90	(revision 13925)
+++ NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/TRA/traadv.F90	(revision 13946)
@@ -25,5 +25,7 @@
    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)
@@ -177,4 +179,9 @@
          !
          CASE ( np_CEN )                                 ! Centered scheme : 2nd / 4th order
+#if defined key_mpi3
+            IF (nn_hls.EQ.2) CALL lbc_lnk_nc_multi( 'traadv', pts(:,:,:,:,Kmm), 'T', 1. )
+#else
+            IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kmm), 'T', 1. )
+#endif
             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
@@ -187,19 +194,33 @@
                CALL lbc_lnk_multi( 'traadv', zuu(:,:,:), 'U', -1., zvv(:,:,:), 'V', -1., zww(:,:,:), 'W', 1.)
 #endif
+#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
-            CALL tra_adv_fct    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v )
          CASE ( np_MUS )                                 ! MUSCL
             ! NOTE: [tiling-comms-merge] I added this lbc_lnk as it did not validate against the trunk when using ln_zco
+            IF (nn_hls.EQ.2) THEN 
+#if defined key_mpi3
+                CALL lbc_lnk_nc_multi( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)
+#else
+                CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)
+#endif
+#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 defined key_mpi3
             IF (nn_hls.EQ.2) CALL lbc_lnk_nc_multi( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)
 #else
-            IF (nn_hls.EQ.2) CALL lbc_lnk_multi( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)
-#endif
-            CALL tra_adv_mus    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )
-         CASE ( np_UBS )                                 ! UBS
-#if defined key_mpi3
-            IF (nn_hls.EQ.2) CALL lbc_lnk_nc_multi( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)
-#else
-            IF (nn_hls.EQ.2) CALL lbc_lnk_multi( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)
+            IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)
 #endif
             CALL tra_adv_ubs    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v   )
Index: NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/TRA/traadv_fct.F90
===================================================================
--- NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/TRA/traadv_fct.F90	(revision 13925)
+++ NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/TRA/traadv_fct.F90	(revision 13946)
@@ -34,4 +34,6 @@
    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
@@ -244,5 +246,5 @@
 #endif
             !
-            DO_3D( 1, 0, 1, 0, 1, jpkm1 )    ! Horizontal advective fluxes
+            DO_3D( 1, 0, 1, 0, 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)
@@ -269,4 +271,6 @@
             IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
 #endif
+            !
+            IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
             !
             DO_3D( 0, 0, 0, 0, 1, jpkm1 )    ! Horizontal advective fluxes
@@ -319,4 +323,10 @@
             CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp)
 #endif
+         END IF
+         !
+         IF (nn_hls.EQ.1) THEN
+            CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1.0_wp, zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp, zwz, 'T', 1.0_wp )
+         ELSE
+            CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp)
          END IF
          !
Index: NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/TRA/traadv_fct_lf.F90
===================================================================
--- NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/TRA/traadv_fct_lf.F90	(revision 13946)
+++ NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/TRA/traadv_fct_lf.F90	(revision 13946)
@@ -0,0 +1,553 @@
+MODULE traadv_fct_lf
+   !!==============================================================================
+   !!                       ***  MODULE  traadv_fct  ***
+   !! Ocean  tracers:  horizontal & vertical advective trend (2nd/4th order Flux Corrected Transport method)
+   !!==============================================================================
+   !! History :  3.7  !  2015-09  (L. Debreu, G. Madec)  original code (inspired from traadv_tvd.F90)
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!  tra_adv_fct-lf : update the tracer trend with a 3D advective trends using a 2nd or 4th order FCT scheme
+   !!                   with sub-time-stepping in the vertical direction - loop fusion version
+   !!  nonosc_lf      : compute monotonic tracer fluxes by a non-oscillatory algorithm - loop fusion version 
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and active tracers
+   USE dom_oce        ! ocean space and time domain
+   USE trc_oce        ! share passive tracers/Ocean variables
+   USE trd_oce        ! trends: ocean variables
+   USE trdtra         ! tracers trends
+   USE diaptr         ! poleward transport diagnostics
+   USE diaar5         ! AR5 diagnostics
+   USE phycst  , ONLY : rho0_rcp
+   USE zdf_oce , ONLY : ln_zad_Aimp
+   !
+   USE in_out_manager ! I/O manager
+   USE iom            ! 
+   USE lib_mpp        ! MPP library
+   USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE traadv_fct
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   tra_adv_fct_lf     ! called by traadv.F90
+
+   LOGICAL  ::   l_trd   ! flag to compute trends
+   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
+   LOGICAL  ::   l_hst   ! flag to compute heat/salt transport
+   REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
+
+   !                                        ! tridiag solver associated indices:
+   INTEGER, PARAMETER ::   np_NH   = 0   ! Neumann homogeneous boundary condition
+   INTEGER, PARAMETER ::   np_CEN2 = 1   ! 2nd order centered  boundary condition
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+
+#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) )
+
+#define search_in_neighbour(out,OP,vec,ji,jj,jk) \
+        out = OP(vec(ji,jj,jk),vec(ji-1,jj,jk),vec(ji+1,jj,jk), \
+        vec(ji,jj-1,jk),vec(ji,jj+1,jk),vec(ji,jj,MAX(jk-1,1)),vec(ji,jj,jk+1))
+
+#define pos_part_of_flux(ji,jj,jk,out) \
+        out = MAX(0.,paa_in(ji-1,jj,jk)) - MIN(0.,paa_in(ji,jj,jk)) \
+        + MAX(0.,pbb_in(ji,jj-1,jk)) - MIN(0.,pbb_in(ji,jj,jk)) \
+        + MAX(0.,pcc_in(ji,jj,jk+1)) - MIN(0.,pcc_in(ji,jj,jk)) 
+
+#define neg_part_of_flux(ji,jj,jk,out) \
+        out = MAX( 0.,paa_in(ji,jj,jk) ) - MIN( 0., paa_in(ji-1,jj,jk)) \
+        + MAX( 0.,pbb_in(ji,jj,jk) ) - MIN( 0., pbb_in(ji,jj-1,jk)) \
+        + MAX( 0.,pcc_in(ji,jj,jk) ) - MIN( 0., pcc_in(ji,jj,jk+1))
+
+#define beta_terms(bt,betup,betdo,up,pos,do,neg,ji,jj,jk) \
+        bt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt ; \
+        betup = ( up            - paft(ji,jj,jk) ) / ( pos + zrtrn ) * bt ; \
+        betdo = ( paft(ji,jj,jk) - do            ) / ( neg + zrtrn ) * bt
+
+#define monotonic_flux(a,b,c,betup_p1,betdo_p1,vec,vec_in,jk) \
+        a = MIN( 1._wp, zbetdo(ji,jj), betup_p1 ) ; \
+        b = MIN( 1._wp, zbetup(ji,jj), betdo_p1 ) ; \
+        c = ( 0.5_wp  + SIGN( 0.5_wp , vec_in(ji,jj,jk) ) ) ; \
+        vec(ji,jj,jk) = vec_in(ji,jj,jk) * ( c * a + ( 1._wp - c) * b )
+
+#define monotonic_flux_k(a,b,c,betup,betdo,vec,vec_in,jk) \
+        a = MIN( 1._wp, betdo, zbetup_ptr(ji,jj) ) ; \
+        b = MIN( 1._wp, betup, zbetdo_ptr(ji,jj) ) ; \
+        c = ( 0.5 + SIGN( 0.5_wp , vec_in(ji,jj,jk) ) ) ; \
+        vec(ji,jj,jk) = vec_in(ji,jj,jk) * ( c * a + ( 1._wp - c) * b )
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: traadv_fct.F90 13660 2020-10-22 10:47:32Z francesca $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   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(inout) ::   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_im1, ztu_ip1, ztv_jm1, ztv_jp1  
+      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwx_3d, zwy_3d, zwz, ztw, zltu_3d, zltv_3d, ztu, ztv
+      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
+            DO jk = 1, jpkm1                 ! Laplacian
+               DO_2D( 1, 0, 1, 0 )                 ! 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_2D
+               DO_2D( 0, 0, 0, 0 )                 ! 2nd derivative * 1/ 6
+                  zltu_3d(ji,jj,jk) = (  ztu(ji,jj,jk) + ztu(ji-1,jj,jk)  ) * r1_6
+                  zltv_3d(ji,jj,jk) = (  ztv(ji,jj,jk) + ztv(ji,jj-1,jk)  ) * r1_6
+               END_2D
+            END DO
+            CALL lbc_lnk_multi( 'traadv_fct', zltu_3d, 'T', 1.0_wp , zltv_3d, 'T', 1.0_wp )   ! 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)
+               !                                                        ! C4 minus upstream advective fluxes 
+               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) ) - 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) ) - zwy_3d(ji,jj,jk)
+            END_3D
+            !
+            CALL lbc_lnk_multi( 'traadv_fct', zwx_3d, 'U', -1.0_wp , zwy_3d, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
+         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_multi( '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  ==!
+         !
+#if defined key_agrif
+         CALL nonosc( Kmm, pt(:,:,:,jn,Kbb), zwx_3d, zwy_3d, zwz, zwi, p2dt )
+#else
+         CALL nonosc_lf( Kmm, pt(:,:,:,jn,Kbb), zwx_3d, zwy_3d, zwz, zwi, p2dt )
+#endif
+         !
+         !        !==  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
+
+   SUBROUTINE nonosc_lf( Kmm, pbef, paa, pbb, pcc, paft, p2dt )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE nonosc  ***
+      !!     
+      !! **  Purpose :   compute monotonic tracer fluxes from the upstream 
+      !!       scheme and the before field by a nonoscillatory algorithm 
+      !!
+      !! **  Method  :   ... ???
+      !!       warning : pbef and paft must be masked, but the boundaries
+      !!       conditions on the fluxes are not necessary zalezak (1979)
+      !!       drange (1995) multi-dimensional forward-in-time and upstream-
+      !!       in-space based differencing for fluid
+      !!----------------------------------------------------------------------
+      INTEGER                          , INTENT(in   ) ::   Kmm             ! time level index 
+      REAL(wp)                         , INTENT(in   ) ::   p2dt            ! tracer time-step
+      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field
+      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) ::   paa, pbb, pcc   ! monotonic fluxes in the 3 directions
+      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   paa_in, pbb_in, pcc_in   ! monotonic fluxes in the 3 directions
+      !
+      REAL(dp), DIMENSION (jpi,jpj,jpk) :: zbup, zbdo
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(dp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn    ! local scalars
+      REAL(dp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
+      REAL(dp) ::   zbt_ip1, zpos_ip1, zneg_ip1, zup_ip1, zdo_ip1, zbetup_ip1, zbetdo_ip1 
+      REAL(dp) ::   zbt_jp1, zpos_jp1, zneg_jp1, zup_jp1, zdo_jp1, zbetup_jp1, zbetdo_jp1
+      REAL(dp), TARGET, DIMENSION(jpi,jpj) :: zbetup_buf, zbetdo_buf, zbetup_ptr_buf, zbetdo_ptr_buf
+      REAL(dp), POINTER, DIMENSION(:,:) :: tmp, zbetup, zbetdo, zbetup_ptr, zbetdo_ptr
+      !!----------------------------------------------------------------------
+      !
+      zbig  = 1.e+40_dp
+      zrtrn = 1.e-15_dp
+
+      paa_in(:,:,:) = paa(:,:,:)
+      pbb_in(:,:,:) = pbb(:,:,:)
+      pcc_in(:,:,:) = pcc(:,:,:)
+      ! Search local extrema
+      ! --------------------
+      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
+      zbup = MAX( pbef * tmask - zbig * ( 1._wp - tmask ),   &
+         &        paft * tmask - zbig * ( 1._wp - tmask )  )
+      zbdo = MIN( pbef * tmask + zbig * ( 1._wp - tmask ),   &
+         &        paft * tmask + zbig * ( 1._wp - tmask )  )
+
+      zbetup => zbetup_buf
+      zbetdo => zbetdo_buf
+      zbetup_ptr => zbetup_ptr_buf
+      zbetdo_ptr => zbetdo_ptr_buf
+
+      DO_2D( 1, 0, 1, 0 )
+
+        ! search maximum in neighbourhood
+        search_in_neighbour(zup,MAX,zbup,ji,jj,1)
+        search_in_neighbour(zup_ip1,MAX,zbup,ji+1,jj,1)
+        search_in_neighbour(zup_jp1,MAX,zbup,ji,jj+1,1)
+
+        ! search minimum in neighbourhood
+        search_in_neighbour(zdo,MIN,zbdo,ji,jj,1)
+        search_in_neighbour(zdo_ip1,MIN,zbdo,ji+1,jj,1)
+        search_in_neighbour(zdo_jp1,MIN,zbdo,ji,jj+1,1)
+
+        ! positive part of the flux
+        pos_part_of_flux(ji,jj,1,zpos)
+        pos_part_of_flux(ji+1,jj,1,zpos_ip1)
+        pos_part_of_flux(ji,jj+1,1,zpos_jp1)
+
+        ! negative part of the flux
+        neg_part_of_flux(ji,jj,1,zneg)
+        neg_part_of_flux(ji+1,jj,1,zneg_ip1)
+        neg_part_of_flux(ji,jj+1,1,zneg_jp1)
+
+        ! up & down beta terms
+        beta_terms(zbt,zbetup(ji,jj),zbetdo(ji,jj),zup,zpos,zdo,zneg,ji,jj,1)
+        beta_terms(zbt_ip1,zbetup_ip1,zbetdo_ip1,zup_ip1,zpos_ip1,zdo_ip1,zneg_ip1,ji+1,jj,1)
+        beta_terms(zbt_jp1,zbetup_jp1,zbetdo_jp1,zup_jp1,zpos_jp1,zdo_jp1,zneg_jp1,ji,jj+1,1)
+
+        ! 3. monotonic flux in the i & j (paa & pbb)
+        ! ----------------------------------------
+        monotonic_flux(zau,zbu,zcu,zbetup_ip1,zbetdo_ip1,paa,paa_in,1)
+        monotonic_flux(zav,zbv,zcv,zbetup_jp1,zbetdo_jp1,pbb,pbb_in,1)
+
+      END_2D
+      tmp => zbetup_ptr
+      zbetup_ptr => zbetup
+      zbetup => tmp
+
+      tmp => zbetdo_ptr
+      zbetdo_ptr => zbetdo
+      zbetdo => tmp
+
+      DO jk = 2, jpk-1
+         DO_2D( 1, 0, 1, 0 )
+
+            ! search maximum in neighbourhood
+            search_in_neighbour(zup,MAX,zbup,ji,jj,jk)
+            search_in_neighbour(zup_ip1,MAX,zbup,ji+1,jj,jk)
+            search_in_neighbour(zup_jp1,MAX,zbup,ji,jj+1,jk)
+
+            ! search minimum in neighbourhood
+            search_in_neighbour(zdo,MIN,zbdo,ji,jj,jk)
+            search_in_neighbour(zdo_ip1,MIN,zbdo,ji+1,jj,jk)
+            search_in_neighbour(zdo_jp1,MIN,zbdo,ji,jj+1,jk)
+
+            ! positive part of the flux
+            pos_part_of_flux(ji,jj,jk,zpos)
+            pos_part_of_flux(ji+1,jj,jk,zpos_ip1)
+            pos_part_of_flux(ji,jj+1,jk,zpos_jp1)
+
+            ! negative part of the flux
+            neg_part_of_flux(ji,jj,jk,zneg)
+            neg_part_of_flux(ji+1,jj,jk,zneg_ip1)
+            neg_part_of_flux(ji,jj+1,jk,zneg_jp1)
+
+            ! up & down beta terms
+            beta_terms(zbt,zbetup(ji,jj),zbetdo(ji,jj),zup,zpos,zdo,zneg,ji,jj,jk)
+            beta_terms(zbt_ip1,zbetup_ip1,zbetdo_ip1,zup_ip1,zpos_ip1,zdo_ip1,zneg_ip1,ji+1,jj,jk)
+            beta_terms(zbt_jp1,zbetup_jp1,zbetdo_jp1,zup_jp1,zpos_jp1,zdo_jp1,zneg_jp1,ji,jj+1,jk)
+
+            ! 3. monotonic flux in the i & j (paa & pbb)
+            ! ----------------------------------------
+            monotonic_flux(zau,zbu,zcu,zbetup_ip1,zbetdo_ip1,paa,paa_in,jk)
+            monotonic_flux(zav,zbv,zcv,zbetup_jp1,zbetdo_jp1,pbb,pbb_in,jk)
+
+            ! monotonic flux in the k direction, i.e. pcc
+            ! -------------------------------------------
+            monotonic_flux_k(za,zb,zc,zbetup(ji,jj),zbetdo(ji,jj),pcc,pcc_in,jk)
+         END_2D
+         tmp => zbetup_ptr
+         zbetup_ptr => zbetup
+         zbetup => tmp
+
+         tmp => zbetdo_ptr
+         zbetdo_ptr => zbetdo
+         zbetdo => tmp
+      END DO
+      !
+      DO_2D( 1, 0, 1, 0 )
+        ! monotonic flux in the k direction, i.e. pcc
+        ! -------------------------------------------
+        monotonic_flux_k(za,zb,zc,0._dp,0._dp,pcc,pcc_in,jpk)
+      END_2D
+   END SUBROUTINE nonosc_lf
+
+END MODULE traadv_fct_lf
Index: NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/TRA/traadv_mus_lf.F90
===================================================================
--- NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/TRA/traadv_mus_lf.F90	(revision 13946)
+++ NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/TRA/traadv_mus_lf.F90	(revision 13946)
@@ -0,0 +1,289 @@
+MODULE traadv_mus_lf
+   !!======================================================================
+   !!                       ***  MODULE  traadv_mus  ***
+   !! Ocean  tracers:  horizontal & vertical advective trend
+   !!======================================================================
+   !! History :       !  2000-06  (A.Estublier)  for passive tracers
+   !!                 !  2001-08  (E.Durand, G.Madec)  adapted for T & S
+   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
+   !!            3.2  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
+   !!            3.4  !  2012-06  (P. Oddo, M. Vichi) include the upstream where needed
+   !!            3.7  !  2015-09  (G. Madec) add the ice-shelf cavities boundary condition
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   tra_adv_mus   : update the tracer trend with the horizontal
+   !!                   and vertical advection trends using MUSCL scheme
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and active tracers
+   USE trc_oce        ! share passive tracers/Ocean variables
+   USE dom_oce        ! ocean space and time domain
+   USE trd_oce        ! trends: ocean variables
+   USE trdtra         ! tracers trends manager
+   USE sbcrnf         ! river runoffs
+   USE diaptr         ! poleward transport diagnostics
+   USE diaar5         ! AR5 diagnostics
+
+   !
+   USE iom            ! XIOS library
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! distribued memory computing
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   tra_adv_mus_lf   ! routine called by traadv.F90
+   
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   upsmsk   !: mixed upstream/centered scheme near some straits
+   !                                                           !  and in closed seas (orca 2 and 1 configurations)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xind     !: mixed upstream/centered index
+   
+   LOGICAL  ::   l_trd   ! flag to compute trends
+   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
+   LOGICAL  ::   l_hst   ! flag to compute heat/salt transport
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+
+#define initial_slop_i(out, ji) \
+        out = umask(ji,jj,jk) * ( pt(ji+1,jj,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
+
+#define initial_slop_j(out, jj) \
+        out = vmask(ji,jj,jk) * ( pt(ji,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
+
+#define initial_slop_k(out, jk) \
+        out = tmask(ji,jj,jk) * ( pt(ji,jj,jk-1,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
+
+#define tracer_slop(out, zzw, zzwm1) \
+        out = ( zzw + zzwm1 ) * ( 0.25 + SIGN( 0.25_wp, zzw * zzwm1 ) )
+
+#define limitation_slop(out, zzslp, zzwm1, zzw) \
+        out = SIGN( 1.0_wp, zzslp ) * MIN( ABS( zzslp ), 2.*ABS( zzwm1 ), 2.*ABS( zzw ) )
+
+#define vertical_adv_flux_i(out, jk, slp, slp1) \
+        z0u = SIGN( 0.5_wp, pU(ji,jj,jk) ) ; \
+        zalpha = 0.5 - z0u ; \
+        zu  = z0u - 0.5 * pU(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) ; \
+        zzwx = pt(ji+1,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * slp1 ; \
+        zzwy = pt(ji  ,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * slp ; \
+        out = pU(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
+
+#define vertical_adv_flux_j(out, jk, slp, slp1) \
+        z0v = SIGN( 0.5_wp, pV(ji,jj,jk) ) ; \
+        zalpha = 0.5 - z0v ; \
+        zv  = z0v - 0.5 * pV(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) ; \
+        zzwx = pt(ji,jj+1,jk,jn,Kbb) + xind(ji,jj,jk) * zv * slp1 ; \
+        zzwy = pt(ji,jj  ,jk,jn,Kbb) + xind(ji,jj,jk) * zv * slp ; \
+        out = pV(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
+
+#define vertical_adv_flux(out, jk, slp, slp1) \
+        z0w = SIGN( 0.5_wp, pW(ji,jj,jk+1) ) ; \
+        zalpha = 0.5 + z0w ; \
+        zw  = z0w - 0.5 * pW(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,Kmm) ; \
+        zzwx = pt(ji,jj,jk+1,jn,Kbb) + xind(ji,jj,jk) * zw * slp1 ; \
+        zzwy = pt(ji,jj,jk  ,jn,Kbb) + xind(ji,jj,jk) * zw * slp ; \
+        out = pW(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk) 
+
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: traadv_mus.F90 13619 2020-10-16 08:41:21Z francesca $ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE tra_adv_mus_lf( kt, kit000, cdtype, p2dt, pU, pV, pW,             &
+      &                    Kbb, Kmm, pt, kjpt, Krhs, ld_msc_ups )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE tra_adv_mus  ***
+      !!
+      !! ** Purpose :   Compute the now trend due to total advection of tracers
+      !!              using a MUSCL scheme (Monotone Upstream-centered Scheme for
+      !!              Conservation Laws) and add it to the general tracer trend.
+      !!
+      !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries
+      !!              ld_msc_ups=T : 
+      !!
+      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
+      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
+      !!             - poleward advective heat and salt transport (ln_diaptr=T)
+      !!
+      !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
+      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
+      !!----------------------------------------------------------------------
+      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
+      LOGICAL                                  , INTENT(in   ) ::   ld_msc_ups      ! use upstream scheme within muscl
+      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
+      INTEGER  ::   ierr             ! local integer
+      REAL(wp) ::   zu, z0u, zw , zalpha   ! local scalars
+      REAL(wp) ::   zv, z0v, z0w           !   -      -
+      REAL(wp) ::   zzwx, zzwxm1, zzwxp1, zzwy, zzwym1, zzwyp1
+      REAL(wp) ::   zzslpx, zzslpxp1, zzslpy, zzslpyp1
+      REAL(wp), TARGET, DIMENSION(jpi,jpj) :: zzwz_buf, zzwzp1_buf, zzwzp2_buf
+      REAL(wp), TARGET, DIMENSION(jpi,jpj) :: zzslpz_buf, zzslpzp1_buf
+      REAL(wp), POINTER, DIMENSION(:,:)    :: tmp, zzwz_ptr, zzwzp1_ptr, zzwzp2_ptr
+      REAL(wp), POINTER, DIMENSION(:,:)    :: zzslpz_ptr, zzslpzp1_ptr 
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz, zwx, zwy   ! 3D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == kit000 )  THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
+         IF(lwp) WRITE(numout,*) '        : mixed up-stream           ', ld_msc_ups
+         IF(lwp) WRITE(numout,*) '~~~~~~~'
+         IF(lwp) WRITE(numout,*)
+         !
+         ! Upstream / MUSCL scheme indicator
+         !
+         ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr )
+         xind(:,:,:) = 1._wp              ! set equal to 1 where up-stream is not needed
+         !
+         IF( ld_msc_ups ) THEN            ! define the upstream indicator (if asked)
+            ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
+            upsmsk(:,:) = 0._wp                             ! not upstream by default
+            !
+            DO jk = 1, jpkm1
+               xind(:,:,jk) = 1._wp                              &                 ! =>1 where up-stream is not needed
+                  &         - MAX ( rnfmsk(:,:) * rnfmsk_z(jk),  &                 ! =>0 near runoff mouths (& closed sea outflows)
+                  &                 upsmsk(:,:)                ) * tmask(:,:,jk)   ! =>0 in some user defined area
+            END DO
+         ENDIF 
+         !
+      ENDIF 
+      !      
+      l_trd = .FALSE.
+      l_hst = .FALSE.
+      l_ptr = .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.
+      !
+         zzwz_ptr => zzwz_buf
+         zzwzp1_ptr => zzwzp1_buf
+         zzwzp2_ptr => zzwzp2_buf
+         zzslpz_ptr => zzslpz_buf
+         zzslpzp1_ptr => zzslpzp1_buf
+         !
+      DO jn = 1, kjpt            !==  loop over the tracers  ==!
+         !
+         zwx(:,:,jpk) = 0._wp                   ! bottom values
+         zwy(:,:,jpk) = 0._wp
+         zwz(:,:, 1 ) = 0._wp                   ! surface & bottom boundary conditions
+         zwz(:,:,jpk) = 0._wp
+         !                          !* Horizontal advective fluxes
+         !
+         !!----------------------------------------------------------------------
+         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+            !-- first guess of the slopes
+            initial_slop_i(zzwxm1, ji-1)
+            initial_slop_i(zzwx, ji)
+            initial_slop_i(zzwxp1, ji+1)
+
+            initial_slop_j(zzwym1, jj-1)
+            initial_slop_j(zzwy, jj)
+            initial_slop_j(zzwyp1, jj+1)
+            !-- Slopes of tracer
+            tracer_slop(zzslpx, zzwx, zzwxm1) 
+            tracer_slop(zzslpxp1, zzwxp1, zzwx) 
+            tracer_slop(zzslpy, zzwy, zzwym1) 
+            tracer_slop(zzslpyp1, zzwyp1, zzwy) 
+            !-- Slopes limitation
+            limitation_slop(zzslpx, zzslpx, zzwxm1, zzwx)
+            limitation_slop(zzslpxp1, zzslpxp1, zzwx, zzwxp1)
+            limitation_slop(zzslpy, zzslpy, zzwym1, zzwy)
+            limitation_slop(zzslpyp1, zzslpyp1, zzwy, zzwyp1)
+            !-- MUSCL horizontal advective fluxes
+            vertical_adv_flux_i(zwx(ji,jj,jk), jk, zzslpx, zzslpxp1) 
+            vertical_adv_flux_j(zwy(ji,jj,jk), jk, zzslpy, zzslpyp1) 
+         END_3D
+         !
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )    !-- Tracer advective trend
+            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk)       &
+            &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     &
+            &                                     * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+         END_3D
+         !                          !* Vertical advective fluxes
+         DO_2D( 0, 0, 0, 0 ) 
+            !-- first guess of the slopes
+            initial_slop_k(zzwz_ptr(ji,jj), 2)
+            initial_slop_k(zzwzp1_ptr(ji,jj), 3)
+            !-- Slopes of tracer
+            tracer_slop(zzslpz_ptr(ji,jj), zzwz_ptr(ji,jj), zzwzp1_ptr(ji,jj))
+            !-- Slopes limitation
+            limitation_slop(zzslpz_ptr(ji,jj), zzslpz_ptr(ji,jj), zzwzp1_ptr(ji,jj), zzwz_ptr(ji,jj))
+            !-- vertical advective flux
+            vertical_adv_flux(zwz(ji,jj,2), 1, 0, zzslpz_ptr(ji,jj))
+         END_2D
+         
+         DO jk = 2, jpk-3    
+            DO_2D( 0, 0, 0, 0 )
+               !-- first guess of the slopes
+               initial_slop_k(zzwzp2_ptr(ji,jj), jk+2)
+               !-- Slopes of tracer
+               tracer_slop(zzslpzp1_ptr(ji,jj), zzwzp1_ptr(ji,jj), zzwzp2_ptr(ji,jj))
+               !-- Slopes limitation
+               limitation_slop(zzslpzp1_ptr(ji,jj), zzslpzp1_ptr(ji,jj), zzwzp2_ptr(ji,jj), zzwzp1_ptr(ji,jj))
+               !-- vertical advective flux
+               vertical_adv_flux(zwz(ji,jj,jk+1), jk, zzslpz_ptr(ji,jj), zzslpzp1_ptr(ji,jj))
+            END_2D
+            tmp => zzwzp1_ptr
+            zzwzp1_ptr => zzwzp2_ptr
+            zzwzp2_ptr => tmp
+
+            tmp => zzslpz_ptr
+            zzslpz_ptr => zzslpzp1_ptr
+            zzslpzp1_ptr => tmp 
+         END DO
+         DO_2D( 0, 0, 0, 0 ) 
+            !-- Slopes of tracer
+            tracer_slop(zzslpzp1_ptr(ji,jj), zzwzp1_ptr(ji,jj), 0)
+            !-- Slopes limitation
+            limitation_slop(zzslpzp1_ptr(ji,jj), zzslpzp1_ptr(ji,jj), 0, zzwzp1_ptr(ji,jj))
+            !-- vertical advective flux
+            vertical_adv_flux(zwz(ji,jj,jpk-1), jpk-2, zzslpz_ptr(ji,jj), zzslpzp1_ptr(ji,jj))
+         END_2D
+
+         IF( ln_linssh ) THEN                   ! top values, linear free surface only
+            IF( ln_isfcav ) THEN                      ! ice-shelf cavities (top of the ocean)
+               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)
+               END_2D
+            ELSE                                      ! no cavities: only at the ocean surface
+               zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
+            ENDIF
+         ENDIF
+         !
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )     !-- vertical advective trend
+            pt(ji,jj,jk,jn,Krhs) =  pt(ji,jj,jk,jn,Krhs) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
+               &                                      * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+         END_3D
+         !                                ! trend horizontal diagnostics
+         IF( l_trd )  THEN
+            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kbb) )
+            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kbb) )
+         END IF
+         !                                 ! "Poleward" heat and salt transports 
+         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
+         !                                 !  heat transport
+         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
+         !
+         !                                ! send vertical trends for diagnostic
+         IF( l_trd )  CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kbb) )
+         !
+      END DO                     ! end of tracer loop
+      !
+   END SUBROUTINE tra_adv_mus_lf
+
+   !!======================================================================
+END MODULE traadv_mus_lf
Index: NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/TRA/traldf.F90
===================================================================
--- NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/TRA/traldf.F90	(revision 13925)
+++ NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/TRA/traldf.F90	(revision 13946)
@@ -93,5 +93,9 @@
          CASE ( np_blp , np_blp_i , np_blp_it )             ! bilaplacian: iso-level & iso-neutral operators
             ! NOTE: [tiling-comms-merge] This lbc_lnk is still needed in the zco case, because zps_hde is not called in step
+#if defined key_mpi3
+            IF(nn_hls.EQ.2) CALL lbc_lnk_nc_multi( 'tra_ldf', pts(:,:,:,:,Kbb), 'T',1.)
+#else
             IF(nn_hls.EQ.2) CALL lbc_lnk( 'tra_ldf', pts(:,:,:,:,Kbb), 'T',1.)
+#endif
             CALL tra_ldf_blp  ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs),             jpts, nldf_tra )
          END SELECT
Index: NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/TRA/zpshde.F90
===================================================================
--- NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/TRA/zpshde.F90	(revision 13925)
+++ NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/OCE/TRA/zpshde.F90	(revision 13946)
@@ -125,6 +125,16 @@
       IF( ln_timing )   CALL timing_start( 'zps_hde')
       ! NOTE: [tiling-comms-merge] Some lbc_lnks in tra_adv and tra_ldf can be taken out in the zps case, because this lbc_lnk is called when zps_hde is called in the stp routine. In the zco case they are still needed.
-      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
+      IF (nn_hls.EQ.2) THEN
+#if defined key_mpi3
+         CALL lbc_lnk_nc_multi( 'zpshde', pta, 'T', 1.0_wp)
+#else
+         CALL lbc_lnk( 'zpshde', pta, 'T', 1.0_wp)
+#endif
+#if defined key_mpi3
+         IF(PRESENT(prd)) CALL lbc_lnk_nc_multi( 'zpshde', prd, 'T', 1.0_wp)
+#else
+         IF(PRESENT(prd)) CALL lbc_lnk( 'zpshde', prd, 'T', 1.0_wp)
+#endif
+      END IF
       !
       pgtu(:,:,:) = 0._wp   ;   zti (:,:,:) = 0._wp   ;   zhi (:,:) = 0._wp
@@ -314,6 +324,16 @@
       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
+      IF (nn_hls.EQ.2) THEN
+#if defined key_mpi3
+         CALL lbc_lnk_nc_multi( 'zpshde', pta, 'T', 1.0_wp)
+#else
+         CALL lbc_lnk( 'zpshde', pta, 'T', 1.0_wp)
+#endif
+#if defined key_mpi3
+         IF (PRESENT(prd)) CALL lbc_lnk_nc_multi( 'zpshde', prd, 'T', 1.0_wp)
+#else
+         IF (PRESENT(prd)) CALL lbc_lnk( 'zpshde', prd, 'T', 1.0_wp)
+#endif
+      END IF
 
       pgtu (:,:,:) = 0._wp   ;   pgtv (:,:,:) =0._wp
Index: NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/TOP/TRP/trcadv.F90
===================================================================
--- NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/TOP/TRP/trcadv.F90	(revision 13925)
+++ NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/TOP/TRP/trcadv.F90	(revision 13946)
@@ -22,5 +22,7 @@
    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)
@@ -127,5 +129,5 @@
          IF (nn_hls.EQ.2) CALL lbc_lnk_nc_multi( 'trcadv', ptr(:,:,:,:,Kmm), 'T', 1.)
 #else
-         IF (nn_hls.EQ.2) CALL lbc_lnk_multi( 'trcadv', ptr(:,:,:,:,Kmm), 'T', 1.)
+         IF (nn_hls.EQ.2) CALL lbc_lnk( 'trcadv', ptr(:,:,:,:,Kmm), 'T', 1.)
 #endif
          CALL tra_adv_cen( kt, nittrc000,'TRC',          zuu, zvv, zww,      Kmm, ptr, jptra, Krhs, nn_cen_h, nn_cen_v )
@@ -139,18 +141,32 @@
             CALL lbc_lnk_multi( 'traadv', zuu(:,:,:), 'U', -1., zvv(:,:,:), 'V', -1., zww(:,:,:), 'W', 1.)
 #endif
+#if defined key_loop_fusion
+            CALL tra_adv_fct_lf( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, nn_fct_h, nn_fct_v )
+#else
+            CALL tra_adv_fct( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, nn_fct_h, nn_fct_v )
+#endif
+         ELSE
+            CALL tra_adv_fct( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, nn_fct_h, nn_fct_v )
          END IF
-         CALL tra_adv_fct( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, nn_fct_h, nn_fct_v )
       CASE ( np_MUS )                                 ! MUSCL
+         IF (nn_hls.EQ.2) THEN
 #if defined key_mpi3
             IF (nn_hls.EQ.2) CALL lbc_lnk_nc_multi( 'trcadv', ptr(:,:,:,:,Kbb), 'T', 1.)
 #else
-            IF (nn_hls.EQ.2) CALL lbc_lnk_multi( 'trcadv', ptr(:,:,:,:,Kbb), 'T', 1.)
-#endif
-         CALL tra_adv_mus( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, ln_mus_ups         ) 
+            IF (nn_hls.EQ.2) CALL lbc_lnk( 'trcadv', ptr(:,:,:,:,Kbb), 'T', 1.)
+#endif
+#if defined key_loop_fusion
+            CALL tra_adv_mus_lf( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, ln_mus_ups ) 
+#else
+            CALL tra_adv_mus( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, ln_mus_ups ) 
+#endif
+         ELSE
+            CALL tra_adv_mus( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, ln_mus_ups ) 
+         END IF
       CASE ( np_UBS )                                 ! UBS
 #if defined key_mpi3
          IF (nn_hls.EQ.2) CALL lbc_lnk_nc_multi( 'trcadv', ptr(:,:,:,:,Kbb), 'T', 1.)
 #else
-         IF (nn_hls.EQ.2) CALL lbc_lnk_multi( 'trcadv', ptr(:,:,:,:,Kbb), 'T', 1.)
+         IF (nn_hls.EQ.2) CALL lbc_lnk( 'trcadv', ptr(:,:,:,:,Kbb), 'T', 1.)
 #endif
          CALL tra_adv_ubs( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, nn_ubs_v           )
Index: NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/TOP/TRP/trcldf.F90
===================================================================
--- NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/TOP/TRP/trcldf.F90	(revision 13925)
+++ NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/TOP/TRP/trcldf.F90	(revision 13946)
@@ -101,5 +101,9 @@
            &                     ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Krhs), jptra, 1 )
       CASE ( np_blp , np_blp_i , np_blp_it )                                                               ! bilaplacian: all operator (iso-level, -neutral)
+#if defined key_mpi3
+         IF(nn_hls.EQ.2) CALL lbc_lnk_nc_multi( 'trc_ldf', ptr(:,:,:,:,Kbb), 'T',1.)
+#else
          IF(nn_hls.EQ.2) CALL lbc_lnk( 'trc_ldf', ptr(:,:,:,:,Kbb), 'T',1.)
+#endif
          CALL tra_ldf_blp  ( kt, Kmm, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi,            &
            &                     ptr(:,:,:,:,Kbb) , ptr(:,:,:,:,Krhs),                 jptra, nldf_trc )
