Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/domzgr_substitute.h90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/domzgr_substitute.h90	(revision 13707)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/domzgr_substitute.h90	(revision 13708)
@@ -16,4 +16,5 @@
 #   define  e3v(i,j,k,t)   (e3v_0(i,j,k)*(1._wp+r3v(i,j,t)))
 #   define  e3f(i,j,k)     (e3f_0(i,j,k)*(1._wp+r3f(i,j)))
+#   define  e3f_vor(i,j,k) (e3f_0vor(i,j,k)*(1._wp+r3f(i,j)))
 #   define  e3w(i,j,k,t)   (e3w_0(i,j,k)*(1._wp+r3t(i,j,t)))
 #   define  e3uw(i,j,k,t)  (e3uw_0(i,j,k)*(1._wp+r3u(i,j,t)))
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/dynvor.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/dynvor.F90	(revision 13707)
+++ 	(revision )
@@ -1,880 +1,0 @@
-MODULE dynvor
-   !!======================================================================
-   !!                       ***  MODULE  dynvor  ***
-   !! Ocean dynamics: Update the momentum trend with the relative and
-   !!                 planetary vorticity trends
-   !!======================================================================
-   !! History :  OPA  ! 1989-12  (P. Andrich)  vor_ens: Original code
-   !!            5.0  ! 1991-11  (G. Madec)  vor_ene, vor_mix: Original code
-   !!            6.0  ! 1996-01  (G. Madec)  s-coord, suppress work arrays
-   !!   NEMO     0.5  ! 2002-08  (G. Madec)  F90: Free form and module
-   !!            1.0  ! 2004-02  (G. Madec)  vor_een: Original code
-   !!             -   ! 2003-08  (G. Madec)  add vor_ctl
-   !!             -   ! 2005-11  (G. Madec)  add dyn_vor (new step architecture)
-   !!            2.0  ! 2006-11  (G. Madec)  flux form advection: add metric term
-   !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme
-   !!            3.3  ! 2010-10  (C. Ethe, G. Madec)  reorganisation of initialisation phase
-   !!            3.7  ! 2014-04  (G. Madec)  trend simplification: suppress jpdyn_trd_dat vorticity 
-   !!             -   ! 2014-06  (G. Madec)  suppression of velocity curl from in-core memory
-   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T)
-   !!            4.0  ! 2017-07  (G. Madec)  linear dynamics + trends diag. with Stokes-Coriolis
-   !!             -   ! 2018-03  (G. Madec)  add two new schemes (ln_dynvor_enT and ln_dynvor_eet)
-   !!             -   ! 2018-04  (G. Madec)  add pre-computed gradient for metric term calculation
-   !!            4.x  ! 2020-03  (G. Madec, A. Nasser)  make ln_dynvor_msk truly efficient on relative vorticity
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   dyn_vor       : Update the momentum trend with the vorticity trend
-   !!       vor_ens   : enstrophy conserving scheme       (ln_dynvor_ens=T)
-   !!       vor_ene   : energy conserving scheme          (ln_dynvor_ene=T)
-   !!       vor_een   : energy and enstrophy conserving   (ln_dynvor_een=T)
-   !!   dyn_vor_init  : set and control of the different vorticity option
-   !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and tracers
-   USE dom_oce        ! ocean space and time domain
-   USE dommsk         ! ocean mask
-   USE dynadv         ! momentum advection
-   USE trd_oce        ! trends: ocean variables
-   USE trddyn         ! trend manager: dynamics
-   USE sbcwave        ! Surface Waves (add Stokes-Coriolis force)
-   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force
-   !
-   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   USE prtctl         ! Print control
-   USE in_out_manager ! I/O manager
-   USE lib_mpp        ! MPP library
-   USE timing         ! Timing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   dyn_vor        ! routine called by step.F90
-   PUBLIC   dyn_vor_init   ! routine called by nemogcm.F90
-
-   !                                   !!* Namelist namdyn_vor: vorticity term
-   LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme          (ENS)
-   LOGICAL, PUBLIC ::   ln_dynvor_ene   !: f-point energy conserving scheme     (ENE)
-   LOGICAL, PUBLIC ::   ln_dynvor_enT   !: t-point energy conserving scheme     (ENT)
-   LOGICAL, PUBLIC ::   ln_dynvor_eeT   !: t-point energy conserving scheme     (EET)
-   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy & enstrophy conserving scheme (EEN)
-   INTEGER, PUBLIC ::      nn_een_e3f      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
-   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                         (MIX)
-   LOGICAL, PUBLIC ::   ln_dynvor_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes)
-
-   INTEGER, PUBLIC ::   nvor_scheme     !: choice of the type of advection scheme
-   !                                       ! associated indices:
-   INTEGER, PUBLIC, PARAMETER ::   np_ENS = 0   ! ENS scheme
-   INTEGER, PUBLIC, PARAMETER ::   np_ENE = 1   ! ENE scheme
-   INTEGER, PUBLIC, PARAMETER ::   np_ENT = 2   ! ENT scheme (t-point vorticity)
-   INTEGER, PUBLIC, PARAMETER ::   np_EET = 3   ! EET scheme (EEN using e3t)
-   INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme
-   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme
-
-   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity 
-   !                               ! associated indices:
-   INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary)
-   INTEGER, PUBLIC, PARAMETER ::   np_RVO = 2         ! relative vorticity
-   INTEGER, PUBLIC, PARAMETER ::   np_MET = 3         ! metric term
-   INTEGER, PUBLIC, PARAMETER ::   np_CRV = 4         ! relative + planetary (total vorticity)
-   INTEGER, PUBLIC, PARAMETER ::   np_CME = 5         ! Coriolis + metric term
-
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2u_2        ! = di(e2u)/2          used in T-point metric term calculation
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       - 
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2u)/(2*e1e2f)  used in F-point metric term calculation
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       - 
-   
-   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4
-   REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8
-   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12
-   
-   !! * Substitutions
-#  include "do_loop_substitute.h90"
-!!st23
-#  include "domzgr_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id: dynvor.F90 12667 2020-04-03 14:22:29Z gm $
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs )
-      !!----------------------------------------------------------------------
-      !!
-      !! ** Purpose :   compute the lateral ocean tracer physics.
-      !!
-      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
-      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
-      !!               and planetary vorticity trends) and send them to trd_dyn 
-      !!               for futher diagnostics (l_trddyn=T)
-      !!----------------------------------------------------------------------
-      INTEGER                             , INTENT( in  ) ::   kt          ! ocean time-step index
-      INTEGER                             , INTENT( in  ) ::   Kmm, Krhs   ! ocean time level indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv    ! ocean velocity field and RHS of momentum equation
-      !
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdu, ztrdv
-      !!----------------------------------------------------------------------
-      !
-      IF( ln_timing )   CALL timing_start('dyn_vor')
-      !
-      IF( l_trddyn ) THEN     !==  trend diagnostics case : split the added trend in two parts  ==!
-         !
-         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )
-         !
-         ztrdu(:,:,:) = puu(:,:,:,Krhs)            !* planetary vorticity trend (including Stokes-Coriolis force)
-         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
-         SELECT CASE( nvor_scheme )
-         CASE( np_ENS )           ;   CALL vor_ens( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! enstrophy conserving scheme
-            IF( ln_stcor )            CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
-         CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme
-            IF( ln_stcor )            CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
-         CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (T-pts)
-            IF( ln_stcor )            CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
-         CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (een with e3t)
-            IF( ln_stcor )            CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
-         CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy & enstrophy scheme
-            IF( ln_stcor )            CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
-         END SELECT
-         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
-         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
-         CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt, Kmm )
-         !
-         IF( n_dynadv /= np_LIN_dyn ) THEN   !* relative vorticity or metric trend (only in non-linear case)
-            ztrdu(:,:,:) = puu(:,:,:,Krhs)
-            ztrdv(:,:,:) = pvv(:,:,:,Krhs)
-            SELECT CASE( nvor_scheme )
-            CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme (T-pts)
-            CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme (een with e3t)
-            CASE( np_ENE )           ;   CALL vor_ene( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme
-            CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! enstrophy conserving scheme
-            CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy & enstrophy scheme
-            END SELECT
-            ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
-            ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
-            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt, Kmm )
-         ENDIF
-         !
-         DEALLOCATE( ztrdu, ztrdv )
-         !
-      ELSE              !==  total vorticity trend added to the general trend  ==!
-         !
-         SELECT CASE ( nvor_scheme )      !==  vorticity trend added to the general trend  ==!
-         CASE( np_ENT )                        !* energy conserving scheme  (T-pts)
-                             CALL vor_enT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
-         CASE( np_EET )                        !* energy conserving scheme (een scheme using e3t)
-                             CALL vor_eeT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
-         CASE( np_ENE )                        !* energy conserving scheme
-                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
-         CASE( np_ENS )                        !* enstrophy conserving scheme
-                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend
-            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend
-         CASE( np_MIX )                        !* mixed ene-ens scheme
-                             CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! relative vorticity or metric trend (ens)
-                             CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! planetary vorticity trend (ene)
-            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
-         CASE( np_EEN )                        !* energy and enstrophy conserving scheme
-                             CALL vor_een( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
-         END SELECT
-         !
-      ENDIF
-      !
-      !                       ! print sum trends (used for debugging)
-      IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' vor  - Ua: ', mask1=umask,               &
-         &                                tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-      !
-      IF( ln_timing )   CALL timing_stop('dyn_vor')
-      !
-   END SUBROUTINE dyn_vor
-
-
-   SUBROUTINE vor_enT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE vor_enT  ***
-      !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to 
-      !!      the general trend of the momentum equation.
-      !!
-      !! ** Method  :   Trend evaluated using now fields (centered in time) 
-      !!       and t-point evaluation of vorticity (planetary and relative).
-      !!       conserves the horizontal kinetic energy.
-      !!         The general trend of momentum is increased due to the vorticity 
-      !!       term which is given by:
-      !!          voru = 1/bu  mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t  mj[vn] ]
-      !!          vorv = 1/bv  mi[ ( mi(mj(bf*rvor))+bt*f_t)/e3f  mj[un] ]
-      !!       where rvor is the relative vorticity at f-point
-      !!
-      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
-      !!----------------------------------------------------------------------
-      INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index
-      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
-      INTEGER                         , INTENT(in   ) ::   kvor             ! total, planetary, relative, or metric
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv           ! now velocities
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! total v-trend
-      !
-      INTEGER  ::   ji, jj, jk           ! dummy loop indices
-      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx, zwy, zwt   ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz             ! 3D workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-      !
-      !
-      SELECT CASE( kvor )                 !== relative vorticity considered  ==!
-      CASE ( np_RVO , np_CRV )                  !* relative vorticity at f-point is used
-         DO jk = 1, jpkm1                                ! Horizontal slab
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
-                  &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
-            END_2D
-            IF( ln_dynvor_msk ) THEN                     ! mask relative vorticity 
-               DO_2D( 1, 0, 1, 0 )
-                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
-               END_2D
-            ENDIF
-         END DO
-         CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
-         !
-      END SELECT
-
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         !
-         SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
-         !
-         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm)
-         CASE ( np_RVO )                           !* relative vorticity
-            DO_2D( 0, 1, 0, 1 )
-               zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   &
-                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) &
-                  &                * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
-            END_2D
-         CASE ( np_MET )                           !* metric term
-            DO_2D( 0, 1, 0, 1 )
-               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)     &
-                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) &
-                  &         * e3t(ji,jj,jk,Kmm)
-            END_2D
-         CASE ( np_CRV )                           !* Coriolis + relative vorticity
-            DO_2D( 0, 1, 0, 1 )
-               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)      &
-                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) &
-                  &         * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
-            END_2D
-         CASE ( np_CME )                           !* Coriolis + metric
-            DO_2D( 0, 1, 0, 1 )
-               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                             &
-                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)    &
-                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) &
-                    &       * e3t(ji,jj,jk,Kmm)
-            END_2D
-         CASE DEFAULT                                             ! error
-            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
-         END SELECT
-         !
-         !                                   !==  compute and add the vorticity term trend  =!
-         DO_2D( 0, 0, 0, 0 )
-            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                    &
-               &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   &
-               &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   )
-               !
-            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                    &
-               &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
-               &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
-         END_2D
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-   END SUBROUTINE vor_enT
-
-
-   SUBROUTINE vor_ene( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE vor_ene  ***
-      !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to 
-      !!      the general trend of the momentum equation.
-      !!
-      !! ** Method  :   Trend evaluated using now fields (centered in time) 
-      !!       and the Sadourny (1975) flux form formulation : conserves the
-      !!       horizontal kinetic energy.
-      !!         The general trend of momentum is increased due to the vorticity 
-      !!       term which is given by:
-      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v pvv(:,:,:,Kmm)) ]
-      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u puu(:,:,:,Kmm)) ]
-      !!       where rvor is the relative vorticity
-      !!
-      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
-      !!
-      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
-      !!----------------------------------------------------------------------
-      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
-      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
-      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
-      !
-      INTEGER  ::   ji, jj, jk           ! dummy loop indices
-      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz   ! 2D workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-      !
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         !
-         SELECT CASE( kvor )                 !==  vorticity considered  ==!
-         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            zwz(:,:) = ff_f(:,:) 
-         CASE ( np_RVO )                           !* relative vorticity
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
-                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
-            END_2D
-            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
-               DO_2D( 1, 0, 1, 0 )
-                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
-               END_2D
-            ENDIF
-         CASE ( np_MET )                           !* metric term
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                  &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
-            END_2D
-         CASE ( np_CRV )                           !* Coriolis + relative vorticity
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
-                  &                        - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
-            END_2D
-            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity (NOT the Coriolis term)
-               DO_2D( 1, 0, 1, 0 )
-                  zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
-               END_2D
-            ENDIF
-         CASE ( np_CME )                           !* Coriolis + metric
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
-            END_2D
-         CASE DEFAULT                                             ! error
-            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
-         END SELECT
-         !
-         !                                   !==  horizontal fluxes and potential vorticity ==!
-         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
-         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
-         zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
-         !
-         !                                   !==  compute and add the vorticity term trend  =!
-         DO_2D( 0, 0, 0, 0 )
-            zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
-            zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
-            zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
-            zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
-            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
-            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
-         END_2D
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-   END SUBROUTINE vor_ene
-
-
-   SUBROUTINE vor_ens( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE vor_ens  ***
-      !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to
-      !!      the general trend of the momentum equation.
-      !!
-      !! ** Method  :   Trend evaluated using now fields (centered in time)
-      !!      and the Sadourny (1975) flux FORM formulation : conserves the
-      !!      potential enstrophy of a horizontally non-divergent flow. the
-      !!      trend of the vorticity term is given by:
-      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v pvv(:,:,:,Kmm)) ]
-      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u puu(:,:,:,Kmm)) ]
-      !!      Add this trend to the general momentum trend:
-      !!          (u(rhs),v(Krhs)) = (u(rhs),v(Krhs)) + ( voru , vorv )
-      !!
-      !! ** Action : - Update (pu_rhs,pv_rhs)) arrays with the now vorticity term trend
-      !!
-      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
-      !!----------------------------------------------------------------------
-      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
-      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
-      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zuav, zvau   ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz, zww   ! 2D workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         !
-         SELECT CASE( kvor )                 !==  vorticity considered  ==!
-         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            zwz(:,:) = ff_f(:,:) 
-         CASE ( np_RVO )                           !* relative vorticity
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
-                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
-            END_2D
-            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
-               DO_2D( 1, 0, 1, 0 )
-                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
-               END_2D
-            ENDIF
-         CASE ( np_MET )                           !* metric term
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                  &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
-            END_2D
-         CASE ( np_CRV )                           !* Coriolis + relative vorticity
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
-                  &                        - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
-            END_2D
-            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity (NOT the Coriolis term)
-               DO_2D( 1, 0, 1, 0 )
-                  zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
-               END_2D
-            ENDIF
-         CASE ( np_CME )                           !* Coriolis + metric
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
-            END_2D
-         CASE DEFAULT                                             ! error
-            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
-         END SELECT
-         !
-         !
-         !                                   !==  horizontal fluxes and potential vorticity ==!
-         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
-         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
-         zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
-         !
-         !                                   !==  compute and add the vorticity term trend  =!
-         DO_2D( 0, 0, 0, 0 )
-            zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  &
-               &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  )
-            zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  &
-               &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  )
-            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
-            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
-         END_2D
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-   END SUBROUTINE vor_ens
-
-
-   SUBROUTINE vor_een( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE vor_een  ***
-      !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to 
-      !!      the general trend of the momentum equation.
-      !!
-      !! ** Method  :   Trend evaluated using now fields (centered in time) 
-      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves 
-      !!      both the horizontal kinetic energy and the potential enstrophy
-      !!      when horizontal divergence is zero (see the NEMO documentation)
-      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
-      !!
-      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
-      !!
-      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
-      !!----------------------------------------------------------------------
-      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
-      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
-      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      INTEGER  ::   ierr         ! local integer
-      REAL(wp) ::   zua, zva     ! local scalars
-      REAL(wp) ::   zmsk, ze3f   ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx , zwy , z1_e3f
-      REAL(wp), DIMENSION(jpi,jpj)     ::   ztnw, ztne, ztsw, ztse
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-      !
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         !
-         SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point
-         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
-            DO_2D( 1, 0, 1, 0 )
-               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
-                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
-                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
-                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
-               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f
-               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
-               ENDIF
-            END_2D
-         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
-            DO_2D( 1, 0, 1, 0 )
-               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
-                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
-                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
-                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
-               zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
-                  &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  )
-               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3f
-               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
-               ENDIF
-            END_2D
-         END SELECT
-         !
-         SELECT CASE( kvor )                 !==  vorticity considered  ==!
-         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj)
-            END_2D
-         CASE ( np_RVO )                           !* relative vorticity
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj,jk) = ( e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
-                  &            - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
-            END_2D
-         CASE ( np_MET )                           !* metric term
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj,jk) = (   ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                  &              - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
-            END_2D
-         CASE ( np_CRV )                           !* Coriolis + relative vorticity
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
-                  &                              - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  )   &
-                  &                           * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
-            END_2D
-         CASE ( np_CME )                           !* Coriolis + metric
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj,jk) = (   ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                  &                            - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
-            END_2D
-         CASE DEFAULT                                             ! error
-            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
-         END SELECT
-         !
-         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
-            END_2D
-         ENDIF
-      END DO                                           !   End of slab
-         !
-      CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
-
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !
-         !                                   !==  horizontal fluxes  ==!
-         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
-         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
-
-         !                                   !==  compute and add the vorticity term trend  =!
-         jj = 2
-         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
-         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
-               ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
-               ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
-               ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
-               ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
-         END DO
-         DO jj = 3, jpj
-            DO ji = 2, jpi   ! vector opt. ok because we start at jj = 3
-               ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
-               ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
-               ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
-               ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
-            END DO
-         END DO
-         DO_2D( 0, 0, 0, 0 )
-            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
-               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
-            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
-               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
-            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
-            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
-         END_2D
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-   END SUBROUTINE vor_een
-
-
-
-   SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE vor_eeT  ***
-      !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to 
-      !!      the general trend of the momentum equation.
-      !!
-      !! ** Method  :   Trend evaluated using now fields (centered in time) 
-      !!      and the Arakawa and Lamb (1980) vector form formulation using 
-      !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
-      !!      The change consists in 
-      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
-      !!
-      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
-      !!
-      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
-      !!----------------------------------------------------------------------
-      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
-      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
-      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
-      !
-      INTEGER  ::   ji, jj, jk     ! dummy loop indices
-      INTEGER  ::   ierr           ! local integer
-      REAL(wp) ::   zua, zva       ! local scalars
-      REAL(wp) ::   zmsk, z1_e3t   ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx , zwy 
-      REAL(wp), DIMENSION(jpi,jpj)     ::   ztnw, ztne, ztsw, ztse
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-      !
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         !
-         !
-         SELECT CASE( kvor )                 !==  vorticity considered  ==!
-         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj,jk) = ff_f(ji,jj)
-            END_2D
-         CASE ( np_RVO )                           !* relative vorticity
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj,jk) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
-                  &             - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
-                  &          * r1_e1e2f(ji,jj)
-            END_2D
-         CASE ( np_MET )                           !* metric term
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj,jk) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                  &          - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
-            END_2D
-         CASE ( np_CRV )                           !* Coriolis + relative vorticity
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
-                  &                              - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
-                  &                         * r1_e1e2f(ji,jj)    )
-            END_2D
-         CASE ( np_CME )                           !* Coriolis + metric
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                  &                        - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
-            END_2D
-         CASE DEFAULT                                             ! error
-            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
-         END SELECT
-         !
-         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
-            END_2D
-         ENDIF
-      END DO
-      !
-      CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
-      !
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !
-         !                                   !==  horizontal fluxes  ==!
-         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
-         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
-
-         !                                   !==  compute and add the vorticity term trend  =!
-         jj = 2
-         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
-         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
-               z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
-               ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
-               ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
-               ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
-               ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
-         END DO
-         DO jj = 3, jpj
-            DO ji = 2, jpi   ! vector opt. ok because we start at jj = 3
-               z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
-               ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
-               ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
-               ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
-               ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
-            END DO
-         END DO
-         DO_2D( 0, 0, 0, 0 )
-            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
-               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
-            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
-               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
-            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
-            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
-         END_2D
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-   END SUBROUTINE vor_eeT
-
-
-   SUBROUTINE dyn_vor_init
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE dyn_vor_init  ***
-      !!
-      !! ** Purpose :   Control the consistency between cpp options for
-      !!              tracer advection schemes
-      !!----------------------------------------------------------------------
-      INTEGER ::   ji, jj, jk    ! dummy loop indices
-      INTEGER ::   ioptio, ios   ! local integer
-      !!
-      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   &
-         &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk
-      !!----------------------------------------------------------------------
-      !
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
-         WRITE(numout,*) '~~~~~~~~~~~~'
-      ENDIF
-      !
-      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_vor in reference namelist' )
-      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist' )
-      IF(lwm) WRITE ( numond, namdyn_vor )
-      !
-      IF(lwp) THEN                    ! Namelist print
-         WRITE(numout,*) '   Namelist namdyn_vor : choice of the vorticity term scheme'
-         WRITE(numout,*) '      enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
-         WRITE(numout,*) '      f-point energy conserving scheme               ln_dynvor_ene = ', ln_dynvor_ene
-         WRITE(numout,*) '      t-point energy conserving scheme               ln_dynvor_enT = ', ln_dynvor_enT
-         WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT
-         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
-         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f
-         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
-         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
-      ENDIF
-
-!!an      IF( ln_dynvor_msk )   CALL ctl_stop( 'dyn_vor_init:   masked vorticity is not currently not available')
-
-!!gm  this should be removed when choosing a unique strategy for fmask at the coast
-      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
-      ! at angles with three ocean points and one land point
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '      change fmask value in the angles (T)           ln_vorlat = ', ln_vorlat
-      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
-         DO_3D( 1, 0, 1, 0, 1, jpk )
-            IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              &
-               & + tmask(ji,jj  ,jk) + tmask(ji+1,jj+1,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
-         END_3D
-         !
-         CALL lbc_lnk( 'dynvor', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
-         !
-      ENDIF
-!!gm end
-
-      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme)
-      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS   ;   ENDIF
-      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE   ;   ENDIF
-      IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF
-      IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF
-      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEN   ;   ENDIF
-      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_MIX   ;   ENDIF
-      !
-      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
-      !                      
-      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
-      ncor = np_COR                       ! planetary vorticity
-      SELECT CASE( n_dynadv )
-      CASE( np_LIN_dyn )
-         IF(lwp) WRITE(numout,*) '   ==>>>   linear dynamics : total vorticity = Coriolis'
-         nrvm = np_COR        ! planetary vorticity
-         ntot = np_COR        !     -         -
-      CASE( np_VEC_c2  )
-         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity' 
-         nrvm = np_RVO        ! relative vorticity
-         ntot = np_CRV        ! relative + planetary vorticity         
-      CASE( np_FLX_c2 , np_FLX_ubs  )
-         IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term'
-         nrvm = np_MET        ! metric term
-         ntot = np_CME        ! Coriolis + metric term
-         !
-         SELECT CASE( nvor_scheme )    ! pre-computed gradients for the metric term:
-         CASE( np_ENT )                      !* T-point metric term :   pre-compute di(e2u)/2 and dj(e1v)/2
-            ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
-            DO_2D( 0, 0, 0, 0 )
-               di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj  ) ) * 0.5_wp
-               dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji  ,jj-1) ) * 0.5_wp
-            END_2D
-            CALL lbc_lnk_multi( 'dynvor', di_e2u_2, 'T', -1. , dj_e1v_2, 'T', -1. )   ! Lateral boundary conditions
-            !
-         CASE DEFAULT                        !* F-point metric term :   pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f)
-            ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
-            DO_2D( 1, 0, 1, 0 )
-               di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj  ) - e2v(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
-               dj_e1u_2e1e2f(ji,jj) = ( e1u(ji  ,jj+1) - e1u(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
-            END_2D
-            CALL lbc_lnk_multi( 'dynvor', di_e2v_2e1e2f, 'F', -1. , dj_e1u_2e1e2f, 'F', -1. )   ! Lateral boundary conditions
-         END SELECT
-         !
-      END SELECT
-      
-      IF(lwp) THEN                   ! Print the choice
-         WRITE(numout,*)
-         SELECT CASE( nvor_scheme )
-         CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)'
-         CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)'
-         CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)'
-         CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)'
-         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)'
-         CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)'
-         END SELECT         
-      ENDIF
-      !
-   END SUBROUTINE dyn_vor_init
-
-   !!==============================================================================
-END MODULE dynvor
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/nemogcm.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/nemogcm.F90	(revision 13707)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/nemogcm.F90	(revision 13708)
@@ -26,5 +26,5 @@
    USE stpRK3         ! NEMO time-stepping               (stp_RK3   routine)
 #else
