Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/ice_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/ice_2.F90	(revision 1857)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/ice_2.F90	(revision 1858)
@@ -4,5 +4,6 @@
    !! Sea Ice physics:  diagnostics variables of ice defined in memory
    !!=====================================================================
-   !! History :  2.0  !  03-08  (C. Ethe)  F90: Free form and module
+   !! History :  LIM  ! 2003-08  (C. Ethe)  F90: Free form and module
+   !!            2.0  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -68,6 +69,8 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qstoif        !: Energy stored in the brine pockets
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fbif          !: Heat flux at the ice base
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdmsnif       !: Variation of snow mass
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdmicif       !: Variation of ice mass
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdm_snw       !: Variation of snow mass over 1 time step           [Kg/m2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdq_snw       !: heat content associated to rdm_snw                [J/m2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdq_ice       !: Variation of ice  mass over 1 time step           [Kg/m2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdq_ice       !: heat content associated to rdm_snw                [J/m2]
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qldif         !: heat balance of the lead (or of the open ocean)
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qcmif         !: Energy needed to bring the ocean surface layer until its freezing 
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 1857)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 1858)
@@ -4,12 +4,12 @@
    !!           computation of the flux at the sea ice/ocean interface
    !!======================================================================
-   !! History : 00-01 (H. Goosse) Original code
-   !!           02-07 (C. Ethe, G. Madec) re-writing F90
-   !!           06-07 (G. Madec) surface module
+   !! History : LIM  ! 2000-01 (H. Goosse) Original code
+   !!           2.0  ! 2002-07 (C. Ethe, G. Madec) re-writing F90
+   !!            -   ! 2006-07 (G. Madec) surface module
+   !!           2.1  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges
    !!----------------------------------------------------------------------
 #if defined key_lim2
    !!----------------------------------------------------------------------
    !!   'key_lim2'                                    LIM 2.0 sea-ice model
-   !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
    !!   lim_sbc_2  : flux at the ice / ocean interface
@@ -211,5 +211,5 @@
 #if defined key_coupled
           zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   ! 
-             &   + rdmsnif(ji,jj) * zrdtir                                      !  freshwaterflux due to snow melting 
+             &   + rdm_snw(ji,jj) * zrdtir                                      !  freshwaterflux due to snow melting 
 #else
 !!$            !  computing freshwater exchanges at the ice/ocean interface
@@ -217,15 +217,15 @@
 !!$               &   + tprecip(ji,jj)                           &   !  total precipitation
 !!$               &   - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )   &   !  remov. snow precip over ice
-!!$               &   - rdmsnif(ji,jj) / rdt_ice                     !  freshwaterflux due to snow melting 
+!!$               &   - rdm_snw(ji,jj) / rdt_ice                     !  freshwaterflux due to snow melting 
             !  computing freshwater exchanges at the ice/ocean interface
             zemp = + emp(ji,jj)     *         frld(ji,jj)      &   !  e-p budget over open ocean fraction 
                &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &   !  liquid precipitation reaches directly the ocean
                &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  taking into account change in ice cover within the time step
-               &   + rdmsnif(ji,jj) * zrdtir                       !  freshwaterflux due to snow melting 
+               &   + rdm_snw(ji,jj) * zrdtir                       !  freshwaterflux due to snow melting 
                !                                                   !  ice-covered fraction:
 #endif            
 
             !  computing salt exchanges at the ice/ocean interface
-            zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdmicif(ji,jj) * zrdtir ) 
+            zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdm_ice(ji,jj) * zrdtir ) 
             
             !  converting the salt flux from ice to a freshwater flux from ocean
@@ -239,7 +239,7 @@
 
       IF( lk_diaar5 ) THEN
