Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/ice_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/ice_2.F90	(revision 1858)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/ice_2.F90	(revision 1859)
@@ -15,5 +15,4 @@
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE par_ice_2          ! LIM sea-ice parameters
 
@@ -34,4 +33,5 @@
    INTEGER , PUBLIC ::   nbitdr = 250       !: maximum number of iterations for relaxation
    REAL(wp), PUBLIC ::   rdt_ice            !: ice time step
+   REAL(wp), PUBLIC ::   r1_rdt_ice         !: =1/rdt_ice
    REAL(wp), PUBLIC ::   epsd   = 1.0e-20   !: tolerance parameter for dynamic
    REAL(wp), PUBLIC ::   alpha  = 0.5       !: coefficient for semi-implicit coriolis
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/iceini_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/iceini_2.F90	(revision 1858)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/iceini_2.F90	(revision 1859)
@@ -11,28 +11,26 @@
    !!   'key_lim2' :                                  LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
+   !!   ice_init_2     : LIM-2 sea-ice model initialization
    !!----------------------------------------------------------------------
-   !!   ice_init_2       : sea-ice model initialization
-   !!   ice_run_2        : Definition some run parameter for ice model
-   !!----------------------------------------------------------------------
-   USE dom_oce
-   USE dom_ice_2
-   USE sbc_oce         ! surface boundary condition: ocean
-   USE sbc_ice         ! surface boundary condition: ice
-   USE phycst          ! Define parameters for the routines
-   USE ice_2
-   USE limmsh_2
-   USE limistate_2
-   USE limrst_2   
-   USE in_out_manager
+   USE dom_oce         ! ocean domain
+   USE dom_ice_2       ! LIM-2 : ice domain
+   USE sbc_oce         ! ocean surface boundary condition
+   USE sbc_ice         ! ice   surface boundary condition
+   USE phycst          ! physical constant
+   USE ice_2           ! LIM-2 variables
+   USE limmsh_2        ! LIM-2 mesh
+   USE limistate_2     ! LIM-2 inital state
+   USE limrst_2        ! LIM-2 restart
+   USE in_out_manager  ! I/O manager
       
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   ice_init_2               ! called by sbcice_lim_2.F90
+   PUBLIC   ice_init_2   ! called by sbcice_lim_2.F90
 
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !! NEMO/LIM 3.3,  UCL-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)
    !!----------------------------------------------------------------------
 
@@ -43,47 +41,19 @@
       !!                  ***  ROUTINE ice_init_2  ***
       !!
-      !! ** purpose :   
+      !! ** purpose :   LIM-2 sea-ice model initialisation
+      !!
+      !! ** input   :   namelist_ice file
+      !!                namelist namicerun (inside namelist_ice file)
       !!----------------------------------------------------------------------
-      !
-      ! Open the namelist file 
-      CALL ctl_opn( numnam_ice, 'namelist_ice', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )      
-      CALL ice_run_2                    !  read in namelist some run parameters
-                 
-      ! Louvain la Neuve Ice model
-      rdt_ice = nn_fsbc * rdttra(1)
-
-      CALL lim_msh_2                  ! ice mesh initialization
-     
-      ! Initial sea-ice state
-      IF( .NOT.ln_rstart ) THEN
-         CALL lim_istate_2            ! start from rest: sea-ice deduced from sst
-      ELSE
-         CALL lim_rst_read_2          ! start from a restart file
-      ENDIF
-      
-      tn_ice(:,:,1) = sist(:,:)         ! initialisation of ice temperature   
-      fr_i  (:,:) = 1.0 - frld(:,:)   ! initialisation of sea-ice fraction    
-      !
-   END SUBROUTINE ice_init_2
-
-
-   SUBROUTINE ice_run_2
-      !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE ice_run_2 ***
-      !!                 
-      !! ** Purpose :   Definition some run parameter for ice model
-      !!
-      !! ** Method  :   Read the namicerun namelist and check the parameter 
-      !!       values called at the first timestep (nit000)
-      !!
-      !! ** input   :   Namelist namicerun
-      !!-------------------------------------------------------------------
       NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, ln_limdmp, acrit, hsndif, hicdif
       !!-------------------------------------------------------------------
-      !                    
-      REWIND ( numnam_ice )                       ! Read Namelist namicerun 
+      !
+      !                                    ! Open the ice namelist file 
+      CALL ctl_opn( numnam_ice, 'namelist_ice', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )   
+      !
+      REWIND ( numnam_ice )                ! Read Namelist namicerun 
       READ   ( numnam_ice , namicerun )
-
-      IF(lwp) THEN
+      !
+      IF(lwp) THEN                         ! control print
          WRITE(numout,*)
          WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice'
@@ -95,6 +65,20 @@
          WRITE(numout,*) '   computation of temp. in ice  (=0) or not (=9999) hicdif = ', hicdif
       ENDIF
+      !         
+      rdt_ice    = nn_fsbc * rdttra(1)     ! set ice time step to surface tracer time step
+      r1_rdt_ice = 1.e0 / rdt_ice
       !
-   END SUBROUTINE ice_run_2
+      CALL lim_msh_2                       ! ice mesh initialization
+      !
+      !                                    ! Initial sea-ice state
+      IF( .NOT.ln_rstart ) THEN   ;   CALL lim_istate_2        ! start from rest: sea-ice deduced from sst
+      ELSE                        ;   CALL lim_rst_read_2      ! start from a restart file
+      ENDIF
+      !
+      tn_ice(:,:,1) = sist(:,:)            ! set initial ice temperature  
+      !
+      fr_i  (:,:) = 1.0 - frld(:,:)        ! set initial ice fraction    
+      !
+   END SUBROUTINE ice_init_2
 
 #else
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 1858)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 1859)
@@ -4,8 +4,9 @@
    !!           computation of the flux at the sea ice/ocean interface
    !!======================================================================
-   !! History : LIM  ! 2000-01 (H. Goosse) Original code
-   !!           2.0  ! 2002-07 (C. Ethe, G. Madec) re-writing F90
-   !!            -   ! 2006-07 (G. Madec) surface module
-   !!           2.1  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges
+   !! History :  1.0  !  2000-01  (H. Goosse) Original code
+   !!            2.0  !  2002-07  (C. Ethe, G. Madec) re-writing F90
+   !!             -   !  2006-07  (G. Madec) surface module
+   !!             -   !  2008-07  (C. Talandier,G.  Madec) 2D fields for soce and sice
+   !!            2.1  !  2010-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -17,15 +18,15 @@
    USE par_oce          ! ocean parameters
    USE dom_oce          ! ocean domain
-   USE sbc_ice          ! surface boundary condition
-   USE sbc_oce          ! surface boundary condition
+   USE sbc_ice          ! ice   surface boundary condition
+   USE sbc_oce          ! ocean surface boundary condition
    USE phycst           ! physical constants
-   USE ice_2            ! LIM sea-ice variables
+   USE albedo           ! albedo parameters
+   USE ice_2            ! LIM-2 sea-ice variables
 
    USE lbclnk           ! ocean lateral boundary condition
    USE in_out_manager   ! I/O manager
+   USE iom              ! 
+   USE prtctl           ! Print control
    USE diaar5, ONLY :   lk_diaar5
-   USE iom              ! 
-   USE albedo           ! albedo parameters
-   USE prtctl           ! Print control
    USE cpl_oasis3, ONLY : lk_cpl
 
@@ -33,16 +34,15 @@
    PRIVATE
 
-   PUBLIC lim_sbc_2     ! called by sbc_ice_lim_2
+   PUBLIC   lim_sbc_2   ! called by sbc_ice_lim_2
 
    REAL(wp)  ::   epsi16 = 1.e-16  ! constant values
    REAL(wp)  ::   rzero  = 0.e0    
    REAL(wp)  ::   rone   = 1.e0
-   REAL(wp), DIMENSION(jpi,jpj)  ::   soce_r
-   REAL(wp), DIMENSION(jpi,jpj)  ::   sice_r
+   REAL(wp), DIMENSION(jpi,jpj)  ::   soce_r, sice_r   ! ocean and ice 2D constant salinity fields (used if lk_vvl=F)
 
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -62,8 +62,8 @@
       !!              - Update 
       !!     
-      !! ** Outputs : - qsr     : sea heat flux:     solar 
-      !!              - qns     : sea heat flux: non solar
-      !!              - emp     : freshwater budget: volume flux 
-      !!              - emps    : freshwater budget: concentration/dillution 
+      !! ** Outputs : - qsr     : solar heat flux
+      !!              - qns     : non-solar heat flux including heat content of mass flux
+      !!              - emp     : mass flux
+      !!              - emps    : salt flux due to Freezing/Melting 
       !!              - utau    : sea surface i-stress (ocean referential)
       !!              - vtau    : sea surface j-stress (ocean referential)
@@ -75,60 +75,62 @@
       !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
       !!---------------------------------------------------------------------
-      INTEGER ::   kt    ! number of iteration
+      INTEGER, INTENT(in) ::   kt    ! number of iteration
       !!
-      INTEGER  ::   ji, jj           ! dummy loop indices
-      INTEGER  ::   ifvt, i1mfr, idfr               ! some switches
-      INTEGER  ::   iflt, ial, iadv, ifral, ifrdv
-      INTEGER  ::   ii0, ii1, ij0, ij1  ! temporary integers
-      REAL(wp) ::   zrdtir           ! 1. / rdt_ice
-      REAL(wp) ::   zqsr  , zqns     ! solar & non solar heat flux
-      REAL(wp) ::   zinda            ! switch for testing the values of ice concentration
-      REAL(wp) ::   zfons            ! salt exchanges at the ice/ocean interface
-      REAL(wp) ::   zemp             ! freshwater exchanges at the ice/ocean interface
-      REAL(wp) ::   zfrldu, zfrldv   ! lead fraction at U- & V-points
-      REAL(wp) ::   zutau , zvtau    ! lead fraction at U- & V-points
-      REAL(wp) ::   zu_io , zv_io    ! 2 components of the ice-ocean velocity
-! interface 2D --> 3D
-      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalb     ! albedo of ice under overcast sky
-      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalbp    ! albedo of ice under clear sky
+      INTEGER  ::   ji, jj                     ! dummy loop indices
+      INTEGER  ::   ifvt, idfr , iadv, i1mfr   ! local integers
+      INTEGER  ::   iflt, ifrdv, ial , ifral   !   -      -
+      INTEGER  ::   ii0, ii1, ij0, ij1         !   -      -
+      REAL(wp) ::   zqsr, zqns, zqhc, zemp     ! local scalars
+      REAL(wp) ::   zinda, zswitch, zcd        !   -      -
+      REAL(wp) ::   zfrldu, zutau, zu_io       !   -      -
+      REAL(wp) ::   zfrldv, zvtau, zv_io       !   -      -
+      REAL(wp) ::   zemp_snw, zfmm, zfsalt     !   -      -
       REAL(wp) ::   zsang, zmod, zztmp, zfm
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! component of ocean stress below sea-ice at I-point
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztiomi           ! module    of ocean stress below sea-ice at I-point
-      REAL(wp), DIMENSION(jpi,jpj) ::   zqnsoce          ! save qns before its modification by ice model
-
+      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalb, zalbp    ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztiomi, zqnsoce  !  -      -
       !!---------------------------------------------------------------------
-     
-      zrdtir = 1. / rdt_ice
-      
+           
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'lim_sbc_2 : LIM 2.0 sea-ice - surface boundary condition'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~   '
-
-         soce_r(:,:) = soce
-         sice_r(:,:) = sice
-         !
-         IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN
-            !                                        ! =======================
-            !                                        !  ORCA_R2 configuration
-            !                                        ! =======================
-            ii0 = 145   ;   ii1 = 180        ! Baltic Sea
+         !                              ! 2D fields for constant ice and ocean salinities
+         soce_r(:,:) = soce             !    in order to use different value in the Baltic sea
+         sice_r(:,:) = sice             !    which is much less salty than polar regions
+         !
+         IF( cp_cfg == "orca" ) THEN    ! ORCA configuration
+            IF( jp_cfg == 2       ) THEN     !  ORCA_R2 configuration
+            ii0 = 145   ;   ii1 = 180              ! Baltic Sea
             ij0 = 113   ;   ij1 = 130   ;   soce_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0
