Index: /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/EXPREF/namelist_cfg	(revision 13004)
+++ /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/EXPREF/namelist_cfg	(revision 13005)
@@ -23,8 +23,7 @@
    nn_it000    =       1   !  first time step
    nn_itend    =  1        !  10 ans -  2 min - last  time step
-!   nn_itend    =  2592000  !  10 ans -  2 min - last  time step
 !   nn_itend    =  1036800  !  10 ans -  5 min
-!   nn_itend    =    34560  !  1  an  - 15 min
-!   nn_itend    =   172800  !  10 ans  - 30 min
+!   nn_itend    =   345600  !  10 ans  - 15 min - AM98 1/4deg + P.Marchand
+!   nn_itend    =   172800  !  10 ans  - 30 min - AM98 1/4deg
 !   nn_itend    =    34560  !  10 ans  - 2h
 !   nn_itend    =    86400  !  10 ans  - 1h
@@ -43,7 +42,11 @@
 &namusr_def    !   AM98 user defined namelist  
 !-----------------------------------------------------------------------
-   nn_AM98     =     4     !  AM98 resolution [1/degrees]
+   nn_AM98     =     4     !  AM98 resolution [1/nn_AM98]
    jpkglo      =     2     !  number of model levels
    rn_theta    =    0.    !  rotation angle fo the grid [deg]
+   !
+   nn_gc       =        2        ! number of ghostcells 
+   rn_domsiz   =  2000000        ! size of the domain (default 2000km)  [m]
+   rn_dx       =   100000        ! gridspacing (default 100km)          [m]
    !
    rn_tau           = 0.20    ! wind stress on the surface    [N/m2]
@@ -53,4 +56,11 @@
    rn_beta = 2.e-11   ! beta-plan coriolis parameter  [1/m.s] 
    rn_f0   = 0.7e-4   ! f-plan coriolis parameter     [1/s] 
+   !
+   nn_dynldf_lap_typ = 1        ! choose type of laplacian (ideally from namelist)
+   !                            !       = 1   divrot    laplacian
+   !                            !       = 2   symmetric laplacian (Griffies&Hallberg 2000)
+   !                            !       = 3   symmetric laplacian (cartesian)
+   ln_dynldf_lap_PM  =.false.   ! if True - apply the P.Marchand boundary condition on the laplacian viscosity
+   !
 /
 !-----------------------------------------------------------------------
@@ -60,7 +70,7 @@
    !
 !   rn_Dt      = 7200.     !  2h    - 1 deg cfg - time step for the dynamics [s]
-!   rn_Dt      = 3600.     !  1h    - 1 deg cfg 
-   rn_Dt      = 1800.     !  30min    - 1/4 deg cfg
-!   rn_Dt      = 900.      !  15min - 1/4 deg cfg
+!   rn_Dt      = 3600.     !  1h     - AM98 1 deg
+   rn_Dt      = 1800.     !  30min    - AM98 1/4deg
+!   rn_Dt      = 900.      !  15min - AM98 1/4deg + P.Marchand
 !   rn_Dt      = 600.      !  10min - 1/4 deg cfg
 !   rn_Dt      = 300.      !  5min - 1/4 deg cfg 
@@ -109,5 +119,6 @@
 &namlbc        !   lateral momentum boundary condition                  (default: NO selection)
 !-----------------------------------------------------------------------
-   rn_shlat    =    0.     !  free slip
+   !                       !  free slip  !   partial slip  !   no slip   ! strong slip
+   rn_shlat    =    0.     !  shlat = 0  !  0 < shlat < 2  !  shlat = 2  !  2 < shlat
 /
 !!======================================================================
@@ -141,17 +152,29 @@
 &namdyn_adv    !   formulation of the momentum advection                (default: NO selection)
 !-----------------------------------------------------------------------
-   ln_dynadv_vec = .true.  !  vector form - 2nd centered scheme
+   ln_dynadv_OFF = .false. !  linear dynamics (no momentum advection)
+   ln_dynadv_vec = .false. !  vector form - 2nd centered scheme
      nn_dynkeg     = 0        ! grad(KE) scheme: =0   C2  ;  =1   Hollingsworth correction
