Index: NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/DOM/domain.F90
===================================================================
--- NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/DOM/domain.F90	(revision 14037)
+++ NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/DOM/domain.F90	(revision 14062)
@@ -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
