Index: branches/dev_004_VVL/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 1407)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 1408)
@@ -127,5 +127,6 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !:
       hur, hvr,          &  !: inverse of u and v-points ocean depth (1/m)
-      hu , hv               !: depth at u- and v-points (meters)
+      hu , hv,           &  !: depth at u- and v-points (meters)
+      hu_0 , hv_0           !: refernce depth at u- and v-points (meters)
 
    !! z-coordinate with full steps (also used in the other cases as reference z-coordinate)
Index: branches/dev_004_VVL/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/DOM/domain.F90	(revision 1407)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/DOM/domain.F90	(revision 1408)
@@ -4,5 +4,14 @@
    !! Ocean initialization : domain initialization
    !!==============================================================================
-
+   !! History :  OPA  !  1990-10  (C. Levy - G. Madec)  Original code
+   !!                 !  1991-11  (G. Madec)
+   !!                 !  1992-01  (M. Imbard) insert time step initialization
+   !!                 !  1996-06  (G. Madec) generalized vertical coordinate 
+   !!                 !  1997-02  (G. Madec) creation of domwri.F
+   !!                 !  2001-05  (E.Durand - G. Madec) insert closed sea
+   !!   NEMO     1.0  !  2002-08  (G. Madec)  F90: Free form and module
+   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization
+   !!----------------------------------------------------------------------
+   
    !!----------------------------------------------------------------------
    !!   dom_init       : initialize the space and time domain
@@ -10,5 +19,4 @@
    !!   dom_ctl        : control print for the ocean domain
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! 
    USE dom_oce         ! ocean space and time domain
@@ -30,14 +38,13 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC dom_init       ! called by opa.F90
+   PUBLIC   dom_init   ! called by opa.F90
 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
-   !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!-------------------------------------------------------------------------
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
-   !!----------------------------------------------------------------------
+   !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!-------------------------------------------------------------------------
 
 CONTAINS
@@ -58,16 +65,5 @@
       !!      - dom_stp: defined the model time step
       !!      - dom_wri: create the meshmask file if nmsh=1
-      !!
-      !! History :
-      !!        !  90-10  (C. Levy - G. Madec)  Original code
-      !!        !  91-11  (G. Madec)
-      !!        !  92-01  (M. Imbard) insert time step initialization
-      !!        !  96-06  (G. Madec) generalized vertical coordinate 
-      !!        !  97-02  (G. Madec) creation of domwri.F
-      !!        !  01-05  (E.Durand - G. Madec) insert closed sea
-      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
-      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
-      !!----------------------------------------------------------------------
-      !! * Local declarations
+      !!----------------------------------------------------------------------
       INTEGER ::   jk                ! dummy loop argument
       INTEGER ::   iconf = 0         ! temporary integers
@@ -90,12 +86,11 @@
       CALL dom_msk                        ! Masks
 
-      IF( lk_vvl )   CALL dom_vvl_ini     ! Vertical variable mesh
+      IF( lk_vvl )   CALL dom_vvl         ! Vertical variable mesh
 
       ! Local depth or Inverse of the local depth of the water column at u- and v-points
       ! ------------------------------
       ! Ocean depth at U- and V-points
-      hu(:,:) = 0.
-      hv(:,:) = 0.
-
+      hu(:,:) = 0.e0
+      hv(:,:) = 0.e0
       DO jk = 1, jpk
          hu(:,:) = hu(:,:) + fse3u(:,:,jk) * umask(:,:,jk)
@@ -105,15 +100,12 @@
       hur(:,:) = fse3u(:,:,1)             ! Lower bound : thickness of the first model level
       hvr(:,:) = fse3v(:,:,1)
-
       DO jk = 2, jpk                      ! Sum of the vertical scale factors
          hur(:,:) = hur(:,:) + fse3u(:,:,jk) * umask(:,:,jk)
          hvr(:,:) = hvr(:,:) + fse3v(:,:,jk) * vmask(:,:,jk)
       END DO
-
       ! Compute and mask the inverse of the local depth
       hur(:,:) = 1. / hur(:,:) * umask(:,:,1)
       hvr(:,:) = 1. / hvr(:,:) * vmask(:,:,1)
 
-
       CALL dom_stp                        ! Time step
 
@@ -121,5 +113,5 @@
 
       IF( .NOT.ln_rstart )   CALL dom_ctl    ! Domain control
-
+      !
    END SUBROUTINE dom_init
 
@@ -134,9 +126,5 @@
       !!              - namdom namelist
       !!              - namcla namelist
-      !!
-      !! History :
-      !!   9.0  !  03-08  (G. Madec)  Original code
-      !!----------------------------------------------------------------------
-      !! * Modules used
+      !!----------------------------------------------------------------------
       USE ioipsl
       NAMELIST/namrun/ no    , cexper, cn_ocerst_in, cn_ocerst_out, ln_rstart, nrstdt,   &
@@ -156,8 +144,6 @@
       ENDIF
 
-      ! Namelist namrun : parameters of the run
-      REWIND( numnam )
+      REWIND( numnam )              ! Namelist namrun : parameters of the run
       READ  ( numnam, namrun )
-
       IF(lwp) THEN
          WRITE(numout,*) '        Namelist namrun'
@@ -228,6 +214,5 @@
       ENDIF
 
-      ! Namelist namdom : space/time domain (bathymetry, mesh, timestep)
-      REWIND( numnam )
+      REWIND( numnam )              ! Namelist namdom : space/time domain (bathymetry, mesh, timestep)
       READ  ( numnam, namdom )
 
@@ -252,9 +237,6 @@
       ENDIF
 
-      ! Default values
       n_cla = 0
-
-      ! Namelist cross land advection
-      REWIND( numnam )
+      REWIND( numnam )              ! Namelist cross land advection
       READ  ( numnam, namcla )
       IF(lwp) THEN
@@ -264,8 +246,6 @@
       ENDIF
 
-      IF( nbit_cmp == 1 .AND. n_cla /= 0 ) THEN
-         CALL ctl_stop( ' Reproductibility tests (nbit_cmp=1) require n_cla = 0' )
-      END IF
-
+      IF( nbit_cmp == 1 .AND. n_cla /= 0 )   CALL ctl_stop( ' Reproductibility tests (nbit_cmp=1) require n_cla = 0' )
+      !
    END SUBROUTINE dom_nam
 
@@ -278,15 +258,9 @@
       !!
       !! ** Method  :   compute and print extrema of masked scale factors
-      !!
-      !! History :
-      !!   8.5  !  02-08  (G. Madec)    Original code
-      !!----------------------------------------------------------------------
-      !! * Local declarations
+      !!----------------------------------------------------------------------
       INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
       INTEGER, DIMENSION(2) ::   iloc      ! 
       REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
       !!----------------------------------------------------------------------
-
-      ! Extrema of the scale factors
 
       IF(lwp)WRITE(numout,*)
@@ -325,5 +299,5 @@
          WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
       ENDIF
-
+      !
    END SUBROUTINE dom_ctl
 
Index: branches/dev_004_VVL/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 1407)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 1408)
@@ -4,39 +4,29 @@
    !! Ocean : 
    !!======================================================================
-   !! History :   9.0  !  06-06  (B. Levier, L. Marie)  original code
-   !!              "   !  07-07  (D. Storkey) Bug fixes and code for BDY option.
+   !! 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
    !!----------------------------------------------------------------------
    !!   'key_vvl'                              variable volume
    !!----------------------------------------------------------------------
+   !!   dom_vvl     : defined coefficients to distribute ssh on each layers
    !!----------------------------------------------------------------------
-   !!   dom_vvl         : empty routine
-   !!   dom_vvl_ini     : defined coefficients to distribute ssh on each layers
-   !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
    USE sbc_oce         ! surface boundary condition: ocean
-   USE dynspg_oce      ! surface pressure gradient variables
    USE phycst          ! physical constants
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! distributed memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE bdy_oce         ! unstructured open boundary conditions
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC dom_vvl_ini    ! called by dom_init.F90
-   PUBLIC dom_vvl        ! called by istate.F90 and step.F90
+   PUBLIC dom_vvl        ! called by domain.F90
 
-   !! * Module variables
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hu_0, hv_0
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ee_t, ee_u, ee_v, ee_f   !:
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ee_t, ee_u, ee_v, ee_f   !: ???
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !:
-      mut, muu, muv, muf                            !:
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   mut, muu, muv, muf   !: ??? 
 
    REAL(wp), DIMENSION(jpk) ::   r2dt               ! vertical profile time-step, = 2 rdttra 
@@ -47,5 +37,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -54,9 +44,7 @@
 CONTAINS       
 
-#if defined key_vvl
-
-   SUBROUTINE dom_vvl_ini
+   SUBROUTINE dom_vvl
       !!----------------------------------------------------------------------
-      !!                ***  ROUTINE dom_vvl_ini  ***
+      !!                ***  ROUTINE dom_vvl  ***
       !!                   
       !! ** Purpose :  compute coefficients muX at T-U-V-F points to spread
@@ -70,6 +58,6 @@
       IF(lwp)   THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'dom_vvl_ini : Variable volume activated'
-         WRITE(numout,*) '~~~~~~~~~~~   compute coef. used to spread ssh over each layers'
+         WRITE(numout,*) 'dom_vvl : Variable volume activated'
+         WRITE(numout,*) '~~~~~~~~  compute coef. used to spread ssh over each layers'
       ENDIF
 