-         CALL iom_put( 'isnwmlt_cea'  ,                 rdmsnif(:,:) * zrdtir )
-         CALL iom_put( 'fsal_virt_cea',   soce_r(:,:) * rdmicif(:,:) * zrdtir )
-         CALL iom_put( 'fsal_real_cea', - sice_r(:,:) * rdmicif(:,:) * zrdtir )
+         CALL iom_put( 'isnwmlt_cea'  ,                 rdm_snw(:,:) * zrdtir )
+         CALL iom_put( 'fsal_virt_cea',   soce_r(:,:) * rdm_ice(:,:) * zrdtir )
+         CALL iom_put( 'fsal_real_cea', - sice_r(:,:) * rdm_ice(:,:) * zrdtir )
       ENDIF
 
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_2.F90	(revision 1857)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_2.F90	(revision 1858)
@@ -4,8 +4,9 @@
    !!              LIM thermo ice model : ice thermodynamic
    !!======================================================================
-   !! History :  1.0  ! 2000-01 (LIM)
-   !!            2.0  ! 2002-07 (C. Ethe, G. Madec) F90
-   !!            2.0  ! 2003-08 (C. Ethe)  add lim_thd_init
-   !!             -   ! 2008-2008  (A. Caubel, G. Madec, E. Maisonnave, S. Masson ) generic coupled interface
+   !! History :  1.0  !  2000-01  (LIM) Original code
+   !!            2.0  !  2002-07  (C. Ethe, G. Madec) F90
+   !!            2.0  !  2003-08  (C. Ethe)  add lim_thd_init
+   !!             -   !  2008-08  (A. Caubel, G. Madec, E. Maisonnave, S. Masson ) generic coupled interface
+   !!            2.1  !  2010-05  (Y. Aksenov G. Madec) add heat associated with mass exchanges
    !!---------------------------------------------------------------------
 #if defined key_lim2
@@ -88,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
       REAL(wp), DIMENSION(jpi,jpj)     ::   ztmp      ! 2D workspace
@@ -132,6 +131,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
 
@@ -193,10 +194,7 @@
       !      Partial computation of forcing for the thermodynamic sea ice model.
       !--------------------------------------------------------------------------
-
-      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)
@@ -211,25 +209,19 @@
             !  net downward heat flux from the ice to the ocean, expressed as a function of ocean 
             !  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) ) 
+            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) + 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) 
@@ -241,6 +233,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.
@@ -251,7 +242,5 @@
          END DO
       END DO
-      
-      sst_m(:,:) = sst_m(:,:) - rt0
-               
+                     
       !         Select icy points and fulfill arrays for the vectorial grid.
       !----------------------------------------------------------------------
@@ -307,5 +296,11 @@
          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) )
+!!gm bug ??  rdm_snw, rdq_snw should be transformed into 1D fields ?  like rdm_ice and rdq_ice ???
+!!gm         I think that  in fact the transformation is useless for both snow and ice as they are already set to zero
+!!gm         at the begining of the routine and only the value over icy point is updated in 1D and converted into 2D
+!!gm         the two lines below are thus useless.
+         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) )
+!!gm bug end.
          CALL tab_2d_1d_2( nbpb, dmgwi_1d   (1:nbpb)     , dmgwi      , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d_2( nbpb, qlbbq_1d   (1:nbpb)     , zqlbsbq    , jpi, jpj, npb(1:nbpb) )
@@ -327,7 +322,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, 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, 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, zdvosif    , npb, dvsbq_1d  (1:nbpb)     , jpi, jpj )
          CALL tab_1d_2d_2( nbpb, zdvobif    , npb, dvbbq_1d  (1:nbpb)     , jpi, jpj )
@@ -387,5 +384,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) )
@@ -398,5 +395,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) )
@@ -412,5 +410,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 )
          !
@@ -443,10 +442,10 @@
       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(:,:) )
+         zlicegr(:,:) = MAX( 0.e0, rdm_ice(:,:)-zlicegr(:,:) )
          CALL iom_put( 'licepro_cea' , zlicegr(:,:) * zztmp     )   ! Latereal sea ice growth    [kg/m2/s]
       ENDIF
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_lac_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_lac_2.F90	(revision 1857)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_lac_2.F90	(revision 1858)
@@ -1,12 +1,16 @@
 MODULE limthd_lac_2
