Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DIA/diawri.F90
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DIA/diawri.F90	(revision 10115)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DIA/diawri.F90	(revision 10116)
@@ -125,6 +125,6 @@
       ! Output of initial vertical scale factor
       CALL iom_put("e3t_0", e3t_0(:,:,:) )
-      CALL iom_put("e3u_0", e3t_0(:,:,:) )
-      CALL iom_put("e3v_0", e3t_0(:,:,:) )
+      CALL iom_put("e3u_0", e3u_0(:,:,:) )
+      CALL iom_put("e3v_0", e3v_0(:,:,:) )
       !
       CALL iom_put( "e3t" , e3t_n(:,:,:) )
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DOM/domvvl.F90
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DOM/domvvl.F90	(revision 10115)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DOM/domvvl.F90	(revision 10116)
@@ -8,4 +8,5 @@
    !!            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.0  !  2018-09  (J. Chanut) improve z_tilde robustness
    !!----------------------------------------------------------------------
 
@@ -31,4 +32,7 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE timing          ! Timing
+   USE bdy_oce         ! ocean open boundary conditions
+   USE sbcrnf          ! river runoff 
+   USE dynspg_ts, ONLY: un_adv, vn_adv
 
    IMPLICIT NONE
@@ -44,20 +48,33 @@
    LOGICAL , PUBLIC :: ln_vvl_ztilde          = .FALSE.    ! ztilde vertical coordinate
    LOGICAL , PUBLIC :: ln_vvl_layer           = .FALSE.    ! level  vertical coordinate
-   LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate
-   LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate
-   LOGICAL , PUBLIC :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer
+   LOGICAL          :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate
+   LOGICAL          :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate
+   LOGICAL          :: ln_vvl_zstar_on_shelf  = .FALSE.    ! revert to zstar on shelves
+   LOGICAL          :: ln_vvl_adv_fct         = .FALSE.    ! Centred thickness advection
+   LOGICAL          :: ln_vvl_adv_cn2         = .TRUE.     ! FCT thickness advection
+   LOGICAL          :: ln_vvl_dbg             = .FALSE.    ! debug control prints
+   LOGICAL          :: ln_vvl_ramp            = .FALSE.    ! Ramp on interfaces displacement
+   LOGICAL          :: ln_vvl_lap             = .FALSE.    ! Laplacian thickness diffusion
+   LOGICAL          :: ln_vvl_blp             = .FALSE.    ! Bilaplacian thickness diffusion 
+   LOGICAL          :: ln_vvl_regrid          = .FALSE.    ! ensure layer separation 
+   LOGICAL          :: ll_shorizd             = .FALSE.    ! Use "shelf horizon depths" 
+   LOGICAL          :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer
    !                                                       ! conservation: not used yet
-   REAL(wp)         :: rn_ahe3                             ! thickness diffusion coefficient
-   REAL(wp)         :: rn_rst_e3t                          ! ztilde to zstar restoration timescale [days]
-   REAL(wp)         :: rn_lf_cutoff                        ! cutoff frequency for low-pass filter  [days]
-   REAL(wp)         :: rn_zdef_max                         ! maximum fractional e3t deformation
-   LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE.                ! debug control prints
+   REAL(wp)         :: rn_ahe3_lap               ! thickness diffusion coefficient (Laplacian)
+   REAL(wp)         :: rn_ahe3_blp               ! thickness diffusion coefficient (Bilaplacian)
+   REAL(wp)         :: rn_rst_e3t                ! ztilde to zstar restoration timescale [days]
+   REAL(wp)         :: rn_lf_cutoff              ! cutoff frequency for low-pass filter  [days]
+   REAL(wp)         :: rn_day_ramp               ! Duration of linear ramp  [days]
+   REAL(wp)         :: hsmall=0.01_wp            ! small thickness [m]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                ! thickness diffusion transport
-   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                     ! low frequency part of hz divergence
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_lf, vn_lf, hdivn_lf      ! low frequency fluxes and divergence
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n    ! baroclinic scale factors
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a   ! baroclinic scale factors
-   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                 ! retoring period for scale factors
-   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: tildemask                   ! mask tilde tendency
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                 ! restoring period for scale factors
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! restoring period for low freq. divergence
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: hsm, dsm                    ! 
+   INTEGER         , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: i_int_bot
 
    !! * Substitutions
@@ -78,5 +95,5 @@
          ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
             &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   &
-            &      STAT = dom_vvl_alloc        )
+            &      tildemask(jpi,jpj) , hsm(jpi,jpj) , dsm(jpi,jpj) , i_int_bot(jpi,jpj), STAT = dom_vvl_alloc )
          IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
          IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
@@ -85,5 +102,6 @@
       ENDIF
       IF( ln_vvl_ztilde ) THEN
-         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
+         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdivn_lf(jpi,jpj,jpk)             ,    &
+            &      un_lf(jpi,jpj,jpk), vn_lf(jpi,jpj,jpk)      , STAT= dom_vvl_alloc )
          IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
          IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
@@ -117,5 +135,5 @@
       INTEGER ::   ji, jj, jk
       INTEGER ::   ii0, ii1, ij0, ij1
-      REAL(wp)::   zcoef
+      REAL(wp)::   zcoef, zwgt, ztmp, zhmin, zhmax
       !!----------------------------------------------------------------------
       !
@@ -129,5 +147,5 @@
       IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
       !
-      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
+      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdivn_lf
       CALL dom_vvl_rst( nit000, 'READ' )
       e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
@@ -199,46 +217,71 @@
 
       !                    !==   z_tilde coordinate case  ==!   (Restoring frequencies)
+      tildemask(:,:) = 1._wp
+
       IF( ln_vvl_ztilde ) THEN
 !!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
+         ! 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
+            tildemask(:,:) = 0._wp
          ENDIF
-         IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator
+        
+         IF ( ln_vvl_zstar_at_eqtor ) THEN
             DO jj = 1, jpj
                DO ji = 1, jpi
-!!gm  case |gphi| >= 6 degrees is useless   initialized just above by default
+!!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    ! Equator strip ==> z-star
+!                     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
                      ! 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                                      ! transition band (2.5 to 6 degrees N/S)
-                     !                                      ! (linearly transition from z-tilde to z-star)
+!                     frq_rst_hdv(ji,jj) =  1.0_wp / rdt
+                     tildemask(ji,jj) = 0._wp
+                  ELSE
+                     ! values in the transition band (linearly vary from ztilde to ztilde as zstar values)
                      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)  &
                         &                                          * 180._wp / 3.5_wp ) )
-                     frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                &
-                        &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   &
-                        &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
-                        &                                          * 180._wp / 3.5_wp ) )
+!                     frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                &
+!                        &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   &
+!                        &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
+!                        &                                          * 180._wp / 3.5_wp ) )
+                     tildemask(ji,jj) = 0.5_wp * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
+                        &                                                 * 180._wp / 3.5_wp ) )
                   ENDIF
                END DO
             END DO
