Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/dom_oce.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/dom_oce.F90	(revision 13426)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/dom_oce.F90	(revision 13427)
@@ -287,5 +287,9 @@
          &       e3w_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk)                      ,  STAT=ierr(ii) )
          !
-#if ! defined key_qco
+#if defined key_qco
+      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) )             
+#else
       ii = ii+1
       ALLOCATE( e3t(jpi,jpj,jpk,jpt) , e3u (jpi,jpj,jpk,jpt) , e3v (jpi,jpj,jpk,jpt) , e3f(jpi,jpj,jpk) ,      &
@@ -294,8 +298,4 @@
          !
       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)  )
@@ -304,8 +304,4 @@
       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)  )
-#else
-      ii = ii+1
-      ALLOCATE(                    hu  (jpi,jpj,jpt),    hv  (jpi,jpj,jpt)                 ,       &
          &                      r1_hu  (jpi,jpj,jpt), r1_hv  (jpi,jpj,jpt)                 ,   STAT=ierr(ii)  )
 #endif
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/domain.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/domain.F90	(revision 13426)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/domain.F90	(revision 13427)
@@ -177,5 +177,5 @@
       !
       IF( ln_linssh ) THEN       != Fix in time : set to the reference one for all
-      !
+         !
          DO jt = 1, jpt                         ! depth of t- and w-grid-points
             gdept(:,:,:,jt) = gdept_0(:,:,:)
@@ -185,20 +185,20 @@
          !
          DO jt = 1, jpt                         ! vertical scale factors
-            e3t(:,:,:,jt) =  e3t_0(:,:,:)
-            e3u(:,:,:,jt) =  e3u_0(:,:,:)
-            e3v(:,:,:,jt) =  e3v_0(:,:,:)
-            e3w(:,:,:,jt) =  e3w_0(:,:,:)
+            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(:,:,:)
+            e3f (:,:,:)    =  e3f_0(:,:,:)
          !
          DO jt = 1, jpt                         ! water column thickness and its inverse
-            hu(:,:,jt)    =    hu_0(:,:)
-            hv(:,:,jt)    =    hv_0(:,:)
+               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(:,:)
+               ht   (:,:) =    ht_0(:,:)
          !
       ELSE                       != time varying : initialize before/now/after variables
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/domqco.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/domqco.F90	(revision 13426)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/domqco.F90	(revision 13427)
@@ -8,13 +8,14 @@
    !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates
    !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
-   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping
-   !!            4.x  !  2020-02  (G. Madec, S. Techene) pure z* (quasi-eulerian) coordinate
+   !!            4.1  !  2019-08  (A. Coward, D. Storkey) add time level indices for prognostic variables
+   !!            4.x  !  2020-02  (S. Techene, G. Madec) quasi-eulerian coordinate (z* or s*) from domvvl
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   dom_qe_init   : define initial vertical scale factors, depths and column thickness
-   !!   dom_qe_r3c    : Compute ssh/h_0 ratioat t-, u-, v-, and optionally f-points
-   !!       qe_rst_read : read/write restart file
-   !!   dom_qe_ctl    : Check the vvl options
+   !!   dom_qco_init   : define initial vertical scale factors, depths and column thickness
+   !!   dom_qco_zgr    : Set ssh/h_0 ratio at t
+   !!   dom_qco_r3c    : Compute ssh/h_0 ratio at t-, u-, v-, and optionally f-points
+   !!       qco_rst_read : read/write restart file
+   !!       qco_ctl    : Check the vvl options
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers
@@ -79,5 +80,6 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
+      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! time level indices
+      !!----------------------------------------------------------------------
       !
       IF(lwp) WRITE(numout,*)
@@ -85,33 +87,33 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
       !
-      CALL dom_qco_ctl     ! choose vertical coordinate (z_star, z_tilde or layer)
-      !
-      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
-      CALL qe_rst_read( nit000, Kbb, Kmm )
-      !
-      CALL dom_qco_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column
-      !
-      ! 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
+      CALL qco_ctl                            ! choose vertical coordinate (z_star, z_tilde or layer)
+      !              ! CAUTION COM A METTRE !!!
+!!st      CALL qco_rst_read2( nit000, Kbb, Kmm )  ! Read or initialize ssh_(Kbb/Kmm) and r3
+!!st CAUTION if read2 removed change restart.F90 ! 
+      !
+      CALL qco_rst_read( nit000, Kbb, Kmm )   ! Read or initialize ssh_(Kbb/Kmm)
+      !
+      CALL dom_qco_zgr( Kbb, Kmm, Kaa )       ! interpolation scale factor, depth and water column
+      !
+      IF(lwxios) THEN   ! define variables in restart file when writing with XIOS
+         CALL iom_set_rstw_var_active('sshb')
+         CALL iom_set_rstw_var_active('sshn')
+      ENDIF
       !
    END SUBROUTINE dom_qco_init
 
 
-   SUBROUTINE dom_qco_zgr(Kbb, Kmm, Kaa)
+   SUBROUTINE dom_qco_zgr( Kbb, Kmm, Kaa )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE dom_qco_init  ***
       !!
-      !! ** Purpose :  Initialization of all ssh. to h._0 ratio
-      !!
-      !! ** Method  :  - interpolate scale factors
+      !! ** Purpose :  Initialization of all ssh./h._0 ratio
+      !!
+      !! ** Method  :  - call domqco using Kbb and Kmm
       !!
       !! ** Action  : - r3(t/u/v)_b
       !!              - r3(t/u/v/f)_n
-      !!
-      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! time level indices
       !!----------------------------------------------------------------------
       !
@@ -121,4 +123,11 @@
       CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) )
       !
+#if defined key_agrif
+      ! We need to define r3[tuv](Kaa) for AGRIF initialisation (should not be a problem for the restartability...)
+      r3t(:,:,Kaa) = r3t(:,:,Kmm)
+      r3u(:,:,Kaa) = r3u(:,:,Kmm)
+      r3v(:,:,Kaa) = r3v(:,:,Kmm)
+#endif
+      !
    END SUBROUTINE dom_qco_zgr
 
@@ -148,5 +157,7 @@
       !                                      !==  ratio at u-,v-point  ==!
       !
-      IF( ln_dynadv_vec ) THEN                     !- Vector Form   (thickness weighted averaging)
+!!st      IF( ln_dynadv_vec ) THEN                     !- Vector Form   (thickness weighted averaging)
+#if ! defined key_qcoTest_FluxForm
+      !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
          DO_2D( 0, 0, 0, 0 )
             pr3u(ji,jj) = 0.5_wp * (  e1e2t(ji  ,jj) * pssh(ji  ,jj)  &
@@ -155,10 +166,12 @@
                &                    + e1e2t(ji,jj+1) * pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj)
          END_2D
-      ELSE                                         !- Flux Form   (simple averaging)
+!!st      ELSE                                         !- Flux Form   (simple averaging)
+#else
          DO_2D( 0, 0, 0, 0 )
-            pr3u(ji,jj) = 0.5_wp * (  pssh(ji  ,jj) + pssh(ji+1,jj)  ) * r1_hu_0(ji,jj)
-            pr3v(ji,jj) = 0.5_wp * (  pssh(ji,jj  ) + pssh(ji,jj+1)  ) * r1_hv_0(ji,jj)
+            pr3u(ji,jj) = 0.5_wp * (  pssh(ji,jj) + pssh(ji+1,jj  )  ) * r1_hu_0(ji,jj)
+            pr3v(ji,jj) = 0.5_wp * (  pssh(ji,jj) + pssh(ji  ,jj+1)  ) * r1_hv_0(ji,jj)
          END_2D
