Index: branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 7349)
+++ branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 7351)
@@ -20,9 +20,12 @@
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
+   USE phycst          ! physical constant
    USE dom_oce         ! ocean space and time domain
    USE sbc_oce         ! ocean surface boundary condition
+   USE wet_dry         ! wetting and drying
+   USE restart         ! ocean restart
+   !
    USE in_out_manager  ! I/O manager
    USE iom             ! I/O manager library
-   USE restart         ! ocean restart
    USE lib_mpp         ! distributed memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -60,12 +63,10 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO-Consortium (2010) 
+   !! NEMO/OPA 3.7 , NEMO-Consortium (2015) 
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -74,5 +75,5 @@
       !!                ***  FUNCTION dom_vvl_alloc  ***
       !!----------------------------------------------------------------------
-      IF( ln_vvl_zstar ) dom_vvl_alloc = 0
+      IF( ln_vvl_zstar )   dom_vvl_alloc = 0
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
          ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
@@ -81,6 +82,6 @@
          IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
          IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
-         un_td = 0.0_wp
-         vn_td = 0.0_wp
+         un_td = 0._wp
+         vn_td = 0._wp
       ENDIF
       IF( ln_vvl_ztilde ) THEN
@@ -89,5 +90,5 @@
          IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
       ENDIF
-
+      !
    END FUNCTION dom_vvl_alloc
 
@@ -103,11 +104,11 @@
       !!               - interpolate scale factors
       !!
-      !! ** Action  : - fse3t_(n/b) and tilde_e3t_(n/b)
-      !!              - Regrid: fse3(u/v)_n
-      !!                        fse3(u/v)_b       
-      !!                        fse3w_n           
-      !!                        fse3(u/v)w_b      
-      !!                        fse3(u/v)w_n      
-      !!                        fsdept_n, fsdepw_n and fsde3w_n
+      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b)
+      !!              - Regrid: e3(u/v)_n
+      !!                        e3(u/v)_b       
+      !!                        e3w_n           
+      !!                        e3(u/v)w_b      
+      !!                        e3(u/v)w_n      
+      !!                        gdept_n, gdepw_n and gde3w_n
       !!              - h(t/u/v)_0
       !!              - frq_rst_e3t and frq_rst_hdv
@@ -115,108 +116,110 @@
       !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
       !!----------------------------------------------------------------------
-      USE phycst,  ONLY : rpi, rsmall, rad
-      !! * Local declarations
-      INTEGER ::   ji,jj,jk
+      INTEGER ::   ji, jj, jk
       INTEGER ::   ii0, ii1, ij0, ij1
       REAL(wp)::   zcoef
       !!----------------------------------------------------------------------
-      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_init')
-
+      !
+      IF( nn_timing == 1 )   CALL timing_start('dom_vvl_init')
+      !
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
-
-      ! choose vertical coordinate (z_star, z_tilde or layer)
-      ! ==========================
-      CALL dom_vvl_ctl
-
-      ! Allocate module arrays
-      ! ======================
+      !
+      CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer)
+      !
+      !                    ! Allocate module arrays
       IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
-
-      ! Read or initialize fse3t_(b/n), tilde_e3t_(b/n) and hdiv_lf (and e3t_a(jpk))
-      ! ============================================================================
+      !
+      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
       CALL dom_vvl_rst( nit000, 'READ' )
