Index: NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/dom_oce.F90
===================================================================
--- NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/dom_oce.F90	(revision 14046)
+++ NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/dom_oce.F90	(revision 14058)
@@ -131,6 +131,6 @@
    !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:) ::   e1e2t , r1_e1e2t                !: associated metrics at t-point
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:) ::   e1e2u , e2_e1u, r1_e1e2u        !: associated metrics at u-point
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:) ::   e1e2v , e1_e2v, r1_e1e2v        !: associated metrics at v-point
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:) ::   e1e2u , r1_e1e2u , e2_e1u       !: associated metrics at u-point
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:) ::   e1e2v , r1_e1e2v , e1_e2v       !: associated metrics at v-point
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:) ::   e1e2f , r1_e1e2f                !: associated metrics at f-point
    !
@@ -162,7 +162,7 @@
 
    !                                                        !  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  
@@ -205,7 +205,9 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)           ::   ssmask, ssumask, ssvmask, ssfmask   !: surface mask at T-,U-, V- and F-pts
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, wmask, fmask   !: land/ocean mask at T-, U-, V-, W- and F-pts
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wumask, wvmask        !: land/ocean mask at WT-, WU- and WV-pts
-
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET ::   tmask, umask, vmask, wmask, fmask   !: land/ocean mask at T-, U-, V-, W- and F-pts
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET ::   wumask, wvmask                      !: land/ocean mask at WU- and WV-pts
+#if defined key_qco   
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET ::   fe3mask                             !: land/ocean mask at F-pts for qco
+#endif
    !!----------------------------------------------------------------------
    !! calendar variables
@@ -306,5 +308,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) ,      &
@@ -313,8 +319,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)  )
@@ -323,8 +325,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
@@ -350,4 +348,9 @@
       ii = ii+1
       ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(ii) )
+#if defined key_qco
+         !
+      ii = ii+1
+      ALLOCATE( fe3mask(jpi,jpj,jpk) , STAT=ierr(ii) )
+#endif
       !
       dom_oce_alloc = MAXVAL(ierr)
Index: NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/domain.F90
===================================================================
--- NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/domain.F90	(revision 14046)
+++ NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/domain.F90	(revision 14058)
@@ -15,5 +15,5 @@
    !!            3.7  !  2015-11  (G. Madec, A. Coward)  time varying zgr by default
    !!            4.0  !  2016-10  (G. Madec, S. Flavoni)  domain configuration / user defined interface
-   !!            4.x  ! 2020-02  (G. Madec, S. Techene) introduce ssh to h0 ratio
+   !!            4.1  !  2020-02  (G. Madec, S. Techene)  introduce ssh to h0 ratio
    !!----------------------------------------------------------------------
    
@@ -28,4 +28,10 @@
    USE oce            ! ocean variables
    USE dom_oce        ! domain: ocean
+#if defined key_qco
+   USE domqco         ! quasi-eulerian
+#else
+   USE domvvl         ! variable volume
+#endif
+   USE sshwzv  , ONLY : ssh_init_rst   ! set initial ssh 
    USE sbc_oce        ! surface boundary condition: ocean
    USE trc_oce        ! shared ocean & passive tracers variab
@@ -35,13 +41,8 @@
    USE dommsk         ! domain: set the mask system
    USE domwri         ! domain: write the meshmask file
-#if ! defined key_qco
-   USE domvvl         ! variable volume
-#else
-   USE domqco          ! variable volume
-#endif
    USE c1d            ! 1D configuration
    USE dyncor_c1d     ! 1D configuration: Coriolis term    (cor_c1d routine)
-   USE wet_dry, ONLY : ll_wd
-   USE closea , ONLY : dom_clo ! closed seas
+   USE wet_dry , ONLY : ll_wd     ! wet & drying flag
+   USE closea  , ONLY : dom_clo   ! closed seas routine
    !
    USE prtctl         ! Print control (prt_ctl_info routine)
@@ -50,4 +51,5 @@
    USE lbclnk         ! ocean lateral boundary condition (or mpp link)
    USE lib_mpp        ! distributed memory computing library
+   USE restart        ! only for lrst_oce
 
    IMPLICIT NONE
@@ -58,4 +60,6 @@
    PUBLIC   dom_tile     ! called by step.F90
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!-------------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -84,4 +88,5 @@
       INTEGER ::   ji, jj, jk, jt   ! dummy loop indices
       INTEGER ::   iconf = 0    ! local integers