-      ENDIF
+!!st      ENDIF
+#endif         
       !
       IF( .NOT.PRESENT( pr3f ) ) THEN              !- lbc on ratio at u-, v-points only
@@ -168,5 +181,8 @@
       ELSE                                   !==  ratio at f-point  ==!
          !
-         IF( ln_dynadv_vec )   THEN                !- Vector Form   (thickness weighted averaging)
+!!st         IF( ln_dynadv_vec )   THEN                !- Vector Form   (thickness weighted averaging)
+#if ! defined key_qcoTest_FluxForm
+         !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
+
             DO_2D( 1, 0, 1, 0 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line
                pr3f(ji,jj) = 0.25_wp * (  e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )  &
@@ -175,10 +191,12 @@
                   &                     + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj)
             END_2D
-         ELSE                                      !- Flux Form   (simple averaging)
+!!st         ELSE                                      !- Flux Form   (simple averaging)
+#else
             DO_2D( 1, 0, 1, 0 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line
-               pr3f(ji,jj) = 0.25_wp * (  pssh(ji  ,jj  ) + pssh(ji+1,jj  )  &
-                  &                     + pssh(ji  ,jj+1) + pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj)
+               pr3f(ji,jj) = 0.25_wp * (  pssh(ji,jj  ) + pssh(ji+1,jj  )  &
+                  &                     + pssh(ji,jj+1) + pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj)
             END_2D
-         ENDIF
+!!st         ENDIF
+#endif
          !                                                 ! lbc on ratio at u-,v-,f-points
          CALL lbc_lnk_multi( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp, pr3f, 'F', 1._wp )
@@ -189,7 +207,7 @@
 
 
-   SUBROUTINE qe_rst_read( kt, Kbb, Kmm )
+   SUBROUTINE qco_rst_read( kt, Kbb, Kmm )
       !!---------------------------------------------------------------------
-      !!                   ***  ROUTINE qe_rst_read  ***
+      !!                   ***  ROUTINE qco_rst_read  ***
       !!
       !! ** Purpose :   Read ssh in restart file
@@ -199,6 +217,6 @@
       !!                it is set to the _0 values.
       !!----------------------------------------------------------------------
-      INTEGER         , INTENT(in) ::   kt        ! ocean time-step
-      INTEGER         , INTENT(in) ::   Kbb, Kmm  ! ocean time level indices
+      INTEGER, INTENT(in) ::   kt         ! ocean time-step
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
       !
       INTEGER ::   ji, jj, jk
@@ -206,96 +224,214 @@
       !!----------------------------------------------------------------------
       !
-         IF( ln_rstart ) THEN                   !* Read the restart file
-            CALL rst_read_open                  !  open the restart file if necessary
-            !
-            id1 = iom_varid( numror, 'sshb', ldstop = .FALSE. )
-            id2 = iom_varid( numror, 'sshn', ldstop = .FALSE. )
-            !
-            !                             ! --------- !
-            !                             ! all cases !
-            !                             ! --------- !
-            !
-            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
-               CALL iom_get( numror, jpdom_auto, 'sshb'   , ssh(:,:,Kbb), ldxios = lrxios    )
-               CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    )
-               ! needed to restart if land processor not computed
-               IF(lwp) write(numout,*) 'qe_rst_read : ssh(:,:,Kbb) and ssh(:,:,Kmm) found in restart files'
-               WHERE ( ssmask(:,:) == 0.0_wp )   !!gm/st ==> sm should not be necessary on ssh when it was required on e3
-                  ssh(:,:,Kmm) = 0._wp
-                  ssh(:,:,Kbb) = 0._wp
-               END WHERE
-               IF( l_1st_euler ) THEN
-                  ssh(:,:,Kbb) = ssh(:,:,Kmm)
-               ENDIF
-            ELSE IF( id1 > 0 ) THEN
-               IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kmm) not found in restart files'
-               IF(lwp) write(numout,*) 'sshn set equal to sshb.'
-               IF(lwp) write(numout,*) 'neuler is forced to 0'
-               CALL iom_get( numror, jpdom_auto, 'sshb', ssh(:,:,Kbb), ldxios = lrxios )
-               ssh(:,:,Kmm) = ssh(:,:,Kbb)
-               l_1st_euler = .TRUE.
-            ELSE IF( id2 > 0 ) THEN
-               IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kbb) not found in restart files'
-               IF(lwp) write(numout,*) 'sshb set equal to sshn.'
-               IF(lwp) write(numout,*) 'neuler is forced to 0'
-               CALL iom_get( numror, jpdom_auto, 'sshn', ssh(:,:,Kmm), ldxios = lrxios )
+      IF( ln_rstart ) THEN                   !* Read the restart file
+         CALL rst_read_open                  !  open the restart file if necessary
+         !
+         id1 = iom_varid( numror, 'sshb', ldstop = .FALSE. )
+         id2 = iom_varid( numror, 'sshn', ldstop = .FALSE. )
+         !
+         !                             ! --------- !
+         !                             ! all cases !
+         !                             ! --------- !
+         !
+         IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
+            CALL iom_get( numror, jpdom_auto, 'sshb'   , ssh(:,:,Kbb), ldxios = lrxios    )
+            CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    )
+            ! needed to restart if land processor not computed
+            IF(lwp) write(numout,*) 'qco_rst_read : ssh(:,:,Kbb) and ssh(:,:,Kmm) found in restart files'
+            !!WHERE ( ssmask(:,:) == 0.0_wp )   !!gm/st ==> sm should not be necessary on ssh while it was required on e3
+            !!   ssh(:,:,Kmm) = 0._wp
+            !!   ssh(:,:,Kbb) = 0._wp
+            !!END WHERE
+            IF( l_1st_euler ) THEN
                ssh(:,:,Kbb) = ssh(:,:,Kmm)
-               l_1st_euler = .TRUE.
-            ELSE
-               IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kmm) not found in restart file'
-               IF(lwp) write(numout,*) 'ssh_b and ssh_n set to zero'
-               IF(lwp) write(numout,*) 'neuler is forced to 0'
-               ssh(:,:,:) = 0._wp
-               l_1st_euler = .TRUE.
             ENDIF
-            !
-         ELSE                                   !* Initialize at "rest"
-            !
-            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
+         ELSE IF( id1 > 0 ) THEN
+            IF(lwp) write(numout,*) 'qco_rst_read WARNING : ssh(:,:,Kmm) not found in restart files'
+            IF(lwp) write(numout,*) 'sshn set equal to sshb.'
+            IF(lwp) write(numout,*) 'neuler is forced to 0'
+            CALL iom_get( numror, jpdom_auto, 'sshb', ssh(:,:,Kbb), ldxios = lrxios )
+            ssh(:,:,Kmm) = ssh(:,:,Kbb)
+            l_1st_euler = .TRUE.
+         ELSE IF( id2 > 0 ) THEN
+            IF(lwp) write(numout,*) 'qco_rst_read WARNING : ssh(:,:,Kbb) not found in restart files'
+            IF(lwp) write(numout,*) 'sshb set equal to sshn.'
+            IF(lwp) write(numout,*) 'neuler is forced to 0'
+            CALL iom_get( numror, jpdom_auto, 'sshn', ssh(:,:,Kmm), ldxios = lrxios )
+            ssh(:,:,Kbb) = ssh(:,:,Kmm)
+            l_1st_euler = .TRUE.
+         ELSE
+            IF(lwp) write(numout,*) 'qco_rst_read WARNING : ssh(:,:,Kmm) not found in restart file'
+            IF(lwp) write(numout,*) 'ssh_b and ssh_n set to zero'
+            IF(lwp) write(numout,*) 'neuler is forced to 0'
+            ssh(:,:,:) = 0._wp
+            l_1st_euler = .TRUE.
+         ENDIF
+         !
+      ELSE                                   !* Initialize at "rest"
+         !
+         IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
+            !
+            IF( cn_cfg == 'wad' ) THEN            ! Wetting and drying test case
+               CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
+               ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
+               ssh(:,:    ,Kmm) = ssh(:,:    ,Kbb)
+               uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb)
+               vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb)
+            ELSE                                  ! if not test case
+               ssh(:,:,Kmm) = -ssh_ref
+               ssh(:,:,Kbb) = -ssh_ref
                !