-            IF( cn_cfg == "orca" .AND. nn_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
-               frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rdt
-            ENDIF
          ENDIF
+         !
+         IF ( ln_vvl_zstar_on_shelf ) THEN
+            zhmin =  50._wp
+            zhmax = 100._wp
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zwgt = 1._wp
+                  IF(( ht_0(ji,jj)>zhmin).AND.(ht_0(ji,jj) <=zhmax)) THEN
+                     zwgt = (ht_0(ji,jj)-zhmin)/(zhmax-zhmin)
+                  ELSEIF ( ht_0(ji,jj)<=zhmin) THEN
+                     zwgt = 0._wp
+                  ENDIF
+                  frq_rst_e3t(ji,jj) = MIN(frq_rst_e3t(ji,jj), frq_rst_e3t(ji,jj)*zwgt)
+                  tildemask(ji,jj)   = MIN(tildemask(ji,jj), zwgt)
+               END DO
+            END DO
+         ENDIF
+         !
+         ztmp = MAXVAL( frq_rst_hdv(:,:) )
+         IF( lk_mpp )   CALL mpp_max( ztmp )                 ! max over the global domain
+         !
+         IF ( (ztmp*rdt) > 1._wp) CALL ctl_stop( 'dom_vvl_init: rn_lf_cuttoff is too small' )
+         !
       ENDIF
       !
@@ -256,5 +299,7 @@
          IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
             !                                        ! ------------ !
-            CALL iom_set_rstw_var_active('hdiv_lf')
+            CALL iom_set_rstw_var_active('un_lf')
+            CALL iom_set_rstw_var_active('vn_lf')
+            CALL iom_set_rstw_var_active('hdivn_lf')
          ENDIF
          !
@@ -279,5 +324,5 @@
       !!                               to the "baroclinic" level thickness.
       !!
-      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
+      !! ** Action  :  - hdivn_lf   : restoring towards full baroclinic divergence in z_tilde case
       !!               - tilde_e3t_a: after increment of vertical scale factor 
       !!                              in z_tilde case
@@ -289,10 +334,15 @@
       INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence
       !
+      LOGICAL                ::   ll_do_bclinic         ! local logical
       INTEGER                ::   ji, jj, jk            ! dummy loop indices
+      INTEGER                ::   ib, ib_bdy, ip, jp    !   "     "     "
       INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
-      REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars
-      LOGICAL                ::   ll_do_bclinic         ! local logical
+      INTEGER                ::   ncall
+      REAL(wp)               ::   z2dt  , z_tmin, z_tmax! local scalars                
+      REAL(wp)               ::   zalpha, zwgt          ! temporary scalars
+      REAL(wp)               ::   zdu, zdv, zramp
+      REAL(wp)               ::   ztra, zbtr, ztout, ztin, zfac, zmsku, zmskv
       REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t, ztu, ztv
       !!----------------------------------------------------------------------
       !
@@ -308,10 +358,21 @@
 
       ll_do_bclinic = .TRUE.
+      ncall = 1
+
       IF( PRESENT(kcall) ) THEN
-         IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
+         ! comment line below if tilda coordinate has to be computed at each call
+         ! This is not working yet
+!         IF (kcall == 2 .AND. (ln_vvl_ztilde.OR.ln_vvl_layer) ) ll_do_bclinic = .FALSE.
+         ncall = kcall  
+      ENDIF
+
+      IF( neuler == 0 .AND. kt == nit000 ) THEN
+         z2dt =  rdt
+      ELSE
+         z2dt = 2.0_wp * rdt
       ENDIF
 
       ! ******************************* !
-      ! After acale factors at t-points !
+      ! After scale factors at t-points !
       ! ******************************* !
       !                                                ! --------------------------------------------- !
@@ -325,175 +386,364 @@
       END DO
       !
-      IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
-         !                                                            ! ------baroclinic part------ !
-         ! I - initialization
-         ! ==================
-
-         ! 1 - barotropic divergence
-         ! -------------------------
-         zhdiv(:,:) = 0._wp
-         zht(:,:)   = 0._wp
+      IF( (ln_vvl_ztilde.OR.ln_vvl_layer) .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
+         !                                                               ! ------baroclinic part------ !
+         tilde_e3t_a(:,:,:) = 0.0_wp  ! tilde_e3t_a used to store tendency terms
+         un_td(:,:,:) = 0.0_wp        ! Transport corrections
+         vn_td(:,:,:) = 0.0_wp
+
+         zhdiv(:,:) = 0.
          DO jk = 1, jpkm1
             zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
-            zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
-         END DO
-         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
-
-         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
-         ! --------------------------------------------------
-         IF( ln_vvl_ztilde ) THEN
-            IF( kt > nit000 ) THEN
+         END DO
+         zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - tmask_i(:,:) ) 
+
+         ze3t(:,:,:) = 0._wp
+         IF( ln_rnf ) THEN
+            CALL sbc_rnf_div( ze3t )          ! runoffs 
+            DO jk=1,jpkm1
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ze3t(:,:,jk) * e3t_n(:,:,jk)
+            END DO   
+         ENDIF 
+
+         ! Thickness advection:
+         ! --------------------
+         ! Set advection velocities and source term
+         IF ( ln_vvl_ztilde ) THEN
+            IF ( ncall==1 ) THEN
+               zalpha  = rdt * 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
                DO jk = 1, jpkm1
-                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
-                     &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
-               END DO
+                  ztu(:,:,jk) = un(:,:,jk) * e3u_n(:,:,jk) * e2u(:,:)
+                  ztv(:,:,jk) = vn(:,:,jk) * e3v_n(:,:,jk) * e1v(:,:)
+                  ze3t(:,:,jk) = -tilde_e3t_a(:,:,jk) - e3t_n(:,:,jk) * zhdiv(:,:)
+               END DO
+               !
+               un_lf(:,:,:)  =    un_lf(:,:,:) * (1._wp - zalpha) + zalpha * ztu(:,:,:)
+               vn_lf(:,:,:)  =    vn_lf(:,:,:) * (1._wp - zalpha) + zalpha * ztv(:,:,:)
+            hdivn_lf(:,:,:)  = hdivn_lf(:,:,:) * (1._wp - zalpha) + zalpha * ze3t(:,:,:)
+! 2nd order filtering:
+!               un_lf2(:,:,:)  =    un_lf2(:,:,:) * (1._wp - zalpha) + zalpha * ztu(:,:,:)
+!               vn_lf2(:,:,:)  =    vn_lf2(:,:,:) * (1._wp - zalpha) + zalpha * ztv(:,:,:)
+!            hdivn_lf2(:,:,:)  = hdivn_lf2(:,:,:) * (1._wp - zalpha) + zalpha * ze3t(:,:,:)
+!               un_lf(:,:,:)   =     un_lf(:,:,:) * (1._wp - zalpha) + zalpha * un_lf2(:,:,:)
+!               vn_lf(:,:,:)   =     vn_lf(:,:,:) * (1._wp - zalpha) + zalpha * vn_lf2(:,:,:)
+!            hdivn_lf(:,:,:)   =  hdivn_lf(:,:,:) * (1._wp - zalpha) + zalpha * hdivn_lf2(:,:,:)
             ENDIF
+            !
+            DO jk = 1, jpkm1
+               ztu(:,:,jk) = (un(:,:,jk)-un_lf(:,:,jk)/e3u_n(:,:,jk)*r1_e2u(:,:))*umask(:,:,jk)
+               ztv(:,:,jk) = (vn(:,:,jk)-vn_lf(:,:,jk)/e3v_n(:,:,jk)*r1_e1v(:,:))*vmask(:,:,jk)
+               tilde_e3t_a(:,:,jk) =  tilde_e3t_a(:,:,jk) + hdivn_lf(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
+            END DO
+
+            !
+         ELSEIF ( ln_vvl_layer ) THEN
+            !
+            DO jk = 1, jpkm1
+               ztu(:,:,jk) = un(:,:,jk)
+               ztv(:,:,jk) = vn(:,:,jk)
+            END DO
+            !
          ENDIF
-
-         ! II - after z_tilde increments of vertical scale factors
-         ! =======================================================
-         tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms
-
-         ! 1 - High frequency divergence term
-         ! ----------------------------------
-         IF( ln_vvl_ztilde ) THEN     ! z_tilde case
+         !
+         ! Block fluxes through small layers:
+         DO jk=1,jpkm1
+            DO ji = 1, jpi
+               DO jj= 1, jpj
+                  zmsku = 0.5_wp * ( 1._wp + SIGN(1._wp, e3u_n(ji,jj,jk) - hsmall) )
+                  un_td(ji,jj,jk) = un_td(ji,jj,jk) - (1. - zmsku) * un(ji,jj,jk) * e3u_n(ji,jj,jk) * e2u(ji,jj)
+                  ztu(ji,jj,jk) = zmsku * ztu(ji,jj,jk)
+                  IF ( ln_vvl_ztilde ) un_lf(ji,jj,jk) = zmsku * un_lf(ji,jj,jk)
+                  !
+                  zmskv = 0.5_wp * ( 1._wp + SIGN(1._wp, e3v_n(ji,jj,jk) - hsmall) )
+                  vn_td(ji,jj,jk) = vn_td(ji,jj,jk) - (1. - zmskv) * vn(ji,jj,jk) * e3v_n(ji,jj,jk) * e1v(ji,jj)
+                  ztv(ji,jj,jk) = zmskv * ztv(ji,jj,jk)
+                  IF ( ln_vvl_ztilde ) vn_lf(ji,jj,jk) = zmskv * vn_lf(ji,jj,jk)
+               END DO
+            END DO
+         END DO      
+         !
+         ! Do advection
+         IF     (ln_vvl_adv_fct) THEN
+            CALL dom_vvl_adv_fct( kt, tilde_e3t_a, ztu, ztv )
+            !
+         ELSEIF (ln_vvl_adv_cn2) THEN
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
-            END DO
-         ELSE                         ! layer case
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     tilde_e3t_a(ji,jj,jk) =   tilde_e3t_a(ji,jj,jk) &
+                     & -(  e2u(ji,jj)*e3u_n(ji,jj,jk) * ztu(ji,jj,jk) - e2u(ji-1,jj  )*e3u_n(ji-1,jj  ,jk) * ztu(ji-1,jj  ,jk)       &
+                     &   + e1v(ji,jj)*e3v_n(ji,jj,jk) * ztv(ji,jj,jk) - e1v(ji  ,jj-1)*e3v_n(ji  ,jj-1,jk) * ztv(ji  ,jj-1,jk)  )    &
+                     & / ( e1t(ji,jj) * e2t(ji,jj) )
+                  END DO
+               END DO
+            END DO
+         ENDIF
+         ! 
+         ! Thickness anomaly diffusion:
+         ! -----------------------------
+         ztu(:,:,:) = 0.0_wp
+         ztv(:,:,:) = 0.0_wp
+
+         ze3t(:,:,1) = 0.e0
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               DO jk = 2, jpk
+                  ze3t(ji,jj,jk) =  ze3t(ji,jj,jk-1) + tilde_e3t_b(ji,jj,jk-1) * tmask(ji,jj,jk-1) 
+               END DO
+            END DO
+         END DO
+
+         IF ( ln_vvl_blp ) THEN  ! Bilaplacian
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
+               DO jj = 1, jpjm1                 ! First derivative (gradient)
+                  DO ji = 1, fs_jpim1   ! vector opt.                  
+!                     ztu(ji,jj,jk) =  umask(ji,jj,jk) * e2_e1u(ji,jj) &
+!                                     &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
+!                     ztv(ji,jj,jk) =  vmask(ji,jj,jk) * e1_e2v(ji,jj) & 
+!                                     &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
+                     ztu(ji,jj,jk) =  umask(ji,jj,jk) * e2_e1u(ji,jj) &
+                                     &  * ( ze3t(ji,jj,jk) - ze3t(ji+1,jj  ,jk) )
+                     ztv(ji,jj,jk) =  vmask(ji,jj,jk) * e1_e2v(ji,jj) & 
+                                     &  * ( ze3t(ji,jj,jk) - ze3t(ji  ,jj+1,jk) )
+                  END DO
+               END DO
+
+               DO jj = 2, jpjm1                 ! Second derivative (divergence) time the eddy diffusivity coefficient
+                  DO ji = fs_2, fs_jpim1 ! vector opt.
+                     zht(ji,jj) =  rn_ahe3_blp * r1_e1e2t(ji,jj) * (   ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   &
+                        &                                            + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)   )
+                  END DO
+               END DO
+
+               ! Open boundary conditions:
+               IF ( ln_bdy ) THEN
+                  DO ib_bdy=1, nb_bdy
+                     DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1)
+                        ji = idx_bdy(ib_bdy)%nbi(ib,1)
+                        jj = idx_bdy(ib_bdy)%nbj(ib,1)
+                        zht(ji,jj) = 0._wp
+                     END DO 
+                  END DO
+               END IF
+
+               CALL lbc_lnk( zht, 'T', 1. )     ! Lateral boundary conditions (unchanged sgn)
+
+               DO jj = 1, jpjm1                 ! third derivative (gradient)
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     ztu(ji,jj,jk) = umask(ji,jj,jk) * e2_e1u(ji,jj) * ( zht(ji+1,jj  ) - zht(ji,jj) )
+                     ztv(ji,jj,jk) = vmask(ji,jj,jk) * e1_e2v(ji,jj) * ( zht(ji  ,jj+1) - zht(ji,jj) )
+                  END DO
+               END DO
             END DO
          ENDIF
 
-         ! 2 - Restoring term (z-tilde case only)
-         ! ------------------
-         IF( ln_vvl_ztilde ) THEN
-            DO jk = 1, jpk
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
-            END DO
-         ENDIF
-
-         ! 3 - Thickness diffusion term
-         ! ----------------------------
-         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.
-                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           &
-                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
-                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           & 
-                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
-                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
-                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         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)
-               vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
-            END DO
-         END DO
-         DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes
+         IF ( ln_vvl_lap ) THEN    ! Laplacian
+            DO jk = 1, jpkm1                    ! First derivative (gradient)
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.                  
+!                     zdu =  rn_ahe3_lap * umask(ji,jj,jk) * e2_e1u(ji,jj) &
+!                         &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
+!                     zdv =  rn_ahe3_lap * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 
+!                         &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
+                     zdu =  rn_ahe3_lap * umask(ji,jj,jk) * e2_e1u(ji,jj) &
+                         &  * ( ze3t(ji,jj,jk) - ze3t(ji+1,jj  ,jk) )
+                     zdv =  rn_ahe3_lap * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 
+                         &  * ( ze3t(ji,jj,jk) - ze3t(ji  ,jj+1,jk) )
+                     ztu(ji,jj,jk) = ztu(ji,jj,jk) + zdu
+                     ztv(ji,jj,jk) = ztv(ji,jj,jk) + zdv
+                  END DO
+               END DO
+            END DO
+         ENDIF 
+
+         ! divergence of diffusive fluxes
+         DO jk = 1, jpkm1
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
-                     &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
+!                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   ztu(ji-1,jj  ,jk) - ztu(ji,jj,jk)    &
+!                     &                                          +     ztv(ji  ,jj-1,jk) - ztv(ji,jj,jk)    &
+!                     &                                            ) * r1_e1e2t(ji,jj)
+                  un_td(ji,jj,jk) = un_td(ji,jj,jk) + ztu(ji,jj,jk+1) - ztu(ji,jj,jk  ) 
+                  vn_td(ji,jj,jk) = vn_td(ji,jj,jk) + ztv(ji,jj,jk+1) - ztv(ji,jj,jk  ) 
+                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   ztu(ji-1,jj  ,jk+1) - ztu(ji,jj,jk+1)    &
+                     &                                               +ztv(ji  ,jj-1,jk+1) - ztv(ji,jj,jk+1)    &
+                     &                                               -ztu(ji-1,jj  ,jk  ) + ztu(ji,jj,jk  )    &
+                     &                                               -ztv(ji  ,jj-1,jk  ) + ztv(ji,jj,jk  )    &
                      &                                            ) * r1_e1e2t(ji,jj)
                END DO
             END DO
          END DO
