Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 3865)
@@ -638,4 +638,13 @@
    ln_dynadv_cen2= .false. !  flux form - 2nd order centered scheme
    ln_dynadv_ubs = .false. !  flux form - 3rd order UBS      scheme
+/
+!-----------------------------------------------------------------------
+&nam_vvl    !   vertical coordinate options
+!-----------------------------------------------------------------------
+   ln_vvl_zstar  = .false. !  zstar vertical coordinate                    
+   ln_vvl_ztilde = .true.  !  hybrid verticalcoordinate: only high frequency variations
+   ln_vvl_layer  = .false. !  full layer vertical coordinate
+   rn_ahe3       = 0.e0    !  thickness diffusion coefficient 
+   ln_vvl_dbg    = .true.  !  debug prints    (T/F)
 /
 !-----------------------------------------------------------------------
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90	(revision 3865)
@@ -237,5 +237,5 @@
             
             !  energy needed to bring ocean surface layer until its freezing
-            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1) * ( tfu(ji,jj) - sst_m(ji,jj) - rt0 ) * ( 1 - zinda )
+            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj) * ( tfu(ji,jj) - sst_m(ji,jj) - rt0 ) * ( 1 - zinda )
             
             !  calculate oceanic heat flux.
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 3865)
@@ -196,6 +196,6 @@
       thick0(:,:) = 0._wp
       DO jk = 1, jpkm1
-         vol0        = vol0        + SUM( area (:,:) * tmask(:,:,jk) * fse3t_0(:,:,jk) )
-         thick0(:,:) = thick0(:,:) +    tmask_i(:,:) * tmask(:,:,jk) * fse3t_0(:,:,jk)
+         vol0        = vol0        + SUM( area (:,:) * tmask(:,:,jk) * e3t_0(:,:,jk) )
+         thick0(:,:) = thick0(:,:) +    tmask_i(:,:) * tmask(:,:,jk) * e3t_0(:,:,jk)
       END DO
       IF( lk_mpp )   CALL mpp_sum( vol0 )
@@ -212,5 +212,5 @@
                ik = mbkt(ji,jj)
                IF( ik > 1 ) THEN
-                  zztmp = ( gdept_1d(ik) - fsdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
+                  zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
                   sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1)
                ENDIF
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90	(revision 3865)
@@ -326,6 +326,4 @@
             htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem) * MIN( fse3t(ji,jj,ilevel+1), zthick(ji,jj) )  &
                                                                    * tmask(ji,jj,ilevel+1)
-            htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem) * MIN( fse3t(ji,jj,ilevel+1), zthick(ji,jj) )   &
-               &                                                   * tmask(ji,jj,ilevel+1)
          END DO
       END DO
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 3865)
@@ -259,4 +259,5 @@
       !
 #if defined key_mpp_mpi
+      ijpjjpk = jpj*jpk
       ish(1) = ijpjjpk  ;   ish2(1) = jpj   ;   ish2(2) = jpk
       zwork(1:ijpjjpk) = RESHAPE( p_fval, ish )
@@ -314,4 +315,5 @@
       END DO
 #if defined key_mpp_mpi
+      ijpjjpk = jpj*jpk
       ish(1) = jpj*jpk   ;   ish2(1) = jpj   ;   ish2(2) = jpk
       zwork(1:ijpjjpk)= RESHAPE( p_fval, ish )
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 3865)
@@ -142,17 +142,30 @@
       ENDIF
 
-      CALL iom_put( "toce"   , tsn(:,:,:,jp_tem)                     )    ! temperature
-      CALL iom_put( "soce"   , tsn(:,:,:,jp_sal)                     )    ! salinity
-      CALL iom_put( "sst"    , tsn(:,:,1,jp_tem)                     )    ! sea surface temperature
-      CALL iom_put( "sst2"   , tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) )    ! square of sea surface temperature
-      CALL iom_put( "sss"    , tsn(:,:,1,jp_sal)                     )    ! sea surface salinity
-      CALL iom_put( "sss2"   , tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) )    ! square of sea surface salinity
-      CALL iom_put( "uoce"   , un                                    )    ! i-current      
-      CALL iom_put( "suoce"  , un(:,:,1)                             )    ! surface i-current      
-      CALL iom_put( "voce"   , vn                                    )    ! j-current
-      CALL iom_put( "svoce"  , vn(:,:,1)                             )    ! surface j-current
- 
-      CALL iom_put( "avt"    , avt                                   )    ! T vert. eddy diff. coef.
-      CALL iom_put( "avm"    , avmu                                  )    ! T vert. eddy visc. coef.
+      IF( lk_vvl ) THEN
+         z3d(:,:,:) = tsn(:,:,:,jp_tem) * fse3t_n(:,:,:)
+         CALL iom_put( "toce" , z3d                        )   ! heat content
+         CALL iom_put( "sst"  , z3d(:,:,1)                 )   ! sea surface heat content
+         CALL iom_put( "sst2" , z3d(:,:,1) * z3d(:,:,1)    )   ! sea surface content of squared temperature
+         z3d(:,:,:) = tsn(:,:,:,jp_sal) * fse3t_n(:,:,:)            
+         CALL iom_put( "soce" , z3d                        )   ! salinity content
+         CALL iom_put( "sss"  , z3d(:,:,1)                 )   ! sea surface salinity content
+         CALL iom_put( "sss2" , z3d(:,:,1) * z3d(:,:,1)    )   ! sea surface content of squared salinity
+      ELSE
+         CALL iom_put( "toce" , tsn(:,:,:,jp_tem)          )   ! temperature
+         CALL iom_put( "sst"  , tsn(:,:,1,jp_tem)          )   ! sea surface temperature
+         CALL iom_put( "sst2" , tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) ) ! square of sea surface temperature
+         CALL iom_put( "soce" , tsn(:,:,:,jp_sal)          )   ! salinity
+         CALL iom_put( "sss"  , tsn(:,:,1,jp_sal)          )   ! sea surface salinity
+         CALL iom_put( "sss2" , tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) ) ! square of sea surface salinity
+      END IF
+      IF( lk_vvl .AND. (.NOT. ln_dynadv_vec) ) THEN
+         CALL iom_put( "uoce" , un(:,:,:) * fse3u_n(:,:,:) )    ! i-transport
+         CALL iom_put( "voce" , vn(:,:,:) * fse3v_n(:,:,:) )    ! j-transport
+      ELSE
+         CALL iom_put( "uoce" , un                         )    ! i-current
+         CALL iom_put( "voce" , vn                         )    ! j-current
+      END IF
+      CALL iom_put(    "avt"  , avt                        )    ! T vert. eddy diff. coef.
+      CALL iom_put(    "avm"  , avmu                       )    ! T vert. eddy visc. coef.
       IF( lk_zdfddm ) THEN
          CALL iom_put( "avs" , fsavs(:,:,:)                          )    ! S vert. eddy diff. coef.
@@ -250,5 +263,5 @@
       !
       CALL wrk_alloc( jpi , jpj      , zw2d )
-      IF ( ln_traldf_gdia )  call wrk_alloc( jpi , jpj , jpk  , zw3d )
+      IF ( ln_traldf_gdia .OR. lk_vvl )  call wrk_alloc( jpi , jpj , jpk  , zw3d )
       !
       ! Output the initial state and forcings
@@ -395,4 +408,12 @@
          CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
             &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
+         IF(  lk_vvl  ) THEN
+            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n
+            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
+            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t_n
+            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
+            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t_n
+            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
+         ENDIF
          !                                                                                      !!! nid_T : 2D
          CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
@@ -406,12 +427,12 @@
          CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-#if ! defined key_vvl
-         CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"        &  ! emp * tsn(:,:,1,jp_tem)
+         IF(  .NOT. lk_vvl  ) THEN
+            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem)
             &                                                                  , "KgC/m2/s",  &  ! sosst_cd
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"           &  ! emp * tsn(:,:,1,jp_sal)
+            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * tsn(:,:,1,jp_sal)
             &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-#endif
+            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+         ENDIF
          CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
@@ -585,5 +606,5 @@
       ! ---------------------
 
-      ! ndex(1) est utilise ssi l'avant dernier argument est diffferent de 
+      ! ndex(1) est utilise ssi l'avant dernier argument est different de 
       ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
       ! donne le nombre d'elements, et ndex la liste des indices a sortir
@@ -595,8 +616,22 @@
 
       ! Write fields on T grid
-      CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem), ndim_T , ndex_T  )   ! temperature
-      CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal), ndim_T , ndex_T  )   ! salinity
-      CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem), ndim_hT, ndex_hT )   ! sea surface temperature
-      CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal), ndim_hT, ndex_hT )   ! sea surface salinity
+      IF( lk_vvl ) THEN
+         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
+         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
+         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
+         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
+      ELSE
+         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature
+         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity
+         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature
+         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
+
+      ENDIF
+      IF( lk_vvl ) THEN
+         z3d(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
+         CALL histwrite( nid_T, "vovvle3t", it, fse3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
+         CALL histwrite( nid_T, "vovvldep", it, fsdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
+         CALL histwrite( nid_T, "vovvldef", it, z3d             , ndim_T , ndex_T  )   ! level thickness deformation
+      ENDIF
       CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
       CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
@@ -604,10 +639,10 @@
                                                                                   ! (includes virtual salt flux beneath ice 
                                                                                   ! in linear free surface case)
-#if ! defined key_vvl
-      zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
-      CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )             ! c/d term on sst
-      zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
-      CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )             ! c/d term on sss
-#endif
+      IF( .NOT. lk_vvl ) THEN
+         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
+         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
+         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
+         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
+      ENDIF
       CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
       CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
@@ -750,5 +785,5 @@
       !
       CALL wrk_dealloc( jpi , jpj      , zw2d )
-      IF ( ln_traldf_gdia )  call wrk_dealloc( jpi , jpj , jpk  , zw3d )
+      IF ( ln_traldf_gdia .OR. lk_vvl )  call wrk_dealloc( jpi , jpj , jpk  , zw3d )
       !
       IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
@@ -839,4 +874,8 @@
       CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2"   ,   &   ! j-wind stress
          &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+      IF( lk_vvl ) THEN
+         CALL histdef( id_i, "vovvldep", "T point depth"         , "m"      ,   &   ! t-point depth
+            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
+      END IF
 
 #if defined key_lim2
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 3865)
@@ -134,10 +134,10 @@
    !! All coordinates
    !! ---------------
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdep3w          !: depth of T-points (sum of e3w) (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept , gdepw   !: analytical depth at T-W  points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3v   , e3f     !: analytical vertical scale factors at  V--F
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t   , e3u     !:                                       T--U  points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3vw            !: analytical vertical scale factors at  VW--
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3w   , e3uw    !:                                        W--UW  points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdep3w_0           !: depth of t-points (sum of e3w) (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_0, gdepw_0   !: analytical (time invariant) depth at t-w  points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3v_0  , e3f_0     !: analytical (time invariant) vertical scale factors at  v-f
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_0  , e3u_0     !:                                      t-u  points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3vw_0             !: analytical (time invariant) vertical scale factors at  vw
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3w_0  , e3uw_0    !:                                      w-uw points (m)
 #if defined key_vvl
    LOGICAL, PUBLIC, PARAMETER ::   lk_vvl = .TRUE.    !: variable grid flag
@@ -145,18 +145,28 @@
    !! 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(:,:), TARGET ::   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(:,:), TARGET ::   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           !: reference depth at t-       points (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   hu_0 , hv_0    !: reference depth at u- and v-points (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   re2u_e1u       !: scale factor coeffs at u points (e2u/e1u)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   re1v_e2v       !: scale factor coeffs at v points (e1v/e2v)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   e12t , r1_e12t !: horizontal cell surface and inverse at t points
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   e12u , r1_e12u !: horizontal cell surface and inverse at u points
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   e12v , r1_e12v !: horizontal cell surface and inverse at v points
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   e12f , r1_e12f !: horizontal cell surface and inverse at f points
 
    INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1)
@@ -165,39 +175,39 @@
    !! z-coordinate with full steps (also used in the other cases as reference z-coordinate)
    !! =-----------------====------
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   :: gdept_1d, gdepw_1d !: reference depth of t- and w-points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   :: e3t_1d  , e3w_1d   !: reference vertical scale factors at T- and W-pts (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e3tp    , e3wp     !: ocean bottom level thickness at T and W points
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   gdept_1d, gdepw_1d !: reference depth of t- and w-points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   e3t_1d  , e3w_1d   !: reference vertical scale factors at T- and W-pts (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3tp    , e3wp     !: ocean bottom level thickness at T and W points
 
    !! s-coordinate and hybrid z-s-coordinate
    !! =----------------======---------------
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   gsigt, gsigw   !: model level depth coefficient at t-, w-levels (analytic)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   gsi3w          !: model level depth coefficient at w-level (sum of gsigw)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   esigt, esigw   !: vertical scale factor coef. at t-, w-levels
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hbatv , hbatf    !: ocean depth at the vertical of  V--F
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hbatt , hbatu    !:                                 T--U  points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   scosrf, scobot   !: ocean surface and bottom topographies 
-   !                                        !  (if deviating from coordinate surfaces in HYBRID)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hifv  , hiff     !: interface depth between stretching at  V--F
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hift  , hifu     !: and quasi-uniform spacing              T--U  points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rx1              !: Maximum grid stiffness ratio
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   gsigt, gsigw       !: model level depth coefficient at t-, w-levels (analytic)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   gsi3w              !: model level depth coefficient at w-level (sum of gsigw)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   esigt, esigw       !: vertical scale factor coef. at t-, w-levels
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hbatv , hbatf      !: ocean depth at the vertical of  v--f
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hbatt , hbatu      !:                                 t--u points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   scosrf, scobot     !: ocean surface and bottom topographies 
+   !                                                                           !  (if deviating from coordinate surfaces in HYBRID)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hifv  , hiff       !: interface depth between stretching at v--f
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hift  , hifu       !: and quasi-uniform spacing             t--u points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rx1                !: Maximum grid stiffness ratio
 
    !!----------------------------------------------------------------------
    !! masks, bathymetry
    !! ---------------------------------------------------------------------
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbathy       !: number of ocean level (=0, 1, ... , jpk-1)
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbkt         !: vertical index of the bottom last T- ocean level
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbku, mbkv   !: vertical index of the bottom last U- and W- ocean level
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy        !: ocean depth (meters)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i      !: interior domain T-point mask
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask        !: land/ocean mask of barotropic stream function
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbathy             !: number of ocean level (=0, 1, ... , jpk-1)
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbkt               !: vertical index of the bottom last T- ocean level
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbku, mbkv         !: vertical index of the bottom last U- and W- ocean level
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy              !: ocean depth (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i            !: interior domain T-point mask
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask              !: land/ocean mask of barotropic stream function
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts
 
-   REAL(wp), PUBLIC, DIMENSION(jpiglo) ::   tpol, fpol          !: north fold mask (jperio= 3 or 4)
+   REAL(wp), PUBLIC, DIMENSION(jpiglo) ::   tpol, fpol                         !: north fold mask (jperio= 3 or 4)
 
 #if defined key_noslip_accurate
-   INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:  ) :: npcoa        !: ???
-   INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: nicoa, njcoa !: ???
+   INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:  )  :: npcoa              !: ???
+   INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: nicoa, njcoa       !: ???
 #endif
 
@@ -279,22 +289,29 @@
          &      glamf(jpi,jpj) , gphif(jpi,jpj) , e1f(jpi,jpj) , e2f(jpi,jpj) , ff   (jpi,jpj) , STAT=ierr(3) )     
          !
-      ALLOCATE( gdep3w(jpi,jpj,jpk) , e3v(jpi,jpj,jpk) , e3f (jpi,jpj,jpk) ,                         &
-         &      gdept (jpi,jpj,jpk) , e3t(jpi,jpj,jpk) , e3u (jpi,jpj,jpk) ,                         &
-         &      gdepw (jpi,jpj,jpk) , e3w(jpi,jpj,jpk) , e3vw(jpi,jpj,jpk) , e3uw(jpi,jpj,jpk) , STAT=ierr(4) )
+      ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 (jpi,jpj,jpk) ,                         &
+         &      gdept_0 (jpi,jpj,jpk) , e3t_0(jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) ,                         &
+         &      gdepw_0 (jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , STAT=ierr(4) )
          !
 #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( gdept_1d(jpk) , gdepw_1d(jpk) ,                                   &
-         &      e3t_1d  (jpk) , e3w_1d  (jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) ,   &
-         &      gsigt  (jpk) , gsigw  (jpk) , gsi3w(jpk)    ,                     &
-         &      esigt  (jpk) , esigw  (jpk)                                 , STAT=ierr(7) )
+      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) ,                       &
+         &      re2u_e1u(jpi,jpj) , re1v_e2v(jpi,jpj) ,                                       &
+         &      e12t    (jpi,jpj) , r1_e12t (jpi,jpj) ,                                       &
+         &      e12u    (jpi,jpj) , r1_e12u (jpi,jpj) ,                                       &
+         &      e12v    (jpi,jpj) , r1_e12v (jpi,jpj) ,                                       &
+         &      e12f    (jpi,jpj) , r1_e12f (jpi,jpj) ,                                   STAT=ierr(6)  )
+         !
+      ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) ,                                     &
+         &      e3t_1d  (jpk) , e3w_1d  (jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) ,     &
+         &      gsigt   (jpk) , gsigw   (jpk) , gsi3w(jpk)    ,                     &
+         &      esigt   (jpk) , esigw   (jpk)                                 , STAT=ierr(7) )
          !
       ALLOCATE( hbatv (jpi,jpj) , hbatf (jpi,jpj) ,     &
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 3865)
@@ -86,13 +86,26 @@
                              CALL dom_msk      ! Masks
       IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency
-      IF( lk_vvl         )   CALL dom_vvl      ! Vertical variable mesh
+      IF( lk_vvl )           CALL dom_vvl_init ! Vertical variable mesh
       !
       IF( lk_c1d         )   CALL cor_c1d      ! 1D configuration: Coriolis set at T-point
+      !
+      ! - ML - Used in dom_vvl_sf_nxt and lateral diffusion routines
+      !        but could be usefull in many other routines
+      e12t    (:,:) = e1t(:,:) * e2t(:,:)
+      e12u    (:,:) = e1u(:,:) * e2u(:,:)
+      e12v    (:,:) = e1v(:,:) * e2v(:,:)
+      e12f    (:,:) = e1f(:,:) * e2f(:,:)
+      r1_e12t (:,:) = 1._wp    / e12t(:,:)
+      r1_e12u (:,:) = 1._wp    / e12u(:,:)
+      r1_e12v (:,:) = 1._wp    / e12v(:,:)
+      r1_e12f (:,:) = 1._wp    / e12f(:,:)
+      re2u_e1u(:,:) = e2u(:,:) / e1u(:,:)
+      re1v_e2v(:,:) = e1v(:,:) / e2v(:,:)
       !
       hu(:,:) = 0._wp                          ! Ocean depth at U- and V-points
       hv(:,:) = 0._wp
       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
@@ -345,20 +358,20 @@
          DO jj = 2, jpjm1
             DO jk = 1, jpkm1
-               zr1(1) = umask(ji-1,jj  ,jk) *abs( (gdepw(ji  ,jj  ,jk  )-gdepw(ji-1,jj  ,jk  )  & 
-                    &                         +gdepw(ji  ,jj  ,jk+1)-gdepw(ji-1,jj  ,jk+1)) &
-                    &                        /(gdepw(ji  ,jj  ,jk  )+gdepw(ji-1,jj  ,jk  )  &
-                    &                         -gdepw(ji  ,jj  ,jk+1)-gdepw(ji-1,jj  ,jk+1) + rsmall) )
-               zr1(2) = umask(ji  ,jj  ,jk) *abs( (gdepw(ji+1,jj  ,jk  )-gdepw(ji  ,jj  ,jk  )  &
-                    &                         +gdepw(ji+1,jj  ,jk+1)-gdepw(ji  ,jj  ,jk+1)) &
-                    &                        /(gdepw(ji+1,jj  ,jk  )+gdepw(ji  ,jj  ,jk  )  &
-                    &                         -gdepw(ji+1,jj  ,jk+1)-gdepw(ji  ,jj  ,jk+1) + rsmall) )
-               zr1(3) = vmask(ji  ,jj  ,jk) *abs( (gdepw(ji  ,jj+1,jk  )-gdepw(ji  ,jj  ,jk  )  &
-                    &                         +gdepw(ji  ,jj+1,jk+1)-gdepw(ji  ,jj  ,jk+1)) &
-                    &                        /(gdepw(ji  ,jj+1,jk  )+gdepw(ji  ,jj  ,jk  )  &
-                    &                         -gdepw(ji  ,jj+1,jk+1)-gdepw(ji  ,jj  ,jk+1) + rsmall) )
-               zr1(4) = vmask(ji  ,jj-1,jk) *abs( (gdepw(ji  ,jj  ,jk  )-gdepw(ji  ,jj-1,jk  )  &
-                    &                         +gdepw(ji  ,jj  ,jk+1)-gdepw(ji  ,jj-1,jk+1)) &
-                    &                        /(gdepw(ji  ,jj  ,jk  )+gdepw(ji  ,jj-1,jk  )  &
-                    &                         -gdepw(ji,  jj  ,jk+1)-gdepw(ji  ,jj-1,jk+1) + rsmall) )
+               zr1(1) = umask(ji-1,jj  ,jk) *abs( (gdepw_0(ji  ,jj  ,jk  )-gdepw_0(ji-1,jj  ,jk  )  & 
+                    &                         +gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji-1,jj  ,jk+1)) &
+                    &                        /(gdepw_0(ji  ,jj  ,jk  )+gdepw_0(ji-1,jj  ,jk  )  &
+                    &                         -gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji-1,jj  ,jk+1) + rsmall) )
+               zr1(2) = umask(ji  ,jj  ,jk) *abs( (gdepw_0(ji+1,jj  ,jk  )-gdepw_0(ji  ,jj  ,jk  )  &
+                    &                         +gdepw_0(ji+1,jj  ,jk+1)-gdepw_0(ji  ,jj  ,jk+1)) &
+                    &                        /(gdepw_0(ji+1,jj  ,jk  )+gdepw_0(ji  ,jj  ,jk  )  &
+                    &                         -gdepw_0(ji+1,jj  ,jk+1)-gdepw_0(ji  ,jj  ,jk+1) + rsmall) )
+               zr1(3) = vmask(ji  ,jj  ,jk) *abs( (gdepw_0(ji  ,jj+1,jk  )-gdepw_0(ji  ,jj  ,jk  )  &
+                    &                         +gdepw_0(ji  ,jj+1,jk+1)-gdepw_0(ji  ,jj  ,jk+1)) &
+                    &                        /(gdepw_0(ji  ,jj+1,jk  )+gdepw_0(ji  ,jj  ,jk  )  &
+                    &                         -gdepw_0(ji  ,jj+1,jk+1)-gdepw_0(ji  ,jj  ,jk+1) + rsmall) )
+               zr1(4) = vmask(ji  ,jj-1,jk) *abs( (gdepw_0(ji  ,jj  ,jk  )-gdepw_0(ji  ,jj-1,jk  )  &
+                    &                         +gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji  ,jj-1,jk+1)) &
+                    &                        /(gdepw_0(ji  ,jj  ,jk  )+gdepw_0(ji  ,jj-1,jk  )  &
+                    &                         -gdepw_0(ji,  jj  ,jk+1)-gdepw_0(ji  ,jj-1,jk+1) + rsmall) )
                zrxmax = MAXVAL(zr1(1:4))
                rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax)
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 3865)
@@ -6,16 +6,24 @@
    !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code
    !!            3.1  !  2009-02  (G. Madec, M. Leclair, R. Benshila)  pure z* coordinate
-   !!----------------------------------------------------------------------
-#if defined key_vvl
+   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl:
+   !!                                          vvl option includes z_star and z_tilde coordinates
    !!----------------------------------------------------------------------
    !!   '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_sf_nxt   : Compute next vertical scale factors
+   !!   dom_vvl_sf_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)
@@ -26,9 +34,25 @@
    PRIVATE
 
-   PUBLIC   dom_vvl         ! called by domain.F90
-   PUBLIC   dom_vvl_2       ! called by domain.F90
-   PUBLIC   dom_vvl_alloc   ! called by nemogcm.F90
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mut , muu , muv , muf    !: 1/H_0 at t-,u-,v-,f-points 
+   !! * 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  = .FALSE.   ! 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)                                              :: rn_ahe3       =  0.e0     ! thickness diffusion coefficient
+   LOGICAL , PUBLIC                                      :: ln_vvl_dbg    = .FALSE.   ! debug control prints
+
+   !! * Module variables
+   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_rst_e3t               ! retoring period for scale factors
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv               ! retoring period for low freq. divergence
 
    REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:)     ::   r2dt   ! vertical profile time-step, = 2 rdttra 
@@ -39,443 +63,734 @@
 #  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) ,     &
-         &      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_zstar ) dom_vvl_alloc = 0
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
+         ALLOCATE( e3t_t_b(jpi,jpj,jpk) , e3t_t_n(jpi,jpj,jpk) , e3t_t_a(jpi,jpj,jpk) ,   &
+            &      un_td  (jpi,jpj,jpk) , vn_td  (jpi,jpj,jpk) , STAT = dom_vvl_alloc        )
+         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
+         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
+      ENDIF
+      IF( ln_vvl_ztilde ) THEN
+         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
+         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 mu coefficients at t-, u-, v- and f-points to 
-      !!              spread ssh over the whole water column (scale factors)
-      !!                set the before and now ssh at u- and v-points 
-      !!              (also f-point in now case)
-      !!----------------------------------------------------------------------
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zcoefu, zcoefv , zcoeff                ! local scalars
-      REAL(wp) ::   zvt   , zvt_ip1, zvt_jp1, zvt_ip1jp1   !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:) ::  zee_t, zee_u, zee_v, zee_f   ! 2D workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('dom_vvl')
-      !
-      CALL wrk_alloc( jpi, jpj, zee_t, zee_u, zee_v, zee_f )
-      !
-      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  ==!
-      zee_t(:,:) = fse3t_0(:,:,1)                ! Lower bound : thickness of the first model level
-      zee_u(:,:) = fse3u_0(:,:,1)
-      zee_v(:,:) = fse3v_0(:,:,1)
-      zee_f(:,:) = fse3f_0(:,:,1)
-      DO jk = 2, jpkm1                          ! Sum of the masked vertical scale factors
-         zee_t(:,:) = zee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk)
-         zee_u(:,:) = zee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk)
-         zee_v(:,:) = zee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk)
-         DO jj = 1, jpjm1                      ! f-point : fmask=shlat at coasts, use the product of umask
-            zee_f(:,jj) = zee_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
-      zee_t(:,:) = 1._wp / zee_t(:,:) * tmask(:,:,1)
-      zee_u(:,:) = 1._wp / zee_u(:,:) * umask(:,:,1)
-      zee_v(:,:) = 1._wp / zee_v(:,:) * vmask(:,:,1)
-      DO jj = 1, jpjm1                               ! f-point case fmask cannot be used 
-         zee_f(:,jj) = 1._wp / zee_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(t/u/v)_0
+      !!              - frq_rst_e3t and frq_rst_hdv
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      USE phycst,  ONLY : rpi
+      !! * Local declarations
+      INTEGER ::   jk
+      !!----------------------------------------------------------------------
+      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_init')
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
+
+      ! 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 fse3t_(b/n), e3t_t_(b/n) and hdiv_lf (and e3t_a(jpk))
+      ! ========================================================================
+      CALL dom_vvl_rst( nit000, 'READ' )
+      fse3t_a(:,:,jpk) = e3t_0(:,:,jpk)
+
+      ! 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( zee_f, 'F', 1. )                 ! lateral boundary condition on ee_f
-      !
-      DO jk = 1, jpk                            ! mu coefficients
-         mut(:,:,jk) = zee_t(:,:) * tmask(:,:,jk)     ! T-point at T levels
-         muu(:,:,jk) = zee_u(:,:) * umask(:,:,jk)     ! U-point at T levels
-         muv(:,:,jk) = zee_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) = zee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk)   ! at T levels
-         END DO
-         muf(:,jpj,jk) = 0._wp
-      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
-         hu_0(:,:) = hu_0(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk)
-         hv_0(:,:) = hv_0(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk)
+         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
+         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
+         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
       END DO