@@ -127,19 +115,10 @@
 
 
-!!debug print
-!         ii=50   ;   ij = 50
-!         do jk= 1, jpk 
-!         WRITE(numout,*) 'domvvl GM  : h0=', SUM( fse3t_0(ii,ij,1:jk) * tmask(ii,ij,1:jk) ), 'e3t0=', fse3t_0(ii,ij,jk),   &
-!            &            'e3t =', fse3t_0(ii,ij,jk) * ( 1 +  mut(ii,ij,jk) ), 'mut', mut(ii,ij,jk),   &
-!            &            'h =', SUM( fse3t_0(ii,ij,1:jk) * ( 1 +  mut(ii,ij,1:jk) ) * tmask(ii,ij,1:jk) )
-!         end do
-!!end debug print
-
       ! Reference ocean depth at U- and V-points
       hu_0(:,:) = 0.e0    
       hv_0(:,:) = 0.e0
       DO jk = 1, jpk
-         hu_0(:,:) = hu_0(:,:) + fse3u(:,:,jk) * umask(:,:,jk)
-         hv_0(:,:) = hv_0(:,:) + fse3v(:,:,jk) * vmask(:,:,jk)
+         hu_0(:,:) = hu_0(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk)
+         hv_0(:,:) = hv_0(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk)
       END DO
 
@@ -149,6 +128,5 @@
             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 * fmask(ji,jj,1)
-!!gm bug used of fmask, even if thereafter multiplied by muf  which is correctly masked)
+            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) )   
@@ -170,28 +148,4 @@
       CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. )
       !
-
-
-   END SUBROUTINE dom_vvl_ini
-
-
-   SUBROUTINE dom_vvl
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE dom_vvl  ***
-      !!                   
-      !! ** Purpose :  compute ssh at U-V-F points, T-W scale factors and local
-      !!               depths at each time step.
-      !!----------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
-      !!----------------------------------------------------------------------
-
- !     IF( kt == nit000 ) THEN
- !        IF(lwp) WRITE(numout,*)
- !        IF(lwp) WRITE(numout,*) 'dom_vvl : '
- !        IF(lwp) WRITE(numout,*) '~~~~~~~ '
- !     ENDIF
-
-       WRITE(*,*) 'dom_vvl  : empty routine, you should not be here'
-
    END SUBROUTINE dom_vvl
 
@@ -200,6 +154,5 @@
    !!   Default option :                                      Empty routine
    !!----------------------------------------------------------------------
-   SUBROUTINE dom_vvl_ini
-   END SUBROUTINE dom_vvl_ini
+CONTAINS
    SUBROUTINE dom_vvl
    END SUBROUTINE dom_vvl
Index: branches/dev_004_VVL/NEMO/OPA_SRC/DOM/domzgr_substitute.h90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/DOM/domzgr_substitute.h90	(revision 1407)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/DOM/domzgr_substitute.h90	(revision 1408)
@@ -78,4 +78,5 @@
 #   define  fse3uw_a(i,j,k)  (fse3uw_0(i,j,k)*(1+sshu_a(i,j)*muu(i,j,k)))
 #   define  fse3vw_a(i,j,k)  (fse3vw_0(i,j,k)*(1+sshv_a(i,j)*muv(i,j,k)))
+
 #else
 ! z- or s-coordinate (1D or 3D + no time dependency) use reference in all cases
Index: branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 1407)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 1408)
@@ -15,5 +15,5 @@
    !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization
    !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines. 
-   !!            3.1  !  2009-02  (G. Madec)  re-introduce the vvl option
+   !!            3.2  !  2009-04  (G. Madec, R.Benshila))  re-introduce the vvl option
    !!----------------------------------------------------------------------
   
@@ -21,5 +21,4 @@
    !!   dyn_nxt      : update the horizontal velocity from the momentum trend
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
@@ -42,13 +41,13 @@
    PRIVATE
 
-   !! * Accessibility
-   PUBLIC dyn_nxt                ! routine called by step.F90
+   PUBLIC    dyn_nxt   ! routine called by step.F90
+
    !! * Substitutions
 #  include "domzgr_substitute.h90"
-   !!----------------------------------------------------------------------
-   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!-------------------------------------------------------------------------
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
    !! $Id$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
-   !!----------------------------------------------------------------------
+   !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
+   !!-------------------------------------------------------------------------
 
 CONTAINS
@@ -85,5 +84,6 @@
       REAL(wp) ::   ze3u_b, ze3u_n, ze3u_a   !    -         -
       REAL(wp) ::   ze3v_b, ze3v_n, ze3v_a   !    -         - 
-      !!----------------------------------------------------------------------
+      REAL(wp) ::   zuf   , zvf              !    -         - 
+     !!----------------------------------------------------------------------
 
       IF( kt == nit000 ) THEN
@@ -100,6 +100,5 @@
 
       ! Lateral boundary conditions on ( ua, va )
-      CALL lbc_lnk( ua, 'U', -1. )
-      CALL lbc_lnk( va, 'V', -1. )
+      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. ) 
 
       ! Next velocity
@@ -110,8 +109,8 @@
       IF( lk_vvl ) THEN                                  ! Varying levels
          DO jk = 1, jpkm1
-            ua(:,:,jk) = (        ub(:,:,jk) * fse3u_b(:,:,jk)     &
+            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)     &
                &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) )   &
                &         / fse3u_a(:,:,jk) * umask(:,:,jk)
-            va(:,:,jk) = (        vb(:,:,jk) * fse3v_b(:,:,jk)     &
+            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)     &
                &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) )   &
                &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
@@ -130,5 +129,5 @@
 
       IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
-         !Flather boundary condition :
+         ! Flather boundary condition :
          !        - Update sea surface height on each open boundary
          !                 sshn (= after ssh) for explicit case
@@ -136,55 +135,44 @@
          !        - Correct the barotropic velocities
          CALL obc_dyn_bt( kt )
-
-         !Boundary conditions on sshn ( after ssh)
+         !
+         ! Boundary conditions on sshn ( after ssh)
          CALL lbc_lnk( sshn, 'T', 1. )
-
-         IF(ln_ctl) THEN         ! print sum trends (used for debugging)
-            CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask)
-         ENDIF
-
-         IF ( ln_vol_cst ) CALL obc_vol( kt )
-
+         !
+         IF( ln_vol_cst )   CALL obc_vol( kt )
+         !
+         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
       ENDIF
 
 # elif defined key_bdy 
       ! Update (ua,va) along open boundaries (for exp or ts options).
-      IF ( lk_dynspg_exp .or. lk_dynspg_ts ) THEN
-  
+      IF( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
+         !
          CALL bdy_dyn_frs( kt )
-
-         IF ( ln_bdy_fla ) THEN
-
-            ua_e(:,:)=0.0
-            va_e(:,:)=0.0
-
+         !
+         IF( ln_bdy_fla ) THEN
+            ua_e(:,:) = 0.e0
+            va_e(:,:) = 0.e0
             ! Set these variables for use in bdy_dyn_fla
             hu_e(:,:) = hu(:,:)
             hv_e(:,:) = hv(:,:)
-
-            DO jk = 1, jpkm1
-               !! Vertically integrated momentum trends
+            DO jk = 1, jpkm1   !! Vertically integrated momentum trends
                ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
                va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
             END DO
-
             DO jk = 1 , jpkm1
                ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:) * hur(:,:)
                va(:,:,jk) = va(:,:,jk) - va_e(:,:) * hvr(:,:)
             END DO
-
             CALL bdy_dta_bt( kt+1, 0)
             CALL bdy_dyn_fla
-
          ENDIF
-
+         !
          DO jk = 1 , jpkm1
             ua(:,:,jk) = ua(:,:,jk) + ua_e(:,:) * hur(:,:)
             va(:,:,jk) = va(:,:,jk) + va_e(:,:) * hvr(:,:)
          END DO
-
-      ENDIF
-
+      ENDIF
 # endif
+
 # if defined key_agrif
       CALL Agrif_dyn( kt )
@@ -196,6 +184,4 @@
       IF( neuler == 0 .AND. kt == nit000 ) THEN
          DO jk = 1, jpkm1                                 
-            ub(:,:,jk) = un(:,:,jk) 
-            vb(:,:,jk) = vn(:,:,jk)
             un(:,:,jk) = ua(:,:,jk)
             vn(:,:,jk) = va(:,:,jk)
@@ -230,17 +216,17 @@
             END DO
          ELSE                                             ! Fixed levels
-!RB_vvl : should be done as in tranxt ?
             DO jk = 1, jpkm1                              ! filter applied on velocities
-               ub(:,:,jk) = atfp * ( ub(:,:,jk) + ua(:,:,jk) ) + atfp1 * un(:,:,jk)
-               vb(:,:,jk) = atfp * ( vb(:,:,jk) + va(:,:,jk) ) + atfp1 * vn(:,:,jk)
-               un(:,:,jk) = ua(:,:,jk)
-               vn(:,:,jk) = va(:,:,jk)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zuf = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk)
+                     zvf = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk)
+                     ub(ji,jj,jk) = zuf
+                     vb(ji,jj,jk) = zvf
+                     un(ji,jj,jk) = ua(ji,jj,jk)
+                     vn(ji,jj,jk) = va(ji,jj,jk)
+                  END DO
+               END DO
             END DO
          ENDIF
-      ENDIF
-
-      IF(ln_ctl)   THEN
-         CALL prt_ctl(tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask, &
-            &         tab3d_2=vn, clinfo2=' Vn: ', mask2=vmask)
       ENDIF
 
@@ -249,4 +235,7 @@
 #endif      
 
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
+         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
+      ! 
    END SUBROUTINE dyn_nxt
 
Index: branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_exp.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_exp.F90	(revision 1407)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_exp.F90	(revision 1408)
@@ -4,4 +4,6 @@
    !! Ocean dynamics:  surface pressure gradient trend
    !!======================================================================
+   !! History :  9.0  !  2005-11  (V. Garnier, G. Madec, L. Bessieres) Original code
+   !!----------------------------------------------------------------------
 #if defined key_dynspg_exp   ||   defined key_esopa
    !!----------------------------------------------------------------------
@@ -12,5 +14,4 @@
    !!                      volume case with vector optimization
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers 
    USE dom_oce         ! ocean space and time domain 