+      REAL(wp)::   zrdt
       CHARACTER (len=64) ::   cform = "(A12, 3(A13, I7))" 
       INTEGER , DIMENSION(jpi,jpj) ::   ik_top , ik_bot       ! top and bottom ocean level
@@ -121,6 +126,5 @@
          WRITE(numout,*)     '         cn_cfg = ', TRIM( cn_cfg ), '   nn_cfg = ', nn_cfg
       ENDIF
-      nn_wxios = 0
-      ln_xios_read = .FALSE.
+      
       !
       !           !==  Reference coordinate system  ==!
@@ -143,10 +147,18 @@
       hv_0(:,:) = 0._wp
       hf_0(:,:) = 0._wp
-      DO jk = 1, jpk
+      DO jk = 1, jpkm1
          ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
          hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
          hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
-         hf_0(:,:) = hf_0(:,:) + e3f_0(:,:,jk) * fmask(:,:,jk)
       END DO
+      !
+      DO jk = 1, jpkm1
+         hf_0(1:jpim1,:) = hf_0(1:jpim1,:) + e3f_0(1:jpim1,:,jk)*vmask(1:jpim1,:,jk)*vmask(2:jpi,:,jk)
+      END DO
+      CALL lbc_lnk('domain', hf_0, 'F', 1._wp)
+      !
+      IF( lk_SWE ) THEN      ! SWE case redefine hf_0
+         hf_0(:,:) = hf_0(:,:) + e3f_0(:,:,1) * ssfmask(:,:)
+      ENDIF
       !
       r1_ht_0(:,:) = ssmask (:,:) / ( ht_0(:,:) + 1._wp -  ssmask (:,:) )
@@ -154,17 +166,29 @@
       r1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp -  ssvmask(:,:) )
       r1_hf_0(:,:) = ssfmask(:,:) / ( hf_0(:,:) + 1._wp -  ssfmask(:,:) )
-
+      !
+      IF( ll_wd ) THEN       ! wet and drying (check ht_0 >= 0)
+         DO_2D( 1, 1, 1, 1 )
+            IF( ht_0(ji,jj) < 0._wp .AND. ssmask(ji,jj) == 1._wp ) THEN
+               CALL ctl_stop( 'ssh_init_rst : ht_0 must be positive at potentially wet points' )
+            ENDIF
+         END_2D
+      ENDIF
+      !
+      !           !==  initialisation of time varying coordinate  ==!
+      !
+      !                                 != ssh initialization
+      IF( .NOT.l_offline .AND. .NOT.l_SAS ) THEN
+         CALL ssh_init_rst( Kbb, Kmm, Kaa )
+      ELSE
+         ssh(:,:,:) = 0._wp
+      ENDIF
       !
 #if defined key_qco
-      !           !==  initialisation of time varying coordinate  ==!   Quasi-Euerian coordinate case
+      !                                 != 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
+      !
+      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
@@ -175,22 +199,22 @@
          !
          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(:,:)
-         !
-      ELSE                       != time varying : initialize before/now/after variables
+               ht   (:,:) =    ht_0(:,:)
+         !
+      ELSE                              != Time varying : initialize before/now/after variables
          !
          IF( .NOT.l_offline )   CALL dom_vvl_init( Kbb, Kmm, Kaa )
@@ -373,5 +397,7 @@
       USE ioipsl
       !!
-      INTEGER  ::   ios   ! Local integer
+      INTEGER ::   ios   ! Local integer
+      REAL(wp)::   zrdt
+      !!----------------------------------------------------------------------
       !
       NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                 &
@@ -393,4 +419,46 @@
       ENDIF
       !
