Index: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 4505)
+++ branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 4506)
@@ -266,7 +266,7 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fbif        !: Heat flux at the ice base
    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(:,:) ::   rdq_snw     !: Heat per grid cell area released by snow melt    [J/m2], >0 if to the ocean
    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(:,:) ::   rdq_ice     !: Heat per grid cell area taken or released by ice growth and melt [J/m2], >0 if to the ocean
    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 
@@ -284,5 +284,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_res     !: residual salt flux due to correction of ice thickness [PSU/m2/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhbri       !: heat flux due to brine rejection
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fheat_mec   !: heat flux associated with porous ridged ice formation [???]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fheat_res   !: residual heat flux due to correction of ice thickness
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fmmec       !: mass flux due to snow loss during compression         [Kg/m2/s]
@@ -455,6 +454,6 @@
          &      fdtcn    (jpi,jpj) , qdtcn  (jpi,jpj) , fstric (jpi,jpj) , fscmbq   (jpi,jpj) ,     &
          &      ffltbif  (jpi,jpj) , fsbbq  (jpi,jpj) , qfvbq  (jpi,jpj) , dmgwi    (jpi,jpj) ,     &
-         &      sfx_res  (jpi,jpj) , sfx_bri(jpi,jpj) , sfx_mec(jpi,jpj) , fheat_mec(jpi,jpj) ,     &
-         &      fhbri    (jpi,jpj) , fmmec  (jpi,jpj) , sfx_thd(jpi,jpj) , fhmec    (jpi,jpj) ,     &
+         &      sfx_res  (jpi,jpj) , sfx_bri(jpi,jpj) , sfx_mec(jpi,jpj) , fhbri    (jpi,jpj) ,     &
+         &      fmmec  (jpi,jpj) , sfx_thd(jpi,jpj) , fhmec    (jpi,jpj) ,                          &
          &      fheat_res(jpi,jpj)                                                            , STAT=ierr(ii) )
 
Index: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 4505)
+++ branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 4506)
@@ -363,4 +363,5 @@
             !-----------------------------------------------------------------------------!
             fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean
+            ! there is a unit problem for e_snow_mlt that we should fix
             fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice     ! heat sink for ocean
 
@@ -908,6 +909,7 @@
       INTEGER ::   ij                ! horizontal index, combines i and j loops
       INTEGER ::   icells            ! number of cells with aicen > puny
-      REAL(wp) ::   zindb, zsrdg2   ! local scalar
+      REAL(wp) ::   zindb, zsrdg2    ! local scalar
       REAL(wp) ::   hL, hR, farea, zdummy, zdummy0, ztmelts    ! left and right limits of integration
+      REAL(wp) ::   zsstK            ! SST in Kelvin
 
       INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices
@@ -1098,5 +1100,5 @@
             esrdg(ji,jj) = esnon_init(ji,jj,jl1) * afrac(ji,jj)
             srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por )
-            srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj)
+            srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) !! MV HC 2014 this line seems useless
 
             ! rafting volumes, heat contents ...
@@ -1120,6 +1122,5 @@
             ! Salinity
             !-------------
-            smsw(ji,jj)  = sss_m(ji,jj) * vsw(ji,jj) * rhoic / rau0       ! salt content of seawater frozen in voids
-
+            smsw(ji,jj)  = vsw(ji,jj) * sss_m(ji,jj)                      ! salt content of seawater frozen in voids !! MV HC2014
             zsrdg2       = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge
 
@@ -1128,6 +1129,8 @@
             !                                                             ! excess of salt is flushed into the ocean
             !sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice
+            ! MV HC 2014 this previous line seems ok too... Maybe units are not that good
 
             !rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic    ! gurvan: increase in ice volume du to seawater frozen in voids             
+            ! MV HC 2014 this previous line seems ok, i'm not sure at this moment of the sign convention
 
             !------------------------------------            
@@ -1158,4 +1161,6 @@
                &                                + rhosn*vsrft(ji,jj)*(1.0-fsnowrft)
 
+            ! MV HC 2014 - I think the units are wrong here... shit there is a 10-9 in the snow...
+            ! we should therefore do something
             esnow_mlt(ji,jj) = esnow_mlt(ji,jj) + esrdg(ji,jj)*(1.0-fsnowrdg)         &   !rafting included
                &                                + esrft(ji,jj)*(1.0-fsnowrft)          
@@ -1187,18 +1192,12 @@
                eirft(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj)
                e_i  (ji,jj,jk,jl1)  = e_i(ji,jj,jk,jl1) - erdg1(ji,jj,jk) - eirft(ji,jj,jk)
-               ! sea water heat content
-               ztmelts          = - tmut * sss_m(ji,jj) + rtt
-               ! heat content per unit volume
-               zdummy0          = - rcp * ( sst_m(ji,jj) + rt0 - rtt ) * vsw(ji,jj)
-
-               ! corrected sea water salinity
-               zindb  = MAX( 0._wp , SIGN( 1._wp , vsw(ji,jj) - epsi20 ) )
-               zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / MAX( ridge_por * vsw(ji,jj), epsi20 )
-
-               ztmelts          = - tmut * zdummy + rtt
-               ersw(ji,jj,jk)   = - rcp * ( ztmelts - rtt ) * vsw(ji,jj)
-
-               ! heat flux
-               fheat_mec(ji,jj) = fheat_mec(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) * r1_rdtice
+               
+               
+               ! enthalpy of the trapped seawater (J/m2, >0)
+               zsstK  = sst_m(ji,jj) + rt0
+               ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 )
+
+               ! heat flux to the ocean
+               rdq_ice(ji,jj) = rdq_ice(ji,jj) + ersw(ji,jj,jk) !! MV HC 2014
 
                ! Correct dimensions to avoid big values
