Index: branches/DEV_r1837_MLF/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- branches/DEV_r1837_MLF/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 2066)
+++ branches/DEV_r1837_MLF/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 2068)
@@ -145,4 +145,7 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3vw_1             !: analytical vertical scale factors at  VW--
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3w_1  , e3uw_1    !:                                       W--UW  points (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3t_b              !: before         -      -      -    -   T      points (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3u_b  , e3v_b     !:   -            -      -      -    -   U--V   points (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3t_d              !: "diffused"     -      -      -    -   T      points (m)
 #else
    LOGICAL, PUBLIC, PARAMETER ::   lk_vvl = .FALSE.   !: fixed grid flag
Index: branches/DEV_r1837_MLF/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- branches/DEV_r1837_MLF/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 2066)
+++ branches/DEV_r1837_MLF/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 2068)
@@ -130,14 +130,14 @@
             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)
+            zv_t_ij       = zs_t(ji  ,jj  ) * sshb(ji  ,jj  )
+            zv_t_ip1j     = zs_t(ji+1,jj  ) * sshb(ji+1,jj  )
+            zv_t_ijp1     = zs_t(ji  ,jj+1) * sshb(ji  ,jj+1)
             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)
+            zv_t_ij       = zs_t(ji  ,jj  ) * sshn(ji  ,jj  )
+            zv_t_ip1j     = zs_t(ji+1,jj  ) * sshn(ji+1,jj  )
+            zv_t_ijp1     = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1)
+            zv_t_ip1jp1   = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1)
             sshu_n(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j )
             sshv_n(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 )
@@ -162,6 +162,6 @@
          END DO
       END DO
-      CALL lbc_lnk( fse3u_b, 'U', 1. )               ! lateral boundary conditions
-      CALL lbc_lnk( fse3v_b, 'U', 1. )
+      CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )               ! lateral boundary conditions
+      CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
       ! Add initial scale factor to scale factor anomaly
       fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
Index: branches/DEV_r1837_MLF/NEMO/OPA_SRC/DOM/domzgr_substitute.h90
===================================================================
--- branches/DEV_r1837_MLF/NEMO/OPA_SRC/DOM/domzgr_substitute.h90	(revision 2066)
+++ branches/DEV_r1837_MLF/NEMO/OPA_SRC/DOM/domzgr_substitute.h90	(revision 2068)
@@ -63,5 +63,6 @@
 #   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)   e3t_m(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_d(i,j,k)   e3t_d(i,j,k)
 
 #   define  fse3t_a(i,j,k)   (fse3t_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k)))
@@ -100,4 +101,5 @@
 
 #   define  fse3t_m(i,j,k)   fse3t_0(i,j,k)
+#   define  fse3t_d(i,j,k)   fse3t_0(i,j,k)
 
 #   define  fse3t_a(i,j,k)   fse3t_0(i,j,k)
Index: branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 2066)
+++ branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 2068)
@@ -22,4 +22,6 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE phycst          ! physical constants
    USE dynspg_oce      ! type of surface pressure gradient
    USE dynadv          ! dynamics: vector invariant versus flux form
@@ -218,5 +220,7 @@
          END DO
       ELSE                                             !* Leap-Frog : Asselin filter and swap
-         IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN          ! applied on velocity
+         !                                ! =============!
+         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
+            !                             ! =============!
             DO jk = 1, jpkm1                              
                DO jj = 1, jpj
@@ -232,62 +236,108 @@
                END DO
             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(:,:,:)
+            !                             ! ================!
+         ELSE                             ! Variable volume !
+            !                             ! ================!
+            ! Before scale factor at t-points
+            ! -------------------------------
+            DO jk = 1, jpkm1
+               fse3t_b(:,:,jk) = fse3t_n(:,:,jk) + atfp * fse3t_d(:,:,jk)
+            ENDDO
             ! 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, 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)