-#if defined key_lim2
    !!======================================================================
    !!                       ***  MODULE limthd_lac_2   ***
    !!                lateral thermodynamic growth of the ice 
    !!======================================================================
-
+   !! History : LIM  ! 2000-04  (UCL) Original code
+   !!           2.0  ! 2002-07  (C. Ethe, G. Madec) F90 & mpp
+   !!            -   ! 2006-07  (G. Madec) surface module
+   !!           2.1  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges
+   !!----------------------------------------------------------------------
+#if defined key_lim2
    !!----------------------------------------------------------------------
    !!   lim_lat_acr_2    : lateral accretion of ice
-   !! * Modules used
+   !!----------------------------------------------------------------------
    USE par_oce          ! ocean parameters
    USE phycst
@@ -60,11 +64,6 @@
       !!               update h_snow_1d, h_ice_1d and tbif_1d(:,:)      
       !! 
-      !! ** References :
-      !!      M. Maqueda, 1995, PhD Thesis, Univesidad Complutense Madrid
-      !!      Fichefet T. and M. Maqueda 1997, J. Geo. Res., 102(C6), 
-      !!                                                12609 -12646   
-      !! History :
-      !!   1.0  !  01-04 (LIM)  original code
-      !!   2.0  !  02-08 (C. Ethe, G. Madec)  F90, mpp
+      !! References : M. Maqueda, 1995, PhD Thesis, Univesidad Complutense Madrid
+      !!              Fichefet T. and M. Maqueda 1997, JGR, 102(C6), 12609 -12646   
       !!-------------------------------------------------------------------
       !! * Arguments
@@ -144,6 +143,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
  
@@ -195,5 +194,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  & 
@@ -213,8 +212,8 @@
       !           dV = Vnew - Vold
       !-------------------------------------------------------------
-      
       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/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_zdf_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_zdf_2.F90	(revision 1857)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_zdf_2.F90	(revision 1858)
@@ -4,6 +4,7 @@
    !!                thermodynamic growth and decay of the ice 
    !!======================================================================
-   !! History :  1.0  !  01-04 (LIM) Original code
-   !!            2.0  !  02-08 (C. Ethe, G. Madec) F90
+   !! History : LIM  ! 2001-04 (UCL) Original code
+   !!           2.0  ! 2002-08 (C. Ethe, G. Madec) re-writing F90
+   !!           2.1  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -158,4 +159,5 @@
           , zdrfrl1, zdrfrl2 &  ! tempory scalars
           , zihsn, zidhb, zihic, zihe, zihq, ziexp, ziqf, zihnf, zibmlt, ziqr, zihgnew, zind
+       REAL(wp) ::   ztmp   ! temporary scalar
        !!----------------------------------------------------------------------
 
@@ -168,16 +170,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
-          zqcmlt(ji,1) = ( MAX ( zzero ,  &
-             &                   rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn )
-          !--computation of energy due to bottom melting
-          zqcmlt(ji,2) = ( 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. ) ) &
+          zihic = MAX( zzero , SIGN( zone , hicdif - h_ice_1d (ji) ) )
+          !--energy required to bring snow to its melting point (rt0_snow)
+          zqcmlt(ji,1) = ( 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)
+          zqcmlt(ji,2) = ( 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
+          !--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) )
@@ -478,5 +478,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
@@ -621,9 +624,16 @@
           !---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 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) )
@@ -657,5 +667,8 @@
           tbif_1d(ji,3) =  zihgnew * ztb3 + ( 1.0 - zihgnew ) * tfu_1d(ji)
           h_ice_1d(ji)  =  zhicnew