+      !                       !=======================!
+      !                       !==  namelist namdom  ==!
+      !                       !=======================!
+      !
+      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
+903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdom in reference namelist' )
+      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
+904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdom in configuration namelist' )
+      IF(lwm) WRITE( numond, namdom )
+      !
+#if defined key_agrif
+      IF( .NOT. Agrif_Root() ) THEN    ! AGRIF child, subdivide the Parent timestep
+         rn_Dt = Agrif_Parent (rn_Dt ) / Agrif_Rhot()
+      ENDIF
+#endif
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '   Namelist : namdom   ---   space & time domain'
+         WRITE(numout,*) '      linear free surface (=T)                ln_linssh   = ', ln_linssh
+         WRITE(numout,*) '      create mesh/mask file                   ln_meshmask = ', ln_meshmask
+         WRITE(numout,*) '      ocean time step                         rn_Dt       = ', rn_Dt
+         WRITE(numout,*) '      asselin time filter parameter           rn_atfp     = ', rn_atfp
+         WRITE(numout,*) '      online coarsening of dynamical fields   ln_crs      = ', ln_crs
+      ENDIF
+      !
+      ! set current model timestep rDt = 2*rn_Dt if MLF or rDt = rn_Dt if RK3
+      rDt   = 2._wp * rn_Dt
+      r1_Dt = 1._wp / rDt
+      !
+      IF( l_SAS .AND. .NOT.ln_linssh ) THEN
+         CALL ctl_warn( 'SAS requires linear ssh : force ln_linssh = T' )
+         ln_linssh = .TRUE.
+      ENDIF
+      !
+#if defined key_qco
+      IF( ln_linssh )   CALL ctl_stop( 'STOP','domain: key_qco and ln_linssh = T are incompatible' )
+#endif
+      !
+      !                       !=======================!
+      !                       !==  namelist namrun  ==!
+      !                       !=======================!
       !
       READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
@@ -452,13 +520,51 @@
       nleapy = nn_leapy
       ninist = nn_istate
+      !
+      !                                        !==  Set parameters for restart reading using xIOS  ==!
+      !
+      IF( TRIM(Agrif_CFixed()) == '0' ) THEN
+         lrxios = ln_xios_read .AND. ln_rstart
+         IF( nn_wxios > 0 )   lwxios = .TRUE.           !* set output file type for XIOS based on NEMO namelist
+         nxioso = nn_wxios
+      ENDIF
+      !                                        !==  Check consistency between ln_rstart and ln_1st_euler  ==!   (i.e. set l_1st_euler)
       l_1st_euler = ln_1st_euler
-      IF( .NOT. l_1st_euler .AND. .NOT. ln_rstart ) THEN
+      !
+      IF( ln_rstart ) THEN                              !*  Restart case
+         !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   open the restart file'
+         CALL rst_read_open                                              !- Open the restart file
+         !
+         IF( iom_varid( numror, 'rdt', ldstop = .FALSE. ) > 0 ) THEN     !- Check time-step consistency and force Euler restart if changed
+            CALL iom_get( numror, 'rdt', zrdt )
+            IF( zrdt /= rn_Dt ) THEN
+               IF(lwp) WRITE( numout,*)
+               IF(lwp) WRITE( numout,*) '   rn_Dt = ', rn_Dt,' not equal to the READ one rdt = ', zrdt
+               IF(lwp) WRITE( numout,*)
+               IF(lwp) WRITE( numout,*) '      ==>>>   forced euler first time-step'
+               l_1st_euler =  .TRUE.
+            ENDIF
+         ENDIF
+         !
+         IF( .NOT.l_SAS .AND. iom_varid( numror, 'sshb', ldstop = .FALSE. ) <= 0 ) THEN   !- Check absence of one of the Kbb field (here sshb)
+            !                                                                             !  (any Kbb field is missing ==> all Kbb fields are missing) 
+            IF( .NOT.l_1st_euler ) THEN
+               CALL ctl_warn('dom_nam : ssh at Kbb not found in restart files ',   &
+                  &                        'l_1st_euler forced to .true. and ' ,   &
+                  &                        'ssh(Kbb) = ssh(Kmm) '                  )
+               l_1st_euler = .TRUE.
+            ENDIF
+         ENDIF
+      ELSEIF( .NOT.l_1st_euler ) THEN                   !*  Initialization case
          IF(lwp) WRITE(numout,*)  
          IF(lwp) WRITE(numout,*)'   ==>>>   Start from rest (ln_rstart=F)'
          IF(lwp) WRITE(numout,*)'           an Euler initial time step is used : l_1st_euler is forced to .true. '   
-         l_1st_euler = .true.
-      ENDIF
-      !                             ! control of output frequency
-      IF( .NOT. ln_rst_list ) THEN     ! we use nn_stock
+         l_1st_euler = .TRUE.
+      ENDIF
+      !
+      !                                        !==  control of output frequency  ==!
+      !
+      IF( .NOT. ln_rst_list ) THEN   ! we use nn_stock
          IF( nn_stock == -1 )   CALL ctl_warn( 'nn_stock = -1 --> no restart will be done' )
          IF( nn_stock == 0 .OR. nn_stock > nitend ) THEN
