Index: /branches/dev_r2091_ZTC/NEMO/LIM_SRC_2/limthd_2.F90
===================================================================
--- /branches/dev_r2091_ZTC/NEMO/LIM_SRC_2/limthd_2.F90	(revision 2115)
+++ /branches/dev_r2091_ZTC/NEMO/LIM_SRC_2/limthd_2.F90	(revision 2116)
@@ -241,5 +241,5 @@
             
             !  energy needed to bring ocean surface layer until its freezing
-            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1)   &
+            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj)   &
                 &          * ( tfu(ji,jj) - sst_m(ji,jj) ) * ( 1 - zinda )
             
Index: /branches/dev_r2091_ZTC/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- /branches/dev_r2091_ZTC/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 2116)
+++ /branches/dev_r2091_ZTC/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 2116)
@@ -0,0 +1,249 @@
+MODULE diahsb
+   !!======================================================================
+   !!                       ***  MODULE  diahsb  ***
+   !! Ocean diagnostics: Heat salt and volume budgets
+   !!======================================================================
+   !! History :  NEMO 3.3  !  2010-09  (M. Leclair)  Original code 
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE sbc_oce         ! surface thermohaline fluxes
+   USE in_out_manager  ! I/O manager
+   USE domvvl          ! vertical scale factors
+   USE traqsr          ! penetrative solar radiation
+   USE lib_mpp         ! distributed memory computing library
+   USE trabbc          ! bottom boundary condition
+   USE bdy_par         ! (for lk_bdy)
+   USE obc_par         ! (for lk_obc)
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC dia_hsb        ! routine called by step.F90
+   PUBLIC dia_hsb_init   ! routine called by opa.F90
+
+   !! * Module variables
+   INTEGER				   ::   num			          !
+   REAL(dp)			           ::   surf_tot  , vol_tot               !
+   REAL(dp)			           ::   frc_t	  , frc_s     , frc_v     ! global forcing trends
+   REAL(dp)				   ::   fact1			          ! conversion factors
+   REAL(dp)				   ::   fact21	  , fact22	          !     -         -
+   REAL(dp)              		   ::   fact31	  , fact32	          !     -         -
+   REAL(dp), DIMENSION(:,:)  , ALLOCATABLE ::   surf	  , ssh_ini		  !
+   REAL(dp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini   !
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
+   !! $Id: trasbc.F90 2091 2010-09-14 20:29:38Z mlelod $
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE dia_hsb( kt )
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_hsb  ***
+      !!     
+      !! ** Purpose: Compute the ocean global heat content, salt content and volume
+      !!             non conservation
+      !!	
+      !! ** Method : - Compute the deviation of heat content, salt content and volume
+      !!	       at the current time step from their values at nit000
+      !!	     - Compute the contribution of forcing and remove it from these
+      !!                deviations
+      !! ** Action : Write the results in the 'heat_salt_volume_budgets.txt' ASCII file
+      !!---------------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !!
+      INTEGER	 ::   jk		          ! dummy loop indice
+      REAL(dp)	 ::   zdiff_hc	  , zdiff_sc	  ! heat and salt content variations
+      REAL(dp)	 ::   zdiff_v1	  , zdiff_v2	  ! volume variation
+      REAL(dp)	 ::   z1_rau0           	  ! local scalars
+      REAL(dp)	 ::   zdeltat			  !    -     -
+      REAL(dp)	 ::   z_frc_trd_t , z_frc_trd_s   !    -     -
+      REAL(dp)	 ::   z_frc_trd_v	          !    -     -
+      !!---------------------------------------------------------------------------
+
+      ! ------------------------- !
+      ! 1 - Trends due to forcing !
+      ! ------------------------- !
+      z1_rau0 = 1.e0 / rau0
+      z_frc_trd_v = z1_rau0 * SUM(    - emp(:,:) * surf(:,:) )     ! volume fluxes
+      z_frc_trd_t =           SUM( sbc_hc_n(:,:) * surf(:,:) )     ! heat fluxes
+      z_frc_trd_s =           SUM( sbc_sc_n(:,:) * surf(:,:) )     ! salt fluxes
+      ! Add penetrative solar radiation
+      IF( ln_traqsr ) z_frc_trd_t = z_frc_trd_t + ro0cpr * SUM( qsr     (:,:) * surf(:,:) )
+      ! Add geothermal heat flux
+      IF( lk_trabbc ) z_frc_trd_t = z_frc_trd_t + ro0cpr * SUM( qgh_trd0(:,:) * surf(:,:) )
+      IF( lk_mpp ) THEN
+         CALL mpp_sum( z_frc_trd_v )
+         CALL mpp_sum( z_frc_trd_t )
+      ENDIF
+      frc_v = frc_v + z_frc_trd_v * rdt
+      frc_t = frc_t + z_frc_trd_t * rdt
+      frc_s = frc_s + z_frc_trd_s * rdt
+
+      ! ----------------------- !
+      ! 2 -  Content variations !
+      ! ----------------------- !
+      zdiff_v2 = 0.d0
+      zdiff_hc = 0.d0
+      zdiff_sc = 0.d0
+      ! volume variation (calculated with ssh)
+      zdiff_v1 = SUM( surf(:,:) * tmask(:,:,1) * ( sshn(:,:) - ssh_ini(:,:) ) )
+      DO jk = 1, jpkm1
+         ! volume variation (calculated with scale factors)
+         zdiff_v2 = zdiff_v2 + SUM( surf(:,:) * tmask(:,:,jk)   &
+            &                       * ( fse3t_n(:,:,jk)         &
+            &                           - e3t_ini(:,:,jk) ) )
+         ! heat content variation
+         zdiff_hc = zdiff_hc + SUM( surf(:,:) * tmask(:,:,jk)          &
+            &                       * ( fse3t_n(:,:,jk) * tn(:,:,jk)   &
+            &                           - hc_loc_ini(:,:,jk) ) )
+         ! salt content variation
+         zdiff_sc = zdiff_sc + SUM( surf(:,:) * tmask(:,:,jk)          &
+            &                       * ( fse3t_n(:,:,jk) * sn(:,:,jk)   &
+            &                           - sc_loc_ini(:,:,jk) ) )
+      ENDDO
+
+      IF( lk_mpp ) THEN
+         CALL mpp_sum( zdiff_hc )
+         CALL mpp_sum( zdiff_sc )
+         CALL mpp_sum( zdiff_v1 )
+         CALL mpp_sum( zdiff_v2 )
+      ENDIF
+
+      ! Substract forcing from heat content, salt content and volume variations
+      zdiff_v1 = zdiff_v1 - frc_v
+      zdiff_v2 = zdiff_v2 - frc_v
+      zdiff_hc = zdiff_hc - frc_t
+      zdiff_sc = zdiff_sc - frc_s
+      
+      ! ----------------------- !
+      ! 3 - Diagnostics writing !
+      ! ----------------------- !
+      zdeltat  = 1.e0 / ( ( kt - nit000 + 1 ) * rdt )
+      WRITE(num , 9020) kt , zdiff_hc / vol_tot , zdiff_hc * fact1  * zdeltat,                                &
+         &                   zdiff_sc / vol_tot , zdiff_sc * fact21 * zdeltat, zdiff_sc * fact22 * zdeltat,   &
+         &                   zdiff_v1           , zdiff_v1 * fact31 * zdeltat, zdiff_v1 * fact32 * zdeltat,   &
+         &                   zdiff_v2           , zdiff_v2 * fact31 * zdeltat, zdiff_v2 * fact32 * zdeltat
+
+      IF ( kt == nitend ) CLOSE( num )
+
+9020  FORMAT(I5,11D15.7)
+
+   END SUBROUTINE dia_hsb
+
+
+   SUBROUTINE dia_hsb_init
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_hsb  ***
+      !!     
+      !! ** Purpose: Initialization for the heat salt volume budgets
+      !!	
+      !! ** Method : Compute initial heat content, salt content and volume
+      !!
+      !! ** Action : - Compute initial heat content, salt content and volume
+      !!             - Initialize forcing trends
+      !!             - Compute coefficients for conversion
+      !!---------------------------------------------------------------------------
+      CHARACTER (len=32) ::   s_name		  ! output file name
+      INTEGER		 ::   jk		  ! dummy loop indice
+      INTEGER		 ::   ierror              ! local integer
+      !!---------------------------------------------------------------------------
+
+      ! ------------------- !
+      ! 1 - Allocate memory !
+      ! ------------------- !
+      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
+      ENDIF
+      ALLOCATE( sc_loc_ini(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to allocate sc_loc_ini' )   ;   RETURN
+      ENDIF
+      ALLOCATE( e3t_ini(jpi,jpj,jpk)   , STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to allocate e3t_ini' )      ;   RETURN
+      ENDIF
+      ALLOCATE( surf(jpi,jpj)          , STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to allocate surf' )         ;   RETURN
+      ENDIF
+      ALLOCATE( ssh_ini(jpi,jpj)       , STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to allocate ssh_ini' )      ;   RETURN
+      ENDIF
+
+      ! ----------------------------------------------- !
+      ! 2 - Time independant variables and file opening !
+      ! ----------------------------------------------- !
+      WRITE(numout,*) "dia_hsb: heat salt volume budgets activated"
+      WRITE(numout,*) "~~~~~~~  output written in the 'heat_salt_volume_budgets.txt' ASCII file"
+      IF( lk_obc .OR. lk_bdy) THEN
+         CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' )         
+      ENDIF
+      s_name    = 'heat_salt_volume_budgets.txt'                         ! name of output file
+      surf(:,:) = e1t(:,:) * e2t(:,:) * tmask(:,:,1) * tmask_i(:,:)      ! masked surface grid cell area
+      surf_tot  = SUM( surf(:,:) )                                       ! total ocean surface area
+      vol_tot   = 0.d0                                                   ! total ocean volume
+      DO jk = 1, jpkm1
+         vol_tot  = vol_tot + SUM( surf(:,:) * tmask(:,:,jk)     &
+            &                      * fse3t_n(:,:,jk)         )
+      ENDDO
+      IF( lk_mpp ) THEN 
+         CALL mpp_sum( vol_tot )
+         CALL mpp_sum( surf_tot )
+      ENDIF
+
+      CALL ctl_opn( num , s_name , 'UNKNOWN' , 'FORMATTED' , 'SEQUENTIAL' , 1 , numout , lwp , 1 )
+      !                   12345678901234567890123456789012345678901234567890123456789012345678901234567890 -> 80
+      WRITE( num, 9010 ) "kt   |     heat content budget     |            salt content budget             ",   &
+         !                                                   123456789012345678901234567890123456789012345 -> 45
+         &                                                  "|            volume budget (ssh)             ",   &
+         !                                                   678901234567890123456789012345678901234567890 -> 45
+         &                                                  "|            volume budget (e3t)             "
+      WRITE( num, 9010 ) "     |      [C]         [W/m2]     |     [psu]        [mmm/s]          [SV]     ",   &
+         &                                                  "|     [m3]         [mmm/s]          [SV]     ",   &
+         &                                                  "|     [m3]         [mmm/s]          [SV]     "
+
+      ! --------------- !
+      ! 3 - Conversions ! (factors will be multiplied by duration afterwards)
+      ! --------------- !
+
+      ! heat content variation   =>   equivalent heat flux:
+      fact1  = rau0 * rcp / surf_tot                                         ! [C*m3]   ->  [W/m2]
+      ! salt content variation   =>   equivalent EMP and equivalent "flow": 
+      fact21 = 1.e3  / ( soce * surf_tot )                                   ! [psu*m3] ->  [mm/s]
+      fact22 = 1.e-6 / soce                                                  ! [psu*m3] ->  [Sv]
+      ! volume variation         =>   equivalent EMP and equivalent "flow":
+      fact31 = 1.e3  / surf_tot                                              ! [m3]     ->  [mm/s]
+      fact32 = 1.e-6                                                         ! [m3]     ->  [SV]
+
+      ! ---------------------------------- !
+      ! 4 - initial conservation variables !
+      ! ---------------------------------- !
+      ssh_ini(:,:) = sshn(:,:)                               ! initial ssh
+      DO jk = 1, jpk
+         e3t_ini   (:,:,jk) = fse3t_n(:,:,jk)                ! initial vertical scale factors
+         hc_loc_ini(:,:,jk) = tn(:,:,jk) * fse3t_n(:,:,jk)   ! initial heat content
+         sc_loc_ini(:,:,jk) = sn(:,:,jk) * fse3t_n(:,:,jk)   ! initial salt content
+      END DO
+      frc_v = 0.d0                                           ! volume       trend due to forcing
+      frc_t = 0.d0                                           ! heat content   -    -   -    -   
+      frc_s = 0.d0                                           ! salt content   -    -   -    -         
+      !
+9010  FORMAT(A80,A45,A45)
+      !
+   END SUBROUTINE dia_hsb_init
+
+   !!======================================================================
+END MODULE diahsb
Index: /branches/dev_r2091_ZTC/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- /branches/dev_r2091_ZTC/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 2115)
+++ /branches/dev_r2091_ZTC/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 2116)
@@ -147,5 +147,4 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3t_b              !: before         -      -      -    -   T      points (m)
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3u_b  , e3v_b     !:   -            -      -      -    -   U--V   points (m)
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3t_d              !: "diffused"     -      -      -    -   T      points (m)
 #else
    LOGICAL, PUBLIC, PARAMETER ::   lk_vvl = .FALSE.   !: fixed grid flag
Index: /branches/dev_r2091_ZTC/NEMO/OPA_SRC/DOM/domzgr_substitute.h90
===================================================================
--- /branches/dev_r2091_ZTC/NEMO/OPA_SRC/DOM/domzgr_substitute.h90	(revision 2115)
+++ /branches/dev_r2091_ZTC/NEMO/OPA_SRC/DOM/domzgr_substitute.h90	(revision 2116)
@@ -7,4 +7,5 @@
    !! History :  1.0  !  2005-10  (A. Beckmann, G. Madec) generalisation to all coord.
    !!            3.1  !  2009-02  (G. Madec, M. Leclair)  pure z* coordinate
+   !!            3.3  !  2010-09  (G. Madec, M. Leclair)  z_tilde coordinate
    !!----------------------------------------------------------------------
 #if defined key_zco
@@ -33,5 +34,44 @@
 #   define  fse3vw_0(i,j,k)  e3vw(i,j,k)
 #endif
-#if defined key_vvl
+#if defined key_ztc
+! z_tilde coordinate 3D arrays at t,u,v points at b,n,a time steps
+#   define  fsdept(i,j,k)  gdept_1(i,j,k)
+#   define  fsdepw(i,j,k)  gdepw_1(i,j,k)
+#   define  fsde3w(i,j,k)  gdep3w_1(i,j,k)
+#   define  fse3t(i,j,k)   e3t_1(i,j,k)
+#   define  fse3u(i,j,k)   e3u_1(i,j,k)
+#   define  fse3v(i,j,k)   e3v_1(i,j,k)
+#   define  fse3f(i,j,k)   e3f_1(i,j,k)
+#   define  fse3w(i,j,k)   e3w_1(i,j,k)
+#   define  fse3uw(i,j,k)  e3uw_1(i,j,k)
+#   define  fse3vw(i,j,k)  e3vw_1(i,j,k)
+
+#   define  fse3t_b(i,j,k)   e3t_b(i,j,k)
+#   define  fse3u_b(i,j,k)   e3u_b(i,j,k)
+#   define  fse3v_b(i,j,k)   e3v_b(i,j,k)
+#   define  fse3uw_b(i,j,k)  e3uw_b(i,j,k)
+#   define  fse3vw_b(i,j,k)  e3vw_b(i,j,k)
+
+#   define  fsdept_n(i,j,k)  gdept_1(i,j,k)
+#   define  fsdepw_n(i,j,k)  gdepw_1(i,j,k)
+#   define  fsde3w_n(i,j,k)  gdep3w_1(i,j,k)
+#   define  fse3t_n(i,j,k)   e3t_n(i,j,k)
+#   define  fse3u_n(i,j,k)   e3u_n(i,j,k)
+#   define  fse3v_n(i,j,k)   e3v_n(i,j,k)
+#   define  fse3f_n(i,j,k)   e3f_1(i,j,k)
+#   define  fse3w_n(i,j,k)   e3w_1(i,j,k)
+#   define  fse3uw_n(i,j,k)  e3uw(i,j,k)
+#   define  fse3vw_n(i,j,k)  e3vw(i,j,k)
+
+! - ML - Caution: only the first level of e3t_m is used (in limthd_2.F90)
+!                 It is replaced b y a 2d array
+#   define  fse3t_m(i,j)     e3t_m(i,j)
+#   define  fse3t_d(i,j,k)   e3t_d(i,j,k)
+
+#   define  fse3t_a(i,j,k)   e3t_a(i,j,k)
+#   define  fse3u_a(i,j,k)   e3u_a(i,j,k)
+#   define  fse3v_a(i,j,k)   e3v_a(i,j,k)
+
+#elif defined key_vvl
 ! s* or z*-coordinate (3D + time dependency) + use of additional now arrays (..._1)
 #   define  fsdept(i,j,k)  gdept_1(i,j,k)
@@ -63,5 +103,7 @@
 #   define  fse3vw_n(i,j,k)  (fse3vw_0(i,j,k)*(1.+sshv_n(i,j)*muv(i,j,k)))
 
-#   define  fse3t_m(i,j,k)   (fse3t_0(i,j,k)*(1.+ssh_m(i,j)*mut(i,j,k)))
+! - ML - Caution: only the first level of e3t_m is used (in limthd_2.F90)
+!                 It is replaced b y a 2d array
+#   define  fse3t_m(i,j)     (fse3t_0(i,j,1)*(1.+ssh_m(i,j)*mut(i,j,1)))
 #   define  fse3t_d(i,j,k)   e3t_d(i,j,k)
 
Index: /branches/dev_r2091_ZTC/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- /branches/dev_r2091_ZTC/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 2115)
+++ /branches/dev_r2091_ZTC/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 2116)
@@ -242,5 +242,7 @@
             ! -------------------------------
             DO jk = 1, jpkm1
-               fse3t_b(:,:,jk) = fse3t_n(:,:,jk) + atfp * fse3t_d(:,:,jk)
+               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
+                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
+                  &                         - 2.e0 * fse3t_n(:,:,jk)            )
             ENDDO
             ! Add volume filter correction only at the first level of t-point scale factors
Index: /branches/dev_r2091_ZTC/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- /branches/dev_r2091_ZTC/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 2115)
+++ /branches/dev_r2091_ZTC/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 2116)
@@ -1,2 +1,3 @@
+
 MODULE sshwzv   
    !!==============================================================================
