Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/asminc.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/asminc.F90	(revision 13059)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/asminc.F90	(revision 13151)
@@ -95,4 +95,7 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+!!st10
+#  include "domzgr_substitute.h90"
+!!st10
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -417,5 +420,6 @@
                      &            - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * u_bkginc(ji-1,jj,jk)    &
                      &            + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * v_bkginc(ji,jj  ,jk)    &
-                     &            - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk)  ) / e3t(ji,jj,jk,Kmm)
+                     &            - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk)  ) &
+                     &            / e3t(ji,jj,jk,Kmm)
                END_2D
                CALL lbc_lnk( 'asminc', zhdiv, 'T', 1. )   ! lateral boundary cond. (no sign change)
@@ -758,6 +762,10 @@
             !
             ssh(:,:,Kbb) = ssh(:,:,Kmm)                        ! Update before fields
+!!st11
+#if ! defined key_qco
             e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
-!!gm why not e3u(:,:,:,Kbb), e3v(:,:,:,Kbb), gdept(:,:,:,Kbb) ????
+#endif
+!!st11
+!!gm why not e3u(:,:,:,Kbb), e3v(:,:,:,Kbb), gdept(:,:,jk,Kbb) ????
             !
             DEALLOCATE( ssh_bkg    )
@@ -970,6 +978,6 @@
 !           ! set to bottom of a level 
 !                 DO jk = jpk-1, 2, -1
-!                   IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN 
-!                     mld=gdepw(ji,jj,jk+1)
+!                   IF ((mld > gdepw(ji,jj,jk,Kmm)) .and. (mld < gdepw(ji,jj,jk+1,Kmm))) THEN 
+!                     mld=gdepw(ji,jj,jk+1,Kmm)
 !                     jkmax=jk
 !                   ENDIF
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/diawri.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/diawri.F90	(revision 13059)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/diawri.F90	(revision 13151)
@@ -85,4 +85,6 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+!!st12
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -137,11 +139,33 @@
       CALL iom_put("e3v_0", e3v_0(:,:,:) )
       !
-      CALL iom_put( "e3t" , e3t(:,:,:,Kmm) )
-      CALL iom_put( "e3u" , e3u(:,:,:,Kmm) )
-      CALL iom_put( "e3v" , e3v(:,:,:,Kmm) )
-      CALL iom_put( "e3w" , e3w(:,:,:,Kmm) )
-      IF( iom_use("e3tdef") )   &
-         CALL iom_put( "e3tdef"  , ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
-
+!!st13
+#if ! defined key_qco
+      IF ( iom_use("e3t") .OR. iom_use("e3tdef") ) THEN  ! time-varying e3t
+         DO jk = 1, jpk
+            z3d(:,:,jk) =  e3t(:,:,jk,Kmm)
+         END DO
+         CALL iom_put( "e3t"     ,     z3d(:,:,:) )
+         CALL iom_put( "e3tdef"  , ( ( z3d(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
+      ENDIF 
+      IF ( iom_use("e3u") ) THEN                         ! time-varying e3u
+         DO jk = 1, jpk
+            z3d(:,:,jk) =  e3u(:,:,jk,Kmm)
+         END DO 
+         CALL iom_put( "e3u" , z3d(:,:,:) )
+      ENDIF
+      IF ( iom_use("e3v") ) THEN                         ! time-varying e3v
+         DO jk = 1, jpk
+            z3d(:,:,jk) =  e3v(:,:,jk,Kmm)
+         END DO 
+         CALL iom_put( "e3v" , z3d(:,:,:) )
+      ENDIF
+      IF ( iom_use("e3w") ) THEN                         ! time-varying e3w
+         DO jk = 1, jpk
+            z3d(:,:,jk) =  e3w(:,:,jk,Kmm)
+         END DO 
+         CALL iom_put( "e3w" , z3d(:,:,:) )
+      ENDIF
+#endif 
+!!st13
       IF( ll_wd ) THEN
          CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying)
@@ -351,5 +375,6 @@
       !
       REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace
+!!st14
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d, ze3t, zgdept       ! 3D workspace
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl   ! ABL 3D workspace
       !!----------------------------------------------------------------------
@@ -391,5 +416,11 @@
       it = kt
       itmod = kt - nit000 + 1
-
+!!st15
+      ! store e3t for subsitute
+      DO jk = 1, jpk
+         ze3t  (:,:,jk) =  e3t  (:,:,jk,Kmm)
+         zgdept(:,:,jk) =  gdept(:,:,jk,Kmm)
+      END DO
+!!st15
 
       ! 1. Define NETCDF files and fields at beginning of first time step
@@ -514,9 +545,9 @@
             &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
          IF(  .NOT.ln_linssh  ) THEN
-            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t(:,:,:,Kmm)
+            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! ze3t(:,:,:,Kmm)
             &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
-            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t(:,:,:,Kmm)
+            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! ze3t(:,:,:,Kmm)
             &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
-            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t(:,:,:,Kmm)
+            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! ze3t(:,:,:,Kmm)
             &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
          ENDIF
@@ -700,10 +731,11 @@
          WRITE(numout,*) '~~~~~~ '
       ENDIF
-
+!!st16
       IF( .NOT.ln_linssh ) THEN
-         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T  )   ! heat content
-         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T  )   ! salt content
-         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface heat content
-         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity content
+         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T  )   ! heat content
+         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T  )   ! salt content
+         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
+         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
+!!st16
       ELSE
          CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T  )   ! temperature
@@ -713,7 +745,9 @@
       ENDIF
       IF( .NOT.ln_linssh ) THEN
-         zw3d(:,:,:) = ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
-         CALL histwrite( nid_T, "vovvle3t", it, e3t (:,:,:,Kmm) , ndim_T , ndex_T  )   ! level thickness
-         CALL histwrite( nid_T, "vovvldep", it, gdept(:,:,:,Kmm) , ndim_T , ndex_T  )   ! t-point depth
+!!st17 if ! defined key_qco 
+         zw3d(:,:,:) = ( ( ze3t(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
+         CALL histwrite( nid_T, "vovvle3t", it, ze3t (:,:,:)     , ndim_T , ndex_T  )   ! level thickness
+         CALL histwrite( nid_T, "vovvldep", it, zgdept , ndim_T , ndex_T  )   ! t-point depth 
+!!st17
          CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
       ENDIF
@@ -854,4 +888,6 @@
       !!
       INTEGER :: inum, jk
+!!st18  TBR
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, zgdept      ! 3D workspace !!st patch to use substitution
       !!----------------------------------------------------------------------
       ! 
@@ -860,5 +896,10 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
       IF(lwp) WRITE(numout,*) '                and named :', cdfile_name, '...nc'
-
+!!st19 TBR
+      DO jk = 1, jpk
+         ze3t(:,:,jk) =  e3t(:,:,jk,Kmm)
+         zgdept(:,:,jk) =  gdept(:,:,jk,Kmm)
+      END DO
+!!st19
 #if defined key_si3
      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
@@ -878,5 +919,5 @@
       ENDIF
       CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep            )    ! now k-velocity
-      CALL iom_rstput( 0, 0, inum, 'ht'     , ht                 )    ! now water column height
+      CALL iom_rstput( 0, 0, inum, 'ht'     , ht(:,:)            )    ! now water column height
 
       IF ( ln_isf ) THEN
@@ -885,19 +926,19 @@
             CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav    )    ! now k-velocity
             CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav    )    ! now k-velocity
-            CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,8)    )    ! now k-velocity
-            CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,8)    )    ! now k-velocity
-            CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,8), ktype = jp_i1 )
+            CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,wp) )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,wp) )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,wp), ktype = jp_i1 )
          END IF
          IF (ln_isfpar_mlt) THEN
-            CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,8)  )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,wp) )    ! now k-velocity
             CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par          )    ! now k-velocity
             CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par    )    ! now k-velocity
             CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par    )    ! now k-velocity
-            CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,8)    )    ! now k-velocity
-            CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,8)    )    ! now k-velocity
-            CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,8), ktype = jp_i1 )
+            CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,wp) )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,wp) )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,wp), ktype = jp_i1 )
          END IF
       END IF
-
+      !
       IF( ALLOCATED(ahtu) ) THEN
          CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point
@@ -914,7 +955,8 @@
       CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
       CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
-      IF(  .NOT.ln_linssh  ) THEN             
-         CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept(:,:,:,Kmm)        )    !  T-cell depth 
-         CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t(:,:,:,Kmm)          )    !  T-cell thickness  
+!!st20 TBR
+      IF(  .NOT.ln_linssh  ) THEN
+         CALL iom_rstput( 0, 0, inum, 'vovvldep', zgdept        )    !  T-cell depth 
+         CALL iom_rstput( 0, 0, inum, 'vovvle3t', ze3t          )    !  T-cell thickness  
       END IF
       IF( ln_wave .AND. ln_sdw ) THEN
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dom_oce.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dom_oce.F90	(revision 13059)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dom_oce.F90	(revision 13151)
@@ -71,7 +71,7 @@
    !                                !  = 6 cyclic East-West AND North fold F-point pivot
    !                                !  = 7 bi-cyclic East-West AND North-South
-   LOGICAL, PUBLIC ::   l_Iperio, l_Jperio   !   should we explicitely take care I/J periodicity 
-
-   !                                 !  domain MPP decomposition parameters
+   LOGICAL, PUBLIC ::   l_Iperio, l_Jperio   !   should we explicitely take care I/J periodicity
+
+   !                             !: domain MPP decomposition parameters
    INTEGER             , PUBLIC ::   nimpp, njmpp     !: i- & j-indexes for mpp-subdomain left bottom
    INTEGER             , PUBLIC ::   nreci, nrecj     !: overlap region in i and j
@@ -136,30 +136,35 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3vw_0   !: vw-vert. scale factor [m]
    !                                                        !  time-dependent scale factors
