Index: /NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/OCE/TRA/traadv_fct_lf.F90
===================================================================
--- /NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/OCE/TRA/traadv_fct_lf.F90	(revision 13881)
+++ /NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/OCE/TRA/traadv_fct_lf.F90	(revision 13882)
@@ -47,33 +47,44 @@
 #  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))
+        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)) 
+        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))
+        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
+        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, betup_p1 ) ; \
-b = MIN( 1._wp, zbetup, 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_p1,betdo_p1,vec,vec_in,jk) \
-a = MIN( 1._wp, betdo_p1, zbetup ) ; \
-b = MIN( 1._wp, betup_p1, zbetdo ) ; \
-c = ( 0.5 + SIGN( 0.5_wp , vec_in(ji,jj,jk+1) ) ) ; \
-vec(ji,jj,jk+1) = vec_in(ji,jj,jk+1) * ( c * a + ( 1._wp - c) * b )
+        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 )
 
    !!----------------------------------------------------------------------
@@ -114,7 +125,7 @@
       INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices  
       REAL(wp) ::   ztra                                     ! local scalar
-      REAL(wp) ::   zwx, zwx_im1, zfp_ui, zfp_ui_m1, zfp_vj, zfp_vj_m1, zfp_wk, zC2t_u, zC4t_u   !   -      -
-      REAL(wp) ::   zwy, 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, zltu_ip1, zltv_jp1, zltu, zltv
+      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
@@ -129,6 +140,4 @@
       ENDIF
       !! -- init to 0
-      !! -- init to 0
-      zwi(:,:,:) = 0._wp
       zwx_3d(:,:,:) = 0._wp
       zwy_3d(:,:,:) = 0._wp
@@ -173,5 +182,4 @@
          !
          !        !==  upstream advection with initial mass fluxes & intermediate update  ==!
-         !                    !* upstream tracer flux in the i and j direction 
          !                               !* upstream tracer flux in the k direction *!
          DO_3D( 1, 1, 1, 1, 2, jpkm1 )      ! Interior value ( multiplied by wmask)
@@ -192,47 +200,20 @@
          ENDIF
          !               
+         !                    !* upstream tracer flux in the i and j direction 
          DO jk = 1, jpkm1
             DO jj = 1, jpj-1
