Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90	(revision 3396)
@@ -19,6 +19,6 @@
    PUBLIC    ice_alloc_2  !  Called in iceini_2.F90
 
-   INTEGER , PUBLIC ::   numit     !: ice iteration index
-   REAL(wp), PUBLIC ::   rdt_ice   !: ice time step
+   INTEGER , PUBLIC ::   numit        !: ice iteration index
+   REAL(wp), PUBLIC ::   rdt_ice      !: ice time step
 
    !                                                                     !!* namicerun read in iceini  *
@@ -98,6 +98,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qstoif        !: Energy stored in the brine pockets
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fbif          !: Heat flux at the ice base
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdmsnif       !: Variation of snow mass
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdmicif       !: Variation of ice mass
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdm_snw       !: Variation of snow mass over 1 time step           [Kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdq_snw       !: Heat content associated with rdm_snw              [J/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdm_ice       !: Variation of ice  mass over 1 time step           [Kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdq_ice       !: Heat content associated with rdm_ice              [J/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qldif         !: heat balance of the lead (or of the open ocean)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qcmif         !: Energy needed to freeze the ocean surface layer
@@ -153,5 +155,7 @@
 
       ALLOCATE(phicif(jpi,jpj) , pfrld  (jpi,jpj) , qstoif (jpi,jpj) ,     &
-         &     fbif  (jpi,jpj) , rdmsnif(jpi,jpj) , rdmicif(jpi,jpj) ,     &
+         &     fbif  (jpi,jpj) , rdm_snw(jpi,jpj) , rdq_snw(jpi,jpj) ,     &
+         &                       rdm_ice(jpi,jpj) , rdq_ice(jpi,jpj) ,     &
+         &     qldif (jpi,jpj) , qcmif  (jpi,jpj) , fdtcn  (jpi,jpj) ,     &
          &     qldif (jpi,jpj) , qcmif  (jpi,jpj) , fdtcn  (jpi,jpj) ,     &
          &     qdtcn (jpi,jpj) , thcm   (jpi,jpj)                    , STAT=ierr(4) )
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 3396)
@@ -9,5 +9,6 @@
    !!            3.3  ! 2009-05 (G. Garric, C. Bricaud) addition of the lim2_evp case
    !!             -   ! 2010-11 (G. Madec) ice-ocean stress computed at each ocean time-step
-   !!            4.0  ! 2011-01  (A. R. Porter, STFC Daresbury) dynamical allocation
+   !!           3.3.1 ! 2011-01 (A. R. Porter, STFC Daresbury) dynamical allocation
+   !!            3.5  ! 2012-11 ((G. Madec, Y. Aksenov, A. Coward) salt and heat fluxes associated with e-p
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -88,8 +89,8 @@
       !!              - Update the fluxes provided to the ocean
       !!     
-      !! ** Outputs : - qsr     : sea heat flux:     solar 
-      !!              - qns     : sea heat flux: non solar
-      !!              - emp     : freshwater budget: volume flux 
-      !!              - emps    : freshwater budget: concentration/dillution 
+      !! ** Outputs : - qsr     : sea heat flux    : solar 
+      !!              - qns     : sea heat flux    : non solar (including heat content of the mass flux)
+      !!              - emp     : freshwater budget: mass flux 
+      !!              - emps    : freshwater budget: salt flux due to Freezing/Melting
       !!              - utau    : sea surface i-stress (ocean referential)
       !!              - vtau    : sea surface j-stress (ocean referential)
@@ -107,6 +108,8 @@
       INTEGER  ::   ifvt, i1mfr, idfr, iflt    !   -       -
       INTEGER  ::   ial, iadv, ifral, ifrdv    !   -       -
-      REAL(wp) ::   zqsr, zqns, zfm            ! local scalars
-      REAL(wp) ::   zinda, zfons, zemp         !   -      -
+      REAL(wp) ::   zqsr,     zqns,   zfmm     ! local scalars
+      REAL(wp) ::   zinda,    zfsalt, zemp     !   -      -
+      REAL(wp) ::   zemp_snw, zqhc,   zcd      !   -      -
+      REAL(wp) ::   zswitch                    !   -      -
       REAL(wp), POINTER, DIMENSION(:,:)   ::   zqnsoce       ! 2D workspace
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp   ! 2D/3D workspace
@@ -116,4 +119,11 @@
       CALL wrk_alloc( jpi, jpj, 1, zalb, zalbp )
 
+      zswitch = 1                                ! Default standard levitating sea-ice (salt exchanges only)
+!!gm ice embedment
+!     SELECT CASE( nn_ice_embd )                 ! levitating/embedded sea-ice option (not yet activated)
+!       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
       !------------------------------------------!
       !      heat flux at the ocean surface      !
@@ -133,22 +143,24 @@
             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.
+!!$            attempt to explain the tricky flags set above....
+!!$            zinda   = 1.0 - AINT( pfrld(ji,jj) )                   ! = 0. if ice-free ocean else 1. (after ice adv, but before ice thermo)
+!!$            i1mfr   = 1.0 - AINT(  frld(ji,jj) )                   ! = 0. if ice-free ocean else 1. (after ice thermo)
+!!$
+!!$            IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda      ! = 1. if there was snow and ice before the ice thermo. which has been completely melted (possibly overmelted)
+!!$            ELSE                             ;   ifvt = 0.         ! 
 !!$            ENDIF
 !!$
-!!$            IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases from previous to current
+!!$            IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases due to ice thermodynamics
 !!$            ELSE                                     ;   idfr = 1.   
 !!$            ENDIF
 !!$
-!!$            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )    !   = 1. if ice (not only snow) at previous and pure ocean at current
+!!$            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )    !   = 1. if ice (not only snow) at previous time and ice-free ocean currently
 !!$
 !!$            ial     = ifvt   * i1mfr    +    ( 1 - ifvt ) * idfr
+!!$                    = i1mfr if ifvt = 1 i.e. 
+!!$                    = idfr  if ifvt = 0
 !!$!                 snow no ice   ice         ice or nothing  lead fraction increases
 !!$!                 at previous   now           at previous
-!!$!                -> ice aera increases  ???         -> ice aera decreases ???
+!!$!                -> ice area increases  ???         -> ice area decreases ???
 !!$
 !!$            iadv    = ( 1  - i1mfr ) * zinda
@@ -174,13 +186,43 @@
 #endif            
             !  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) ) * r1_rdtice    &
-               &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) )                   * r1_rdtice 
-
-            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! ???
+            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) ) * r1_rdtice  &
+               &       + ifrdv   * (       qfvbq(ji,jj) +             qdtcn(ji,jj) ) * r1_rdtice 
+
+            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! store residual heat flux (to put into the ocean at the next time-step)
+            zqhc = ( rdq_snw(ji,jj) + rdq_ice(ji,jj) ) * r1_rdtice       ! heat flux due to snow & ice heat content exchange
             !
             qsr  (ji,jj) = zqsr                                          ! solar heat flux 
-            qns  (ji,jj) = zqns - fdtcn(ji,jj)                           ! non 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_rdtice                    ! snow melting = pure water that enters the ocean
+            zfmm     = rdm_ice(ji,jj) * r1_rdtice                    ! Freezing minus Melting (F-M)
+
+            ! salt flux at the ice/ocean interface (sea ice fraction) [PSU*kg/m2/s]
+            zfsalt = - sice_0(ji,jj) * zfmm                          ! F-M salt exchange
+            zcd    =   soce_0(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
@@ -190,38 +232,8 @@
       CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1.e0 - pfrld(:,:)) )
 
-      !------------------------------------------!
-      !      mass flux at the ocean surface      !
-      !------------------------------------------!
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-#if defined key_coupled
-            ! freshwater exchanges at the ice-atmosphere / ocean interface (coupled mode)
-            zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   ! 
-               &   + rdmsnif(ji,jj) * r1_rdtice                                   !  freshwaterflux due to snow melting 
-#else
-            !  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) )    &   !  change in ice cover within the time step
-               &   + rdmsnif(ji,jj) * r1_rdtice                    !  freshwater flux due to snow melting 
-#endif            
-            !
-            !  computing salt exchanges at the ice/ocean interface
-            zfons = ( soce_0(ji,jj) - sice_0(ji,jj) ) * ( rdmicif(ji,jj) * r1_rdtice ) 
-            !
-            !  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       ! AR5 diagnostics
-         CALL iom_put( 'isnwmlt_cea'  ,                 rdmsnif(:,:) * r1_rdtice )
-         CALL iom_put( 'fsal_virt_cea',   soce_0(:,:) * rdmicif(:,:) * r1_rdtice )
-         CALL iom_put( 'fsal_real_cea', - sice_0(:,:) * rdmicif(:,:) * r1_rdtice )
+         CALL iom_put( 'isnwmlt_cea'  ,                 rdm_snw(:,:) * r1_rdtice )
+         CALL iom_put( 'fsal_virt_cea',   soce_0(:,:) * rdm_ice(:,:) * r1_rdtice )
+         CALL iom_put( 'fsal_real_cea', - sice_0(:,:) * rdm_ice(:,:) * r1_rdtice )
       ENDIF
 
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90	(revision 3396)
@@ -89,6 +89,4 @@
       REAL(wp) ::   za , zh, zthsnice    !
       REAL(wp) ::   zfric_u              ! friction velocity 
-      REAL(wp) ::   zfnsol               ! total non solar heat
-      REAL(wp) ::   zfontn               ! heat flux from snow thickness
       REAL(wp) ::   zfntlat, zpareff     ! test. the val. of lead heat budget
 
@@ -129,5 +127,4 @@
       zdvolif(:,:) = 0.e0   ! total variation of ice volume
       zdvonif(:,:) = 0.e0   ! transformation of snow to sea-ice volume
-!      zdvonif(:,:) = 0.e0   ! lateral variation of ice volume
       zlicegr(:,:) = 0.e0   ! lateral variation of ice volume
       zdvomif(:,:) = 0.e0   ! variation of ice volume at bottom due to melting only
@@ -137,6 +134,8 @@
       ffltbif(:,:) = 0.e0   ! linked with fstric
       qfvbq  (:,:) = 0.e0   ! linked with fstric
-      rdmsnif(:,:) = 0.e0   ! variation of snow mass per unit area
-      rdmicif(:,:) = 0.e0   ! variation of ice mass per unit area
+      rdm_snw(:,:) = 0.e0   ! variation of snow mass over 1 time step
+      rdq_snw(:,:) = 0.e0   ! heat content associated with rdm_snw
+      rdm_ice(:,:) = 0.e0   ! variation of ice mass over 1 time step
+      rdq_ice(:,:) = 0.e0   ! heat content associated with rdm_ice
       zmsk (:,:,:) = 0.e0
 
@@ -199,9 +198,7 @@
       !--------------------------------------------------------------------------
 