-      
-      DO jj = 1, jpjm1                          ! initialise before and now Sea Surface Height at u-, v-, f-points
-         DO ji = 1, jpim1   ! NO vector opt.
-            zcoefu = 0.50_wp / ( e1u(ji,jj) * e2u(ji,jj) ) * umask(ji,jj,1)
-            zcoefv = 0.50_wp / ( e1v(ji,jj) * e2v(ji,jj) ) * vmask(ji,jj,1)
-            zcoeff = 0.25_wp / ( e1f(ji,jj) * e2f(ji,jj) ) * umask(ji,jj,1) * umask(ji,jj+1,1)
+
+      ! Restoring frequencies for z_tilde coordinate
+      ! ============================================
+      IF( ln_vvl_ztilde ) THEN
+         ! - ML - In the future, this should be tunable by the user (namelist)
+         frq_rst_e3t(:,:) = 2.e0_wp * rpi / ( 30.e0_wp * 86400.e0_wp )
+         frq_rst_hdv(:,:) = 2.e0_wp * rpi / (  5.e0_wp * 86400.e0_wp )
+      ENDIF
+
+      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_init')
+
+   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_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
+      !!               - z_tilde_case: after scale factor increment = 
+      !!                                    high frequency part of horizontal divergence
+      !!                                  + retsoring towards the background grid
+      !!                                  + thickness difusion
+      !!                               Then repartition of ssh INCREMENT proportionnaly
+      !!                               to the "baroclinic" 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 Madec, G. 2011, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zht, z_scale, zwu, zwv, zhdiv
+      !! * 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
+      !!----------------------------------------------------------------------
+      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt')
+      CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
+      CALL wrk_alloc( jpi, jpj, jpk, ze3t                     )
+
+      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_zstar ) THEN                          ! z_star 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(:,:) * tmask(:,:,jk)
+         END DO
+
+      !                                                ! --------------------------- !
+      ELSEIF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN   ! z_tilde or layer coordinate !
+         !                                             ! --------------------------- !
+
+         ! I - initialization
+         ! ==================
+
+         ! 1 - barotropic divergence
+         ! -------------------------
+         zhdiv(:,:) = 0.
+         zht(:,:)   = 0.
+         DO jk = 1, jpkm1
+            zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
+            zht  (:,:) = zht  (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+         END DO
+         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask(:,:,1) )
+
+         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
+         ! --------------------------------------------------
+         IF( ln_vvl_ztilde ) THEN
+            IF( kt .GT. nit000 ) THEN
+               DO jk = 1, jpkm1
+                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
+                     &          * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
+               END DO
+            ENDIF
+         END IF
+
+         ! II - 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
+         ! ----------------------------------
+         IF( ln_vvl_ztilde ) THEN   ! z_tilde case
+            DO jk = 1, jpkm1
+               e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
+            END DO
+         ELSE                       ! layer case
+            DO jk = 1, jpkm1
+               e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) )
+            END DO
+         END IF
+
+         ! 2 - Restoring term (z-tilde case only)
+         ! ------------------
+         IF( ln_vvl_ztilde ) THEN
+            DO jk = 1, jpk
+               e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - frq_rst_e3t(:,:) * e3t_t_b(:,:,jk)
+            END DO
+         END IF
+
+         ! 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) = rn_ahe3 * umask(ji,jj,jk) * re2u_e1u(ji,jj) * ( e3t_t_b(ji,jj,jk) - e3t_t_b(ji+1,jj  ,jk) )
+                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(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) )  &
+                     &                                    * r1_e12t(ji,jj)
+               END DO
+            END DO
+         END DO
+         ! d - thickness diffusion transport: boundary conditions
+         !     (stored for tracer advction 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. )
+         e3t_t_a(:,:,:) = e3t_t_b(:,:,:) + z2dt * tmask(:,:,:) * 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) / e3t_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(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
+            WRITE(numout, *) 'at i, j, k=', ijk_max
+            WRITE(numout, *) 'MIN( e3t_t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
+            WRITE(numout, *) 'at i, j, k=', ijk_min            
+            CALL ctl_stop('MAX( ABS( e3t_t_a(:,:,:) ) / e3t_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(:,:,:),   z_def_max * e3t_0(:,:,:) )
+         e3t_t_a(:,:,:) = MAX( e3t_t_a(:,:,:), - z_def_max * e3t_0(:,:,:) )
+
+         ! Add "tilda" part to the after scale factor
+         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         fse3t_a(:,:,:) = e3t_0(:,:,:) + e3t_t_a(:,:,:)
+
+         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
+         ! ==================================================================================
+         ! add e3t(n-1) "star" Asselin-filtered
+         DO jk = 1, jpkm1
+            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_b(:,:,jk) - e3t_0(:,:,jk) - e3t_t_b(:,:,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.
+         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
+         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
+
+      ENDIF
+
+      IF( ln_vvl_dbg ) THEN   ! - ML - test: control prints for debuging
+         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
+            WRITE(numout, *) 'kt =', kt
+            WRITE(numout, *) 'MAXVAL(abs(SUM(e3t_t_a))) =',   &
+               &              MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( 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+ssha-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' )
+
+      CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
+      CALL wrk_dealloc( jpi, jpj, jpk, ze3t                     )
+
+      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt')
+
+   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.
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( in )               :: kt       ! time step
+      !! * Local declarations
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_e3t_def
+      INTEGER                             :: jk       ! dummy loop indices
+      !!----------------------------------------------------------------------
+
+      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_swp')
+      !
+      CALL wrk_alloc( jpi, jpj, jpk, z_e3t_def                )
+      !
+      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
+      ! 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. - vmask(:,:,1) )
+
+      ! Write outputs
+      ! =============
+      z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
+      CALL iom_put( "e3t_n"  , fse3t_n  (:,:,:) )
+      CALL iom_put( "dept_n" , fsde3w_n (:,:,:) )
+      CALL iom_put( "e3tdef" , z_e3t_def(:,:,:) )
+
+      ! write restart file
+      ! ==================
+      IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
+      !
+      CALL wrk_dealloc( jpi, jpj, jpk, z_e3t_def )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_swp')
+
+   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
+      !!----------------------------------------------------------------------
+      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_interpol')
+      SELECT CASE ( pout )
+         !               ! ------------------------------------- !
+      CASE( 'U' )        ! interpolation from T-point to U-point !
+         !               ! ------------------------------------- !
+         ! horizontal surface weighted interpolation
+         DO jk = 1, jpk
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  pe3_out(ji,jj,jk) = 0.5 * umask(ji,jj,jk) * r1_e12u(ji,jj)                                      &
+                     &                    * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
+                     &                        + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
+               END DO
+            END DO
+         END DO
+         ! boundary conditions
+         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1. )
+         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
+         !               ! ------------------------------------- !
+      CASE( 'V' )        ! interpolation from T-point to V-point !
+         !               ! ------------------------------------- !
+         ! horizontal surface weighted interpolation
+         DO jk = 1, jpk
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  pe3_out(ji,jj,jk) = 0.5 * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                      &
+                     &                    * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
+                     &                        + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
+               END DO
+            END DO
+         END DO
+         ! boundary conditions
+         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1. )
+         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
+         !               ! ------------------------------------- !
+      CASE( 'F' )        ! interpolation from U-point to F-point !
+         !               ! ------------------------------------- !
+         ! horizontal surface weighted interpolation
+         DO jk = 1, jpk
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  pe3_out(ji,jj,jk) = 0.5 * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)                  &
+                     &                    * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
+                     &                        + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
+               END DO
+            END DO
+         END DO
+         ! boundary conditions
+         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1. )
+         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
+         !               ! ------------------------------------- !
+      CASE( 'W' )        ! interpolation from T-point to W-point !
+         !               ! ------------------------------------- !
+         ! vertical simple interpolation
+         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
+         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
+         DO jk = 2, jpk
+            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1. - 0.5 * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   &
+               &                              +        0.5 * tmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
+         END DO
+         !               ! -------------------------------------- !
+      CASE( 'UW' )       ! interpolation from U-point to UW-point !
+         !               ! -------------------------------------- !
+         ! vertical simple interpolation
+         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
+         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
+         DO jk = 2, jpk
+            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1. - 0.5 * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   &
+               &                               +        0.5 * umask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
+         END DO
+         !               ! -------------------------------------- !
+      CASE( 'VW' )       ! interpolation from V-point to VW-point !
+         !               ! -------------------------------------- !
+         ! vertical simple interpolation
+         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
+         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
+         DO jk = 2, jpk
+            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1. - 0.5 * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   &
+               &                               +        0.5 * vmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
+         END DO
+      END SELECT
+
+      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_interpol')
+
+   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( nn_timing == 1 )  CALL timing_start('dom_vvl_rst')
+      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. )
+            !                             ! --------- !
+            !                             ! all cases !
+            !                             ! --------- !
+            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( 'dom_vvl_rst: vvl cannot restart from a non vvl run' )
+            ENDIF
+            !                             ! ----------- !
+            IF( ln_vvl_zstar ) THEN       ! z_star case !
+               !                          ! ----------- !
+               IF( MIN( id3, id4 ) > 0 ) THEN
+                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
+               ENDIF
+               !                          ! ----------------------- !
+            ELSE                          ! z_tilde and layer cases !
+               !                          ! ----------------------- !
+               IF( MIN( id3, id4 ) > 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(:,:,:) )
+               ELSE                            ! one at least array is missing
+                  e3t_t_b(:,:,:) = 0.e0
+                  e3t_t_n(:,:,:) = 0.e0
+               ENDIF
+               !                          ! ------------ !
+               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
+                  !                       ! ------------ !
+                  IF( id5 > 0 ) THEN  ! required array exists
+                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
+                  ELSE                ! array is missing
+                     hdiv_lf(:,:,:) = 0.e0
+                  ENDIF
+               ENDIF
+            ENDIF
             !