@@ -156,11 +157,6 @@
       ! 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
-      IF( neuler == 0 .AND. kt == nit000 ) THEN
-         z1_rau0 = 1.e0 / rau0
-         ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 *   emp(:,:)                + zhdiv(:,:) )  ) * tmask(:,:,1)
-      ELSE
-         z1_rau0 = 0.5 / rau0
-         ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
-      ENDIF
+      z1_rau0 = 0.5 / rau0
+      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
 
 #if defined key_obc
@@ -186,13 +182,4 @@
          CALL lbc_lnk( sshv_a, 'V', 1. )
       ENDIF
-      !                                           !----------------------------------------!
-      !                                           !     vertical scale factor laplacian    !
-      !                                           !----------------------------------------!
-      ! Needed for Robert-Asselin time filter and for Brown & Campana semi implicit hydrostatic presure gradient
-      DO jk = 1, jpk
-         fse3t_d(:,:,jk) =          fse3t_b(:,:,jk)   &
-            &              - 2.e0 * fse3t_n(:,:,jk)   &
-            &              +        fse3t_a(:,:,jk)
-      ENDDO
       !                                           !------------------------------!
       !                                           !     Now Vertical Velocity    !
Index: /branches/dev_r2091_ZTC/NEMO/OPA_SRC/IOM/in_out_manager.F90
===================================================================
--- /branches/dev_r2091_ZTC/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 2115)
+++ /branches/dev_r2091_ZTC/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 2116)
@@ -78,4 +78,5 @@
    INTEGER            ::   nn_bench     =    0    !: benchmark parameter (0/1)
    INTEGER            ::   nn_bit_cmp   =    0    !: bit reproducibility  (0/1)