-                     !
-                     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)
-                  END DO
-               END DO
-            END DO
-            CALL lbc_lnk( ub, 'U', -1. )         ! local domain boundaries
-            CALL lbc_lnk( vb, 'V', -1. ) 
+            !
+            IF( ln_dynadv_vec ) THEN
+               ! Before scale factor 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
+               ! lateral boundary conditions
+               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
+               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
+               ! Add initial scale factor to scale factor anomaly
+               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
+               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
+               ! Leap-Frog - Asselin filter and swap: applied on velocity
+               ! -----------------------------------
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpj
+                     DO ji = 1, jpi
+                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
+                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
+                        !
+                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
+                        vb(ji,jj,jk) = zvf
+                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
+                        vn(ji,jj,jk) = va(ji,jj,jk)
+                     END DO
+                  END DO
+               END DO
+               !
+            ELSE
+               ! Temporary filered scale factor at (u/v)-points (will become before scale factor)
+               !-----------------------------------------------
+               ! 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
+               ! lateral boundary conditions
+               CALL lbc_lnk( ze3u_f, 'U', 1. )
+               CALL lbc_lnk( ze3v_f, 'V', 1. )
+               ! Add initial scale factor to scale factor anomaly
+               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
+               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
+               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
+               ! -----------------------------------             ===========================
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpj
+                     DO ji = 1, jpim1
+                        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)
+                        !
+                        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
+                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
+                        vn(ji,jj,jk) = va(ji,jj,jk)
+                     END DO
+                  END DO
+               END DO
+               fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor
+               fse3v_b(:,:,:) = ze3v_f(:,:,:)
+               CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions
+               CALL lbc_lnk( vb, 'V', -1. )
+            ENDIF
+            !
          ENDIF
+         !
       ENDIF
 
Index: branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 2066)
+++ branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 2068)
@@ -104,12 +104,4 @@
             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
-            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  ) )                     &
@@ -118,5 +110,5 @@
                END DO
             END DO
-            CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. )
+            CALL lbc_lnk( sshf_n, 'F', 1. )
          ENDIF
          !
@@ -157,6 +149,4 @@
       !                                           !   After Sea Surface Height   !
       !                                           !------------------------------!
-      z1_rau0 = 0.5 / rau0
-      !
       zhdiv(:,:) = 0.e0
       DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
@@ -166,5 +156,11 @@
       ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
       ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
-      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
+      IF( neuler == 0 .AND. kt == nit000 ) THEN
+         z1_rau0 = 1.e0 / rau0
+         ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 *   emp(:,:)                + zhdiv(:,:) )  ) * tmask(:,:,1)
+      ELSE
+         z1_rau0 = 0.5 / rau0
+         ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
+      ENDIF
 
 #if defined key_obc
@@ -194,7 +190,9 @@
       !                                           !----------------------------------------!
       ! 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(:,:,:)
+      DO jk = 1, jpk
+         fse3t_d(:,:,jk) =          fse3t_b(:,:,jk)   &
+            &              - 2.e0 * fse3t_n(:,:,jk)   &
+            &              +        fse3t_a(:,:,jk)
+      ENDDO
       !                                           !------------------------------!
       !                                           !     Now Vertical Velocity    !
@@ -251,5 +249,5 @@
       !!
       INTEGER  ::   ji, jj   ! dummy loop indices
-      REAL(wp) ::   zec      ! temporary scalare
+      REAL(wp) ::   zec      ! temporary scalar
       !!----------------------------------------------------------------------
 
@@ -281,4 +279,5 @@
             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
Index: branches/DEV_r1837_MLF/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- branches/DEV_r1837_MLF/NEMO/OPA_SRC/IOM/restart.F90	(revision 2066)
+++ branches/DEV_r1837_MLF/NEMO/OPA_SRC/IOM/restart.F90	(revision 2068)
@@ -26,4 +26,5 @@
    USE zdfmxl          ! mixed layer depth
    USE trdmld_oce      ! ocean active mixed layer tracers trends variables
+   USE domvvl          ! variable volume
 
    IMPLICIT NONE
@@ -38,4 +39,5 @@
 
    !! * Substitutions
