Index: branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_flt.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 1362)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 1368)
@@ -427,32 +427,4 @@
       END DO
 
-      ! Sea surface elevation time stepping
-      ! -----------------------------------
-      IF( lk_vvl ) THEN   ! after free surface elevation
-         zssha(:,:) = ssha(:,:)
-      ELSE
-         zssha(:,:) = sshb(:,:) + z2dt * ( wn(:,:,1) - emp(:,:) * zraur ) * tmask(:,:,1)
-      ENDIF
-#if defined key_obc
-# if defined key_agrif
-      IF ( Agrif_Root() ) THEN
-# endif
-         zssha(:,:)=zssha(:,:)*obctmsk(:,:)
-         CALL lbc_lnk(zssha,'T',1.)  ! absolutly compulsory !! (jmm)
-# if defined key_agrif
-      ENDIF
-# endif
-#endif
-
-      IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler (forward) time stepping, no time filter
-         ! swap of arrays
-         sshb(:,:) = sshn (:,:)
-         sshn(:,:) = zssha(:,:)
-      ELSE                                           ! Leap-frog time stepping and time filter
-         ! time filter and array swap
-         sshb(:,:) = atfp * ( sshb(:,:) + zssha(:,:) ) + atfp1 * sshn(:,:)
-         sshn(:,:) = zssha(:,:)
-      ENDIF
-
       ! write filtered free surface arrays in restart file
       ! --------------------------------------------------
Index: branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 1362)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 1368)
@@ -566,49 +566,5 @@
       ! Phase 3 : Update sea surface height from time averaged barotropic variables
       ! ---------------------------------------------------------------------------
-
- 
-      ! Horizontal divergence of time averaged barotropic transports
-      !-------------------------------------------------------------
-      IF( .NOT. lk_vvl ) THEN
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj  ) * un_b(ji-1,jj  )     &
-                  &            +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji  ,jj-1) * vn_b(ji  ,jj-1) )   &
-                  &           / ( e1t(ji,jj) * e2t(ji,jj) )
-            END DO
-         END DO
-      ENDIF
-
-#if defined key_obc && ! defined key_vvl
-      ! open boundaries (div must be zero behind the open boundary)
-      !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
-      IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east
-      IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west
-      IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
-      IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south
-#endif
-
-#if defined key_bdy
-         DO jj = 1, jpj
-           DO ji = 1, jpi
-             zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj)
-           END DO
-         END DO
-#endif
-
-      ! sea surface height
-      !-------------------
-      IF( lk_vvl ) THEN
-         sshbb(:,:) = sshb(:,:)
-         sshb (:,:) = sshn(:,:)
-         sshn (:,:) = ssha(:,:)
-      ELSE
-         sshb (:,:) = sshn(:,:)
-         sshn(:,:) = (  sshb_b(:,:) - z2dt_b * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
-      ENDIF
-
-      ! ... Boundary conditions on sshn
-      IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. )
-
+!RB_vvl now done in ssh_wzv and ssh_nxt
 
       ! -----------------------------------------------------------------------------
Index: branches/dev_004_VVL/NEMO/OPA_SRC/DYN/wzvmod.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/DYN/wzvmod.F90	(revision 1362)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/DYN/wzvmod.F90	(revision 1368)
@@ -29,4 +29,6 @@
    !! * Routine accessibility
    PUBLIC wzv          ! routine called by step.F90 and inidtr.F90
+   PUBLIC ssh_wzv
+   PUBLIC ssh_nxt
 
    !! * Substitutions
@@ -54,100 +56,144 @@
       !! ** action  :   wn array : the now vertical velocity
       !!----------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-