+   LOGICAL            ::   ln_hsb     = .FALSE.   !: check the heat and salt budgets
 
    !                                              !: OLD namelist names
Index: /branches/dev_r2091_ZTC/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- /branches/dev_r2091_ZTC/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 2115)
+++ /branches/dev_r2091_ZTC/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 2116)
@@ -45,5 +45,5 @@
    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) ::   qsr               !: 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]
@@ -53,9 +53,9 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp_tot           !: total E-P-R over ocean and ice               [Kg/m2/s]
    ! - ML - begin
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_hc_n          !: sbc heat content trend now                   [K/m/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_hc_n          !: sbc heat content trend now                   [K.m/s]
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_hc_b          !:  "   "      "      "   before                   "
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_sc_n          !: sbc salt content trend now                   [psu/m/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_sc_n          !: sbc salt content trend now                   [psu.m/s]
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_sc_b          !:  "   "      "      "   before                   "
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   qsr_hc_n      !: heat content trend due to qsr flux now       [K/m/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   qsr_hc_n      !: heat content trend due to qsr flux now       [K.m/s]
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   qsr_hc_b      !:  "      "      "    "  "   "   "   before       "
    ! - ML - end
@@ -78,4 +78,7 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sss_m     !: mean (nn_fsbc time-step) surface sea salinity            [psu]
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ssh_m     !: mean (nn_fsbc time-step) sea surface height                [m]
+#if defined key_ztc
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   e3t_m     !: mean (nn_fsbc time-step) first ocean level thickness       [m]
+#endif
 
    !!----------------------------------------------------------------------
Index: /branches/dev_r2091_ZTC/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/dev_r2091_ZTC/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 2115)
+++ /branches/dev_r2091_ZTC/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 2116)
@@ -195,10 +195,9 @@
          vtau_b(:,:) = vtau(:,:)                         ! (except at nit000 where before fields
          qns_b (:,:) = qns (:,:)                         !  are set at the end of the routine)
-         qsr_b (:,:) = qsr (:,:)
+         ! The 3D heat content due to qsr forcing is treated in traqsr
+         ! qsr_b (:,:) = qsr (:,:)
          emp_b (:,:) = emp (:,:)
          emps_b(:,:) = emps(:,:)
-         
-      ENDIF
-
+      ENDIF
       !                                            ! ---------------------------------------- !
       !                                            !        forcing field computation         !
@@ -259,5 +258,6 @@
             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)