-                                            sice_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 2.e0
+                                            sice_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
+!!gm add here the R1 R05 and R025  cases
+!!          ELSEIF( jp_cfg == 1   ) THEN           !  ORCA_R1   configuration
+!!          ELSEIF( jp_cfg == 05  ) THEN           !  ORCA_R05  configuration
+!!          ELSEIF( jp_cfg == 025 ) THEN           !  ORCA_R025 configuration
+!!
+!!gm or better introduce the baltic change as a function of lat/lon of the baltic sea
+!!end gm
+            ENDIF
          ENDIF
          !
       ENDIF
 
-      !------------------------------------------!
-      !      heat flux at the ocean surface      !
-      !------------------------------------------!
-
-!!gm
-!!gm CAUTION   
-!!gm re-verifies the non solar expression, especially over open ocen
-!!gm
-      zqnsoce(:,:) = qns(:,:)
+      zqnsoce(:,:) = qns(:,:)      ! save non-solar flux prior to its modification by ice-ocean fluxes (diag.)
+      !
+      zswitch = 1       ! standard levitating sea-ice : salt exchange only
+      !
+!!gm ice embedment
+!      SELECT CASE( nn_ice_embd )       ! levitating/embedded sea-ice
+!      CASE( 0    )   ;   zswitch = 1       ! standard levitating sea-ice : salt exchange only
+!      CASE( 1, 2 )   ;   zswitch = 0       ! other levitating sea-ice or embedded sea-ice : salt and volume fluxes
+!      END SELECT                           !    
+!!gm end embedment
+
       DO jj = 1, jpj
          DO ji = 1, jpi
+            !                          !------------------------------------------!
+            !                          !      heat flux at the ocean surface      !
+            !                          !------------------------------------------!
             zinda   = 1.0   - MAX( rzero , SIGN( rone, - ( 1.0 - pfrld(ji,jj) )   ) )
             ifvt    = zinda * MAX( rzero , SIGN( rone,  - phicif(ji,jj)           ) )
@@ -138,42 +140,42 @@
             ial     = ifvt   * i1mfr + ( 1 - ifvt ) * idfr
             iadv    = ( 1  - i1mfr ) * zinda
-            ifral   = ( 1  - i1mfr * ( 1 - ial ) )   
-            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv 
-
-!!$            zinda   = 1.0 - AINT( pfrld(ji,jj) )                   !   = 0. if pure ocean else 1. (at previous time)
-!!$
-!!$            i1mfr   = 1.0 - AINT(  frld(ji,jj) )                   !   = 0. if pure ocean else 1. (at current  time)
-!!$
-!!$            IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda      !   = 1. if (snow and no ice at previous time) else 0. ???
-!!$            ELSE                             ;   ifvt = 0.
-!!$            ENDIF
-!!$
-!!$            IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases from previous to current
-!!$            ELSE                                     ;   idfr = 1.   
-!!$            ENDIF
-!!$
-!!$            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )    !   = 1. if ice (not only snow) at previous and pure ocean at current
-!!$
+            ifral   = ( 1  - i1mfr * ( 1 - ial ) )
+            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv
+
+!!gm  attempt to understand and comment the tricky flags used here....
+!
+!gm      zinda   = 1.0 - AINT( pfrld(ji,jj) )    ! = 0. free-ice ocean else 1. (after ice adv, but before ice thermo)
+!gm      i1mfr   = 1.0 - AINT(  frld(ji,jj) )    ! = 0. free-ice ocean else 1. (after ice thermo)
+!
+!gm      IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda   ! = 1. if (snow and no ice at previous time) else 0. ???
+!gm      ELSE                             ;   ifvt = 0.      ! correspond to a overmelting of snow in surface ablation
+!gm      ENDIF                                               ! 
+!
+!gm      IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases due to ice thermo
+!gm      ELSE                                     ;   idfr = 1.   
+!gm      ENDIF
+!
+!!$      iflt    = zinda  * (1 - i1mfr) * (1 - ifvt ) ! = 1. if ice (not only snow) at previous and pure ocean at current
+!
 !!$            ial     = ifvt   * i1mfr    +    ( 1 - ifvt ) * idfr
 !!$!                 snow no ice   ice         ice or nothing  lead fraction increases
 !!$!                 at previous   now           at previous
 !!$!                -> ice aera increases  ???         -> ice aera decreases ???
-!!$
+!
 !!$            iadv    = ( 1  - i1mfr ) * zinda  
 !!$!                     pure ocean      ice at
 !!$!                     at current      previous
 !!$!                        -> = 1. if ice disapear between previous and current
-!!$
+!
 !!$            ifral   = ( 1  - i1mfr * ( 1 - ial ) )  
 !!$!                            ice at     ???
 !!$!                            current         
 !!$!                         -> ???
-!!$ 
+!
 !!$            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv 
 !!$!                                                    ice disapear                           
-!!$
-!!$
-
-            !   computation the solar flux at ocean surface
+!
+            !
+            ! - computation the solar flux at ocean surface
 #if defined key_coupled 
             zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) )
@@ -181,60 +183,56 @@
             zqsr = pfrld(ji,jj) * qsr(ji,jj)  + ( 1.  - pfrld(ji,jj) ) * fstric(ji,jj)
 #endif            
-            !  computation the non solar heat flux at ocean surface
+            !
+            ! - computation the non solar heat flux at ocean surface
             zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr   &   ! part of the solar energy used in leads
-               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                            &
-               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * zrdtir    &
-               &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) * zrdtir
-
-            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! ???
-            
-            qsr  (ji,jj) = zqsr                                          ! solar heat flux 
-            qns  (ji,jj) = zqns - fdtcn(ji,jj)                           ! non solar heat flux
+               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                                &
+               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdt_ice    &
+               &       + ifrdv   * (       qfvbq(ji,jj) +             qdtcn(ji,jj) ) * r1_rdt_ice
+
+            ! - store residual heat flux (put in the ocean at the next time-step)
+            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)   ! ???
+            !
+            ! - heat content of mass exchanged between ocean and sea-ice
+            zqhc = ( rdq_snw(ji,jj) + rdq_ice(ji,jj) ) * r1_rdt_ice    ! heat flux due to sown & ice heat content exchanges
+            !            
+            qsr(ji,jj) = zqsr                                          ! solar heat flux 
+            qns(ji,jj) = zqns - fdtcn(ji,jj) + zqhc                    ! non solar heat flux
+  
+            !                          !------------------------------------------!
+            !                          !      mass flux at the ocean surface      !
+            !                          !------------------------------------------!
+            !
+            ! mass flux at the ocean-atmosphere interface (open ocean fraction = leads area)
+#if defined key_coupled
+            !                                                       ! coupled mode: 
+            zemp = + emp_tot(ji,jj)                              &       ! net mass flux over the grid cell (ice+ocean area)
+               &   - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )              ! minus the mass flux intercepted by sea-ice
+#else
+            !                                                       ! forced  mode: 
+            zemp = + emp(ji,jj)     *         frld(ji,jj)      &         ! mass flux over open ocean fraction 
+               &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &         ! liquid precip. over ice reaches directly the ocean
+               &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &         ! snow is intercepted by sea-ice (previous frld)
+#endif            
+            !
+            ! mass flux at the ocean/ice interface (sea ice fraction)
+            zemp_snw = rdm_snw(ji,jj) * r1_rdt_ice                  ! snow melting = pure water that enters the ocean
+            zfmm     = rdm_ice(ji,jj) * r1_rdt_ice                  ! Freezing minus Melting (F-M)
+
+            ! salt flux at the ice/ocean interface (sea ice fraction) [PSU*kg/m2/s]
+            zfsalt = - sice_r(ji,jj) * zfmm                         ! F-M salt exchange
+            zcd    =   soce_r(ji,jj) * zfmm                         ! concentration/dilution term due to F-M
+            !
+            ! salt flux only       : add concentration dilution term in salt flux  and no  F-M term in volume flux
+            ! salt and mass fluxes : non concentartion dilution term in salt flux  and add F-M term in volume flux
+            emps(ji,jj) = zfsalt +                  zswitch  * zcd   ! salt flux (+ C/D if no ice/ocean mass exchange)
+            emp (ji,jj) = zemp   + zemp_snw + ( 1.- zswitch) * zfmm  ! mass flux (- F/M mass flux if no ice/ocean mass exchange)
+            !
          END DO
       END DO
+
 
       CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )      
       CALL iom_put( 'qns_io_cea', qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )      
       CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1. - pfrld(:,:)) )
-
-      !------------------------------------------!
-      !      mass flux at the ocean surface      !
-      !------------------------------------------!
-
-!!gm
-!!gm CAUTION   
-!!gm re-verifies the emp & emps expression, especially the absence of 1-frld on zfm
-!!gm
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            
-#if defined key_coupled
-          zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   ! 
-             &   + rdm_snw(ji,jj) * zrdtir                                      !  freshwaterflux due to snow melting 
-#else
-!!$            !  computing freshwater exchanges at the ice/ocean interface
-!!$            zpme = - evap(ji,jj)    *   frld(ji,jj)           &   !  evaporation over oceanic fraction
-!!$               &   + tprecip(ji,jj)                           &   !  total precipitation
-!!$               &   - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )   &   !  remov. snow precip over ice
-!!$               &   - rdm_snw(ji,jj) / rdt_ice                     !  freshwaterflux due to snow melting 
-            !  computing freshwater exchanges at the ice/ocean interface
-            zemp = + emp(ji,jj)     *         frld(ji,jj)      &   !  e-p budget over open ocean fraction 
-               &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &   !  liquid precipitation reaches directly the ocean
-               &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  taking into account change in ice cover within the time step
-               &   + rdm_snw(ji,jj) * zrdtir                       !  freshwaterflux due to snow melting 
-               !                                                   !  ice-covered fraction:
-#endif            
-
-            !  computing salt exchanges at the ice/ocean interface
-            zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdm_ice(ji,jj) * zrdtir ) 
-            
-            !  converting the salt flux from ice to a freshwater flux from ocean
-            zfm  = zfons / ( sss_m(ji,jj) + epsi16 )
-            
-            emps(ji,jj) = zemp + zfm      ! surface ocean concentration/dilution effect (use on SSS evolution)
-            emp (ji,jj) = zemp            ! surface ocean volume flux (use on sea-surface height evolution)
-
-         END DO
-      END DO
 
       IF( lk_diaar5 ) THEN
@@ -250,5 +248,5 @@
       IF ( ln_limdyn ) THEN                        ! Update the stress over ice-over area (only in ice-dynamic case)
          !                                         ! otherwise the atmosphere-ocean stress is used everywhere
-
+         !
          ! ... ice stress over ocean with a ice-ocean rotation angle (at I-point)
 !CDIR NOVERRCHK
@@ -290,22 +288,22 @@
             END DO
          END DO
-
-         ! boundary condition on the stress (utau,vtau,taum)
-         CALL lbc_lnk( utau, 'U', -1. )
-         CALL lbc_lnk( vtau, 'V', -1. )
+         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )         ! lateral boundary conditions
          CALL lbc_lnk( taum, 'T',  1. )
-
+         !
       ENDIF
 
       !-----------------------------------------------!
-      !   Coupling variables                          !
+      !   Storing the transmitted variables           !
       !-----------------------------------------------!
-
-      IF ( lk_cpl ) THEN           
-         ! Ice surface temperature 
+!!gm where this is done ?????   ==>>> limthd_2   not logic ???
+!!gm      fr_i(:,:) = 1.0 - frld(:,:)       ! sea-ice fraction
+!!gm
+
+      IF ( lk_cpl ) THEN      ! coupled mode :
          tn_ice(:,:,1) = sist(:,:)          ! sea-ice surface temperature       
-         ! Computation of snow/ice and ocean albedo
+         !                                  ! snow/ice and ocean albedo
          CALL albedo_ice( tn_ice, reshape( hicif, (/jpi,jpj,1/) ), reshape( hsnif, (/jpi,jpj,1/) ), zalbp, zalb )
          alb_ice(:,:,1) =  0.5 * ( zalbp(:,:,1) + zalb (:,:,1) )   ! Ice albedo (mean clear and overcast skys)
+         !
          CALL iom_put( "icealb_cea", alb_ice(:,:,1) * fr_i(:,:) )  ! ice albedo
       ENDIF
