Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90	(revision 12580)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90	(revision 12581)
@@ -45,4 +45,5 @@
    PUBLIC  dom_qe_sf_update  ! called by step.F90
    PUBLIC  dom_qe_interpol   ! called by dynnxt.F90
+   PUBLIC  dom_qe_r3c        ! called by steplf.F90
 
    !                                                      !!* Namelist nam_vvl
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatf.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatf.F90	(revision 12580)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatf.F90	(revision 12581)
@@ -13,5 +13,5 @@
    !!             -   !  2002-10  (C. Talandier, A-M. Treguier) Open boundary cond.
    !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization
-   !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines. 
+   !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines.
    !!            3.2  !  2009-06  (G. Madec, R.Benshila)  re-introduce the vvl option
    !!            3.3  !  2010-09  (D. Storkey, E.O'Dea) Bug fix for BDY module
@@ -22,5 +22,5 @@
    !!            4.1  !  2019-08  (A. Coward, D. Storkey) Rename dynnxt.F90 -> dynatf.F90. Now just does time filtering.
    !!-------------------------------------------------------------------------
-  
+
    !!----------------------------------------------------------------------------------------------
    !!   dyn_atf       : apply Asselin time filtering to "now" velocities and vertical scale factors
@@ -42,5 +42,5 @@
    USE trdken         ! trend manager: kinetic energy
    USE isf_oce   , ONLY: ln_isf     ! ice shelf
-   USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine 
+   USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine
    !
    USE in_out_manager ! I/O manager
@@ -63,5 +63,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$ 
+   !! $Id$
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
@@ -71,6 +71,6 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_atf  ***
-      !!                   
-      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary 
+      !!
+      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary
       !!             condition on the after velocity and apply the Asselin time
       !!             filter to the now fields.
@@ -79,5 +79,5 @@
       !!             estimate (ln_dynspg_ts=T)
       !!
-      !!              * Apply lateral boundary conditions on after velocity 
+      !!              * Apply lateral boundary conditions on after velocity
       !!             at the local domain boundaries through lbc_lnk call,
       !!             at the one-way open boundaries (ln_bdy=T),
@@ -86,5 +86,5 @@
       !!              * Apply the Asselin time filter to the now fields
       !!             arrays to start the next time step:
-      !!                (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm)) 
+      !!                (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm))
       !!                                    + atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Kaa),pvv(Kaa)) - 2 (puu(Kmm),pvv(Kmm)) ]
       !!             Note that with flux form advection and non linear free surface,
@@ -92,5 +92,5 @@
       !!             As a result, dyn_atf MUST be called after tra_atf.
       !!
-      !! ** Action :   puu(Kmm),pvv(Kmm)   filtered now horizontal velocity 
+      !! ** Action :   puu(Kmm),pvv(Kmm)   filtered now horizontal velocity
       !!----------------------------------------------------------------------
       INTEGER                             , INTENT(in   ) :: kt               ! ocean time-step index
@@ -103,5 +103,5 @@
       REAL(wp) ::   zve3a, zve3n, zve3b, z1_2dt   !   -      -
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve, zwfld
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3t_f, ze3u_f, ze3v_f, zua, zva 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3t_f, ze3u_f, ze3v_f, zua, zva
       !!----------------------------------------------------------------------
       !
@@ -131,17 +131,17 @@
          !
          IF( .NOT.ln_bt_fw ) THEN
-            ! Remove advective velocity from "now velocities" 
-            ! prior to asselin filtering     
-            ! In the forward case, this is done below after asselin filtering   
-            ! so that asselin contribution is removed at the same time 
+            ! Remove advective velocity from "now velocities"
+            ! prior to asselin filtering
+            ! In the forward case, this is done below after asselin filtering
+            ! so that asselin contribution is removed at the same time
             DO jk = 1, jpkm1
                puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk)
                pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk)
-            END DO  
+            END DO
          ENDIF
       ENDIF
 
       ! Update after velocity on domain lateral boundaries
-      ! --------------------------------------------------      
+      ! --------------------------------------------------
 # if defined key_agrif
       CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
@@ -157,5 +157,5 @@
       !
       IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
-         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step 
+         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
          IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
          !
@@ -177,6 +177,6 @@
       ! Time filter and swap of dynamics arrays
       ! ------------------------------------------
-         
-      IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN    !* Leap-Frog : Asselin time filter 
+
+      IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN    !* Leap-Frog : Asselin time filter
          !                                ! =============!
          IF( ln_linssh ) THEN             ! Fixed volume !
@@ -202,5 +202,5 @@
             DO jk = 1, jpkm1
                ze3t_f(:,:,jk) = ze3t_f(:,:,jk) - zcoef * zwfld(:,:) * tmask(:,:,jk) &
-                              &                        * pe3t(:,:,jk,Kmm) / ( ht(:,:) + 1._wp - ssmask(:,:) ) 
+                              &                        * pe3t(:,:,jk,Kmm) / ( ht(:,:) + 1._wp - ssmask(:,:) )
             END DO
             !
@@ -239,5 +239,5 @@
                   pvv(ji,jj,jk,Kmm) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
                END_3D
-               pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1)  
+               pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1)
                pe3v(:,:,1:jpkm1,Kmm) = ze3v_f(:,:,1:jpkm1)
                !
@@ -250,10 +250,10 @@
          IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
             ! Revert filtered "now" velocities to time split estimate
-            ! Doing it here also means that asselin filter contribution is removed  
+            ! Doing it here also means that asselin filter contribution is removed
             zue(:,:) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1)
-            zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)    
+            zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)
             DO jk = 2, jpkm1
                zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk)
-               zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)    
+               zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)
             END DO
             DO jk = 1, jpkm1
@@ -307,5 +307,5 @@
       IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' nxt  - puu(:,:,:,Kaa): ', mask1=umask,   &
          &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=' pvv(:,:,:,Kaa): '       , mask2=vmask )