-            zvt           = e1e2t(ji  ,jj  ) * sshb(ji  ,jj  )    ! before fields
-            zvt_ip1       = e1e2t(ji+1,jj  ) * sshb(ji+1,jj  )
-            zvt_jp1       = e1e2t(ji  ,jj+1) * sshb(ji  ,jj+1)
-            sshu_b(ji,jj) = zcoefu * ( zvt + zvt_ip1 )
-            sshv_b(ji,jj) = zcoefv * ( zvt + zvt_jp1 )
-            !
-            zvt           = e1e2t(ji  ,jj  ) * sshn(ji  ,jj  )    ! now fields
-            zvt_ip1       = e1e2t(ji+1,jj  ) * sshn(ji+1,jj  )
-            zvt_jp1       = e1e2t(ji  ,jj+1) * sshn(ji  ,jj+1)
-            zvt_ip1jp1    = e1e2t(ji+1,jj+1) * sshn(ji+1,jj+1)
-            sshu_n(ji,jj) = zcoefu * ( zvt + zvt_ip1 )
-            sshv_n(ji,jj) = zcoefv * ( zvt + zvt_jp1 )
-            sshf_n(ji,jj) = zcoeff * ( zvt + zvt_ip1 + zvt_jp1 + zvt_ip1jp1 )
-         END DO
-      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. )
-      !
-      CALL wrk_dealloc( jpi, jpj, zee_t, zee_u, zee_v, zee_f )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl')
-      !
-   END SUBROUTINE dom_vvl
-
-
-   SUBROUTINE dom_vvl_2( kt, pe3u_b, pe3v_b )
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE dom_vvl_2  ***
-      !!                   
-      !! ** Purpose :   compute the vertical scale factors at u- and v-points
-      !!              in variable volume case.
-      !!
-      !! ** Method  :   In variable volume case (non linear sea surface) the 
-      !!              the vertical scale factor at velocity points is computed
-      !!              as the average of the cell surface weighted e3t.
-      !!                It uses the sea surface heigth so it have to be initialized
-      !!              after ssh is read/set
-      !!----------------------------------------------------------------------
-      INTEGER                   , INTENT(in   ) ::   kt               ! ocean time-step index
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe3u_b, pe3v_b   ! before vertical scale factor at u- & v-pts
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      INTEGER  ::   iku, ikv     ! local integers    
-      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
-      REAL(wp) ::   zvt          ! local scalars
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_2')
-      !
-      IF( lwp .AND. kt == nit000 ) THEN
+         ELSE                                   !* Initialize at "rest"
+            fse3t_b(:,:,:) = e3t_0(:,:,:)
+            fse3t_n(:,:,:) = e3t_0(:,:,:)
+            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
+               e3t_t_b(:,:,:) = 0.e0
+               e3t_t_n(:,:,:) = 0.e0
+               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.e0
+            END IF
+         ENDIF
+
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+         !                                   ! ===================
+         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
+         !                                           ! --------- !
+         !                                           ! all 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 and layer cases !
+            !                                        ! ----------------------- !
+            CALL iom_rstput( kt, nitrst, numrow, 'e3t_t_b', e3t_t_b(:,:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'e3t_t_n', e3t_t_n(:,:,:) )
+         END IF
+         !                                           ! -------------!    
+         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
+            !                                        ! ------------ !
+            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )
+         ENDIF
+
+      ENDIF
+      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_rst')
+
+   END SUBROUTINE dom_vvl_rst
+
+
+   SUBROUTINE dom_vvl_ctl
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_vvl_ctl  ***
+      !!                
+      !! ** Purpose :   Control the consistency between namelist options
+      !!                for vertical coordinate
+      !!----------------------------------------------------------------------
+      INTEGER ::   ioptio
+
+      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, rn_ahe3, ln_vvl_dbg! , 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_2 : Variable volume, fse3t_b initialization'
-         WRITE(numout,*) '~~~~~~~~~ '
-         pe3u_b(:,:,jpk) = fse3u_0(:,:,jpk)
-         pe3v_b(:,:,jpk) = fse3u_0(:,:,jpk)
-      ENDIF
-      
-      DO jk = 1, jpkm1           ! set the before scale factors at u- & v-points
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zvt = fse3t_b(ji,jj,jk) * e1e2t(ji,jj)
-               pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1e2t(ji+1,jj) ) / ( e1u(ji,jj) * e2u(ji,jj) )
-               pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e1e2t(ji,jj+1) ) / ( e1v(ji,jj) * e2v(ji,jj) )
-            END DO
-         END DO
-      END DO
-
-      ! Correct scale factors at locations that have been individually modified in domhgr
-      ! Such modifications break the relationship between e1e2t and e1u*e2u etc. Recompute
-      ! scale factors ignoring the modified metric.
-      !                                                ! =====================
-      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
-         !                                             ! =====================
-         IF( nn_cla == 0 ) THEN
-            !
-            ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u was modified)
-            ij0 = 102   ;   ij1 = 102   
-            DO jk = 1, jpkm1                 ! set the before scale factors at u-points
-               DO jj = mj0(ij0), mj1(ij1)
-                  DO ji = mi0(ii0), mi1(ii1)
-                     zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
-                     pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
-                  END DO
-               END DO
-            END DO
-            !
-            ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u and e1v were modified)
-            ij0 =  88   ;   ij1 =  88   
-            DO jk = 1, jpkm1                 ! set the before scale factors at u-points
-               DO jj = mj0(ij0), mj1(ij1)
-                  DO ji = mi0(ii0), mi1(ii1)
-                     zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
-                     pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
-                  END DO
-               END DO
-            END DO
-            DO jk = 1, jpkm1                 ! set the before scale factors at v-points
-               DO jj = mj0(ij0), mj1(ij1)
-                  DO ji = mi0(ii0), mi1(ii1)
-                     zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
-                     pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
-                  END DO
-               END DO
-            END DO
-         ENDIF
-
-         ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u was modified)
-         ij0 = 116   ;   ij1 = 116   
-         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
-            DO jj = mj0(ij0), mj1(ij1)
-               DO ji = mi0(ii0), mi1(ii1)
-                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
-                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
-               END DO
-            END DO
-         END DO
-         !
-      ENDIF
-         !                                             ! =====================
-      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
-         !                                             ! =====================
-
-         ii0 = 281   ;   ii1 = 282        ! Gibraltar Strait (e2u was modified)
-         ij0 = 200   ;   ij1 = 200   
-         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
-            DO jj = mj0(ij0), mj1(ij1)
-               DO ji = mi0(ii0), mi1(ii1)
-                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
-                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
-               END DO
-            END DO
-         END DO
-
-         ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait (e2u was modified)
-         ij0 = 208   ;   ij1 = 208   
-         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
-            DO jj = mj0(ij0), mj1(ij1)
-               DO ji = mi0(ii0), mi1(ii1)
-                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
-                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
-               END DO
-            END DO
-         END DO
-
-         ii0 =  44   ;   ii1 =  44        ! Lombok Strait (e1v was modified)
-         ij0 = 124   ;   ij1 = 125   
-         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
-            DO jj = mj0(ij0), mj1(ij1)
-               DO ji = mi0(ii0), mi1(ii1)
-                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
-                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
-               END DO
-            END DO
-         END DO
-
-         ii0 =  48   ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
-         ij0 = 124   ;   ij1 = 125   
-         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
-            DO jj = mj0(ij0), mj1(ij1)
-               DO ji = mi0(ii0), mi1(ii1)
-                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
-                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
-               END DO
-            END DO
-         END DO
-
-         ii0 =  53   ;   ii1 =  53        ! Ombai Strait (e1v was modified)
-         ij0 = 124   ;   ij1 = 125   
-         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
-            DO jj = mj0(ij0), mj1(ij1)
-               DO ji = mi0(ii0), mi1(ii1)
-                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
-                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
-               END DO
-            END DO
-         END DO
-
-         ii0 =  56   ;   ii1 =  56        ! Timor Passage (e1v was modified)
-         ij0 = 124   ;   ij1 = 125   
-         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
-            DO jj = mj0(ij0), mj1(ij1)
-               DO ji = mi0(ii0), mi1(ii1)
-                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
-                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
-               END DO
-            END DO
-         END DO
-
-         ii0 =  55   ;   ii1 =  55        ! West Halmahera Strait (e1v was modified)
-         ij0 = 141   ;   ij1 = 142   
-         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
-            DO jj = mj0(ij0), mj1(ij1)
-               DO ji = mi0(ii0), mi1(ii1)
-                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
-                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
-               END DO
-            END DO
-         END DO
-
-         ii0 =  58   ;   ii1 =  58        ! East Halmahera Strait (e1v was modified)
-         ij0 = 141   ;   ij1 = 142   
-         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
-            DO jj = mj0(ij0), mj1(ij1)
-               DO ji = mi0(ii0), mi1(ii1)
-                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
-                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
-               END DO
-            END DO
-         END DO
-
-         !
-      ENDIF
-      !                                                ! ======================
-      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
-         !                                             ! ======================
-         ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u was modified)
-         ij0 = 327   ;   ij1 = 327   
-         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
-            DO jj = mj0(ij0), mj1(ij1)
-               DO ji = mi0(ii0), mi1(ii1)
-                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
-                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
-               END DO
-            END DO
-         END DO
-         !
-         ii0 = 627   ;   ii1 = 628        ! Bosphore Strait (e2u was modified)
-         ij0 = 343   ;   ij1 = 343   
-         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
-            DO jj = mj0(ij0), mj1(ij1)
-               DO ji = mi0(ii0), mi1(ii1)
-                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
-                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
-               END DO
-            END DO
-         END DO
-         !
-         ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u was modified)
-         ij0 = 232   ;   ij1 = 232   
-         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
-            DO jj = mj0(ij0), mj1(ij1)
-               DO ji = mi0(ii0), mi1(ii1)
-                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
-                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
-               END DO
-            END DO
-         END DO
-         !
-         ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u was modified)
-         ij0 = 232   ;   ij1 = 232   
-         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
-            DO jj = mj0(ij0), mj1(ij1)
-               DO ji = mi0(ii0), mi1(ii1)
-                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
-                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
-               END DO
-            END DO
-         END DO
-         !
-         ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u was modified)
-         ij0 = 270   ;   ij1 = 270   
-         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
-            DO jj = mj0(ij0), mj1(ij1)
-               DO ji = mi0(ii0), mi1(ii1)
-                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
-                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
-               END DO
-            END DO
-         END DO
-         !
-         ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v was modified)
-         ij0 = 232   ;   ij1 = 233   
-         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
-            DO jj = mj0(ij0), mj1(ij1)
-               DO ji = mi0(ii0), mi1(ii1)
-                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
-                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
-               END DO
-            END DO
-         END DO
-         !
-         ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v was modified)
-         ij0 = 276   ;   ij1 = 276   
-         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
-            DO jj = mj0(ij0), mj1(ij1)
-               DO ji = mi0(ii0), mi1(ii1)
-                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
-                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
-               END DO
-            END DO
-         END DO
-         !
-      ENDIF
-      ! End of individual corrections to scale factors
-
-      IF( ln_zps ) THEN          ! minimum of the e3t at partial cell level
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               iku = mbku(ji,jj)
-               ikv = mbkv(ji,jj)
-               pe3u_b(ji,jj,iku) = MIN( fse3t_b(ji,jj,iku), fse3t_b(ji+1,jj  ,iku) ) 
-               pe3v_b(ji,jj,ikv) = MIN( fse3t_b(ji,jj,ikv), fse3t_b(ji  ,jj+1,ikv) ) 
-            END DO
-         END DO
-      ENDIF
-
-      pe3u_b(:,:,:) = pe3u_b(:,:,:) - fse3u_0(:,:,:)      ! anomaly to avoid zero along closed boundary/extra halos
-      pe3v_b(:,:,:) = pe3v_b(:,:,:) - fse3v_0(:,:,:)
-      CALL lbc_lnk( pe3u_b(:,:,:), 'U', 1. )               ! lateral boundary conditions
-      CALL lbc_lnk( pe3v_b(:,:,:), 'V', 1. )
-      pe3u_b(:,:,:) = pe3u_b(:,:,:) + fse3u_0(:,:,:)      ! recover the full scale factor
-      pe3v_b(:,:,:) = pe3v_b(:,:,:) + fse3v_0(:,:,:)
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_2')
-      !
-   END SUBROUTINE dom_vvl_2
-   
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :                                      Empty routine
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE dom_vvl
-   END SUBROUTINE dom_vvl
-   SUBROUTINE dom_vvl_2(kdum, pudum, pvdum )
-      USE par_kind
-      INTEGER                   , INTENT(in   ) ::   kdum       
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pudum, pvdum
-   END SUBROUTINE dom_vvl_2
-#endif
+         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,*) '                                         rn_ahe3        = ', rn_ahe3
+         WRITE(numout,*) '           Namelist nam_vvl : debug prints'
+         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg
+      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(lwp) THEN                   ! Print the choice
+         WRITE(numout,*)
+         IF( ln_vvl_zstar      ) WRITE(numout,*) '              zstar vertical coordinate is used'
+         IF( ln_vvl_ztilde     ) WRITE(numout,*) '              ztilde vertical coordinate is used'
+         IF( ln_vvl_layer      ) 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
 
    !!======================================================================
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 3865)
@@ -183,12 +183,12 @@
          CALL iom_rstput( 0, 0, inum4, 'esigw', esigw )
          !
-         CALL iom_rstput( 0, 0, inum4, 'e3t', e3t )             !    ! scale factors
-         CALL iom_rstput( 0, 0, inum4, 'e3u', e3u )
-         CALL iom_rstput( 0, 0, inum4, 'e3v', e3v )
-         CALL iom_rstput( 0, 0, inum4, 'e3w', e3w )
+         CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 )         !    ! scale factors
+         CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 )
+         CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 )
+         CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 )
          CALL iom_rstput( 0, 0, inum4, 'rx1', rx1 )             !    ! Max. grid stiffness ratio
          !
-         CALL iom_rstput( 0, 0, inum4, 'gdept' , gdept )    !    ! stretched system
-         CALL iom_rstput( 0, 0, inum4, 'gdepw' , gdepw )
+         CALL iom_rstput( 0, 0, inum4, 'gdept_1d' , gdept_1d )  !    ! stretched system
+         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d' , gdepw_1d )
       ENDIF
       
@@ -196,13 +196,13 @@
          !
          IF( nmsh <= 6 ) THEN                                   !    ! 3D vertical scale factors
-            CALL iom_rstput( 0, 0, inum4, 'e3t', e3t )         
-            CALL iom_rstput( 0, 0, inum4, 'e3u', e3u )
-            CALL iom_rstput( 0, 0, inum4, 'e3v', e3v )
-            CALL iom_rstput( 0, 0, inum4, 'e3w', e3w )
+            CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 )         
+            CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 )
+            CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 )
+            CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 )
          ELSE                                                   !    ! 2D masked bottom ocean scale factors
             DO jj = 1,jpj   
                DO ji = 1,jpi
-                  e3tp(ji,jj) = e3t(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)
-                  e3wp(ji,jj) = e3w(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)
+                  e3tp(ji,jj) = e3t_0(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)
+                  e3wp(ji,jj) = e3w_0(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)
                END DO
             END DO
@@ -212,10 +212,10 @@
          !
          IF( nmsh <= 3 ) THEN                                   !    ! 3D depth
-            CALL iom_rstput( 0, 0, inum4, 'gdept', gdept, ktype = jp_r4 )     
+            CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r4 )     
             DO jk = 1,jpk   
                DO jj = 1, jpjm1   
                   DO ji = 1, fs_jpim1   ! vector opt.
-                     zdepu(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji+1,jj  ,jk) )
-                     zdepv(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji  ,jj+1,jk) )
+                     zdepu(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk) )
+                     zdepv(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk) )
                   END DO   
                END DO   
@@ -224,10 +224,10 @@
             CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 )
             CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r4 )
-            CALL iom_rstput( 0, 0, inum4, 'gdepw', gdepw, ktype = jp_r4 )
+            CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r4 )
          ELSE                                                   !    ! 2D bottom depth
             DO jj = 1,jpj   
                DO ji = 1,jpi
-                  zprt(ji,jj) = gdept(ji,jj,mbkt(ji,jj)  ) * tmask(ji,jj,1)
-                  zprw(ji,jj) = gdepw(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1)
+                  zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj)  ) * tmask(ji,jj,1)
+                  zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1)
                END DO
             END DO
@@ -246,5 +246,5 @@
          CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )   !    ! depth
          CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
-         CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )     !    ! scale factors
+         CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )   !    ! scale factors
          CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   )
       ENDIF
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 3865)
@@ -146,17 +146,17 @@
       IF( nprint == 1 .AND. lwp )   THEN
          WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
-         WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ),   &
-            &                   ' w ',   MINVAL( fsdepw(:,:,:) ), '3w ', MINVAL( fsde3w(:,:,:) )
-         WRITE(numout,*) ' MIN val e3    t ', MINVAL( fse3t(:,:,:) ), ' f ', MINVAL( fse3f(:,:,:) ),  &
-            &                   ' u ',   MINVAL( fse3u(:,:,:) ), ' u ', MINVAL( fse3v(:,:,:) ),  &
-            &                   ' uw',   MINVAL( fse3uw(:,:,:)), ' vw', MINVAL( fse3vw(:,:,:)),   &
-            &                   ' w ',   MINVAL( fse3w(:,:,:) )
-
-         WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ),   &
-            &                   ' w ',   MAXVAL( fsdepw(:,:,:) ), '3w ', MAXVAL( fsde3w(:,:,:) )
-         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( fse3t(:,:,:) ), ' f ', MAXVAL( fse3f(:,:,:) ),  &
-            &                   ' u ',   MAXVAL( fse3u(:,:,:) ), ' u ', MAXVAL( fse3v(:,:,:) ),  &
-            &                   ' uw',   MAXVAL( fse3uw(:,:,:)), ' vw', MAXVAL( fse3vw(:,:,:)),   &
-            &                   ' w ',   MAXVAL( fse3w(:,:,:) )
+         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
+            &                   ' w ',   MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gdep3w_0(:,:,:) )
+         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ),  &
+            &                   ' u ',   MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ),  &
+            &                   ' uw',   MINVAL( e3uw_0(:,:,:)), ' vw', MINVAL( e3vw_0(:,:,:)),   &
+            &                   ' w ',   MINVAL( e3w_0(:,:,:) )
+
+         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
+            &                   ' w ',   MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gdep3w_0(:,:,:) )
+         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ),  &
+            &                   ' u ',   MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ),  &
+            &                   ' uw',   MAXVAL( e3uw_0(:,:,:)), ' vw', MAXVAL( e3vw_0(:,:,:)),   &
+            &                   ' w ',   MAXVAL( e3w_0(:,:,:) )
       ENDIF
       !
@@ -177,8 +177,8 @@
       !!      function the derivative of which gives the scale factors.
       !!        both depth and scale factors only depend on k (1d arrays).
-      !!              w-level: gdepw_1d  = fsdep(k)
-      !!                       e3w_1d(k) = dk(fsdep)(k)     = fse3(k)
-      !!              t-level: gdept_1d  = fsdep(k+0.5)
-      !!                       e3t_1d(k) = dk(fsdep)(k+0.5) = fse3(k+0.5)
+      !!              w-level: gdepw_1d  = gdep(k)
+      !!                       e3w_1d(k) = dk(gdep)(k)     = e3(k)
+      !!              t-level: gdept_1d  = gdep(k+0.5)
+      !!                       e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5)
       !!
       !! ** Action  : - gdept_1d, gdepw_1d : depth of T- and W-point (m)
@@ -298,10 +298,10 @@
          WRITE(numout,*)
          WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
-         WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" )
+         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
          WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
       ENDIF
       DO jk = 1, jpk                      ! control positivity
-         IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( 'dom:zgr_z: e3w or e3t =< 0 '    )
-         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( 'dom:zgr_z: gdepw or gdept < 0 ' )
+         IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 '    )
+         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' )
       END DO
       !
@@ -369,5 +369,5 @@
             IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin'
             idta(:,:) = jpkm1                            ! before last level
-            zdta(:,:) = gdepw_1d(jpk)                     ! last w-point depth
+            zdta(:,:) = gdepw_1d(jpk)                    ! last w-point depth
             h_oce     = gdepw_1d(jpk)
          ELSE                                         ! bump centered in the basin
@@ -378,5 +378,5 @@
             r_bump  = 50000._wp                            ! bump radius (meters)       
             h_bump  =  2700._wp                            ! bump height (meters)
-            h_oce   = gdepw_1d(jpk)                         ! background ocean depth (meters)
+            h_oce   = gdepw_1d(jpk)                        ! background ocean depth (meters)
             IF(lwp) WRITE(numout,*) '            bump characteristics: '
             IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ii_bump
@@ -440,7 +440,7 @@
             CALL iom_close( inum )
             mbathy(:,:) = INT( bathy(:,:) )
-            !
+            !                                                ! =====================
             IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
-               !
+               !                                             ! =====================
                IF( nn_cla == 0 ) THEN
                   ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open 
@@ -780,14 +780,14 @@
       !
       DO jk = 1, jpk
