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 9528)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 9529)
@@ -96,4 +96,5 @@
          <field id="tpt_dep"      long_name="T-point depth"                  standard_name="depth_below_geoid"   unit="m"   grid_ref="grid_T_3D" />
          <field id="e3tdef"       long_name="T-cell thickness deformation"                                       unit="%"   grid_ref="grid_T_3D" />
+         <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="depw_tilde"   long_name="Interface displacement"                                             unit="m"   grid_ref="grid_W_3D" />
@@ -408,5 +409,4 @@
          <field id="uoces"        long_name="ocean transport along i-axis times salinity (CRS)"                                                  unit="0.001*m/s"   grid_ref="grid_U_3D" />
          <field id="un_lf_tilde"  long_name="Sea Water X Velocity (low pass)"                        standard_name="low_pass_sea_water_x_velocity"     unit="m/s"        grid_ref="grid_U_3D" />
-
 
          <!-- variables available with MLE -->
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 9528)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/cfg.txt	(revision 9529)
@@ -12,4 +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 9528)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 9529)
@@ -157,6 +157,11 @@
          CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
       CALL iom_put("tpt_dep", fsdept_n(:,:,:) )
-
-
+      IF( iom_use("gdepwt") )  THEN
+         z3d(:,:,1) = (gdepw_n(:,:,1)-sshn(:,:))*tmask(:,:,1)
+         DO jk=2,jpk
+            z3d(:,:,jk) = (gdepw_n(:,:,jk)-sshn(:,:))*tmask(:,:,jk-1)
+         END DO
+         CALL iom_put("gdepwt" , z3d(:,:,:) )
+      END IF
 
       CALL iom_put( "ssh" , sshn )                 ! sea surface height
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 9528)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 9529)
@@ -60,4 +60,5 @@
    LOGICAL                                               :: ln_vvl_blp                          ! Bilaplacian thickness diffusion 
    LOGICAL                                               :: ln_vvl_regrid                       ! ensure layer separation 
+   LOGICAL                                               :: ll_shorizd=.FALSE.                  ! Use "shelf horizon depths"         
    !                                                                                            ! conservation: not used yet
    REAL(wp)                                              :: rn_ahe3_lap               ! thickness diffusion coefficient (Laplacian)
@@ -66,8 +67,9 @@
    REAL(wp)                                              :: rn_lf_cutoff              ! cutoff frequency for low-pass filter  [days]
    REAL(wp)                                              :: rn_zdef_max               ! maximum fractional e3t deformation
+   REAL(wp)                                              :: hsmall=0.01               ! small thickness
 
    !! * Module variables
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td              ! thickness diffusion transport
-   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdivn_lf,  un_lf,  vn_lf  ! 1st order filtered arrays         
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_lf, vn_lf, hdivn_lf    ! 1st order filtered arrays   
    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
@@ -226,5 +228,4 @@
          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 )
-         ! tendency mask:
          !
          IF( ln_vvl_ztilde_as_zstar ) THEN
@@ -331,9 +332,11 @@
       INTEGER                                :: ib, ib_bdy, ip, jp    !   "     "     "
       INTEGER , DIMENSION(3)                 :: ijk_max, ijk_min      ! temporary integers
+      INTEGER                                :: ncall
       REAL(wp)                               :: z2dt                  ! temporary scalars
       REAL(wp)                               :: z_tmin, z_tmax        ! temporary scalars
       REAL(wp)                               :: zalpha, zwgt          ! temporary scalars
-      REAL(wp)                               :: zdu, zdv
+      REAL(wp)                               :: zdu, zdv, zramp
       LOGICAL                                :: ll_do_bclinic         ! temporary logical
+      REAL(wp)                               :: ztra, zbtr, ztout, ztin, zfac, zmsku, zmskv
       !!----------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt')
@@ -348,10 +351,19 @@
 
       ll_do_bclinic = .TRUE.
+      ncall=1
       IF( PRESENT(kcall) ) THEN
-         IF ( kcall == 2 .AND. ln_vvl_ztilde.OR.ln_vvl_layer ) ll_do_bclinic = .FALSE.
+         ! 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) ncall=2    
       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 !
       ! ******************************* !
 
@@ -383,16 +395,10 @@
          IF ( ln_vvl_ztilde ) THEN
             !
-            IF ((kt==nit000).AND.(neuler==0)) THEN
-               DO jk = 1, jpkm1
-                  ztu(:,:,jk) = un(:,:,jk)
-                  ztv(:,:,jk) = vn(:,:,jk)
-               END DO
-            ELSE
-               DO jk = 1, jpkm1
-                  ztu(:,:,jk) = (un(:,:,jk)-un_lf(:,:,jk)/fse3u_n(:,:,jk)*r1_e2u(:,:))*umask(:,:,jk)
-                  ztv(:,:,jk) = (vn(:,:,jk)-vn_lf(:,:,jk)/fse3v_n(:,:,jk)*r1_e1v(:,:))*vmask(:,:,jk)
-                  tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) + hdivn_lf(:,:,jk)
-               END DO
-            ENDIF
+            DO jk = 1, jpkm1
+               ztu(:,:,jk) = (un(:,:,jk)-un_lf(:,:,jk)/fse3u_n(:,:,jk)*r1_e2u(:,:))*umask(:,:,jk)
+               ztv(:,:,jk) = (vn(:,:,jk)-vn_lf(:,:,jk)/fse3v_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
@@ -405,4 +411,21 @@
          ENDIF
          !
+         ! 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, fse3u_n(ji,jj,jk) - hsmall) )
+!!                  un_td(ji,jj,jk) = un_td(ji,jj,jk) - (1. - zmsku) * un(ji,jj,jk) * fse3u_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, fse3v_n(ji,jj,jk) - hsmall) )
+!!                  vn_td(ji,jj,jk) = vn_td(ji,jj,jk) - (1. - zmskv) * vn(ji,jj,jk) * fse3v_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
@@ -411,14 +434,28 @@
          ELSEIF (ln_vvl_adv_cn2) THEN
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * hdivn(:,:,jk)
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     tilde_e3t_a(ji,jj,jk) =   &
+                    -(  e2u(ji,jj)*fse3u(ji,jj,jk) * ztu(ji,jj,jk) - e2u(ji-1,jj  )*fse3u(ji-1,jj  ,jk) * ztu(ji-1,jj  ,jk)       &
+                      + e1v(ji,jj)*fse3v(ji,jj,jk) * ztv(ji,jj,jk) - e1v(ji  ,jj-1)*fse3v(ji  ,jj-1,jk) * ztv(ji  ,jj-1,jk)  )    &
+                     / ( e1t(ji,jj) * e2t(ji,jj) )
+                  END DO
+               END DO
             END DO
          ENDIF
          ! 
