Index: /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 9544)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 9545)
@@ -98,4 +98,6 @@
          <field id="gdepwt"       long_name="T-cell interface depth"                                             unit="m"   grid_ref="grid_T_3D" />
          <field id="stiff_tilde"  long_name="Stiffness"                                                          unit=" "   grid_ref="grid_T_3D" />
+         <field id="gdepwu"       long_name="U-cell interface depth"                                             unit="m"   grid_ref="grid_U_3D" />
+         <field id="gdepwv"       long_name="V-cell interface depth"                                             unit="m"   grid_ref="grid_V_3D" />
          <field id="depw_tilde"   long_name="Interface displacement"                                             unit="m"   grid_ref="grid_W_3D" />
          <field id="dz_tilde"     long_name="Vertical grid deformation"                                          unit=" "   grid_ref="grid_W_3D" />
Index: /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 9544)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 9545)
@@ -833,4 +833,6 @@
    rn_zdef_max   = 0.9e0            !  maximum fractional e3t deformation
    ln_vvl_dbg    = .true.           !  debug prints    (T/F)
+   ln_vvl_ramp   = .false.          !  linear ramp at startup
+      rn_day_ramp = 3.e0                ! ramp duration [days]
 /
 !-----------------------------------------------------------------------
Index: /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/cfg.txt
===================================================================
--- /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/cfg.txt	(revision 9544)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/cfg.txt	(revision 9545)
@@ -12,5 +12,5 @@
 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC
 ORCA2LIM3_16 OPA_SRC LIM_SRC_3 NST_SRC
-ORCA2LIM3_LONG OPA_SRC LIM_SRC_3 NST_SRC
 RIDGE OPA_SRC
 COMODO_IW OPA_SRC
+ORCA2LIM3_LONG OPA_SRC LIM_SRC_3 NST_SRC
Index: /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 9544)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 9545)
@@ -158,9 +158,60 @@
       CALL iom_put("tpt_dep", fsdept_n(:,:,:) )
       IF( iom_use("gdepwt") )  THEN
-         z3d(:,:,1) = (gdepw_n(:,:,1)-sshn(:,:))*tmask(:,:,1)
+#if defined key_vvl
+         z3d(:,:,1) = (fsdepw(:,:,1)-sshn(:,:))*tmask(:,:,1)
          DO jk=2,jpk
-            z3d(:,:,jk) = (gdepw_n(:,:,jk)-sshn(:,:))*tmask(:,:,jk-1)
-         END DO
+            z3d(:,:,jk) = (fsdepw(:,:,jk)-sshn(:,:))*tmask(:,:,jk-1)
+         END DO
+#else
+         z3d(:,:,1) = 0._wp
+         DO jk=2,jpk
+            z3d(:,:,jk) = fsdepw(:,:,jk)*tmask(:,:,jk-1)
+         END DO
+#endif
          CALL iom_put("gdepwt" , z3d(:,:,:) )
+      END IF
+
+      IF( iom_use("gdepwu") )  THEN
+         z2d(:,:) = 0._wp
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               z2d(ji,jj) = 0.5_wp * umask(ji,jj,1) * r1_e12u(ji,jj)                             &
+                        &          * ( e12t(ji,jj) * sshn(ji,jj) + e12t(ji+1,jj) * sshn(ji+1,jj) )
+            END DO
+         END DO
+         CALL lbc_lnk( z2d(:,:), 'U', 1._wp )
+
+         z3d(:,:,1)=0._wp
+         DO jk=2,jpk
+            z3d(:,:,jk) = z3d(:,:,jk-1) + fse3u_n(:,:,jk-1) * umask(:,:,jk-1)
+         END DO
+
+         z3d(:,:,1) = (z3d(:,:,1)-z2d(:,:))*umask(:,:,1)
+         DO jk=2,jpk
+            z3d(:,:,jk) = (z3d(:,:,jk)-z2d(:,:))*umask(:,:,jk-1)
+         END DO
+         CALL iom_put("gdepwu" , z3d(:,:,:) )
+      END IF
+
+      IF( iom_use("gdepwv") )  THEN
+         z2d(:,:) = 0._wp
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               z2d(ji,jj) = 0.5_wp * vmask(ji,jj,1) * r1_e12v(ji,jj)                             &
+                        &          * ( e12t(ji,jj) * sshn(ji,jj) + e12t(ji,jj+1) * sshn(ji,jj+1) )
+            END DO
+         END DO
+         CALL lbc_lnk( z2d(:,:), 'V', 1._wp )
+
+         z3d(:,:,1)=0._wp
+         DO jk=2,jpk
+            z3d(:,:,jk) = z3d(:,:,jk-1) + fse3v_n(:,:,jk-1) * vmask(:,:,jk-1)
+         END DO
+
+         z3d(:,:,1) = (z3d(:,:,1)-z2d(:,:))*vmask(:,:,1)
+         DO jk=2,jpk
+            z3d(:,:,jk) = (z3d(:,:,jk)-z2d(:,:))*vmask(:,:,jk-1)
+         END DO
+         CALL iom_put("gdepwv" , z3d(:,:,:) )
       END IF
 