@@ -1206,7 +1205,10 @@
 
                ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
+               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean 
+               !! MV HC 2014
                ersw (ji,jj,jk)  = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / REAL( nlay_i )
 
                erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk)
+
             END DO ! ij
          END DO !jk
Index: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4505)
+++ branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4506)
@@ -106,5 +106,6 @@
       INTEGER  ::   iflt, ial , iadv , ifral, ifrdv   !   -      -
       REAL(wp) ::   zinda, zemp, zemp_snow, zfmm      ! local scalars
-      REAL(wp) ::   zemp_snw                          !   -      -
+      REAL(wp) ::   zemp_snw, zqmass,   zcd           !   -      -  
+      REAL(wp) ::   zswitch                           !   -      -  
       REAL(wp) ::   zfcm1 , zfcm2                     !   -      -
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp     ! 2D/3D workspace
@@ -113,4 +114,10 @@
       
       IF( lk_cpl )   CALL wrk_alloc( jpi, jpj, jpl, zalb, zalbp )
+
+      SELECT CASE( nn_ice_embd )                 ! levitating or embedded sea-ice option                              
+        CASE( 0    )   ;   zswitch = 1           ! (0) standard levitating sea-ice : salt exchange only              
+        CASE( 1, 2 )   ;   zswitch = 0           ! (1) levitating sea-ice: salt and volume exchange but no pressure effect
+                                                 ! (2) embedded sea-ice : salt and volume fluxes and pressure            
+      END SELECT                                 !                                                                      
 
       !------------------------------------------!
@@ -172,5 +179,4 @@
                &    + ifrdv   * (       qfvbq(ji,jj) +             qdtcn(ji,jj) ) * r1_rdtice   &
                &    + fhmec(ji,jj)                                                              & ! snow melt when ridging
-               &    + fheat_mec(ji,jj)                                                          & ! ridge formation
                &    + fheat_res(ji,jj)                                                            ! residual heat flux
             ! qcmif   Energy needed to bring the ocean surface layer until its freezing (ok)
@@ -184,8 +190,10 @@
             fsbbq(ji,jj) = ( 1._wp - ( ifvt + iflt ) ) * fscmbq(ji,jj)     
 
-            ! used to compute the oceanic heat flux at the next time step
+            ! heat flux associated with ice-ocean mass exchange
+            zqmass = ( rdq_snw(ji,jj)                                         & 
+                 & + rdq_ice(ji,jj) * ( 1.- zswitch) ) * r1_rdtice   !  heat flux due to snow & ice heat content 
             qsr(ji,jj) = zfcm1                                       ! solar heat flux 
-            qns(ji,jj) = zfcm2 - fdtcn(ji,jj)                        ! non solar heat flux
-            !                           ! fdtcn : turbulent oceanic heat flux
+            qns(ji,jj) = zfcm2 - fdtcn(ji,jj) + zqmass               ! non solar heat flux   
+                               ! fdtcn : turbulent oceanic heat flux
          END DO
       END DO
Index: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4505)
+++ branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4506)
@@ -184,5 +184,8 @@
                &                            + fdtcn(ji,jj)                        &   ! turbulent ice-ocean heat
                &                            + fsbbq(ji,jj) * ( 1.0 - zinda )  )   &   ! residual heat from previous step
+               ! MV HC 2014 check that this is good.
+               ! We should remove the heat content of precip that has fallen on sea ice
                & - pfrld(ji,jj)**betas * sprecip(ji,jj) * lfus                    )   ! latent heat of sprecip melting
+               ! MV HC 2014 partie heat content manque
             !
             ! Positive heat budget is used for bottom ablation
@@ -282,5 +285,7 @@
             CALL tab_2d_1d( nbpb, qldif_1d   (1:nbpb), qldif           , jpi, jpj, npb(1:nbpb) )
             CALL tab_2d_1d( nbpb, rdm_ice_1d (1:nbpb), rdm_ice         , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, rdq_ice_1d (1:nbpb), rdq_ice         , jpi, jpj, npb(1:nbpb) )
             CALL tab_2d_1d( nbpb, rdm_snw_1d (1:nbpb), rdm_snw         , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, rdq_snw_1d (1:nbpb), rdq_snw         , jpi, jpj, npb(1:nbpb) )
             CALL tab_2d_1d( nbpb, dmgwi_1d   (1:nbpb), dmgwi           , jpi, jpj, npb(1:nbpb) )
             CALL tab_2d_1d( nbpb, qlbbq_1d   (1:nbpb), zqlbsbq         , jpi, jpj, npb(1:nbpb) )
@@ -349,5 +354,7 @@
                CALL tab_1d_2d( nbpb, qfvbq         , npb, qfvbq_1d  (1:nbpb)   , jpi, jpj )
                CALL tab_1d_2d( nbpb, rdm_ice       , npb, rdm_ice_1d(1:nbpb)   , jpi, jpj )