-         ! Thickness anomlaly diffusion:
+         ! Thickness anomaly diffusion:
          ! -----------------------------
-         zwu(:,:) = 0.0_wp
-         zwv(:,:) = 0.0_wp
          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
@@ -426,8 +463,12 @@
                DO jj = 1, jpjm1                 ! First derivative (gradient)
                   DO ji = 1, fs_jpim1   ! vector opt.                  
+!                     ztu(ji,jj,jk) =  umask(ji,jj,jk) * re2u_e1u(ji,jj) &
+!                                     &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
+!                     ztv(ji,jj,jk) =  vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
+!                                     &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
                      ztu(ji,jj,jk) =  umask(ji,jj,jk) * re2u_e1u(ji,jj) &
-                                     &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
+                                     &  * ( ze3t(ji,jj,jk) - ze3t(ji+1,jj  ,jk) )
                      ztv(ji,jj,jk) =  vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
-                                     &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
+                                     &  * ( ze3t(ji,jj,jk) - ze3t(ji  ,jj+1,jk) )
                   END DO
                END DO
@@ -455,6 +496,4 @@
                      ztu(ji,jj,jk) = umask(ji,jj,jk) * re2u_e1u(ji,jj) * ( zht(ji+1,jj  ) - zht(ji,jj) )
                      ztv(ji,jj,jk) = vmask(ji,jj,jk) * re1v_e2v(ji,jj) * ( zht(ji  ,jj+1) - zht(ji,jj) )
-                     zwu(ji,jj) = zwu(ji,jj) + ztu(ji,jj,jk)
-                     zwv(ji,jj) = zwv(ji,jj) + ztv(ji,jj,jk)
                   END DO
                END DO
@@ -466,10 +505,12 @@
                DO jj = 1, jpjm1
                   DO ji = 1, fs_jpim1   ! vector opt.                  
+!                     zdu =  rn_ahe3_lap * umask(ji,jj,jk) * re2u_e1u(ji,jj) &
+!                         &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
+!                     zdv =  rn_ahe3_lap * vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
+!                         &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
                      zdu =  rn_ahe3_lap * umask(ji,jj,jk) * re2u_e1u(ji,jj) &
-                         &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
+                         &  * ( ze3t(ji,jj,jk) - ze3t(ji+1,jj  ,jk) )
                      zdv =  rn_ahe3_lap * vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
-                         &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
-                     zwu(ji,jj) = zwu(ji,jj) + zdu
-                     zwv(ji,jj) = zwv(ji,jj) + zdv
+                         &  * ( 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
@@ -478,20 +519,4 @@
             END DO
          ENDIF 
-
-         ! Ensure barotropic fluxes are null:
-!         DO jj = 1, jpj
-!            DO ji = 1, jpi
-!               DO jk = 1, jpkm1
-!                  ztu(ji,jj,jk) = ztu(ji,jj,jk) - zwu(ji,jj)*fse3u_b(ji,jj,jk)*hur_b(ji,jj)*umask(ji,jj,jk)
-!                  ztv(ji,jj,jk) = ztv(ji,jj,jk) - zwv(ji,jj)*fse3v_b(ji,jj,jk)*hvr_b(ji,jj)*vmask(ji,jj,jk)
-!               END DO
-!            END DO
-!         END DO
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ztu(ji,jj,mbku(ji,jj)) = ztu(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
-               ztv(ji,jj,mbkv(ji,jj)) = ztv(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
-            END DO
-         END DO
 
          ! divergence of diffusive fluxes
@@ -499,24 +524,33 @@
             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)    &
+!                  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_e12t(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_e12t(ji,jj)
                END DO
             END DO
          END DO
+
  
-         un_td(:,:,:) = un_td(:,:,:) + ztu(:,:,:)
-         vn_td(:,:,:) = vn_td(:,:,:) + ztv(:,:,:)
+!         un_td(:,:,:) = un_td(:,:,:) + ztu(:,:,:)
+!         vn_td(:,:,:) = vn_td(:,:,:) + ztv(:,:,:)
+
          CALL lbc_lnk( un_td , 'U' , -1.)
-         CALL lbc_lnk( vn_td , 'V' , -1.)  
-         !
-         !
-         ! Restoring:
-         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
-
+         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
@@ -526,16 +560,10 @@
          ! Leapfrog time stepping
          ! ~~~~~~~~~~~~~~~~~~~~~~
-         IF( neuler == 0 .AND. kt == nit000 ) THEN
-            z2dt =  rdt
-         ELSE
-            z2dt = 2.0_wp * rdt
-         ENDIF
-
-         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
-
-         ! Revert to zstar locally:
-         ! ~~~~~~~~~~~~~~~~~~~~~~~~
+         zramp = 1._wp
+!         IF (.NOT.ln_rstart) zramp = MIN(MAX( ((kt-nit000)*rdt)/(3._wp*rday),0._wp),1._wp)
+
          DO jk=1,jpkm1
-            tilde_e3t_a(:,:,jk) =  tilde_e3t_a(:,:,jk) * tildemask(:,:) 
+            tilde_e3t_a(:,:,jk) = tilde_e3t_b(:,:,jk) + z2dt * tmask(:,:,jk) * tilde_e3t_a(:,:,jk) &
+                               & * tildemask(:,:) * zramp
          END DO
 
@@ -595,7 +623,9 @@
                CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')
             ENDIF
-         ENDIF
-
-         IF ( ln_vvl_ztilde ) THEN
+         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
@@ -603,21 +633,20 @@
                ztv(:,:,jk) = vn(:,:,jk) * fse3v_n(:,:,jk) * e1v(:,:) + vn_td(:,:,jk)
                ze3t(:,:,jk) = -fse3t_n(:,:,jk) * zhdiv(:,:)
-!
-!               DO jj = 2, jpjm1
-!                  DO ji = fs_2, fs_jpim1   ! vector opt.
-!                     
-!                     ze3t(ji,jj,jk) =   -fse3t_n(ji,jj,jk) * zhdiv(ji,jj)   &
-!                        & + (  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
             !
-            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 * 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 !
+      !                                             ! ---baroclinic part--------- !
+         DO jk = 1, jpkm1
+            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + (tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk))                     
+         END DO
       ENDIF
 
@@ -639,14 +668,12 @@
                                   & (hdivn(:,:,jk) - hdivb(:,:,jk) - zhdiv(:,:))
          END DO
+         DO jk = 1, jpkm1
+            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) - z2dt * fse3t_n(:,:,jk) * & 
+                                  & (hdivn(:,:,jk) - hdivb(:,:,jk) - zhdiv(:,:))
+         END DO
       ENDIF
 
-      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
-      !                                             ! ---baroclinic part--------- !
-         DO jk = 1, jpkm1
-            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + (tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk))                     
-         END DO
-      ENDIF
-
-      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
+
+      IF( ln_vvl_dbg .AND. ( ncall==2 ) ) THEN   ! - ML - test: control prints for debuging
          !
          zht(:,:) = 0.0_wp
@@ -663,17 +690,73 @@
          z_tmin = MINVAL( fse3t_n(:,:,:), mask = tmask(:,:,:) == 1.e0  )
          IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3t_n) =', z_tmin         