-      ! 
+      !
       IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
       IF( l_trddyn     )   DEALLOCATE( zua, zva )
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatfLF.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatfLF.F90	(revision 12581)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatfLF.F90	(revision 12581)
@@ -0,0 +1,298 @@
+MODULE dynatfLF
+   !!=========================================================================
+   !!                       ***  MODULE  dynatfLF  ***
+   !! Ocean dynamics: time filtering
+   !!=========================================================================
+   !! History :  OPA  !  1987-02  (P. Andrich, D. L Hostis)  Original code
+   !!                 !  1990-10  (C. Levy, G. Madec)
+   !!            7.0  !  1993-03  (M. Guyon)  symetrical conditions
+   !!            8.0  !  1997-02  (G. Madec & M. Imbard)  opa, release 8.0
+   !!            8.2  !  1997-04  (A. Weaver)  Euler forward step
+   !!             -   !  1997-06  (G. Madec)  lateral boudary cond., lbc routine
+   !!    NEMO    1.0  !  2002-08  (G. Madec)  F90: Free form and module
+   !!             -   !  2002-10  (C. Talandier, A-M. Treguier) Open boundary cond.
+   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization
+   !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines.
+   !!            3.2  !  2009-06  (G. Madec, R.Benshila)  re-introduce the vvl option
+   !!            3.3  !  2010-09  (D. Storkey, E.O'Dea) Bug fix for BDY module
+   !!            3.3  !  2011-03  (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL
+   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
+   !!            3.6  !  2014-04  (G. Madec) add the diagnostic of the time filter trends
+   !!            3.7  !  2015-11  (J. Chanut) Free surface simplification
+   !!            4.1  !  2019-08  (A. Coward, D. Storkey) Rename dynnxt.F90 -> dynatfLF.F90. Now just does time filtering.
+   !!-------------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------------------------------
+   !!   dyn_atfLF       : apply Asselin time filtering to "now" velocities and vertical scale factors
+   !!----------------------------------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE sbc_oce        ! Surface boundary condition: ocean fields
+   USE sbcrnf         ! river runoffs
+   USE phycst         ! physical constants
+   USE dynadv         ! dynamics: vector invariant versus flux form
+   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
+   USE domvvl         ! variable volume
+   USE bdy_oce   , ONLY: ln_bdy
+   USE bdydta         ! ocean open boundary conditions
+   USE bdydyn         ! ocean open boundary conditions
+   USE bdyvol         ! ocean open boundary condition (bdy_vol routines)
+   USE trd_oce        ! trends: ocean variables
+   USE trddyn         ! trend manager: dynamics
+   USE trdken         ! trend manager: kinetic energy
+   USE isf_oce   , ONLY: ln_isf     ! ice shelf
+   USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine
+   !
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
+   USE lbclnk         ! lateral boundary condition (or mpp link)
+   USE lib_mpp        ! MPP library
+   USE prtctl         ! Print control
+   USE timing         ! Timing
+#if defined key_agrif
+   USE agrif_oce_interp
+#endif
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC    dyn_atf_lf   ! routine called by step.F90
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dyn_atf_lf ( kt, Kbb, Kmm, Kaa, puu, pvv, pe3t, pe3u, pe3v )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_atf_lf  ***
+      !!
+      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary
+      !!             condition on the after velocity and apply the Asselin time
+      !!             filter to the now fields.
+      !!
+      !! ** Method  : * Ensure after velocities transport matches time splitting
+      !!             estimate (ln_dynspg_ts=T)
+      !!
+      !!              * Apply lateral boundary conditions on after velocity
+      !!             at the local domain boundaries through lbc_lnk call,
+      !!             at the one-way open boundaries (ln_bdy=T),
+      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
+      !!
+      !!              * Apply the Asselin time filter to the now fields
+      !!             arrays to start the next time step:
+      !!                (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm))
+      !!                                    + atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Kaa),pvv(Kaa)) - 2 (puu(Kmm),pvv(Kmm)) ]
+      !!             Note that with flux form advection and non linear free surface,
+      !!             the time filter is applied on thickness weighted velocity.
+      !!             As a result, dyn_atf_lf MUST be called after tra_atf.
+      !!
+      !! ** Action :   puu(Kmm),pvv(Kmm)   filtered now horizontal velocity
+      !!----------------------------------------------------------------------
+      INTEGER                             , INTENT(in   ) :: kt               ! ocean time-step index
+      INTEGER                             , INTENT(in   ) :: Kbb, Kmm, Kaa    ! before and after time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv         ! velocities to be time filtered
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: pe3t, pe3u, pe3v ! scale factors to be time filtered
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zue3a, zue3n, zue3b, zcoef    ! local scalars
+      REAL(wp) ::   zve3a, zve3n, zve3b, z1_2dt   !   -      -
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zua, zva
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing    )   CALL timing_start('dyn_atf_lf')
+      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
+      IF( l_trddyn     )   ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) )
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dyn_atf_lf : Asselin time filtering'
+         IF(lwp) WRITE(numout,*) '~~~~~~~'
+      ENDIF
+
+      IF ( ln_dynspg_ts ) THEN
+         ! Ensure below that barotropic velocities match time splitting estimate
+         ! Compute actual transport and replace it with ts estimate at "after" time step
+         zue(:,:) = pe3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
+         zve(:,:) = pe3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
+         DO jk = 2, jpkm1
+            zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
+            zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
+         END DO
+         DO jk = 1, jpkm1
+            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk)
+            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk)
+         END DO
+         !
+         IF( .NOT.ln_bt_fw ) THEN
+            ! Remove advective velocity from "now velocities"
+            ! prior to asselin filtering
+            ! In the forward case, this is done below after asselin filtering
+            ! so that asselin contribution is removed at the same time
+            DO jk = 1, jpkm1
+               puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk)
+               pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk)
+            END DO
+         ENDIF
+      ENDIF
+
+      ! Update after velocity on domain lateral boundaries
+      ! --------------------------------------------------
+# if defined key_agrif
+      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
+# endif
+      !
+      CALL lbc_lnk_multi( 'dynatfLF', puu(:,:,:,Kaa), 'U', -1., pvv(:,:,:,Kaa), 'V', -1. )     !* local domain boundaries
+      !
+      !                                !* BDY open boundaries
+      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa )
+      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only=.true. )
+
+!!$   Do we need a call to bdy_vol here??
+      !
+      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
+         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
+         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
+         !
+         !                                  ! Kinetic energy and Conversion
+         IF( ln_KE_trd  )   CALL trd_dyn( puu(:,:,:,Kaa), pvv(:,:,:,Kaa), jpdyn_ken, kt, Kmm )
+         !
+         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
+            zua(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) * z1_2dt
+            zva(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) * z1_2dt
+            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
+            CALL iom_put( "vtrd_tot", zva )
+         ENDIF
+         !
+         zua(:,:,:) = puu(:,:,:,Kmm)             ! save the now velocity before the asselin filter
+         zva(:,:,:) = pvv(:,:,:,Kmm)             ! (caution: there will be a shift by 1 timestep in the
+         !                                  !  computation of the asselin filter trends)
+      ENDIF
+
+      ! Time filter and swap of dynamics arrays
+      ! ------------------------------------------
+
+      IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN    !* Leap-Frog : Asselin time filter
+         !                                ! =============!
+         IF( ln_linssh ) THEN             ! Fixed volume !
+            !                             ! =============!
+            DO_3D_11_11( 1, jpkm1 )
+               puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )
+               pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )
+            END_3D
+            !                             ! ================!
+         ELSE                             ! Variable volume !
+            !                             ! ================!
+            ! Time-filtered scale factor at t-points
+            ! ----------------------------------------------------
+            DO jk = 1, jpk                                          ! filtered scale factor at T-points
+               pe3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + r3t_f(:,:) * tmask(:,:,jk) )
+            END DO
+            !
+            !
+            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
+               ! Before filtered scale factor at (u/v)-points
+               DO jk = 1, jpk
+                  pe3u(:,:,jk,Kmm) = e3u_0(:,:,jk) * ( 1._wp + r3u_f(:,:) * umask(:,:,jk) )
+                  pe3v(:,:,jk,Kmm) = e3v_0(:,:,jk) * ( 1._wp + r3v_f(:,:) * vmask(:,:,jk) )
+               END DO
+               !
+               DO_3D_11_11( 1, jpkm1 )
+                  puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )
+                  pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )
+               END_3D
+               !
+            ELSE                          ! Asselin filter applied on thickness weighted velocity
+               !
+               DO_3D_11_11( 1, jpkm1 )
+                  zue3a = pe3u(ji,jj,jk,Kaa) * puu(ji,jj,jk,Kaa)
+                  zve3a = pe3v(ji,jj,jk,Kaa) * pvv(ji,jj,jk,Kaa)
+                  zue3n = pe3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm)
+                  zve3n = pe3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm)
+                  zue3b = pe3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb)
+                  zve3b = pe3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb)
+                  !                                                 ! filtered scale factor at U-,V-points
+                  pe3u(ji,jj,jk,Kmm) = e3u_0(ji,jj,jk) * ( 1._wp + r3u_f(ji,jj) * umask(ji,jj,jk) )
+                  pe3v(ji,jj,jk,Kmm) = e3v_0(ji,jj,jk) * ( 1._wp + r3v_f(ji,jj) * vmask(ji,jj,jk) )
+                  !
+                  puu(ji,jj,jk,Kmm) = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / pe3u(ji,jj,jk,Kmm) !!st ze3u_f(ji,jj,jk)
+                  pvv(ji,jj,jk,Kmm) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / pe3v(ji,jj,jk,Kmm) !!st ze3v_f(ji,jj,jk)
+               END_3D
+               !
+            ENDIF
+            !
+         ENDIF
+         !
+         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
+            ! Revert filtered "now" velocities to time split estimate
+            ! Doing it here also means that asselin filter contribution is removed
+            zue(:,:) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1)
+            zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)
+            DO jk = 2, jpkm1
+               zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk)
+               zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)
+            END DO
+            DO jk = 1, jpkm1
+               puu(:,:,jk,Kmm) = puu(:,:,jk,Kmm) - (zue(:,:) * r1_hu(:,:,Kmm) - uu_b(:,:,Kmm)) * umask(:,:,jk)
+               pvv(:,:,jk,Kmm) = pvv(:,:,jk,Kmm) - (zve(:,:) * r1_hv(:,:,Kmm) - vv_b(:,:,Kmm)) * vmask(:,:,jk)
+            END DO
+         ENDIF
+         !
+      ENDIF ! neuler /= 0
+      !
+      ! Set "now" and "before" barotropic velocities for next time step:
+      ! JC: Would be more clever to swap variables than to make a full vertical
+      ! integration
+      !
+      IF(.NOT.ln_linssh ) THEN
+         hu(:,:,Kmm) = pe3u(:,:,1,Kmm ) * umask(:,:,1)
+         hv(:,:,Kmm) = pe3v(:,:,1,Kmm ) * vmask(:,:,1)
+         DO jk = 2, jpkm1
+            hu(:,:,Kmm) = hu(:,:,Kmm) + pe3u(:,:,jk,Kmm ) * umask(:,:,jk)
+            hv(:,:,Kmm) = hv(:,:,Kmm) + pe3v(:,:,jk,Kmm ) * vmask(:,:,jk)
+         END DO
+         r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
+         r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) )
+      ENDIF
+      !
+      uu_b(:,:,Kaa) = pe3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
+      uu_b(:,:,Kmm) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1)
+      vv_b(:,:,Kaa) = pe3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
+      vv_b(:,:,Kmm) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)
+      DO jk = 2, jpkm1
+         uu_b(:,:,Kaa) = uu_b(:,:,Kaa) + pe3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
+         uu_b(:,:,Kmm) = uu_b(:,:,Kmm) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk)
+         vv_b(:,:,Kaa) = vv_b(:,:,Kaa) + pe3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
+         vv_b(:,:,Kmm) = vv_b(:,:,Kmm) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)
+      END DO
+      uu_b(:,:,Kaa) = uu_b(:,:,Kaa) * r1_hu(:,:,Kaa)
+      vv_b(:,:,Kaa) = vv_b(:,:,Kaa) * r1_hv(:,:,Kaa)
+      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm)
+      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm)
+      !
+      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
+         CALL iom_put(  "ubar", uu_b(:,:,Kmm) )
+         CALL iom_put(  "vbar", vv_b(:,:,Kmm) )
+      ENDIF
+      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
+         zua(:,:,:) = ( puu(:,:,:,Kmm) - zua(:,:,:) ) * z1_2dt
+         zva(:,:,:) = ( pvv(:,:,:,Kmm) - zva(:,:,:) ) * z1_2dt
+         CALL trd_dyn( zua, zva, jpdyn_atf, kt, Kmm )
+      ENDIF
+      !
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' nxt  - puu(:,:,:,Kaa): ', mask1=umask,   &
+         &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=' pvv(:,:,:,Kaa): '       , mask2=vmask )
+      !
+      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
+      IF( l_trddyn     )   DEALLOCATE( zua, zva )
+      IF( ln_timing    )   CALL timing_stop('dyn_atf_lf')
+      !
+   END SUBROUTINE dyn_atf_lf
+
+   !!=========================================================================
+END MODULE dynatfLF
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatf.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatf.F90	(revision 12580)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatf.F90	(revision 12581)
@@ -26,5 +26,5 @@
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain variables 
+   USE dom_oce         ! ocean space and time domain variables
    USE sbc_oce         ! surface boundary condition: ocean
    USE sbcrnf          ! river runoffs