@@ -30,6 +31,5 @@
    PRIVATE
 
-   !! * Accessibility
-   PUBLIC dyn_spg_exp  ! routine called by step.F90
+   PUBLIC   dyn_spg_exp   ! routine called by step.F90
 
    !! * Substitutions
@@ -37,9 +37,8 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    !!----------------------------------------------------------------------
-
 
 CONTAINS
@@ -62,15 +61,7 @@
       !!
       !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT( in )  ::   kt         ! ocean time-step index
       !!
-      !! References :
-      !!
-      !! History :
-      !!   9.0  !  05-11  (V. Garnier, G. Madec, L. Bessieres) Original code
-      !!
-      !!---------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in )  ::   kt         ! ocean time-step index
-
-      !! * Local declarations
       INTEGER  ::   ji, jj, jk               ! dummy loop indices
       !!----------------------------------------------------------------------
@@ -86,11 +77,8 @@
       ENDIF
 
-      ! 0. Initialization
-      ! -----------------
       ! read or estimate sea surface height and vertically integrated velocities
       IF( lk_obc )   CALL obc_dta_bt( kt, 0 )
 
-      ! 1. Surface pressure gradient (now)
-      ! ----------------------------
+      ! Surface pressure gradient (now)
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -100,6 +88,5 @@
       END DO
 
-      ! 2. Add the surface pressure trend to the general trend
-      ! ------------------------------------------------------
+      ! Add the surface pressure trend to the general trend
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1
Index: branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_flt.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 1407)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 1408)
@@ -4,17 +4,17 @@
    !! Ocean dynamics:  surface pressure gradient trend
    !!======================================================================
-   !! History    8.0  !  98-05  (G. Roullet)  free surface
-   !!                 !  98-10  (G. Madec, M. Imbard)  release 8.2
-   !!            8.5  !  02-08  (G. Madec)  F90: Free form and module
-   !!            " "  !  02-11  (C. Talandier, A-M Treguier) Open boundaries
-   !!            9.0  !  04-08  (C. Talandier) New trends organization
-   !!            " "  !  05-11  (V. Garnier) Surface pressure gradient organization
-   !!            " "  !  06-07  (S. Masson)  distributed restart using iom
-   !!            " "  !  05-01  (J.Chanut, A.Sellar) Calls to BDY routines. 
+   !! History    OPA  !  1998-05  (G. Roullet)  free surface
+   !!                 !  1998-10  (G. Madec, M. Imbard)  release 8.2
+   !!   NEMO     O.1  !  2002-08  (G. Madec)  F90: Free form and module
+   !!             -   !  2002-11  (C. Talandier, A-M Treguier) Open boundaries
+   !!            1.0  !  2004-08  (C. Talandier) New trends organization
+   !!             -   !  2005-11  (V. Garnier) Surface pressure gradient organization
+   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
+   !!             -   !  2006-08  (J.Chanut, A.Sellar) Calls to BDY routines. 
+   !!            3.2  !  2009-03  (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module
    !!----------------------------------------------------------------------
 #if defined key_dynspg_flt   ||   defined key_esopa  
    !!----------------------------------------------------------------------
    !!   'key_dynspg_flt'                              filtered free surface
-   !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
    !!   dyn_spg_flt  : update the momentum trend with the surface pressure
@@ -53,6 +53,6 @@
    PRIVATE
 
-   PUBLIC dyn_spg_flt  ! routine called by step.F90
-   PUBLIC flt_rst      ! routine called by istate.F90
+   PUBLIC   dyn_spg_flt  ! routine called by step.F90
+   PUBLIC   flt_rst      ! routine called by istate.F90
 
    !! * Substitutions
@@ -60,5 +60,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
@@ -80,5 +80,5 @@
       !!         spgv = 1/rau0 d/dy(ps) =  1/e2v dj( sshn + rnu btda )
       !!      where sshn is the free surface elevation and btda is the after
-      !!      of the free surface elevation
+      !!      time derivative of the free surface elevation
       !!       -1- evaluate the surface presure trend (including the addi-
       !!      tional force) in three steps:
@@ -104,15 +104,14 @@
       !! References : Roullet and Madec 1999, JGR.
       !!---------------------------------------------------------------------
-      !! * Modules used
-      USE oce    , ONLY :   zub   => ta,             &  ! ta used as workspace
-                            zvb   => sa                 ! sa   "          "
-
-      INTEGER, INTENT( in )  ::   kt         ! ocean time-step index
-      INTEGER, INTENT( out ) ::   kindic     ! solver convergence flag (<0 if not converge)
+      USE oce, ONLY :   zub   => ta   ! ta used as workspace
+      USE oce, ONLY :   zvb   => sa   ! ta used as workspace
+      !!
+      INTEGER, INTENT( in )  ::   kt       ! ocean time-step index
+      INTEGER, INTENT( out ) ::   kindic   ! solver convergence flag (<0 if not converge)
       !!                                   
       INTEGER  ::   ji, jj, jk                          ! dummy loop indices
-      REAL(wp) ::   z2dt, z2dtg, zraur, znugdt,      &  ! temporary scalars
-         &          znurau, zgcb, zbtd,              &  !   "          "
-         &          ztdgu, ztdgv                        !   "          "
+      REAL(wp) ::   z2dt, z2dtg, zraur, znugdt          ! temporary scalars
+      REAL(wp) ::   znurau, zgcb, zbtd                  !   "          "
+      REAL(wp) ::   ztdgu, ztdgv                        !   "          "
       !!----------------------------------------------------------------------
       !
@@ -130,5 +129,5 @@
          ! when using agrif, sshn, gcx have to be read in istate
          IF (.NOT. lk_agrif) CALL flt_rst( nit000, 'READ' )       ! read or initialize the following fields:
-         !                                                        ! gcx, gcxb, sshb, sshn
+         !                                                        ! gcx, gcxb
       ENDIF
 
@@ -354,6 +353,6 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ua(ji,jj,jk) = (ua(ji,jj,jk) + spgu(ji,jj)) * umask(ji,jj,jk)
-               va(ji,jj,jk) = (va(ji,jj,jk) + spgv(ji,jj)) * vmask(ji,jj,jk)
+               ua(ji,jj,jk) = ( ua(ji,jj,jk) + spgu(ji,jj) ) * umask(ji,jj,jk)
+               va(ji,jj,jk) = ( va(ji,jj,jk) + spgv(ji,jj) ) * vmask(ji,jj,jk)
             END DO
          END DO
@@ -363,7 +362,4 @@
       ! --------------------------------------------------
       IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' )
-
-      ! print sum trends (used for debugging)
-      IF(ln_ctl)     CALL prt_ctl( tab2d_1=sshn, clinfo1=' spg  - ssh: ', mask1=tmask )
       !
    END SUBROUTINE dyn_spg_flt
@@ -386,7 +382,5 @@
            CALL iom_get( numror, jpdom_autoglo, 'gcx' , gcx (1:jpi,1:jpj) )
            CALL iom_get( numror, jpdom_autoglo, 'gcxb', gcxb(1:jpi,1:jpj) )
-           IF( neuler == 0 ) THEN
-              gcxb(:,:) = gcx (:,:)
-           ENDIF
+           IF( neuler == 0 )   gcxb(:,:) = gcx (:,:)
         ELSE
            gcx (:,:) = 0.e0
@@ -396,5 +390,5 @@
 ! Caution : extra-hallow
 ! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
-        CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx( 1:jpi,1:jpj) )
+        CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) )
         CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
      ENDIF
Index: branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 1407)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 1408)
@@ -1,10 +1,10 @@
 MODULE dynspg_ts
    !!======================================================================
-   !! History :   9.0  !  04-12  (L. Bessieres, G. Madec)  Original code
-   !!             " "  !  05-11  (V. Garnier, G. Madec)  optimization
-   !!             9.0  !  06-08  (S. Masson)  distributed restart using iom
-   !!              "   !  08-01  (R. Benshila)  change averaging method
-   !!              "   !  07-07  (D. Storkey) calls to BDY routines
-   !!---------------------------------------------------------------------
+   !! History :   1.0  !  04-12  (L. Bessieres, G. Madec)  Original code
+   !!              -   !  05-11  (V. Garnier, G. Madec)  optimization
+   !!              -   !  06-08  (S. Masson)  distributed restart using iom
+   !!             2.0  !  07-07  (D. Storkey) calls to BDY routines
+   !!              -   !  08-01  (R. Benshila)  change averaging method
+  !!---------------------------------------------------------------------
 #if defined key_dynspg_ts   ||   defined key_esopa
    !!----------------------------------------------------------------------
@@ -16,5 +16,4 @@
    !!   ts_rst      : read/write the time-splitting restart fields in the ocean restart file
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
@@ -46,16 +45,15 @@
    PUBLIC ts_rst      ! routine called by istate.F90
 
-   REAL(wp), DIMENSION(jpi,jpj) ::  ftnw, ftne,   &  ! triad of coriolis parameter
-      &                             ftsw, ftse       ! (only used with een vorticity scheme)
-
+   REAL(wp), DIMENSION(jpi,jpj) ::  ftnw, ftne   ! triad of coriolis parameter
+   REAL(wp), DIMENSION(jpi,jpj) ::  ftsw, ftse   ! (only used with een vorticity scheme)
 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!-------------------------------------------------------------------------
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
-   !!----------------------------------------------------------------------
+   !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
+   !!-------------------------------------------------------------------------
 
 CONTAINS
@@ -90,6 +88,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT( in )  ::   kt           ! ocean time-step index
-
-      !! * Local declarations
+      !!
       INTEGER  ::  ji, jj, jk, jit             ! dummy loop indices
       INTEGER  ::  icycle                      ! temporary scalar
Index: branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynvor.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 1407)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 1408)
@@ -5,12 +5,13 @@
    !!                 planetary vorticity trends
    !!======================================================================