-               IF( cn_cfg == 'wad' ) THEN            ! Wetting and drying test case
-                  CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
-                  ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
-                  ssh(:,:    ,Kmm) = ssh(:,:    ,Kbb)
-                  uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb)
-                  vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb)
-               ELSE                                  ! if not test case
-                  ssh(:,:,Kmm) = -ssh_ref
-                  ssh(:,:,Kbb) = -ssh_ref
-                  !
-                  DO_2D( 1, 1, 1, 1 )
-                     IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
-                        ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
-                        ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
-                     ENDIF
-                  END_2D
-               ENDIF
-
-               DO ji = 1, jpi
-                  DO jj = 1, jpj
-                     IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
-                       CALL ctl_stop( 'qe_rst_read: ht_0 must be positive at potentially wet points' )
-                     ENDIF
-                  END DO
+               DO_2D( 1, 1, 1, 1 )
+                  IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
+                     ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
+                     ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
+                  ENDIF
+               END_2D
+            ENDIF
+            !
+            DO ji = 1, jpi
+               DO jj = 1, jpj
+                  IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
+                    CALL ctl_stop( 'qco_rst_read: ht_0 must be positive at potentially wet points' )
+                  ENDIF
                END DO
+            END DO
+            !
+         ELSE
+            !
+            ! Just to read set ssh in fact, called latter once vertical grid
+            ! is set up:
+!           CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
+!           !
+            ssh(:,:,:) = 0._wp
+            !
+         ENDIF           ! end of ll_wd edits
+         !
+      ENDIF
+      !
+   END SUBROUTINE qco_rst_read
+
+   
+   SUBROUTINE qco_rst_read2( kt, Kbb, Kmm )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE qco_rst_read  ***
+      !!
+      !! ** Purpose :   Read ssh in restart file
+      !!
+      !! ** Method  :   use of IOM library
+      !!                if the restart does not contain ssh,
+      !!                it is set to the _0 values.
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt         ! ocean time-step
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
+      !
+      INTEGER ::   ji, jj, jk
+      INTEGER ::   id1, id2     ! local integers
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_rstart ) THEN                   !* Read the restart file
+         CALL rst_read_open                  !  open the restart file if necessary
+         !
+         id1 = iom_varid( numror, 'sshb', ldstop = .FALSE. )
+         id2 = iom_varid( numror, 'sshn', ldstop = .FALSE. )
+         !
+         !                             ! --------- !
+         !                             ! all cases !
+         !                             ! --------- !
+         !
+         IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
+            CALL iom_get( numror, jpdom_auto, 'sshb'   , ssh(:,:,Kbb), ldxios = lrxios    )
+            CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    )
+            CALL iom_get( numror, jpdom_auto, 'r3tb'   , r3t(:,:,Kbb), ldxios = lrxios    )
+            CALL iom_get( numror, jpdom_auto, 'r3tn'   , r3t(:,:,Kmm), ldxios = lrxios    )
+            CALL iom_get( numror, jpdom_auto, 'r3ub'   , r3u(:,:,Kbb), ldxios = lrxios, cd_type = 'U'    )
+            CALL iom_get( numror, jpdom_auto, 'r3un'   , r3u(:,:,Kmm), ldxios = lrxios, cd_type = 'U'    )
+            CALL iom_get( numror, jpdom_auto, 'r3vb'   , r3v(:,:,Kbb), ldxios = lrxios, cd_type = 'V'    )
+            CALL iom_get( numror, jpdom_auto, 'r3vn'   , r3v(:,:,Kmm), ldxios = lrxios, cd_type = 'V'    )
+            CALL iom_get( numror, jpdom_auto, 'r3f'    , r3f(:,:)    , ldxios = lrxios, cd_type = 'F'    )
+            
+            ! needed to restart if land processor not computed
+            IF(lwp) write(numout,*) 'qco_rst_read : ssh(:,:,Kbb) and ssh(:,:,Kmm) found in restart files'
+            !!WHERE ( ssmask(:,:) == 0.0_wp )   !!gm/st ==> sm should not be necessary on ssh while it was required on e3
+            !!   ssh(:,:,Kmm) = 0._wp
+            !!   ssh(:,:,Kbb) = 0._wp
+            !!END WHERE
+            IF( l_1st_euler ) THEN
+               ssh(:,:,Kbb) = ssh(:,:,Kmm)
+            ENDIF
+         ELSE IF( id1 > 0 ) THEN
+            IF(lwp) write(numout,*) 'qco_rst_read WARNING : ssh(:,:,Kmm) not found in restart files'
+            IF(lwp) write(numout,*) 'sshn set equal to sshb.'
+            IF(lwp) write(numout,*) 'neuler is forced to 0'
+            CALL iom_get( numror, jpdom_auto, 'sshb', ssh(:,:,Kbb), ldxios = lrxios )
+            ssh(:,:,Kmm) = ssh(:,:,Kbb)
+            l_1st_euler = .TRUE.
+         ELSE IF( id2 > 0 ) THEN
+            IF(lwp) write(numout,*) 'qco_rst_read WARNING : ssh(:,:,Kbb) not found in restart files'
+            IF(lwp) write(numout,*) 'sshb set equal to sshn.'
+            IF(lwp) write(numout,*) 'neuler is forced to 0'
+            CALL iom_get( numror, jpdom_auto, 'sshn', ssh(:,:,Kmm), ldxios = lrxios )
+            ssh(:,:,Kbb) = ssh(:,:,Kmm)
+            l_1st_euler = .TRUE.
+         ELSE
+            IF(lwp) write(numout,*) 'qco_rst_read WARNING : ssh(:,:,Kmm) not found in restart file'
+            IF(lwp) write(numout,*) 'ssh_b and ssh_n set to zero'
+            IF(lwp) write(numout,*) 'neuler is forced to 0'
+            ssh(:,:,:) = 0._wp
+            l_1st_euler = .TRUE.
+         ENDIF
+         !
+      ELSE                                   !* Initialize at "rest"
+         !
+         IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
+            !
+            IF( cn_cfg == 'wad' ) THEN            ! Wetting and drying test case
+               CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
+               ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
+               ssh(:,:    ,Kmm) = ssh(:,:    ,Kbb)
+               uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb)
+               vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb)
+            ELSE                                  ! if not test case
+               ssh(:,:,Kmm) = -ssh_ref
+               ssh(:,:,Kbb) = -ssh_ref
                !
-            ELSE
-               !
-               ! Just to read set ssh in fact, called latter once vertical grid
-               ! is set up:
-!               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
-!               !
-                ssh(:,:,:) = 0._wp
-               !
-            ENDIF           ! end of ll_wd edits
-            !
-         ENDIF
-      !
-   END SUBROUTINE qe_rst_read
-
-
-   SUBROUTINE dom_qco_ctl
+               DO_2D( 1, 1, 1, 1 )
+                  IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
+                     ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
+                     ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
+                  ENDIF
+               END_2D
+            ENDIF
+            !
+            DO ji = 1, jpi
+               DO jj = 1, jpj
+                  IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
+                    CALL ctl_stop( 'qco_rst_read: ht_0 must be positive at potentially wet points' )
+                  ENDIF
+               END DO
+            END DO
+            !
+         ELSE
+            !
+            ! Just to read set ssh in fact, called latter once vertical grid
+            ! is set up:
+!           CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
+!           !
+            ssh(:,:,:) = 0._wp
+            r3t(:,:,:) = 0._wp
+            r3u(:,:,:) = 0._wp
+            r3v(:,:,:) = 0._wp
+            r3f(:,:  ) = 0._wp
+            !
+         ENDIF           ! end of ll_wd edits
+         !
+      ENDIF
+      !
+   END SUBROUTINE qco_rst_read2
+
+
+   SUBROUTINE qco_ctl
       !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE dom_qco_ctl  ***