@@ -33,5 +33,5 @@
    USE domvvl          ! variable volume
    USE trd_oce         ! trends: ocean variables
-   USE trdtra          ! trends manager: tracers 
+   USE trdtra          ! trends manager: tracers
    USE traqsr          ! penetrative solar radiation (needed for nksr)
    USE phycst          ! physical constant
@@ -69,14 +69,14 @@
       !!                   ***  ROUTINE traatf  ***
       !!
-      !! ** Purpose :   Apply the boundary condition on the after temperature  
+      !! ** Purpose :   Apply the boundary condition on the after temperature
       !!             and salinity fields and add the Asselin time filter on now fields.
-      !! 
-      !! ** Method  :   At this stage of the computation, ta and sa are the 
+      !!
+      !! ** Method  :   At this stage of the computation, ta and sa are the
       !!             after temperature and salinity as the time stepping has
       !!             been performed in trazdf_imp or trazdf_exp module.
       !!
-      !!              - Apply lateral boundary conditions on (ta,sa) 
-      !!             at the local domain   boundaries through lbc_lnk call, 
-      !!             at the one-way open boundaries (ln_bdy=T), 
+      !!              - Apply lateral boundary conditions on (ta,sa)
+      !!             at the local domain   boundaries through lbc_lnk call,
+      !!             at the one-way open boundaries (ln_bdy=T),
       !!             at the AGRIF zoom   boundaries (lk_agrif=T)
       !!