-      fse3t_a(:,:,jpk) = e3t_0(:,:,jpk)
-
-      ! Reconstruction of all vertical scale factors at now and before time steps
-      ! =============================================================================
-      ! Horizontal scale factor interpolations
-      ! --------------------------------------
-      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
-      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
-      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
-      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
-      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
-      ! Vertical scale factor interpolations
-      ! ------------------------------------
-      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
-      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
-      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W'  )
-      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
-      ! t- and w- points depth
-      ! ----------------------
-      ! set the isf depth as it is in the initial step
-      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
-      fsdepw_n(:,:,1) = 0.0_wp
-      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
-      fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1)
-      fsdepw_b(:,:,1) = 0.0_wp
-
-      DO jk = 2, jpk
+      e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
+      !
+      !                    !== Set of all other vertical scale factors  ==!  (now and before)
+      !                                ! Horizontal interpolation of e3t
+      CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )    ! from T to U
+      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
+      CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )    ! from T to V 
+      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
+      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )    ! from U to F
+      !                                ! Vertical interpolation of e3t,u,v 
+      CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )  ! from T to W
+      CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W'  )
+      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )  ! from U to UW
+      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
+      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )  ! from V to UW
+      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
+      !
+      !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
+      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)       ! reference to the ocean surface (used for MLD and light penetration)
+      gdepw_n(:,:,1) = 0.0_wp
+      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)  ! reference to a common level z=0 for hpg
+      gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1)
+      gdepw_b(:,:,1) = 0.0_wp
+      DO jk = 2, jpk                               ! vertical sum
          DO jj = 1,jpj
             DO ji = 1,jpi
-              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
-                                                     ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
-                                                     ! 0.5 where jk = mikt  
-               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
-               fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  &
-                   &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) +       fse3w_n(ji,jj,jk)) 
-               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj)
-               fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1)
-               fsdept_b(ji,jj,jk) =      zcoef  * ( fsdepw_b(ji,jj,jk  ) + 0.5 * fse3w_b(ji,jj,jk))  &
-                   &                + (1-zcoef) * ( fsdept_b(ji,jj,jk-1) +       fse3w_b(ji,jj,jk)) 
+               !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
+               !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
+               !                             ! 0.5 where jk = mikt     
+!!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ??
+               zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
+               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
+               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
+                  &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
+               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
+               gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
+               gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
+                  &                + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
             END DO
          END DO
       END DO
-
-      ! Before depth and Inverse of the local depth of the water column at u- and v- points
-      ! -----------------------------------------------------------------------------------
-      hu_b(:,:) = 0.
-      hv_b(:,:) = 0.
-      DO jk = 1, jpkm1
-         hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
-         hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
+      !
+      !                    !==  thickness of the water column  !!   (ocean portion only)
+      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
+      hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
+      hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1)
+      hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
+      hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1)
+      DO jk = 2, jpkm1
+         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
+         hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
+         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
+         hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
+         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
       END DO
-      hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) )
-      hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) )
-
-      ! Restoring frequencies for z_tilde coordinate
-      ! ============================================
+      !
+      !                    !==  inverse of water column thickness   ==!   (u- and v- points)
+      r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
+      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
+      r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
+      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
+
+      !                    !==   z_tilde coordinate case  ==!   (Restoring frequencies)
       IF( ln_vvl_ztilde ) THEN
-         ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings
-         frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
-         frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
-         IF( ln_vvl_ztilde_as_zstar ) THEN
-            ! Ignore namelist settings and use these next two to emulate z-star using z-tilde
-            frq_rst_e3t(:,:) = 0.0_wp 
-            frq_rst_hdv(:,:) = 1.0_wp / rdt
+!!gm : idea: add here a READ in a file of custumized restoring frequency
+         !                                   ! Values in days provided via the namelist
+         !                                   ! use rsmall to avoid possible division by zero errors with faulty settings
+         frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
+         frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
+         !
+         IF( ln_vvl_ztilde_as_zstar ) THEN   ! z-star emulation using z-tile
+            frq_rst_e3t(:,:) = 0._wp               !Ignore namelist settings
+            frq_rst_hdv(:,:) = 1._wp / rdt
          ENDIF
-         IF ( ln_vvl_zstar_at_eqtor ) THEN
+         IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator
             DO jj = 1, jpj
                DO ji = 1, jpi
+!!gm  case |gphi| >= 6 degrees is useless   initialized just above by default
                   IF( ABS(gphit(ji,jj)) >= 6.) THEN
                      ! values outside the equatorial band and transition zone (ztilde)
                      frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp )
                      frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
