Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynadv.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynadv.F90	(revision 14353)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynadv.F90	(revision 14381)
@@ -51,5 +51,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs )
+   SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_adv  ***
@@ -64,7 +64,8 @@
       !!      (see dynvor module).
       !!----------------------------------------------------------------------
-      INTEGER                             , INTENT( in )  ::  kt               ! ocean time-step index
-      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs   ! ocean time level indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv         ! ocean velocities and RHS of momentum equation
+      INTEGER                                       , INTENT(in   ) ::  kt               ! ocean time-step index
+      INTEGER                                       , INTENT(in   ) ::  Kbb, Kmm, Krhs   ! ocean time level indices
+      REAL(wp), DIMENSION(:,:,:)          , OPTIONAL, TARGET, INTENT(in   ) ::  pau, pav, paw    ! advective velocity
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt)          , TARGET, INTENT(inout) ::  puu, pvv         ! ocean velocities and RHS of momentum Eq.
       !!----------------------------------------------------------------------
       !
@@ -73,10 +74,10 @@
       SELECT CASE( n_dynadv )    !==  compute advection trend and add it to general trend  ==!
       CASE( np_VEC_c2  )     
-         CALL dyn_keg     ( kt, nn_dynkeg,      Kmm, puu, pvv, Krhs )    ! vector form : horizontal gradient of kinetic energy
-         CALL dyn_zad     ( kt,                 Kmm, puu, pvv, Krhs )    ! vector form : vertical advection
+         CALL dyn_keg     ( kt, nn_dynkeg     , Kmm, puu, pvv, Krhs )                  ! vector form : horizontal gradient of kinetic energy
+         CALL dyn_zad     ( kt                , Kmm, puu, pvv, Krhs )                  ! vector form : vertical advection
       CASE( np_FLX_c2  ) 
-         CALL dyn_adv_cen2( kt,                 Kmm, puu, pvv, Krhs )    ! 2nd order centered scheme
+         CALL dyn_adv_cen2( kt                , Kmm, puu, pvv, Krhs, pau, pav, paw )   ! 2nd order centered scheme
       CASE( np_FLX_ubs )   
-         CALL dyn_adv_ubs ( kt,            Kbb, Kmm, puu, pvv, Krhs )    ! 3rd order UBS      scheme (UP3)
+         CALL dyn_adv_ubs ( kt           , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )   ! 3rd order UBS      scheme (UP3)
       END SELECT
       !
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynadv_cen2.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynadv_cen2.F90	(revision 14353)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynadv_cen2.F90	(revision 14381)
@@ -36,5 +36,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs )
+   SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs, pau, pav, paw )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_adv_cen2  ***
@@ -47,11 +47,13 @@
       !! ** Action  :   (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the now vorticity term trend
       !!----------------------------------------------------------------------
-      INTEGER                             , INTENT( in )  ::  kt           ! ocean time-step index
-      INTEGER                             , INTENT( in )  ::  Kmm, Krhs    ! ocean time level indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv     ! ocean velocities and RHS of momentum equation
+      INTEGER                                               , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                                               , INTENT(in   ) ::   Kmm, Krhs       ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt)          , TARGET, INTENT(inout) ::   puu, pvv        ! ocean velocities and RHS of momentum equation
+      REAL(wp), DIMENSION(:,:,:)          , OPTIONAL, TARGET, INTENT(in   ) ::   pau, pav, paw   ! advective velocity
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zfu_t, zfu_f, zfu_uw, zfu
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zfv_t, zfv_f, zfv_vw, zfv, zfw
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zfu_t, zfu_f, zfu_uw, zfu
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zfv_t, zfv_f, zfv_vw, zfv, zfw
+      REAL(wp), DIMENSION(:,:,:), POINTER ::   zptu, zptv, zptw
       !!----------------------------------------------------------------------
       !
@@ -67,9 +69,19 @@
       ENDIF
       !