+#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -123,5 +125,5 @@
       CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb    )
       IF( lk_vvl ) THEN
-         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b )
+         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
       ENDIF
       !
@@ -157,4 +159,5 @@
       REAL(wp) ::   zrdt, zrdttra1
       INTEGER  ::   jlibalt = jprstlib
+      INTEGER  ::   jk                    ! dummy loop indices
       LOGICAL  ::   llok
       !!----------------------------------------------------------------------
@@ -195,5 +198,5 @@
       CALL iom_get( numror, jpdom_autoglo, 'sshb' , sshb  )
       IF( lk_vvl ) THEN
-         CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b )
+         CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
       ENDIF
       !
@@ -224,5 +227,7 @@
          sshb (:,:)   = sshn (:,:)
          IF( lk_vvl ) THEN
-            fse3t_b(:,:,:) = fse3t_n(:,:,:)
+            DO jk = 1, jpk
+               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
+            ENDDO
          ENDIF
       ENDIF
Index: branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 2066)
+++ branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 2068)
@@ -52,5 +52,5 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emps   , emps_b   !: freshwater budget: concentration/dillution   [Kg/m2/s]
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp_tot           !: total E-P-R over ocean and ice               [Kg/m2/s]
-   ! - ML - beginning
+   ! - ML - begin
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_trd_hc_n      !: sbc heat content trend now                   [K/m/s]
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_trd_hc_b      !:  "   "      "      "   before                   "
Index: branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 2066)
+++ branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 2068)
@@ -194,12 +194,12 @@
          utau_b(:,:) = utau(:,:)                         ! Swap the ocean forcing fields
          vtau_b(:,:) = vtau(:,:)                         ! (except at nit000 where before fields
-         qns_b (:,:) = qns (:,:)                         !  are set the end of the routine)
+         qns_b (:,:) = qns (:,:)                         !  are set at the end of the routine)
          qsr_b (:,:) = qsr (:,:)
          emp_b (:,:) = emp (:,:)
          emps_b(:,:) = emps(:,:)
          ! - ML -
-         sbc_trd_hc_b(:,:) = sbc_trd_hc_n(:,:)
-         qsr_trd_hc_b(:,:,:) = qsr_trd_hc_n(:,:,:)
-         IF ( .NOT. lk_vvl )  sbc_trd_sc_b(:,:) = sbc_trd_sc_n(:,:)
+         ! sbc_trd_hc_b(:,:) = sbc_trd_hc_n(:,:)
+         ! qsr_trd_hc_b(:,:,:) = qsr_trd_hc_n(:,:,:)
+         ! IF ( .NOT. lk_vvl )  sbc_trd_sc_b(:,:) = sbc_trd_sc_n(:,:)
          
       ENDIF
@@ -265,11 +265,11 @@
             CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b  )   ! before     solar heat flux (T-point)
             CALL iom_get( numror, jpdom_autoglo, 'emp_b' , emp_b  )   ! before     freshwater flux (T-point)
-            CALL iom_get( numror, jpdom_autoglo, 'emps_b', emp_b  )   ! before C/D freshwater flux (T-point)
+            CALL iom_get( numror, jpdom_autoglo, 'emps_b', emps_b )   ! before C/D freshwater flux (T-point)
             ! - ML -
-            CALL iom_get( numror, jpdom_autoglo, 'sbc_trd_hc_b', sbc_trd_hc_b )   ! before heat content sbc trend
-            CALL iom_get( numror, jpdom_autoglo, 'qsr_trd_hc_b', qsr_trd_hc_b )   ! before heat content trend due to Qsr flux
-            IF ( .NOT. lk_vvl ) THEN
-               CALL iom_get( numror, jpdom_autoglo, 'sbc_trd_sc_b', sbc_trd_sc_b )   ! before salt content sbc trend
-            ENDIF
+            ! CALL iom_get( numror, jpdom_autoglo, 'sbc_trd_hc_b', sbc_trd_hc_b )   ! before heat content sbc trend
+            ! CALL iom_get( numror, jpdom_autoglo, 'qsr_trd_hc_b', qsr_trd_hc_b )   ! before heat content trend due to Qsr flux
+            ! IF ( .NOT. lk_vvl ) THEN
+            !    CALL iom_get( numror, jpdom_autoglo, 'sbc_trd_sc_b', sbc_trd_sc_b )   ! before salt content sbc trend
+            ! ENDIF
             !
          ELSE                                                   !* no restart: set from nit000 values
