Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 2004)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 2005)
@@ -51,5 +51,7 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zcoefu, zcoefv, zcoeff
+      REAL(wp) ::   zcoefu , zcoefv   , zcoeff                   ! temporary scalars
+      REAL(wp) ::   zv_t_ij, zv_t_ip1j, zv_t_ijp1, zv_t_ip1jp1   !     -        -
+      REAL(wp), DIMENSION(jpi,jpj) ::  zs_t, zs_u_1, zs_v_1      !     -     2D workspace
       !!----------------------------------------------------------------------
 
@@ -115,27 +117,54 @@
          hv_0(:,:) = hv_0(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk)
       END DO
+      
+      ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
+      ! for ssh and scale factors
+      zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
+      zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:)
+      zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:)
 
       DO jj = 1, jpjm1                          ! initialise before and now Sea Surface Height at u-, v-, f-points
          DO ji = 1, jpim1   ! NO vector opt.
-            zcoefu = 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) )   
-            sshf_b(ji,jj) = zcoeff * ( sshb(ji  ,jj) + sshb(ji  ,jj+1)                 &
-               &                     + sshb(ji+1,jj) + sshb(ji+1,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) )     
-            sshf_n(ji,jj) = zcoeff * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)                 &
-               &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) )               
+            zcoefu = umask(ji,jj,1) * zs_u_1(ji,jj)
+            zcoefv = vmask(ji,jj,1) * zs_v_1(ji,jj)
+            zcoeff = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) / ( e1f(ji,jj) * e2f(ji,jj) )
+            ! before fields
+            zv_t_ij       = zs_t(ji  ,jj  ) * sshb(ji  ,jj  ,jk)
+            zv_t_ip1j     = zs_t(ji+1,jj  ) * sshb(ji+1,jj  ,jk)
+            zv_t_ijp1     = zs_t(ji  ,jj+1) * sshb(ji  ,jj+1,jk)
+            sshu_b(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j )
+            sshv_b(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 )
+            ! now fields
+            zv_t_ij       = zs_t(ji  ,jj  ) * sshn(ji  ,jj  ,jk)
+            zv_t_ip1j     = zs_t(ji+1,jj  ) * sshn(ji+1,jj  ,jk)
+            zv_t_ijp1     = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1,jk)
+            zv_t_ip1jp1   = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1,jk)
+            sshu_n(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j )
+            sshv_n(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 )
+            sshf_n(ji,jj) = zcoeff * ( zv_t_ij + zv_t_ip1j + zv_t_ijp1 + zv_t_ip1jp1 )
          END DO
       END DO
-      CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )      ! lateral boundary conditions
-      CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
-      CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. )
+      CALL lbc_lnk( sshu_n, 'U', 1. )   ;   CALL lbc_lnk( sshu_b, 'U', 1. )      ! lateral boundary conditions
+      CALL lbc_lnk( sshv_n, 'V', 1. )   ;   CALL lbc_lnk( sshv_b, 'V', 1. )
+      CALL lbc_lnk( sshf_n, 'F', 1. )
+
+                                                ! initialise before scale factors at (u/v)-points
+      ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
+      DO jk = 1, jpkm1
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
+               zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
+               zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
+               fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
+               fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( fse3u_b, 'U', 1. )               ! lateral boundary conditions
+      CALL lbc_lnk( fse3v_b, 'U', 1. )
+      ! Add initial scale factor to scale factor anomaly
+      fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
+      fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
       !
    END SUBROUTINE dom_vvl
Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DOM/domzgr_substitute.h90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DOM/domzgr_substitute.h90	(revision 2004)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DOM/domzgr_substitute.h90	(revision 2005)
@@ -46,12 +46,7 @@
 #   define  fse3vw(i,j,k)  e3vw_1(i,j,k)
 
