Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 1974)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 1975)
@@ -244,8 +244,8 @@
                      zve3b = vb(ji,jj,jk) * ze3v_b
                      !
-                     zuf  = ( atfp  * ( zue3b  + zue3a  ) + atfp1 * zue3n  )   &
-                        & / ( atfp  * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk)
-                     zvf  = ( atfp  * ( zve3b  + zve3a  ) + atfp1 * zve3n  )   &
-                        & / ( atfp  * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk)
+                     zuf  = ( atfp * ( zue3b  + zue3a  ) + atfp1 * zue3n  )   &
+                        & / ( atfp * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk)
+                     zvf  = ( atfp * ( zve3b  + zve3a  ) + atfp1 * zve3n  )   &
+                        & / ( atfp * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk)
                      !
                      ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 1974)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 1975)
@@ -71,5 +71,5 @@
       INTEGER  ::   ji, jj, jk      ! dummy loop indices
       REAL(wp) ::   zcoefu, zcoefv, zcoeff      ! temporary scalars
-      REAL(wp) ::   z2dt, zraur     ! temporary scalars
+      REAL(wp) ::   z2dt, z1_2dt, z1_rau0       ! temporary scalars
       REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv       ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj) ::   z2d         ! 2D workspace
@@ -94,16 +94,28 @@
                   sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
                      &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
-                  sshf_b(ji,jj) = zcoeff * ( sshb(ji  ,jj) + sshb(ji  ,jj+1)                 &
-                     &                     + sshb(ji+1,jj) + sshb(ji+1,jj+1) )
                   sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
                      &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
                   sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
                      &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
-                  sshf_n(ji,jj) = zcoeff * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)                 &
-                     &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) )
                END DO
             END DO
             CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
             CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1      ! NO Vector Opt.
+                  sshf_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  ) )                     &
+                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
+                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
+               END DO
+            END DO
             CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. )
          ENDIF
@@ -145,5 +157,5 @@
       !                                           !   After Sea Surface Height   !
       !                                           !------------------------------!
-      zraur = 0.5 / rau0
+      z1_rau0 = 0.5 / rau0
       !
       zhdiv(:,:) = 0.e0
@@ -152,5 +164,7 @@
       END DO
       !                                                ! Sea surface elevation time stepping
-      ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
+      ! 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 defined key_obc
@@ -170,11 +184,18 @@
                   &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
                   &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
-               sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1)                                 & 
-                  &                                  * ( ssha(ji  ,jj) + ssha(ji  ,jj+1)                 &
-                  &                                    + ssha(ji+1,jj) + ssha(ji+1,jj+1) )
             END DO
          END DO
-         CALL lbc_lnk( sshu_a, 'U', 1. )               ! Boundaries conditions
+         ! Boundaries conditions
+         CALL lbc_lnk( sshu_a, 'U', 1. )
          CALL lbc_lnk( sshv_a, 'V', 1. )
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1      ! NO Vector Opt.
+               sshf_a(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
+                    &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
+                    &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_a(ji,jj  )     &
+                    &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_a(ji,jj+1) )
+            END DO
+         END DO
+         ! Boundaries conditions
          CALL lbc_lnk( sshf_a, 'F', 1. )
       ENDIF
@@ -183,7 +204,10 @@
       !                                           !     Now Vertical Velocity    !
       !                                           !------------------------------!
+      z1_2dt = 1.e0 / 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
+         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
+         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
+            &                      - (  fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
+            &                         * tmask(:,:,jk) * z1_2dt
       END DO
 
@@ -248,10 +272,9 @@
             sshf_n(:,:) = sshf_a(:,:)
          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) )
