Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/stpMLF.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/stepLF.F90	(revision 13426)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/stpMLF.F90	(revision 13603)
@@ -1,29 +1,29 @@
-MODULE stpLF
+MODULE stpMLF
    !!======================================================================
-   !!                       ***  MODULE step  ***
+   !!                       ***  MODULE stpMLF  ***
    !! Time-stepping   : manager of the shallow water equation time stepping
+   !!                   Modified Leap Frog scheme
    !!======================================================================
    !! History :  NEMO !  2020-03  (A. Nasser, G. Madec)  Original code from  4.0.2
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   stp             : Shallow Water time-stepping
-   !!----------------------------------------------------------------------
-   USE step_oce         ! time stepping definition modules
-   USE phycst           ! physical constants
-   USE usrdef_nam
+   !!             -   !  2020-10  (S. Techene, G. Madec)  cleanning
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   stp_MLF       : MLF Shallow Water Eq. time-stepping
+   !!----------------------------------------------------------------------
+   USE stp_oce        ! modules used in nemo_init and stp_MLF
    !
-   USE iom              ! xIOs server 
-   USE domqco
-
+   USE domqco         ! quasi-eulerian coordinate
+   USE phycst         ! physical constants
+   USE usrdef_nam     ! user defined namelist parameters
+!!st   USE usrdef_sbc     ! user defined surface boundary cond
+   
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   stp_LF   ! called by nemogcm.F90
+   PUBLIC   stp_MLF   ! called by nemogcm.F90
    
-   !!----------------------------------------------------------------------
-   !! time level indices
-   !!----------------------------------------------------------------------
-   INTEGER, PUBLIC ::   Nbb, Nnn, Naa, Nrhs   !! used by nemo_init
+   !                                          !**  time level indices  **!
+   INTEGER, PUBLIC ::   Nbb, Nnn, Naa, Nrhs   !: used by nemo_init
       
    !! * Substitutions
@@ -37,13 +37,7 @@
 CONTAINS
 
-#if defined key_agrif
-   RECURSIVE SUBROUTINE stp_LF( )
-      INTEGER             ::   kstp   ! ocean time-step index
-#else
-   SUBROUTINE stp_LF( kstp )
-      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
-#endif
+   SUBROUTINE stp_MLF( kstp )
       !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE stp  ***
+      !!                     ***  ROUTINE stp_MLF  ***
       !!
       !! ** Purpose : - Time stepping of shallow water (SHW) (momentum and ssh eqs.)
@@ -56,28 +50,14 @@
       !!              -6- 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)
-      REAL(wp)::   z1_2rho0     ! local scalars
-      
-      REAL(wp) ::   zue3a, zue3n, zue3b    ! local scalars
-      REAL(wp) ::   zve3a, zve3n, zve3b    !   -      -
-      REAL(wp) ::   ze3t_tf, ze3u_tf, ze3v_tf, zua, zva
+      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
+      !
+      INTEGER ::   ji, jj, jk             ! dummy loop indices
+      INTEGER ::   indic                  ! error indicator if < 0
+      REAL(wp)::   z1_2rho0               ! local scalars
+      REAL(wp)::   zrhs_u, zue3a, zue3n, zue3b, zua   ! local scalars
+      REAL(wp)::   zrhs_v, zve3a, zve3n, zve3b, zva   !   -      -
       !! ---------------------------------------------------------------------
-#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('stp_LF')
+      !
+      IF( ln_timing )   CALL timing_start('stp_MLF')
       !
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -85,10 +65,10 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !
-      IF( l_1st_euler ) THEN     ! start or restart with Euler 1st time-step
-         rDt   =  rn_Dt   
+      IF( l_1st_euler ) THEN   ! start or restart with Euler 1st time-step
+         rDt   = rn_Dt   
          r1_Dt = 1._wp / rDt
       ENDIF
 
-      IF ( kstp == nit000 )   ww(:,:,:) = 0._wp   ! initialize vertical velocity one for all to zero
+      IF( kstp == nit000 )   ww(:,:,:) = 0._wp   ! initialize vertical velocity one for all to zero
 
       !
@@ -97,89 +77,75 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
                              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