+!!st1
+#if ! defined key_qco
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e3t, e3u, e3v, e3w, e3uw, e3vw  !: vert. scale factor [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   e3f                             !: F-point vert. scale factor [m]
+#endif
+   !                                                        !  time-dependent ratio ssh / h_0
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   r3t, r3u, r3v                   !: time-dependent    ratio at t-, u- and v-point [-]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   r3f                             !: mid-time-level    ratio at f-point            [-]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   r3t_f, r3u_f, r3v_f             !: now time-filtered ratio at t-, u- and v-point [-]
 
    !                                                        !  reference depths of cells
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_0  !: t- depth              [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdepw_0  !: w- depth              [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gde3w_0  !: w- depth (sum of e3w) [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   gdept_0  !: t- depth              [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   gdepw_0  !: w- depth              [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   gde3w_0  !: w- depth (sum of e3w) [m]
    !                                                        !  time-dependent depths of cells
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  gdept, gdepw  
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  gde3w  
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   gdept, gdepw
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   gde3w
+!!st2
+   !                                                        !  reference heights of ocean water column and its inverse
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   ht_0, r1_ht_0   !: t-depth        [m] and [1/m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   hu_0, r1_hu_0   !: u-depth        [m] and [1/m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   hv_0, r1_hv_0   !: v-depth        [m] and [1/m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   hf_0, r1_hf_0   !: f-depth        [m] and [1/m]
+   !                                                        ! time-dependent heights of ocean water column
+#if ! defined key_qco
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   ht          !: t-points           [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hu, r1_hu   !: u-depth            [m] and [1/m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hv, r1_hv   !: v-depth            [m] and [1/m]
+#endif
+!!st2
    
-   !                                                      !  reference heights of water column
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0  !: t-depth              [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0  !: u-depth              [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hv_0  !: v-depth              [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hf_0  !: f-depth              [m]
-   !                                                      !  inverse of reference heights of water column
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   r1_ht_0  !: t-depth              [1/m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   r1_hu_0  !: u-depth              [1/m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   r1_hv_0  !: v-depth              [1/m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   r1_hf_0  !: f-depth              [1/m]
-   
-                                                          ! time-dependent heights of water column
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ht                     !: height of water column at T-points [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hu, hv, r1_hu, r1_hv   !: height of water column [m] and reciprocal [1/m]
-
    INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1)
    INTEGER, PUBLIC ::   nlb10              !: shallowest W level Bellow ~10m (nla10 + 1) 
@@ -176,6 +181,6 @@
    !! ---------------------------------------------------------------------
 !!gm Proposition of new name for top/bottom vertical indices
-!   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mtk_t, mtk_u, mtk_v   !: top first wet T-, U-, V-, F-level (ISF)
-!   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbk_t, mbk_u, mbk_v   !: bottom last wet T-, U- and V-level
+!   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mtk_t, mtk_u, mtk_v   !: top    first wet T-, U-, and V-level (ISF)
+!   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbk_t, mbk_u, mbk_v   !: bottom last  wet T-, U-, and V-level
 !!gm
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbkt, mbku, mbkv   !: bottom last wet T-, U- and V-level
@@ -244,16 +249,20 @@
    END FUNCTION Agrif_CFixed
 #endif
-
+!!st3: dom_oce_alloc modified to ease the ifdef if necessary (gm stuff)
    INTEGER FUNCTION dom_oce_alloc()
       !!----------------------------------------------------------------------
-      INTEGER, DIMENSION(12) :: ierr
+      INTEGER                ::   ii
+      INTEGER, DIMENSION(30) :: ierr
       !!----------------------------------------------------------------------
-      ierr(:) = 0
+      ii = 0   ;   ierr(:) = 0
       !
-      ALLOCATE( mig(jpi), mjg(jpj), STAT=ierr(1) )
-         !
-      ALLOCATE( mi0(jpiglo)   , mi1 (jpiglo),  mj0(jpjglo)   , mj1 (jpjglo) ,     &
-         &      tpol(jpiglo)  , fpol(jpiglo)                                , STAT=ierr(2) )
-         !
+      ii = ii+1
+      ALLOCATE( mig(jpi), mjg(jpj), STAT=ierr(ii) )
+         !
+      ii = ii+1
+      ALLOCATE( mi0 (jpiglo) , mi1 (jpiglo),  mj0(jpjglo) , mj1 (jpjglo) ,     &
+         &      tpol(jpiglo) , fpol(jpiglo)                              , STAT=ierr(ii) )
+         !
+      ii = ii+1
       ALLOCATE( glamt(jpi,jpj) ,    glamu(jpi,jpj) ,  glamv(jpi,jpj) ,  glamf(jpi,jpj) ,     &
          &      gphit(jpi,jpj) ,    gphiu(jpi,jpj) ,  gphiv(jpi,jpj) ,  gphif(jpi,jpj) ,     &
@@ -266,33 +275,56 @@
          &      e1e2v(jpi,jpj) , r1_e1e2v(jpi,jpj) , e1_e2v(jpi,jpj)                   ,     &
          &      e1e2f(jpi,jpj) , r1_e1e2f(jpi,jpj)                                     ,     &
-         &      ff_f (jpi,jpj) ,    ff_t (jpi,jpj)                                     , STAT=ierr(3) )
-         !
+         &      ff_f (jpi,jpj) ,    ff_t (jpi,jpj)                                     , STAT=ierr(ii) )
+         !
+      ii = ii+1
+      ALLOCATE(  e3t_0(jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) , e3v_0 (jpi,jpj,jpk) , e3f_0(jpi,jpj,jpk) ,      &
+         &       e3w_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk)                      ,  STAT=ierr(ii) )
+         !
+      ii = ii+1
       ALLOCATE( gdept_0(jpi,jpj,jpk)     , gdepw_0(jpi,jpj,jpk)     , gde3w_0(jpi,jpj,jpk) ,      &
-         &      gdept  (jpi,jpj,jpk,jpt) , gdepw  (jpi,jpj,jpk,jpt) , gde3w  (jpi,jpj,jpk) , STAT=ierr(4) )
-         !
-      ALLOCATE( e3t_0(jpi,jpj,jpk)     , e3u_0(jpi,jpj,jpk)     , e3v_0(jpi,jpj,jpk)     , e3f_0(jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk)     ,   &
-         &      e3t  (jpi,jpj,jpk,jpt) , e3u  (jpi,jpj,jpk,jpt) , e3v  (jpi,jpj,jpk,jpt) , e3f  (jpi,jpj,jpk) , e3w  (jpi,jpj,jpk,jpt) ,   & 
-         &      e3uw_0(jpi,jpj,jpk)     , e3vw_0(jpi,jpj,jpk)     ,         &
-         &      e3uw  (jpi,jpj,jpk,jpt) , e3vw  (jpi,jpj,jpk,jpt) ,    STAT=ierr(5) )                       
-         !    
-      ALLOCATE(    ht_0(jpi,jpj) ,    hu_0(jpi,jpj)     ,    hv_0(jpi,jpj)     ,    hf_0(jpi,jpj) ,     &
-         &      r1_ht_0(jpi,jpj) , r1_hu_0(jpi,jpj)     , r1_hv_0(jpi,jpj)     , r1_hf_0(jpi,jpj) ,     &
-         &         ht  (jpi,jpj) ,    hu  (jpi,jpj,jpt) ,    hv  (jpi,jpj,jpt)                    ,     &
-         &                         r1_hu  (jpi,jpj,jpt) , r1_hv  (jpi,jpj,jpt)                    , STAT=ierr(6)  )
-         !
-      ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(7)  ) 
-         !
-      ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(8) )
-         !
-      ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                                           & 
+         &      gdept  (jpi,jpj,jpk,jpt) , gdepw  (jpi,jpj,jpk,jpt) , gde3w  (jpi,jpj,jpk) , STAT=ierr(ii) )
+         !
+!!st4
+#if ! defined key_qco
+      ii = ii+1
+      ALLOCATE( e3t(jpi,jpj,jpk,jpt) , e3u (jpi,jpj,jpk,jpt) , e3v (jpi,jpj,jpk,jpt) , e3f(jpi,jpj,jpk) ,      &
+         &      e3w(jpi,jpj,jpk,jpt) , e3uw(jpi,jpj,jpk,jpt) , e3vw(jpi,jpj,jpk,jpt)                    ,  STAT=ierr(ii) )
+#endif
+!!st4
+         !
+      ii = ii+1
+      ALLOCATE( r3t  (jpi,jpj,jpt)   , r3u  (jpi,jpj,jpt)    , r3v  (jpi,jpj,jpt)    , r3f  (jpi,jpj) ,  &
+         &      r3t_f(jpi,jpj)       , r3u_f(jpi,jpj)        , r3v_f(jpi,jpj)                         ,  STAT=ierr(ii) )       
+         !
+      ii = ii+1
+      ALLOCATE( ht_0(jpi,jpj) ,    hu_0(jpi,jpj)    ,    hv_0(jpi,jpj)     , hf_0(jpi,jpj) ,       &
+         &   r1_ht_0(jpi,jpj) , r1_hu_0(jpi,jpj) ,    r1_hv_0(jpi,jpj),   r1_hf_0(jpi,jpj) ,   STAT=ierr(ii)  )
+         !
+#if ! defined key_qco
+      ii = ii+1
+      ALLOCATE( ht  (jpi,jpj) ,    hu  (jpi,jpj,jpt),    hv  (jpi,jpj,jpt)                 ,       &
+         &                      r1_hu  (jpi,jpj,jpt), r1_hv  (jpi,jpj,jpt)                 ,   STAT=ierr(ii)  )
+#endif
+         !
+      ii = ii+1
+      ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(ii)  ) 
+         !
+      ii = ii+1
+      ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(ii) )
+         !
+      ii = ii+1
+      ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                                           &
          &      ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) ,     &
-         &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj)                    , STAT=ierr(9) )
-         !
-      ALLOCATE( mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj), mikf(jpi,jpj), STAT=ierr(10) )
-         !
-      ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) ,     & 
-         &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk) , STAT=ierr(11) )
-         !
-      ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(12) )
+         &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) ,                    STAT=ierr(ii) )
+         !
+      ii = ii+1
+      ALLOCATE( mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj), mikf(jpi,jpj), STAT=ierr(ii) )
+         !
+      ii = ii+1
+      ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) ,     &
+         &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk) , STAT=ierr(ii) )
+         !
+      ii = ii+1
+      ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(ii) )
       !
       dom_oce_alloc = MAXVAL(ierr)
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domain.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domain.F90	(revision 13059)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domain.F90	(revision 13151)
@@ -34,5 +34,11 @@
    USE dommsk         ! domain: set the mask system
    USE domwri         ! domain: write the meshmask file
+!!st5
+#if ! defined key_qco
    USE domvvl         ! variable volume
+#else
+   USE domqco          ! variable volume
+#endif
+!!st5
    USE c1d            ! 1D configuration
    USE dyncor_c1d     ! 1D configuration: Coriolis term    (cor_c1d routine)
@@ -78,5 +84,7 @@
       CHARACTER (len=*), INTENT(in) :: cdstr                  ! model: NEMO or SAS. Determines core restart variables
       !
-      INTEGER ::   ji, jj, jk, ik   ! dummy loop indices
+!!st6
+      INTEGER ::   ji, jj, jk, jt   ! dummy loop indices
+!!st6
       INTEGER ::   iconf = 0    ! local integers
       CHARACTER (len=64) ::   cform = "(A12, 3(A13, I7))" 
@@ -114,5 +122,5 @@
          CASE( 7 )   ;   WRITE(numout,*) '         (i.e. cyclic east-west and north-south)'
          CASE DEFAULT
-            CALL ctl_stop( 'jperio is out of range' )
+            CALL ctl_stop( 'dom_init:   jperio is out of range' )
          END SELECT
          WRITE(numout,*)     '      Ocean model configuration used:'
@@ -144,9 +152,9 @@
       IF( ln_closea ) CALL dom_clo      ! Read in masks to define closed seas and lakes
 
-      CALL dom_zgr( ik_top, ik_bot )    ! Vertical mesh and bathymetry
+      CALL dom_zgr( ik_top, ik_bot )    ! Vertical mesh and bathymetry (return top and bottom ocean t-level indices)
 
       CALL dom_msk( ik_top, ik_bot )    ! Masks
-
-      ht_0(:,:) = 0._wp                 ! Reference ocean thickness
+      !
+      ht_0(:,:) = 0._wp  ! Reference ocean thickness
       hu_0(:,:) = 0._wp
       hv_0(:,:) = 0._wp
@@ -190,37 +198,48 @@
 !            r1_e1e2t(ji,jj) = r1_e1e2t(ji,jj) / zcoeff
 !!an45
+!!st7 : make it easier to use key_qco condition (gm stuff)
+#if defined key_qco
+      !           !==  initialisation of time varying coordinate  ==!   Quasi-Euerian coordinate case
+      !
+      IF( .NOT.l_offline )   CALL dom_qco_init( Kbb, Kmm, Kaa )
+      !
+      IF( ln_linssh )        CALL ctl_stop('STOP','domain: key_qco and ln_linssh = T are incompatible')
+      !
+#else
       !           !==  time varying part of coordinate system  ==!
       !
       IF( ln_linssh ) THEN       != Fix in time : set to the reference one for all