-       END DO
+          ! update the ice heat content given to the ocean in freezing case (part from rt0_ice to tfu_1d)
+         ztmp = ( 1. - zidhb ) * rhoic * dvbbq_1d(ji)
+         rdqicif_1d(ji) = rdqicif_1d(ji) + cpic * ztmp * ( tfu_1d(ji) - rt0_ice )
+      END DO
 
 
@@ -698,6 +711,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.
@@ -746,8 +764,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 at 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 at 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 rdqicif_1d
+!!gm                 cpic * ztmp * ( rt0_snow - rt0 )
+!!gm   Currently we put all the heat in rdqicif_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/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/thd_ice_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/thd_ice_2.F90	(revision 1857)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/thd_ice_2.F90	(revision 1858)
@@ -5,6 +5,6 @@
    !! LIM 2.0 sea-ice :   Ice thermodynamics in 1D
    !!=====================================================================
-   !! History :
-   !!   2.0  !  02-11  (C. Ethe)  F90: Free form and module
+   !! History : 2.0  ! 2002-11  (C. Ethe)  F90: Free form and module
+   !!           2.1  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges
    !!----------------------------------------------------------------------
    !!   LIM 2.0, UCL-LOCEAN-IPSL (2005)
@@ -67,6 +67,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
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/DOM/phycst.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/DOM/phycst.F90	(revision 1857)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/DOM/phycst.F90	(revision 1858)
@@ -4,9 +4,9 @@
    !!     Definition of of both ocean and ice parameters used in the code
    !!=====================================================================
-   !! History :        !  90-10  (C. Levy - G. Madec)  Original code
-   !!                  !  91-11  (G. Madec)
-   !!                  !  91-12  (M. Imbard)
-   !!             8.5  !  02-08  (G. Madec, C. Ethe)  F90, add ice constants
-   !!             9.0  !  06-08  (G. Madec)  style 
+   !! History :  OPA  !  1990-10  (C. Levy - G. Madec)  Original code
+   !!                 !  1991-11  (G. Madec, M. Imbard)
+   !!  NEMO      1.0  !  2002-08  (G. Madec, C. Ethe)  F90, add ice constants
+   !!             -   !  2006-08  (G. Madec)  style 
+   !!            3.0  !  2008-08  (G. Madec)  add LIM-3  
    !!----------------------------------------------------------------------
 
@@ -26,66 +26,66 @@
    REAL(wp), PUBLIC ::   rsmall = 0.5 * EPSILON( 1. )           !: smallest real computer value
    
