Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 3859)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 3865)
@@ -8,9 +8,11 @@
    !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
    !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   ssh_wzv        : after ssh & now vertical velocity
-   !!   ssh_nxt        : filter ans swap the ssh arrays
+   !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   ssh_nxt        : after ssh
+   !!   ssh_swp        : filter ans swap the ssh arrays
+   !!   wzv            : compute now vertical velocity
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers variables
@@ -20,4 +22,5 @@
    USE divcur          ! hor. divergence and curl      (div & cur routines)
    USE iom             ! I/O library
+   USE restart         ! only for lrst_oce
    USE in_out_manager  ! I/O manager
    USE prtctl          ! Print control
@@ -44,6 +47,7 @@
    PRIVATE
 
-   PUBLIC   ssh_wzv    ! called by step.F90
    PUBLIC   ssh_nxt    ! called by step.F90
+   PUBLIC   wzv        ! called by step.F90
+   PUBLIC   ssh_swp    ! called by step.F90
 
    !! * Substitutions
@@ -57,97 +61,36 @@
 CONTAINS
 
-   SUBROUTINE ssh_wzv( kt ) 
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE ssh_wzv  ***
+   SUBROUTINE ssh_nxt( kt )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE ssh_nxt  ***
       !!                   
-      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity
-      !!              and update the now vertical coordinate (lk_vvl=T).
-      !!
-      !! ** Method  : - Using the incompressibility hypothesis, the vertical 
-      !!      velocity is computed by integrating the horizontal divergence  
-      !!      from the bottom to the surface minus the scale factor evolution.
-      !!        The boundary conditions are w=0 at the bottom (no flux) and.
+      !! ** Purpose :   compute the after ssh (ssha)
+      !!
+      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
+      !!      is computed by integrating the horizontal divergence and multiply by
+      !!      by the time step.
       !!
       !! ** action  :   ssha    : after sea surface height
-      !!                wn      : now vertical velocity
-      !!                sshu_a, sshv_a, sshf_a  : after sea surface height (lk_vvl=T)
-      !!                hu, hv, hur, hvr        : ocean depth and its inverse at u-,v-points
       !!
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! time step
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d, zhdiv
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('ssh_wzv')
-      !
-      CALL wrk_alloc( jpi, jpj, z2d, zhdiv ) 
+      !
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv
+      INTEGER, INTENT(in) ::   kt                      ! time step
+      ! 
+      INTEGER             ::   jk                      ! dummy loop indice
+      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
+      !
+      CALL wrk_alloc( jpi, jpj, zhdiv ) 
       !
       IF( kt == nit000 ) THEN
          !
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity '
+         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
          !
-         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
-         !
-         IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only)
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
-                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
-                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
-                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
-                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
-                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
-                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
-                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
-                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
-                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
-                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
-                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
-               END DO
-            END DO
-            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
-            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1      ! NO Vector Opt.
-                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
-                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
-                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
-                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
-               END DO
-            END DO
-            CALL lbc_lnk( sshf_n, 'F', 1. )
-         ENDIF
-         !
-      ENDIF
-
-      !                                           !------------------------------------------!
-      IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case)
-         !                                        !------------------------------------------!
-         DO jk = 1, jpkm1
-            fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays
-            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
-            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
-            !
-            fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays
-            fse3u (:,:,jk) = fse3u_n (:,:,jk)
-            fse3v (:,:,jk) = fse3v_n (:,:,jk)
-            fse3f (:,:,jk) = fse3f_n (:,:,jk)
-            fse3w (:,:,jk) = fse3w_n (:,:,jk)
-            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
-            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
-         END DO
-         !
-         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
-         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
-         !                                            ! now masked inverse of the ocean depth (at u- and v-points)
-         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) )
-         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1._wp - vmask(:,:,1) )
-         ! 
       ENDIF
       !
@@ -162,5 +105,5 @@
       zhdiv(:,:) = 0._wp
       DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