-   !! History :  1.0  !  89-12  (P. Andrich)  vor_ens: Original code
-   !!            5.0  !  91-11  (G. Madec) vor_ene, vor_mix: Original code
-   !!            6.0  !  96-01  (G. Madec)  s-coord, suppress work arrays
-   !!            8.5  !  02-08  (G. Madec)  F90: Free form and module
-   !!            8.5  !  04-02  (G. Madec)  vor_een: Original code
-   !!            9.0  !  03-08  (G. Madec)  vor_ctl: Original code
-   !!            9.0  !  05-11  (G. Madec)  dyn_vor: Original code (new step architecture)
-   !!            9.0  !  06-11  (G. Madec)  flux form advection: add metric term
+   !! History :  OPA  !  1989-12  (P. Andrich)  vor_ens: Original code
+   !!            5.0  !  1991-11  (G. Madec) vor_ene, vor_mix: Original code
+   !!            6.0  !  1996-01  (G. Madec)  s-coord, suppress work arrays
+   !!            8.5  !  2002-08  (G. Madec)  F90: Free form and module
+   !!   NEMO     1.0  !  2004-02  (G. Madec)  vor_een: Original code
+   !!             -   !  2003-08  (G. Madec)  add vor_ctl
+   !!             -   !  2005-11  (G. Madec)  add dyn_vor (new step architecture)
+   !!            2.0  !  2006-11  (G. Madec)  flux form advection: add metric term
+   !!            3.2  !  2009-04  (R. Benshila)  vvl: correction of een scheme
    !!----------------------------------------------------------------------
 
@@ -37,5 +38,5 @@
    PUBLIC   dyn_vor   ! routine called by step.F90
 
-   !!* Namelist nam_dynvor: vorticity term
+   !                                             !!* Namelist nam_dynvor: vorticity term
    LOGICAL, PUBLIC ::   ln_dynvor_ene = .FALSE.   !: energy conserving scheme
    LOGICAL, PUBLIC ::   ln_dynvor_ens = .TRUE.    !: enstrophy conserving scheme
@@ -52,5 +53,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2006) 
+   !! NEMO/OPA 3,2 , LOCEAN-IPSL (2009) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -174,5 +175,5 @@
       IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask, &
          &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-
+      !
    END SUBROUTINE dyn_vor
 
@@ -206,5 +207,5 @@
       INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
       INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
-         !                                                        ! =nrvm (relative vorticity or metric)
+      !                                                           ! =nrvm (relative vorticity or metric)
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
@@ -222,6 +223,5 @@
       ENDIF
 
-      ! Local constant initialization
-      zfact2 = 0.5 * 0.5
+      zfact2 = 0.5 * 0.5      ! Local constant initialization
 
 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
@@ -229,4 +229,5 @@
       DO jk = 1, jpkm1                                 ! Horizontal slab
          !                                             ! ===============
+         !
          ! Potential vorticity and horizontal fluxes
          ! -----------------------------------------
@@ -315,5 +316,5 @@
       INTEGER, INTENT(in) ::   kt   ! ocean timestep index
       !!
-      INTEGER ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   zfact1, zua, zcua, zx1, zy1   ! temporary scalars
       REAL(wp) ::   zfact2, zva, zcva, zx2, zy2   !    "         "
@@ -327,6 +328,5 @@
       ENDIF
 
-      ! Local constant initialization
-      zfact1 = 0.5 * 0.25
+      zfact1 = 0.5 * 0.25      ! Local constant initialization
       zfact2 = 0.5 * 0.5
 
@@ -335,5 +335,5 @@
       DO jk = 1, jpkm1                                 ! Horizontal slab
          !                                             ! ===============
-
+         !
          ! Relative and planetary potential vorticity and horizontal fluxes
          ! ----------------------------------------------------------------
@@ -438,6 +438,5 @@
       ENDIF
 
-      ! Local constant initialization
-      zfact1 = 0.5 * 0.25
+      zfact1 = 0.5 * 0.25      ! Local constant initialization
 
 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
@@ -445,4 +444,5 @@
       DO jk = 1, jpkm1                                 ! Horizontal slab
          !                                             ! ===============
+         !
          ! Potential vorticity and horizontal fluxes
          ! -----------------------------------------
@@ -465,10 +465,10 @@
                        &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               &
+                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                &
                        &       )
                END DO
             END DO
          END SELECT
-
+         !
          IF( ln_sco ) THEN
             DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop 
@@ -487,5 +487,5 @@
             END DO
          ENDIF
-
+         !
          ! Compute and add the vorticity term trend
          ! ----------------------------------------
@@ -514,12 +514,8 @@
       !!
       !! ** Method  :   Trend evaluated using now fields (centered in time) 
-      !!      and the Arakawa and Lamb (19XX) flux form formulation : conserves 
+      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves 
       !!      both the horizontal kinetic energy and the potential enstrophy
-      !!      when horizontal divergence is zero.
-      !!      The trend of the vorticity term is given by:
-      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
-      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
-      !!      Add this trend to the general momentum trend (ua,va):
-      !!          (ua,va) = (ua,va) + ( voru , vorv )
+      !!      when horizontal divergence is zero (see the NEMO documentation)
+      !!      Add this trend to the general momentum trend (ua,va).
       !!
       !! ** Action : - Update (ua,va) with the now vorticity term trend
@@ -531,13 +527,17 @@
       INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
       INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
-         !                                                        ! =nrvm (relative vorticity or metric)
+      !                                                           ! =nrvm (relative vorticity or metric)
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
       !!
-      INTEGER ::   ji, jj, jk          ! dummy loop indices
+      INTEGER  ::   ji, jj, jk         ! dummy loop indices
       REAL(wp) ::   zfac12, zua, zva   ! temporary scalars
       REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz            ! temporary 2D workspace
       REAL(wp), DIMENSION(jpi,jpj) ::   ztnw, ztne, ztsw, ztse   ! temporary 3D workspace
+#if defined key_vvl
+      REAL(wp), DIMENSION(jpi,jpj,jpk)       ::   ze3f
+#else
       REAL(wp), DIMENSION(jpi,jpj,jpk), SAVE ::   ze3f
+#endif
       !!----------------------------------------------------------------------
 
@@ -548,5 +548,5 @@
       ENDIF
 
-      IF( kt == nit000 .OR. lk_vvl ) THEN
+      IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t)
          DO jk = 1, jpk
             DO jj = 1, jpjm1
@@ -561,6 +561,5 @@
       ENDIF
 
-      ! Local constant initialization
-      zfac12 = 1.e0 / 12.e0
+      zfac12 = 1.e0 / 12.e0      ! Local constant initialization
 
       
@@ -573,6 +572,8 @@
          ! -----------------------------------------
          SELECT CASE( kvor )      ! vorticity considered
-         CASE ( 1 )   ;   zwz(:,:) = ff(:,:)      * ze3f(:,:,jk)   ! planetary vorticity (Coriolis)
-         CASE ( 2 )   ;   zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk)   ! relative  vorticity
+         CASE ( 1 )                                                ! planetary vorticity (Coriolis)
+            zwz(:,:) = ff(:,:)      * ze3f(:,:,jk)
+         CASE ( 2 )                                                ! relative  vorticity
+            zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk)
          CASE ( 3 )                                                ! metric term
             DO jj = 1, jpjm1
@@ -583,5 +584,6 @@
                END DO
             END DO
-         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) ! total (relative + planetary vorticity)
+         CASE ( 4 )                                                ! total (relative + planetary vorticity)
+            zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk)
          CASE ( 5 )                                                ! total (coriolis + metric)
             DO jj = 1, jpjm1
@@ -590,5 +592,5 @@
                        &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               &
+                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                &
                        &       ) * ze3f(ji,jj,jk)
                END DO
@@ -601,6 +603,6 @@
          ! Compute and add the vorticity term trend
          ! ----------------------------------------
-         jj=2
-         ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0
+         jj = 2
+         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
          DO ji = 2, jpi   
                ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
@@ -638,5 +640,5 @@
       !!
       !! ** Purpose :   Control the consistency between cpp options for
-      !!      tracer advection schemes
+      !!              tracer advection schemes
       !!----------------------------------------------------------------------
       INTEGER ::   ioptio          ! temporary integer
Index: branches/dev_004_VVL/NEMO/OPA_SRC/DYN/wzvmod.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/DYN/wzvmod.F90	(revision 1407)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/DYN/wzvmod.F90	(revision 1408)
@@ -1,18 +1,16 @@
 MODULE wzvmod
+!! MODULE sshwzv   
    !!==============================================================================
-   !!                       ***  MODULE  wzvmod  ***
-   !! Ocean diagnostic variable : vertical velocity
+   !!                       ***  MODULE  sshwzv  ***
+   !! Ocean dynamics : sea surface height and vertical velocity
    !!==============================================================================
-   !! History :   5.0  !  90-10  (C. Levy, G. Madec)  Original code
-   !!             7.0  !  96-01  (G. Madec)  Statement function for e3
-   !!             8.5  !  02-07  (G. Madec)  Free form, F90
-   !!              "   !  07-07  (D. Storkey) Zero zhdiv at open boundary (BDY) 
-   !!----------------------------------------------------------------------
-   !!   wzv            : empty routine 
+   !! History :  3.1  !  2009-02  (G. Madec, M. Leclair)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
    !!   ssh_wzv        : after ssh & now vertical velocity
    !!   ssh_nxt        : filter ans swap the ssh arrays
    !!   ssh_rst        : read/write ssh restart fields in the ocean restart file
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers variables
    USE dom_oce         ! ocean space and time domain variables 
@@ -24,5 +22,4 @@
    USE prtctl          ! Print control
    USE phycst
-   USE bdy_oce         ! unstructured open boundaries
    USE lbclnk          ! ocean lateral boundary condition (or mpp link)
    USE obc_par         ! open boundary cond. parameter