-      !
-         !       before        !          now          !       after         !
-            gdept(:,:,:,Kbb) = gdept_0  ;   gdept(:,:,:,Kmm) = gdept_0   ;   gdept(:,:,:,Kaa) = gdept_0   ! depth of grid-points
-            gdepw(:,:,:,Kbb) = gdepw_0  ;   gdepw(:,:,:,Kmm) = gdepw_0   ;   gdepw(:,:,:,Kaa) = gdepw_0   !
-                                   gde3w = gde3w_0   !        ---          !
-         !                                                                  
-              e3t(:,:,:,Kbb) =   e3t_0  ;     e3t(:,:,:,Kmm) =   e3t_0   ;   e3t(:,:,:,Kaa) =  e3t_0    ! scale factors
-              e3u(:,:,:,Kbb) =   e3u_0  ;     e3u(:,:,:,Kmm) =   e3u_0   ;   e3u(:,:,:,Kaa) =  e3u_0    !
-              e3v(:,:,:,Kbb) =   e3v_0  ;     e3v(:,:,:,Kmm) =   e3v_0   ;   e3v(:,:,:,Kaa) =  e3v_0    !
-                                     e3f =   e3f_0   !        ---          !
-              e3w(:,:,:,Kbb) =   e3w_0  ;     e3w(:,:,:,Kmm) =   e3w_0   ;    e3w(:,:,:,Kaa) =   e3w_0   ! 
-             e3uw(:,:,:,Kbb) =  e3uw_0  ;    e3uw(:,:,:,Kmm) =  e3uw_0   ;   e3uw(:,:,:,Kaa) =  e3uw_0   !  
-             e3vw(:,:,:,Kbb) =  e3vw_0  ;    e3vw(:,:,:,Kmm) =  e3vw_0   ;   e3vw(:,:,:,Kaa) =  e3vw_0   !
-         !
-         z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )     ! _i mask due to ISF
-         z1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) )
-         !
-         !        before       !          now          !       after         !
-                                      ht =    ht_0   !                     ! water column thickness
-               hu(:,:,Kbb) =    hu_0  ;      hu(:,:,Kmm) =    hu_0   ;    hu(:,:,Kaa) =    hu_0   ! 
-               hv(:,:,Kbb) =    hv_0  ;      hv(:,:,Kmm) =    hv_0   ;    hv(:,:,Kaa) =    hv_0   !
-            r1_hu(:,:,Kbb) = z1_hu_0  ;   r1_hu(:,:,Kmm) = z1_hu_0   ; r1_hu(:,:,Kaa) = z1_hu_0   ! inverse of water column thickness
-            r1_hv(:,:,Kbb) = z1_hv_0  ;   r1_hv(:,:,Kmm) = z1_hv_0   ; r1_hv(:,:,Kaa) = z1_hv_0   !
-         !
+         !
+         DO jt = 1, jpt                         ! depth of t- and w-grid-points
+            gdept(:,:,:,jt) = gdept_0(:,:,:)
+            gdepw(:,:,:,jt) = gdepw_0(:,:,:)
+         END DO
+            gde3w(:,:,:)    = gde3w_0(:,:,:)    ! = gdept as the sum of e3t
+         !
+         DO jt = 1, jpt                         ! vertical scale factors
+            e3t(:,:,:,jt) =  e3t_0(:,:,:)
+            e3u(:,:,:,jt) =  e3u_0(:,:,:)
+            e3v(:,:,:,jt) =  e3v_0(:,:,:)
+            e3w(:,:,:,jt) =  e3w_0(:,:,:)
+            e3uw(:,:,:,jt) = e3uw_0(:,:,:)
+            e3vw(:,:,:,jt) = e3vw_0(:,:,:)
+         END DO
+            e3f(:,:,:)    =  e3f_0(:,:,:)
+         !
+         DO jt = 1, jpt                         ! water column thickness and its inverse
+            hu(:,:,jt)    =    hu_0(:,:)
+            hv(:,:,jt)    =    hv_0(:,:)
+            r1_hu(:,:,jt) = r1_hu_0(:,:)
+            r1_hv(:,:,jt) = r1_hv_0(:,:)
+         END DO
+            ht(:,:) =    ht_0(:,:)
          !
       ELSE                       != time varying : initialize before/now/after variables
          !
-         IF( .NOT.l_offline )  CALL dom_vvl_init( Kbb, Kmm, Kaa )
-         !
-      ENDIF
+         IF( .NOT.l_offline )   CALL dom_vvl_init( Kbb, Kmm, Kaa )
+         !
+      ENDIF
+#endif
+!!st7
       !
       IF( lk_c1d         )   CALL cor_c1d       ! 1D configuration: Coriolis set at T-point
@@ -238,5 +257,4 @@
          WRITE(numout,*) 
       ENDIF
-      
       !
    END SUBROUTINE dom_init
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domvvl.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domvvl.F90	(revision 13059)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domvvl.F90	(revision 13151)
@@ -1,2 +1,3 @@
+
 MODULE domvvl
    !!======================================================================
@@ -11,12 +12,4 @@
    !!----------------------------------------------------------------------
 
-   !!----------------------------------------------------------------------
-   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
-   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
-   !!   dom_vvl_sf_update   : Swap vertical scale factors and update the vertical grid
-   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
-   !!   dom_vvl_rst      : read/write restart file
-   !!   dom_vvl_ctl      : Check the vvl options
-   !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
    USE phycst          ! physical constant
@@ -36,10 +29,4 @@
    PRIVATE
 
-   PUBLIC  dom_vvl_init       ! called by domain.F90
-   PUBLIC  dom_vvl_zgr        ! called by isfcpl.F90
-   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90
-   PUBLIC  dom_vvl_sf_update  ! called by step.F90
-   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90
-
    !                                                      !!* Namelist nam_vvl
    LOGICAL , PUBLIC :: ln_vvl_zstar           = .FALSE.    ! zstar  vertical coordinate
@@ -62,5 +49,33 @@
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                 ! retoring period for scale factors
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence
-
+!!stoops
+#if defined key_qco
+   !!----------------------------------------------------------------------
+   !!   'key_qco'      EMPTY MODULE      Quasi-Eulerian vertical coordonate
+   !!----------------------------------------------------------------------
+#else
+   !!----------------------------------------------------------------------
+   !!   Default key      Old management of time varying vertical coordinate
+   !!----------------------------------------------------------------------
+!!st
+   !!----------------------------------------------------------------------
+   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
+   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
+   !!   dom_vvl_sf_update   : Swap vertical scale factors and update the vertical grid
+   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
+   !!   dom_vvl_rst      : read/write restart file
+   !!   dom_vvl_ctl      : Check the vvl options
+   !!----------------------------------------------------------------------
+  
+   PUBLIC  dom_vvl_init       ! called by domain.F90
+   PUBLIC  dom_vvl_zgr        ! called by isfcpl.F90
+   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90
+   PUBLIC  dom_vvl_sf_update  ! called by step.F90
+   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90
+   PUBLIC  dom_vvl_interpol_st! called by dynnxt.F90
+   PUBLIC  dom_vvl_sf_nxt_st  ! called by step.F90
+   PUBLIC  dom_vvl_sf_update_st
+!!st
+   
    !! * Substitutions
 #  include "do_loop_substitute.h90"
@@ -132,5 +147,5 @@
       e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
       !
-      CALL dom_vvl_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column
+      CALL dom_vvl_zgr_st(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column
       !
    END SUBROUTINE dom_vvl_init
@@ -290,4 +305,82 @@
    END SUBROUTINE dom_vvl_zgr
 
+   
+   SUBROUTINE dom_vvl_zgr_st(Kbb, Kmm, Kaa)
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_init  ***
+      !!                   
+      !! ** Purpose :  Interpolation of all scale factors, 
+      !!               depths and water column heights
+      !!
+      !! ** Method  :  - interpolate scale factors
+      !!
+      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b)
+      !!              - Regrid: e3(u/v)_n
+      !!                        e3(u/v)_b       
+      !!                        e3w_n           
+      !!                        e3(u/v)w_b      
+      !!                        e3(u/v)w_n      
+      !!                        gdept_n, gdepw_n and gde3w_n
+      !!              - h(t/u/v)_0
+      !!              - frq_rst_e3t and frq_rst_hdv
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
+      !!----------------------------------------------------------------------
+      INTEGER ::   ji, jj, jk
+      INTEGER ::   ii0, ii1, ij0, ij1
+      REAL(wp)::   zcoef
+      !!----------------------------------------------------------------------
+      !
+      !                    !== Set of all other vertical scale factors  ==!  (now and before)
+      !                                ! Horizontal interpolation of e3t
+      CALL dom_vvl_interpol_st( r3u(:,:,Kbb), e3u(:,:,:,Kbb), 'U' )    ! from T to U
+      CALL dom_vvl_interpol_st( r3u(:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
+      CALL dom_vvl_interpol_st( r3v(:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V 
+      CALL dom_vvl_interpol_st( r3v(:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
+      CALL dom_vvl_interpol_st( r3f(:,:), e3f(:,:,:), 'F' )    ! from U to F
+      !                                ! Vertical interpolation of e3t,u,v 
+      CALL dom_vvl_interpol_st( r3t(:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )  ! from T to W
+      CALL dom_vvl_interpol_st( r3t(:,:,Kbb), e3w (:,:,:,Kbb), 'W'  )
+      CALL dom_vvl_interpol_st( r3u(:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )  ! from U to UW
+      CALL dom_vvl_interpol_st( r3u(:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
+      CALL dom_vvl_interpol_st( r3v(:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )  ! from V to UW
+      CALL dom_vvl_interpol_st( r3v(:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
+
+      ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
+      e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm)
+      e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm)
+      e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm)
+      !
+      DO_3D_11_11( 1, jpk )
+         gdepw(ji,jj,jk,Kmm) =  gdepw_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kmm))
+         gdept(ji,jj,jk,Kmm) =  gdept_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kmm)) 
+         gde3w(ji,jj,jk    ) =  gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
+         gdepw(ji,jj,jk,Kbb) =  gdepw_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kbb))
+         gdept(ji,jj,jk,Kbb) =  gdept_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kbb)) 
+      END_3D     
+      !
+      !                    !==  thickness of the water column  !!   (ocean portion only)
+      ht(:,:) = ht_0(:,:) + ssh(:,:,Kmm)
+      hu(:,:,Kbb) = (hu_0(:,:)*(1._wp+r3u(:,:,Kbb)))
+      hv(:,:,Kbb) = (hv_0(:,:)*(1._wp+r3v(:,:,Kbb)))
+      hu(:,:,Kbb) = (hu_0(:,:)*(1._wp+r3u(:,:,Kmm)))
+      hv(:,:,Kbb) = (hv_0(:,:)*(1._wp+r3v(:,:,Kmm)))
+      !                    !==  inverse of water column thickness   ==!   (u- and v- points)
+      r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
+      r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
+      r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) )
+      r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) ) 
+      !
+      IF(lwxios) THEN
+! define variables in restart file when writing with XIOS
+         CALL iom_set_rstw_var_active('e3t_b')
+         CALL iom_set_rstw_var_active('e3t_n')
+         !
+      ENDIF
+      !
+   END SUBROUTINE dom_vvl_zgr_st
+   
 
    SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall ) 