@@ -318,6 +316,6 @@
          CALL prt_ctl(tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice  : ')
       ENDIF 
-   
-    END SUBROUTINE lim_sbc_2
+      !
+   END SUBROUTINE lim_sbc_2
 
 #else
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/thd_ice_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/thd_ice_2.F90	(revision 1858)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/thd_ice_2.F90	(revision 1859)
@@ -8,9 +8,8 @@
    !!           2.1  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges
    !!----------------------------------------------------------------------
-   !!   LIM 2.0, UCL-LOCEAN-IPSL (2005)
+   !! NEMO/LIM 3.3, UCL-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)
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE par_ice_2
 
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_3/ice.F90	(revision 1858)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_3/ice.F90	(revision 1859)
@@ -1,39 +1,22 @@
 MODULE ice
+   !!======================================================================
+   !!                        ***  MODULE ice  ***
+   !! LIM-3 sea ice :   variables defined in memory
+   !!=====================================================================
+   !! History :  3.0  ! 2008-03  (M. Vancoppenolle)  LIM3 original code
+   !!            3.1  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges
+   !!----------------------------------------------------------------------
 #if defined key_lim3
    !!----------------------------------------------------------------------
    !!   'key_lim3' :                                   LIM3 sea-ice model
    !!----------------------------------------------------------------------
-   !! History :
-   !!   2.0  !  03-08  (C. Ethe)  F90: Free form and module
-   !!   3.0  !  08-03  (M. Vancoppenolle) : LIM3 !
-   !!----------------------------------------------------------------------
-   !!  LIM 3.0, UCL-LOCEAN-IPSL (2005)
-   !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
-   !!----------------------------------------------------------------------
-   !! * Modules used
    USE par_ice          ! LIM sea-ice parameters
 
    IMPLICIT NONE
    PRIVATE
-   !!
-   !!======================================================================
-   !!                        ***  MODULE ice  ***
-   !!
-   !!                           **************
-   !!                           * L I M  3.0 *
-   !!                           **************
-   !!
-   !!                         ''in ice we trust''
-   !!
-   !!                   This module contains the sea ice 
-   !!                 diagnostics variables of ice defined 
-   !!                             in memory
-   !!
-   !!======================================================================
-   !!
+
+   !!----------------------------------------------------------------------
    !! LIM3 by the use of sweat, agile fingers and sometimes brain juice, 
    !!  was developed in Louvain-la-Neuve by : 
-   !!
    !!    * Martin Vancoppenolle (UCL-ASTR, Belgium)
    !!    * Sylvain Bouillon (UCL-ASTR, Belgium)
@@ -41,28 +24,21 @@
    !! 
    !! Based on extremely valuable earlier work by
-   !!
    !!    * Thierry Fichefet
    !!    * Hugues Goosse
    !!
    !! The following persons also contributed to the code in various ways
-   !!
    !!    * Gurvan Madec, Claude Talandier, Christian Ethe 
    !!      and Rachid Benshila (LOCEAN-IPSL, France)
    !!    * Xavier Fettweis (UCL-ASTR), Ralph Timmermann (AWI, Germany)
-   !!    * Bill Lipscomb (LANL), Cecilia Bitz (UWa) 
-   !!      and Elisabeth Hunke (LANL), USA.
+   !!    * Bill Lipscomb (LANL), Cecilia Bitz (UWa) and Elisabeth Hunke (LANL), USA.
    !! 
    !! (c) UCL-ASTR, 2005-2008
    !!
-   !! For more info, the interested user is kindly invited to consult the
-   !! following references
+   !! For more info, the interested user is kindly invited to consult the following references
    !!    For model description and validation :
    !!    * Vancoppenolle et al., Ocean Modelling, 2008a.
    !!    * Vancoppenolle et al., Ocean Modelling, 2008b.
-   !!
    !!    For a specific description of EVP :
    !!    * Bouillon et al., in prep for 2008.
-   !!
-   !!    Or the reference manual, that should be available by 2009
    !!
    !!======================================================================
@@ -183,6 +159,5 @@
    !!=====================================================================
 
-   LOGICAL, PUBLIC ::    &
-      con_i = .false.           ! switch for conservation test
+   LOGICAL, PUBLIC ::   con_i = .false.   ! switch for conservation test
 
    !!--------------------------------------------------------------------------
@@ -310,6 +285,8 @@
       phicif ,   &  !: Old ice thickness
       fbif   ,   &  !: Heat flux at the ice base
-      rdmsnif,   &  !: Variation of snow mass
-      rdmicif,   &  !: Variation of ice mass
+      rdm_snw,   &  !: Variation of snow mass over 1 time step           [Kg/m2]
+      rdq_snw,   &  !: heat content associated to rdm_snw                [J/m2]
+      rdm_ice,   &  !: Variation of ice  mass over 1 time step           [Kg/m2]
+      rdq_ice,   &  !: heat content associated to rdm_ice                [J/m2]
       qldif  ,   &  !: heat balance of the lead (or of the open ocean)
       qcmif  ,   &  !: Energy needed to bring the ocean surface layer until its freezing 
@@ -506,4 +483,8 @@
 #endif
 
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM 3.0, UCL-LOCEAN-IPSL (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!======================================================================
 END MODULE ice
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_3/thd_ice.F90	(revision 1858)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_3/thd_ice.F90	(revision 1859)
@@ -2,14 +2,9 @@
    !!======================================================================
    !!                       ***  MODULE thd_ice  ***
-   !! LIM sea-ice :   Ice thermodynamics in 1D
+   !! LIM-3 sea-ice :   Ice thermodynamics in 1D
    !!=====================================================================
-   !! History :
-   !!   2.0  !  02-11  (C. Ethe)  F90: Free form and module
+   !! History :  3.0  ! 2008-03  (M. Vancoppenolle)  LIM3 original code
+   !!            3.1  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges
    !!----------------------------------------------------------------------
-   !!   LIM 3.0, UCL-LOCEAN-IPSL (2008)
-   !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-   !! * Modules used
    USE par_ice
 
@@ -73,6 +68,8 @@
       at_i_b      ,     &  !:    "                  "      frld
       fbif_1d     ,     &  !:    "                  "      fbif
-      rdmicif_1d  ,     &  !:    "                  "      rdmicif
-      rdmsnif_1d  ,     &  !:    "                  "      rdmsnif
+      rdm_ice_1d  ,     &  !:    "                  "      rdm_ice
+      rdq_ice_1d  ,     &  !:    "                  "      rdm_ice
+      rdm_snw_1d  ,     &  !:    "                  "      rdm_snw
+      rdq_snw_1d  ,     &  !:    "                  "      rdm_snw
       qlbbq_1d    ,     &  !:    "                  "      qlbsbq
       dmgwi_1d    ,     &  !:    "                  "      dmgwi
@@ -149,15 +146,14 @@
       ftotal_fin               !: final total heat flux
 
-   REAL(wp), PUBLIC, DIMENSION(jpij,0:nlay_s) ::   &  !:
-      fc_s
-   REAL(wp), PUBLIC, DIMENSION(jpij,0:jkmax)  ::   &  !:
-      fc_i
-   REAL(wp), PUBLIC, DIMENSION(jpij,nlay_s) ::   &  !:
-      de_s_lay
-   REAL(wp), PUBLIC, DIMENSION(jpij,jkmax)  ::   &  !:
-      de_i_lay
-   INTEGER , PUBLIC ::                           &
-      jiindex_1d   ! 1D index of debugging point
+   REAL(wp), PUBLIC, DIMENSION(jpij,0:nlay_s) ::   fc_s
+   REAL(wp), PUBLIC, DIMENSION(jpij,0:jkmax)  ::   fc_i
+   REAL(wp), PUBLIC, DIMENSION(jpij,nlay_s)   ::   de_s_lay
+   REAL(wp), PUBLIC, DIMENSION(jpij,jkmax)    ::   de_i_lay
+   INTEGER , PUBLIC                           ::   jiindex_1d   ! 1D index of debugging point
 
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM 3.0, UCL-LOCEAN-IPSL (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!======================================================================
 END MODULE thd_ice
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcana.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcana.F90	(revision 1858)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcana.F90	(revision 1859)
@@ -6,4 +6,5 @@
    !! History :  3.0   ! 2006-06  (G. Madec)  Original code
    !!            3.2   ! 2009-07  (G. Madec)  Style only
+   !!            3.3  !  2010-07  (Y. Aksenov G. Madec) salt flux + heat associated with emp
    !!----------------------------------------------------------------------
 
@@ -39,5 +40,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)
@@ -60,5 +61,5 @@
       !!
       !! ** Action  : - set the ocean surface boundary condition, i.e.  
-      !!                   utau, vtau, taum, wndm, qns, qsr, emp, emps
+      !!                   utau, vtau, taum, wndm, qns, qsr, emp
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt       ! ocean time step
@@ -88,8 +89,7 @@
          !
          nn_tau000 = MAX( nn_tau000, 1 )   ! must be >= 1
-         qns   (:,:) = rn_qns0
-         qsr   (:,:) = rn_qsr0
-         emp   (:,:) = rn_emp0
-         emps  (:,:) = rn_emp0
+         emp(:,:) = rn_emp0
+         qns(:,:) = rn_qns0 - emp(:,:) * sst_m(:,:) * rcp      ! including heat content associated with mass flux at SST
+         qsr(:,:) = rn_qsr0
          !
       ENDIF
@@ -123,5 +123,5 @@
       !!
       !! ** Action  : - set the ocean surface boundary condition, i.e.   
-      !!                   utau, vtau, taum, wndm, qns, qsr, emp, emps
+      !!                   utau, vtau, taum, wndm, qns, qsr, emp
       !!
       !! Reference : Hazeleger, W., and S. Drijfhout, JPO, 30, 677-695, 2000.
@@ -204,5 +204,4 @@
          END DO
       END DO
-      emps(:,:) = emp(:,:)
 
       ! Compute the emp flux such as its integration on the whole domain at each time is zero
@@ -226,7 +225,8 @@
       ENDIF
 
-      !salinity terms
-      emp (:,:) = emp(:,:) - zsumemp * tmask(:,:,1)
-      emps(:,:) = emp(:,:)
+      
+      ! freshwater (mass flux) and update of qns with heat content of emp
+      emp (:,:) = emp(:,:) - zsumemp * tmask(:,:,1)        ! freshwater flux (=0 in domain average)
+      qns (:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp   ! evap and precip are at SST
 
 
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcblk_clio.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 1858)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 1859)
@@ -9,13 +9,14 @@
    !!            3.0  !  2008-03 (C. Talandier, G. Madec) surface module + LIM3
    !!            3.2  !  2009-04 (B. Lemaire) Introduce iom_put
+   !!            3.3  !  2010-05 (Y. Aksenov G. Madec) salt flux + heat associated with emp
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   sbc_blk_clio   : CLIO bulk formulation: read and update required input fields
-   !!   blk_clio_oce   : ocean CLIO bulk formulea: compute momentum, heat and freswater fluxes for the ocean
-   !!   blk_ice_clio   : ice   CLIO bulk formulea: compute momentum, heat and freswater fluxes for the sea-ice
+   !!   sbc_blk_clio     : CLIO bulk formulation: read and update required input fields
+   !!   blk_clio_oce     : ocean CLIO bulk formulea: compute momentum, heat and freswater fluxes for the ocean
+   !!   blk_ice_clio     : ice   CLIO bulk formulea: compute momentum, heat and freswater fluxes for the sea-ice
    !!   blk_clio_qsr_oce : shortwave radiation for ocean computed from the cloud cover
    !!   blk_clio_qsr_ice : shortwave radiation for ice   computed from the cloud cover
-   !!   flx_blk_declin : solar declinaison
+   !!   flx_blk_declin   : solar declinaison
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
@@ -47,6 +48,6 @@
    INTEGER , PARAMETER ::   jp_vtau = 2           ! index of wind stress (j-component)      (N/m2)    at V-point
    INTEGER , PARAMETER ::   jp_wndm = 3           ! index of 10m wind module                 (m/s)    at T-point