+         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3t_n) =', z_tmin
+         IF ( z_tmin .LE. 0._wp ) THEN
+            IF( lk_mpp ) THEN
+               CALL mpp_minloc(fse3t_n(:,:,:), tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
+            ELSE
+               ijk_min = MINLOC( fse3t_n(:,:,:), mask = tmask(:,:,:) == 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, fse3t_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( fse3u_n(:,:,:), mask = umask(:,:,:) == 1.e0 )
          IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
          IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3u_n) =', z_tmin
+         IF ( z_tmin .LE. 0._wp ) THEN
+            IF( lk_mpp ) THEN
+               CALL mpp_minloc(fse3u_n(:,:,:), umask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
+            ELSE
+               ijk_min = MINLOC( fse3u_n(:,:,:), mask = umask(:,:,:) == 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, fse3u_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( fse3v_n(:,:,:), mask = vmask(:,:,:) == 1.e0 )
          IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
          IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3v_n) =', z_tmin
+         IF ( z_tmin .LE. 0._wp ) THEN
+            IF( lk_mpp ) THEN
+               CALL mpp_minloc(fse3v_n(:,:,:), vmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
+            ELSE
+               ijk_min = MINLOC( fse3v_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, fse3v_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( fse3f_n(:,:,:), mask = fmask(:,:,:) == 1.e0 )
          IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
          IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3f_n) =', z_tmin
+         IF ( z_tmin .LE. 0._wp ) THEN
+            IF( lk_mpp ) THEN
+               CALL mpp_minloc(fse3f_n(:,:,:), fmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
+            ELSE
+               ijk_min = MINLOC( fse3f_n(:,:,:), mask = fmask(:,:,:) == 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, fse3f_n =', z_tmin
+               WRITE(numout, *) 'at i, j, k=', ijk_min            
+               CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor')
+            ENDIF
+         ENDIF
          !
          zht(:,:) = 0.0_wp
@@ -956,40 +1039,23 @@
          zout(:,:,jpk)=0._wp
          DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zufim1 = umask(ji-1,jj,jk) * re2u_e1u(ji-1,jj) * ( zwdw(ji-1,jj,jk) - zwdw(ji  ,jj  ,jk) )
-                  zufi   = umask(ji  ,jj,jk) * re2u_e1u(ji  ,jj) * ( zwdw(ji  ,jj,jk) - zwdw(ji+1,jj  ,jk) )
-                  zvfjm1 = vmask(ji,jj-1,jk) * re1v_e2v(ji,jj-1) * ( zwdw(ji,jj-1,jk) - zwdw(ji  ,jj  ,jk) )
-                  zvfj   = vmask(ji,jj  ,jk) * re1v_e2v(ji,jj  ) * ( zwdw(ji,jj  ,jk) - zwdw(ji  ,jj+1,jk) )
-!                  ztmp1  = (zufim1-zufi+zvfjm1-zvfj) * SQRT(r1_e12t(ji,jj))  
-                  ztmp1  = (zufim1-zufi+zvfjm1-zvfj) * r1_e12t(ji,jj) 
-!                  zout(ji,jj,jk) = ABS(ztmp1)*tmask(ji,jj,jk)
-                  zout(ji,jj,jk) = ztmp1
-               END DO
-            END DO            
-         END DO
-         ! Mask open boundaries:
-#if defined key_bdy
-         IF (lk_bdy) THEN
-            DO jk = 1, jpkm1
-               zout(:,:,jk) = zout(:,:,jk) * bdytmask(:,:)
-            END DO
-         ENDIF
-#endif
-         CALL lbc_lnk( zout(:,:,:), 'T', 1. )
-         zwdw(:,:,:) = zout(:,:,:)
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.                  
+                  ua(ji,jj,jk) =  umask(ji,jj,jk) * re2u_e1u(ji,jj) &
+                                  &  * ( zwdw(ji,jj,jk) - zwdw(ji+1,jj  ,jk) )
+                  va(ji,jj,jk) =  vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
+                                  &  * ( zwdw(ji,jj,jk) - zwdw(ji  ,jj+1,jk) )
+               END DO
+            END DO
+         END DO
+            
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zufim1 = umask(ji-1,jj,jk) * re2u_e1u(ji-1,jj) * ( zwdw(ji-1,jj,jk) - zwdw(ji  ,jj  ,jk) )
-                  zufi   = umask(ji  ,jj,jk) * re2u_e1u(ji  ,jj) * ( zwdw(ji  ,jj,jk) - zwdw(ji+1,jj  ,jk) )
-                  zvfjm1 = vmask(ji,jj-1,jk) * re1v_e2v(ji,jj-1) * ( zwdw(ji,jj-1,jk) - zwdw(ji  ,jj  ,jk) )
-                  zvfj   = vmask(ji,jj  ,jk) * re1v_e2v(ji,jj  ) * ( zwdw(ji,jj  ,jk) - zwdw(ji  ,jj+1,jk) )
-!                  ztmp1  = (zufim1-zufi+zvfjm1-zvfj) * SQRT(r1_e12t(ji,jj))  
-                  ztmp1  = (zufim1-zufi+zvfjm1-zvfj) * r1_e12t(ji,jj) 
-                  zout(ji,jj,jk) = ABS(ztmp1)*tmask(ji,jj,jk)
-               END DO
-            END DO            
-         END DO
+                  ztmp1 = ( (ua(ji-1,jj  ,jk) - ua(ji,jj,jk))    &
+                     &  +   (va(ji  ,jj-1,jk) - va(ji,jj,jk)) ) * SQRT(r1_e12t(ji,jj))
+                  zout(ji,jj,jk) = ABS(ztmp1)*tmask(ji,jj,jk)           
+               END DO
+            END DO
+         END DO 
          ! Mask open boundaries:
 #if defined key_bdy
@@ -1065,5 +1131,5 @@
       INTEGER ::   jkbot                                               ! bottom level
       LOGICAL ::   l_is_orca                                           ! local logical
-      REAL(wp) :: zmin, zdo, zup, ztap
+      REAL(wp) :: zmin, zdo, zup, ztap, zsmall
       REAL(wp), POINTER, DIMENSION(:,:)   :: zs                        ! Surface interface depth anomaly
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zw                        ! Interface depth anomaly
@@ -1075,9 +1141,10 @@
 
  
-      nmet=1 ! Original method (Surely wrong)
-!      nmet=1 ! Interface interpolation
-!      nmet=2 ! Internal interfaces interpolation only, spread barotropic increment
-!      nmet=1
-      ztap = 0.1_wp ! Minimum fraction of T-point thickness at cell interfaces
+!      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
 
       IF ( (nmet==1).OR.(nmet==2) ) THEN
@@ -1091,25 +1158,31 @@
             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
+!            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
-                     jkbot = mbkt(ji,jj)
-                     DO jk=jkbot,1,-1
-                        zw(ji,jj,jk) = zw(ji,jj,jk+1) - ( pe3_in(ji,jj,jk)                     &
+                     DO jk=1,jpk
+                        zw(ji,jj,jk) = (zw(ji,jj,jk) + zs(ji,jj))                                          &
                                      & * ht_0(ji,jj) / (ht_0(ji,jj) + zs(ji,jj) + 1._wp - tmask(ji,jj,1))  &
-                                     & - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk)
+                                     & * tmask(ji,jj,jk)
                      END DO
                   END DO
@@ -1144,28 +1217,24 @@
                   ! Correction at last level:
                   jkbot = mbku(ji,jj)
-                  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))
+                  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)
                      ENDIF
-                     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)
+                     zup = MIN(zup, zdo-zsmall)
+                     pe3_out(ji,jj,jk) = zdo - zup - e3u_0(ji,jj,jk)
+                     zdo = zup
                   END DO
                END DO
@@ -1217,29 +1286,24 @@
                   ! Correction at last level:
                   jkbot = mbkv(ji,jj)
-                  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))
+                  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)   )
+                     !
+                     ! 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
-                     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) )
-                     !
-                     pe3_out(ji,jj,jk) = zdo - zup 
-                     zdo =  zdo - pe3_out(ji,jj,jk)
+                     zup = MIN(zup, zdo-zsmall)
+                     pe3_out(ji,jj,jk) = zdo - zup - e3v_0(ji,jj,jk)
+                     zdo = zup
                   END DO
                END DO