-      sst_m(:,:) = sst_m(:,:) + rt0
-
-!CDIR NOVERRCHK
-      DO jj = 1, jpj
-!CDIR NOVERRCHK
+      !CDIR NOVERRCHK
+      DO jj = 1, jpj
+         !CDIR NOVERRCHK
          DO ji = 1, jpi
             zthsnice       = hsnif(ji,jj) + hicif(ji,jj)
@@ -217,24 +214,18 @@
             !  temperature and turbulent mixing (McPhee, 1992)
             zfric_u        = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin )  ! friction velocity
-            fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( sst_m(ji,jj) - tfu(ji,jj) ) 
+            fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( sst_m(ji,jj) + rt0 - tfu(ji,jj) ) 
             qdtcn(ji,jj)  = zindb * fdtcn(ji,jj) * frld(ji,jj) * rdt_ice
                         
             !  partial computation of the lead energy budget (qldif)
 #if defined key_coupled 
-            qldif(ji,jj)   = tms(ji,jj) * rdt_ice                                             &
+            qldif(ji,jj)   = tms(ji,jj) * rdt_ice                                                  &
                &    * (   ( qsr_tot(ji,jj) - qsr_ice(ji,jj,1) * zfricp ) * ( 1.0 - thcm(ji,jj) )   &
                &        + ( qns_tot(ji,jj) - qns_ice(ji,jj,1) * zfricp )                           &
                &        + frld(ji,jj) * ( fdtcn(ji,jj) + ( 1.0 - zindb ) * fsbbq(ji,jj) )   )
 #else
-            zfontn         = ( sprecip(ji,jj) / rhosn ) * xlsn  !   energy for melting solid precipitation
-            zfnsol         = qns(ji,jj)                         !  total non solar flux over the ocean
-            qldif(ji,jj)   = tms(ji,jj) * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )   &
-               &                               + zfnsol + fdtcn(ji,jj) - zfontn     &
-               &                               + ( 1.0 - zindb ) * fsbbq(ji,jj) )   &
-               &                        * frld(ji,jj) * rdt_ice    
-!!$            qldif(ji,jj)   = tms(ji,jj) * rdt_ice * frld(ji,jj) 
-!!$               &           * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )      &
-!!$               &             + qns(ji,jj)  + fdtcn(ji,jj) - zfontn     &
-!!$               &             + ( 1.0 - zindb ) * fsbbq(ji,jj)      )   &
+            qldif(ji,jj)   = tms(ji,jj) * rdt_ice * frld(ji,jj)                    &
+               &                        * (  qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )    &
+               &                           + qns(ji,jj)  +  fdtcn(ji,jj)           &
+               &                           + ( 1.0 - zindb ) * fsbbq(ji,jj)      )
 #endif
             !  parlat : percentage of energy used for lateral ablation (0.0) 
@@ -246,6 +237,5 @@
             
             !  energy needed to bring ocean surface layer until its freezing
-            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1)   &
-                &          * ( tfu(ji,jj) - sst_m(ji,jj) ) * ( 1 - zinda )
+            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1) * ( tfu(ji,jj) - sst_m(ji,jj) - rt0 ) * ( 1 - zinda )
             
             !  calculate oceanic heat flux.
@@ -257,6 +247,4 @@
       END DO
       
-      sst_m(:,:) = sst_m(:,:) - rt0
-               
       !         Select icy points and fulfill arrays for the vectorial grid.
       !----------------------------------------------------------------------