@@ -291,16 +291,16 @@
             &                    'at it= ', kt,' date= ', ndastp
          IF(lwp) WRITE(numout,*) '~~~~'
-         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau_b )    ! 
-         CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau_b )
-         CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns_b  )
-         CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr_b  )
-         CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp_b  )
-         CALL iom_rstput( kt, nitrst, numrow, 'emps_b' , emps_b )
+         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau )    ! 
+         CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau )
+         CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns  )
+         CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  )
+         CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp  )
+         CALL iom_rstput( kt, nitrst, numrow, 'emps_b' , emps )
          ! - ML -
-         CALL iom_rstput( kt, nitrst, numrow, 'sbc_trd_hc_b', sbc_trd_hc_b )
-         CALL iom_rstput( kt, nitrst, numrow, 'qsr_trd_hc_b', qsr_trd_hc_b )
-         IF ( .NOT. lk_vvl ) THEN
-            CALL iom_rstput( kt, nitrst, numrow, 'sbc_trd_sc_b', sbc_trd_sc_b )
-         ENDIF
+         ! CALL iom_rstput( kt, nitrst, numrow, 'sbc_trd_hc_b', sbc_trd_hc_n )
+         ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_trd_hc_b', qsr_trd_hc_n )
+         ! IF ( .NOT. lk_vvl ) THEN
+         !    CALL iom_rstput( kt, nitrst, numrow, 'sbc_trd_sc_b', sbc_trd_sc_n )
+         ! ENDIF
          !
       ENDIF
Index: branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/trabbc.F90
===================================================================
--- branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 2066)
+++ branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 2068)
@@ -35,6 +35,6 @@
    REAL(wp) ::   rn_geoflx_cst = 86.4e-3      ! Constant value of geothermal heat flux
 
-   INTEGER , DIMENSION(jpi,jpj) ::   nbotlevt   ! ocean bottom level index at T-pt
-   REAL(wp), DIMENSION(jpi,jpj) ::   qgh_trd0   ! geothermal heating trend
+   INTEGER , PUBLIC, DIMENSION(jpi,jpj) ::   nbotlevt   ! ocean bottom level index at T-pt
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qgh_trd0   ! geothermal heating trend
  
    !! * Substitutions
Index: branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 2066)
+++ branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 2068)
@@ -178,5 +178,5 @@
       !!
       INTEGER  ::   ji, jj, jk      ! dummy loop indices
-      REAL(wp) ::   zt_m, zs_m      ! temporary scalars
+      REAL(wp) ::   zt_d, zs_d      ! temporary scalars
       REAL(wp) ::   ztn, zsn        !    -         -
       !!----------------------------------------------------------------------
@@ -203,16 +203,16 @@
                   !                                         ! time laplacian on tracers
                   !                                         ! used for both Asselin and Brown & Campana filters
-                  zt_m = ta(ji,jj,jk) - 2. * tn(ji,jj,jk) + tb(ji,jj,jk)
-                  zs_m = sa(ji,jj,jk) - 2. * sn(ji,jj,jk) + sb(ji,jj,jk)
+                  zt_d = ta(ji,jj,jk) - 2. * tn(ji,jj,jk) + tb(ji,jj,jk)
+                  zs_d = sa(ji,jj,jk) - 2. * sn(ji,jj,jk) + sb(ji,jj,jk)
                   !
                   !                                         ! swap of arrays