-#   define  fsdept_b(i,j,k)  (fsdept_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k)))
-#   define  fsdepw_b(i,j,k)  (fsdepw_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k)))
-#   define  fsde3w_b(i,j,k)  (fsde3w_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k))-sshb(i,j))
-#   define  fse3t_b(i,j,k)   (fse3t_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k)))
-#   define  fse3u_b(i,j,k)   (fse3u_0(i,j,k)*(1.+sshu_b(i,j)*muu(i,j,k)))
-#   define  fse3v_b(i,j,k)   (fse3v_0(i,j,k)*(1.+sshv_b(i,j)*muv(i,j,k)))
-#   define  fse3f_b(i,j,k)   (fse3f_0(i,j,k)*(1.+sshf_b(i,j)*muf(i,j,k)))
-#   define  fse3w_b(i,j,k)   (fse3w_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k)))
+#   define  fse3t_b(i,j,k)   e3t_b(i,j,k)
+#   define  fse3u_b(i,j,k)   e3u_b(i,j,k)
+#   define  fse3v_b(i,j,k)   e3v_b(i,j,k)
 #   define  fse3uw_b(i,j,k)  (fse3uw_0(i,j,k)*(1.+sshu_b(i,j)*muu(i,j,k)))
 #   define  fse3vw_b(i,j,k)  (fse3vw_0(i,j,k)*(1.+sshv_b(i,j)*muv(i,j,k)))
@@ -68,16 +63,9 @@
 #   define  fse3vw_n(i,j,k)  (fse3vw_0(i,j,k)*(1.+sshv_n(i,j)*muv(i,j,k)))
 
-#   define  fse3t_m(i,j,k)   (fse3t_0(i,j,k)*(1.+ssh_m(i,j)*mut(i,j,k)))
+#   define  fse3t_m(i,j,k)   e3t_m(i,j,k)
 
-#   define  fsdept_a(i,j,k)  (fsdept_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k)))
-#   define  fsdepw_a(i,j,k)  (fsdepw_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k)))
-#   define  fsde3w_a(i,j,k)  (fsde3w_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k))-ssha(i,j))
 #   define  fse3t_a(i,j,k)   (fse3t_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k)))
 #   define  fse3u_a(i,j,k)   (fse3u_0(i,j,k)*(1.+sshu_a(i,j)*muu(i,j,k)))
 #   define  fse3v_a(i,j,k)   (fse3v_0(i,j,k)*(1.+sshv_a(i,j)*muv(i,j,k)))
-#   define  fse3f_a(i,j,k)   (fse3f_0(i,j,k)*(1.+sshf_a(i,j)*muf(i,j,k)))
-#   define  fse3w_a(i,j,k)   (fse3w_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k)))
-#   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
@@ -94,12 +82,7 @@
 #   define  fse3vw(i,j,k)  fse3vw_0(i,j,k)
 
-#   define  fsdept_b(i,j,k)  fsdept_0(i,j,k)
-#   define  fsdepw_b(i,j,k)  fsdepw_0(i,j,k)
-#   define  fsde3w_b(i,j,k)  fsde3w_0(i,j,k)
 #   define  fse3t_b(i,j,k)   fse3t_0(i,j,k)
 #   define  fse3u_b(i,j,k)   fse3u_0(i,j,k)
 #   define  fse3v_b(i,j,k)   fse3v_0(i,j,k)
-#   define  fse3f_b(i,j,k)   fse3f_0(i,j,k)
-#   define  fse3w_b(i,j,k)   fse3w_0(i,j,k)
 #   define  fse3uw_b(i,j,k)  fse3uw_0(i,j,k)
 #   define  fse3vw_b(i,j,k)  fse3vw_0(i,j,k)
@@ -118,14 +101,7 @@
 #   define  fse3t_m(i,j,k)   fse3t_0(i,j,k)
 
-#   define  fsdept_a(i,j,k)  fsdept_0(i,j,k)
-#   define  fsdepw_a(i,j,k)  fsdepw_0(i,j,k)
-#   define  fsde3w_a(i,j,k)  fsde3w_0(i,j,k)
 #   define  fse3t_a(i,j,k)   fse3t_0(i,j,k)
 #   define  fse3u_a(i,j,k)   fse3u_0(i,j,k)
 #   define  fse3v_a(i,j,k)   fse3v_0(i,j,k)
