Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ASM/asminc.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ASM/asminc.F90	(revision 12623)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ASM/asminc.F90	(revision 12624)
@@ -760,5 +760,5 @@
             !
             ssh(:,:,Kbb) = ssh(:,:,Kmm)                        ! Update before fields
-#if ! defined key_LF
+#if ! defined key_QCO
             e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
 #endif
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/CRS/crsini.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/CRS/crsini.F90	(revision 12623)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/CRS/crsini.F90	(revision 12624)
@@ -176,10 +176,4 @@
      
      !
-#if ! defined key_LF
-     ze3t(:,:,:) = e3t(:,:,:,Kmm)
-     ze3u(:,:,:) = e3u(:,:,:,Kmm)
-     ze3v(:,:,:) = e3v(:,:,:,Kmm)
-     ze3w(:,:,:) = e3w(:,:,:,Kmm)
-#else
      DO jk = 1, jpk
         ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
@@ -187,6 +181,5 @@
         ze3v(:,:,jk) = e3v(:,:,jk,Kmm)
         ze3w(:,:,jk) = e3w(:,:,jk,Kmm)
-     END DO 
-#endif 
+     END DO  
 
      !    3.d.2   Surfaces 
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 12623)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90	(revision 12624)
@@ -44,5 +44,5 @@
    PUBLIC  dom_qe_sf_nxt     ! called by steplf.F90
    PUBLIC  dom_qe_sf_update  ! called by steplf.F90
-   PUBLIC  dom_qe_interpol   ! called by dynnxt.F90
+   PUBLIC  dom_h_nxt         ! called by steplf.F90
    PUBLIC  dom_qe_r3c        ! called by steplf.F90
 
@@ -292,4 +292,46 @@
 
 
+   SUBROUTINE dom_h_nxt( kt, Kbb, Kmm, Kaa, kcall )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_qe_sf_nxt  ***
+      !!
+      !! ** Purpose :  - compute the after water heigh used in tra_zdf, dynnxt,
+      !!                 tranxt and dynspg routines
+      !!
+      !! ** Method  :  - z_star case:  Proportionnaly to the water column thickness.
+      !!
+      !! ** Action  :  - h(u/v) update wrt ssh/h(u/v)_0
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in )           ::   kt             ! time step
+      INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step
+      INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence
+      !
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
+      !
+      IF( ln_timing )   CALL timing_start('dom_h_nxt')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dom_h_nxt : compute after scale factors'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
+      ENDIF
+      !
+      ! *********************************** !
+      ! After depths at u- v points         !
+      ! *********************************** !
+      hu(:,:,Kaa) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kaa) )
+      hv(:,:,Kaa) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kaa) )
+      !                                        ! Inverse of the local depth
+      r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
+      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
+      !
+      IF( ln_timing )   CALL timing_stop('dom_h_nxt')
+      !
+   END SUBROUTINE dom_h_nxt
+
+
    SUBROUTINE dom_qe_sf_update( kt, Kbb, Kmm, Kaa )
       !!----------------------------------------------------------------------
@@ -398,102 +440,4 @@
       !
    END SUBROUTINE dom_qe_sf_update
-
-
-   SUBROUTINE dom_qe_interpol( pe3_in, pe3_out, pout )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE dom_qe_interpol  ***
-      !!
-      !! ** Purpose :   interpolate scale factors from one grid point to another
-      !!
-      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
-      !!                - horizontal interpolation: grid cell surface averaging
-      !!                - vertical interpolation: simple averaging
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::  pe3_in    ! input e3 to be interpolated
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  pe3_out   ! output interpolated e3
-      CHARACTER(LEN=*)                , INTENT(in   ) ::  pout      ! grid point of out scale factors
-      !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
-      !
-      INTEGER ::   ji, jj, jk                                       ! dummy loop indices
-      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F
-      !!----------------------------------------------------------------------
-      !
-      IF(ln_wd_il) THEN
-        zlnwd = 1.0_wp
-      ELSE
-        zlnwd = 0.0_wp
-      END IF
-      !
-      SELECT CASE ( pout )    !==  type of interpolation  ==!
-         !
-      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
-         DO_3D_10_10( 1, jpk )
-            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
-               &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
-               &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
-         END_3D
-         CALL lbc_lnk( 'domqe', pe3_out(:,:,:), 'U', 1._wp )
-         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
-         !
-      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
-         DO_3D_10_10( 1, jpk )
-            pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
-               &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
-               &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
-         END_3D
-         CALL lbc_lnk( 'domqe', pe3_out(:,:,:), 'V', 1._wp )
-         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
-         !
-      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
-         DO_3D_10_10( 1, jpk )
-            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
-               &                       *    r1_e1e2f(ji,jj)                                                  &
-               &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
-               &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
-         END_3D
-         CALL lbc_lnk( 'domqe', pe3_out(:,:,:), 'F', 1._wp )
-         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
-         !
-      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
-         !
-         !zlnwd = 1.0_wp
-         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
-         ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
-!!gm BUG? use here wmask in case of ISF ?  to be checked
-         DO jk = 2, jpk
-            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   &
-               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               &
-               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     &
-               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
-         END DO
-         !
-      CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean
-         !
-         !zlnwd = 1.0_wp
-         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
-         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
-!!gm BUG? use here wumask in case of ISF ?  to be checked
-         DO jk = 2, jpk
-            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
-               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              &
-               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
-               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
-         END DO
-         !
-      CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean
-         !
-         !zlnwd = 1.0_wp
-         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
-         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
-!!gm BUG? use here wvmask in case of ISF ?  to be checked
-         DO jk = 2, jpk
-            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
-               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              &
-               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
-               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
-         END DO
-      END SELECT
-      !
-   END SUBROUTINE dom_qe_interpol
 
 
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domzgr_substitute.h90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domzgr_substitute.h90	(revision 12623)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domzgr_substitute.h90	(revision 12624)
@@ -11,5 +11,5 @@
 !! Software governed by the CeCILL license (see ./LICENSE)
 !!----------------------------------------------------------------------
-#if defined key_LF
+#if defined key_QCO
 #   define  e3t(i,j,k,t)   (e3t_0(i,j,k)*(1.+r3t(i,j,t)*tmask(i,j,k)))
 #   define  e3u(i,j,k,t)   (e3u_0(i,j,k)*(1.+r3u(i,j,t)*umask(i,j,k)))