-   REAL(wp), PUBLIC ::          & !:
-      rday = 24.*60.*60.  ,     & !: day (s)
-      rsiyea              ,     & !: sideral year (s)
-      rsiday              ,     & !: sideral day (s)
-      raamo =  12._wp     ,     & !: number of months in one year
-      rjjhh =  24._wp     ,     & !: number of hours in one day
-      rhhmm =  60._wp     ,     & !: number of minutes in one hour
-      rmmss =  60._wp     ,     & !: number of seconds in one minute
-!!!   omega = 7.292115083046061e-5_wp ,  &  !: change the last digit!
-      omega               ,    &  !: earth rotation parameter
-      ra    = 6371229._wp ,    &  !: earth radius (meter)
-      grav  = 9.80665_wp          !: gravity (m/s2)
+   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)
-      rt0      = 273.15_wp  ,  &  !: freezing point of water (Kelvin)
-#if defined key_lim3
-      rt0_snow = 273.16_wp  ,  &  !: melting point of snow  (Kelvin)
-      rt0_ice  = 273.16_wp  ,  &  !: melting point of ice   (Kelvin)
-#else
-      rt0_snow = 273.15_wp  ,  &  !: melting point of snow  (Kelvin)
-      rt0_ice  = 273.05_wp  ,  &  !: melting point of ice   (Kelvin)
-#endif
-      rau0     = 1035._wp   ,  &  !: volumic mass of reference (kg/m3)
-      rauw     = 1000._wp   ,  &  !: volumic mass of pure water (kg/m3)
-      rcp      =    4.e+3_wp,  &  !: ocean specific heat
-      ro0cpr                      !: = 1. / ( rau0 * rcp )
-
-   REAL(wp), PUBLIC ::            &  !:
-#if defined key_lim3
-      rcdsn   =   0.31_wp     ,   &  !: thermal conductivity of snow
-      rcdic   =   2.034396_wp ,   &  !: thermal conductivity of fresh ice
-      cpic    = 2067.0        ,   &
-      ! add the following lines
-      lsub    = 2.834e+6      ,   &  !: pure ice latent heat of sublimation (J.kg-1)
-      lfus    = 0.334e+6      ,   &  !: latent heat of fusion of fresh ice   (J.kg-1)
-      rhoic   = 917._wp       ,   &  !: volumic mass of sea ice (kg/m3)
-      tmut    =   0.054       ,   &  !: decrease of seawater meltpoint with salinity
-#else
-      rcdsn   =   0.22_wp     ,   &  !: conductivity of the snow
-      rcdic   =   2.034396_wp ,   &  !: conductivity of the ice
-      rcpsn   =   6.9069e+5_wp,   &  !: density times specific heat for snow
-      rcpic   =   1.8837e+6_wp,   &  !: volumetric latent heat fusion of sea ice
-      xlsn    = 110.121e+6_wp ,   &  !: volumetric latent heat fusion of snow
-      xlic    = 300.33e+6_wp  ,   &  !: volumetric latent heat fusion of ice
-      xsn     =   2.8e+6      ,   &  !: latent heat of sublimation of snow
-      rhoic   = 900._wp       ,   &  !: volumic mass of sea ice (kg/m3)
-#endif
-      rhosn   = 330._wp       ,   &  !: volumic mass of snow (kg/m3)
-      emic    =   0.97_wp     ,   &  !: emissivity of snow or ice
-      sice    =   6.0_wp      ,   &  !: salinity of ice (psu)
-      soce    =  34.7_wp      ,   &  !: salinity of sea (psu)
-      cevap   =   2.5e+6_wp   ,   &  !: latent heat of evaporation (water)
-      srgamma =   0.9_wp      ,   &  !: correction factor for solar radiation (Oberhuber, 1974)
-      vkarmn  =   0.4_wp      ,   &  !: von Karman constant
-      stefan  =   5.67e-8_wp         !: Stefan-Boltzmann constant 
-      !!----------------------------------------------------------------------
-      !!  OPA 9.0 , LOCEAN-IPSL (2005) 
-      !! $Id$ 
-      !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
-      !!----------------------------------------------------------------------
-   
+   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 ::   rau0     = 1035._wp         !: volumic mass of reference     [kg/m3]
+   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)
+
+#if defined key_lim3
+   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
+   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 
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
+   !! $Id$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
    
@@ -153,10 +153,21 @@
       IF(lwp) WRITE(numout,*) '          melting point of ice             rt0_ice  = ', rt0_ice , ' K'
 
-      ro0cpr = 1. / ( 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,*) '          ocean specific heat                rcp    = ', rcp
-      IF(lwp) WRITE(numout,*) '                       1. / ( rau0 * rcp ) = ro0cpr = ', ro0cpr
+      r1_rau0     = 1.e0 / rau0
+      r1_rcp      = 1.e0 / rcp
+      r1_rau0_rcp = 1.e0 / ( rau0 * rcp )
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '          volumic mass of pure water          rauw  = ', rauw   , ' kg/m3'
+      IF(lwp) WRITE(numout,*) '          volumic mass of reference           rau0  = ', rau0   , ' kg/m3'
+      IF(lwp) WRITE(numout,*) '          1. / rau0                        r1_rau0  = ', r1_rau0, ' m3/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
+      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
@@ -164,21 +175,21 @@
          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'
 #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,*) '          density times specific heat for snow      = ', rcpsn   , ' J/m3/K' 