@@ -32,8 +29,6 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC wzv          ! routine called by step.F90 and inidtr.F90
-   PUBLIC ssh_wzv
-   PUBLIC ssh_nxt
+   PUBLIC   ssh_wzv    ! called by step.F90
+   PUBLIC   ssh_nxt    ! called by step.F90
 
    !! * Substitutions
@@ -42,5 +37,5 @@
 
    !!----------------------------------------------------------------------
-   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -48,27 +43,4 @@
 
 CONTAINS
-
-   SUBROUTINE wzv( kt )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE wzv  ***
-      !!
-      !! ** Purpose :   Compute the now vertical velocity after the array swap
-      !!
-      !! ** Method  :   Using the incompressibility hypothesis, the vertical
-      !!      velocity is computed by integrating the horizontal divergence 
-      !!      from the bottom to the surface.
-      !!        The boundary conditions are w=0 at the bottom (no flux) and,
-      !!      in regid-lid case, w=0 at the sea surface.
-      !!
-      !! ** action  :   wn array : the now vertical velocity
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) :: kt
-      
-      ! Empty routine
-
-      WRITE(*,*) 'wzv : you should not be here : error ?'  
-
-   END SUBROUTINE wzv
-
 
    SUBROUTINE ssh_wzv( kt ) 
@@ -90,5 +62,4 @@
       !! if lk_vvl=T:   sshu_a, sshv_a, sshf_a  : after sea surface height
       !!                          at u-, v-, f-point s
-      !!                .._1    : now vertical coordinate arrays
       !!                hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points
       !!----------------------------------------------------------------------
@@ -103,5 +74,5 @@
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'ssh_wzv : vertical velocity from continuity eq. (vvl option)'
+         IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity '
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
          !
@@ -152,5 +123,5 @@
 
       !                                                ! Sea surface elevation time stepping
-      ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * emp(:,:)                       + zhdiv(:,:) )  ) * tmask(:,:,1)
+      ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
 
 #if defined key_obc
@@ -158,6 +129,6 @@
       IF ( Agrif_Root() ) THEN 
 # endif
-         ssha(:,:) = ssha(:,:)*obctmsk(:,:)
-         CALL lbc_lnk(ssha,'T',1.)  ! absolutly compulsory !! (jmm)
+         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
+         CALL lbc_lnk( ssha, 'T', 1. )  ! absolutly compulsory !! (jmm)
 # if defined key_agrif
       ENDIF
@@ -169,12 +140,12 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! 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) )
-   !!gm bug used of fmask, even if thereafter multiplied by muf  which is correctly masked)
-               sshf_a(ji,jj) = 0.25 * fmask(ji,jj,1) * ( ssha(ji  ,jj) + ssha(ji  ,jj+1)                 &
+               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) )
+               sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1)   &   ! Caution : fmask not used
+                  &                                  * ( ssha(ji  ,jj) + ssha(ji  ,jj+1)                 &
                   &                                    + ssha(ji+1,jj) + ssha(ji+1,jj+1) )
             END DO
@@ -199,5 +170,5 @@
       !                                           !------------------------------!
       IF( lk_vvl ) THEN                           ! only in vvl case)
-      !                                                ! now local depth and scale factors (stored in fse3. arrays)
+         !                                             ! now local depth and scale factors (stored in fse3. arrays)
          DO jk = 1, jpkm1
             fsdept(:,:,jk) = fsdept_n(:,:,jk)               ! depths
@@ -217,11 +188,6 @@
          hv(:,:) = hv_0(:,:) + sshv_n(:,:)
          !                                             ! masked inverse of the ocean depth (at u- and v-points)
-         hur(:,:) = 1. / MAX( hu(:,:), fse3u_0(:,:,1) ) * umask(:,:,1)
-         hvr(:,:) = 1. / MAX( hv(:,:), fse3v_0(:,:,1) ) * vmask(:,:,1)
-!!gm to be corrected (the above case does not work properly with 1 ocean level only)
-!        hur(:,:) = 1. / MAX( hu(:,:), 1.e-15 ) * umask(:,:,1)
-!        hvr(:,:) = 1. / MAX( hv(:,:), 1.e-15 ) * vmask(:,:,1)
-!!gm
-!!add end
+         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) )
+         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) )
 
       ENDIF
@@ -247,5 +213,5 @@
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
       !!
-      INTEGER  ::   ji, jj, jk           ! dummy loop indices
+      INTEGER  ::   ji, jj               ! dummy loop indices
       !!----------------------------------------------------------------------
 
@@ -330,6 +296,6 @@
          ENDIF
       ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
-         CALL iom_rstput( kt, nitrst, numrow, 'sshb'  , sshb  (:,:) )
-         CALL iom_rstput( kt, nitrst, numrow, 'sshn'  , sshn  (:,:) )
+         CALL iom_rstput( kt, nitrst, numrow, 'sshb'  , sshb(:,:) )
+         CALL iom_rstput( kt, nitrst, numrow, 'sshn'  , sshn(:,:) )
       ENDIF
       !
Index: branches/dev_004_VVL/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 1407)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 1408)
@@ -14,8 +14,11 @@
    !!            2.0  !  2006-02  (L. Debreu, C. Mazauric) Agrif implementation
    !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf
+   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
    !!   tra_nxt       : time stepping on temperature and salinity
+   !!   tra_nxt_fix   : time stepping on temperature and salinity : fixed    volume case
+   !!   tra_nxt_vvl   : time stepping on temperature and salinity : variable volume case
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers variables
@@ -127,6 +130,6 @@
 #endif      
 
-      ! trends diagnostics :  Asselin filter trend : (tb filtered - tb)/2dt
-      IF( l_trdtra ) THEN     
+      ! trends computation
+      IF( l_trdtra ) THEN              ! trend of the Asselin filter (tb filtered - tb)/dt     
          DO jk = 1, jpkm1
             zfact = 1.e0 / r2dt_t(jk)             
@@ -136,4 +139,5 @@
          CALL trd_mod( ztrdt, ztrds, jptra_trd_atf, 'TRA', kt )
       END IF
+
       !                        ! control print
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=tn, clinfo1=' nxt  - Tn: ', mask1=tmask,   &
@@ -182,16 +186,10 @@
          !                                           ! ----------------------- !
          !
-         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! case of Euler time-stepping at first time-step
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     ztm = 0.25 * ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) )      ! mean t
-                     zsm = 0.25 * ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) )
-                     tb(ji,jj,jk) = tn(ji,jj,jk)                                           ! tb <-- tn
-                     sb(ji,jj,jk) = sn(ji,jj,jk)
+         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
+            DO jk = 1, jpkm1                              ! (only swap)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
                      tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tb <-- tn
                      sn(ji,jj,jk) = sa(ji,jj,jk)
-                     ta(ji,jj,jk) = ztm                                                    ! ta <-- mean t
-                     sa(ji,jj,jk) = zsm
                   END DO
                END DO
@@ -219,10 +217,8 @@
          !                                           ! ----------------------- !
          !
-         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! case of Euler time-stepping at first time-step
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     tb(ji,jj,jk) = tn(ji,jj,jk)                                           ! tb <-- tn
-                     sb(ji,jj,jk) = sn(ji,jj,jk)
+         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
+            DO jk = 1, jpkm1
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
                      tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta
                      sn(ji,jj,jk) = sa(ji,jj,jk)
@@ -234,11 +230,8 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-!RBvvl for reproducibility
-!                     ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )       ! Asselin filter on t 
-!                     zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) )
-!                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                     ! tb <-- filtered tn 
-!                     sb(ji,jj,jk) = sn(ji,jj,jk) + zsf
-                     tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk)
-                     sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk)
+                     ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )       ! Asselin filter on t 
+                     zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) )
+                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                     ! tb <-- filtered tn 
+                     sb(ji,jj,jk) = sn(ji,jj,jk) + zsf
                      tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta
                      sn(ji,jj,jk) = sa(ji,jj,jk)
@@ -295,29 +288,10 @@
          !                                           ! ----------------------- !
          !
-         IF( neuler == 0 .AND. kt == nit000 ) THEN
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     ze3t_b = fse3t_b(ji,jj,jk)
-                     ze3t_n = fse3t_n(ji,jj,jk)
-                     ze3t_a = fse3t_a(ji,jj,jk)
-                     !                                         ! tracer content at Before, now and after
-                     ztcb = tb(ji,jj,jk) *  ze3t_b   ;   zscb = sb(ji,jj,jk) * ze3t_b
-                     ztcn = tn(ji,jj,jk) *  ze3t_n   ;   zscn = sn(ji,jj,jk) * ze3t_n
-                     ztca = ta(ji,jj,jk) *  ze3t_a   ;   zsca = sa(ji,jj,jk) * ze3t_a
-                     !
-                     !                                         ! mean thickness and tracer
-                     ze3mr= 1.e0  / ( ze3t_a + 2.* ze3t_n + ze3t_b )
-                     ztm  = ze3mr * ( ztca   + 2.* ztcn   + ztcb   )
-                     zsm  = ze3mr * ( zsca   + 2.* zscn   + zscb   )
-!!gm mean e3t have to be saved and used in dynhpg  or it can be recomputed in dynhpg !!
-!!gm                 e3t_m(ji,jj,jk) = 0.25 / ze3mr
-                     !                                         ! swap of arrays
-                     tb(ji,jj,jk) = tn(ji,jj,jk)                    ! tb <-- tn
-                     sb(ji,jj,jk) = sn(ji,jj,jk)
+         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
+            DO jk = 1, jpkm1                              ! (only swap)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
                      tn(ji,jj,jk) = ta(ji,jj,jk)                    ! tn <-- ta
                      sn(ji,jj,jk) = sa(ji,jj,jk)
