Index: /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynadv.F90
===================================================================
--- /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynadv.F90	(revision 13153)
+++ /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynadv.F90	(revision 13154)
@@ -74,4 +74,5 @@
       CASE( np_VEC_c2  )     
          CALL dyn_keg     ( kt, nn_dynkeg,      Kmm, puu, pvv, Krhs )    ! vector form : horizontal gradient of kinetic energy
+!!an SWE : w = 0
          CALL dyn_zad     ( kt,                 Kmm, puu, pvv, Krhs )    ! vector form : vertical advection
       CASE( np_FLX_c2  ) 
Index: /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynkeg.F90
===================================================================
--- /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynkeg.F90	(revision 13153)
+++ /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynkeg.F90	(revision 13154)
@@ -29,4 +29,5 @@
 
    PUBLIC   dyn_keg    ! routine called by step module
+   PUBLIC   dyn_kegAD   ! routine called by step module
    
    INTEGER, PARAMETER, PUBLIC  ::   nkeg_C2     = 0   !: 2nd order centered scheme (standard scheme)
@@ -155,5 +156,116 @@
       !
    END SUBROUTINE dyn_keg
-
+   
+   
+   SUBROUTINE dyn_kegAD( kt, kscheme, puu, pvv, pu_rhs, pv_rhs )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_kegAD  ***
+      !!
+      !! ** Purpose :   Compute the now momentum trend due to the horizontal
+      !!      gradient of the horizontal kinetic energy and add it to the 
+      !!      general momentum trend.
+      !!
+      !! ** Method  : * kscheme = nkeg_C2 : 2nd order centered scheme that 
+      !!      conserve kinetic energy. Compute the now horizontal kinetic energy 
+      !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
+      !!              * kscheme = nkeg_HW : Hollingsworth correction following
+      !!      Arakawa (2001). The now horizontal kinetic energy is given by:
+      !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((u(j+1)+u(j-1))/2)^2  )
+      !!                    + mj-1(  2 * vn^2 + ((v(i+1)+v(i-1))/2)^2  ) ]
+      !!      
+      !!      Take its horizontal gradient and add it to the general momentum
+      !!      trend.
+      !!         u(rhs) = u(rhs) - 1/e1u di[ zhke ]
+      !!         v(rhs) = v(rhs) - 1/e2v dj[ zhke ]
+      !!
+      !! ** Action : - Update the (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) with the hor. ke gradient trend
+      !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing
+      !!
+      !! ** References : Arakawa, A., International Geophysics 2001.
+      !!                 Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983.
+      !!----------------------------------------------------------------------
+      !
+      INTEGER                                  , INTENT( in )  ::  kt               ! ocean time-step index
+      INTEGER                                  , INTENT( in )  ::  kscheme          ! =0/1/2   type of KEG scheme 
+      REAL(wp), DIMENSION(jpi,jpj,jpk)         , INTENT(inout) ::  puu, pvv         ! ocean velocities at Kmm
+      REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) ::  pu_rhs, pv_rhs   ! RHS 
+      !
+      INTEGER  ::   ji, jj, jk             ! dummy loop indices
+      REAL(wp) ::   zu, zv                   ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zhke
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('dyn_keg')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dyn_kegAD : kinetic energy gradient trend, scheme number=', kscheme
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~'
+      ENDIF
+      
+      zhke(:,:,jpk) = 0._wp
+
+      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
+!!an45 to be ADDED : que cas C2 - "wet points only" il suffit de x2 le terme quadratic a la coast (nn_dynkeg_adv = 2)
+      CASE ( nkeg_C2_wpo )                          !--  Standard scheme  --!
+         DO_3D_01_01( 1, jpkm1 )
+            zu =  (   puu(ji-1,jj  ,jk) * puu(ji-1,jj  ,jk)   &
+               &    + puu(ji  ,jj  ,jk) * puu(ji  ,jj  ,jk)   ) * ( 2._wp - umask(ji-1,jj,jk) * umask(ji,jj,jk) )
+            zv =  (   pvv(ji  ,jj-1,jk) * pvv(ji  ,jj-1,jk)   &
+               &    + pvv(ji  ,jj  ,jk) * pvv(ji  ,jj  ,jk)   ) * ( 2._wp - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )
+            zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
+         END_3D
+!!an45         
+      !
+      CASE ( nkeg_C2 )                          !--  Standard scheme  --!
+         DO_3D_01_01( 1, jpkm1 )
+            zu =    puu(ji-1,jj  ,jk) * puu(ji-1,jj  ,jk)   &
+               &  + puu(ji  ,jj  ,jk) * puu(ji  ,jj  ,jk)
+            zv =    pvv(ji  ,jj-1,jk) * pvv(ji  ,jj-1,jk)   &
+               &  + pvv(ji  ,jj  ,jk) * pvv(ji  ,jj  ,jk)
+            zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
+         END_3D
+!!an 00_00 ?
+      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
+         DO_3D_00_00( 1, jpkm1 )
+            zu = 8._wp * ( puu(ji-1,jj  ,jk) * puu(ji-1,jj  ,jk)    &
+               &         + puu(ji  ,jj  ,jk) * puu(ji  ,jj  ,jk) )  &
+               &   +     ( puu(ji-1,jj-1,jk) + puu(ji-1,jj+1,jk) ) * ( puu(ji-1,jj-1,jk) + puu(ji-1,jj+1,jk) )   &
+               &   +     ( puu(ji  ,jj-1,jk) + puu(ji  ,jj+1,jk) ) * ( puu(ji  ,jj-1,jk) + puu(ji  ,jj+1,jk) )
+               !
+            zv = 8._wp * ( pvv(ji  ,jj-1,jk) * pvv(ji  ,jj-1,jk)    &
+               &         + pvv(ji  ,jj  ,jk) * pvv(ji  ,jj  ,jk) )  &
+               &  +      ( pvv(ji-1,jj-1,jk) + pvv(ji+1,jj-1,jk) ) * ( pvv(ji-1,jj-1,jk) + pvv(ji+1,jj-1,jk) )   &
+               &  +      ( pvv(ji-1,jj  ,jk) + pvv(ji+1,jj  ,jk) ) * ( pvv(ji-1,jj  ,jk) + pvv(ji+1,jj  ,jk) )
+            zhke(ji,jj,jk) = r1_48 * ( zv + zu )
+         END_3D
+         CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. )
+         !
+      END SELECT 
+      !
+         IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN     !***  NO alternating direction  ***!
+            !
+            DO_3D_00_00( 1, jpkm1 )
+               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
+               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
+            END_3D
+            !
+         ELSEIF(       PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : i-component  ***!
+            !
+            DO_3D_00_00( 1, jpkm1 )
+               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
+            END_3D
+            !
+         ELSEIF( .NOT. PRESENT( pu_rhs ) .AND.       PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : j-component  ***!
+            !
+            DO_3D_00_00( 1, jpkm1 )
+               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
+            END_3D
+            !
+         ENDIF
+      IF( ln_timing )   CALL timing_stop('dyn_kegAD')
+      !
+   END SUBROUTINE dyn_kegAD
    !!======================================================================
 END MODULE dynkeg
Index: /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90
===================================================================
--- /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90	(revision 13153)
+++ /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90	(revision 13154)
@@ -22,4 +22,6 @@
    !!             -   ! 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
+   !!            4.x  ! 2020-03  (G. Madec, A. Nasser)  alternate direction computation of vorticity tendancy
+   !!                 !                                 for ENS, ENE
    !!----------------------------------------------------------------------
 
@@ -45,4 +47,7 @@
    USE lib_mpp        ! MPP library
    USE timing         ! Timing
+!!anAD only
+   USE dynkeg, ONLY : dyn_kegAD
+   USE dynadv, ONLY : nn_dynkeg
 
    IMPLICIT NONE
@@ -53,12 +58,18 @@
 
    !                                   !!* 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)
+   LOGICAL, PUBLIC ::   ln_dynvor_ens      !: enstrophy conserving scheme          (ENS)
+   LOGICAL, PUBLIC ::   ln_dynvor_ens_adVO   = .FALSE.   !: AD enstrophy conserving scheme       (ENS_ad)
+   LOGICAL, PUBLIC ::   ln_dynvor_ens_adKE   = .FALSE.   !: AD enstrophy conserving scheme       (ENS_ad)
+   LOGICAL, PUBLIC ::   ln_dynvor_ens_adKEVO = .FALSE.   !: AD enstrophy conserving scheme       (ENS_ad)
+   LOGICAL, PUBLIC ::   ln_dynvor_ene      !: f-point energy conserving scheme     (ENE)
+   LOGICAL, PUBLIC ::   ln_dynvor_ene_adVO   = .FALSE.   !: f-point AD energy conserving scheme  (ENE_ad)
+   LOGICAL, PUBLIC ::   ln_dynvor_ene_adKE   = .FALSE.   !: f-point AD energy conserving scheme  (ENE_ad)
+   LOGICAL, PUBLIC ::   ln_dynvor_ene_adKEVO = .FALSE.   !: f-point AD energy conserving scheme  (ENE_ad)
+   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
@@ -70,6 +81,15 @@
    INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme
    INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme
-
-   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity 
+!!an
+   INTEGER, PUBLIC, PARAMETER ::   np_ENS_adKE   = 11   ! ENS scheme - AD scheme (KE only)
+   INTEGER, PUBLIC, PARAMETER ::   np_ENS_adVO   = 12   ! ENS scheme - AD scheme (VOR only)
+   INTEGER, PUBLIC, PARAMETER ::   np_ENS_adKEVO = 13   ! ENS scheme - AD scheme (KE+VOR)
+   INTEGER, PUBLIC, PARAMETER ::   np_ENE_adKE   = 21   ! ENE scheme - AD scheme (KE only)
+   INTEGER, PUBLIC, PARAMETER ::   np_ENE_adVO   = 22   ! ENE scheme - AD scheme (VOR only)
+   INTEGER, PUBLIC, PARAMETER ::   np_ENE_adKEVO = 23   ! ENE scheme - AD scheme (KE+VOR)
+!!an      
+   
+!!an ds step on pourra spécifier la valeur de ntot = np_COR ou np_COR + np_RVO
+   INTEGER, PUBLIC ::   ncor, nrvm, ntot   ! choice of calculated vorticity 
    !                               ! associated indices:
    INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary)
