Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/cfgs/SHARED/field_def_nemo-oce.xml	(revision 10163)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/cfgs/SHARED/field_def_nemo-oce.xml	(revision 10164)
@@ -21,4 +21,7 @@
         <field id="e3t"          long_name="T-cell thickness"                    standard_name="cell_thickness"     unit="m"   grid_ref="grid_T_3D" />
         <field id="e3t_surf"     long_name="T-cell thickness"   field_ref="e3t"  standard_name="cell_thickness"     unit="m"   grid_ref="grid_T_SFC"/>
+        <field id="e3t_0"        long_name="Initial T-cell thickness"            standard_name="ref_cell_thickness" unit="m"   grid_ref="grid_T_3D" />
+        <field id="e3t_star"     long_name="Barotropic T-cell thickness anomaly" standard_name="barotropic_cell_thickness_anomaly" unit="m"   grid_ref="grid_T_3D" />
+        <field id="e3t_tilde"    long_name="Baroclinic T-cell thickness anomaly" standard_name="baroclinic_cell_thickness_anomaly" unit="m"   grid_ref="grid_T_3D" />
         <field id="e3t_0"        long_name="Initial T-cell thickness"            standard_name="ref_cell_thickness" unit="m"   grid_ref="grid_T_3D" />
 
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 10163)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DIA/diawri.F90	(revision 10164)
@@ -58,4 +58,5 @@
    USE diurnal_bulk    ! diurnal warm layer
    USE cool_skin       ! Cool skin
+   USE domvvl
 
    IMPLICIT NONE
@@ -124,5 +125,5 @@
 
       ! Output of initial vertical scale factor
-      CALL iom_put("e3t_0", e3t_0(:,:,:) )
+      CALL iom_put("e3t_0", e3t_0(:,:,:))
       CALL iom_put("e3u_0", e3u_0(:,:,:) )
       CALL iom_put("e3v_0", e3v_0(:,:,:) )
@@ -132,4 +133,15 @@
       CALL iom_put( "e3v" , e3v_n(:,:,:) )
       CALL iom_put( "e3w" , e3w_n(:,:,:) )
+      !
+      IF ( (ln_vvl_ztilde).OR.(ln_vvl_layer) ) THEN
+         IF( iom_use("e3t_star") ) THEN    ! Barotropic cell thickness anomaly
+            z3d(:,:,:) = (e3t_n(:,:,:)-tilde_e3t_n(:,:,:)-e3t_0(:,:,:))*tmask(:,:,:) 
+            CALL iom_put( "e3t_star" , z3d(:,:,:) )
+         ENDIF
+         IF( iom_use("e3t_tilde") )  THEN  ! Baroclinic cell thickness anomaly
+            CALL iom_put( "e3t_tilde" , tilde_e3t_n(:,:,:) )
+         ENDIF
+      ENDIF
+
       IF( iom_use("e3tdef") )   &
          CALL iom_put( "e3tdef"  , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
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 10163)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DOM/domvvl.F90	(revision 10164)
@@ -61,4 +61,5 @@
    LOGICAL          :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer
    !                                                       ! conservation: not used yet
+   INTEGER          :: nn_filt_order=1
    REAL(wp)         :: rn_ahe3_lap               ! thickness diffusion coefficient (Laplacian)
    REAL(wp)         :: rn_ahe3_blp               ! thickness diffusion coefficient (Bilaplacian)
@@ -69,7 +70,7 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                ! thickness diffusion transport
-   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_lf, vn_lf, hdivn_lf      ! low frequency fluxes and 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), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a                 ! baroclinic scale factors
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: tildemask                   ! mask tilde tendency
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                 ! restoring period for scale factors
@@ -94,5 +95,5 @@
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
          ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
-            &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   &
+            &      un_td  (jpi,jpj,jpk)      , vn_td  (jpi,jpj,jpk)     ,                              &
             &      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 )
@@ -102,6 +103,6 @@
       ENDIF
       IF( ln_vvl_ztilde ) THEN
