Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/GYRE/EXP00/namelist
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/GYRE/EXP00/namelist	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/GYRE/EXP00/namelist	(revision 2905)
@@ -500,4 +500,5 @@
 !!                      ***  Dynamics namelists  ***
 !!======================================================================
+!!   nam_vvl       vertical coordinate options
 !!   namdyn_adv    formulation of the momentum advection
 !!   namdyn_vor    advection scheme
@@ -507,4 +508,13 @@
 !!======================================================================
 !
+!-----------------------------------------------------------------------
+&nam_vvl    !   vertical coordinate options
+!-----------------------------------------------------------------------
+&nam_vvl
+   ln_vvl_zstar  = .false. !  zstar vertical coordinate                    
+   ln_vvl_ztilde = .false. !  hybrid verticalcoordinate: only high frequency variations
+   ln_vvl_layer  = .false. !  full layer vertical coordinate
+   ahe3          = 0.e0    !  thickness diffusion coefficient 
+/
 !-----------------------------------------------------------------------
 &namdyn_adv    !   formulation of the momentum advection
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist	(revision 2905)
@@ -500,4 +500,5 @@
 !!                      ***  Dynamics namelists  ***
 !!======================================================================
+!!   nam_vvl       vertical coordinate options
 !!   namdyn_adv    formulation of the momentum advection
 !!   namdyn_vor    advection scheme
@@ -507,4 +508,13 @@
 !!======================================================================
 !
+!-----------------------------------------------------------------------
+&nam_vvl    !   vertical coordinate options
+!-----------------------------------------------------------------------
+&nam_vvl
+   ln_vvl_zstar  = .false. !  zstar vertical coordinate                    
+   ln_vvl_ztilde = .false. !  hybrid verticalcoordinate: only high frequency variations
+   ln_vvl_layer  = .false. !  full layer vertical coordinate
+   ahe3          = 0.e0    !  thickness diffusion coefficient 
+/
 !-----------------------------------------------------------------------
 &namdyn_adv    !   formulation of the momentum advection
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef.xml
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef.xml	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef.xml	(revision 2905)
@@ -48,4 +48,8 @@
 	<field id="botpres"      description="Pressure at sea floor"                     unit="dbar"                     />
 	<field id="cellthc"      description="Cell thickness"                            unit="m"     axis_ref="deptht"  />
+	<!-- variables available with key_vvl -->
+	<field id="e3t_n"        description="vertical scale factor"                     unit="m"                        />
+	<field id="e3tdef"       description="vertical scale factor variance"            unit="%^2"                      />
+	<field id="dept_n"       description="depth"                                     unit="m"                        />
      </group>
 
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef_ar5.xml
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef_ar5.xml	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef_ar5.xml	(revision 2905)
@@ -48,4 +48,8 @@
 	<field id="botpres"      description="Pressure at sea floor"                     unit="dbar"                     />
 	<field id="cellthc"      description="Cell thickness"                            unit="m"     axis_ref="deptht"  />
+	<!-- variables available with key_vvl -->
+	<field id="e3t_n"        description="vertical scale factor"                     unit="m"                        />
+	<field id="e3tdef"       description="vertical scale factor variance"            unit="%^2"                      />
+	<field id="dept_n"       description="depth"                                     unit="m"                        />
      </group>
 
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 2905)
@@ -500,4 +500,5 @@
 !!                      ***  Dynamics namelists  ***
 !!======================================================================
+!!   nam_vvl       vertical coordinate options
 !!   namdyn_adv    formulation of the momentum advection
 !!   namdyn_vor    advection scheme
@@ -507,4 +508,13 @@
 !!======================================================================
 !
+!-----------------------------------------------------------------------
+&nam_vvl    !   vertical coordinate options
+!-----------------------------------------------------------------------
+&nam_vvl
+   ln_vvl_zstar  = .true.  !  zstar vertical coordinate                    
+   ln_vvl_ztilde = .false. !  hybrid verticalcoordinate: only high frequency variations
+   ln_vvl_layer  = .false. !  full layer vertical coordinate
+   ahe3          = 0.e0    !  thickness diffusion coefficient 
+/
 !-----------------------------------------------------------------------
 &namdyn_adv    !   formulation of the momentum advection
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/cpp_ORCA2_LIM.fcm
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/cpp_ORCA2_LIM.fcm	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/cpp_ORCA2_LIM.fcm	(revision 2905)
@@ -1,1 +1,1 @@
- bld::tool::fppkeys key_trabbl key_vectopt_loop key_orca_r2 key_lim2 key_dynspg_flt key_diaeiv key_ldfslp key_traldf_c2d key_traldf_eiv key_dynldf_c3d key_dtatem key_dtasal key_tradmp key_zdftke key_zdfddm key_zdftmx key_iomput 
+ bld::tool::fppkeys key_trabbl key_vectopt_loop key_orca_r2 key_lim2 key_dynspg_flt key_diaeiv key_ldfslp key_traldf_c2d key_traldf_eiv key_dynldf_c3d key_dtatem key_dtasal key_tradmp key_zdftke key_zdfddm key_zdftmx key_iomput key_vvl
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist	(revision 2905)
@@ -500,4 +500,5 @@
 !!                      ***  Dynamics namelists  ***
 !!======================================================================
+!!   nam_vvl       vertical coordinate options
 !!   namdyn_adv    formulation of the momentum advection
 !!   namdyn_vor    advection scheme
@@ -508,4 +509,13 @@
 !!======================================================================
 !
+!-----------------------------------------------------------------------
+&nam_vvl    !   vertical coordinate options
+!-----------------------------------------------------------------------
+&nam_vvl
+   ln_vvl_zstar  = .false. !  zstar vertical coordinate                    
+   ln_vvl_ztilde = .false. !  hybrid verticalcoordinate: only high frequency variations
+   ln_vvl_layer  = .false. !  full layer vertical coordinate
+   ahe3          = 0.e0    !  thickness diffusion coefficient 
+/
 !-----------------------------------------------------------------------
 &namdyn_adv    !   formulation of the momentum advection
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/POMME/EXP00/namelist
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/POMME/EXP00/namelist	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/POMME/EXP00/namelist	(revision 2905)
@@ -500,4 +500,5 @@
 !!                      ***  Dynamics namelists  ***
 !!======================================================================
+!!   nam_vvl       vertical coordinate options
 !!   namdyn_adv    formulation of the momentum advection
 !!   namdyn_vor    advection scheme
@@ -507,4 +508,13 @@
 !!======================================================================
 !
+!-----------------------------------------------------------------------
+&nam_vvl    !   vertical coordinate options
+!-----------------------------------------------------------------------
+&nam_vvl
+   ln_vvl_zstar  = .false. !  zstar vertical coordinate                    
+   ln_vvl_ztilde = .false. !  hybrid verticalcoordinate: only high frequency variations
+   ln_vvl_layer  = .false. !  full layer vertical coordinate
+   ahe3          = 0.e0    !  thickness diffusion coefficient 
+/
 !-----------------------------------------------------------------------
 &namdyn_adv    !   formulation of the momentum advection
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90	(revision 2905)
@@ -249,5 +249,5 @@
             
             !  energy needed to bring ocean surface layer until its freezing