@@ -479,5 +585,5 @@
       IF( Agrif_Root() ) THEN
          IF(lwp) WRITE(numout,*)
-         SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
+         SELECT CASE ( nleapy )                !==  Choose calendar for IOIPSL  ==!
          CASE (  1 ) 
             CALL ioconf_calendar('gregorian')
@@ -491,31 +597,9 @@
          END SELECT
       ENDIF
-
-      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
-903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdom in reference namelist' )
-      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
-904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdom in configuration namelist' )
-      IF(lwm) WRITE( numond, namdom )
-      !
-#if defined key_agrif
-      IF( .NOT. Agrif_Root() ) THEN
-            rn_Dt = Agrif_Parent(rn_Dt) / Agrif_Rhot()
-      ENDIF
-#endif
-      !
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) '   Namelist : namdom   ---   space & time domain'
-         WRITE(numout,*) '      linear free surface (=T)                ln_linssh   = ', ln_linssh
-         WRITE(numout,*) '      create mesh/mask file                   ln_meshmask = ', ln_meshmask
-         WRITE(numout,*) '      ocean time step                         rn_Dt       = ', rn_Dt
-         WRITE(numout,*) '      asselin time filter parameter           rn_atfp     = ', rn_atfp
-         WRITE(numout,*) '      online coarsening of dynamical fields   ln_crs      = ', ln_crs
-      ENDIF
-      !
-      !! Initialise current model timestep rDt = 2*rn_Dt if MLF or rDt = rn_Dt if RK3
-      rDt  = 2._wp * rn_Dt
-      r1_Dt = 1._wp / rDt
-
+      !
+      !                       !========================!
+      !                       !==  namelist namtile  ==!
+      !                       !========================!
+      !
       READ  ( numnam_ref, namtile, IOSTAT = ios, ERR = 905 )
 905   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtile in reference namelist' )
@@ -537,14 +621,10 @@
          ENDIF
       ENDIF
-
-      IF( TRIM(Agrif_CFixed()) == '0' ) THEN
-         lrxios = ln_xios_read.AND.ln_rstart
-!set output file type for XIOS based on NEMO namelist 
-         IF (nn_wxios > 0) lwxios = .TRUE. 
-         nxioso = nn_wxios
-      ENDIF
-
+      !
 #if defined key_netcdf4
-      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
+      !                       !=======================!
+      !                       !==  namelist namnc4  ==!   NetCDF 4 case   ("key_netcdf4" defined)
+      !                       !=======================!
+      !
       READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
 907   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namnc4 in reference namelist' )
@@ -555,5 +635,5 @@
       IF(lwp) THEN                        ! control print
          WRITE(numout,*)
-         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
+         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters ("key_netcdf4" defined)'
          WRITE(numout,*) '      number of chunks in i-dimension             nn_nchunks_i = ', nn_nchunks_i
          WRITE(numout,*) '      number of chunks in j-dimension             nn_nchunks_j = ', nn_nchunks_j
@@ -618,5 +698,5 @@
    SUBROUTINE domain_cfg( cd_cfg, kk_cfg, kpi, kpj, kpk, kperio )
       !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE dom_nam  ***
+      !!                     ***  ROUTINE domain_cfg  ***
       !!                    
       !! ** Purpose :   read the domain size in domain configuration file
Index: NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/dommsk.F90
===================================================================
--- NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/dommsk.F90	(revision 14046)
+++ NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/dommsk.F90	(revision 14058)
@@ -181,4 +181,14 @@
       ssvmask(:,:) = MAXVAL( vmask(:,:,:), DIM=3 )
       ssfmask(:,:) = MAXVAL( fmask(:,:,:), DIM=3 )
+      IF( lk_SWE ) THEN      ! Shallow Water Eq. case : redefine ssfmask
+         DO_2D( 0,0 , 0,0 )
+            ssfmask(ji,jj) = MAX(  ssmask(ji,jj+1), ssmask(ji+1,jj+1),  & 
+               &                   ssmask(ji,jj  ), ssmask(ji+1,jj  )   )
+         END_2D
+         CALL lbc_lnk( 'dommsk', ssfmask, 'F', 1.0_wp )
+      ENDIF
+#if defined key_qco
+      fe3mask(:,:,:) = fmask(:,:,:)
+#endif
 
       ! Interior domain mask  (used for global sum)