-         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 )
+         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdivn_lf(jpi,jpj,jpk,nn_filt_order),  &
+            &      un_lf(jpi,jpj,jpk,nn_filt_order), vn_lf(jpi,jpj,jpk,nn_filt_order), 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')
@@ -263,5 +264,5 @@
          !
          IF ( ln_vvl_zstar_on_shelf ) THEN
-            zhmin =  50._wp
+            zhmin = 50._wp
             zhmax = 100._wp
             DO jj = 1, jpj
@@ -284,4 +285,39 @@
          IF ( (ztmp*rdt) > 1._wp) CALL ctl_stop( 'dom_vvl_init: rn_lf_cuttoff is too small' )
          !
+      ENDIF
+
+      IF( ln_vvl_layer ) THEN
+         IF ( ln_vvl_zstar_on_shelf ) THEN
+            zhmin = 100._wp
+            zhmax = 150._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
+                  tildemask(ji,jj)   = MIN(tildemask(ji,jj), zwgt)
+               END DO
+            END DO
+         ENDIF
+         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  
+                  IF( ABS(gphit(ji,jj)) >= 6.) THEN
+                     ! values outside the equatorial band and transition zone (ztilde)
+              
+                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
+                     ! values inside the equatorial band (ztilde as zstar)
+                     tildemask(ji,jj) = 0._wp
+                  ELSE
+                     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
+         ENDIF
       ENDIF
       !
@@ -335,5 +371,5 @@
       !
       LOGICAL                ::   ll_do_bclinic         ! local logical
-      INTEGER                ::   ji, jj, jk            ! dummy loop indices
+      INTEGER                ::   ji, jj, jk, jo        ! dummy loop indices
       INTEGER                ::   ib, ib_bdy, ip, jp    !   "     "     "
       INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
@@ -341,5 +377,5 @@
       REAL(wp)               ::   z2dt  , z_tmin, z_tmax! local scalars                
       REAL(wp)               ::   zalpha, zwgt          ! temporary scalars
-      REAL(wp)               ::   zdu, zdv, zramp
+      REAL(wp)               ::   zdu, zdv, zramp, zmet
       REAL(wp)               ::   ztra, zbtr, ztout, ztin, zfac, zmsku, zmskv
       REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
@@ -357,4 +393,6 @@
       ENDIF
 
+      zmet = 1._wp
+
       ll_do_bclinic = .TRUE.
       ncall = 1
@@ -362,6 +400,5 @@
       IF( PRESENT(kcall) ) THEN
          ! 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.
+         IF (kcall == 2 .AND. (ln_vvl_ztilde.OR.ln_vvl_layer) ) ll_do_bclinic = .FALSE.
          ncall = kcall  
       ENDIF
@@ -386,6 +423,12 @@
       END DO
       !