-            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1)   &
+            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj)   &
                 &          * ( tfu(ji,jj) - sst_m(ji,jj) ) * ( 1 - zinda )
             
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 2905)
@@ -139,18 +139,26 @@
    !! All coordinates
    !! ---------------
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdep3w_1           !: depth of T-points (sum of e3w) (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_1, gdepw_1   !: analytical depth at T-W  points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3v_1  , e3f_1     !: analytical vertical scale factors at  V--F
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_1  , e3u_1     !:                                       T--U  points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3vw_1             !: analytical vertical scale factors at  VW--
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3w_1  , e3uw_1    !:                                       W--UW  points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_b              !: before         -      -      -    -   T      points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3u_b  , e3v_b     !:   -            -      -      -    -   U--V   points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdep3w_n           !: now depth of T-points (sum of e3w) (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_n, gdepw_n   !: now depth at T-W  points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_n              !: now    vertical scale factors at  T       point  (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3u_n  , e3v_n     !:            -      -      -    -   U --V   points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3w_n  , e3f_n     !:            -      -      -    -   w --F   points  (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3uw_n , e3vw_n    !:            -      -      -    -   UW--VW  points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_b              !: before     -      -      -    -   T       points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3u_b  , e3v_b     !:   -        -      -      -    -   U --V   points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3uw_b , e3vw_b    !:   -        -      -      -    -   UW--VW  points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_a              !: after      -      -      -    -   T       point  (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3u_a  , e3v_a     !:   -        -      -      -    -   U --V   points (m)
 #else
    LOGICAL, PUBLIC, PARAMETER ::   lk_vvl = .FALSE.   !: fixed grid flag
 #endif
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hur  , hvr    !: inverse of u and v-points ocean depth (1/m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu   , hv     !: depth at u- and v-points (meters)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0 , hv_0   !: refernce depth at u- and v-points (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hur    , hvr       !: inverse of u and v-points ocean depth (1/m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu     , hv        !: depth at u- and v-points (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0               !: refernce depth at t-       points (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0   , hv_0      !: refernce depth at u- and v-points (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e1ur   , e2vr      !: scale factor coeffs at U--V points
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e12t   , e12t_1    !: horizontal cell surface and inverse  at T    points
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e12u   , e12u_1    !: horizontal cell surface and inverse  at U    points
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e12v_1 , e12f_1    !: inverse horizontal cell surface      at V--F   points
 
    INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1)
@@ -277,12 +285,16 @@
          !
 #if defined key_vvl
-      ALLOCATE( gdep3w_1(jpi,jpj,jpk) , e3v_1(jpi,jpj,jpk) , e3f_1 (jpi,jpj,jpk) ,                           &
-         &      gdept_1 (jpi,jpj,jpk) , e3t_1(jpi,jpj,jpk) , e3u_1 (jpi,jpj,jpk) ,                           &
-         &      gdepw_1 (jpi,jpj,jpk) , e3w_1(jpi,jpj,jpk) , e3vw_1(jpi,jpj,jpk) , e3uw_1(jpi,jpj,jpk) ,     &
-         &      e3t_b   (jpi,jpj,jpk) , e3u_b(jpi,jpj,jpk) , e3v_b (jpi,jpj,jpk)                       , STAT=ierr(5) )
-#endif
-         !
-      ALLOCATE( hu(jpi,jpj) , hur(jpi,jpj) , hu_0(jpi,jpj) ,     &
-         &      hv(jpi,jpj) , hvr(jpi,jpj) , hv_0(jpi,jpj) , STAT=ierr(6) )
+      ALLOCATE( gdep3w_n(jpi,jpj,jpk) , e3t_n (jpi,jpj,jpk) , e3u_n (jpi,jpj,jpk) ,                           &
+         &      gdept_n (jpi,jpj,jpk) , e3v_n (jpi,jpj,jpk) , e3w_n (jpi,jpj,jpk) ,                           &
+         &      gdepw_n (jpi,jpj,jpk) , e3f_n (jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) , e3uw_n(jpi,jpj,jpk) ,     &
+         &      e3t_b   (jpi,jpj,jpk) , e3u_b (jpi,jpj,jpk) , e3v_b (jpi,jpj,jpk) ,                           &
+         &      e3uw_b  (jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) ,                                                 &
+         &      e3t_a   (jpi,jpj,jpk) , e3u_a (jpi,jpj,jpk) , e3v_a (jpi,jpj,jpk)                       , STAT=ierr(5) )
+#endif
+         !
+      ALLOCATE( hu  (jpi,jpj) , hur   (jpi,jpj) , hu_0  (jpi,jpj) , ht_0  (jpi,jpj) ,     &
+         &      hv  (jpi,jpj) , hvr   (jpi,jpj) , hv_0  (jpi,jpj) ,                       &
+         &      e1ur(jpi,jpj) , e2vr  (jpi,jpj) , e12t  (jpi,jpj) , e12t_1(jpi,jpj) ,     &
+         &      e12u(jpi,jpj) , e12u_1(jpi,jpj) , e12v_1(jpi,jpj) , e12f_1(jpi,jpj) , STAT=ierr(6)  )
          !
       ALLOCATE( gdept_0(jpk) , gdepw_0(jpk) ,                                     &
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 2905)
@@ -1,2 +1,3 @@
+
 MODULE domain
    !!==============================================================================
@@ -81,5 +82,5 @@
                              CALL dom_zgr      ! Vertical mesh and bathymetry
                              CALL dom_msk      ! Masks
-      IF( lk_vvl         )   CALL dom_vvl      ! Vertical variable mesh
+      IF( lk_vvl )           CALL dom_vvl_init ! Vertical variable mesh
       !
       IF( lk_c1d ) THEN                        ! 1D configuration 
@@ -89,9 +90,19 @@
       END IF
       !
+      ! - ML - only used in domvvl but could be usefull in many routines
+      e12t  (:,:) = e1t(:,:) * e2t(:,:)
+      e12u  (:,:) = e1u(:,:) * e2u(:,:)
+      e12u_1(:,:) = 0.5 / e12u(:,:)
+      e12v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) )
+      e12f_1(:,:) = 0.5 / ( e1f(:,:) * e2f(:,:) )
+      ! - ML - used in domvvl and traldf_(lab/bilap/iso)
+      e1ur  (:,:) = e2u(:,:) / e1u(:,:)
+      e2vr  (:,:) = e1v(:,:) / e2v(:,:)
+      !
       hu(:,:) = 0.e0                           ! Ocean depth at U- and V-points
       hv(:,:) = 0.e0
       DO jk = 1, jpk
-         hu(:,:) = hu(:,:) + fse3u(:,:,jk) * umask(:,:,jk)
-         hv(:,:) = hv(:,:) + fse3v(:,:,jk) * vmask(:,:,jk)
+         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
+         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
       END DO
       !                                        ! Inverse of the local depth
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 2905)
@@ -6,16 +6,25 @@
    !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code
    !!            3.1  !  2009-02  (G. Madec, M. Leclair, R. Benshila)  pure z* coordinate
-   !!----------------------------------------------------------------------
+   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl:
+   !!                                          vvl option includes z_star and z_tilde coordinates
 #if defined key_vvl
    !!----------------------------------------------------------------------
    !!   'key_vvl'                              variable volume
    !!----------------------------------------------------------------------
-   !!   dom_vvl     : defined coefficients to distribute ssh on each layers
    !!----------------------------------------------------------------------
+   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
+   !!   dom_vvl_nxt      : Compute next vertical scale factors
+   !!   dom_vvl_swp      : Swap vertical scale factors and update the vertical grid
+   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
+   !!   dom_vvl_rst      : read/write restart file
+   !!   dom_vvl_ctl      : Check the vvl options
+   !!----------------------------------------------------------------------
+   !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
-   USE sbc_oce         ! surface boundary condition: ocean
-   USE phycst          ! physical constants
+   USE sbc_oce         ! ocean surface boundary condition
    USE in_out_manager  ! I/O manager
+   USE iom             ! I/O manager library
+   USE restart         ! ocean restart
    USE lib_mpp         ! distributed memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -24,12 +33,26 @@
    PRIVATE
 
-   PUBLIC   dom_vvl       ! called by domain.F90
-   PUBLIC   dom_vvl_alloc ! called by nemogcm.F90
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ee_t, ee_u, ee_v, ee_f   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mut , muu , muv , muf    !: ??? 
-
-   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:)     ::   r2dt   ! vertical profile time-step, = 2 rdttra 
-      !                                                              ! except at nit000 (=rdttra) if neuler=0
+   !! * Routine accessibility
+   PUBLIC dom_vvl_init       ! called by domain.F90
+   PUBLIC dom_vvl_sf_nxt     ! called by step.F90
+   PUBLIC dom_vvl_sf_swp     ! called by step.F90
+   PUBLIC dom_vvl_interpol   ! called by dynnxt.F90
+
+   !!* Namelist nam_vvl
+   LOGICAL , PUBLIC                                      :: ln_vvl_zstar  = .TRUE.    ! zstar  vertical coordinate
+   LOGICAL , PUBLIC                                      :: ln_vvl_ztilde = .FALSE.   ! ztilde vertical coordinate
+   LOGICAL , PUBLIC                                      :: ln_vvl_layer  = .FALSE.   ! level  vertical coordinate
+   LOGICAL , PUBLIC                                      :: ln_vvl_kepe   = .FALSE.   ! kinetic/potential energy transfer
+   !                                                                                  ! conservation: not used yet
+   REAL(wp), PUBLIC                                      :: ahe3          =  0.e0     ! thickness diffusion coefficient
+
+   !! * Module variables
+   INTEGER                                               :: nvvl                      ! choice of vertical coordinate
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td              ! thickness diffusion transport
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                   ! low frequency part of hz divergence
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_t_b, e3t_t_n, e3t_t_a ! baroclinic scale factors
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_restore_e3t           ! retoring period for scale factors
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_restore_hdv           ! retoring period for low freq. divergence
+
 
    !! * Substitutions
@@ -37,167 +60,750 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! NEMO/OPA 3.3 , NEMO-Consortium (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-CONTAINS       
+
+CONTAINS
 
    INTEGER FUNCTION dom_vvl_alloc()
       !!----------------------------------------------------------------------
-      !!                ***  ROUTINE dom_vvl_alloc  ***
-      !!----------------------------------------------------------------------
-      !
-      ALLOCATE( mut (jpi,jpj,jpk) , muu (jpi,jpj,jpk) , muv (jpi,jpj,jpk) , muf (jpi,jpj,jpk) ,     &
-         &      ee_t(jpi,jpj)     , ee_u(jpi,jpj)     , ee_v(jpi,jpj)     , ee_f(jpi,jpj)     ,     &
-         &      r2dt        (jpk)                                                             , STAT=dom_vvl_alloc )
-         !
-      IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
-      IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
-      !
+      !!                ***  FUNCTION dom_vvl_alloc  ***
+      !!----------------------------------------------------------------------
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
+         ALLOCATE(                                                                                             &
+            &      un_td  (jpi,jpj,jpk) , vn_td  (jpi,jpj,jpk) , hdiv_lf(jpi,jpj,jpk) ,                        &
+            &      e3t_t_b(jpi,jpj,jpk) , e3t_t_n(jpi,jpj,jpk) , e3t_t_a(jpi,jpj,jpk) ,                        &
+            &      frq_restore_e3t(jpi,jpj) , frq_restore_hdv(jpi,jpj)                , STAT= dom_vvl_alloc )
+         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
+         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
+      ENDIF
+
    END FUNCTION dom_vvl_alloc
 
 
-   SUBROUTINE dom_vvl
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE dom_vvl  ***
+   SUBROUTINE dom_vvl_init
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_init  ***
       !!                   
-      !! ** Purpose :  compute coefficients muX at T-U-V-F points to spread
-      !!               ssh over the whole water column (scale factors)
-      !!----------------------------------------------------------------------
-      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3     ! 2D workspace
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zcoefu , zcoefv   , zcoeff                   ! local scalars
-      REAL(wp) ::   zv_t_ij, zv_t_ip1j, zv_t_ijp1, zv_t_ip1jp1   !   -      -
-      !!----------------------------------------------------------------------
-
-      IF( wrk_in_use(2, 1,2,3) ) THEN
-         CALL ctl_stop('dom_vvl: requested workspace arrays unavailable')   ;   RETURN
-      ENDIF
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'dom_vvl : Variable volume initialization'
-         WRITE(numout,*) '~~~~~~~~  compute coef. used to spread ssh over each layers'
-      ENDIF
-      
-      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl : unable to allocate arrays' )
-
-      fsdept(:,:,:) = gdept (:,:,:)
-      fsdepw(:,:,:) = gdepw (:,:,:)
-      fsde3w(:,:,:) = gdep3w(:,:,:)
-      fse3t (:,:,:) = e3t   (:,:,:)
-      fse3u (:,:,:) = e3u   (:,:,:)
-      fse3v (:,:,:) = e3v   (:,:,:)
-      fse3f (:,:,:) = e3f   (:,:,:)
-      fse3w (:,:,:) = e3w   (:,:,:)
-      fse3uw(:,:,:) = e3uw  (:,:,:)
-      fse3vw(:,:,:) = e3vw  (:,:,:)
-
-      !                                 !==  mu computation  ==!
-      ee_t(:,:) = fse3t_0(:,:,1)                ! Lower bound : thickness of the first model level
-      ee_u(:,:) = fse3u_0(:,:,1)
-      ee_v(:,:) = fse3v_0(:,:,1)
-      ee_f(:,:) = fse3f_0(:,:,1)
-      DO jk = 2, jpkm1                          ! Sum of the masked vertical scale factors
-         ee_t(:,:) = ee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk)
-         ee_u(:,:) = ee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk)
-         ee_v(:,:) = ee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk)
-         DO jj = 1, jpjm1                      ! f-point : fmask=shlat at coasts, use the product of umask
-            ee_f(:,jj) = ee_f(:,jj) + fse3f_0(:,jj,jk) *  umask(:,jj,jk) * umask(:,jj+1,jk)
-         END DO
-      END DO  
-      !                                         ! Compute and mask the inverse of the local depth at T, U, V and F points
-      ee_t(:,:) = 1. / ee_t(:,:) * tmask(:,:,1)
-      ee_u(:,:) = 1. / ee_u(:,:) * umask(:,:,1)
-      ee_v(:,:) = 1. / ee_v(:,:) * vmask(:,:,1)
-      DO jj = 1, jpjm1                               ! f-point case fmask cannot be used 
-         ee_f(:,jj) = 1. / ee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1)
+      !! ** Purpose :  Initialization of all scale factors, depths
+      !!               and water column heights
+      !!
+      !! ** Method  :  - use restart file and/or initialize
+      !!               - interpolate scale factors
+      !!
+      !! ** Action  : - fse3t_(n/b) and e3t_t_(n/b)
+      !!              - Regrid: fse3(u/v)_n
+      !!                        fse3(u/v)_b       
+      !!                        fse3w_n           
+      !!                        fse3(u/v)w_b      
+      !!                        fse3(u/v)w_n      
+      !!                        fsdept_n, fsdepw_n and fsde3w_n
+      !!                        h(u/v) and h(u/v)r
+      !!              - ht_0 and ht1_0
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      !! * Local declarations
+      INTEGER ::   jk
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
+
+#if defined key_zco
+      CALL ctl_stop( 'dom_vvl_init : options key_zco/key_dynspg_rl are incompatible with variable volume option key_vvl')
+#endif
+
+      ! choose vertical coordinate (z_star, z_tilde or layer)
+      ! ==========================
+      CALL dom_vvl_ctl
+
+      ! Allocate module arrays
+      ! ======================
+      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
+
+      ! read or initialize e3t_t_(b/n) and fse3t_(b/n) 
+      ! ==============================================
+      CALL dom_vvl_rst( nit000, 'READ' )
+
+      ! Reconstruction of all vertical scale factors at now and before time steps
+      ! =========================================================================
+      ! Horizontal scale factor interpolations
+      ! --------------------------------------
+      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
+      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
+      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
+      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
+      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
+      ! Vertical scale factor interpolations
+      ! ------------------------------------
+      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
+      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
+      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
+      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
+      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
+      ! t- and w- points depth
+      ! ----------------------
+      fsdept_n(:,:,1) = 0.5 * fse3w_n(:,:,1)
+      fsdepw_n(:,:,1) = 0.e0
+      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
+      DO jk = 2, jpk
+         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
+         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
+         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
       END DO
-      CALL lbc_lnk( ee_f, 'F', 1. )                  ! lateral boundary condition on ee_f
-      !
-      DO jk = 1, jpk                            ! mu coefficients
-         mut(:,:,jk) = ee_t(:,:) * tmask(:,:,jk)     ! T-point at T levels
-         muu(:,:,jk) = ee_u(:,:) * umask(:,:,jk)     ! U-point at T levels
-         muv(:,:,jk) = ee_v(:,:) * vmask(:,:,jk)     ! V-point at T levels
-      END DO
-      DO jk = 1, jpk                                 ! F-point : fmask=shlat at coasts, use the product of umask
-         DO jj = 1, jpjm1
-               muf(:,jj,jk) = ee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk)   ! at T levels
-         END DO
-         muf(:,jpj,jk) = 0.e0
-      END DO
-      CALL lbc_lnk( muf, 'F', 1. )                   ! lateral boundary condition
-
-
-      hu_0(:,:) = 0.e0                          ! Reference ocean depth at U- and V-points
+      ! Reference water column height at t-, u- and v- point
+      ! ----------------------------------------------------
+      ht_0(:,:) = 0.e0
+      hu_0(:,:) = 0.e0
       hv_0(:,:) = 0.e0
       DO jk = 1, jpk
+         ht_0(:,:) = ht_0(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk)
          hu_0(:,:) = hu_0(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk)
          hv_0(:,:) = hv_0(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk)
       END DO
-      
-      ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
-      ! for ssh and scale factors
-      zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
-      zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:)
-      zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:)
-
-      DO jj = 1, jpjm1                          ! initialise before and now Sea Surface Height at u-, v-, f-points
-         DO ji = 1, jpim1   ! NO vector opt.
-            zcoefu = umask(ji,jj,1) * zs_u_1(ji,jj)
-            zcoefv = vmask(ji,jj,1) * zs_v_1(ji,jj)
-            zcoeff = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) / ( e1f(ji,jj) * e2f(ji,jj) )
-            ! before fields
-            zv_t_ij       = zs_t(ji  ,jj  ) * sshb(ji  ,jj  )
-            zv_t_ip1j     = zs_t(ji+1,jj  ) * sshb(ji+1,jj  )
-            zv_t_ijp1     = zs_t(ji  ,jj+1) * sshb(ji  ,jj+1)
-            sshu_b(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j )
-            sshv_b(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 )
-            ! now fields
-            zv_t_ij       = zs_t(ji  ,jj  ) * sshn(ji  ,jj  )
-            zv_t_ip1j     = zs_t(ji+1,jj  ) * sshn(ji+1,jj  )
-            zv_t_ijp1     = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1)
-            zv_t_ip1jp1   = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1)
-            sshu_n(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j )
-            sshv_n(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 )
-            sshf_n(ji,jj) = zcoeff * ( zv_t_ij + zv_t_ip1j + zv_t_ijp1 + zv_t_ip1jp1 )
-         END DO
+
+   END SUBROUTINE dom_vvl_init
+
+
+   SUBROUTINE dom_vvl_sf_nxt( kt ) 
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
+      !!                   
+      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
+      !!                 tranxt and dynspg routines
+      !!
+      !! ** Method  :  - z_tilde_case: after scale factor increment computed with
+      !!                               high frequency part of horizontal divergence + retsoring to
+      !!                               towards the background grid + thickness difusion.
+      !!                               Then repartition of ssh INCREMENT proportionnaly
+      !!                               to the "baroclinic" level thickness.
+      !!               - z_star case:  Repartition of ssh proportionnaly to the level thickness.
+      !!
+      !! ** Action  :  - hdiv_lf: restoring towards full baroclinic divergence in z_tilde case
+      !!               - e3t_t_a: after increment of vertical scale factor 
+      !!                          in z_tilde case
+      !!               - fse3(t/u/v)_a
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released, iwrk_in_use, iwrk_not_released
+      USE oce     , ONLY:   ze3t    => ta         ! ua used as workspace
+      USE wrk_nemo, ONLY:   zht     => wrk_2d_1   ! 2D REAL workspace
+      USE wrk_nemo, ONLY:   z_scale => wrk_2d_2   ! 2D REAL workspace
+      USE wrk_nemo, ONLY:   zwu     => wrk_2d_3   ! 2D REAL workspace
+      USE wrk_nemo, ONLY:   zwv     => wrk_2d_4   ! 2D REAL workspace
+      USE wrk_nemo, ONLY:   zhdiv   => wrk_2d_5   ! 2D REAL workspace
+      !! * Arguments
+      INTEGER, INTENT( in )                  :: kt                    ! time step
+      !! * Local declarations
+      INTEGER                                :: ji, jj, jk            ! dummy loop indices
+      INTEGER , DIMENSION(3)                 :: ijk_max, ijk_min      ! temporary integers
+      REAL(wp)                               :: z2dt                  ! temporary scalars
+      REAL(wp)                               :: z_def_max             ! temporary scalar
+      REAL(wp)                               :: z_tmin, z_tmax        ! temporary scalars
+      LOGICAL                                :: ln_debug = .FALSE.    ! local logical for debug prints
+      !!----------------------------------------------------------------------
+      IF( wrk_in_use(2, 1, 2, 3, 4, 5) ) THEN
+         CALL ctl_stop('dom_vvl_sf_nxt: requested workspace arrays unavailable')   ;   RETURN
+      END IF
+
+      IF(kt == nit000)   THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
+      ENDIF
+
+      ! ******************************* !
+      ! After acale factors at t-points !
+      ! ******************************* !
+      !                                             !----------------------------!
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! ZTILDE or LAYER coordinate !
+         !                                          !----------------------------!
+
+         ! I - Initialisations
+         ! ===================
+         IF( kt == nit000 ) THEN
+            ! - ML - In the future, this should be tunable by the user
+            IF( ln_vvl_ztilde ) THEN       ! ZTILDE CASE
+               ! DO jj = 1, jpj
+               !    DO ji = 1, jpi
+               !       frq_restore_hdv(ji,jj) = 2.e0 * rpi / 86400.e0 / 5.e0   &
+               !          &                     * MAX( SIN( ABS( gphit(ji,jj) ) / .5e0, 1.e0 / 6.e0) )
+               !    END DO
+               ! END DO
+               ! frq_restore_e3t(:,:) = ( 1.e0 / 6.e0 ) * frq_restore_hdv(:,:)
+               frq_restore_e3t(:,:) = 2.e0 * rpi / ( 30.e0 * 86400.e0 )
+               frq_restore_hdv(:,:) = 2.e0 * rpi / (  5.e0 * 86400.e0 )
+               ! frq_restore_hdv(:,:) = 2.e0 * rpi / (  2.e0 * 86400.e0 )
+            ELSE                           ! LAYER CASE
+               frq_restore_e3t(:,:) = 0.e0
+               frq_restore_hdv(:,:) = 0.e0
+            ENDIF
+
+         ENDIF
+
+         ! II - Low frequency horizontal divergence
+         ! ========================================
+
+         ! 1 - barotropic divergence
+         ! -------------------------
+         zhdiv(:,:) = 0.
+         zht(:,:)   = 0.
+         DO jk = 1, jpkm1
+            zhdiv(:,:) = zhdiv(:,:) + hdivn(:,:,jk) * fse3t_n(:,:,jk)
+            zht  (:,:) = zht  (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+         END DO
+         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask(:,:,1) )
+
+         ! 2 - restoring equation
+         ! ----------------------
+         IF( kt .GT. nit000 ) THEN
+            DO jk = 1, jpkm1
+               hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_restore_hdv(:,:)   &
+                  &          * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
+            END DO
+         ENDIF
+
+         ! III - after z_tilde increments of vertical scale factors
+         ! =========================================================
+         e3t_t_a(:,:,:) = 0.e0   ! e3t_t_a used to store tendency terms
+
+         ! 1 - High frequency divergence term
+         ! ----------------------------------
+         DO jk = 1, jpkm1
+            e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
+         END DO
+
+         ! 2 - Restoring term
+         ! ------------------
+         DO jk = 1, jpk
+            e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - frq_restore_e3t(:,:) * e3t_t_b(:,:,jk)
+         END DO
+
+         ! 3 - Thickness diffusion term
+         ! ----------------------------
+         zwu(:,:) = 0.e0
+         zwv(:,:) = 0.e0
+
+         ! a - first derivative: diffusive fluxes
+         DO jk = 1, jpkm1
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  un_td(ji,jj,jk) = ahe3 * umask(ji,jj,jk) * e1ur(ji,jj) * ( e3t_t_b(ji,jj,jk) - e3t_t_b(ji+1,jj  ,jk) )
+                  vn_td(ji,jj,jk) = ahe3 * vmask(ji,jj,jk) * e2vr(ji,jj) * ( e3t_t_b(ji,jj,jk) - e3t_t_b(ji  ,jj+1,jk) )
+                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
+                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         ! b - correction for last oceanic u-v points
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
+               vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
+            END DO
+         END DO
+         ! c - second derivative: divergence of diffusive fluxes
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  e3t_t_a(ji,jj,jk) = e3t_t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
+                     &                                      + vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk) )  &
+                     &                                    * e12t_1(ji,jj)
+               END DO
+            END DO
+         END DO
+         ! d - thickness diffusion equivalent transport: boundary conditions
+         !     (stored for tracer advaction and continuity equation)
+         CALL lbc_lnk( un_td , 'U' , -1.)
+         CALL lbc_lnk( vn_td , 'V' , -1.)
+
+         ! 4 - Time stepping of baroclinic scale factors
+         ! ---------------------------------------------
+         ! Leapfrog time stepping
+         ! ~~~~~~~~~~~~~~~~~~~~~~
+         IF( neuler == 0 .AND. kt == nit000 ) THEN
+            z2dt =  rdt
+         ELSE
+            z2dt = 2.e0 * rdt
+         ENDIF
+         CALL lbc_lnk( e3t_t_a(:,:,:), 'T', 1. )    ! - ML - Do not use lbc_lnk_e3: e3t_t_a is a tendency term at this stage
+         e3t_t_a(:,:,:) = e3t_t_b(:,:,:) + z2dt * tmask(:,:,:) * e3t_t_a(:,:,:)
+         fse3t_a(:,:,:) = fse3t_0(:,:,:) + e3t_t_a(:,:,:)
+
+         ! Maximum deformation control
+         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         ! - ML - Should perhaps be put in the namelist
+         z_def_max = 0.9e0
+         ze3t(:,:,jpk) = 0.e0
+         DO jk = 1, jpkm1
+            ze3t(:,:,jk) = e3t_t_a(:,:,jk) / fse3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
+         END DO
+         z_tmax = MAXVAL( ze3t(:,:,:) )
+         z_tmin = MINVAL( ze3t(:,:,:) )
+         ! - ML - test: for the moment, stop simulation for too large e3_t variations
+         IF( ( z_tmax .GT. z_def_max ) .OR. ( z_tmin .LT. - z_def_max ) ) THEN
+            ijk_max = MAXLOC( ze3t(:,:,:) )
+            ijk_min = MINLOC( ze3t(:,:,:) )
+            WRITE(numout, *) 'MAX( e3t_t_a(:,:,:) / fse3t_0(:,:,:) ) =', z_tmax
+            WRITE(numout, *) 'at i, j, k=', ijk_max
+            WRITE(numout, *) 'MIN( e3t_t_a(:,:,:) / fse3t_0(:,:,:) ) =', z_tmin
+            WRITE(numout, *) 'at i, j, k=', ijk_min            
+            CALL ctl_stop('MAX( ABS( e3t_t_a(:,:,:) ) / fse3t_0(:,:,:) ) too high')
+         ENDIF
+         ! - ML - end test
+         ! - ML - This will cause a baroclinicity error if the ctl_stop above is not used
+         e3t_t_a(:,:,:) = MIN( e3t_t_a(:,:,:), ( 1.e0 + z_def_max ) * fse3t_0(:,:,:) )
+         e3t_t_a(:,:,:) = MAX( e3t_t_a(:,:,:), ( 1.e0 - z_def_max ) * fse3t_0(:,:,:) )
+
+         ! Boundary conditions
+         ! ~~~~~~~~~~~~~~~~~~~
+         ! - ML - think it's not necessary
+         ! CALL lbc_lnk( e3t_t_a(:,:,:), 'T', 1. )    ! - ML - Do not use lbc_lnk_e3: e3t_t_a is a level thickness ANOMALY
+
+         ! IV - Barotropic repartition of the sea surface height over the baroclinic profile
+         ! =================================================================================
+         ! add e3t(n-1) "star" Asselin-filtered
+         DO jk = 1, jpkm1
+            ! - ML - : multiplication by tmask not necessary a priori. Just to be sure for the moment.
+            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + ( fse3t_b(:,:,jk) - fse3t_0(:,:,jk) - e3t_t_b(:,:,jk) ) * tmask(:,:,jk)
+         END DO
+         ! add ( ssh increment + "baroclinicity error" ) proportionnaly to e3t(n)
+         ! - ML - baroclinicity error should be better treated in the future
+         !        i.e. locally and not spread over the water column.
+         zht(:,:) = 0.
+         DO jk = 1, jpkm1
+            zht(:,:) = zht(:,:) + e3t_t_a(:,:,jk) * tmask(:,:,jk)
+         END DO
+         z_scale(:,:) = ( ssha(:,:) - sshb(:,:) - zht(:,:) ) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
+         DO jk = 1, jpkm1
+            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+         END DO
+
+         !                                !------------------!
+      ELSEIF( ln_vvl_zstar ) THEN         ! Zstar coordinate !
+         !                                !------------------!
+         z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
+         DO jk = 1, jpkm1
+            fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:)
+         END DO
+
+      ENDIF
+
+      IF( ln_debug ) THEN   ! - ML - test: control prints for debuging
+         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
+            WRITE(numout, *) 'kt  =', kt
+            WRITE(numout, *) 'MAXVAL(abs(ht_0-SUM(e3t_t_a))) =',   &
+               &              MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) - zht(:,:) ) )
+         END IF
+         zht(:,:) = 0.e0
+         DO jk = 1, jpkm1
+            zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+         END DO
+         WRITE(numout, *) 'MAXVAL(abs(ht_0+sshn-SUM(fse3t_n))) =',   &
+            &              MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
+         zht(:,:) = 0.e0
+         DO jk = 1, jpkm1
+            zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)
+         END DO
+         WRITE(numout, *) 'MAXVAL(abs(ht_0+sshn-SUM(fse3t_a))) =',   &
+            &              MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
+      END IF
+
+      ! *********************************** !
+      ! After scale factors at u- v- points !
+      ! *********************************** !
+
+      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' )
+      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' )
+
+      IF( wrk_not_released(2, 1, 2, 3, 4, 5) ) THEN
+         CALL ctl_stop( 'dom_vvl_sf_nxt: failed to release workspace arrays' )
+      ENDIF
+
+   END SUBROUTINE dom_vvl_sf_nxt
+
+
+   SUBROUTINE dom_vvl_sf_swp( kt )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_sf_swp  ***
+      !!                   
+      !! ** Purpose :  compute time filter and swap of scale factors 
+      !!               compute all depths and related variables for next time step
+      !!               write outputs and restart file
+      !!
+      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
+      !!               - reconstruct scale factor at other grid points (interpolate)
+      !!               - recompute depths and water height fields
+      !!
+      !! ** Action  :  - fse3t_(b/n), e3t_t_(b/n) and fse3(u/v)_n ready for next time step
+      !!               - Recompute:
+      !!                    fse3(u/v)_b       
+      !!                    fse3w_n           
+      !!                    fse3(u/v)w_b      
+      !!                    fse3(u/v)w_n      
+      !!                    fsdept_n, fsdepw_n  and fsde3w_n
+      !!                    h(u/v) and h(u/v)r
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
+      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      USE oce, ONLY:   z_e3t_def => ta                ! square of baroclinic scale factors deformation (in %)
+      !! * Arguments
+      INTEGER, INTENT( in )               :: kt       ! time step
+      !! * Local declarations
+      INTEGER                             :: jk       ! dummy loop indices
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 )   THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
+      ENDIF
+
+      ! Time filter and swap of scale factors
+      ! =====================================
+      ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
+         IF( neuler == 0 .AND. kt == nit000 ) THEN
+            e3t_t_b(:,:,:) = e3t_t_n(:,:,:)
+         ELSE
+            e3t_t_b(:,:,:) = e3t_t_n(:,:,:) + atfp * (e3t_t_b(:,:,:) - 2.e0 * e3t_t_n(:,:,:) + e3t_t_a(:,:,:) )
+         ENDIF
+         e3t_t_n(:,:,:) = e3t_t_a(:,:,:)
+      ENDIF
+      fse3t_n(:,:,:) = fse3t_a(:,:,:)
+      fse3u_n(:,:,:) = fse3u_a(:,:,:)
+      fse3v_n(:,:,:) = fse3v_a(:,:,:)
+
+      ! Compute all missing vertical scale factor and depths
+      ! ====================================================
+      ! Horizontal scale factor interpolations
+      ! --------------------------------------
+      ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt
+      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F'  )
+      ! Vertical scale factor interpolations
+      ! ------------------------------------
+      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
+      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
+      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
+      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
+      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
+      ! t- and w- points depth
+      ! ----------------------
+      fsdept_n(:,:,1) = 0.5 * fse3w_n(:,:,1)
+      fsdepw_n(:,:,1) = 0.e0
+      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
+      DO jk = 2, jpk
+         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
+         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
+         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
       END DO