+            ! The 3D heat content due to qsr forcing is treated in traqsr
+            ! 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', emps_b )   ! before C/D freshwater flux (T-point)
@@ -267,5 +267,5 @@
             vtau_b(:,:) = vtau(:,:)
             qns_b (:,:) = qns (:,:)
-            qsr_b (:,:) = qsr (:,:)
+            ! qsr_b (:,:) = qsr (:,:)
             emp_b (:,:) = emp (:,:)
             emps_b(:,:) = emps(:,:)
@@ -283,5 +283,6 @@
          CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau )
          CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns  )
-         CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  )
+         ! The 3D heat content due to qsr forcing is treated in traqsr
+         ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  )
          CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp  )
          CALL iom_rstput( kt, nitrst, numrow, 'emps_b' , emps )
Index: /branches/dev_r2091_ZTC/NEMO/OPA_SRC/SBC/sbcssm.F90
===================================================================
--- /branches/dev_r2091_ZTC/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 2115)
+++ /branches/dev_r2091_ZTC/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 2116)
@@ -64,4 +64,5 @@
          sss_m(:,:) = sn(:,:,1)
          ssh_m(:,:) = sshn(:,:)
+         IF( lk_vvl ) fse3t_m(:,:) = fse3t_n(:,:,1)
          !
       ELSE