Index: NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/domqco.F90
===================================================================
--- NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/domqco.F90	(revision 14046)
+++ NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/domqco.F90	(revision 14058)
@@ -8,13 +8,13 @@
    !!            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
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   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
+   !!            4.1  !  2019-08  (A. Coward, D. Storkey) add time level indices for prognostic variables
+   !!             -   !  2020-02  (S. Techene, G. Madec) quasi-eulerian coordinate (z* or s*)
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   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_ctl   : Check the vvl options
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers
@@ -55,6 +55,4 @@
    LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE.                ! debug control prints
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                ! thickness diffusion transport
-
    !! * Substitutions
 #  include "do_loop_substitute.h90"
@@ -79,5 +77,6 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
+      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! time level indices
+      !!----------------------------------------------------------------------
       !
       IF(lwp) WRITE(numout,*)
@@ -85,33 +84,37 @@
       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
+      CALL qco_ctl                            ! choose vertical coordinate (z_star, z_tilde or layer)
+      !
+      CALL dom_qco_zgr( Kbb, Kmm )            ! interpolation scale factor, depth and water column
+      !
+#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_init
 
 
-   SUBROUTINE dom_qco_zgr(Kbb, Kmm, Kaa)
+   SUBROUTINE dom_qco_zgr( Kbb, Kmm )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE dom_qco_init  ***
       !!
-      !! ** Purpose :  Initialization of all ssh. to h._0 ratio
-      !!
-      !! ** Method  :  - interpolate scale factors
-      !!
-      !! ** 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
+      !! ** Purpose :  Initialization of all r3. = ssh./h._0 ratios
+      !!
+      !! ** Method  :  Call domqco using Kbb and Kmm
+      !!               NB: dom_qco_zgr is called by dom_qco_init it uses ssh from ssh_init 
+      !!
+      !! ** Action  : - r3(t/u/v)(Kbb)
+      !!              - r3(t/u/v/f)(Kmm)
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! time level indices
       !!----------------------------------------------------------------------
       !
       !                    !== Set of all other vertical scale factors  ==!  (now and before)
       !                                ! Horizontal interpolation of e3t
-      CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb) )
+      CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb)           )
       CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) )
       !
@@ -143,5 +146,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)  &
@@ -150,10 +155,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
@@ -163,6 +170,9 @@
       ELSE                                   !==  ratio at f-point  ==!
          !
-         IF( ln_dynadv_vec )   THEN                !- Vector Form   (thickness weighted averaging)
-            DO_2D( 1, 0, 1, 0 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line
+!!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 )                               ! 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  )  &
                   &                     + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )  &
@@ -170,10 +180,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)
-            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)
+!!st         ELSE                                      !- Flux Form   (simple averaging)
+#else
+            DO_2D( 0, 0, 0, 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)
             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 )
@@ -184,113 +196,7 @@
 
 
-   SUBROUTINE qe_rst_read( kt, Kbb, Kmm )
+   SUBROUTINE qco_ctl
       !!---------------------------------------------------------------------
-      !!                   ***  ROUTINE qe_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)    )
-               CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm)    )
-               ! 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) )
-               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) )
-               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
-               !
-               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
-               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 qe_rst_read
-
-
-   SUBROUTINE dom_qco_ctl
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE dom_qco_ctl  ***
+      !!                  ***  ROUTINE qco_ctl  ***
       !!
       !! ** Purpose :   Control the consistency between namelist options
@@ -312,6 +218,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
@@ -357,5 +263,5 @@
 #endif
       !
-   END SUBROUTINE dom_qco_ctl
+   END SUBROUTINE qco_ctl
 
    !!======================================================================
Index: NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/domvvl.F90
===================================================================
--- NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/domvvl.F90	(revision 14046)
+++ NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/domvvl.F90	(revision 14058)
@@ -9,5 +9,5 @@
    !!            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) introduce ssh to h0 ratio
+   !!             -   !  2020-02  (G. Madec, S. Techene) introduce ssh to h0 ratio
    !!----------------------------------------------------------------------
 
@@ -766,9 +766,18 @@
       !! ** Purpose :   Read or write VVL file in restart file
       !!