-        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
+        zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
       END DO
       !                                                ! Sea surface elevation time stepping
@@ -181,5 +124,5 @@
 #if defined key_bdy
       ssha(:,:) = ssha(:,:) * bdytmask(:,:)
-      CALL lbc_lnk( ssha, 'T', 1. )                    ! absolutly compulsory !! (jmm)
+      CALL lbc_lnk( ssha, 'T', 1. ) 
 #endif
 #if defined key_asminc
@@ -190,43 +133,102 @@
       ENDIF
 #endif
-      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
-      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1      ! NO Vector Opt.
-               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
-                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
-                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
-               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
-                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
-                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
+
+      !                                           !------------------------------!
+      !                                           !           outputs            !
+      !                                           !------------------------------!
+      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
+      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
+      !
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
+      !
+      CALL wrk_dealloc( jpi, jpj, zhdiv ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
+      !
+   END SUBROUTINE ssh_nxt
+
+   
+   SUBROUTINE wzv( kt )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE wzv  ***
+      !!                   
+      !! ** Purpose :   compute the now vertical velocity
+      !!
+      !! ** Method  : - Using the incompressibility hypothesis, the vertical 
+      !!      velocity is computed by integrating the horizontal divergence  
+      !!      from the bottom to the surface minus the scale factor evolution.
+      !!        The boundary conditions are w=0 at the bottom (no flux) and.
+      !!
+      !! ** action  :   wn      : now vertical velocity
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      !
+      INTEGER, INTENT(in) ::   kt           ! time step
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
+      !
+      INTEGER             ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp)            ::   z1_2dt       ! local scalars
+      !!----------------------------------------------------------------------
+      
+      IF( nn_timing == 1 )  CALL timing_start('wzv')
+      !
+      CALL wrk_alloc( jpi, jpj, z2d ) 
+      CALL wrk_alloc( jpi, jpj, jpk, z3d, zhdiv ) 
+      !
+      IF( kt == nit000 ) THEN
+         !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
+         IF(lwp) WRITE(numout,*) '~~~ '
+         !
+         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
+         !
+      ENDIF
+      !                                           !------------------------------!
+      !                                           !     Now Vertical Velocity    !
+      !                                           !------------------------------!
+      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
+      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
+      !
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
+         DO jk = 1, jpkm1
+            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
+            ! - ML - note: computation allready done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zhdiv(ji,jj,jk) = r1_e12t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
+               END DO
             END DO
          END DO
-         CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions
-      ENDIF
-
-      !                                           !------------------------------!
-      !                                           !     Now Vertical Velocity    !
-      !                                           !------------------------------!
-      z1_2dt = 1.e0 / z2dt
-      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
-         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
-         wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
-            &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
-            &                         * tmask(:,:,jk) * z1_2dt
+         CALL lbc_lnk(zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
+         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
+         !                             ! Same question holds for hdivn. Perhaps just for security
+         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
+            ! computation of w
+            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)                    &
+               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
+         END DO
+         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
+      ELSE   ! z_star and linear free surface cases
+         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
+            ! computation of w
+            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk)                                   &
+               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
+         END DO
+      ENDIF
+
 #if defined key_bdy
          wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
 #endif
-      END DO
 
       !                                           !------------------------------!
       !                                           !           outputs            !
       !                                           !------------------------------!
-      CALL iom_put( "woce", wn                    )   ! vertical velocity
-      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
-      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
+      CALL iom_put( "woce", wn )   ! vertical velocity
       IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
          ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
-         CALL wrk_alloc( jpi,jpj,jpk, z3d )
-         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
+         z2d(:,:) = rau0 * e12t(:,:)
          DO jk = 1, jpk
             z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
@@ -234,17 +236,16 @@
          CALL iom_put( "w_masstr" , z3d                     )  
          CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