@@ -79,4 +80,7 @@
                CALL iom_get( numror, jpdom_autoglo, 'sss_m'  , sss_m  )   !   "         "    salinity    (T-point)
                CALL iom_get( numror, jpdom_autoglo, 'ssh_m'  , ssh_m  )   !   "         "    height      (T-point)
+               IF( lk_vvl ) THEN
+                  CALL iom_get( numror, jpdom_autoglo, 'e3t_m'  , fse3t_m  )   !   first ocean level thickness (T-point)
+               ENDIF
                !
                IF( zf_sbc /= REAL( nn_fsbc, wp ) ) THEN      ! nn_fsbc has changed between 2 runs
@@ -89,4 +93,5 @@
                   sss_m(:,:) = zcoef * sss_m(:,:)
                   ssh_m(:,:) = zcoef * ssh_m(:,:)
+                  IF( lk_vvl ) fse3t_m(:,:) = zcoef * fse3t_m(:,:)
                ELSE
                   IF(lwp) WRITE(numout,*) '~~~~~~~   mean fields read in the ocean restart file'
@@ -100,4 +105,5 @@
                sss_m(:,:) = zcoef * sn(:,:,1)
                ssh_m(:,:) = zcoef * sshn(:,:)
+               IF( lk_vvl ) fse3t_m(:,:) = zcoef * fse3t_n(:,:,1)
             ENDIF
             !                                             ! ---------------------------------------- !