+   ln_dynadv_cen2 = .false. !  flux form - 2nd order centered scheme
+   ln_dynadv_ubs = .false. !  flux form - 3rd order UBS      scheme
 /
 !-----------------------------------------------------------------------
 &namdyn_vor    !   Vorticity / Coriolis scheme                          (default: NO selection)
 !-----------------------------------------------------------------------
-   ln_dynvor_ene = .true.   !  energy conserving scheme
-   ln_dynvor_ens = .false.  !  enstrophy conserving scheme
-   ln_dynvor_een = .false.  !  energy & enstrophy scheme
+   ln_dynvor_ene = .false. !  energy    conserving scheme
+   ln_dynvor_ens = .false. !  enstrophy conserving scheme
+   ln_dynvor_enT = .false. !  energy conserving scheme (T-point)
+   ln_dynvor_eeT = .false. !  energy conserving scheme (een using e3t)
+   ln_dynvor_een = .false. !  energy & enstrophy scheme
       nn_een_e3f = 0          ! =0  e3f = mi(mj(e3t))/4
       !                       ! =1  e3f = mi(mj(e3t))/mi(mj( tmask))
-   ln_dynvor_msk = .false.  !  vorticity multiplied by fmask (=T)
-      !                     !  (f-point vorticity schemes only)
+   ln_dynvor_msk = .false. !  vorticity multiplied by fmask (=T)        
+      !                    !  (f-point vorticity schemes only)
+      !
+   ln_dynvor_ens_adVO     = .false.   ! Alternative Direction
+   ln_dynvor_ens_adKE     = .false.
+   ln_dynvor_ens_adKEVO   = .false.
+   ln_dynvor_ene_adVO     = .false.
+   ln_dynvor_ene_adKE     = .false.
+   ln_dynvor_ene_adKEVO   = .false.
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_hgr.F90
===================================================================
--- /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_hgr.F90	(revision 13004)
+++ /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_hgr.F90	(revision 13005)
@@ -68,6 +68,7 @@
       INTEGER  ::   ji, jj               ! dummy loop indices
       REAL(wp) ::   zlam1, zlam0, zcos_theta, zim1 , zjm1 , ze1  , ze1deg ! local scalars
-      REAL(wp) ::   zphi1, zphi0, zsin_theta, zim05, zjm05, znorme      !   -      -
-      REAL(wp) ::   zl, zgl, zbl      !   -      -
+      REAL(wp) ::   zphi1, zphi0, zsin_theta, zim05, zjm05, znorme        !   -      -
+      REAL(wp) ::   zgl, zbl      !   -      -
+
       !!-------------------------------------------------------------------------------
       !
@@ -81,20 +82,6 @@
       !                       !==  grid point position  ==!
       !
-      ze1 =  100000._wp / REAL(nn_AM98, wp)       ! (100km) gridspacing              [m]
-       zl = 2000000._wp                           ! lentgh side square domain        [m]  
-      zgl = 2000000._wp + 4._wp*ze1               ! length side square + ghost cells [m]
-                                                  !    2000km x 2000km + 2 cells around sides
-      ! 
-      ! biggest L = 3500 km (worst case)
-      !zbl = 3500000._wp                           ! length side bigger domain [m]
-      ! non rotated (0deg) - worst case (centered in the 3500kmx3500km)
-      !zlam0 = 875000._wp    ! m
-      !zphi0 = 875000._wp    ! m
-      ! rotated case (45eg)
-      ! origin (in meters)  of the 2000km x 2000km domain 
-      !zlam0 = - 1000000._wp            ! - 1000km
-      !zphi0 =   1474873.7341529163_wp  ! 3500*sin(45deg) - 1000km
-
-
+      ze1 =  rn_dx / REAL(nn_AM98, wp)                   ! [m] gridspacing used
+      zgl =  rn_domsiz + 2._wp * REAL(nn_gc, wp) * ze1   ! [m] length of the square with ghostcells
       ! fit the best square around the square + ghost cells 
       zbl = zgl * ( COS( rn_theta * rad ) + SIN( rn_theta * rad ) )   ! length side bigger domain [m]                                 