-               zfp_ui = pU(1,jj,jk) + ABS( pU(1,jj,jk) )
-               zfm_ui = pU(1,jj,jk) - ABS( pU(1,jj,jk) )
-               zwx = 0.5 * ( zfp_ui * pt(1,jj,jk,jn,Kbb) + zfm_ui * pt(2,jj  ,jk,jn,Kbb) )
-               zwx_3d(1,jj,jk) = zwx
-
-               zfp_vj = pV(1,jj,jk) + ABS( pV(1,jj,jk) )
-               zfm_vj = pV(1,jj,jk) - ABS( pV(1,jj,jk) )
-               zwy = 0.5 * ( zfp_vj * pt(1,jj,jk,jn,Kbb) + zfm_vj * pt(1  ,jj+1,jk,jn,Kbb) )
-               zwy_3d(1,jj,jk) = zwy
+               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
-               zfp_ui = pU(ji,1,jk) + ABS( pU(ji,1,jk) )
-               zfm_ui = pU(ji,1,jk) - ABS( pU(ji,1,jk) )
-               zwx = 0.5 * ( zfp_ui * pt(ji,1,jk,jn,Kbb) + zfm_ui * pt(ji+1,1  ,jk,jn,Kbb) )
-               zwx_3d(ji,1,jk) = zwx
-
-               zfp_vj = pV(ji,1,jk) + ABS( pV(ji,1,jk) )
-               zfm_vj = pV(ji,1,jk) - ABS( pV(ji,1,jk) )
-               zwy = 0.5 * ( zfp_vj * pt(ji,1,jk,jn,Kbb) + zfm_vj * pt(ji  ,2,jk,jn,Kbb) )
-               zwy_3d(ji,1,jk) = zwy
+               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 )
-               zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )
-               zfm_ui = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) )
-               zwx = 0.5 * ( zfp_ui * pt(ji,jj,jk,jn,Kbb) + zfm_ui * pt(ji+1,jj  ,jk,jn,Kbb) )
-               zwx_3d(ji,jj,jk) = zwx
-
-               zfp_ui_m1 = pU(ji-1,jj,jk) + ABS( pU(ji-1,jj,jk) )
-               zfm_ui_m1 = pU(ji-1,jj,jk) - ABS( pU(ji-1,jj,jk) )
-               zwx_im1 = 0.5 * ( zfp_ui_m1 * pt(ji-1,jj,jk,jn,Kbb) + zfm_ui_m1 * pt(ji,jj  ,jk,jn,Kbb) )
-
-               zfp_vj = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) )
-               zfm_vj = pV(ji,jj,jk) - ABS( pV(ji,jj,jk) )
-               zwy = 0.5 * ( zfp_vj * pt(ji,jj,jk,jn,Kbb) + zfm_vj * pt(ji  ,jj+1,jk,jn,Kbb) )
-               zwy_3d(ji,jj,jk) = zwy
-
-               zfp_vj_m1 = pV(ji,jj-1,jk) + ABS( pV(ji,jj-1,jk) )
-               zfm_vj_m1 = pV(ji,jj-1,jk) - ABS( pV(ji,jj-1,jk) )
-               zwy_jm1 = 0.5 * ( zfp_vj_m1 * pt(ji,jj-1,jk,jn,Kbb) + zfm_vj_m1 * pt(ji,jj,jk,jn,Kbb) )
-               !                               ! total intermediate advective trends
-               ztra = - (  zwx - zwx_im1 + zwy - zwy_jm1 + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
+               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   &
@@ -449,11 +430,10 @@
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      INTEGER  ::   ikm1         ! local integer
       REAL(dp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn    ! local scalars
       REAL(dp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
-      REAL(dp) ::   zbetup, zbetdo
       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) ::   zbt_kp1, zpos_kp1, zneg_kp1, zup_kp1, zdo_kp1, zbetup_kp1, zbetdo_kp1
+      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
       !!----------------------------------------------------------------------
       !
@@ -472,84 +452,100 @@
          &        paft * tmask + zbig * ( 1._wp - tmask )  )
 
-      DO_3D( 1, 0, 1, 0, 1, jpk-2 )
+      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,jk)
-        search_in_neighbour(zup_ip1,MAX,zbup,ji+1,jj,jk)
-        search_in_neighbour(zup_jp1,MAX,zbup,ji,jj+1,jk)
-        search_in_neighbour(zup_kp1,MAX,zbup,ji,jj,jk+1)
+        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,jk)
-        search_in_neighbour(zdo_ip1,MIN,zbdo,ji+1,jj,jk)
-        search_in_neighbour(zdo_jp1,MIN,zbdo,ji,jj+1,jk)
-        search_in_neighbour(zdo_kp1,MIN,zbdo,ji,jj,jk+1)
+        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,jk,zpos)
-        pos_part_of_flux(ji+1,jj,jk,zpos_ip1)
-        pos_part_of_flux(ji,jj+1,jk,zpos_jp1)
-        pos_part_of_flux(ji,jj,jk+1,zpos_kp1)
+        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,jk,zneg)
-        neg_part_of_flux(ji+1,jj,jk,zneg_ip1)
-        neg_part_of_flux(ji,jj+1,jk,zneg_jp1)
-        neg_part_of_flux(ji,jj,jk+1,zneg_kp1)
+        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,zbetdo,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)
-        beta_terms(zbt_kp1,zbetup_kp1,zbetdo_kp1,zup_kp1,zpos_kp1,zdo_kp1,zneg_kp1,ji,jj,jk+1)
+        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,jk)
-        monotonic_flux(zav,zbv,zcv,zbetup_jp1,zbetdo_jp1,pbb,pbb_in,jk)
-
+        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,zbetup_kp1,zbetdo_kp1,pcc,pcc_in,jk)
-      END_3D
-      !
-      DO_2D( 1, 0, 1, 0 )
-
-        ! search maximum in neighbourhood
-        search_in_neighbour(zup,MAX,zbup,ji,jj,jpk-1)
-        search_in_neighbour(zup_ip1,MAX,zbup,ji+1,jj,jpk-1)
-        search_in_neighbour(zup_jp1,MAX,zbup,ji,jj+1,jpk-1)
-
-        ! search minimum in neighbourhood
-        search_in_neighbour(zdo,MIN,zbdo,ji,jj,jk)
-        search_in_neighbour(zdo_ip1,MIN,zbdo,ji+1,jj,jpk-1)
-        search_in_neighbour(zdo_jp1,MIN,zbdo,ji,jj+1,jpk-1)
-
-        ! positive part of the flux
-        pos_part_of_flux(ji,jj,jpk-1,zpos)
-        pos_part_of_flux(ji+1,jj,jpk-1,zpos_ip1)
-        pos_part_of_flux(ji,jj+1,jpk-1,zpos_jp1)
-
-        ! negative part of the flux
-        neg_part_of_flux(ji,jj,jpk-1,zneg)
-        neg_part_of_flux(ji+1,jj,jpk-1,zneg_ip1)
-        neg_part_of_flux(ji,jj+1,jpk-1,zneg_jp1)
-
-        ! up & down beta terms
-        beta_terms(zbt,zbetup,zbetdo,zup,zpos,zdo,zneg,ji,jj,jpk-1)
-        beta_terms(zbt_ip1,zbetup_ip1,zbetdo_ip1,zup_ip1,zpos_ip1,zdo_ip1,zneg_ip1,ji+1,jj,jpk-1)
-        beta_terms(zbt_jp1,zbetup_jp1,zbetdo_jp1,zup_jp1,zpos_jp1,zdo_jp1,zneg_jp1,ji,jj+1,jpk-1)
-
-        zbetup_kp1 = 0._dp 
-        zbetdo_kp1 = 0._dp
-
-        ! 3. monotonic flux in the i & j direction (paa & pbb)
-        ! ----------------------------------------
-        monotonic_flux(zau,zbu,zcu,zbetup_ip1,zbetdo_ip1,paa,paa_in,jpk-1)
-        monotonic_flux(zav,zbv,zcv,zbetup_jp1,zbetdo_jp1,pbb,pbb_in,jpk-1)
-
-        ! monotonic flux in the k direction, i.e. pcc
-        ! -------------------------------------------
-        monotonic_flux_k(za,zb,zc,zbetup_kp1,zbetdo_kp1,pcc,pcc_in,jpk-1)
+        monotonic_flux_k(za,zb,zc,0._dp,0._dp,pcc,pcc_in,jpk)
       END_2D
    END SUBROUTINE nonosc_lf