-         !                       ! d - thickness diffusion transport: boundary conditions
-         !                             (stored for tracer advction and continuity equation)
-         CALL lbc_lnk_multi( un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
-
-         ! 4 - Time stepping of baroclinic scale factors
-         ! ---------------------------------------------
+
+ 
+!         un_td(:,:,:) = un_td(:,:,:) + ztu(:,:,:)
+!         vn_td(:,:,:) = vn_td(:,:,:) + ztv(:,:,:)
+
+         CALL lbc_lnk( un_td , 'U' , -1.)
+         CALL lbc_lnk( vn_td , 'V' , -1.) 
+         !
+         CALL dom_vvl_ups_cor( kt, tilde_e3t_a, un_td, vn_td )
+
+!         IF ( ln_vvl_ztilde ) THEN
+!            ztu(:,:,:) = -un_lf(:,:,:)
+!            ztv(:,:,:) = -vn_lf(:,:,:)
+!            CALL dom_vvl_ups_cor( kt, tilde_e3t_a, ztu, ztv )
+!         ENDIF
+         !
+         ! Remove "external thickness" tendency:
+         DO jk = 1, jpkm1
+            tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) +   e3t_n(:,:,jk) * zhdiv(:,:)
+         END DO 
+                   
          ! Leapfrog time stepping
          ! ~~~~~~~~~~~~~~~~~~~~~~
-         IF( neuler == 0 .AND. kt == nit000 ) THEN
-            z2dt =  rdt
-         ELSE
-            z2dt = 2.0_wp * rdt
+         zramp = 1._wp
+         IF ((.NOT.ln_rstart).AND.ln_vvl_ramp) zramp = MIN(MAX( ((kt-nit000)*rdt)/(rn_day_ramp*rday),0._wp),1._wp)
+
+         DO jk=1,jpkm1
+            tilde_e3t_a(:,:,jk) = tilde_e3t_b(:,:,jk) + z2dt * tmask(:,:,jk) * tilde_e3t_a(:,:,jk) &
+                               & * tildemask(:,:) * zramp
+         END DO
+
+         ! Ensure layer separation:
+         ! ~~~~~~~~~~~~~~~~~~~~~~~~
+         IF ( ln_vvl_regrid ) CALL dom_vvl_regrid( kt )
+  
+         ! Boundary conditions:
+         ! ~~~~~~~~~~~~~~~~~~~~
+         IF ( ln_bdy ) THEN
+            DO ib_bdy = 1, nb_bdy
+               DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1)
+!!               DO ib = 1, idx_bdy(ib_bdy)%nblen(1)
+                  ji = idx_bdy(ib_bdy)%nbi(ib,1)
+                  jj = idx_bdy(ib_bdy)%nbj(ib,1)
+                  zwgt = idx_bdy(ib_bdy)%nbw(ib,1)
+                  ip = bdytmask(ji+1,jj  ) - bdytmask(ji-1,jj  )
+                  jp = bdytmask(ji  ,jj+1) - bdytmask(ji  ,jj-1)
+                  DO jk = 1, jpkm1  
+                     tilde_e3t_a(ji,jj,jk) = 0.e0
+!!                     tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) * (1._wp - zwgt)
+!!                     tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji+ip,jj+jp,jk) * tmask(ji+ip,jj+jp,jk)
+                  END DO
+               END DO 
+            END DO
          ENDIF
-         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp )
-         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
-
-         ! Maximum deformation control
-         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         ze3t(:,:,jpk) = 0._wp
+
+         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )        
+        
+      ENDIF
+
+      IF( ln_vvl_ztilde.AND.( ncall==1))  THEN
+         zalpha  = rdt * 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
+         !
+         ! divergence of diffusive fluxes
          DO jk = 1, jpkm1
-            ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
-         END DO
-         z_tmax = MAXVAL( ze3t(:,:,:) )
-         IF( lk_mpp )   CALL mpp_max( z_tmax )                 ! max over the global domain
-         z_tmin = MINVAL( ze3t(:,:,:) )
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ze3t(ji,jj,jk) = (  un_td(ji,jj,jk) - un_td(ji-1,jj  ,jk) &
+                                 &  + vn_td(ji,jj,jk) - vn_td(ji  ,jj-1,jk) &
+                                 & ) / ( e1t(ji,jj) * e2t(ji,jj) )
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk( ze3t(:,:,:), 'T', 1. )
+         hdivn_lf(:,:,:)  =  hdivn_lf(:,:,:)  + zalpha * ze3t(:,:,:) 
+! 2nd order filtering:
+!         hdivn_lf2(:,:,:)  = hdivn_lf2(:,:,:) + zalpha * ze3t(:,:,:)
+!         hdivn_lf(:,:,:)  =  hdivn_lf(:,:,:)  + zalpha * zalpha * ze3t(:,:,:)
+      ENDIF
+
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
+      !                                             ! ---baroclinic part--------- !
+         DO jk = 1, jpkm1
+            e3t_a(:,:,jk) = e3t_a(:,:,jk) + (tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk))                     
+         END DO
+      ENDIF
+
+      IF( ln_vvl_dbg.AND.(ln_vvl_ztilde .OR. ln_vvl_layer) ) THEN   ! - ML - test: control prints for debuging
+         !
+         zht(:,:) = 0.0_wp
+         DO jk = 1, jpkm1
+            zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
+         END DO
+         IF( lwp ) WRITE(numout, *) 'kt = ', kt
+         IF( lwp ) WRITE(numout, *) 'ncall = ', ncall
+         IF( lwp ) WRITE(numout, *) 'll_do_bclinic', ll_do_bclinic 
+         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
+            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ), mask = tmask(:,:,1) == 1.e0 )
+            IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain
+            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
+         END IF
+         !
+         z_tmin = MINVAL( e3t_n(:,:,:)  )
          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 >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN
+         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(e3t_n) =', z_tmin
+         IF ( z_tmin .LE. 0._wp ) THEN
             IF( lk_mpp ) THEN
-               CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) )
-               CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
+               CALL mpp_minloc(e3t_n(:,:,:), tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
             ELSE
-               ijk_max = MAXLOC( ze3t(:,:,:) )
-               ijk_max(1) = ijk_max(1) + nimpp - 1
-               ijk_max(2) = ijk_max(2) + njmpp - 1
-               ijk_min = MINLOC( ze3t(:,:,:) )
+               ijk_min = MINLOC( e3t_n(:,:,:)  )
                ijk_min(1) = ijk_min(1) + nimpp - 1
                ijk_min(2) = ijk_min(2) + njmpp - 1
             ENDIF
             IF (lwp) THEN
-               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
-               WRITE(numout, *) 'at i, j, k=', ijk_max
-               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
+               WRITE(numout, *) 'Negative scale factor, e3t_n =', z_tmin
                WRITE(numout, *) 'at i, j, k=', ijk_min            
-               CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')
+               CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor')
+            ENDIF
+         ENDIF         
+         !
+         z_tmin = MINVAL( e3u_n(:,:,:))
+         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
+         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(e3u_n) =', z_tmin
+         IF ( z_tmin .LE. 0._wp ) THEN
+            IF( lk_mpp ) THEN
+               CALL mpp_minloc(e3u_n(:,:,:), umask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
+            ELSE
+               ijk_min = MINLOC( e3u_n(:,:,:)  )
+               ijk_min(1) = ijk_min(1) + nimpp - 1
+               ijk_min(2) = ijk_min(2) + njmpp - 1
+            ENDIF
+            IF (lwp) THEN
+               WRITE(numout, *) 'Negative scale factor, e3u_n =', z_tmin
+               WRITE(numout, *) 'at i, j, k=', ijk_min            
+               CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor')
+            ENDIF
+         ENDIF 
+         !
+         z_tmin = MINVAL( e3v_n(:,:,:) )
+         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
+         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(e3v_n) =', z_tmin
+         IF ( z_tmin .LE. 0._wp ) THEN
+            IF( lk_mpp ) THEN
+               CALL mpp_minloc(e3v_n(:,:,:), vmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
+            ELSE
+               ijk_min = MINLOC( e3v_n(:,:,:), mask = vmask(:,:,:) == 1.e0   )
+               ijk_min(1) = ijk_min(1) + nimpp - 1
+               ijk_min(2) = ijk_min(2) + njmpp - 1
+            ENDIF
+            IF (lwp) THEN
+               WRITE(numout, *) 'Negative scale factor, e3v_n =', z_tmin
+               WRITE(numout, *) 'at i, j, k=', ijk_min            
+               CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor')
+            ENDIF
+         ENDIF 
+         !
+         z_tmin = MINVAL( e3f_n(:,:,:))
+         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
+         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(e3f_n) =', z_tmin
+         IF ( z_tmin .LE. 0._wp ) THEN
+            IF( lk_mpp ) THEN
+               CALL mpp_minloc(e3f_n(:,:,:), fmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
+            ELSE
+               ijk_min = MINLOC( e3f_n(:,:,:) )
+               ijk_min(1) = ijk_min(1) + nimpp - 1
+               ijk_min(2) = ijk_min(2) + njmpp - 1
+            ENDIF
+            IF (lwp) THEN
+               WRITE(numout, *) 'Negative scale factor, e3f_n =', z_tmin
+               WRITE(numout, *) 'at i, j, k=', ijk_min            
+               CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor')
             ENDIF
          ENDIF
-         ! - ML - end test
-         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
-         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
-         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
-
-         !
-         ! "tilda" change in the after scale factor
-         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         DO jk = 1, jpkm1
-            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
-         END DO
-         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
-         ! ==================================================================================
-         ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
-         ! - ML - baroclinicity error should be better treated in the future
-         !        i.e. locally and not spread over the water column.
-         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
-         zht(:,:) = 0.
-         DO jk = 1, jpkm1
-            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
-         DO jk = 1, jpkm1
-            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
-         END DO
-
-      ENDIF
-
-      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
-      !                                           ! ---baroclinic part--------- !
-         DO jk = 1, jpkm1
-            e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
-         END DO
-      ENDIF
-
-      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
-         !
-         IF( lwp ) WRITE(numout, *) 'kt =', kt
-         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
-            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
-            IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain
-            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
-         END IF
          !
          zht(:,:) = 0.0_wp
@@ -503,33 +753,58 @@
          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(e3t_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(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
+            zht(:,:) = zht(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
+         END DO
+         zwu(:,:) = 0._wp
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zwu(ji,jj) = 0.5_wp * umask(ji,jj,1) * r1_e1e2u(ji,jj)                             &
+                        &          * ( e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj) )
+            END DO
+         END DO
+         CALL lbc_lnk( zwu(:,:), 'U', 1._wp )
+         z_tmax = MAXVAL( ssumask(:,:) * ssumask(:,:) * ABS( hu_0(:,:) + zwu(:,:) - 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(e3t_a))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(hu_0+sshu_n-SUM(e3u_n))) =', z_tmax
          !
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
+            zht(:,:) = zht(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
+         END DO
+         zwv(:,:) = 0._wp
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zwv(ji,jj) = 0.5_wp * vmask(ji,jj,1) * r1_e1e2v(ji,jj)                             &
+                        &          * ( e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1) )
+            END DO
+         END DO
+         CALL lbc_lnk( zwv(:,:), 'V', 1._wp )
+         z_tmax = MAXVAL( ssvmask(:,:) * ssvmask(:,:) * ABS( hv_0(:,:) + zwv(:,:) - 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(e3t_b))) =', z_tmax
-         !
-         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) )
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(hv_0+sshv_n-SUM(e3v_n))) =', z_tmax
+         !
+         zht(:,:) = 0.0_wp
+         DO jk = 1, jpkm1
+            DO jj = 1, jpjm1
+               zht(:,jj) = zht(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk)*umask(:,jj+1,jk)
+            END DO
+         END DO
+         zwu(:,:) = 0._wp
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zwu(ji,jj) = 0.25_wp * umask(ji,jj,1) * umask(ji,jj+1,1) * r1_e1e2f(ji,jj)  &
+                        &          * (  e1e2t(ji  ,jj) * sshn(ji  ,jj) + e1e2t(ji  ,jj+1) * sshn(ji  ,jj+1) &
+                        &             + e1e2t(ji+1,jj) * sshn(ji+1,jj) + e1e2t(ji+1,jj+1) * sshn(ji+1,jj+1) )
+            END DO
+         END DO
+         CALL lbc_lnk( zht(:,:), 'F', 1._wp )
+         CALL lbc_lnk( zwu(:,:), 'F', 1._wp )
+         z_tmax = MAXVAL( fmask(:,:,1) * fmask(:,:,1) * ABS( hf_0(:,:) + zwu(:,:) - zht(:,:) ) )
          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax
-         !
-         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) )
-         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax
-         !
-         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) )
-         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(hf_0+sshf_n-SUM(e3f_n))) =', z_tmax
+         !
       END IF
 