-                  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) )
-                  sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(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
                   sshu_n(ji,jj) = sshu_a(ji,jj)
@@ -260,4 +283,27 @@
                END DO
             END DO
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1      ! NO Vector Opt.
+                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
+                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
+                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
+                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
+                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
+                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
+               END DO
+            END DO
+            ! Boundaries conditions
+            CALL lbc_lnk( sshu_b, 'U', 1. )
+            CALL lbc_lnk( sshv_b, 'V', 1. )
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1      ! NO Vector Opt.
+                  sshf_b(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
+                     &                 / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
+                     &                 * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_b(ji,jj  )     &
+                     &                   + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_b(ji,jj+1) )
+               END DO
+            END DO
+            ! Boundaries conditions
+            CALL lbc_lnk( sshf_b, 'F', 1. )
          ENDIF
          !                    !--------------------------!
@@ -268,9 +314,7 @@
             !
          ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
-            zec = atfp * rdt / rau0
             DO jj = 1, jpj
                DO ji = 1, jpi                                ! before <-- now filtered
-                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   &
-                     &                      - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
+                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
                   sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after
                END DO
@@ -279,13 +323,4 @@
          !
       ENDIF
-
-      ! 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 )
Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/IOM/restart.F90	(revision 1974)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/IOM/restart.F90	(revision 1975)
@@ -217,4 +217,5 @@
          hdivb(:,:,:) = hdivn(:,:,:)
          sshb (:,:)   = sshn (:,:)
+         ! - ML - sshbnc 
       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 1974)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 1975)
@@ -52,4 +52,12 @@
    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
+   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                   "
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_trd_sc_n      !: sbc salt content trend now                   [psu/m/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_trd_sc_b      !:  "   "      "      "   before                   "
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   qsr_trd_hc_n  !: heat content trend due to qsr flux now       [K/m/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   qsr_trd_hc_b  !:  "      "      "    "  "   "   "   before       "
+   ! - ML - end
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tprecip           !: total precipitation                          [Kg/m2/s]
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sprecip           !: solid precipitation                          [Kg/m2/s]
Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 1974)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 1975)
@@ -193,9 +193,14 @@
          !                                         ! ---------------------------------------- !
          utau_b(:,:) = utau(:,:)                         ! Swap the ocean forcing fields
-         utau_b(:,:) = utau(:,:)                         ! (except at nitOOO where before fields
+         vtau_b(:,:) = vtau(:,:)                         ! (except at nit000 where before fields
          qns_b (:,:) = qns (:,:)                         !  are set 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(:,:)
+         
       ENDIF
 
@@ -256,14 +261,20 @@
             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, 'vtau_b', vtau_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)
+            ! - 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
             !
          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(:,:)
+            vtau_b(:,:) = vtau(:,:)
             qns_b (:,:) = qns (:,:)
             qsr_b (:,:) = qsr (:,:)
@@ -280,10 +291,16 @@
             &                    '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  )
+         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 )
+         ! - 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
          !
       ENDIF
Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 1974)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 1975)
@@ -25,4 +25,5 @@
    USE oce             ! ocean dynamics and tracers variables
    USE dom_oce         ! ocean space and time domain variables 
+   USE sbc_oce         ! surface boundary condition: ocean
    USE zdf_oce         ! ???
    USE domvvl          ! variable volume
@@ -37,4 +38,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl          ! Print control
+   USE traqsr          ! penetrative solar radiation (needed for nksr)
    USE agrif_opa_update
    USE agrif_opa_interp
@@ -45,5 +47,5 @@
    PUBLIC   tra_nxt    ! routine called by step.F90
 
-   REAL(wp)                 ::   rbcp, r1_2bcp   ! Brown & Campana parameters for semi-implicit hpg
+   REAL(wp)                 ::   rbcp            ! Brown & Campana parameters for semi-implicit hpg
    REAL(wp), DIMENSION(jpk) ::   r2dt_t          ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler)
 
@@ -98,6 +100,5 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~'
          !
-         rbcp    = 0.25 * (1 + atfp) * (1 + atfp * atfp)       ! Brown & Campana parameter for semi-implicit hpg
-         r1_2bcp = 1.e0 - 2.e0 * rbcp
+         rbcp    = 0.25 * (1. + atfp) * (1. + atfp) * ( 1. - atfp)       ! Brown & Campana parameter for semi-implicit hpg
       ENDIF
       !                                      ! set time step size (Euler/Leapfrog)
@@ -164,7 +165,7 @@
       !!              - swap tracer fields to prepare the next time_step.
       !!                This can be summurized for temperature as:
-      !!             ztm = (rbcp*ta+(2-rbcp)*tn+rbcp*tb)       ln_dynhpg_imp = T
+      !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T
       !!             ztm = 0                                   otherwise