@@ -88,5 +88,5 @@
       INTEGER                                  , INTENT(in   ) :: kt             ! ocean time-step index
       INTEGER                                  , INTENT(in   ) :: Kbb, Kmm, Kaa  ! time level indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers 
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers
       !!
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -104,5 +104,5 @@
 
       ! Update after tracer on domain lateral boundaries
-      ! 
+      !
 #if defined key_agrif
       CALL Agrif_tra                     ! AGRIF zoom boundaries
@@ -112,5 +112,5 @@
       !
       IF( ln_bdy )   CALL bdy_tra( kt, Kbb, pts, Kaa )  ! BDY open boundaries
- 
+
       ! set time step size (Euler/Leapfrog)
       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =        rdt   ! at nit000             (Euler)
@@ -119,13 +119,13 @@
 
       ! trends computation initialisation
-      IF( l_trdtra )   THEN                    
+      IF( l_trdtra )   THEN
          ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
          ztrdt(:,:,jpk) = 0._wp
          ztrds(:,:,jpk) = 0._wp
-         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend 
+         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend
             CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_zdfp, ztrdt )
             CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_zdfp, ztrds )
          ENDIF
-         ! total trend for the non-time-filtered variables. 
+         ! total trend for the non-time-filtered variables.
          zfact = 1.0 / rdt
          ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from pts(Kmm) terms
@@ -137,12 +137,12 @@
          CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_tot, ztrds )
          IF( ln_linssh ) THEN       ! linear sea surface height only
-            ! Store now fields before applying the Asselin filter 
+            ! Store now fields before applying the Asselin filter
             ! in order to calculate Asselin filter trend later.
-            ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kmm) 
+            ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kmm)
             ztrds(:,:,:) = pts(:,:,:,jp_sal,Kmm)
          ENDIF
       ENDIF
 
-      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping 
+      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping
          !
          IF (l_trdtra .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
@@ -156,5 +156,5 @@
       ELSE                                            ! Leap-Frog + Asselin filter time stepping
          !
-         IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000,      'TRA', pts, jpts )  ! linear free surface 
+         IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000,      'TRA', pts, jpts )  ! linear free surface
          ELSE                   ;   CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nit000, rdt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface
          ENDIF
@@ -164,8 +164,8 @@
                   &                    pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1.  )
          !
-      ENDIF     
-      !
-      IF( l_trdtra .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
-         zfact = 1._wp / r2dt             
+      ENDIF
+      !
+      IF( l_trdtra .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt
+         zfact = 1._wp / r2dt
          DO jk = 1, jpkm1
             ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kmm) - ztrdt(:,:,jk) ) * zfact