+               CALL tab_1d_2d( nbpb, rdq_ice       , npb, rdq_ice_1d(1:nbpb)   , jpi, jpj )
                CALL tab_1d_2d( nbpb, rdm_snw       , npb, rdm_snw_1d(1:nbpb)   , jpi, jpj )
+               CALL tab_1d_2d( nbpb, rdq_snw       , npb, rdq_snw_1d(1:nbpb)   , jpi, jpj )
                CALL tab_1d_2d( nbpb, dmgwi         , npb, dmgwi_1d  (1:nbpb)   , jpi, jpj )
                CALL tab_1d_2d( nbpb, rdvosif       , npb, dvsbq_1d  (1:nbpb)   , jpi, jpj )
Index: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 4505)
+++ branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 4506)
@@ -74,4 +74,5 @@
       INTEGER  ::   ji , jk        ! dummy loop indices
       INTEGER  ::   ii, ij         ! 2D corresponding indices to ji
+      INTEGER  ::   i_use          ! debugging flag
       INTEGER  ::   isnow          ! switch for presence (1) or absence (0) of snow
       INTEGER  ::   isnowic        ! snow ice formation not
@@ -85,11 +86,18 @@
       REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
       REAL(wp) ::   zcoeff       ! dummy argument for snowfall partitioning over ice and leads
-      REAL(wp) ::   zsm_snowice  ! snow-ice salinity
+      REAL(wp) ::   zs_snic  ! snow-ice salinity
       REAL(wp) ::   zswi1        ! switch for computation of bottom salinity
       REAL(wp) ::   zswi12       ! switch for computation of bottom salinity
       REAL(wp) ::   zswi2        ! switch for computation of bottom salinity
       REAL(wp) ::   zgrr         ! bottom growth rate
-      REAL(wp) ::   ztform       ! bottom formation temperature
-      !
+      REAL(wp) ::   zt_i_new     ! bottom formation temperature
+
+      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
+      REAL(wp) ::   zEi          ! specific enthalpy of sea ice (J/kg)
+      REAL(wp) ::   zEw          ! specific enthalpy of exchanged water (J/kg)
+      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg)
+      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean
+      REAL(wp) ::   zsstK        ! SST in Kelvin
+
       REAL(wp), POINTER, DIMENSION(:) ::   zh_i        ! ice layer thickness
       REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness
@@ -119,5 +127,5 @@
 
       ! mass and salt flux (clem)
-      REAL(wp) :: zdvres, zdvsur, zdvbot
+      REAL(wp) :: zdvres, zdvsur, zdvbot, zswitch_sal
       REAL(wp), POINTER, DIMENSION(:) ::   zviold, zvsold   ! old ice volume...
 
@@ -129,4 +137,10 @@
       REAL(wp), POINTER, DIMENSION(:) ::   zfbase, zdq_i
       !!------------------------------------------------------------------
+
+      ! Discriminate between varying salinity (num_sal=2) and prescribed cases (other values)
+      SELECT CASE( num_sal )                       ! varying salinity or not
+         CASE( 1, 3, 4 ) ;   zswitch_sal = 0       ! prescribed salinity profile
+         CASE( 2 )       ;   zswitch_sal = 1       ! varying salinity profile
+      END SELECT
 
       CALL wrk_alloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i )
@@ -222,5 +236,5 @@
       DO ji = kideb, kiut
          ! tatm_ice is now in K
-         zqprec   (ji)   =  rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus )  
+         zqprec   (ji)   =  rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )  
          zqfont_su(ji)   =  z_f_surf(ji) * rdt_ice
          zdeltah  (ji,1) =  MIN( 0.e0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) )
@@ -272,18 +286,34 @@
       DO jk = 1, nlay_i
          DO ji = kideb, kiut 
-            !                                                    ! melt of layer jk
-            zdeltah  (ji,jk) = - zqfont_su(ji) / q_i_b(ji,jk)
-            !                                                    ! recompute heat available
-            zqfont_su(ji   ) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk) 
-            !                                                    ! melt of layer jk cannot be higher than its thickness
-            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) )
-            !                                                    ! update surface melt
-            dh_i_surf(ji   ) = dh_i_surf(ji) + zdeltah(ji,jk) 
-            !                                                    ! for energy conservation
-            zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice
-            !
-            ! clem
-            sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    &
+            zEi            = - q_i_b(ji,jk) / rhoic                ! Specific enthalpy of layer k [J/kg, <0]
+
+            ztmelts        = - tmut * s_i_b(ji,jk) + rtt           ! Melting point of layer k [K]
+
+            zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0]
+
+            zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0
+
+            zfmdt          = - zqfont_su(ji) / zdE                 ! Mass flux to the ocean [kg/m2, >0]
+
+            zdeltah(ji,jk) = - zfmdt / rhoic                       ! Melt of layer jk [m, <0]
+
+            zqfont_su(ji)  = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * rhoic * ( - zdE )
+                                                                   ! Energy remaining in case of melting of the full layer [J/m2, >0]
+            zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0]
+
+            dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt
+
+            zfmdt          = - rhoic*zdeltah(ji,jk)                ! Recompute mass flux [kg/m2, >0]
+
+            zQm            = MAX ( zfmdt, 0._wp ) * zEw            ! Energy of the melt water sent to the ocean [J/m2, <0]
+
+            ! Contribution to salt flux 
+            sfx_thd_1d(ji)   = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    &
                &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice
+            ! Contribution to heat flux
+            rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm * a_i_b(ji)
+            
+            ! Conservation test
+            zdq_i(ji)      = zdq_i(ji) + zdeltah(ji,jk) * rhoic * zdE * r1_rdtice  ! ? still don't know
          END DO
       END DO
@@ -364,94 +394,84 @@
       !------------------------------------------------------------------------------!
       !
-      ! Ice basal growth / melt is given by the ratio of heat budget over basal
-      ! ice heat content.  Basal heat budget is given by the difference between
-      ! the inner conductive flux  (fc_bo_i), from the open water heat flux 
+      !------------------
+      ! 4.1 Basal growth 
+      !------------------
+      ! Basal growth is driven by heat imbalance at the ice-ocean interface,
+      ! between the inner conductive flux  (fc_bo_i), from the open water heat flux 
       ! (qlbbqb) and the turbulent ocean flux (fbif). 
       ! fc_bo_i is positive downwards. fbif and qlbbq are positive to the ice 
 
-      !-----------------------------------------------------
-      ! 4.1 Basal growth - (a) salinity not varying in time 
-      !-----------------------------------------------------
-      IF(  num_sal /= 2  ) THEN   ! ice salinity constant in time
+      ! If salinity varies in time, an iterative procedure is required, because
+      ! the involved quantities are inter-dependent.
+      ! Basal growth (dh_i_bott) depends upon new ice specific enthalpy (zEi),
+      ! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bott)
+      ! -> need for an iterative procedure, which converges quickly
+
+      IF ( num_sal == 2 ) THEN
+         num_iter_max = 5
+      ELSE
+         num_iter_max = 1
+      ENDIF
+
+      ! Initialize dh_i_bott
+      dh_i_bott(:) = 0.0e0
+
+      ! Iterative procedure
+      DO iter = 1, num_iter_max
          DO ji = kideb, kiut
-            IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0._wp  ) THEN
-               s_i_new(ji)         =  sm_i_b(ji)
-               ! Melting point in K
-               ztmelts             =  - tmut * s_i_new(ji) + rtt 
-               ! New ice heat content (Bitz and Lipscomb, 1999)
-               ztform              =  t_i_b(ji,nlay_i)  ! t_bo_b crashes in the
-               ! Baltic
-               q_i_b(ji,nlay_i+1)  = rhoic * (  cpic * ( ztmelts - ztform )                                &
-                  &                           + lfus * (  1.0 - ( ztmelts - rtt ) / ( ztform - rtt )  )    &
-                  &                           - rcp  * ( ztmelts - rtt )                                 )
-               ! Basal growth rate = - F*dt / q
-               dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
-               sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice
-            ENDIF
-         END DO
-      ENDIF
-
-      !-------------------------------------------------
-      ! 4.1 Basal growth - (b) salinity varying in time 
-      !-------------------------------------------------
-      IF(  num_sal == 2  ) THEN
-         ! the growth rate (dh_i_bott) is function of the new ice heat content (q_i_b(nlay_i+1)). 
-         ! q_i_b depends on the new ice salinity (snewice). 
-         ! snewice depends on dh_i_bott ; it converges quickly, so, no problem
-         ! See Vancoppenolle et al., OM08 for more info on this
-
-         ! Initial value (tested 1D, can be anything between 1 and 20)
-         num_iter_max = 4
-         s_i_new(:)   = 4.0
-
-         ! Iterative procedure
-         DO iter = 1, num_iter_max
-            DO ji = kideb, kiut
-               IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0  ) THEN
-                  ii = MOD( npb(ji) - 1, jpi ) + 1
-                  ij = ( npb(ji) - 1 ) / jpi + 1
-                  ! Melting point in K
-                  ztmelts             =   - tmut * s_i_new(ji) + rtt 
-                  ! New ice heat content (Bitz and Lipscomb, 1999)
-                  q_i_b(ji,nlay_i+1)  =  rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             &
-                     &                            + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )   &
-                     &                            - rcp * ( ztmelts-rtt )                                     )
-                  ! Bottom growth rate = - F*dt / q
-                  dh_i_bott(ji) =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)
-                  ! New ice salinity ( Cox and Weeks, JGR, 1988 )
-                  ! zswi2  (1) if dh_i_bott/rdt .GT. 3.6e-7
-                  ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8
-                  ! zswi1  (1) if dh_i_bott/rdt .LT. 2.0e-8
-                  zgrr   = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi13 ) )
-                  zswi2  = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) ) 
-                  zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
-                  zswi1  = 1. - zswi2 * zswi12 
-                  zfracs = zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
-                     &                   + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 
-                  zfracs = MIN( 0.5 , zfracs )
-                  s_i_new(ji) = zfracs * sss_m(ii,ij)
-               ENDIF ! fc_bo_i
-            END DO ! ji
-         END DO ! iter
-
-         ! Final values
-         DO ji = kideb, kiut
-            IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN
-               ! New ice salinity must not exceed 20 psu
-               s_i_new(ji) = MIN( s_i_new(ji), s_i_max )
-               ! Metling point in K
-               ztmelts     =   - tmut * s_i_new(ji) + rtt 
-               ! New ice heat content (Bitz and Lipscomb, 1999)
-               q_i_b(ji,nlay_i+1)  =  rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                              &
-                  &                            + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )    &
-                  &                            - rcp * ( ztmelts - rtt )                                    )
-               ! Basal growth rate = - F*dt / q
-               dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)
-               ! Salinity update
-               ! entrapment during bottom growth
-               sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice
-            ENDIF ! heat budget
-         END DO
-      ENDIF
+
+            IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0  ) THEN
+
+               ! New bottom ice salinity (Cox & Weeks, JGR88 )
+               !--- zswi1  if dh/dt < 2.0e-8
+               !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7 
+               !--- zswi2  if dh/dt > 3.6e-7
+               zgrr               = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi13 ) )
+               zswi2              = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) )
+               zswi12             = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
+               zswi1              = 1. - zswi2 * zswi12
+               zfracs             = MIN ( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
+                  &               + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 )
+
+               ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
+
+               s_i_new(ji)        = zswitch_sal * zfracs * sss_m(ii,ij)  &  ! New ice salinity
+                                  + ( 1. - zswitch_sal ) * sm_i_b(ji) 
+               ! New ice growth
+               ztmelts            = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K)
+
+               zt_i_new           = zswitch_sal * t_bo_b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)
+               
+               zEi                = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)      
+                  &               - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) )   &
+                  &               + rcp  * ( ztmelts-rtt )          
+
+               zEw                = rcp  * ( t_bo_b(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
+
+               zdE                = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
+
+               dh_i_bott(ji)      = rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / ( zdE * rhoic )
+
+               !!! not sure we still need the next line... useful to keep this in memory ? check limthd_ent...
+               q_i_b(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0)
+
+            ENDIF ! fc_bo_i
+         END DO ! ji
+      END DO ! iter
+
+      ! Contribution to Energy and Salt Fluxes
+      DO ji = kideb, kiut
+
+         zEw            = rcp  * ( t_bo_b(ji) - rt0 )                   ! Specific enthalpy of seawater (J/kg, < 0)
+
+         zfmdt          = - rhoic * dh_i_bott(ji)                       ! Mass flux x time step (kg/m2, < 0)
+
+         ! Contribution to energy flux to the ocean (J/m2)
+         rdq_ice_1d(ji) = rdq_ice_1d(ji) + zEw         * a_i_b(ji) * zfmdt
+
+         ! Contribution to salt flux  ()
+         sfx_thd_1d(ji) = sfx_thd_1d(ji) + s_i_new(ji) * a_i_b(ji) * zfmdt * r1_rdtice
+
+      END DO
 
       !----------------
@@ -465,5 +485,5 @@
          ! heat convergence at the surface > 0
          IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0._wp  ) THEN
-            s_i_new(ji)   =  s_i_b(ji,nlay_i)
+            s_i_new(ji) = s_i_b(ji,nlay_i) ! is this line still useful ?
             zqfont_bo(ji) =  rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) )
             zfbase(ji)    =  zqfont_bo(ji) * r1_rdtice     ! heat conservation test
