Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynzdf_exp.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynzdf_exp.F90	(revision 1869)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynzdf_exp.F90	(revision 1870)
@@ -4,7 +4,8 @@
    !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
    !!==============================================================================
-   !! History :      !  90-10  (B. Blanke)  Original code
-   !!                !  97-05  (G. Madec)  vertical component of isopycnal
-   !!           8.5  !  02-08  (G. Madec)  F90: Free form and module
+   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
+   !!            8.0  !  1997-05  (G. Madec)  vertical component of isopycnal
+   !!   NEMO     1.0  !  1002-08  (G. Madec)  F90: Free form and module
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
    !!----------------------------------------------------------------------
 
@@ -13,5 +14,4 @@
    !!                  sion using an explicit time-stepping scheme.
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
@@ -24,6 +24,5 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC dyn_zdf_exp    ! called by step.F90
+   PUBLIC   dyn_zdf_exp   ! called by step.F90
 
    !! * Substitutions
@@ -31,5 +30,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -47,5 +46,5 @@
       !!      technique. The vertical diffusion of momentum is given by:
       !!         diffu = dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ub) )
-      !!      Surface boundary conditions: wind stress input
+      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
       !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F90)
       !!      Add this trend to the general trend ua :
@@ -54,37 +53,29 @@
       !! ** Action : - Update (ua,va) with the vertical diffusive trend
       !!---------------------------------------------------------------------
-      !! * Arguments
-      INTEGER , INTENT( in ) ::   kt                           ! ocean time-step index
-      REAL(wp), INTENT( in ) ::   p2dt                         ! time-step 
-
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk, jl                              ! dummy loop indices
-      REAL(wp) ::   zrau0r, zlavmr, zua, zva                   ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpk) ::   zwx, zwy, zwz, zww     ! temporary workspace arrays
+      INTEGER , INTENT(in) ::   kt     ! ocean time-step index
+      REAL(wp), INTENT(in) ::   p2dt   ! time-step 
+      !!
+      INTEGER ::   ji, jj, jk, jl                            ! dummy loop indices
+      REAL(wp) ::   zrau0r, zlavmr, zua, zva                 ! temporary scalars
+      REAL(wp), DIMENSION(jpi,jpk) ::   zwx, zwy, zwz, zww   ! 2D workspace
       !!----------------------------------------------------------------------
 
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn_zdf_exp : vertical momentum diffusion explicit operator'
+         IF(lwp) WRITE(numout,*) 'dyn_zdf_exp : vertical momentum diffusion - explicit operator'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
       ENDIF
 
-      ! Local constant initialization
-      ! -----------------------------
-      zrau0r = 1. / rau0                                   ! inverse of the reference density
-      zlavmr = 1. / float( nn_zdfexp )                      ! inverse of the number of sub time step
+      zrau0r = 1. / rau0                              ! Local constant initialization
+      zlavmr = 1. / REAL( nn_zdfexp )
 
-      !                                                ! ===============
-      DO jj = 2, jpjm1                                 !  Vertical slab
-         !                                             ! ===============
-
-         ! Surface boundary condition
-         DO ji = 2, jpim1
-            zwy(ji,1) = utau(ji,jj) * zrau0r
-            zww(ji,1) = vtau(ji,jj) * zrau0r
+      !                                               ! ===============
+      DO jj = 2, jpjm1                                !  Vertical slab
+         !                                            ! ===============
+         DO ji = 2, jpim1         ! Surface boundary condition
+            zwy(ji,1) = ( utau_b(ji,jj) + utau(ji,jj) ) * zrau0r
+            zww(ji,1) = ( vtau_b(ji,jj) + vtau(ji,jj) ) * zrau0r
          END DO  
-
-         ! Initialization of x, z and contingently trends array
-         DO jk = 1, jpk
+         DO jk = 1, jpk         ! Initialization of x, z and contingently trends array
             DO ji = 2, jpim1
                zwx(ji,jk) = ub(ji,jj,jk)
@@ -92,10 +83,8 @@
             END DO  
          END DO  
-
-         ! Time splitting loop
-         DO jl = 1, nn_zdfexp
-
-            ! First vertical derivative
-            DO jk = 2, jpk
+         !
+         DO jl = 1, nn_zdfexp     ! Time splitting loop
+            !
+            DO jk = 2, jpk            ! First vertical derivative
                DO ji = 2, jpim1
                   zwy(ji,jk) = avmu(ji,jj,jk) * ( zwx(ji,jk-1) - zwx(ji,jk) ) / fse3uw(ji,jj,jk) 
@@ -103,24 +92,20 @@
                END DO  
             END DO  
-
-            ! Second vertical derivative and trend estimation at kt+l*rdt/nn_zdfexp
-            DO jk = 1, jpkm1
+            DO jk = 1, jpkm1          ! Second vertical derivative and trend estimation at kt+l*rdt/nn_zdfexp
                DO ji = 2, jpim1
-                  zua = zlavmr*( zwy(ji,jk) - zwy(ji,jk+1) ) / fse3u(ji,jj,jk)
-                  zva = zlavmr*( zww(ji,jk) - zww(ji,jk+1) ) / fse3v(ji,jj,jk)
+                  zua = zlavmr * ( zwy(ji,jk) - zwy(ji,jk+1) ) / fse3u(ji,jj,jk)
+                  zva = zlavmr * ( zww(ji,jk) - zww(ji,jk+1) ) / fse3v(ji,jj,jk)
                   ua(ji,jj,jk) = ua(ji,jj,jk) + zua
                   va(ji,jj,jk) = va(ji,jj,jk) + zva
-
-                  zwx(ji,jk) = zwx(ji,jk) + p2dt*zua*umask(ji,jj,jk)
-                  zwz(ji,jk) = zwz(ji,jk) + p2dt*zva*vmask(ji,jj,jk)
+                  !
+                  zwx(ji,jk) = zwx(ji,jk) + p2dt * zua * umask(ji,jj,jk)
+                  zwz(ji,jk) = zwz(ji,jk) + p2dt * zva * vmask(ji,jj,jk)
                END DO  
             END DO  
-
+            !
          END DO  
-
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-
+         !                                            ! ===============
+      END DO                                          !   End of slab
+      !                                               ! ===============
    END SUBROUTINE dyn_zdf_exp
 
Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 1869)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 1870)
@@ -4,7 +4,8 @@
    !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
    !!==============================================================================