-#   define  fse3f_a(i,j,k)   fse3f_0(i,j,k)
-#   define  fse3w_a(i,j,k)   fse3w_0(i,j,k)
-#   define  fse3uw_a(i,j,k)  fse3uw_0(i,j,k)
-#   define  fse3vw_a(i,j,k)  fse3vw_0(i,j,k)
 #endif
    !!----------------------------------------------------------------------
Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 2004)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 2005)
@@ -87,4 +87,8 @@
       !!               un,vn   now horizontal velocity of next time-step
       !!----------------------------------------------------------------------
+#if defined key_vvl
+      USE oce, ONLY :   ze3u_f => ta   ! use ta as 3D workspace
+      USE oce, ONLY :   ze3v_f => sa   ! use sa as 3D workspace
+#endif
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
       !!
@@ -95,7 +99,9 @@
       REAL(wp) ::   zue3a , zue3n , zue3b    ! temporary scalar
       REAL(wp) ::   zve3a , zve3n , zve3b    !    -         -
-      REAL(wp) ::   ze3u_b, ze3u_n, ze3u_a   !    -         -
-      REAL(wp) ::   ze3v_b, ze3v_n, ze3v_a   !    -         - 
       REAL(wp) ::   zuf   , zvf              !    -         - 
+      REAL(wp) ::   zec                      !    -         - 
+      REAL(wp) ::   zv_t_ij  , zv_t_ip1j     !     -        -
+      REAL(wp) ::   zv_t_ijp1                !     -        -
+      REAL(wp), DIMENSION(jpi,jpj) ::  zs_t, zs_u_1, zs_v_1      ! temporary 2D workspace
       !!----------------------------------------------------------------------
 
@@ -216,6 +222,6 @@
                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)
+                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
+                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
                      !
                      ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
@@ -227,28 +233,53 @@
             END DO
          ELSE                                                ! applied on thickness weighted velocity
+            ! Before scale factors at (t/u/v)-points (actually "now filtered" and futur "before")
+            ! ======================================
+            ! Scale factor at t-points
+            ! ------------------------
+            fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * fse3t_m(:,:,:)
+            ! Add volume filter correction only at the first level of t-point scale factors
+            zec = atfp * rdt / rau0
+            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
+            ! Scale factor at (u/v)-points
+            ! ------------------------
+            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
+            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
+            zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:)
+            zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:)
+            ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
+            DO jk = 1, jpkm1
+               DO jj = 1, jpjm1
+                  DO ji = 1, jpim1
+                     zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
+                     zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
+                     zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
+                     ze3u_f(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
+                     ze3v_f(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
+                  END DO
+               END DO
+            END DO
+            CALL lbc_lnk( ze3u_f, 'U', 1. )               ! lateral boundary conditions
+            CALL lbc_lnk( ze3v_f, 'U', 1. )
+            ! Add initial scale factor to scale factor anomaly
+            ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
+            ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
+            
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     ze3u_a = fse3u_a(ji,jj,jk)
-                     ze3v_a = fse3v_a(ji,jj,jk)
-                     ze3u_n = fse3u_n(ji,jj,jk)
-                     ze3v_n = fse3v_n(ji,jj,jk)
-                     ze3u_b = fse3u_b(ji,jj,jk)
-                     ze3v_b = fse3v_b(ji,jj,jk)
+                  DO ji = 1, jpim
+                     zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
+                     zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
+                     zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
+                     zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
+                     zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
+                     zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
                      !
-                     zue3a = ua(ji,jj,jk) * ze3u_a
-                     zve3a = va(ji,jj,jk) * ze3v_a
-                     zue3n = un(ji,jj,jk) * ze3u_n
-                     zve3n = vn(ji,jj,jk) * ze3v_n
-                     zue3b = ub(ji,jj,jk) * ze3u_b
-                     zve3b = vb(ji,jj,jk) * ze3v_b
-                     !
-                     zuf  = ( atfp * ( zue3b  + zue3a  ) + atfp1 * zue3n  )   &
-                        & / ( atfp * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk)
-                     zvf  = ( atfp * ( zve3b  + zve3a  ) + atfp1 * zve3n  )   &
-                        & / ( atfp * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk)
+                     zuf  = ( zue3n + atfp * ( zue3b  - 2.e0 * zue3n  + zue3a  ) ) / ze3u_f(ji,jj,jk)
+                     zvf  = ( zve3n + atfp * ( zve3b  - 2.e0 * zve3n  + zve3a  ) ) / ze3v_f(ji,jj,jk)
                      !
                      ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
                      vb(ji,jj,jk) = zvf
+                     fse3u_b(ji,jj,jk) = ze3u_f(ji,jj,jk)    ! e3u_b <-- filtered scale factor
+                     fse3v_b(ji,jj,jk) = ze3v_f(ji,jj,jk)
                      un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
                      vn(ji,jj,jk) = va(ji,jj,jk)
@@ -256,4 +287,6 @@
                END DO
             END DO
+            CALL lbc_lnk( ub, 'U', -1. )         ! local domain boundaries
+            CALL lbc_lnk( vb, 'V', -1. ) 
          ENDIF
       ENDIF
Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 2004)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 2005)
@@ -189,16 +189,12 @@
          CALL lbc_lnk( sshu_a, 'U', 1. )
          CALL lbc_lnk( sshv_a, 'V', 1. )
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1      ! NO Vector Opt.
-               sshf_a(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_a(ji,jj  )     &
-                    &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_a(ji,jj+1) )
-            END DO
-         END DO
-         ! Boundaries conditions
-         CALL lbc_lnk( sshf_a, 'F', 1. )
-      ENDIF
-
+      ENDIF
+      !                                           !----------------------------------------!
+      !                                           !     vertical scale factor laplacian    !
+      !                                           !----------------------------------------!
+      ! Needed for Robert-Asselin time filter and for Brown & Campana semi implicit hydrostatic presure gradient
+      fse3t_m(:,:,:) =          fse3t_b(:,:,:)   &
+         &             - 2.e0 * fse3t_n(:,:,:)   &
+         &             +        fse3t_a(:,:,:)
       !                                           !------------------------------!
       !                                           !     Now Vertical Velocity    !