-            gdept(:,:,jk) = gdept_1d(jk)
-            gdepw(:,:,jk) = gdepw_1d(jk)
-            gdep3w(:,:,jk) = gdepw_1d(jk)
-            e3t (:,:,jk) = e3t_1d(jk)
-            e3u (:,:,jk) = e3t_1d(jk)
-            e3v (:,:,jk) = e3t_1d(jk)
-            e3f (:,:,jk) = e3t_1d(jk)
-            e3w (:,:,jk) = e3w_1d(jk)
-            e3uw(:,:,jk) = e3w_1d(jk)
-            e3vw(:,:,jk) = e3w_1d(jk)
+         gdept_0 (:,:,jk) = gdept_1d(jk)
+         gdepw_0 (:,:,jk) = gdepw_1d(jk)
+         gdep3w_0(:,:,jk) = gdepw_1d(jk)
+         e3t_0   (:,:,jk) = e3t_1d  (jk)
+         e3u_0   (:,:,jk) = e3t_1d  (jk)
+         e3v_0   (:,:,jk) = e3t_1d  (jk)
+         e3f_0   (:,:,jk) = e3t_1d  (jk)
+         e3w_0   (:,:,jk) = e3w_1d  (jk)
+         e3uw_0  (:,:,jk) = e3w_1d  (jk)
+         e3vw_0  (:,:,jk) = e3w_1d  (jk)
       END DO
       !
@@ -814,8 +814,8 @@
       !!      with partial steps  on 3d arrays ( i, j, k ).
       !!
-      !!              w-level: gdepw(i,j,k)  = fsdep(k)
-      !!                       e3w(i,j,k) = dk(fsdep)(k)     = fse3(i,j,k)
-      !!              t-level: gdept(i,j,k)  = fsdep(k+0.5)
-      !!                       e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5)
+      !!              w-level: gdepw_0(i,j,k)  = gdep(k)
+      !!                       e3w_0(i,j,k) = dk(gdep)(k)     = e3(i,j,k)
+      !!              t-level: gdept_0(i,j,k)  = gdep(k+0.5)
+      !!                       e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
       !!
       !!        With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
@@ -825,5 +825,5 @@
       !!              - bathy = 0 => mbathy = 0
       !!              - 1 < mbathy < jpkm1    
-      !!              - bathy > gdepw(jpk) => mbathy = jpkm1  
+      !!              - bathy > gdepw_0(jpk) => mbathy = jpkm1  
       !!
       !!        Then, for each case, we find the new depth at t- and w- levels
@@ -838,5 +838,5 @@
       !!
       !!         c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives
-      !!         - - - - - - -   gdept, gdepw and e3. are positives
+      !!         - - - - - - -   gdept_0, gdepw_0 and e3. are positives
       !!      
       !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
@@ -891,8 +891,8 @@
       ! Scale factors and depth at T- and W-points
       DO jk = 1, jpk                        ! intitialization to the reference z-coordinate
-         gdept(:,:,jk) = gdept_1d(jk)
-         gdepw(:,:,jk) = gdepw_1d(jk)
-         e3t  (:,:,jk) = e3t_1d  (jk)
-         e3w  (:,:,jk) = e3w_1d  (jk)
+         gdept_0(:,:,jk) = gdept_1d(jk)
+         gdepw_0(:,:,jk) = gdepw_1d(jk)
+         e3t_0  (:,:,jk) = e3t_1d  (jk)
+         e3w_0  (:,:,jk) = e3w_1d  (jk)
       END DO
       ! 
@@ -906,29 +906,29 @@
                   ze3tp  = bathy(ji,jj) - gdepw_1d(ik)
                   ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
-                  e3t(ji,jj,ik  ) = ze3tp
-                  e3t(ji,jj,ik+1) = ze3tp
-                  e3w(ji,jj,ik  ) = ze3wp
-                  e3w(ji,jj,ik+1) = ze3tp
-                  gdepw(ji,jj,ik+1) = zdepwp
-                  gdept(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
-                  gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
+                  e3t_0(ji,jj,ik  ) = ze3tp
+                  e3t_0(ji,jj,ik+1) = ze3tp
+                  e3w_0(ji,jj,ik  ) = ze3wp
+                  e3w_0(ji,jj,ik+1) = ze3tp
+                  gdepw_0(ji,jj,ik+1) = zdepwp
+                  gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
+                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
                   !
                ELSE                         ! standard case
-                  IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw(ji,jj,ik+1) = bathy(ji,jj)
-                  ELSE                                       ;   gdepw(ji,jj,ik+1) = gdepw_1d(ik+1)
+                  IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
+                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
                   ENDIF
 !gm Bug?  check the gdepw_1d
                   !       ... on ik
-                  gdept(ji,jj,ik) = gdepw_1d(ik) + ( gdepw  (ji,jj,ik+1) - gdepw_1d(ik) )   &
-                     &                           * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
-                     &                           / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
-                  e3t  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw  (ji,jj,ik+1) - gdepw_1d(ik)  )  & 
-                     &                           / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
-                  e3w  (ji,jj,ik) = 0.5_wp * ( gdepw(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   &
+                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0   (ji,jj,ik+1) - gdepw_1d(ik) )   &
+                     &                           * ((gdept_1d(      ik  ) - gdepw_1d(ik) )   &
+                     &                           / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ))
+                  e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
+                     &                          / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
+                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   &
                      &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
                   !       ... on ik+1
-                  e3w  (ji,jj,ik+1) = e3t  (ji,jj,ik)
-                  e3t  (ji,jj,ik+1) = e3t  (ji,jj,ik)
-                  gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
+                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
+                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
+                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
                ENDIF
             ENDIF
@@ -941,14 +941,14 @@
             ik = mbathy(ji,jj)
             IF( ik > 0 ) THEN               ! ocean point only
-               e3tp (ji,jj) = e3t(ji,jj,ik  )
-               e3wp (ji,jj) = e3w(ji,jj,ik  )
+               e3tp (ji,jj) = e3t_0(ji,jj,ik)
+               e3wp (ji,jj) = e3w_0(ji,jj,ik)
                ! test
-               zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik  )
+               zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  )
                IF( zdiff <= 0._wp .AND. lwp ) THEN 
                   it = it + 1
                   WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
                   WRITE(numout,*) ' bathy = ', bathy(ji,jj)
-                  WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zdiff
-                  WRITE(numout,*) ' e3tp  = ', e3t  (ji,jj,ik), ' e3wp  = ', e3w  (ji,jj,ik  )
+                  WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
+                  WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  )
                ENDIF
             ENDIF
@@ -958,63 +958,63 @@
       ! Scale factors and depth at U-, V-, UW and VW-points
       DO jk = 1, jpk                        ! initialisation to z-scale factors
-         e3u (:,:,jk) = e3t_1d(jk)
-         e3v (:,:,jk) = e3t_1d(jk)
-         e3uw(:,:,jk) = e3w_1d(jk)
-         e3vw(:,:,jk) = e3w_1d(jk)
+         e3u_0 (:,:,jk) = e3t_1d(jk)
+         e3v_0 (:,:,jk) = e3t_1d(jk)
+         e3uw_0(:,:,jk) = e3w_1d(jk)
+         e3vw_0(:,:,jk) = e3w_1d(jk)
       END DO
       DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk) )
-               e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk) )
-               e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
-               e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( e3u , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw, 'U', 1._wp )   ! lateral boundary conditions
-      CALL lbc_lnk( e3v , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw, 'V', 1._wp )
+               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
+               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
+               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
+               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
+      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp )
       !
       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
-         WHERE( e3u (:,:,jk) == 0._wp )   e3u (:,:,jk) = e3t_1d(jk)
-         WHERE( e3v (:,:,jk) == 0._wp )   e3v (:,:,jk) = e3t_1d(jk)
-         WHERE( e3uw(:,:,jk) == 0._wp )   e3uw(:,:,jk) = e3w_1d(jk)
-         WHERE( e3vw(:,:,jk) == 0._wp )   e3vw(:,:,jk) = e3w_1d(jk)
+         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk)
+         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk)
+         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk)
+         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk)
       END DO
       
       ! Scale factor at F-point
       DO jk = 1, jpk                        ! initialisation to z-scale factors
-         e3f(:,:,jk) = e3t_1d(jk)
+         e3f_0(:,:,jk) = e3t_1d(jk)
       END DO
       DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( e3f, 'F', 1._wp )       ! Lateral boundary conditions
+               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions
       !
       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
-         WHERE( e3f(:,:,jk) == 0._wp )   e3f(:,:,jk) = e3t_1d(jk)
+         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk)
       END DO
 !!gm  bug ? :  must be a do loop with mj0,mj1
       ! 
-      e3t(:,mj0(1),:) = e3t(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
-      e3w(:,mj0(1),:) = e3w(:,mj0(2),:) 
-      e3u(:,mj0(1),:) = e3u(:,mj0(2),:) 
-      e3v(:,mj0(1),:) = e3v(:,mj0(2),:) 
-      e3f(:,mj0(1),:) = e3f(:,mj0(2),:) 
+      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
+      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 
+      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 
+      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 
+      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 
 
       ! Control of the sign
-      IF( MINVAL( e3t  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t   <= 0' )
-      IF( MINVAL( e3w  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w   <= 0' )
-      IF( MINVAL( gdept(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw <  0' )
-      IF( MINVAL( gdepw(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw <  0' )
+      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' )
+      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' )
+      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' )
+      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' )
      
-      ! Compute gdep3w (vertical sum of e3w)
-      gdep3w(:,:,1) = 0.5_wp * e3w(:,:,1)
+      ! Compute gdep3w_0 (vertical sum of e3w)
+      gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
       DO jk = 2, jpk
-         gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk) 
+         gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
       END DO
         
@@ -1025,10 +1025,10 @@
             DO ji = 1, jpi
                ik = MAX( mbathy(ji,jj), 1 )
-               zprt(ji,jj,1) = e3t   (ji,jj,ik)
-               zprt(ji,jj,2) = e3w   (ji,jj,ik)
-               zprt(ji,jj,3) = e3u   (ji,jj,ik)
-               zprt(ji,jj,4) = e3v   (ji,jj,ik)
-               zprt(ji,jj,5) = e3f   (ji,jj,ik)
-               zprt(ji,jj,6) = gdep3w(ji,jj,ik)
+               zprt(ji,jj,1) = e3t_0   (ji,jj,ik)
+               zprt(ji,jj,2) = e3w_0   (ji,jj,ik)
+               zprt(ji,jj,3) = e3u_0   (ji,jj,ik)
+               zprt(ji,jj,4) = e3v_0   (ji,jj,ik)
+               zprt(ji,jj,5) = e3f_0   (ji,jj,ik)
+               zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)
             END DO
          END DO
@@ -1356,34 +1356,34 @@
       ENDIF 
 
-      CALL lbc_lnk( e3t , 'T', 1._wp )
-      CALL lbc_lnk( e3u , 'U', 1._wp )
-      CALL lbc_lnk( e3v , 'V', 1._wp )
-      CALL lbc_lnk( e3f , 'F', 1._wp )
-      CALL lbc_lnk( e3w , 'W', 1._wp )
-      CALL lbc_lnk( e3uw, 'U', 1._wp )
-      CALL lbc_lnk( e3vw, 'V', 1._wp )
-
-      fsdepw(:,:,:) = gdepw (:,:,:)
-      fsde3w(:,:,:) = gdep3w(:,:,:)
-      !
-      where (e3t   (:,:,:).eq.0.0)  e3t(:,:,:) = 1.0
-      where (e3u   (:,:,:).eq.0.0)  e3u(:,:,:) = 1.0
-      where (e3v   (:,:,:).eq.0.0)  e3v(:,:,:) = 1.0
-      where (e3f   (:,:,:).eq.0.0)  e3f(:,:,:) = 1.0
-      where (e3w   (:,:,:).eq.0.0)  e3w(:,:,:) = 1.0
-      where (e3uw  (:,:,:).eq.0.0)  e3uw(:,:,:) = 1.0
-      where (e3vw  (:,:,:).eq.0.0)  e3vw(:,:,:) = 1.0
-
-
-      fsdept(:,:,:) = gdept (:,:,:)
-      fsdepw(:,:,:) = gdepw (:,:,:)
-      fsde3w(:,:,:) = gdep3w(:,:,:)
-      fse3t (:,:,:) = e3t   (:,:,:)
-      fse3u (:,:,:) = e3u   (:,:,:)
-      fse3v (:,:,:) = e3v   (:,:,:)
-      fse3f (:,:,:) = e3f   (:,:,:)
-      fse3w (:,:,:) = e3w   (:,:,:)
-      fse3uw(:,:,:) = e3uw  (:,:,:)
-      fse3vw(:,:,:) = e3vw  (:,:,:)
+      CALL lbc_lnk( e3t_0 , 'T', 1._wp )
+      CALL lbc_lnk( e3u_0 , 'U', 1._wp )
+      CALL lbc_lnk( e3v_0 , 'V', 1._wp )
+      CALL lbc_lnk( e3f_0 , 'F', 1._wp )
+      CALL lbc_lnk( e3w_0 , 'W', 1._wp )
+      CALL lbc_lnk( e3uw_0, 'U', 1._wp )
+      CALL lbc_lnk( e3vw_0, 'V', 1._wp )
+
+      fsdepw(:,:,:) = gdepw_0 (:,:,:)
+      fsde3w(:,:,:) = gdep3w_0(:,:,:)
+      !
+      where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
+      where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
+      where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
+      where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
+      where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
+      where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
+      where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
+
+
+      fsdept(:,:,:) = gdept_0 (:,:,:)
+      fsdepw(:,:,:) = gdepw_0 (:,:,:)
+      fsde3w(:,:,:) = gdep3w_0(:,:,:)
+      fse3t (:,:,:) = e3t_0   (:,:,:)
+      fse3u (:,:,:) = e3u_0   (:,:,:)
+      fse3v (:,:,:) = e3v_0   (:,:,:)
+      fse3f (:,:,:) = e3f_0   (:,:,:)
+      fse3w (:,:,:) = e3w_0   (:,:,:)
+      fse3uw(:,:,:) = e3uw_0  (:,:,:)
+      fse3vw(:,:,:) = e3vw_0  (:,:,:)
 !!
       ! HYBRID : 
@@ -1400,18 +1400,18 @@
 
       IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
-         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)   ), ' MAX ', MAXVAL( mbathy(:,:) )
-         WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ),   &
-            &                          ' w ', MINVAL( fsdepw(:,:,:) ), '3w '  , MINVAL( fsde3w(:,:,:) )
-         WRITE(numout,*) ' MIN val e3    t ', MINVAL( fse3t (:,:,:) ), ' f '  , MINVAL( fse3f (:,:,:) ),   &
-            &                          ' u ', MINVAL( fse3u (:,:,:) ), ' u '  , MINVAL( fse3v (:,:,:) ),   &
-            &                          ' uw', MINVAL( fse3uw(:,:,:) ), ' vw'  , MINVAL( fse3vw(:,:,:) ),   &
-            &                          ' w ', MINVAL( fse3w (:,:,:) )
-
-         WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ),   &
-            &                          ' w ', MAXVAL( fsdepw(:,:,:) ), '3w '  , MAXVAL( fsde3w(:,:,:) )
-         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( fse3t (:,:,:) ), ' f '  , MAXVAL( fse3f (:,:,:) ),   &
-            &                          ' u ', MAXVAL( fse3u (:,:,:) ), ' u '  , MAXVAL( fse3v (:,:,:) ),   &
-            &                          ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw'  , MAXVAL( fse3vw(:,:,:) ),   &
-            &                          ' w ', MAXVAL( fse3w (:,:,:) )
+         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
+         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
+            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gdep3w_0(:,:,:) )
+         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0   (:,:,:) ),   &
+            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0   (:,:,:) ),   &
+            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0  (:,:,:) ),   &
+            &                          ' w ', MINVAL( e3w_0  (:,:,:) )
+
+         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
+            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gdep3w_0(:,:,:) )
+         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0   (:,:,:) ),   &
+            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0   (:,:,:) ),   &
+            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0  (:,:,:) ),   &
+            &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
       ENDIF
       !  END DO