-      CALL lbc_lnk( sshu_n, 'U', 1. )   ;   CALL lbc_lnk( sshu_b, 'U', 1. )      ! lateral boundary conditions
-      CALL lbc_lnk( sshv_n, 'V', 1. )   ;   CALL lbc_lnk( sshv_b, 'V', 1. )
-      CALL lbc_lnk( sshf_n, 'F', 1. )
-
-                                                ! initialise before scale factors at (u/v)-points
-      ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
-      DO jk = 1, jpkm1
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
-               zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
-               zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
-               fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
-               fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
+      ! Local depth and Inverse of the local depth of the water column at u- and v- points
+      ! ----------------------------------------------------------------------------------
+      hu(:,:) = 0.
+      hv(:,:) = 0.
+      DO jk = 1, jpk
+         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
+         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
+      END DO
+      ! Inverse of the local depth
+      hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1. - umask(:,:,1) )
+      hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1. - umask(:,:,1) )
+
+      ! Write outputs
+      ! =============
+      ! - ML - add output variables in xml file for all configurations
+      z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - fse3t_0(:,:,:) ) / fse3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
+      CALL iom_put( "e3t_n"  , fse3t_n  (:,:,:) )
+      CALL iom_put( "dept"   , fsdept_n (:,:,:) )
+      CALL iom_put( "e3tdef" , z_e3t_def(:,:,:) )
+
+      ! write restart file
+      ! ==================
+      IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
+
+   END SUBROUTINE dom_vvl_sf_swp
+
+
+   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_vvl__interpol  ***
+      !!
+      !! ** Purpose :   interpolate scale factors from one grid point to another
+      !!
+      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
+      !!                - horizontal interpolation: grid cell surface averaging
+      !!                - vertical interpolation: simple averaging
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
+      CHARACTER(len=1), INTENT( in )                    ::  pout       ! grid point of out scale factors
+      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
+      !! * Local declarations
+      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
+      !!----------------------------------------------------------------------
+      SELECT CASE ( pout )
+         !               ! ------------------------------------- !
+      CASE( 'U' )        ! interpolation from T-point to U-point !
+         !               ! ------------------------------------- !
+         ! horizontal surface weighted interpolation
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  pe3_out(ji,jj,jk) = umask(ji,jj,jk) * e12u_1(ji,jj)                                          &
+                     &                  * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - fse3t_0(ji  ,jj,jk) )     &
+                     &                      + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - fse3t_0(ji+1,jj,jk) ) )
+               END DO
             END DO
          END DO