@@ -219,4 +215,5 @@
       CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
       IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
+         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
          z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
          DO jk = 1, jpk
@@ -264,13 +261,24 @@
 
       !                       !--------------------------!
-      IF( lk_vvl ) THEN       !  Variable volume levels  !   ssh at t-, u-, v, f-points
+      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(:,:)
-            sshf_n(:,:) = sshf_a(:,:)
+            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
+            ! Boundaries conditions
+            CALL lbc_lnk( sshf_n, 'F', 1. )
          ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
-            zec = atfp * rdt / rau0
             DO jj = 1, jpj
                DO ji = 1, jpi                                ! before <-- now filtered
@@ -280,7 +288,16 @@
                   sshu_n(ji,jj) = sshu_a(ji,jj)
                   sshv_n(ji,jj) = sshv_a(ji,jj)
-                  sshf_n(ji,jj) = sshf_a(ji,jj)
-               END DO
-            END DO
+               END DO
+            END DO
+            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
+            ! Boundaries conditions
+            CALL lbc_lnk( sshf_n, 'F', 1. )
             DO jj = 1, jpjm1
                DO ji = 1, jpim1      ! NO Vector Opt.
@@ -296,18 +313,11 @@
             CALL lbc_lnk( sshu_b, 'U', 1. )
             CALL lbc_lnk( sshv_b, 'V', 1. )
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1      ! NO Vector Opt.
-                  sshf_b(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_b(ji,jj  )     &
-                     &                   + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_b(ji,jj+1) )
-               END DO
-            END DO
-            ! Boundaries conditions
-            CALL lbc_lnk( sshf_b, 'F', 1. )
          ENDIF
          !                    !--------------------------!
-      ELSE                    !        fixed levels      !   ssh at t-point only
+      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)
Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/IOM/restart.F90	(revision 2004)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/IOM/restart.F90	(revision 2005)
@@ -122,4 +122,7 @@
       CALL iom_rstput( kt, nitrst, numrow, 'hdivb'  , hdivb   )
       CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb    )