@@ -312,6 +300,9 @@
          CALL tab_2d_1d_2( nbpb, qldif_1d   (1:nbpb)     , qldif      , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d_2( nbpb, qstbif_1d  (1:nbpb)     , qstoif     , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d_2( nbpb, rdmicif_1d (1:nbpb)     , rdmicif    , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d_2( nbpb, rdm_ice_1d (1:nbpb)     , rdm_ice    , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d_2( nbpb, rdq_ice_1d (1:nbpb)     , rdq_ice    , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d_2( nbpb, dmgwi_1d   (1:nbpb)     , dmgwi      , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d_2( nbpb, rdm_snw_1d (1:nbpb)     , rdm_snw    , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d_2( nbpb, rdq_snw_1d (1:nbpb)     , rdq_snw    , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d_2( nbpb, qlbbq_1d   (1:nbpb)     , zqlbsbq    , jpi, jpj, npb(1:nbpb) )
          !
@@ -332,7 +323,9 @@
          CALL tab_1d_2d_2( nbpb, qfvbq      , npb, qfvbq_1d  (1:nbpb)     , jpi, jpj )
          CALL tab_1d_2d_2( nbpb, qstoif     , npb, qstbif_1d (1:nbpb)     , jpi, jpj )
-         CALL tab_1d_2d_2( nbpb, rdmicif    , npb, rdmicif_1d(1:nbpb)     , jpi, jpj )
+         CALL tab_1d_2d_2( nbpb, rdm_ice    , npb, rdm_ice_1d(1:nbpb)     , jpi, jpj )
+         CALL tab_1d_2d_2( nbpb, rdq_ice    , npb, rdq_ice_1d(1:nbpb)     , jpi, jpj )
          CALL tab_1d_2d_2( nbpb, dmgwi      , npb, dmgwi_1d  (1:nbpb)     , jpi, jpj )
-         CALL tab_1d_2d_2( nbpb, rdmsnif    , npb, rdmsnif_1d(1:nbpb)     , jpi, jpj )
+         CALL tab_1d_2d_2( nbpb, rdm_snw    , npb, rdm_snw_1d(1:nbpb)     , jpi, jpj )
+         CALL tab_1d_2d_2( nbpb, rdq_snw    , npb, rdq_snw_1d(1:nbpb)     , jpi, jpj )
          CALL tab_1d_2d_2( nbpb, zdvosif    , npb, dvsbq_1d  (1:nbpb)     , jpi, jpj )
          CALL tab_1d_2d_2( nbpb, zdvobif    , npb, dvbbq_1d  (1:nbpb)     , jpi, jpj )
@@ -393,5 +386,5 @@
       IF( nbpac > 0 ) THEN
          !
-         zlicegr(:,:) = rdmicif(:,:)		! to output the lateral sea-ice growth 
+         zlicegr(:,:) = rdm_ice(:,:)		! to output the lateral sea-ice growth 
          !...Put the variable in a 1-D array for lateral accretion
          CALL tab_2d_1d_2( nbpac, frld_1d   (1:nbpac)     , frld       , jpi, jpj, npac(1:nbpac) )
@@ -404,5 +397,6 @@
          CALL tab_2d_1d_2( nbpac, qcmif_1d  (1:nbpac)     , qcmif      , jpi, jpj, npac(1:nbpac) )
          CALL tab_2d_1d_2( nbpac, qstbif_1d (1:nbpac)     , qstoif     , jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d_2( nbpac, rdmicif_1d(1:nbpac)     , rdmicif    , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d_2( nbpac, rdm_ice_1d(1:nbpac)     , rdm_ice    , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d_2( nbpac, rdq_ice_1d(1:nbpac)     , rdq_ice    , jpi, jpj, npac(1:nbpac) )
          CALL tab_2d_1d_2( nbpac, dvlbq_1d  (1:nbpac)     , zdvolif    , jpi, jpj, npac(1:nbpac) )
          CALL tab_2d_1d_2( nbpac, tfu_1d    (1:nbpac)     , tfu        , jpi, jpj, npac(1:nbpac) )
@@ -418,5 +412,6 @@
          CALL tab_1d_2d_2( nbpac, tbif(:,:,3), npac(1:nbpac), tbif_1d   (1:nbpac , 3 ), jpi, jpj )
          CALL tab_1d_2d_2( nbpac, qstoif     , npac(1:nbpac), qstbif_1d (1:nbpac)     , jpi, jpj )
-         CALL tab_1d_2d_2( nbpac, rdmicif    , npac(1:nbpac), rdmicif_1d(1:nbpac)     , jpi, jpj )
+         CALL tab_1d_2d_2( nbpac, rdm_ice    , npac(1:nbpac), rdm_ice_1d(1:nbpac)     , jpi, jpj )
+         CALL tab_1d_2d_2( nbpac, rdq_ice    , npac(1:nbpac), rdq_ice_1d(1:nbpac)     , jpi, jpj )
          CALL tab_1d_2d_2( nbpac, zdvolif    , npac(1:nbpac), dvlbq_1d  (1:nbpac)     , jpi, jpj )
          !
@@ -449,11 +444,11 @@
       CALL iom_put( 'iceprod_cea' , hicifp (:,:) * zztmp     )   ! Ice produced               [m/s]
       IF( lk_diaar5 ) THEN
-         CALL iom_put( 'snowmel_cea' , rdmsnif(:,:) * zztmp     )   ! Snow melt                  [kg/m2/s]
+         CALL iom_put( 'snowmel_cea' , rdm_snw(:,:) * zztmp     )   ! Snow melt                  [kg/m2/s]
          zztmp = rhoic / rdt_ice
          CALL iom_put( 'sntoice_cea' , zdvonif(:,:) * zztmp     )   ! Snow to Ice transformation [kg/m2/s]
          CALL iom_put( 'ticemel_cea' , zdvosif(:,:) * zztmp     )   ! Melt at Sea Ice top        [kg/m2/s]
          CALL iom_put( 'bicemel_cea' , zdvomif(:,:) * zztmp     )   ! Melt at Sea Ice bottom     [kg/m2/s]
-         zlicegr(:,:) = MAX( 0.e0, rdmicif(:,:)-zlicegr(:,:) )
-         CALL iom_put( 'licepro_cea' , zlicegr(:,:) * zztmp     )   ! Latereal sea ice growth    [kg/m2/s]
+         zlicegr(:,:) = MAX( 0.e0, rdm_ice(:,:)-zlicegr(:,:) )
+         CALL iom_put( 'licepro_cea' , zlicegr(:,:) * zztmp     )   ! Lateral sea ice growth     [kg/m2/s]
       ENDIF
       !
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limthd_lac_2.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limthd_lac_2.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limthd_lac_2.F90	(revision 3396)
@@ -145,6 +145,6 @@
          frld_1d   (ji) = MAX( zfrlnew , zfrlmin(ji) )
          !--computation of the remaining part of ice thickness which has been already used
-         zdhicbot(ji) =  ( frld_1d(ji) - zfrlnew ) * zhice0(ji) / ( 1.0 - zfrlmin(ji) ) & 
-                      -  (  ( 1.0 - zfrrate ) / ( 1.0 - frld_1d(ji) ) )  * ( zqbgow(ji) / xlic ) 
+         zdhicbot(ji) =  ( frld_1d(ji) - zfrlnew ) * zhice0(ji) / ( 1.0 - zfrlmin(ji) )   & 
+            &         -  (  ( 1.0 - zfrrate ) / ( 1.0 - frld_1d(ji) ) )  * ( zqbgow(ji) / xlic ) 
       END DO
  
@@ -196,5 +196,5 @@
             &          ) / zah
          
-         tbif_1d(ji,3) =     (  iiceform * ( zhnews2 - zdh3 )                                          * zta1  &
+         tbif_1d(ji,3) =     ( iiceform * ( zhnews2 - zdh3 )                                           * zta1  &
             &              + ( iiceform * zdh3 + ( 1 - iiceform ) * zdh1 )                             * zta2  &
             &              + ( iiceform * ( zhnews2 - zdh5 ) + ( 1 - iiceform ) * ( zhnews2 - zdh1 ) ) * zta3  & 
@@ -217,5 +217,6 @@
       DO ji = kideb , kiut
          dvlbq_1d  (ji) = ( 1. - frld_1d(ji) ) * h_ice_1d(ji) - ( 1. - zfrl_old(ji) ) * zhice_old(ji)
-         rdmicif_1d(ji) = rdmicif_1d(ji) + rhoic * dvlbq_1d(ji)
+         rdm_ice_1d(ji) = rdm_ice_1d(ji) + rhoic * dvlbq_1d(ji)
+         rdq_ice_1d(ji) = rdq_ice_1d(ji) + rcpic * dvlbq_1d(ji) * ( tfu_1d(ji) - rt0 )      ! heat content relative to rt0
       END DO
       
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limthd_zdf_2.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limthd_zdf_2.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limthd_zdf_2.F90	(revision 3396)
@@ -86,5 +86,5 @@
       REAL(wp), POINTER, DIMENSION(:) ::   zrcpdt         ! h_su*rho_su*cp_su/dt(h_su being the thick. of surf. layer)
       REAL(wp), POINTER, DIMENSION(:) ::   zts_old        ! previous surface temperature
-      REAL(wp), POINTER, DIMENSION(:) ::   zidsn , z1midsn , zidsnic ! tempory variables
+      REAL(wp), POINTER, DIMENSION(:) ::   zidsn , z1midsn , zidsnic ! temporary variables
       REAL(wp), POINTER, DIMENSION(:) ::   zfnet          ! net heat flux at the top surface( incl. conductive heat flux)
       REAL(wp), POINTER, DIMENSION(:) ::   zsprecip       ! snow accumulation
@@ -98,62 +98,68 @@
       REAL(wp), POINTER, DIMENSION(:) ::   zep            ! internal temperature of the 2nd layer of the snow/ice system
       REAL(wp), DIMENSION(3) :: & 
-          zplediag  &    ! principle diagonal, subdiag. and supdiag. of the 
+            zplediag  &    ! principle diagonal, subdiag. and supdiag. of the 
           , zsubdiag  &    ! tri-diagonal matrix coming from the computation
           , zsupdiag  &    ! of the temperatures inside the snow-ice system
           , zsmbr          ! second member
-       REAL(wp) :: & 
-          zhsu     &     ! thickness of surface layer
-          , zhe      &     ! effective thickness for compu. of equ. thermal conductivity
-          , zheshth  &     ! = zhe / thth
-          , zghe     &     ! correction factor of the thermal conductivity
-          , zumsb    &     ! parameter for numerical method to solve heat-diffusion eq.
-          , zkhsn    &     ! conductivity at the snow layer
-          , zkhic    &     ! conductivity at the ice layers
-          , zkint    &     ! equivalent conductivity at the snow-ice interface
-          , zkhsnint &     ! = zkint*dt / (hsn*rhosn*cpsn)  
-          , zkhicint &     ! = 2*zkint*dt / (hic*rhoic*cpic)
-          , zpiv1 , zpiv2  &       ! tempory scalars used to solve the tri-diagonal system
-          , zb2 , zd2 , zb3 , zd3 &
+       REAL(wp) ::    &
+            zhsu      &    ! thickness of surface layer
+          , zhe       &    ! effective thickness for compu. of equ. thermal conductivity
+          , zheshth   &    ! = zhe / thth
+          , zghe      &    ! correction factor of the thermal conductivity
+          , zumsb     &    ! parameter for numerical method to solve heat-diffusion eq.
+          , zkhsn     &    ! conductivity at the snow layer
+          , zkhic     &    ! conductivity at the ice layers
+          , zkint     &    ! equivalent conductivity at the snow-ice interface
+          , zkhsnint  &    ! = zkint*dt / (hsn*rhosn*cpsn)  
+          , zkhicint  &    ! = 2*zkint*dt / (hic*rhoic*cpic)
+          , zpiv1, zpiv2 & ! temporary scalars used to solve the tri-diagonal system
+          , zb2, zd2  &    ! temporary scalars used to solve the tri-diagonal system
+          , zb3, zd3  &    ! temporary scalars used to solve the tri-diagonal system
           , ztint          ! equivalent temperature at the snow-ice interface
-       REAL(wp) :: & 
-          zexp      &     ! exponential function of the ice thickness
-          , zfsab     &     ! part of solar radiation stored in brine pockets
-          , zfts      &     ! value of energy balance function when the temp. equal surf. temp.
-          , zdfts     &     ! value of derivative of ztfs when the temp. equal surf. temp.
-          , zdts      &     ! surface temperature increment
-          , zqsnw_mlt &     ! energy needed to melt snow
-          , zdhsmlt   &     ! change in snow thickness due to melt
-          , zhsn      &     ! snow thickness (previous+accumulation-melt)
-          , zqsn_mlt_rem &  ! remaining heat coming from snow melting
-          , zqice_top_mlt & ! energy used to melt ice at top surface
-          , zdhssub      &  ! change in snow thick. due to sublimation or evaporation
-          , zdhisub      &  ! change in ice thick. due to sublimation or evaporation    
-          , zdhsn        &  ! snow ice thickness increment
-          , zdtsn        &  ! snow internal temp. increment
-          , zdtic        &  ! ice internal temp. increment
+       REAL(wp) ::    & 
+            zexp      &    ! exponential function of the ice thickness
+          , zfsab     &    ! part of solar radiation stored in brine pockets
+          , zfts      &    ! value of energy balance function when the temp. equal surf. temp.
+          , zdfts     &    ! value of derivative of ztfs when the temp. equal surf. temp.
+          , zdts      &    ! surface temperature increment
+          , zqsnw_mlt &    ! energy needed to melt snow
+          , zdhsmlt   &    ! change in snow thickness due to melt
+          , zhsn      &    ! snow thickness (previous+accumulation-melt)
+          , zqsn_mlt_rem & ! remaining heat coming from snow melting
+          , zqice_top_mlt &! energy used to melt ice at top surface
+          , zdhssub   &    ! change in snow thick. due to sublimation or evaporation
+          , zdhisub   &    ! change in ice thick. due to sublimation or evaporation    
+          , zdhsn     &    ! snow ice thickness increment
+          , zdtsn     &    ! snow internal temp. increment
+          , zdtic     &    ! ice internal temp. increment
           , zqnes          ! conductive energy due to ice melting in the first ice layer
-       REAL(wp) :: & 
-          ztbot     &      ! temperature at the bottom surface
-          , zfcbot    &      ! conductive heat flux at bottom surface
-          , zqice_bot &      ! energy used for bottom melting/growing
-          , zqice_bot_mlt &  ! energy used for bottom melting
-          , zqstbif_bot  &  ! part of energy stored in brine pockets used for bottom melting
-          , zqstbif_old  &  ! tempory var. for zqstbif_bot
-          , zdhicmlt      &  ! change in ice thickness due to bottom melting
-          , zdhicm        &  ! change in ice thickness var. 
-          , zdhsnm        &  ! change in snow thickness var. 
-          , zhsnfi        &  ! snow thickness var. 
-          , zc1, zpc1, zc2, zpc2, zp1, zp2 & ! tempory variables
-          , ztb2, ztb3
-       REAL(wp) :: & 
-          zdrmh         &   ! change in snow/ice thick. after snow-ice formation
-          , zhicnew       &   ! new ice thickness
-          , zhsnnew       &   ! new snow thickness
-          , zquot , ztneq &   ! tempory temp. variables
-          , zqice, zqicetot & ! total heat inside the snow/ice system
-          , zdfrl         &   ! change in ice concentration
-          , zdvsnvol      &   ! change in snow volume
-          , zdrfrl1, zdrfrl2 &  ! tempory scalars
-          , zihsn, zidhb, zihic, zihe, zihq, ziexp, ziqf, zihnf, zibmlt, ziqr, zihgnew, zind
+       REAL(wp) ::    & 
+            ztbot     &    ! temperature at the bottom surface
+          , zfcbot    &    ! conductive heat flux at bottom surface
+          , zqice_bot &    ! energy used for bottom melting/growing
+          , zqice_bot_mlt &! energy used for bottom melting
+          , zqstbif_bot  & ! part of energy stored in brine pockets used for bottom melting
+          , zqstbif_old  & ! temporary var. for zqstbif_bot
+          , zdhicmlt  &    ! change in ice thickness due to bottom melting
+          , zdhicm    &    ! change in ice thickness var. 
+          , zdhsnm    &    ! change in snow thickness var. 
+          , zhsnfi    &    ! snow thickness var. 
+          , zc1, zpc1 &    ! temporary variables
+          , zc2, zpc2 &    ! temporary variables
+          , zp1, zp2  &    ! temporary variables
+          , ztb2, ztb3     ! temporary variables
+       REAL(wp) ::    & 
+            zdrmh     &    ! change in snow/ice thick. after snow-ice formation
+          , zhicnew   &    ! new ice thickness
+          , zhsnnew   &    ! new snow thickness
+          , zquot     &
+          , ztneq     &    ! temporary temp. variables
+          , zqice     &
+          , zqicetot  &    ! total heat inside the snow/ice system
+          , zdfrl     &    ! change in ice concentration
+          , zdvsnvol  &    ! change in snow volume
+          , zdrfrl1, zdrfrl2, zihsn, zidhb, zihic &  ! temporary scalars
+          , zihe, zihq, ziexp, ziqf, zihnf        &  ! temporary scalars
+          , zibmlt, ziqr, zihgnew, zind, ztmp        ! temporary scalars
        !!----------------------------------------------------------------------
        CALL wrk_alloc( jpij, ztsmlt, ztbif  , zksn    , zkic    , zksndh , zfcsu  , zfcsudt , zi0      , z1mi0   , zqmax    )
@@ -169,16 +175,14 @@
        
        DO ji = kideb , kiut
+          ! do nothing if the snow (ice) thickness falls below its minimum thickness
           zihsn = MAX( zzero , SIGN( zone , hsndif - h_snow_1d(ji) ) )
           zihic = MAX( zzero , SIGN( zone , hicdif - h_ice_1d(ji) ) )
-          !--computation of energy due to surface melting
-          zqcmlts(ji) = ( MAX ( zzero ,  &
-             &                   rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn )
-          !--computation of energy due to bottom melting
-          zqcmltb(ji) = ( MAX( zzero , &
-             &                  rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) &
-             &           + MAX( zzero , &
-             &                  rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) &
-             &           ) * ( 1.0 - zihic  )
-          !--limitation of  snow/ice system internal temperature
+          !--energy required to bring snow to its melting point (rt0_snow)
+          zqcmlts(ji) = ( MAX ( zzero , rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn )
+          !--energy required to bring ice to its melting point (rt0_ice)
+          zqcmltb(ji) = ( MAX( zzero , rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) )  &
+             &          + MAX( zzero , rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) )  &
+             &          ) * ( 1.0 - zihic  )
+          !--limitation of snow/ice system internal temperature
           tbif_1d(ji,1)   = MIN( rt0_snow, tbif_1d(ji,1) )
           tbif_1d(ji,2)   = MIN( rt0_ice , tbif_1d(ji,2) )
@@ -480,5 +484,8 @@
           dvsbq_1d(ji) =  ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnw_old(ji) - zsprecip(ji) )
           dvsbq_1d(ji) =  MIN( zzero , dvsbq_1d(ji) )
-          rdmsnif_1d(ji) =  rhosn * dvsbq_1d(ji)
+          ztmp = rhosn * dvsbq_1d(ji)
+          rdm_snw_1d(ji) =  ztmp
+          !--heat content of the water provided to the ocean (referenced to rt0)
+          rdq_snw_1d(ji) =  cpic * ztmp * ( rt0_snow - rt0 )
           !-- If the snow is completely melted the remaining heat is used to melt ice
           zqsn_mlt_rem  = MAX( zzero , -zhsn ) * xlsn
@@ -623,9 +630,17 @@
           !---updating new ice thickness and computing the newly formed ice mass
           zhicnew   =  zihgnew * zhicnew
-          rdmicif_1d(ji) =  rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) ) * rhoic
+          ztmp    =  ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) ) * rhoic
+          rdm_ice_1d(ji) =  rdm_ice_1d(ji) + ztmp
+          !---heat content of the water provided to the ocean (referenced to rt0)
+          !   use of rt0_ice is OK for melting ice; in the case of freezing, tfu_1d should be used. 
+          !   This is done in 9.5 section (see below)
+          rdq_ice_1d(ji) =  cpic * ztmp * ( rt0_ice - rt0 )
           !---updating new snow thickness and computing the newly formed snow mass
           zhsnfi   = zhsn + zdhsnm
           h_snow_1d(ji) = MAX( zzero , zhsnfi )
-          rdmsnif_1d(ji) =  rdmsnif_1d(ji) + ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsn ) * rhosn
+          ztmp = ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsn ) * rhosn
+          rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp
+          !---updating the heat content of the water provided to the ocean (referenced to rt0)
+          rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( rt0_snow - rt0 )
           !--remaining energy in case of total ablation
           zqocea(ji) = - ( zihsn * xlic * zdhicm + xlsn * ( zhsnfi - h_snow_1d(ji) ) ) * ( 1.0 - frld_1d(ji) )
@@ -659,4 +674,8 @@
           tbif_1d(ji,3) =  zihgnew * ztb3 + ( 1.0 - zihgnew ) * tfu_1d(ji)
           h_ice_1d(ji)  =  zhicnew
+          ! update the ice heat content given to the ocean in freezing case 
+          ! (part due to difference between rt0_ice and tfu_1d)
+          ztmp = ( 1. - zidhb ) * rhoic * dvbbq_1d(ji)
+          rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * ( tfu_1d(ji) - rt0_ice )
        END DO
 
@@ -700,6 +719,11 @@
           dmgwi_1d(ji) = dmgwi_1d(ji) + ( 1.0 -frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnnew ) * rhosn
           !---  volume change of ice and snow (used for ocean-ice freshwater flux computation)
-          rdmicif_1d(ji) = rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhicnew - h_ice_1d (ji) ) * rhoic
-          rdmsnif_1d(ji) = rdmsnif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhsnnew - h_snow_1d(ji) ) * rhosn
+          ztmp = ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d (ji) ) * rhoic
+          rdm_ice_1d(ji) = rdm_ice_1d(ji) + ztmp
+          rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * ( tfu_1d(ji) - rt0 )
+          !!gm BUG ??   snow ==>  only needed for nn_ice_embd == 0  (standard levitating sea-ice)
+          ztmp = ( 1.0 - frld_1d(ji) )   * ( zhsnnew - h_snow_1d(ji) ) * rhosn
+          rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp
+          rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( rt0_snow - rt0 )
 
           !---  Actualize new snow and ice thickness.
@@ -748,8 +772,28 @@
           !--variation of ice volume and ice mass 
           dvlbq_1d(ji)   = zihic * ( zfrl_old(ji) - frld_1d(ji) ) * h_ice_1d(ji)
-          rdmicif_1d(ji) = rdmicif_1d(ji) + dvlbq_1d(ji) * rhoic
+          ztmp = dvlbq_1d(ji) * rhoic
+          rdm_ice_1d(ji) = rdm_ice_1d(ji) + ztmp
+!!gm
+!!gm   This should be split in two parts:
+!!gm         1-  heat required to bring sea-ice to tfu  : this part should be added to the heat flux taken from the ocean
+!!gm                 cpic * ztmp * 0.5 * ( tbif_1d(ji,2) + tbif_1d(ji,3) - 2.* rt0_ice )
+!!gm         2-  heat content of lateral ablation referenced to rt0 : this part only put in rdq_ice_1d
+!!gm                 cpic * ztmp * ( rt0_ice - rt0 )
+!!gm   Currently we put all the heat in rdq_ice_1d
+          rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * 0.5 * ( tbif_1d(ji,2) + tbif_1d(ji,3) - 2.* rt0 )
+          !
           !--variation of snow volume and snow mass 
-          zdvsnvol    = zihsn * ( zfrl_old(ji) - frld_1d(ji) ) * h_snow_1d(ji)
-          rdmsnif_1d(ji) = rdmsnif_1d(ji) + zdvsnvol * rhosn
+          zdvsnvol = zihsn * ( zfrl_old(ji) - frld_1d(ji) ) * h_snow_1d(ji)
+          ztmp     = zdvsnvol * rhosn
+          rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp
+!!gm
+!!gm   This should be split in two parts:
+!!gm         1-  heat required to bring snow to tfu  : this part should be added to the heat flux taken from the ocean
+!!gm                 cpic * ztmp * ( tbif_1d(ji,1) - rt0_snow )
+!!gm         2-  heat content of lateral ablation referenced to rt0 : this part only put in rdq_snw_1d
+!!gm                 cpic * ztmp * ( rt0_snow - rt0 )
+!!gm   Currently we put all the heat in rdq_snw_1d
+          rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( tbif_1d(ji,1) - rt0 )
+
           h_snow_1d(ji)  = ziqf * h_snow_1d(ji)
 
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/thd_ice_2.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/thd_ice_2.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/thd_ice_2.F90	(revision 3396)
@@ -68,6 +68,8 @@
       qstbif_1d   ,     &  !:    "                  "      qstoif
       fbif_1d     ,     &  !:    "                  "      fbif
-      rdmicif_1d  ,     &  !:    "                  "      rdmicif
-      rdmsnif_1d  ,     &  !:    "                  "      rdmsnif
+      rdm_ice_1d  ,     &  !:    "                  "      rdm_ice
+      rdq_ice_1d  ,     &  !:    "                  "      rdq_ice
+      rdm_snw_1d  ,     &  !:    "                  "      rdm_snw
+      rdq_snw_1d  ,     &  !:    "                  "      rdq_snw
       qlbbq_1d    ,     &  !:    "                  "      qlbsbq
       dmgwi_1d    ,     &  !:    "                  "      dmgwi
@@ -108,5 +110,6 @@
          &      qstbif_1d(jpij),  fbif_1d(jpij),  Stat=ierr(2))
          !
-      ALLOCATE( rdmicif_1d(jpij), rdmsnif_1d(jpij), qlbbq_1d(jpij),   &
+      ALLOCATE( rdm_ice_1d(jpij), rdq_ice_1d(jpij)                  , &
+         &      rdm_snw_1d(jpij), rdq_snw_1d(jpij), qlbbq_1d(jpij)  , &
          &      dmgwi_1d(jpij)  , dvsbq_1d(jpij)  , rdvomif_1d(jpij), &
          &      dvbbq_1d(jpij)  , dvlbq_1d(jpij)  , dvnbq_1d(jpij)  , &
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 3396)
@@ -264,6 +264,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   phicif      !: Old ice thickness
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fbif        !: Heat flux at the ice base
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdmsnif     !: Variation of snow mass
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdmicif     !: Variation of ice mass
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdm_snw     !: Variation of snow mass over 1 time step     [Kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdq_snw     !: Heat content associated with rdm_snw        [J/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdm_ice     !: Variation of ice mass over 1 time step      [Kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdq_ice     !: Heat content associated with rdm_ice        [J/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qldif       !: heat balance of the lead (or of the open ocean)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qcmif       !: Energy needed to bring the ocean to freezing 
@@ -448,5 +450,6 @@
          &      sist     (jpi,jpj) , icethi (jpi,jpj) , t_bo     (jpi,jpj) , hicifp   (jpi,jpj) ,     &
          &      frld     (jpi,jpj) , pfrld  (jpi,jpj) , phicif   (jpi,jpj) , fbif     (jpi,jpj) ,     &
-         &      rdmsnif  (jpi,jpj) , rdmicif(jpi,jpj) , qldif    (jpi,jpj) , qcmif    (jpi,jpj) ,     &
+         &      rdm_snw  (jpi,jpj) , rdq_snw(jpi,jpj) , rdm_ice  (jpi,jpj) , rdq_ice  (jpi,jpj) ,     &
+         &                                              qldif    (jpi,jpj) , qcmif    (jpi,jpj) ,     &
          &      fdtcn    (jpi,jpj) , qdtcn  (jpi,jpj) , fstric   (jpi,jpj) , fscmbq   (jpi,jpj) ,     &
          &      ffltbif  (jpi,jpj) , fsbbq  (jpi,jpj) , qfvbq    (jpi,jpj) , dmgwi    (jpi,jpj) ,     &
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 3396)
@@ -83,7 +83,7 @@
       z_frc_trd_s =           SUM( sbc_tsc(:,:,jp_sal) * surf(:,:) )     ! salt fluxes
       ! Add penetrative solar radiation
-      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + ro0cpr * SUM( qsr     (:,:) * surf(:,:) )
+      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * SUM( qsr     (:,:) * surf(:,:) )
       ! Add geothermal heat flux
-      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t + ro0cpr * SUM( qgh_trd0(:,:) * surf(:,:) )
+      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * SUM( qgh_trd0(:,:) * surf(:,:) )
       IF( lk_mpp ) THEN
          CALL mpp_sum( z_frc_trd_v )
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 3396)
@@ -27,62 +27,65 @@
    REAL(wp), PUBLIC ::   rsmall = 0.5 * EPSILON( 1.e0 )         !: smallest real computer value
    
-   REAL(wp), PUBLIC ::   rday = 24.*60.*60.       !: day (s)
-   REAL(wp), PUBLIC ::   rsiyea                   !: sideral year (s)
-   REAL(wp), PUBLIC ::   rsiday                   !: sideral day (s)
-   REAL(wp), PUBLIC ::   raamo =  12._wp          !: number of months in one year
-   REAL(wp), PUBLIC ::   rjjhh =  24._wp          !: number of hours in one day
-   REAL(wp), PUBLIC ::   rhhmm =  60._wp          !: number of minutes in one hour
-   REAL(wp), PUBLIC ::   rmmss =  60._wp          !: number of seconds in one minute
-!! REAL(wp), PUBLIC ::   omega = 7.292115083046061e-5_wp ,  &  !: change the last digit!
-   REAL(wp), PUBLIC ::   omega                    !: earth rotation parameter
-   REAL(wp), PUBLIC ::   ra    = 6371229._wp      !: earth radius (meter)
-   REAL(wp), PUBLIC ::   grav  = 9.80665_wp       !: gravity (m/s2)
-   
-   REAL(wp), PUBLIC ::   rtt      = 273.16_wp     !: triple point of temperature (Kelvin)
-   REAL(wp), PUBLIC ::   rt0      = 273.15_wp     !: freezing point of water (Kelvin)
+   REAL(wp), PUBLIC ::   rday = 24.*60.*60.     !: day                                [s]
+   REAL(wp), PUBLIC ::   rsiyea                 !: sideral year                       [s]
+   REAL(wp), PUBLIC ::   rsiday                 !: sideral day                        [s]
+   REAL(wp), PUBLIC ::   raamo =  12._wp        !: number of months in one year
+   REAL(wp), PUBLIC ::   rjjhh =  24._wp        !: number of hours in one day
+   REAL(wp), PUBLIC ::   rhhmm =  60._wp        !: number of minutes in one hour
+   REAL(wp), PUBLIC ::   rmmss =  60._wp        !: number of seconds in one minute
+   REAL(wp), PUBLIC ::   omega                  !: earth rotation parameter           [s-1]
+   REAL(wp), PUBLIC ::   ra    = 6371229._wp    !: earth radius                       [m]
+   REAL(wp), PUBLIC ::   grav  = 9.80665_wp     !: gravity                            [m/s2]
+   
+   REAL(wp), PUBLIC ::   rtt      = 273.16_wp        !: triple point of temperature   [Kelvin]
+   REAL(wp), PUBLIC ::   rt0      = 273.15_wp        !: freezing point of fresh water [Kelvin]
 #if defined key_lim3
-   REAL(wp), PUBLIC ::   rt0_snow = 273.16_wp     !: melting point of snow  (Kelvin)
-   REAL(wp), PUBLIC ::   rt0_ice  = 273.16_wp     !: melting point of ice   (Kelvin)
-#else
-   REAL(wp), PUBLIC ::   rt0_snow = 273.15_wp     !: melting point of snow  (Kelvin)
-   REAL(wp), PUBLIC ::   rt0_ice  = 273.05_wp     !: melting point of ice   (Kelvin)
-#endif
-
+   REAL(wp), PUBLIC ::   rt0_snow = 273.16_wp        !: melting point of snow         [Kelvin]
+   REAL(wp), PUBLIC ::   rt0_ice  = 273.16_wp        !: melting point of ice          [Kelvin]
+#else
+   REAL(wp), PUBLIC ::   rt0_snow = 273.15_wp        !: melting point of snow         [Kelvin]
+   REAL(wp), PUBLIC ::   rt0_ice  = 273.05_wp        !: melting point of ice          [Kelvin]
+#endif
 #if defined key_cice
-   REAL(wp), PUBLIC ::   rau0     = 1026._wp      !: reference volumic mass (density)  (kg/m3)
-#else
-   REAL(wp), PUBLIC ::   rau0     = 1035._wp      !: reference volumic mass (density)  (kg/m3)
-#endif
-   REAL(wp), PUBLIC ::   rau0r                    !: reference specific volume         (m3/kg)
-   REAL(wp), PUBLIC ::   rcp      =    4.e+3_wp   !: ocean specific heat
-   REAL(wp), PUBLIC ::   ro0cpr                   !: = 1. / ( rau0 * rcp )
+   REAL(wp), PUBLIC ::   rau0     = 1026._wp         !: volumic mass of reference     [kg/m3]
+#else
+   REAL(wp), PUBLIC ::   rau0     = 1035._wp         !: volumic mass of reference     [kg/m3]
+#endif
+   REAL(wp), PUBLIC ::   r1_rau0                     !: = 1. / rau0                   [m3/kg]
+   REAL(wp), PUBLIC ::   rauw     = 1000._wp         !: volumic mass of pure water    [m3/kg]
+   REAL(wp), PUBLIC ::   rcp      =    4.e3_wp       !: ocean specific heat           [J/Kelvin]
+   REAL(wp), PUBLIC ::   r1_rcp                      !: = 1. / rcp                    [Kelvin/J]
+   REAL(wp), PUBLIC ::   r1_rau0_rcp                 !: = 1. / ( rau0 * rcp )
+
+   REAL(wp), PUBLIC ::   rhosn    =  330._wp         !: volumic mass of snow          [kg/m3]
+   REAL(wp), PUBLIC ::   emic     =    0.97_wp       !: emissivity of snow or ice
+   REAL(wp), PUBLIC ::   sice     =    6.0_wp        !: salinity of ice               [psu]
+   REAL(wp), PUBLIC ::   soce     =   34.7_wp        !: salinity of sea               [psu]
+   REAL(wp), PUBLIC ::   cevap    =    2.5e+6_wp     !: latent heat of evaporation (water)
+   REAL(wp), PUBLIC ::   srgamma  =    0.9_wp        !: correction factor for solar radiation (Oberhuber, 1974)
+   REAL(wp), PUBLIC ::   vkarmn   =    0.4_wp        !: von Karman constant
+   REAL(wp), PUBLIC ::   stefan   =    5.67e-8_wp    !: Stefan-Boltzmann constant 
 
 #if defined key_lim3 || defined key_cice
-   REAL(wp), PUBLIC ::   rcdsn   =   0.31_wp      !: thermal conductivity of snow
-   REAL(wp), PUBLIC ::   rcdic   =   2.034396_wp  !: thermal conductivity of fresh ice
-   REAL(wp), PUBLIC ::   cpic    = 2067.0         !: specific heat of sea ice
-   REAL(wp), PUBLIC ::   lsub    = 2.834e+6       !: pure ice latent heat of sublimation (J.kg-1)
-   REAL(wp), PUBLIC ::   lfus    = 0.334e+6       !: latent heat of fusion of fresh ice   (J.kg-1)
-   REAL(wp), PUBLIC ::   rhoic   = 917._wp        !: volumic mass of sea ice (kg/m3)
-   REAL(wp), PUBLIC ::   tmut    =   0.054        !: decrease of seawater meltpoint with salinity
-#else
-   REAL(wp), PUBLIC ::   rcdsn   =   0.22_wp      !: conductivity of the snow
-   REAL(wp), PUBLIC ::   rcdic   =   2.034396_wp  !: conductivity of the ice
-   REAL(wp), PUBLIC ::   rcpsn   =   6.9069e+5_wp !: density times specific heat for snow
-   REAL(wp), PUBLIC ::   rcpic   =   1.8837e+6_wp !: volumetric latent heat fusion of sea ice
-   REAL(wp), PUBLIC ::   lfus    =   0.3337e+6    !: latent heat of fusion of fresh ice   (J.kg-1) 	
-   REAL(wp), PUBLIC ::   xlsn    = 110.121e+6_wp  !: volumetric latent heat fusion of snow
-   REAL(wp), PUBLIC ::   xlic    = 300.33e+6_wp   !: volumetric latent heat fusion of ice
-   REAL(wp), PUBLIC ::   xsn     =   2.8e+6       !: latent heat of sublimation of snow
-   REAL(wp), PUBLIC ::   rhoic   = 900._wp        !: volumic mass of sea ice (kg/m3)
-#endif
-   REAL(wp), PUBLIC ::   rhosn   = 330._wp        !: volumic mass of snow (kg/m3)
-   REAL(wp), PUBLIC ::   emic    =   0.97_wp      !: emissivity of snow or ice
-   REAL(wp), PUBLIC ::   sice    =   6.0_wp       !: reference salinity of ice (psu)
-   REAL(wp), PUBLIC ::   soce    =  34.7_wp       !: reference salinity of sea (psu)
-   REAL(wp), PUBLIC ::   cevap   =   2.5e+6_wp    !: latent heat of evaporation (water)
-   REAL(wp), PUBLIC ::   srgamma =   0.9_wp       !: correction factor for solar radiation (Oberhuber, 1974)
-   REAL(wp), PUBLIC ::   vkarmn  =   0.4_wp       !: von Karman constant
-   REAL(wp), PUBLIC ::   stefan  =   5.67e-8_wp   !: Stefan-Boltzmann constant 
+   REAL(wp), PUBLIC ::   rhoic    =  917._wp         !: volumic mass of sea ice                               [kg/m3]
+   REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: thermal conductivity of fresh ice
+   REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow
+   REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat for ice 
+   REAL(wp), PUBLIC ::   lsub     =    2.834e+6_wp   !: pure ice latent heat of sublimation                   [J/kg]
+   REAL(wp), PUBLIC ::   lfus     =    0.334e+6_wp   !: latent heat of fusion of fresh ice                    [J/kg]
+   REAL(wp), PUBLIC ::   tmut     =    0.054_wp      !: decrease of seawater meltpoint with salinity
+   REAL(wp), PUBLIC ::   xlsn                        !: = lfus*rhosn (volumetric latent heat fusion of snow)  [J/m3]
+#else
+   REAL(wp), PUBLIC ::   rhoic    =  900._wp         !: volumic mass of sea ice                               [kg/m3]
+   REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: conductivity of the ice                               [W/m/K]
+   REAL(wp), PUBLIC ::   rcpic    =    1.8837e+6_wp  !: volumetric specific heat for ice                      [J/m3/K]
+   REAL(wp), PUBLIC ::   cpic                        !: = rcpic / rhoic  (specific heat for ice)              [J/Kg/K]
+   REAL(wp), PUBLIC ::   rcdsn    =    0.22_wp       !: conductivity of the snow                              [W/m/K]
+   REAL(wp), PUBLIC ::   rcpsn    =    6.9069e+5_wp  !: volumetric specific heat for snow                     [J/m3/K]
+   REAL(wp), PUBLIC ::   xlsn     =  110.121e+6_wp   !: volumetric latent heat fusion of snow                 [J/m3]
+   REAL(wp), PUBLIC ::   lfus                        !: = xlsn / rhosn   (latent heat of fusion of fresh ice) [J/Kg]
+   REAL(wp), PUBLIC ::   xlic     =  300.33e+6_wp    !: volumetric latent heat fusion of ice                  [J/m3]
+   REAL(wp), PUBLIC ::   xsn      =    2.8e+6_wp     !: volumetric latent heat of sublimation of snow         [J/m3]
+#endif
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -102,21 +105,11 @@
       !!----------------------------------------------------------------------
 
-      !                                   ! Define additional parameters
-      rsiyea = 365.25 * rday * 2. * rpi / 6.283076
-      rsiday = rday / ( 1. + rday / rsiyea )
-#if defined key_cice
-      omega =  7.292116e-05
-#else
-      omega  = 2. * rpi / rsiday 
-#endif
-
-      rau0r  = 1. /   rau0  
-      ro0cpr = 1. / ( rau0 * rcp )
-
-
-      IF(lwp) THEN                        ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) ' phy_cst : initialization of ocean parameters and constants'
-         WRITE(numout,*) ' ~~~~~~~'
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' phy_cst : initialization of ocean parameters and constants'
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~'
+
+      ! Ocean Parameters
+      ! ----------------
+      IF(lwp) THEN
          WRITE(numout,*) '       Domain info'
          WRITE(numout,*) '          dimension of model'
@@ -131,46 +124,76 @@
          WRITE(numout,*) '             jpnij   : ', jpnij
          WRITE(numout,*) '          lateral domain boundary condition type : jperio  = ', jperio
-         WRITE(numout,*)
-         WRITE(numout,*) '       Constants'
-         WRITE(numout,*)
-         WRITE(numout,*) '          mathematical constant                 rpi = ', rpi
-         WRITE(numout,*) '          day                                rday   = ', rday,   ' s'
-         WRITE(numout,*) '          sideral year                       rsiyea = ', rsiyea, ' s'
-         WRITE(numout,*) '          sideral day                        rsiday = ', rsiday, ' s'
-         WRITE(numout,*) '          omega                              omega  = ', omega,  ' s-1'
-         WRITE(numout,*)
-         WRITE(numout,*) '          nb of months per year               raamo = ', raamo, ' months'
-         WRITE(numout,*) '          nb of hours per day                 rjjhh = ', rjjhh, ' hours'
-         WRITE(numout,*) '          nb of minutes per hour              rhhmm = ', rhhmm, ' mn'
-         WRITE(numout,*) '          nb of seconds per minute            rmmss = ', rmmss, ' s'
-         WRITE(numout,*)
-         WRITE(numout,*) '          earth radius                         ra   = ', ra, ' m'
-         WRITE(numout,*) '          gravity                              grav = ', grav , ' m/s^2'
-         WRITE(numout,*)
-         WRITE(numout,*) '          triple point of temperature      rtt      = ', rtt     , ' K'
-         WRITE(numout,*) '          freezing point of water          rt0      = ', rt0     , ' K'
-         WRITE(numout,*) '          melting point of snow            rt0_snow = ', rt0_snow, ' K'
-         WRITE(numout,*) '          melting point of ice             rt0_ice  = ', rt0_ice , ' K'
-         WRITE(numout,*)
-         WRITE(numout,*) '          ocean reference volumic mass       rau0   = ', rau0 , ' kg/m^3'
-         WRITE(numout,*) '          ocean reference specific volume    rau0r  = ', rau0r, ' m^3/Kg'
-         WRITE(numout,*) '          ocean specific heat                rcp    = ', rcp
-         WRITE(numout,*) '                       1. / ( rau0 * rcp ) = ro0cpr = ', ro0cpr
+      ENDIF
+
+      ! Define constants
+      ! ----------------
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '       Constants'
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '          mathematical constant                 rpi = ', rpi
+
+      rsiyea = 365.25_wp * rday * 2._wp * rpi / 6.283076_wp
+      rsiday = rday / ( 1._wp + rday / rsiyea )
+#if defined key_cice
+      omega  = 7.292116e-05
+#else
+      omega  = 2._wp * rpi / rsiday 
+#endif
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '          day                                rday   = ', rday,   ' s'
+      IF(lwp) WRITE(numout,*) '          sideral year                       rsiyea = ', rsiyea, ' s'
+      IF(lwp) WRITE(numout,*) '          sideral day                        rsiday = ', rsiday, ' s'
+      IF(lwp) WRITE(numout,*) '          omega                              omega  = ', omega,  ' s^-1'
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '          nb of months per year               raamo = ', raamo, ' months'
+      IF(lwp) WRITE(numout,*) '          nb of hours per day                 rjjhh = ', rjjhh, ' hours'
+      IF(lwp) WRITE(numout,*) '          nb of minutes per hour              rhhmm = ', rhhmm, ' mn'
+      IF(lwp) WRITE(numout,*) '          nb of seconds per minute            rmmss = ', rmmss, ' s'
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '          earth radius                         ra   = ', ra, ' m'
+      IF(lwp) WRITE(numout,*) '          gravity                              grav = ', grav , ' m/s^2'
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '          triple point of temperature      rtt      = ', rtt     , ' K'
+      IF(lwp) WRITE(numout,*) '          freezing point of water          rt0      = ', rt0     , ' K'
+      IF(lwp) WRITE(numout,*) '          melting point of snow            rt0_snow = ', rt0_snow, ' K'
+      IF(lwp) WRITE(numout,*) '          melting point of ice             rt0_ice  = ', rt0_ice , ' K'
+
+      r1_rau0     = 1._wp / rau0
+      r1_rcp      = 1._wp / rcp
+      r1_rau0_rcp = 1._wp / ( rau0 * rcp )
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '          volumic mass of pure water          rauw  = ', rauw   , ' kg/m^3'
+      IF(lwp) WRITE(numout,*) '          volumic mass of reference           rau0  = ', rau0   , ' kg/m^3'
+      IF(lwp) WRITE(numout,*) '          1. / rau0                        r1_rau0  = ', r1_rau0, ' m^3/kg'
+      IF(lwp) WRITE(numout,*) '          ocean specific heat                 rcp   = ', rcp    , ' J/Kelvin'
+      IF(lwp) WRITE(numout,*) '          1. / ( rau0 * rcp )           r1_rau0_rcp = ', r1_rau0_rcp
+
+
+#if defined key_lim3 || defined key_cice
+      xlsn = lfus * rhosn        ! volumetric latent heat fusion of snow [J/m3]
+#else
+      cpic = rcpic / rhoic       ! specific heat for ice   [J/Kg/K]
+      lfus = xlsn / rhosn        ! latent heat of fusion of fresh ice
+#endif
+
+      IF(lwp) THEN
          WRITE(numout,*)
          WRITE(numout,*) '          thermal conductivity of the snow          = ', rcdsn   , ' J/s/m/K'
          WRITE(numout,*) '          thermal conductivity of the ice           = ', rcdic   , ' J/s/m/K'
-#if defined key_lim3
          WRITE(numout,*) '          fresh ice specific heat                   = ', cpic    , ' J/kg/K'
          WRITE(numout,*) '          latent heat of fusion of fresh ice / snow = ', lfus    , ' J/kg'
+#if defined key_lim3
          WRITE(numout,*) '          latent heat of subl.  of fresh ice / snow = ', lsub    , ' J/kg'
-#elif defined key_cice
-         WRITE(numout,*) '          latent heat of fusion of fresh ice / snow = ', lfus    , ' J/kg'
 #else
          WRITE(numout,*) '          density times specific heat for snow      = ', rcpsn   , ' J/m^3/K' 
          WRITE(numout,*) '          density times specific heat for ice       = ', rcpic   , ' J/m^3/K'
          WRITE(numout,*) '          volumetric latent heat fusion of sea ice  = ', xlic    , ' J/m' 
-         WRITE(numout,*) '          volumetric latent heat fusion of snow     = ', xlsn    , ' J/m' 
          WRITE(numout,*) '          latent heat of sublimation of snow        = ', xsn     , ' J/kg' 
 #endif
+         WRITE(numout,*) '          volumetric latent heat fusion of snow     = ', xlsn    , ' J/m^3' 
          WRITE(numout,*) '          density of sea ice                        = ', rhoic   , ' kg/m^3'
          WRITE(numout,*) '          density of snow                           = ', rhosn   , ' kg/m^3'
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_exp.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_exp.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_exp.F90	(revision 3396)
@@ -61,5 +61,5 @@
       !
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      REAL(wp) ::   zrau0r, zlavmr, zua, zva   ! local scalars
+      REAL(wp) ::   zlavmr, zua, zva   ! local scalars
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwx, zwy, zwz, zww
       !!----------------------------------------------------------------------
@@ -75,5 +75,4 @@
       ENDIF
 
-      zrau0r = 1. / rau0               ! Local constant initialization
       zlavmr = 1. / REAL( nn_zdfexp )
 
@@ -81,6 +80,6 @@
       DO jj = 2, jpjm1                 ! Surface boundary condition
          DO ji = 2, jpim1
-            zwy(ji,jj,1) = ( utau_b(ji,jj) + utau(ji,jj) ) * zrau0r
-            zww(ji,jj,1) = ( vtau_b(ji,jj) + vtau(ji,jj) ) * zrau0r
+            zwy(ji,jj,1) = ( utau_b(ji,jj) + utau(ji,jj) ) * r1_rau0
+            zww(ji,jj,1) = ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_rau0
          END DO  
       END DO  
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 3396)
@@ -161,5 +161,5 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             ua(ji,jj,1) = ub(ji,jj,1) + p2dt * (  ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
-               &                                                       / ( fse3u(ji,jj,1) * rau0       )  )
+               &                                                       * r1_rau0 / fse3u(ji,jj,1)       )
          END DO
       END DO
@@ -247,5 +247,5 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             va(ji,jj,1) = vb(ji,jj,1) + p2dt * (  va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
-               &                                                       / ( fse3v(ji,jj,1) * rau0       )  )
+               &                                                       * r1_rau0 / fse3v(ji,jj,1)       )
          END DO
       END DO
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90	(revision 3396)
@@ -60,5 +60,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
@@ -89,8 +89,7 @@
          nn_tau000 = MAX( nn_tau000, 1 )     ! must be >= 1
          !