@@ -572,4 +665,95 @@
 
 
+
+   SUBROUTINE dom_vvl_sf_nxt_st( kt, Kbb, Kmm, Kaa, kcall ) 
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
+      !!                   
+      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
+      !!                 tranxt and dynspg routines
+      !!
+      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
+      !!               - z_tilde_case: after scale factor increment = 
+      !!                                    high frequency part of horizontal divergence
+      !!                                  + retsoring towards the background grid
+      !!                                  + thickness difusion
+      !!                               Then repartition of ssh INCREMENT proportionnaly
+      !!                               to the "baroclinic" level thickness.
+      !!
+      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
+      !!               - tilde_e3t_a: after increment of vertical scale factor 
+      !!                              in z_tilde case
+      !!               - e3(t/u/v)_a
+      !!
+      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in )           ::   kt             ! time step
+      INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step
+      INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence
+      !
+      INTEGER                ::   ji, jj, jk            ! dummy loop indices
+      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
+      REAL(wp)               ::   z_tmin, z_tmax        ! local scalars
+      LOGICAL                ::   ll_do_bclinic         ! local logical
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
+      !
+      IF( ln_timing )   CALL timing_start('dom_vvl_sf_nxt')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
+      ENDIF
+
+      ll_do_bclinic = .TRUE.
+      IF( PRESENT(kcall) ) THEN
+         IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
+      ENDIF
+
+      ! ******************************* !
+      ! After acale factors at t-points !
+      ! ******************************* !
+      !
+      DO jk = 1, jpkm1
+         e3t(:,:,jk,Kaa) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kaa) )
+         e3u(:,:,jk,Kaa) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kaa) )
+         e3v(:,:,jk,Kaa) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kaa) )
+      END DO
+      !
+      ! *********************************** !
+      ! After scale factors at u- v- points !
+      ! *********************************** !
+      
+      !!st CALL dom_vvl_interpol_st( r3u(:,:,Kaa), e3u(:,:,:,Kaa), 'U' )
+      !!st CALL dom_vvl_interpol_st( r3v(:,:,Kaa), e3v(:,:,:,Kaa), 'V' )
+
+      ! *********************************** !
+      ! After depths at u- v points         !
+      ! *********************************** !
+
+      !!st hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1)
+      !!st hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1)
+      !!st DO jk = 2, jpkm1
+      !!st    hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk)
+      !!st    hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk)
+      !!st    
+      !!st END DO
+      hu(:,:,Kaa) = (hu_0(:,:)*(1._wp+r3u(:,:,Kaa)))
+      hv(:,:,Kaa) = (hv_0(:,:)*(1._wp+r3v(:,:,Kaa)))
+      !                                        ! Inverse of the local depth
+!!gm BUG ?  don't understand the use of umask_i here .....
+      r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
+      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
+      !
+      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
+      !
+   END SUBROUTINE dom_vvl_sf_nxt_st
+   
+
+
    SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa )
       !!----------------------------------------------------------------------
@@ -672,5 +856,140 @@
       !
    END SUBROUTINE dom_vvl_sf_update
-
+   
+
+   SUBROUTINE dom_vvl_sf_update_st( kt, Kbb, Kmm, Kaa )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_sf_update  ***
+      !!                   
+      !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors 
+      !!               compute all depths and related variables for next time step
+      !!               write outputs and restart file
+      !!
+      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE)
+      !!               - reconstruct scale factor at other grid points (interpolate)
+      !!               - recompute depths and water height fields
+      !!
+      !! ** Action  :  - tilde_e3t_(b/n) ready for next time step
+      !!               - Recompute:
+      !!                    e3(u/v)_b       
+      !!                    e3w(:,:,:,Kmm)           
+      !!                    e3(u/v)w_b      
+      !!                    e3(u/v)w_n      
+      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
+      !!                    h(u/v) and h(u/v)r
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
+      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt              ! time step
+      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zcoef        ! local scalar
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
+      !
+      IF( ln_timing )   CALL timing_start('dom_vvl_sf_update')
+      !
+      IF( kt == nit000 )   THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
+      ENDIF
+      !
+
+      ! Compute all missing vertical scale factor and depths
+      ! ====================================================
+      ! Horizontal scale factor interpolations
+      ! --------------------------------------
+      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
+      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
+      
+      CALL dom_vvl_interpol_st( r3f(:,:), e3f(:,:,:), 'F'  )
+
+      ! Vertical scale factor interpolations
+      CALL dom_vvl_interpol_st( r3t(:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  )
+      CALL dom_vvl_interpol_st( r3u(:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
+      CALL dom_vvl_interpol_st( r3v(:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
+      CALL dom_vvl_interpol_st( r3t(:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  )
+      CALL dom_vvl_interpol_st( r3u(:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
+      CALL dom_vvl_interpol_st( r3v(:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
+
+      ! t- and w- points depth (set the isf depth as it is in the initial step)
+      DO_3D_11_11( 1, jpk )
+         gdepw(ji,jj,jk,Kmm) =  gdepw_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kmm))
+         gdept(ji,jj,jk,Kmm) =  gdept_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kmm)) 
+         gde3w(ji,jj,jk    ) =  gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
+      END_3D
+
+      ! Local depth and Inverse of the local depth of the water
+      ! -------------------------------------------------------
+      !
+      ht(:,:) = ht_0(:,:) + ssh(:,:,Kmm)
+
+      ! write restart file
+      ! ==================
+      IF( lrst_oce  )   CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' )
+      !
+      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_update')
+      !
+   END SUBROUTINE dom_vvl_sf_update_st
+   
+
+    
+   SUBROUTINE dom_vvl_interpol_st( rc3, pe3, cdp )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_vvl__interpol  ***
+      !!
+      !! ** Purpose :   interpolate scale factors from one grid point to another
+      !!
+      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
+      !!                - horizontal interpolation: grid cell surface averaging
+      !!                - vertical interpolation: simple averaging
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   rc3     ! input e3   NOT used here (ssh is used instead)
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe3     ! scale factor e3 to be updated   [m]
+      CHARACTER(LEN=*)          , INTENT(in   ) ::   cdp     ! grid point of the scale factor ( 'U', 'V', 'W, 'F', 'UW' or 'VW' )
+      !
+      INTEGER ::   ji, jj, jk                 ! dummy loop indices
+      REAL(wp), DIMENSION(jpi,jpj) ::   zc3   ! 2D workspace
+      !!----------------------------------------------------------------------
+      !
+      SELECT CASE ( cdp )     !==  type of interpolation  ==!
+         !
+      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
+         DO jk = 1, jpkm1
+            pe3(:,:,jk) = e3u_0(:,:,jk) * ( 1.0_wp + rc3(:,:) )
+         END DO
+         !
+      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
+         DO jk = 1, jpkm1
+            pe3(:,:,jk) = e3v_0(:,:,jk) * ( 1.0_wp + rc3(:,:) )
+         END DO
+         !
+      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
+         DO jk = 1, jpkm1                    ! Horizontal interpolation of e3f from ssh
+            e3f(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + rc3(:,:) )
+         END DO
+         !
+      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
+         DO jk = 1, jpk
+            pe3(:,:,jk) = e3w_0(:,:,jk) * ( 1.0_wp + rc3(:,:) )
+         END DO
+         !
+      CASE( 'UW' )                  !* from U- to UW-point
+         DO jk = 1, jpk
+            pe3(:,:,jk) = e3uw_0(:,:,jk) * ( 1.0_wp + rc3(:,:) )
+         END DO
+      CASE( 'VW' )                  !* from U- to UW-point : vertical simple mean
+         DO jk = 1, jpk
+            pe3(:,:,jk) = e3vw_0(:,:,jk) * ( 1.0_wp + rc3(:,:) )
+         END DO
+         !
+      END SELECT
+      !
+   END SUBROUTINE dom_vvl_interpol_st
+   
 
    SUBROUTINE dom_vvl_interpol( pssh, pe3, cdp )
@@ -723,13 +1042,4 @@
                &                    + 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 )
          !
@@ -771,5 +1081,5 @@
       !
    END SUBROUTINE dom_vvl_interpol
-
+   
 
    SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw )
@@ -1036,4 +1346,7 @@
    END SUBROUTINE dom_vvl_ctl
 
+#endif
+!!stoops
+
    !!======================================================================
 END MODULE domvvl
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domzgr_substitute.h90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domzgr_substitute.h90	(revision 13151)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domzgr_substitute.h90	(revision 13151)
@@ -0,0 +1,30 @@
+!!----------------------------------------------------------------------
+!!                    ***  domzgr_substitute.h90   ***
+!!----------------------------------------------------------------------
+!! ** purpose :   substitute fsdep. and fse.., the vert. depth and scale
+!!      factors depending on the vertical coord. used, using CPP macro.
+!!----------------------------------------------------------------------
+!! History :  4.2  !  2020-02  (S. Techene, G. Madec)  star coordinate
+!!----------------------------------------------------------------------
+!! NEMO/OCE 4.2 , NEMO Consortium (2020)
+!! $Id$
+!! Software governed by the CeCILL license (see ./LICENSE)
+!!----------------------------------------------------------------------
+#if defined key_qco
+#   define  e3t(i,j,k,t)   (e3t_0(i,j,k)*(1._wp+r3t(i,j,t)))
+#   define  e3u(i,j,k,t)   (e3u_0(i,j,k)*(1._wp+r3u(i,j,t)))
+#   define  e3v(i,j,k,t)   (e3v_0(i,j,k)*(1._wp+r3v(i,j,t)))
+#   define  e3f(i,j,k)     (e3f_0(i,j,k)*(1._wp+r3f(i,j)))
+#   define  e3w(i,j,k,t)   (e3w_0(i,j,k)*(1._wp+r3t(i,j,t)))
+#   define  e3uw(i,j,k,t)  (e3uw_0(i,j,k)*(1._wp+r3u(i,j,t)))
+#   define  e3vw(i,j,k,t)  (e3vw_0(i,j,k)*(1._wp+r3v(i,j,t)))
+#   define  ht(i,j)        (ht_0(i,j)+ssh(i,j,Kmm))
+#   define  hu(i,j,t)      (hu_0(i,j)*(1._wp+r3u(i,j,t)))
+#   define  hv(i,j,t)      (hv_0(i,j)*(1._wp+r3v(i,j,t)))
+#   define  r1_hu(i,j,t)   (r1_hu_0(i,j)/(1._wp+r3u(i,j,t)))
+#   define  r1_hv(i,j,t)   (r1_hv_0(i,j)/(1._wp+r3v(i,j,t)))
+#   define  gdept(i,j,k,t) (gdept_0(i,j,k)*(1._wp+r3t(i,j,t)))
+#   define  gdepw(i,j,k,t) (gdepw_0(i,j,k)*(1._wp+r3t(i,j,t)))
+#   define  gde3w(i,j,k)   (gdept_0(i,j,k)*(1._wp+r3t(i,j,Kmm))-ssh(i,j,Kmm))
+#endif
+!!----------------------------------------------------------------------
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynadv.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynadv.F90	(revision 13059)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynadv.F90	(revision 13151)
@@ -74,5 +74,4 @@
       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/dynatf.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynatf.F90	(revision 13059)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynatf.F90	(revision 13151)
@@ -58,4 +58,21 @@
 
    PUBLIC    dyn_atf   ! routine called by step.F90
+!!st22
+#if defined key_qco
+   !!----------------------------------------------------------------------
+   !!   'key_qco'      EMPTY ROUTINE     Quasi-Eulerian vertical coordonate
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dyn_atf ( kt, Kbb, Kmm, Kaa, puu, pvv, pe3t, pe3u, pe3v )
+      INTEGER                             , INTENT(in   ) :: kt               ! ocean time-step index
+      INTEGER                             , INTENT(in   ) :: Kbb, Kmm, Kaa    ! before and after time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv         ! velocities to be time filtered
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: pe3t, pe3u, pe3v ! scale factors to be time filtered
+
+      WRITE(*,*) 'dyn_atf: You should not have seen this print! error?', kt
+   END SUBROUTINE dyn_atf
+
+#else
 
    !! * Substitutions
@@ -312,4 +329,6 @@
    END SUBROUTINE dyn_atf
 
+#endif
+!!st22
    !!=========================================================================
 END MODULE dynatf
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynkeg.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynkeg.F90	(revision 13059)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynkeg.F90	(revision 13151)
@@ -29,5 +29,4 @@
 
    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)
@@ -156,116 +155,5 @@
       !
    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 13059)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynldf_lap_blp.F90	(revision 13151)
@@ -20,8 +20,5 @@
    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
@@ -34,8 +31,10 @@
    INTEGER, PUBLIC, PARAMETER ::   np_dynldf_lap_symN = 3         ! symmetric laplacian (cartesian)
    
-   !INTEGER, PUBLIC, PARAMETER ::   nn_dynldf_lap_typ = 1         ! choose type of laplacian (ideally from namelist)
+   INTEGER, PUBLIC, PARAMETER ::   ln_dynldf_lap_typ = 1         ! choose type of laplacian (ideally from namelist)
 !!anSYM
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+!!st21
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -80,9 +79,9 @@
          WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass
          WRITE(numout,*) '~~~~~~~ '
-         WRITE(numout,*) '                  nn_dynldf_lap_typ = ', nn_dynldf_lap_typ
-         SELECT CASE( nn_dynldf_lap_typ )             ! print the choice of operator
+         WRITE(numout,*) '                  ln_dynldf_lap_typ = ', ln_dynldf_lap_typ
+         SELECT CASE( ln_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 (cartesian form)'
+         CASE( np_dynldf_lap_symN)   ;   WRITE(numout,*) '   ==>>>   symmetric laplacian (simple form)'
          END SELECT
       ENDIF
@@ -92,5 +91,5 @@
       ENDIF
       !
-      SELECT CASE( nn_dynldf_lap_typ )  
+      SELECT CASE( ln_dynldf_lap_typ )  
          !              
          CASE ( np_dynldf_lap_rot )       !==  Vorticity-Divergence form  ==!
@@ -102,8 +101,9 @@
 !!gm open question here : e3f  at before or now ?    probably now... 
 !!gm note that ahmf has already been multiplied by fmask
-                  zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &
-                     &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  &
-                     &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  )
-                  !                                      ! ahm * div        (computed from 2 to jpi/jpj)
+            zcur(ji-1,jj-1) =  &
+               &      ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)      &
+               &  * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  &
+               &     - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  )
+            !                                      ! ahm * div        (computed from 2 to jpi/jpj)
 !!gm note that ahmt has already been multiplied by tmask
                   zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)                                         &