+      IF( lk_vvl ) THEN
+         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b )
+      ENDIF
       !
       CALL iom_rstput( kt, nitrst, numrow, 'un'     , un      )     ! now fields
@@ -191,4 +194,7 @@
       CALL iom_get( numror, jpdom_autoglo, 'hdivb', hdivb )
       CALL iom_get( numror, jpdom_autoglo, 'sshb' , sshb  )
+      IF( lk_vvl ) THEN
+         CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b )
+      ENDIF
       !
       CALL iom_get( numror, jpdom_autoglo, 'un'   , un    )        ! now    fields
@@ -217,5 +223,7 @@
          hdivb(:,:,:) = hdivn(:,:,:)
          sshb (:,:)   = sshn (:,:)
-         ! - ML - sshbnc 
+         IF( lk_vvl ) THEN
+            fse3t_b(:,:,:) = fse3t_n(:,:,:)
+         ENDIF
       ENDIF
       !
Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 2004)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 2005)
@@ -282,4 +282,5 @@
                   ze3t_n = fse3t_n(ji,jj,jk)
                   ze3t_a = fse3t_a(ji,jj,jk)
+                  ze3t_m = fse3t_m(ji,jj,jk)
                   !                                         ! tracer content at Before, now and after
                   ztc_b  = tb(ji,jj,jk) * ze3t_b   ;   zsc_b = sb(ji,jj,jk) * ze3t_b
@@ -287,14 +288,13 @@
                   ztc_a  = ta(ji,jj,jk) * ze3t_a   ;   zsc_a = sa(ji,jj,jk) * ze3t_a
                   !
-                  !                                         ! Time laplacian on thickness and tracer content
+                  !                                         ! Time laplacian on tracer contents
                   !                                         ! used for both Asselin and Brown & Campana filters
-                  ze3t_m = ze3t_a - 2. * ze3t_n + ze3t_b
                   ztc_m  = ztc_a  - 2. * ztc_n  + ztc_b
                   zsc_m  = zsc_a  - 2. * zsc_n  + zsc_b
-                  !                                         ! Asselin Filter + correction
+                  !                                         ! Asselin Filter on thicknesses and tracer contents
                   ze3t_f = ze3t_n + atfp * ze3t_m
                   ztc_f  = ztc_n  + atfp * ztc_m
                   zsc_f  = zsc_n  + atfp * zsc_m
-
+                  !                                         ! Filter correction
                   IF( jk == 1 ) THEN
                      ze3t_f = ze3t_f - zfact2 * ( emp_b       (ji,jj)    - emp         (ji,jj)    )
@@ -312,7 +312,7 @@
                   !                                         ! semi imlicit hpg computation (Brown & Campana)
                   IF( ln_dynhpg_imp ) THEN
-                     ze3t_m = 1.e0 / ( ze3t_n + rbcp * ze3t_m )
-                     ta(ji,jj,jk) =  ( ztc_n  + rbcp * ztc_m  ) * ze3t_m  ! ta <-- Brown & Campana average
-                     sa(ji,jj,jk) =  ( zsc_n  + rbcp * zsc_m  ) * ze3t_m  ! sa <-- Brown & Campana average
+                     ze3t_m       = 1.e0   / ( ze3t_n + rbcp * ze3t_m )
+                     ta(ji,jj,jk) = ze3t_m * ( ztc_n  + rbcp * ztc_m  )   ! ta <-- Brown & Campana average
+                     sa(ji,jj,jk) = ze3t_m * ( zsc_n  + rbcp * zsc_m  )   ! sa <-- Brown & Campana average
                   ENDIF
                END DO
Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/oce.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/oce.F90	(revision 2004)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/oce.F90	(revision 2005)
@@ -47,5 +47,5 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshu_b ,  sshu_n  ,  sshu_a  !: sea surface height at u-point [m]
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshv_b ,  sshv_n  ,  sshv_a  !: sea surface height at u-point [m]
-   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) ::             sshf_n             !: sea surface height at f-point [m]
 
    !!----------------------------------------------------------------------