Index: /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 9544)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 9545)
@@ -35,4 +35,5 @@
    USE timing          ! Timing
    USE bdy_oce         ! ocean open boundary conditions
+   USE sbcrnf          ! river runoff 
 
    IMPLICIT NONE
@@ -57,4 +58,5 @@
    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                          ! Laplacian thickness diffusion
    LOGICAL                                               :: ln_vvl_blp                          ! Bilaplacian thickness diffusion 
@@ -66,4 +68,5 @@
    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)                                              :: rn_zdef_max               ! maximum fractional e3t deformation
    REAL(wp)                                              :: hsmall=0.01               ! small thickness
@@ -354,5 +357,5 @@
       IF( PRESENT(kcall) ) THEN
          ! comment line below if tilda coordinate has to be computed at each call
-         IF (kcall == 2 .AND. ln_vvl_ztilde.OR.ln_vvl_layer ) ll_do_bclinic = .FALSE.
+         IF (kcall == 2 .AND. (ln_vvl_ztilde.OR.ln_vvl_layer) ) ll_do_bclinic = .FALSE.
          IF (kcall == 2) ncall=2    
       ENDIF
@@ -388,5 +391,13 @@
             zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
          END DO
-         zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - tmask_i(:,:) )      
+         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) * fse3t_n(:,:,jk)
+            END DO   
+         ENDIF 
 
          ! Thickness advection:
@@ -394,4 +405,19 @@
          ! 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
+                  ztu(:,:,jk) = un(:,:,jk) * fse3u_n(:,:,jk) * e2u(:,:)
+                  ztv(:,:,jk) = vn(:,:,jk) * fse3v_n(:,:,jk) * e1v(:,:)
+                  ze3t(:,:,jk) = -tilde_e3t_a(:,:,jk) - fse3t_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(:,:,:)
+!               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
@@ -561,5 +587,5 @@
          ! ~~~~~~~~~~~~~~~~~~~~~~
          zramp = 1._wp