-      !! ** Method  :   use of IOM library
-      !!                if the restart does not contain vertical scale factors,
-      !!                they are set to the _0 values
-      !!                if the restart does not contain vertical scale factors increments (z_tilde),
-      !!                they are set to 0.
+      !! ** Method  : * restart comes from a linear ssh simulation :
+      !!                   an attempt to read e3t_n stops simulation
+      !!              * restart comes from a z-star, z-tilde, or layer :
+      !!                   read e3t_n and e3t_b
+      !!              * restart comes from a z-star :
+      !!                   set tilde_e3t_n, tilde_e3t_n, and hdiv_lf to 0
+      !!              * restart comes from layer :
+      !!                   read tilde_e3t_n and tilde_e3t_b
+      !!                   set hdiv_lf to 0
+      !!              * restart comes from a z-tilde:
+      !!                   read tilde_e3t_n, tilde_e3t_b, and hdiv_lf
+      !!
+      !!              NB: if l_1st_euler = T (ln_1st_euler or ssh_b not found)
+      !!                   Kbb fields set to Kmm ones
       !!----------------------------------------------------------------------
       INTEGER         , INTENT(in) ::   kt        ! ocean time-step
@@ -776,154 +785,91 @@
       CHARACTER(len=*), INTENT(in) ::   cdrw      ! "READ"/"WRITE" flag
       !
-      INTEGER ::   ji, jj, jk
-      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
-      !!----------------------------------------------------------------------
-      !
-      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
-         !                                   ! ===============
-         IF( ln_rstart ) THEN                   !* Read the restart file
-            CALL rst_read_open                  !  open the restart file if necessary
-            CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm)    )
+      INTEGER ::   ji, jj, jk      ! dummy loop indices
+      INTEGER ::   id3, id4, id5   ! local integers
+      !!----------------------------------------------------------------------
+      !
+      !                                      !=====================!
+      IF( TRIM(cdrw) == 'READ' ) THEN        !  Read / initialise  !
+         !                                   !=====================!
+         !
+         IF( ln_rstart ) THEN                   !==  Read the restart file  ==!
             !
-            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
-            id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
-            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
+            CALL rst_read_open                                          !*  open the restart file if necessary
+            !                                         ! --------- !
+            !                                         ! all cases !
+            !                                         ! --------- !
+            !
+            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )  !*  check presence
             id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
-            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
+            id5 = iom_varid( numror, 'hdiv_lf'    , ldstop = .FALSE. )
             !
-            !                             ! --------- !
-            !                             ! all cases !
-            !                             ! --------- !
-            !
-            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
+            !                                                           !*  scale factors
+            IF(lwp) WRITE(numout,*)    '          Kmm scale factor read in the restart file'
+            CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm) )
+            WHERE ( tmask(:,:,:) == 0.0_wp ) 
+               e3t(:,:,:,Kmm) = e3t_0(:,:,:)
+            END WHERE
+            IF( l_1st_euler ) THEN                       ! euler
+               IF(lwp) WRITE(numout,*) '          Euler first time step : e3t(Kbb) = e3t(Kmm)'
+               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
+            ELSE                                         ! leap frog
+               IF(lwp) WRITE(numout,*) '          Kbb scale factor read in the restart file'
                CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb) )
-               CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm) )
-               ! needed to restart if land processor not computed 
-               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files'
                WHERE ( tmask(:,:,:) == 0.0_wp ) 
-                  e3t(:,:,:,Kmm) = e3t_0(:,:,:)
                   e3t(:,:,:,Kbb) = e3t_0(:,:,:)
                END WHERE
-               IF( l_1st_euler ) THEN
-                  e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
-               ENDIF
-            ELSE IF( id1 > 0 ) THEN
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files'
-               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
-               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
-               CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb) )
-               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
-               l_1st_euler = .true.
-            ELSE IF( id2 > 0 ) THEN
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files'
-               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
-               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
-               CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm) )
-               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
-               l_1st_euler = .true.
-            ELSE
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file'
-               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
-               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
-               DO jk = 1, jpk
-                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
-                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
-                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
-               END DO
-               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
-               l_1st_euler = .true.
             ENDIF
-            !                             ! ----------- !
-            IF( ln_vvl_zstar ) THEN       ! z_star case !
-               !                          ! ----------- !
+            !                                         ! ------------ !
+            IF( ln_vvl_zstar ) THEN                   !  z_star case !
+               !                                      ! ------------ !
                IF( MIN( id3, id4 ) > 0 ) THEN
                   CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
                ENDIF