-   INTEGER , PARAMETER ::   jp_humi = 4           ! index of specific humidity               ( - )
-   INTEGER , PARAMETER ::   jp_ccov = 5           ! index of cloud cover                     ( - )
+   INTEGER , PARAMETER ::   jp_humi = 4           ! index of specific humidity               ( % )
+   INTEGER , PARAMETER ::   jp_ccov = 5           ! index of cloud cover                     ( % )
    INTEGER , PARAMETER ::   jp_tair = 6           ! index of 10m air temperature             (Kelvin)
    INTEGER , PARAMETER ::   jp_prec = 7           ! index of total precipitation (rain+snow) (Kg/m2/s)
@@ -81,5 +82,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)
@@ -98,14 +99,13 @@
       !!      the i-component of the stress                (N/m2)
       !!      the j-component of the stress                (N/m2)
-      !!      the 10m wind pseed module                    (m/s)
+      !!      the 10m wind speed module                    (m/s)
       !!      the 10m air temperature                      (Kelvin)
-      !!      the 10m specific humidity                    (-)
-      !!      the cloud cover                              (-)
+      !!      the 10m specific humidity                    (%)
+      !!      the cloud cover                              (%)
       !!      the total precipitation (rain+snow)          (Kg/m2/s)
       !!              (2) CALL blk_oce_clio
       !!
       !!      C A U T I O N : never mask the surface stress fields
-      !!                      the stress is assumed to be in the mesh referential
-      !!                      i.e. the (i,j) referential
+      !!                      the stress is assumed to be in the (i,j) mesh referential
       !!
       !! ** Action  :   defined at each time-step at the air-sea interface
@@ -113,6 +113,8 @@
       !!              - taum        wind stress module at T-point
       !!              - wndm        10m wind module at T-point
-      !!              - qns, qsr    non-slor and solar heat flux
-      !!              - emp, emps   evaporation minus precipitation
+      !!              - qns         non-solar heat flux including latent heat of solid 
+      !!                            precip. melting and emp heat content
+      !!              - qsr         solar heat flux
+      !!              - emp         upward mass flux (evap. - precip)
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in  ) ::   kt   ! ocean time step
@@ -175,10 +177,12 @@
       !                                         ! ====================== !
       !
-      CALL fld_read( kt, nn_fsbc, sf )                ! input fields provided at the current time-step
+      CALL fld_read( kt, nn_fsbc, sf )                                        ! input fields  at the current time-step
       !
 #if defined key_lim3      
-      tatm_ice(:,:) = sf(jp_tair)%fnow(:,:)     !RB ugly patch
+      tatm_ice(:,:) = sf(jp_tair)%fnow(:,:)                                   ! Tair needed in LIM-3 (!RB ugly patch)
 #endif
-      !
+      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce_clio( sf, sst_m )      ! surface ocean fluxes using CLIO bulk formulea
+      ENDIF                                               ! 
+      
       IF(lwp .AND. nitend-nit000 <= 100 ) THEN
          IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
@@ -210,8 +214,4 @@
          ENDIF
       ENDIF
-
-      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
-          CALL blk_oce_clio( sf, sst_m )                  ! compute the surface ocean fluxes using CLIO bulk formulea
-      ENDIF                                               ! 
       !
    END SUBROUTINE sbc_blk_clio
@@ -239,6 +239,8 @@
       !!               - taum        wind stress module at T-point
       !!               - wndm        10m wind module at T-point
-      !!               - qns, qsr    non-slor and solar heat flux
-      !!               - emp, emps   evaporation minus precipitation
+      !!               - qns         non-solar heat flux including latent heat of solid 
+      !!                             precip. melting and emp heat content
+      !!               - qsr         solar heat flux
+      !!               - emp         suface mass flux (evap.-precip.)
       !!  ** Nota    :   sf has to be a dummy argument for AGRIF on NEC
       !!----------------------------------------------------------------------
@@ -257,5 +259,5 @@
       REAL(wp) ::   zsst, ztatm, zcco1, zpatm, zcmax, zrmax     !    -         -
       REAL(wp) ::   zrhoa, zev, zes, zeso, zqatm, zevsqr        !    -         -
-      REAL(wp) ::   ztx2, zty2                                  !    -         -
+      REAL(wp) ::   ztx2, zty2, zcevap, zcprec                  !    -         -
       !!
       REAL(wp), DIMENSION(jpi,jpj) ::   zqlw        ! long-wave heat flux over ocean
@@ -270,10 +272,13 @@
       !------------------------------------!
 !CDIR COLLAPSE
-      DO jj = 1 , jpj
-         DO ji = 1, jpi
-            utau(ji,jj) = sf(jp_utau)%fnow(ji,jj)
-            vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj)
-         END DO
-      END DO
+      utau(:,:) = sf(jp_utau)%fnow(:,:)
+!CDIR COLLAPSE
+      vtau(:,:) = sf(jp_vtau)%fnow(:,:)
+
+      !------------------------------------!
+      !   store the wind speed  (wndm )    !
+      !------------------------------------!
+!CDIR COLLAPSE
+      wndm(:,:) = sf(jp_wndm)%fnow(:,:)
 
       !------------------------------------!
@@ -291,18 +296,7 @@
       CALL lbc_lnk( taum, 'T', 1. )
 
-      !------------------------------------!
-      !   store the wind speed  (wndm )    !
-      !------------------------------------!
-!CDIR COLLAPSE
-      DO jj = 1 , jpj
-         DO ji = 1, jpi
-            wndm(ji,jj) = sf(jp_wndm)%fnow(ji,jj)
-         END DO
-      END DO
-
       !------------------------------------------------!
       !   Shortwave radiation for ocean and snow/ice   !
       !------------------------------------------------!
-      
       CALL blk_clio_qsr_oce( qsr )
 
@@ -401,16 +395,17 @@
       !     III    Total FLUXES                                                       !
       ! ----------------------------------------------------------------------------- !
-
-!CDIR COLLAPSE
-!CDIR NOVERRCHK
-      DO jj = 1, jpj
-!CDIR NOVERRCHK
-         DO ji = 1, jpi
-            qns (ji,jj) = zqlw(ji,jj) - zqsb(ji,jj) - zqla(ji,jj)      ! Downward Non Solar flux
-            emp (ji,jj) = zqla(ji,jj) / cevap - sf(jp_prec)%fnow(ji,jj) / rday * tmask(ji,jj,1)
-         END DO
-      END DO
-      emps(:,:) = emp(:,:)
-      !
+      zcevap = rcp /  cevap    ! convert zqla ==> evap (Kg/m2/s) ==> m/s ==> W/m2
+      zcprec = rcp /  rday     ! convert prec ( mm/day ==> m/s)  ==> W/m2
+
+!CDIR COLLAPSE
+      emp(:,:) = zqla(:,:) / cevap                                        &   ! freshwater flux
+         &     - sf(jp_prec)%fnow(:,:) / rday * tmask(:,:,1)
+      !
+!CDIR COLLAPSE
+      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                        &   ! Downward Non Solar flux
+         &     - zqla(:,:)             * pst(:,:)              * zcevap   &   ! remove evap.   heat content at SST in Celcius
+         &     + sf(jp_prec)%fnow(:,:) * sf(jp_tair)%fnow(:,:) * zcprec       ! add    precip. heat content at Tair in Celcius
+      ! NB: if sea-ice model, the snow precip are computed and the associated heat is added to qns (see blk_ice_clio)
+
       CALL iom_put( "qlw_oce",   zqlw )   ! output downward longwave  heat over the ocean
       CALL iom_put( "qsb_oce", - zqsb )   ! output downward sensible  heat over the ocean
@@ -425,5 +420,5 @@
             &         tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask )
       ENDIF
-
+      !
    END SUBROUTINE blk_oce_clio
 
@@ -447,12 +442,13 @@
       !!
       !!  ** Action  :   call albedo_oce/albedo_ice to compute ocean/ice albedo 
-      !!          computation of snow precipitation
-      !!          computation of solar flux at the ocean and ice surfaces
-      !!          computation of the long-wave radiation for the ocean and sea/ice
-      !!          computation of turbulent heat fluxes over water and ice
-      !!          computation of evaporation over water
-      !!          computation of total heat fluxes sensitivity over ice (dQ/dT)
-      !!          computation of latent heat flux sensitivity over ice (dQla/dT)
-      !!
+      !!               - snow precipitation
+      !!               - solar flux at the ocean and ice surfaces
+      !!               - the long-wave radiation for the ocean and sea/ice
+      !!               - turbulent heat fluxes over water and ice
+      !!               - evaporation over water
+      !!               - total heat fluxes sensitivity over ice (dQ/dT)
+      !!               - latent heat flux sensitivity over ice (dQla/dT)
+      !!               - qns  :  modified the non solar heat flux over the ocean
+      !!                         to take into account solid precip latent heat flux
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   pst      ! ice surface temperature                   [Kelvin]
@@ -633,11 +629,20 @@
       !
       ! ----------------------------------------------------------------------------- !
-      !    Total FLUXES                                                       !
+      !    Total FLUXES                                                               !
       ! ----------------------------------------------------------------------------- !
       !
 !CDIR COLLAPSE
-      p_qns(:,:,:) = z_qlw (:,:,:) - z_qsb (:,:,:) - p_qla (:,:,:)      ! Downward Non Solar flux
+      p_qns(:,:,:) = z_qlw(:,:,:) - z_qsb(:,:,:) - p_qla(:,:,:)         ! Downward Non Solar flux
 !CDIR COLLAPSE
       p_tpr(:,:)   = sf(jp_prec)%fnow(:,:) / rday                       ! total precipitation [kg/m2/s]
+      !
+      ! ----------------------------------------------------------------------------- !
+      !    Correct the OCEAN non solar flux with the existence of solid precipitation !
+      ! ---------------=====--------------------------------------------------------- !
+!CDIR COLLAPSE
+      qns(:,:) = qns(:,:)                                           &   ! update the non-solar heat flux with:
+         &     - p_spr(:,:) * lfus                                                  &   ! remove melting solid precip
+         &     + p_spr(:,:) * MIN( sf(jp_tair)%fnow(:,:), rt0_snow - rt0 ) * cpic   &   ! add solid P at least below melting
+         &     - p_spr(:,:) * sf(jp_tair)%fnow(:,:)                        * rcp        ! remove solid precip. at Tair
       !
 !!gm : not necessary as all input data are lbc_lnk...
@@ -667,6 +672,5 @@
          CALL prt_ctl(tab2d_1=p_taui , clinfo1=' blk_ice_clio: p_taui : ', tab2d_2=p_tauj , clinfo2=' p_tauj : ')
       ENDIF
-
-
+      !
    END SUBROUTINE blk_ice_clio
 
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 1858)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 1859)
@@ -12,4 +12,5 @@
    !!            3.0  !  2006-06  (G. Madec) sbc rewritting   
    !!            3.2  !  2009-04  (B. Lemaire)  Introduce iom_put
+   !!            3.3  !  2010-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp
    !!----------------------------------------------------------------------
 
@@ -45,5 +46,5 @@
    INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point
    INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point
-   INTEGER , PARAMETER ::   jp_humi = 3           ! index of specific humidity               ( - )
+   INTEGER , PARAMETER ::   jp_humi = 3           ! index of specific humidity               ( % )
    INTEGER , PARAMETER ::   jp_qsr  = 4           ! index of solar heat                      (W/m2)
    INTEGER , PARAMETER ::   jp_qlw  = 5           ! index of Long wave                       (W/m2)
@@ -62,8 +63,8 @@
    REAL(wp), PARAMETER ::   Cice =    1.63e-3     ! transfer coefficient over ice
 
-   !                                !!* Namelist namsbc_core : CORE bulk parameters
-   LOGICAL  ::   ln_2m     = .FALSE.     ! logical flag for height of air temp. and hum
-   LOGICAL  ::   ln_taudif = .FALSE.     ! logical flag to use the "mean of stress module - module of mean stress" data
-   REAL(wp) ::   rn_pfac   = 1.          ! multiplication factor for precipitation
+   !                                    !!* Namelist namsbc_core : CORE bulk parameters
+   LOGICAL  ::   ln_2m     = .FALSE.     ! air temperature and humidity given at 2m (T) or 10m (F)
+   LOGICAL  ::   ln_taudif = .FALSE.     ! (T) use the "mean of stress module - module of mean stress" data or (F) not
+   REAL(wp) ::   rn_pfac   = 1.          ! multiplicative factor for precipitation
 
    !! * Substitutions