@@ -1294,35 +1358,6 @@
                   ! bottom correction:
                   jkbot = MIN(mbku(ji,jj), mbku(ji,jj+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
+                  zdo = hf_0(ji,jj)
+                  DO jk=jkbot,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)  & 
@@ -1330,11 +1365,33 @@
                            &           + 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 
                      !
-                     zdo =  zdo - pe3_out(ji,jj,jk)                                    
+                     ! 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
                   END DO
-                  !
-               END DO
-            END DO 
+               END DO
+            END DO
          ENDIF
          !
@@ -1432,8 +1489,7 @@
             id(3) = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
             id(4) = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
-            id(5) = iom_varid( numror, 'hdivn_lf', ldstop = .FALSE. )
-            id(6) = iom_varid( numror, 'un_lf'   , ldstop = .FALSE. )
-            id(7) = iom_varid( numror, 'vn_lf'   , ldstop = .FALSE. )
-
+            id(5) = iom_varid( numror, 'un_lf'   , ldstop = .FALSE. )
+            id(6) = iom_varid( numror, 'vn_lf'   , ldstop = .FALSE. )
+            id(7) = iom_varid( numror, 'hdivn_lf', ldstop = .FALSE. )
             !                             ! --------- !
             !                             ! all cases !
@@ -1498,11 +1554,11 @@
                   !                       ! ------------ !
                   IF( MINVAL(id(5:7) ) > 0 ) THEN  ! all required arrays exist
-                     CALL iom_get( numror, jpdom_autoglo, 'hdivn_lf', hdivn_lf(:,:,:) )
                      CALL iom_get( numror, jpdom_autoglo, 'un_lf'   ,    un_lf(:,:,:) )
                      CALL iom_get( numror, jpdom_autoglo, 'vn_lf'   ,    vn_lf(:,:,:) )
+                     CALL iom_get( numror, jpdom_autoglo, 'hdivn_lf', hdivn_lf(:,:,:) )
                   ELSE                            ! one at least array is missing
-                     hdivn_lf(:,:,:) = 0.0_wp
                      un_lf(:,:,:)    = 0.0_wp
                      vn_lf(:,:,:)    = 0.0_wp
+                     hdivn_lf(:,:,:) = 0.0_wp
                      neuler = 0
                   ENDIF
@@ -1520,7 +1576,7 @@
             END IF
             IF( ln_vvl_ztilde ) THEN
-               hdivn_lf(:,:,:) = 0.0_wp
                un_lf(:,:,:)    = 0.0_wp
                vn_lf(:,:,:)    = 0.0_wp
+               hdivn_lf(:,:,:) = 0.0_wp
             ENDIF
          ENDIF
@@ -1543,5 +1599,4 @@
          IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
             !                                        ! ------------ !
-            CALL iom_rstput( kt, nitrst, numrow, 'hdivn_lf', hdivn_lf(:,:,:) )
             CALL iom_rstput( kt, nitrst, numrow, 'un_lf'   , un_lf(:,:,:)    )
             CALL iom_rstput( kt, nitrst, numrow, 'vn_lf'   , vn_lf(:,:,:)    )
@@ -1651,4 +1706,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()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' )
@@ -2099,11 +2162,11 @@
                 za3 = 0.5_wp * (za1 + za2) 
                 zdiff = ABS(za3-za4)/za4
-!                IF (zdiff>=0.8) THEN
-!                   zwt(ji,jj,jk) =  zr_tscale * MIN(zdiff,1._wp) * za3 / p2dt * tmask(ji,jj,jk)
-                   zwt(ji,jj,jk) =  dsm(ji,jj)/ht_0(ji,jj)*(1._wp-tanh((mbkt(ji,jj)+1-jk)*0.2))*tmask(ji,jj,jk)
-
-!                ELSE
-!                   zwt(ji,jj,jk) = 0.e0*tmask(ji,jj,jk)
-!                ENDIF
+                IF (zdiff>=0.8) THEN
+                   zwt(ji,jj,jk) =  zr_tscale * MIN(zdiff,1._wp) * za3 / p2dt * tmask(ji,jj,jk)
+!                   zwt(ji,jj,jk) =  dsm(ji,jj)/ht_0(ji,jj)*(1._wp-tanh((mbkt(ji,jj)+1-jk)*0.2))*tmask(ji,jj,jk)
+
+                ELSE
+                   zwt(ji,jj,jk) = 0.e0*tmask(ji,jj,jk)
+                ENDIF
             END DO
          END DO
@@ -2300,5 +2363,5 @@
       !!
       !! ** 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.      
@@ -2318,8 +2381,8 @@
       INTEGER                             :: ji, jj, jk ! dummy loop indices
       LOGICAL                             :: ll_chk_bot2top, ll_chk_top2bot, ll_lapdiff_cond
-      LOGICAL                             :: ll_e3tsurf_const, ll_zdiff_cond, ll_blpdiff_cond
+      LOGICAL                             :: ll_zdiff_cond, ll_blpdiff_cond
       INTEGER                             :: jkbot
-      REAL(wp)                            :: zh_min, zh_0, zh2, zdiff, zh_max, ztmph, ztmpd, dzmin_vvl
-      REAL(wp)                            :: zufim1, zufi, zvfjm1, zvfj
+      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
@@ -2336,19 +2399,20 @@
       ll_chk_bot2top   = .TRUE.
       ll_chk_top2bot   = .TRUE.
-      ll_e3tsurf_const = .FALSE.
-      dzmin_vvl  = 1._wp   ! Absolute minimum depth (in meters)
-      zfrch_stp  = 0.5_wp   ! Maximum fractionnal thickness change in one time step (<= 1.)
-      zfrch_rel  = 0.0_wp   ! Maximum relative thickness change in the vertical (<= 1.)
-      zfrac_bot  = 0.01_wp  ! Fraction of bottom level allowed to change
-      zscal_bot  = 2._wp    ! Depth lengthscale
-      ll_zdiff_cond    = .FALSE.  ! Conditionnal vertical diffusion of interfaces
-         zvdiff  = 0.1_wp   ! m/s
-         zvlim   = 0.4_wp   ! max d2h/dh
-      ll_lapdiff_cond  = .FALSE.  ! Conditionnal Laplacian diffusion of interfaces
-         zhdiff  = 0.1_wp   ! ad.
-         zhlim   = 0.01_wp  ! max lap(z)/e1
-      ll_blpdiff_cond  = .FALSE.  ! Conditionnal Bilaplacian diffusion of interfaces
-         zhdiff2 = 0.8_wp   ! ad.
-         zhlim2  = 5.e-11_wp  ! max bilap(z)/e1**3
+      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_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.1_wp   ! m
+         zvlim   = 0.5_wp   ! max d2h/dh
+      ll_lapdiff_cond  = .TRUE.  ! Conditionnal Laplacian diffusion of interfaces
+         zhdiff  = 0.2_wp   ! ad.
+         zhlim   = 0.03_wp  ! ad. max lap(z)*e1
+      ll_blpdiff_cond  = .TRUE.  ! Conditionnal Bilaplacian diffusion of interfaces
+         zhdiff2 = 0.2_wp   ! ad.
+!         zhlim2  = 3.e-11_wp  !  max bilap(z)
+         zhlim2  = 0.03_wp  ! ad. max bilap(z)*e1**3
       ! --------------------------------------------------------------------------------------- 
       !
@@ -2372,5 +2436,5 @@
                           &     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 ) 