-      !!                   with rbcp=(1+atfp)*(1+atfp*atfp)/4
+      !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp)
       !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
       !!             tn  = ta 
@@ -177,6 +178,6 @@
       !!
       INTEGER  ::   ji, jj, jk      ! dummy loop indices
-      REAL(wp) ::   ztm, ztf, zac   ! temporary scalars
-      REAL(wp) ::   zsm, zsf        !    -         -
+      REAL(wp) ::   zt_m, zs_m      ! temporary scalars
+      REAL(wp) ::   ztn, zsn        !    -         -
       !!----------------------------------------------------------------------
 
@@ -187,75 +188,35 @@
       ENDIF
       !
-      !                                              ! ----------------------- !
-      IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case !
-         !                                           ! ----------------------- !
-         !
-         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
-            DO jk = 1, jpkm1                              ! (only swap)
-               tn(:,:,jk) = ta(:,:,jk)                                                          ! tb <-- tn
-               sn(:,:,jk) = sa(:,:,jk)
-            END DO
-         ELSE                                             ! general case (Leapfrog + Asselin filter)
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     ztm = rbcp * ( ta(ji,jj,jk) + tb(ji,jj,jk) ) + r1_2bcp * tn(ji,jj,jk)      ! mean t
-                     zsm = rbcp * ( sa(ji,jj,jk) + sb(ji,jj,jk) ) + r1_2bcp * sn(ji,jj,jk) 
-                     ztf = atfp * ( ta(ji,jj,jk) + tb(ji,jj,jk)   - 2.      * tn(ji,jj,jk) )    ! Asselin filter on t 
-                     zsf = atfp * ( sa(ji,jj,jk) + sb(ji,jj,jk)   - 2.      * sn(ji,jj,jk) )
-                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                          ! tb <-- filtered tn 
-                     sb(ji,jj,jk) = sn(ji,jj,jk) + zsf
-                     tn(ji,jj,jk) = ta(ji,jj,jk)                                                ! tn <-- ta
-                     sn(ji,jj,jk) = sa(ji,jj,jk)
-                     ta(ji,jj,jk) = ztm                                                         ! ta <-- mean t
-                     sa(ji,jj,jk) = zsm
-                  END DO
+      IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
+         DO jk = 1, jpkm1                                                 ! (only swap)
+            tn(:,:,jk) = ta(:,:,jk)                                       ! tn <-- ta
+            sn(:,:,jk) = sa(:,:,jk)                                       ! sn <-- sa
+         END DO
+      ELSE                                             ! General case (Leapfrog + Asselin filter)
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF( ln_dynhpg_imp ) THEN                  ! implicit hpg: keep tn, sn in memory
+                     ztn = tn(ji,jj,jk)
+                     zsn = sn(ji,jj,jk)
+                  ENDIF
+                  !                                         ! 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)
+                  !
+                  !                                         ! 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
+                  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
+                  ENDIF
                END DO
             END DO
-         ENDIF
-         !                                           ! ----------------------- !
-      ELSE                                           !    explicit hpg case    !
-         !                                           ! ----------------------- !
-         !
-         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
-            DO jk = 1, jpkm1
-               tn(:,:,jk) = ta(:,:,jk)                                                    ! tn <-- ta
-               sn(:,:,jk) = sa(:,:,jk)
-            END DO
-         ELSE                                             ! general case (Leapfrog + Asselin filter)
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )      ! Asselin filter on t 
-                     zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) )
-                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                    ! tb <-- filtered tn 
-                     sb(ji,jj,jk) = sn(ji,jj,jk) + zsf
-                     tn(ji,jj,jk) = ta(ji,jj,jk)                                          ! tn <-- ta
-                     sn(ji,jj,jk) = sa(ji,jj,jk)
-                  END DO
-               END DO
-            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
+         END DO
       ENDIF
       !
@@ -276,6 +237,6 @@
       !!              - swap tracer fields to prepare the next time_step.
       !!                This can be summurized for temperature as:
-      !!             ztm = ( rbcp*e3t_a*ta + (2-rbtp)*e3t_n*tn + rbcp*e3t_b*tb )   ln_dynhpg_imp = T
-      !!                  /( rbcp*e3t_a    + (2-rbcp)*e3t_n    + rbcp*e3t_b    )   
+      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
+      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
       !!             ztm = 0                                                       otherwise
       !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