-!         IF (.NOT.ln_rstart) zramp = MIN(MAX( ((kt-nit000)*rdt)/(3._wp*rday),0._wp),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
@@ -626,21 +652,21 @@
       ENDIF
 
-      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
-               ztu(:,:,jk) = un(:,:,jk) * fse3u_n(:,:,jk) * e2u(:,:) + un_td(:,:,jk)
-               ztv(:,:,jk) = vn(:,:,jk) * fse3v_n(:,:,jk) * e1v(:,:) + vn_td(:,:,jk)
-               ze3t(:,:,jk) = -fse3t_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(:,:,:)
-!               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
-      ENDIF
+!      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
+!               ztu(:,:,jk) = un(:,:,jk) * fse3u_n(:,:,jk) * e2u(:,:) + un_td(:,:,jk)
+!               ztv(:,:,jk) = vn(:,:,jk) * fse3v_n(:,:,jk) * e1v(:,:) + vn_td(:,:,jk)
+!               ze3t(:,:,jk) = -fse3t_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(:,:,:)
+!!               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
+!      ENDIF
 
       IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
@@ -675,5 +701,5 @@
 
 
-      IF( ln_vvl_dbg .AND. ( ncall==2 ) ) THEN   ! - ML - test: control prints for debuging
+      IF( (ln_vvl_dbg .AND. ( ncall==2 )).AND.(ln_vvl_ztilde .OR. ln_vvl_layer) ) THEN   ! - ML - test: control prints for debuging
          !
          zht(:,:) = 0.0_wp
@@ -1131,5 +1157,5 @@
       INTEGER ::   jkbot                                               ! bottom level
       LOGICAL ::   l_is_orca                                           ! local logical
-      REAL(wp) :: zmin, zdo, zup, ztap, zsmall
+      REAL(wp) :: zmin, zdo, zup, ztap
       REAL(wp), POINTER, DIMENSION(:,:)   :: zs                        ! Surface interface depth anomaly
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zw                        ! Interface depth anomaly
@@ -1141,10 +1167,9 @@
 
  
+      nmet=1 
 !      nmet=0 ! Original method (Surely wrong)
-!      nmet= 1 ! Interface interpolation
-      nmet=1 ! Internal interfaces interpolation only, spread barotropic increment
-
-      ztap = 0._wp ! Minimum fraction of T-point thickness at cell interfaces
-      zsmall = 1e-3_wp
+!      nmet=1 ! Interface interpolation
+!      nmet=2 ! Internal interfaces interpolation only, spread barotropic increment
+      ztap = 0.1_wp ! Minimum fraction of T-point thickness at cell interfaces
 
       IF ( (nmet==1).OR.(nmet==2) ) THEN
@@ -1158,31 +1183,25 @@
             zw(:,:,:) =  0._wp    
             !
-!            DO jj = 1, jpj
-!               DO ji = 1, jpi
-!                  jkbot = mbkt(ji,jj)
-!                  DO jk=jkbot,1,-1
-!                     zw(ji,jj,jk) = zw(ji,jj,jk+1) - pe3_in(ji,jj,jk) + e3t_0(ji,jj,jk)
-!                  END DO
-!               END DO
-!            END DO 
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  jkbot = mbkt(ji,jj)
+                  DO jk=jkbot,1,-1
+                     zw(ji,jj,jk) = zw(ji,jj,jk+1) - pe3_in(ji,jj,jk) + e3t_0(ji,jj,jk)
+                  END DO
+               END DO
+            END DO 
             !
-            DO jk=2,jpk
-               zw(:,:,jk) = zw(:,:,jk-1) + pe3_in(:,:,jk-1)*tmask(:,:,jk-1)
-            END DO 
-            ! Interface depth anomalies:
-            DO jk=1,jpkm1
-               zw(:,:,jk) = zw(:,:,jk) - zw(:,:,jpk) + ht_0(:,:)
-            END DO
-            zw(:,:,jpk) = ht_0(:,:)
             !
             IF (nmet==2) THEN        ! Consider "internal" interfaces only
                zs(:,:) = - zw(:,:,1) ! Save surface anomaly (ssh)
                !
+               zw(:,:,:) = 0._wp
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     DO jk=1,jpk
-                        zw(ji,jj,jk) = (zw(ji,jj,jk) + zs(ji,jj))                                          &
+                     jkbot = mbkt(ji,jj)
+                     DO jk=jkbot,1,-1
+                        zw(ji,jj,jk) = zw(ji,jj,jk+1) - ( pe3_in(ji,jj,jk)                     &
                                      & * ht_0(ji,jj) / (ht_0(ji,jj) + zs(ji,jj) + 1._wp - tmask(ji,jj,1))  &
-                                     & * tmask(ji,jj,jk)
+                                     & - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk)
                      END DO
                   END DO
@@ -1217,24 +1236,28 @@
                   ! Correction at last level:
                   jkbot = mbku(ji,jj)
-                  zdo   = hu_0(ji,jj)
-                  DO jk=jkbot,1,-1
-                     zup = 0.5_wp * umask(ji,jj,jk)   * r1_e12u(ji,jj)   &
-                           &      * (   e12t(ji  ,jj) * zw(ji  ,jj,jk)   &
-                           &          + e12t(ji+1,jj) * zw(ji+1,jj,jk)   )
-                     !
-                     ! 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 * pe3_in(ji+1,jj,jk)
-                           zmin = ztap * (zw(ji+1,jj,jk+1)-zw(ji+1,jj,jk))
-                        ELSE
-!                           zmin = ztap * pe3_in(ji  ,jj,jk)
-                           zmin = ztap * (zw(ji  ,jj,jk+1)-zw(ji  ,jj,jk))
-                        ENDIF
-                        zup = MIN(zup, zdo-zmin)
+                  pe3_out(ji,jj,jkbot) = -0.5_wp * umask(ji,jj,jkbot) * r1_e12u(ji,jj) &
+                        &                    * (   e12t(ji  ,jj) * zw(ji  ,jj,jkbot)   &
+                        &                        + e12t(ji+1,jj) * zw(ji+1,jj,jkbot)   )
+                  !
+                  ! 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) ) ) THEN
+                     IF ( ht_0(ji+1,jj) < ht_0(ji,jj) ) THEN
+!                        zmin = ztap * pe3_in(ji+1,jj,jkbot)
+                        zmin = ztap * (-zw(ji+1,jj,jkbot)+e3t_0(ji+1,jj,jkbot))
+                     ELSE
+!                        zmin = ztap * pe3_in(ji  ,jj,jkbot)
+                        zmin = ztap * (-zw(ji,jj,jkbot)+e3t_0(ji,jj,jkbot))
                      ENDIF