-      END DO
-      CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )               ! lateral boundary conditions
-      CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
-      ! Add initial scale factor to scale factor anomaly
-      fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
-      fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
+         ! boundary conditions
+         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1. )
+         pe3_out(:,:,:) = pe3_out(:,:,:) + fse3u_0(:,:,:)
+         !               ! ------------------------------------- !
+      CASE( 'V' )        ! interpolation from T-point to V-point !
+         !               ! ------------------------------------- !
+         ! horizontal surface weighted interpolation
+         DO jk = 1, jpkm1
+            DO jj = 1, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  pe3_out(ji,jj,jk) = umask(ji,jj,jk) * e12v_1(ji,jj)                                          &
+                     &                  * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - fse3t_0(ji,jj  ,jk) )     &
+                     &                      + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - fse3t_0(ji,jj+1,jk) ) )
+               END DO
+            END DO
+         END DO
+         ! boundary conditions
+         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1. )
+         pe3_out(:,:,:) = pe3_out(:,:,:) + fse3v_0(:,:,:)
+         !               ! ------------------------------------- !
+      CASE( 'F' )        ! interpolation from U-point to F-point !
+         !               ! ------------------------------------- !
+         ! horizontal surface weighted interpolation
+         DO jk = 1, jpkm1
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  pe3_out(ji,jj,jk) = umask(ji,jj,jk) * umask(ji,jj+1,jk) * e12f_1(ji,jj)                      &
+                     &                  * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - fse3u_0(ji,jj  ,jk) )     &
+                     &                      + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - fse3u_0(ji,jj+1,jk) ) )
+               END DO
+            END DO
+         END DO
+         ! boundary conditions
+         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1. )
+         pe3_out(:,:,:) = pe3_out(:,:,:) + fse3f_0(:,:,:)
+         !               ! ------------------------------------- !
+      CASE( 'W' )        ! interpolation from T-point to W-point !
+         !               ! ------------------------------------- !
+         ! vertical simple interpolation
+         pe3_out(:,:,1) = fse3w_0(:,:,1) + pe3_in(:,:,1) - fse3t_0(:,:,1)
+         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without undirect adressing
+         DO jk = 2, jpk
+            pe3_out(:,:,jk) = fse3w_0(:,:,jk) + ( 1. - 0.5 * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - fse3t_0(:,:,jk-1) )   &
+               &                              + (      0.5 * tmask(:,:,jk) ) * ( pe3_in(:,:,jk  ) - fse3t_0(:,:,jk  ) )
+         END DO
+         !               ! -------------------------------------- !
+      CASE( 'UW' )       ! interpolation from U-point to UW-point !
+         !               ! -------------------------------------- !
+         ! vertical simple interpolation
+         pe3_out(:,:,1) = fse3uw_0(:,:,1) + pe3_in(:,:,1) - fse3u_0(:,:,1)
+         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without undirect adressing
+         DO jk = 2, jpk
+            pe3_out(:,:,jk) = fse3uw_0(:,:,jk) + ( 1. - 0.5 * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - fse3u_0(:,:,jk-1) )   &
+               &                               + (      0.5 * umask(:,:,jk) ) * ( pe3_in(:,:,jk  ) - fse3u_0(:,:,jk  ) )
+         END DO
+         !               ! -------------------------------------- !
+      CASE( 'VW' )       ! interpolation from V-point to VW-point !
+         !               ! -------------------------------------- !
+         ! vertical simple interpolation
+         pe3_out(:,:,1) = fse3vw_0(:,:,1) + pe3_in(:,:,1) - fse3v_0(:,:,1)
+         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without undirect adressing
+         DO jk = 2, jpk
+            pe3_out(:,:,jk) = fse3vw_0(:,:,jk) + ( 1. - 0.5 * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - fse3v_0(:,:,jk-1) )   &
+               &                               + (      0.5 * vmask(:,:,jk) ) * ( pe3_in(:,:,jk  ) - fse3v_0(:,:,jk  ) )
+         END DO
+      END SELECT
+
+   END SUBROUTINE dom_vvl_interpol
+
+
+   SUBROUTINE dom_vvl_rst( kt, cdrw )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE dom_vvl_rst  ***
+      !!                     
+      !! ** Purpose :   Read or write VVL file in restart file
+      !!
+      !! ** Method  :   use of IOM library
+      !!                if the restart does not contain vertical scale factors,
+      !!                they are set to the _0 values
+      !!                if the restart does not contain vertical scale factors increments (z_tilde),
+      !!                they are set to 0.
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
+      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      !! * Local declarations
+      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
+      !!----------------------------------------------------------------------
       !