-                  tb(ji,jj,jk) = tn(ji,jj,jk) + atfp * zt_m               ! tb <-- tn filtered
-                  sb(ji,jj,jk) = sn(ji,jj,jk) + atfp * zs_m               ! sb <-- sn filtered
+                  tb(ji,jj,jk) = tn(ji,jj,jk) + atfp * zt_d               ! tb <-- tn filtered
+                  sb(ji,jj,jk) = sn(ji,jj,jk) + atfp * zs_d               ! sb <-- sn filtered
                   tn(ji,jj,jk) = ta(ji,jj,jk)                             ! tn <-- ta
                   sn(ji,jj,jk) = sa(ji,jj,jk)                             ! sn <-- sa
                   !                                         ! semi imlicit hpg computation (Brown & Campana)
                   IF( ln_dynhpg_imp ) THEN
-                     ta(ji,jj,jk) = ztn + rbcp * zt_m                     ! ta <-- Brown & Campana average
-                     sa(ji,jj,jk) = zsn + rbcp * zs_m                     ! sa <-- Brown & Campana average
+                     ta(ji,jj,jk) = ztn + rbcp * zt_d                     ! ta <-- Brown & Campana average
+                     sa(ji,jj,jk) = zsn + rbcp * zs_d                     ! sa <-- Brown & Campana average
                   ENDIF
                END DO
@@ -254,6 +254,6 @@
       REAL     ::   ztc_a, ztc_n, ztc_b            !    -         -
       REAL     ::   zsc_a, zsc_n, zsc_b            !    -         -
-      REAL     ::   ztc_f, zsc_f, ztc_m, zsc_m     !    -         -
-      REAL     ::   ze3t_f, ze3t_m                 !    -         -
+      REAL     ::   ztc_f, zsc_f, ztc_d, zsc_d     !    -         -
+      REAL     ::   ze3t_f, ze3t_d                 !    -         -
       REAL     ::   zfact1, zfact2                 !    -         -
       !!----------------------------------------------------------------------
@@ -274,5 +274,5 @@
       ELSE                                             ! apply filter on thickness weighted tracer and swap
          DO jk = 1, jpkm1
-            zfact1 = atfp * r2dt_t(jk)
+            zfact1 = atfp * rdttra(jk)
             zfact2 = zfact1 / rau0
             DO jj = 1, jpj
@@ -282,5 +282,5 @@
                   ze3t_n = fse3t_n(ji,jj,jk)
                   ze3t_a = fse3t_a(ji,jj,jk)
-                  ze3t_m = fse3t_m(ji,jj,jk)
+                  ze3t_d = fse3t_d(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
@@ -290,19 +290,22 @@
                   !                                         ! Time laplacian on tracer contents
                   !                                         ! used for both Asselin and Brown & Campana filters
-                  ztc_m  = ztc_a  - 2. * ztc_n  + ztc_b
-                  zsc_m  = zsc_a  - 2. * zsc_n  + zsc_b
+                  ztc_d  = ztc_a - 2. * ztc_n + ztc_b
+                  zsc_d  = zsc_a - 2. * zsc_n + zsc_b
                   !                                         ! 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
+                  ze3t_f = ze3t_n + atfp * ze3t_d
+                  ztc_f  = ztc_n  + atfp * ztc_d
+                  zsc_f  = zsc_n  + atfp * zsc_d
                   !                                         ! Filter correction
                   IF( jk == 1 ) THEN
-                     ze3t_f = ze3t_f - zfact2 * ( emp_b       (ji,jj)    - emp         (ji,jj)    )
-                     ztc_f  = ztc_f  - zfact1 * ( sbc_trd_hc_n(ji,jj)    - sbc_trd_hc_b(ji,jj)    )
+                     ! WRITE(numout,*) 'filter correction: sbc_trd_hc_n'
+                     ze3t_f = ze3t_f - zfact2 * ( emp_b       (ji,jj) - emp         (ji,jj) )
+                     ztc_f  = ztc_f  - zfact1 * ( sbc_trd_hc_n(ji,jj) - sbc_trd_hc_b(ji,jj) )
                   ENDIF
                   IF( ln_traqsr .AND. ( jk .LE. nksr ) ) THEN
+                     ! WRITE(numout,*) 'jk =', jk
+                     ! WRITE(numout,*) 'filter correction: qsr_trd_hc_n'
                      ztc_f  = ztc_f  - zfact1 * ( qsr_trd_hc_n(ji,jj,jk) - qsr_trd_hc_b(ji,jj,jk) )
                   ENDIF
-                  !                                         ! swap of arrays
+                                                          ! swap of arrays
                   ze3t_f = 1.e0 / ze3t_f
                   tb(ji,jj,jk) = ztc_f * ze3t_f                           ! tb <-- tn filtered
@@ -312,7 +315,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) = 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
+                     ze3t_d       = 1.e0   / ( ze3t_n + rbcp * ze3t_d )
+                     ta(ji,jj,jk) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
+                     sa(ji,jj,jk) = ze3t_d * ( zsc_n  + rbcp * zsc_d  )   ! sa <-- Brown & Campana average
                   ENDIF
                END DO