@@ -111,10 +98,10 @@
 
       ! exact origin in meters
-      zlam1 =  zbl * COS((rn_theta + 45 )* rad ) / SQRT( 2._wp )  - zl/2._wp 
-      zphi1 =  zbl * SIN((rn_theta + 45 )* rad ) / SQRT( 2._wp )  - zl/2._wp 
+      zlam1 =  zbl * COS((rn_theta + 45 )* rad ) / SQRT( 2._wp )  - rn_domsiz/2._wp 
+      zphi1 =  zbl * SIN((rn_theta + 45 )* rad ) / SQRT( 2._wp )  - rn_domsiz/2._wp 
       ! origin put in the true corner of a cell so there will be no cropping 
       ! of the edge cells
-      zlam0 = REAL( anint( zlam1 / ze1 ) + 1, wp ) * ze1
-      zphi0 = REAl( anint( zphi1 / ze1 ) + 1, wp ) * ze1
+      zlam0 = REAL( anint( zlam1 / ze1 ), wp ) * ze1
+      zphi0 = REAl( anint( zphi1 / ze1 ), wp ) * ze1
 
       IF(lwp) WRITE(numout,*) '                  origin position    zlam0   = ', zlam0/1000,   ' km'
Index: /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_nam.F90
===================================================================
--- /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_nam.F90	(revision 13004)
+++ /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_nam.F90	(revision 13005)
@@ -28,9 +28,18 @@
    INTEGER , PUBLIC ::   nn_AM98            ! 1/nn_AM98 = the resolution chosen in degrees and thus defining the horizontal domain size
    REAL(wp), PUBLIC ::   rn_theta           ! rotation angle (in degree) of the grid 
-   REAL(wp), PUBLIC ::   rn_tau             ! wind stress on the surface    [N/m2]
-   REAL(wp), PUBLIC ::   rn_f0              !    f-plan coriolis parameter  [1/s] 
-   REAL(wp), PUBLIC ::   rn_beta            ! beta-plan coriolis parameter  [1/m.s] 
-   REAL(wp), PUBLIC ::   rn_modified_grav   ! modified gravity              [m/s2] 
-   REAL(wp), PUBLIC ::   rn_rfr             ! layer friction                [1/s]
+   INTEGER , PUBLIC ::   nn_gc              ! number of ghostcells
+   REAL(wp), PUBLIC ::   rn_domsiz          ! size of the domain (default 2000km)  [m]
+   REAL(wp), PUBLIC ::   rn_dx              ! gridspacing (default 100km)          [m]
+   REAL(wp), PUBLIC ::   rn_tau             ! wind stress on the surface           [N/m2]
+   REAL(wp), PUBLIC ::   rn_f0              !    f-plan coriolis parameter         [1/s] 
+   REAL(wp), PUBLIC ::   rn_beta            ! beta-plan coriolis parameter         [1/m.s] 
+   REAL(wp), PUBLIC ::   rn_modified_grav   ! modified gravity                     [m/s2] 
+   REAL(wp), PUBLIC ::   rn_rfr             ! layer friction                       [1/s]
+   !                              !!* temporary *!!
+   INTEGER , PUBLIC ::   nn_dynldf_lap_typ            ! choose type of laplacian (ideally from namelist)
+   !                                                  ! = 1   divrot    laplacian
+   !                                                  ! = 2   symmetric laplacian (Griffies&Hallberg 2000)
+   !                                                  ! = 3   symmetric laplacian (cartesian)
+   LOGICAL , PUBLIC ::   ln_dynldf_lap_PM             ! if True - apply the P.Marchand boundary condition on the laplacian
    !
    !!----------------------------------------------------------------------
@@ -61,7 +70,9 @@
       REAL(wp) ::   ze1, zgl, zbl   ! gridspacing, length of the biggest square
       !!