-                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
+                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star
                      ! values inside the equatorial band (ztilde as zstar)
                      frq_rst_e3t(ji,jj) =  0.0_wp
                      frq_rst_hdv(ji,jj) =  1.0_wp / rdt
-                  ELSE
-                     ! values in the transition band (linearly vary from ztilde to ztilde as zstar values)
+                  ELSE                                      ! transition band (2.5 to 6 degrees N/S)
+                     !                                      ! (linearly transition from z-tilde to z-star)
                      frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   &
                         &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
@@ -229,6 +232,6 @@
                END DO
             END DO
-            IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN
-               ii0 = 103   ;   ii1 = 111        ! Suppress ztilde in the Foxe Basin for ORCA2
+            IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN   ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2
+               ii0 = 103   ;   ii1 = 111       
                ij0 = 128   ;   ij1 = 135   ;   
                frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0.0_wp
@@ -237,7 +240,7 @@
          ENDIF
       ENDIF
-
+      !
       IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_init')
-
+      !
    END SUBROUTINE dom_vvl_init
 
@@ -261,25 +264,27 @@
       !!               - tilde_e3t_a: after increment of vertical scale factor 
       !!                              in z_tilde case
-      !!               - fse3(t/u/v)_a
+      !!               - e3(t/u/v)_a
       !!
       !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
       !!----------------------------------------------------------------------
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zht, z_scale, zwu, zwv, zhdiv
-      !! * Arguments
-      INTEGER, INTENT( in )                  :: kt                    ! time step
-      INTEGER, INTENT( in ), OPTIONAL        :: kcall                 ! optional argument indicating call sequence
-      !! * Local declarations
-      INTEGER                                :: ji, jj, jk            ! dummy loop indices
-      INTEGER , DIMENSION(3)                 :: ijk_max, ijk_min      ! temporary integers
-      REAL(wp)                               :: z2dt                  ! temporary scalars
-      REAL(wp)                               :: z_tmin, z_tmax        ! temporary scalars
-      LOGICAL                                :: ll_do_bclinic         ! temporary logical
-      !!----------------------------------------------------------------------
-      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt')
-      CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
-      CALL wrk_alloc( jpi, jpj, jpk, ze3t                     )
-
-      IF(kt == nit000)   THEN
+      INTEGER, INTENT( in )           ::   kt      ! time step
+      INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence
+      !
+      INTEGER                ::   ji, jj, jk            ! dummy loop indices
+      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
+      REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars
+      LOGICAL                ::   ll_do_bclinic         ! local logical
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze3t
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zht, z_scale, zwu, zwv, zhdiv
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
+      !
+      IF( nn_timing == 1 )   CALL timing_start('dom_vvl_sf_nxt')
+      !
+      CALL wrk_alloc( jpi,jpj,zht,   z_scale, zwu, zwv, zhdiv )
+      CALL wrk_alloc( jpi,jpj,jpk,   ze3t )
+
+      IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
@@ -289,5 +294,5 @@
       ll_do_bclinic = .TRUE.
       IF( PRESENT(kcall) ) THEN
-         IF ( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE.
+         IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
       ENDIF
 
@@ -295,18 +300,16 @@
       ! After acale factors at t-points !
       ! ******************************* !
-
       !                                                ! --------------------------------------------- !
-                                                       ! z_star coordinate and barotropic z-tilde part !
+      !                                                ! z_star coordinate and barotropic z-tilde part !
       !                                                ! --------------------------------------------- !
-
+      !
       z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
       DO jk = 1, jpkm1
-         ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)
-         fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+         ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0)
+         e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
       END DO
-
+      !
       IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
          !                                                            ! ------baroclinic part------ !
-
          ! I - initialization
          ! ==================
@@ -314,9 +317,9 @@
          ! 1 - barotropic divergence
          ! -------------------------
-         zhdiv(:,:) = 0.
-         zht(:,:)   = 0.
+         zhdiv(:,:) = 0._wp
+         zht(:,:)   = 0._wp
          DO jk = 1, jpkm1