Index: branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 2066)
+++ branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 2068)
@@ -112,4 +112,18 @@
          ztrdt(:,:,:) = ta(:,:,:) 
          ztrds(:,:,:) = 0.e0
+      ENDIF
+
+      !                                           ! ---------------------------------------- !
+      !                                           !          Swap of forcing field           !
+      !                                           ! ---------------------------------------- !
+      IF( kt /= nit000 ) qsr_trd_hc_b(:,:,:) = qsr_trd_hc_n(:,:,:)
+      !                                           ! ---------------------------------------- !
+      IF( kt == nit000 ) THEN                     !   set the forcing field at nit000 - 1    !
+         !                                        ! ---------------------------------------- !
+         IF( ln_rstart .AND.    &                               !* Restart: read in restart file
+            & iom_varid( numror, 'qsr_trd_hc_b', ldstop = .FALSE. ) > 0 ) THEN 
+            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields red in the restart file'
+            CALL iom_get( numror, jpdom_autoglo, 'qsr_trd_hc_b', qsr_trd_hc_b )   ! before heat content trend due to Qsr flux
+         ENDIF
       ENDIF
 
@@ -225,4 +239,15 @@
             END DO
          END DO
+      ENDIF
+
+      !                                            ! ---------------------------------------- !
+      IF( lrst_oce ) THEN                          !      Write in the ocean restart file     !
+         !                                         ! ---------------------------------------- !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'qsr : penetrative solar radiation forcing field written in ocean restart file ',   &
+            &                    'at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*) '~~~~'
+         CALL iom_rstput( kt, nitrst, numrow, 'qsr_trd_hc_b', qsr_trd_hc_n )
+         !
       ENDIF
 
Index: branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 2066)
+++ branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 2068)
@@ -20,5 +20,7 @@
    USE trdmod          ! ocean trends 
    USE trdmod_oce      ! ocean variables trends
+   USE iom
    USE in_out_manager  ! I/O manager
+   USE restart         ! ocean restart
    USE prtctl          ! Print control
 
@@ -132,8 +134,26 @@
          ENDIF
       ENDIF