@@ -699,7 +974,7 @@
       !!----------------------------------------------------------------------
       !
-      nmet = 0  ! Original method (Surely wrong)
-!     nmet = 1  ! Interface interpolation
-!     nmet = 2  ! Internal interfaces interpolation only, spread barotropic increment
+!     nmet = 0  ! Original method (Surely wrong)
+     nmet = 2  ! Interface interpolation
+!      nmet = 2  ! Internal interfaces interpolation only, spread barotropic increment
 !     Note that we kept surface weighted interpolation for barotropic increment to be compliant
 !     with what is done in surface pressure module.
@@ -711,6 +986,6 @@
       END IF
       !
-      ztap   = 0.0_wp   ! Minimum fraction of T-point thickness at cell interfaces
-      zsmall = 1.e-3_wp ! Minimum thickness at U or V points (m)
+      ztap   = 0._wp   ! Minimum fraction of T-point thickness at cell interfaces
+      zsmall = 1.e-8_wp ! Minimum thickness at U or V points (m)
       !
       IF ( (nmet==1).OR.(nmet==2) ) THEN
@@ -744,4 +1019,5 @@
                END DO 
             ENDIF 
+            zw(:,:,:) = (zw(:,:,:) - gdepw_0(:,:,:))*tmask(:,:,:)
             !
          END SELECT
@@ -768,19 +1044,20 @@
                   ! Correction at last level:
                   jkbot = mbku(ji,jj)
-                  zdo   = hu_0(ji,jj)
+                  zdo = 0._wp
                   DO jk=jkbot,1,-1
-                     zup = 0.5_wp * ( zw(ji  ,jj,jk) + zw(ji+1,jj,jk) ) 
+                     zup = 0.5_wp * ( e1e2t(ji  ,jj)*zw(ji  ,jj,jk) &
+                         &          + e1e2t(ji+1,jj)*zw(ji+1,jj,jk) ) * r1_e1e2u(ji,jj)
                      !
                      ! If there is a step, taper bottom interface:
-                     IF ((hu_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji+1,jj) ) ).AND.(zup>zdo)) THEN
-                        IF ( ht_0(ji+1,jj) < ht_0(ji,jj) ) THEN
-                           zmin = ztap * (zw(ji+1,jj,jk+1)-zw(ji+1,jj,jk))
-                        ELSE
-                           zmin = ztap * (zw(ji  ,jj,jk+1)-zw(ji  ,jj,jk))
-                        ENDIF
-                        zup = MIN(zup, zdo-zmin)
-                     ENDIF
-                     zup = MIN(zup, zdo-zsmall)
-                     pe3_out(ji,jj,jk) = (zdo - zup - e3u_0(ji,jj,jk)) * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )
+!                     IF ((hu_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji+1,jj) ) ).AND.(jk==jkbot)) THEN
+!                        IF ( ht_0(ji+1,jj) < ht_0(ji,jj) ) THEN
+!                           zmin = ztap * (zw(ji+1,jj,jk+1)-zw(ji+1,jj,jk))
+!                        ELSE
+!                           zmin = ztap * (zw(ji  ,jj,jk+1)-zw(ji  ,jj,jk))
+!                        ENDIF
+!                        zup = MIN(zup, zdo-zmin)
+!                     ENDIF
+                     zup = MIN(zup, zdo+e3u_0(ji,jj,jk)-zsmall)
+                     pe3_out(ji,jj,jk) = (zdo - zup) * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )
                      zdo = zup
                   END DO
@@ -795,5 +1072,5 @@
                      pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                   &
                            &               + ( pe3_out(ji,jj,jk) + e3u_0(ji,jj,jk) )               & 
-                           &               / ( hu_0(ji,jj) + 1._wp - umask(ji,jj,1) )              &
+                           &               / ( hu_0(ji,jj) + 1._wp - ssumask(ji,jj) )              &
                            &               * 0.5_wp * r1_e1e2u(ji,jj)                              &
                            &               * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )        &
@@ -824,19 +1101,20 @@
                   ! Correction at last level:
                   jkbot = mbkv(ji,jj)
-                  zdo   = hv_0(ji,jj)
+                  zdo = 0._wp
                   DO jk=jkbot,1,-1
-                     zup = 0.5_wp * ( zw(ji,jj  ,jk) + zw(ji,jj+1,jk) ) 
+                     zup = 0.5_wp * ( e1e2t(ji,jj  ) * zw(ji,jj  ,jk) & 
+                         &          + e1e2t(ji,jj+1) * zw(ji,jj+1,jk) ) * r1_e1e2v(ji,jj)
                      !
                      ! If there is a step, taper bottom interface:
-                     IF ((hv_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji,jj+1) ) ).AND.(zup>zdo)) THEN
-                        IF ( ht_0(ji,jj+1) < ht_0(ji,jj) ) THEN
-                           zmin = ztap * (zw(ji,jj+1,jk+1)-zw(ji,jj+1,jk))
-                        ELSE
-                           zmin = ztap * (zw(ji  ,jj,jk+1)-zw(ji  ,jj,jk))
-                        ENDIF
-                        zup = MIN(zup, zdo-zmin)
-                     ENDIF
-                     zup = MIN(zup, zdo-zsmall)
-                     pe3_out(ji,jj,jk) = (zdo - zup - e3v_0(ji,jj,jk)) * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )
+!                     IF ((hv_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji,jj+1) ) ).AND.(jk==jkbot)) THEN
+!                        IF ( ht_0(ji,jj+1) < ht_0(ji,jj) ) THEN
+!                           zmin = ztap * (zw(ji,jj+1,jk+1)-zw(ji,jj+1,jk))
+!                        ELSE
+!                           zmin = ztap * (zw(ji  ,jj,jk+1)-zw(ji  ,jj,jk))
+!                        ENDIF
+!                        zup = MIN(zup, zdo-zmin)
+!                     ENDIF
+                     zup = MIN(zup, zdo + e3v_0(ji,jj,jk) - zsmall)
+                     pe3_out(ji,jj,jk) = (zdo - zup) * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )
                      zdo = zup
                   END DO
@@ -851,5 +1129,5 @@
                      pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                                       &
                            &               + ( pe3_out(ji,jj,jk) + e3v_0(ji,jj,jk) )                                   & 
-                           &               / ( hv_0(ji,jj) + 1._wp - vmask(ji,jj,1) )                                  &
+                           &               / ( hv_0(ji,jj) + 1._wp - ssvmask(ji,jj) )                                  &
                            &               * 0.5_wp * r1_e1e2v(ji,jj)                                                  &
                                            * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )                            &
@@ -883,31 +1161,31 @@
                   ! bottom correction:
                   jkbot = MIN(mbku(ji,jj), mbku(ji,jj+1)) 
-                  zdo = hf_0(ji,jj)
+                  zdo = 0._wp
                   DO jk=jkbot,1,-1
-                     zup =  0.25_wp * (  zw(ji  ,jj  ,jk)  & 
-                           &           + zw(ji+1,jj  ,jk)  & 
-                           &           + zw(ji  ,jj+1,jk)  & 
-                           &           + zw(ji+1,jj+1,jk)  )
+                     zup =  0.25_wp * (  e1e2t(ji  ,jj  ) * zw(ji  ,jj  ,jk)  & 
+                           &           + e1e2t(ji+1,jj  ) * zw(ji+1,jj  ,jk)  & 
+                           &           + e1e2t(ji  ,jj+1) * zw(ji  ,jj+1,jk)  & 
+                           &           + e1e2t(ji+1,jj+1) * zw(ji+1,jj+1,jk)  ) *    r1_e1e2f(ji,jj)
                      !
                      ! If there is a step, taper bottom interface:
-                     IF ((hf_0(ji,jj) < 0.5_wp * ( hu_0(ji,jj  ) + hu_0(ji,jj+1) ) ).AND.(zup>zdo)) THEN
-                        IF ( hu_0(ji,jj+1) < hu_0(ji,jj) ) THEN
-                           IF ( ht_0(ji+1,jj+1) < ht_0(ji  ,jj+1) ) THEN
-                              zmin = ztap * (zw(ji+1,jj+1,jk+1)-zw(ji+1,jj+1,jk))
-                           ELSE
-                              zmin = ztap * (zw(ji  ,jj+1,jk+1)-zw(ji  ,jj+1,jk))
-                           ENDIF
-                        ELSE
-                           IF ( ht_0(ji+1,jj  ) < ht_0(ji  ,jj  ) ) THEN
-                              zmin = ztap * (zw(ji+1,jj  ,jk+1)-zw(ji+1,jj  ,jk))
-                           ELSE
-                              zmin = ztap * (zw(ji  ,jj  ,jk+1)-zw(ji  ,jj  ,jk))
-                           ENDIF
-                        ENDIF
-                        zup = MIN(zup, zdo-zmin)
-                     ENDIF
-                     zup = MIN(zup, zdo-zsmall)
+!                     IF ((hf_0(ji,jj) < 0.5_wp * ( hu_0(ji,jj  ) + hu_0(ji,jj+1) ) ).AND.(jk==jkbot)) THEN
+!                        IF ( hu_0(ji,jj+1) < hu_0(ji,jj) ) THEN
+!                           IF ( ht_0(ji+1,jj+1) < ht_0(ji  ,jj+1) ) THEN
+!                              zmin = ztap * (zw(ji+1,jj+1,jk+1)-zw(ji+1,jj+1,jk))
+!                           ELSE
+!                              zmin = ztap * (zw(ji  ,jj+1,jk+1)-zw(ji  ,jj+1,jk))
+!                           ENDIF
+!                        ELSE
+!                           IF ( ht_0(ji+1,jj  ) < ht_0(ji  ,jj  ) ) THEN
+!                              zmin = ztap * (zw(ji+1,jj  ,jk+1)-zw(ji+1,jj  ,jk))
+!                           ELSE
+!                              zmin = ztap * (zw(ji  ,jj  ,jk+1)-zw(ji  ,jj  ,jk))
+!                           ENDIF
+!                        ENDIF
+!                        zup = MIN(zup, zdo-zmin)
+!                     ENDIF
+                     zup = MIN(zup, zdo+e3f_0(ji,jj,jk)-zsmall)
                      !