@@ -476,19 +496,53 @@
          DO ji = kideb, kiut
             IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  >=  0._wp  ) THEN
-               ztmelts = - tmut * s_i_b(ji,jk) + rtt 
-               IF( t_i_b(ji,jk) >= ztmelts ) THEN   !!gm : a comment is needed
-                  zdeltah   (ji,jk) = - zh_i(ji)
-                  dh_i_bott (ji   ) = dh_i_bott(ji) + zdeltah(ji,jk)
-                  zinnermelt(ji   ) = 1._wp
-               ELSE                                  ! normal ablation
-                  zdeltah  (ji,jk) = - zqfont_bo(ji) / q_i_b(ji,jk)
-                  zqfont_bo(ji   ) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk)
-                  zdeltah  (ji,jk) = MAX(zdeltah(ji,jk), - zh_i(ji) )
-                  dh_i_bott(ji   ) = dh_i_bott(ji) + zdeltah(ji,jk)
-                  zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice
+
+               ztmelts = - tmut * s_i_b(ji,jk) + rtt  ! Melting point of layer jk (K)
+
+               IF( t_i_b(ji,jk) >= ztmelts ) THEN !!! Internal melting
+
+                  zdeltah   (ji,jk) = - zh_i(ji)  ! internal melting occurs when the internal temperature is above freezing     
+                                                  ! this should normally not happen, but sometimes, heat diffusion leads to this
+
+                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk)
+
+                  zinnermelt(ji)    = 1._wp
+
+                  zQm               = 0.          ! Not sure which specific enthalpy we should use here (MV HC 2014)
+                                                  ! If that happens, heat is probably not well counted
+                                                  ! Put zero by default, but more bug analysis should be done to investigate this case
+
+               ELSE                               !!! Basal melting
+
+                  zEi               = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)
+
+                  zEw               = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0)
+
+                  zdE               = zEi - zEw              ! Specific enthalpy difference   (J/kg, <0)
+
+                  zfmdt             = - zqfont_bo(ji) / zdE  ! Mass flux x time step (kg/m2, >0)
+
+                  zdeltah(ji,jk)    = - zfmdt / rhoic        ! Gross thickness change
+
+                  zqfont_bo(ji)     = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * rhoic * ( - zdE ) ! Update heat available
+
+                  zdeltah(ji,jk)    = MAX( zdeltah(ji,jk), - zh_i(ji) ) ! Update thickness change
+
+                  dh_i_bott(ji)     = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt
+
+                  zfmdt             = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step
+
+                  zQm               = MAX ( zfmdt , 0.0 ) * zEw         ! Heat exchanged with ocean
+
+                  zdq_i(ji)         = zdq_i(ji) + zdeltah(ji,jk) * rhoic * zdE * r1_rdtice ! for heat conservation 
+
                ENDIF