@@ -1420,7 +1420,7 @@
          WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
          WRITE(numout,*) ' ~~~~~~  --------------------'
-         WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
-         WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk),     &
-            &                                 fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk )
+         WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
+         WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
+            &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
          DO jj = mj0(20), mj1(20)
             DO ji = mi0(20), mi1(20)
@@ -1428,7 +1428,7 @@
                WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
                WRITE(numout,*) ' ~~~~~~  --------------------'
-               WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
-               WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk),     &
-                  &                                 fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
+               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
+               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
+                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
             END DO
          END DO
@@ -1438,7 +1438,7 @@
                WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
                WRITE(numout,*) ' ~~~~~~  --------------------'
-               WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
-               WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk),     &
-                  &                                 fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
+               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
+               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
+                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
             END DO
          END DO
@@ -1558,7 +1558,7 @@
                zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
                zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
-               gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
-               gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
-               gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
+               gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
+               gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
+               gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
             END DO
            !
@@ -1581,12 +1581,12 @@
                   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
                !
-               e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
-               e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
-               e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
-               e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
                !
-               e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
-               e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
-               e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
             END DO
         END DO
@@ -1686,7 +1686,7 @@
 
           DO jk = 1, jpk
-             gdept (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
-             gdepw (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
-             gdep3w(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
+             gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
+             gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
+             gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
           END DO
 
@@ -1710,12 +1710,12 @@
                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) )
 
-             e3t(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
-             e3u(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
-             e3v(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
-             e3f(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
+             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
+             e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
+             e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
+             e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
              !
-             e3w(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
-             e3uw(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
-             e3vw(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
+             e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
+             e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
+             e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
           END DO
 
@@ -1723,8 +1723,8 @@
       ENDDO
       !
-      CALL lbc_lnk(e3t ,'T',1.) ; CALL lbc_lnk(e3u ,'T',1.)
-      CALL lbc_lnk(e3v ,'T',1.) ; CALL lbc_lnk(e3f ,'T',1.)
-      CALL lbc_lnk(e3w ,'T',1.)
-      CALL lbc_lnk(e3uw,'T',1.) ; CALL lbc_lnk(e3vw,'T',1.)
+      CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.)
+      CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.)
+      CALL lbc_lnk(e3w_0 ,'T',1.)
+      CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
       !
       !                                               ! =============
@@ -1784,7 +1784,7 @@
          zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
          zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
-         gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
-         gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
-         gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
+         gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
+         gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
+         gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
       END DO
 !!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
@@ -1792,12 +1792,12 @@
          DO ji = 1, jpi
             DO jk = 1, jpk
-              e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
-              e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
-              e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
-              e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
+              e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
+              e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
+              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
+              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
               !
-              e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
-              e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
-              e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
+              e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
+              e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
+              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
             END DO
          END DO
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr_substitute.h90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr_substitute.h90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr_substitute.h90	(revision 3865)
@@ -8,86 +8,85 @@
    !!            3.1  !  2009-02  (G. Madec, M. Leclair)  pure z* coordinate
    !!----------------------------------------------------------------------
-! reference for s- or zps-coordinate (3D no time dependency)
-#   define  fsdept_0(i,j,k)  gdept(i,j,k)
-#   define  fsdepw_0(i,j,k)  gdepw(i,j,k)
-#   define  fsde3w_0(i,j,k)  gdep3w(i,j,k)
-#   define  fse3t_0(i,j,k)   e3t(i,j,k)
-#   define  fse3u_0(i,j,k)   e3u(i,j,k)
-#   define  fse3v_0(i,j,k)   e3v(i,j,k)
-#   define  fse3f_0(i,j,k)   e3f(i,j,k)
-#   define  fse3w_0(i,j,k)   e3w(i,j,k)
-#   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 ...
+! ... In that case all occurence of the above statement functions
+!     have to be replaced in the code by xxx_n
+#   define  fsdept(i,j,k)    gdept_n(i,j,k)
+#   define  fsdepw(i,j,k)    gdepw_n(i,j,k)
+#   define  fsde3w(i,j,k)    gdep3w_n(i,j,k)
+#   define  fse3t(i,j,k)     e3t_n(i,j,k) 
+#   define  fse3u(i,j,k)     e3u_n(i,j,k) 
+#   define  fse3v(i,j,k)     e3v_n(i,j,k) 
+#   define  fse3f(i,j,k)     e3f_n(i,j,k) 
+#   define  fse3w(i,j,k)     e3w_n(i,j,k) 
+#   define  fse3uw(i,j,k)    e3uw_n(i,j,k)
+#   define  fse3vw(i,j,k)    e3vw_n(i,j,k)
 
 #else
 ! z- or s-coordinate (1D or 3D + no time dependency) use reference in all cases
-#   define  fsdept(i,j,k)  fsdept_0(i,j,k)
-#   define  fsdepw(i,j,k)  fsdepw_0(i,j,k)
-#   define  fsde3w(i,j,k)  fsde3w_0(i,j,k)
-#   define  fse3t(i,j,k)   fse3t_0(i,j,k)
-#   define  fse3u(i,j,k)   fse3u_0(i,j,k)
-#   define  fse3v(i,j,k)   fse3v_0(i,j,k)
-#   define  fse3f(i,j,k)   fse3f_0(i,j,k)
-#   define  fse3w(i,j,k)   fse3w_0(i,j,k)
-#   define  fse3uw(i,j,k)  fse3uw_0(i,j,k)
-#   define  fse3vw(i,j,k)  fse3vw_0(i,j,k)
 
-#   define  fse3t_b(i,j,k)   fse3t_0(i,j,k)
-#   define  fse3u_b(i,j,k)   fse3u_0(i,j,k)
-#   define  fse3v_b(i,j,k)   fse3v_0(i,j,k)
-#   define  fse3uw_b(i,j,k)  fse3uw_0(i,j,k)
-#   define  fse3vw_b(i,j,k)  fse3vw_0(i,j,k)
+#   define  fse3t_b(i,j,k)   e3t_0(i,j,k)
+#   define  fse3u_b(i,j,k)   e3u_0(i,j,k)
+#   define  fse3v_b(i,j,k)   e3v_0(i,j,k)
+#   define  fse3uw_b(i,j,k)  e3uw_0(i,j,k)
+#   define  fse3vw_b(i,j,k)  e3vw_0(i,j,k)
 
-#   define  fsdept_n(i,j,k)  fsdept_0(i,j,k)
-#   define  fsdepw_n(i,j,k)  fsdepw_0(i,j,k)
-#   define  fsde3w_n(i,j,k)  fsde3w_0(i,j,k)
-#   define  fse3t_n(i,j,k)   fse3t_0(i,j,k)
-#   define  fse3u_n(i,j,k)   fse3u_0(i,j,k)
-#   define  fse3v_n(i,j,k)   fse3v_0(i,j,k)
-#   define  fse3f_n(i,j,k)   fse3f_0(i,j,k)
-#   define  fse3w_n(i,j,k)   fse3w_0(i,j,k)
-#   define  fse3uw_n(i,j,k)  fse3uw_0(i,j,k)
-#   define  fse3vw_n(i,j,k)  fse3vw_0(i,j,k)
+#   define  fsdept_n(i,j,k)  gdept_0(i,j,k)
+#   define  fsdepw_n(i,j,k)  gdepw_0(i,j,k)
+#   define  fsde3w_n(i,j,k)  gdep3w_0(i,j,k)
+#   define  fse3t_n(i,j,k)   e3t_0(i,j,k)
+#   define  fse3u_n(i,j,k)   e3u_0(i,j,k)
+#   define  fse3v_n(i,j,k)   e3v_0(i,j,k)
+#   define  fse3f_n(i,j,k)   e3f_0(i,j,k)
+#   define  fse3w_n(i,j,k)   e3w_0(i,j,k)
+#   define  fse3uw_n(i,j,k)  e3uw_0(i,j,k)
+#   define  fse3vw_n(i,j,k)  e3vw_0(i,j,k)
 
-#   define  fse3t_m(i,j,k)   fse3t_0(i,j,k)
+#   define  fse3t_a(i,j,k)   e3t_0(i,j,k)
+#   define  fse3u_a(i,j,k)   e3u_0(i,j,k)
+#   define  fse3v_a(i,j,k)   e3v_0(i,j,k)
 
-#   define  fse3t_a(i,j,k)   fse3t_0(i,j,k)
-#   define  fse3u_a(i,j,k)   fse3u_0(i,j,k)
-#   define  fse3v_a(i,j,k)   fse3v_0(i,j,k)
+#   define  fse3t_m(i,j)     e3t_0(i,j,1)
+
+! This part should be removed one day ...
+! ... In that case all occurence of the above statement functions
+!     have to be replaced in the code by xxx_n
+#   define  fsdept(i,j,k)    gdept_0(i,j,k)
+#   define  fsdepw(i,j,k)    gdepw_0(i,j,k)
+#   define  fsde3w(i,j,k)    gdep3w_0(i,j,k)
+#   define  fse3t(i,j,k)     e3t_0(i,j,k)
+#   define  fse3u(i,j,k)     e3u_0(i,j,k)
+#   define  fse3v(i,j,k)     e3v_0(i,j,k)
+#   define  fse3f(i,j,k)     e3f_0(i,j,k)
+#   define  fse3w(i,j,k)     e3w_0(i,j,k)
+#   define  fse3uw(i,j,k)    e3uw_0(i,j,k)
+#   define  fse3vw(i,j,k)    e3vw_0(i,j,k)
+
 #endif
    !!----------------------------------------------------------------------
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90	(revision 3865)
@@ -222,5 +222,5 @@
             DO ji = 1, jpi
                DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
-                  zl = fsdept_0(ji,jj,jk)
+                  zl = gdept_0(ji,jj,jk)
                   IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data
                      ztp(jk) =  ptsd(ji,jj,1    ,jp_tem)
@@ -260,5 +260,5 @@
                   ik = mbkt(ji,jj) 
                   IF( ik > 1 ) THEN
-                     zl = ( gdept_1d(ik) - fsdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
+                     zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
                      ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem)
                      ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal)
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 3865)
@@ -90,6 +90,4 @@
          neuler = 1                              ! Set time-step indicator at nit000 (leap-frog)
          CALL rst_read                           ! Read the restart file
-         !                                       ! define e3u_b, e3v_b from e3t_b read in restart file
-         CALL dom_vvl_2( nit000, fse3u_b(:,:,:), fse3v_b(:,:,:) )
          CALL day_init                           ! model calendar (using both namelist and restart infos)
       ELSE
@@ -131,6 +129,4 @@
             ENDDO
          ENDIF
-         !                                       ! define e3u_b, e3v_b from e3t_b initialized in domzgr
-         CALL dom_vvl_2( nit000, fse3u_b(:,:,:), fse3v_b(:,:,:) )
          ! 
       ENDIF
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 3865)
@@ -760,6 +760,6 @@
       DO jj = 2, jpjm1
         DO ji = 2, jpim1
-          zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshu_n(ji,jj) * znad)
-          zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshv_n(ji,jj) * znad)
+          zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshn(ji,jj) * znad)    ! probable bug: changed from sshu_n for ztilde compilation
+          zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshn(ji,jj) * znad)    ! probable bug: changed from sshv_n for ztilde compilation
         END DO
       END DO
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 3865)
@@ -201,6 +201,6 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi    
-                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
-                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
+                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0_wp * un(ji,jj,jk) + ua(ji,jj,jk) )
+                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0_wp * vn(ji,jj,jk) + va(ji,jj,jk) )
                      !
                      ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
@@ -214,23 +214,24 @@
          ELSE                             ! Variable volume !
             !                             ! ================!
+            ! Before scale factor at t-points
+            ! (used as a now filtered scale factor until the swap)
+            ! ----------------------------------------------------
+            fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) )
+            ! Add volume filter correction: compatibility with tracer advection scheme
+            ! => time filter + conservation correction (only at the first level)
+            fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
             !
-            DO jk = 1, jpkm1                 ! Before scale factor at t-points
-               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
-                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
-                  &                         - 2._wp * fse3t_n(:,:,jk)            )
-            END DO
-            zec = atfp * rdt / rau0          ! Add filter correction only at the 1st level of t-point scale factors
-            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
-            !
-            IF( ln_dynadv_vec ) THEN         ! vector invariant form (no thickness weighted calulation)
-               !
-               !                                      ! before scale factors at u- & v-pts (computed from fse3t_b)
-               CALL dom_vvl_2( kt, fse3u_b(:,:,:), fse3v_b(:,:,:) )
-               !
-               DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap: applied on velocity
-                  DO jj = 1, jpj                      !                                                 --------
+            IF( ln_dynadv_vec ) THEN
+               ! Before scale factor at (u/v)-points
+               ! -----------------------------------
+               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
+               ! -----------------------------------
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpj
                      DO ji = 1, jpi
-                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
-                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
+                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
+                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
                         !
                         ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
@@ -242,10 +243,13 @@
                END DO
                !
-            ELSE                             ! flux form (thickness weighted calulation)
-               !
-               CALL dom_vvl_2( kt, ze3u_f, ze3v_f )   ! before scale factors at u- & v-pts (computed from fse3t_b)
-               !
-               DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap: 
-                  DO jj = 1, jpj                      !                   applied on thickness weighted velocity
+            ELSE
+               ! 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
+               ! -----------------------------------             ===========================
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpj
                      DO ji = 1, jpi                   !                              ---------------------------
                         zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 3865)
@@ -126,4 +126,5 @@
       REAL(wp), POINTER, DIMENSION(:,:) :: zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e 
       REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum
+      REAL(wp), POINTER, DIMENSION(:,:) :: zhu_b, zhv_b
       !!----------------------------------------------------------------------
       !
@@ -133,4 +134,5 @@
       CALL wrk_alloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e   )
       CALL wrk_alloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum )
+      CALL wrk_alloc( jpi, jpj, zhu_b, zhv_b                                     )
       !
       IF( kt == nit000 ) THEN             !* initialisation