-                     pe3_out(ji,jj,jk) = ( zdo - zup - e3f_0(ji,jj,jk) ) & 
+                     pe3_out(ji,jj,jk) = ( zdo - zup ) & 
                                       & *( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd )
                      zdo = zup
@@ -1006,5 +1284,5 @@
             id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
             id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
-            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
+            id5 = iom_varid( numror, 'hdivn_lf', ldstop = .FALSE. )
             !                             ! --------- !
             !                             ! all cases !
@@ -1068,7 +1346,7 @@
                   !                       ! ------------ !
                   IF( id5 > 0 ) THEN  ! required array exists
-                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios )
+                     CALL iom_get( numror, jpdom_autoglo, 'hdivn_lf', hdivn_lf(:,:,:), ldxios = lrxios )
                   ELSE                ! array is missing
-                     hdiv_lf(:,:,:) = 0.0_wp
+                     hdivn_lf(:,:,:) = 0.0_wp
                   ENDIF
                ENDIF
@@ -1140,5 +1418,9 @@
                tilde_e3t_b(:,:,:) = 0._wp
                tilde_e3t_n(:,:,:) = 0._wp
-               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
+               IF( ln_vvl_ztilde ) THEN 
+                  hdivn_lf(:,:,:) = 0._wp 
+                  un_lf(:,:,:) = 0._wp 
+                  vn_lf(:,:,:) = 0._wp 
+               ENDIF
             END IF
          ENDIF
@@ -1162,5 +1444,5 @@
          IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
             !                                        ! ------------ !
-            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lwxios)
+            CALL iom_rstput( kt, nitrst, numrow, 'hdivn_lf', hdivn_lf(:,:,:), ldxios = lwxios)
          ENDIF
          !
@@ -1179,9 +1461,16 @@
       !!----------------------------------------------------------------------
       INTEGER ::   ioptio, ios
-      !!
-      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
-         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
-         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
-      !!---------------------------------------------------------------------- 
+
+      NAMELIST/nam_vvl/ ln_vvl_zstar               , ln_vvl_ztilde                       , &
+                      & ln_vvl_layer               , ln_vvl_ztilde_as_zstar              , &
+                      & ln_vvl_zstar_at_eqtor      , ln_vvl_zstar_on_shelf               , &
+                      & ln_vvl_adv_cn2             , ln_vvl_adv_fct                      , &
+                      & ln_vvl_lap                 , ln_vvl_blp                          , &
+                      & rn_ahe3_lap                , rn_ahe3_blp                         , &
+                      & rn_rst_e3t                 , rn_lf_cutoff                        , &
+                      & ln_vvl_regrid                                                    , &
+                      & ln_vvl_ramp                , rn_day_ramp                         , &
+                      & ln_vvl_dbg   ! not yet implemented: ln_vvl_kepe
+      !!----------------------------------------------------------------------  
       !
       REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist : 
@@ -1197,25 +1486,48 @@
          WRITE(numout,*) 'dom_vvl_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
-         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
-         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
-         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
+         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate'
+         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
+         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
+         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
+         WRITE(numout,*) '              ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
+         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
+         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
          WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
-         WRITE(numout,*) '      !'
-         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
-         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
+         WRITE(numout,*) '      ztilde on shelves          ln_vvl_zstar_on_shelf  = ', ln_vvl_zstar_on_shelf
+         WRITE(numout,*) '           Namelist nam_vvl : thickness advection scheme'
+         WRITE(numout,*) '              2nd order                  ln_vvl_adv_cn2 = ', ln_vvl_adv_cn2
+         WRITE(numout,*) '              2nd order FCT              ln_vvl_adv_fct = ', ln_vvl_adv_fct                    
+         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion scheme'
+         WRITE(numout,*) '              Laplacian                  ln_vvl_lap     = ', ln_vvl_lap
+         WRITE(numout,*) '              Bilaplacian                ln_vvl_blp     = ', ln_vvl_blp
+         WRITE(numout,*) '              Laplacian   coefficient    rn_ahe3_lap    = ', rn_ahe3_lap
+         WRITE(numout,*) '              Bilaplacian coefficient    rn_ahe3_blp    = ', rn_ahe3_blp
+         WRITE(numout,*) '           Namelist nam_vvl : layers regriding'
+         WRITE(numout,*) '              ln_vvl_regrid                             = ', ln_vvl_regrid
+         WRITE(numout,*) '           Namelist nam_vvl : linear ramp at startup'
+         WRITE(numout,*) '              ln_vvl_ramp                               = ', ln_vvl_ramp
+         WRITE(numout,*) '              rn_day_ramp                               = ', rn_day_ramp
          IF( ln_vvl_ztilde_as_zstar ) THEN
-            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
-            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
-            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
-            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
-            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
-            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rdt'
+            WRITE(numout,*) '           ztilde running in zstar emulation mode; '
+            WRITE(numout,*) '           ignoring namelist timescale parameters and using:'
+            WRITE(numout,*) '                 hard-wired : z-tilde to zstar restoration timescale (days)'
+            WRITE(numout,*) '                                         rn_rst_e3t     =    0.0'
+            WRITE(numout,*) '                 hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
+            WRITE(numout,*) '                                         rn_lf_cutoff   =    1.0/rdt'
          ELSE
-            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
-            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
+            WRITE(numout,*) '           Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
+            WRITE(numout,*) '                                         rn_rst_e3t     = ', rn_rst_e3t
+            WRITE(numout,*) '           Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
+            WRITE(numout,*) '                                         rn_lf_cutoff   = ', rn_lf_cutoff
          ENDIF
-         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
+         WRITE(numout,*) '           Namelist nam_vvl : debug prints'
+         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg
+      ENDIF
+      !
+      IF ( ln_vvl_ztilde.OR.ln_vvl_layer ) THEN 
+         ioptio = 0                      ! Choose one advection scheme at most
+         IF( ln_vvl_adv_cn2         )        ioptio = ioptio + 1
+         IF( ln_vvl_adv_fct         )        ioptio = ioptio + 1
+         IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE thickness advection scheme in namelist nam_vvl' )
       ENDIF
       !
@@ -1237,4 +1549,12 @@
       ENDIF
       !
+      ! Use of "shelf horizon depths" should be allowed with s-z coordinates, but we restrict it to zco and zps 
+      ! for the time being
+      IF ( ln_sco ) THEN 
+        ll_shorizd=.FALSE.
+      ELSE
+        ll_shorizd=.TRUE.
+      ENDIF
+      !
 #if defined key_agrif
       IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
@@ -1243,4 +1563,872 @@
    END SUBROUTINE dom_vvl_ctl
 