@@ -71,5 +72,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)
@@ -88,36 +89,32 @@
       !!      the 10m wind velocity (i-component) (m/s)    at T-point
       !!      the 10m wind velocity (j-component) (m/s)    at T-point
-      !!      the specific humidity               ( - )
+      !!      the 10m or 2m specific humidity     ( % )
       !!      the solar heat                      (W/m2)
       !!      the Long wave                       (W/m2)
-      !!      the 10m air temperature             (Kelvin)
+      !!      the 10m or 2m air temperature       (Kelvin)
       !!      the total precipitation (rain+snow) (Kg/m2/s)
       !!      the snow (solid prcipitation)       (kg/m2/s)
-      !!   OPTIONAL parameter (see ln_taudif namelist flag):
-      !!      the tau diff associated to HF tau   (N/m2)   at T-point 
+      !!      the tau diff associated to HF tau   (N/m2)   at T-point   (ln_taudif=T)
       !!              (2) CALL blk_oce_core
       !!
       !!      C A U T I O N : never mask the surface stress fields
-      !!                      the stress is assumed to be in the mesh referential
-      !!                      i.e. the (i,j) referential
+      !!                      the stress is assumed to be in the (i,j) mesh referential
       !!
       !! ** Action  :   defined at each time-step at the air-sea interface
       !!              - utau, vtau  i- and j-component of the wind stress
-      !!              - taum        wind stress module at T-point
-      !!              - wndm        10m wind module at T-point
-      !!              - qns, qsr    non-slor and solar heat flux
-      !!              - emp, emps   evaporation minus precipitation
+      !!              - taum, wndm  wind stress and 10m wind modules at T-point
+      !!              - qns, qsr    non-solar and solar heat flux
+      !!              - emp         upward mass flux (evapo. - precip.)
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in  ) ::   kt   ! ocean time step
       !!
+      INTEGER  ::   jf       ! dummy loop indice
+      INTEGER  ::   ifld     ! number of files to be read
       INTEGER  ::   ierror   ! return error code
-      INTEGER  ::   ifpr     ! dummy loop indice
-      INTEGER  ::   jfld     ! dummy loop arguments
       !!
       CHARACTER(len=100) ::  cn_dir   !   Root directory for location of core files
       TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read
-      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read
-      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !   "                                 "
-      TYPE(FLD_N) ::   sn_tdif                                 !   "                                 "
+      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr             ! informations about the fields to be read
+      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow, sn_tdif   !       -                       -
       NAMELIST/namsbc_core/ cn_dir , ln_2m  , ln_taudif, rn_pfac,           &
          &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           &
@@ -156,14 +153,14 @@
          ! do we use HF tau information?
          lhftau = ln_taudif
-         jfld = jpfld - COUNT( (/.NOT. lhftau/) )
+         ifld = jpfld - COUNT( (/.NOT. lhftau/) )
          !
          ! set sf structure
-         ALLOCATE( sf(jfld), STAT=ierror )
+         ALLOCATE( sf(ifld), STAT=ierror )
          IF( ierror > 0 ) THEN
             CALL ctl_stop( 'sbc_blk_core: unable to allocate sf structure' )   ;   RETURN
          ENDIF
-         DO ifpr= 1, jfld
-            ALLOCATE( sf(ifpr)%fnow(jpi,jpj) )
-            ALLOCATE( sf(ifpr)%fdta(jpi,jpj,2) )
+         DO jf = 1, ifld
+            ALLOCATE( sf(jf)%fnow(jpi,jpj) )
+            ALLOCATE( sf(jf)%fdta(jpi,jpj,2) )
          END DO
          !
@@ -173,14 +170,14 @@
       ENDIF
 
+!!gm    all the below lines should be executed only at nn_fbc frequency, no???   check fldread capability
+
       CALL fld_read( kt, nn_fsbc, sf )                   ! input fields provided at the current time-step
-
+      !
 #if defined key_lim3
-      tatm_ice(:,:) = sf(jp_tair)%fnow(:,:)
+      tatm_ice(:,:) = sf(jp_tair)%fnow(:,:)              ! air temperature over ice (LIM3 only)
 #endif
-
-      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
-          CALL blk_oce_core( sf, sst_m, ssu_m, ssv_m )   ! compute the surface ocean fluxes using CLIO bulk formulea
-      ENDIF
-      !                                                  ! using CORE bulk formulea
+      !                                                  ! surface ocean fluxes using CORE bulk formulea
+      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce_core( sf, sst_m, ssu_m, ssv_m )
+      !
    END SUBROUTINE sbc_blk_core
    
@@ -196,23 +193,23 @@
       !!      fields read in sbc_read
       !! 
-      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2)
-      !!              - vtau    : j-component of the stress at V-point  (N/m2)
-      !!              - taum    : Wind stress module at T-point         (N/m2)
-      !!              - wndm    : Wind speed module at T-point          (m/s)
-      !!              - qsr     : Solar heat flux over the ocean        (W/m2)
-      !!              - qns     : Non Solar heat flux over the ocean    (W/m2)
-      !!              - evap    : Evaporation over the ocean            (kg/m2/s)
-      !!              - emp(s)  : evaporation minus precipitation       (kg/m2/s)
+      !! ** Action  : - utau  : i-component of the stress at U-point  (N/m2)
+      !!              - vtau  : j-component of the stress at V-point  (N/m2)
+      !!              - taum  : Wind stress module at T-point         (N/m2)
+      !!              - wndm  : 10m Wind speed module at T-point      (m/s)
+      !!              - qsr   : Solar heat flux over the ocean        (W/m2)
+      !!              - qns   : Non Solar heat flux over the ocean    (W/m2)
+      !!                        including the latent heat of solid 
+      !!                        precip. melting and emp heat content
+      !!              - emp   : upward mass flux (evap. - precip.)    (kg/m2/s)
       !!
       !!  ** Nota  :   sf has to be a dummy argument for AGRIF on NEC
       !!---------------------------------------------------------------------
-      TYPE(fld), INTENT(in), DIMENSION(:)       ::   sf    ! input data
-      REAL(wp),  INTENT(in), DIMENSION(jpi,jpj) ::   pst   ! surface temperature                      [Celcius]
-      REAL(wp),  INTENT(in), DIMENSION(jpi,jpj) ::   pu    ! surface current at U-point (i-component) [m/s]
-      REAL(wp),  INTENT(in), DIMENSION(jpi,jpj) ::   pv    ! surface current at V-point (j-component) [m/s]
-
+      TYPE(fld), INTENT(in), DIMENSION(:)       ::   sf    ! input data (forcing field structure)
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pst   ! surface temperature                      [Celcius]
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pu    ! surface current at U-point (i-component) [m/s]
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pv    ! surface current at V-point (j-component) [m/s]
+      !!
       INTEGER  ::   ji, jj     ! dummy loop indices
-      REAL(wp) ::   zcoef_qsatw
-      REAL(wp) ::   zztmp                                 ! temporary variable
+      REAL(wp) ::   zcoef_qsatw, zztmp                    ! temporary scalar
       REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point
       REAL(wp), DIMENSION(jpi,jpj) ::   zqsatw            ! specific humidity at pst
@@ -230,5 +227,5 @@
       zcoef_qsatw = 0.98 * 640380. / rhoa
       
-      zst(:,:) = pst(:,:) + rt0      ! converte Celcius to Kelvin (and set minimum value far above 0 K)
+      zst(:,:) = pst(:,:) + rt0      ! converte SST from Celcius to Kelvin (and set minimum value far above 0 K)
 
       ! ----------------------------------------------------------------------------- !
@@ -262,5 +259,5 @@
       ! ocean albedo assumed to be 0.066
 !CDIR COLLAPSE
-      qsr (:,:) = ( 1. - 0.066 ) * sf(jp_qsr)%fnow(:,:) * tmask(:,:,1)                                 ! Short Wave
+      qsr (:,:) = ( 1. - 0.066 ) * sf(jp_qsr)%fnow(:,:) * tmask(:,:,1)                                     ! Short Wave
 !CDIR COLLAPSE
       zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave
@@ -353,9 +350,14 @@
      
 !CDIR COLLAPSE
-      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)      ! Downward Non Solar flux
-!CDIR COLLAPSE
-      emp (:,:) = zevap(:,:) - sf(jp_prec)%fnow(:,:) * rn_pfac * tmask(:,:,1)
-!CDIR COLLAPSE
-      emps(:,:) = emp(:,:)
+      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.)
+         &         - sf(jp_prec)%fnow(:,:) * rn_pfac  ) * tmask(:,:,1)
+!CDIR COLLAPSE
+      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                       &   ! Downward Non Solar flux
+         &     - sf(jp_snow)%fnow(:,:) * lfus                            &   ! remove latent melting heat for solid precip
+         &     - zevap(:,:) * pst(ji,jj) * rcp                           &   ! remove evap heat content at SST
+         &     + ( sf(jp_prec)%fnow(:,:) - sf(jp_snow)%fnow(:,:) )       &   ! add liquid precip heat content at Tair
+         &     * ( sf(jp_tair)%fnow(:,:) - rt0 ) * rcp                   &   
+         &     + sf(jp_snow)%fnow(:,:)                                   &   ! add solid  precip heat content at min(Tair,Tsnow)
+         &     * ( MIN( sf(jp_tair)%fnow(:,:), rt0_snow ) - rt0 ) * cpic 
       !
       CALL iom_put( "qlw_oce",   zqlw )                 ! output downward longwave heat over the ocean