@@ -287,13 +248,15 @@
       !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      INTEGER, INTENT(in) ::   kt                  ! ocean time-step index
       !!     
-      INTEGER  ::   ji, jj, jk             ! dummy loop indices
-      REAL(wp) ::   ztm , ztc_f , ztf , ztca, ztcn, ztcb   ! temporary scalar
-      REAL(wp) ::   zsm , zsc_f , zsf , zsca, zscn, zscb   !    -         -
-      REAL(wp) ::   ze3mr, ze3fr                           !    -         -
-      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f         !    -         -
-      !!----------------------------------------------------------------------
-
+      INTEGER  ::   ji, jj, jk                     ! dummy loop indices
+      REAL     ::   ze3t_a, ze3t_n, ze3t_b         ! temporary scalars
+      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     ::   zfact1, zfact2                 !    -         -
+      !!----------------------------------------------------------------------
+      
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
@@ -302,130 +265,59 @@
       ENDIF
 
-      !                                              ! ----------------------- !
-      IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case !  (mean tracer computation)
-         !                                           ! ----------------------- !
-         !
-         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
-            DO jk = 1, jpkm1                                   ! (only swap)
-               tn(:,:,jk) = ta(:,:,jk)                         ! tn <-- ta
-               sn(:,:,jk) = sa(:,:,jk)
-            END DO
-            !                                             ! general case (Leapfrog + Asselin filter)
-         ELSE                                             ! apply filter on thickness weighted tracer and swap
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     ze3t_b = fse3t_b(ji,jj,jk)
-                     ze3t_n = fse3t_n(ji,jj,jk)
-                     ze3t_a = fse3t_a(ji,jj,jk)
-                     !                                         ! tracer content at Before, now and after
-                     ztcb = tb(ji,jj,jk) *  ze3t_b   ;   zscb = sb(ji,jj,jk) * ze3t_b
-                     ztcn = tn(ji,jj,jk) *  ze3t_n   ;   zscn = sn(ji,jj,jk) * ze3t_n
-                     ztca = ta(ji,jj,jk) *  ze3t_a   ;   zsca = sa(ji,jj,jk) * ze3t_a
-                     !
-                     !                                         ! Asselin filter on thickness and tracer content
-                     ze3t_f = atfp * ( ze3t_a + ze3t_b - 2.* ze3t_n )
-                     ztc_f  = atfp * ( ztca   + ztcb   - 2.* ztcn   ) 
-                     zsc_f  = atfp * ( zsca   + zscb   - 2.* zscn   ) 
-                     !
-                     !                                         ! filtered tracer including the correction 
-                     ze3fr = 1.e0  / ( ze3t_n + ze3t_f )
-                     ztf   = ze3fr * ( ztcn   + ztc_f  )
-                     zsf   = ze3fr * ( zscn   + zsc_f  )
-                     !                                         ! mean thickness and tracer
-                     ze3mr = 1.e0  / (  rbcp * ( ze3t_a + ze3t_b ) + r1_2bcp * ze3t_n )
-                     ztm   = ze3mr * (  rbcp * ( ztca   + ztcb   ) + r1_2bcp * ztcn   )
-                     zsm   = ze3mr * (  rbcp * ( zsca   + zscb   ) + r1_2bcp * zscn   )
-!!gm mean e3t have to be saved and used in dynhpg  or it can be recomputed in dynhpg !!
-!!gm                 e3t_m(ji,jj,jk) = 1.e0 / ze3mr
-                     !                                         ! swap of arrays
-                     tb(ji,jj,jk) = ztf                             ! tb <-- tn + filter
-                     sb(ji,jj,jk) = zsf
-                     tn(ji,jj,jk) = ta(ji,jj,jk)                    ! tn <-- ta
-                     sn(ji,jj,jk) = sa(ji,jj,jk)
-                     ta(ji,jj,jk) = ztm                             ! ta <-- mean t
-                     sa(ji,jj,jk) = zsm
-                  END DO
+      IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
+         !                                                  ! (only swap)
+         DO jk = 1, jpkm1                                   ! tn <-- ta
+            tn(:,:,jk) = ta(:,:,jk)                         ! sn <-- sa
+            sn(:,:,jk) = sa(:,:,jk)
+         END DO
+         !                                             ! General case (Leapfrog + Asselin filter)
+      ELSE                                             ! apply filter on thickness weighted tracer and swap
+         DO jk = 1, jpkm1
+            zfact1 = atfp * r2dt_t(jk)
+            zfact2 = zfact1 / rau0
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  !                                         ! scale factors at Before, now and after
+                  ze3t_b = fse3t_b(ji,jj,jk)
+                  ze3t_n = fse3t_n(ji,jj,jk)
+                  ze3t_a = fse3t_a(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
+                  ztc_n  = tn(ji,jj,jk) * ze3t_n   ;   zsc_n = sn(ji,jj,jk) * ze3t_n
+                  ztc_a  = ta(ji,jj,jk) * ze3t_a   ;   zsc_a = sa(ji,jj,jk) * ze3t_a
+                  !
+                  !                                         ! Time laplacian on thickness and tracer content
+                  !                                         ! used for both Asselin and Brown & Campana filters
+                  ze3t_m = ze3t_a - 2. * ze3t_n + ze3t_b
+                  ztc_m  = ztc_a  - 2. * ztc_n  + ztc_b
+                  zsc_m  = zsc_a  - 2. * zsc_n  + zsc_b
+                  !                                         ! Asselin Filter + correction
+                  ze3t_f = ze3t_n + atfp * ze3t_m
+                  ztc_f  = ztc_n  + atfp * ztc_m
+                  zsc_f  = zsc_n  + atfp * zsc_m
+
+                  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)    )
+                  ENDIF
+                  IF( ln_traqsr .AND. ( jk .LE. nksr ) ) THEN
+                     ztc_f  = ztc_f  - zfact1 * ( qsr_trd_hc_n(ji,jj,jk) - qsr_trd_hc_b(ji,jj,jk) )
+                  ENDIF
+                  !                                         ! swap of arrays
+                  ze3t_f = 1.e0 / ze3t_f
+                  tb(ji,jj,jk) = ztc_f * ze3t_f                           ! tb <-- tn filtered
+                  sb(ji,jj,jk) = zsc_f * ze3t_f                           ! 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
+                     ze3t_m = 1.e0 / ( ze3t_n + rbcp * ze3t_m )
+                     ta(ji,jj,jk) =  ( ztc_n  + rbcp * ztc_m  ) * ze3t_m  ! ta <-- Brown & Campana average
+                     sa(ji,jj,jk) =  ( zsc_n  + rbcp * zsc_m  ) * ze3t_m  ! sa <-- Brown & Campana average
+                  ENDIF
                END DO
             END DO