-                         CALL sbc     ( kstp, Nbb, Nnn )              ! Sea Boundary Condition
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Ocean physics update
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !  LATERAL  PHYSICS
-      !                                                                        ! eddy diffusivity coeff.
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Update external forcing   (SWE: surface boundary condition only)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+                             CALL sbc     ( kstp, Nbb, Nnn )                   ! Sea Boundary Condition
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Ocean physics update   (SWE: eddy viscosity only)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
       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)
+      !  RHS of horizontal velocity
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
                          uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
                          vv(:,:,:,Nrhs) = 0._wp
 
-      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( .NOT.ln_linssh )  CALL dom_qco_r3c   ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! "after" ssh./h._0 ratio explicit
-      !IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt_st( kstp, Nbb, Nnn,      Naa )    ! after vertical scale factors 
-!!an - calcul du gradient de pression horizontal (explicit)
+                         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                     ==> RHS
+
+      z1_2rho0 = 0.5_wp * r1_rho0
+      !
       DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
-         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
+         !                                                              ! horizontal pressure gradient
+         zrhs_u =        - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
+         zrhs_v =        - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
+         !                                                              ! wind stress and layer friction
+         zrhs_u = zrhs_u + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn)   &
+            &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
+         zrhs_v = zrhs_v + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   &
+            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
+         !                                                              ! ==> RHS
+         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + zrhs_u
+         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + zrhs_v
       END_3D
-      !
-      ! add wind stress forcing and layer linear friction to the RHS 
-      z1_2rho0 = 0.5_wp * r1_rho0
-      DO_3D( 0, 0, 0, 0,1,jpkm1)
-         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn)   &
-            &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
-         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   &
-            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
-      END_3D   
-!!an         
-  
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Leap-Frog time splitting + Robert-Asselin time filter on u,v,e3 
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-!!st ssh_atf : add ssh filtering up there
-      IF ( .NOT.( l_1st_euler ) ) THEN   ! Only do time filtering for leapfrog timesteps
-         !                                                  ! filtering "now" field
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Time stepping of ssh Eq. (and update r3_Naa)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !        ! Leap Frog time stepping ==> ssh_Naa and r3_Naa
+                         CALL ssh_nxt( kstp, Nbb, Nnn, ssh   , Naa )    ! after ssh 
+      !                                                                 ! after ssh/h_0 ratio explicit
+                         CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )
+      !        ! Asselin filter ==> ssh_Nnn filtered
+      IF ( .NOT.( l_1st_euler ) ) THEN   ! Time filtering of now ssh
          ssh(:,:,Nnn) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2._wp * ssh(:,:,Nnn) + ssh(:,:,Naa) )
       ENDIF
-!!st ssh_atf
-          
-!! what about  IF( .NOT.ln_linssh )  ?
-!!an futur module dyn_nxt (a la place de dyn_atf)
-      
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Time stepping of dynamics (u,v)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
       IF( ln_dynadv_vec ) THEN      ! vector invariant form : applied on velocity
-         IF( l_1st_euler ) THEN           ! Euler time stepping (no Asselin filter)
-            DO_3D( 0, 0, 0, 0,1,jpkm1)
+         IF( l_1st_euler ) THEN          ! Euler time stepping (no Asselin filter)
+            DO_3D( 0,0, 0,0, 1,jpkm1)
                uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
                vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
              END_3D          
-         ELSE                             ! Leap Frog time stepping + Asselin filter         
-            DO_3D( 1, 1, 1, 1,1,jpkm1)
+         ELSE                            ! Leap Frog time stepping + Asselin filter         
+            DO_3D( 0,0, 0,0, 1,jpkm1)
                zua = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
                zva = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
@@ -191,22 +157,16 @@
                vv(ji,jj,jk,Naa) = zva
             END_3D
-            CALL dom_qco_r3c( ssh(:,:,Nnn), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn) )   ! "now" ssh/h_0 ratio from filtrered ssh
-#if ! defined key_qco
-            DO jk = 1, jpkm1
-               e3t(:,:,jk,Nnn) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Nnn) )
-               e3u(:,:,jk,Nnn) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Nnn) )
-               e3v(:,:,jk,Nnn) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Nnn) )
-            END DO
-#endif
+            !                            ! Update r3_Nnn
+            CALL dom_qco_r3c( ssh(:,:,Nnn), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn) )   ! now ssh/h_0 ratio from filtrered ssh
          ENDIF
          !
       ELSE                          ! flux form : applied on thickness weighted velocity