@@ -99,5 +119,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs )
+   SUBROUTINE dyn_vor( kt, Kbb, Kmm, puu, pvv, Krhs)
       !!----------------------------------------------------------------------
       !!
@@ -109,9 +129,11 @@
       !!               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
+      INTEGER ::   ji, jj, jk   ! dummy loop indice
+      INTEGER                             , INTENT( in  ) ::   kt               ! ocean time-step index
+      INTEGER                             , INTENT( in  ) ::   Kmm, Krhs, Kbb   ! 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
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zuu, zvv
       !!----------------------------------------------------------------------
       !
@@ -167,8 +189,162 @@
             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 dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )  
+                             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_ENE_adVO )                     !* energy conserving scheme with alternating direction
+            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components
+                              
+                              !==  Alternative direction - VOR only  ==!
+
+                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
+
+                             ALLOCATE( zuu(jpi,jpj,jpk) )
+                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
+                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)
+                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )
+                             CALL vor_ene( kt, Kmm, ntot, zuu(:,:,:)     , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
+                             DEALLOCATE( zuu )
+            ELSE
+                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
+
+                             ALLOCATE( zvv(jpi,jpj,jpk) ) 
+                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
+                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)
+                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )
+                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
+                             DEALLOCATE( zvv )
+            ENDIF
+            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+         CASE( np_ENE_adKE )                     !* energy conserving scheme with alternating direction
+            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components
+                                 
+                                 !==  Alternative direction - KEG only  ==!
+
                              CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