-      NAMELIST/namusr_def/ nn_AM98, rn_theta, jpkglo, rn_tau,   &   ! 
+      NAMELIST/namusr_def/ nn_AM98, rn_theta, jpkglo,           &   ! 
+         &                 nn_gc ,rn_domsiz, rn_dx,             &   ! domain parameters
          &                 rn_f0 ,rn_beta,                      &   ! coriolis parameter
-         &                 rn_modified_grav, rn_rfr                 !
+         &                 rn_modified_grav, rn_rfr , rn_tau,   &   ! reduced gravity, friction, wind
+         &                 nn_dynldf_lap_typ, ln_dynldf_lap_PM      ! temporary parameter
       !!----------------------------------------------------------------------
       !
@@ -78,16 +89,10 @@
 
       kk_cfg = nn_AM98
-      
-      ! square domain of 2000km x 2000km with dx = dy = 100km 
-      ! kpi = 20 * nn_AM98 + 2        ! Global Domain size
-      ! kpj = 20 * nn_AM98 + 2
-     
-      ! Number of cells
-      ! lenght of the biggest square + 2 ghost cells on each side 
-      ze1 =  100000._wp / REAL(nn_AM98, wp)       ! (100km) gridspacing [m]
-      zgl = 2000000._wp + 4._wp*ze1               ! length side square + ghostcells [m]
-      !zbl = zgl / COS( rn_theta * rad )           ! length side bigger domain [m]  
+       !
+      ze1 =  rn_dx / REAL(nn_AM98, wp)                   ! [m] gridspacing used
+      zgl =  rn_domsiz + 2._wp * REAL(nn_gc, wp) * ze1   ! [m] length of the square with ghostcells
+      ! rotation
       zbl = zgl * ( COS( rn_theta * rad ) + SIN( rn_theta * rad ) )   ! length side bigger domain [m]  
-        
+      !
       kpi = ceiling(zbl / ze1 )    ! Global Domain size + ghost cells                
       kpj = ceiling(zbl / ze1 )    ! Global Domain size + ghost cells              
@@ -104,4 +109,8 @@
       IF( rn_modified_grav /= 0._wp) grav = rn_modified_grav   ! update the gravity 
       !
+      !        !== Temporary namelist parameter ==!
+      !
+      ! ll_dynldf_lap_PM = ln_dynldf_lap_PM
+      !
       kpk = jpkglo
       !                             ! Set the lateral boundary condition of the global domain
@@ -114,10 +123,17 @@
          WRITE(numout,*) '~~~~~~~~~~~ '
          WRITE(numout,*) '   Namelist namusr_def : AM98 case'
+         WRITE(numout,*) '                                   domain size       rn_domsiz   = ', rn_domsiz, 'm'
+         WRITE(numout,*) '                                   gridspacing           rn_dx   = ', rn_dx, 'm'
          WRITE(numout,*) '      inverse resolution & implied domain size         nn_AM98   = ', nn_AM98
+         WRITE(numout,*) '                           implied gridspacing           rn_dx   = ', rn_dx, 'm'
+         WRITE(numout,*) '                          number of ghostcells           nn_gc   = ', nn_gc
+         WRITE(numout,*) '   '
+         WRITE(numout,*) '                    rotation angle chosen             rn_theta   = ', rn_theta, 'deg'
+         WRITE(numout,*) '                    modified gravity          rn_modified_grav   = ', rn_modified_grav, 'm2/s'
 #if defined key_agrif
          IF( Agrif_Root() ) THEN
 #endif
-         WRITE(numout,*) '         jpiglo = Nx*nn_AM98+4                            jpiglo = ', kpi
-         WRITE(numout,*) '         jpjglo = Ny*nn_AM98+4                            jpjglo = ', kpj
+         WRITE(numout,*) '         jpiglo = Nx*nn_AM98 + 2*n_ghost                    jpiglo = ', kpi
+         WRITE(numout,*) '         jpjglo = Nx*nn_AM98 + 2*n_ghost                    jpjglo = ', kpj
 #if defined key_agrif
          ENDIF