-         ENDIF
-         !                                           ! ----------------------- !
-      ELSE                                           !    explicit hpg case    !  (NO mean tracer computation)
-         !                                           ! ----------------------- !
-         !
-         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! case of Euler time-stepping at first time-step
-            DO jk = 1, jpkm1                                   ! (only swap)
-               tn(:,:,jk) = ta(:,:,jk)                         ! tn <-- ta
-               sn(:,:,jk) = sa(:,:,jk)
-            END DO
-            !                                             ! general case (Leapfrog + Asselin filter)
-         ELSE                                             ! apply filter on thickness weighted tracer and swap
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     ze3t_b = fse3t_b(ji,jj,jk)
-                     ze3t_n = fse3t_n(ji,jj,jk)
-                     ze3t_a = fse3t_a(ji,jj,jk)
-                     !                                         ! tracer content at Before, now and after
-                     ztcb = tb(ji,jj,jk) *  ze3t_b   ;   zscb = sb(ji,jj,jk) * ze3t_b
-                     ztcn = tn(ji,jj,jk) *  ze3t_n   ;   zscn = sn(ji,jj,jk) * ze3t_n
-                     ztca = ta(ji,jj,jk) *  ze3t_a   ;   zsca = sa(ji,jj,jk) * ze3t_a
-                     !
-                     !                                         ! Asselin filter on thickness and tracer content
-                     ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b )
-                     ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   ) 
-                     zsc_f  = atfp * ( zsca   - 2.* zscn   + zscb   ) 
-                     !
-                     !                                         ! filtered tracer including the correction 
-                     ze3fr = 1.e0 / ( ze3t_n + ze3t_f )
-                     ztf   =  ( ztcn  + ztc_f ) * ze3fr
-                     zsf   =  ( zscn  + zsc_f ) * ze3fr
-                     !                                         ! swap of arrays
-                     tb(ji,jj,jk) = ztf                             ! tb <-- tn filtered
-                     sb(ji,jj,jk) = zsf
-                     tn(ji,jj,jk) = ta(ji,jj,jk)                    ! tn <-- ta
-                     sn(ji,jj,jk) = sa(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-         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 DO
+      ENDIF
       !
    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 1974)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 1975)