-            zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
-            zht  (:,:) = zht  (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+            zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
+            zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
          END DO
          zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
@@ -325,15 +328,15 @@
          ! --------------------------------------------------
          IF( ln_vvl_ztilde ) THEN
-            IF( kt .GT. nit000 ) THEN
+            IF( kt > nit000 ) THEN
                DO jk = 1, jpkm1
                   hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
-                     &          * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
+                     &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
                END DO
             ENDIF
-         END IF
+         ENDIF
 
          ! II - after z_tilde increments of vertical scale factors
          ! =======================================================
-         tilde_e3t_a(:,:,:) = 0.0_wp  ! tilde_e3t_a used to store tendency terms
+         tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms
 
          ! 1 - High frequency divergence term
@@ -341,11 +344,11 @@
          IF( ln_vvl_ztilde ) THEN     ! z_tilde case
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
             END DO
          ELSE                         ! layer case
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
-            END DO
-         END IF
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
+            END DO
+         ENDIF
 
          ! 2 - Restoring term (z-tilde case only)
@@ -355,12 +358,11 @@
                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
             END DO
-         END IF
+         ENDIF
 
          ! 3 - Thickness diffusion term
          ! ----------------------------
-         zwu(:,:) = 0.0_wp
-         zwv(:,:) = 0.0_wp
-         ! a - first derivative: diffusive fluxes
-         DO jk = 1, jpkm1
+         zwu(:,:) = 0._wp
+         zwv(:,:) = 0._wp
+         DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
@@ -374,6 +376,5 @@
             END DO
          END DO
-         ! b - correction for last oceanic u-v points
-         DO jj = 1, jpj
+         DO jj = 1, jpj          ! b - correction for last oceanic u-v points
             DO ji = 1, jpi
                un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
@@ -381,6 +382,5 @@
             END DO
          END DO
-         ! c - second derivative: divergence of diffusive fluxes
-         DO jk = 1, jpkm1
+         DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -391,6 +391,6 @@
             END DO
          END DO
-         ! d - thickness diffusion transport: boundary conditions
-         !     (stored for tracer advction and continuity equation)
+         !                       ! d - thickness diffusion transport: boundary conditions
+         !                             (stored for tracer advction and continuity equation)
          CALL lbc_lnk( un_td , 'U' , -1._wp)
          CALL lbc_lnk( vn_td , 'V' , -1._wp)
@@ -410,5 +410,5 @@
          ! Maximum deformation control
          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         ze3t(:,:,jpk) = 0.0_wp
+         ze3t(:,:,jpk) = 0._wp
          DO jk = 1, jpkm1
             ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
@@ -419,5 +419,5 @@
          IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
          ! - ML - test: for the moment, stop simulation for too large e3_t variations
-         IF( ( z_tmax .GT. rn_zdef_max ) .OR. ( z_tmin .LT. - rn_zdef_max ) ) THEN
+         IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN
             IF( lk_mpp ) THEN
                CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) )
@@ -462,5 +462,5 @@
          z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
          DO jk = 1, jpkm1
-            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
          END DO
 
@@ -470,5 +470,5 @@
       !                                           ! ---baroclinic part--------- !
          DO jk = 1, jpkm1
-            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
+            e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
       ENDIF
@@ -485,25 +485,25 @@
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+            zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
          END DO
          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(fse3t_n))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax
          !
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)
+            zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(fse3t_a))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax
          !
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk)
+            zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk)
          END DO
          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(fse3t_b))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax
          !
          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) )
@@ -524,6 +524,6 @@
       ! *********************************** !
 
-      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' )
-      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' )
+      CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' )
+      CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' )
 
       ! *********************************** !
@@ -531,19 +531,20 @@
       ! *********************************** !
 