Index: /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_sbc.F90
===================================================================
--- /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_sbc.F90	(revision 13004)
+++ /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_sbc.F90	(revision 13005)
@@ -92,5 +92,4 @@
       END DO
      
-     
       ! module of wind stress and wind speed at T-point
       zcoef = 1. / ( zrhoa * zcdrag ) 
Index: /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_zgr.F90
===================================================================
--- /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_zgr.F90	(revision 13004)
+++ /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_zgr.F90	(revision 13005)
@@ -3,5 +3,5 @@
    !!                       ***  MODULE  usrdef_zgr  ***
    !!
-   !!                       ===  GYRE configuration  ===
+   !!                       ===  AM98 configuration  ===
    !!
    !! User defined : vertical coordinate system of a user configuration
@@ -64,5 +64,5 @@
       !
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'usr_def_zgr : GYRE configuration (z-coordinate closed flat box ocean without cavities)'
+      IF(lwp) WRITE(numout,*) 'usr_def_zgr : AM98 configuration (z-coordinate closed flat box ocean without cavities)'
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
       !
@@ -70,5 +70,5 @@
       ! type of vertical coordinate
       ! ---------------------------
-      ld_zco    = .FALSE.         ! GYRE case:  z-coordinate without ocean cavities
+      ld_zco    = .FALSE.         ! AM98 case:  z-coordinate without ocean cavities
       ld_zps    = .FALSE.
       ld_sco    = .TRUE.
@@ -170,5 +170,5 @@
       !! ** Purpose :   set the masked top and bottom ocean t-levels
       !!
-      !! ** Method  :   GYRE case = closed flat box ocean without ocean cavities
+      !! ** Method  :   AM98 case = closed flat box ocean without ocean cavities
       !!                   k_top = 1     except along north, south, east and west boundaries
       !!                   k_bot = jpk-1 except along north, south, east and west boundaries
@@ -182,4 +182,5 @@
       INTEGER  ::   ji, jj                    ! dummy loop indices
       REAL(wp) ::   zylim0, zylim1, zxlim0, zxlim1 ! limit of the domain [m]
+      REAL(WP) ::   zcoeff    ! local scalar
       !!----------------------------------------------------------------------
       !
@@ -187,5 +188,5 @@
       IF(lwp) WRITE(numout,*) '    zgr_top_bot : defines the top and bottom wet ocean levels.'
       IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
-      IF(lwp) WRITE(numout,*) '       GYRE case : closed flat box ocean without ocean cavities'
+      IF(lwp) WRITE(numout,*) '       AM98 case : closed flat box ocean without ocean cavities'
       !
       z2d(:,:) = REAL( jpkm1 , wp )          ! flat bottom
@@ -193,5 +194,4 @@
       CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )           ! set surrounding land to zero (here jperio=0 ==>> closed)
       !
-
       ! rotated domain (45deg)
       !zylim0 =  100000._wp * SQRT( 2._wp ) / 3    !  dy*sqrt(2)/2 * 2/3 
@@ -214,5 +214,4 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-         
          ! if T point in the 2000 km x 2000 km domain
          ! IF ( gphit(ji,jj) > zylim0 .AND. gphit(ji,jj) < zylim1 .AND. & 
@@ -227,8 +226,17 @@
          k_bot(ji,jj) = 0
          END IF
-
          END DO
       END DO
-
+      ! mask the lonely corners
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+         zcoeff = k_top(ji+1,jj) + k_top(ji,jj+1)   &
+            +     k_top(ji-1,jj) + k_top(ji,jj-1)
+         IF ( zcoeff <= 1._wp )   THEN
+            k_top(ji,jj) = 0    ! = land
+            k_bot(ji,jj) = 0
+         END IF
+         END DO
+      END DO
       ! k_bot(:,:) = NINT( z2d(:,:) )           ! =jpkm1 over the ocean point, =0 elsewhere
       !
Index: /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domvvl.F90
===================================================================
--- /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domvvl.F90	(revision 13004)
+++ /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domvvl.F90	(revision 13005)
@@ -723,4 +723,13 @@
                &                    + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj)
          END_2D