@@ -47,5 +47,5 @@
    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)
+   INTEGER , PUBLIC          ::   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
 
@@ -97,5 +97,4 @@
       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
@@ -123,5 +122,5 @@
                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) 
+                  qsr_trd_hc_n(ji,jj,jk) = ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) )
                END DO
             END DO
@@ -155,11 +154,10 @@
                zsi0r = 1.e0 / rn_si0
                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
+
+               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                                     ! deeper values
@@ -183,5 +181,5 @@
                !
                DO jk = 1, nksr                                       ! compute and add qsr trend to ta
-                  ta(:,:,jk) = ta(:,:,jk) + ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk)
+                  qsr_trd_hc_n(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) )
                END DO
                zea(:,:,nksr+1:jpk) = 0.e0     ! below 400m set to zero
@@ -190,5 +188,5 @@
             ELSE                                                 !*  Constant Chlorophyll
                DO jk = 1, nksr
-                  ta(:,:,jk) = ta(:,:,jk) + etot3(:,:,jk) * 0.5 * ( qsr_b(:,:) + qsr(:,:) )
+                  qsr_trd_hc_n(:,:,jk) = etot3(:,:,jk) * qsr(:,:)
                END DO
             ENDIF
@@ -201,5 +199,5 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     ta(ji,jj,jk) = ta(ji,jj,jk) + etot3(ji,jj,jk) * 0.5 * ( qsr_b(:,:) + qsr(:,:) )
+                     qsr_trd_hc_n(ji,jj,jk) = etot3(ji,jj,jk) * qsr(ji,jj)
                   END DO
                END DO
@@ -208,4 +206,23 @@
          ENDIF
          !
+      ENDIF
+
+      ! Add qsr trend to ta in all cases
+      IF( neuler == 0 .AND. kt == nit000 ) THEN
+         DO jk = 1, nksr
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ta(ji,jj,jk) = ta(ji,jj,jk) + qsr_trd_hc_n(ji,jj,jk) / fse3t(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ELSE
+         DO jk = 1, nksr
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ta(ji,jj,jk) = ta(ji,jj,jk) + 0.5 * ( qsr_trd_hc_b(ji,jj,jk) + qsr_trd_hc_n(ji,jj,jk) ) / fse3t(ji,jj,jk)
+               END DO
+            END DO
+         END DO
       ENDIF
 
@@ -382,5 +399,5 @@
                !
                DO jk = 1, nksr
-                  etot3(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk)
+                  etot3(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) )
                END DO
                etot3(:,:,nksr+1:jpk) = 0.e0                ! below 400m set to zero
@@ -404,5 +421,5 @@
                      zc0 = rn_abs * EXP( -fsdepw(ji,jj,jk  )*zsi0r ) + (1.-rn_abs) * EXP( -fsdepw(ji,jj,jk  )*zsi1r )
                      zc1 = rn_abs * EXP( -fsdepw(ji,jj,jk+1)*zsi0r ) + (1.-rn_abs) * EXP( -fsdepw(ji,jj,jk+1)*zsi1r )
-                     etot3(ji,jj,jk) = ro0cpr * (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) / fse3t(ji,jj,jk)
+                     etot3(ji,jj,jk) = ro0cpr * (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  )
                   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 1974)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 1975)
@@ -104,6 +104,6 @@
       INTEGER, INTENT(in) ::   kt     ! ocean time-step index
       !!