+         WRITE(numout,*) '          density times specific heat for ice       = ', rcpic   , ' J/m3/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,*) '          density of sea ice                        = ', rhoic   , ' kg/m^3'
-         WRITE(numout,*) '          density of snow                           = ', rhosn   , ' kg/m^3'
+         WRITE(numout,*) '          volumetric latent heat fusion of snow     = ', xlsn    , ' J/m^3' 
+         WRITE(numout,*) '          density of sea ice                        = ', rhoic   , ' kg/m3'
+         WRITE(numout,*) '          density of snow                           = ', rhosn   , ' kg/m3'
          WRITE(numout,*) '          emissivity of snow or ice                 = ', emic  
          WRITE(numout,*) '          salinity of ice                           = ', sice    , ' psu'
          WRITE(numout,*) '          salinity of sea                           = ', soce    , ' psu'
-         WRITE(numout,*) '          latent heat of evaporation (water)        = ', cevap   , ' J/m^3' 
+         WRITE(numout,*) '          latent heat of evaporation (water)        = ', cevap   , ' J/m3' 
          WRITE(numout,*) '          correction factor for solar radiation     = ', srgamma 
          WRITE(numout,*) '          von Karman constant                       = ', vkarmn 
@@ -191,5 +202,5 @@
          WRITE(numout,*) '          smallest real computer value       rsmall = ', rsmall
       ENDIF
-
+      !
    END SUBROUTINE phy_cst
 
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/DYN/dynzdf_exp.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/DYN/dynzdf_exp.F90	(revision 1857)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/DYN/dynzdf_exp.F90	(revision 1858)
@@ -54,12 +54,10 @@
       !! ** Action : - Update (ua,va) with the vertical diffusive trend
       !!---------------------------------------------------------------------
-      !! * Arguments
-      INTEGER , INTENT( in ) ::   kt                           ! ocean time-step index
-      REAL(wp), INTENT( in ) ::   p2dt                         ! time-step 
-
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk, jl                              ! dummy loop indices
-      REAL(wp) ::   zrau0r, zlavmr, zua, zva                   ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpk) ::   zwx, zwy, zwz, zww     ! temporary workspace arrays
+      INTEGER , INTENT( in ) ::   kt     ! ocean time-step index
+      REAL(wp), INTENT( in ) ::   p2dt   ! time-step 
+      !! 
+      INTEGER  ::   ji, jj, jk, jl     ! dummy loop indices
+      REAL(wp) ::   zlavmr, zua, zva   ! temporary scalars
+      REAL(wp), DIMENSION(jpi,jpk) ::   zwx, zwy, zwz, zww   ! temporary workspace arrays
       !!----------------------------------------------------------------------
 
@@ -72,5 +70,4 @@
       ! Local constant initialization
       ! -----------------------------
-      zrau0r = 1. / rau0                                   ! inverse of the reference density
       zlavmr = 1. / float( nn_zdfexp )                      ! inverse of the number of sub time step
 
@@ -81,6 +78,6 @@
          ! Surface boundary condition
          DO ji = 2, jpim1
-            zwy(ji,1) = utau(ji,jj) * zrau0r
-            zww(ji,1) = vtau(ji,jj) * zrau0r
+            zwy(ji,1) = utau(ji,jj) * r1_rau0
+            zww(ji,1) = vtau(ji,jj) * r1_rau0
          END DO  
 
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 1857)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 1858)
@@ -59,20 +59,16 @@
       !!         ua = ua + dz( avmu dz(u) )
       !!
-      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive
-      !!               mixing trend.
+      !! ** Action : -  (ua,va) updated with the after vertical diffusive mixing trend
       !!---------------------------------------------------------------------