-      hu_a(:,:) = 0._wp                        ! Ocean depth at U-points
-      hv_a(:,:) = 0._wp                        ! Ocean depth at V-points
-      DO jk = 1, jpkm1
-         hu_a(:,:) = hu_a(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)
-         hv_a(:,:) = hv_a(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk)
+      hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1)
+      hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)
+      DO jk = 2, jpkm1
+         hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)
+         hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)
       END DO
       !                                        ! Inverse of the local depth
-      hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
-      hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
-
-      CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
-      CALL wrk_dealloc( jpi, jpj, jpk, ze3t                     )
-
+!!gm BUG ?  don't understand the use of umask_i here .....
+      r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
+      r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )
+      !
+      CALL wrk_dealloc( jpi,jpj,       zht, z_scale, zwu, zwv, zhdiv )
+      CALL wrk_dealloc( jpi,jpj,jpk,   ze3t )
+      !
       IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_nxt')
-
+      !
    END SUBROUTINE dom_vvl_sf_nxt
 
@@ -561,11 +562,11 @@
       !!               - recompute depths and water height fields
       !!
-      !! ** Action  :  - fse3t_(b/n), tilde_e3t_(b/n) and fse3(u/v)_n ready for next time step
+      !! ** Action  :  - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step
       !!               - Recompute:
-      !!                    fse3(u/v)_b       
-      !!                    fse3w_n           
-      !!                    fse3(u/v)w_b      
-      !!                    fse3(u/v)w_n      
-      !!                    fsdept_n, fsdepw_n  and fsde3w_n
+      !!                    e3(u/v)_b       
+      !!                    e3w_n           
+      !!                    e3(u/v)w_b      
+      !!                    e3(u/v)w_n      
+      !!                    gdept_n, gdepw_n  and gde3w_n
       !!                    h(u/v) and h(u/v)r
       !!
@@ -573,11 +574,12 @@
       !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
       !!----------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in )               :: kt       ! time step
-      !! * Local declarations
-      INTEGER                             :: ji,jj,jk       ! dummy loop indices
-      REAL(wp)                            :: zcoef
-      !!----------------------------------------------------------------------
-
+      INTEGER, INTENT( in ) ::   kt   ! time step
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zcoef        ! local scalar
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
+      !
       IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_swp')
       !
@@ -590,5 +592,5 @@
       ! Time filter and swap of scale factors
       ! =====================================
-      ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.
+      ! - ML - e3(t/u/v)_b are allready computed in dynnxt.
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
          IF( neuler == 0 .AND. kt == nit000 ) THEN
@@ -600,10 +602,10 @@
          tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
       ENDIF
-      fsdept_b(:,:,:) = fsdept_n(:,:,:)
-      fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
-
-      fse3t_n(:,:,:) = fse3t_a(:,:,:)
-      fse3u_n(:,:,:) = fse3u_a(:,:,:)
-      fse3v_n(:,:,:) = fse3v_a(:,:,:)
+      gdept_b(:,:,:) = gdept_n(:,:,:)
+      gdepw_b(:,:,:) = gdepw_n(:,:,:)
+
+      e3t_n(:,:,:) = e3t_a(:,:,:)
+      e3u_n(:,:,:) = e3u_a(:,:,:)
+      e3v_n(:,:,:) = e3v_a(:,:,:)
 
       ! Compute all missing vertical scale factor and depths
@@ -611,22 +613,21 @@
       ! Horizontal scale factor interpolations
       ! --------------------------------------
-      ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt
+      ! - ML - e3u_b and e3v_b are allready computed in dynnxt
       ! - JC - hu_b, hv_b, hur_b, hvr_b also
-      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F'  )
+      
+      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  )
+      
       ! Vertical scale factor interpolations
-      ! ------------------------------------
-      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
-      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
-      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W'  )
-      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
-      ! t- and w- points depth
-      ! ----------------------
-      ! set the isf depth as it is in the initial step
-      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
-      fsdepw_n(:,:,1) = 0.0_wp
-      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
-
+      CALL dom_vvl_interpol( e3t_n(:,:,:),  e3w_n(:,:,:), 'W'  )
+      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
+      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
+      CALL dom_vvl_interpol( e3t_b(:,:,:),  e3w_b(:,:,:), 'W'  )
+      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
+      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
+
+      ! t- and w- points depth (set the isf depth as it is in the initial step)
+      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
+      gdepw_n(:,:,1) = 0.0_wp
+      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
       DO jk = 2, jpk
          DO jj = 1,jpj