+
+                             ALLOCATE( zuu(jpi,jpj,jpk) )
+                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
+                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)
+                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )
+                             CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:)     , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
+                             DEALLOCATE( zuu )
+            ELSE
+
+                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend                      
+
+                             ALLOCATE( zvv(jpi,jpj,jpk) )
+                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
+                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)
+                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )
+                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm)  , zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
+                             DEALLOCATE( zvv )
+            ENDIF
             IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            
+         CASE( np_ENE_adKEVO )                     !* energy conserving scheme with alternating direction
+         ! equivalent ADE pas sur coriolis
+         ! appel vor_ene(ncor)
+            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components
+                              
+                              !==  Alternative direction - KE + VOR  ==!
+
+                             ALLOCATE( zuu(jpi,jpj,jpk) )
+                             !puis call que sur nrvm 
+                             CALL vor_ene(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
+                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! 
+                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)
+                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )
+                             !puis call que sur nrvm 
+                             CALL vor_ene(   kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
+                             CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   
+                             DEALLOCATE( zuu )
+            ELSE
+
+                             ALLOCATE( zvv(jpi,jpj,jpk) ) 
+                             CALL vor_ene(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
+                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   
+                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)
+                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )
+                             CALL vor_ene(   kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
+                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )   
+                             DEALLOCATE( zvv )
+            ENDIF   
          CASE( np_ENS )                        !* enstrophy conserving scheme