+      !!                  ***  ROUTINE qco_ctl  ***
       !!
       !! ** Purpose :   Control the consistency between namelist options
@@ -317,6 +453,6 @@
       IF(lwp) THEN                    ! Namelist print
          WRITE(numout,*)
-         WRITE(numout,*) 'dom_qco_ctl : choice/control of the variable vertical coordinate'
-         WRITE(numout,*) '~~~~~~~~~~~'
+         WRITE(numout,*) 'qco_ctl : choice/control of the variable vertical coordinate'
+         WRITE(numout,*) '~~~~~~~~'
          WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
          WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
@@ -362,5 +498,5 @@
 #endif
       !
-   END SUBROUTINE dom_qco_ctl
+   END SUBROUTINE qco_ctl
 
    !!======================================================================
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/domzgr_substitute.h90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/domzgr_substitute.h90	(revision 13426)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/domzgr_substitute.h90	(revision 13427)
@@ -19,5 +19,5 @@
 #   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  ht(i,j)        (ht_0(i,j)*(1._wp+r3t(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)))
@@ -29,2 +29,5 @@
 #endif
 !!----------------------------------------------------------------------
+!!#   define  e3t_f(i,j,k)   (e3t_0(i,j,k)*(1._wp+r3t_f(i,j)*tmask(i,j,k)))
+!!#   define  e3u_f(i,j,k)   (e3u_0(i,j,k)*(1._wp+r3u_f(i,j)*umask(i,j,k)))
+!!#   define  e3v_f(i,j,k)   (e3v_0(i,j,k)*(1._wp+r3v_f(i,j)*vmask(i,j,k)))
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/istate.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/istate.F90	(revision 13426)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/istate.F90	(revision 13427)
@@ -42,5 +42,5 @@
    PRIVATE
 
-   PUBLIC   istate_init   ! routine called by step.F90
+   PUBLIC   istate_init   ! routine called by nemogcm.F90
 
    !! * Substitutions
@@ -80,8 +80,8 @@
 !!gm
 
-      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
-      rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
-      ts  (:,:,:,:,Kaa) = 0._wp                                   ! set one for all to 0 at level jpk
-      rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk
+      rhd  (:,:,:      ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
+      rn2b (:,:,:      ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
+      ts   (:,:,:,:,Kaa) = 0._wp                                   ! set one for all to 0 at level jpk
+      rab_b(:,:,:,:    ) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk
 #if defined key_agrif
       uu   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization
@@ -96,8 +96,8 @@
          CALL agrif_istate( Kbb, Kmm, Kaa )   ! Interp from parent
          !
-         ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) 
-         ssh (:,:,Kmm)     = ssh(:,:,Kbb)
-         uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
-         vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
+         ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) 
+         ssh(:,:,    Kmm) = ssh(:,:    ,Kbb)
+         uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb)
+         vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb)
       ELSE
 #endif
@@ -117,7 +117,7 @@
             CALL dta_tsd( nit000, ts(:,:,:,:,Kbb) )       ! read 3D T and S data at nit000
             !
-            ssh(:,:,Kbb)   = 0._wp               ! set the ocean at rest
-            uu  (:,:,:,Kbb) = 0._wp
-            vv  (:,:,:,Kbb) = 0._wp  
+            ssh(:,:  ,Kbb) = 0._wp               ! set the ocean at rest
+            uu (:,:,:,Kbb) = 0._wp
+            vv (:,:,:,Kbb) = 0._wp  
             !
             IF( ll_wd ) THEN
@@ -139,8 +139,8 @@
             CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )         
          ENDIF
-         ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
-         ssh (:,:,Kmm)     = ssh(:,:,Kbb)   
-         uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
-         vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
+         ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
+         ssh(:,:    ,Kmm) = ssh(:,:    ,Kbb)   
+         uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb)
+         vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb)
 
 !!gm POTENTIAL BUG :
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynadv.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynadv.F90	(revision 13426)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynadv.F90	(revision 13427)
@@ -127,5 +127,7 @@
       IF( ioptio /= 1 )   CALL ctl_stop( 'choose ONE and only ONE advection scheme' )
       IF( nn_dynkeg /= nkeg_C2 .AND. nn_dynkeg /= nkeg_HW )   CALL ctl_stop( 'KEG scheme wrong value of nn_dynkeg' )
-
+#if defined key_qcoTest_FluxForm
+      IF( ln_dynadv_vec  ) THEN CALL ctl_stop( 'STOP', 'key_qcoTest_FluxForm requires flux form advection' )
+#endif
 
       IF(lwp) THEN                    ! Print the choice
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynatf_qco.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynatf_qco.F90	(revision 13426)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynatf_qco.F90	(revision 13427)
@@ -199,19 +199,19 @@
       ! JC: Would be more clever to swap variables than to make a full vertical
       ! integration
-      !
+      ! CAUTION : calculation need to be done in the same way than see GM  
       uu_b(:,:,Kaa) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
-      uu_b(:,:,Kmm) = e3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1)
+      uu_b(:,:,Kmm) = (e3u_0(:,:,1) * ( 1._wp + r3u_f(:,:) * umask(:,:,1) )) * puu(:,:,1,Kmm) * umask(:,:,1)
       vv_b(:,:,Kaa) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
-      vv_b(:,:,Kmm) = e3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)
+      vv_b(:,:,Kmm) = (e3v_0(:,:,1) * ( 1._wp + r3v_f(:,:) * vmask(:,:,1))) * pvv(:,:,1,Kmm) * vmask(:,:,1)
       DO jk = 2, jpkm1
          uu_b(:,:,Kaa) = uu_b(:,:,Kaa) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
-         uu_b(:,:,Kmm) = uu_b(:,:,Kmm) + e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk)
+         uu_b(:,:,Kmm) = uu_b(:,:,Kmm) + (e3u_0(:,:,jk) * ( 1._wp + r3u_f(:,:) * umask(:,:,jk) )) * puu(:,:,jk,Kmm) * umask(:,:,jk)
          vv_b(:,:,Kaa) = vv_b(:,:,Kaa) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
-         vv_b(:,:,Kmm) = vv_b(:,:,Kmm) + e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)
+         vv_b(:,:,Kmm) = vv_b(:,:,Kmm) + (e3v_0(:,:,jk) * ( 1._wp + r3v_f(:,:) * vmask(:,:,jk) )) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)
       END DO
       uu_b(:,:,Kaa) = uu_b(:,:,Kaa) * r1_hu(:,:,Kaa)
       vv_b(:,:,Kaa) = vv_b(:,:,Kaa) * r1_hv(:,:,Kaa)
-      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm)
-      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm)
+      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * (r1_hu_0(:,:)/( 1._wp + r3u_f(:,:) ))
+      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * (r1_hv_0(:,:)/( 1._wp + r3v_f(:,:) ))
       !
       IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynspg_ts.F90	(revision 13426)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynspg_ts.F90	(revision 13427)
@@ -306,7 +306,7 @@
       ENDIF
       !