Index: /NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/OCE/TRA/traadv_mus_lf.F90
===================================================================
--- /NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/OCE/TRA/traadv_mus_lf.F90	(revision 13881)
+++ /NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/OCE/TRA/traadv_mus_lf.F90	(revision 13882)
@@ -129,7 +129,10 @@
       REAL(wp) ::   zv, z0v, z0w           !   -      -
       REAL(wp) ::   zzwx, zzwxm1, zzwxp1, zzwy, zzwym1, zzwyp1
-      REAL(wp) ::   zzwz, zzwzp1, zzwzp2, zzslpz, zzslpzp1
-      REAL(wp) ::   zzslpx, zzslpx_ip1, zzslpy, zzslpy_jp1
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zwy, zwz 
+      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
       !!----------------------------------------------------------------------
       !
@@ -167,36 +170,41 @@
          &                          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  ==!
          !
-         !                          !* Horizontal advective fluxes
-         !
-         !!----------------------------------------------------------------------
          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(zzwxm1, ji-1)
             initial_slop_i(zzwxp1, ji+1)
 
+            initial_slop_j(zzwym1, jj-1)
             initial_slop_j(zzwy, jj)
-            initial_slop_j(zzwym1, jj-1)
             initial_slop_j(zzwyp1, jj+1)
             !-- Slopes of tracer
             tracer_slop(zzslpx, zzwx, zzwxm1) 
-            tracer_slop(zzslpx_ip1, zzwxp1, zzwx) 
+            tracer_slop(zzslpxp1, zzwxp1, zzwx) 
             tracer_slop(zzslpy, zzwy, zzwym1) 
-            tracer_slop(zzslpy_jp1, zzwyp1, zzwy) 
+            tracer_slop(zzslpyp1, zzwyp1, zzwy) 
             !-- Slopes limitation
             limitation_slop(zzslpx, zzslpx, zzwxm1, zzwx)