@@ -635,44 +636,28 @@
                                                                  ! 1 for jk = mikt
                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
-               fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  &
-                   &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) +       fse3w_n(ji,jj,jk)) 
-               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj)
+               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
+               gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  &
+                   &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) ) 
+               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
             END DO
          END DO
       END DO
 
-      ! Local depth and Inverse of the local depth of the water column at u- and v- points
-      ! ----------------------------------------------------------------------------------
-      hu (:,:) = hu_a (:,:)
-      hv (:,:) = hv_a (:,:)
-
-      ! Inverse of the local depth
-      hur(:,:) = hur_a(:,:)
-      hvr(:,:) = hvr_a(:,:)
-
-      ! Local depth of the water column at t- points
-      ! --------------------------------------------
-      ht(:,:) = 0.
-      DO jk = 1, jpkm1
-         ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+      ! Local depth and Inverse of the local depth of the water
+      ! -------------------------------------------------------
+      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
+      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
+      !
+      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)
+      DO jk = 2, jpkm1
+         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
       END DO
-
-      ! Write outputs
-      ! =============
-      CALL iom_put(     "e3t" , fse3t_n  (:,:,:) )
-      CALL iom_put(     "e3u" , fse3u_n  (:,:,:) )
-      CALL iom_put(     "e3v" , fse3v_n  (:,:,:) )
-      CALL iom_put(     "e3w" , fse3w_n  (:,:,:) )
-      CALL iom_put( "tpt_dep" , fsde3w_n (:,:,:) )
-      IF( iom_use("e3tdef") )   &
-         CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
 
       ! write restart file
       ! ==================
-      IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_swp')
-
+      IF( lrst_oce )   CALL dom_vvl_rst( kt, 'WRITE' )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('dom_vvl_sf_swp')
+      !
    END SUBROUTINE dom_vvl_sf_swp
 
@@ -693,8 +678,15 @@
       !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
       !
-      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_interpol')
+      INTEGER ::   ji, jj, jk                                       ! dummy loop indices
+      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd = T/F
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('dom_vvl_interpol')
+      !
+      IF(ln_wd) THEN
+        zlnwd = 1.0_wp
+      ELSE
+        zlnwd = 0.0_wp
+      END IF
       !
       SELECT CASE ( pout )    !==  type of interpolation  ==!
@@ -704,5 +696,5 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e1e2u(ji,jj)                                   &
+                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
                      &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
                      &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
@@ -717,5 +709,5 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e1e2v(ji,jj)                                   &
+                  pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
                      &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
                      &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
@@ -730,5 +722,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e1e2f(ji,jj)               &
+                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
+                     &                       *    r1_e1e2f(ji,jj)                                                  &
                      &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
                      &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
@@ -745,6 +738,8 @@
 !!gm BUG? use here wmask in case of ISF ?  to be checked
          DO jk = 2, jpk
-            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   &
-               &                            +            0.5_wp * tmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
+            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   &
+               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               &
+               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     &
+               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
          END DO
          !
@@ -755,6 +750,8 @@
 !!gm BUG? use here wumask in case of ISF ?  to be checked
          DO jk = 2, jpk
-            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   &
-               &                             +            0.5_wp * umask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
+            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
+               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              &
+               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
+               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
          END DO
          !
@@ -765,10 +762,12 @@
 !!gm BUG? use here wvmask in case of ISF ?  to be checked
          DO jk = 2, jpk
-            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   &
-               &                             +            0.5_wp * vmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
+            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
+               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              &
+               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
+               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
          END DO
       END SELECT
       !