+                zdw(ji,jj) = MAX(zscal_bot * zdw(ji,jj), rsmall )
             END DO
          END DO
@@ -2385,5 +2449,5 @@
                   &                           + 4._wp*  zdw(ji  ,jj  )                       )
             END DO
-         END DO
+         END DO         
 
          CALL lbc_lnk( dsm(:,:), 'T', 1. )
@@ -2395,4 +2459,5 @@
                   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
@@ -2402,4 +2467,5 @@
                   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
@@ -2515,5 +2581,5 @@
                   ztmp1 = ( (ua(ji-1,jj  ,jk) - ua(ji,jj,jk))    &
                      &  +   (va(ji  ,jj-1,jk) - va(ji,jj,jk)) ) * r1_e12t(ji,jj)
-                  zh2 = MAX(abs(ztmp1)-zhlim2, 0._wp)
+                  zh2 = MAX(abs(ztmp1)-zhlim2*SQRT(r1_e12t(ji,jj))*r1_e12t(ji,jj), 0._wp)
                   ztmp = SIGN(zh2, ztmp1)
                   zeu2 = zhdiff2 * e12t(ji,jj)*e12t(ji,jj) / 16._wp
@@ -2555,13 +2621,13 @@
                   !
                   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 = tilde_e3t_a(ji,jj,jk  ) + zh_0