@@ -109,4 +115,5 @@
             sss_m(:,:) = 0.e0
             ssh_m(:,:) = 0.e0
+            IF( lk_vvl ) fse3t_m(:,:) = 0.e0
          ENDIF
          !                                                ! ---------------------------------------- !
@@ -118,8 +125,10 @@
          sss_m(:,:) = sss_m(:,:) + sn(:,:,1)
          ssh_m(:,:) = ssh_m(:,:) + sshn(:,:)
-
+         IF( lk_vvl ) fse3t_m(:,:) = fse3t_m(:,:) + fse3t_n(:,:,1)
          !                                                ! ---------------------------------------- !
          IF( MOD( kt - 1 , nn_fsbc ) == 0 ) THEN          !   Mean value at each nn_fsbc time-step   !
             !                                             ! ---------------------------------------- !
+            ! - ML - In a pure rigourous implementation, the averaging of sst and sss should be weighted
+            !        by the level thickness in vvl case. We suppose here that the effect is negligible.
             zcoef = 1. / REAL( nn_fsbc, wp )
             sst_m(:,:) = sst_m(:,:) * zcoef           ! mean SST             [Celcius]
@@ -128,4 +137,5 @@
             ssv_m(:,:) = ssv_m(:,:) * zcoef           !
             ssh_m(:,:) = ssh_m(:,:) * zcoef           ! mean SSH             [m]