-      !! * Modules used
-      USE oce, ONLY :  zwd   => ta,   &                ! use ta as workspace
-                       zws   => sa                     ! use sa as workspace
-
-      !! * Arguments
-      INTEGER , INTENT( in ) ::   kt                   ! ocean time-step index
-      REAL(wp), INTENT( in ) ::  p2dt                  ! vertical profile of tracer time-step
-
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk                          ! dummy loop indices
-      REAL(wp) ::   zrau0r, z2dtf, zcoef, zzws, zrhs   ! temporary scalars
-      REAL(wp) ::   zzwi                               ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! temporary workspace arrays
+      USE oce, ONLY :  zwd   => ta   ! use ta as workspace
+      USE oce, ONLY :  zws   => sa   !  -  sa         -
+      !!
+      INTEGER , INTENT( in ) ::   kt    ! ocean time-step index
+      REAL(wp), INTENT( in ) ::  p2dt   ! vertical profile of tracer time-step
+      !!
+      INTEGER ::   ji, jj, jk                     ! dummy loop indices
+      REAL(wp) ::   z2dtf, zcoef, zzws, zrhs      ! temporary scalars
+      REAL(wp) ::   zzwi                          !    -         -
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi   ! 3D workspace 
       !!----------------------------------------------------------------------
 
@@ -82,8 +78,4 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
       ENDIF
-
-      ! 0. Local constant initialization
-      ! --------------------------------
-      zrau0r = 1. / rau0      ! inverse of the reference density
 
       ! 1. Vertical diffusion on u
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 1857)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 1858)
@@ -114,5 +114,4 @@
       REAL(wp) ::   zd , zc , zaw, za    !    -         -
       REAL(wp) ::   zb1, za1, zkw, zk0   !    -         -
-      REAL(wp) ::   zrau0r               !    -         -
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zws   ! temporary workspace
       !!----------------------------------------------------------------------
@@ -121,5 +120,4 @@
       !
       CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==!
-         zrau0r = 1.e0 / rau0
 !CDIR NOVERRCHK
          zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) )
@@ -162,5 +160,5 @@
                   ! masked in situ density anomaly
                   prd(ji,jj,jk) = (  zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  )    &
-                     &             - rau0  ) * zrau0r * tmask(ji,jj,jk)
+                     &             - rau0  ) * r1_rau0 * tmask(ji,jj,jk)
                END DO
             END DO
@@ -237,5 +235,5 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw   ! temporary 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), DIMENSION(jpi,jpj,jpk) ::   zws   ! 3D workspace
       !!----------------------------------------------------------------------
@@ -244,5 +242,4 @@
       !
       CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==!
-         zrau0r = 1.e0 / rau0
 !CDIR NOVERRCHK
          zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) )
@@ -288,5 +285,5 @@
                   ! masked in situ density anomaly
                   prd(ji,jj,jk) = (  zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  )    &
-                     &             - rau0  ) * zrau0r * tmask(ji,jj,jk)
+                     &             - rau0  ) * r1_rau0 * tmask(ji,jj,jk)
                END DO
             END DO
@@ -296,5 +293,5 @@
          DO jk = 1, jpkm1
             prd  (:,:,jk) = ( 0.0285 - rn_alpha * ptem(:,:,jk) )        * tmask(:,:,jk)
-            prhop(:,:,jk) = ( 1.e0   +          prd (:,:,jk) ) * rau0 * tmask(:,:,jk)
+            prhop(:,:,jk) = ( 1.e0   +            prd (:,:,jk) ) * rau0 * tmask(:,:,jk)
          END DO
          !
@@ -302,5 +299,5 @@
          DO jk = 1, jpkm1
             prd  (:,:,jk) = ( rn_beta  * psal(:,:,jk) - rn_alpha * ptem(:,:,jk) )        * tmask(:,:,jk)
-            prhop(:,:,jk) = ( 1.e0   + prd (:,:,jk) )                         * rau0 * tmask(:,:,jk)
+            prhop(:,:,jk) = ( 1.e0   + prd (:,:,jk) )                             * rau0 * tmask(:,:,jk)
          END DO
          !
@@ -409,5 +406,5 @@
                !
                ! masked in situ density anomaly
-               prd(ji,jj) = ( zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 ) / rau0 * zmask
+               prd(ji,jj) = ( zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 ) * r1_rau0 * zmask
             END DO
          END DO
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/trabbc.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 1857)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 1858)
@@ -71,5 +71,5 @@
       !!----------------------------------------------------------------------
       USE oce, ONLY :   ztrdt => ua   ! use ua as 3D workspace   