-      !! * Local declarations
-      INTEGER  ::           jk           ! dummy loop indices
-      !! Variable volume
-      INTEGER  ::   ji, jj               ! dummy loop indices
-      REAL(wp) ::   z2dt, zraur          ! temporary scalar
-      REAL(wp), DIMENSION (jpi,jpj) ::   zssha, zun, zvn, zhdiv
-#if defined key_bdy
-      INTEGER  ::     jgrd, jb           ! temporary scalars
-#endif
+      INTEGER, INTENT(in) :: kt
+      
+      ! Empty routine
+
+      WRITE(*,*) 'wzv : you should not be here : error ?'  
+
+   END SUBROUTINE wzv
+
+   SUBROUTINE ssh_wzv( kt ) 
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_ssh  ***
+      !!                   
+      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity
+      !!              Cand update the now vertical coordinate (lk_vvl=T).
+      !!
+      !! ** Method  : -  
+
+      !!              - Using the incompressibility hypothesis, the vertical
+      !!      velocity is computed by integrating the horizontal divergence 
+      !!      from the bottom to the surface minus the scale factor evolution.
+      !!        The boundary conditions are w=0 at the bottom (no flux) and.
+      !!
+      !! ** action  :   wn array : the now vertical velocity
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! time step
+      !!
+      INTEGER  ::   ji, jj, jk      ! dummy loop indices
+      REAL(wp) ::   z2dt, zraur     ! temporary scalars
+      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv       ! 2D workspace
       !!----------------------------------------------------------------------
 
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'wzv     : vertical velocity from continuity eq.'
-         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
-
-         ! bottom boundary condition: w=0 (set once for all)
-         wn(:,:,jpk) = 0.e0
-      ENDIF
-
-      IF( lk_vvl ) THEN                ! Variable volume
+         IF(lwp) WRITE(numout,*) 'ssh_wzv : vertical velocity from continuity eq. (vvl option)'
+         IF(lwp) WRITE(numout,*) '~~~~~~~ '
          !
-         z2dt = 2. * rdt                                       ! time step: leap-frog
-         IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt       ! time step: Euler if restart from rest
-         zraur  = 1. / rauw
-
-         ! Vertically integrated quantities
-         ! --------------------------------
-         zun(:,:) = 0.e0
-         zvn(:,:) = 0.e0
-         !
-         DO jk = 1, jpkm1             ! Vertically integrated transports (now)
-            zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un(:,:,jk)
-            zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
+         wn(:,:,jpk) = 0.e0      ! bottom boundary condition: w=0 (set once for all)
+      ENDIF
+
+      ! set time step size (Euler/Leapfrog)
+      z2dt = 2. * rdt 
+      IF( neuler == 0 .AND. kt == nit000 ) z2dt =rdt
+
+      zraur = 1. / rauw
+
+      !                                           !------------------------------!
+      !                                           !   After Sea Surface Height   !
+      !                                           !------------------------------!
+      zhdiv(:,:) = 0.e0
+      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
+        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
+      END DO
+
+      !                                                ! Sea surface elevation time stepping
+      ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * emp(:,:)                       + zhdiv(:,:) )  ) * tmask(:,:,1)
+
+      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
+      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
+!RB_vvl
+! Compute ssh at u,v, f points and update vertical coordinate
+! Currently done in dom_vvl
+      ENDIF
+
+      !                                           !------------------------------!
+      !                                           !     Now Vertical Velocity    !
+      !                                           !------------------------------!
+      !                                                ! integrate from the bottom the hor. divergence
+         DO jk = jpkm1, 1, -1
+            wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)   &
+                 &                    - (  fse3t_a(:,:,jk)                   &
+                 &                       - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt
          END DO
 
-         ! Horizontal divergence of barotropic transports
-         !--------------------------------------------------
-         zhdiv(:,:) = 0.e0
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1   ! vector opt.
-               zhdiv(ji,jj) = (  e2u(ji  ,jj  ) * zun(ji  ,jj  )     &
-                  &            - e2u(ji-1,jj  ) * zun(ji-1,jj  )     &
-                  &            + e1v(ji  ,jj  ) * zvn(ji  ,jj  )     &
-                  &            - e1v(ji  ,jj-1) * zvn(ji  ,jj-1) )   &
-                  &           / ( e1t(ji,jj) * e2t(ji,jj) )
+      !                                           !------------------------------!
+      !                                           !  Update Now Vertical coord.  !
+      !                                           !------------------------------!
+      IF( lk_vvl ) THEN                           ! only in vvl case)
+      !                                                ! now local depth and scale factors (stored in fse3. arrays)
+!RB_vvl to be done
+
+
+      ENDIF
+      !
+   END SUBROUTINE ssh_wzv
+
+
+   SUBROUTINE ssh_nxt( kt )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE ssh_nxt  ***
+      !!
+      !! ** Purpose :   achieve the sea surface  height time stepping by 
+      !!              applying Asselin time filter and swapping the arrays
+      !!              ssha  already computed in ssh_wzv  
+      !!
+      !! ** Method  : - apply Asselin time fiter to now ssh and swap :
+      !!             sshn = ssha + atfp * ( sshb -2 sshn + ssha )
+      !!             sshn = ssha
+      !!
+      !! ** action  : - sshb, sshn   : before & now sea surface height
+      !!                               ready for the next time step
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk           ! dummy loop indices
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
+         IF(lwp) WRITE(numout,*) '~~~~~~~ '
+      ENDIF
+
+      ! Time filter and swap of the ssh
+      ! -------------------------------
+
+!RB_vvl ssh at u, f,v point to be added
+
+     IF( neuler == 0 .AND. kt == nit000 ) THEN   ! Euler time stepping
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ! before <-- now
+               sshb  (ji,jj) = sshn(ji,jj) 
+               ! now <-- after
+               sshn  (ji,jj) = ssha  (ji,jj)
             END DO
          END DO