-      INTEGER  ::   ji, jj                   ! dummy loop indices
-      REAL(wp) ::   zta, zsa, zsrau, zse3t   ! temporary scalars
+      INTEGER  ::   ji, jj            ! dummy loop indices
+      REAL(wp) ::   zsrau, zse3t      ! temporary scalars
       !!----------------------------------------------------------------------
 
@@ -129,66 +129,78 @@
          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
-!!         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
+!!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
+         IF ( neuler == 0 .AND. kt == nit000 ) THEN
+            DO jj = 2, jpj
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ! temperature : heat flux + cooling/heating effet of EMP flux
+                  sbc_trd_hc_n(ji,jj) = ro0cpr * qns(ji,jj) - zsrau * emp(ji,jj) * tn(ji,jj,1)
+#if ! defined key_zco
+                  zse3t = 1. / fse3t(ji,jj,1)
+#endif
+                  ta(ji,jj,1) = ta(ji,jj,1) + zse3t * sbc_trd_hc_n(ji,jj)
+                END DO
+            END DO
+         ELSE
+            DO jj = 2, jpj
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ! temperature : heat flux + cooling/heating effet of EMP flux
+                  sbc_trd_hc_n(ji,jj) = ro0cpr * qns(ji,jj) - zsrau * emp(ji,jj) * tn(ji,jj,1)
+#if ! defined key_zco
+                  zse3t = 1. / fse3t(ji,jj,1)
+#endif
+                  ta(ji,jj,1) = ta(ji,jj,1) + 0.5 * ( sbc_trd_hc_b(ji,jj) + sbc_trd_hc_n(ji,jj) ) * zse3t
+               END DO
+            END DO
+         ENDIF
          !                          ! ---------------------- !
       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
+         IF ( neuler == 0 .AND. kt == nit000 ) THEN
+            DO jj = 2, jpj
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ! temperature : heat flux
+                  sbc_trd_hc_n(ji,jj) = ro0cpr * qns(ji,jj)
+                  ! salinity    : salt flux + concent./dilut. effect (both in emps)
+                  sbc_trd_sc_n(ji,jj) = zsrau * emps(ji,jj) * sn(ji,jj,1)
+#if ! defined key_zco
+                  zse3t = 1. / fse3t(ji,jj,1)
+#endif
+                  ta(ji,jj,1) = ta(ji,jj,1) + sbc_trd_hc_n(ji,jj) * zse3t
+                  sa(ji,jj,1) = sa(ji,jj,1) + sbc_trd_sc_n(ji,jj) * zse3t
+               END DO
+            END DO
+         ELSE
+            DO jj = 2, jpj
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ! temperature : heat flux
+                  sbc_trd_hc_n(ji,jj) = ro0cpr * qns(ji,jj)
+                  ! salinity    : salt flux + concent./dilut. effect (both in emps)
+                  sbc_trd_sc_n(ji,jj) = zsrau * emps(ji,jj) * sn(ji,jj,1)
+#if ! defined key_zco
+                  zse3t = 1. / fse3t(ji,jj,1)
+#endif
+                  ! temperature : heat flux
+                  ta(ji,jj,1) = ta(ji,jj,1) + 0.5 * ( sbc_trd_hc_b(ji,jj) + sbc_trd_hc_n(ji,jj) ) * zse3t 
+                  sa(ji,jj,1) = sa(ji,jj,1) + 0.5 * ( sbc_trd_sc_b(ji,jj) + sbc_trd_sc_n(ji,jj) ) * zse3t
+               END DO
+            END DO
+         ENDIF
          !
       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
+      !                             ! 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' )
Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/step.F90	(revision 1974)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/step.F90	(revision 1975)
@@ -116,4 +116,6 @@
    USE restart         ! ocean restart                    (rst_wri routine)
    USE prtctl          ! Print control                    (prt_ctl routine)
+   ! - ML - 
+   USE diatrb          ! global tracer conservation       (dia_trb  routine)
 
 #if defined key_agrif
@@ -317,4 +319,5 @@
 
                                CALL ssh_nxt( kstp )         ! sea surface height at next time step
+                               CALL dia_trb( kstp )         ! - ML - global conservation diagnostics
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