-   !! History :      !  90-10  (B. Blanke)  Original code
-   !!                !  97-05  (G. Madec)  vertical component of isopycnal
-   !!           8.5  !  02-08  (G. Madec)  F90: Free form and module
+   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
+   !!            8.0  !  1997-05  (G. Madec)  vertical component of isopycnal
+   !!   NEMO     1.0  !  2002-08  (G. Madec)  F90: Free form and module
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
    !!----------------------------------------------------------------------
 
@@ -13,9 +14,4 @@
    !!                  sion using a implicit time-stepping.
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
-   !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
@@ -28,6 +24,5 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC dyn_zdf_imp    ! called by step.F90
+   PUBLIC   dyn_zdf_imp   ! called by step.F90
 
    !! * Substitutions
@@ -35,11 +30,10 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
 CONTAINS
-
 
    SUBROUTINE dyn_zdf_imp( kt, p2dt )
@@ -54,25 +48,21 @@
       !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
       !!      backward time stepping
-      !!      Surface boundary conditions: wind stress input
+      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
       !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
       !!      Add this trend to the general trend ua :
       !!         ua = ua + dz( avmu dz(u) )
       !!
-      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive
-      !!               mixing trend.
+      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
       !!---------------------------------------------------------------------
-      !! * Modules used
-      USE oce, ONLY :  zwd   => ta,   &                ! use ta as workspace
-                       zws   => sa                     ! use sa as workspace
-
-      !! * Arguments
-      INTEGER , INTENT( in ) ::   kt                   ! ocean time-step index
-      REAL(wp), INTENT( in ) ::  p2dt                  ! vertical profile of tracer time-step
-
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk                          ! dummy loop indices
-      REAL(wp) ::   zrau0r, z2dtf, zcoef, zzws, zrhs   ! temporary scalars
-      REAL(wp) ::   zzwi                               ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! temporary workspace arrays
+      USE oce, ONLY :  zwd   => ta      ! use ta as workspace
+      USE oce, ONLY :  zws   => sa      ! use sa as workspace
+      !!
+      INTEGER , INTENT( in ) ::   kt    ! ocean time-step index
+      REAL(wp), INTENT( in ) ::  p2dt   ! vertical profile of tracer time-step
+      !!
+      INTEGER  ::   ji, jj, jk             ! dummy loop indices
+      REAL(wp) ::   zrau0r, z2dtf, zcoef   ! temporary scalars
+      REAL(wp) ::   zzwi, zzws, zrhs       ! temporary scalars
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! 3D workspace
       !!----------------------------------------------------------------------
 
@@ -95,19 +85,18 @@
       ! friction velocity in dyn_bfr using values from the previous timestep. There
       ! is no need to include these in the implicit calculation.
-      DO jk = 1, jpkm1
+      !
+      DO jk = 1, jpkm1        ! Matrix
          DO jj = 2, jpjm1 
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zcoef = - p2dt / fse3u(ji,jj,jk)
-               zzwi          = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  )
-               zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk)
-               zzws          = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
-               zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1)
+               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
+               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
+               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
+               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
                zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
             END DO
          END DO
       END DO
-
-      ! Surface boudary conditions
-      DO jj = 2, jpjm1 
+      DO jj = 2, jpjm1        ! Surface boudary conditions
          DO ji = fs_2, fs_jpim1   ! vector opt.
             zwi(ji,jj,1) = 0.
@@ -130,7 +119,6 @@
       !   The solution (the after velocity) is in ua
       !-----------------------------------------------------------------------
-
-      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
-      DO jk = 2, jpkm1
+      !
+      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -139,14 +127,9 @@
          END DO
       END DO
-
-      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
-      DO jj = 2, jpjm1   
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-!!! change les resultats (derniers digit, pas significativement + rapide 1* de moins)
-!!!         ua(ji,jj,1) = ub(ji,jj,1)  &
-!!!                      + p2dt * ( ua(ji,jj,1) + utau(ji,jj) / ( fse3u(ji,jj,1)*rau0 ) )
-            z2dtf = p2dt / ( fse3u(ji,jj,1)*rau0 )
-            ua(ji,jj,1) = ub(ji,jj,1)  &
-                         + p2dt *  ua(ji,jj,1) + z2dtf * utau(ji,jj)
+      !
+      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            ua(ji,jj,1) = ub(ji,jj,1)   &
+               &        + p2dt * (  ua(ji,jj,1) + 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(ji,jj,1) * rau0 )  )
          END DO
       END DO
@@ -159,7 +142,6 @@
          END DO
       END DO
-
-      ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk
-      DO jj = 2, jpjm1   
+      !
+      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
          DO ji = fs_2, fs_jpim1   ! vector opt.
             ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
@@ -173,7 +155,6 @@
          END DO
       END DO
-
-      ! Normalization to obtain the general momentum trend ua
-      DO jk = 1, jpkm1
+      !
+      DO jk = 1, jpkm1        !==  Normalization to obtain the general momentum trend ua  ==
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -192,11 +173,12 @@
       ! friction velocity in dyn_bfr using values from the previous timestep. There
       ! is no need to include these in the implicit calculation.
-      DO jk = 1, jpkm1
+      !
+      DO jk = 1, jpkm1        ! Matrix
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zcoef = -p2dt / fse3v(ji,jj,jk)
-               zzwi          = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  )
+               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
                zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk)
-               zzws       = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
+               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
                zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
                zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
@@ -204,7 +186,5 @@
          END DO
       END DO
-
-      ! Surface boudary conditions
-      DO jj = 2, jpjm1   
+      DO jj = 2, jpjm1        ! Surface boudary conditions
          DO ji = fs_2, fs_jpim1   ! vector opt.
             zwi(ji,jj,1) = 0.e0
@@ -223,12 +203,10 @@
       !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
       !
-      !   m is decomposed in the product of an upper and lower triangular
-      !   matrix
+      !   m is decomposed in the product of an upper and lower triangular matrix
       !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
       !   The solution (after velocity) is in 2d array va
       !-----------------------------------------------------------------------
-
-      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
-      DO jk = 2, jpkm1
+      !
+      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -237,14 +215,9 @@
          END DO
       END DO