-               ! clem: contribution to salt flux
+
+               ! Contribution to salt flux
                sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    &
                     &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice
+
+               ! Contribution to energy flux to the ocean (J/m2)
+               rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm * a_i_b(ji)
+
             ENDIF
          END DO ! ji
@@ -571,5 +625,5 @@
          dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji)
 
-        ! remaining heat
+         ! remaining heat
          zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) +  zqfont_bo(ji) ) 
 
@@ -661,5 +715,6 @@
       DO ji = kideb, kiut
          ht_i_b(ji) = zhgnew(ji)
-      END DO  ! ji
+      END DO ! ji 
+      
       !
       !------------------------------------------------------------------------------|
@@ -681,12 +736,10 @@
          !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew     - ht_s_b(ji) ) * rhosn 
 
-         !        Equivalent salt flux (1) Snow-ice formation component
-         !        -----------------------------------------------------
-         ii = MOD( npb(ji) - 1, jpi ) + 1
-         ij =    ( npb(ji) - 1 ) / jpi + 1
-
-         IF( num_sal == 2 ) THEN   ;   zsm_snowice = sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic
-         ELSE                      ;   zsm_snowice = sm_i_b(ji)   
-         ENDIF
+         ! Salinity of snow ice
+         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
+
+         zs_snic = zswitch_sal          * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic     &
+                 + ( 1. - zswitch_sal ) * sm_i_b(ji)
+
          ! entrapment during snow ice formation
          ! clem: new salinity difference stored (to be used in limthd_ent.F90)
@@ -694,5 +747,5 @@
             i_ice_switch = MAX( 0._wp , SIGN( 1._wp , zhgnew(ji) - epsi10 ) )
             ! salinity dif due to snow-ice formation
-            dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch     
+            dsm_i_si_1d(ji) = ( zs_snic - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch     
             ! salinity dif due to bottom growth 
             IF (  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  < 0._wp ) THEN
@@ -709,13 +762,19 @@
 
          ! diagnostic ( snow ice growth )
-         ii = MOD( npb(ji) - 1, jpi ) + 1
-         ij =    ( npb(ji) - 1 ) / jpi + 1
+         ii = MOD( npb(ji) - 1, jpi ) + 1 !MV HC 2014 useless
+         ij =    ( npb(ji) - 1 ) / jpi + 1 ! MV HC 2014 useless
          diag_sni_gr(ii,ij)  = diag_sni_gr(ii,ij) + dh_snowice(ji)*a_i_b(ji) * r1_rdtice
-         !
-         ! salt flux
-         sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice
-         !--------------------------------
-         ! Update mass fluxes (clem)
-         !--------------------------------
+
+         ! Contribution to energy flux to the ocean [J/m2], >0
+         zfmdt          = ( rhosn - rhoic ) * MAX( dh_snowice(ji), 0.0 )    ! <0
+         zsstK          = sst_m(ii,ij) + rt0                                
+         zEw            = rcp * ( zsstK - rt0 )
+         zQm            = zfmdt * zEw 
+         rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm * a_i_b(ji)
+
+         ! Contribution to salt flux
+         sfx_thd_1d(ji) = sfx_thd_1d(ji) + sss_m(ii,ij) * a_i_b(ji) * zfmdt * r1_rdtice 
+          
+         ! Contribution to mass fluxes
          rdm_ice_1d(ji) = rdm_ice_1d(ji) + ( a_i_b(ji) * ht_i_b(ji) - zviold(ji) ) * rhoic 
          rdm_snw_1d(ji) = rdm_snw_1d(ji) + ( a_i_b(ji) * ht_s_b(ji) - zvsold(ji) ) * rhosn 
Index: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
===================================================================
--- branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 4505)
+++ branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 4506)
@@ -22,4 +22,5 @@
    USE domain         !
    USE phycst         ! physical constants
+   USE sbc_oce        ! Surface boundary condition: ocean fields
    USE ice            ! LIM variables
    USE par_ice        ! LIM parameters
@@ -86,4 +87,5 @@
          ztmelts        ,   &  ! ice melting point
          zqsnic         ,   &  ! enthalpy of snow ice layer
+         zsstK          ,   &  ! SST in Kelvin
          zhsnow         ,   &  ! temporary snow thickness variable
          zswitch        ,   &  ! dummy switch argument
@@ -540,10 +542,11 @@
       DO ji = kideb, kiut
          ! energy of the flooding seawater