-      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dom_vvl: failed to release workspace arrays')
-      !
-   END SUBROUTINE dom_vvl
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+         !                                   ! ===============
+         IF( ln_rstart ) THEN                   !* Read the restart file
+            id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. )
+            id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. )
+            id3 = iom_varid( numror, 'e3t_t_b', ldstop = .FALSE. )
+            id4 = iom_varid( numror, 'e3t_t_n', ldstop = .FALSE. )
+            id5 = iom_varid( numror, 'hdif_lf', ldstop = .FALSE. )
+
+            !                         ! ----------- !
+            IF( ln_vvl_zstar ) THEN   ! z_star case !
+               !                      ! ----------- !
+               IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
+                  CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
+                  CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
+                  IF( neuler == 0 ) THEN
+                     fse3t_b(:,:,:) = fse3t_n(:,:,:)
+                  ENDIF
+               ELSE                                 ! one at least array is missing
+                  CALL ctl_stop( 'vvl cannot restart from a non vvl run' )
+               ENDIF
+               IF( MIN( id3, id4, id5 ) > 0 ) CALL ctl_stop( 'z_star coordinate cannot restart from a z_tilde run' )
+               !                      ! ------------ !
+            ELSE                      ! z_tilde case !
+               !                      ! ------------ !
+               IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
+                  CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
+                  CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
+                  IF( neuler == 0 ) THEN
+                     fse3t_b(:,:,:) = fse3t_n(:,:,:)
+                  ENDIF
+               ELSE                                 ! one at least array is missing
+                  CALL ctl_stop( 'vvl cannot restart from a non vvl run' )
+               ENDIF
+               IF( MIN( id3, id4, id5 ) > 0 ) THEN  ! all required arrays exist
+                  CALL iom_get( numror, jpdom_autoglo, 'e3t_t_b', e3t_t_b(:,:,:) )
+                  CALL iom_get( numror, jpdom_autoglo, 'e3t_t_n', e3t_t_n(:,:,:) )
+                  CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
+               ELSE                                 ! one at least array is missing
+                  e3t_t_b(:,:,:) = 0.e0
+                  e3t_t_n(:,:,:) = 0.e0
+                  hdiv_lf(:,:,:) = 0.e0
+               ENDIF
+
+            ENDIF
+
+         ELSE                                   !* Initialize at "rest"
+            fse3t_b(:,:,:) = fse3t_0(:,:,:)
+            fse3t_n(:,:,:) = fse3t_0(:,:,:)
+            e3t_t_b(:,:,:) = 0.e0
+            e3t_t_n(:,:,:) = 0.e0
+            hdiv_lf(:,:,:) = 0.e0
+         ENDIF
+
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+         !                                   ! ===================
+         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
+         !                                         ! ---------------------- !
+         !                                         ! z_star & z_tilde cases !
+         !                                         ! ---------------------- !
+         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
+         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )
+         !                                           ! ----------------- !
+         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde case only !
+            !                                        ! ----------------- !
+            CALL iom_rstput( kt, nitrst, numrow, 'e3t_t_b', e3t_t_b(:,:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'e3t_t_n', e3t_t_n(:,:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )
+         ENDIF
+
+      ENDIF
+
+   END SUBROUTINE dom_vvl_rst
+
+
+   SUBROUTINE dom_vvl_ctl
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_vvl_ctl  ***
+      !!                
+      !! ** Purpose :   Control the consistency between namelist options for 
+      !!                vertical coordinate and set nvvl
+      !!----------------------------------------------------------------------
+      INTEGER ::   ioptio
+
+      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ahe3! , ln_vvl_kepe
+      !!---------------------------------------------------------------------- 
+
+      REWIND ( numnam )               ! Read Namelist nam_vvl : vertical coordinate
+      READ   ( numnam, nam_vvl )
+
+      IF(lwp) THEN                    ! Namelist print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
+         WRITE(numout,*) '~~~~~~~~~~~'
+         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate'
+         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
+         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
+         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
+         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
+         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
+         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion coefficient'
+         WRITE(numout,*) '                                         ahe3           = ', ahe3
+      ENDIF
+
+      ioptio = 0                      ! Parameter control
+      IF( ln_vvl_zstar     )   ioptio = ioptio + 1
+      IF( ln_vvl_ztilde    )   ioptio = ioptio + 1
+      IF( ln_vvl_layer     )   ioptio = ioptio + 1
+
+      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
+
+      IF( ln_vvl_zstar  )   nvvl = 1
+      IF( ln_vvl_ztilde )   nvvl = 2
+      IF( ln_vvl_layer  )   nvvl = 3
+
+      IF(lwp) THEN                   ! Print the choice
+         WRITE(numout,*)
+         IF( nvvl ==  1        ) WRITE(numout,*) '              zstar vertical coordinate is used'
+         IF( nvvl ==  2        ) WRITE(numout,*) '              ztilde vertical coordinate is used'
+         IF( nvvl ==  3        ) WRITE(numout,*) '              layer vertical coordinate is used'
+         ! - ML - Option not developed yet
+         ! IF(       ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option used'
+         ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option not used'
+      ENDIF
+
+   END SUBROUTINE dom_vvl_ctl
+
 
 #else
+
+
    !!----------------------------------------------------------------------
    !!   Default option :                                      Empty routine
    !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE dom_vvl
-   END SUBROUTINE dom_vvl
+   SUBROUTINE dom_vvl_init
+      WRITE(*,*) 'dom_vvl_init: You should not have seen this print! error?'
+   END SUBROUTINE dom_vvl_init
+   SUBROUTINE dom_vvl_sf_nxt( kt ) 
+      !! * Arguments
+      INTEGER, INTENT( in )  ::    kt
+      WRITE(*,*) 'dom_vvl_sf_nxt: You should not have seen this print! error?', kt
+   END SUBROUTINE dom_vvl_sf_nxt
+   SUBROUTINE dom_vvl_sf_swp( kt ) 
+      !! * Arguments
+      INTEGER, INTENT( in )  ::    kt
+      WRITE(*,*) 'dom_vvl_sf_swp: You should not have seen this print! error?', kt
+   END SUBROUTINE dom_vvl_sf_swp
+   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out
+      CHARACTER(len=1), INTENT( in )                    ::  pout
+      WRITE(*,*) 'dom_vvl_interpol: You should not have seen this print! error?'
+   END SUBROUTINE dom_vvl_interpol
+
+
 #endif
 
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr_substitute.h90
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr_substitute.h90	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr_substitute.h90	(revision 2905)
@@ -19,39 +19,42 @@
 #   define  fse3uw_0(i,j,k)  e3uw(i,j,k)
 #   define  fse3vw_0(i,j,k)  e3vw(i,j,k)
+
 #if defined key_vvl
-! s* or z*-coordinate (3D + time dependency) + use of additional now arrays (..._1)
-#   define  fsdept(i,j,k)  gdept_1(i,j,k)
-#   define  fsdepw(i,j,k)  gdepw_1(i,j,k)
-#   define  fsde3w(i,j,k)  gdep3w_1(i,j,k)
-#   define  fse3t(i,j,k)   e3t_1(i,j,k)
-#   define  fse3u(i,j,k)   e3u_1(i,j,k)
-#   define  fse3v(i,j,k)   e3v_1(i,j,k)
-#   define  fse3f(i,j,k)   e3f_1(i,j,k)
-#   define  fse3w(i,j,k)   e3w_1(i,j,k)
-#   define  fse3uw(i,j,k)  e3uw_1(i,j,k)
-#   define  fse3vw(i,j,k)  e3vw_1(i,j,k)
+! s* or z*-coordinate (3D + time dependency) + use of additional now arrays (..._n)
 
 #   define  fse3t_b(i,j,k)   e3t_b(i,j,k)
 #   define  fse3u_b(i,j,k)   e3u_b(i,j,k)
 #   define  fse3v_b(i,j,k)   e3v_b(i,j,k)
-#   define  fse3uw_b(i,j,k)  (fse3uw_0(i,j,k)*(1.+sshu_b(i,j)*muu(i,j,k)))
-#   define  fse3vw_b(i,j,k)  (fse3vw_0(i,j,k)*(1.+sshv_b(i,j)*muv(i,j,k)))
+#   define  fse3uw_b(i,j,k)  e3uw_b(i,j,k)
+#   define  fse3vw_b(i,j,k)  e3vw_b(i,j,k)
 
-#   define  fsdept_n(i,j,k)  (fsdept_0(i,j,k)*(1.+sshn(i,j)*mut(i,j,k)))
-#   define  fsdepw_n(i,j,k)  (fsdepw_0(i,j,k)*(1.+sshn(i,j)*mut(i,j,k)))
-#   define  fsde3w_n(i,j,k)  (fsde3w_0(i,j,k)*(1.+sshn(i,j)*mut(i,j,k))-sshn(i,j))
-#   define  fse3t_n(i,j,k)   (fse3t_0(i,j,k)*(1.+sshn(i,j)*mut(i,j,k)))
-#   define  fse3u_n(i,j,k)   (fse3u_0(i,j,k)*(1.+sshu_n(i,j)*muu(i,j,k)))
-#   define  fse3v_n(i,j,k)   (fse3v_0(i,j,k)*(1.+sshv_n(i,j)*muv(i,j,k)))
-#   define  fse3f_n(i,j,k)   (fse3f_0(i,j,k)*(1.+sshf_n(i,j)*muf(i,j,k)))
-#   define  fse3w_n(i,j,k)   (fse3w_0(i,j,k)*(1.+sshn(i,j)*mut(i,j,k)))
-#   define  fse3uw_n(i,j,k)  (fse3uw_0(i,j,k)*(1.+sshu_n(i,j)*muu(i,j,k)))
-#   define  fse3vw_n(i,j,k)  (fse3vw_0(i,j,k)*(1.+sshv_n(i,j)*muv(i,j,k)))
+#   define  fsdept_n(i,j,k)  gdept_n(i,j,k)
+#   define  fsdepw_n(i,j,k)  gdepw_n(i,j,k)
+#   define  fsde3w_n(i,j,k)  gdep3w_n(i,j,k)
+#   define  fse3t_n(i,j,k)   e3t_n(i,j,k)
+#   define  fse3u_n(i,j,k)   e3u_n(i,j,k)
+#   define  fse3v_n(i,j,k)   e3v_n(i,j,k)
+#   define  fse3f_n(i,j,k)   e3f_n(i,j,k)
+#   define  fse3w_n(i,j,k)   e3w_n(i,j,k)
+#   define  fse3uw_n(i,j,k)  e3uw_n(i,j,k)
+#   define  fse3vw_n(i,j,k)  e3vw_n(i,j,k)
 
-#   define  fse3t_m(i,j,k)   (fse3t_0(i,j,k)*(1.+ssh_m(i,j)*mut(i,j,k)))
+#   define  fse3t_a(i,j,k)   e3t_a(i,j,k)
+#   define  fse3u_a(i,j,k)   e3u_a(i,j,k)
+#   define  fse3v_a(i,j,k)   e3v_a(i,j,k)
 
-#   define  fse3t_a(i,j,k)   (fse3t_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k)))
-#   define  fse3u_a(i,j,k)   (fse3u_0(i,j,k)*(1.+sshu_a(i,j)*muu(i,j,k)))
-#   define  fse3v_a(i,j,k)   (fse3v_0(i,j,k)*(1.+sshv_a(i,j)*muv(i,j,k)))
+#   define  fse3t_m(i,j)   e3t_m(i,j)
+
+! This part should be removed one day ...
+#   define  fsdept(i,j,k)  fsdept_n(i,j,k)
+#   define  fsdepw(i,j,k)  fsdepw_n(i,j,k)
+#   define  fsde3w(i,j,k)  fsde3w_n(i,j,k)
+#   define  fse3t(i,j,k)   fse3t_n(i,j,k) 
+#   define  fse3u(i,j,k)   fse3u_n(i,j,k) 
+#   define  fse3v(i,j,k)   fse3v_n(i,j,k) 
+#   define  fse3f(i,j,k)   fse3f_n(i,j,k) 
+#   define  fse3w(i,j,k)   fse3w_n(i,j,k) 
+#   define  fse3uw(i,j,k)  fse3uw_n(i,j,k)
+#   define  fse3vw(i,j,k)  fse3vw_n(i,j,k)
 
 #else
@@ -85,5 +88,5 @@
 #   define  fse3vw_n(i,j,k)  fse3vw_0(i,j,k)
 
-#   define  fse3t_m(i,j,k)   fse3t_0(i,j,k)
+#   define  fse3t_m(i,j)     fse3t_0(i,j,1)
 
 #   define  fse3t_a(i,j,k)   fse3t_0(i,j,k)
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 2905)
@@ -92,7 +92,5 @@
       !!               un,vn   now horizontal velocity of next time-step
       !!----------------------------------------------------------------------
-      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
       USE oce     , ONLY:   ze3u_f => ta       , ze3v_f => sa       ! (ta,sa) used as 3D workspace
-      USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3
       !
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
@@ -104,10 +102,6 @@
       REAL(wp) ::   zue3a, zue3n, zue3b, zuf    ! local scalars
       REAL(wp) ::   zve3a, zve3n, zve3b, zvf    !   -      -
-      REAL(wp) ::   zec, zv_t_ij, zv_t_ip1j, zv_t_ijp1
+      REAL(wp) ::   zec                         !   -      -
       !!----------------------------------------------------------------------
-
-      IF( wrk_in_use(2, 1,2,3) ) THEN
-         CALL ctl_stop('dyn_nxt: requested workspace arrays unavailable')   ;   RETURN
-      ENDIF
 
       IF( kt == nit000 ) THEN
@@ -185,6 +179,6 @@
          hvr_e(:,:) = hvr(:,:)
          DO jk = 1, jpkm1   !! Vertically integrated momentum trends
-            ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
-            va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
+            ua_e(:,:) = ua_e(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
+            va_e(:,:) = va_e(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
          END DO
          ua_e(:,:) = ua_e(:,:) * hur(:,:)
@@ -239,5 +233,6 @@
             !                             ! ================!
             ! Before scale factor at t-points
-            ! -------------------------------
+            ! (used as a now filtered scale factor until the swap)
+            ! ----------------------------------------------------
             DO jk = 1, jpkm1
                fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
@@ -245,33 +240,16 @@
                   &                         - 2.e0 * fse3t_n(:,:,jk)            )
             ENDDO
-            ! Add volume filter correction only at the first level of t-point scale factors
+            ! Add volume filter correction: comatibility with tracer advection scheme
+            ! => time filter + conservation correction (only at the first level)
             zec = atfp * rdt / rau0
             fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
-            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
-            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
-            zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:)
-            zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:)
+            ! swap of emp fields
+            emp_b(:,:) = emp(:,:)
             !
             IF( ln_dynadv_vec ) THEN
                ! Before scale factor at (u/v)-points
                ! -----------------------------------
-               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
-               DO jk = 1, jpkm1
-                  DO jj = 1, jpjm1
-                     DO ji = 1, jpim1
-                        zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
-                        zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
-                        zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
-                        fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
-                        fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
-                     END DO
-                  END DO
-               END DO
-               ! lateral boundary conditions
-               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
-               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
-               ! Add initial scale factor to scale factor anomaly
-               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
-               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
+               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
+               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
                ! Leap-Frog - Asselin filter and swap: applied on velocity
                ! -----------------------------------
@@ -291,24 +269,8 @@
                !
             ELSE
-               ! Temporary filered scale factor at (u/v)-points (will become before scale factor)
-               !-----------------------------------------------
-               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
-               DO jk = 1, jpkm1
-                  DO jj = 1, jpjm1
-                     DO ji = 1, jpim1
-                        zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
-                        zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
-                        zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
-                        ze3u_f(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
-                        ze3v_f(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
-                     END DO
-                  END DO
-               END DO
-               ! lateral boundary conditions
-               CALL lbc_lnk( ze3u_f, 'U', 1. )
-               CALL lbc_lnk( ze3v_f, 'V', 1. )
-               ! Add initial scale factor to scale factor anomaly
-               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
-               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
+               ! Temporary filtered scale factor at (u/v)-points (will become before scale factor)
+               !------------------------------------------------
+               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3u_f, 'U' )
+               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3v_f, 'V' )
                ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
                ! -----------------------------------             ===========================
@@ -345,6 +307,5 @@
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
          &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
-      ! 
-      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays')
+      !
       !
    END SUBROUTINE dyn_nxt
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 2905)
@@ -114,4 +114,5 @@
       USE wrk_nemo, ONLY:   zun => wrk_2d_14 , zun_e => wrk_2d_15 , zub_e => wrk_2d_16 , zu_sum => wrk_2d_17
       USE wrk_nemo, ONLY:   zvn => wrk_2d_18 , zvn_e => wrk_2d_19 , zvb_e => wrk_2d_20 , zv_sum => wrk_2d_21