@@ -160,5 +160,5 @@
             END DO                                           !   End of slab
             !
-         CASE ( np_dynldf_lap_symN )       !==  Symmetric form  ==!   (cartesian way)
+         CASE ( np_dynldf_lap_symN )       !==  Symmetric form  ==!   (naive way)
             !
             DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -193,5 +193,5 @@
             !  
          CASE DEFAULT                                     ! error
-            CALL ctl_stop('STOP','dyn_ldf_lap: wrong value for nn_dynldf_lap_typ'  )
+            CALL ctl_stop('STOP','dyn_ldf_lap: wrong value for ln_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 13059)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90	(revision 13151)
@@ -22,6 +22,4 @@
    !!             -   ! 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
    !!----------------------------------------------------------------------
 
@@ -47,7 +45,4 @@
    USE lib_mpp        ! MPP library
    USE timing         ! Timing
-!!anAD only
-   USE dynkeg, ONLY : dyn_kegAD
-   USE dynadv, ONLY : nn_dynkeg
 
    IMPLICIT NONE
@@ -58,18 +53,12 @@
 
    !                                   !!* Namelist namdyn_vor: vorticity term
-   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)
+   LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme          (ENS)
+   LOGICAL, PUBLIC ::   ln_dynvor_ene   !: f-point energy conserving scheme     (ENE)
+   LOGICAL, PUBLIC ::   ln_dynvor_enT   !: t-point energy conserving scheme     (ENT)
+   LOGICAL, PUBLIC ::   ln_dynvor_eeT   !: t-point energy conserving scheme     (EET)
+   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy & enstrophy conserving scheme (EEN)
+   INTEGER, PUBLIC ::      nn_een_e3f      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
+   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                         (MIX)
+   LOGICAL, PUBLIC ::   ln_dynvor_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes)
 
    INTEGER, PUBLIC ::   nvor_scheme     !: choice of the type of advection scheme
@@ -81,15 +70,6 @@
    INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme
    INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme
-!!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 
+
+   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity 
    !                               ! associated indices:
    INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary)
@@ -110,4 +90,6 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+!!st23
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -117,5 +99,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_vor( kt, Kbb, Kmm, puu, pvv, Krhs)
+   SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs )
       !!----------------------------------------------------------------------
       !!
@@ -127,11 +109,9 @@
       !!               for futher diagnostics (l_trddyn=T)
       !!----------------------------------------------------------------------
-      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
+      INTEGER                             , INTENT( in  ) ::   kt          ! ocean time-step index
+      INTEGER                             , INTENT( in  ) ::   Kmm, Krhs   ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv    ! ocean velocity field and RHS of momentum equation
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdu, ztrdv
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zuu, zvv
       !!----------------------------------------------------------------------
       !
@@ -187,158 +167,8 @@
             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
+                             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
@@ -427,21 +257,25 @@
             DO_2D_01_01
                zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   &
-                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
+                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) &
+                  &                * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
             END_2D
          CASE ( np_MET )                           !* metric term
             DO_2D_01_01
-               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   &
-                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t(ji,jj,jk,Kmm)
+               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)     &
+                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) &
+                  &         * e3t(ji,jj,jk,Kmm)
             END_2D
          CASE ( np_CRV )                           !* Coriolis + relative vorticity
             DO_2D_01_01
-               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    &
-                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
+               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)      &
+                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) &
+                  &         * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
             END_2D
          CASE ( np_CME )                           !* Coriolis + metric
             DO_2D_01_01
-               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           &
-                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  &
-                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) * e3t(ji,jj,jk,Kmm)
+               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                             &
+                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)    &
+                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) &
+                    &       * e3t(ji,jj,jk,Kmm)
             END_2D
          CASE DEFAULT                                             ! error
@@ -485,9 +319,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),OPTIONAL, 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), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
       !
       INTEGER  ::   ji, jj, jk           ! dummy loop indices
@@ -543,51 +377,18 @@
          END SELECT
          !
-         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
+         !                                   !==  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
          !                                             ! ===============
       END DO                                           !   End of slab
@@ -616,9 +417,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),OPTIONAL, 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), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -674,51 +475,18 @@
          !
          !
-!!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
+         !                                   !==  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
          !                                             ! ===============
       END DO                                           !   End of slab
@@ -772,6 +540,8 @@
          CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
             DO_2D_10_10
-               ze3f = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
-                  &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
+               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
+                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
+                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
+                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
                IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f
                ELSE                       ;   z1_e3f(ji,jj) = 0._wp
@@ -780,6 +550,8 @@
          CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
             DO_2D_10_10
-               ze3f = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
-                  &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
+               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
+                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
+                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
+                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
                zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
                   &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  )
@@ -1000,7 +772,5 @@
       !!
       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_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
+         &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk
       !!----------------------------------------------------------------------
       !
@@ -1019,12 +789,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
 
@@ -1039,5 +809,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  ,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
+               & + tmask(ji,jj  ,jk) + tmask(ji+1,jj+1,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
          END_3D
          !
@@ -1049,11 +819,5 @@
       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
@@ -1072,5 +836,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  )
@@ -1102,19 +866,10 @@
          WRITE(numout,*)
          SELECT CASE( nvor_scheme )
-         
-         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)'
+         CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)'
+         CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)'
+         CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)'
+         CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)'
+         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)'
+         CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)'
          END SELECT         
       ENDIF
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/ldfdyn.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/ldfdyn.F90	(revision 13059)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/ldfdyn.F90	(revision 13151)
@@ -25,7 +25,5 @@
    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
@@ -62,7 +60,5 @@
    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)
@@ -327,38 +323,33 @@
          IF( .NOT.l_ldfdyn_time ) THEN             !* No time variation 
             IF(     ln_dynldf_lap ) THEN                 !   laplacian operator (mask only)
-!!an          !
-            WRITE(numout,*) '   ln_dynldf_lap_PM = ',ln_dynldf_lap_PM 
-               IF(     ln_dynldf_lap_PM ) THEN                 !   laplacian operator (mask only)
+               ahmt(:,:,1:jpkm1) =       ahmt(:,:,1:jpkm1)   * tmask(:,:,1:jpkm1)
+               WRITE(numout,*) ' ahmt tmask '
 !! 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,*) '  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         !                 
+!              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'
+
             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/nemogcm.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/nemogcm.F90	(revision 13059)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/nemogcm.F90	(revision 13151)
@@ -60,5 +60,11 @@
    USE diacfl         ! CFL diagnostics               (dia_cfl_init routine)
    USE diamlr         ! IOM context management for multiple-linear-regression analysis
+#if defined key_RK3
+   USE stpRK3
+#elif defined key_qco
+   USE stpLF
+#else
    USE step           ! NEMO time-stepping                 (stp     routine)
+#endif
    USE isfstp         ! ice shelf                     (isf_stp_init routine)
    USE icbini         ! handle bergs, initialisation
@@ -175,6 +181,11 @@
             IF ( istp ==         nitend ) elapsed_time = zstptiming - elapsed_time
          ENDIF
-         
-         CALL stp        ( istp ) 
+#if defined key_RK3
+         CALL stp_RK3    ( istp )
+#elif defined key_qco
+         CALL stp_LF     ( istp )
+#else
+         CALL stp        ( istp )
+#endif
          istp = istp + 1
 
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/sbcice_cice.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/sbcice_cice.F90	(revision 13059)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/sbcice_cice.F90	(revision 13151)
@@ -12,5 +12,11 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
+!!st8
+# if ! defined key_qco
    USE domvvl
+# else
+   USE domqco
+# endif
+!!st8
    USE phycst, only : rcp, rho0, r1_rho0, rhos, rhoi
    USE in_out_manager  ! I/O manager
@@ -233,10 +239,14 @@
 !!gm This should be put elsewhere....   (same remark for limsbc)
 !!gm especially here it is assumed zstar coordinate, but it can be ztilde....
+!!st9
+#if defined key_qco
+            IF( .NOT.ln_linssh )   CALL dom_qco_zgr( Kbb, Kmm, Kaa )   ! interpolation scale factor, depth and water column
+#else
             IF( .NOT.ln_linssh ) THEN
                !
                DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
-                  e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kmm)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
-                  e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kbb)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
-               ENDDO
+                  e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kmm)*r1_ht_0(:,:)*tmask(:,:,jk) )
+                  e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kbb)*r1_ht_0(:,:)*tmask(:,:,jk) )
+               END DO
                e3t(:,:,:,Krhs) = e3t(:,:,:,Kbb)
                ! Reconstruction of all vertical scale factors at now and before time-steps
@@ -267,4 +277,6 @@
                END DO
             ENDIF
+#endif
+!!st9
          ENDIF
       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 13059)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90	(revision 13151)
@@ -6,5 +6,9 @@
    !! History :  NEMO !  2020-03  (A. Nasser, G. Madec)  Original code from  4.0.2
    !!----------------------------------------------------------------------
-
+#if defined key_qco
+   !!----------------------------------------------------------------------
+   !!   'key_qco'      EMPTY MODULE      Quasi-Eulerian vertical coordonate
+   !!----------------------------------------------------------------------
+#else
    !!----------------------------------------------------------------------
    !!   stp             : Shallow Water time-stepping