-         IF( l_1st_euler ) THEN           ! Euler time stepping (no Asselin filter)
-            DO_3D( 0, 0, 0, 0,1,jpkm1)
+         IF( l_1st_euler ) THEN          ! Euler time stepping (no Asselin filter)
+            DO_3D( 0,0, 0,0, 1,jpkm1)
                zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
                zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
-               !                                                ! LF time stepping
-               zue3a = zue3b + rDt * e3u(ji,jj,jk,Nrhs) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
-               zve3a = zve3b + rDt * e3v(ji,jj,jk,Nrhs) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
+               !                                                ! Euler time stepping
+               zue3a = zue3b + rDt * e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
+               zve3a = zve3b + rDt * e3v(ji,jj,jk,Nnn) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
                !
                uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)    
@@ -214,59 +174,33 @@
             END_3D
          ELSE                             ! Leap Frog time stepping + Asselin filter
-            CALL dom_qco_r3c( ssh(:,:,Nnn), r3t_f(:,:), r3u_f(:,:), r3v_f(:,:) )   ! "now" ssh/h_0 ratio from filtrered ssh
-            DO_3D( 1, 1, 1, 1,1,jpkm1)
-               zue3n = e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nnn)
-               zve3n = e3v(ji,jj,jk,Nnn) * vv(ji,jj,jk,Nnn)
-               zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
-               zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
+            CALL dom_qco_r3c( ssh(:,:,Nnn), r3t_f(:,:), r3u_f(:,:), r3v_f(:,:) )   ! now ssh/h_0 ratio from filtrered ssh
+            DO_3D( 0,0, 0,0, 1,jpkm1)
+               zue3n = ( 1._wp + r3u(ji,jj,Nnn) ) * uu(ji,jj,jk,Nnn)
+               zve3n = ( 1._wp + r3v(ji,jj,Nnn) ) * vv(ji,jj,jk,Nnn)
+               zue3b = ( 1._wp + r3u(ji,jj,Nbb) ) * uu(ji,jj,jk,Nbb)
+               zve3b = ( 1._wp + r3v(ji,jj,Nbb) ) * vv(ji,jj,jk,Nbb)
                !                                                ! LF time stepping
-               zue3a = zue3b + rDt * e3u(ji,jj,jk,Nrhs) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
-               zve3a = zve3b + rDt * e3v(ji,jj,jk,Nrhs) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
-               !                                                ! Asselin time filter on e3u/v/t
-               ze3u_tf = e3u(ji,jj,jk,Nnn) + rn_atfp * ( e3u(ji,jj,jk,Nbb) - 2._wp * e3u(ji,jj,jk,Nnn)  + e3u(ji,jj,jk,Naa) ) 
-               ze3v_tf = e3v(ji,jj,jk,Nnn) + rn_atfp * ( e3v(ji,jj,jk,Nbb) - 2._wp * e3v(ji,jj,jk,Nnn)  + e3v(ji,jj,jk,Naa) )
-               ze3t_tf = e3t(ji,jj,jk,Nnn) + rn_atfp * ( e3t(ji,jj,jk,Nbb) - 2._wp * e3t(ji,jj,jk,Nnn)  + e3t(ji,jj,jk,Naa) ) 
+               zue3a = zue3b + rDt * ( 1._wp + r3u(ji,jj,Nnn) ) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
+               zve3a = zve3b + rDt * ( 1._wp + r3v(ji,jj,Nnn) ) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
                !                                                ! Asselin time filter on u,v (Nnn)
-               uu(ji,jj,jk,Nnn) = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / (e3u_0(ji,jj,jk) * ( 1._wp + r3u_f(ji,jj) ))
-               vv(ji,jj,jk,Nnn) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / (e3v_0(ji,jj,jk) * ( 1._wp + r3v_f(ji,jj) ))
+               uu(ji,jj,jk,Nnn) = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ( 1._wp + r3u_f(ji,jj) )
+               vv(ji,jj,jk,Nnn) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ( 1._wp + r3v_f(ji,jj) )
                !