+      USE wrk_nemo, ONLY:   zhu_b => wrk_2d_22 , zhv_b => wrk_2d_23
       !
       INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
@@ -127,5 +128,5 @@
 
       IF( wrk_in_use(2,  1, 2, 3, 4, 5, 6, 7, 8, 9,10,     &
-         &              11,12,13,14,15,16,17,18,19,20,21 ) ) THEN
+         &              11,12,13,14,15,16,17,18,19,20,21,22,23 ) ) THEN
          CALL ctl_stop( 'dyn_spg_ts: requested workspace arrays unavailable' )   ;   RETURN
       ENDIF
@@ -188,20 +189,30 @@
 #endif
                !                                                                              ! now trend
-               zua(ji,jj) = zua(ji,jj) + fse3u  (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
-               zva(ji,jj) = zva(ji,jj) + fse3v  (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
+               zua(ji,jj) = zua(ji,jj) + fse3u_n(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
+               zva(ji,jj) = zva(ji,jj) + fse3v_n(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
                !                                                                              ! now velocity 
-               zun(ji,jj) = zun(ji,jj) + fse3u  (ji,jj,jk) * un(ji,jj,jk)
-               zvn(ji,jj) = zvn(ji,jj) + fse3v  (ji,jj,jk) * vn(ji,jj,jk)               
+               zun(ji,jj) = zun(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk)
+               zvn(ji,jj) = zvn(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk)               
                !
-#if defined key_vvl
-               ub_b(ji,jj) = ub_b(ji,jj) + (fse3u_0(ji,jj,jk)*(1.+sshu_b(ji,jj)*muu(ji,jj,jk)))* ub(ji,jj,jk) 
-               vb_b(ji,jj) = vb_b(ji,jj) + (fse3v_0(ji,jj,jk)*(1.+sshv_b(ji,jj)*muv(ji,jj,jk)))* vb(ji,jj,jk)   
-#else
-               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_0(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk)
-               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_0(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk)
-#endif
+               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk)
+               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk)
             END DO
          END DO
       END DO
+
+      ! before inverse water column height at u- and v- points
+      IF( lk_vvl ) THEN
+         zhu_b(:,:) = 0.
+         zhv_b(:,:) = 0.
+         DO jk = 1, jpk
+            zhu_b(:,:) = zhu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
+            zhv_b(:,:) = zhv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
+         END DO
+         zhu_b(:,:) = umask(:,:,1) / ( zhu_b(:,:) + 1. - umask(:,:,1) )
+         zhv_b(:,:) = vmask(:,:,1) / ( zhv_b(:,:) + 1. - vmask(:,:,1) )
+      ELSE
+         zhu_b(:,:) = hur(:,:)
+         zhv_b(:,:) = hvr(:,:)
+      ENDIF
 
       !                                   !* baroclinic momentum trend (remove the vertical mean trend)
@@ -289,26 +300,15 @@
             !RBbug for vvl and external mode we may need to use varying fse3
             !!gm  Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx.
-            zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  )
-            zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  )
+            zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u_n(ji,jj,mbku(ji,jj)) * zcoef  )
+            zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v_n(ji,jj,mbkv(ji,jj)) * zcoef  )
          END DO
       END DO
-
-      IF( lk_vvl ) THEN
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   &
-                  &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) )
-               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj)   &
-                  &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) )
-            END DO
-         END DO
-      ELSE
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj)
-               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj)
-            END DO
-         END DO
-      ENDIF
+      
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * zhu_b(ji,jj)
+            zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * zhv_b(ji,jj)
+         END DO
+      END DO
 
       !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity)
@@ -317,11 +317,6 @@
       zva(:,:) = zva(:,:) * hvr(:,:)
       !
-      IF( lk_vvl ) THEN
-         ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) )
-         vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) )
-      ELSE
-         ub_b(:,:) = ub_b(:,:) * hur(:,:)
-         vb_b(:,:) = vb_b(:,:) * hvr(:,:)
-      ENDIF
+      ub_b(:,:) = ub_b(:,:) * zhu_b(:,:)
+      vb_b(:,:) = vb_b(:,:) * zhv_b(:,:)
 
       ! -----------------------------------------------------------------------
@@ -640,15 +635,6 @@
             ENDIF
 
-            IF( lk_vvl ) THEN
-               ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) )
-               vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) )
-            ELSE
-               ub_b(:,:) = ub_b(:,:) * hur(:,:)
-               vb_b(:,:) = vb_b(:,:) * hvr(:,:)
-            ENDIF
-         ELSE                                 ! neuler==0
-            ub_b (:,:) = un_b (:,:)
-            vb_b (:,:) = vn_b (:,:)
-         ENDIF
+            ub_b(:,:) = ub_b(:,:) * zhu_b(:,:)
+            vb_b(:,:) = vb_b(:,:) * zhv_b(:,:)
 
          IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 2905)
@@ -8,4 +8,5 @@
    !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
    !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
+   !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work
    !!----------------------------------------------------------------------
 
@@ -43,6 +44,7 @@
    PRIVATE
 
-   PUBLIC   ssh_wzv    ! called by step.F90
    PUBLIC   ssh_nxt    ! called by step.F90