+            IF( lk_vvl ) fse3t_m(:,:) = fse3t_m(:,:) * zcoef  !   first ocean level thickness [m]
             !
          ENDIF
@@ -144,4 +154,7 @@
             CALL iom_rstput( kt, nitrst, numrow, 'sss_m'  , sss_m  )
             CALL iom_rstput( kt, nitrst, numrow, 'ssh_m'  , ssh_m  )
+            IF( lk_vvl ) THEN
+               CALL iom_rstput( kt, nitrst, numrow, 'e3t_m'  , fse3t_m  )
+            ENDIF
             !
          ENDIF
Index: /branches/dev_r2091_ZTC/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- /branches/dev_r2091_ZTC/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 2115)
+++ /branches/dev_r2091_ZTC/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 2116)
@@ -282,5 +282,5 @@
                   ze3t_n = fse3t_n(ji,jj,jk)
                   ze3t_a = fse3t_a(ji,jj,jk)
-                  ze3t_d = fse3t_d(ji,jj,jk)
+                  ze3t_d = ze3t_b - 2.e0 * ze3t_n + ze3t_a
                   !                                         ! tracer content at Before, now and after
                   ztc_b  = tb(ji,jj,jk) * ze3t_b   ;   zsc_b = sb(ji,jj,jk) * ze3t_b