@@ -191,5 +191,5 @@
       !!
       !! ** Purpose :   fixed volume: apply the Asselin time filter to the "now" field
-      !! 
+      !!
       !! ** Method  : - Apply a Asselin time filter on now fields.
       !!
@@ -216,5 +216,5 @@
          !
          DO_3D_00_00( 1, jpkm1 )
-            ztn = pt(ji,jj,jk,jn,Kmm)                                    
+            ztn = pt(ji,jj,jk,jn,Kmm)
             ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb)  ! time laplacian on tracers
             !
@@ -231,6 +231,6 @@
       !!                   ***  ROUTINE tra_atf_vvl  ***
       !!
-      !! ** Purpose :   Time varying volume: apply the Asselin time filter  
-      !! 
+      !! ** Purpose :   Time varying volume: apply the Asselin time filter
+      !!
       !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
       !!             pt(Kmm)  = ( e3t(Kmm)*pt(Kmm) + atfp*[ e3t(Kbb)*pt(Kbb) - 2 e3t(Kmm)*pt(Kmm) + e3t_a*pt(Kaa) ] )
@@ -262,5 +262,5 @@
       ENDIF
       !
-      IF( cdtype == 'TRA' )  THEN   
+      IF( cdtype == 'TRA' )  THEN
          ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
          ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
@@ -268,6 +268,6 @@
       ELSE                          ! passive tracers case
          ll_traqsr  = .FALSE.          ! NO solar penetration
-         ll_rnf     = .FALSE.          ! NO river runoffs ????          !!gm BUG ?  
-         ll_isf     = .FALSE.          ! NO ice shelf melting/freezing  !!gm BUG ?? 
+         ll_rnf     = .FALSE.          ! NO river runoffs ????          !!gm BUG ?
+         ll_isf     = .FALSE.          ! NO ice shelf melting/freezing  !!gm BUG ??
       ENDIF
       !
@@ -279,5 +279,5 @@
       zfact1 = atfp * p2dt
       zfact2 = zfact1 * r1_rau0
-      DO jn = 1, kjpt      
+      DO jn = 1, kjpt
          DO_3D_00_00( 1, jpkm1 )
             ze3t_b = e3t(ji,jj,jk,Kbb)
@@ -296,7 +296,7 @@
             !
             ! Add asselin correction on scale factors:
-            zscale = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) 
-            ze3t_f = ze3t_f - zfact2 * zscale * ( emp_b(ji,jj) - emp(ji,jj) ) 
-            IF ( ll_rnf ) ze3t_f = ze3t_f + zfact2 * zscale * (    rnf_b(ji,jj) -    rnf(ji,jj) ) 
+            zscale = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )
+            ze3t_f = ze3t_f - zfact2 * zscale * ( emp_b(ji,jj) - emp(ji,jj) )
+            IF ( ll_rnf ) ze3t_f = ze3t_f + zfact2 * zscale * (    rnf_b(ji,jj) -    rnf(ji,jj) )
             IF ( ll_isf ) THEN
                IF ( ln_isfcav_mlt ) ze3t_f = ze3t_f - zfact2 * zscale * ( fwfisf_cav_b(ji,jj) - fwfisf_cav(ji,jj) )
@@ -304,15 +304,15 @@
             ENDIF
             !
-            IF( jk == mikt(ji,jj) ) THEN           ! first level 
+            IF( jk == mikt(ji,jj) ) THEN           ! first level
                ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
             ENDIF
             !
             ! solar penetration (temperature only)
-            IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            & 
-               &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
+            IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            &
+               &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) )
                !
             !
             IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
-               &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
+               &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) &
                &                              * e3t(ji,jj,jk,Kmm) / h_rnf(ji,jj)
 
@@ -328,5 +328,5 @@
                         &                     * e3t(ji,jj,jk,Kmm) / rhisf_tbl_cav(ji,jj)
                   END IF
-                  ! level partially include in Losch_2008 ice shelf boundary layer 
+                  ! level partially include in Losch_2008 ice shelf boundary layer
                   IF ( jk == misfkb_cav(ji,jj) ) THEN
                      ztc_f  = ztc_f  - zfact1 * ( risf_cav_tsc(ji,jj,jn) - risf_cav_tsc_b(ji,jj,jn) )  &
@@ -342,5 +342,5 @@
                             &                 * e3t(ji,jj,jk,Kmm) / rhisf_tbl_par(ji,jj)
                   END IF
-                  ! level partially include in Losch_2008 ice shelf boundary layer 
+                  ! level partially include in Losch_2008 ice shelf boundary layer
                   IF ( jk == misfkb_par(ji,jj) ) THEN
                      ztc_f  = ztc_f  - zfact1 * ( risf_par_tsc(ji,jj,jn) - risf_par_tsc_b(ji,jj,jn) )  &
@@ -371,9 +371,9 @@
             !
          END_3D
-         ! 
+         !
       END DO
       !
       IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) )   THEN