-      USE oce, ONLY :   ztrds => va   ! use va as 3D workspace   
+      USE oce, ONLY :   ztrds => va   !  -  va     -      -  
       !!
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
@@ -98,5 +98,5 @@
             DO ji = 2, jpim1
 #endif
-               zqgh_trd = ro0cpr * qgh_trd0(ji,jj) / fse3t(ji,jj,nbotlevt(ji,jj))
+               zqgh_trd = r1_rau0_rcp * qgh_trd0(ji,jj) / fse3t(ji,jj,nbotlevt(ji,jj))
                ta(ji,jj,nbotlevt(ji,jj)) = ta(ji,jj,nbotlevt(ji,jj)) + zqgh_trd
             END DO
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 1857)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 1858)
@@ -120,5 +120,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ta(ji,jj,jk) = ta(ji,jj,jk) + ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk) 
+                  ta(ji,jj,jk) = ta(ji,jj,jk) + r1_rau0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
                END DO
             END DO
@@ -177,5 +177,5 @@
                !
                DO jk = 1, nksr                                        ! compute and add qsr trend to ta
-                  ta(:,:,jk) = ta(:,:,jk) + ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk)
+                  ta(:,:,jk) = ta(:,:,jk) + r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk)
                END DO
                zea(:,:,nksr+1:jpk) = 0.e0     ! below 400m set to zero
@@ -377,5 +377,5 @@
                !
                DO jk = 1, nksr
-                  etot3(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk)
+                  etot3(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk)
                END DO
                etot3(:,:,nksr+1:jpk) = 0.e0                ! below 400m set to zero
@@ -399,5 +399,5 @@
                      zc0 = rn_abs * EXP( -fsdepw(ji,jj,jk  )*zsi0r ) + (1.-rn_abs) * EXP( -fsdepw(ji,jj,jk  )*zsi1r )
                      zc1 = rn_abs * EXP( -fsdepw(ji,jj,jk+1)*zsi0r ) + (1.-rn_abs) * EXP( -fsdepw(ji,jj,jk+1)*zsi1r )
-                     etot3(ji,jj,jk) = ro0cpr * (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) / fse3t(ji,jj,jk)
+                     etot3(ji,jj,jk) = r1_rau0_rcp * ( zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
                   END DO
                END DO
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 1857)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 1858)
@@ -99,5 +99,5 @@
       !!----------------------------------------------------------------------
       USE oce, ONLY :   ztrdt => ua   ! use ua as 3D workspace   
-      USE oce, ONLY :   ztrds => va   ! use va as 3D workspace   
+      USE oce, ONLY :   ztrds => va   !  -  va     -      -   
       !!
       INTEGER, INTENT(in) ::   kt     ! ocean time-step index
@@ -132,9 +132,9 @@
 #endif
             IF( lk_vvl) THEN
-               zta = ro0cpr * qns(ji,jj) * zse3t &                   ! temperature : heat flux
+               zta = r1_rau0_rcp * qns(ji,jj) * zse3t &              ! temperature : heat flux
                 &    - emp(ji,jj) * zsrau * tn(ji,jj,1)  * zse3t     ! & cooling/heating effet of EMP flux
                zsa = 0.e0                                            ! No salinity concent./dilut. effect
             ELSE
-               zta = ro0cpr * qns(ji,jj) * zse3t     ! temperature : heat flux
+               zta = r1_rau0_rcp * qns(ji,jj) * zse3t                ! temperature : heat flux
                zsa = emps(ji,jj) * zsrau * sn(ji,jj,1)   * zse3t     ! salinity :  concent./dilut. effect
             ENDIF
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/ZDF/zdfkpp.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/ZDF/zdfkpp.F90	(revision 1857)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/ZDF/zdfkpp.F90	(revision 1858)
@@ -455,5 +455,5 @@
             zBo   (ji,jj) = grav * zthermal * qns(ji,jj) -  grav * zhalin * emps(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) * sn(ji,jj,1) * rcs ) * tmask(ji,jj,1)