-   USE stpMLF         ! NEMO time-stepping               (stp_MLF   routine)
+   USE stpmlf         ! NEMO time-stepping               (stp_MLF   routine)
 #endif
    !
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/stpMLF.F90	(revision 13707)
+++ 	(revision )
@@ -1,253 +1,0 @@
-MODULE stpMLF
-   !!======================================================================
-   !!                       ***  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
-   !!             -   !  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 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_MLF   ! called by nemogcm.F90
-   
-   !                                          !**  time level indices  **!
-   INTEGER, PUBLIC ::   Nbb, Nnn, Naa, Nrhs   !: used by nemo_init
-      
-   !! * Substitutions
-#  include "do_loop_substitute.h90"
-#  include "domzgr_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id: step.F90 12614 2020-03-26 14:59:52Z gm $
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE stp_MLF( kstp )
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE stp_MLF  ***
-      !!
-      !! ** Purpose : - Time stepping of shallow water (SHW) (momentum and ssh eqs.)
-      !!
-      !! ** Method  : -1- Update forcings
-      !!              -2- Update the ssh at Naa
-      !!              -3- Compute the momentum trends (Nrhs)
-      !!              -4- Update the horizontal velocity
-      !!              -5- Apply Asselin time filter to uu,vv,ssh
-      !!              -6- Outputs and diagnostics
-      !!----------------------------------------------------------------------
-      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( ln_timing )   CALL timing_start('stp_MLF')
-      !
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! model timestep
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !
-      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
-
-      !
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! 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)
-                             CALL iom_init_closedef
-      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
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! 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. 
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      !  RHS of horizontal velocity
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
-                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
-                         vv(:,:,:,Nrhs) = 0._wp
-
-                         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 )
-         !                                                              ! 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
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! 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
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! 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)
-               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( 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)
-               !                                                                  ! Asselin time filter on u,v (Nnn)
-               uu(ji,jj,jk,Nnn) = uu(ji,jj,jk,Nnn) + rn_atfp * (uu(ji,jj,jk,Nbb) - 2._wp * uu(ji,jj,jk,Nnn) + zua)
-               vv(ji,jj,jk,Nnn) = vv(ji,jj,jk,Nnn) + rn_atfp * (vv(ji,jj,jk,Nbb) - 2._wp * vv(ji,jj,jk,Nnn) + zva)
-               !              
-               uu(ji,jj,jk,Naa) = zua
-               vv(ji,jj,jk,Naa) = zva
-            END_3D
-            !                            ! 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)
-               zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
-               zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
-               !                                                ! 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)    
-               vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)
-            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( 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 * ( 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 ) ) / ( 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 / ( 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(:,:)
-         ENDIF
-      ENDIF
-
-      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
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
-      ! Swap time levels
-      Nrhs = Nbb
-      Nbb = Nnn
-      Nnn = Naa
-      Naa = Nrhs
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! diagnostics and outputs
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      
-      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, Nnn )
-
-      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
-
-      !
-#if defined key_iomput
-      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   
-         r1_Dt = 1._wp / rDt
-         l_1st_euler = .FALSE.      
-      ENDIF
-      !
-      IF( ln_timing )   CALL timing_stop('stp_MLF')
-      !
-   END SUBROUTINE stp_MLF
-
-   !!======================================================================
-END MODULE stpMLF
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/stpmlf.F90	(revision 13708)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/stpmlf.F90	(revision 13708)
@@ -0,0 +1,253 @@
+MODULE stpmlf
+   !!======================================================================
+   !!                       ***  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
+   !!             -   !  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 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_MLF   ! called by nemogcm.F90
+   
+   !                                          !**  time level indices  **!
+   INTEGER, PUBLIC ::   Nbb, Nnn, Naa, Nrhs   !: used by nemo_init
+      
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: step.F90 12614 2020-03-26 14:59:52Z gm $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE stp_MLF( kstp )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE stp_MLF  ***
+      !!
+      !! ** Purpose : - Time stepping of shallow water (SHW) (momentum and ssh eqs.)
+      !!
+      !! ** Method  : -1- Update forcings
+      !!              -2- Update the ssh at Naa
+      !!              -3- Compute the momentum trends (Nrhs)
+      !!              -4- Update the horizontal velocity
+      !!              -5- Apply Asselin time filter to uu,vv,ssh
+      !!              -6- Outputs and diagnostics
+      !!----------------------------------------------------------------------
+      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( ln_timing )   CALL timing_start('stp_MLF')
+      !
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! model timestep
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      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
+
+      !
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! 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)
+                             CALL iom_init_closedef
+      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
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! 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. 
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      !  RHS of horizontal velocity
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
+                         vv(:,:,:,Nrhs) = 0._wp
+
+                         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 )
+         !                                                              ! 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
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! 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
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! 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)
+               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( 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)
+               !                                                                  ! Asselin time filter on u,v (Nnn)
+               uu(ji,jj,jk,Nnn) = uu(ji,jj,jk,Nnn) + rn_atfp * (uu(ji,jj,jk,Nbb) - 2._wp * uu(ji,jj,jk,Nnn) + zua)
+               vv(ji,jj,jk,Nnn) = vv(ji,jj,jk,Nnn) + rn_atfp * (vv(ji,jj,jk,Nbb) - 2._wp * vv(ji,jj,jk,Nnn) + zva)
+               !              
+               uu(ji,jj,jk,Naa) = zua
+               vv(ji,jj,jk,Naa) = zva
+            END_3D
+            !                            ! 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)
+               zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
+               zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
+               !                                                ! 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)    
+               vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)
+            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( 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 * ( 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 ) ) / ( 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 / ( 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(:,:)
+         ENDIF
+      ENDIF
+
+      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
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+      ! Swap time levels
+      Nrhs = Nbb
+      Nbb = Nnn
+      Nnn = Naa
+      Naa = Nrhs
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! diagnostics and outputs
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      
+      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, Nnn )
+
+      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
+
+      !
+#if defined key_iomput
+      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   
+         r1_Dt = 1._wp / rDt
+         l_1st_euler = .FALSE.      
+      ENDIF
+      !
+      IF( ln_timing )   CALL timing_stop('stp_MLF')
+      !
+   END SUBROUTINE stp_MLF
+
+   !!======================================================================
+END MODULE stpmlf