-      !                                   !=  Add atmospheric pressure forcing  =!
-      !                                   !  ----------------------------------  !
-      IF( ln_bt_fw ) THEN                        ! Add wind forcing
+      !                                   !=  Add wind forcing  =!
+      !                                   !  ------------------  !
+      IF( ln_bt_fw ) THEN
          DO_2D( 0, 0, 0, 0 )
             zu_frc(ji,jj) =  zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm)
@@ -386,11 +386,11 @@
       !
       IF( ln_linssh ) THEN    ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0)
-         zhup2_e(:,:) = hu(:,:,Kmm)
-         zhvp2_e(:,:) = hv(:,:,Kmm)
-         zhtp2_e(:,:) = ht(:,:)
-      ENDIF
-      !
-      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                    
-         sshn_e(:,:) =    pssh(:,:,Kmm)            
+         zhup2_e(:,:) = hu_0(:,:)
+         zhvp2_e(:,:) = hv_0(:,:)
+         zhtp2_e(:,:) = ht_0(:,:)
+      ENDIF
+      !
+      IF( ln_bt_fw ) THEN                 ! FORWARD integration: start from NOW fields                    
+         sshn_e(:,:) =    pssh (:,:,Kmm)            
          un_e  (:,:) =    puu_b(:,:,Kmm)            
          vn_e  (:,:) =    pvv_b(:,:,Kmm)
@@ -401,5 +401,5 @@
          hvr_e (:,:) = r1_hv(:,:,Kmm)
       ELSE                                ! CENTRED integration: start from BEFORE fields
-         sshn_e(:,:) =    pssh(:,:,Kbb)
+         sshn_e(:,:) =    pssh (:,:,Kbb)
          un_e  (:,:) =    puu_b(:,:,Kbb)         
          vn_e  (:,:) =    pvv_b(:,:,Kbb)
@@ -412,9 +412,9 @@
       !
       ! Initialize sums:
-      puu_b  (:,:,Kaa) = 0._wp       ! After barotropic velocities (or transport if flux form)          
-      pvv_b  (:,:,Kaa) = 0._wp
+      puu_b (:,:,Kaa) = 0._wp       ! After barotropic velocities (or transport if flux form)          
+      pvv_b (:,:,Kaa) = 0._wp
       pssh  (:,:,Kaa) = 0._wp       ! Sum for after averaged sea level
-      un_adv(:,:) = 0._wp       ! Sum for now transport issued from ts loop
-      vn_adv(:,:) = 0._wp
+      un_adv(:,:)     = 0._wp       ! Sum for now transport issued from ts loop
+      vn_adv(:,:)     = 0._wp
       !
       IF( ln_wd_dl ) THEN
@@ -475,4 +475,14 @@
             !
             !                          ! ocean u- and v-depth at mid-step   (separate DO-loops remove the need of a lbc_lnk)
+#if defined key_qcoTest_FluxForm
+            !                                ! 'key_qcoTest_FluxForm' : simple ssh average
+            DO_2D( 1, 1, 1, 0 )
+               zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * (  zsshp2_e(ji,jj) + zsshp2_e(ji+1,jj  )  ) * ssumask(ji,jj)
+            END_2D
+            DO_2D( 1, 0, 1, 1 )
+               zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * (  zsshp2_e(ji,jj) + zsshp2_e(ji  ,jj+1)  ) * ssvmask(ji,jj)
+            END_2D
+#else
+            !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
             DO_2D( 1, 1, 1, 0 )
                zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj)                        &
@@ -485,4 +495,5 @@
                     &                                 + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj)
             END_2D
+#endif               
             !
          ENDIF
@@ -540,14 +551,22 @@
          !  
          ! Sea Surface Height at u-,v-points (vvl case only)
-         IF( .NOT.ln_linssh ) THEN                                
+         IF( .NOT.ln_linssh ) THEN
+#if defined key_qcoTest_FluxForm
+            !                                ! 'key_qcoTest_FluxForm' : simple ssh average
+            DO_2D( 1, 1, 1, 0 )
+               zsshu_a(ji,jj) = r1_2 * (  ssha_e(ji,jj) + ssha_e(ji+1,jj  )  ) * ssumask(ji,jj)
+            END_2D
+            DO_2D( 1, 0, 1, 1 )
+               zsshv_a(ji,jj) = r1_2 * (  ssha_e(ji,jj) + ssha_e(ji  ,jj+1)  ) * ssvmask(ji,jj)
+            END_2D
+#else
             DO_2D( 0, 0, 0, 0 )
-               zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    &
-                  &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
-                  &              +   e1e2t(ji+1,jj  )  * ssha_e(ji+1,jj  ) )
-               zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj)    &
-                  &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
-                  &              +   e1e2t(ji  ,jj+1)  * ssha_e(ji  ,jj+1) )
-            END_2D
-         ENDIF   
+               zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji  ,jj  ) * ssha_e(ji  ,jj  )   &
+                  &                                      + e1e2t(ji+1,jj  ) * ssha_e(ji+1,jj  ) ) * ssumask(ji,jj)
+               zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji  ,jj  ) * ssha_e(ji  ,jj  )   &
+                  &                                      + e1e2t(ji  ,jj+1) * ssha_e(ji  ,jj+1) ) * ssvmask(ji,jj)
+            END_2D
+#endif
+         ENDIF
          !         
          ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2
@@ -624,8 +643,14 @@
                !                    ! hu_e, hv_e hold depth at jn,  zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2
                !                    ! backward interpolated depth used in spg terms at jn+1/2
+#if defined key_qcoTest_FluxForm
+            !                                ! 'key_qcoTest_FluxForm' : simple ssh average
+               zhu_bck = hu_0(ji,jj) + r1_2 * (  zsshp2_e(ji,jj) + zsshp2_e(ji+1,jj  )  ) * ssumask(ji,jj)
+               zhv_bck = hv_0(ji,jj) + r1_2 * (  zsshp2_e(ji,jj) + zsshp2_e(ji  ,jj+1)  ) * ssvmask(ji,jj)
+#else
                zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)    &
                     &                                          + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj)
                zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )    &
                     &                                          + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj)
+#endif
                !                    ! inverse depth at jn+1
                z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) )
@@ -646,14 +671,14 @@
          IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
             DO_2D( 0, 0, 0, 0 )
-                  ua_e(ji,jj) =  ua_e(ji,jj) /(1.0 -   rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj))
-                  va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj))
+               ua_e(ji,jj) =  ua_e(ji,jj) / ( 1._wp - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj) )
+               va_e(ji,jj) =  va_e(ji,jj) / ( 1._wp - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj) )
             END_2D
          ENDIF
        
          IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
-            hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
-            hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) )
-            hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)
-            hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) )
+            hu_e (2:jpim1,2:jpjm1) =    hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
+            hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / (  hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1)  )
+            hv_e (2:jpim1,2:jpjm1) =    hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)
+            hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / (  hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1)  )
          ENDIF
          !
@@ -743,17 +768,25 @@
       ELSE
          ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points
+#if defined key_qcoTest_FluxForm
+         !                                ! 'key_qcoTest_FluxForm' : simple ssh average
          DO_2D( 1, 0, 1, 0 )
-            zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj) &
-               &              * ( e1e2t(ji  ,jj) * pssh(ji  ,jj,Kaa)      &
-               &              +   e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) )
-            zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj) &
-               &              * ( e1e2t(ji,jj  ) * pssh(ji,jj  ,Kaa)      &
-               &              +   e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) )
-         END_2D
+            zsshu_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji+1,jj  ,Kaa) ) * ssumask(ji,jj)
+            zsshv_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji  ,jj+1,Kaa) ) * ssvmask(ji,jj)
+         END_2D
+#else
+         DO_2D( 1, 0, 1, 0 )
+            zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji  ,jj) * pssh(ji  ,jj,Kaa)   &
+               &                                      + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) * ssumask(ji,jj)
+            zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj  ) * pssh(ji,jj  ,Kaa)   &
+               &                                      + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj)
+         END_2D
+#endif   
          CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
          !
          DO jk=1,jpkm1