+   PUBLIC   wzv        ! called by step.F90
+   PUBLIC   ssh_swp    ! called by step.F90
 
    !! * Substitutions
@@ -56,35 +58,29 @@
 CONTAINS
 
-   SUBROUTINE ssh_wzv( kt ) 
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE ssh_wzv  ***
+   SUBROUTINE ssh_nxt( kt )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE ssh_nxt  ***
       !!                   
-      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity
-      !!              and update the now vertical coordinate (lk_vvl=T).
-      !!
-      !! ** Method  : - Using the incompressibility hypothesis, the vertical 
-      !!      velocity is computed by integrating the horizontal divergence  
-      !!      from the bottom to the surface minus the scale factor evolution.
-      !!        The boundary conditions are w=0 at the bottom (no flux) and.
+      !! ** Purpose :   compute the after ssh (ssha)
+      !!
+      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
+      !!      is computed by integrating the horizontal divergence and multiply by
+      !!      by the time step.
       !!
       !! ** action  :   ssha    : after sea surface height
-      !!                wn      : now vertical velocity
-      !!                sshu_a, sshv_a, sshf_a  : after sea surface height (lk_vvl=T)
-      !!                hu, hv, hur, hvr        : ocean depth and its inverse at u-,v-points
       !!
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
       !!----------------------------------------------------------------------
       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE oce     , ONLY:   z3d   => ta                           ! ta used as 3D workspace
-      USE wrk_nemo, ONLY:   zhdiv => wrk_2d_1 , z2d => wrk_2d_2   ! 2D workspace
-      !
-      INTEGER, INTENT(in) ::   kt   ! time step
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars
-      !!----------------------------------------------------------------------
-
-      IF( wrk_in_use(2, 1,2) ) THEN
-         CALL ctl_stop('ssh_wzv: requested workspace arrays unavailable')   ;   RETURN
+      USE wrk_nemo, ONLY:   zhdiv => wrk_2d_1                     ! 2D workspace
+      !
+      INTEGER, INTENT(in) ::   kt                      ! time step
+      ! 
+      INTEGER             ::   jk                      ! dummy loop indice
+      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars
+      !!----------------------------------------------------------------------
+
+      IF( wrk_in_use(2, 1) ) THEN
+         CALL ctl_stop('ssh_nxt: requested workspace arrays unavailable')   ;   RETURN
       ENDIF
 
@@ -92,63 +88,7 @@
          !
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity '
+         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
          !
-         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
-         !
-         IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only)
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
-                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
-                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
-                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
-                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
-                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
-                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
-                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
-                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
-                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
-                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
-                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
-               END DO
-            END DO
-            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
-            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1      ! NO Vector Opt.
-                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
-                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
-                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
-                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
-               END DO
-            END DO
-            CALL lbc_lnk( sshf_n, 'F', 1. )
-         ENDIF
-         !
-      ENDIF
-
-      !                                           !------------------------------------------!
-      IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case)
-         !                                        !------------------------------------------!
-         DO jk = 1, jpkm1
-            fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays
-            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
-            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
-            !
-            fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays
-            fse3u (:,:,jk) = fse3u_n (:,:,jk)
-            fse3v (:,:,jk) = fse3v_n (:,:,jk)
-            fse3f (:,:,jk) = fse3f_n (:,:,jk)
-            fse3w (:,:,jk) = fse3w_n (:,:,jk)
-            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
-            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
-         END DO
-         !
-         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
-         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
-         !                                            ! now masked inverse of the ocean depth (at u- and v-points)
-         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) )
-         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1._wp - vmask(:,:,1) )
-         ! 
       ENDIF
       !
@@ -184,20 +124,4 @@
       CALL lbc_lnk( ssha, 'T', 1. ) 
 #endif
-
-      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
-      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1      ! NO Vector Opt.
-               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
-                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
-                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
-               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
-                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
-                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
-            END DO
-         END DO
-         CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions
-      ENDIF
-      
 #if defined key_asminc
       !                                                ! Include the IAU weighted SSH increment
@@ -209,23 +133,95 @@
 
       !                                           !------------------------------!
+      !                                           !           outputs            !
+      !                                           !------------------------------!
+      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
+      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
+      !
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
+      !
+      IF( wrk_not_released(2, 1) )   CALL ctl_stop('ssh_nxt: failed to release workspace arrays')
+      !
+   END SUBROUTINE ssh_nxt
+
+   
+   SUBROUTINE wzv( kt )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE wzv  ***
+      !!                   
+      !! ** Purpose :   compute the now vertical velocity
+      !!
+      !! ** Method  : - Using the incompressibility hypothesis, the vertical 
+      !!      velocity is computed by integrating the horizontal divergence  
+      !!      from the bottom to the surface minus the scale factor evolution.
+      !!        The boundary conditions are w=0 at the bottom (no flux) and.
+      !!
+      !! ** action  :   wn      : now vertical velocity
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
+      USE oce     , ONLY:   z3d   => ta                           ! ta used as 3D workspace
+      USE wrk_nemo, ONLY:   zhdiv => wrk_3d_1 , z2d => wrk_2d_1   ! 2D workspace
+      !
+      INTEGER, INTENT(in) ::   kt           ! time step
+      !
+      INTEGER             ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp)            ::   z1_2dt       ! local scalars
+      !!----------------------------------------------------------------------
+      
+      IF( wrk_in_use(2, 1) .OR. wrk_in_use(3, 1) ) THEN
+         CALL ctl_stop('wzv: requested workspace arrays unavailable')   ;   RETURN
+      ENDIF
+
+      IF( kt == nit000 ) THEN
+         !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
+         IF(lwp) WRITE(numout,*) '~~~ '
+         !
+         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
+         !
+      ENDIF
+      !                                           !------------------------------!
       !                                           !     Now Vertical Velocity    !
       !                                           !------------------------------!
-      z1_2dt = 1.e0 / z2dt
-      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
-         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
-         wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
-            &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
-            &                         * tmask(:,:,jk) * z1_2dt
+      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
+      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
+      !
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
+         DO jk = 1, jpkm1
+            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
+            ! - ML - note: computation allready done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zhdiv(ji,jj,jk) = e12t_1(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk(zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
+         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
+         !                             ! Same question holds for hdivn. Perhaps just for security
+         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
+            ! computation of w
+            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)                    &
+               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
+         END DO
+         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
+      ELSE   ! z_star and linear free surface cases
+         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
+            ! computation of w
+            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk)                                   &
+               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
+         END DO
+      ENDIF
+
 #if defined key_bdy
          wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
 #endif
-      END DO
 
       !                                           !------------------------------!
       !                                           !           outputs            !
       !                                           !------------------------------!
-      CALL iom_put( "woce", wn                    )   ! vertical velocity
-      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
-      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
+      CALL iom_put( "woce", wn )   ! vertical velocity
       IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
          ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
@@ -237,13 +233,11 @@
          CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
       ENDIF
-      !
-      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
-      !
-      IF( wrk_not_released(2, 1,2) )   CALL ctl_stop('ssh_wzv: failed to release workspace arrays')
-      !
-   END SUBROUTINE ssh_wzv
-
-
-   SUBROUTINE ssh_nxt( kt )
+
+      IF( wrk_not_released(2, 1) .OR. wrk_not_released(3, 1) )   CALL ctl_stop('wzv: failed to release workspace arrays')
+
+   END SUBROUTINE wzv
+
+
+   SUBROUTINE ssh_swp( kt )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE ssh_nxt  ***
@@ -251,5 +245,5 @@
       !! ** Purpose :   achieve the sea surface  height time stepping by 
       !!              applying Asselin time filter and swapping the arrays
-      !!              ssha  already computed in ssh_wzv  
+      !!              ssha  already computed in ssh_nxt  
       !!
       !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
@@ -266,83 +260,24 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !!
-      INTEGER  ::   ji, jj   ! dummy loop indices
-      REAL(wp) ::   zec      ! temporary scalar
+      REAL(wp)            ::   zec  ! temporary scalar
       !!----------------------------------------------------------------------
 
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
+         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
       ENDIF
 
-      !                       !--------------------------!
-      IF( lk_vvl ) THEN       !  Variable volume levels  !     (ssh at t-, u-, v, f-points)
-         !                    !--------------------------!
-         !
-         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
-            sshn  (:,:) = ssha  (:,:)                       ! now <-- after  (before already = now)
-            sshu_n(:,:) = sshu_a(:,:)
-            sshv_n(:,:) = sshv_a(:,:)
-            DO jj = 1, jpjm1                                ! ssh now at f-point
-               DO ji = 1, jpim1      ! NO Vector Opt.
-                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
-                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
-                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
-                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
-               END DO
-            END DO
-            CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions
-            !
-         ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
+      IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
+         sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now)
+      ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
+         IF( lk_vvl ) THEN                               ! before <-- now filtered
             zec = atfp * rdt / rau0
-            DO jj = 1, jpj
-               DO ji = 1, jpi                               ! before <-- now filtered
-                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   &
-                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
-                  sshn  (ji,jj) = ssha  (ji,jj)             ! now <-- after
-                  sshu_n(ji,jj) = sshu_a(ji,jj)
-                  sshv_n(ji,jj) = sshv_a(ji,jj)
-               END DO
-            END DO
-            DO jj = 1, jpjm1                                ! ssh now at f-point
-               DO ji = 1, jpim1      ! NO Vector Opt.
-                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
-                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
-                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
-                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
-               END DO
-            END DO
-            CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions
-            !
-            DO jj = 1, jpjm1                                ! ssh before at u- & v-points
-               DO ji = 1, jpim1      ! NO Vector Opt.
-                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
-                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
-                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
-                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
-                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
-                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
-               END DO
-            END DO
-            CALL lbc_lnk( sshu_b, 'U', 1. )
-            CALL lbc_lnk( sshv_b, 'V', 1. )            !  Boundaries conditions
-            !
+            sshb  (:,:) = sshn  (:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )   &
+               &                          - zec  * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
+         ELSE
+            sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )
          ENDIF
-         !                    !--------------------------!
-      ELSE                    !        fixed levels      !     (ssh at t-point only)
-         !                    !--------------------------!
-         !
-         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
-            sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now)
-            !
-         ELSE                                               ! Leap-Frog time-stepping: Asselin filter + swap
-            DO jj = 1, jpj
-               DO ji = 1, jpi                               ! before <-- now filtered
-                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
-                  sshn(ji,jj) = ssha(ji,jj)                 ! now <-- after
-               END DO
-            END DO
-         ENDIF
-         !
+         sshn(:,:) = ssha(:,:)                           ! now <-- after
       ENDIF
       !
@@ -354,5 +289,5 @@
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
       !
-   END SUBROUTINE ssh_nxt
+   END SUBROUTINE ssh_swp
 
    !!======================================================================
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 2905)
@@ -23,5 +23,4 @@
    USE eosbn2          ! equation of state            (eos bn2 routine)
    USE trdmld_oce      ! ocean active mixed layer tracers trends variables
-   USE domvvl          ! variable volume
    USE traswp          ! swap from 4D T-S to 3D T & S and vice versa
 
@@ -122,5 +121,4 @@
                      CALL iom_rstput( kt, nitrst, numrow, 'hdivb'  , hdivb     )
                      CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb      )
-      IF( lk_vvl )   CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
                      !
                      CALL iom_rstput( kt, nitrst, numrow, 'un'     , un        )     ! now fields
@@ -191,5 +189,4 @@
                      CALL iom_get( numror, jpdom_autoglo, 'hdivb'  , hdivb   )
                      CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb    )
-      IF( lk_vvl )   CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
                      !
                      CALL iom_get( numror, jpdom_autoglo, 'un'     , un      )   ! now    fields