-
-      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-!!! change les resultats (derniers digit, pas significativement + rapide 1* de moins)
-!!!         va(ji,jj,1) = vb(ji,jj,1)  &
-!!!                      + p2dt * ( va(ji,jj,1) + vtau(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) )
-            z2dtf = p2dt / ( fse3v(ji,jj,1)*rau0 )
-            va(ji,jj,1) = vb(ji,jj,1)  &
-                         + p2dt * va(ji,jj,1) + z2dtf * vtau(ji,jj)
+      !
+      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            va(ji,jj,1) = vb(ji,jj,1)   &
+               &        + p2dt * (  va(ji,jj,1) + 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(ji,jj,1) * rau0 )  )
          END DO
       END DO
@@ -257,7 +230,6 @@
          END DO
       END DO
-
-      ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk
-      DO jj = 2, jpjm1   
+      !
+      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
          DO ji = fs_2, fs_jpim1   ! vector opt.
             va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
@@ -271,7 +243,6 @@
          END DO
       END DO
-
-      ! Normalization to obtain the general momentum trend va
-      DO jk = 1, jpkm1
+      !
+      DO jk = 1, jpkm1     !==  Normalization to obtain the general momentum trend va  ==
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -280,5 +251,5 @@
          END DO
       END DO
-
+      !
    END SUBROUTINE dyn_zdf_imp
 
Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 1869)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 1870)
@@ -5,4 +5,5 @@
    !!==============================================================================
    !! History :  3.1  !  2009-02  (G. Madec, M. Leclair)  Original code
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA 
    !!----------------------------------------------------------------------
 
@@ -37,7 +38,6 @@
 #  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
-
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -53,7 +53,5 @@
       !!              and update the now vertical coordinate (lk_vvl=T).
       !!
-      !! ** Method  : -
-      !!
-      !!              - Using the incompressibility hypothesis, the vertical 
+      !! ** 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.
@@ -62,7 +60,8 @@
       !! ** action  :   ssha    : after sea surface height
       !!                wn      : now vertical velocity
-      !! if lk_vvl=T:   sshu_a, sshv_a, sshf_a  : after sea surface height
-      !!                          at u-, v-, f-point s
-      !!                hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points
+      !!                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.
       !!----------------------------------------------------------------------
       USE oce, ONLY :   z3d => ta   ! use ta as 3D workspace
@@ -78,4 +77,5 @@
 
       IF( kt == nit000 ) THEN
+         !
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity '
@@ -111,13 +111,13 @@
       ENDIF
 
-      !                                           !------------------------------!
-      IF( lk_vvl ) THEN                           !  Update Now Vertical coord.  !   (only in vvl case)
-      !                                           !------------------------------!
+      !                                           !------------------------------------------!
+      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
+            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
+            fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays
             fse3u (:,:,jk) = fse3u_n (:,:,jk)
             fse3v (:,:,jk) = fse3v_n (:,:,jk)
@@ -127,8 +127,8 @@
             fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
          END DO
-         !                                             ! now ocean depth (at u- and v-points)
-         hu(:,:) = hu_0(:,:) + sshu_n(:,:)
+         !
+         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)
+         !                                            ! now masked inverse of the ocean depth (at u- and v-points)
          hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) )
          hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) )
@@ -136,31 +136,28 @@
       ENDIF
 
-                         CALL div_cur( kt )            ! Horizontal divergence & Relative vorticity
-      IF( n_cla == 1 )   CALL div_cla( kt )            ! Cross Land Advection (Update Hor. divergence)
-
-      ! set time step size (Euler/Leapfrog)
-      z2dt = 2. * rdt 
+                         CALL div_cur( kt )           ! Horizontal divergence & Relative vorticity
+      IF( n_cla == 1 )   CALL div_cla( kt )           ! Cross Land Advection (Update Hor. divergence)
+
+      z2dt = 2. * rdt                                 ! set time step size (Euler/Leapfrog)
       IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt
 
-      zraur = 1. / rau0
-
       !                                           !------------------------------!
       !                                           !   After Sea Surface Height   !
       !                                           !------------------------------!
+      zraur = 0.5 / rau0
+      !
       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)
+      ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
 
 #if defined key_obc
-      IF ( Agrif_Root() ) THEN 
+      IF( Agrif_Root() ) THEN 
          ssha(:,:) = ssha(:,:) * obctmsk(:,:)
-         CALL lbc_lnk( ssha, 'T', 1. )  ! absolutly compulsory !! (jmm)
+         CALL lbc_lnk( ssha, 'T', 1. )                ! absolutly compulsory !! (jmm)
       ENDIF
 #endif
-
       !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
       IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
@@ -186,22 +183,25 @@
       !                                           !     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
+      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
+         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)     &
+            &                      - (  fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt
       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
-      IF( lk_diaar5 ) THEN
+      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
          z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
          DO jk = 1, jpk
             z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
          END DO
-         CALL iom_put( "w_masstr" , z3d                     )   !           vertical mass transport
-         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )   ! square of vertical mass transport
-      ENDIF
+         CALL iom_put( "w_masstr" , z3d                     )  
+         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
+      ENDIF
+      !
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
       !
    END SUBROUTINE ssh_wzv
@@ -216,14 +216,19 @@
       !!              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
+      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
+      !!              from the filter, see Leclair and Madec 2010) and swap :
+      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
+      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
+      !!                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               ! dummy loop indices
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj   ! dummy loop indices
+      REAL(wp) ::   zec      ! temporary scalare
       !!----------------------------------------------------------------------
 
@@ -234,9 +239,7 @@
       ENDIF
 
-      ! Time filter and swap of the ssh
-      ! -------------------------------
-
-      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)
@@ -247,5 +250,5 @@
             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) )
+                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) )
                   sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) )
                   sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) )
@@ -258,14 +261,16 @@
             END DO
          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)
             !
          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) )    
+                  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
                END DO
@@ -274,6 +279,15 @@
          !
       ENDIF
-      !
-      IF(ln_ctl)   CALL prt_ctl(tab2d_1=sshb    , clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
+
+      ! time filter with global conservation correction and array swap
+      sshbb(:,:) = sshb(:,:)
+      sshb (:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + zssha(:,:) )    &
+           &       - zfact  * 
+      sshn (:,:) = zssha(:,:)
+      empb (:,:) = emp(:,:)
+
+
+      !
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
       !
    END SUBROUTINE ssh_nxt
Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 1869)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 1870)
@@ -4,6 +4,7 @@
    !! Surface module :   variables defined in core memory 
    !!======================================================================