-         qns (:,:) = rn_qns0
+         emp (:,:) = rn_emp0
+         qns (:,:) = rn_qns0 - emp(:,:) * sst_m(:,:) * rcp      ! including heat content associated with mass flux at SST
          qsr (:,:) = rn_qsr0
-         emp (:,:) = rn_emp0
-         emps(:,:) = rn_emp0
          !
          utau(:,:) = rn_utau0
@@ -130,5 +129,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.
@@ -211,5 +210,4 @@
          END DO
       END DO
-      emps(:,:) = emp(:,:)
 
       ! Compute the emp flux such as its integration on the whole domain at each time is zero
@@ -224,7 +222,7 @@
       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/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 3396)
@@ -12,10 +12,10 @@
 
    !!----------------------------------------------------------------------
-   !!   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 declination
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
@@ -50,6 +50,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)
@@ -100,14 +100,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
@@ -115,8 +114,10 @@
       !!              - 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
+      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
       !!
       INTEGER  ::   ifpr, jfpr   ! dummy indices
@@ -205,6 +206,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
       !!----------------------------------------------------------------------
@@ -223,5 +226,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), POINTER, DIMENSION(:,:) ::   zqlw        ! long-wave heat flux over ocean
       REAL(wp), POINTER, DIMENSION(:,:) ::   zqla        ! latent heat flux over ocean