@@ -392,34 +394,34 @@
       !! caution : the net upward water flux has with mm/day unit
       !!---------------------------------------------------------------------
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)      ::   pst      ! ice surface temperature (>0, =rt0 over land) [Kelvin]
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)    ::   pui      ! ice surface velocity (i- and i- components      [m/s]
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)    ::   pvi      !    at I-point (B-grid) or U & V-point (C-grid)
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)      ::   palb     ! ice albedo (clear sky) (alb_ice_cs)               [%]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)    ::   p_taui   ! i- & j-components of surface ice stress        [N/m2]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)    ::   p_tauj   !   at I-point (B-grid) or U & V-point (C-grid)
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:)      ::   p_qns    ! non solar heat flux over ice (T-point)         [W/m2]
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:)      ::   p_qsr    !     solar heat flux over ice (T-point)         [W/m2]
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:)      ::   p_qla    ! latent    heat flux over ice (T-point)         [W/m2]
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:)      ::   p_dqns   ! non solar heat sensistivity  (T-point)         [W/m2]
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:)      ::   p_dqla   ! latent    heat sensistivity  (T-point)         [W/m2]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)    ::   p_tpr    ! total precipitation          (T-point)      [Kg/m2/s]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)    ::   p_spr    ! solid precipitation          (T-point)      [Kg/m2/s]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)    ::   p_fr1    ! 1sr fraction of qsr penetration in ice (T-point)  [%]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)    ::   p_fr2    ! 2nd fraction of qsr penetration in ice (T-point)  [%]
-      CHARACTER(len=1), INTENT(in   )                ::   cd_grid  ! ice grid ( C or B-grid)
-      INTEGER, INTENT(in   )                         ::   pdim     ! number of ice categories
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   pst      ! ice surface temperature (>0, =rt0 over land) [Kelvin]
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   pui      ! ice surface velocity (i- and i- components      [m/s]
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   pvi      !    at I-point (B-grid) or U & V-point (C-grid)
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   palb     ! ice albedo (clear sky) (alb_ice_cs)               [%]
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_taui   ! i- & j-components of surface ice stress        [N/m2]
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_tauj   !   at I-point (B-grid) or U & V-point (C-grid)
+      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_qns    ! non solar heat flux over ice (T-point)         [W/m2]
+      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_qsr    !     solar heat flux over ice (T-point)         [W/m2]
+      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_qla    ! latent    heat flux over ice (T-point)         [W/m2]
+      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_dqns   ! non solar heat sensistivity  (T-point)         [W/m2]
+      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_dqla   ! latent    heat sensistivity  (T-point)         [W/m2]
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_tpr    ! total precipitation          (T-point)      [Kg/m2/s]
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_spr    ! solid precipitation          (T-point)      [Kg/m2/s]
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_fr1    ! 1sr fraction of qsr penetration in ice (T-point)  [%]
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_fr2    ! 2nd fraction of qsr penetration in ice (T-point)  [%]
+      CHARACTER(len=1), INTENT(in   )             ::   cd_grid  ! ice grid ( C or B-grid)
+      INTEGER, INTENT(in   )                      ::   pdim     ! number of ice categories
       !!
       INTEGER  ::   ji, jj, jl    ! dummy loop indices
       INTEGER  ::   ijpl          ! number of ice categories (size of 3rd dim of input arrays)
-      REAL(wp) ::   zst2, zst3
-      REAL(wp) ::   zcoef_wnorm, zcoef_wnorm2, zcoef_dqlw, zcoef_dqla, zcoef_dqsb
-      REAL(wp) ::   zcoef_frca                                   ! fractional cloud amount
-      REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f                  ! relative wind module and components at F-point
-      REAL(wp) ::             zwndi_t , zwndj_t                  ! relative wind components at T-point
-      REAL(wp), DIMENSION(jpi,jpj) ::   z_wnds_t                 ! wind speed ( = | U10m - U_ice | ) at T-point
-      REAL(wp), DIMENSION(jpi,jpj,pdim) ::   z_qlw               ! long wave heat flux over ice
-      REAL(wp), DIMENSION(jpi,jpj,pdim) ::   z_qsb               ! sensible  heat flux over ice
-      REAL(wp), DIMENSION(jpi,jpj,pdim) ::   z_dqlw              ! long wave heat sensitivity over ice
-      REAL(wp), DIMENSION(jpi,jpj,pdim) ::   z_dqsb              ! sensible  heat sensitivity over ice
+      REAL(wp) ::   zst2, zcoef_wnorm , zcoef_dqlw              !
+      REAL(wp) ::   zst3, zcoef_wnorm2, zcoef_dqla, zcoef_dqsb  !
+      REAL(wp) ::   zcoef_frca                                  ! fractional cloud amount
+      REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f                 ! relative wind module and components at F-point
+      REAL(wp) ::             zwndi_t , zwndj_t                 ! relative wind components at T-point
+      REAL(wp), DIMENSION(jpi,jpj)      ::   z_wnds_t           ! wind speed ( = | U10m - U_ice | ) at T-point
+      REAL(wp), DIMENSION(jpi,jpj,pdim) ::   z_qlw              ! long wave heat flux over ice
+      REAL(wp), DIMENSION(jpi,jpj,pdim) ::   z_qsb              ! sensible  heat flux over ice
+      REAL(wp), DIMENSION(jpi,jpj,pdim) ::   z_dqlw             ! long wave heat sensitivity over ice
+      REAL(wp), DIMENSION(jpi,jpj,pdim) ::   z_dqsb             ! sensible  heat sensitivity over ice
       !!---------------------------------------------------------------------
 
@@ -576,10 +578,10 @@
          CALL prt_ctl(tab2d_1=z_wnds_t, clinfo1=' blk_ice_core: z_wnds_t : ')
       ENDIF
-
+      !
    END SUBROUTINE blk_ice_core
   
 
    SUBROUTINE TURB_CORE_1Z(zu, sst, T_a, q_sat, q_a,   &
-      &                        dU, Cd, Ch, Ce   )
+      &                        dU , Cd , Ch   , Ce   )
       !!----------------------------------------------------------------------
       !!                      ***  ROUTINE  turb_core  ***
@@ -704,5 +706,6 @@
 
 
-    SUBROUTINE TURB_CORE_2Z(zt, zu, sst, T_zt, q_sat, q_zt, dU, Cd, Ch, Ce, T_zu, q_zu)
+    SUBROUTINE TURB_CORE_2Z( zt  , zu, sst, T_zt, q_sat,   &
+      &                      q_zt, dU, Cd , Ch  , Ce   , T_zu, q_zu)
       !!----------------------------------------------------------------------
       !!                      ***  ROUTINE  turb_core  ***
@@ -838,7 +841,6 @@
          Ce  = Ce_n10*sqrt_Cd/sqrt_Cd_n10/xct
          !!
-         !!
       END DO
-      !!
+      !
     END SUBROUTINE TURB_CORE_2Z
 
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 1858)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 1859)
@@ -4,7 +4,8 @@
    !! Surface Boundary Condition :  momentum, heat and freshwater fluxes in coupled mode
    !!======================================================================
-   !! History :  2.0  !  06-2007  (R. Redler, N. Keenlyside, W. Park) Original code split into flxmod & taumod
-   !!            3.0  !  02-2008  (G. Madec, C Talandier)  surface module
-   !!            3.1  !  02-2009  (G. Madec, S. Masson, E. Maisonave, A. Caubel) generic coupled interface
+   !! History :  2.0  !  2007-06  (R. Redler, N. Keenlyside, W. Park) Original code split into flxmod & taumod
+   !!            3.0  !  2008-02  (G. Madec, C Talandier)  surface module
+   !!            3.1  !  2009-02  (G. Madec, S. Masson, E. Maisonave, A. Caubel) generic coupled interface
+   !!            3.3  !  2010-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp
    !!----------------------------------------------------------------------
 #if defined key_oasis3 || defined key_oasis4
@@ -156,5 +157,5 @@
 #  include "vectopt_loop_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)
@@ -562,6 +563,8 @@
       !! ** Action  :   update  utau, vtau   ocean stress at U,V grid 
       !!                        taum, wndm   wind stres and wind speed module at T-point
-      !!                        qns , qsr    non solar and solar ocean heat fluxes   ('ocean only case)
-      !!                        emp = emps   evap. - precip. (- runoffs) (- calving) ('ocean only case)
+      !!                        qns          non solar heat fluxes including emp heat content    (ocean only case)
+      !!                                     and the latent heat flux of solid precip. melting
+      !!                        qsr          solar ocean heat fluxes   (ocean only case)
+      !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt       ! ocean model time step index
@@ -697,5 +700,5 @@
       ENDIF
 
-      ! u(v)tau and taum will be modified by ice model (wndm will be changed by PISCES)
+      ! u(v)tau and taum will be modified by ice model (and wndm will be changed by PISCES)
       ! -> need to be reset before each call of the ice/fsbc      
       IF( MOD( kt-1, k_fsbc ) == 0 ) THEN
@@ -712,15 +715,5 @@
          !                                                   ! ========================= !
          !
-         !                                                       ! non solar heat flux over the ocean (qns)
-         IF( srcv(jpr_qnsoce)%laction )   qns(:,:) = frcv(:,:,jpr_qnsoce)
-         IF( srcv(jpr_qnsmix)%laction )   qns(:,:) = frcv(:,:,jpr_qnsmix)        
-         !   energy for melting solid precipitation over free ocean
-         zcoef = xlsn / rhosn
-         qns(:,:) = qns(:,:) - frcv(:,:,jpr_snow) * zcoef
-         !                                                       ! solar flux over the ocean          (qsr)
-         IF( srcv(jpr_qsroce)%laction )   qsr(:,:) = frcv(:,:,jpr_qsroce) 
-         IF( srcv(jpr_qsrmix)%laction )   qsr(:,:) = frcv(:,:,jpr_qsrmix)
-         !
-         !                                                       ! total freshwater fluxes over the ocean (emp, emps)
+         !                                                       ! total freshwater fluxes over the ocean (emp)
          SELECT CASE( TRIM( cn_rcv_emp ) )                                    ! evaporation - precipitation
          CASE( 'conservative' )
@@ -752,15 +745,22 @@
 !!         ENDIF
 !!gm  end of internal cooking
-         !
-         emps(:,:) = emp(:,:)                                        ! concentration/dilution = emp
-  
-         !                                                           ! 10 m wind speed
-         IF( srcv(jpr_w10m)%laction )   wndm(:,:) = frcv(:,:,jpr_w10m)
+         !  
+         !                                                       ! non solar heat flux over the ocean (qns)
+         IF( srcv(jpr_qnsoce)%laction )   qns(:,:) = frcv(:,:,jpr_qnsoce)
+         IF( srcv(jpr_qnsmix)%laction )   qns(:,:) = frcv(:,:,jpr_qnsmix)
+         !
+         zcoef = xlsn / rhosn                                    ! qns update over free ocean with:
+         qns(:,:) = qns(:,:) - frcv(:,:,jpr_snow) * zcoef            ! energy for melting solid precipitation over free ocean
+            &                - emp(:,:) * sst_m(:,:) * rcp           ! remove heat content due to mass flux (assumed to be at SST)
+         !
+         !                                                       ! solar flux over the ocean          (qsr)
+         IF( srcv(jpr_qsroce)%laction )   qsr(:,:) = frcv(:,:,jpr_qsroce) 
+         IF( srcv(jpr_qsrmix)%laction )   qsr(:,:) = frcv(:,:,jpr_qsrmix)
          !
 #if defined  key_cpl_carbon_cycle
-         !                                                              ! atmosph. CO2 (ppm)
+         !                                                       ! atmosph. CO2 (ppm)
          IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(:,:,jpr_co2)
 #endif
-
+         !
       ENDIF
       !
@@ -1046,5 +1046,5 @@
       !!----------------------------------------------------------------------
       zicefr(:,:,1) = 1.- p_frld(:,:,1)
-      IF( lk_diaar5 )   zcptn(:,:) = rcp * tn(:,:,1)
+      zcptn(:,:) = rcp * sst_m(:,:)
       !
       !                                                      ! ========================= !
@@ -1118,23 +1118,22 @@
             &                                                   +          pist(:,:,1)   * zicefr(:,:,1) ) )
       END SELECT
-      !                                                           ! snow melting heat flux ....
-      !   energy for melting solid precipitation over ice-free ocean
-      zcoef = xlsn / rhosn
+      !
+      zcoef = xlsn / rhosn                                        ! qns_tot update over free ocean with:
       ztmp(:,:) = p_frld(:,:,1) * zsnow(:,:) * zcoef
-      pqns_tot(:,:) = pqns_tot(:,:) - ztmp(:,:)
-      IF( lk_diaar5 )   CALL iom_put( 'hflx_snow_cea', ztmp + zsnow(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
-!!gm
-!!    currently it is taken into account in leads budget but not in the qns_tot, and thus not in 
-!!    the flux that enter the ocean....
-!!    moreover 1 - it is not diagnose anywhere.... 
-!!             2 - it is unclear for me whether this heat lost is taken into account in the atmosphere or not...
-!!
-!! similar job should be done for snow and precipitation temperature
-      !                                                           ! Iceberg melting heat flux ....
-      !   energy for iceberg melting 
-      IF( srcv(jpr_cal)%laction ) THEN 
+      pqns_tot(:,:) = pqns_tot(:,:)                       &
+         &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting
+         &          + (  pemp_tot(:,:)                    &            ! remove the heat content of mass flux (assumed to be at SST)
+         &             - pemp_ice(:,:) * p_frld(:,:,1)  ) * zcptn(:,:) 
+         !
+      IF( lk_diaar5 )   CALL iom_put( 'hflx_snow_cea', ztmp(:,:) + zsnow(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
+
+!!gm BUG ???   just above the right value should be : ztmp + zsnow*p_frld*zcptn
+
+      !                               
+      IF( srcv(jpr_cal)%laction ) THEN                                 ! remove the latent heat flux of iceberg melting
          zcoef = xlic / rhoic
          ztmp(:,:) = frcv(:,:,jpr_cal) * zcoef
          pqns_tot(:,:) = pqns_tot(:,:) - ztmp(:,:)
+         !
          IF( lk_diaar5 )   CALL iom_put( 'hflx_cal_cea', ztmp + frcv(:,:,jpr_cal) * zcptn(:,:) )   ! heat flux from calving
       ENDIF
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcflx.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 1858)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 1859)
@@ -4,5 +4,6 @@
    !! Ocean forcing:  momentum, heat and freshwater flux formulation
    !!=====================================================================
-   !! History :  9.0   !  06-06  (G. Madec)  Original code
+   !! History :  1.0  !  2006-06  (G. Madec)  Original code
+   !!            3.3  !  2010-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp
    !!----------------------------------------------------------------------
 
@@ -52,5 +53,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2006) 
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -74,6 +75,5 @@
       !!
       !!      CAUTION :  - never mask the surface stress fields
-      !!                 - the stress is assumed to be in the mesh referential
-      !!                   i.e. the (i,j) referential
+      !!                 - the stress is assumed to be in the (i,j) mesh referential
       !!
       !! ** Action  :   update at each time-step
@@ -81,6 +81,7 @@
       !!              - taum        wind stress module at T-point
       !!              - wndm        10m wind module at T-point
-      !!              - qns, qsr    non-slor and solar heat flux
-      !!              - emp, emps   evaporation minus precipitation
+      !!              - qns         non solar heat flux including heat flux due to emp
+      !!              - qsr         solar heat flux
+      !!              - emp         upward mass flux (evap. - precip.)
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
@@ -136,10 +137,7 @@
       ENDIF
 
-      CALL fld_read( kt, nn_fsbc, sf )           ! Read input fields and provides the
-      !                                          ! input fields at the current time-step
+      CALL fld_read( kt, nn_fsbc, sf )           ! input fields at the current time-step
 
-      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
-         !
-         ! set the ocean fluxes from read fields
+      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN       ! set the ocean fluxes from read fields
 !CDIR COLLAPSE
          DO jj = 1, jpj
@@ -151,8 +149,9 @@
                emp (ji,jj) = sf(jp_emp )%fnow(ji,jj)
             END DO
-         END DO
-         
-         ! module of wind stress and wind speed at T-point
-         zcoef = 1. / ( zrhoa * zcdrag ) 
+         END DO         
+         !                                       ! add to qns the heat due to e-p
+         qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp      ! mass flux are at SST
+
+         zcoef = 1. / ( zrhoa * zcdrag )         ! module of wind stress and wind speed at T-point
 !CDIR NOVERRCHK
          DO jj = 2, jpjm1
@@ -167,10 +166,6 @@
          END DO
          CALL lbc_lnk( taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( wndm(:,:), 'T', 1. )
-
-         ! Initialization of emps (when no ice model)
-         emps(:,:) = emp (:,:) 
                   
-         ! control print (if less than 100 time-step asked)
-         IF( nitend-nit000 <= 100 .AND. lwp ) THEN
+         IF( nitend-nit000 <= 100 .AND. lwp ) THEN      ! control print (if less than 100 time-step asked)
             WRITE(numout,*) 
             WRITE(numout,*) '        read daily momentum, heat and freshwater fluxes OK'
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcfwb.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 1858)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 1859)
@@ -4,8 +4,9 @@
    !! Ocean fluxes   : domain averaged freshwater budget
    !!======================================================================