-   !! History :  3.0   !  2006-06  (G. Madec)  Original code
-   !!             -    !  2008-08  (G. Madec)  namsbc moved from sbcmod
+   !! History :  3.0  !  2006-06  (G. Madec)  Original code
+   !!             -   !  2008-08  (G. Madec)  namsbc moved from sbcmod
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
    !!----------------------------------------------------------------------
    USE par_oce          ! ocean parameters
@@ -25,9 +26,9 @@
    LOGICAL , PUBLIC ::   ln_ssr      = .FALSE.   !: Sea Surface restoring on SST and/or SSS      
    INTEGER , PUBLIC ::   nn_ice      = 0         !: flag on ice in the surface boundary condition (=0/1/2/3)
-   INTEGER , PUBLIC ::   nn_fwb      = 0         !: FreshWater Budget: 
+   INTEGER , PUBLIC ::   nn_fwb      = 0         !: FreshWater Budget control : 
    !                                             !:  = 0 unchecked 
    !                                             !:  = 1 global mean of e-p-r set to zero at each nn_fsbc time step
    !                                             !:  = 2 annual global mean of e-p-r set to zero
-   INTEGER , PUBLIC ::   nn_ico_cpl  = 0          !: ice-ocean coupling indicator
+   INTEGER , PUBLIC ::   nn_ico_cpl  = 0         !: ice-ocean coupling indicator
    !                                             !:  = 0   LIM-3 old case
    !                                             !:  = 1   stresses computed using now ocean velocity
@@ -37,25 +38,26 @@
    !!              Ocean Surface Boundary Condition fields
    !!----------------------------------------------------------------------