@@ -363,10 +366,17 @@
       !     III    Total FLUXES                                                       !
       ! ----------------------------------------------------------------------------- !
-
-!CDIR COLLAPSE
-      emp (:,:) = zqla(:,:) / cevap - sf(jp_prec)%fnow(:,:,1) / rday * tmask(:,:,1)
-      qns (:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)         ! Downward Non Solar flux
-      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(:,:,1) / 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(:,:,1) * sf(jp_tair)%fnow(:,:,1) * 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
@@ -407,12 +417,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]
@@ -594,5 +605,5 @@
       !
       ! ----------------------------------------------------------------------------- !
-      !    Total FLUXES                                                       !
+      !    Total FLUXES                                                               !
       ! ----------------------------------------------------------------------------- !
       !
@@ -601,4 +612,13 @@
 !CDIR COLLAPSE
       p_tpr(:,:)   = sf(jp_prec)%fnow(:,:,1) / 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(:,:,1), rt0_snow - rt0 ) * cpic &   ! add solid P at least below melting
+         &     - p_spr(:,:) * sf(jp_tair)%fnow(:,:,1)                        * rcp      ! remove solid precip. at Tair
       !
 !!gm : not necessary as all input data are lbc_lnk...
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 3396)
@@ -52,5 +52,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)
@@ -69,5 +69,5 @@
    REAL(wp), PARAMETER ::   Stef =    5.67e-8     ! Stefan Boltzmann constant
    REAL(wp), PARAMETER ::   Cice =    1.63e-3     ! transfer coefficient over ice