-   !! History :  8.2  !  01-02  (E. Durand)  Original code
-   !!            8.5  !  02-06  (G. Madec)  F90: Free form and module
-   !!            9.0  !  06-08  (G. Madec)  Surface module
-   !!            9.2  !  09-07  (C. Talandier) emp mean s spread over erp area 
+   !! History :  OPA  !  2001-02  (E. Durand)  Original code
+   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
+   !!            3.0  !  2006-08  (G. Madec)  Surface module
+   !!            3.2  !  2009-07  (C. Talandier) emp mean s spread over erp area 
+   !!            3.3  !  2010-05  (Y. Aksenov G. Madec) embedded sea-ice case
    !!----------------------------------------------------------------------
 
@@ -27,10 +28,10 @@
    PRIVATE
 
-   PUBLIC   sbc_fwb      ! routine called by step
-
-   REAL(wp) ::   a_fwb_b            ! annual domain averaged freshwater budget
-   REAL(wp) ::   a_fwb              ! for 2 year before (_b) and before year.
-   REAL(wp) ::   empold             ! empold to be suppressed
-   REAL(wp) ::   area               ! global mean ocean surface (interior domain)
+   PUBLIC   sbc_fwb    ! routine called by step
+
+   REAL(wp) ::   a_fwb_b   ! annual domain averaged freshwater budget
+   REAL(wp) ::   a_fwb     ! for 2 year before (_b) and before year.
+   REAL(wp) ::   empold    ! empold to be suppressed
+   REAL(wp) ::   area      ! global mean ocean surface (interior domain)
 
    REAL(wp), DIMENSION(jpi,jpj) ::   e1e2_i    ! area of the interior domain (e1t*e2t*tmask_i)
@@ -40,5 +41,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!  OPA 9.0 , LOCEAN-IPSL (2006) 
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -57,5 +58,5 @@
       !!                =2 annual global mean corrected from previous year
       !!                =3 global mean of emp set to zero at each nn_fsbc time step
-      !!                   & spread out over erp area depending its sign
+      !!                   and spread out over erp area depending its sign
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
@@ -63,10 +64,10 @@
       INTEGER, INTENT( in ) ::   kn_fwb   ! ocean time-step index
       !!
-      INTEGER  ::   inum                  ! temporary logical unit
-      INTEGER  ::   ikty, iyear           ! 
-      REAL(wp) ::   z_emp, z_emp_nsrf, zsum_emp, zsum_erp       ! temporary scalars
-      REAL(wp) ::   zsurf_neg, zsurf_pos, zsurf_tospread
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztmsk_neg, ztmsk_pos, ztmsk_tospread
-      REAL(wp), DIMENSION(jpi,jpj) ::   z_wgt, zerp_cor
+      INTEGER  ::   inum          ! temporary logical unit
+      INTEGER  ::   ikty, iyear   ! 
+      REAL(wp) ::   z_emp, z_emp_nsrf, zsum_emp, zsum_erp, zcoef     ! temporary scalars
+      REAL(wp) ::   zsurf_neg, zsurf_pos, zsurf_tospread             !    -         -
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztmsk_neg, ztmsk_tospread    ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztmsk_pos, z_wgt, zerp_cor   !  -      -
       !!----------------------------------------------------------------------
       !
@@ -83,5 +84,4 @@
             IF( kn_fwb == 3 .AND. nn_sssr /= 2 )   &
                &   CALL ctl_stop( 'The option nn_fwb = 3 must be associated to nn_sssr = 2 ' )
-            
          ENDIF
          !
@@ -99,12 +99,11 @@
          CALL ctl_stop( ctmp1 )
          !
-         
-      !
       CASE ( 1 )                               ! global mean emp set to zero
          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
             z_emp = SUM( e1e2_i(:,:) * emp(:,:) ) / area
             IF( lk_mpp )   CALL  mpp_sum( z_emp    )   ! sum over the global domain
-            emp (:,:) = emp (:,:) - z_emp
-            emps(:,:) = emps(:,:) - z_emp
+            zcoef = z_emp * rcp
+            emp(:,:) = emp(:,:) - z_emp
+            qns(:,:) = qns(:,:) + zcoef * sst_m(:,:)
          ENDIF
          !
@@ -138,6 +137,7 @@
          ! correct the freshwater fluxes
          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
-            emp (:,:) = emp (:,:) + empold
-            emps(:,:) = emps(:,:) + empold
+            zcoef = z_emp * rcp
+            emp(:,:) = emp (:,:) + empold
+            qns(:,:) = qns(:,:) - zcoef * sst_m(:,:)
          ENDIF
          !
@@ -152,26 +152,24 @@
          !
          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
-            ! Select <0 and >0 area of erp
+            !                                       ! Select <0 and >0 area of erp
             ztmsk_pos(:,:) = tmask_i(:,:)
-            WHERE( erp < 0.e0 ) ztmsk_pos = 0.e0
+            WHERE( erp(:,:) < 0.e0 )   ztmsk_pos(:,:) = 0.e0
             ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:)
 
-            ! Area filled by <0 and >0 erp 
-            zsurf_neg = SUM( e1e2_i(:,:)*ztmsk_neg(:,:) )
-            zsurf_pos = SUM( e1e2_i(:,:)*ztmsk_pos(:,:) )
+            !                                       ! Area filled by <0 and >0 erp 
+            zsurf_neg = SUM(  e1e2_i(:,:) * ztmsk_neg(:,:)  )
+            zsurf_pos = SUM(  e1e2_i(:,:) * ztmsk_pos(:,:)  )
         
-            ! emp global mean 
+            !                                       ! emp global mean 
             z_emp = SUM( e1e2_i(:,:) * emp(:,:) ) / area
             !
-            IF( lk_mpp )   CALL  mpp_sum( z_emp )
+            IF( lk_mpp )   CALL  mpp_sum( z_emp     )
             IF( lk_mpp )   CALL  mpp_sum( zsurf_neg )
             IF( lk_mpp )   CALL  mpp_sum( zsurf_pos )
             
-            IF( z_emp < 0.e0 ) THEN
-                ! to spread out over >0 erp area to increase evaporation damping process
+            IF( z_emp < 0.e0 ) THEN                 ! spread out over >0 erp area to increase evaporation damping process
                 zsurf_tospread = zsurf_pos
                 ztmsk_tospread(:,:) = ztmsk_pos(:,:)
-            ELSE
-                ! to spread out over <0 erp area to increase precipitation damping process
+            ELSE                                    ! spread out over <0 erp area to increase precipitation damping process
                 zsurf_tospread = zsurf_neg
                 ztmsk_tospread(:,:) = ztmsk_neg(:,:)
@@ -192,7 +190,7 @@
             CALL lbc_lnk( zerp_cor, 'T', 1. )
 
-            emp (:,:) = emp (:,:) + zerp_cor(:,:)
-            emps(:,:) = emps(:,:) + zerp_cor(:,:)
-            erp (:,:) = erp (:,:) + zerp_cor(:,:)
+            emp(:,:) = emp(:,:) + zerp_cor(:,:)
+            qns(:,:) = qns(:,:) - zerp_cor(:,:) * rcp * sst_m(:,:)
+            erp(:,:) = erp(:,:) + zerp_cor(:,:)
             
             IF( nprint == 1 .AND. lwp ) THEN
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcice_if.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcice_if.F90	(revision 1858)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcice_if.F90	(revision 1859)
@@ -30,5 +30,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)
@@ -72,9 +72,9 @@
          !             !   file    ! frequency !  variable  ! time intep !  clim  ! 'yearly' or ! weights  ! rotation   !
          !             !   name    !  (hours)  !   name     !   (T/F)    !  (T/F) !  'monthly'  ! filename ! pairs      ! 
-         sn_ice = FLD_N('ice_cover',    -1    ,  'ice_cov' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
-
-         REWIND ( numnam )               ! ... read in namlist namiif
+         sn_ice = FLD_N('ice_cover',    -1     ,  'ice_cov' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
+         !
+         REWIND ( numnam )               ! read in namlist namiif
          READ   ( numnam, namsbc_iif )
-
+         !
          ALLOCATE( sf_ice(1), STAT=ierror )
          IF( ierror > 0 ) THEN
@@ -83,13 +83,11 @@
          ALLOCATE( sf_ice(1)%fnow(jpi,jpj) )
          ALLOCATE( sf_ice(1)%fdta(jpi,jpj,2) )
-
-
-         ! fill sf_ice with sn_ice and control print
+         !
+         !                               ! fill sf_ice with sn_ice and control print
          CALL fld_fill( sf_ice, (/ sn_ice /), cn_dir, 'sbc_ice_if', 'ice-if sea-ice model', 'namsbc_iif' )
          !
       ENDIF
 
-      CALL fld_read( kt, nn_fsbc, sf_ice )           ! Read input fields and provides the
-      !                                              ! input fields at the current time-step
+      CALL fld_read( kt, nn_fsbc, sf_ice )           ! Read input fields at the current time-step
       
       IF( MOD( kt-1, nn_fsbc) == 0 ) THEN
@@ -101,7 +99,6 @@
          fr_i(:,:) = tfreez( sss_m ) * tmask(:,:,1)      ! sea surface freezing temperature [Celcius]
 
-         ! Flux and ice fraction computation
 !CDIR COLLAPSE
-         DO jj = 1, jpj
+         DO jj = 1, jpj          ! Flux and ice fraction computation
             DO ji = 1, jpi
                !
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 1858)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 1859)
@@ -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  !  2007-06  (G. Madec)  Original code
+   !!            3.1  !  2008-08  (S. Masson, E. .... ) coupled interface
+   !!            3.3  !  2010-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp
    !!----------------------------------------------------------------------
 
@@ -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
@@ -123,4 +124,7 @@
       ENDIF
       IF( nn_ice == 0  )   fr_i(:,:) = 0.e0        ! no ice in the domain, ice fraction is always zero
+
+      emps(:,:) = 0.e0                             ! the salt flux will be computed (i.e. will be non-zero) only if 
+      !                                            ! sea-ice is present, or lk_vvl=F, or surface salt restoring is used.
 
       !                                            ! restartability   
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 1858)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 1859)
@@ -8,4 +8,5 @@
    !!            3.0  !  2006-07  (G. Madec)  Surface module 
    !!            3.2  !  2009-04  (B. Lemaire)  Introduce iom_put