@@ -199,9 +201,9 @@
 #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
@@ -215,4 +217,19 @@
          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)
@@ -355,4 +372,6 @@
          vb_b(:,:) = vb_b(:,:) * hvr(:,:)
       ENDIF
+      ub_b(:,:) = ub_b(:,:) * zhu_b(:,:)
+      vb_b(:,:) = vb_b(:,:) * zhv_b(:,:)
 
       ! -----------------------------------------------------------------------
@@ -683,4 +702,5 @@
       CALL wrk_dealloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e   )
       CALL wrk_dealloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum )
+      CALL wrk_dealloc( jpi, jpj, zhu_b, zhv_b                                     )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
@@ -698,4 +718,5 @@
       CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
       !
+      REAL(wp), POINTER, DIMENSION(:,:) :: zzhu_b, zzhv_b
       INTEGER ::  ji, jk        ! dummy loop indices
       !!----------------------------------------------------------------------
@@ -706,4 +727,5 @@
             CALL iom_get( numror, jpdom_autoglo, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
          ELSE
+            CALL wrk_alloc( jpi, jpj, zzhu_b, zzhv_b )
             un_b (:,:) = 0._wp
             vn_b (:,:) = 0._wp
@@ -712,12 +734,12 @@
                DO jk = 1, jpkm1
                   DO ji = 1, jpij
-                     un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
-                     vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
+                     un_b(ji,1) = un_b(ji,1) + fse3u_n(ji,1,jk) * un(ji,1,jk)
+                     vn_b(ji,1) = vn_b(ji,1) + fse3v_n(ji,1,jk) * vn(ji,1,jk)
                   END DO
                END DO
             ELSE                             ! No  vector opt.
                DO jk = 1, jpkm1
-                  un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
-                  vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
+                  un_b(:,:) = un_b(:,:) + fse3u_n(:,:,jk) * un(:,:,jk)
+                  vn_b(:,:) = vn_b(:,:) + fse3v_n(:,:,jk) * vn(:,:,jk)
                END DO
             ENDIF
@@ -747,7 +769,47 @@
 
             IF( lk_vvl ) THEN
-               ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) )
-               vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) )
-            ELSE
+               CALL wrk_alloc( jpi, jpj, zzhu_b, zzhv_b )
+               ub_b (:,:) = 0.
+               vb_b (:,:) = 0.
+               zzhu_b(:,:) = 0.
+               zzhv_b(:,:) = 0.
+               ! vertical sum
+               IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
+                  DO jk = 1, jpkm1
+                     DO ji = 1, jpij
+                        ub_b  (ji,1) = ub_b (ji,1) + fse3u_b(ji,1,jk) * ub   (ji,1,jk)
+                        vb_b  (ji,1) = vb_b (ji,1) + fse3v_b(ji,1,jk) * vb   (ji,1,jk)
+                        zzhu_b(ji,1) = zhu_b(ji,1) + fse3u_b(ji,1,jk) * umask(ji,1,jk)
+                        zzhv_b(ji,1) = zhv_b(ji,1) + fse3v_b(ji,1,jk) * vmask(ji,1,jk)
+                     END DO
+                  END DO
+               ELSE                             ! No  vector opt.
+                  DO jk = 1, jpkm1
+                     ub_b  (:,:) = ub_b  (:,:) + fse3u_b(:,:,jk) * ub   (:,:,jk)
+                     vb_b  (:,:) = vb_b  (:,:) + fse3v_b(:,:,jk) * vb   (:,:,jk)
+                     zzhu_b(:,:) = zzhu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
+                     zzhv_b(:,:) = zzhv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
+                  END DO
+               ENDIF
+               ub_b(:,:) = ub_b(:,:) / ( zzhu_b(:,:) + 1. - umask(:,:,1) )
+               vb_b(:,:) = vb_b(:,:) / ( zzhv_b(:,:) + 1. - vmask(:,:,1) )
+               CALL wrk_dealloc( jpi, jpj, zzhu_b, zzhv_b )
+            ELSE             
+               ub_b (:,:) = 0.e0
+               vb_b (:,:) = 0.e0
+               ! vertical sum
+               IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
+                  DO jk = 1, jpkm1
+                     DO ji = 1, jpij
+                        ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk)
+                        vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk)
+                     END DO
+                  END DO
+               ELSE                             ! No  vector opt.
+                  DO jk = 1, jpkm1
+                     ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk)
+                     vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk)
+                  END DO
+               ENDIF
                ub_b(:,:) = ub_b(:,:) * hur(:,:)
                vb_b(:,:) = vb_b(:,:) * hvr(:,:)
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 3865)
@@ -8,9 +8,11 @@
    !!             -   !  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
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   ssh_wzv        : after ssh & now vertical velocity
-   !!   ssh_nxt        : filter ans swap the ssh arrays
+   !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   ssh_nxt        : after ssh
+   !!   ssh_swp        : filter ans swap the ssh arrays
+   !!   wzv            : compute now vertical velocity
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers variables
@@ -20,4 +22,5 @@
    USE divcur          ! hor. divergence and curl      (div & cur routines)
    USE iom             ! I/O library
+   USE restart         ! only for lrst_oce
    USE in_out_manager  ! I/O manager
    USE prtctl          ! Print control
@@ -44,6 +47,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
@@ -57,97 +61,36 @@
 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.
       !!----------------------------------------------------------------------
-      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
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d, zhdiv
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('ssh_wzv')
-      !
-      CALL wrk_alloc( jpi, jpj, z2d, zhdiv ) 
+      !
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv
+      INTEGER, INTENT(in) ::   kt                      ! time step
+      ! 
+      INTEGER             ::   jk                      ! dummy loop indice
+      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
+      !
+      CALL wrk_alloc( jpi, jpj, zhdiv ) 
       !
       IF( kt == nit000 ) THEN
          !
          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
       !
@@ -162,5 +105,5 @@
       zhdiv(:,:) = 0._wp
       DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
-        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
+        zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
       END DO
       !                                                ! Sea surface elevation time stepping
@@ -181,5 +124,5 @@
 #if defined key_bdy
       ssha(:,:) = ssha(:,:) * bdytmask(:,:)
-      CALL lbc_lnk( ssha, 'T', 1. )                    ! absolutly compulsory !! (jmm)
+      CALL lbc_lnk( ssha, 'T', 1. ) 
 #endif
 #if defined key_asminc
@@ -190,43 +133,102 @@
       ENDIF
 #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) )
+
+      !                                           !------------------------------!
+      !                                           !           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 )
+      !
+      CALL wrk_dealloc( jpi, jpj, zhdiv ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
+      !
+   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.
+      !!----------------------------------------------------------------------
+      !
+      INTEGER, INTENT(in) ::   kt           ! time step
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
+      !
+      INTEGER             ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp)            ::   z1_2dt       ! local scalars
+      !!----------------------------------------------------------------------
+      
+      IF( nn_timing == 1 )  CALL timing_start('wzv')
+      !
+      CALL wrk_alloc( jpi, jpj, z2d ) 
+      CALL wrk_alloc( jpi, jpj, jpk, z3d, zhdiv ) 
+      !
+      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. / ( 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) = r1_e12t(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( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions
-      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
+         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.
-         CALL wrk_alloc( jpi,jpj,jpk, z3d )
-         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
+         z2d(:,:) = rau0 * e12t(:,:)
          DO jk = 1, jpk
             z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
@@ -234,17 +236,16 @@
          CALL iom_put( "w_masstr" , z3d                     )  
          CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
-         CALL wrk_dealloc( jpi,jpj,jpk, z3d )
-      ENDIF
-      !
-      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
-      !
-      CALL wrk_dealloc( jpi, jpj, z2d, zhdiv ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('ssh_wzv')
-      !
-   END SUBROUTINE ssh_wzv
-
-
-   SUBROUTINE ssh_nxt( kt )
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi, jpj, z2d  ) 
+      CALL wrk_dealloc( jpi, jpj, jpk, z3d, zhdiv ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('wzv')
+
+
+   END SUBROUTINE wzv
+
+
+   SUBROUTINE ssh_swp( kt )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE ssh_nxt  ***
@@ -252,5 +253,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,86 +267,20 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      !!
-      INTEGER  ::   ji, jj   ! dummy loop indices
-      REAL(wp) ::   zec      ! temporary scalar
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('ssh_swp')
       !
       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
-            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
-            !
-         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
-         !
+      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
+         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered
+         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
+         sshn(:,:) = ssha(:,:)                           ! now <-- after
       ENDIF
       !
@@ -357,7 +292,7 @@
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
       !
-      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
-      !
-   END SUBROUTINE ssh_nxt
+      IF( nn_timing == 1 )  CALL timing_stop('ssh_swp')
+      !
+   END SUBROUTINE ssh_swp
 
    !!======================================================================
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 3865)
@@ -91,4 +91,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
@@ -105,5 +106,5 @@
       !!                  ***  FUNCTION sbc_oce_alloc  ***
       !!---------------------------------------------------------------------
-      INTEGER :: ierr(4)
+      INTEGER :: ierr(5)
       !!---------------------------------------------------------------------
       ierr(:) = 0
@@ -126,4 +127,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/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 3865)
@@ -412,5 +412,5 @@
 ! Freezing/melting potential
 ! Calculated over NEMO leapfrog timestep (hence 2*dt)
-      nfrzmlt(:,:)=rau0*rcp*fse3t_m(:,:,1)*(Tocnfrz-sst_m(:,:))/(2.0*dt)
+      nfrzmlt(:,:)=rau0*rcp*fse3t_m(:,:)*(Tocnfrz-sst_m(:,:))/(2.0*dt)
 
       ztmp(:,:) = nfrzmlt(:,:)
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 3865)
@@ -24,10 +24,9 @@
    PRIVATE
 
-   PUBLIC   sbc_ssm         ! routine called by step.F90
-   PUBLIC   sbc_ssm_init    ! routine called by sbcmod.F90
-   
+   PUBLIC   sbc_ssm    ! routine called by step.F90
+
    LOGICAL, SAVE  ::   l_ssm_mean = .FALSE.       ! keep track of whether means have been read
                                                   ! from restart file
-
+   
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -67,5 +66,6 @@
          ELSE                    ;   ssh_m(:,:) = sshn(:,:)
          ENDIF
-
+         !
+         IF( lk_vvl )   fse3t_m(:,:) = fse3t_n(:,:,1)
          !
       ELSE
@@ -84,4 +84,5 @@
             ELSE                    ;   ssh_m(:,:) = zcoef *   sshn(:,:)
             ENDIF
+            IF( lk_vvl )   fse3t_m(:,:) = zcoef * fse3t_n(:,:,1)
             !                                             ! ---------------------------------------- !
          ELSEIF( MOD( kt - 2 , nn_fsbc ) == 0 ) THEN      !   Initialisation: New mean computation   !
@@ -92,4 +93,5 @@
             sss_m(:,:) = 0.e0
             ssh_m(:,:) = 0.e0
+            IF( lk_vvl )   fse3t_m(:,:) = 0.e0
          ENDIF
          !                                                ! ---------------------------------------- !
@@ -104,4 +106,5 @@
          ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:)
          ENDIF
+         IF( lk_vvl )   fse3t_m(:,:) = fse3t_m(:,:) + fse3t_n(:,:,1)
 
          !                                                ! ---------------------------------------- !
@@ -114,4 +117,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
@@ -130,4 +134,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
@@ -168,4 +175,5 @@
             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 ) CALL iom_get( numror, jpdom_autoglo, 'fse3t_m', fse3t_m(:,:) )
             !
             IF( zf_sbc /= REAL( nn_fsbc, wp ) ) THEN      ! nn_fsbc has changed between 2 runs
@@ -178,4 +186,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'
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 3865)
@@ -155,10 +155,10 @@
              DO jj = 1, jpj-1
                 DO ji = 1, jpi-1
-                   usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept(ji,jj,jk) , gdept(ji+1,jj  ,jk))))
-                   vsd3d(ji,jj,jk) = vsd2d(ji,jj)*exp(2.0*vwavenum(ji,jj)*(-MIN( gdept(ji,jj,jk) , gdept(ji  ,jj+1,jk))))
+                   usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk))))
+                   vsd3d(ji,jj,jk) = vsd2d(ji,jj)*exp(2.0*vwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk))))
                 END DO
              END DO
-             usd3d(jpi,:,jk) = usd2d(jpi,:)*exp( 2.0*uwavenum(jpi,:)*(-gdept(jpi,:,jk)) )
-             vsd3d(:,jpj,jk) = vsd2d(:,jpj)*exp( 2.0*vwavenum(:,jpj)*(-gdept(:,jpj,jk)) )
+             usd3d(jpi,:,jk) = usd2d(jpi,:)*exp( 2.0*uwavenum(jpi,:)*(-gdept_0(jpi,:,jk)) )
+             vsd3d(:,jpj,jk) = vsd2d(:,jpj)*exp( 2.0*vwavenum(:,jpj)*(-gdept_0(:,jpj,jk)) )
           END DO
 
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 3865)
@@ -74,5 +74,5 @@
 CONTAINS
 
-   SUBROUTINE eos_insitu( pts, prd )
+   SUBROUTINE eos_insitu( pts, prd, pdep )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE eos_insitu  ***
@@ -114,4 +114,5 @@
       !                                                      ! 2 : salinity               [psu]
       REAL(wp), DIMENSION(:,:,:)  , INTENT(  out) ::   prd   ! in situ density            [-]
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pdep  ! depth                      [m]
       !!
       INTEGER  ::   ji, jj, jk           ! dummy loop indices
@@ -140,5 +141,5 @@
                   zt = pts   (ji,jj,jk,jp_tem)
                   zs = pts   (ji,jj,jk,jp_sal)
-                  zh = fsdept(ji,jj,jk)        ! depth
+                  zh = pdep(ji,jj,jk)        ! depth
                   zsr= zws   (ji,jj,jk)        ! square root salinity
                   !
@@ -198,5 +199,5 @@
 
 
-   SUBROUTINE eos_insitu_pot( pts, prd, prhop )
+   SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE eos_insitu_pot  ***
@@ -249,4 +250,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-]
       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m]
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -271,5 +273,5 @@
                   zt = pts   (ji,jj,jk,jp_tem)
                   zs = pts   (ji,jj,jk,jp_sal)
-                  zh = fsdept(ji,jj,jk)        ! depth
+                  zh = pdep(ji,jj,jk)        ! depth
                   zsr= zws   (ji,jj,jk)        ! square root salinity
                   !
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 3865)
@@ -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)
@@ -92,4 +93,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/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 3865)
@@ -66,5 +66,4 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ahu_bbl_0, ahv_bbl_0   ! diffusive bbl flux coefficients at u and v-points
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   e3u_bbl_0, e3v_bbl_0   ! thichness of the bbl (e3) at u and v-points (PUBLIC for TAM)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   r1_e1e2t               ! inverse of the cell surface at t-point      [1/m2] (PUBLIC for TAM)
 
    !! * Substitutions