+                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )  
+                             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_ENS_adVO )                     !* enstrophy conserving scheme with alternating direction
+            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components
+   
+                                 !==  Alternative direction - VOR only  ==!
+
+                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )  
+
+                             ALLOCATE( zuu(jpi,jpj,jpk) )
+                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
+                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)
+                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )
+                             CALL vor_ens( kt, Kmm, ntot, zuu(:,:,:)     , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
+                             DEALLOCATE( zuu )
+            ELSE
+                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )
+
+                             ALLOCATE( zvv(jpi,jpj,jpk) )
+                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
+                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)
+                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )
+                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , zvv(:,:,:), pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
+                             DEALLOCATE( zvv )
+            ENDIF
+            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend
+         CASE( np_ENS_adKE )                     !* enstrophy conserving scheme with alternating direction
+            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components
+            
+                                !==  Alternative direction - KEG only  ==!
+
                              CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend
+
+                             ALLOCATE( zuu(jpi,jpj,jpk) )
+                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
+                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)
+                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )
+                             CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:)     , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
+                             DEALLOCATE( zuu )
+            ELSE
+
+                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend                                                          
+
+                             ALLOCATE( zvv(jpi,jpj,jpk) )
+                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
+                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)
+                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )
+                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , zvv(:,:,:)  , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
+                             DEALLOCATE( zvv )
+            ENDIF
+         CASE( np_ENS_adKEVO )                     !* enstrophy conserving scheme with alternating direction
+            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components
+                              
+                              !==  Alternative direction - KE + VOR  ==!
+
+                             ALLOCATE( zuu(jpi,jpj,jpk) )
+                             CALL vor_ens(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
+                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! 
+                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)
+                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )
+                             CALL vor_ens(   kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
+                             CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   
+                             DEALLOCATE( zuu )
+            ELSE
+
+                             ALLOCATE( zvv(jpi,jpj,jpk) ) 
+                             CALL vor_ens(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend
+                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   
+                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)
+                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )
+                             CALL vor_ens(   kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend
+                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )   
+                             DEALLOCATE( zvv )
+            ENDIF   
             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
@@ -319,13 +495,14 @@
       !! 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                                  , 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),OPTIONAL, 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
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t        ! unpenalised scale factors
       !!----------------------------------------------------------------------
       !
@@ -377,18 +554,56 @@
          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_00_00
-            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
+         IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN     !***  NO alternating direction  ***!
+            !
+            !                                   !==  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_00_00
+               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
+            !            
+            !
+         ELSEIF(       PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : i-component  ***!
+            !
+            !
+            !                                   !==  horizontal fluxes and potential vorticity ==!
+!!anAD remplacer pu par uu(Nnn) et on fait correctement la direction alternée 
+!      on transportera (en AD) de la bonne manière
+            !zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
+            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * vv(:,:,jk,Kmm)
+            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
+            !
+            !                                   !==  compute and add the vorticity term trend  =!
+            DO_2D_00_00
+               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
+               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
+               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 )
+            END_2D
+            !
+         ELSEIF( .NOT. PRESENT( pu_rhs ) .AND.       PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : j-component  ***!
+            !
+            !
+            !                                   !==  horizontal fluxes and potential vorticity ==!
+!!anAD remplacer pu par uu(Nnn) et on fait correctement la direction alternée 
+            !zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
+            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * uu(:,:,jk,Kmm)
+            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
+            !
+            !                                   !==  compute and add the vorticity term trend  =!
+            DO_2D_00_00
+               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
+               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
+               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
+            !
+         ENDIF
          !                                             ! ===============
       END DO                                           !   End of slab
@@ -417,9 +632,9 @@
       !! 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                                  , 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),OPTIONAL, INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -475,18 +690,51 @@
          !
          !