+   !!             -   !  2010-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp
    !!----------------------------------------------------------------------
 
@@ -41,8 +42,9 @@
    REAL(wp), PUBLIC, DIMENSION(jpk)     ::   rnfmsk_z    !: river mouth mask (vert.)
 
+   REAL(wp) ::   rfact_rcp   ! = rn_rfact * rcp
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnf   ! structure of input SST (file information, fields read)
 
    !!----------------------------------------------------------------------
-   !! 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)
@@ -64,5 +66,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt          ! ocean time step
-      !!
+      !
       INTEGER  ::   ji, jj   ! dummy loop indices
       INTEGER  ::   ierror   ! temporary integer
@@ -78,5 +80,7 @@
             ALLOCATE( sf_rnf(1)%fdta(jpi,jpj,2) )
          ENDIF
-         CALL sbc_rnf_init(sf_rnf)
+         CALL sbc_rnf_init( sf_rnf )
+         !
+         rfact_rcp = rn_rfact * rcp
       ENDIF
 
@@ -85,11 +89,9 @@
          !                                                !-------------------!
          !
-         CALL fld_read( kt, nn_fsbc, sf_rnf )   ! Read Runoffs data and provides it
-         !                                      ! at the current time-step
-
-         ! Runoff reduction only associated to the ORCA2_LIM configuration
-         ! when reading the NetCDF file runoff_1m_nomask.nc
-         IF( cp_cfg == 'orca' .AND. jp_cfg == 2 )   THEN
-            DO jj = 1, jpj
+         CALL fld_read( kt, nn_fsbc, sf_rnf )                    ! Read Runoffs data at the current time-step
+         !
+!!gm CAUTION this is ugly  ===>>> to be removed!
+         IF( cp_cfg == 'orca' .AND. jp_cfg == 2 )   THEN         ! Runoff reduction only associated to the ORCA2_LIM configuration
+            DO jj = 1, jpj                                       ! when reading the NetCDF file runoff_1m_nomask.nc
                DO ji = 1, jpi
                   IF( gphit(ji,jj) > 40 .AND. gphit(ji,jj) < 65 )   sf_rnf(1)%fnow(ji,jj) = 0.85 * sf_rnf(1)%fnow(ji,jj)
@@ -97,10 +99,9 @@
             END DO
          ENDIF
-
-         ! C a u t i o n : runoff is negative and in kg/m2/s 
-
-         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
-            emp (:,:) = emp (:,:) - rn_rfact * ABS( sf_rnf(1)%fnow(:,:) )
-            emps(:,:) = emps(:,:) - rn_rfact * ABS( sf_rnf(1)%fnow(:,:) )
+         !
+         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN                    ! C a u t i o n : runoff is negative and in kg/m2/s 
+            emp(:,:) = emp(:,:) - rn_rfact  * ABS( sf_rnf(1)%fnow(:,:) )                ! mass flux
+            qns(:,:) = qns(:,:) + rfact_rcp * ABS( sf_rnf(1)%fnow(:,:) ) * sst_m(:,:)   ! its associated heat content (at SST)
+            !
             CALL iom_put( "runoffs", sf_rnf(1)%fnow )         ! runoffs
          ENDIF
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcssr.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 1858)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 1859)
@@ -6,4 +6,6 @@
    !! History :  3.0  !  2006-06  (G. Madec)  Original code
    !!            3.2  !  2009-04  (B. Lemaire)  Introduce iom_put
+   !!             -   !  2009-07  (C. Talandier, G. Madec)  Add a bound to the Erp 
+   !!            3.3  !  2010-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp
    !!----------------------------------------------------------------------
 
@@ -37,5 +39,5 @@
    REAL(wp)        ::   rn_deds     = -27.70    ! restoring factor on SST and SSS
    LOGICAL         ::   ln_sssr_bnd = .false.   ! flag to bound erp term 
-   REAL(wp)        ::   rn_sssr_bnd =   0.e0    ! ABS(Max./Min.) value of erp term [mm/day]
+   REAL(wp)        ::   rn_sssr_bnd =   4.e0    ! ABS(Max./Min.) value of erp term [mm/day]
 
    REAL(wp) , ALLOCATABLE, DIMENSION(:) ::   buffer   ! Temporary buffer for exchange
@@ -46,5 +48,5 @@
 #  include "domzgr_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)
@@ -70,11 +72,8 @@
       !!
       INTEGER  ::   ji, jj   ! dummy loop indices
-      REAL(wp) ::   zerp     ! local scalar for evaporation damping
-      REAL(wp) ::   zqrp     ! local scalar for heat flux damping
-      REAL(wp) ::   zsrp     ! local scalar for unit conversion of rn_deds factor
-      REAL(wp) ::   zerp_bnd ! local scalar for unit conversion of rn_epr_max factor
       INTEGER  ::   ierror   ! return error code
-      !!
-      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
+      REAL(wp) ::   zerp, zqrp, zsrp, zerp_bnd    ! local scalar
+      !!
+      CHARACTER(len=100) ::  cn_dir = './'   ! Root directory for location of ssr files
       TYPE(FLD_N) ::   sn_sst, sn_sss        ! informations about the fields to be read
       NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, rn_dqdt, rn_deds, sn_sst, sn_sss, ln_sssr_bnd, rn_sssr_bnd
@@ -85,5 +84,4 @@
          !                                            ! -------------------- !
          !                            !* set file information
-         cn_dir  = './'            ! directory in which the model is executed
          ! ... default values (NB: frequency positive => hours, negative => months)
          !            !   file    ! frequency !  variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   !
@@ -158,8 +156,8 @@
                   END DO
                END DO
-               CALL iom_put( "qrp", qrp )                             ! heat flux damping
-            ENDIF
-            !
-            IF( nn_sssr == 1 ) THEN                   !* Salinity damping term (salt flux, emps only)
+               CALL iom_put( "qrp", qrp )                             ! heat flux damping 
+            ENDIF
+            !
+            IF( nn_sssr == 1 ) THEN                   !* Salinity damping term (salt flux only (emps))
                zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s]
 !CDIR COLLAPSE
@@ -167,13 +165,12 @@
                   DO ji = 1, jpi
                      zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
-                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj) )   &
-                        &        / ( sss_m(ji,jj) + 1.e-20   )
+                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj) )
                      emps(ji,jj) = emps(ji,jj) + zerp
-                     erp( ji,jj) = zerp
+                     erp( ji,jj) = zerp / MAX( sss_m(ji,jj), 1.e-20  )  ! converted into an equivalent emp (diag. only)
                   END DO
                END DO
                CALL iom_put( "erp", erp )                             ! freshwater flux damping
                !
-            ELSEIF( nn_sssr == 2 ) THEN               !* Salinity damping term (volume flux, emp and emps)
+            ELSEIF( nn_sssr == 2 ) THEN               !* Salinity damping term (volume flux (emp) and qns)
                zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s]
                zerp_bnd = rn_sssr_bnd / rday                          !       -              -    
@@ -183,9 +180,10 @@
                      zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
                         &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj) )   &
-                        &        / ( sss_m(ji,jj) + 1.e-20   )
+                        &        / MAX(  sss_m(ji,jj), 1.e-20  )
                      IF( ln_sssr_bnd )   zerp = SIGN( 1., zerp ) * MIN( zerp_bnd, ABS(zerp) )
-                     emp (ji,jj) = emp (ji,jj) + zerp
-                     emps(ji,jj) = emps(ji,jj) + zerp
-                     erp (ji,jj) = zerp
+!!gm better coding   IF( ln_sssr_bnd )   zerp = MAX( -zerp_bnd, MIN( zerp, zerp_bnd )  )
+                     emp(ji,jj) = emp(ji,jj) + zerp
+                     qns(ji,jj) = qns(ji,jj) - zerp * rcp * sst_m(ji,jj)
+                     erp(ji,jj) = zerp
                   END DO
                END DO
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 1858)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 1859)
@@ -1,10 +1,11 @@
 MODULE trasbc
-   !!==============================================================================
+   !!======================================================================
    !!                       ***  MODULE  trasbc  ***
    !! 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  !  1998-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-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp
    !!----------------------------------------------------------------------
 
@@ -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)
@@ -52,4 +53,6 @@
       !!         at the surface by evaporation, precipitations and runoff (E-P-R); 
       !!      (3) Fwe, tracer carried with the water that is exchanged. 
+      !!            - salinity    : salt flux only due to freezing/melting
+      !!            sa = sa +  fsalt / rau0 / e3t  for k=1
       !!
       !!      Fext, flux through the air-sea interface for temperature and salt: 
@@ -76,9 +79,11 @@
       !!            (Tp P - Te E) + SST (P-E) = 0 when Tp=Te=SST
       !!            - salinity    : evaporation, precipitation and runoff
-      !!         water has a zero salinity (Fwe=0), thus only Fwi remains:
-      !!            sa = sa + emp * sn / e3t   for k=1
+      !!         water has a zero salinity  but there is a salt flux due to 
+      !!         freezing/melting, thus:
+      !!            sa = sa + emp * sn / rau0 / e3t   for k=1
+      !!                    + fsalt    / rau0 / e3t
       !!         where emp, the surface freshwater budget (evaporation minus
       !!         precipitation minus runoff) given in kg/m2/s is divided
-      !!         by 1035 kg/m3 (density of ocena water) to obtain m/s.    
+      !!         by rau0 = 1020 kg/m3 (density of sea water) to obtain m/s.    
       !!         Note: even though Fwe does not appear explicitly for 
       !!         temperature in this routine, the heat carried by the water
@@ -92,5 +97,5 @@
       !!         deal with it in this routine.
       !!           - temperature: Fwe=SST (P-E+R) is added to Fext.
-      !!           - salinity:  Fwe = 0, there is no surface flux of salt.
+      !!           - salinity:  Fwe sa = sa +  fsalt / rau0 / e3t.
       !!
       !! ** Action  : - Update the 1st level of (ta,sa) with the trend associated
@@ -103,6 +108,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) ::   z1_e3t_rau0   ! local scalars
       !!----------------------------------------------------------------------
 
@@ -113,9 +118,4 @@
       ENDIF
 
-      zsrau = 1. / rau0             ! initialization
-#if defined key_zco
-      zse3t = 1. / e3t_0(1)
-#endif
-
       IF( l_trdtra ) THEN           ! Save ta and sa trends
          ztrdt(:,:,:) = ta(:,:,:) 
@@ -123,24 +123,29 @@
       ENDIF
 
+!!gm useless staff ???
       IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration
+!!gm
 
-      ! 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 = r1_rau0_rcp * 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 = r1_rau0_rcp * 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
+      
+      IF( lk_vvl ) THEN          ! Variable Volume Layers case   ===>> heat content of mass flux in qns
+         DO jj = 2, jpj
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               z1_e3t_rau0 = 1./ ( fse3t(ji,jj,1) * rau0 )
+               ta(ji,jj,1) = ta(ji,jj,1) + z1_e3t_rau0 * qns (ji,jj)  * r1_rcp             ! non solar heat flux
+               sa(ji,jj,1) = sa(ji,jj,1) + z1_e3t_rau0 * emps(ji,jj)                       ! salt flux (freezing/melting)
+            END DO
          END DO
-      END DO
+         !
+      ELSE                       ! Constant Volume layers case   ===>> Concentration dillution effect
+         DO jj = 2, jpj
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               z1_e3t_rau0 = 1./ ( fse3t(ji,jj,1) * rau0 )
+               ta(ji,jj,1) = ta(ji,jj,1) + z1_e3t_rau0 * (  qns (ji,jj) * r1_rcp       &   ! non solar heat flux
+                  &                                       + emp (ji,jj) * tn(ji,jj,1)  )   ! concent./dilut. effect
+               sa(ji,jj,1) = sa(ji,jj,1) + z1_e3t_rau0 * (  emps(ji,jj)                &   ! salt flux (freezing/melting)
+                  &                                       + emp (ji,jj) * sn(ji,jj,1)  )   ! concent./dilut. effect
+            END DO
+         END DO
+      ENDIF
 
       IF( l_trdtra ) THEN      ! save the sbc trends for diagnostic