+   SUBROUTINE dom_vvl_regrid( kt )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_regrid  ***
+      !!                   
+      !! ** Purpose :  Ensure "well-behaved" vertical grid
+      !!
+      !! ** Method  :  More or less adapted from references below.
+      !!regrid
+      !! ** Action  :  Ensure that thickness are above a given value, spaced enough
+      !!               and revert to Eulerian coordinates near the bottom.      
+      !!
+      !! References :  Bleck, R. and S. Benjamin, 1993: Regional Weather Prediction
+      !!               with a Model Combining Terrain-following and Isentropic
+      !!               coordinates. Part I: Model Description. Monthly Weather Rev.,
+      !!               121, 1770-1785.
+      !!               Toy, M., 2011: Incorporating Condensational Heating into a
+      !!               Nonhydrostatic Atmospheric Model Based on a Hybrid Isentropic-
+      !!               Sigma Vertical Coordinate. Monthly Weather Rev., 139, 2940-2954.
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( in )               :: kt         ! time step
+
+      !! * Local declarations
+      INTEGER                             :: ji, jj, jk ! dummy loop indices
+      LOGICAL                             :: ll_chk_bot2top, ll_chk_top2bot, ll_lapdiff_cond
+      LOGICAL                             :: ll_zdiff_cond, ll_blpdiff_cond
+      INTEGER                             :: jkbot
+      REAL(wp)                            :: zh_min, zh_0, zh2, zdiff, zh_max, ztmph, ztmpd
+      REAL(wp)                            :: zufim1, zufi, zvfjm1, zvfj, dzmin_int, dzmin_surf
+      REAL(wp)                            :: zh_new, zh_old, zh_bef, ztmp, ztmp1, z2dt, zh_up, zh_dwn
+      REAL(wp)                            :: zeu2, zev2, zfrch_stp, zfrch_rel, zfrac_bot, zscal_bot
+      REAL(wp)                            :: zhdiff, zhdiff2, zvdiff, zhlim, zhlim2, zvlim
+      REAL(wp), DIMENSION(jpi,jpj)        :: zdw, zwu, zwv
+      REAL(wp), DIMENSION(jpi,jpj,jpk)    :: zwdw, zwdw_b
+      !!----------------------------------------------------------------------
+
+      IF( ln_timing )  CALL timing_start('dom_vvl_regrid')
+      !
+      !
+      ! Some user defined parameters below:
+      ll_chk_bot2top   = .TRUE.
+      ll_chk_top2bot   = .TRUE.
+      dzmin_int  = 0.1_wp   ! Absolute minimum depth in the interior (in meters)
+      dzmin_surf = 1.0_wp   ! Absolute minimum depth at the surface (in meters)
+      zfrch_stp  = 0.5_wp   ! Maximum fractionnal thickness change in one time step (<= 1.)
+      zfrch_rel  = 0.4_wp   ! Maximum relative thickness change in the vertical (<= 1.)
+      zfrac_bot  = 0.05_wp  ! Fraction of bottom level allowed to change
+      zscal_bot  = 2.0_wp   ! Depth lengthscale
+      ll_zdiff_cond    = .TRUE.  ! Conditionnal vertical diffusion of interfaces
+         zvdiff  = 0.2_wp   ! m
+         zvlim   = 0.5_wp   ! max d2h/dh
+      ll_lapdiff_cond  = .TRUE.  ! Conditionnal Laplacian diffusion of interfaces
+         zhdiff  = 0.01_wp   ! ad.
+         zhlim   = 0.03_wp  ! ad. max lap(z)*e1
+      ll_blpdiff_cond  = .FALSE.  ! Conditionnal Bilaplacian diffusion of interfaces
+         zhdiff2 = 0.2_wp   ! ad.
+!         zhlim2  = 3.e-11_wp  !  max bilap(z)
+         zhlim2  = 0.01_wp  ! ad. max bilap(z)*e1**3
+      ! --------------------------------------------------------------------------------------- 
+      !
+      ! Set arrays determining maximum vertical displacement at the bottom:
+      !--------------------------------------------------------------------
+      IF ( kt==nit000 ) THEN
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               jk = MIN(mbkt(ji,jj), mbkt(ji+1,jj), mbkt(ji-1,jj), mbkt(ji,jj+1), mbkt(ji,jj-1))
+               jk = MIN(jk,mbkt(ji-1,jj-1), mbkt(ji-1,jj+1), mbkt(ji+1,jj+1), mbkt(ji+1,jj-1))
+               i_int_bot(ji,jj) = jk
+            END DO
+         END DO
+         dsm(:,:) = REAL( i_int_bot(:,:), wp )   ;   CALL lbc_lnk(dsm(:,:),'T',1.)
+         i_int_bot(:,:) = MAX( INT( dsm(:,:) ), 1 )
+
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zdw(ji,jj) = MAX(ABS(ht_0(ji,jj)-ht_0(ji+1,jj))*umask(ji  ,jj,1), &
+                          &     ABS(ht_0(ji,jj)-ht_0(ji-1,jj))*umask(ji-1,jj,1), &
+                          &     ABS(ht_0(ji,jj)-ht_0(ji,jj+1))*vmask(ji,jj  ,1), &
+                          &     ABS(ht_0(ji,jj)-ht_0(ji,jj-1))*vmask(ji,jj-1,1)  )
+                zdw(ji,jj) = MAX(zscal_bot * zdw(ji,jj), rsmall )
+            END DO
+         END DO
+         CALL lbc_lnk( zdw(:,:), 'T', 1. )
+
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               dsm(ji,jj) = 1._wp/16._wp * (            zdw(ji-1,jj-1) + zdw(ji+1,jj-1)      &
+                  &                           +         zdw(ji-1,jj+1) + zdw(ji+1,jj+1)      &
+                  &                           + 2._wp*( zdw(ji  ,jj-1) + zdw(ji-1,jj  )      &
+                  &                           +         zdw(ji+1,jj  ) + zdw(ji  ,jj+1) )    &
+                  &                           + 4._wp*  zdw(ji  ,jj  )                       )
+            END DO
+         END DO         
+
+         CALL lbc_lnk( dsm(:,:), 'T', 1. )   
+     
+         IF (ln_zps) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  jk = i_int_bot(ji,jj)
+                  hsm(ji,jj) = zfrac_bot * e3w_1d(jk)
+!                  dsm(ji,jj) = MAX(dsm(ji,jj), 0.1_wp*ht_0(ji,jj))
+               END DO
+            END DO
+         ELSE
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  jk = i_int_bot(ji,jj)
+                  hsm(ji,jj) = zfrac_bot * e3w_0(ji,jj,jk)
+!                  dsm(ji,jj) = MAX(dsm(ji,jj), 0.1_wp*ht_0(ji,jj))
+               END DO
+            END DO
+         ENDIF
+      END IF
+
+      ! Provisionnal interface depths:
+      !-------------------------------
+      zwdw(:,:,1) = 0.e0
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            DO jk = 2, jpk
+               zwdw(ji,jj,jk) =  zwdw(ji,jj,jk-1) + & 
+                              & (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      !
+      ! Conditionnal horizontal Laplacian diffusion:
+      !---------------------------------------------
+      IF ( ll_lapdiff_cond ) THEN
+         !
+         zwdw_b(:,:,1) = 0._wp
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               DO jk=2,jpk
+                  zwdw_b(ji,jj,jk) =  zwdw_b(ji,jj,jk-1) + & 
+                                   & (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
+               END DO
+            END DO
+         END DO
+         !
+         DO jk = 2, jpkm1
+            zwu(:,:) = 0._wp
+            zwv(:,:) = 0._wp
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.                  
+                  zwu(ji,jj) =  umask(ji,jj,jk) * e2_e1u(ji,jj) &
+                             &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj  ,jk) )
+                  zwv(ji,jj) =  vmask(ji,jj,jk) * e1_e2v(ji,jj) & 
+                             &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji  ,jj+1,jk) )
+               END DO
+            END DO
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ztmp1 = ( zwu(ji-1,jj  ) - zwu(ji,jj) &
+                     &  +   zwv(ji  ,jj-1) - zwv(ji,jj) ) * r1_e1e2t(ji,jj)
+                  zh2 = MAX(abs(ztmp1)-zhlim*SQRT(r1_e1e2t(ji,jj)), 0._wp)
+                  ztmp = SIGN(zh2, ztmp1)
+                  zeu2 = zhdiff * e1e2t(ji,jj)*e1e2t(ji,jj)/(e1t(ji,jj)*e1t(ji,jj) + e2t(ji,jj)*e2t(ji,jj))
+                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp *  tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO         
+         !
+      ENDIF
+
+      ! Conditionnal horizontal Bilaplacian diffusion:
+      !-----------------------------------------------
+      IF ( ll_blpdiff_cond ) THEN
+         !
+         zwdw_b(:,:,1) = 0._wp
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               DO jk = 2,jpkm1
+                  zwdw_b(ji,jj,jk) =  zwdw_b(ji,jj,jk-1) + & 
+                                   & (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
+               END DO
+            END DO
+         END DO
+         !
+         DO jk = 2, jpkm1
+            zwu(:,:) = 0._wp
+            zwv(:,:) = 0._wp
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.                  
+                  zwu(ji,jj) =  umask(ji,jj,jk) * e2_e1u(ji,jj) &
+                             &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj  ,jk) )
+                  zwv(ji,jj) =  vmask(ji,jj,jk) * e1_e2v(ji,jj) & 
+                             &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji  ,jj+1,jk) )
+               END DO
+            END DO
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zwdw_b(ji,jj,jk) = -( (zwu(ji-1,jj  ) - zwu(ji,jj)) &
+                                &  +    (zwv(ji  ,jj-1) - zwv(ji,jj)) ) * r1_e1e2t(ji,jj)
+               END DO
+            END DO
+         END DO   
+         !
+         CALL lbc_lnk( zwdw_b(:,:,:), 'T', 1. )      
+         !
+         DO jk = 2, jpkm1
+            zwu(:,:) = 0._wp
+            zwv(:,:) = 0._wp
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.                  
+                  zwu(ji,jj) =  umask(ji,jj,jk) * e2_e1u(ji,jj) &
+                             &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj  ,jk) )
+                  zwv(ji,jj) =  vmask(ji,jj,jk) * e1_e2v(ji,jj) & 
+                             &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji  ,jj+1,jk) )
+               END DO
+            END DO
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ztmp1 = ( (zwu(ji-1,jj  ) - zwv(ji,jj)) &
+                     &  +   (zwv(ji  ,jj-1) - zwv(ji,jj)) ) * r1_e1e2t(ji,jj)
+                  zh2 = MAX(abs(ztmp1)-zhlim2*SQRT(r1_e1e2t(ji,jj))*r1_e1e2t(ji,jj), 0._wp)
+                  ztmp = SIGN(zh2, ztmp1)
+                  zeu2 = zhdiff2 * e1e2t(ji,jj)*e1e2t(ji,jj) / 16._wp
+                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp *  tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO   
+         !
+      ENDIF
+
+      ! Conditionnal vertical diffusion:
+      !---------------------------------
+      IF ( ll_zdiff_cond ) THEN
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.    
+                  ztmp  = -( (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1))*tmask(ji,jj,jk-1) & 
+                            -(tilde_e3t_b(ji,jj,jk  )+e3t_0(ji,jj,jk  ))*tmask(ji,jj,jk  ) ) 
+                  ztmp1 = 0.5_wp * ( tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1) & 
+                        &           +tilde_e3t_b(ji,jj,jk  ) + e3t_0(ji,jj,jk  ) )
+                  zh2   = MAX(abs(ztmp)-zvlim*ztmp1, 0._wp)
+                  ztmp  = SIGN(zh2, ztmp)
+                  IF ((jk==mbkt(ji,jj)).AND.(ln_zps)) ztmp=0.e0
+                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zvdiff * ztmp * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ENDIF
+      !
+      ! Check grid from the bottom to the surface
+      !------------------------------------------
+      IF ( ll_chk_bot2top ) THEN
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               jkbot = mbkt(ji,jj)                   
+               DO jk = jkbot,2,-1
+                  !
+                  zh_0   = e3t_0(ji,jj,jk)
+                  zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1))
+                  zh_old = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
+                  zh_min = MIN(zh_0/3._wp, dzmin_int)
+                  ! 
+                  ! Set maximum and minimum vertical excursions   
+                  ztmph = hsm(ji,jj)
+                  ztmpd = dsm(ji,jj)
+                  zh2   = ztmph * exp(-(gdepw_0(ji,jj,jk)-gdepw_0(ji,jj,i_int_bot(ji,jj)))/ztmpd)
+!                  zh2   = ztmph * exp(-(gdepw_0(ji,jj,jk)-gdepw_0(ji,jj,i_int_bot(ji,jj)+1))/ztmpd)
+                  zh2 = MAX(zh2,0.001_wp) ! Extend tolerance a bit for stability reasons (to be explored)
+                  zdiff = cush_max(gdepw_0(ji,jj,jk)-zwdw(ji,jj,jk), zh2 )
+                  zwdw(ji,jj,jk) = MAX(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) - zdiff)
+                  zdiff = cush_max(zwdw(ji,jj,jk)-gdepw_0(ji,jj,jk), zh2 )
+                  zwdw(ji,jj,jk) = MIN(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) + zdiff)
+                  !
+                  ! New layer thickness:
+                  zh_new =  zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
+                  !
+                  ! Ensure minimum layer thickness:                  
+!                  zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new)
+                  zh_new = cush(zh_new, zh_min)
+                  !
+                  ! Final flux:
+                  zdiff = (zh_new - zh_old) * tmask(ji,jj,jk)
+                  !
+                  ! Limit thickness change in 1 time step: 
+                  ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef )
+                  zdiff = SIGN(ztmp, zh_new - zh_old)
+                  zh_new = zdiff + zh_old
+                  !
+                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk+1) - zh_new       
+               END DO    
+            END DO 
+         END DO
+      END IF    
+      !
+      ! Check grid from the surface to the bottom 
+      !------------------------------------------ 
+      IF ( ll_chk_top2bot ) THEN      
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               jkbot = mbkt(ji,jj)   
+               DO jk = 1, jkbot-1
+                  !
+                  zh_0   = e3t_0(ji,jj,jk)
+                  zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk+1) + e3t_0(ji,jj,jk+1))
+                  zh_old = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
+                  zh_min = MIN(zh_0/3._wp, dzmin_int)
+                  !
+                  zwdw(ji,jj,jk+1) = MAX(zwdw(ji,jj,jk+1), REAL(jk)*dzmin_surf)
+                  !
+                  ! New layer thickness:
+                  zh_new =  zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
+                  !
+                  ! Ensure minimum layer thickness:
+!                  zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new)
+                  zh_new = cush(zh_new, zh_min)
+                  !
+                  ! Final flux:
+                  zdiff = (zh_new -zh_old) * tmask(ji,jj,jk)
+                  ! 
+                  ! Limit flux:                 
+!                  ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef )
+!                  zdiff = SIGN(ztmp, zh_new - zh_old)
+                  zh_new = zdiff + zh_old
+                  !
+                  zwdw(ji,jj,jk+1) = zwdw(ji,jj,jk) + zh_new
+               END DO
+               !               
+            END DO
+         END DO
+      ENDIF
+      !
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            DO jk = 1, jpkm1
+               tilde_e3t_a(ji,jj,jk) =  (zwdw(ji,jj,jk+1)-zwdw(ji,jj,jk)-e3t_0(ji,jj,jk)) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      ! 
+      !
+      IF( ln_timing )  CALL timing_stop('dom_vvl_regrid')
+      !
+   END SUBROUTINE dom_vvl_regrid
+
+   FUNCTION cush(hin, hmin)  RESULT(hout)
+      !!----------------------------------------------------------------------
+      !!                 ***  FUNCTION cush  ***
+      !!
+      !! ** Purpose :   
+      !!
+      !! ** Method  :
+      !!
+      !!----------------------------------------------------------------------
+      IMPLICIT NONE
+      REAL(wp), INTENT(in) ::  hin, hmin
+      REAL(wp)             ::  hout, zx, zh_cri
+      !!----------------------------------------------------------------------
+      zh_cri = 3._wp * hmin
+      !
+      IF ( hin<=0._wp ) THEN
+         hout = hmin
+      !
+      ELSEIF ( (hin>0._wp).AND.(hin<=zh_cri) ) THEN
+         zx = hin/zh_cri
+         hout = hmin * (1._wp + zx + zx*zx)      
+      !
+      ELSEIF ( hin>zh_cri ) THEN
+         hout = hin
+      !
+      ENDIF
+      !
+   END FUNCTION cush
+
+   FUNCTION cush_max(hin, hmax)  RESULT(hout)
+      !!----------------------------------------------------------------------
+      !!                 ***  FUNCTION cush  ***
+      !!
+      !! ** Purpose :   
+      !!
+      !! ** Method  :
+      !!
+      !!----------------------------------------------------------------------
+      IMPLICIT NONE
+      REAL(wp), INTENT(in) ::  hin, hmax
+      REAL(wp)             ::  hout, hmin, zx, zh_cri
+      !!----------------------------------------------------------------------
+      hmin = 0.1_wp * hmax
+      zh_cri = 3._wp * hmin
+      !
+      IF ( (hin>=(hmax-zh_cri)).AND.(hin<=(hmax-hmin))) THEN
+         zx = (hmax-hin)/zh_cri
+         hout = hmax - hmin * (1._wp + zx + zx*zx)
+         !
+      ELSEIF ( hin>(hmax-zh_cri) ) THEN
+         hout = hmax - hmin
+         !
+      ELSE
+         hout = hin
+         !
+      ENDIF
+      !
+   END FUNCTION cush_max
+
+   SUBROUTINE dom_vvl_adv_fct( kt, pta, uin, vin )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_vvl_adv_fct  ***
+      !! 
+      !! **  Purpose :  Do thickness advection
+      !!
+      !! **  Method  :  FCT scheme to ensure positivity 
+      !!
+      !! **  Action  : - Update pta thickness tendency and diffusive fluxes
+      !!               - this is the total trend, hence it does include sea level motions
+      !!----------------------------------------------------------------------
+      !
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta    ! thickness baroclinic trend 
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   uin, vin  ! input velocities
+      !
+      INTEGER  ::   ji, jj, jk, ib, ib_bdy               ! dummy loop indices  
+      INTEGER  ::   ikbu, ikbv, ibot
+      REAL(wp) ::   z2dtt, zbtr, ztra        ! local scalar
+      REAL(wp) ::   zdi, zdj, zmin           !   -      -
+      REAL(wp) ::   zfp_ui, zfp_vj           !   -      -
+      REAL(wp) ::   zfm_ui, zfm_vj           !   -      -
+      REAL(wp) ::   zfp_hi, zfp_hj           !   -      -
+      REAL(wp) ::   zfm_hi, zfm_hj           !   -      -
+      REAL(wp) ::   ztout,  ztin, zfac       !   -      -
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zwy, zwi
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )  CALL timing_start('dom_vvl_adv_fct')
+      !
+      !
+      ! 1. Initializations
+      ! ------------------
+      !
+      IF( neuler == 0 .AND. kt == nit000 ) THEN
+         z2dtt =  rdt
+      ELSE
+         z2dtt = 2.0_wp * rdt
+      ENDIF
+      !
+      zwi(:,:,:) = 0.e0
+      zwx(:,:,:) = 0.e0 
+      zwy(:,:,:) = 0.e0
+      !
+      !
+      ! 2. upstream advection with initial mass fluxes & intermediate update
+      ! --------------------------------------------------------------------
+      IF ( ll_shorizd ) THEN
+         DO jk = 1, jpkm1
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  !
+                  zfp_hi = MAX(hu_b(ji,jj) - gdepw_b(ji  ,jj  ,jk), 0._wp)
+                  zfp_hi = MIN(zfp_hi, e3t_b(ji  ,jj  ,jk))
+                  zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) ) 
+                  !
+                  zfm_hi = MAX(hu_b(ji,jj) - gdepw_b(ji+1,jj  ,jk), 0._wp)
+                  zfm_hi = MIN(zfm_hi, e3t_b(ji+1,jj  ,jk))
+                  zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) ) 
+                  ! 
+                  zfp_hj = MAX(hv_b(ji,jj) - gdepw_b(ji  ,jj  ,jk), 0._wp)
+                  zfp_hj = MIN(zfp_hj, e3t_b(ji  ,jj  ,jk))
+                  zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) ) 
+                  !
+                  zfm_hj = MAX(hv_b(ji,jj) - gdepw_b(ji  ,jj+1,jk), 0._wp)
+                  zfm_hj = MIN(zfm_hj, e3t_b(ji  ,jj+1,jk))
+                  zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) )
+                  !
+                  zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) )
+                  zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) )
+                  zfp_vj = vin(ji,jj,jk) + ABS( vin(ji,jj,jk) )
+                  zfm_vj = vin(ji,jj,jk) - ABS( vin(ji,jj,jk) )
+                  zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
+                  zwy(ji,jj,jk) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ELSE
+         DO jk = 1, jpkm1
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  !
+                  zfp_hi = e3t_b(ji  ,jj  ,jk)
+                  zfm_hi = e3t_b(ji+1,jj  ,jk)             
+                  zfp_hj = e3t_b(ji  ,jj  ,jk)
+                  zfm_hj = e3t_b(ji  ,jj+1,jk)
+                  !
+                  zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) )
+                  zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) )
+                  zfp_vj = vin(ji,jj,jk) + ABS( vin(ji,jj,jk) )
+                  zfm_vj = vin(ji,jj,jk) - ABS( vin(ji,jj,jk) )
+                  zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
+                  zwy(ji,jj,jk) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ENDIF
+
+      ! total advective trend
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zbtr = r1_e1e2t(ji,jj)
+               ! total intermediate advective trends
+               ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )  &
+                  &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )  )
+               !
+               ! update and guess with monotonic sheme
+               pta(ji,jj,jk) =   pta(ji,jj,jk) + ztra
+               zwi(ji,jj,jk) = (e3t_b(ji,jj,jk) + z2dtt * ztra ) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      CALL lbc_lnk( zwi, 'T', 1. )  
+
+      IF ( ln_bdy ) THEN
+         DO ib_bdy=1, nb_bdy
+            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1)
+               ji = idx_bdy(ib_bdy)%nbi(ib,1)
+               jj = idx_bdy(ib_bdy)%nbj(ib,1)
+               DO jk = 1, jpkm1  
+                  zwi(ji,jj,jk) = e3t_a(ji,jj,jk)
+               END DO
+            END DO 
+         END DO
+      ENDIF
+
+      IF ( ln_vvl_dbg ) THEN
+         zmin = MINVAL( zwi(:,:,:), mask = tmask(:,:,:) == 1.e0 ) 
+         IF( lk_mpp )   CALL mpp_min( zmin )
+         IF( zmin < 0._wp) THEN
+            IF(lwp) CALL ctl_warn('vvl_adv: CFL issue here')
+            IF(lwp) WRITE(numout,*) zmin
+         ENDIF
+      ENDIF
+
+      ! 3. antidiffusive flux : high order minus low order
+      ! --------------------------------------------------
+      ! antidiffusive flux on i and j
+      DO jk = 1, jpkm1
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zwx(ji,jj,jk) = (e2u(ji,jj) * uin(ji,jj,jk) * e3u_n(ji,jj,jk) &
+                                & - zwx(ji,jj,jk)) * umask(ji,jj,jk)
+               zwy(ji,jj,jk) = (e1v(ji,jj) * vin(ji,jj,jk) * e3v_n(ji,jj,jk) & 
+                                & - zwy(ji,jj,jk)) * vmask(ji,jj,jk)
+               !
+               ! Update advective fluxes
+               un_td(ji,jj,jk) = un_td(ji,jj,jk) - zwx(ji,jj,jk)
+               vn_td(ji,jj,jk) = vn_td(ji,jj,jk) - zwy(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      
+      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral boundary conditions
+
+      ! 4. monotonicity algorithm
+      ! -------------------------
+      CALL nonosc_2d( e3t_b(:,:,:), zwx, zwy, zwi, z2dtt )
+
+      ! 5. final trend with corrected fluxes
+      ! ------------------------------------
+      !
+      ! Update advective fluxes
+      un_td(:,:,:) = (un_td(:,:,:) + zwx(:,:,:))*umask(:,:,:)
+      vn_td(:,:,:) = (vn_td(:,:,:) + zwy(:,:,:))*vmask(:,:,:)
+      !
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.  
+               !
+               zbtr = r1_e1e2t(ji,jj)
+               ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                  &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )
+               ! add them to the general tracer trends
+               pta(ji,jj,jk) = pta(ji,jj,jk) + ztra
+            END DO
+         END DO
+      END DO
+      !
+      IF( ln_timing )  CALL timing_stop('dom_vvl_adv_fct')
+      !
+   END SUBROUTINE dom_vvl_adv_fct
+
+   SUBROUTINE dom_vvl_ups_cor( kt, pta, uin, vin ) 
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_vvl_adv_fct  ***
+      !! 
+      !! **  Purpose :  Correct for addionnal barotropic fluxes 
+      !!                in the upstream direction
+      !!
+      !! **  Method  :   
+      !!
+      !! **  Action  : - Update diffusive fluxes uin, vin
+      !!               - Remove divergence from thickness tendency
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta       ! thickness baroclinic trend 
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   uin, vin  ! input fluxes
+      INTEGER  ::   ji, jj, jk               ! dummy loop indices  
+      INTEGER  ::   ikbu, ikbv, ibot
+      REAL(wp) ::   zbtr, ztra               ! local scalar
+      REAL(wp) ::   zdi, zdj                 !   -      -
+      REAL(wp) ::   zfp_hi, zfp_hj           !   -      -
+      REAL(wp) ::   zfm_hi, zfm_hj           !   -      -
+      REAL(wp) ::   zfp_ui, zfp_vj           !   -      -
+      REAL(wp) ::   zfm_ui, zfm_vj           !   -      -
+      REAL(wp), DIMENSION(jpi,jpj) :: zbu, zbv, zhu_b, zhv_b
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zwy
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )  CALL timing_start('dom_vvl_ups_cor')
+      !
+      ! Compute barotropic flux difference:
+      zbu(:,:) = 0.e0
+      zbv(:,:) = 0.e0
+      DO jj = 1, jpj
+         DO ji = 1, jpi   ! vector opt.
+            DO jk = 1, jpkm1
+               zbu(ji,jj) = zbu(ji,jj) - uin(ji,jj,jk) * umask(ji,jj,jk)
+               zbv(ji,jj) = zbv(ji,jj) - vin(ji,jj,jk) * vmask(ji,jj,jk)
+            END DO
+         END DO
+      ENDDO
+
+      ! Compute upstream depths:
+      zhu_b(:,:) = 0.e0
+      zhv_b(:,:) = 0.e0
+
+      IF ( ll_shorizd ) THEN
+         ! Correct bottom value
+         ! considering "shelf horizon depth"
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1   ! vector opt.
+               zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj))
+               zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj))
+               DO jk=1, jpkm1
+                  zfp_hi = MAX(hu_b(ji,jj) - gdepw_b(ji  ,jj  ,jk), 0._wp)
+                  zfp_hi = MIN(zfp_hi, e3t_b(ji  ,jj  ,jk))
+                  zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) ) 
+                  !
+                  zfm_hi = MAX(hu_b(ji,jj) - gdepw_b(ji+1,jj  ,jk), 0._wp)
+                  zfm_hi = MIN(zfm_hi, e3t_b(ji+1,jj  ,jk))
+                  zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) )
+                  ! 
+                  zfp_hj = MAX(hv_b(ji,jj) - gdepw_b(ji  ,jj  ,jk), 0._wp)
+                  zfp_hj = MIN(zfp_hj, e3t_b(ji  ,jj  ,jk))
+                  zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) ) 
+                  !
+                  zfm_hj = MAX(hv_b(ji,jj) - gdepw_b(ji  ,jj+1,jk), 0._wp)
+                  zfm_hj = MIN(zfm_hj, e3t_b(ji  ,jj+1,jk))
+                  zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) )
+                  !
+                  zhu_b(ji,jj) = zhu_b(ji,jj) + (         zdi  * zfp_hi &
+                               &                 + (1._wp-zdi) * zfm_hi &
+                               &                ) * umask(ji,jj,jk)
+                  zhv_b(ji,jj) = zhv_b(ji,jj) + (         zdj  * zfp_hj &
+                               &                 + (1._wp-zdj) * zfm_hj &
+                               &                ) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ELSE
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1   ! vector opt.
+               zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj)) 
+               zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj))
+               DO jk = 1, jpkm1
+                  zfp_hi = e3t_b(ji  ,jj  ,jk)
+                  zfm_hi = e3t_b(ji+1,jj  ,jk)
+                  zfp_hj = e3t_b(ji  ,jj  ,jk)
+                  zfm_hj = e3t_b(ji  ,jj+1,jk)
+                  !
+                  zhu_b(ji,jj) = zhu_b(ji,jj) + (          zdi  * zfp_hi &
+                               &                  + (1._wp-zdi) * zfm_hi &
+                               &                ) * umask(ji,jj,jk)
+                  !
+                  zhv_b(ji,jj) = zhv_b(ji,jj) + (          zdj  * zfp_hj &
+                               &                  + (1._wp-zdj) * zfm_hj &
+                               &                ) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ENDIF
+
+      ! Corrective barotropic velocity (times hor. scale factor)
+      zbu(:,:) = zbu(:,:)/ (zhu_b(:,:)*umask(:,:,1)+1._wp-umask(:,:,1))
+      zbv(:,:) = zbv(:,:)/ (zhv_b(:,:)*vmask(:,:,1)+1._wp-vmask(:,:,1))
+
+      CALL lbc_lnk( zbu(:,:), 'U', -1. )
+      CALL lbc_lnk( zbv(:,:), 'V', -1. )
+      
+      ! Set corrective fluxes in upstream direction:
+      !
+      zwx(:,:,:) = 0.e0
+      zwy(:,:,:) = 0.e0
+      IF ( ll_shorizd ) THEN
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               ! upstream scheme
+               zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) )
+               zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) )
+               zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) )
+               zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) )
+               DO jk = 1, jpkm1
+                  zfp_hi = MAX(hu_b(ji,jj) - gdepw_b(ji  ,jj  ,jk), 0._wp)
+                  zfp_hi = MIN(e3t_b(ji  ,jj  ,jk), zfp_hi)
+                  zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) ) 
+                  !
+                  zfm_hi = MAX(hu_b(ji,jj) - gdepw_b(ji+1,jj  ,jk), 0._wp)
+                  zfm_hi = MIN(e3t_b(ji+1,jj  ,jk), zfm_hi)
+                  zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) ) 
+                  !
+                  zfp_hj = MAX(hv_b(ji,jj) - gdepw_b(ji  ,jj  ,jk), 0._wp)
+                  zfp_hj = MIN(e3t_b(ji  ,jj  ,jk), zfp_hj) 
+                  zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) ) 
+                  !
+                  zfm_hj = MAX(hv_b(ji,jj) - gdepw_b(ji  ,jj+1,jk), 0._wp)
+                  zfm_hj = MIN(e3t_b(ji  ,jj+1,jk), zfm_hj)
+                  zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) ) 
+                  !
+                  zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
+                  zwy(ji,jj,jk) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ELSE
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               ! upstream scheme
+               zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) )
+               zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) )
+               zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) )
+               zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) )
+               DO jk = 1, jpkm1
+                  zfp_hi = e3t_b(ji  ,jj  ,jk)
+                  zfm_hi = e3t_b(ji+1,jj  ,jk)
+                  zfp_hj = e3t_b(ji  ,jj  ,jk)
+                  zfm_hj = e3t_b(ji  ,jj+1,jk)
+                  !
+                  zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
+                  zwy(ji,jj,jk) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ENDIF
+      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral boundary conditions
+
+      uin(:,:,:) = uin(:,:,:) + zwx(:,:,:)
+      vin(:,:,:) = vin(:,:,:) + zwy(:,:,:)
+      !
+      ! Update trend with corrective fluxes:
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.  
+               !
+               zbtr = r1_e1e2t(ji,jj)
+               ! total advective trends
+               ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                  &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )
+               ! add them to the general tracer trends
+               pta(ji,jj,jk) = pta(ji,jj,jk) + ztra
+            END DO
+         END DO
+      END DO
+      !
+      IF( ln_timing )  CALL timing_stop('dom_vvl_ups_cor')
+      !
+   END SUBROUTINE dom_vvl_ups_cor
+
+   SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, p2dt )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE nonosc_2d  ***
+      !!     
+      !! **  Purpose :   compute monotonic thickness fluxes from the upstream 
+      !!       scheme and the before field by a nonoscillatory algorithm 
+      !!
+      !! **  Method  :   ... ???
+      !!       warning : pbef and paft must be masked, but the boundaries
+      !!       conditions on the fluxes are not necessary zalezak (1979)
+      !!       drange (1995) multi-dimensional forward-in-time and upstream-
+      !!       in-space based differencing for fluid
+      !!----------------------------------------------------------------------
+      !
+      !!----------------------------------------------------------------------
+      REAL(wp)                         , INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
+      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field
+      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) ::   paa, pbb        ! monotonic fluxes in the 3 directions
+      !
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt   ! local scalars
+      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
+      REAL(wp) ::   zupip1, zupim1, zupjp1, zupjm1, zupb, zupa
+      REAL(wp) ::   zdoip1, zdoim1, zdojp1, zdojm1, zdob, zdoa
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )  CALL timing_start('nonosc2')
+      !
+      zbig  = 1.e+40_wp
+      zrtrn = 1.e-15_wp
+      zbetup(:,:,jpk) = 0._wp   ;   zbetdo(:,:,jpk) = 0._wp
+
+
+      ! Search local extrema
+      ! --------------------
+      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
+      zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ),   &
+         &        paft * tmask - zbig * ( 1.e0 - tmask )  )
+      zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ),   &
+         &        paft * tmask + zbig * ( 1.e0 - tmask )  )
+
+      DO jk = 1, jpkm1
+         z2dtt = p2dt
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+
+               ! search maximum in neighbourhood
+               zup = MAX(  zbup(ji  ,jj  ,jk  ),   &
+                  &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   &
+                  &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  ))
+
+               ! search minimum in neighbourhood
+               zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   &
+                  &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   &
+                  &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  ))
+
+               ! positive part of the flux
+               zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   &
+                  & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   
+
+               ! negative part of the flux
+               zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   &
+                  & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )
+
+               ! up & down beta terms
+               zbt = e1t(ji,jj) * e2t(ji,jj) / z2dtt
+               zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
+               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt
+            END DO
+         END DO
+      END DO
+
+      CALL lbc_lnk( zbetup, 'T', 1. )   ;   CALL lbc_lnk( zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
+
+      ! 3. monotonic flux in the i & j direction (paa & pbb)
+      ! ----------------------------------------
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
+               zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
+               zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) )
+               paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu )
+
+               zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
+               zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
+               zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) )
+               pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( paa, 'U', -1. )   ;   CALL lbc_lnk( pbb, 'V', -1. )   ! lateral boundary condition (changed sign)
+      !
+      IF( ln_timing )  CALL timing_stop('nonosc2')
+      !
+   END SUBROUTINE nonosc_2d
+
    !!======================================================================
 END MODULE domvvl
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DYN/dynnxt.F90
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DYN/dynnxt.F90	(revision 10115)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DYN/dynnxt.F90	(revision 10116)
@@ -95,7 +95,8 @@
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      INTEGER  ::   ikt          ! local integers
+      INTEGER  ::   ikt, jkbot   ! local integers
       REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars
       REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