-         !                                   !==  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_00_00
-            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
+!!an wut ? v et u 
+         IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN     !***  NO alternating direction  ***!
+            !
+            !                                   !==  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_00_00
+               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
+            !
+         ELSEIF(       PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : i-component  ***!
+            !
+            !
+            !                                   !==  horizontal fluxes and potential vorticity ==!
+            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
+            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
+            !
+            !                                   !==  compute and add the vorticity term trend  =!
+            DO_2D_00_00
+               zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  &
+                  &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  )
+               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
+            END_2D
+            !
+         ELSEIF( .NOT. PRESENT( pu_rhs ) .AND.       PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : j-component  ***!
+            !
+            !
+            !                                   !==  horizontal fluxes and potential vorticity ==!
+            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
+            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
+            !
+            !                                   !==  compute and add the vorticity term trend  =!
+            DO_2D_00_00
+               zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  &
+                  &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  )
+               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
+            END_2D
+            !
+         ENDIF
          !                                             ! ===============
       END DO                                           !   End of slab
@@ -772,5 +1020,7 @@
       !!
       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
+         &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk,   &
+         &                 ln_dynvor_ens_adVO, ln_dynvor_ens_adKE, ln_dynvor_ens_adKEVO, &   ! Alternative direction parameters
+         &                 ln_dynvor_ene_adVO, ln_dynvor_ene_adKE, ln_dynvor_ene_adKEVO
       !!----------------------------------------------------------------------
       !
@@ -809,5 +1059,5 @@
          DO_3D_10_10( 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
+               & + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
          END_3D
          !
@@ -819,5 +1069,11 @@
       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_ens_adVO ) THEN   ;   ioptio = ioptio + 1     ;   nvor_scheme = np_ENS_adVO   ;   ENDIF
+      IF( ln_dynvor_ens_adKE ) THEN   ;   ioptio = ioptio + 1     ;   nvor_scheme = np_ENS_adKE   ;   ENDIF
+      IF( ln_dynvor_ens_adKEVO ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS_adKEVO   ;   ENDIF
       IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE   ;   ENDIF
+      IF( ln_dynvor_ene_adVO ) THEN   ;   ioptio = ioptio + 1     ;   nvor_scheme = np_ENE_adVO   ;   ENDIF
+      IF( ln_dynvor_ene_adKE ) THEN   ;   ioptio = ioptio + 1     ;   nvor_scheme = np_ENE_adKE   ;   ENDIF
+      IF( ln_dynvor_ene_adKEVO ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE_adKEVO   ;   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
@@ -866,10 +1122,19 @@
          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)'
+         
+         CASE( np_ENS    )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)'
+         CASE( np_ENS_adVO ) ;   WRITE(numout,*) '   ==>>>   AD enstrophy conserving scheme (ENS_adVO) on vorticity only'
+         CASE( np_ENS_adKE ) ;   WRITE(numout,*) '   ==>>>   AD enstrophy conserving scheme (ENS_adKE) on kinetic energy only'
+         CASE( np_ENS_adKEVO ) ;   WRITE(numout,*) '   ==>>>   AD enstrophy conserving scheme (ENS_adKEVO) on kinetic energy and vorticity'
+         
+         CASE( np_ENE    )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)'
+         CASE( np_ENE_adVO ) ;   WRITE(numout,*) '   ==>>>   AD energy conserving scheme (Coriolis at F-points) (ENE_adVO) on vorticity only'
+         CASE( np_ENE_adKE ) ;   WRITE(numout,*) '   ==>>>   AD energy conserving scheme (Coriolis at F-points) (ENE_adKE) on kinetic energy only'
+         CASE( np_ENE_adKEVO ) ;   WRITE(numout,*) '   ==>>>   AD energy conserving scheme (Coriolis at F-points) (ENE_adKEVO) on kinetic energy and vorticity'
+         
+         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
Index: /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90
===================================================================
--- /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90	(revision 13153)
+++ /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90	(revision 13154)
@@ -17,5 +17,5 @@
    USE phycst           ! physical constants
    USE usrdef_nam