-            puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b
-            pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b
+            puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm)   &
+               &             * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b
+            pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm)   &
+               &             * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b
          END DO
          ! Save barotropic velocities not transport:
@@ -1101,14 +1134,14 @@
          CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
             DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                    &
-                    &           ht(ji  ,jj  ) + ht(ji+1,jj  )   ) * 0.25_wp  
+               zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1)   &
+                    &       + ht(ji,jj  ) + ht(ji+1,jj  ) ) * 0.25_wp  
                IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
             END_2D
          CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
             DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj) =             (  ht  (ji  ,jj+1) + ht  (ji+1,jj+1)      &
-                    &                    + ht  (ji  ,jj  ) + ht  (ji+1,jj  )  )   &
-                    &       / ( MAX( 1._wp,  ssmask(ji  ,jj+1) + ssmask(ji+1,jj+1)      &
-                    &                      + ssmask(ji  ,jj  ) + ssmask(ji+1,jj  )  )   )
+               zwz(ji,jj) =     (    ht(ji,jj+1) +     ht(ji+1,jj+1)      &
+                    &            +   ht(ji,jj  ) +     ht(ji+1,jj  )  )   &
+                    &    / ( MAX(ssmask(ji,jj+1) + ssmask(ji+1,jj+1)      &
+                    &          + ssmask(ji,jj  ) + ssmask(ji+1,jj  ) , 1._wp  )   )
                IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
             END_2D
@@ -1116,5 +1149,5 @@
          CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
          !
-         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
+         ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;   ftsw(1,:) = 0._wp
          DO_2D( 0, 1, 0, 1 )
             ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
@@ -1125,5 +1158,5 @@
          !
       CASE( np_EET )                  != EEN scheme using e3t energy conserving scheme
-         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
+         ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;   ftsw(1,:) = 0._wp
          DO_2D( 0, 1, 0, 1 )
             z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )
@@ -1167,9 +1200,9 @@
          ENDIF
          !
-         DO jj = 1, jpjm1
+         DO jj = 1, jpjm1   ! keep only the value at the coast if sco
             zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
          END DO
          !
-         DO jk = 1, jpkm1
+         DO jk = 1, jpkm1   ! ocean point : sum of masked e3f
             DO jj = 1, jpjm1
                zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/sshwzv.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/sshwzv.F90	(revision 13426)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/sshwzv.F90	(revision 13427)
@@ -299,4 +299,5 @@
          !                                                  ! filtered "now" field
          pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )
+         !
          IF( .NOT.ln_linssh ) THEN                          ! "now" <-- with forcing removed
             zcoef = rn_atfp * rn_Dt * r1_rho0
@@ -307,10 +308,11 @@
 
             ! ice sheet coupling
-            IF ( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1) pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:)
+            IF( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1 )   &
+               &   pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:)
 
          ENDIF
       ENDIF
       !
-      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' pssh(:,:,Kmm)  - : ', mask1=tmask )
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' atf  - pssh(:,:,Kmm): ', mask1=tmask )
       !
       IF( ln_timing )   CALL timing_stop('ssh_atf')
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/IOM/restart.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/IOM/restart.F90	(revision 13426)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/IOM/restart.F90	(revision 13427)
@@ -160,4 +160,14 @@
                      CALL iom_rstput( kt, nitrst, numrow, 'sshn'   ,ssh(:,:         ,Kmm), ldxios = lwxios      )
                      CALL iom_rstput( kt, nitrst, numrow, 'rhop'   , rhop, ldxios = lwxios      )
+!!st uncomment in case of dom_qco_read2 usage
+!!#if defined key_qco
+!!                     CALL iom_rstput( kt, nitrst, numrow, 'r3tb'   , r3t(:,:,Kbb) , ldxios = lwxios      )
+!!                     CALL iom_rstput( kt, nitrst, numrow, 'r3ub'   , r3u(:,:,Kbb) , ldxios = lwxios      )
+!!                     CALL iom_rstput( kt, nitrst, numrow, 'r3vb'   , r3v(:,:,Kbb) , ldxios = lwxios      )                     
+!!                     CALL iom_rstput( kt, nitrst, numrow, 'r3tn'   , r3t(:,:,Kmm) , ldxios = lwxios      )
+!!                     CALL iom_rstput( kt, nitrst, numrow, 'r3un'   , r3u(:,:,Kmm) , ldxios = lwxios      )
+!!                     CALL iom_rstput( kt, nitrst, numrow, 'r3vn'   , r3v(:,:,Kmm) , ldxios = lwxios      )
+!!                     CALL iom_rstput( kt, nitrst, numrow, 'r3f'   , r3f(:,:) , ldxios = lwxios      )
+!!#endif 
       ENDIF
       
@@ -275,5 +285,7 @@
          CALL iom_get( numror, jpdom_auto, 'tb'     , ts(:,:,:,jp_tem,Kbb), ldxios = lrxios )
          CALL iom_get( numror, jpdom_auto, 'sb'     , ts(:,:,:,jp_sal,Kbb), ldxios = lrxios )
+#if ! defined key_qco
          CALL iom_get( numror, jpdom_auto, 'sshb'   ,ssh(:,:         ,Kbb), ldxios = lrxios )
+#endif
       ELSE
          l_1st_euler =  .TRUE.      ! before field not found, forced euler 1st time-step
@@ -285,5 +297,7 @@
       CALL iom_get( numror, jpdom_auto, 'tn'     , ts(:,:,:,jp_tem,Kmm), ldxios = lrxios )
       CALL iom_get( numror, jpdom_auto, 'sn'     , ts(:,:,:,jp_sal,Kmm), ldxios = lrxios )
+#if ! defined key_qco
       CALL iom_get( numror, jpdom_auto, 'sshn'   ,ssh(:,:         ,Kmm), ldxios = lrxios )
+#endif
       IF( iom_varid( numror, 'rhop', ldstop = .FALSE. ) > 0 ) THEN
          CALL iom_get( numror, jpdom_auto, 'rhop'   , rhop, ldxios = lrxios )   ! now    potential density
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/SBC/sbcice_cice.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/SBC/sbcice_cice.F90	(revision 13426)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/SBC/sbcice_cice.F90	(revision 13427)
@@ -12,8 +12,8 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
-# if ! defined key_qco
-   USE domvvl
+# if defined key_qco
+   USE domqco         ! Variable volume
 # else
-   USE domqco
+   USE domvvl         ! Variable volume
 # endif
    USE phycst, only : rcp, rho0, r1_rho0, rhos, rhoi
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/ZDF/zdfddm.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/ZDF/zdfddm.F90	(revision 13426)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/ZDF/zdfddm.F90	(revision 13427)
@@ -31,4 +31,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/nemogcm.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/nemogcm.F90	(revision 13426)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/nemogcm.F90	(revision 13427)
@@ -183,9 +183,9 @@
       !
       DO WHILE( istp <= nitend .AND. nstop == 0 )
-#if defined key_qco
+#  if defined key_qco
          CALL stp_MLF
-#else
+#  else
          CALL stp
-#endif
+#  endif
          istp = istp + 1
       END DO
@@ -204,9 +204,9 @@
             ENDIF
             
-#if defined key_qco
+#  if defined key_qco
             CALL stp_MLF      ( istp )
-#else
+#  else
             CALL stp        ( istp ) 
-#endif
+#  endif
             istp = istp + 1
 
@@ -416,7 +416,7 @@
 
       ! Initialise time level indices