@@ -13,7 +17,4 @@
    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 
@@ -122,5 +123,5 @@
       !  LATERAL  PHYSICS
       !                                                                        ! eddy diffusivity coeff.
-      IF( l_ldfdyn_time                    )   CALL ldf_dyn( kstp, Nbb )       ! eddy viscosity coeff. 
+      IF( l_ldfdyn_time )   CALL ldf_dyn( kstp, Nbb )                          ! eddy viscosity coeff. 
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -129,16 +130,20 @@
 
                             CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor)
-
-      IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn,      Naa )    ! after vertical scale factors 
-                                                       
                          uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
                          vv(:,:,:,Nrhs) = 0._wp
 
-      IF( ln_bdy     )   CALL bdy_dyn3d_dmp ( kstp, Nbb,      uu, vv, Nrhs )  ! bdy damping trends
+      IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn,      Naa )    ! after vertical scale factors 
+
+      IF( ln_bdy     )      CALL bdy_dyn3d_dmp ( kstp, Nbb,      uu, vv, Nrhs )  ! bdy damping trends
 
 #if defined key_agrif
       IF(.NOT. Agrif_Root())  & 
-               &         CALL Agrif_Sponge_dyn        ! momentum sponge
-#endif
+               &            CALL 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)
@@ -148,25 +153,4 @@
       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
@@ -175,5 +159,5 @@
             &                                  - 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)  
+            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
       END_3D   
 !!an         
@@ -182,5 +166,6 @@
       ! Leap-Frog time splitting + Robert-Asselin time filter on u,v,e3 
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      
+          
+!! what about  IF( .NOT.ln_linssh )  ?
 !!an futur module dyn_nxt (a la place de dyn_atf)
       
@@ -209,5 +194,5 @@
                uu(ji,jj,jk,Naa) = zua
                vv(ji,jj,jk,Naa) = zva
-             END_3D   
+            END_3D
          ENDIF
          !
@@ -220,9 +205,9 @@
                zue3a = zue3b + rDt * e3u(ji,jj,jk,Nrhs) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
                zve3a = zve3b + rDt * e3v(ji,jj,jk,Nrhs) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
-               !                                                
+               !
                uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)    
-               vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)     
-             END_3D          
-         ELSE                             ! Leap Frog time stepping + Asselin filter         
+               vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)
+            END_3D
+         ELSE                             ! Leap Frog time stepping + Asselin filter
             DO_3D_11_11(1,jpkm1)
                zue3n = e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nnn)
@@ -239,5 +224,5 @@
                !                                                ! Asselin time filter on u,v (Nnn)
                uu(ji,jj,jk,Nnn) = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_tf
-               vv(ji,jj,jk,Nnn) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_tf           
+               vv(ji,jj,jk,Nnn) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_tf
                !
                e3u(ji,jj,jk,Nnn) = ze3u_tf
@@ -246,9 +231,10 @@
                !
                uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)    
-               vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)     
-             END_3D   
+               vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)
+            END_3D       
          ENDIF
       ENDIF
-      
+
+
       CALL lbc_lnk_multi( 'stp', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.,   &   !* local domain boundaries
          &                       uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.    )     
@@ -263,11 +249,9 @@
 !!                         CALL dyn_atf       ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v  )  ! time filtering of "now" velocities and scale factors
 !!an TO BE ADDED : a simplifier
-!                         CALL ssh_atf       ( kstp, Nbb, Nnn, Naa, ssh )                     ! time filtering of "now" sea surface height
- 
+!!                         CALL ssh_atf       ( kstp, Nbb, Nnn, Naa, ssh )                     ! time filtering of "now" sea surface height
       IF ( .NOT.( l_1st_euler ) ) THEN   ! Only do time filtering for leapfrog timesteps
          !                                                  ! filtering "now" field
          ssh(:,:,Nnn) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2 * ssh(:,:,Nnn) + ssh(:,:,Naa) )
       ENDIF
-      
 !!an 
 
@@ -280,5 +264,4 @@
       !
                          CALL dom_vvl_sf_update( kstp, Nbb, Nnn, Naa )  ! recompute vertical scale factors
-
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! diagnostics and outputs
@@ -287,6 +270,5 @@
       IF( ln_diacfl  )   CALL dia_cfl   ( kstp,      Nnn )      ! Courant number diagnostics
     
-                        CALL dia_wri   ( kstp,      Nnn )      ! ocean model: outputs
-
+                         CALL dia_wri   ( kstp,      Nnn )      ! ocean model: outputs
       !
       IF( lrst_oce   )   CALL rst_write    ( kstp, Nbb, Nnn )   ! write output ocean restart file
@@ -335,4 +317,5 @@
       !
    END SUBROUTINE stp