-   LOGICAL , PUBLIC ::   lhftau = .FALSE.              !: HF tau contribution: mean of stress module - module of the mean stress
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   utau      !: sea surface i-stress (ocean referential)     [N/m2]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   vtau      !: sea surface j-stress (ocean referential)     [N/m2]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   taum      !: module of sea surface stress (at T-point)    [N/m2] 
-   !! wndm is used only in PISCES to compute gases exchanges at the surface of the free ocean or in the leads in sea-ice parts
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   wndm      !: wind speed module at T-point (=|U10m-Uoce|)  [m/s] 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qsr       !: sea heat flux:     solar                     [W/m2]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qns       !: sea heat flux: non solar                     [W/m2]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qsr_tot   !: total     solar heat flux (over sea and ice) [W/m2]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qns_tot   !: total non solar heat flux (over sea and ice) [W/m2]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp       !: freshwater budget: volume flux               [Kg/m2/s]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emps      !: freshwater budget: concentration/dillution   [Kg/m2/s]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp_tot   !: total evaporation - (liquid + solid) precpitation over oce and ice
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tprecip   !: total precipitation           [Kg/m2/s]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sprecip   !: solid precipitation           [Kg/m2/s]
-!!$   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rrunoff       !: runoff
-!!$   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   calving       !: calving
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fr_i      !: ice fraction  (between 0 to 1)               -
+   LOGICAL , PUBLIC ::   lhftau = .FALSE.              !: HF tau used in TKE: mean(stress module) - module(mean stress)
+   !!                                   !!   now    ! before   !!
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   utau   , utau_b   !: sea surface i-stress (ocean referential)     [N/m2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   vtau   , vtau_b   !: sea surface j-stress (ocean referential)     [N/m2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   taum              !: module of sea surface stress (at T-point)    [N/m2] 
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   wndm              !: wind speed module at T-point (=|U10m-Uoce|)  [m/s]
+   !! wndm is used only in PISCES to compute surface gases exchanges in ice-free ocean or leads
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qsr    , qsr_b    !: sea heat flux:     solar                     [W/m2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qns    , qns_b    !: sea heat flux: non solar                     [W/m2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qsr_tot           !: total     solar heat flux (over sea and ice) [W/m2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qns_tot           !: total non solar heat flux (over sea and ice) [W/m2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp    , emp_b    !: freshwater budget: volume flux               [Kg/m2/s]
+   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]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tprecip           !: total precipitation                          [Kg/m2/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sprecip           !: solid precipitation                          [Kg/m2/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fr_i              !: ice fraction = 1 - lead fraction      (between 0 to 1)
 #if defined key_cpl_carbon_cycle
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   atm_co2   !: atmospheric pCO2                             [ppm]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   atm_co2           !: atmospheric pCO2                             [ppm]
 #endif
+!!$   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rrunoff        !: runoff
+!!$   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   calving        !: calving
 
    !!----------------------------------------------------------------------
@@ -70,8 +72,7 @@
 
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2006) 
-   !! $ Id: $
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
+   !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!======================================================================
-
 END MODULE sbc_oce
Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 1869)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 1870)
@@ -4,6 +4,7 @@
    !! Surface module :  provide to the ocean its surface boundary condition
    !!======================================================================
-   !! History :  3.0   !  07-2006  (G. Madec)  Original code
-   !!             -    !  08-2008  (S. Masson, E. .... ) coupled interface
+   !! History :  3.0  !  2006-07  (G. Madec)  Original code
+   !!            3.1  !  2008-08  (S. Masson, E. Maisonnave, G. Madec) coupled interface
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
    !!----------------------------------------------------------------------
 
@@ -49,5 +50,5 @@
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -86,5 +87,5 @@
 !!gm  IF( lk_sbc_cpl       ) THEN   ;   ln_cpl      = .TRUE.   ;   ELSE   ;   ln_cpl      = .FALSE.   ;   ENDIF
 
-      IF ( Agrif_Root() ) THEN
+      IF( Agrif_Root() ) THEN
         IF( lk_lim2 )            nn_ice      = 2
         IF( lk_lim3 )            nn_ice      = 3
@@ -179,6 +180,8 @@
       !!                CAUTION : never mask the surface stress field (tke sbc)
       !!
-      !! ** Action  : - set the ocean surface boundary condition, i.e.  
-      !!                utau, vtau, qns, qsr, emp, emps, qrp, erp
+      !! ** Action  : - set the ocean surface boundary condition at before and now 
+      !!                time step, i.e.  
+      !!                utau_b, vtau_b, qns_b, qsr_b, emp_n, emps_b, qrp_b, erp_b
+      !!                utau  , vtau  , qns  , qsr  , emp  , emps  , qrp  , erp
       !!              - updte the ice fraction : fr_i
       !!----------------------------------------------------------------------
@@ -186,17 +189,28 @@
       !!---------------------------------------------------------------------
 
-      CALL iom_setkt( kt + nn_fsbc - 1 )         !  in sbc, iom_put is called every nn_fsbc time step
-      !
-      ! ocean to sbc mean sea surface variables (ss._m)
-      ! ---------------------------------------
-      CALL sbc_ssm( kt )                         ! sea surface mean currents (at U- and V-points), 
-      !                                          ! temperature and salinity (at T-point) over nf_sbc time-step
-      !                                          ! (i.e. sst_m, sss_m, ssu_m, ssv_m)
-
-      ! sbc formulation
-      ! ---------------
-         
-      SELECT CASE( nsbc )                        ! Compute ocean surface boundary condition
-      !                                          ! (i.e. utau,vtau, qns, qsr, emp, emps)
+      !                                            ! ---------------------------------------- !
+      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
+         !                                         ! ---------------------------------------- !
+         utau_b(:,:) = utau(:,:)                         ! Swap the ocean forcing fields
+         utau_b(:,:) = utau(:,:)                         ! (except at nitOOO where before fields
+         qns_b (:,:) = qns (:,:)                         !  are set the end of the routine)
+         qsr_b (:,:) = qsr (:,:)
+         emp_b (:,:) = emp (:,:)
+         emps_b(:,:) = emps(:,:)
+      ENDIF
+
+      !                                            ! ---------------------------------------- !
+      !                                            !        forcing field computation         !
+      !                                            ! ---------------------------------------- !
+
+      CALL iom_setkt( kt + nn_fsbc - 1 )                 ! in sbc, iom_put is called every nn_fsbc time step
+      !
+      CALL sbc_ssm( kt )                                 ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
+      !                                                  ! averaged over nf_sbc time-step
+
+                                                   !==  sbc formulation  ==!
+                                                            
+      SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition
+      !                                                  ! (i.e. utau,vtau, qns, qsr, emp, emps)
       CASE(  0 )   ;   CALL sbc_gyre    ( kt )                    ! analytical formulation : GYRE configuration
       CASE(  1 )   ;   CALL sbc_ana     ( kt )                    ! analytical formulation : uniform sbc
@@ -214,6 +228,5 @@
       END SELECT
 
-      ! Misc. Options
-      ! -------------
+      !                                            !==  Misc. Options  ==!
 
 !!gm  IF( ln_dm2dc       )   CALL sbc_dcy( kt )                 ! Daily mean qsr distributed over the Diurnal Cycle
@@ -236,4 +249,47 @@
       !                                                         ! (update freshwater fluxes)
       !
+      !                                                ! ---------------------------------------- !
+      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
+         !                                             ! ---------------------------------------- !
+         IF( ln_rstart .AND.    &                               !* Restart: read in restart file
+            & iom_varid( numror, 'utau_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, 'utau_b', utau_b )   ! before i-stress  (U-point)
+            CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b )   ! before j-stress  (V-point)
+            CALL iom_get( numror, jpdom_autoglo, 'qns_b' , qns_b  )   ! before non solar heat flux (T-point)
+            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)
+            !
+         ELSE                                                   !* no restart: set from nit000 values
+            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields set to nit000'
+            utau_b(:,:) = utau(:,:) 
+            utau_b(:,:) = utau(:,:)
+            qns_b (:,:) = qns (:,:)
+            qsr_b (:,:) = qsr (:,:)
+            emp_b (:,:) = emp (:,:)
+            emps_b(:,:) = emps(:,:)
+         ENDIF
+      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, 'utau_b' , utau )    ! 
+         CALL iom_rstput( kt, nitrst, numrow, 'utau_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' , emp  )
+         !
+      ENDIF
+
+      !                                                ! ---------------------------------------- !
+      !                                                !        Outputs and control print         !
+      !                                                ! ---------------------------------------- !
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
          CALL iom_put( "emp"    , emp       )                   ! upward water flux
@@ -242,5 +298,5 @@
          CALL iom_put( "qns"    , qns       )                   ! solar heat flux    moved after the call to iom_setkt)
          CALL iom_put( "qsr"    ,       qsr )                   ! solar heat flux    moved after the call to iom_setkt)
-         IF(  nn_ice > 0 )   CALL iom_put( "ice_cover", fr_i )  ! ice fraction 
+         IF( nn_ice > 0 )   CALL iom_put( "ice_cover", fr_i )   ! ice fraction 
       ENDIF
       !
Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 1869)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 1870)
@@ -15,5 +15,5 @@
    !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf
    !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option
-   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  semi-implicit hpg with asselin filter
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  semi-implicit hpg with asselin filter + modified LF-RA
    !!----------------------------------------------------------------------
 
@@ -81,4 +81,6 @@
       !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
       !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
       !!----------------------------------------------------------------------
       USE oce, ONLY :    ztrdt => ua   ! use ua as 3D workspace   
@@ -91,5 +93,5 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( kt == nit000 ) THEN          !==  initialisation  ==!
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
@@ -99,32 +101,30 @@
          r1_2bcp = 1.e0 - 2.e0 * rbcp
       ENDIF