-               !                          ! ----------------------- !
-            ELSE                          ! z_tilde and layer cases !
-               !                          ! ----------------------- !
-               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
-                  CALL iom_get( numror, jpdom_auto, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
+               !                                      ! ------------------------ !
+            ELSE                                      !  z_tilde and layer cases !
+               !                                      ! ------------------------ !
+               !
+               IF( id4 > 0 ) THEN                                       !*  scale factor increments
+                  IF(lwp) WRITE(numout,*)    '          Kmm scale factor increments read in the restart file'
                   CALL iom_get( numror, jpdom_auto, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
-               ELSE                            ! one at least array is missing
+                  IF( l_1st_euler ) THEN                 ! euler
+                     IF(lwp) WRITE(numout,*) '          Euler first time step : tilde_e3t(Kbb) = tilde_e3t(Kmm)'
+                     tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
+                  ELSE                                   ! leap frog
+                     IF(lwp) WRITE(numout,*) '          Kbb scale factor increments read in the restart file'
+                     CALL iom_get( numror, jpdom_auto, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
+                  ENDIF
+               ELSE 
                   tilde_e3t_b(:,:,:) = 0.0_wp
                   tilde_e3t_n(:,:,:) = 0.0_wp
                ENDIF
-               !                          ! ------------ !
-               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
-                  !                       ! ------------ !
+               !                                      ! ------------ !
+               IF( ln_vvl_ztilde ) THEN               ! z_tilde case !
+                  !                                   ! ------------ !
                   IF( id5 > 0 ) THEN  ! required array exists
                      CALL iom_get( numror, jpdom_auto, 'hdiv_lf', hdiv_lf(:,:,:) )
                   ELSE                ! array is missing
-                     hdiv_lf(:,:,:) = 0.0_wp
+                     hdiv_lf(:,:,:) = 0.0_wp 
                   ENDIF
                ENDIF
             ENDIF
             !
-         ELSE                                   !* Initialize at "rest"
+         ELSE                                   !==  Initialize at "rest" with ssh  ==!
             !
-
-            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
-
-                  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 !If test case else
-
-               ! Adjust vertical metrics for all wad
-               DO jk = 1, jpk
-                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm)  ) &
-                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
-                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
-               END DO
-               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
-
-               DO_2D( 1, 1, 1, 1 )
-                  IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
-                     CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
-                  ENDIF
-               END_2D
-               !
-            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)  )
-!               !
-!               DO jk=1,jpk
-!                  e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) &
-!                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk)
-!               END DO
-!               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
-                ssh(:,:,Kmm)=0._wp
-                e3t(:,:,:,Kmm)=e3t_0(:,:,:)
-                e3t(:,:,:,Kbb)=e3t_0(:,:,:)
-               !
-            END IF           ! end of ll_wd edits
-
+            DO jk = 1, jpk
+               e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( 1._wp + ssh(:,:,Kmm) * r1_ht_0(:,:) * tmask(:,:,jk) )
+            END DO
+            e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
+            !
             IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
                tilde_e3t_b(:,:,:) = 0._wp
                tilde_e3t_n(:,:,:) = 0._wp
                IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
-            END IF
+            ENDIF
          ENDIF
-         !
-      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
-         !                                   ! ===================
+         !                                       !=======================!
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN       !  Create restart file  !
+         !                                       !=======================!
+         !
          IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
          !                                           ! --------- !
Index: NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/domzgr_substitute.h90
===================================================================
--- NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/domzgr_substitute.h90	(revision 14046)
+++ NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/domzgr_substitute.h90	(revision 14058)
@@ -15,9 +15,10 @@
 #   define  e3u(i,j,k,t)   (e3u_0(i,j,k)*(1._wp+r3u(i,j,t)*umask(i,j,k)))
 #   define  e3v(i,j,k,t)   (e3v_0(i,j,k)*(1._wp+r3v(i,j,t)*vmask(i,j,k)))
-#   define  e3f(i,j,k)     (e3f_0(i,j,k)*(1._wp+r3f(i,j)*fmask(i,j,k)))
+#   define  e3f(i,j,k)     (e3f_0(i,j,k)*(1._wp+r3f(i,j)*fe3mask(i,j,k)))
+#   define  e3f_vor(i,j,k) (e3f_0vor(i,j,k)*(1._wp+r3f(i,j)*fe3mask(i,j,k)))
 #   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  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 +30,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_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/istate.F90