-   REAL(wp), PARAMETER ::   albo =    0.066       ! ocean albedo assumed to be contant
+   REAL(wp), PARAMETER ::   albo =    0.066       ! ocean albedo assumed to be constant
 
    !                                  !!* Namelist namsbc_core : CORE bulk parameters
@@ -96,24 +96,21 @@
       !!      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
@@ -125,7 +122,6 @@
       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 ,           &
@@ -221,5 +217,5 @@
       !!              - 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)
+      !!              - emp     : evaporation minus precipitation       (kg/m2/s)
       !!
       !!  ** Nota  :   sf has to be a dummy argument for AGRIF on NEC
@@ -378,9 +374,14 @@
      
 !CDIR COLLAPSE
-      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)      ! Downward Non Solar flux
-!CDIR COLLAPSE
-      emp(:,:) = zevap(:,:) - sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)
-!CDIR COLLAPSE
-      emps(:,:) = emp(:,:)
+      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.)
+         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1)
+!CDIR COLLAPSE
+      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                       &   ! Downward Non Solar flux
+         &     - sf(jp_snow)%fnow(:,:,1) * lfus                          &   ! remove latent melting heat for solid precip
+         &     - zevap(:,:) * pst(:,:) * rcp                             &   ! remove evap heat content at SST
+         &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) )   &   ! add liquid precip heat content at Tair
+         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                 &   
+         &     + sf(jp_snow)%fnow(:,:,1)                                 &   ! add solid  precip heat content at min(Tair,Tsnow)
+         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic 
       !
       CALL iom_put( "qlw_oce",   zqlw )                 ! output downward longwave heat over the ocean
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 3396)
@@ -664,6 +664,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
@@ -821,15 +823,4 @@
          !                                                   ! ========================= !
          !