-
-      ! Update after tracer on domain lateral boundaries
+      !                                      ! set time step size (Euler/Leapfrog)
+      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt_t(:) =     rdttra(:)      ! at nit000             (Euler)
+      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt_t(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
+      ENDIF
+
+      !                                !==  Update after tracer on domain lateral boundaries  ==!
       ! 
-      CALL lbc_lnk( ta, 'T', 1. )      ! local domain boundaries  (T-point, unchanged sign)
+      CALL lbc_lnk( ta, 'T', 1. )            ! local domain boundaries  (T-point, unchanged sign)
       CALL lbc_lnk( sa, 'T', 1. )
       !
 #if defined key_obc
-      CALL obc_tra( kt )               ! OBC open boundaries
+      CALL obc_tra( kt )                     ! OBC open boundaries
 #endif
 #if defined key_bdy
-      CALL bdy_tra( kt )               ! BDY open boundaries
+      CALL bdy_tra( kt )                     ! BDY open boundaries
 #endif
 #if defined key_agrif
-      CALL Agrif_tra                   ! AGRIF zoom boundaries
+      CALL Agrif_tra                         ! AGRIF zoom boundaries
 #endif
  
-      ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt_t(:) =     rdttra(:)      ! at nit000             (Euler)
-      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt_t(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
-      ENDIF
-
-      ! trends computation initialisation
-      IF( l_trdtra ) THEN              ! store now fields before applying the Asselin filter
+      IF( l_trdtra ) THEN              ! trends computation: store now fields before applying the Asselin filter
          ztrdt(:,:,:) = tn(:,:,:)      
          ztrds(:,:,:) = sn(:,:,:)
       ENDIF
 
-      ! Leap-Frog + Asselin filter time stepping
+      !                                !==  modifed Leap-Frog + Asselin filter (modified LF-RA) scheme  ==!
       IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl( kt )      ! variable volume level (vvl)
       ELSE                  ;   CALL tra_nxt_fix( kt )      ! fixed    volume level
@@ -132,10 +132,8 @@
 
 #if defined key_agrif
-      ! Update tracer at AGRIF zoom boundaries
-      IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )      ! children only
+      IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )   ! Update tracer at AGRIF zoom boundaries (children only)
 #endif      
 
-      ! trends computation
-      IF( l_trdtra ) THEN              ! trend of the Asselin filter (tb filtered - tb)/dt     
+      IF( l_trdtra ) THEN              ! trends computation: trend of the Asselin filter (tb filtered - tb)/dt     
          DO jk = 1, jpkm1
             zfact = 1.e0 / r2dt_t(jk)             
@@ -176,9 +174,9 @@
       !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
-      !!
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   ztm, ztf     ! temporary scalars
-      REAL(wp) ::   zsm, zsf     !    -         -
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk      ! dummy loop indices
+      REAL(wp) ::   ztm, ztf, zac   ! temporary scalars
+      REAL(wp) ::   zsm, zsf        !    -         -
       !!----------------------------------------------------------------------
 
@@ -225,5 +223,5 @@
                sn(:,:,jk) = sa(:,:,jk)
             END DO
-         ELSE                                             ! general case (Leapfrog + Asselin filter
+         ELSE                                             ! general case (Leapfrog + Asselin filter)
             DO jk = 1, jpkm1
                DO jj = 1, jpj
@@ -239,4 +237,25 @@
             END DO
          ENDIF
+      ENDIF
+      !
+!!gm from Matthieu : unchecked
+      IF( neuler /= 0 .OR. kt /= nit000 ) THEN           ! remove the forcing from the Asselin filter
+         zac = atfp * rdttra(1)
+         tb(:,:,1) = tb(:,:,1) - zac * ( qns (:,:) - qns_b (:,:) )         ! non solar surface flux
+         sb(:,:,1) = sn(:,:,1) - zac * ( emps(:,:) - emps_b(:,:) )         ! surface salt flux
+         !
+         IF( ln_traqsr )                                                   ! solar penetrating flux
+            DO jk = 1, nksr
+                  DO jj = 1, jpj
+                     DO ji = 1, jpi
+                        IF( ln_sco ) THEN
+                           z_cor_qsr = etot3(ji,jj,jk) * ( qsr(ji,jj) - qsrb(ji,jj) )
+                        ELSEIF( ln_zps .OR. ln_zco ) THEN
+                           z_cor_qsr = ( qsr(ji,jj) - qsrb(ji,jj) ) *   &
+                                &   ( gdsr(jk) * tmask(ji,jj,jk) - gdsr(jk+1) * tmask(ji,jj,jk+1) )
+                        ENDIF
+                        zt = zt - zfact1 * z_cor_qsr
+                     END DO
+                  END DO
       ENDIF
       !
@@ -372,4 +391,41 @@
          ENDIF
       ENDIF
+!!gm from Matthieu : unchecked
+      !
+      IF( neuler /= 0 .OR. kt /= nit000 ) THEN           ! remove the forcing from the Asselin filter
+               IF( lk_vvl ) THEN                                        ! Varying levels
+                  DO jj = 1, jpj
+                     DO ji = 1, jpi
+                        ! - ML - modification for varaiance diagnostics
+                        zssh1 = tmask(ji,jj,jk) / ( fse3t(ji,jj,jk) + atfp * ( zfse3ta(ji,jj,jk) - 2*fse3t(ji,jj,jk) &
+                           &                                                 + fse3tb(ji,jj,jk) ) )
+                        zt = zssh1 * ( tn(ji,jj,jk) +  atfp * ( ta(ji,jj,jk) - 2 * tn(ji,jj,jk) + tb(ji,jj,jk) ) )
+                        zs = zssh1 * ( sn(ji,jj,jk) +  atfp * ( sa(ji,jj,jk) - 2 * sn(ji,jj,jk) + sb(ji,jj,jk) ) )
+                        ! filter correction for global conservation
+                        !------------------------------------------
+                        zfact1 = atfp * rdttra(1) * zssh1
+                        IF (jk == 1) THEN                        ! remove sbc trend from time filter
+                           zt = zt - zfact1 * ( qn(ji,jj) - qb(ji,jj) )
+!!???                           zs = zs - zfact1 * ( - emp( ji,jj) + empb(ji,jj) ) * zsrau * 35.5e0
+                        ENDIF
+                        ! remove qsr trend from time filter
+                        IF (jk <= nksr) THEN
+                           IF( ln_sco ) THEN
+                              z_cor_qsr = etot3(ji,jj,jk) * ( qsr(ji,jj) - qsrb(ji,jj) )
+                           ELSEIF( ln_zps .OR. ln_zco ) THEN
+                              z_cor_qsr = ( qsr(ji,jj) - qsrb(ji,jj) ) *   &
+                                   &   ( gdsr(jk) * tmask(ji,jj,jk) - gdsr(jk+1) * tmask(ji,jj,jk+1) )
+                           ENDIF
+                           zt = zt - zfact1 * z_cor_qsr
+                           IF (jk == nksr) qsrb(ji,jj) = qsr(ji,jj)
+                        ENDIF
+                        ! - ML - end of modification
+                        zssh1 = tmask(ji,jj,1) / zfse3ta(ji,jj,jk)
+                        tn(ji,jj,jk) = ta(ji,jj,jk) * zssh1
+                        sn(ji,jj,jk) = sa(ji,jj,jk) * zssh1
+                     END DO
+                  END DO
+      ENDIF
+!!gm end
       !
    END SUBROUTINE tra_nxt_vvl
Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 1869)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 1870)
@@ -10,4 +10,5 @@
    !!             -   !  2005-11  (G. Madec) zco, zps, sco coordinate
    !!            3.2  !  2009-04  (G. Madec & NEMO team) 
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
    !!----------------------------------------------------------------------
 