-                     zup = MIN(zup, zdo-zsmall)
-                     pe3_out(ji,jj,jk) = zdo - zup - e3u_0(ji,jj,jk)
-                     zdo = zup
+                     zmin = -e3u_0(ji,jj,jkbot) + zmin
+                     pe3_out(ji,jj,jkbot) = MAX(pe3_out(ji,jj,jkbot), zmin)
+                  ENDIF
+                  !
+                  zdo =  -pe3_out(ji,jj,jkbot)
+                  DO jk=jkbot-1,1,-1
+                     zup = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji  ,jj) &
+                              &                     *(   e12t(ji  ,jj) * zw(ji  ,jj,jk)  &
+                              &                         +e12t(ji+1,jj) * zw(ji+1,jj,jk)  )
+                     pe3_out(ji,jj,jk) = zdo - zup  
+                     zdo =  zdo - pe3_out(ji,jj,jk)
                   END DO
                END DO
@@ -1286,24 +1309,29 @@
                   ! Correction at last level:
                   jkbot = mbkv(ji,jj)
-                  zdo   = hv_0(ji,jj)
-                  DO jk=jkbot,1,-1
-                     zup = 0.5_wp * vmask(ji,jj,jk)   * r1_e12v(ji,jj)  &
-                           &      * (   e12t(ji,jj  ) * zw(ji,jj  ,jk)   &
-                           &          + e12t(ji,jj+1) * zw(ji,jj+1,jk)   )
+                  pe3_out(ji,jj,jkbot) = -0.5_wp * vmask(ji,jj,jkbot) * r1_e12v(ji,jj) &
+                        &                    * (   e12t(ji,jj  ) * zw(ji,jj  ,jkbot)   &
+                        &                        + e12t(ji,jj+1) * zw(ji,jj+1,jkbot)   )
+                  !
+                  ! 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) ) ) THEN
+                     IF ( ht_0(ji,jj+1) < ht_0(ji,jj) ) THEN
+!                        zmin = ztap * pe3_in(ji,jj+1,jkbot)
+                        zmin = ztap * (-zw(ji,jj+1,jkbot)+e3t_0(ji,jj+1,jkbot))
+                     ELSE
+!                        zmin = ztap * pe3_in(ji,jj  ,jkbot)
+                        zmin = ztap * (-zw(ji,jj,jkbot)+e3t_0(ji,jj,jkbot))
+                     ENDIF
+                     zmin = -e3v_0(ji,jj,jkbot) + zmin
+                     pe3_out(ji,jj,jkbot) = MAX(pe3_out(ji,jj,jkbot), zmin)
+                  ENDIF
+                  !
+                  zdo =  -pe3_out(ji,jj,jkbot)
+                  DO jk=jkbot-1,1,-1
+                     zup = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj  ) &
+                              &                     *  ( e12t(ji,jj  ) * zw(ji,jj  ,jk) &
+                              &                         +e12t(ji,jj+1) * zw(ji,jj+1,jk) )
                      !