+      IF ((ln_vvl_ztilde.OR.ln_vvl_layer).AND.(zmet==1._wp)) THEN
+         DO jk = 1, jpkm1
+            e3t_a(:,:,jk) = e3t_a(:,:,jk) - tilde_e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+         END DO  
+      ENDIF
+
       IF( (ln_vvl_ztilde.OR.ln_vvl_layer) .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
-         !                                                               ! ------baroclinic part------ !
+         !                                                             ! ------baroclinic part------ !
          tilde_e3t_a(:,:,:) = 0.0_wp  ! tilde_e3t_a used to store tendency terms
          un_td(:,:,:) = 0.0_wp        ! Transport corrections
@@ -415,25 +458,23 @@
                   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(:,:)
+                  ze3t(:,:,jk) = -tilde_e3t_a(:,:,jk) - (e3t_n(:,:,jk)-zmet*tilde_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(:,:,:)
+                  un_lf(:,:,:,nn_filt_order)  =    un_lf(:,:,:,nn_filt_order) * (1._wp - zalpha) + zalpha * ztu(:,:,:)
+                  vn_lf(:,:,:,nn_filt_order)  =    vn_lf(:,:,:,nn_filt_order) * (1._wp - zalpha) + zalpha * ztv(:,:,:)
+               hdivn_lf(:,:,:,nn_filt_order)  = hdivn_lf(:,:,:,nn_filt_order) * (1._wp - zalpha) + zalpha * ze3t(:,:,:)
+
+               DO jo = nn_filt_order-1,1,-1
+                     un_lf(:,:,:,jo)  =    un_lf(:,:,:,jo) * (1._wp - zalpha) + zalpha *    un_lf(:,:,:,jo+1)
+                     vn_lf(:,:,:,jo)  =    vn_lf(:,:,:,jo) * (1._wp - zalpha) + zalpha *    vn_lf(:,:,:,jo+1)
+                  hdivn_lf(:,:,:,jo)  = hdivn_lf(:,:,:,jo) * (1._wp - zalpha) + zalpha * hdivn_lf(:,:,:,jo+1)
+               END DO
             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
-
+               ztu(:,:,jk) = (un(:,:,jk)-un_lf(:,:,jk,1)/e3u_n(:,:,jk)*r1_e2u(:,:))*umask(:,:,jk)
+               ztv(:,:,jk) = (vn(:,:,jk)-vn_lf(:,:,jk,1)/e3v_n(:,:,jk)*r1_e1v(:,:))*vmask(:,:,jk)
+               tilde_e3t_a(:,:,jk) =  tilde_e3t_a(:,:,jk) + hdivn_lf(:,:,jk,1) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
+            END DO
             !
          ELSEIF ( ln_vvl_layer ) THEN
@@ -447,19 +488,19 @@
          !
          ! 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 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
@@ -481,5 +522,5 @@
          ! 
          ! Thickness anomaly diffusion:
-         ! -----------------------------
+         ! ----------------------------
          ztu(:,:,:) = 0.0_wp
          ztv(:,:,:) = 0.0_wp
@@ -498,8 +539,4 @@
                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) )
@@ -542,8 +579,4 @@
                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) )
@@ -561,7 +594,4 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-!                  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  ) 
@@ -575,8 +605,4 @@
          END DO
 
- 
-!         un_td(:,:,:) = un_td(:,:,:) + ztu(:,:,:)
-!         vn_td(:,:,:) = vn_td(:,:,:) + ztv(:,:,:)
-
          CALL lbc_lnk_multi( un_td, 'U', -1., vn_td, 'V', -1. )     !* local domain boundaries
          !
@@ -591,5 +617,5 @@
          ! Remove "external thickness" tendency:
          DO jk = 1, jpkm1
-            tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) +   e3t_n(:,:,jk) * zhdiv(:,:)
+            tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) +   (e3t_n(:,:,jk)-zmet*tilde_e3t_n(:,:,jk)) * zhdiv(:,:) * tmask(:,:,jk)
          END DO 
                    
@@ -646,15 +672,50 @@
          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(:,:,:)
+
+         DO jo = nn_filt_order,1,-1
+            hdivn_lf(:,:,:,jo) =  hdivn_lf(:,:,:,jo)  + zalpha**(nn_filt_order-jo+1) * ze3t(:,:,:)
+         END DO
       ENDIF
 
       IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
       !                                             ! ---baroclinic part--------- !
+
+         IF ( (ncall==2).AND.(.NOT.ll_do_bclinic) ) THEN
+
+            DO jk = 1, jpkm1
+               ztu(:,:,jk) = (un_adv(:,:)*r1_hu_n(:,:) - un_b(:,:) ) * e3u_n(:,:,jk) * e2u(:,:) * umask(:,:,jk)
+               ztv(:,:,jk) = (vn_adv(:,:)*r1_hv_n(:,:) - vn_b(:,:) ) * e3v_n(:,:,jk) * e1v(:,:) * vmask(:,:,jk)
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     ze3t(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk) &
+                                    &  + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk) &
+                                    & ) / ( e1t(ji,jj) * e2t(ji,jj) )
+                  END DO
+               END DO
+            END DO
+            !
+            zhdiv(:,:) = 0.
+            DO jk = 1, jpkm1
+               zhdiv(:,:) = zhdiv(:,:) + ze3t(:,:,jk) * tmask(:,:,jk)
+            END DO
+            zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - tmask_i(:,:) ) 
+            !
+            DO jk = 1, jpkm1
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - z2dt * (ze3t(:,:,jk) & 
+                                   & - zhdiv(:,:)*(e3t_n(:,:,jk)-zmet*tilde_e3t_n(:,:,jk))*tmask(:,:,jk))
+            END DO
+            CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )
+         ENDIF
+         !
          DO jk = 1, jpkm1