-                     ta(ji,jj,jk) = ztm                             ! ta <-- mean t
-                     sa(ji,jj,jk) = zsm
                   END DO
                END DO
@@ -369,6 +343,4 @@
                DO jj = 1, jpj                             ! ONLY swap
                   DO ji = 1, jpi
-                     tb(ji,jj,jk) = tn(ji,jj,jk)                                 ! tb <-- tn
-                     sb(ji,jj,jk) = sn(ji,jj,jk)
                      tn(ji,jj,jk) = ta(ji,jj,jk)                                 ! tn <-- ta
                      sn(ji,jj,jk) = sa(ji,jj,jk)
@@ -394,7 +366,6 @@
                      zsc_f  = atfp * ( zsca   - 2.* zscn   + zscb   ) 
                      !
-!!gm tmask useless below
                      !                                         ! filtered tracer including the correction 
-                     ze3fr = tmask(ji,jj,jk) / ( ze3t_n + ze3t_f )
+                     ze3fr = 1.e0 / ( ze3t_n + ze3t_f )
                      ztf   =  ( ztcn  + ztc_f ) * ze3fr
                      zsf   =  ( zscn  + zsc_f ) * ze3fr
Index: branches/dev_004_VVL/NEMO/OPA_SRC/TRA/trazdf_exp.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/TRA/trazdf_exp.F90	(revision 1407)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/TRA/trazdf_exp.F90	(revision 1408)
@@ -42,5 +42,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA  3.0 , LOCEAN-IPSL (2008) 
+   !! NEMO/OPA  3.2 , LOCEAN-IPSL (2009) 
    !! $Id$ 
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -78,5 +78,5 @@
       INTEGER  ::   ji, jj, jk, jl            ! dummy loop indices
       REAL(wp) ::   zlavmr, zave3r, ze3tr     ! temporary scalars
-      REAL(wp) ::   zta, zsa, ze3tb, zcoef    ! temporary scalars
+      REAL(wp) ::   zta, zsa, ze3tb           ! temporary scalars
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zwy, zwz, zww   ! 3D workspace
       !!---------------------------------------------------------------------
@@ -130,11 +130,9 @@
             DO jj = 2, jpjm1 
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ze3tb = fse3t_b(ji,jj,jk)                                            ! before e3t
+                  ze3tb = fse3t_b(ji,jj,jk) / fse3t(ji,jj,jk)                          ! before e3t
                   zta   = zwx(ji,jj,jk) - tb(ji,jj,jk) + p2dt(jk) * ta(ji,jj,jk)       ! total trends * 2*rdt
                   zsa   = zwz(ji,jj,jk) - sb(ji,jj,jk) + p2dt(jk) * sa(ji,jj,jk)     
-                  zcoef = 1.e0 / fse3t_n(ji,jj,jk) * tmask(ji,jj,jk)
-!RBvvl : which e3t ?
-                  ta(ji,jj,jk) = (  ze3tb * tb(ji,jj,jk) + fse3t(ji,jj,jk) * zta  ) * zcoef
-                  sa(ji,jj,jk) = (  ze3tb * sb(ji,jj,jk) + fse3t(ji,jj,jk) * zsa  ) * zcoef
+                  ta(ji,jj,jk) = (  ze3tb * tb(ji,jj,jk) + zta  ) * tmask(ji,jj,jk)
+                  sa(ji,jj,jk) = (  ze3tb * sb(ji,jj,jk) + zsa  ) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -144,6 +142,6 @@
             DO jj = 2, jpjm1 
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ta(ji,jj,jk) = ( zwx(ji,jj,jk) + p2dt(jk) * ta(ji,jj,jk) ) *tmask(ji,jj,jk)
-                  sa(ji,jj,jk) = ( zwz(ji,jj,jk) + p2dt(jk) * sa(ji,jj,jk) ) *tmask(ji,jj,jk)
+                  ta(ji,jj,jk) = ( zwx(ji,jj,jk) + p2dt(jk) * ta(ji,jj,jk) ) * tmask(ji,jj,jk)
+                  sa(ji,jj,jk) = ( zwz(ji,jj,jk) + p2dt(jk) * sa(ji,jj,jk) ) * tmask(ji,jj,jk)
                END DO
             END DO
Index: branches/dev_004_VVL/NEMO/OPA_SRC/TRA/trazdf_imp.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 1407)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 1408)
@@ -1,22 +1,23 @@
 MODULE trazdf_imp
-   !!==============================================================================
+   !!======================================================================
    !!                 ***  MODULE  trazdf_imp  ***
    !! Ocean active tracers:  vertical component of the tracer mixing trend
-   !!==============================================================================
-   !! History :
-   !!   6.0  !  90-10  (B. Blanke)  Original code
-   !!   7.0  !  91-11 (G. Madec)
-   !!        !  92-06 (M. Imbard) correction on tracer trend loops
-   !!        !  96-01 (G. Madec) statement function for e3
-   !!        !  97-05 (G. Madec) vertical component of isopycnal
-   !!        !  97-07 (G. Madec) geopotential diffusion in s-coord
-   !!        !  00-08 (G. Madec) double diffusive mixing
-   !!   8.5  !  02-08 (G. Madec)  F90: Free form and module
-   !!   9.0  !  06-11 (G. Madec)  New step reorganisation
+   !!======================================================================
+   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
+   !!            7.0  !  1991-11  (G. Madec)
+   !!                 !  1992-06  (M. Imbard) correction on tracer trend loops
+   !!                 !  1996-01  (G. Madec) statement function for e3
+   !!                 !  1997-05  (G. Madec) vertical component of isopycnal
+   !!                 !  1997-07  (G. Madec) geopotential diffusion in s-coord
+   !!                 !  2000-08  (G. Madec) double diffusive mixing
+   !!   NEMO     1.0  !  2002-08  (G. Madec) F90: Free form and module
+   !!            2.0  !  2006-11  (G. Madec) New step reorganisation
+   !!            3.2  !  2009-03  (G. Madec)  heat and salt content trends
+   !!----------------------------------------------------------------------
+  
    !!----------------------------------------------------------------------
    !!   tra_zdf_imp : Update the tracer trend with the diagonal vertical  
    !!                 part of the mixing tensor.
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers variables
    USE dom_oce         ! ocean space and time domain variables 
@@ -36,6 +37,5 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC tra_zdf_imp   !  routine called by step.F90
+   PUBLIC   tra_zdf_imp   !  routine called by step.F90
 
    !! * Substitutions
@@ -45,6 +45,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
-   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -90,18 +89,14 @@
       !!
       !!---------------------------------------------------------------------
-      !! * Modules used
-      USE oce    , ONLY :   zwd   => ua,  &  ! ua used as workspace
-                            zws   => va      ! va   "          "
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
-      REAL(wp), DIMENSION(jpk), INTENT( in ) ::   &
-         p2dt                                ! vertical profile of tracer time-step
-
-      !! * Local declarations
-      INTEGER  ::   ji, jj, jk               ! dummy loop indices
-      REAL(wp) ::   zavi, zrhs, znvvl,     & ! temporary scalars
-         ze3tb, ze3tn, ze3ta                 ! variable vertical scale factors
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-         zwi, zwt, zavsi                     ! workspace arrays
+      USE oce    , ONLY :   zwd   => ua   ! ua used as workspace
+      USE oce    , ONLY :   zws   => va   ! va  -          -
+      !!
+      INTEGER                 , INTENT(in) ::   kt     ! ocean time-step index
+      REAL(wp), DIMENSION(jpk), INTENT(in) ::   p2dt   ! vertical profile of tracer time-step
+      !!
+      INTEGER  ::   ji, jj, jk            ! dummy loop indices
+      REAL(wp) ::   zavi, zrhs, znvvl     ! temporary scalars
+      REAL(wp) ::   ze3tb, ze3tn, ze3ta   ! variable vertical scale factors
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi, zwt, zavsi   ! workspace arrays
       !!---------------------------------------------------------------------
 
@@ -169,6 +164,8 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ze3ta = ( 1. - znvvl ) + znvvl*fse3t_a(ji,jj,jk)                        ! after scale factor at T-point
-               ze3tn = ( 1. - znvvl )*fse3t_n(ji,jj,jk) + znvvl                        ! now   scale factor at T-point
+               ze3ta =  ( 1. - znvvl )   &                            ! after scale factor at T-point
+                  &   +        znvvl    * fse3t_a(ji,jj,jk) 
+               ze3tn =    znvvl          &                            ! now   scale factor at T-point
+                  &   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)
                zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
                zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
@@ -181,5 +178,6 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            ze3ta = ( 1. - znvvl ) + znvvl*fse3t_a(ji,jj,1)                             ! after scale factor at T-point
+            ze3ta = ( 1. - znvvl )        &                           ! after scale factor at T-point
+               &   +       znvvl      * fse3t_a(ji,jj,1) 
             zwi(ji,jj,1) = 0.e0
             zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
@@ -225,7 +223,6 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1
-            ze3tb = ( 1. - znvvl ) + znvvl*fse3t_b(ji,jj,1)
-!RBvvl which 
-            ze3tn = ( 1. - znvvl ) + znvvl*fse3t(ji,jj,1)
+            ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1)
+            ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1)
             ta(ji,jj,1) = ze3tb * tb(ji,jj,1) + p2dt(1) * ze3tn * ta(ji,jj,1)
          END DO
@@ -234,9 +231,8 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
-               ze3tb = ( 1. - znvvl ) + znvvl*fse3t_b(ji,jj,jk)
-!RB_vvl
-               ze3tn = ( 1. - znvvl ) + znvvl*fse3t (ji,jj,jk)
+               ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk)
+               ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk)
                zrhs = ze3tb * tb(ji,jj,jk) + p2dt(jk) * ze3tn * ta(ji,jj,jk)   ! zrhs=right hand side 
-               ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *ta(ji,jj,jk-1)
+               ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * ta(ji,jj,jk-1)
             END DO
          END DO