-      Nbb = 1; Nnn = 2; Naa = 3; Nrhs = Naa
+      Nbb = 1   ;   Nnn = 2   ;   Naa = 3   ;   Nrhs = Naa
 #if defined key_agrif
-      Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
+      Kbb_a = Nbb   ;   Kmm_a = Nnn   ;   Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
 #endif 
       !                             !-------------------------------!
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/stpMLF.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/stpMLF.F90	(revision 13426)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/stpMLF.F90	(revision 13427)
@@ -32,8 +32,14 @@
    !!            4.0  !  2017-05  (G. Madec)  introduction of the vertical physics manager (zdfphy)
    !!            4.1  !  2019-08  (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   stp_MLF       : OPA system time-stepping
+   !!            4.x  !  2020-08  (S. Techene, G. Madec)  quasi eulerian coordinate time stepping 
+   !!----------------------------------------------------------------------
+
+#if defined key_qco
+   !!----------------------------------------------------------------------
+   !!   'key_qco'       Quasi-Eulerian vertical coordonate
+   !!----------------------------------------------------------------------
+   
+   !!----------------------------------------------------------------------
+   !!   stp_MLF       : NEMO modified Leap Frog time-stepping with qco
    !!----------------------------------------------------------------------
    USE step_oce       ! time stepping definition modules
@@ -43,5 +49,5 @@
    USE traatfqco      ! time filtering                   (tra_atf_qco routine)
    USE dynatfqco      ! time filtering                   (dyn_atf_qco routine)
-   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme (define un_adv)
+   
    USE bdydyn         ! ocean open boundary conditions (define bdy_dyn)
 
@@ -49,9 +55,4 @@
    PRIVATE
    
-#if ! defined key_qco
-   !!----------------------------------------------------------------------
-   !!   'key_qco'      EMPTY MODULE      Quasi-Eulerian vertical coordonate
-   !!----------------------------------------------------------------------
-#else
    PUBLIC   stp_MLF   ! called by nemogcm.F90
 
@@ -195,9 +196,12 @@
          zgdept(:,:,jk) = gdept(:,:,jk,Nnn)
       END DO
-                            CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh,  Naa )   ! after ssh (includes call to div_hor)
-      IF( .NOT.ln_linssh )  CALL dom_qco_r3c    ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa) )   ! "after" ssh./h._0 ratio
-                            CALL wzv           ( kstp, Nbb, Nnn, Naa, ww  )    ! Nnn cross-level velocity
-      IF( ln_zad_Aimp )     CALL wAimp         ( kstp,      Nnn           )    ! Adaptive-implicit vertical advection partitioning
-                            CALL eos    ( ts(:,:,:,:,Nnn), rhd, rhop, zgdept ) ! now in situ density for hpg computation
+                            CALL ssh_nxt    ( kstp, Nbb, Nnn, ssh,  Naa )   ! after ssh (includes call to div_hor)
+      IF( .NOT.ln_linssh )  THEN
+                             CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa)           )   ! "after" ssh/h_0 ratio at t,u,v pts
+         IF( ln_dynspg_exp ) CALL dom_qco_r3c( ssh(:,:,Nnn), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn), r3f(:,:) )   ! spg_exp : needed only for "now" ssh/h_0 ratio at f point
+      ENDIF
+                            CALL wzv        ( kstp, Nbb, Nnn, Naa, ww  )    ! Nnn cross-level velocity
+      IF( ln_zad_Aimp )     CALL wAimp      ( kstp,      Nnn           )    ! Adaptive-implicit vertical advection partitioning
+                            CALL eos        ( ts(:,:,:,:,Nnn), rhd, rhop, zgdept ) ! now in situ density for hpg computation
 
 
@@ -218,9 +222,9 @@
                          CALL dyn_hpg( kstp,      Nnn      , uu, vv, Nrhs )  ! horizontal gradient of Hydrostatic pressure
                          CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa )  ! surface pressure gradient
-
                                                       ! With split-explicit free surface, since now transports have been updated and ssh(:,:,Nrhs) as well
+
       IF( ln_dynspg_ts ) THEN                         ! vertical scale factors and vertical velocity need to be updated
                             CALL div_hor    ( kstp, Nbb, Nnn )                ! Horizontal divergence  (2nd call in time-split case)
-         IF(.NOT.ln_linssh) CALL dom_qco_r3c ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! "after" ssh./h._0 ratio
+         IF(.NOT.ln_linssh) CALL dom_qco_r3c ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! update ssh/h_0 ratio at t,u,v,f pts 
       ENDIF
                             CALL dyn_zdf    ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion
@@ -305,9 +309,9 @@
 !!    place.
 !!
-                         CALL mlf_baro_corr ( Nnn, Naa, uu, vv )                    ! barotrope ajustment
-                         CALL finalize_sbc  ( kstp, Nbb, Naa, uu, vv, ts )          ! boundary condifions
-                         CALL tra_atf_qco   ( kstp, Nbb, Nnn, Naa, ts )             ! time filtering of "now" tracer arrays
-                         CALL dyn_atf_qco   ( kstp, Nbb, Nnn, Naa, uu, vv  )        ! time filtering of "now" velocities and scale factors
-                         r3t(:,:,Nnn) = r3t_f(:,:)
+      IF( ln_dynspg_ts ) CALL mlf_baro_corr (            Nnn, Naa, uu, vv     )   ! barotrope adjustment
+                         CALL finalize_lbc  ( kstp, Nbb     , Naa, uu, vv, ts )   ! boundary conditions
+                         CALL tra_atf_qco   ( kstp, Nbb, Nnn, Naa        , ts )   ! time filtering of "now" tracer arrays
+                         CALL dyn_atf_qco   ( kstp, Nbb, Nnn, Naa, uu, vv     )   ! time filtering of "now" velocities 
+                         r3t(:,:,Nnn) = r3t_f(:,:)                                ! update now ssh/h_0 with time filtered values
                          r3u(:,:,Nnn) = r3u_f(:,:)
                          r3v(:,:,Nnn) = r3v_f(:,:)
@@ -347,5 +351,5 @@
       ! Control
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL stp_ctl      ( kstp, Nbb, Nnn, indic )
+                         CALL stp_ctl      ( kstp, Nnn )
 
       IF( kstp == nit000 ) THEN                          ! 1st time step only
@@ -364,5 +368,5 @@
       IF( kstp == nitend .OR. indic < 0 ) THEN
                       CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
-         IF(lrxios) CALL iom_context_finalize(      crxios_context          )
+         IF( lrxios ) CALL iom_context_finalize(      crxios_context          )
          IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !
       ENDIF
@@ -380,5 +384,5 @@
 
 
-   SUBROUTINE mlf_baro_corr (Kmm, Kaa, puu, pvv)
+   SUBROUTINE mlf_baro_corr( Kmm, Kaa, puu, pvv )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE mlf_baro_corr  ***
@@ -389,42 +393,38 @@
       !!             estimate (ln_dynspg_ts=T)
       !!
-      !! ** Action :   puu(Kmm),pvv(Kmm),puu(Kaa),pvv(Kaa)   now and after horizontal velocity
-      !!----------------------------------------------------------------------
-      INTEGER                             , INTENT(in   ) :: Kmm, Kaa    ! before and after time level indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv         ! velocities
-      !
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
+      !! ** Action :   puu(Kmm),pvv(Kmm)   updated now horizontal velocity (ln_bt_fw=F)
+      !!               puu(Kaa),pvv(Kaa)   after horizontal velocity
+      !!----------------------------------------------------------------------
+      USE dynspg_ts, ONLY : un_adv, vn_adv   ! updated Kmm barotropic transport 
+      !!
+      INTEGER                             , INTENT(in   ) ::   Kmm, Kaa   ! before and after time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv   ! velocities
       !
       INTEGER  ::   jk   ! dummy loop indices