-      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_interpol')
+      IF( nn_timing == 1 )   CALL timing_stop('dom_vvl_interpol')
       !
    END SUBROUTINE dom_vvl_interpol
@@ -790,5 +789,5 @@
       CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
       !
-      INTEGER ::   jk
+      INTEGER ::   ji, jj, jk
       INTEGER ::   id1, id2, id3, id4, id5     ! local integers
       !!----------------------------------------------------------------------
@@ -801,6 +800,6 @@
             CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn    )
             !
-            id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. )
-            id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. )
+            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. )
             id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
@@ -810,39 +809,39 @@
             !                             ! --------- !
             IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
-               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
-               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) )
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) )
                ! needed to restart if land processor not computed 
-               IF(lwp) write(numout,*) 'dom_vvl_rst : fse3t_b and fse3t_n found in restart files'
+               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files'
                WHERE ( tmask(:,:,:) == 0.0_wp ) 
-                  fse3t_n(:,:,:) = e3t_0(:,:,:)
-                  fse3t_b(:,:,:) = e3t_0(:,:,:)
+                  e3t_n(:,:,:) = e3t_0(:,:,:)
+                  e3t_b(:,:,:) = e3t_0(:,:,:)
                END WHERE
                IF( neuler == 0 ) THEN
-                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
+                  e3t_b(:,:,:) = e3t_n(:,:,:)
                ENDIF
             ELSE IF( id1 > 0 ) THEN
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files'
-               IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.'
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'
+               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
-               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
-               fse3t_n(:,:,:) = fse3t_b(:,:,:)
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) )
+               e3t_n(:,:,:) = e3t_b(:,:,:)
                neuler = 0
             ELSE IF( id2 > 0 ) THEN
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files'
-               IF(lwp) write(numout,*) 'fse3t_b set equal to fse3t_n.'
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files'
+               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
-               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
-               fse3t_b(:,:,:) = fse3t_n(:,:,:)
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) )
+               e3t_b(:,:,:) = e3t_n(:,:,:)
                neuler = 0
             ELSE
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart file'
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file'
                IF(lwp) write(numout,*) 'Compute scale factor from sshn'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
-               DO jk=1,jpk
-                  fse3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
-                      &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
-                      &            + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
+               DO jk = 1, jpk
+                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
+                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
+                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
                END DO
-               fse3t_b(:,:,:) = fse3t_n(:,:,:)
+               e3t_b(:,:,:) = e3t_n(:,:,:)
                neuler = 0
             ENDIF
@@ -875,7 +874,23 @@
             !
          ELSE                                   !* Initialize at "rest"
-            fse3t_b(:,:,:) = e3t_0(:,:,:)
-            fse3t_n(:,:,:) = e3t_0(:,:,:)
+            e3t_b(:,:,:) = e3t_0(:,:,:)
+            e3t_n(:,:,:) = e3t_0(:,:,:)
             sshn(:,:) = 0.0_wp
+
+            IF( ln_wd ) THEN
+              DO jj = 1, jpj
+                DO ji = 1, jpi
+                  IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1 ) THEN
+                     e3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1 
+                     e3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1 
+                     e3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1 
+                     sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj)
+                     sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj)
+                     ssha(ji,jj) = rn_wdmin1 - bathy(ji,jj)
+                  ENDIF
+                ENDDO
+              ENDDO
+            END IF
+
             IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
                tilde_e3t_b(:,:,:) = 0.0_wp
@@ -891,6 +906,6 @@
          !                                           ! all cases !
          !                                           ! --------- !
-         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
-         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )
+         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:) )
+         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:) )
          !                                           ! ----------------------- !
          IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
@@ -975,5 +990,5 @@
       !
       IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
-      IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
+      IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
       !
       IF(lwp) THEN                   ! Print the choice
@@ -989,5 +1004,5 @@
       !
 #if defined key_agrif
-      IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' )
+      IF(.NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented with non-linear free surface' )
 #endif
       !
@@ -996,5 +1011,2 @@
    !!======================================================================
 END MODULE domvvl
-
-
-