===================================================================
--- NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/istate.F90	(revision 14046)
+++ NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/istate.F90	(revision 14058)
@@ -42,5 +42,5 @@
    PRIVATE
 
-   PUBLIC   istate_init   ! routine called by step.F90
+   PUBLIC   istate_init   ! routine called by nemogcm.F90
 
    !! * Substitutions
@@ -59,9 +59,11 @@
       !! 
       !! ** Purpose :   Initialization of the dynamics and tracer fields.
+      !!
+      !! ** Method  :   
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in )  ::  Kbb, Kmm, Kaa   ! ocean time level indices
       !
       INTEGER ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zgdept     ! 3D table  !!st patch to use gdept subtitute
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zgdept     ! 3D table for qco substitute
 !!gm see comment further down
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace
@@ -73,15 +75,12 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
 
-!!gm  Why not include in the first call of dta_tsd ?  
-!!gm  probably associated with the use of internal damping...
        CALL dta_tsd_init        ! Initialisation of T & S input data
-!!gm to be moved in usrdef of C1D case
+
 !      IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data
-!!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 +95,11 @@
          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)
+!!st
+!!st need for a recent agrif version to be displaced toward ssh_init_rst with agrif_istate_ssh
+         ssh(:,:,    Kmm) = ssh(:,:    ,Kbb)
+!!st end
+         uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb)
+         vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb)
       ELSE
 #endif
@@ -117,55 +119,27 @@
             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  
+            uu (:,:,:,Kbb) = 0._wp
+            vv (:,:,:,Kbb) = 0._wp  
             !
-            IF( ll_wd ) THEN
-               ssh(:,:,Kbb) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD 
-               !
-               ! Apply minimum wetdepth criterion
-               !
-               DO_2D( 1, 1, 1, 1 )
-                  IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN
-                     ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )
-                  ENDIF
-               END_2D
-            ENDIF 
-             !
          ELSE                                 ! user defined initial T and S
             DO jk = 1, jpk
                zgdept(:,:,jk) = gdept(:,:,jk,Kbb)
             END DO
-            CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )         
+            CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,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)
-
-!!gm POTENTIAL BUG :
-!!gm  ISSUE :  if ssh(:,:,Kbb) /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
-!!             as well as gdept_ and gdepw_....   !!!!! 
-!!      ===>>>>   probably a call to domvvl initialisation here....
-
+         ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
+         uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb)
+         vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb)
 
          !
-!!gm to be moved in usrdef of C1D case
-!         IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
-!            ALLOCATE( zuvd(jpi,jpj,jpk,2) )
-!            CALL dta_uvd( nit000, zuvd )
-!            uu(:,:,:,Kbb) = zuvd(:,:,:,1) ;  uu(:,:,:,Kmm) = uu(:,:,:,Kbb)
-!            vv(:,:,:,Kbb) = zuvd(:,:,:,2) ;  vv(:,:,:,Kmm) = vv(:,:,:,Kbb)
-!            DEALLOCATE( zuvd )
-!         ENDIF
+!!gm ==>>>  to be moved in usrdef_istate of C1D case 
+         IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
+            ALLOCATE( zuvd(jpi,jpj,jpk,2) )
+            CALL dta_uvd( nit000, Kbb, zuvd )
+            uu(:,:,:,Kbb) = zuvd(:,:,:,1)  ;  uu(:,:,:,Kmm) = uu(:,:,:,Kbb)
+            vv(:,:,:,Kbb) = zuvd(:,:,:,2)  ;  vv(:,:,:,Kmm) = vv(:,:,:,Kbb)
+            DEALLOCATE( zuvd )
+         ENDIF
          !
-!!gm This is to be changed !!!!
-!         ! - ML - ssh(:,:,Kmm) could be modified by istate_eel, so that initialization of e3t(:,:,:,Kbb) is done here
-!         IF( .NOT.ln_linssh ) THEN
-!            DO jk = 1, jpk
-!               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm)
-!            END DO
-!         ENDIF
-!!gm 
          ! 
       ENDIF 
Index: NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/phycst.F90
===================================================================
--- NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/phycst.F90	(revision 14046)
+++ NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/phycst.F90	(revision 14058)
@@ -66,4 +66,5 @@
    REAL(wp), PUBLIC ::   r1_rhos                     !: 1 / rhos
    REAL(wp), PUBLIC ::   r1_rcpi                     !: 1 / rcpi
+   
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