@@ -269,6 +265,8 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ze3ta = ( 1. - znvvl ) + znvvl*fse3t_a(ji,jj,jk)                            ! after scale factor at T-point
-               ze3tn = ( 1. - znvvl )*fse3t_n(ji,jj,jk) + znvvl                            ! now   scale factor at T-point
+               ze3ta =  ( 1. - znvvl )                        &         ! after scale factor at T-point
+                  &   +        znvvl   * fse3t_a(ji,jj,jk)           
+               ze3tn =    znvvl                               &         ! now   scale factor at T-point
+                  &   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)
                zwi(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
                zws(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
@@ -281,5 +279,5 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            ze3ta = ( 1. - znvvl ) + znvvl*fse3t_a(ji,jj,1)
+            ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,1)
             zwi(ji,jj,1) = 0.e0
             zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
@@ -324,6 +322,7 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1
-            ze3tb = ( 1. - znvvl ) + znvvl*fse3t_b(ji,jj,1)                      ! before scale factor at T-point
-            ze3tn = ( 1. - znvvl ) + znvvl*fse3t(ji,jj,1)                        ! now    scale factor at T-point
+            ze3tb = ( 1. - znvvl )   &                                           ! before scale factor at T-point
+               &   +  znvvl       * fse3t_b(ji,jj,1)
+            ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,1)                    ! now    scale factor at T-point
             sa(ji,jj,1) = ze3tb * sb(ji,jj,1) + p2dt(1) * ze3tn * sa(ji,jj,1)
          END DO