-   !
+   USE lib_mpp        ! MPP library
    USE iom              ! xIOs server 
 
@@ -66,5 +66,7 @@
       REAL(wp) ::   zue3a, zue3n, zue3b    ! local scalars
       REAL(wp) ::   zve3a, zve3n, zve3b    !   -      -
-      REAL(wp) ::   ze3t_tf, ze3u_tf, ze3v_tf, zua, zva
+      REAL(wp) ::   ze3t_tf, ze3u_tf, ze3v_tf, zua, zva     
+      REAL(wp), DIMENSION(jpi,jpj)    ::   zssh
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3u, ze3v
       !! ---------------------------------------------------------------------
 #if defined key_agrif
@@ -141,25 +143,61 @@
                &            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
 
 !!an - calcul du gradient de pression horizontal (explicit)
+      zssh(:,:) = ssh(:,:,Nnn)
+# if defined key_bvp
+      !   gradient and divergence not penalised
+      zssh(:,:) = zssh(:,:) * r1_rpo(:,:,1)  
+#endif 
       DO_3D_00_00( 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)
+         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( zssh(ji+1,jj) - zssh(ji,jj) ) * r1_e1u(ji,jj)
+         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( zssh(ji,jj+1) - zssh(ji,jj) ) * r1_e2v(ji,jj)
       END_3D
       !
-      ! add wind stress forcing and layer linear friction to the RHS 
+      
+!      IF( kstp == nit000 .AND. lwp ) THEN
+!         WRITE(numout,*)
+!         WRITE(numout,*) 'step.F90 : classic script used'
+!         WRITE(numout,*) '~~~~~~~'
+!         IF(       ln_dynvor_ens_adVO .OR. ln_dynvor_ens_adKE .OR. ln_dynvor_ens_adKEVO   &
+!         &    .OR. ln_dynvor_ene_adVO .OR. ln_dynvor_ene_adKE .OR. ln_dynvor_ene_adKEVO   ) THEN
+!            CALL ctl_stop('STOP','step : alternative direction asked but classis step'  )
+!         ENDIF
+!      ENDIF
+!!an     
+!                         CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)	==> RHS
+! 
+!                         CALL dyn_vor( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! vorticity           	==> RHS
+! 
+!!an     In dynvor, dynkegAD is called even if not AD, so we keep the same step.F90
+  
+                         CALL dyn_vor( kstp, Nbb, Nnn      , uu, vv, Nrhs)  ! vorticity           	==> RHS
+  
+                         CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing
+
+      ! add wind stress forcing and layer linear friction to the RHS
+      ze3u(:,:,:) = e3u(:,:,:,Nnn)
+      ze3v(:,:,:) = e3v(:,:,:,Nnn)
+# if defined key_bvp
+      !ze3u(:,:,:) = ze3u(:,:,:) * r1_rpo(:,:,:)  
+      !ze3v(:,:,:) = ze3v(:,:,:) * r1_rpo(:,:,:)  
+#endif  
       z1_2rho0 = 0.5_wp * r1_rho0
       DO_3D_00_00(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)   &
+         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / ze3u(ji,jj,jk)   &
             &                                  - 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)
+         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ze3v(ji,jj,jk)   &
+            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)  
       END_3D   
-!!an         
+      !
+# if defined key_bvp
+      !  Add frictionnal term   - sigma * u
+      ! can be done via sigmaT (mean_ij)
+      DO_3D_00_00(1,jpkm1)
+         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - bmpu(ji,jj,jk) * uu(ji,jj,jk,Nbb)
+         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - bmpv(ji,jj,jk) * vv(ji,jj,jk,Nbb)
+      END_3D 
+#endif 
+      !         
   
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