-         CALL wrk_dealloc( jpi,jpj,jpk, z3d )
-      ENDIF
-      !
-      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
-      !
-      CALL wrk_dealloc( jpi, jpj, z2d, zhdiv ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('ssh_wzv')
-      !
-   END SUBROUTINE ssh_wzv
-
-
-   SUBROUTINE ssh_nxt( kt )
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi, jpj, z2d  ) 
+      CALL wrk_dealloc( jpi, jpj, jpk, z3d, zhdiv ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('wzv')
+
+
+   END SUBROUTINE wzv
+
+
+   SUBROUTINE ssh_swp( kt )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE ssh_nxt  ***
@@ -252,5 +253,5 @@
       !! ** Purpose :   achieve the sea surface  height time stepping by 
       !!              applying Asselin time filter and swapping the arrays
-      !!              ssha  already computed in ssh_wzv  
+      !!              ssha  already computed in ssh_nxt  
       !!
       !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
@@ -266,86 +267,20 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      !!
-      INTEGER  ::   ji, jj   ! dummy loop indices
-      REAL(wp) ::   zec      ! temporary scalar
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('ssh_swp')
       !
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
+         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
       ENDIF
 
-      !                       !--------------------------!
-      IF( lk_vvl ) THEN       !  Variable volume levels  !     (ssh at t-, u-, v, f-points)
-         !                    !--------------------------!
-         !
-         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
-            sshn  (:,:) = ssha  (:,:)                       ! now <-- after  (before already = now)
-            sshu_n(:,:) = sshu_a(:,:)
-            sshv_n(:,:) = sshv_a(:,:)
-            DO jj = 1, jpjm1                                ! ssh now at f-point
-               DO ji = 1, jpim1      ! NO Vector Opt.
-                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
-                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
-                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
-                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
-               END DO
-            END DO
-            CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions
-            !
-         ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
-            zec = atfp * rdt / rau0
-            DO jj = 1, jpj
-               DO ji = 1, jpi                               ! before <-- now filtered
-                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   &
-                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
-                  sshn  (ji,jj) = ssha  (ji,jj)             ! now <-- after
-                  sshu_n(ji,jj) = sshu_a(ji,jj)
-                  sshv_n(ji,jj) = sshv_a(ji,jj)
-               END DO
-            END DO
-            DO jj = 1, jpjm1                                ! ssh now at f-point
-               DO ji = 1, jpim1      ! NO Vector Opt.
-                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
-                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
-                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
-                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
-               END DO
-            END DO
-            CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions
-            !
-            DO jj = 1, jpjm1                                ! ssh before at u- & v-points
-               DO ji = 1, jpim1      ! NO Vector Opt.
-                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
-                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
-                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
-                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
-                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
-                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
-               END DO
-            END DO
-            CALL lbc_lnk( sshu_b, 'U', 1. )
-            CALL lbc_lnk( sshv_b, 'V', 1. )            !  Boundaries conditions
-            !
-         ENDIF
-         !                    !--------------------------!
-      ELSE                    !        fixed levels      !     (ssh at t-point only)
-         !                    !--------------------------!
-         !
-         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
-            sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now)
-            !
-         ELSE                                               ! Leap-Frog time-stepping: Asselin filter + swap
-            DO jj = 1, jpj
-               DO ji = 1, jpi                               ! before <-- now filtered
-                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
-                  sshn(ji,jj) = ssha(ji,jj)                 ! now <-- after
-               END DO
-            END DO
-         ENDIF
-         !
+      IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
+         sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now)
+      ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
+         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered
+         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
+         sshn(:,:) = ssha(:,:)                           ! now <-- after
       ENDIF
       !
@@ -357,7 +292,7 @@
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
       !
-      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
-      !
-   END SUBROUTINE ssh_nxt
+      IF( nn_timing == 1 )  CALL timing_stop('ssh_swp')
+      !
+   END SUBROUTINE ssh_swp
 
    !!======================================================================