+#endif
    !
    !!======================================================================
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/stepLF.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/stepLF.F90	(revision 13151)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/stepLF.F90	(revision 13151)
@@ -0,0 +1,335 @@
+MODULE stpLF
+   !!======================================================================
+   !!                       ***  MODULE step  ***
+   !! Time-stepping   : manager of the shallow water equation time stepping
+   !!======================================================================
+   !! History :  NEMO !  2020-03  (A. Nasser, G. Madec)  Original code from  4.0.2
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   stp             : Shallow Water time-stepping
+   !!----------------------------------------------------------------------
+   USE step_oce         ! time stepping definition modules
+   USE phycst           ! physical constants
+   USE usrdef_nam
+   !
+   USE iom              ! xIOs server 
+   USE domqco
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   stp_LF   ! called by nemogcm.F90
+   
+   !!----------------------------------------------------------------------
+   !! time level indices
+   !!----------------------------------------------------------------------
+   INTEGER, PUBLIC ::   Nbb, Nnn, Naa, Nrhs   !! used by nemo_init
+      
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: step.F90 12614 2020-03-26 14:59:52Z gm $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+#if defined key_agrif
+   RECURSIVE SUBROUTINE stp_LF( )
+      INTEGER             ::   kstp   ! ocean time-step index
+#else
+   SUBROUTINE stp_LF( kstp )
+      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
+#endif
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE stp  ***
+      !!
+      !! ** Purpose : - Time stepping of shallow water (SHW) (momentum and ssh eqs.)
+      !!
+      !! ** Method  : -1- Update forcings
+      !!              -2- Update the ssh at Naa
+      !!              -3- Compute the momentum trends (Nrhs)
+      !!              -4- Update the horizontal velocity
+      !!              -5- Apply Asselin time filter to uu,vv,ssh
+      !!              -6- Outputs and diagnostics
+      !!----------------------------------------------------------------------
+      INTEGER ::   ji, jj, jk   ! dummy loop indice
+      INTEGER ::   indic        ! error indicator if < 0
+!!gm kcall can be removed, I guess
+      INTEGER ::   kcall        ! optional integer argument (dom_vvl_sf_nxt)
+      REAL(wp)::   z1_2rho0     ! local scalars
+      
+      REAL(wp) ::   zue3a, zue3n, zue3b    ! local scalars
+      REAL(wp) ::   zve3a, zve3n, zve3b    !   -      -
+      REAL(wp) ::   ze3t_tf, ze3u_tf, ze3v_tf, zua, zva
+      !! ---------------------------------------------------------------------
+#if defined key_agrif
+      kstp = nit000 + Agrif_Nb_Step()
+      Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
+      IF( lk_agrif_debug ) THEN
+         IF( Agrif_Root() .and. lwp)   WRITE(*,*) '---'
+         IF(lwp)   WRITE(*,*) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint()
+      ENDIF
+      IF( kstp == nit000 + 1 )   lk_agrif_fstep = .FALSE.
+# if defined key_iomput
+      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context )
+# endif
+#endif
+      !
+      IF( ln_timing )   CALL timing_start('stp_LF')
+      !
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! model timestep
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      IF( l_1st_euler ) THEN     ! start or restart with Euler 1st time-step
+         rDt   =  rn_Dt   
+         r1_Dt = 1._wp / rDt
+      ENDIF
+
+      IF ( kstp == nit000 )   ww(:,:,:) = 0._wp   ! initialize vertical velocity one for all to zero
+
+      !
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! update I/O and calendar 
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                             indic = 0                ! reset to no error condition
+                             
+      IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
+                             CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including passible AGRIF zoom)
+         IF( lk_diamlr   )   CALL dia_mlr_iom_init    ! with additional setup for multiple-linear-regression analysis
+                             CALL iom_init_closedef
+         IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid
+      ENDIF
+      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
+                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell IOM we are at time step kstp
+      IF( ln_crs         )   CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell IOM we are at time step kstp
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Update external forcing (tides, open boundaries, ice shelf interaction and surface boundary condition (including sea-ice)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( ln_tide    )   CALL tide_update( kstp )                     ! update tide potential
+      IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                        ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib) 
+      IF( ln_bdy     )   CALL bdy_dta ( kstp, Nnn )                   ! update dynamic & tracer data at open boundaries
+                         CALL sbc     ( kstp, Nbb, Nnn )              ! Sea Boundary Condition
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Ocean physics update
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !  LATERAL  PHYSICS
+      !                                                                        ! eddy diffusivity coeff.
+      IF( l_ldfdyn_time )   CALL ldf_dyn( kstp, Nbb )                          ! eddy viscosity coeff. 
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      !  Ocean dynamics : hdiv, ssh, e3, u, v, w
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+                            CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor)
+                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
+                         vv(:,:,:,Nrhs) = 0._wp
+
+      IF( ln_bdy     )      CALL bdy_dyn3d_dmp ( kstp, Nbb,      uu, vv, Nrhs )  ! bdy damping trends
+
+#if defined key_agrif
+      IF(.NOT. Agrif_Root())  & 
+               &            CALL Agrif_Sponge_dyn        ! momentum sponge
+#endif
+                            CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)	==> RHS
+ 
+                            CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity           	==> RHS
+ 
+                            CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing
+
+      IF( .NOT.ln_linssh )  CALL dom_qco_r3c   ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! "after" ssh./h._0 ratio explicit
+      !IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt_st( kstp, Nbb, Nnn,      Naa )    ! after vertical scale factors 
+!!an - calcul du gradient de pression horizontal (explicit)
+      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)
+      END_3D
+      !
+      ! add wind stress forcing and layer linear friction to the RHS 
+      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)   &
+            &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
+         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   &
+            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
+      END_3D   
+!!an         
+  
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Leap-Frog time splitting + Robert-Asselin time filter on u,v,e3 
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+!!st ssh_atf : add ssh filtering up there
+      IF ( .NOT.( l_1st_euler ) ) THEN   ! Only do time filtering for leapfrog timesteps
+         !                                                  ! filtering "now" field
+         ssh(:,:,Nnn) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2._wp * ssh(:,:,Nnn) + ssh(:,:,Naa) )
+      ENDIF
+!!st ssh_atf
+          
+!! what about  IF( .NOT.ln_linssh )  ?
+!!an futur module dyn_nxt (a la place de dyn_atf)
+      
+      IF( ln_dynadv_vec ) THEN      ! vector invariant form : applied on velocity
+         IF( l_1st_euler ) THEN           ! Euler time stepping (no Asselin filter)
+            DO_3D_00_00(1,jpkm1)
+               uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
+               vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
+             END_3D          
+         ELSE                             ! Leap Frog time stepping + Asselin filter         
+            DO_3D_11_11(1,jpkm1)
+               zua = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
+               zva = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
+               !                                                                  ! Asselin time filter on u,v (Nnn)
+               uu(ji,jj,jk,Nnn) = uu(ji,jj,jk,Nnn) + rn_atfp * (uu(ji,jj,jk,Nbb) - 2._wp * uu(ji,jj,jk,Nnn) + zua)
+               vv(ji,jj,jk,Nnn) = vv(ji,jj,jk,Nnn) + rn_atfp * (vv(ji,jj,jk,Nbb) - 2._wp * vv(ji,jj,jk,Nnn) + zva)
+               !              
+               uu(ji,jj,jk,Naa) = zua
+               vv(ji,jj,jk,Naa) = zva
+            END_3D
+            CALL dom_qco_r3c( ssh(:,:,Nnn), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn) )   ! "now" ssh/h_0 ratio from filtrered ssh
+#if ! defined key_qco
+            DO jk = 1, jpkm1
+               e3t(:,:,jk,Nnn) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Nnn) )
+               e3u(:,:,jk,Nnn) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Nnn) )
+               e3v(:,:,jk,Nnn) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Nnn) )
+            END DO
+#endif
+         ENDIF
+         !
+      ELSE                          ! flux form : applied on thickness weighted velocity
+         IF( l_1st_euler ) THEN           ! Euler time stepping (no Asselin filter)
+            DO_3D_00_00(1,jpkm1)
+               zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
+               zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
+               !                                                ! LF time stepping
+               zue3a = zue3b + rDt * e3u(ji,jj,jk,Nrhs) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
+               zve3a = zve3b + rDt * e3v(ji,jj,jk,Nrhs) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
+               !
+               uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)    
+               vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)
+            END_3D
+         ELSE                             ! Leap Frog time stepping + Asselin filter
+            CALL dom_qco_r3c( ssh(:,:,Nnn), r3t_f(:,:), r3u_f(:,:), r3v_f(:,:) )   ! "now" ssh/h_0 ratio from filtrered ssh
+            DO_3D_11_11(1,jpkm1)
+               zue3n = e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nnn)
+               zve3n = e3v(ji,jj,jk,Nnn) * vv(ji,jj,jk,Nnn)
+               zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
+               zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
+               !                                                ! LF time stepping
+               zue3a = zue3b + rDt * e3u(ji,jj,jk,Nrhs) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
+               zve3a = zve3b + rDt * e3v(ji,jj,jk,Nrhs) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
+               !                                                ! Asselin time filter on e3u/v/t
+               ze3u_tf = e3u(ji,jj,jk,Nnn) + rn_atfp * ( e3u(ji,jj,jk,Nbb) - 2._wp * e3u(ji,jj,jk,Nnn)  + e3u(ji,jj,jk,Naa) ) 
+               ze3v_tf = e3v(ji,jj,jk,Nnn) + rn_atfp * ( e3v(ji,jj,jk,Nbb) - 2._wp * e3v(ji,jj,jk,Nnn)  + e3v(ji,jj,jk,Naa) )
+               ze3t_tf = e3t(ji,jj,jk,Nnn) + rn_atfp * ( e3t(ji,jj,jk,Nbb) - 2._wp * e3t(ji,jj,jk,Nnn)  + e3t(ji,jj,jk,Naa) ) 
+               !                                                ! Asselin time filter on u,v (Nnn)
+               uu(ji,jj,jk,Nnn) = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / (e3u_0(ji,jj,jk) * ( 1._wp + r3u_f(ji,jj) ))
+               vv(ji,jj,jk,Nnn) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / (e3v_0(ji,jj,jk) * ( 1._wp + r3v_f(ji,jj) ))
+               !
+               uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)    
+               vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)
+            END_3D
+            r3t(:,:,Nnn) = r3t_f(:,:)
+            r3u(:,:,Nnn) = r3u_f(:,:)
+            r3v(:,:,Nnn) = r3v_f(:,:)
+#if ! defined key_qco
+            DO jk = 1, jpkm1
+               e3t(:,:,jk,Nnn) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Nnn) )
+               e3u(:,:,jk,Nnn) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Nnn) )
+               e3v(:,:,jk,Nnn) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Nnn) )
+            END DO
+#endif         
+         ENDIF
+      ENDIF
+
+
+      CALL lbc_lnk_multi( 'stp', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.,   &   !* local domain boundaries
+         &                       uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.    )     
+
+!!an         
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Set boundary conditions, time filter and swap time levels
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+!!an TO BE ADDED : dyn_nxt
+!!                         CALL dyn_atf       ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v  )  ! time filtering of "now" velocities and scale factors
+!!an TO BE ADDED : a simplifier
+!                         CALL ssh_atf       ( kstp, Nbb, Nnn, Naa, ssh )                     ! time filtering of "now" sea surface height
+!!st copied above 
+!!      IF ( .NOT.( l_1st_euler ) ) THEN   ! Only do time filtering for leapfrog timesteps
+!!         !                                                  ! filtering "now" field
+!!         ssh(:,:,Nnn) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2 * ssh(:,:,Nnn) + ssh(:,:,Naa) )
+!!      ENDIF
+!!st      
+!!an 
+
+
+      ! Swap time levels
+      Nrhs = Nbb
+      Nbb = Nnn
+      Nnn = Naa
+      Naa = Nrhs
+      !
+!                         CALL dom_vvl_sf_update_st( kstp, Nbb, Nnn, Naa )  ! recompute vertical scale factors
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! diagnostics and outputs
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( ln_floats  )   CALL flo_stp   ( kstp, Nbb, Nnn )      ! drifting Floats
+      IF( ln_diacfl  )   CALL dia_cfl   ( kstp,      Nnn )      ! Courant number diagnostics
+    
+                         CALL dia_wri   ( kstp,      Nnn )      ! ocean model: outputs
+
+      !
+      IF( lrst_oce   )   CALL rst_write    ( kstp, Nbb, Nnn )   ! write output ocean restart file
+          
+
+#if defined key_agrif
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! AGRIF
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<      
+                         Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs      ! agrif_oce module copies of time level indices
+                         CALL Agrif_Integrate_ChildGrids( stp_LF )       ! allows to finish all the Child Grids before updating
+
+                         IF( Agrif_NbStepint() == 0 ) THEN
+                            CALL Agrif_update_all( )                  ! Update all components
+                         ENDIF
+#endif
+      IF( ln_diaobs  )   CALL dia_obs      ( kstp, Nnn )      ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Control
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         CALL stp_ctl      ( kstp, Nbb, Nnn, indic )
+               
+                         
+      IF( kstp == nit000 ) THEN                          ! 1st time step only
+                                        CALL iom_close( numror )   ! close input  ocean restart file
+         IF(lwm)                        CALL FLUSH    ( numond )   ! flush output namelist oce
+         IF(lwm .AND. numoni /= -1 )    CALL FLUSH    ( numoni )   ! flush output namelist ice (if exist)
+      ENDIF
+
+      !
+#if defined key_iomput
+      IF( kstp == nitend .OR. indic < 0 ) THEN 
+                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
+                      IF(lrxios) CALL iom_context_finalize(      crxios_context          )
+      ENDIF
+#endif
+      !
+      IF( l_1st_euler ) THEN         ! recover Leap-frog timestep
+         rDt = 2._wp * rn_Dt   
+         r1_Dt = 1._wp / rDt
+         l_1st_euler = .FALSE.      
+      ENDIF
+      !
+      IF( ln_timing )   CALL timing_stop('stp')
+      !
+   END SUBROUTINE stp_LF
+   !
+   !!======================================================================
+END MODULE stpLF
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/stpRK3.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/stpRK3.F90	(revision 13151)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/stpRK3.F90	(revision 13151)
@@ -0,0 +1,378 @@
+MODULE stpRK3
+   !!======================================================================
+   !!                       ***  MODULE step  ***
+   !! Time-stepping   : manager of the shallow water equation time stepping
+   !!======================================================================
+   !! History :  NEMO !  2020-03  (A. Nasser, G. Madec)  Original code from  4.0.2
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   stpRK3             : Shallow Water time-stepping
+   !!----------------------------------------------------------------------
+   USE step_oce         ! time stepping definition modules
+   USE phycst           ! physical constants
+   USE usrdef_nam
+   !
+   USE iom              ! xIOs server 
+   USE domqco
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   stp_RK3   ! called by nemogcm.F90
+   
+   !!----------------------------------------------------------------------
+   !! time level indices
+   !!----------------------------------------------------------------------
+   INTEGER, PUBLIC ::   Nbb, Nnn, Naa, Nrhs   !! used by nemo_init
+      
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: step.F90 12614 2020-03-26 14:59:52Z gm $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+#if defined key_agrif
+   RECURSIVE SUBROUTINE stp_RK3( )
+      INTEGER             ::   kstp   ! ocean time-step index
+#else
+   SUBROUTINE stp_RK3( kstp )
+      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
+#endif
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE stp_RK3  ***
+      !!
+      !! ** Purpose : - Time stepping of shallow water (SHW) (momentum and ssh eqs.)
+      !!
+      !! ** Method  : -1- Update forcings
+      !!              -2- Update the ssh at Naa
+      !!              -3- Compute the momentum trends (Nrhs)
+      !!              -4- Update the horizontal velocity
+      !!              -5- Apply Asselin time filter to uu,vv,ssh
+      !!              -6- Outputs and diagnostics
+      !!----------------------------------------------------------------------
+      INTEGER ::   ji, jj, jk   ! dummy loop indice
+      INTEGER ::   indic        ! error indicator if < 0
+!!gm kcall can be removed, I guess
+      INTEGER ::   kcall        ! optional integer argument (dom_vvl_sf_nxt)
+      REAL(wp)::   z1_2rho0,  z5_6,  z3_4  ! local scalars
+      
+      REAL(wp) ::   zue3a, zue3n, zue3b    ! local scalars
+      REAL(wp) ::   zve3a, zve3n, zve3b    !   -      -
+      REAL(wp) ::   ze3t_tf, ze3u_tf, ze3v_tf, zua, zva
+      !! ---------------------------------------------------------------------
+#if defined key_agrif
+      kstp = nit000 + Agrif_Nb_Step()
+      Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
+      IF( lk_agrif_debug ) THEN
+         IF( Agrif_Root() .and. lwp)   WRITE(*,*) '---'
+         IF(lwp)   WRITE(*,*) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint()
+      ENDIF
+      IF( kstp == nit000 + 1 )   lk_agrif_fstep = .FALSE.
+# if defined key_iomput
+      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context )
+# endif
+#endif
+      !
+      IF( ln_timing )   CALL timing_start('stp_RK3')
+      !
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! model timestep
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      IF ( kstp == nit000 )   ww(:,:,:) = 0._wp   ! initialize vertical velocity one for all to zero
+
+      !
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! update I/O and calendar 
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                             indic = 0                ! reset to no error condition
+                             
+      IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
+                             CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including passible AGRIF zoom)
+         IF( lk_diamlr   )   CALL dia_mlr_iom_init    ! with additional setup for multiple-linear-regression analysis
+                             CALL iom_init_closedef
+         IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid
+      ENDIF
+      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
+                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell IOM we are at time step kstp
+      IF( ln_crs         )   CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell IOM we are at time step kstp
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Update external forcing (tides, open boundaries, ice shelf interaction and surface boundary condition (including sea-ice)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( ln_tide    )   CALL tide_update( kstp )                     ! update tide potential
+      IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                        ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib) 
+      IF( ln_bdy     )   CALL bdy_dta ( kstp, Nnn )                   ! update dynamic & tracer data at open boundaries
+                         CALL sbc     ( kstp, Nbb, Nnn )              ! Sea Boundary Condition
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Ocean physics update
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !  LATERAL  PHYSICS
+      !                                                                        ! eddy diffusivity coeff.
+      IF( l_ldfdyn_time )   CALL ldf_dyn( kstp, Nbb )                          ! eddy viscosity coeff. 
+
+
+      !======================================================================
+      !======================================================================
+      !                     =====       RK3       =====
+      !======================================================================
+      !======================================================================
+
+      
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      !  RK3 1st stage Ocean dynamics : hdiv, ssh, e3, u, v, w
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      rDt   = rn_Dt / 3._wp  
+      r1_Dt = 1._wp / rDt
+      
+                            CALL ssh_nxt       ( kstp, Nbb, Nbb, ssh, Naa )    ! after ssh (includes call to div_hor)
+
+                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
+                         vv(:,:,:,Nrhs) = 0._wp
+
+                            CALL dyn_adv( kstp, Nbb, Nbb      , uu, vv, Nrhs )  ! advection (VF or FF)	==> RHS
+ 
+                            CALL dyn_vor( kstp,      Nbb      , uu, vv, Nrhs )  ! vorticity           	==> RHS
+#if defined key_RK3all 
+                            CALL dyn_ldf( kstp, Nbb, Nbb      , uu, vv, Nrhs )  ! lateral mixing
+#endif
+      !
+!!an - calcul du gradient de pression horizontal (explicit)
+      DO_3D_00_00( 1, jpkm1 )
+         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nbb) - ssh(ji,jj,Nbb) ) * r1_e1u(ji,jj)
+         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nbb) - ssh(ji,jj,Nbb) ) * r1_e2v(ji,jj)
+      END_3D
+      !
+#if defined key_RK3all 
+      ! add wind stress forcing and layer linear friction to the RHS 
+      z5_6 = 5._wp/6._wp
+      DO_3D_00_00(1,jpkm1)
+         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + r1_rho0 * ( z5_6*utau_b(ji,jj) + (1._wp - z5_6)*utau(ji,jj) ) / e3u(ji,jj,jk,Nbb)   &
+            &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
+         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + r1_rho0 * ( z5_6*vtau_b(ji,jj) + (1._wp - z5_6)*vtau(ji,jj) ) / e3v(ji,jj,jk,Nbb)   &
+            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
+      END_3D
+#endif
+!!an
+                            CALL dom_qco_r3c   ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! "after" ssh./h._0 ratio explicit
+      IF( ln_dynadv_vec ) THEN      ! vector invariant form : applied on velocity
+         DO_3D_00_00(1,jpkm1)
+            uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
+            vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
+         END_3D          
+      ELSE
+         DO_3D_00_00(1,jpkm1)       ! flux form : applied on thickness weighted velocity
+            uu(ji,jj,jk,Naa) = (         uu(ji,jj,jk,Nbb )*e3u(ji,jj,jk,Nbb)                              &
+               &                 + rDt * uu(ji,jj,jk,Nrhs)*e3t(ji,jj,jk,Nbb) * umask(ji,jj,jk)        )   &
+               &               /                           e3t(ji,jj,jk,Naa)
+            vv(ji,jj,jk,Naa) = (         vv(ji,jj,jk,Nbb )*e3v(ji,jj,jk,Nbb)                              &
+               &                 + rDt * vv(ji,jj,jk,Nrhs)*e3t(ji,jj,jk,Nbb) * vmask(ji,jj,jk)        )   &
+               &               /                           e3t(ji,jj,jk,Naa)
+         END_3D
+      ENDIF
+      ! Swap time levels
+      Nrhs= Nnn
+      Nnn = Naa
+      Naa = Nrhs
+
+      
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      !  RK3 2nd stage Ocean dynamics : hdiv, ssh, e3, u, v, w
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      rDt   = rn_Dt / 2._wp  
+      r1_Dt = 1._wp / rDt
+      
+                            CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor)
+
+                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
+                         vv(:,:,:,Nrhs) = 0._wp
+!!st TBC for dyn_adv
+                            CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)	==> RHS 
+ 
+                            CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity           	==> RHS
+#if defined key_RK3all  
+                            CALL dyn_ldf( kstp, Nbb, Nbb      , uu, vv, Nrhs )  ! lateral mixing
+#endif
+                            
+      !
+!!an - calcul du gradient de pression horizontal (explicit)
+      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)
+      END_3D
+      !
+      ! add wind stress forcing and layer linear friction to the RHS 
+#if defined key_RK3all
+      z3_4 = 3._wp/4._wp
+      DO_3D_00_00(1,jpkm1)
+         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + r1_rho0 * ( z3_4*utau_b(ji,jj) + (1._wp - z3_4)*utau(ji,jj) ) / e3u(ji,jj,jk,Nbb)   &
+            &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
+         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + r1_rho0 * ( z3_4*vtau_b(ji,jj) + (1._wp - z3_4)*vtau(ji,jj) ) / e3v(ji,jj,jk,Nbb)   &
+            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
+      END_3D
+#endif
+!!an
+                           CALL dom_qco_r3c   ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! "after" ssh./h._0 ratio explicit
+      IF( ln_dynadv_vec ) THEN      ! vector invariant form : applied on velocity
+         DO_3D_00_00(1,jpkm1)
+            uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
+            vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
+         END_3D          
+      ELSE
+         DO_3D_00_00(1,jpkm1)       ! flux form : applied on thickness weighted velocity
+            uu(ji,jj,jk,Naa) = (         uu(ji,jj,jk,Nbb )*e3u(ji,jj,jk,Nbb)                              &
+               &                 + rDt * uu(ji,jj,jk,Nrhs)*e3t(ji,jj,jk,Nnn) * umask(ji,jj,jk)        )   &
+               &               /                           e3t(ji,jj,jk,Naa)
+            vv(ji,jj,jk,Naa) = (         vv(ji,jj,jk,Nbb )*e3v(ji,jj,jk,Nbb)                              &
+               &                 + rDt * vv(ji,jj,jk,Nrhs)*e3t(ji,jj,jk,Nnn) * vmask(ji,jj,jk)        )   &
+               &               /                           e3t(ji,jj,jk,Naa)
+         END_3D
+      ENDIF
+      ! Swap time levels
+      Nrhs= Nnn
+      Nnn = Naa
+      Naa = Nrhs
+       
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      !  RK3 3rd stage Ocean dynamics : hdiv, ssh, e3, u, v, w
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      rDt   = rn_Dt
+      r1_Dt = 1._wp / rDt
+
+                            CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor)
+      
+                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
+                         vv(:,:,:,Nrhs) = 0._wp
+
+      IF( ln_bdy     )      CALL bdy_dyn3d_dmp ( kstp, Nbb,      uu, vv, Nrhs )  ! bdy damping trends
+
+#if defined key_agrif
+      IF(.NOT. Agrif_Root())  & 
+               &            CALL Agrif_Sponge_dyn        ! momentum sponge
+#endif
+                            CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)	==> RHS
+ 
+                            CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity           	==> RHS
+ 
+                            CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing
+
+!!an - calcul du gradient de pression horizontal (explicit)
+      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)
+      END_3D
+      !
+      ! add wind stress forcing and layer linear friction to the RHS 
+      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)   &
+            &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
+         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   &
+            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
+      END_3D   
+!!an         
+                            CALL dom_qco_r3c   ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! "after" ssh./h._0 ratio explicit      
+      IF( ln_dynadv_vec ) THEN      ! vector invariant form : applied on velocity
+         DO_3D_11_11(1,jpkm1)
+            zua = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
+            zva = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
+            !                                                                  ! Asselin time filter on u,v (Nnn)
+            uu(ji,jj,jk,Nnn) = uu(ji,jj,jk,Nnn) + rn_atfp * (uu(ji,jj,jk,Nbb) - 2._wp * uu(ji,jj,jk,Nnn) + zua)
+            vv(ji,jj,jk,Nnn) = vv(ji,jj,jk,Nnn) + rn_atfp * (vv(ji,jj,jk,Nbb) - 2._wp * vv(ji,jj,jk,Nnn) + zva)
+            !              
+            uu(ji,jj,jk,Naa) = zua
+            vv(ji,jj,jk,Naa) = zva
+         END_3D
+         !
+      ELSE                          ! flux form : applied on thickness weighted velocity
+         DO_3D_11_11(1,jpkm1)
+            zue3n = e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nnn)
+            zve3n = e3v(ji,jj,jk,Nnn) * vv(ji,jj,jk,Nnn)
+            zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
+            zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
+            !                                                ! LF time stepping
+            zue3a = zue3b + rDt * e3t(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
+            zve3a = zve3b + rDt * e3t(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
+            !
+            uu(ji,jj,jk,Naa) = zue3a / e3t(ji,jj,jk,Naa)    
+            vv(ji,jj,jk,Naa) = zve3a / e3t(ji,jj,jk,Naa)
+         END_3D
+!!st je ne comprends pas l'histoire des e3t et du Nbb et pas du Nnn pour le rhs ?        
+      ENDIF
+
+
+      CALL lbc_lnk_multi( 'stp_RK3', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.,   &   !* local domain boundaries
+         &                       uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.    )     
+
+!!an         
+
+
+      ! Swap time levels
+      Nrhs = Nbb
+      Nbb = Naa
+      Naa = Nrhs
+      !
+!                         CALL dom_vvl_sf_update_st( kstp, Nbb, Nnn, Naa )  ! recompute vertical scale factors
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! diagnostics and outputs
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( ln_floats  )   CALL flo_stp   ( kstp, Nbb, Nnn )      ! drifting Floats
+      IF( ln_diacfl  )   CALL dia_cfl   ( kstp,      Nnn )      ! Courant number diagnostics
+    
+                         CALL dia_wri   ( kstp,      Nnn )      ! ocean model: outputs
+
+      !
+      IF( lrst_oce   )   CALL rst_write    ( kstp, Nbb, Nnn )   ! write output ocean restart file
+          
+
+#if defined key_agrif
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! AGRIF
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<      
+                         Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs      ! agrif_oce module copies of time level indices
+                         CALL Agrif_Integrate_ChildGrids( stp_RK3 )       ! allows to finish all the Child Grids before updating
+
+                         IF( Agrif_NbStepint() == 0 ) THEN
+                            CALL Agrif_update_all( )                  ! Update all components
+                         ENDIF
+#endif
+      IF( ln_diaobs  )   CALL dia_obs      ( kstp, Nnn )      ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Control
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         CALL stp_ctl      ( kstp, Nbb, Nnn, indic )
+               
+                         
+      IF( kstp == nit000 ) THEN                          ! 1st time step only
+                                        CALL iom_close( numror )   ! close input  ocean restart file
+         IF(lwm)                        CALL FLUSH    ( numond )   ! flush output namelist oce
+         IF(lwm .AND. numoni /= -1 )    CALL FLUSH    ( numoni )   ! flush output namelist ice (if exist)
+      ENDIF
+
+      !
+#if defined key_iomput
+      IF( kstp == nitend .OR. indic < 0 ) THEN 
+                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
+                      IF(lrxios) CALL iom_context_finalize(      crxios_context          )
+      ENDIF
+#endif
+      !
+      IF( l_1st_euler ) THEN         ! recover Leap-frog timestep
+         rDt = 2._wp * rn_Dt   
+         r1_Dt = 1._wp / rDt
+         l_1st_euler = .FALSE.      
+      ENDIF
+      !
+      IF( ln_timing )   CALL timing_stop('stp_RK3')
+      !
+   END SUBROUTINE stp_RK3
+   !
+   !!======================================================================
+END MODULE stpRK3
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/stpctl.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/stpctl.F90	(revision 13059)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/stpctl.F90	(revision 13151)
@@ -35,4 +35,6 @@
    INTEGER  ::   idrun, idtime, idssh, idu, ids1, ids2, idt1, idt2, idc1, idw1, istatus
    LOGICAL  ::   lsomeoce
+!!stoops
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