@@ -44,8 +45,8 @@
    REAL(wp), PUBLIC ::   rn_si2     = 61.8_wp   !: deepest depth of extinction (blue & 0.01 mg.m-3)     (RGB)
    
-   ! Module variables
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
-   INTEGER ::   nksr   ! levels below which the light cannot penetrate ( depth larger than 391 m)
-   REAL(wp), DIMENSION(3,61) ::   rkrgb   !: tabulated attenuation coefficients for RGB absorption
+   
+   INTEGER                   ::   nksr    ! levels below which the light cannot penetrate (depth larger than 391 m)
+   REAL(wp), DIMENSION(3,61) ::   rkrgb   ! tabulated attenuation coefficients for RGB absorption
 
    !! * Substitutions
@@ -53,5 +54,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -96,4 +97,5 @@
       REAL(wp) ::   zchl, zcoef, zsi0r   ! temporary scalars
       REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         -
+      REAL(wp) ::   zqsr                 !    -         -
       REAL(wp), DIMENSION(jpi,jpj)     ::   zekb, zekg, zekr            ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze0, ze1 , ze2, ze3, zea    ! 3D workspace
@@ -105,5 +107,5 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~'
          CALL tra_qsr_init
-         IF( .NOT.ln_traqsr )   RETURN
+         IF( .NOT.ln_traqsr )   RETURN             !!gm not authirized by Coding rules
       ENDIF
 
@@ -120,4 +122,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
+!!gm  how to stecify the mean of time step here : TOP versus OPA time stepping strategy not obvious
                   ta(ji,jj,jk) = ta(ji,jj,jk) + ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk) 
                END DO
@@ -135,9 +138,9 @@
             IF( nn_chldta ==1 ) THEN                             !*  Variable Chlorophyll
                !
-               CALL fld_read( kt, 1, sf_chl )                         ! Read Chl data and provides it at the current time step
+               CALL fld_read( kt, 1, sf_chl )                        ! Read Chl data and provides it at the current time step
                !         
 !CDIR COLLAPSE
 !CDIR NOVERRCHK
-               DO jj = 1, jpj                                         ! Separation in R-G-B depending of the surface Chl
+               DO jj = 1, jpj                                        ! Separation in R-G-B depending of the surface Chl
 !CDIR NOVERRCHK
                   DO ji = 1, jpi
@@ -149,19 +152,21 @@
                   END DO
                END DO
-               !
+               !                                                     ! surface value
                zsi0r = 1.e0 / rn_si0
-               zcoef  = ( 1. - rn_abs ) / 3.e0                        ! equi-partition in R-G-B
-               ze0(:,:,1) = rn_abs  * qsr(:,:)
-               ze1(:,:,1) = zcoef * qsr(:,:)
-               ze2(:,:,1) = zcoef * qsr(:,:)
-               ze3(:,:,1) = zcoef * qsr(:,:)
-               zea(:,:,1) =         qsr(:,:)
-               !
-               DO jk = 2, nksr+1
+               zcoef = ( 1. - rn_abs ) / 3.e0                           ! equi-partition in R-G-B
+!!gm bug !!! zqsr only a constant not an array
+               zqsr  = 0.5 * ( qsr_b(:,:) + qsr(:,:) )                  ! mean over 2 time steps
+               ze0(:,:,1) = rn_abs  * zqsr
+               ze1(:,:,1) = zcoef   * zqsr
+               ze2(:,:,1) = zcoef   * zqsr
+               ze3(:,:,1) = zcoef   * zqsr
+               zea(:,:,1) =           zqsr
+               !
+               DO jk = 2, nksr+1                                     ! deeper values
 !CDIR NOVERRCHK
                   DO jj = 1, jpj
 !CDIR NOVERRCHK   
                      DO ji = 1, jpi
-                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zsi0r     )
+                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zsi0r       )
                         zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) )
                         zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) )
@@ -175,6 +180,7 @@
                   END DO
                END DO
-               !
-               DO jk = 1, nksr                                        ! compute and add qsr trend to ta
+!!gm add here  etot3(:,:,jk) = ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk)
+               !
+               DO jk = 1, nksr                                       ! compute and add qsr trend to ta
                   ta(:,:,jk) = ta(:,:,jk) + ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk)
                END DO
@@ -184,17 +190,16 @@
             ELSE                                                 !*  Constant Chlorophyll
                DO jk = 1, nksr
-                  ta(:,:,jk) = ta(:,:,jk) + etot3(:,:,jk) * qsr(:,:)
+                  ta(:,:,jk) = ta(:,:,jk) + etot3(:,:,jk) * 0.5 * ( qsr_b(:,:) + qsr(:,:) )
                END DO
             ENDIF
-
+            !
          ENDIF
          !                                                ! ------------------------- !
          IF( ln_qsr_2bd ) THEN                            !  2 band light penetration !
             !                                             ! ------------------------- !
-            !
             DO jk = 1, nksr
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     ta(ji,jj,jk) = ta(ji,jj,jk) + etot3(ji,jj,jk) * qsr(ji,jj)
+                     ta(ji,jj,jk) = ta(ji,jj,jk) + etot3(ji,jj,jk) * 0.5 * ( qsr_b(:,:) + qsr(:,:) )
                   END DO
                END DO
Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 1869)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 1870)
@@ -4,7 +4,8 @@
    !! Ocean active tracers:  surface boundary condition
    !!==============================================================================