-         !                                                       ! non solar heat flux over the ocean (qns)
-         IF( srcv(jpr_qnsoce)%laction )   qns(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
-         IF( srcv(jpr_qnsmix)%laction )   qns(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
-         ! add the latent heat of solid precip. melting
-         IF( srcv(jpr_snow  )%laction )   qns(:,:) = qns(:,:) - frcv(jpr_snow)%z3(:,:,1) * lfus              
-
-         !                                                       ! solar flux over the ocean          (qsr)
-         IF( srcv(jpr_qsroce)%laction )   qsr(:,:) = frcv(jpr_qsroce)%z3(:,:,1)
-         IF( srcv(jpr_qsrmix)%laction )   qsr(:,:) = frcv(jpr_qsrmix)%z3(:,:,1)
-         IF( ln_dm2dc )   qsr(:,:) = sbc_dcy( qsr )                           ! modify qsr to include the diurnal cycle
-         !
          !                                                       ! total freshwater fluxes over the ocean (emp, emps)
          SELECT CASE( TRIM( sn_rcv_emp%cldes ) )                                    ! evaporation - precipitation
@@ -863,5 +854,18 @@
 !!gm  end of internal cooking
          !
-         emps(:,:) = emp(:,:)                                        ! concentration/dilution = emp
+         !                                                       ! non solar heat flux over the ocean (qns)
+         IF( srcv(jpr_qnsoce)%laction )   qns(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
+         IF( srcv(jpr_qnsmix)%laction )   qns(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
+         ! add the latent heat of solid precip. melting
+         IF( srcv(jpr_snow  )%laction )   THEN                         ! update qns over the free ocean with:
+              qns(:,:) = qns(:,:) - frcv(jpr_snow)%z3(:,:,1) * lfus  & ! energy for melting solid precipitation over the free ocean
+           &           - emp(:,:) * sst_m(:,:) * rcp                   ! remove heat content due to mass flux (assumed to be at SST)
+         ENDIF
+
+         !                                                       ! solar flux over the ocean          (qsr)
+         IF( srcv(jpr_qsroce)%laction )   qsr(:,:) = frcv(jpr_qsroce)%z3(:,:,1)
+         IF( srcv(jpr_qsrmix)%laction )   qsr(:,:) = frcv(jpr_qsrmix)%z3(:,:,1)
+         IF( ln_dm2dc )   qsr(:,:) = sbc_dcy( qsr )                           ! modify qsr to include the diurnal cycle
+         !
   
       ENDIF
@@ -1141,5 +1145,6 @@
 
       zicefr(:,:) = 1.- p_frld(:,:)
-      IF( lk_diaar5 )   zcptn(:,:) = rcp * tsn(:,:,1,jp_tem)
+      !zcptn(:,:) = rcp * tsn(:,:,1,jp_tem)
+      zcptn(:,:) = rcp * sst_m(:,:)
       !
       !                                                      ! ========================= !
@@ -1233,6 +1238,9 @@
             &                                                   +          pist(:,:,1)   * zicefr(:,:) ) )
       END SELECT
-      ztmp(:,:) = p_frld(:,:) * sprecip(:,:) * lfus               ! add the latent heat of solid precip. melting
-      qns_tot(:,:) = qns_tot(:,:) - ztmp(:,:)                     ! over free ocean 
+      ztmp(:,:) = p_frld(:,:) * sprecip(:,:) * lfus
+      qns_tot(:,:) = qns_tot(:,:)                         &            ! qns_tot update over free ocean with:
+         &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting
+         &          + (  emp_tot(:,:)                     &            ! remove the heat content of mass flux (assumed to be at SST)
+         &             - emp_ice(:,:) * p_frld(:,:,1)  ) * zcptn(:,:) 
       IF( lk_diaar5 )   CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
 !!gm
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 3396)
@@ -61,6 +61,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
@@ -68,6 +67,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
@@ -139,4 +139,7 @@
             END DO
          END DO
+         !                                                        ! add to qns the heat due to e-p
+         qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
+         !
          !                                                        ! module of wind stress and wind speed at T-point
          zcoef = 1. / ( zrhoa * zcdrag )
@@ -154,6 +157,4 @@
          CALL lbc_lnk( taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( wndm(:,:), 'T', 1. )
 
-         emps(:,:) = emp (:,:)                                    ! Initialization of emps (needed when no ice model)
-                  
          IF( nitend-nit000 <= 100 .AND. lwp ) THEN                ! control print (if less than 100 time-step asked)
             WRITE(numout,*) 
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 3396)
@@ -64,8 +64,9 @@
       INTEGER, INTENT( in ) ::   kn_fwb   ! ocean time-step index
       !
-      INTEGER  ::   inum, ikty, iyear   ! local integers
-      REAL(wp) ::   z_fwf, z_fwf_nsrf, zsum_fwf, zsum_erp   ! local scalars
-      REAL(wp) ::   zsurf_neg, zsurf_pos, zsurf_tospread    !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:) ::   ztmsk_neg, ztmsk_pos, ztmsk_tospread, z_wgt, zerp_cor
+      INTEGER  ::   inum, ikty, iyear     ! local integers
+      REAL(wp) ::   z_fwf, z_fwf_nsrf, zsum_fwf, zsum_erp                ! local scalars
+      REAL(wp) ::   zsurf_neg, zsurf_pos, zsurf_tospread, zcoef          !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztmsk_neg, ztmsk_pos, z_wgt ! 2D workspaces
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztmsk_tospread, zerp_cor    !   -      -
       !!----------------------------------------------------------------------
       !
@@ -96,6 +97,7 @@
          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
             z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) ) ) / area   ! sum over the global domain
-            emp (:,:) = emp (:,:) - z_fwf 
-            emps(:,:) = emps(:,:) - z_fwf 
+            zcoef = z_fwf * rcp
+            emp(:,:) = emp(:,:) - z_fwf 
+            qns(:,:) = qns(:,:) + zcoef * sst_m(:,:)                          ! ensure fw correction does not change the heat budget
          ENDIF
          !
@@ -103,5 +105,5 @@
          !
          IF( kt == nit000 ) THEN                   ! initialisation
-            !                                         ! Read the corrective factor on precipitations (fwfold)
+            !                                      ! Read the corrective factor on precipitations (fwfold)
             CALL ctl_opn( inum, 'EMPave_old.dat', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
             READ ( inum, "(24X,I8,2ES24.16)" ) iyear, a_fwb_b, a_fwb
@@ -125,10 +127,11 @@
          ENDIF
          ! 
-         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN      ! correct the freshwater fluxes
-            emp (:,:) = emp (:,:) + fwfold
-            emps(:,:) = emps(:,:) + fwfold
-         ENDIF
-         !
-         IF( kt == nitend .AND. lwp ) THEN         ! save fwfold value in a file
+         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN         ! correct the freshwater fluxes
+            zcoef = fwfold * rcp
+            emp(:,:) = emp(:,:) + fwfold
+            qns(:,:) = qns(:,:) - zcoef * sst_m(:,:)  ! ensure fw correction does not change the heat budget
+         ENDIF
+         !
+         IF( kt == nitend .AND. lwp ) THEN            ! save fwfold value in a file
             CALL ctl_opn( inum, 'EMPave.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
             WRITE( inum, "(24X,I8,2ES24.16)" ) nyear, a_fwb_b, a_fwb
@@ -143,8 +146,8 @@
             ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:)
             !
-            zsurf_neg = glob_sum( e1e2t(:,:)*ztmsk_neg(:,:) )   ! Area filled by <0 and >0 erp 
+            zsurf_neg = glob_sum( e1e2t(:,:)*ztmsk_neg(:,:) )  ! Area filled by <0 and >0 erp 
             zsurf_pos = glob_sum( e1e2t(:,:)*ztmsk_pos(:,:) )
             !                                                  ! fwf global mean 
-            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) ) ) / area
+            z_fwf     = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) ) ) / area
             !            
             IF( z_fwf < 0._wp ) THEN         ! spread out over >0 erp area to increase evaporation