@@ -218,9 +215,4 @@
          hdivb(:,:,:) = hdivn(:,:,:)
          sshb (:,:)   = sshn (:,:)
-         IF( lk_vvl ) THEN
-            DO jk = 1, jpk
-               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
-            END DO
-         ENDIF
       ENDIF
       !
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 2905)
@@ -82,4 +82,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sss_m     !: mean (nn_fsbc time-step) surface sea salinity            [psu]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssh_m     !: mean (nn_fsbc time-step) sea surface height                [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3t_m     !: mean (nn_fsbc time-step) sea surface height                [m]
 
    !! * Substitutions
@@ -96,5 +97,5 @@
       !!                  ***  FUNCTION sbc_oce_alloc  ***
       !!---------------------------------------------------------------------
-      INTEGER :: ierr(4)
+      INTEGER :: ierr(5)
       !!---------------------------------------------------------------------
       ierr(:) = 0
@@ -117,4 +118,8 @@
          &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) ,                       &
          &      ssv_m  (jpi,jpj) , sss_m  (jpi,jpj), ssh_m(jpi,jpj) , STAT=ierr(4) )
+         !
+#if defined key_vvl
+      ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) )
+#endif
          !
       sbc_oce_alloc = MAXVAL( ierr )
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 2905)
@@ -70,5 +70,6 @@
          ELSE                    ;   ssh_m(:,:) = sshn(:,:)
          ENDIF
-
+         !
+         IF( lk_vvl )   fse3t_m(:,:) = fse3t_n(:,:,1)
          !
       ELSE
@@ -86,4 +87,8 @@
                CALL iom_get( numror, jpdom_autoglo, 'sss_m'  , sss_m  )   !   "         "    salinity    (T-point)
                CALL iom_get( numror, jpdom_autoglo, 'ssh_m'  , ssh_m  )   !   "         "    height      (T-point)
+               IF( lk_vvl ) THEN
+                  CALL iom_get( numror, jpdom_autoglo, 'fse3t_m', fse3t_m(:,:) )
+                  !                                                       !   "         "    vertical scale factor (T-point)
+               ENDIF
                !
                IF( zf_sbc /= REAL( nn_fsbc, wp ) ) THEN      ! nn_fsbc has changed between 2 runs
@@ -96,4 +101,5 @@
                   sss_m(:,:) = zcoef * sss_m(:,:)
                   ssh_m(:,:) = zcoef * ssh_m(:,:)
+                  IF( lk_vvl )   fse3t_m(:,:) = zcoef * fse3t_m(:,:)
                ELSE
                   IF(lwp) WRITE(numout,*) '~~~~~~~   mean fields read in the ocean restart file'
@@ -110,5 +116,6 @@
                ELSE                    ;   ssh_m(:,:) = zcoef *   sshn(:,:)
                ENDIF
-
+               IF( lk_vvl )   fse3t_m(:,:) = zcoef * fse3t_n(:,:,1)
+               !
             ENDIF
             !                                             ! ---------------------------------------- !
@@ -120,4 +127,5 @@
             sss_m(:,:) = 0.e0
             ssh_m(:,:) = 0.e0
+            IF( lk_vvl )   fse3t_m(:,:) = 0.e0
          ENDIF
          !                                                ! ---------------------------------------- !
@@ -132,4 +140,5 @@
          ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:)
          ENDIF
+         IF( lk_vvl )   fse3t_m(:,:) = fse3t_m(:,:) + fse3t_n(:,:,1)
 
          !                                                ! ---------------------------------------- !
@@ -142,4 +151,5 @@
             ssv_m(:,:) = ssv_m(:,:) * zcoef           !
             ssh_m(:,:) = ssh_m(:,:) * zcoef           ! mean SSH             [m]
+            IF( lk_vvl )   fse3t_m(:,:) = fse3t_m(:,:) * zcoef   ! mean vertical scale factor [m]
             !
          ENDIF
@@ -158,4 +168,7 @@
             CALL iom_rstput( kt, nitrst, numrow, 'sss_m'  , sss_m  )
             CALL iom_rstput( kt, nitrst, numrow, 'ssh_m'  , ssh_m  )
+            IF( lk_vvl ) THEN
+               CALL iom_rstput( kt, nitrst, numrow, 'fse3t_m'  , fse3t_m(:,:)  )
+            END IF
             !
          ENDIF
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 2905)
@@ -14,4 +14,5 @@
    USE oce             ! ocean dynamics and active tracers
    USE dom_oce         ! ocean space and time domain
+   USE domvvl          ! variable vertical scale factors
    USE traadv_cen2     ! 2nd order centered scheme (tra_adv_cen2   routine)
    USE traadv_tvd      ! TVD      scheme           (tra_adv_tvd    routine)
@@ -88,4 +89,10 @@
          zwn(:,:,jk) = e1t(:,:) * e2t(:,:)      * wn(:,:,jk)
       END DO
+      !
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
+         zun(:,:,:) = zun(:,:,:) + un_td(:,:,:)
+         zvn(:,:,:) = zvn(:,:,:) + vn_td(:,:,:)
+      ENDIF
+      !
       zun(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
       zvn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90	(revision 2905)
@@ -107,6 +107,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
-                  zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
+                  zeeu(ji,jj) = e1ur(ji,jj) * fse3u_n(ji,jj,jk) * umask(ji,jj,jk)
+                  zeev(ji,jj) = e2vr(ji,jj) * fse3v_n(ji,jj,jk) * vmask(ji,jj,jk)
                END DO
             END DO
@@ -130,5 +130,5 @@
             DO jj = 2, jpjm1                 ! Second derivative (divergence) time the eddy diffusivity coefficient
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t_n(ji,jj,jk) )
                   zlt(ji,jj) = fsahtt(ji,jj,jk) * zbtr * (   ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   &
                      &                                     + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)   )
@@ -148,5 +148,5 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! horizontal diffusive trends
-                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t_n(ji,jj,jk) )
                   ztra = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
                   ! add it to the general tracer trends
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 2905)
@@ -172,6 +172,6 @@
             DO jj = 1 , jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)
-                  zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)
+                  zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * e1ur(ji,jj) * fse3u_n(ji,jj,jk)
+                  zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * e2vr(ji,jj) * fse3v_n(ji,jj,jk)
                   !
                   zmsku = 1. / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   &
@@ -197,5 +197,5 @@
             DO jj = 2 , jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t_n(ji,jj,jk) )
                   ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )
                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
@@ -282,5 +282,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t_n(ji,jj,jk) )
                   ztra = (  ztfw(ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * zbtr
                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90	(revision 2905)
@@ -30,6 +30,4 @@
 
    PUBLIC   tra_ldf_lap   ! routine called by step.F90
-
-   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:) ::   e1ur, e2vr   ! scale factor coefficients
 
    !! * Substitutions
@@ -81,11 +79,4 @@
          IF(lwp) WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
-         !
-         ALLOCATE( e1ur(jpi,jpj), e2vr(jpi,jpj), STAT=ierr )
-         IF( lk_mpp    )   CALL mpp_sum( ierr )
-         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'tra_ldf_lap : unable to allocate arrays' )
-         !
-         e1ur(:,:) = e2u(:,:) / e1u(:,:)
-         e2vr(:,:) = e1v(:,:) / e2v(:,:)
       ENDIF
 
@@ -99,6 +90,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zabe1 = fsahtu(ji,jj,jk) * umask(ji,jj,jk) * e1ur(ji,jj) * fse3u(ji,jj,jk)
-                  zabe2 = fsahtv(ji,jj,jk) * vmask(ji,jj,jk) * e2vr(ji,jj) * fse3v(ji,jj,jk)
+                  zabe1 = fsahtu(ji,jj,jk) * umask(ji,jj,jk) * e1ur(ji,jj) * fse3u_n(ji,jj,jk)
+                  zabe2 = fsahtv(ji,jj,jk) * vmask(ji,jj,jk) * e2vr(ji,jj) * fse3v_n(ji,jj,jk)
                   ztu(ji,jj,jk) = zabe1 * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
                   ztv(ji,jj,jk) = zabe2 * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
@@ -112,9 +103,9 @@
                      ikv = mbkv(ji,jj)
                      IF( iku == jk ) THEN
-                        zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * e1ur(ji,jj) * fse3u(ji,jj,iku)
+                        zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * e1ur(ji,jj) * fse3u_n(ji,jj,iku)
                         ztu(ji,jj,jk) = zabe1 * pgu(ji,jj,jn)
                      ENDIF
                      IF( ikv == jk ) THEN
-                        zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * e2vr(ji,jj) * fse3v(ji,jj,ikv)
+                        zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * e2vr(ji,jj) * fse3v_n(ji,jj,ikv)
                         ztv(ji,jj,jk) = zabe2 * pgv(ji,jj,jn)
                      ENDIF
@@ -128,5 +119,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1._wp / ( e1t(ji,jj) *e2t(ji,jj) * fse3t(ji,jj,jk) )
+                  zbtr = 1._wp / ( e1t(ji,jj) *e2t(ji,jj) * fse3t_n(ji,jj,jk) )
                   ! horizontal diffusive trends added to the general tracer trends
                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   &
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/oce.F90
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 2905)
@@ -36,7 +36,4 @@
    !! ------------                                      !  fields  ! fields ! trends !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshb   , sshn   , ssha   !: sea surface height at t-point [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshu_b , sshu_n , sshu_a !: sea surface height at u-point [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshv_b , sshv_n , sshv_a !: sea surface height at u-point [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::            sshf_n          !: sea surface height at f-point [m]
    !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   spgu, spgv               !: horizontal surface pressure gradient
@@ -74,7 +71,4 @@
          &     rhop(jpi,jpj,jpk) ,                                         &
          &     sshb  (jpi,jpj)   , sshn  (jpi,jpj) , ssha  (jpi,jpj) ,     &
-         &     sshu_b(jpi,jpj)   , sshu_n(jpi,jpj) , sshu_a(jpi,jpj) ,     &
-         &     sshv_b(jpi,jpj)   , sshv_n(jpi,jpj) , sshv_a(jpi,jpj) ,     &
-         &                         sshf_n(jpi,jpj) ,                       &
          &     spgu  (jpi,jpj)   , spgv(jpi,jpj)   ,                       &
          &     gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts),                     &
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 2905)
@@ -102,8 +102,9 @@
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      !  Ocean dynamics : ssh, wn, hdiv, rot                                 !
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL ssh_wzv( kstp )         ! after ssh & vertical velocity
-
+      !  Ocean dynamics : hdiv, rot, ssh, e3, wn
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         CALL ssh_nxt       ( kstp )  ! after ssh
+      IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors
+                         CALL wzv           ( kstp )  ! now cross-level velocity
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Ocean physics update                (ua, va, ta, sa used as workspace)
@@ -232,8 +233,9 @@
                                CALL dyn_nxt( kstp )         ! lateral velocity at next time step
 
-                               CALL ssh_nxt( kstp )         ! sea surface height at next time step
+                               CALL ssh_swp( kstp )         ! swap of sea surface height
+      IF( lk_vvl           )   CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors
 
       IF( ln_diahsb        )   CALL dia_hsb( kstp )         ! - ML - global conservation diagnostics
-      IF( lk_diaobs  )         CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
+      IF( lk_diaobs        )   CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Index: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 2863)
+++ branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 2905)
@@ -57,4 +57,6 @@
 
    USE sshwzv           ! vertical velocity and ssh        (ssh_wzv routine)
+   USE domvvl           ! variable vertical scale factors  (dom_vvl_sf_nxt routine)
+   !                                                       (dom_vvl_sf_swp routine)
 
    USE ldfslp           ! iso-neutral slopes               (ldf_slp routine)