-   !! History :  8.2  !  98-10  (G. Madec, G. Roullet, M. Imbard)  Original code
-   !!            8.2  !  01-02  (D. Ludicone)  sea ice and free surface
-   !!            8.5  !  02-06  (G. Madec)  F90: Free form and module
+   !! History :  OPA  !  2998-10  (G. Madec, G. Roullet, M. Imbard)  Original code
+   !!            8.2  !  2001-02  (D. Ludicone)  sea ice and free surface
+   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
    !!----------------------------------------------------------------------
 
@@ -31,5 +32,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -46,7 +47,7 @@
       !!      and add it to the general trend of tracer equations.
       !!
-      !! ** Method :
-      !!      Following Roullet and Madec (2000), the air-sea flux can be divided 
-      !!      into three effects: (1) Fext, external forcing; 
+      !! ** Method  :   Following Roullet and Madec (2000), the air-sea flux 
+      !!      can be divided into three effects: 
+      !!      (1) Fext, external forcing; 
       !!      (2) Fwi, concentration/dilution effect due to water exchanged 
       !!         at the surface by evaporation, precipitations and runoff (E-P-R); 
@@ -58,5 +59,5 @@
       !!         solar part is added in traqsr.F routine.
       !!            ta = ta + q /(rau0 rcp e3t)  for k=1
-      !!            - salinity    : no salt flux
+      !!            - salinity    : only salt flux exchanged with sea-ice
       !!
       !!      The formulation for Fwb and Fwi vary according to the free 
@@ -123,31 +124,71 @@
       ENDIF
 
-      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration
+!!gm      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration
+      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
+         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
+         qsr(:,:) = 0.e0                     ! qsr set to zero
+         IF( kt == nit000 ) THEN             ! idem on before field at nit000
+            qns_b(:,:) = qns_b(:,:) + qsr_b(:,:)  
+            qsr_b(:,:) = 0.e0                     
+         ENDIF
+      ENDIF
 
       ! Concentration dillution effect on (t,s)
-      DO jj = 2, jpj
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-#if ! defined key_zco
-            zse3t = 1. / fse3t(ji,jj,1)
-#endif
-            IF( lk_vvl) THEN
-               zta = ro0cpr * qns(ji,jj) * zse3t &                   ! temperature : heat flux
-                &    - emp(ji,jj) * zsrau * tn(ji,jj,1)  * zse3t     ! & cooling/heating effet of EMP flux
-               zsa = 0.e0                                            ! No salinity concent./dilut. effect
-            ELSE
-               zta = ro0cpr * qns(ji,jj) * zse3t     ! temperature : heat flux
-               zsa = emps(ji,jj) * zsrau * sn(ji,jj,1)   * zse3t     ! salinity :  concent./dilut. effect
-            ENDIF
-            ta(ji,jj,1) = ta(ji,jj,1) + zta                          ! add the trend to the general tracer trend
-            sa(ji,jj,1) = sa(ji,jj,1) + zsa
+!!      DO jj = 2, jpj
+!!         DO ji = fs_2, fs_jpim1   ! vector opt.
+!!#if ! defined key_zco
+!!            zse3t = 1. / fse3t(ji,jj,1)
+!!#endif
+!!            IF( lk_vvl) THEN
+!!               zta = ro0cpr * qns(ji,jj) * zse3t &                   ! temperature : heat flux
+!!                &    - emp(ji,jj) * zsrau * tn(ji,jj,1)  * zse3t     ! & cooling/heating effet of EMP flux
+!!               zsa = 0.e0                                            ! No salinity concent./dilut. effect
+!!            ELSE
+!!               zta = ro0cpr * qns(ji,jj) * zse3t     ! temperature : heat flux
+!!               zsa = emps(ji,jj) * zsrau * sn(ji,jj,1)   * zse3t     ! salinity :  concent./dilut. effect
+!!            ENDIF
+!!            ta(ji,jj,1) = ta(ji,jj,1) + zta                          ! add the trend to the general tracer trend
+!!            sa(ji,jj,1) = sa(ji,jj,1) + zsa
+!!         END DO
+!!      END DO
+
+      !                             ! ---------------------- !
+      IF( lk_vvl ) THEN             !  Variable Volume case  !
+         !                          ! ---------------------- !
+         DO jj = 2, jpj
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zse3t = 0.5 / fse3t(ji,jj,1)
+               !                           ! temperature: heat flux + cooling/heating effet of EMP flux
+               ta(ji,jj,1) = ta(ji,jj,1) + (  ro0cpr * ( qns(ji,jj) + qns_b(ji,jj) )                & 
+                  &                         - zsrau  * ( emp(ji,jj) + emp_b(ji,jj) ) * tn(ji,jj,1)  ) * zse3t
+               !                           ! salinity: salt flux
+               sa(ji,jj,1) = sa(ji,jj,1) + ( emps(ji,jj) + emps_b(ji,jj) ) * zse3t
+
+!!gm BUG : in key_vvl emps must be modified to only include the salt flux due to with sea-ice freezing/melting
+!!gm       otherwise this flux will be missing  ==> modification required in limsbc,  limsbc_2 and CICE interface.
+
+            END DO
          END DO
-      END DO
+         !                          ! ---------------------- !
+      ELSE                          !  Constant Volume case  !
+         !                          ! ---------------------- !
+         DO jj = 2, jpj
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zse3t = 0.5 / fse3t(ji,jj,1)
+               !                           ! temperature: heat flux 
+               ta(ji,jj,1) = ta(ji,jj,1) + ro0cpr * ( qns (ji,jj) + qns_b (ji,jj) )               * zse3t 
+               !                           ! salinity: salt flux + concent./dilut. effect (both in emps)
+               sa(ji,jj,1) = sa(ji,jj,1) + zsrau  * ( emps(ji,jj) + emps_b(ji,jj) ) * sn(ji,jj,1) * zse3t
+            END DO
+         END DO
+         !
+      ENDIF
 
-      IF( l_trdtra ) THEN      ! save the sbc trends for diagnostic
+      IF( l_trdtra ) THEN        ! save the sbc trends for diagnostic
          ztrdt(:,:,:) = ta(:,:,:) - ztrdt(:,:,:)
          ztrds(:,:,:) = sa(:,:,:) - ztrds(:,:,:)
          CALL trd_mod(ztrdt, ztrds, jptra_trd_nsr, 'TRA', kt)
       ENDIF
-      !
+      !                          ! control print
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' sbc  - Ta: ', mask1=tmask,   &
          &                       tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