+!!an   dans le cas tourné, hf augmente et trend VOR diminue
+!         DO_2D_10_10
+!            zc3(ji,jj) =           (  e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )  &
+!               &                    + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )  &
+!               &                    + e1e2t(ji  ,jj+1) * pssh(ji  ,jj+1)  &
+!               &                    + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj) &
+!               &           / MAX( tmask(ji,jj) + tmask(ji+1,jj) + tmask(ji,jj+1) + tmask(ji+1,jj+1), 1._wp ) 
+!         END_2D
+         
          CALL lbc_lnk( 'domvvl', zc3(:,:), 'F', 1._wp )
          !
Index: /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynadv.F90
===================================================================
--- /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynadv.F90	(revision 13004)
+++ /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynadv.F90	(revision 13005)
@@ -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 13004)
+++ /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynkeg.F90	(revision 13005)
@@ -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/dynldf_lap_blp.F90
===================================================================
--- /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynldf_lap_blp.F90	(revision 13004)
+++ /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynldf_lap_blp.F90	(revision 13005)
@@ -20,5 +20,8 @@
    USE in_out_manager ! I/O manager
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-
+   USE lib_mpp        ! MPP library
+   !
+   USE usrdef_nam , ONLY : nn_dynldf_lap_typ    ! use laplacian parameter
+   !
    IMPLICIT NONE
    PRIVATE
@@ -31,5 +34,5 @@
    INTEGER, PUBLIC, PARAMETER ::   np_dynldf_lap_symN = 3         ! symmetric laplacian (cartesian)
    
-   INTEGER, PUBLIC, PARAMETER ::   ln_dynldf_lap_typ = 1         ! choose type of laplacian (ideally from namelist)
+   !INTEGER, PUBLIC, PARAMETER ::   nn_dynldf_lap_typ = 1         ! choose type of laplacian (ideally from namelist)
 !!anSYM
    !! * Substitutions
@@ -77,9 +80,9 @@
          WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass
          WRITE(numout,*) '~~~~~~~ '
-         WRITE(numout,*) '                  ln_dynldf_lap_typ = ', ln_dynldf_lap_typ
-         SELECT CASE( ln_dynldf_lap_typ )             ! print the choice of operator
+         WRITE(numout,*) '                  nn_dynldf_lap_typ = ', nn_dynldf_lap_typ
+         SELECT CASE( nn_dynldf_lap_typ )             ! print the choice of operator
          CASE( np_dynldf_lap_rot )   ;   WRITE(numout,*) '   ==>>>   div-rot   laplacian'
          CASE( np_dynldf_lap_sym )   ;   WRITE(numout,*) '   ==>>>   symmetric laplacian (covariant form)'
-         CASE( np_dynldf_lap_symN)   ;   WRITE(numout,*) '   ==>>>   symmetric laplacian (simple form)'
+         CASE( np_dynldf_lap_symN)   ;   WRITE(numout,*) '   ==>>>   symmetric laplacian (cartesian form)'
          END SELECT
       ENDIF
@@ -89,5 +92,5 @@
       ENDIF
       !
-      SELECT CASE( ln_dynldf_lap_typ )  
+      SELECT CASE( nn_dynldf_lap_typ )  
          !              
          CASE ( np_dynldf_lap_rot )       !==  Vorticity-Divergence form  ==!
@@ -157,5 +160,5 @@
             END DO                                           !   End of slab
             !
-         CASE ( np_dynldf_lap_symN )       !==  Symmetric form  ==!   (naive way)
+         CASE ( np_dynldf_lap_symN )       !==  Symmetric form  ==!   (cartesian way)
             !
             DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -190,5 +193,5 @@
             !  
          CASE DEFAULT                                     ! error
-            CALL ctl_stop('STOP','dyn_ldf_lap: wrong value for ln_dynldf_lap_typ'  )
+            CALL ctl_stop('STOP','dyn_ldf_lap: wrong value for nn_dynldf_lap_typ'  )
          END SELECT
          !