-         IF( l_trdtra .AND. cdtype == 'TRA' ) THEN 
+         IF( l_trdtra .AND. cdtype == 'TRA' ) THEN
             CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) )
             CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) )
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatfLF.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatfLF.F90	(revision 12581)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatfLF.F90	(revision 12581)
@@ -0,0 +1,385 @@
+MODULE traatfLF
+   !!======================================================================
+   !!                       ***  MODULE  traatfLF  ***
+   !! Ocean active tracers:  Asselin time filtering for temperature and salinity
+   !!======================================================================
+   !! History :  OPA  !  1991-11  (G. Madec)  Original code
+   !!            7.0  !  1993-03  (M. Guyon)  symetrical conditions
+   !!            8.0  !  1996-02  (G. Madec & M. Imbard)  opa release 8.0
+   !!             -   !  1996-04  (A. Weaver)  Euler forward step
+   !!            8.2  !  1999-02  (G. Madec, N. Grima)  semi-implicit pressure grad.
+   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module
+   !!             -   !  2002-11  (C. Talandier, A-M Treguier) Open boundaries
+   !!             -   !  2005-04  (C. Deltel) Add Asselin trend in the ML budget
+   !!            2.0  !  2006-02  (L. Debreu, C. Mazauric) Agrif implementation
+   !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf
+   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  semi-implicit hpg with asselin filter + modified LF-RA
+   !!             -   !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
+   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rename tranxt.F90 -> traatfLF.F90. Now only does time filtering.
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   tra_atf       : time filtering on tracers
+   !!   tra_atf_fix   : time filtering on tracers : fixed    volume case
+   !!   tra_atf_vvl   : time filtering on tracers : variable volume case
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables
+   USE sbc_oce         ! surface boundary condition: ocean
+   USE sbcrnf          ! river runoffs
+   USE isf_oce         ! ice shelf melting
+   USE zdf_oce         ! ocean vertical mixing
+   USE domvvl          ! variable volume
+   USE trd_oce         ! trends: ocean variables
+   USE trdtra          ! trends manager: tracers
+   USE traqsr          ! penetrative solar radiation (needed for nksr)
+   USE phycst          ! physical constant
+   USE ldftra          ! lateral physics : tracers
+   USE ldfslp          ! lateral physics : slopes
+   USE bdy_oce  , ONLY : ln_bdy
+   USE bdytra          ! open boundary condition (bdy_tra routine)
+   !
+   USE in_out_manager  ! I/O manager
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl          ! Print control
+   USE timing          ! Timing
+#if defined key_agrif
+   USE agrif_oce_interp
+#endif
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   tra_atf_lf       ! routine called by step.F90
+   PUBLIC   tra_atf_fix_lf   ! to be used in trcnxt !!st WARNING discrepancy here interpol is used
+   PUBLIC   tra_atf_qe_lf    ! to be used in trcnxt !!st WARNING discrepancy here interpol is used
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE tra_atf_lf( kt, Kbb, Kmm, Kaa, pts )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE traatfLF  ***
+      !!
+      !! ** Purpose :   Apply the boundary condition on the after temperature
+      !!             and salinity fields and add the Asselin time filter on now fields.
+      !!
+      !! ** Method  :   At this stage of the computation, ta and sa are the
+      !!             after temperature and salinity as the time stepping has
+      !!             been performed in trazdf_imp or trazdf_exp module.
+      !!
+      !!              - Apply lateral boundary conditions on (ta,sa)
+      !!             at the local domain   boundaries through lbc_lnk call,
+      !!             at the one-way open boundaries (ln_bdy=T),
+      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
+      !!
+      !!              - Update lateral boundary conditions on AGRIF children
+      !!             domains (lk_agrif=T)
+      !!
+      !! ** Action  : - ts(Kmm) time filtered
+      !!----------------------------------------------------------------------
+      INTEGER                                  , INTENT(in   ) :: kt             ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) :: Kbb, Kmm, Kaa  ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers
+      !!
+      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
+      REAL(wp) ::   zfact            ! local scalars
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start( 'tra_atf_lf')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_atf_lf : apply Asselin time filter to "now" fields'
+         IF(lwp) WRITE(numout,*) '~~~~~~~'
+      ENDIF
+
+      ! Update after tracer on domain lateral boundaries
+      !
+#if defined key_agrif
+      CALL Agrif_tra                     ! AGRIF zoom boundaries
+#endif
+      !                                              ! local domain boundaries  (T-point, unchanged sign)
+      CALL lbc_lnk_multi( 'traatfLF', pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. )
+      !
+      IF( ln_bdy )   CALL bdy_tra( kt, Kbb, pts, Kaa )  ! BDY open boundaries
+
+      ! set time step size (Euler/Leapfrog)
+      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =        rdt   ! at nit000             (Euler)
+      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp* rdt   ! at nit000 or nit000+1 (Leapfrog)
+      ENDIF
+
+      ! trends computation initialisation
+      IF( l_trdtra )   THEN
+         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
+         ztrdt(:,:,jpk) = 0._wp
+         ztrds(:,:,jpk) = 0._wp
+         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend
+            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_zdfp, ztrdt )
+            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_zdfp, ztrds )
+         ENDIF
+         ! total trend for the non-time-filtered variables.
+         zfact = 1.0 / rdt
+         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from pts(Kmm) terms
+         DO jk = 1, jpkm1
+            ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - pts(:,:,jk,jp_tem,Kmm)) * zfact
+            ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - pts(:,:,jk,jp_sal,Kmm)) * zfact
+         END DO
+         CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_tot, ztrdt )
+         CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_tot, ztrds )
+         IF( ln_linssh ) THEN       ! linear sea surface height only
+            ! Store now fields before applying the Asselin filter
+            ! in order to calculate Asselin filter trend later.
+            ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kmm)
+            ztrds(:,:,:) = pts(:,:,:,jp_sal,Kmm)
+         ENDIF
+      ENDIF
+
+      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping
+         !
+         IF (l_trdtra .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
+            !                                        ! Asselin filter is output by tra_atf_vvl that is not called on this time step
+            ztrdt(:,:,:) = 0._wp
+            ztrds(:,:,:) = 0._wp
+            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_atf, ztrdt )
+            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_atf, ztrds )
+         END IF
+         !
+      ELSE                                            ! Leap-Frog + Asselin filter time stepping
+         !
+         IF( ln_linssh ) THEN   ;   CALL tra_atf_fix_lf( kt, Kbb, Kmm, Kaa, nit000,      'TRA', pts, jpts )  ! linear free surface
+         ELSE                   ;   CALL tra_atf_qe_lf( kt, Kbb, Kmm, Kaa, nit000, rdt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface
+         ENDIF
+         !
+         CALL lbc_lnk_multi( 'traatfLF', pts(:,:,:,jp_tem,Kbb) , 'T', 1., pts(:,:,:,jp_sal,Kbb) , 'T', 1., &
+                  &                    pts(:,:,:,jp_tem,Kmm) , 'T', 1., pts(:,:,:,jp_sal,Kmm) , 'T', 1., &
+                  &                    pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1.  )
+         !
+      ENDIF
+      !
+      IF( l_trdtra .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt
+         zfact = 1._wp / r2dt
+         DO jk = 1, jpkm1
+            ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kmm) - ztrdt(:,:,jk) ) * zfact
+            ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kmm) - ztrds(:,:,jk) ) * zfact
+         END DO
+         CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_atf, ztrdt )
+         CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_atf, ztrds )
+      END IF
+      IF( l_trdtra )   DEALLOCATE( ztrdt , ztrds )
+      !
+      !                        ! control print
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kmm), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
+         &                                  tab3d_2=pts(:,:,:,jp_sal,Kmm), clinfo2=       ' Sn: ', mask2=tmask )
+      !
+      IF( ln_timing )   CALL timing_stop('tra_atf_lf')
+      !
+   END SUBROUTINE tra_atf_lf
+
+
+   SUBROUTINE tra_atf_fix_lf( kt, Kbb, Kmm, Kaa, kit000, cdtype, pt, kjpt )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE tra_atf_fix  ***
+      !!
+      !! ** Purpose :   fixed volume: apply the Asselin time filter to the "now" field
+      !!
+      !! ** Method  : - Apply a Asselin time filter on now fields.
+      !!
+      !! ** Action  : - pt(Kmm) ready for the next time step
+      !!----------------------------------------------------------------------
+      INTEGER                                  , INTENT(in   ) ::  kt            ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) ::  Kbb, Kmm, Kaa ! 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
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::  pt            ! tracer fields
+      !
+      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
+      REAL(wp) ::   ztn, ztd         ! local scalars
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == kit000 )  THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_atf_fix_lf : time filtering', cdtype
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      ENDIF
+      !
+      DO jn = 1, kjpt
+         !
+         DO_3D_00_00( 1, jpkm1 )
+            ztn = pt(ji,jj,jk,jn,Kmm)
+            ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb)  ! time laplacian on tracers
+            !
+            pt(ji,jj,jk,jn,Kmm) = ztn + atfp * ztd                      ! pt <-- filtered pt
+         END_3D
+         !
+      END DO
+      !
+   END SUBROUTINE tra_atf_fix_lf
+
+
+   SUBROUTINE tra_atf_qe_lf( kt, Kbb, Kmm, Kaa, kit000, p2dt, cdtype, pt, psbc_tc, psbc_tc_b, kjpt )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE tra_atf_vvl  ***
+      !!
+      !! ** Purpose :   Time varying volume: apply the Asselin time filter
+      !!
+      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
+      !!             pt(Kmm)  = ( e3t(Kmm)*pt(Kmm) + atfp*[ e3t(Kbb)*pt(Kbb) - 2 e3t(Kmm)*pt(Kmm) + e3t_a*pt(Kaa) ] )
+      !!                       /( e3t(Kmm)         + atfp*[ e3t(Kbb)         - 2 e3t(Kmm)         + e3t(Kaa)    ] )
+      !!
+      !! ** Action  : - pt(Kmm) ready for the next time step
+      !!----------------------------------------------------------------------
+      INTEGER                                  , INTENT(in   ) ::  kt        ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices
+      INTEGER                                  , INTENT(in   ) ::  kit000    ! first time step index
+      REAL(wp)                                 , INTENT(in   ) ::  p2dt      ! time-step
+      CHARACTER(len=3)                         , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator)
+      INTEGER                                  , INTENT(in   ) ::  kjpt      ! number of tracers
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::  pt        ! tracer fields
+      REAL(wp), DIMENSION(jpi,jpj    ,kjpt)    , INTENT(in   ) ::  psbc_tc   ! surface tracer content
+      REAL(wp), DIMENSION(jpi,jpj    ,kjpt)    , INTENT(in   ) ::  psbc_tc_b ! before surface tracer content
+      !
+      LOGICAL  ::   ll_traqsr, ll_rnf, ll_isf   ! local logical
+      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
+      REAL(wp) ::   zfact, zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
+      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d  !   -      -
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrd_atf
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == kit000 )  THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_atf_vvl : time filtering', cdtype
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      ENDIF
+      !
+      IF( cdtype == 'TRA' )  THEN
+         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
+         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
+         ll_isf     = ln_isf           ! active  tracers case  and  ice shelf melting
+      ELSE                          ! passive tracers case
+         ll_traqsr  = .FALSE.          ! NO solar penetration
+         ll_rnf     = .FALSE.          ! NO river runoffs ????          !!gm BUG ?
+         ll_isf     = .FALSE.          ! NO ice shelf melting/freezing  !!gm BUG ??
+      ENDIF
+      !
+      IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) )   THEN
+         ALLOCATE( ztrd_atf(jpi,jpj,jpk,kjpt) )
+         ztrd_atf(:,:,:,:) = 0.0_wp
+      ENDIF
+      zfact = 1._wp / p2dt
+      zfact1 = atfp * p2dt
+      zfact2 = zfact1 * r1_rau0
+      DO jn = 1, kjpt
+         DO_3D_00_00( 1, jpkm1 )
+            ze3t_b = e3t(ji,jj,jk,Kbb)
+            ze3t_n = e3t(ji,jj,jk,Kmm)
+            ze3t_a = e3t(ji,jj,jk,Kaa)
+            !                                         ! tracer content at Before, now and after
+            ztc_b  = pt(ji,jj,jk,jn,Kbb) * ze3t_b
+            ztc_n  = pt(ji,jj,jk,jn,Kmm) * ze3t_n
+            ztc_a  = pt(ji,jj,jk,jn,Kaa) * ze3t_a
+            !
+            ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
+            ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
+            !
+            ztc_f  = ztc_n  + atfp * ztc_d
+            !
+            ! Asselin correction on scale factors is done via ssh in r3t_f
+            ze3t_f = e3t_0(ji,jj,jk) * ( 1._wp + r3t_f(ji,jj) * tmask(ji,jj,jk) )
+
+            !
+            IF( jk == mikt(ji,jj) ) THEN           ! first level
+               ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
+            ENDIF
+            !
+            ! solar penetration (temperature only)
+            IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            &
+               &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) )
+               !
+            !
+            IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
+               &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) &
+               &                              * e3t(ji,jj,jk,Kmm) / h_rnf(ji,jj)
+
+            !
+            ! ice shelf
+            IF( ll_isf ) THEN
+               !
+               ! melt in the cavity
+               IF ( ln_isfcav_mlt ) THEN
+                  ! level fully include in the Losch_2008 ice shelf boundary layer
+                  IF ( jk >= misfkt_cav(ji,jj) .AND. jk < misfkb_cav(ji,jj) ) THEN
+                     ztc_f  = ztc_f  - zfact1 * ( risf_cav_tsc(ji,jj,jn) - risf_cav_tsc_b(ji,jj,jn) ) &
+                        &                     * e3t(ji,jj,jk,Kmm) / rhisf_tbl_cav(ji,jj)
+                  END IF
+                  ! level partially include in Losch_2008 ice shelf boundary layer
+                  IF ( jk == misfkb_cav(ji,jj) ) THEN
+                     ztc_f  = ztc_f  - zfact1 * ( risf_cav_tsc(ji,jj,jn) - risf_cav_tsc_b(ji,jj,jn) )  &
+                            &                 * e3t(ji,jj,jk,Kmm) / rhisf_tbl_cav(ji,jj) * rfrac_tbl_cav(ji,jj)
+                  END IF
+               END IF
+               !
+               ! parametrised melt (cavity closed)
+               IF ( ln_isfpar_mlt ) THEN
+                  ! level fully include in the Losch_2008 ice shelf boundary layer
+                  IF ( jk >= misfkt_par(ji,jj) .AND. jk < misfkb_par(ji,jj) ) THEN
+                     ztc_f  = ztc_f  - zfact1 * ( risf_par_tsc(ji,jj,jn) - risf_par_tsc_b(ji,jj,jn) )  &
+                            &                 * e3t(ji,jj,jk,Kmm) / rhisf_tbl_par(ji,jj)
+                  END IF
+                  ! level partially include in Losch_2008 ice shelf boundary layer
+                  IF ( jk == misfkb_par(ji,jj) ) THEN
+                     ztc_f  = ztc_f  - zfact1 * ( risf_par_tsc(ji,jj,jn) - risf_par_tsc_b(ji,jj,jn) )  &
+                            &                 * e3t(ji,jj,jk,Kmm) / rhisf_tbl_par(ji,jj) * rfrac_tbl_par(ji,jj)
+                  END IF
+               END IF
+               !
+               ! ice sheet coupling correction
+               IF ( ln_isfcpl ) THEN
+                  !
+                  ! at kt = nit000,  risfcpl_vol_n = 0 and risfcpl_vol_b = risfcpl_vol so contribution nul
+                  IF ( ln_rstart .AND. kt == nit000+1 ) THEN
+                     ztc_f  = ztc_f  + zfact1 * risfcpl_tsc(ji,jj,jk,jn) * r1_e1e2t(ji,jj)
+                     ! Shouldn't volume increment be spread according thanks to zscale  ?
+                  END IF
+                  !
+               END IF
+               !
+            END IF
+            !
+            ze3t_f = 1.e0 / ze3t_f
+            pt(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f    ! time filtered "now" field
+            !
+            IF( ( l_trdtra .and. cdtype == 'TRA' ) .OR. ( l_trdtrc .and. cdtype == 'TRC' ) ) THEN
+               ztrd_atf(ji,jj,jk,jn) = (ztc_f - ztc_n) * zfact/ze3t_n
+            ENDIF
+            !
+         END_3D
+         !
+      END DO
+      !
+      IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) )   THEN
+         IF( l_trdtra .AND. cdtype == 'TRA' ) THEN
+            CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) )
+            CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) )
+         ENDIF
+         IF( l_trdtrc .AND. cdtype == 'TRC' ) THEN
+            DO jn = 1, kjpt
+               CALL trd_tra( kt, Kmm, Kaa, cdtype, jn, jptra_atf, ztrd_atf(:,:,:,jn) )
+            END DO
+         ENDIF
+         DEALLOCATE( ztrd_atf )
+      ENDIF
+      !
+   END SUBROUTINE tra_atf_qe_lf
+
+   !!======================================================================
+END MODULE traatfLF
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/steplf.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/steplf.F90	(revision 12580)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/steplf.F90	(revision 12581)
@@ -41,4 +41,6 @@
    USE iom              ! xIOs server
    USE domqe
+   USE traatfLF          ! time filtering                   (tra_atf_lf routine)
+   USE dynatfLF          ! time filtering                   (dyn_atf_lf routine)
 
    IMPLICIT NONE
@@ -286,7 +288,8 @@
 !!
 !!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine
-                         CALL tra_atf       ( kstp, Nbb, Nnn, Naa, ts )                      ! time filtering of "now" tracer arrays
-                         CALL dyn_atf       ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v  )  ! time filtering of "now" velocities and scale factors
                          CALL ssh_atf       ( kstp, Nbb, Nnn, Naa, ssh )                     ! time filtering of "now" sea surface height
+                         CALL dom_qe_r3c    ( ssh(:,:,Nnn), r3t_f, r3u_f, r3v_f )            ! "now" ssh/h_0 ratio from filtrered ssh
+                         CALL tra_atf_lf    ( kstp, Nbb, Nnn, Naa, ts )                      ! time filtering of "now" tracer arrays
+                         CALL dyn_atf_lf    ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v  )  ! time filtering of "now" velocities and scale factors
       !
       ! Swap time levels