@@ -85,5 +84,5 @@
          &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d  (jpi,jpj) , mgrhv(jpi,jpj) ,     &
          &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                          &
-         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) , r1_e1e2t(jpi,jpj)                  , STAT= tra_bbl_alloc )
+         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) , STAT= tra_bbl_alloc                )
          !
       IF( lk_mpp            )   CALL mpp_sum ( tra_bbl_alloc )
@@ -217,5 +216,5 @@
 #  endif
                ik = mbkt(ji,jj)                            ! bottom T-level index
-               zbtr = r1_e1e2t(ji,jj)  / fse3t(ji,jj,ik)
+               zbtr = r1_e12t(ji,jj)  / fse3t(ji,jj,ik)
                pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         &
                   &               + (   ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   &
@@ -279,15 +278,15 @@
                   !
                   !                                               ! up  -slope T-point (shelf bottom point)
-                  zbtr = r1_e1e2t(iis,jj) / fse3t(iis,jj,ikus)
+                  zbtr = r1_e12t(iis,jj) / fse3t(iis,jj,ikus)
                   ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr
                   pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra
                   !
                   DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
-                     zbtr = r1_e1e2t(iid,jj) / fse3t(iid,jj,jk)
+                     zbtr = r1_e12t(iid,jj) / fse3t(iid,jj,jk)
                      ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr
                      pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra
                   END DO
                   !
-                  zbtr = r1_e1e2t(iid,jj) / fse3t(iid,jj,ikud)
+                  zbtr = r1_e12t(iid,jj) / fse3t(iid,jj,ikud)
                   ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr
                   pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra
@@ -301,15 +300,15 @@
                   !
                   ! up  -slope T-point (shelf bottom point)
-                  zbtr = r1_e1e2t(ji,ijs) / fse3t(ji,ijs,ikvs)
+                  zbtr = r1_e12t(ji,ijs) / fse3t(ji,ijs,ikvs)
                   ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
                   pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra
                   !
                   DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
-                     zbtr = r1_e1e2t(ji,ijd) / fse3t(ji,ijd,jk)
+                     zbtr = r1_e12t(ji,ijd) / fse3t(ji,ijd,jk)
                      ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
                      pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra
                   END DO
                   !                                               ! down-slope T-point (deep bottom point)
-                  zbtr = r1_e1e2t(ji,ijd) / fse3t(ji,ijd,ikvd)
+                  zbtr = r1_e12t(ji,ijd) / fse3t(ji,ijd,ikvd)
                   ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
                   pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra
@@ -423,5 +422,5 @@
             ztb (ji,jj) = tsb(ji,jj,ik,jp_tem) * tmask(ji,jj,1)      ! bottom before T and S
             zsb (ji,jj) = tsb(ji,jj,ik,jp_sal) * tmask(ji,jj,1)
-            zdep(ji,jj) = fsdept_0(ji,jj,ik)        ! bottom T-level reference depth
+            zdep(ji,jj) = gdept_0(ji,jj,ik)         ! bottom T-level reference depth
             !
             zub(ji,jj) = un(ji,jj,mbku(ji,jj))      ! bottom velocity
@@ -594,8 +593,4 @@
       IF( nn_eos /= 0 )   CALL ctl_stop ( ' bbl parameterisation requires eos = 0. We stop.' )
 
-
-      !                             !* inverse of surface of T-cells
-      r1_e1e2t(:,:) = 1._wp / ( e1t(:,:) * e2t(:,:) )
-
       !                             !* vertical index of  "deep" bottom u- and v-points
       DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
@@ -605,21 +600,21 @@
          END DO
       END DO
-      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
+      ! convert into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
       zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
       zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
 
-                                        !* sign of grad(H) at u- and v-points
+                                     !* sign of grad(H) at u- and v-points
       mgrhu(jpi,:) = 0.    ;    mgrhu(:,jpj) = 0.   ;    mgrhv(jpi,:) = 0.    ;    mgrhv(:,jpj) = 0.
       DO jj = 1, jpjm1
          DO ji = 1, jpim1
-            mgrhu(ji,jj) = INT(  SIGN( 1.e0, fsdept_0(ji+1,jj,mbkt(ji+1,jj)) - fsdept_0(ji,jj,mbkt(ji,jj)) )  )
-            mgrhv(ji,jj) = INT(  SIGN( 1.e0, fsdept_0(ji,jj+1,mbkt(ji,jj+1)) - fsdept_0(ji,jj,mbkt(ji,jj)) )  )
+            mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
+            mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
          END DO
       END DO
 
       DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
-         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
-            e3u_bbl_0(ji,jj) = MIN( fse3u_0(ji,jj,mbkt(ji+1,jj  )), fse3u_0(ji,jj,mbkt(ji,jj)) )
-            e3v_bbl_0(ji,jj) = MIN( fse3v_0(ji,jj,mbkt(ji  ,jj+1)), fse3v_0(ji,jj,mbkt(ji,jj)) )
+         DO ji = 1, jpim1           ! minimum of top & bottom e3u_0 (e3v_0)
+            e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
+            e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
          END DO
       END DO
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90	(revision 3865)
@@ -110,6 +110,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) = re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk) * umask(ji,jj,jk)
+                  zeev(ji,jj) = re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk) * vmask(ji,jj,jk)
                END DO
             END DO
@@ -133,5 +133,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 / ( e12t(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)   )
@@ -151,5 +151,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 / ( e12t(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/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90	(revision 3865)
@@ -210,6 +210,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, jpim1
-                  zabe1 = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)
-                  zabe2 = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)
+                  zabe1 = re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk)
+                  zabe2 = re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk)
                   
                   zmku = 1./MAX( tmask(ji+1,jj,jk  )+tmask(ji,jj,jk+1)   &
@@ -279,5 +279,5 @@
             DO jk = 2, jpkm1
                DO ji = 2, jpim1
-                  zcof0 = e1t(ji,jj) * e2t(ji,jj) / fse3w(ji,jj,jk)   &
+                  zcof0 = e12t(ji,jj) / fse3w_n(ji,jj,jk)   &
                      &     * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)        &
                      &        + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
@@ -310,5 +310,5 @@
                   DO ji = 2, jpim1
                      ! eddy coef. divided by the volume element
-                     zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                     zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) )
                      ! vertical divergence
                      ztav = fsahtt(ji,jj,jk) * ( zftw(ji,jk) - zftw(ji,jk+1) )
@@ -322,5 +322,5 @@
                   DO ji = 2, jpim1
                      ! inverse of the volume element
-                     zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                     zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) )
                      ! vertical divergence
                      ztav = zftw(ji,jk) - zftw(ji,jk+1)
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 3865)
@@ -176,6 +176,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 ) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk)
+                  zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk)
                   !
                   zmsku = 1. / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   &
@@ -201,5 +201,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 / ( e12t(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
@@ -288,5 +288,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 / ( e12t(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/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90	(revision 3865)
@@ -31,6 +31,4 @@
 
    PUBLIC   tra_ldf_lap   ! routine called by step.F90
-
-   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:) ::   e1ur, e2vr   ! scale factor coefficients
 
    !! * Substitutions
@@ -85,15 +83,4 @@
          IF(lwp) WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
-         !
-         IF( .NOT. ALLOCATED( e1ur ) ) THEN
-            ! This routine may be called for both active and passive tracers. 
-            ! Allocate and set saved arrays on first call only.
-            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
       ENDIF
 
@@ -107,6 +94,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) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk)
+                  zabe2 = fsahtv(ji,jj,jk) * vmask(ji,jj,jk) * re1v_e2v(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) )
@@ -120,9 +107,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) * re2u_e1u(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) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,ikv)
                         ztv(ji,jj,jk) = zabe2 * pgv(ji,jj,jn)
                      ENDIF
@@ -136,5 +123,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 / ( e12t(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/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 3865)
@@ -448,8 +448,8 @@
 !CDIR NOVERRCHK   
                         DO ji = 1, jpi
-                           zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * xsi0r     )
-                           zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * zekb(ji,jj) )
-                           zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * zekg(ji,jj) )
-                           zc3 = ze3(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * zekr(ji,jj) )
+                           zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r     )
+                           zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekb(ji,jj) )
+                           zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekg(ji,jj) )
+                           zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekr(ji,jj) )
                            ze0(ji,jj,jk) = zc0
                            ze1(ji,jj,jk) = zc1
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/oce.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 3865)
@@ -36,7 +36,4 @@
    !! ------------                                      !  fields  ! fields ! trends !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   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
@@ -77,13 +74,10 @@
          &      rn2b (jpi,jpj,jpk)      , rn2  (jpi,jpj,jpk)                              , STAT=ierr(1) )
          !
-      ALLOCATE( rhd (jpi,jpj,jpk) ,                                         &
-         &      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),                     &
-         &      gru(jpi,jpj)      , grv(jpi,jpj)                      , STAT=ierr(2) )
+      ALLOCATE(rhd (jpi,jpj,jpk) ,                                         &
+         &     rhop(jpi,jpj,jpk) ,                                         &
+         &     sshb  (jpi,jpj)   , sshn  (jpi,jpj) , ssha  (jpi,jpj) ,     &
+         &     spgu  (jpi,jpj)   , spgv(jpi,jpj)   ,                       &
+         &     gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts),                     &
+         &     gru(jpi,jpj)      , grv(jpi,jpj)                      , STAT=ierr(2) )
          !
       ALLOCATE( snwice_mass(jpi,jpj)  , snwice_mass_b(jpi,jpj),             &
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 3865)
@@ -97,7 +97,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
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -139,5 +141,5 @@
       !
       IF( lk_ldfslp ) THEN                            ! slope of lateral mixing
-                         CALL eos( tsb, rhd )                ! before in situ density
+                         CALL eos( tsb, rhd, gdept_0(:,:,:) )             ! before in situ density
          IF( ln_zps )    CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
             &                                      rhd, gru , grv  )      ! of t, s, rd at the last ocean level
@@ -202,10 +204,10 @@
          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
                              CALL tra_nxt( kstp )                ! tracer fields at next time step
-                             CALL eos    ( tsa, rhd, rhop )      ! Time-filtered in situ density for hpg computation
+                             CALL eos    ( tsa, rhd, rhop, fsdept_n(:,:,:) )  ! Time-filtered in situ density for hpg computation
          IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsa, gtsu, gtsv,  &    ! zps: time filtered hor. derivative
             &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level
 
       ELSE                                                  ! centered hpg  (eos then time stepping)
-                             CALL eos    ( tsn, rhd, rhop )      ! now in situ density for hpg computation
+                             CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation
          IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &    ! zps: now hor. derivative
             &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level
@@ -238,5 +240,6 @@
                                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
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 3865)
@@ -64,5 +64,9 @@
    USE bdydyn3d         ! bdy cond. for baroclinic vel.  (bdy_dyn3d routine)
 
-   USE sshwzv           ! vertical velocity and ssh        (ssh_wzv routine)
+   USE sshwzv           ! vertical velocity and ssh        (ssh_nxt routine)
+   !                                                       (ssh_swp routine)
+   !                                                       (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)
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/C14b/trcwri_c14b.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/C14b/trcwri_c14b.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/C14b/trcwri_c14b.F90	(revision 3865)
@@ -36,5 +36,9 @@
       DO jn = jp_c14b0, jp_c14b1
          cltra = TRIM( ctrcnm(jn) )                  ! short title for tracer
-         CALL iom_put( cltra, trn(:,:,:,jn) )
+         IF( lk_vvl ) THEN
+            CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) )
+         ELSE
+            CALL iom_put( cltra, trn(:,:,:,jn) )
+         ENDIF
       END DO
       !
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/CFC/trcwri_cfc.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/CFC/trcwri_cfc.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/CFC/trcwri_cfc.F90	(revision 3865)
@@ -36,5 +36,9 @@
       DO jn = jp_cfc0, jp_cfc1
          cltra = TRIM( ctrcnm(jn) )                  ! short title for tracer
-         CALL iom_put( cltra, trn(:,:,:,jn) )
+         IF( lk_vvl ) THEN
+            CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) )
+         ELSE
+            CALL iom_put( cltra, trn(:,:,:,jn) )
+         ENDIF
       END DO
       !
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcwri_my_trc.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcwri_my_trc.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcwri_my_trc.F90	(revision 3865)
@@ -36,5 +36,9 @@
       DO jn = jp_myt0, jp_myt1
          cltra = TRIM( ctrcnm(jn) )                  ! short title for tracer
-         CALL iom_put( cltra, trn(:,:,:,jn) )
+         IF( lk_vvl ) THEN
+            CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) )
+         ELSE
+            CALL iom_put( cltra, trn(:,:,:,jn) )
+         ENDIF
       END DO
       !
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/PISCES/trcwri_pisces.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/PISCES/trcwri_pisces.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/PISCES/trcwri_pisces.F90	(revision 3865)
@@ -39,4 +39,9 @@
       DO jn = jp_pcs0, jp_pcs1
          cltra = TRIM( ctrcnm(jn) )                  ! short title for tracer
+         IF( lk_vvl ) THEN
+            CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) )
+         ELSE
+            CALL iom_put( cltra, trn(:,:,:,jn) )
+         ENDIF
          CALL iom_put( cltra, trn(:,:,:,jn) * zrfact )
       END DO
@@ -47,5 +52,9 @@
          IF( jn == jppo4  )                 zrfact = po4r * 1.0e+6
          cltra = TRIM( ctrcnm(jn) )                  ! short title for tracer
-         CALL iom_put( cltra, trn(:,:,:,jn) * zrfact )
+         IF( lk_vvl ) THEN
+            CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) * zrfact )
+         ELSE
+            CALL iom_put( cltra, trn(:,:,:,jn) * zrfact )
+         ENDIF
       END DO
 #endif
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90	(revision 3864)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90	(revision 3865)
@@ -77,11 +77,4 @@
    USE oce , ONLY :   sshb    =>    sshb    !: sea surface height at t-point [m]   
    USE oce , ONLY :   ssha    =>    ssha    !: sea surface height at t-point [m]   
-   USE oce , ONLY :   sshu_n  =>    sshu_n  !: sea surface height at u-point [m]   
-   USE oce , ONLY :   sshu_b  =>    sshu_b  !: sea surface height at u-point [m]   
-   USE oce , ONLY :   sshu_a  =>    sshu_a  !: sea surface height at u-point [m]   
-   USE oce , ONLY :   sshv_n  =>    sshv_n  !: sea surface height at v-point [m]   
-   USE oce , ONLY :   sshv_b  =>    sshv_b  !: sea surface height at v-point [m]   
-   USE oce , ONLY :   sshv_a  =>    sshv_a  !: sea surface height at v-point [m]   
-   USE oce , ONLY :   sshf_n  =>    sshf_n  !: sea surface height at v-point [m]   
    USE oce , ONLY :   l_traldf_rot => l_traldf_rot  !: rotated laplacian operator for lateral diffusion
 #if defined key_offline