-         zqsnic = rau0 * rcp * ( rtt - t_bo_b(ji) ) * dh_snowice(ji) * &
-            (rhoic - rhosn) / rhoic * REAL(snicswi(ji)) ! generally positive
+         ii = MOD( npb(ji) - 1, jpi ) + 1
+         ij =    ( npb(ji) - 1 ) / jpi + 1
+         zsstK  = sst_m(ii,ij) + rt0
+         zqsnic = ( rhosn - rhoic ) * dh_snowice(ji) * rcp * ( zsstK - rt0 ) * REAL ( snicswi(ji) ) ! MV HC 2014
+         
          ! Heat conservation diagnostic
          qt_i_in(ji,jl) = qt_i_in(ji,jl) + zqsnic 
-
-         qldif_1d(ji)   = qldif_1d(ji) + zqsnic * a_i_b(ji)
 
          ! enthalpy of the newly formed snow-ice layer
Index: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4505)
+++ branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4506)
@@ -81,5 +81,10 @@
       LOGICAL  ::   iterate_frazil   ! iterate frazil ice collection thickness
       CHARACTER (len = 15) :: fieldid
-      !
+
+      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean)
+      REAL(wp) ::   zEi          ! sea ice specific enthalpy (J/kg)
+      REAL(wp) ::   zEw          ! seawater specific enthalpy (J/kg)
+      REAL(wp) ::   zfmdt        ! mass flux x time step (kg/m2, >0 towards ocean)
+       
       INTEGER , POINTER, DIMENSION(:) ::   zcatac      ! indexes of categories where new ice grows
       REAL(wp), POINTER, DIMENSION(:) ::   zswinew     ! switch for new ice or not
@@ -323,4 +328,5 @@
          CALL tab_2d_1d( nbpac, sfx_thd_1d(1:nbpac)     , sfx_thd, jpi, jpj, npac(1:nbpac) )
          CALL tab_2d_1d( nbpac, rdm_ice_1d(1:nbpac)     , rdm_ice, jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, rdq_ice_1d(1:nbpac)     , rdq_ice, jpi, jpj, npac(1:nbpac) )
          CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) )
          CALL tab_2d_1d( nbpac, zvrel_ac  (1:nbpac)     , zvrel  , jpi, jpj, npac(1:nbpac) )
@@ -365,4 +371,6 @@
                &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )   &
                &                       - rcp  *         ( ztmelts - rtt )  )
+            ! MV HC 2014 comment I dont see why this line below is here... ?
+            ! This implies that ze_newice gets to rhoic*Lfus if it was negative, but this should never happen
             ze_newice(ji) =   MAX( ze_newice(ji) , 0._wp )    &
                &          +   MAX(  0.0 , SIGN( 1.0 , - ze_newice(ji) )  ) * rhoic * lfus
@@ -386,5 +394,19 @@
          !-------------------
          DO ji = 1, nbpac
-            zv_newice(ji) = - zqbgow(ji) / ze_newice(ji)
+
+            zEi           = - ze_newice(ji) / rhoic                ! specific enthalpy of forming ice [J/kg]
+
+            zEw           = rcp * ( t_bo_b(ji) - rt0 )             ! specific enthalpy of seawater at t_bo_b [J/kg]
+
+            zdE           = zEi - zEw                              ! specific enthalpy difference [J/kg]
+
+            zfmdt         = - zqbgow(ji) / zdE                     ! Fm.dt [kg/m2] (<0) 
+
+            zv_newice(ji) = - zfmdt / rhoic
+
+            zQm           = zfmdt * zEw                            ! heat to the ocean >0 associated with mass flux  
+
+            ! Contribution to energy flux to the ocean [J/m2], >0  
+            rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm 
 
             ! A fraction zfrazb of frazil ice is accreted at the ice bottom
@@ -539,4 +561,5 @@
          END DO
 !!gm ???  why the previous do loop  if ocerwriten by the following one ?
+!! MV HC 2014 just because it'snot the same index and the expression is different
          DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
             DO ji = 1, nbpac
@@ -623,4 +646,5 @@
          CALL tab_1d_2d( nbpac, sfx_thd, npac(1:nbpac), sfx_thd_1d(1:nbpac), jpi, jpj )
          CALL tab_1d_2d( nbpac, rdm_ice, npac(1:nbpac), rdm_ice_1d(1:nbpac), jpi, jpj )
+         CALL tab_1d_2d( nbpac, rdq_ice, npac(1:nbpac), rdq_ice_1d(1:nbpac), jpi, jpj )
          !
       ENDIF ! nbpac > 0
Index: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
===================================================================
--- branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 4505)
+++ branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 4506)
@@ -448,5 +448,5 @@
          CALL prt_ctl(tab2d_1=fmmec  , clinfo1= ' lim_update1 : fmmec : ', tab2d_2=fhmec     , clinfo2= ' fhmec     : ')
          CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update1 : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ')
-         CALL prt_ctl(tab2d_1=fhbri  , clinfo1= ' lim_update1 : fhbri : ', tab2d_2=fheat_mec , clinfo2= ' fheat_mec : ')
+         CALL prt_ctl(tab2d_1=fhbri  , clinfo1= ' lim_update1 : fhbri : ') 
 
          CALL prt_ctl_info(' ')
Index: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
===================================================================
--- branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 4505)
+++ branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 4506)
@@ -598,5 +598,5 @@
          CALL prt_ctl(tab2d_1=fmmec  , clinfo1= ' lim_update2 : fmmec : ', tab2d_2=fhmec     , clinfo2= ' fhmec     : ')
          CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update2 : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ')