-!                  zh_dwn = tilde_e3t_a(ji,jj,jk+1) + e3t_0(ji,jj,jk+1)
-                  zh_min = MIN(zh_0/3._wp, dzmin_vvl)
+                  zh_bef = tilde_e3t_b(ji,jj,jk) + zh_0
+                  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)
                   ! 
                   ! 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)+1))/ztmpd)               
+                  zh2   = ztmph * exp(-(gdepw_0(ji,jj,jk)-gdepw_0(ji,jj,i_int_bot(ji,jj)+1))/ztmpd)
                   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)
@@ -2573,5 +2639,6 @@
                   !
                   ! Ensure minimum layer thickness:                  
-!                  zh_new = MIN(zh_old, zh_dwn * zfrch_rel / (2._wp-zfrch_rel) )
+!                  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 = cush(zh_new, zh_min)
                   !
@@ -2579,12 +2646,12 @@
                   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
+                  ! 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)
+!                  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)            
+!                  tilde_e3t_a(ji,jj,jk-1) = (-zdiff + tilde_e3t_a(ji,jj,jk-1) ) * tmask(ji,jj,jk-1)         
                END DO    
             END DO 
@@ -2601,13 +2668,18 @@
                   !
                   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 = tilde_e3t_a(ji,jj,jk  ) + zh_0   
-!                  zh_up  = tilde_e3t_a(ji,jj,jk-1) + e3t_0(ji,jj,jk-1)
-                  zh_min = MIN(zh_0/3._wp, dzmin_vvl)
-                  IF ((jk<=5).AND.ll_e3tsurf_const) zh_min = MAX(e3t_0(ji,jj,1)/3._wp, dzmin_vvl)
+                  zh_bef = tilde_e3t_b(ji,jj,jk) + zh_0
+                  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)
+                  !
+                  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=MIN(zh_old, zh_up * zfrch_rel / (2._wp-zfrch_rel) )
-                  zh_new = cush(zh_old, zh_min)
+!                  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 = cush(zh_new, zh_min)
                   !
                   ! Final flux:
@@ -2615,11 +2687,11 @@
                   ! 
                   ! Limit flux:                 
-                  ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef )
-                  zdiff = SIGN(ztmp, zdiff)
-                  zh_new = zdiff + zh_old
+!                  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)
+!                  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)
+!                  tilde_e3t_a(ji,jj,jk+1) = (-zdiff + tilde_e3t_a(ji,jj,jk+1) ) * tmask(ji,jj,jk+1)
                END DO
                !               
@@ -2627,4 +2699,12 @@
          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
       !
       CALL wrk_dealloc( jpi, jpj, jpk, zwdw ) 
@@ -2703,13 +2783,11 @@
       !! **  Action  : - Update pta thickness tendency and diffusive fluxes
       !!               - this is the total trend, hence it does include sea level motions
-      !!               - Upstream corrections to antidiffusive fluxes ensure
-      !!                 that barotropic transport matches what is contained in input fluxes
       !!----------------------------------------------------------------------
       !
       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(in)    ::   uin, vin  ! input velocities
-      !
-      INTEGER  ::   ji, jj, jk               ! dummy loop indices  
+      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
@@ -2719,5 +2797,5 @@
       REAL(wp) ::   zfp_hi, zfp_hj           !   -      -
       REAL(wp) ::   zfm_hi, zfm_hj           !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:)   :: zbu, zbv, zhu_b, zhv_b
+      REAL(wp) ::   ztout,  ztin, zfac       !   -      -
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zwy, zwi
       !!----------------------------------------------------------------------
@@ -2725,5 +2803,4 @@
       IF( nn_timing == 1 )  CALL timing_start('dom_vvl_adv_fct')
       !
-      CALL wrk_alloc( jpi, jpj, zhu_b, zhv_b, zbu, zbv)
       CALL wrk_alloc( jpi, jpj, jpk, zwx, zwy, zwi)
       !
@@ -2744,42 +2821,51 @@
       ! 2. upstream advection with initial mass fluxes & intermediate update
       ! --------------------------------------------------------------------