Index: /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90
===================================================================
--- /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90	(revision 13004)
+++ /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90	(revision 13005)
@@ -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)
@@ -97,5 +117,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs )
+   SUBROUTINE dyn_vor( kt, Kbb, Kmm, puu, pvv, Krhs)
       !!----------------------------------------------------------------------
       !!
@@ -107,9 +127,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
       !!----------------------------------------------------------------------
       !
@@ -165,8 +187,158 @@
             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
+            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_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 )
+                             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
@@ -313,9 +485,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
@@ -371,18 +543,51 @@
          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 ==!
+            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  )
+               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 ==!
+            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
+            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
@@ -411,9 +616,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
@@ -469,18 +674,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
@@ -762,5 +1000,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
       !!----------------------------------------------------------------------
       !
@@ -779,12 +1019,12 @@
       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
+         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
 
@@ -799,5 +1039,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
          !
@@ -809,5 +1049,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
@@ -826,5 +1072,5 @@
       CASE( np_VEC_c2  )
          IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity' 
-         nrvm = np_RVO        ! relative vorticity
+         nrvm = np_RVO        ! relative vorticity      
          ntot = np_CRV        ! relative + planetary vorticity         
       CASE( np_FLX_c2 , np_FLX_ubs  )
@@ -856,10 +1102,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/ldfdyn.F90
===================================================================
--- /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/ldfdyn.F90	(revision 13004)
+++ /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/ldfdyn.F90	(revision 13005)
@@ -25,5 +25,7 @@
    USE lib_mpp         ! distribued memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-
+   !
+   USE usrdef_nam , ONLY : ln_dynldf_lap_PM
+   !
    IMPLICIT NONE
    PRIVATE
@@ -60,5 +62,7 @@
    INTEGER           , PUBLIC ::   nldf_dyn         !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals)
    LOGICAL           , PUBLIC ::   l_ldfdyn_time    !: flag for time variation of the lateral eddy viscosity coef.
-
+!!an
+   !LOGICAL           , PUBLIC ::   ll_dynldf_lap_PM !: flag for P.Marchand modification on viscosity
+!!an
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahmt, ahmf   !: eddy viscosity coef. at T- and F-points [m2/s or m4/s]
    REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dtensq       !: horizontal tension squared         (Smagorinsky only)
@@ -323,33 +327,38 @@
          IF( .NOT.l_ldfdyn_time ) THEN             !* No time variation 
             IF(     ln_dynldf_lap ) THEN                 !   laplacian operator (mask only)
-               ahmt(:,:,1:jpkm1) =       ahmt(:,:,1:jpkm1)   * tmask(:,:,1:jpkm1)
-               WRITE(numout,*) ' ahmt tmask '
+!!an          !
+            WRITE(numout,*) '   ln_dynldf_lap_PM = ',ln_dynldf_lap_PM 
+               IF(     ln_dynldf_lap_PM ) THEN                 !   laplacian operator (mask only)
 !! mask tension at the coast (equivalent of ghostpoints at T)