Index: MO/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 12623)
+++ 	(revision )
@@ -1,298 +1,0 @@
-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/DYN/dynatfQCO.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatfQCO.F90	(revision 12624)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatfQCO.F90	(revision 12624)
@@ -0,0 +1,272 @@
+MODULE dynatfqco
+   !!=========================================================================
+   !!                       ***  MODULE  dynatfqco  ***
+   !! 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_atf_qco       : 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_qco   ! routine called by step.F90
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dyn_atf_qco ( kt, Kbb, Kmm, Kaa, puu, pvv, pe3t, pe3u, pe3v )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_atf_qco  ***
+      !!
+      !! ** 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_qco')
+      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_qco : Asselin time filtering'
+         IF(lwp) WRITE(numout,*) '~~~~~~~'
+      ENDIF
+      !
+      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)
+                  !
+                  zue3a = ( 1._wp + r3u(ji,jj,Kaa) * umask(ji,jj,jk) ) * puu(ji,jj,jk,Kaa)
+                  zve3a = ( 1._wp + r3v(ji,jj,Kaa) * vmask(ji,jj,jk) ) * pvv(ji,jj,jk,Kaa)
+                  zue3n = ( 1._wp + r3u(ji,jj,Kmm) * umask(ji,jj,jk) ) * puu(ji,jj,jk,Kmm)
+                  zve3n = ( 1._wp + r3v(ji,jj,Kmm) * vmask(ji,jj,jk) ) * pvv(ji,jj,jk,Kmm)
+                  zue3b = ( 1._wp + r3u(ji,jj,Kbb) * umask(ji,jj,jk) ) * puu(ji,jj,jk,Kbb)
+                  zve3b = ( 1._wp + r3v(ji,jj,Kbb) * vmask(ji,jj,jk) ) * 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 ) ) / ( 1._wp + r3u_f(ji,jj)*umask(ji,jj,jk) ) !!st ze3u_f(ji,jj,jk)
+                  pvv(ji,jj,jk,Kmm) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ( 1._wp + r3v_f(ji,jj)*vmask(ji,jj,jk) ) !!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
+            zue(:,:) = e3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1)
+            zve(:,:) = e3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)
+            DO jk = 2, jpkm1
+               zue(:,:) = zue(:,:) + e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk)
+               zve(:,:) = zve(:,:) + e3v(:,:,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) = e3u(:,:,1,Kmm ) * umask(:,:,1)
+         hv(:,:,Kmm) = e3v(:,:,1,Kmm ) * vmask(:,:,1)
+         DO jk = 2, jpkm1
+            hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm ) * umask(:,:,jk)
+            hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,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) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
+      uu_b(:,:,Kmm) = e3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1)
+      vv_b(:,:,Kaa) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
+      vv_b(:,:,Kmm) = e3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)
+      DO jk = 2, jpkm1
+         uu_b(:,:,Kaa) = uu_b(:,:,Kaa) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
+         uu_b(:,:,Kmm) = uu_b(:,:,Kmm) + e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk)
+         vv_b(:,:,Kaa) = vv_b(:,:,Kaa) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
+         vv_b(:,:,Kmm) = vv_b(:,:,Kmm) + e3v(:,:,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_qco')
+      !
+   END SUBROUTINE dyn_atf_qco
+
+   !!=========================================================================
+END MODULE dynatfqco
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfcpl.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfcpl.F90	(revision 12623)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfcpl.F90	(revision 12624)
@@ -15,5 +15,5 @@
    USE isfutils, ONLY : debug
    USE lib_mpp , ONLY: mpp_sum, mpp_max ! mpp routine
-#if ! defined key_LF
+#if ! defined key_QCO
    USE domvvl  , ONLY: dom_vvl_zgr      ! vertical scale factor interpolation
 #else
@@ -117,5 +117,5 @@
       vv   (:,:,:,Kbb)   = vv   (:,:,:,Kmm)
       ssh (:,:,Kbb)     = ssh (:,:,Kmm)
-#if ! defined key_LF
+#if ! defined key_QCO
       e3t(:,:,:,Kbb)   = e3t(:,:,:,Kmm)
 #endif 
@@ -224,5 +224,5 @@
       IF(lwp) write(numout,*) 'isfcpl_ssh : recompute scale factor from ssh (new wet cell,Kmm)'
       IF(lwp) write(numout,*) '~~~~~~~~~~~'
-#if ! defined key_LF
+#if ! defined key_QCO
       DO jk = 1, jpk
          e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcice_cice.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcice_cice.F90	(revision 12623)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcice_cice.F90	(revision 12624)
@@ -233,5 +233,5 @@
 !!gm This should be put elsewhere....   (same remark for limsbc)
 !!gm especially here it is assumed zstar coordinate, but it can be ztilde....
-#if defined key_LF
+#if defined key_QCO
             IF( .NOT.ln_linssh )   CALL dom_qe_zgr( Kbb, Kmm, Kaa )   ! interpolation scale factor, depth and water column
 #else
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traadv.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traadv.F90	(revision 12623)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traadv.F90	(revision 12624)
@@ -65,5 +65,6 @@
    INTEGER, PARAMETER ::   np_UBS     = 4   ! 3rd order Upstream Biased Scheme
    INTEGER, PARAMETER ::   np_QCK     = 5   ! QUICK scheme
-   
+
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -103,7 +104,10 @@
       IF( ln_wave .AND. ln_sdw )  THEN
          DO jk = 1, jpkm1                                                       ! eulerian transport + Stokes Drift
-            zuu(:,:,jk) = e2u  (:,:) * e3u(:,:,jk,Kmm) * ( uu(:,:,jk,Kmm) + usd(:,:,jk) )
-            zvv(:,:,jk) = e1v  (:,:) * e3v(:,:,jk,Kmm) * ( vv(:,:,jk,Kmm) + vsd(:,:,jk) )
-            zww(:,:,jk) = e1e2t(:,:)                 * ( ww(:,:,jk) + wsd(:,:,jk) )
+            zuu(:,:,jk) =   &
+               &  e2u  (:,:) * e3u(:,:,jk,Kmm) * ( uu(:,:,jk,Kmm) + usd(:,:,jk) )
+            zvv(:,:,jk) =   & 
+               &  e1v  (:,:) * e3v(:,:,jk,Kmm) * ( vv(:,:,jk,Kmm) + vsd(:,:,jk) )
+            zww(:,:,jk) =   & 
+               &  e1e2t(:,:)                 * ( ww(:,:,jk) + wsd(:,:,jk) )
          END DO
       ELSE
Index: MO/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 12623)
+++ 	(revision )
@@ -1,385 +1,0 @@
-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/TRA/traatfQCO.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatfQCO.F90	(revision 12624)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatfQCO.F90	(revision 12624)
@@ -0,0 +1,392 @@
+MODULE traatfqco
+   !!======================================================================
+   !!                       ***  MODULE  traatfqco  ***
+   !! 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_qco       ! 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_qco_lf    ! to be used in trcnxt !!st WARNING discrepancy here interpol is used
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE tra_atf_qco( 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_qco')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_atf_qco : 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
+            ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kaa) * (1._wp + r3t(:,:,Kaa) * tmask(:,:,jk))/(1._wp + r3t(:,:,Kmm) * tmask(:,:,jk))  &
+               &            - pts(:,:,jk,jp_tem,Kmm) ) * zfact
+            ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kaa) * (1._wp + r3t(:,:,Kaa) * tmask(:,:,jk))/(1._wp + r3t(:,:,Kmm) * tmask(:,:,jk))  &
+               &            - 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_qco_lf( kt, Kbb, Kmm, Kaa, nit000, rdt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface
+         ENDIF
+         !
+         CALL lbc_lnk_multi( 'traatfqco', 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_qco')
+      !
+   END SUBROUTINE tra_atf_qco
+
+
+   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_qco_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_m*pt(Kmm) + atfp*[ e3t_b*pt(Kbb) - 2 e3t_m*pt(Kmm) + e3t_a*pt(Kaa) ] )
+      !!                       /( e3t_m         + atfp*[ e3t_b         - 2 e3t_m         + e3t_a    ] )
+      !!
+      !! ** 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_qco : 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_qco_lf
+
+   !!======================================================================
+END MODULE traatfqco
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfmxl.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfmxl.F90	(revision 12623)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfmxl.F90	(revision 12624)
@@ -38,4 +38,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
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 12624)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/stepLF.F90	(revision 12624)
@@ -0,0 +1,460 @@
+MODULE steplf
+   !!======================================================================
+   !!                       ***  MODULE step  ***
+   !! Time-stepping   : manager of the ocean, tracer and ice time stepping
+   !!======================================================================
+   !! History :  OPA  !  1991-03  (G. Madec)  Original code
+   !!             -   !  1991-11  (G. Madec)
+   !!             -   !  1992-06  (M. Imbard)  add a first output record
+   !!             -   !  1996-04  (G. Madec)  introduction of dynspg
+   !!             -   !  1996-04  (M.A. Foujols)  introduction of passive tracer
+   !!            8.0  !  1997-06  (G. Madec)  new architecture of call
+   !!            8.2  !  1997-06  (G. Madec, M. Imbard, G. Roullet)  free surface
+   !!             -   !  1999-02  (G. Madec, N. Grima)  hpg implicit
+   !!             -   !  2000-07  (J-M Molines, M. Imbard)  Open Bondary Conditions
+   !!   NEMO     1.0  !  2002-06  (G. Madec)  free form, suppress macro-tasking
+   !!             -   !  2004-08  (C. Talandier) New trends organization
+   !!             -   !  2005-01  (C. Ethe) Add the KPP closure scheme
+   !!             -   !  2005-11  (G. Madec)  Reorganisation of tra and dyn calls
+   !!             -   !  2006-01  (L. Debreu, C. Mazauric)  Agrif implementation
+   !!             -   !  2006-07  (S. Masson)  restart using iom
+   !!            3.2  !  2009-02  (G. Madec, R. Benshila)  reintroduicing z*-coordinate
+   !!             -   !  2009-06  (S. Masson, G. Madec)  TKE restart compatible with key_cpl
+   !!            3.3  !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
+   !!             -   !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA
+   !!            3.4  !  2011-04  (G. Madec, C. Ethe) Merge of dtatem and dtasal
+   !!            3.6  !  2012-07  (J. Simeon, G. Madec. C. Ethe)  Online coarsening of outputs
+   !!            3.6  !  2014-04  (F. Roquet, G. Madec) New equations of state
+   !!            3.6  !  2014-10  (E. Clementi, P. Oddo) Add Qiao vertical mixing in case of waves
+   !!            3.7  !  2014-10  (G. Madec)  LDF simplication
+   !!             -   !  2014-12  (G. Madec) remove KPP scheme
+   !!             -   !  2015-11  (J. Chanut) free surface simplification (remove filtered free surface)
+   !!            4.0  !  2017-05  (G. Madec)  introduction of the vertical physics manager (zdfphy)
+   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   stplf             : OPA system time-stepping
+   !!----------------------------------------------------------------------
+   USE step_oce         ! time stepping definition modules
+   !
+   USE iom              ! xIOs server
+   USE domqe
+   USE traatfqco         ! time filtering                   (tra_atf_qco routine)
+   USE dynatfqco         ! time filtering                   (dyn_atf_qco routine)
+   USE dynspg_ts         ! surface pressure gradient: split-explicit scheme (define un_adv)
+   USE bdydyn            ! ocean open boundary conditions (define bdy_dyn)
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   stplf   ! called by nemogcm.F90
+
+   !!----------------------------------------------------------------------
+   !! time level indices
+   !!----------------------------------------------------------------------
+   INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs          !! used by nemo_init
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: step.F90 12377 2020-02-12 14:39:06Z acc $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+#if defined key_agrif
+   RECURSIVE SUBROUTINE stplf( )
+      INTEGER             ::   kstp   ! ocean time-step index
+#else
+   SUBROUTINE stplf( kstp )
+      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
+#endif
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE stplf  ***
+      !!
+      !! ** Purpose : - Time stepping of OPA  (momentum and active tracer eqs.)
+      !!              - Time stepping of SI3 (dynamic and thermodynamic eqs.)
+      !!              - Time stepping of TRC  (passive tracer eqs.)
+      !!
+      !! ** Method  : -1- Update forcings and data
+      !!              -2- Update ocean physics
+      !!              -3- Compute the t and s trends
+      !!              -4- Update t and s
+      !!              -5- Compute the momentum trends
+      !!              -6- Update the horizontal velocity
+      !!              -7- Compute the diagnostics variables (rd,N2, hdiv,w)
+      !!              -8- Outputs and diagnostics
+      !!----------------------------------------------------------------------
+      INTEGER ::   ji, jj, jk   ! dummy loop indice
+      INTEGER ::   indic        ! error indicator if < 0
+!!gm kcall can be removed, I guess
+      INTEGER ::   kcall        ! optional integer argument (dom_vvl_sf_nxt)
+      !! ---------------------------------------------------------------------
+#if defined key_agrif
+      kstp = nit000 + Agrif_Nb_Step()
+      Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
+      IF( lk_agrif_debug ) THEN
+         IF( Agrif_Root() .and. lwp)   WRITE(*,*) '---'
+         IF(lwp)   WRITE(*,*) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint()
+      ENDIF
+      IF( kstp == nit000 + 1 )   lk_agrif_fstep = .FALSE.
+# if defined key_iomput
+      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context )
+# endif
+#endif
+      !
+      IF( ln_timing )   CALL timing_start('stplf')
+      !
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! update I/O and calendar
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                             indic = 0                ! reset to no error condition
+
+      IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
+                             CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including passible AGRIF zoom)
+         IF( lk_diamlr   )   CALL dia_mlr_iom_init    ! with additional setup for multiple-linear-regression analysis
+                             CALL iom_init_closedef
+         IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid
+      ENDIF
+      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
+                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell IOM we are at time step kstp
+      IF( ln_crs         )   CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell IOM we are at time step kstp
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Update external forcing (tides, open boundaries, ice shelf interaction and surface boundary condition (including sea-ice)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( ln_tide    )   CALL tide_update( kstp )                     ! update tide potential
+      IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                        ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)
+      IF( ln_bdy     )   CALL bdy_dta ( kstp, Nnn )                   ! update dynamic & tracer data at open boundaries
+      IF( ln_isf     )   CALL isf_stp ( kstp, Nnn )
+                         CALL sbc     ( kstp, Nbb, Nnn )              ! Sea Boundary Condition (including sea-ice)
+
+      !             !!st      !!!!!!!!!!!!!!!!!!!!!!!
+      !
+      ! emp = 0._wp
+      ! emp_b = 0._wp
+      ! qns = 0._wp
+      ! qsr = 0._wp
+      ! qns_b = 0._wp
+      !
+      !             !!st
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Update stochastic parameters and random T/S fluctuations
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      IF( ln_sto_eos ) CALL sto_par( kstp )          ! Stochastic parameters
+      IF( ln_sto_eos ) CALL sto_pts( ts(:,:,:,:,Nnn)  )          ! Random T/S fluctuations
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Ocean physics update
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !  THERMODYNAMICS
+                         CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nnn )       ! before local thermal/haline expension ratio at T-points
+                         CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn )       ! now    local thermal/haline expension ratio at T-points
+                         CALL bn2    ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
+                         CALL bn2    ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn  ) ! now    Brunt-Vaisala frequency
+
+      !  VERTICAL PHYSICS
+                         CALL zdf_phy( kstp, Nbb, Nnn, Nrhs )   ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
+
+      !  LATERAL  PHYSICS
+      !
+      IF( l_ldfslp ) THEN                             ! slope of lateral mixing
+                         CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )               ! before in situ density
+
+         IF( ln_zps .AND. .NOT. ln_isfcav)                                    &
+            &            CALL zps_hde    ( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
+            &                                          rhd, gru , grv    )       ! of t, s, rd at the last ocean level
+
+         IF( ln_zps .AND.       ln_isfcav)                                                &
+            &            CALL zps_hde_isf( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+            &                                          rhd, gru , grv , grui, grvi   )       ! of t, s, rd at the first ocean level
+         IF( ln_traldf_triad ) THEN
+                         CALL ldf_slp_triad( kstp, Nbb, Nnn )             ! before slope for triad operator
+         ELSE
+                         CALL ldf_slp     ( kstp, rhd, rn2b, Nbb, Nnn )   ! before slope for standard operator
+         ENDIF
+      ENDIF
+      !                                                                        ! eddy diffusivity coeff.
+      IF( l_ldftra_time .OR. l_ldfeiv_time )   CALL ldf_tra( kstp, Nbb, Nnn )  !       and/or eiv coeff.
+      IF( l_ldfdyn_time                    )   CALL ldf_dyn( kstp, Nbb )       ! eddy viscosity coeff.
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      !  Ocean dynamics : hdiv, ssh, e3, u, v, w
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+                            CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor)
+      IF( .NOT.ln_linssh )  CALL dom_qe_r3c    ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa) )
+      IF( .NOT.ln_linssh )  CALL dom_h_nxt     ( kstp, Nbb, Nnn,      Naa )    ! after vertical scale factors
+      !IF( .NOT.ln_linssh )  CALL dom_qe_sf_nxt ( kstp, Nbb, Nnn,      Naa )    ! after vertical scale factors
+                            CALL wzv           ( kstp, Nbb, Nnn, ww,  Naa )    ! now cross-level velocity
+      IF( ln_zad_Aimp )     CALL wAimp         ( kstp,      Nnn           )  ! Adaptive-implicit vertical advection partitioning
+                            CALL eos    ( ts(:,:,:,:,Nnn), rhd, rhop, gdept(:,:,:,Nnn) )  ! now in situ density for hpg computation
+
+
+                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
+                         vv(:,:,:,Nrhs) = 0._wp
+
+      IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   &
+               &         CALL dyn_asm_inc   ( kstp, Nbb, Nnn, uu, vv, Nrhs )  ! apply dynamics assimilation increment
+      IF( ln_bdy     )   CALL bdy_dyn3d_dmp ( kstp, Nbb,      uu, vv, Nrhs )  ! bdy damping trends
+#if defined key_agrif
+      IF(.NOT. Agrif_Root())  &
+               &         CALL Agrif_Sponge_dyn        ! momentum sponge
+#endif
+                         CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! 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
+      IF( ln_zdfosm  )   CALL dyn_osm( kstp,      Nnn      , uu, vv, Nrhs )  ! OSMOSIS non-local velocity fluxes ==> RHS
+                         CALL dyn_hpg( kstp,      Nnn      , uu, vv, Nrhs )  ! horizontal gradient of Hydrostatic pressure
+                         CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa )  ! surface pressure gradient
+
+                                                      ! With split-explicit free surface, since now transports have been updated and ssh(:,:,Nrhs) as well
+      IF( ln_dynspg_ts ) THEN                         ! vertical scale factors and vertical velocity need to be updated
+                            CALL div_hor       ( kstp, Nbb, Nnn )                ! Horizontal divergence  (2nd call in time-split case)
+         IF(.NOT.ln_linssh) CALL dom_qe_r3c    ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )
+         IF(.NOT.ln_linssh) CALL dom_qe_sf_nxt ( kstp, Nbb, Nnn, Naa, kcall=2 )  ! after vertical scale factors (update depth average component)
+         !IF(.NOT.ln_linssh) CALL dom_h_nxt ( kstp, Nbb, Nnn, Naa, kcall=2 )  ! after vertical scale factors (update depth average component)
+      ENDIF
+                            CALL dyn_zdf    ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion
+      IF( ln_dynspg_ts ) THEN                                                       ! vertical scale factors and vertical velocity need to be updated
+                            CALL wzv        ( kstp, Nbb, Nnn, ww, Naa )             ! now cross-level velocity
+         IF( ln_zad_Aimp )  CALL wAimp      ( kstp,      Nnn )                      ! Adaptive-implicit vertical advection partitioning
+      ENDIF
+
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! cool skin
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF ( ln_diurnal )  CALL diurnal_layers( kstp )
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! diagnostics and outputs
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( ln_floats  )   CALL flo_stp   ( kstp, Nbb, Nnn )      ! drifting Floats
+      IF( ln_diacfl  )   CALL dia_cfl   ( kstp,      Nnn )      ! Courant number diagnostics
+                         CALL dia_hth   ( kstp,      Nnn )      ! Thermocline depth (20 degres isotherm depth)
+      IF( ln_diadct  )   CALL dia_dct   ( kstp,      Nnn )      ! Transports
+                         CALL dia_ar5   ( kstp,      Nnn )      ! ar5 diag
+                         CALL dia_ptr   ( kstp,      Nnn )      ! Poleward adv/ldf TRansports diagnostics
+                         CALL dia_wri   ( kstp,      Nnn )      ! ocean model: outputs
+      IF( ln_crs     )   CALL crs_fld   ( kstp,      Nnn )      ! ocean model: online field coarsening & output
+      IF( lk_diadetide ) CALL dia_detide( kstp )                ! Weights computation for daily detiding of model diagnostics
+      IF( lk_diamlr  )   CALL dia_mlr                           ! Update time used in multiple-linear-regression analysis
+
+#if defined key_top
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Passive Tracer Model
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         CALL trc_stp       ( kstp, Nbb, Nnn, Nrhs, Naa )  ! time-stepping
+#endif
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Active tracers
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         ts(:,:,:,:,Nrhs) = 0._wp         ! set tracer trends to zero
+
+      IF(  lk_asminc .AND. ln_asmiau .AND. &
+         & ln_trainc )   CALL tra_asm_inc( kstp, Nbb, Nnn, ts, Nrhs )  ! apply tracer assimilation increment
+                         CALL tra_sbc    ( kstp,      Nnn, ts, Nrhs )  ! surface boundary condition
+      IF( ln_traqsr  )   CALL tra_qsr    ( kstp,      Nnn, ts, Nrhs )  ! penetrative solar radiation qsr
+      IF( ln_isf     )   CALL tra_isf    ( kstp,      Nnn, ts, Nrhs )  ! ice shelf heat flux
+      IF( ln_trabbc  )   CALL tra_bbc    ( kstp,      Nnn, ts, Nrhs )  ! bottom heat flux
+      IF( ln_trabbl  )   CALL tra_bbl    ( kstp, Nbb, Nnn, ts, Nrhs )  ! advective (and/or diffusive) bottom boundary layer scheme
+      IF( ln_tradmp  )   CALL tra_dmp    ( kstp, Nbb, Nnn, ts, Nrhs )  ! internal damping trends
+      IF( ln_bdy     )   CALL bdy_tra_dmp( kstp, Nbb,      ts, Nrhs )  ! bdy damping trends
+#if defined key_agrif
+      IF(.NOT. Agrif_Root())  &
+               &         CALL Agrif_Sponge_tra        ! tracers sponge
+#endif
+                         CALL tra_adv    ( kstp, Nbb, Nnn, ts, Nrhs )  ! hor. + vert. advection	==> RHS
+      IF( ln_zdfosm  )   CALL tra_osm    ( kstp,      Nnn, ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS
+      IF( lrst_oce .AND. ln_zdfosm ) &
+           &             CALL osm_rst    ( kstp,      Nnn, 'WRITE'  )  ! write OSMOSIS outputs + ww (so must do here) to restarts
+                         CALL tra_ldf    ( kstp, Nbb, Nnn, ts, Nrhs )  ! lateral mixing
+
+                         CALL tra_zdf    ( kstp, Nbb, Nnn, Nrhs, ts, Naa  )  ! vertical mixing and after tracer fields
+      IF( ln_zdfnpc  )   CALL tra_npc    ( kstp,      Nnn, Nrhs, ts, Naa  )  ! update after fields by non-penetrative convection
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Set boundary conditions, time filter and swap time levels
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+!!jc1: For agrif, it would be much better to finalize tracers/momentum here (e.g. bdy conditions) and move the swap
+!!    (and time filtering) after Agrif update. Then restart would be done after and would contain updated fields.
+!!    If so:
+!!    (i) no need to call agrif update at initialization time
+!!    (ii) no need to update "before" fields
+!!
+!!    Apart from creating new tra_swp/dyn_swp routines, this however:
+!!    (i) makes boundary conditions at initialization time computed from updated fields which is not the case between
+!!    two restarts => restartability issue. One can circumvent this, maybe, by assuming "interface separation",
+!!    e.g. a shift of the feedback interface inside child domain.
+!!    (ii) requires that all restart outputs of updated variables by agrif (e.g. passive tracers/tke/barotropic arrays) are done at the same
+!!    place.
+!!
+!!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine
+                         CALL zdyn_ts       ( Nnn, Naa, e3u, e3v, uu, vv )                   ! barotrope ajustment
+                         CALL finalize_sbc  ( kstp, Nbb, Naa, uu, vv, ts )                   ! boundary condifions
+                         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_qco   ( kstp, Nbb, Nnn, Naa, ts )                      ! time filtering of "now" tracer arrays
+                         CALL dyn_atf_qco   ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v  )  ! time filtering of "now" velocities and scale factors
+                         r3t(:,:,Nnn) = r3t_f(:,:)
+                         r3u(:,:,Nnn) = r3u_f(:,:)
+                         r3v(:,:,Nnn) = r3v_f(:,:)
+
+      !
+      ! Swap time levels
+      Nrhs = Nbb
+      Nbb = Nnn
+      Nnn = Naa
+      Naa = Nrhs
+      !
+      IF(.NOT.ln_linssh) CALL dom_qe_sf_update( kstp, Nbb, Nnn, Naa )  ! recompute vertical scale factors
+      !
+      IF( ln_diahsb  )   CALL dia_hsb       ( kstp, Nbb, Nnn )  ! - ML - global conservation diagnostics
+
+!!gm : This does not only concern the dynamics ==>>> add a new title
+!!gm2: why ouput restart before AGRIF update?
+!!
+!!jc: That would be better, but see comment above
+!!
+      IF( lrst_oce   )   CALL rst_write    ( kstp, Nbb, Nnn )   ! write output ocean restart file
+      IF( ln_sto_eos )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters
+
+#if defined key_agrif
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! AGRIF
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs      ! agrif_oce module copies of time level indices
+                         CALL Agrif_Integrate_ChildGrids( stplf )       ! allows to finish all the Child Grids before updating
+
+                         IF( Agrif_NbStepint() == 0 ) THEN
+                            CALL Agrif_update_all( )                  ! Update all components
+                         ENDIF
+#endif
+      IF( ln_diaobs  )   CALL dia_obs      ( kstp, Nnn )      ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Control
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         CALL stp_ctl      ( kstp, Nbb, Nnn, indic )
+
+      IF( kstp == nit000 ) THEN                          ! 1st time step only
+                                        CALL iom_close( numror )   ! close input  ocean restart file
+         IF(lwm)                        CALL FLUSH    ( numond )   ! flush output namelist oce
+         IF(lwm .AND. numoni /= -1 )    CALL FLUSH    ( numoni )   ! flush output namelist ice (if exist)
+      ENDIF
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Coupled mode
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+!!gm why lk_oasis and not lk_cpl ????
+      IF( lk_oasis   )   CALL sbc_cpl_snd( kstp, Nbb, Nnn )        ! coupled mode : field exchanges
+      !
+#if defined key_iomput
+      IF( kstp == nitend .OR. indic < 0 ) THEN
+                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
+                      IF(lrxios) CALL iom_context_finalize(      crxios_context          )
+         IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !
+      ENDIF
+#endif
+      !
+      IF( ln_timing )   CALL timing_stop('stplf')
+      !
+   END SUBROUTINE stplf
+
+
+   SUBROUTINE zdyn_ts (Kmm, Kaa, pe3u, pe3v, puu, pvv)
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdyn_ts  ***
+      !!
+      !! ** Purpose :   Finalize after horizontal velocity.
+      !!
+      !! ** Method  : * Ensure after velocities transport matches time splitting
+      !!             estimate (ln_dynspg_ts=T)
+      !!
+      !! ** Action :   puu(Kmm),pvv(Kmm),puu(Kaa),pvv(Kaa)   now and after horizontal velocity
+      !!----------------------------------------------------------------------
+      INTEGER                             , INTENT(in   ) :: Kmm, Kaa    ! before and after time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv         ! velocities
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in   ) :: pe3u, pe3v   ! scale factors
+      !
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
+      !
+      INTEGER  ::   jk   ! dummy loop indices
+      !!----------------------------------------------------------------------
+
+      IF ( ln_dynspg_ts ) THEN
+         ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
+         ! 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
+         !
+         DEALLOCATE( zue, zve )
+         !
+      ENDIF
+      !
+   END SUBROUTINE zdyn_ts
+
+
+   SUBROUTINE finalize_sbc (kt, Kbb, Kaa, puu, pvv, pts)
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE finalize_sbc  ***
+      !!
+      !! ** Purpose :   Apply the boundary condition on the after velocity
+      !!
+      !! ** Method  : * 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)
+      !!
+      !! ** Action :   puu(Kaa),pvv(Kaa)   after horizontal velocity and tracers
+      !!----------------------------------------------------------------------
+      INTEGER                             , INTENT(in   ) :: kt               ! ocean time-step index
+      INTEGER                             , INTENT(in   ) :: Kbb, 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,jpts,jpt), INTENT(inout) :: pts            ! active tracers
+      !
+      ! Update after tracer and velocity on domain lateral boundaries
+      !
+#if defined key_agrif
+            CALL Agrif_tra                     ! AGRIF zoom boundaries
+            CALL Agrif_dyn( kt )               !* AGRIF zoom boundaries
+#endif
+      !                                        ! local domain boundaries  (T-point, unchanged sign)
+      CALL lbc_lnk_multi( 'finalize_sbc', puu(:,:,:,       Kaa), 'U', -1., pvv(:,:,:       ,Kaa), 'V', -1.   &
+                       &                , pts(:,:,:,jp_tem,Kaa), 'T',  1., pts(:,:,:,jp_sal,Kaa), 'T',  1. )     !* local domain boundaries
+      !
+      !                                        !* BDY open boundaries
+      IF( ln_bdy )   THEN
+                               CALL bdy_tra( kt, Kbb, pts,      Kaa )
+         IF( ln_dynspg_exp )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa )
+         IF( ln_dynspg_ts  )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only=.true. )
+      ENDIF
+      !
+   END SUBROUTINE finalize_sbc
+
+   !
+   !!======================================================================
+END MODULE steplf
Index: MO/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 12623)
+++ 	(revision )
@@ -1,458 +1,0 @@
-MODULE steplf
-   !!======================================================================
-   !!                       ***  MODULE step  ***
-   !! Time-stepping   : manager of the ocean, tracer and ice time stepping
-   !!======================================================================
-   !! History :  OPA  !  1991-03  (G. Madec)  Original code
-   !!             -   !  1991-11  (G. Madec)
-   !!             -   !  1992-06  (M. Imbard)  add a first output record
-   !!             -   !  1996-04  (G. Madec)  introduction of dynspg
-   !!             -   !  1996-04  (M.A. Foujols)  introduction of passive tracer
-   !!            8.0  !  1997-06  (G. Madec)  new architecture of call
-   !!            8.2  !  1997-06  (G. Madec, M. Imbard, G. Roullet)  free surface
-   !!             -   !  1999-02  (G. Madec, N. Grima)  hpg implicit
-   !!             -   !  2000-07  (J-M Molines, M. Imbard)  Open Bondary Conditions
-   !!   NEMO     1.0  !  2002-06  (G. Madec)  free form, suppress macro-tasking
-   !!             -   !  2004-08  (C. Talandier) New trends organization
-   !!             -   !  2005-01  (C. Ethe) Add the KPP closure scheme
-   !!             -   !  2005-11  (G. Madec)  Reorganisation of tra and dyn calls
-   !!             -   !  2006-01  (L. Debreu, C. Mazauric)  Agrif implementation
-   !!             -   !  2006-07  (S. Masson)  restart using iom
-   !!            3.2  !  2009-02  (G. Madec, R. Benshila)  reintroduicing z*-coordinate
-   !!             -   !  2009-06  (S. Masson, G. Madec)  TKE restart compatible with key_cpl
-   !!            3.3  !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
-   !!             -   !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA
-   !!            3.4  !  2011-04  (G. Madec, C. Ethe) Merge of dtatem and dtasal
-   !!            3.6  !  2012-07  (J. Simeon, G. Madec. C. Ethe)  Online coarsening of outputs
-   !!            3.6  !  2014-04  (F. Roquet, G. Madec) New equations of state
-   !!            3.6  !  2014-10  (E. Clementi, P. Oddo) Add Qiao vertical mixing in case of waves
-   !!            3.7  !  2014-10  (G. Madec)  LDF simplication
-   !!             -   !  2014-12  (G. Madec) remove KPP scheme
-   !!             -   !  2015-11  (J. Chanut) free surface simplification (remove filtered free surface)
-   !!            4.0  !  2017-05  (G. Madec)  introduction of the vertical physics manager (zdfphy)
-   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   stplf             : OPA system time-stepping
-   !!----------------------------------------------------------------------
-   USE step_oce         ! time stepping definition modules
-   !
-   USE iom              ! xIOs server
-   USE domqe
-   USE traatfLF          ! time filtering                   (tra_atf_lf routine)
-   USE dynatfLF          ! time filtering                   (dyn_atf_lf routine)
-   USE dynspg_ts         ! surface pressure gradient: split-explicit scheme (define un_adv)
-   USE bdydyn            ! ocean open boundary conditions (define bdy_dyn)
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   stplf   ! called by nemogcm.F90
-
-   !!----------------------------------------------------------------------
-   !! time level indices
-   !!----------------------------------------------------------------------
-   INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs          !! used by nemo_init
-
-   !!----------------------------------------------------------------------
-   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id: step.F90 12377 2020-02-12 14:39:06Z acc $
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-#if defined key_agrif
-   RECURSIVE SUBROUTINE stplf( )
-      INTEGER             ::   kstp   ! ocean time-step index
-#else
-   SUBROUTINE stplf( kstp )
-      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
-#endif
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE stplf  ***
-      !!
-      !! ** Purpose : - Time stepping of OPA  (momentum and active tracer eqs.)
-      !!              - Time stepping of SI3 (dynamic and thermodynamic eqs.)
-      !!              - Time stepping of TRC  (passive tracer eqs.)
-      !!
-      !! ** Method  : -1- Update forcings and data
-      !!              -2- Update ocean physics
-      !!              -3- Compute the t and s trends
-      !!              -4- Update t and s
-      !!              -5- Compute the momentum trends
-      !!              -6- Update the horizontal velocity
-      !!              -7- Compute the diagnostics variables (rd,N2, hdiv,w)
-      !!              -8- Outputs and diagnostics
-      !!----------------------------------------------------------------------
-      INTEGER ::   ji, jj, jk   ! dummy loop indice
-      INTEGER ::   indic        ! error indicator if < 0
-!!gm kcall can be removed, I guess
-      INTEGER ::   kcall        ! optional integer argument (dom_vvl_sf_nxt)
-      !! ---------------------------------------------------------------------
-#if defined key_agrif
-      kstp = nit000 + Agrif_Nb_Step()
-      Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
-      IF( lk_agrif_debug ) THEN
-         IF( Agrif_Root() .and. lwp)   WRITE(*,*) '---'
-         IF(lwp)   WRITE(*,*) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint()
-      ENDIF
-      IF( kstp == nit000 + 1 )   lk_agrif_fstep = .FALSE.
-# if defined key_iomput
-      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context )
-# endif
-#endif
-      !
-      IF( ln_timing )   CALL timing_start('stplf')
-      !
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! update I/O and calendar
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                             indic = 0                ! reset to no error condition
-
-      IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
-                             CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including passible AGRIF zoom)
-         IF( lk_diamlr   )   CALL dia_mlr_iom_init    ! with additional setup for multiple-linear-regression analysis
-                             CALL iom_init_closedef
-         IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid
-      ENDIF
-      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
-                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell IOM we are at time step kstp
-      IF( ln_crs         )   CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell IOM we are at time step kstp
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Update external forcing (tides, open boundaries, ice shelf interaction and surface boundary condition (including sea-ice)
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( ln_tide    )   CALL tide_update( kstp )                     ! update tide potential
-      IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                        ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)
-      IF( ln_bdy     )   CALL bdy_dta ( kstp, Nnn )                   ! update dynamic & tracer data at open boundaries
-      IF( ln_isf     )   CALL isf_stp ( kstp, Nnn )
-                         CALL sbc     ( kstp, Nbb, Nnn )              ! Sea Boundary Condition (including sea-ice)
-
-      !             !!st      !!!!!!!!!!!!!!!!!!!!!!!
-      !
-      ! emp = 0._wp
-      ! emp_b = 0._wp
-      ! qns = 0._wp
-      ! qsr = 0._wp
-      ! qns_b = 0._wp
-      !
-      !             !!st
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Update stochastic parameters and random T/S fluctuations
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      IF( ln_sto_eos ) CALL sto_par( kstp )          ! Stochastic parameters
-      IF( ln_sto_eos ) CALL sto_pts( ts(:,:,:,:,Nnn)  )          ! Random T/S fluctuations
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Ocean physics update
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !  THERMODYNAMICS
-                         CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nnn )       ! before local thermal/haline expension ratio at T-points
-                         CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn )       ! now    local thermal/haline expension ratio at T-points
-                         CALL bn2    ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
-                         CALL bn2    ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn  ) ! now    Brunt-Vaisala frequency
-
-      !  VERTICAL PHYSICS
-                         CALL zdf_phy( kstp, Nbb, Nnn, Nrhs )   ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
-
-      !  LATERAL  PHYSICS
-      !
-      IF( l_ldfslp ) THEN                             ! slope of lateral mixing
-                         CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )               ! before in situ density
-
-         IF( ln_zps .AND. .NOT. ln_isfcav)                                    &
-            &            CALL zps_hde    ( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
-            &                                          rhd, gru , grv    )       ! of t, s, rd at the last ocean level
-
-         IF( ln_zps .AND.       ln_isfcav)                                                &
-            &            CALL zps_hde_isf( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
-            &                                          rhd, gru , grv , grui, grvi   )       ! of t, s, rd at the first ocean level
-         IF( ln_traldf_triad ) THEN
-                         CALL ldf_slp_triad( kstp, Nbb, Nnn )             ! before slope for triad operator
-         ELSE
-                         CALL ldf_slp     ( kstp, rhd, rn2b, Nbb, Nnn )   ! before slope for standard operator
-         ENDIF
-      ENDIF
-      !                                                                        ! eddy diffusivity coeff.
-      IF( l_ldftra_time .OR. l_ldfeiv_time )   CALL ldf_tra( kstp, Nbb, Nnn )  !       and/or eiv coeff.
-      IF( l_ldfdyn_time                    )   CALL ldf_dyn( kstp, Nbb )       ! eddy viscosity coeff.
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      !  Ocean dynamics : hdiv, ssh, e3, u, v, w
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
-                            CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor)
-      IF( .NOT.ln_linssh )  CALL dom_qe_r3c    ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa) )
-      IF( .NOT.ln_linssh )  CALL dom_qe_sf_nxt ( kstp, Nbb, Nnn,      Naa )    ! after vertical scale factors
-                            CALL wzv           ( kstp, Nbb, Nnn, ww,  Naa )    ! now cross-level velocity
-      IF( ln_zad_Aimp )     CALL wAimp         ( kstp,      Nnn           )  ! Adaptive-implicit vertical advection partitioning
-                            CALL eos    ( ts(:,:,:,:,Nnn), rhd, rhop, gdept(:,:,:,Nnn) )  ! now in situ density for hpg computation
-
-
-                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
-                         vv(:,:,:,Nrhs) = 0._wp
-
-      IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   &
-               &         CALL dyn_asm_inc   ( kstp, Nbb, Nnn, uu, vv, Nrhs )  ! apply dynamics assimilation increment
-      IF( ln_bdy     )   CALL bdy_dyn3d_dmp ( kstp, Nbb,      uu, vv, Nrhs )  ! bdy damping trends
-#if defined key_agrif
-      IF(.NOT. Agrif_Root())  &
-               &         CALL Agrif_Sponge_dyn        ! momentum sponge
-#endif
-                         CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! 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
-      IF( ln_zdfosm  )   CALL dyn_osm( kstp,      Nnn      , uu, vv, Nrhs )  ! OSMOSIS non-local velocity fluxes ==> RHS
-                         CALL dyn_hpg( kstp,      Nnn      , uu, vv, Nrhs )  ! horizontal gradient of Hydrostatic pressure
-                         CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa )  ! surface pressure gradient
-
-                                                      ! With split-explicit free surface, since now transports have been updated and ssh(:,:,Nrhs) as well
-      IF( ln_dynspg_ts ) THEN                         ! vertical scale factors and vertical velocity need to be updated
-                            CALL div_hor       ( kstp, Nbb, Nnn )                ! Horizontal divergence  (2nd call in time-split case)
-         IF(.NOT.ln_linssh) CALL dom_qe_r3c    ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )
-         IF(.NOT.ln_linssh) CALL dom_qe_sf_nxt ( kstp, Nbb, Nnn, Naa, kcall=2 )  ! after vertical scale factors (update depth average component)
-      ENDIF
-                            CALL dyn_zdf    ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion
-      IF( ln_dynspg_ts ) THEN                                                       ! vertical scale factors and vertical velocity need to be updated
-                            CALL wzv        ( kstp, Nbb, Nnn, ww, Naa )             ! now cross-level velocity
-         IF( ln_zad_Aimp )  CALL wAimp      ( kstp,      Nnn )                      ! Adaptive-implicit vertical advection partitioning
-      ENDIF
-
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! cool skin
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF ( ln_diurnal )  CALL diurnal_layers( kstp )
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! diagnostics and outputs
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( ln_floats  )   CALL flo_stp   ( kstp, Nbb, Nnn )      ! drifting Floats
-      IF( ln_diacfl  )   CALL dia_cfl   ( kstp,      Nnn )      ! Courant number diagnostics
-                         CALL dia_hth   ( kstp,      Nnn )      ! Thermocline depth (20 degres isotherm depth)
-      IF( ln_diadct  )   CALL dia_dct   ( kstp,      Nnn )      ! Transports
-                         CALL dia_ar5   ( kstp,      Nnn )      ! ar5 diag
-                         CALL dia_ptr   ( kstp,      Nnn )      ! Poleward adv/ldf TRansports diagnostics
-                         CALL dia_wri   ( kstp,      Nnn )      ! ocean model: outputs
-      IF( ln_crs     )   CALL crs_fld   ( kstp,      Nnn )      ! ocean model: online field coarsening & output
-      IF( lk_diadetide ) CALL dia_detide( kstp )                ! Weights computation for daily detiding of model diagnostics
-      IF( lk_diamlr  )   CALL dia_mlr                           ! Update time used in multiple-linear-regression analysis
-
-#if defined key_top
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Passive Tracer Model
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL trc_stp       ( kstp, Nbb, Nnn, Nrhs, Naa )  ! time-stepping
-#endif
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Active tracers
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         ts(:,:,:,:,Nrhs) = 0._wp         ! set tracer trends to zero
-
-      IF(  lk_asminc .AND. ln_asmiau .AND. &
-         & ln_trainc )   CALL tra_asm_inc( kstp, Nbb, Nnn, ts, Nrhs )  ! apply tracer assimilation increment
-                         CALL tra_sbc    ( kstp,      Nnn, ts, Nrhs )  ! surface boundary condition
-      IF( ln_traqsr  )   CALL tra_qsr    ( kstp,      Nnn, ts, Nrhs )  ! penetrative solar radiation qsr
-      IF( ln_isf     )   CALL tra_isf    ( kstp,      Nnn, ts, Nrhs )  ! ice shelf heat flux
-      IF( ln_trabbc  )   CALL tra_bbc    ( kstp,      Nnn, ts, Nrhs )  ! bottom heat flux
-      IF( ln_trabbl  )   CALL tra_bbl    ( kstp, Nbb, Nnn, ts, Nrhs )  ! advective (and/or diffusive) bottom boundary layer scheme
-      IF( ln_tradmp  )   CALL tra_dmp    ( kstp, Nbb, Nnn, ts, Nrhs )  ! internal damping trends
-      IF( ln_bdy     )   CALL bdy_tra_dmp( kstp, Nbb,      ts, Nrhs )  ! bdy damping trends
-#if defined key_agrif
-      IF(.NOT. Agrif_Root())  &
-               &         CALL Agrif_Sponge_tra        ! tracers sponge
-#endif
-                         CALL tra_adv    ( kstp, Nbb, Nnn, ts, Nrhs )  ! hor. + vert. advection	==> RHS
-      IF( ln_zdfosm  )   CALL tra_osm    ( kstp,      Nnn, ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS
-      IF( lrst_oce .AND. ln_zdfosm ) &
-           &             CALL osm_rst    ( kstp,      Nnn, 'WRITE'  )  ! write OSMOSIS outputs + ww (so must do here) to restarts
-                         CALL tra_ldf    ( kstp, Nbb, Nnn, ts, Nrhs )  ! lateral mixing
-
-                         CALL tra_zdf    ( kstp, Nbb, Nnn, Nrhs, ts, Naa  )  ! vertical mixing and after tracer fields
-      IF( ln_zdfnpc  )   CALL tra_npc    ( kstp,      Nnn, Nrhs, ts, Naa  )  ! update after fields by non-penetrative convection
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Set boundary conditions, time filter and swap time levels
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-!!jc1: For agrif, it would be much better to finalize tracers/momentum here (e.g. bdy conditions) and move the swap
-!!    (and time filtering) after Agrif update. Then restart would be done after and would contain updated fields.
-!!    If so:
-!!    (i) no need to call agrif update at initialization time
-!!    (ii) no need to update "before" fields
-!!
-!!    Apart from creating new tra_swp/dyn_swp routines, this however:
-!!    (i) makes boundary conditions at initialization time computed from updated fields which is not the case between
-!!    two restarts => restartability issue. One can circumvent this, maybe, by assuming "interface separation",
-!!    e.g. a shift of the feedback interface inside child domain.
-!!    (ii) requires that all restart outputs of updated variables by agrif (e.g. passive tracers/tke/barotropic arrays) are done at the same
-!!    place.
-!!
-!!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine
-                         CALL zdyn_ts       ( Nnn, Naa, e3u, e3v, uu, vv )                   ! barotrope ajustment
-                         CALL finalize_sbc  ( kstp, Nbb, Naa, uu, vv, ts )                   ! boundary condifions
-                         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
-                         r3t(:,:,Nnn) = r3t_f(:,:)
-                         r3u(:,:,Nnn) = r3u_f(:,:)
-                         r3v(:,:,Nnn) = r3v_f(:,:)
-
-      !
-      ! Swap time levels
-      Nrhs = Nbb
-      Nbb = Nnn
-      Nnn = Naa
-      Naa = Nrhs
-      !
-      IF(.NOT.ln_linssh) CALL dom_qe_sf_update( kstp, Nbb, Nnn, Naa )  ! recompute vertical scale factors
-      !
-      IF( ln_diahsb  )   CALL dia_hsb       ( kstp, Nbb, Nnn )  ! - ML - global conservation diagnostics
-
-!!gm : This does not only concern the dynamics ==>>> add a new title
-!!gm2: why ouput restart before AGRIF update?
-!!
-!!jc: That would be better, but see comment above
-!!
-      IF( lrst_oce   )   CALL rst_write    ( kstp, Nbb, Nnn )   ! write output ocean restart file
-      IF( ln_sto_eos )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters
-
-#if defined key_agrif
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! AGRIF
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs      ! agrif_oce module copies of time level indices
-                         CALL Agrif_Integrate_ChildGrids( stplf )       ! allows to finish all the Child Grids before updating
-
-                         IF( Agrif_NbStepint() == 0 ) THEN
-                            CALL Agrif_update_all( )                  ! Update all components
-                         ENDIF
-#endif
-      IF( ln_diaobs  )   CALL dia_obs      ( kstp, Nnn )      ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Control
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL stp_ctl      ( kstp, Nbb, Nnn, indic )
-
-      IF( kstp == nit000 ) THEN                          ! 1st time step only
-                                        CALL iom_close( numror )   ! close input  ocean restart file
-         IF(lwm)                        CALL FLUSH    ( numond )   ! flush output namelist oce
-         IF(lwm .AND. numoni /= -1 )    CALL FLUSH    ( numoni )   ! flush output namelist ice (if exist)
-      ENDIF
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Coupled mode
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-!!gm why lk_oasis and not lk_cpl ????
-      IF( lk_oasis   )   CALL sbc_cpl_snd( kstp, Nbb, Nnn )        ! coupled mode : field exchanges
-      !
-#if defined key_iomput
-      IF( kstp == nitend .OR. indic < 0 ) THEN
-                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
-                      IF(lrxios) CALL iom_context_finalize(      crxios_context          )
-         IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !
-      ENDIF
-#endif
-      !
-      IF( ln_timing )   CALL timing_stop('stplf')
-      !
-   END SUBROUTINE stplf
-
-
-   SUBROUTINE zdyn_ts (Kmm, Kaa, pe3u, pe3v, puu, pvv)
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE zdyn_ts  ***
-      !!
-      !! ** Purpose :   Finalize after horizontal velocity.
-      !!
-      !! ** Method  : * Ensure after velocities transport matches time splitting
-      !!             estimate (ln_dynspg_ts=T)
-      !!
-      !! ** Action :   puu(Kmm),pvv(Kmm),puu(Kaa),pvv(Kaa)   now and after horizontal velocity
-      !!----------------------------------------------------------------------
-      INTEGER                             , INTENT(in   ) :: Kmm, Kaa    ! before and after time level indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv         ! velocities
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in   ) :: pe3u, pe3v   ! scale factors
-      !
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
-      !
-      INTEGER  ::   jk   ! dummy loop indices
-      !!----------------------------------------------------------------------
-
-      IF ( ln_dynspg_ts ) THEN
-         ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
-         ! 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
-         !
-         DEALLOCATE( zue, zve )
-         !
-      ENDIF
-      !
-   END SUBROUTINE zdyn_ts
-
-
-   SUBROUTINE finalize_sbc (kt, Kbb, Kaa, puu, pvv, pts)
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE finalize_sbc  ***
-      !!
-      !! ** Purpose :   Apply the boundary condition on the after velocity
-      !!
-      !! ** Method  : * 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)
-      !!
-      !! ** Action :   puu(Kaa),pvv(Kaa)   after horizontal velocity and tracers
-      !!----------------------------------------------------------------------
-      INTEGER                             , INTENT(in   ) :: kt               ! ocean time-step index
-      INTEGER                             , INTENT(in   ) :: Kbb, 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,jpts,jpt), INTENT(inout) :: pts            ! active tracers
-      !
-      ! Update after tracer and velocity on domain lateral boundaries
-      !
-#if defined key_agrif
-            CALL Agrif_tra                     ! AGRIF zoom boundaries
-            CALL Agrif_dyn( kt )               !* AGRIF zoom boundaries
-#endif
-      !                                        ! local domain boundaries  (T-point, unchanged sign)
-      CALL lbc_lnk_multi( 'finalize_sbc', puu(:,:,:,       Kaa), 'U', -1., pvv(:,:,:       ,Kaa), 'V', -1.   &
-                       &                , pts(:,:,:,jp_tem,Kaa), 'T',  1., pts(:,:,:,jp_sal,Kaa), 'T',  1. )     !* local domain boundaries
-      !
-      !                                        !* BDY open boundaries
-      IF( ln_bdy )   THEN
-                               CALL bdy_tra( kt, Kbb, pts,      Kaa )
-         IF( ln_dynspg_exp )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa )
-         IF( ln_dynspg_ts  )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only=.true. )
-      ENDIF
-      !
-   END SUBROUTINE finalize_sbc
-
-   !
-   !!======================================================================
-END MODULE steplf