@@ -160,5 +163,5 @@
             z_fwf_nsrf =  zsum_fwf / ( zsurf_tospread + rsmall )
             !                                                  ! weight to respect erp field 2D structure 
-            zsum_erp = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e1e2t(:,:) )
+            zsum_erp   = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e1e2t(:,:) )
             z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( zsum_erp + rsmall )
             !                                                  ! final correction term to apply
@@ -168,7 +171,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                   ! control print
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 3396)
@@ -140,4 +140,7 @@
       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   
       IF( MOD( nitend - nit000 + 1, nn_fsbc) /= 0 .OR.   &
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 3396)
@@ -56,5 +56,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rnf_tsc_b, rnf_tsc  !: before and now T & S runoff contents   [K.m/s & PSU.m/s]
    
-   REAL(wp) ::   r1_rau0   ! = 1 / rau0 
 
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnf       ! structure: river runoff (file information, fields read)
@@ -129,5 +128,4 @@
             rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) )  
             !
-            r1_rau0 = 1._wp / rau0
             !                                                     ! set temperature & salinity content of runoffs
             IF( ln_rnf_tem ) THEN                                       ! use runoffs temperature data
@@ -199,5 +197,4 @@
       !!
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   r1_rau0   ! local scalar
       REAL(wp) ::   zfact     ! local scalar
       !!----------------------------------------------------------------------
@@ -205,5 +202,4 @@
       zfact = 0.5_wp
       !
-      r1_rau0 = 1._wp / rau0
       IF( ln_rnf_depth ) THEN      !==   runoff distributed over several levels   ==!
          IF( lk_vvl ) THEN             ! variable volume case 
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 3396)
@@ -156,5 +156,5 @@
             !                                      ! ========================= !
             !
-            IF( nn_sstr == 1 ) THEN                   !* Temperature restoring term
+            IF( nn_sstr == 1 ) THEN                                   !* Temperature restoring term
 !CDIR COLLAPSE
                DO jj = 1, jpj
@@ -168,5 +168,5 @@
             ENDIF
             !
-            IF( nn_sssr == 1 ) THEN                   !* Salinity damping term (salt flux, emps only)
+            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
@@ -174,13 +174,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,1) )   &
-                        &        / ( sss_m(ji,jj) + 1.e-20   )
-                     emps(ji,jj) = emps(ji,jj) + zerp
-                     erp( ji,jj) = zerp
+                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) 
+                     emps(ji,jj) = emps(ji,jj) + zerp                 ! salt flux
+                     erp( ji,jj) = zerp / MAX( sss_m(ji,jj), 1.e-20 ) ! converted into an equivalent volume flux (diagnostic 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 associated heat flux (qns)
                zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s]
                zerp_bnd = rn_sssr_bnd / rday                          !       -              -    
@@ -190,9 +189,9 @@
                      zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
                         &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )   &
-                        &        / ( 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
+                     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/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 3396)
@@ -121,5 +121,4 @@
       REAL(wp) ::   zd , zc , zaw, za    !   -      -
       REAL(wp) ::   zb1, za1, zkw, zk0   !   -      -
-      REAL(wp) ::   zrau0r               !   -      -
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
       !!----------------------------------------------------------------------
@@ -133,5 +132,4 @@
       !
       CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==!
-         zrau0r = 1.e0 / rau0
 !CDIR NOVERRCHK
          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
@@ -174,5 +172,5 @@
                   ! masked in situ density anomaly
                   prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    &
-                     &             - rau0  ) * zrau0r * tmask(ji,jj,jk)
+                     &             - rau0  ) * r1_rau0 * tmask(ji,jj,jk)
                END DO
             END DO
@@ -254,5 +252,5 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw   ! local scalars
-      REAL(wp) ::   zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zrau0r       !   -      -
+      REAL(wp) ::   zb, zd, zc, zaw, za, zb1, za1, zkw, zk0               !   -      -
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
       !!----------------------------------------------------------------------
@@ -265,5 +263,4 @@
       !
       CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==!
-         zrau0r = 1.e0 / rau0
 !CDIR NOVERRCHK
          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
@@ -309,5 +306,5 @@
                   ! masked in situ density anomaly
                   prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    &
-                     &             - rau0  ) * zrau0r * tmask(ji,jj,jk)
+                     &             - rau0  ) * r1_rau0 * tmask(ji,jj,jk)
                END DO
             END DO
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 3396)
@@ -155,5 +155,5 @@
          CASE ( 1 )                          !* constant flux
             IF(lwp) WRITE(numout,*) '      *** constant heat flux  =   ', rn_geoflx_cst
-            qgh_trd0(:,:) = ro0cpr * rn_geoflx_cst
+            qgh_trd0(:,:) = r1_rau0_rcp * rn_geoflx_cst
             !
          CASE ( 2 )                          !* variable geothermal heat flux : read the geothermal fluxes in mW/m2
@@ -162,5 +162,5 @@
             CALL iom_get  ( inum, jpdom_data, 'heatflow', qgh_trd0 )
             CALL iom_close( inum )
-            qgh_trd0(:,:) = ro0cpr * qgh_trd0(:,:) * 1.e-3     ! conversion in W/m2
+            qgh_trd0(:,:) = r1_rau0_rcp * qgh_trd0(:,:) * 1.e-3     ! conversion in W/m2
             !
          CASE DEFAULT
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 3396)
@@ -147,5 +147,5 @@
          !                                        ! ============================================== !
          DO jk = 1, jpkm1
-            qsr_hc(:,:,jk) = ro0cpr * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
+            qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
          END DO
          !                                        Add to the general trend
@@ -219,5 +219,5 @@
                !
                DO jk = 1, nksr                                        ! compute and add qsr trend to ta
-                  qsr_hc(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) )
+                  qsr_hc(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) )
                END DO
                zea(:,:,nksr+1:jpk) = 0.e0     ! below 400m set to zero
@@ -236,6 +236,6 @@
             !
             IF( lk_vvl ) THEN                                  !* variable volume
-               zz0   =        rn_abs   * ro0cpr
-               zz1   = ( 1. - rn_abs ) * ro0cpr
+               zz0   =        rn_abs   * r1_rau0_rcp
+               zz1   = ( 1. - rn_abs ) * r1_rau0_rcp
                DO jk = 1, nksr                    ! solar heat absorbed at T-point in the top 400m 
                   DO jj = 1, jpj
@@ -463,5 +463,5 @@
                   !
                   DO jk = 1, nksr
-                     etot3(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) 
+                     etot3(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) ) 
                   END DO
                   etot3(:,:,nksr+1:jpk) = 0.e0                ! below 400m set to zero
@@ -484,6 +484,6 @@
                IF(lwp) WRITE(numout,*) '        key_vvl: light distribution will be computed at each time step'
             ELSE                                ! constant volume: computes one for all
-               zz0 =        rn_abs   * ro0cpr
-               zz1 = ( 1. - rn_abs ) * ro0cpr
+               zz0 =        rn_abs   * r1_rau0_rcp
+               zz1 = ( 1. - rn_abs ) * r1_rau0_rcp
                DO jk = 1, nksr                    !*  solar heat absorbed at T-point computed once for all
                   DO jj = 1, jpj                              ! top 400 meters
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 3396)
@@ -60,4 +60,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: 
@@ -84,6 +86,8 @@
       !!            (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
@@ -109,5 +113,5 @@
       !!
       INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices  
-      REAL(wp) ::   zfact, z1_e3t, zsrau, zdep
+      REAL(wp) ::   zfact, z1_e3t, zdep
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
       !!----------------------------------------------------------------------
@@ -120,6 +124,4 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
       ENDIF
-
-      zsrau = 1. / rau0             ! initialization
 
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
@@ -163,11 +165,9 @@
                                                    ! evaporation, precipitation and qns, but not river runoff 
                                                
-      IF( lk_vvl ) THEN                            ! Variable Volume case
+      IF( lk_vvl ) THEN                            ! Variable Volume case  ==>> heat content of mass flux is in qns
          DO jj = 1, jpj
             DO ji = 1, jpi 
-               ! temperature : heat flux + cooling/heating effet of EMP flux
-               sbc_tsc(ji,jj,jp_tem) = ro0cpr * qns(ji,jj) - zsrau * emp(ji,jj) * tsn(ji,jj,1,jp_tem)
-               ! concent./dilut. effect due to sea-ice melt/formation and (possibly) SSS restoration
-               sbc_tsc(ji,jj,jp_sal) = ( emps(ji,jj) - emp(ji,jj) ) * zsrau * tsn(ji,jj,1,jp_sal)
+               sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)                              ! non solar heat flux
+               sbc_tsc(ji,jj,jp_sal) = r1_rau0     * emps(ji,jj)                             ! salt flux (freezing/melting)
             END DO
          END DO
@@ -176,7 +176,9 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                ! temperature : heat flux
-               sbc_tsc(ji,jj,jp_tem) = ro0cpr * qns(ji,jj)
+               sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)                          &   ! non solar heat flux
+                  &                  + r1_rau0     * emp(ji,jj)  * tsn(ji,jj,1,jp_tem)       ! concent./dilut. effect
                ! salinity    : salt flux + concent./dilut. effect (both in emps)
-               sbc_tsc(ji,jj,jp_sal) = zsrau * emps(ji,jj) * tsn(ji,jj,1,jp_sal)
+               sbc_tsc(ji,jj,jp_sal) = r1_rau0  * (  emps(ji,jj)                         &   ! salt flux (freezing/melting)
+                  &                                + emp (ji,jj) * tsn(ji,jj,1,jp_sal) )     ! concent./dilut. effect
             END DO
          END DO
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 3396)
@@ -167,5 +167,5 @@
             ! 
             ! surface friction 
-            ustars2(ji,jj) = rau0r * taum(ji,jj) * tmask(ji,jj,1)
+            ustars2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1)
             !
             ! bottom friction (explicit before friction)
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90	(revision 3386)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90	(revision 3396)
@@ -428,5 +428,5 @@
             zBo   (ji,jj) = grav * zthermal * qns(ji,jj) -  grav * zhalin * ( emps(ji,jj)-rnf(ji,jj) ) 
             ! Surface Temperature flux for non-local term
-            wt0(ji,jj) = - ( qsr(ji,jj) + qns(ji,jj) )* ro0cpr * tmask(ji,jj,1)
+            wt0(ji,jj) = - ( qsr(ji,jj) + qns(ji,jj) )* r1_rau0_rcp * tmask(ji,jj,1)
             ! Surface salinity flux for non-local term
             ws0(ji,jj) = - ( ( emps(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal) * rcs ) * tmask(ji,jj,1) 