+      IF( PRESENT( pau ) ) THEN     ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww)
+         zptu => pau(:,:,:)
+         zptv => pav(:,:,:)
+         zptw => paw(:,:,:)
+      ELSE                          ! MLF: advective velocity = (puu,pvv,ww)
+         zptu => puu(:,:,:,Kmm)
+         zptv => pvv(:,:,:,Kmm)
+         zptw => ww (:,:,:    )
+      ENDIF
+      !
       !                             !==  Horizontal advection  ==!
       !
       DO jk = 1, jpkm1                    ! horizontal transport
-         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm)
-         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
+         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * zptu(:,:,jk)
+         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * zptv(:,:,jk)
          DO_2D( 1, 0, 1, 0 )              ! horizontal momentum fluxes (at T- and F-point)
             zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj  ,jk,Kmm) )
@@ -104,11 +116,11 @@
       IF( ln_linssh ) THEN                ! linear free surface: advection through the surface
          DO_2D( 0, 0, 0, 0 )
-            zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
-            zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji,jj+1) * ww(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
+            zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zptw(ji,jj,1) + e1e2t(ji+1,jj) * zptw(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
+            zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zptw(ji,jj,1) + e1e2t(ji,jj+1) * zptw(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
          END_2D
       ENDIF
       DO jk = 2, jpkm1                    ! interior advective fluxes
          DO_2D( 0, 1, 0, 1 )                  ! 1/4 * Vertical transport
-            zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
+            zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * zptw(ji,jj,jk)
          END_2D
          DO_2D( 0, 0, 0, 0 )
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynadv_ubs.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynadv_ubs.F90	(revision 14353)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynadv_ubs.F90	(revision 14381)
@@ -42,5 +42,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs )
+   SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_adv_ubs  ***
@@ -69,7 +69,8 @@
       !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling. 
       !!----------------------------------------------------------------------
-      INTEGER                             , INTENT( in )  ::  kt              ! ocean time-step index
-      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs  ! ocean time level indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv        ! ocean velocities and RHS of momentum equation
+      INTEGER                                               , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                                               , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt)          , TARGET, INTENT(inout) ::   puu, pvv        ! ocean velocities and RHS of momentum equation
+      REAL(wp), DIMENSION(:,:,:)          , OPTIONAL, TARGET, INTENT(in   ) ::   pau, pav, paw   ! advective velocity
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -79,4 +80,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlu_uu, zlu_uv
       REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlv_vv, zlv_vu
+      REAL(wp), DIMENSION(:,:,:), POINTER ::   zptu, zptv, zptw
       !!----------------------------------------------------------------------
       !
@@ -101,10 +103,21 @@
          zfv_vw(:,:,:) = pvv(:,:,:,Krhs)
       ENDIF
+      !
+      IF( PRESENT( pau ) ) THEN     ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww)
+         zptu => pau(:,:,:)
+         zptv => pav(:,:,:)
+         zptw => paw(:,:,:)
+      ELSE                          ! MLF: advective velocity = (puu,pvv,ww)
+         zptu => puu(:,:,:,Kmm)
+         zptv => pvv(:,:,:,Kmm)
+         zptw => ww (:,:,:    )
+      ENDIF
+      !
       !                                      ! =========================== !
       DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
          !                                   ! =========================== !
          !                                         ! horizontal volume fluxes
-         zfu(:,:,jk) = e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm)
-         zfv(:,:,jk) = e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
+         zfu(:,:,jk) = e2u(:,:) * e3u(:,:,jk,Kmm) * zptu(:,:,jk)
+         zfv(:,:,jk) = e1v(:,:) * e3v(:,:,jk,Kmm) * zptv(:,:,jk)
          !            
          DO_2D( 0, 0, 0, 0 )                       ! laplacian
@@ -133,6 +146,6 @@
       DO jk = 1, jpkm1                       ! ====================== !
          !                                         ! horizontal volume fluxes
-         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm)
-         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
+         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * zptu(:,:,jk)
+         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * zptv(:,:,jk)
          !
          DO_2D( 1, 0, 1, 0 )                       ! horizontal momentum fluxes at T- and F-point
@@ -195,11 +208,11 @@
       IF( ln_linssh ) THEN                         ! constant volume : advection through the surface
          DO_2D( 0, 0, 0, 0 )