Index: /branches/dev_r2091_ZTC/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- /branches/dev_r2091_ZTC/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 2115)
+++ /branches/dev_r2091_ZTC/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 2116)
@@ -127,7 +127,7 @@
          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
-            qsr_b(:,:) = 0.e0                     
-         ENDIF
+         ! IF( kt == nit000 ) THEN             ! idem on before field at nit000
+         !    qsr_b(:,:) = 0.e0                     
+         ! ENDIF
       ENDIF
       !                                          Set before sbc tracer content fields
@@ -186,10 +186,12 @@
       !                                          Write in the ocean restart file
       !                                          *******************************
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'sbc : ocean surface tracer content forcing fields written in ocean restart file ',   &
-           &                    'at it= ', kt,' date= ', ndastp
-      IF(lwp) WRITE(numout,*) '~~~~'
-      CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_hc_n )
-      CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_sc_n )
+      IF( lrst_oce ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'sbc : ocean surface tracer content forcing fields written in ocean restart file ',   &
+            &                    'at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*) '~~~~'
+         CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_hc_n )
+         CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_sc_n )
+      ENDIF
       
       IF( l_trdtra ) THEN           ! save the sbc trends for diagnostic
Index: /branches/dev_r2091_ZTC/NEMO/OPA_SRC/opa.F90
===================================================================
--- /branches/dev_r2091_ZTC/NEMO/OPA_SRC/opa.F90	(revision 2115)
+++ /branches/dev_r2091_ZTC/NEMO/OPA_SRC/opa.F90	(revision 2116)
@@ -54,4 +54,5 @@
    USE trdmod          ! momentum/tracers trends       (trd_mod_init routine)
    USE diaptr          ! poleward transports           (dia_ptr_init routine)
+   USE diahsb          ! heat, salt and volume budgets (dia_hsb_init routine)
    USE step            ! OPA time-stepping                  (stp     routine)
 #if defined key_oasis3
@@ -180,5 +181,5 @@
       !!
       NAMELIST/namctl/ ln_ctl  , nn_print, nn_ictls, nn_ictle,   &
-         &             nn_isplt, nn_jsplt, nn_jctls, nn_jctle, nn_bench, nn_bit_cmp
+         &             nn_isplt, nn_jsplt, nn_jctls, nn_jctle, nn_bench, nn_bit_cmp, ln_hsb
       !!----------------------------------------------------------------------
       !
@@ -279,4 +280,5 @@
       CALL dia_ptr_init                         ! Poleward TRansports initialization
       CALL trd_mod_init                         ! Mixed-layer/Vorticity/Integral constraints trends
+      IF(ln_hsb)   CALL dia_hsb_init            ! heat content, salt content and volume budgets
       !
    END SUBROUTINE opa_init
@@ -313,4 +315,5 @@
          WRITE(numout,*) '      benchmark parameter (0/1)       nn_bench   = ', nn_bench
          WRITE(numout,*) '      bit comparison mode (0/1)       nn_bit_cmp = ', nn_bit_cmp
+         WRITE(numout,*) '      check the heat and salt budgets ln_hsb     = ', ln_hsb
       ENDIF
 
Index: /branches/dev_r2091_ZTC/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/dev_r2091_ZTC/NEMO/OPA_SRC/step.F90	(revision 2115)
+++ /branches/dev_r2091_ZTC/NEMO/OPA_SRC/step.F90	(revision 2116)
@@ -117,5 +117,5 @@
    USE prtctl          ! Print control                    (prt_ctl routine)
    ! - ML - 
-   USE diatrb          ! global tracer conservation       (dia_trb  routine)
+   USE diahsb          ! global tracer conservation       (dia_hsb  routine)
 
 #if defined key_agrif
@@ -320,5 +320,5 @@
 
                                CALL ssh_nxt( kstp )         ! sea surface height at next time step
-                               CALL dia_trb( kstp )         ! - ML - global conservation diagnostics
+      IF( ln_hsb           )   CALL dia_hsb( kstp )         ! - ML - global conservation diagnostics
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