-!              DO jk = 1, jpk
-!                 DO jj = 1, jpjm1
-!                    DO ji = 1, jpim1      ! NO vector opt.
-!                       ! si sum(fmask)==3 = mouillé (on touche pas)
-!                       ! si sum = 2 alors on met a 0 zsum = fmask + fmask + fmask,.. et si zsum < 2 -> 0 sinon = 1
-!                       zsum =   fmask(ji,jj  ,jk) + fmask(ji+1,jj  ,jk)   &
-!                          &   + fmask(ji,jj+1,jk) + fmask(ji+1,jj+1,jk)
-!                       IF ( zsum < 2._wp )   ahmt(ji,jj,jk) = 0
-!                       !
-!                       !ahmt(ji,jj,jk) = ahmt(ji,jj,jk) * fmask(ji,jj  ,jk) * fmask(ji+1,jj  ,jk)   &
-!                       !   &                            * fmask(ji,jj+1,jk) * fmask(ji+1,jj+1,jk)
-!                    END DO
-!                 END DO
-!              END DO
-!              ahmt(jpi,:,1:jpkm1) =  0._wp
-!              ahmt(:,jpj,1:jpkm1) =  0._wp
-!              WRITE(numout,*) '   an45 ahmt x0'
-
-               ahmf(:,:,1:jpkm1) =       ahmf(:,:,1:jpkm1)   * fmask(:,:,1:jpkm1)
-               WRITE(numout,*) ' ahmf fmask '
-!!an apply no slip at the coast (ssfmask = 1 within the domain and at the coast contrary to fmask in free slip)
-!               DO jk = 1, jpkm1
-!                  ahmf(:,:,jk) =    ahmf(:,:,jk) * ( 2._wp * ssfmask(:,:) - fmask(:,:,jk) )
-!               END DO
-!               WRITE(numout,*) '   an45 ahmf x2'
-
+                  DO jk = 1, jpk
+                     DO jj = 1, jpjm1
+                        DO ji = 1, jpim1      ! NO vector opt.
+                           ! si sum(fmask)==3 = mouillé (on touche pas)
+                           ! si sum = 2 alors on met a 0 zsum = fmask + fmask + fmask,.. et si zsum < 2 -> 0 sinon = 1
+                           zsum =   fmask(ji,jj  ,jk) + fmask(ji+1,jj  ,jk)   &
+                              &   + fmask(ji,jj+1,jk) + fmask(ji+1,jj+1,jk)
+                           IF ( zsum < 2._wp )   ahmt(ji,jj,jk) = 0
+                           !
+                           !ahmt(ji,jj,jk) = ahmt(ji,jj,jk) * fmask(ji,jj  ,jk) * fmask(ji+1,jj  ,jk)   &
+                           !   &                            * fmask(ji,jj+1,jk) * fmask(ji+1,jj+1,jk)
+                        END DO
+                     END DO
+                  END DO
+                  ahmt(jpi,:,1:jpkm1) =  0._wp
+                  ahmt(:,jpj,1:jpkm1) =  0._wp
+                  WRITE(numout,*) '  ahmt x0'
+!! apply no slip at the coast (ssfmask = 1 within the domain and at the coast contrary to fmask in free slip)
+                   DO jk = 1, jpkm1
+                      ahmf(:,:,jk) =    ahmf(:,:,jk) * ( 2._wp * ssfmask(:,:) - fmask(:,:,jk) )
+                   END DO
+                   WRITE(numout,*) '  ahmf x2'
+               ELSE
+               ! classic boundary condition on the viscosity coefficient
+                  ahmt(:,:,1:jpkm1) =       ahmt(:,:,1:jpkm1)   * tmask(:,:,1:jpkm1)
+                  WRITE(numout,*) ' ahmt tmasked '
+                  ahmf(:,:,1:jpkm1) =       ahmf(:,:,1:jpkm1)   * fmask(:,:,1:jpkm1)
+                  WRITE(numout,*) ' ahmf fmasked '
+               ENDIF
+!!an         !                 
             ELSEIF( ln_dynldf_blp ) THEN                 ! bilaplacian operator (square root + mask)
                ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1)
Index: /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90
===================================================================
--- /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90	(revision 13004)
+++ /NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90	(revision 13005)
@@ -13,4 +13,7 @@
    USE phycst           ! physical constants
    USE usrdef_nam
+   USE lib_mpp        ! MPP library
+   USE dynvor , ONLY : ln_dynvor_ens_adVO, ln_dynvor_ens_adKE, ln_dynvor_ens_adKEVO, &
+    &                  ln_dynvor_ene_adVO, ln_dynvor_ene_adKE, ln_dynvor_ene_adKEVO    
    !
    USE iom              ! xIOs server 
@@ -138,9 +141,4 @@
                &         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)
@@ -150,4 +148,25 @@
       END_3D
       !
+      
+!      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 
       z1_2rho0 = 0.5_wp * r1_rho0