-            zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
-            zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji,jj+1) * ww(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
+            zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zptw(ji,jj,1) + e1e2t(ji+1,jj) * zptw(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
+            zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zptw(ji,jj,1) + e1e2t(ji,jj+1) * zptw(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
          END_2D
       ENDIF
       DO jk = 2, jpkm1                          ! interior fluxes
          DO_2D( 0, 1, 0, 1 )
-            zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
+            zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * zptw(ji,jj,jk)
          END_2D
          DO_2D( 0, 0, 0, 0 )
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/traadv.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/traadv.F90	(revision 14353)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/traadv.F90	(revision 14381)
@@ -80,5 +80,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv( kt, Kbb, Kmm, pts, Krhs )
+   SUBROUTINE tra_adv( kt, Kbb, Kmm, pts, Krhs, pau, pav, paw )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv  ***
@@ -88,14 +88,16 @@
       !! ** Method  : - Update (uu(:,:,:,Krhs),vv(:,:,:,Krhs)) with the advection term following nadv
       !!----------------------------------------------------------------------
-      INTEGER                                  , INTENT(in)    :: kt             ! ocean time-step index
-      INTEGER                                  , INTENT(in)    :: Kbb, Kmm, Krhs ! time level indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers and RHS of tracer equation
+      INTEGER                                     , INTENT(in   ) ::   kt             ! ocean time-step index
+      INTEGER                                     , INTENT(in   ) ::   Kbb, Kmm, Krhs ! time level indices
+      REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in   ) ::   pau, pav, paw  ! advective velocity
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt)   , INTENT(inout) ::   pts            ! active tracers and RHS of tracer equation
       !
       INTEGER ::   ji, jj, jk   ! dummy loop index
+      REAL(wp), DIMENSION(:,:,:), POINTER ::   zptu, zptv, zptw
       ! TEMP: [tiling] This change not necessary and can be A2D(nn_hls) if using XIOS (subdomain support)
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zuu, zvv, zww   ! 3D workspace
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds
+      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
-      LOGICAL :: lskip
+      LOGICAL ::   lskip
       !!----------------------------------------------------------------------
       !
@@ -119,17 +121,29 @@
          ENDIF
       ENDIF
+      !
       IF( .NOT. lskip ) THEN
-         !                                         !==  effective transport  ==!
+         !                                         !==  effective advective transport  ==!
+         !
+         IF( PRESENT( pau ) ) THEN     ! RK3: advective velocity (pau,pav,paw) /= advected velocity (uu,vv,ww)
+            zptu => pau(:,:,:)
+            zptv => pav(:,:,:)
+            zptw => paw(:,:,:)
+         ELSE                          ! MLF: advective velocity = (uu,vv,ww)
+            zptu => uu(:,:,:,Kmm)
+            zptv => vv(:,:,:,Kmm)
+            zptw => ww(:,:,:    )
+         ENDIF
+         !
          IF( ln_wave .AND. ln_sdw )  THEN
             DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 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) )
-               zww(ji,jj,jk) = e1e2t(ji,jj)                     * ( ww(ji,jj,jk)     + wsd(ji,jj,jk) )
+               zuu(ji,jj,jk) = e2u  (ji,jj) * e3u(ji,jj,jk,Kmm) * ( zptu(ji,jj,jk) + usd(ji,jj,jk) )
+               zvv(ji,jj,jk) = e1v  (ji,jj) * e3v(ji,jj,jk,Kmm) * ( zptv(ji,jj,jk) + vsd(ji,jj,jk) )
+               zww(ji,jj,jk) = e1e2t(ji,jj)                     * ( zptw(ji,jj,jk) + wsd(ji,jj,jk) )
             END_3D
          ELSE
             DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 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)
-               zww(ji,jj,jk) = e1e2t(ji,jj)                     * ww(ji,jj,jk)
+               zuu(ji,jj,jk) = e2u  (ji,jj) * e3u(ji,jj,jk,Kmm) * zptu(ji,jj,jk)               ! eulerian transport only
+               zvv(ji,jj,jk) = e1v  (ji,jj) * e3v(ji,jj,jk,Kmm) * zptv(ji,jj,jk)
+               zww(ji,jj,jk) = e1e2t(ji,jj)                     * zptw(ji,jj,jk)
             END_3D
          ENDIF