-      !!----------------------------------------------------------------------
-
-      IF ( ln_dynspg_ts ) THEN
-         ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
-         ! Ensure below that barotropic velocities match time splitting estimate
-         ! Compute actual transport and replace it with ts estimate at "after" time step
-         zue(:,:) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
-         zve(:,:) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
-         DO jk = 2, jpkm1
-            zue(:,:) = zue(:,:) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
-            zve(:,:) = zve(:,:) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
+      REAL(wp), DIMENSION(jpi,jpj) ::   zue, zve
+      !!----------------------------------------------------------------------
+
+      ! Ensure below that barotropic velocities match time splitting estimate
+      ! Compute actual transport and replace it with ts estimate at "after" time step
+      zue(:,:) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
+      zve(:,:) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
+      DO jk = 2, jpkm1
+         zue(:,:) = zue(:,:) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
+         zve(:,:) = zve(:,:) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
+      END DO
+      DO jk = 1, jpkm1
+         puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk)
+         pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk)
+      END DO
+      !
+      IF( .NOT.ln_bt_fw ) THEN
+         ! Remove advective velocity from "now velocities"
+         ! prior to asselin filtering
+         ! In the forward case, this is done below after asselin filtering
+         ! so that asselin contribution is removed at the same time
+         DO jk = 1, jpkm1
+            puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk)
+            pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk)
          END DO
-         DO jk = 1, jpkm1
-            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk)
-            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk)
-         END DO
-         !
-         IF( .NOT.ln_bt_fw ) THEN
-            ! Remove advective velocity from "now velocities"
-            ! prior to asselin filtering
-            ! In the forward case, this is done below after asselin filtering
-            ! so that asselin contribution is removed at the same time
-            DO jk = 1, jpkm1
-               puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk)
-               pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk)
-            END DO
-         ENDIF
-         !
-         DEALLOCATE( zue, zve )
-         !
       ENDIF
       !
@@ -432,7 +432,7 @@
 
 
-   SUBROUTINE finalize_sbc (kt, Kbb, Kaa, puu, pvv, pts)
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE finalize_sbc  ***
+   SUBROUTINE finalize_lbc( kt, Kbb, Kaa, puu, pvv, pts )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE finalize_lbc  ***
       !!
       !! ** Purpose :   Apply the boundary condition on the after velocity
@@ -445,18 +445,19 @@
       !! ** Action :   puu(Kaa),pvv(Kaa)   after horizontal velocity and tracers
       !!----------------------------------------------------------------------
-      INTEGER                             , INTENT(in   ) :: kt               ! ocean time-step index
-      INTEGER                             , INTENT(in   ) :: Kbb, 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,jpts,jpt), INTENT(inout) :: pts            ! active tracers
+      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) ::   Kbb, 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,jpts,jpt), INTENT(inout) ::   pts        ! active tracers
+      !!----------------------------------------------------------------------
       !
       ! Update after tracer and velocity on domain lateral boundaries
       !
-#if defined key_agrif
-            CALL Agrif_tra                     ! AGRIF zoom boundaries
-            CALL Agrif_dyn( kt )               !* AGRIF zoom boundaries
-#endif
+# if defined key_agrif
+            CALL Agrif_tra                     !* AGRIF zoom boundaries
+            CALL Agrif_dyn( kt )
+# endif
       !                                        ! local domain boundaries  (T-point, unchanged sign)
-      CALL lbc_lnk_multi( 'finalize_sbc', puu(:,:,:,       Kaa), 'U', -1., pvv(:,:,:       ,Kaa), 'V', -1.   &
-                       &                , pts(:,:,:,jp_tem,Kaa), 'T',  1., pts(:,:,:,jp_sal,Kaa), 'T',  1. )     !* local domain boundaries
+      CALL lbc_lnk_multi( 'finalize_lbc', puu(:,:,:,       Kaa), 'U', -1., pvv(:,:,:       ,Kaa), 'V', -1.   &
+                       &                , pts(:,:,:,jp_tem,Kaa), 'T',  1., pts(:,:,:,jp_sal,Kaa), 'T',  1. )
       !
       !                                        !* BDY open boundaries
@@ -467,7 +468,12 @@
       ENDIF
       !
-   END SUBROUTINE finalize_sbc
-#endif
-   !
+   END SUBROUTINE finalize_lbc
+
+#else
+   !!----------------------------------------------------------------------
+   !!   default option             EMPTY MODULE           qco not activated
+   !!----------------------------------------------------------------------
+#endif
+   
    !!======================================================================
 END MODULE stepMLF
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OFF/dtadyn.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OFF/dtadyn.F90	(revision 13426)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OFF/dtadyn.F90	(revision 13427)
@@ -23,8 +23,8 @@
    USE c1d             ! 1D configuration: lk_c1d
    USE dom_oce         ! ocean domain: variables
-#if ! defined key_qco 
-   USE domvvl          ! variable volume
+#if defined key_qco 
+   USE domqco          ! variable volume
 #else
-   USE domqco
+   USE domvvl
 #endif
    USE zdf_oce         ! ocean vertical physics: variables
@@ -97,4 +97,6 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+   
    !!----------------------------------------------------------------------
    !! NEMO/OFF 4.0 , NEMO Consortium (2018)
@@ -385,6 +387,6 @@
         gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm)
         !
-      ENDIF
 #endif
+      ENDIF
       !
       IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN       ! read depht over which runoffs are distributed
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OFF/nemogcm.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OFF/nemogcm.F90	(revision 13426)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OFF/nemogcm.F90	(revision 13427)
@@ -64,5 +64,9 @@
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
    USE lbcnfd  , ONLY : isendto, nsndto   ! Setup of north fold exchanges
-   USE step, ONLY : Nbb, Nnn, Naa, Nrhs   ! time level indices
+#if defined key_qco
+   USE stepMLF , ONLY : Nbb, Nnn, Naa, Nrhs   ! time level indices
+#else
+   USE step    , ONLY : Nbb, Nnn, Naa, Nrhs   ! time level indices
+#endif
    USE halo_mng
 
@@ -126,12 +130,12 @@
                                 CALL dta_dyn_atf( istp, Nbb, Nnn, Naa )       ! time filter of sea  surface height and vertical scale factors
 # if defined key_qco
-                                CALL dom_qco_r3c( ssh(:,:,Kmm), r3t_f, r3u_f, r3v_f )
+                                CALL dom_qco_r3c( ssh(:,:,Nnn), r3t_f, r3u_f, r3v_f )
 # endif
          ENDIF
                                 CALL trc_stp    ( istp, Nbb, Nnn, Nrhs, Naa ) ! time-stepping
 # if defined key_qco
-                                !r3t(:,:,Kmm) = r3t_f(:,:)                     ! update ssh to h0 ratio
-                                !r3u(:,:,Kmm) = r3u_f(:,:)
-                                !r3v(:,:,Kmm) = r3v_f(:,:)
+                                !r3t(:,:,Nnn) = r3t_f(:,:)                     ! update ssh to h0 ratio
+                                !r3u(:,:,Nnn) = r3u_f(:,:)
+                                !r3v(:,:,Nnn) = r3v_f(:,:)
 # endif
 #endif
@@ -143,7 +147,7 @@
          !
 #if ! defined key_qco
-#if ! defined key_sed_off
+# if ! defined key_sed_off
          IF( .NOT.ln_linssh )   CALL dta_dyn_sf_interp( istp, Nnn )  ! calculate now grid parameters
-#endif
+# endif
 #endif         
                                 CALL stp_ctl    ( istp )             ! Time loop: control and print