-            e3t_a(:,:,jk) = e3t_a(:,:,jk) + (tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk))                     
-         END DO
+            e3t_a(:,:,jk) = e3t_a(:,:,jk) + tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
+         END DO
+!         IF( ln_vvl_ztilde ) THEN !  Relax barotropic component:
+!            DO jk = 1, jpkm1
+!               e3t_a(:,:,jk) = e3t_a(:,:,jk)  &
+!                             & - z2dt * frq_rst_e3t(:,:) * (e3t_b(:,:,jk) - tilde_e3t_b(:,:,jk)   & 
+!                             & - e3t_0(:,:,jk) * (ht_0(:,:) + sshb(:,:))/ (ht_0(:,:)*tmask(:,:,1) + 1._wp - tmask(:,:,1)))                 
+!            END DO
+!         ENDIF
       ENDIF
 
@@ -686,4 +747,5 @@
             ENDIF
             IF (lwp) THEN
+               ji = ijk_min(1) ; jj = ijk_min(2) ; jk = ijk_min(3)
                WRITE(numout, *) 'Negative scale factor, e3t_n =', z_tmin
                WRITE(numout, *) 'at i, j, k=', ijk_min            
@@ -1270,5 +1332,5 @@
       !
       INTEGER ::   ji, jj, jk
-      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
+      INTEGER ::   id1, id2, id3, id4, id5, id6, id7     ! local integers
       !!----------------------------------------------------------------------
       !
@@ -1284,4 +1346,6 @@
             id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
             id5 = iom_varid( numror, 'hdivn_lf', ldstop = .FALSE. )
+            id6 = iom_varid( numror, 'un_lf', ldstop = .FALSE. )
+            id7 = iom_varid( numror, 'vn_lf', ldstop = .FALSE. )
             !                             ! --------- !
             !                             ! all cases !
@@ -1344,8 +1408,12 @@
                IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
                   !                       ! ------------ !
-                  IF( id5 > 0 ) THEN  ! required array exists
-                     CALL iom_get( numror, jpdom_autoglo, 'hdivn_lf', hdivn_lf(:,:,:), ldxios = lrxios )
+                  IF( MIN(id5, id6, id7) > 0 ) THEN  ! required arrays exist
+                     CALL iom_get( numror, jpdom_autoglo, 'hdivn_lf', hdivn_lf(:,:,:,1), ldxios = lrxios )
+                     CALL iom_get( numror, jpdom_autoglo, 'un_lf', un_lf(:,:,:,1), ldxios = lrxios )
+                     CALL iom_get( numror, jpdom_autoglo, 'vn_lf', vn_lf(:,:,:,1), ldxios = lrxios )
                   ELSE                ! array is missing
-                     hdivn_lf(:,:,:) = 0.0_wp
+                     hdivn_lf(:,:,:,:) = 0.0_wp
+                     un_lf(:,:,:,:) = 0.0_wp
+                     vn_lf(:,:,:,:) = 0.0_wp
                   ENDIF
                ENDIF
@@ -1418,7 +1486,7 @@
                tilde_e3t_n(:,:,:) = 0._wp
                IF( ln_vvl_ztilde ) THEN 
-                  hdivn_lf(:,:,:) = 0._wp 
-                  un_lf(:,:,:) = 0._wp 
-                  vn_lf(:,:,:) = 0._wp 
+                  hdivn_lf(:,:,:,:) = 0._wp 
+                  un_lf(:,:,:,:) = 0._wp 
+                  vn_lf(:,:,:,:) = 0._wp 
                ENDIF
             END IF
@@ -1443,5 +1511,7 @@
          IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
             !                                        ! ------------ !