@@ -275,6 +289,6 @@
          WRITE(numout,*) '      No advection on T & S                     ln_traadv_OFF = ', ln_traadv_OFF
          WRITE(numout,*) '      centered scheme                           ln_traadv_cen = ', ln_traadv_cen
-         WRITE(numout,*) '            horizontal 2nd/4th order               nn_cen_h   = ', nn_fct_h
-         WRITE(numout,*) '            vertical   2nd/4th order               nn_cen_v   = ', nn_fct_v
+         WRITE(numout,*) '            horizontal 2nd/4th order               nn_cen_h   = ', nn_cen_h
+         WRITE(numout,*) '            vertical   2nd/4th order               nn_cen_v   = ', nn_cen_v
          WRITE(numout,*) '      Flux Corrected Transport scheme           ln_traadv_fct = ', ln_traadv_fct
          WRITE(numout,*) '            horizontal 2nd/4th order               nn_fct_h   = ', nn_fct_h
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/oce.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/oce.F90	(revision 14353)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/oce.F90	(revision 14381)
@@ -21,6 +21,6 @@
    !! dynamics and tracer fields
    !! --------------------------                            
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:)   ::   uu   ,  vv     !: horizontal velocities        [m/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   ww             !: vertical velocity            [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:), TARGET   ::   uu   ,  vv     !: horizontal velocities        [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  , TARGET   ::   ww             !: vertical velocity            [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   wi             !: vertical vel. (adaptive-implicit) [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   hdiv           !: horizontal divergence        [s-1]
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stpmlf.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stpmlf.F90	(revision 14353)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stpmlf.F90	(revision 14381)
@@ -94,5 +94,6 @@
       !!----------------------------------------------------------------------
       INTEGER ::   ji, jj, jk, jtile   ! dummy loop indice
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zgdept
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)       ::   zgdept
+      REAL(wp), TARGET     , DIMENSION(jpi,jpj,jpk) ::   zau, zav, zaw   ! advective velocity
       !! ---------------------------------------------------------------------
 #if defined key_agrif
@@ -222,5 +223,9 @@
                          uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
                          vv(:,:,:,Nrhs) = 0._wp
-
+!!st                         
+                         zau(:,:,:) = uu(:,:,:,Nnn)        !!st   patch for MLF will be computed in RK3
+                         zav(:,:,:) = vv(:,:,:,Nnn)
+                         zaw(:,:,:) = ww(:,:,:    )
+!!st
       IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   &
                &         CALL dyn_asm_inc   ( kstp, Nbb, Nnn, uu, vv, Nrhs )  ! apply dynamics assimilation increment
@@ -230,5 +235,5 @@
                &         CALL Agrif_Sponge_dyn        ! momentum sponge
 #endif
-                         CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)	==> RHS
+                         CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs, zau, zav, zaw )  ! advection (VF or FF)	==> RHS
                          CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity           	==> RHS
                          CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing
@@ -247,6 +252,9 @@
          IF( ln_zad_Aimp )  CALL wAimp      ( kstp,      Nnn )                      ! Adaptive-implicit vertical advection partitioning
       ENDIF
-
-
+!!st
+                            zau(:,:,:) = uu(:,:,:,Nnn)        !!st   patch for MLF will be computed in RK3
+                            zav(:,:,:) = vv(:,:,:,Nnn)
+                            zaw(:,:,:) = ww(:,:,:    )
+!!st
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! cool skin
@@ -313,5 +321,5 @@
          IF( ln_tile    )   CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile )
 
-                            CALL tra_adv    ( kstp, Nbb, Nnn, ts, Nrhs )  ! hor. + vert. advection	==> RHS
+                            CALL tra_adv    ( kstp, Nbb, Nnn, ts, Nrhs, zau, zav, zaw )  ! hor. + vert. advection	==> RHS
          IF( ln_zdfmfc  )   CALL tra_mfc    ( kstp, Nbb,      ts, Nrhs )  ! Mass Flux Convection
          IF( ln_zdfosm  ) THEN