+      REAL(wp) ::   zhcri, zmsku, zwgtu, zmskv, zwgtv  !   -      -
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva 
@@ -136,4 +137,32 @@
             END DO  
          ENDIF
+      ENDIF
+
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 
+!         ! Duplicate baroclinic velocities from above in massless layers
+!         ! Correction to handle correct barotropic velocities right after
+         zhcri = 0.1_wp
+         DO ji = 2, jpim1
+            DO jj= 2, jpjm1
+               jkbot=jpkm1
+               DO jk=jpkm1,2,-1
+                  IF (e3u_a(ji,jj,jk)>zhcri) jkbot = jk+1
+               END DO
+               DO jk=jkbot,jpkm1
+                  zwgtu = MIN(zhcri, e3u_a(ji,jj,jk))
+                  zmsku = 0.5_wp * ( 1._wp + SIGN(1._wp, e3u_a(ji,jj,jk) - 0.01_wp) )
+                  ua(ji,jj,jk) = zmsku * (zwgtu * ua(ji,jj,jk) + (zhcri - zwgtu) * ua(ji,jj,jk-1)*zwgtu/(zwgtu+0.01_wp)) / zhcri
+               END DO
+               jkbot=jpkm1
+               DO jk=jpkm1,2,-1
+                  IF (e3v_a(ji,jj,jk)>zhcri) jkbot = jk+1
+               END DO
+               DO jk=jkbot,jpkm1
+                  zwgtv = MIN(zhcri, e3v_a(ji,jj,jk))
+                  zmskv = 0.5_wp * ( 1._wp + SIGN(1._wp, e3v_a(ji,jj,jk) - 0.01_wp) )
+                  va(ji,jj,jk) = zmskv * (zwgtv * va(ji,jj,jk) + (zhcri - zwgtv) * va(ji,jj,jk-1)*zwgtv/(zwgtv+0.01_wp)) / zhcri
+               END DO
+            END DO
+         END DO
       ENDIF
 
@@ -291,7 +320,9 @@
                         zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
                         zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
+                        zmsku = 0.5_wp * ( 1._wp + SIGN(1._wp, ze3u_f(ji,jj,jk) - 0.01_wp) )
+                        zmskv = 0.5_wp * ( 1._wp + SIGN(1._wp, ze3v_f(ji,jj,jk) - 0.01_wp) )
                         !
-                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
-                        vb(ji,jj,jk) = zvf
+                        ub(ji,jj,jk) = zmsku * zuf             ! ub <-- filtered velocity
+                        vb(ji,jj,jk) = zmskv * zvf
                         un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
                         vn(ji,jj,jk) = va(ji,jj,jk)
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/SBC/sbcwave.F90
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/SBC/sbcwave.F90	(revision 10115)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/SBC/sbcwave.F90	(revision 10116)
@@ -21,5 +21,4 @@
    USE zdf_oce,  ONLY : ln_zdfswm
    USE bdy_oce        ! open boundary condition variables
-   USE domvvl         ! domain: variable volume layers
    !
    USE iom            ! I/O manager library