-            CALL iom_rstput( kt, nitrst, numrow, 'hdivn_lf', hdivn_lf(:,:,:), ldxios = lwxios)
+            CALL iom_rstput( kt, nitrst, numrow, 'hdivn_lf', hdivn_lf(:,:,:,1), ldxios = lwxios)
+            CALL iom_rstput( kt, nitrst, numrow, 'un_lf', un_lf(:,:,:,1), ldxios = lwxios)
+            CALL iom_rstput( kt, nitrst, numrow, 'vn_lf', vn_lf(:,:,:,1), ldxios = lwxios)
          ENDIF
          !
@@ -1604,7 +1674,7 @@
       ll_chk_bot2top   = .TRUE.
       ll_chk_top2bot   = .TRUE.
-      dzmin_int  = 0.1_wp   ! Absolute minimum depth in the interior (in meters)
+      dzmin_int  = 1.0_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_stp  = 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
@@ -1614,9 +1684,8 @@
          zvlim   = 0.5_wp   ! max d2h/dh
       ll_lapdiff_cond  = .TRUE.  ! Conditionnal Laplacian diffusion of interfaces
-         zhdiff  = 0.01_wp   ! ad.
+         zhdiff  = 0.01_wp  ! ad.
          zhlim   = 0.03_wp  ! ad. max lap(z)*e1
-      ll_blpdiff_cond  = .FALSE.  ! Conditionnal Bilaplacian diffusion of interfaces
+      ll_blpdiff_cond  = .TRUE.  ! 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
       ! --------------------------------------------------------------------------------------- 
@@ -1819,4 +1888,5 @@
                   zh_old = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
                   zh_min = MIN(zh_0/3._wp, dzmin_int)
+!                  zh_min = MAX(zh_min, zh_min-e3t_a(ji,jj,jk)+e3t_0(ji,jj,jk))
                   ! 
                   ! Set maximum and minimum vertical excursions   
@@ -1836,5 +1906,6 @@
                   ! Ensure minimum layer thickness:                  
 !                  zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new)
-                  zh_new = cush(zh_new, zh_min)
+!                  zh_new = cush(zh_new, zh_min)
+                  zh_new = MAX(zh_new, zh_min)           
                   !
                   ! Final flux:
@@ -1842,6 +1913,6 @@
                   !
                   ! Limit thickness change in 1 time step: 
-                  ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef )
-                  zdiff = SIGN(ztmp, zh_new - zh_old)
+!                  ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef )
+!                  zdiff = SIGN(ztmp, zh_new - zh_old)
                   zh_new = zdiff + zh_old
                   !
@@ -1864,6 +1935,7 @@
                   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)
+!                  zh_min = MAX(zh_min, zh_min-e3t_a(ji,jj,jk)+e3t_0(ji,jj,jk))
+                  !
+!                  zwdw(ji,jj,jk+1) = MAX(zwdw(ji,jj,jk+1), REAL(jk)*dzmin_surf)
                   !
                   ! New layer thickness:
@@ -1872,5 +1944,6 @@
                   ! Ensure minimum layer thickness:
 !                  zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new)
-                  zh_new = cush(zh_new, zh_min)
+!                  zh_new = cush(zh_new, zh_min)
+                  zh_new = MAX(zh_new, zh_min)   
                   !
                   ! Final flux:
@@ -2089,12 +2162,12 @@
       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
+!      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
@@ -2138,7 +2211,17 @@
                ! add them to the general tracer trends
                pta(ji,jj,jk) = pta(ji,jj,jk) + ztra
+               zwi(ji,jj,jk) = (e3t_b(ji,jj,jk) + z2dtt * pta(ji,jj,jk)* bdytmask(ji,jj) ) * tmask(ji,jj,jk)
             END DO
          END DO
       END DO
+
+      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
       !
       IF( ln_timing )  CALL timing_stop('dom_vvl_adv_fct')
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/namelist_cfg	(revision 10163)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/namelist_cfg	(revision 10164)
@@ -237,5 +237,5 @@
       rn_ahe3_lap = 30.e0
    ln_vvl_blp = .true.              !  Bilaplacian thickness diffusion
-      rn_ahe3_blp = -2.e8
+      rn_ahe3_blp = -1.e8
    rn_rst_e3t    = 30.e0            !  ztilde to zstar restoration timescale [days]
    rn_lf_cutoff  = 5.0e0            !  cutoff frequency for low-pass filter  [days]