-
-      !                             ! ---------------------- !
-      IF( lk_vvl ) THEN             !  Variable Volume case  !
-         !                          ! ---------------------- !
+      !                                            ! ---------------------------------------- !
+      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
+         !                                         ! ---------------------------------------- !
+         sbc_trd_hc_b(:,:) = sbc_trd_hc_n(:,:)                         ! Swap the ocean forcing fields except at nit000
+         IF ( .NOT. lk_vvl ) sbc_trd_sc_b(:,:)   = sbc_trd_sc_n(:,:)
+      ENDIF
+      !                                            ! ---------------------------------------- !
+      IF( kt == nit000 ) THEN                      !   set the forcing field at nit000 - 1    !
+         !                                         ! ---------------------------------------- !
+         IF( ln_rstart .AND.    &                               !* Restart: read in restart file
+            & iom_varid( numror, 'sbc_trd_hc_b', ldstop = .FALSE. ) > 0 ) THEN 
+            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields red in the restart file'
+            CALL iom_get( numror, jpdom_autoglo, 'sbc_trd_hc_b', sbc_trd_hc_b )   ! before heat content sbc trend
+            CALL iom_get( numror, jpdom_autoglo, 'qsr_trd_hc_b', qsr_trd_hc_b )   ! before heat content trend due to Qsr flux
+            IF ( .NOT. lk_vvl ) THEN
+               CALL iom_get( numror, jpdom_autoglo, 'sbc_trd_sc_b', sbc_trd_sc_b )   ! before salt content sbc trend
+            ENDIF
+         ENDIF
+      ENDIF
+      !                                            ! ---------------------- !
+      IF( lk_vvl ) THEN                            !  Variable Volume case  !
+         !                                         ! ---------------------- !
 !!gm BUG : in key_vvl emps must be modified to only include the salt flux due to sea-ice freezing/melting
 !!gm       otherwise this flux will be missing  ==> modification required in limsbc,  limsbc_2 and CICE interface.s
@@ -161,7 +181,7 @@
             END DO
          ENDIF
-         !                          ! ---------------------- !
-      ELSE                          !  Constant Volume case  !
-         !                          ! ---------------------- !
+         !                                         ! ---------------------- !
+      ELSE                                         !  Constant Volume case  !
+         !                                         ! ---------------------- !
          IF ( neuler == 0 .AND. kt == nit000 ) THEN
             DO jj = 2, jpj
@@ -197,4 +217,16 @@
       ENDIF
 
+      !                                            ! ---------------------------------------- !
+      IF( lrst_oce ) THEN                          !      Write in the ocean restart file     !
+         !                                         ! ---------------------------------------- !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ',   &
+            &                    'at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*) '~~~~'
+         CALL iom_rstput( kt, nitrst, numrow, 'sbc_trd_hc_b', sbc_trd_hc_n )
+         IF ( .NOT. lk_vvl ) CALL iom_rstput( kt, nitrst, numrow, 'sbc_trd_sc_b', sbc_trd_sc_n )
+         !
+      ENDIF
+
       IF( l_trdtra ) THEN           ! save the sbc trends for diagnostic
          ztrdt(:,:,:) = ta(:,:,:) - ztrdt(:,:,:)
@@ -205,5 +237,4 @@
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' sbc  - Ta: ', mask1=tmask,   &
          &                       tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-      !
    END SUBROUTINE tra_sbc
 
Index: branches/DEV_r1837_MLF/NEMO/OPA_SRC/istate.F90
===================================================================
--- branches/DEV_r1837_MLF/NEMO/OPA_SRC/istate.F90	(revision 2066)
+++ branches/DEV_r1837_MLF/NEMO/OPA_SRC/istate.F90	(revision 2068)
@@ -66,4 +66,6 @@
       !!----------------------------------------------------------------------
       USE eosbn2          ! eq. of state, Brunt Vaisala frequency (eos     routine)
+      ! - ML - needed for initialization of e3t_b
+      INTEGER  ::  jk     ! dummy loop indice
 
       IF(lwp) WRITE(numout,*)
@@ -128,5 +130,11 @@
             &                                  gtu, gsu, gru, &  ! of t, s, rd at the bottom ocean level
             &                                  gtv, gsv, grv )
-         
+      
+         ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here
+         IF( lk_vvl ) THEN
+            DO jk = 1, jpk
+               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
+            ENDDO
+         ENDIF
       ENDIF
       !
Index: branches/DEV_r1837_MLF/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/DEV_r1837_MLF/NEMO/OPA_SRC/step.F90	(revision 2066)
+++ branches/DEV_r1837_MLF/NEMO/OPA_SRC/step.F90	(revision 2068)
@@ -195,4 +195,5 @@
       !  Ocean dynamics : ssh, wn, hdiv, rot                                 !
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
                          CALL ssh_wzv( kstp )         ! after ssh & vertical velocity
 