-      DO jk = 1, jpkm1
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               !
-               zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) )
-               zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) )
-               zfp_hi = fse3t_b(ji  ,jj  ,jk)
-               zfm_hi = fse3t_b(ji+1,jj  ,jk)
-               zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(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) )
-               zfp_hj = fse3t_b(ji  ,jj  ,jk)
-               zfm_hj = fse3t_b(ji  ,jj+1,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
-
-      IF ( .NOT.ln_sco ) THEN
-         ! Correct bottom upstream fluxes
-         ! considering "shelf horizon depths"
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               ! upstream scheme
-               ikbu = mbku(ji,jj)
-               ikbv = mbkv(ji,jj)
-               zfp_ui = uin(ji,jj,ikbu) + ABS( uin(ji,jj,ikbu) )
-               zfm_ui = uin(ji,jj,ikbu) - ABS( uin(ji,jj,ikbu) )
-               zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbu), 0._wp)
-               zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj  ,ikbu), 0._wp)
-               zwx(ji,jj,ikbu) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,ikbu)
-               !
-               zfp_vj = vin(ji,jj,ikbv) + ABS( vin(ji,jj,ikbv) )
-               zfm_vj = vin(ji,jj,ikbv) - ABS( vin(ji,jj,ikbv) )
-               zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbv), 0._wp)
-               zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj+1,ikbv), 0._wp)
-               zwy(ji,jj,ikbv) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,ikbv)
+      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) - fsdepw_b(ji  ,jj  ,jk), 0._wp)
+                  zfp_hi = MIN(zfp_hi, fse3t_b(ji  ,jj  ,jk))
+                  zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) ) 
+                  !
+                  zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj  ,jk), 0._wp)
+                  zfm_hi = MIN(zfm_hi, fse3t_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) - fsdepw_b(ji  ,jj  ,jk), 0._wp)
+                  zfp_hj = MIN(zfp_hj, fse3t_b(ji  ,jj  ,jk))
+                  zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) ) 
+                  !
+                  zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj+1,jk), 0._wp)
+                  zfm_hj = MIN(zfm_hj, fse3t_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 = fse3t_b(ji  ,jj  ,jk)
+                  zfm_hi = fse3t_b(ji+1,jj  ,jk)             
+                  zfp_hj = fse3t_b(ji  ,jj  ,jk)
+                  zfm_hj = fse3t_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
@@ -2801,6 +2887,18 @@
          END DO
       END DO
-      !                             ! Lateral boundary conditions on zwi  (unchanged sign)
+
       CALL lbc_lnk( zwi, 'T', 1. )  
+
+#if defined key_bdy
+         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) = fse3t_a(ji,jj,jk)
+               END DO
+            END DO 
+         END DO
+#endif
 
       IF ( ln_vvl_dbg ) THEN
@@ -2808,5 +2906,6 @@
          IF( lk_mpp )   CALL mpp_min( zmin )
          IF( zmin < 0._wp) THEN
-            IF(lwp) CALL ctl_stop('vvl_adv: CFL issue here')
+            IF(lwp) CALL ctl_warn('vvl_adv: CFL issue here')
+            IF(lwp) WRITE(numout,*) zmin
          ENDIF
       ENDIF
@@ -2855,8 +2954,43 @@
          END DO
       END DO
-
-      !
-      ! 6. Correct barotropic flux
-      ! --------------------------
+      !
+      CALL wrk_dealloc( jpi, jpj, jpk, zwx, zwy, zwi)
+      !
+      IF( nn_timing == 1 )  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), POINTER, DIMENSION(:,:)   :: zbu, zbv, zhu_b, zhv_b
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zwy
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_ups_cor')
+      !
+      CALL wrk_alloc( jpi, jpj, zhu_b, zhv_b, zbu, zbv)
+      CALL wrk_alloc( jpi, jpj, jpk, zwx, zwy)
+      !
       ! Compute barotropic flux difference:
       zbu(:,:) = 0.e0
@@ -2865,6 +2999,6 @@
          DO ji = 1, jpi   ! vector opt.
             DO jk = 1, jpkm1
-               zbu(ji,jj) = zbu(ji,jj) - un_td(ji,jj,jk) * umask(ji,jj,jk)
-               zbv(ji,jj) = zbv(ji,jj) - vn_td(ji,jj,jk) * vmask(ji,jj,jk)
+               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
@@ -2874,31 +3008,6 @@
       zhu_b(:,:) = 0.e0
       zhv_b(:,:) = 0.e0
-      IF ( ln_sco ) THEN; ibot=0 ; ELSE ; ibot=1 ; ENDIF
-
-      DO jj = 1, jpjm1
-         DO ji = 1, jpim1   ! vector opt.
-            zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj))
-            ikbu = mbku(ji,jj) - ibot 
-            DO jk = 1, ikbu
-               zfp_hi = fse3t_b(ji  ,jj  ,jk)
-               zfm_hi = fse3t_b(ji+1,jj  ,jk)
-               zhu_b(ji,jj) = zhu_b(ji,jj) + (          zdi  * zfp_hi &
-                            &                  + (1._wp-zdi) * zfm_hi &
-                            &                ) * umask(ji,jj,jk)
-            END DO
-            !
-            zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj))
-            ikbv = mbkv(ji,jj) - ibot
-            DO jk = 1, ikbv
-               zfp_hj = fse3t_b(ji  ,jj  ,jk)
-               zfm_hj = fse3t_b(ji  ,jj+1,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
-
-      IF ( .NOT.ln_sco ) THEN
+
+      IF ( ll_shorizd ) THEN
          ! Correct bottom value
          ! considering "shelf horizon depth"
@@ -2907,16 +3016,49 @@
                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))