-                     ! 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 * pe3_in(ji,jj+1,jk)
-                           zmin = ztap * (zw(ji,jj+1,jk+1)-zw(ji,jj+1,jk))
-                        ELSE
-!                           zmin = ztap * pe3_in(ji  ,jj,jk)
-                           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)
-                     zdo = zup
+                     pe3_out(ji,jj,jk) = zdo - zup 
+                     zdo =  zdo - pe3_out(ji,jj,jk)
                   END DO
                END DO
@@ -1358,6 +1386,35 @@
                   ! bottom correction:
                   jkbot = MIN(mbku(ji,jj), mbku(ji,jj+1)) 
-                  zdo = hf_0(ji,jj)
-                  DO jk=jkbot,1,-1
+                  pe3_out(ji,jj,jkbot) = -0.25_wp * umask(ji,jj,jkbot) * umask(ji,jj+1,jkbot) * r1_e12f(ji,jj) &
+                        &                         * (  e12t(ji  ,jj  ) * zw(ji  ,jj  ,jkbot)   &
+                        &                            + e12t(ji+1,jj  ) * zw(ji+1,jj  ,jkbot)   &
+                        &                            + e12t(ji  ,jj+1) * zw(ji  ,jj+1,jkbot)   &
+                        &                            + e12t(ji+1,jj+1) * zw(ji+1,jj+1,jkbot)   )
+                  !
+                  ! 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) ) ) 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 * pe3_in(ji+1,jj+1,jkbot)
+                           zmin = ztap * (-zw(ji+1,jj+1,jkbot)+e3t_0(ji+1,jj+1,jkbot))
+                        ELSE
+!                           zmin = ztap * pe3_in(ji  ,jj+1,jkbot)
+                           zmin = ztap * (-zw(ji,jj+1,jkbot)+e3t_0(ji,jj+1,jkbot))
+                        ENDIF
+                     ELSE
+                        IF ( ht_0(ji+1,jj  ) < ht_0(ji  ,jj  ) ) THEN
+!                           zmin = ztap * pe3_in(ji+1,jj  ,jkbot)
+                           zmin = ztap * (-zw(ji+1,jj,jkbot)+e3t_0(ji+1,jj,jkbot))
+                        ELSE
+!                           zmin = ztap * pe3_in(ji  ,jj  ,jkbot)
+                           zmin = ztap * (-zw(ji,jj,jkbot)+e3t_0(ji,jj,jkbot))
+                        ENDIF
+                     ENDIF
+                     zmin = -e3f_0(ji,jj,jkbot) + zmin
+                     pe3_out(ji,jj,jkbot) = MAX(pe3_out(ji,jj,jkbot), zmin)
+                  ENDIF
+                  !
+                  zdo =  -pe3_out(ji,jj,jkbot)
+                  DO jk=jkbot-1,1,-1
                      zup =  0.25_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) &
                            &        * (  e12t(ji  ,jj  ) * zw(ji  ,jj  ,jk)  & 
@@ -1365,33 +1422,11 @@
                            &           + e12t(ji  ,jj+1) * zw(ji  ,jj+1,jk)  & 
                            &           + e12t(ji+1,jj+1) * zw(ji+1,jj+1,jk)  )
+                     pe3_out(ji,jj,jk) = zdo - zup 
                      !
-                     ! 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 * pe3_in(ji+1,jj+1,jk)
-                              zmin = ztap * (zw(ji+1,jj+1,jk+1)-zw(ji+1,jj+1,jk))
-                           ELSE
-!                              zmin = ztap * pe3_in(ji  ,jj+1,jk)
-                              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 * pe3_in(ji+1,jj  ,jk)
-                              zmin = ztap * (zw(ji+1,jj  ,jk+1)-zw(ji+1,jj  ,jk))
-                           ELSE
-!                              zmin = ztap * pe3_in(ji  ,jj  ,jk)
-                              zmin = ztap * (zw(ji  ,jj  ,jk+1)-zw(ji  ,jj  ,jk))
-                           ENDIF
-                        ENDIF
-                        zup = MIN(zup, zdo-zmin)
-                     ENDIF
-                     zup = MIN(zup, zdo-zsmall)
-                     !
-                     pe3_out(ji,jj,jk) = zdo - zup - e3f_0(ji,jj,jk)
-                     zdo = zup
+                     zdo =  zdo - pe3_out(ji,jj,jk)                                    
                   END DO