-
-#if defined key_obc && ( defined key_dynspg_exp || defined key_dynspg_ts )
-         ! open boundaries (div must be zero behind the open boundary)
-         !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
-         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east
-         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west
-         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
-         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south
-#endif
-
-#if defined key_bdy
-         jgrd=1 !: tracer grid.
-         DO jb = 1, nblenrim(jgrd)
-           ji = nbi(jb,jgrd)
-           jj = nbj(jb,jgrd)
-           zhdiv(ji,jj) = 0.e0
+     ELSE                                         ! Leap-frog time stepping
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ! before <-- now filtered
+               sshb  (ji,jj) = sshn(ji,jj)   + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) )    !&
+               ! now <-- after
+               sshn  (ji,jj) = ssha  (ji,jj)
+            END DO
          END DO
-#endif
-
-         CALL lbc_lnk( zhdiv, 'T', 1. )
-
-         ! Sea surface elevation time stepping
-         ! -----------------------------------
-         zssha(:,:) = sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) * tmask(:,:,1)
-
-         ! Vertical velocity computed from bottom
-         ! --------------------------------------
-         DO jk = jpkm1, 1, -1
-            wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) &
-              &        - ( zssha(:,:) - sshb(:,:) ) * fse3t_b(:,:,jk) * mut(:,:,jk) / z2dt
-         END DO
-
-      ELSE                             ! Fixed volume 
-
-         ! Vertical velocity computed from bottom
-         ! --------------------------------------
-         DO jk = jpkm1, 1, -1
-            wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk)
-         END DO
-      
-      ENDIF 
-
-      IF(ln_ctl)   CALL prt_ctl(tab3d_1=wn, clinfo1=' w**2 -   : ', mask1=wn)
-
-   END SUBROUTINE wzv
+      ENDIF
+
+      IF(ln_ctl)   CALL prt_ctl(tab2d_1=sshb    , clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
+      !
+   END SUBROUTINE ssh_nxt
 
    !!======================================================================
Index: branches/dev_004_VVL/NEMO/OPA_SRC/istate.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/istate.F90	(revision 1362)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/istate.F90	(revision 1368)
@@ -156,13 +156,6 @@
          CALL bn2( tb, sb, rn2 )    ! before Brunt Vaissala frequency
 
-         IF( .NOT. ln_rstart ) CALL wzv( nit000 ) 
-
       ENDIF
 
-      !                                       ! Vertical velocity
-      !                                       ! -----------------
-
-      IF( .NOT. lk_vvl )    CALL wzv( nit000 )                         ! from horizontal divergence
-      !
    END SUBROUTINE istate_init
 
Index: branches/dev_004_VVL/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/step.F90	(revision 1362)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/step.F90	(revision 1368)
@@ -201,4 +201,14 @@
       ENDIF
 
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Computation of diagnostic variables
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
+      !-----------------------------------------------------------------------
+                       CALL div_cur( kstp )                 ! Horizontal divergence & Relative vorticity
+      IF( n_cla == 1 ) CALL div_cla( kstp )                 ! Cross Land Advection (Update Hor. divergence)
+                       CALL ssh_wzv( kstp       )           ! after ssh & vertical velocity
+      IF( lk_vvl )     CALL dom_vvl                         ! vertical mesh at next time step
+
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -360,15 +370,6 @@
                                CALL dyn_spg( kstp, indic )    ! surface pressure gradient
                                CALL dyn_nxt( kstp )           ! lateral velocity at next time step
-      IF( lk_vvl )             CALL dom_vvl                   ! vertical mesh at next time step
-
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Computation of diagnostic variables
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
-      !-----------------------------------------------------------------------
-                       CALL div_cur( kstp )                 ! Horizontal divergence & Relative vorticity
-      IF( n_cla == 1 ) CALL div_cla( kstp )                 ! Cross Land Advection (Update Hor. divergence)
-                       CALL wzv( kstp )                     ! Vertical velocity
+
+                               CALL ssh_nxt( kstp )           ! sea surface height at next time step
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