-               ikbu = mbku(ji,jj)
-               ikbv = mbkv(ji,jj)
-               zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbu), 0._wp)
-               zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj  ,ikbu), 0._wp)
-               zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbv), 0._wp)
-               zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj+1,ikbv), 0._wp)
-               zhu_b(ji,jj) = zhu_b(ji,jj) + (         zdi  * zfp_hi &
-                            &                 + (1._wp-zdi) * zfm_hi &
-                            &                ) * umask(ji,jj,ikbu)
-               zhv_b(ji,jj) = zhv_b(ji,jj) + (         zdj  * zfp_hj &
-                            &                 + (1._wp-zdj) * zfm_hj &
-                            &                ) * vmask(ji,jj,ikbv)
+               DO jk=1, jpkm1
+                  zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji  ,jj  ,jk), 0._wp)
+                  zfp_hi = MIN(zfp_hi, fse3t_b(ji  ,jj  ,jk))
+                  zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) ) 
+                  !
+                  zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj  ,jk), 0._wp)
+                  zfm_hi = MIN(zfm_hi, fse3t_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) - fsdepw_b(ji  ,jj  ,jk), 0._wp)
+                  zfp_hj = MIN(zfp_hj, fse3t_b(ji  ,jj  ,jk))
+                  zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) ) 
+                  !
+                  zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj+1,jk), 0._wp)
+                  zfm_hj = MIN(zfm_hj, fse3t_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 = fse3t_b(ji  ,jj  ,jk)
+                  zfm_hi = fse3t_b(ji+1,jj  ,jk)
+                  zfp_hj = fse3t_b(ji  ,jj  ,jk)
+                  zfm_hj = fse3t_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
@@ -2929,5 +3071,4 @@
       CALL lbc_lnk( zbu(:,:), 'U', -1. )
       CALL lbc_lnk( zbv(:,:), 'V', -1. )
-      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral boundary conditions
       
       ! Set corrective fluxes in upstream direction:
@@ -2935,55 +3076,58 @@
       zwx(:,:,:) = 0.e0
       zwy(:,:,:) = 0.e0
-      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 = fse3t_b(ji  ,jj  ,jk)
-               zfm_hi = fse3t_b(ji+1,jj  ,jk)
-               !
-               zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
-
-               zfp_hj = fse3t_b(ji  ,jj  ,jk)
-               zfm_hj = fse3t_b(ji  ,jj+1,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
-
-      IF ( .NOT.ln_sco ) THEN
-      ! Bottom correction:
-      DO jj = 1, jpjm1
-         DO ji = 1, fs_jpim1   ! vector opt.
-            ! upstream scheme
-            ikbu = mbku(ji,jj)
-            ikbv = mbkv(ji,jj)
-            !
-            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) )
-            !
-            zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbu), 0._wp)
-            zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj  ,ikbu), 0._wp)
-            !
-            zwx(ji,jj,ikbu) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi )
-            !
-            zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbv), 0._wp)
-            zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj+1,ikbv), 0._wp)
-            !
-            zwy(ji,jj,ikbv) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj )
-         END DO
-      END DO
+      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) - fsdepw_b(ji  ,jj  ,jk), 0._wp)
+                  zfp_hi = MIN(fse3t_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) - fsdepw_b(ji+1,jj  ,jk), 0._wp)
+                  zfm_hi = MIN(fse3t_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) - fsdepw_b(ji  ,jj  ,jk), 0._wp)
+                  zfp_hj = MIN(fse3t_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) - fsdepw_b(ji  ,jj+1,jk), 0._wp)
+                  zfm_hj = MIN(fse3t_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 = fse3t_b(ji  ,jj  ,jk)
+                  zfm_hi = fse3t_b(ji+1,jj  ,jk)
+                  zfp_hj = fse3t_b(ji  ,jj  ,jk)
+                  zfm_hj = fse3t_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
 
-      un_td(:,:,:) = un_td(:,:,:) + zwx(:,:,:)
-      vn_td(:,:,:) = vn_td(:,:,:) + zwy(:,:,:)
+      uin(:,:,:) = uin(:,:,:) + zwx(:,:,:)
+      vin(:,:,:) = vin(:,:,:) + zwy(:,:,:)
       !
       ! Update trend with corrective fluxes:
@@ -3001,11 +3145,11 @@
          END DO
       END DO
-
+      !
       CALL wrk_dealloc( jpi, jpj, zhu_b, zhv_b, zbu, zbv)
-      CALL wrk_dealloc( jpi, jpj, jpk, zwx, zwy, zwi)
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_adv_fct')
-      !
-   END SUBROUTINE dom_vvl_adv_fct
+      CALL wrk_dealloc( jpi, jpj, jpk, zwx, zwy)
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_ups_cor')
+      !
+   END SUBROUTINE dom_vvl_ups_cor
 
    SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, p2dt )
Index: /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 9528)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 9529)
@@ -198,4 +198,5 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   jk                     ! dummy loop indices
+      LOGICAL  ::   ll_e3_dep
       REAL(wp) ::   zt, zw                 ! temporary scalars
       REAL(wp) ::   zsur, za0, za1, zkth   ! Values set from parameters in
@@ -209,4 +210,5 @@
       ! Set variables from parameters
       ! ------------------------------
+       ll_e3_dep = .TRUE.
        zkth = ppkth       ;   zacr = ppacr
        zdzmin = ppdzmin   ;   zhmax = pphmax
@@ -312,5 +314,5 @@
       ENDIF
 
-      IF ( ln_isfcav ) THEN
+      IF ( ln_isfcav .OR. ll_e3_dep ) THEN
 ! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
 ! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively
Index: /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 9528)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 9529)
@@ -144,5 +144,5 @@
 
    
-   SUBROUTINE wzv( kt )
+   SUBROUTINE wzv( kt, kcall )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE wzv  ***
@@ -161,7 +161,9 @@
       !
       INTEGER, INTENT(in) ::   kt           ! time step
+      INTEGER, INTENT( in ), OPTIONAL     :: kcall   ! optional argument
       REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
       !
+      LOGICAL             ::   ll_use_totflx
       INTEGER             ::   ji, jj, jk   ! dummy loop indices
       REAL(wp)            ::   z1_2dt       ! local scalars
@@ -170,4 +172,9 @@
       IF( nn_timing == 1 )  CALL timing_start('wzv')
       !
+      ll_use_totflx=.FALSE.
+      IF (( PRESENT(kcall) ).AND.(ln_vvl_ztilde .OR. ln_vvl_layer)) THEN
+         IF ( kcall==2 ) ll_use_totflx=.TRUE.
+      ENDIF
+
       IF( kt == nit000 ) THEN
          !
@@ -185,5 +192,5 @@
       IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
       !
-      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
+      IF (( ln_vvl_ztilde .OR. ln_vvl_layer ).AND.ll_use_totflx) THEN      ! z_tilde and layer cases
          CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 
          !
Index: /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 9528)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 9529)
@@ -222,5 +222,5 @@
                                   CALL div_cur( kstp )         ! Horizontal divergence & Relative vorticity (2nd call in time-split case)
           IF( lk_vvl     )        CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component)
-                                  CALL wzv           ( kstp )  ! now cross-level velocity 
+                                  CALL wzv           ( kstp, kcall=2 )  ! now cross-level velocity 
       ENDIF
 