@@ -332,6 +331,7 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
-               ze3tb = ( 1. - znvvl ) + znvvl*fse3t_b(ji,jj,jk)                  ! before scale factor at T-point
-               ze3tn = ( 1. - znvvl ) + znvvl*fse3t(ji,jj,jk)                    ! now    scale factor at T-point
+               ze3tb = ( 1. - znvvl )   &                                        ! before scale factor at T-point
+                  &   +  znvvl       * fse3t_b(ji,jj,jk)
+               ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk)                ! now    scale factor at T-point
                zrhs = ze3tb * sb(ji,jj,jk) + p2dt(jk) * ze3tn * sa(ji,jj,jk)     ! zrhs=right hand side
                sa(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *sa(ji,jj,jk-1)
@@ -355,5 +355,5 @@
          END DO
       END DO
-
+      !
    END SUBROUTINE tra_zdf_imp
 
Index: branches/dev_004_VVL/NEMO/OPA_SRC/ZDF/zdfevd.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/ZDF/zdfevd.F90	(revision 1407)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/ZDF/zdfevd.F90	(revision 1408)
@@ -5,11 +5,14 @@
    !!                of vertical eddy mixing coefficient
    !!======================================================================
+   !! History :  OPA  !  1997-06  (G. Madec, A. Lazar)  Original code
+   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
+   !!             -   !  2005-06  (C. Ethe) KPP parameterization
+   !!            3.2  !  2009-03  (M. Leclair, G. Madec, R. Benshila) test on both before & after
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   zdf_evd      : update momentum and tracer Kz at the location of
-   !!                  statically unstable portion of the water column
-   !!                  (called if ln_zdfevd=T)
+   !!   zdf_evd      : increase the momentum and tracer Kz at the location of
+   !!                  statically unstable portion of the water column (ln_zdfevd=T)
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers variables
    USE dom_oce         ! ocean space and time domain variables
@@ -22,13 +25,12 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC zdf_evd      ! called by step.F90
+   PUBLIC   zdf_evd    ! called by step.F90
 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
-   !! $Id$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
+   !! $Id:$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    !!----------------------------------------------------------------------
 
@@ -48,15 +50,8 @@
       !! ** Action  :   Update avt, avmu, avmv in statically instable cases
       !!                and avt_evd which is avt due to convection
-      !! References :
-      !!      Lazar, A., these de l'universite Paris VI, France, 1997
-      !! History :
-      !!   7.0  !  97-06  (G. Madec, A. Lazar)  Original code
-      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
-      !!   9.0  !  05-06  (C. Ethe) KPP parameterization
+      !! References :   Lazar, A., these de l'universite Paris VI, France, 1997
       !!----------------------------------------------------------------------
-      !! * Arguments
       INTEGER, INTENT( in ) ::   kt         ! ocean time-step indexocean time step
-
-      !! * Local declarations
+      !!
       INTEGER ::   ji, jj, jk               ! dummy loop indices
       !!----------------------------------------------------------------------
@@ -80,29 +75,16 @@
             !                                             ! ===============
 #if defined key_vectopt_loop
-!!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL!
-            jj = 1                     ! big loop forced
-            DO ji = jpi+2, jpij   
-# if defined key_zdfkpp
-!! no implicit mixing in the boundary layer with KPP
-               IF( ( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) .AND. ( fsdepw(ji,jj,jk) > hkpp(ji,jj) ) ) THEN
-# else
-               IF(   MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN
-# endif
-                  avt (ji  ,jj  ,jk) = avevd * tmask(ji  ,jj  ,jk)
-                  avmu(ji  ,jj  ,jk) = avevd * umask(ji  ,jj  ,jk)
-                  avmu(ji-1,jj  ,jk) = avevd * umask(ji-1,jj  ,jk)
-                  avmv(ji  ,jj  ,jk) = avevd * vmask(ji  ,jj  ,jk)
-                  avmv(ji  ,jj-1,jk) = avevd * vmask(ji  ,jj-1,jk)
-               ENDIF
-            END DO
+            DO jj = 1, 1                     ! big loop forced
+               DO ji = jpi+2, jpij   
 #else
             DO jj = 2, jpj             ! no vector opt.
                DO ji = 2, jpi
-# if defined key_zdfkpp
+#endif
+#if defined key_zdfkpp
 !! no implicit mixing in the boundary layer with KPP
-               IF( ( MIN( rn2(ji,jj,jk),  rn2b(ji,jj,jk) ) <= -1.e-12 ) .AND. ( fsdepw(ji,jj,jk) > hkpp(ji,jj) ) ) THEN
-# else
-               IF(   MIN( rn2(ji,jj,jk),  rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN
-# endif
+                  IF( ( MIN( rn2(ji,jj,jk),  rn2b(ji,jj,jk) ) <= -1.e-12 ) .AND. ( fsdepw(ji,jj,jk) > hkpp(ji,jj) ) ) THEN
+#else
+                  IF(   MIN( rn2(ji,jj,jk),  rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN
+#endif
                      avt (ji  ,jj  ,jk) = avevd * tmask(ji  ,jj  ,jk)
                      avmu(ji  ,jj  ,jk) = avevd * umask(ji  ,jj  ,jk)
@@ -113,5 +95,4 @@
                END DO
             END DO
-#endif
             !                                             ! ===============
          END DO                                           !   End of slab
@@ -130,27 +111,19 @@
 !!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL! 
 #if defined key_vectopt_loop
-            jj = 1                     ! big loop forced
-            DO ji = 1, jpij   
-# if defined key_zdfkpp
-!! no implicit mixing in the boundary layer with KPP
-               IF( ( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) .AND. ( fsdepw(ji,jj,jk) > hkpp(ji,jj) ) ) &              
-# else
-               IF(   MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) &
-# endif
-                  avt(ji,jj,jk) = avevd * tmask(ji,jj,jk)
-            END DO
+            DO jj = 1, 1                     ! big loop forced
+               DO ji = 1, jpij   
 #else
             DO jj = 1, jpj             ! loop over the whole domain (no lbc_lnk call)
                DO ji = 1, jpi
-# if defined key_zdfkpp
+#endif
+#if defined key_zdfkpp
 !! no implicit mixing in the boundary layer with KPP
                   IF( ( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) .AND. ( fsdepw(ji,jj,jk) > hkpp(ji,jj) ) ) &          
-# else
+#else
                   IF(   MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) &
-# endif
+#endif
                      avt(ji,jj,jk) = avevd * tmask(ji,jj,jk)
                END DO
             END DO
-#endif
             !                                             ! ===============
          END DO                                           !   End of slab
@@ -159,6 +132,6 @@
 
       ! update of avt_evd and avmu_evd
-      avt_evd  (:,:,:) = avt (:,:,:)  - avt_evd  (:,:,:) 
-      avmu_evd (:,:,:) = avmu(:,:,:)  - avmu_evd (:,:,:) 
+      avt_evd (:,:,:) = avt (:,:,:)  - avt_evd (:,:,:) 
+      avmu_evd(:,:,:) = avmu(:,:,:)  - avmu_evd(:,:,:) 
 
    END SUBROUTINE zdf_evd
Index: branches/dev_004_VVL/NEMO/OPA_SRC/oce.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/oce.F90	(revision 1407)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/oce.F90	(revision 1408)
@@ -4,13 +4,8 @@
    !! Ocean        :  dynamics and active tracers defined in memory 
    !!======================================================================
-   !! History :
-   !!   8.5  !  02-11  (G. Madec)  F90: Free form and module
-   !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
+   !! History :  0.1  !  2002-11  (G. Madec)  F90: Free form and module
+   !!            1.0  !  2005-11  (V. Garnier) Surface pressure gradient organization
+   !!            3.1  !  2009-02  (G. Madec, M. Leclair)  pure z* coordinate
    !!----------------------------------------------------------------------
-   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
-   !! $Id$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
-   !!----------------------------------------------------------------------
-   !! * Modules used
    USE par_oce      ! ocean parameters
 
@@ -20,44 +15,34 @@
    !! Physics and algorithm flags
    !! ---------------------------
-   LOGICAL, PUBLIC ::   l_traldf_rot    = .FALSE.  !: rotated laplacian operator for lateral diffusion 
+   LOGICAL, PUBLIC ::   l_traldf_rot    = .FALSE.  !: rotated laplacian operator for lateral diffusion
    LOGICAL, PUBLIC ::   ln_dynhpg_imp   = .FALSE.  !: semi-implicite hpg flag
    INTEGER, PUBLIC ::   nn_dynhpg_rst   = 0        !: add dynhpg implicit variables in restart ot not
 
-   !! dynamics and tracer fields
-   !! --------------------------
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
-      ! before !  now      !  after  !      ! the after trends becomes the fields
-      ! fields !  fields   !  trends !      ! only in dyn(tra)_zdf and dyn(tra)_nxt
-      ub       ,  un       ,  ua     ,   &  !: i-horizontal velocity (m/s)
-      vb       ,  vn       ,  va     ,   &  !: j-horizontal velocity (m/s)
-                  wn       ,             &  !: vertical velocity (m/s)
-      rotb     ,  rotn     ,             &  !: relative vorticity (1/s)
-      hdivb    ,  hdivn    ,             &  !: horizontal divergence (1/s)
-      tb       ,  tn       ,  ta     ,   &  !: potential temperature (celcius)
-      sb       ,  sn       ,  sa            !: salinity (psu)
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
-      rhd ,                              &  !: in situ density anomalie rhd=(rho-rau0)/rau0 (no units)
-      rhop,                              &  !: potential volumic mass (kg/m3)
-      rn2,                               &  !: now    brunt-vaisala frequency (1/s2)
-      rn2b                                  !: before brunt-vaisala frequency (1/s2)
+   !! dynamics and tracer fields             !  before  !  now     !  after   ! the after trends becomes the fields
+   !! --------------------------             !  fields  !  fields  !  trends  ! only after tra_zdf and dyn_spg
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   ub     ,  un      ,  ua      !: i-horizontal velocity      [m/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   vb     ,  vn      ,  va      !: j-horizontal velocity      [m/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::             wn                 !: vertical velocity          [m/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   rotb   ,  rotn               !: relative vorticity         [s-1]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   hdivb  ,  hdivn              !: horizontal divergence      [s-1]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   tb     ,  tn      ,  ta      !: potential temperature      [Celcius]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   sb     ,  sn      ,  sa      !: salinity                   [psu]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   rn2b   ,  rn2                !: brunt-vaisala frequency**2 [s-2]
+   !
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   rhd    !: in situ density anomalie rhd=(rho-rau0)/rau0     [no units]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   rhop   !: potential volumic mass                           [kg/m3]
 
-
-      !! advection scheme choice
-      !! -----------------------
-      CHARACTER(len=3), PUBLIC  ::   l_adv   !: 'ce2' centre scheme used
-         !                                   !: 'tvd' TVD scheme used
-         !                                   !: 'mus' MUSCL scheme used
-         !                                   !: 'mu2' MUSCL2 scheme used
+   !! advection scheme choice
+   !! -----------------------
+   CHARACTER(len=3), PUBLIC  ::   l_adv   !: flag for the advection scheme used (= 'ce2', 'tvd', 'mus' or ...)
 
    !! surface pressure gradient
    !! -------------------------
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !:
-      spgu, spgv             !: horizontal surface pressure gradient
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   spgu, spgv      !: horizontal surface pressure gradient
 
    !! interpolated gradient (only used in zps case)
    !! ---------------------
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !:
-      gtu, gsu, gru,      &  !: t-, s- and rd horizontal gradient at u- and 
-      gtv, gsv, grv          !: v-points at bottom ocean level 
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   gtu, gsu, gru   !: horizontal gradient of T, S and rd at bottom u-point
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   gtv, gsv, grv   !: horizontal gradient of T, S and rd at bottom v-point 
 
    !! free surface                       !  before  !  now     !  after   !
@@ -68,14 +53,14 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshf_b ,  sshf_n  ,  sshf_a  !: sea surface height at f-point [m]
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshbb         !: before before & after sea surface height at t-point
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshn_f        !: filtered now sea surface height at t-point
-
 #if defined key_dynspg_rl   ||   defined key_esopa
    !! rigid-lid formulation
    !! ---------------------
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !:
-      bsfb, bsfn,         &  !: before, now barotropic streamfunction (m3/s)
-      bsfd                   !: now trend of barotropic streamfunction (m3/s2)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   bsfb, bsfn   !: before, now barotropic streamfunction (m3/s)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   bsfd         !: now trend of barotropic streamfunction (m3/s2)
 #endif
    !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2008) 
+   !! $Id:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!======================================================================
 END MODULE oce
Index: branches/dev_004_VVL/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/step.F90	(revision 1407)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/step.F90	(revision 1408)
@@ -4,21 +4,20 @@
    !! Time-stepping    : manager of the ocean, tracer and ice time stepping
    !!======================================================================
-   !! History :        !  91-03  (G. Madec)  Original code
-   !!                  !  92-06  (M. Imbard)  add a first output record
-   !!                  !  96-04  (G. Madec)  introduction of dynspg
-   !!                  !  96-04  (M.A. Foujols)  introduction of passive tracer
-   !!             8.0  !  97-06  (G. Madec)  new architecture of call
-   !!             8.2  !  97-06  (G. Madec, M. Imbard, G. Roullet)  free surface
-   !!             8.2  !  99-02  (G. Madec, N. Grima)  hpg implicit
-   !!             8.2  !  00-07  (J-M Molines, M. Imbard)  Open Bondary Conditions
-   !!             9.0  !  02-06  (G. Madec)  free form, suppress macro-tasking
-   !!             " "  !  04-08  (C. Talandier) New trends organization
-   !!             " "  !  05-01  (C. Ethe) Add the KPP closure scheme
-   !!             " "  !  05-11  (V. Garnier) Surface pressure gradient organization
-   !!             " "  !  05-11  (G. Madec)  Reorganisation of tra and dyn calls
-   !!             " "  !  06-01  (L. Debreu, C. Mazauric)  Agrif implementation
-   !!             " "  !  06-07  (S. Masson)  restart using iom
-   !!             " "  !  06-08  (G. Madec)  surface module 
-   !!             " "  !  07-07  (J. Chanut, A. Sellar) Unstructured open boundaries (BDY)
+   !! History :  OPA  !  1991-03  (G. Madec)  Original code
+   !!                 !  1991-11  (G. Madec)
+   !!                 !  1992-06  (M. Imbard)  add a first output record
+   !!                 !  1996-04  (G. Madec)  introduction of dynspg
+   !!                 !  1996-04  (M.A. Foujols)  introduction of passive tracer
+   !!            8.0  !  1997-06  (G. Madec)  new architecture of call
+   !!            8.2  !  1997-06  (G. Madec, M. Imbard, G. Roullet)  free surface
+   !!             -   !  1999-02  (G. Madec, N. Grima)  hpg implicit
+   !!             -   !  2000-07  (J-M Molines, M. Imbard)  Open Bondary Conditions
+   !!   NEMO     1.0  !  2002-06  (G. Madec)  free form, suppress macro-tasking
+   !!             -   !  2004-08  (C. Talandier) New trends organization
+   !!                 !  2005-01  (C. Ethe) Add the KPP closure scheme
+   !!                 !  2005-11  (G. Madec)  Reorganisation of tra and dyn calls
+   !!                 !  2006-01  (L. Debreu, C. Mazauric)  Agrif implementation
+   !!                 !  2006-07  (S. Masson)  restart using iom
+   !!            3.2  !  2009-02  (G. Madec, R. Benshila)  reintroduicing z*-coordinate
    !!----------------------------------------------------------------------
 
@@ -116,5 +115,4 @@
    USE restart         ! ocean restart                    (rst_wri routine)
    USE prtctl          ! Print control                    (prt_ctl routine)
-   USE domvvl          ! variable volume                  (dom_vvl routine)
 
 #if defined key_agrif
@@ -140,6 +138,8 @@
 #if defined key_agrif
    SUBROUTINE stp( )
+      INTEGER             ::   kstp   ! ocean time-step index
 #else
    SUBROUTINE stp( kstp )
+      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
 #endif
       !!----------------------------------------------------------------------
@@ -159,10 +159,4 @@
       !!              -8- Outputs and diagnostics
       !!----------------------------------------------------------------------
-      !! * Arguments
-#if defined key_agrif   
-      INTEGER             ::   kstp   ! ocean time-step index
-#else
-      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
-#endif      
       INTEGER ::   jk       ! dummy loop indice
       INTEGER ::   indic    ! error indicator if < 0
@@ -202,8 +196,6 @@
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Computation of diagnostic variables
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
-      !-----------------------------------------------------------------------
+      !                Ocean dynamics : ssh, wn, hdiv, rot                   !
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
                        CALL div_cur( kstp )                 ! Horizontal divergence & Relative vorticity
       IF( n_cla == 1 ) CALL div_cla( kstp )                 ! Cross Land Advection (Update Hor. divergence)
@@ -213,10 +205,10 @@
       ! Ocean physics update
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
+      !-----------------------------------------------------------------------
 #if defined key_zdftke2
       IF ( ln_dynhpg_imp ) THEN
       !-----------------------------------------------------------------------
       !  LATERAL PHYSICS
-      !-----------------------------------------------------------------------
-      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
       !-----------------------------------------------------------------------
                                CALL zdf_mxl( kstp )                 ! mixed layer depth
@@ -229,6 +221,4 @@
       !-----------------------------------------------------------------------
       !  VERTICAL PHYSICS
-      !-----------------------------------------------------------------------
-      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
       !-----------------------------------------------------------------------
       IF( neuler == 0 .AND. kstp == nit000 ) THEN
@@ -273,6 +263,4 @@
       !-----------------------------------------------------------------------
       !  LATERAL PHYSICS
-      !-----------------------------------------------------------------------
-      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
       !-----------------------------------------------------------------------
       IF( lk_ldfslp     )   CALL ldf_slp( kstp, rhd, rn2b )      ! before slope of the lateral mixing