-            limitation_slop(zzslpx_ip1, zzslpx_ip1, zzwx, zzwxp1)
+            limitation_slop(zzslpxp1, zzslpxp1, zzwx, zzwxp1)
             limitation_slop(zzslpy, zzslpy, zzwym1, zzwy)
-            limitation_slop(zzslpy_jp1, zzslpy_jp1, zzwy, zzwyp1)
+            limitation_slop(zzslpyp1, zzslpyp1, zzwy, zzwyp1)
             !-- MUSCL horizontal advective fluxes
-            vertical_adv_flux_i(zwx(ji,jj,jk), jk, zzslpx, zzslpx_ip1) 
-            vertical_adv_flux_j(zwy(ji,jj,jk), jk, zzslpy, zzslpy_jp1) 
+            vertical_adv_flux_i(zwx(ji,jj,jk), jk, zzslpx, zzslpxp1) 
+            vertical_adv_flux_j(zwy(ji,jj,jk), jk, zzslpy, zzslpyp1) 
          END_3D
          !
@@ -207,43 +215,42 @@
          END_3D
          !                          !* Vertical advective fluxes
-         !
          DO_2D( 0, 0, 0, 0 ) 
             !-- first guess of the slopes
-            initial_slop_k(zzwzp1, 2)
-            initial_slop_k(zzwzp2, 3)
+            initial_slop_k(zzwz_ptr(ji,jj), 2)
+            initial_slop_k(zzwzp1_ptr(ji,jj), 3)
             !-- Slopes of tracer
-            tracer_slop(zzslpzp1, zzwzp1, zzwzp2)
+            tracer_slop(zzslpz_ptr(ji,jj), zzwz_ptr(ji,jj), zzwzp1_ptr(ji,jj))
             !-- Slopes limitation
-            limitation_slop(zzslpzp1, zzslpzp1, zzwzp2, zzwzp1)
+            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, zzslpzp1)
+            vertical_adv_flux(zwz(ji,jj,2), 1, 0, zzslpz_ptr(ji,jj))
          END_2D
-         DO_3D( 0, 0, 0, 0, 2, jpk-3 )    
-            !-- first guess of the slopes
-            initial_slop_k(zzwz, jk)
-            initial_slop_k(zzwzp1, jk+1)
-            initial_slop_k(zzwzp2, jk+2)
+         
+         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(zzslpz, zzwz, zzwzp1)
-            tracer_slop(zzslpzp1, zzwzp1, zzwzp2)
+            tracer_slop(zzslpzp1_ptr(ji,jj), zzwzp1_ptr(ji,jj), 0)
             !-- Slopes limitation
-            limitation_slop(zzslpz, zzslpz, zzwzp1, zzwz)
-            limitation_slop(zzslpzp1, zzslpzp1, zzwzp2, zzwzp1)
+            limitation_slop(zzslpzp1_ptr(ji,jj), zzslpzp1_ptr(ji,jj), 0, zzwzp1_ptr(ji,jj))
             !-- vertical advective flux
-            vertical_adv_flux(zwz(ji,jj,jk+1), jk, zzslpz, zzslpzp1)
-         END_3D
-         DO_2D( 0, 0, 0, 0 ) 
-            !-- first guess of the slopes
-            initial_slop_k(zzwz, jpk-2)
-            initial_slop_k(zzwzp1, jpk-1)
-            zzwzp2 = 0
-            !-- Slopes of tracer
-            tracer_slop(zzslpz, zzwz, zzwzp1)
-            tracer_slop(zzslpzp1, zzwzp1, zzwzp2)
-            !-- Slopes limitation
-            limitation_slop(zzslpz, zzslpz, zzwzp1, zzwz)
-            limitation_slop(zzslpzp1, zzslpzp1, zzwzp2, zzwzp1)
-            !-- vertical advective flux
-            vertical_adv_flux(zwz(ji,jj,jpk-1), jpk-2, zzslpz, zzslpzp1)
+            vertical_adv_flux(zwz(ji,jj,jpk-1), jpk-2, zzslpz_ptr(ji,jj), zzslpzp1_ptr(ji,jj))
          END_2D
 