-               uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)    
-               vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)
+               uu(ji,jj,jk,Naa) = zue3a / ( 1._wp + r3u(ji,jj,Naa) )    
+               vv(ji,jj,jk,Naa) = zve3a / ( 1._wp + r3v(ji,jj,Naa) )
             END_3D
+            !                             ! Update r3_Nnn with time filtered values
             r3t(:,:,Nnn) = r3t_f(:,:)
             r3u(:,:,Nnn) = r3u_f(:,:)
             r3v(:,:,Nnn) = r3v_f(:,:)
-#if ! defined key_qco
-            DO jk = 1, jpkm1
-               e3t(:,:,jk,Nnn) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Nnn) )
-               e3u(:,:,jk,Nnn) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Nnn) )
-               e3v(:,:,jk,Nnn) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Nnn) )
-            END DO
-#endif         
          ENDIF
       ENDIF
 
-
-      CALL lbc_lnk_multi( 'stp', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.,   &   !* local domain boundaries
-         &                       uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.    )     
-
-!!an         
+      CALL lbc_lnk_multi( 'stp_MLF', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.,   &   !* local domain boundaries
+         &                           uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.    )     
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Set boundary conditions, time filter and swap time levels
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
-!!an TO BE ADDED : dyn_nxt
-!!                         CALL dyn_atf       ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v  )  ! time filtering of "now" velocities and scale factors
-!!an TO BE ADDED : a simplifier
-!                         CALL ssh_atf       ( kstp, Nbb, Nnn, Naa, ssh )                     ! time filtering of "now" sea surface height
-!!st copied above 
-!!      IF ( .NOT.( l_1st_euler ) ) THEN   ! Only do time filtering for leapfrog timesteps
-!!         !                                                  ! filtering "now" field
-!!         ssh(:,:,Nnn) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2 * ssh(:,:,Nnn) + ssh(:,:,Naa) )
-!!      ENDIF
-!!st      
-!!an 
-
 
       ! Swap time levels
@@ -275,37 +209,19 @@
       Nnn = Naa
       Naa = Nrhs
-      !
-!                         CALL dom_vvl_sf_update_st( kstp, Nbb, Nnn, Naa )  ! recompute vertical scale factors
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! 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_wri   ( kstp,      Nnn )      ! ocean model: outputs
-
-      !
-      IF( lrst_oce   )   CALL rst_write    ( kstp, Nbb, Nnn )   ! write output ocean restart file
-          
-
-#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( stp_LF )       ! 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)
+      
+      IF( ln_diacfl  )   CALL dia_cfl  ( kstp,      Nnn )      ! Courant number diagnostics
+                         CALL dia_wri  ( kstp,      Nnn )      ! ocean model: outputs
+      !
+      IF( lrst_oce   )   CALL rst_write( kstp, Nbb, Nnn )      ! write output ocean restart file
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Control
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL stp_ctl      ( kstp, Nbb, Nnn, indic )
-               
-                         
+                         CALL stp_ctl      ( kstp, Nnn )
+
       IF( kstp == nit000 ) THEN                          ! 1st time step only
                                         CALL iom_close( numror )   ! close input  ocean restart file
@@ -316,20 +232,22 @@
       !
 #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( kstp == nitend .OR. indic < 0 ) THEN
+!!st : cxios_context needed ? because opened earlier ???         
+         CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF
+!!st : crxios_context not needed associated with coarsening ! 
+         IF(lrxios) CALL iom_context_finalize( crxios_context )
       ENDIF
 #endif
       !
       IF( l_1st_euler ) THEN         ! recover Leap-frog timestep
-         rDt = 2._wp * rn_Dt   
+         rDt   = 2._wp * rn_Dt   
          r1_Dt = 1._wp / rDt
          l_1st_euler = .FALSE.      
       ENDIF
       !
-      IF( ln_timing )   CALL timing_stop('stp')
-      !
-   END SUBROUTINE stp_LF
-   !
+      IF( ln_timing )   CALL timing_stop('stp_MLF')
+      !
+   END SUBROUTINE stp_MLF
+
    !!======================================================================
-END MODULE stpLF
+END MODULE stpMLF