-               END DO
-            END DO
+                  !
+               END DO
+            END DO 
          ENDIF
          !
@@ -1456,4 +1491,5 @@
 
    END SUBROUTINE dom_vvl_interpol
+
 
    SUBROUTINE dom_vvl_rst( kt, cdrw )
@@ -1626,4 +1662,5 @@
                       & rn_rst_e3t                 , rn_lf_cutoff                        , &
                       & ln_vvl_regrid              , rn_zdef_max                         , &
+                      & ln_vvl_ramp                , rn_day_ramp                         , &
                       & ln_vvl_dbg   ! not yet implemented: ln_vvl_kepe
       !!---------------------------------------------------------------------- 
@@ -1660,5 +1697,8 @@
          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,*) '              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
          WRITE(numout,*) '           Namelist nam_vvl : maximum e3t deformation fractional change'
          WRITE(numout,*) '                                         rn_zdef_max    = ', rn_zdef_max
@@ -2401,5 +2441,5 @@
       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.1_wp   ! Maximum fractionnal thickness change in one time step (<= 1.)
+      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
@@ -2621,7 +2661,6 @@
                   !
                   zh_0   = e3t_0(ji,jj,jk)
-                  zh_bef = tilde_e3t_b(ji,jj,jk) + zh_0
+                  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_dwn = tilde_e3t_a(ji,jj,jk-1) + e3t_0(ji,jj,jk-1)
                   zh_min = MIN(zh_0/3._wp, dzmin_int)
                   ! 
@@ -2639,6 +2678,5 @@
                   !
                   ! Ensure minimum layer thickness:                  
-!                  zh_new = MAX(zh_new, zh_dwn * zfrch_rel / (2._wp-zfrch_rel) )
-                  zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new)
+!                  zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new)
                   zh_new = cush(zh_new, zh_min)
                   !
@@ -2647,11 +2685,9 @@
                   !
                   ! Limit thickness change in 1 time step: 
-!                  zh_new = MIN( ABS(zh_new-zh_bef), (1._wp-zfrch_stp)*zh_bef )
-!                  zdiff = SIGN(ztmp, zh_new - zh_old)
-!                  zh_new = zdiff + zh_old
-                  !
-!                  tilde_e3t_a(ji,jj,jk  ) = (zh_new - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk)
-                  zwdw(ji,jj,jk)          = zwdw(ji,jj,jk+1) - zh_new
-!                  tilde_e3t_a(ji,jj,jk-1) = (-zdiff + tilde_e3t_a(ji,jj,jk-1) ) * tmask(ji,jj,jk-1)         
+                  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 
@@ -2668,7 +2704,6 @@
                   !
                   zh_0   = e3t_0(ji,jj,jk)
-                  zh_bef = tilde_e3t_b(ji,jj,jk) + zh_0
+                  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_up  = tilde_e3t_a(ji,jj,jk+1) + e3t_0(ji,jj,jk+1)
                   zh_min = MIN(zh_0/3._wp, dzmin_int)
                   !
@@ -2679,6 +2714,5 @@
                   !
                   ! Ensure minimum layer thickness:
-!                  zh_new=MAX(zh_new, zh_up * zfrch_rel / (2._wp-zfrch_rel) )
-                  zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new)
+!                  zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new)
                   zh_new = cush(zh_new, zh_min)
                   !
@@ -2687,11 +2721,9 @@
                   ! 
                   ! Limit flux:                 
-!                  ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef )
-!                  zdiff = SIGN(ztmp, zdiff)
-!                  zh_new = zdiff + zh_old
-                  !
-!                  tilde_e3t_a(ji,jj,jk  ) = (zh_new - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk)
-                  zwdw(ji,jj,jk+1)        = zwdw(ji,jj,jk) + zh_new
-!                  tilde_e3t_a(ji,jj,jk+1) = (-zdiff + tilde_e3t_a(ji,jj,jk+1) ) * tmask(ji,jj,jk+1)
+                  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
                !               