-         CALL prt_ctl(tab2d_1=fhbri  , clinfo1= ' lim_update2 : fhbri : ', tab2d_2=fheat_mec , clinfo2= ' fheat_mec : ')
+         CALL prt_ctl(tab2d_1=fhbri  , clinfo1= ' lim_update2 : fhbri : ')
 
          CALL prt_ctl_info(' ')
Index: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 4505)
+++ branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 4506)
@@ -66,5 +66,7 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fbif_1d       !: <==> the 2D  fbif
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdm_ice_1d    !: <==> the 2D  rdm_ice
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdq_ice_1d    !: <==> the 2D  rdq_ice
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdm_snw_1d    !: <==> the 2D  rdm_snw
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdq_snw_1d    !: <==> the 2D  rdq_snw
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qlbbq_1d      !: <==> the 2D  qlbsbq
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dmgwi_1d      !: <==> the 2D  dmgwi
@@ -164,4 +166,5 @@
       ALLOCATE( sprecip_1d (jpij) , frld_1d    (jpij) , at_i_b     (jpij) ,     &
          &      fbif_1d    (jpij) , rdm_ice_1d (jpij) , rdm_snw_1d (jpij) ,     &
+         &                          rdq_ice_1d (jpij) , rdq_snw_1d (jpij) ,     & !MV HC 2014
          &      qlbbq_1d   (jpij) , dmgwi_1d   (jpij) , dvsbq_1d   (jpij) ,     &
          &      dvbbq_1d   (jpij) , dvlbq_1d   (jpij) , dvnbq_1d   (jpij) ,     &
Index: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
===================================================================
--- branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 4505)
+++ branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 4506)
@@ -54,6 +54,6 @@
    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 ::   rcp      =    4.e3_wp       !: ocean specific heat           [J/kg/K]
+   REAL(wp), PUBLIC ::   r1_rcp                      !: = 1. / rcp                    [kg.K/J]
    REAL(wp), PUBLIC ::   r1_rau0_rcp                 !: = 1. / ( rau0 * rcp )
 
@@ -69,10 +69,10 @@
 #if defined key_lim3 || defined key_cice
    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 ::   rcdic    =    2.034396_wp   !: thermal conductivity of fresh ice                     [W/m/K]
+   REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow                          [W/m/K] 
+   REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat for ice                                 [J/kg/K]
    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 ::   tmut     =    0.054_wp      !: decrease of seawater meltpoint with salinity          [degC/ppt]
    REAL(wp), PUBLIC ::   xlsn                        !: = lfus*rhosn (volumetric latent heat fusion of snow)  [J/m3]
 #else
Index: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 4505)
+++ branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 4506)
@@ -136,5 +136,5 @@
       REAL(wp) ::   zcoef   ! local scalar
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice_os, zalb_ice_cs  ! albedo of the ice under overcast/clear sky
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice      ! mean albedo of ice (for coupled)
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice      ! mean albedo of ice 
 
       REAL(wp), POINTER, DIMENSION(:,:) :: zalb_ice_all    ! Mean albedo over all categories
@@ -152,8 +152,7 @@
       IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim')
 
-      CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs )
+      CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs, zalb_ice )
 
 #if defined key_coupled
-      IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_alloc( jpi,jpj,jpl, zalb_ice)
       IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) &
          &   CALL wrk_alloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all)
@@ -168,5 +167,5 @@
          !
          IF( ln_nicep ) THEN      ! control print at a given point
-            jiindx = 177   ;   jjindx = 112
+            jiindx = 14    ;   jjindx =  25
             IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx
          ENDIF
@@ -310,5 +309,6 @@
          sfx    (:,:) = 0._wp   ;   sfx_thd  (:,:) = 0._wp
          sfx_bri(:,:) = 0._wp   ;   sfx_mec  (:,:) = 0._wp   ;   sfx_res  (:,:) = 0._wp
-         fhbri  (:,:) = 0._wp   ;   fheat_mec(:,:) = 0._wp   ;   fheat_res(:,:) = 0._wp
+         fhbri  (:,:) = 0._wp   ;   
+         fheat_res(:,:) = 0._wp 
          fhmec  (:,:) = 0._wp   ;   
          fmmec  (:,:) = 0._wp
@@ -327,4 +327,6 @@
          rdm_snw(:,:) = 0._wp   ! variation of snow mass per unit area
          rdm_ice(:,:) = 0._wp   ! variation of ice mass per unit area
+         rdq_snw(:,:) = 0._wp   ! heat flux associated with mass exchange (ice contribution)
+         rdq_ice(:,:) = 0._wp   ! heat flux associated with mass exchange (snow contribution)
          hicifp (:,:) = 0._wp   ! daily thermodynamic ice production. 
          !
@@ -421,9 +423,8 @@
       
 !!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
-      !
-      CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs )
+      CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs, zalb_ice )
 
 #if defined key_coupled
-      IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice)
+  
       IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) &
          &    CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all)
@@ -807,5 +808,4 @@
                WRITE(numout,*) ' fhmec      : ', fhmec    (ji,jj)
                WRITE(numout,*) ' fhbri      : ', fhbri    (ji,jj)
-               WRITE(numout,*) ' fheat_mec  : ', fheat_mec(ji,jj)
                WRITE(numout,*) 
                WRITE(numout,*) ' sst        : ', sst_m(ji,jj)  
