Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90	(revision 3605)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90	(revision 3625)
@@ -19,6 +19,6 @@
    PUBLIC    ice_alloc_2  !  Called in iceini_2.F90
 
-   INTEGER , PUBLIC ::   numit     !: ice iteration index
-   REAL(wp), PUBLIC ::   rdt_ice   !: ice time step
+   INTEGER , PUBLIC ::   numit        !: ice iteration index
+   REAL(wp), PUBLIC ::   rdt_ice      !: ice time step
 
    !                                                                     !!* namicerun read in iceini  *
@@ -98,6 +98,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qstoif        !: Energy stored in the brine pockets
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fbif          !: Heat flux at the ice base
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdmsnif       !: Variation of snow mass
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdmicif       !: Variation of ice mass
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdm_snw       !: Variation of snow mass over 1 time step           [Kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdq_snw       !: Heat content associated with rdm_snw              [J/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdm_ice       !: Variation of ice  mass over 1 time step           [Kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdq_ice       !: Heat content associated with rdm_ice              [J/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qldif         !: heat balance of the lead (or of the open ocean)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qcmif         !: Energy needed to freeze the ocean surface layer
@@ -153,5 +155,6 @@
 
       ALLOCATE(phicif(jpi,jpj) , pfrld  (jpi,jpj) , qstoif (jpi,jpj) ,     &
-         &     fbif  (jpi,jpj) , rdmsnif(jpi,jpj) , rdmicif(jpi,jpj) ,     &
+         &     fbif  (jpi,jpj) , rdm_snw(jpi,jpj) , rdq_snw(jpi,jpj) ,     &
+         &                       rdm_ice(jpi,jpj) , rdq_ice(jpi,jpj) ,     &
          &     qldif (jpi,jpj) , qcmif  (jpi,jpj) , fdtcn  (jpi,jpj) ,     &
          &     qdtcn (jpi,jpj) , thcm   (jpi,jpj)                    , STAT=ierr(4) )
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/iceini_2.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/iceini_2.F90	(revision 3605)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/iceini_2.F90	(revision 3625)
@@ -13,21 +13,22 @@
    !!   'key_lim2' :                                  LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   !!   ice_init_2       : sea-ice model initialization
-   !!   ice_run_2        : Definition some run parameter for ice model
+   !!   ice_init_2    : sea-ice model initialization
+   !!   ice_run_2     : Definition some run parameter for ice model
    !!----------------------------------------------------------------------
-   USE phycst           ! physical constants
-   USE dom_oce          ! ocean domain
-   USE sbc_oce          ! surface boundary condition: ocean
-   USE sbc_ice          ! LIM2 surface boundary condition
-   USE dom_ice_2        ! LIM2 ice domain
-   USE par_ice_2        ! LIM2 parameters
-   USE thd_ice_2        ! LIM2 thermodynamical variables
-   USE ice_2            ! LIM2 ice variable
-   USE limmsh_2         ! LIM2 mesh
-   USE limistate_2      ! LIM2 initial state
-   USE limrst_2         ! LIM2 restart
-   USE limsbc_2         ! LIM2 surface boundary condition
-   USE in_out_manager   ! I/O manager
-   USE lib_mpp          ! MPP library
+   USE phycst         ! physical constants
+   USE dom_oce        ! ocean domain
+   USE sbc_oce        ! surface boundary condition: ocean
+   USE sbc_ice        ! LIM2 surface boundary condition
+   USE dom_ice_2      ! LIM2 ice domain
+   USE par_ice_2      ! LIM2 parameters
+   USE thd_ice_2      ! LIM2 thermodynamical variables
+   USE ice_2          ! LIM2 ice variable
+   USE limmsh_2       ! LIM2 mesh
+   USE limistate_2    ! LIM2 initial state
+   USE limrst_2       ! LIM2 restart
+   USE limsbc_2       ! LIM2 surface boundary condition
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 
    IMPLICIT NONE
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limadv_2.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limadv_2.F90	(revision 3605)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limadv_2.F90	(revision 3625)
@@ -14,6 +14,6 @@
    !!   'key_lim2'                                    LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   !!   lim_adv_x_2  : advection of sea ice on x axis
-   !!   lim_adv_y_2  : advection of sea ice on y axis
+   !!   lim_adv_x_2   : advection of sea ice on x axis
+   !!   lim_adv_y_2   : advection of sea ice on y axis
    !!----------------------------------------------------------------------
    USE dom_oce
@@ -21,8 +21,9 @@
    USE ice_2
    USE lbclnk
-   USE in_out_manager     ! I/O manager
-   USE lib_mpp            ! MPP library
-   USE wrk_nemo           ! work arrays
-   USE prtctl             ! Print control
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE wrk_nemo       ! work arrays
+   USE prtctl         ! Print control
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 
    IMPLICIT NONE
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limdia_2.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limdia_2.F90	(revision 3605)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limdia_2.F90	(revision 3625)
@@ -24,4 +24,5 @@
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! MPP library
+   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 
    IMPLICIT NONE
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limdmp_2.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limdmp_2.F90	(revision 3605)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limdmp_2.F90	(revision 3625)
@@ -19,4 +19,5 @@
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! MPP library
+   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 
    IMPLICIT NONE
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90	(revision 3605)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90	(revision 3625)
@@ -31,4 +31,5 @@
    USE in_out_manager   ! I/O manager
    USE prtctl           ! Print control
+   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 
    IMPLICIT NONE
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limhdf_2.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limhdf_2.F90	(revision 3605)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limhdf_2.F90	(revision 3625)
@@ -21,4 +21,5 @@
    USE prtctl           ! Print control
    USE in_out_manager   ! I/O manager
+   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 
    IMPLICIT NONE
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limistate_2.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limistate_2.F90	(revision 3605)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limistate_2.F90	(revision 3625)
@@ -27,4 +27,5 @@
    USE iom
    USE in_out_manager
+   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 
    IMPLICIT NONE
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90	(revision 3605)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90	(revision 3625)
@@ -23,4 +23,5 @@
    USE wrk_nemo         ! work arrays
 #endif
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 
    IMPLICIT NONE
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90	(revision 3605)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90	(revision 3625)
@@ -30,4 +30,6 @@
    USE in_out_manager ! I/O manager
    USE prtctl         ! Print control
+   USE oce     , ONLY : snwice_mass, snwice_mass_b
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 
    IMPLICIT NONE
@@ -80,8 +82,9 @@
       REAL(wp) ::   zs21_11, zs21_12, zs21_21, zs21_22
       REAL(wp) ::   zs22_11, zs22_12, zs22_21, zs22_22
+      REAL(wp) ::   zintb, zintn
       REAL(wp), POINTER, DIMENSION(:,:) ::   zfrld, zmass, zcorl
       REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct, za2ct, zresr
       REAL(wp), POINTER, DIMENSION(:,:) ::   zc1u, zc1v, zc2u, zc2v
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zsang
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zsang, zpice
       REAL(wp), POINTER, DIMENSION(:,:) ::   zu0, zv0
       REAL(wp), POINTER, DIMENSION(:,:) ::   zu_n, zv_n
@@ -93,5 +96,5 @@
       
       CALL wrk_alloc( jpi,jpj, zfrld, zmass, zcorl, za1ct, za2ct, zresr )
-      CALL wrk_alloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang )
+      CALL wrk_alloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang, zpice )
       CALL wrk_alloc( jpi,jpj+2, zu0, zv0, zu_n, zv_n, zu_a, zv_a, zviszeta, zviseta, kjstart = 0 )
       CALL wrk_alloc( jpi,jpj+2, zzfrld, zztms, zi1, zi2, zmasst, zpresh, kjstart = 0 )
@@ -129,4 +132,20 @@
 !i    zviszeta(:,jpj+1) = 0._wp    ;    zviseta(:,jpj+1) = 0._wp
 
+      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==!
+          !
+          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
+          !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
+         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
+          !
+          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
+          !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
+         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
+          !
+         zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0
+          !
+         !
+      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
+         zpice(:,:) = ssh_m(:,:)
+      ENDIF
 
       ! Ice mass, ice strength, and wind stress at the center            |
@@ -196,8 +215,8 @@
 
             ! Gradient of the sea surface height
-            zgsshx =  (   (ssh_m(ji  ,jj  ) - ssh_m(ji-1,jj  ))/e1u(ji-1,jj  )   &
-               &       +  (ssh_m(ji  ,jj-1) - ssh_m(ji-1,jj-1))/e1u(ji-1,jj-1)   ) * 0.5_wp
-            zgsshy =  (   (ssh_m(ji  ,jj  ) - ssh_m(ji  ,jj-1))/e2v(ji  ,jj-1)   &
-               &       +  (ssh_m(ji-1,jj  ) - ssh_m(ji-1,jj-1))/e2v(ji-1,jj-1)   ) * 0.5_wp
+            zgsshx =  (   (zpice(ji  ,jj  ) - zpice(ji-1,jj  ))/e1u(ji-1,jj  )   &
+               &       +  (zpice(ji  ,jj-1) - zpice(ji-1,jj-1))/e1u(ji-1,jj-1)   ) * 0.5_wp
+            zgsshy =  (   (zpice(ji  ,jj  ) - zpice(ji  ,jj-1))/e2v(ji  ,jj-1)   &
+               &       +  (zpice(ji-1,jj  ) - zpice(ji-1,jj-1))/e2v(ji-1,jj-1)   ) * 0.5_wp
 
             ! Computation of the velocity field taking into account the ice-ice interaction.                                 
@@ -575,5 +594,5 @@
 
       CALL wrk_dealloc( jpi,jpj, zfrld, zmass, zcorl, za1ct, za2ct, zresr )
-      CALL wrk_dealloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang )
+      CALL wrk_dealloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang, zpice )
       CALL wrk_dealloc( jpi,jpj+2, zu0, zv0, zu_n, zv_n, zu_a, zv_a, zviszeta, zviseta, kjstart = 0 )
       CALL wrk_dealloc( jpi,jpj+2, zzfrld, zztms, zi1, zi2, zmasst, zpresh, kjstart = 0 )
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 3605)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 3625)
@@ -9,5 +9,6 @@
    !!            3.3  ! 2009-05 (G. Garric, C. Bricaud) addition of the lim2_evp case
    !!             -   ! 2010-11 (G. Madec) ice-ocean stress computed at each ocean time-step
-   !!            4.0  ! 2011-01  (A. R. Porter, STFC Daresbury) dynamical allocation
+   !!           3.3.1 ! 2011-01 (A. R. Porter, STFC Daresbury) dynamical allocation
+   !!            3.5  ! 2012-11 ((G. Madec, Y. Aksenov, A. Coward) salt and heat fluxes associated with e-p
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -28,5 +29,6 @@
    USE sbc_oce          ! surface boundary condition: ocean
    USE sbccpl
-
+   USE cpl_oasis3, ONLY : lk_cpl
+   USE oce       , ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass 
    USE albedo           ! albedo parameters
    USE lbclnk           ! ocean lateral boundary condition - MPP exchanges
@@ -37,5 +39,5 @@
    USE iom              ! I/O library
    USE prtctl           ! Print control
-   USE cpl_oasis3, ONLY : lk_cpl
+   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 
    IMPLICIT NONE
@@ -88,8 +90,8 @@
       !!              - Update the fluxes provided to the ocean
       !!     
-      !! ** Outputs : - qsr     : sea heat flux:     solar 
-      !!              - qns     : sea heat flux: non solar
-      !!              - emp     : freshwater budget: volume flux 
-      !!              - emps    : freshwater budget: concentration/dillution 
+      !! ** Outputs : - qsr     : sea heat flux    : solar 
+      !!              - qns     : sea heat flux    : non solar (including heat content of the mass flux)
+      !!              - emp     : freshwater budget: mass flux 
+      !!              - sfx     : freshwater budget: salt flux due to Freezing/Melting
       !!              - utau    : sea surface i-stress (ocean referential)
       !!              - vtau    : sea surface j-stress (ocean referential)
@@ -107,6 +109,8 @@
       INTEGER  ::   ifvt, i1mfr, idfr, iflt    !   -       -
       INTEGER  ::   ial, iadv, ifral, ifrdv    !   -       -
-      REAL(wp) ::   zqsr, zqns, zfm            ! local scalars
-      REAL(wp) ::   zinda, zfons, zemp         !   -      -
+      REAL(wp) ::   zqsr,     zqns,   zfmm     ! local scalars
+      REAL(wp) ::   zinda,    zfsalt, zemp     !   -      -
+      REAL(wp) ::   zemp_snw, zqhc,   zcd      !   -      -
+      REAL(wp) ::   zswitch                    !   -      -
       REAL(wp), POINTER, DIMENSION(:,:)   ::   zqnsoce       ! 2D workspace
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp   ! 2D/3D workspace
@@ -115,4 +119,10 @@
       CALL wrk_alloc( jpi, jpj, zqnsoce )
       CALL wrk_alloc( jpi, jpj, 1, 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                                 !    
 
       !------------------------------------------!
@@ -133,22 +143,24 @@
             ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv 
 
-!!$            zinda   = 1.0 - AINT( pfrld(ji,jj) )                   !   = 0. if pure ocean else 1. (at previous time)
-!!$
-!!$            i1mfr   = 1.0 - AINT(  frld(ji,jj) )                   !   = 0. if pure ocean else 1. (at current  time)
-!!$
-!!$            IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda      !   = 1. if (snow and no ice at previous time) else 0. ???
-!!$            ELSE                             ;   ifvt = 0.
+!!$            attempt to explain the tricky flags set above....
+!!$            zinda   = 1.0 - AINT( pfrld(ji,jj) )                   ! = 0. if ice-free ocean else 1. (after ice adv, but before ice thermo)
+!!$            i1mfr   = 1.0 - AINT(  frld(ji,jj) )                   ! = 0. if ice-free ocean else 1. (after ice thermo)
+!!$
+!!$            IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda      ! = zinda if previous thermodynamic step overmelted the ice???
+!!$            ELSE                             ;   ifvt = 0.         ! 
 !!$            ENDIF
 !!$
-!!$            IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases from previous to current
+!!$            IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases due to ice thermodynamics
 !!$            ELSE                                     ;   idfr = 1.   
 !!$            ENDIF
 !!$
-!!$            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )    !   = 1. if ice (not only snow) at previous and pure ocean at current
+!!$            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )    !   = 1. if ice (not only snow) at previous time and ice-free ocean currently
 !!$
 !!$            ial     = ifvt   * i1mfr    +    ( 1 - ifvt ) * idfr
+!!$                    = i1mfr if ifvt = 1 i.e. 
+!!$                    = idfr  if ifvt = 0
 !!$!                 snow no ice   ice         ice or nothing  lead fraction increases
 !!$!                 at previous   now           at previous
-!!$!                -> ice aera increases  ???         -> ice aera decreases ???
+!!$!                -> ice area increases  ???         -> ice area decreases ???
 !!$
 !!$            iadv    = ( 1  - i1mfr ) * zinda
@@ -174,15 +186,57 @@
 #endif            
             !  computation the non solar heat flux at ocean surface
-            zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr   &   ! part of the solar energy used in leads
-               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                            &
-               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice    &
-               &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) )                   * r1_rdtice 
-
-            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! ???
-            !
+            zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr                                              &   ! part of the solar energy used in leads
+               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                             &
+               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice  &
+               &       + ifrdv   * (       qfvbq(ji,jj) +             qdtcn(ji,jj) ) * r1_rdtice 
+
+            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! store residual heat flux (to put into the ocean at the next time-step)
+            zqhc = ( rdq_snw(ji,jj)                                     &
+                 & + rdq_ice(ji,jj) * ( 1.- zswitch) ) * r1_rdtice       ! heat flux due to snow ( & ice heat content, 
+            !                                                            !           if ice/ocean mass exchange active) 
             qsr  (ji,jj) = zqsr                                          ! solar heat flux 
-            qns  (ji,jj) = zqns - fdtcn(ji,jj)                           ! non solar heat flux
+            qns  (ji,jj) = zqns - fdtcn(ji,jj) + zqhc                    ! non solar heat flux 
+            !
+            !                          !------------------------------------------!
+            !                          !  mass and salt flux at the ocean surface !
+            !                          !------------------------------------------!
+            !
+            ! mass flux at the ocean-atmosphere interface (open ocean fraction = leads area)
+#if defined key_coupled
+            !                                                  ! coupled mode: 
+            zemp = + emp_tot(ji,jj)                            &     ! net mass flux over the grid cell (ice+ocean area)
+               &   - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )          ! minus the mass flux intercepted by sea-ice
+#else
+            !                                                  ! forced  mode: 
+            zemp = + emp(ji,jj)     *         frld(ji,jj)      &     ! mass flux over open ocean fraction 
+               &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &     ! liquid precip. over ice reaches directly the ocean
+               &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )          ! snow is intercepted by sea-ice (previous frld)
+#endif            
+            !
+            ! mass flux at the ocean/ice interface (sea ice fraction)
+            zemp_snw = rdm_snw(ji,jj) * r1_rdtice                    ! snow melting = pure water that enters the ocean
+            zfmm     = rdm_ice(ji,jj) * r1_rdtice                    ! Freezing minus Melting (F-M)
+
+            ! salt flux at the ice/ocean interface (sea ice fraction) [PSU*kg/m2/s]
+            zfsalt = - sice_0(ji,jj) * zfmm                          ! F-M salt exchange
+            zcd    =   soce_0(ji,jj) * zfmm                          ! concentration/dilution term due to F-M
+            !
+            ! salt flux only       : add concentration dilution term in salt flux  and no  F-M term in volume flux
+            ! salt and mass fluxes : non concentration dilution term in salt flux  and add F-M term in volume flux
+            sfx (ji,jj) = zfsalt +                  zswitch  * zcd   ! salt flux (+ C/D if no ice/ocean mass exchange)
+            emp (ji,jj) = zemp   + zemp_snw + ( 1.- zswitch) * zfmm  ! mass flux (+ F/M mass flux if ice/ocean mass exchange)
+            !
          END DO
       END DO
+      !                                !------------------------------------------!
+      !                                !    mass of snow and ice per unit area    !
+      !                                !------------------------------------------!
+      IF( nn_ice_embd /= 0 ) THEN      ! embedded sea-ice (mass required)
+         snwice_mass_b(:,:) = snwice_mass(:,:)                  ! save mass from the previous ice time step
+         !                                                      ! new mass per unit area
+         snwice_mass  (:,:) = tms(:,:) * ( rhosn * hsnif(:,:) + rhoic * hicif(:,:)  ) * ( 1.0 - frld(:,:) )
+         !                                                      ! time evolution of snow+ice mass
+         snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / rdt_ice
+      ENDIF
 
       CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )      
@@ -190,38 +244,8 @@
       CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1.e0 - pfrld(:,:)) )
 
-      !------------------------------------------!
-      !      mass flux at the ocean surface      !
-      !------------------------------------------!
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-#if defined key_coupled
-            ! freshwater exchanges at the ice-atmosphere / ocean interface (coupled mode)
-            zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   ! 
-               &   + rdmsnif(ji,jj) * r1_rdtice                                   !  freshwaterflux due to snow melting 
-#else
-            !  computing freshwater exchanges at the ice/ocean interface
-            zemp = + emp(ji,jj)     *         frld(ji,jj)      &   !  e-p budget over open ocean fraction 
-               &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &   !  liquid precipitation reaches directly the ocean
-               &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  change in ice cover within the time step
-               &   + rdmsnif(ji,jj) * r1_rdtice                    !  freshwater flux due to snow melting 
-#endif            
-            !
-            !  computing salt exchanges at the ice/ocean interface
-            zfons = ( soce_0(ji,jj) - sice_0(ji,jj) ) * ( rdmicif(ji,jj) * r1_rdtice ) 
-            !
-            !  converting the salt flux from ice to a freshwater flux from ocean
-            zfm  = zfons / ( sss_m(ji,jj) + epsi16 )
-            !
-            emps(ji,jj) = zemp + zfm      ! surface ocean concentration/dilution effect (use on SSS evolution)
-            emp (ji,jj) = zemp            ! surface ocean volume flux (use on sea-surface height evolution)
-            !
-         END DO
-      END DO
-
       IF( lk_diaar5 ) THEN       ! AR5 diagnostics
-         CALL iom_put( 'isnwmlt_cea'  ,                 rdmsnif(:,:) * r1_rdtice )
-         CALL iom_put( 'fsal_virt_cea',   soce_0(:,:) * rdmicif(:,:) * r1_rdtice )
-         CALL iom_put( 'fsal_real_cea', - sice_0(:,:) * rdmicif(:,:) * r1_rdtice )
+         CALL iom_put( 'isnwmlt_cea'  ,                 rdm_snw(:,:) * r1_rdtice )
+         CALL iom_put( 'fsal_virt_cea',   soce_0(:,:) * rdm_ice(:,:) * r1_rdtice )
+         CALL iom_put( 'fsal_real_cea', - sice_0(:,:) * rdm_ice(:,:) * r1_rdtice )
       ENDIF
 
@@ -243,5 +267,5 @@
       IF(ln_ctl) THEN            ! control print
          CALL prt_ctl(tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns   , clinfo2=' qns     : ')
-         CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps  , clinfo2=' emps    : ')
+         CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=sfx   , clinfo2=' sfx     : ')
          CALL prt_ctl(tab2d_1=utau  , clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
             &         tab2d_2=vtau  , clinfo2=' vtau    : '        , mask2=vmask )
@@ -439,4 +463,17 @@
          END WHERE
       ENDIF
+      !                                      ! embedded sea ice
+      IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
+         snwice_mass  (:,:) = tms(:,:) * ( rhosn * hsnif(:,:) + rhoic * hicif(:,:)  ) * ( 1.0 - frld(:,:) )
+         snwice_mass_b(:,:) = snwice_mass(:,:)
+      ELSE
+         snwice_mass  (:,:) = 0.e0           ! no mass exchanges
+         snwice_mass_b(:,:) = 0.e0           ! no mass exchanges
+      ENDIF
+      IF( nn_ice_embd == 2 .AND.          &  ! full embedment (case 2) & no restart : 
+         &   .NOT.ln_rstart ) THEN           ! deplete the initial ssh below sea-ice area
+         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
+         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
+      ENDIF
       !
    END SUBROUTINE lim_sbc_init_2
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90	(revision 3605)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90	(revision 3625)
@@ -13,27 +13,28 @@
    !!   'key_lim2' :                                  LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   !!   lim_thd_2      : thermodynamic of sea ice
-   !!   lim_thd_init_2 : initialisation of sea-ice thermodynamic
+   !!   lim_thd_2       : thermodynamic of sea ice
+   !!   lim_thd_init_2  : initialisation of sea-ice thermodynamic
    !!----------------------------------------------------------------------
-   USE phycst          ! physical constants
-   USE dom_oce         ! ocean space and time domain variables
+   USE phycst           ! physical constants
+   USE dom_oce          ! ocean space and time domain variables
    USE domvvl
    USE lbclnk
-   USE in_out_manager  ! I/O manager
+   USE in_out_manager   ! I/O manager
    USE lib_mpp
-   USE wrk_nemo        ! work arrays
-   USE iom             ! IOM library
-   USE ice_2           ! LIM sea-ice variables
-   USE sbc_oce         ! 
-   USE sbc_ice         ! 
-   USE thd_ice_2       ! LIM thermodynamic sea-ice variables
-   USE dom_ice_2       ! LIM sea-ice domain
+   USE wrk_nemo         ! work arrays
+   USE iom              ! IOM library
+   USE ice_2            ! LIM sea-ice variables
+   USE sbc_oce          ! 
+   USE sbc_ice          ! 
+   USE thd_ice_2        ! LIM thermodynamic sea-ice variables
+   USE dom_ice_2        ! LIM sea-ice domain
    USE limthd_zdf_2
    USE limthd_lac_2
    USE limtab_2
-   USE prtctl          ! Print control
-   USE cpl_oasis3, ONLY : lk_cpl
-   USE diaar5, ONLY :   lk_diaar5
-      
+   USE prtctl           ! Print control
+   USE cpl_oasis3, ONLY :   lk_cpl
+   USE diaar5    , ONLY :   lk_diaar5
+   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   
    IMPLICIT NONE
    PRIVATE
@@ -55,5 +56,4 @@
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -89,6 +89,4 @@
       REAL(wp) ::   za , zh, zthsnice    !
       REAL(wp) ::   zfric_u              ! friction velocity 
-      REAL(wp) ::   zfnsol               ! total non solar heat
-      REAL(wp) ::   zfontn               ! heat flux from snow thickness
       REAL(wp) ::   zfntlat, zpareff     ! test. the val. of lead heat budget
 
@@ -129,5 +127,4 @@
       zdvolif(:,:) = 0.e0   ! total variation of ice volume
       zdvonif(:,:) = 0.e0   ! transformation of snow to sea-ice volume
-!      zdvonif(:,:) = 0.e0   ! lateral variation of ice volume
       zlicegr(:,:) = 0.e0   ! lateral variation of ice volume
       zdvomif(:,:) = 0.e0   ! variation of ice volume at bottom due to melting only
@@ -137,6 +134,8 @@
       ffltbif(:,:) = 0.e0   ! linked with fstric
       qfvbq  (:,:) = 0.e0   ! linked with fstric
-      rdmsnif(:,:) = 0.e0   ! variation of snow mass per unit area
-      rdmicif(:,:) = 0.e0   ! variation of ice mass per unit area
+      rdm_snw(:,:) = 0.e0   ! variation of snow mass over 1 time step
+      rdq_snw(:,:) = 0.e0   ! heat content associated with rdm_snw
+      rdm_ice(:,:) = 0.e0   ! variation of ice mass over 1 time step
+      rdq_ice(:,:) = 0.e0   ! heat content associated with rdm_ice
       zmsk (:,:,:) = 0.e0
 
@@ -199,9 +198,7 @@
       !--------------------------------------------------------------------------
 
-      sst_m(:,:) = sst_m(:,:) + rt0
-
-!CDIR NOVERRCHK
-      DO jj = 1, jpj
-!CDIR NOVERRCHK
+      !CDIR NOVERRCHK
+      DO jj = 1, jpj
+         !CDIR NOVERRCHK
          DO ji = 1, jpi
             zthsnice       = hsnif(ji,jj) + hicif(ji,jj)
@@ -217,24 +214,18 @@
             !  temperature and turbulent mixing (McPhee, 1992)
             zfric_u        = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin )  ! friction velocity
-            fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( sst_m(ji,jj) - tfu(ji,jj) ) 
+            fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( sst_m(ji,jj) + rt0 - tfu(ji,jj) ) 
             qdtcn(ji,jj)  = zindb * fdtcn(ji,jj) * frld(ji,jj) * rdt_ice
                         
             !  partial computation of the lead energy budget (qldif)
 #if defined key_coupled 
-            qldif(ji,jj)   = tms(ji,jj) * rdt_ice                                             &
+            qldif(ji,jj)   = tms(ji,jj) * rdt_ice                                                  &
                &    * (   ( qsr_tot(ji,jj) - qsr_ice(ji,jj,1) * zfricp ) * ( 1.0 - thcm(ji,jj) )   &
                &        + ( qns_tot(ji,jj) - qns_ice(ji,jj,1) * zfricp )                           &
                &        + frld(ji,jj) * ( fdtcn(ji,jj) + ( 1.0 - zindb ) * fsbbq(ji,jj) )   )
 #else
-            zfontn         = ( sprecip(ji,jj) / rhosn ) * xlsn  !   energy for melting solid precipitation
-            zfnsol         = qns(ji,jj)                         !  total non solar flux over the ocean
-            qldif(ji,jj)   = tms(ji,jj) * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )   &
-               &                               + zfnsol + fdtcn(ji,jj) - zfontn     &
-               &                               + ( 1.0 - zindb ) * fsbbq(ji,jj) )   &
-               &                        * frld(ji,jj) * rdt_ice    
-!!$            qldif(ji,jj)   = tms(ji,jj) * rdt_ice * frld(ji,jj) 
-!!$               &           * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )      &
-!!$               &             + qns(ji,jj)  + fdtcn(ji,jj) - zfontn     &
-!!$               &             + ( 1.0 - zindb ) * fsbbq(ji,jj)      )   &
+            qldif(ji,jj)   = tms(ji,jj) * rdt_ice * frld(ji,jj)                    &
+               &                        * (  qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )    &
+               &                           + qns(ji,jj)  +  fdtcn(ji,jj)           &
+               &                           + ( 1.0 - zindb ) * fsbbq(ji,jj)      )
 #endif
             !  parlat : percentage of energy used for lateral ablation (0.0) 
@@ -246,6 +237,5 @@
             
             !  energy needed to bring ocean surface layer until its freezing
-            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1)   &
-                &          * ( tfu(ji,jj) - sst_m(ji,jj) ) * ( 1 - zinda )
+            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1) * ( tfu(ji,jj) - sst_m(ji,jj) - rt0 ) * ( 1 - zinda )
             
             !  calculate oceanic heat flux.
@@ -257,6 +247,4 @@
       END DO
       
-      sst_m(:,:) = sst_m(:,:) - rt0
-               
       !         Select icy points and fulfill arrays for the vectorial grid.
       !----------------------------------------------------------------------
@@ -312,6 +300,9 @@
          CALL tab_2d_1d_2( nbpb, qldif_1d   (1:nbpb)     , qldif      , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d_2( nbpb, qstbif_1d  (1:nbpb)     , qstoif     , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d_2( nbpb, rdmicif_1d (1:nbpb)     , rdmicif    , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d_2( nbpb, rdm_ice_1d (1:nbpb)     , rdm_ice    , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d_2( nbpb, rdq_ice_1d (1:nbpb)     , rdq_ice    , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d_2( nbpb, dmgwi_1d   (1:nbpb)     , dmgwi      , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d_2( nbpb, rdm_snw_1d (1:nbpb)     , rdm_snw    , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d_2( nbpb, rdq_snw_1d (1:nbpb)     , rdq_snw    , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d_2( nbpb, qlbbq_1d   (1:nbpb)     , zqlbsbq    , jpi, jpj, npb(1:nbpb) )
          !
@@ -332,7 +323,9 @@
          CALL tab_1d_2d_2( nbpb, qfvbq      , npb, qfvbq_1d  (1:nbpb)     , jpi, jpj )
          CALL tab_1d_2d_2( nbpb, qstoif     , npb, qstbif_1d (1:nbpb)     , jpi, jpj )
-         CALL tab_1d_2d_2( nbpb, rdmicif    , npb, rdmicif_1d(1:nbpb)     , jpi, jpj )
+         CALL tab_1d_2d_2( nbpb, rdm_ice    , npb, rdm_ice_1d(1:nbpb)     , jpi, jpj )
+         CALL tab_1d_2d_2( nbpb, rdq_ice    , npb, rdq_ice_1d(1:nbpb)     , jpi, jpj )
          CALL tab_1d_2d_2( nbpb, dmgwi      , npb, dmgwi_1d  (1:nbpb)     , jpi, jpj )
-         CALL tab_1d_2d_2( nbpb, rdmsnif    , npb, rdmsnif_1d(1:nbpb)     , jpi, jpj )
+         CALL tab_1d_2d_2( nbpb, rdm_snw    , npb, rdm_snw_1d(1:nbpb)     , jpi, jpj )
+         CALL tab_1d_2d_2( nbpb, rdq_snw    , npb, rdq_snw_1d(1:nbpb)     , jpi, jpj )
          CALL tab_1d_2d_2( nbpb, zdvosif    , npb, dvsbq_1d  (1:nbpb)     , jpi, jpj )
          CALL tab_1d_2d_2( nbpb, zdvobif    , npb, dvbbq_1d  (1:nbpb)     , jpi, jpj )
@@ -393,5 +386,5 @@
       IF( nbpac > 0 ) THEN
          !
-         zlicegr(:,:) = rdmicif(:,:)		! to output the lateral sea-ice growth 
+         zlicegr(:,:) = rdm_ice(:,:)		! to output the lateral sea-ice growth 
          !...Put the variable in a 1-D array for lateral accretion
          CALL tab_2d_1d_2( nbpac, frld_1d   (1:nbpac)     , frld       , jpi, jpj, npac(1:nbpac) )
@@ -404,5 +397,6 @@
          CALL tab_2d_1d_2( nbpac, qcmif_1d  (1:nbpac)     , qcmif      , jpi, jpj, npac(1:nbpac) )
          CALL tab_2d_1d_2( nbpac, qstbif_1d (1:nbpac)     , qstoif     , jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d_2( nbpac, rdmicif_1d(1:nbpac)     , rdmicif    , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d_2( nbpac, rdm_ice_1d(1:nbpac)     , rdm_ice    , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d_2( nbpac, rdq_ice_1d(1:nbpac)     , rdq_ice    , jpi, jpj, npac(1:nbpac) )
          CALL tab_2d_1d_2( nbpac, dvlbq_1d  (1:nbpac)     , zdvolif    , jpi, jpj, npac(1:nbpac) )
          CALL tab_2d_1d_2( nbpac, tfu_1d    (1:nbpac)     , tfu        , jpi, jpj, npac(1:nbpac) )
@@ -418,5 +412,6 @@
          CALL tab_1d_2d_2( nbpac, tbif(:,:,3), npac(1:nbpac), tbif_1d   (1:nbpac , 3 ), jpi, jpj )
          CALL tab_1d_2d_2( nbpac, qstoif     , npac(1:nbpac), qstbif_1d (1:nbpac)     , jpi, jpj )
-         CALL tab_1d_2d_2( nbpac, rdmicif    , npac(1:nbpac), rdmicif_1d(1:nbpac)     , jpi, jpj )
+         CALL tab_1d_2d_2( nbpac, rdm_ice    , npac(1:nbpac), rdm_ice_1d(1:nbpac)     , jpi, jpj )
+         CALL tab_1d_2d_2( nbpac, rdq_ice    , npac(1:nbpac), rdq_ice_1d(1:nbpac)     , jpi, jpj )
          CALL tab_1d_2d_2( nbpac, zdvolif    , npac(1:nbpac), dvlbq_1d  (1:nbpac)     , jpi, jpj )
          !
@@ -449,11 +444,11 @@
       CALL iom_put( 'iceprod_cea' , hicifp (:,:) * zztmp     )   ! Ice produced               [m/s]
       IF( lk_diaar5 ) THEN
-         CALL iom_put( 'snowmel_cea' , rdmsnif(:,:) * zztmp     )   ! Snow melt                  [kg/m2/s]
+         CALL iom_put( 'snowmel_cea' , rdm_snw(:,:) * zztmp     )   ! Snow melt                  [kg/m2/s]
          zztmp = rhoic / rdt_ice
          CALL iom_put( 'sntoice_cea' , zdvonif(:,:) * zztmp     )   ! Snow to Ice transformation [kg/m2/s]
          CALL iom_put( 'ticemel_cea' , zdvosif(:,:) * zztmp     )   ! Melt at Sea Ice top        [kg/m2/s]
          CALL iom_put( 'bicemel_cea' , zdvomif(:,:) * zztmp     )   ! Melt at Sea Ice bottom     [kg/m2/s]
-         zlicegr(:,:) = MAX( 0.e0, rdmicif(:,:)-zlicegr(:,:) )
-         CALL iom_put( 'licepro_cea' , zlicegr(:,:) * zztmp     )   ! Latereal sea ice growth    [kg/m2/s]
+         zlicegr(:,:) = MAX( 0.e0, rdm_ice(:,:)-zlicegr(:,:) )
+         CALL iom_put( 'licepro_cea' , zlicegr(:,:) * zztmp     )   ! Lateral sea ice growth     [kg/m2/s]
       ENDIF
       !
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limthd_lac_2.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limthd_lac_2.F90	(revision 3605)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limthd_lac_2.F90	(revision 3625)
@@ -7,13 +7,14 @@
 
    !!----------------------------------------------------------------------
-   !!   lim_lat_acr_2   : lateral accretion of ice
-   !!----------------------------------------------------------------------
-   USE par_oce          ! ocean parameters
+   !!   lim_lat_acr_2 : lateral accretion of ice
+   !!----------------------------------------------------------------------
+   USE par_oce        ! ocean parameters
    USE phycst
    USE thd_ice_2
    USE ice_2
    USE limistate_2 
-   USE lib_mpp          ! MPP library
-   USE wrk_nemo         ! work arrays
+   USE lib_mpp        ! MPP library
+   USE wrk_nemo       ! work arrays
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 
    IMPLICIT NONE
@@ -145,6 +146,6 @@
          frld_1d   (ji) = MAX( zfrlnew , zfrlmin(ji) )
          !--computation of the remaining part of ice thickness which has been already used
-         zdhicbot(ji) =  ( frld_1d(ji) - zfrlnew ) * zhice0(ji) / ( 1.0 - zfrlmin(ji) ) & 
-                      -  (  ( 1.0 - zfrrate ) / ( 1.0 - frld_1d(ji) ) )  * ( zqbgow(ji) / xlic ) 
+         zdhicbot(ji) =  ( frld_1d(ji) - zfrlnew ) * zhice0(ji) / ( 1.0 - zfrlmin(ji) )   & 
+            &         -  (  ( 1.0 - zfrrate ) / ( 1.0 - frld_1d(ji) ) )  * ( zqbgow(ji) / xlic ) 
       END DO
  
@@ -196,5 +197,5 @@
             &          ) / zah
          
-         tbif_1d(ji,3) =     (  iiceform * ( zhnews2 - zdh3 )                                          * zta1  &
+         tbif_1d(ji,3) =     ( iiceform * ( zhnews2 - zdh3 )                                           * zta1  &
             &              + ( iiceform * zdh3 + ( 1 - iiceform ) * zdh1 )                             * zta2  &
             &              + ( iiceform * ( zhnews2 - zdh5 ) + ( 1 - iiceform ) * ( zhnews2 - zdh1 ) ) * zta3  & 
@@ -217,5 +218,6 @@
       DO ji = kideb , kiut
          dvlbq_1d  (ji) = ( 1. - frld_1d(ji) ) * h_ice_1d(ji) - ( 1. - zfrl_old(ji) ) * zhice_old(ji)
-         rdmicif_1d(ji) = rdmicif_1d(ji) + rhoic * dvlbq_1d(ji)
+         rdm_ice_1d(ji) = rdm_ice_1d(ji) + rhoic * dvlbq_1d(ji)
+         rdq_ice_1d(ji) = rdq_ice_1d(ji) + rcpic * dvlbq_1d(ji) * ( tfu_1d(ji) - rt0 )      ! heat content relative to rt0
       END DO
       
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limthd_zdf_2.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limthd_zdf_2.F90	(revision 3605)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limthd_zdf_2.F90	(revision 3625)
@@ -18,9 +18,10 @@
    USE ice_2
    USE limistate_2
+   USE cpl_oasis3, ONLY : lk_cpl
    USE in_out_manager
    USE lib_mpp          ! MPP library
    USE wrk_nemo         ! work arrays
-   USE cpl_oasis3, ONLY : lk_cpl
-      
+   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+    
    IMPLICIT NONE
    PRIVATE
@@ -86,5 +87,5 @@
       REAL(wp), POINTER, DIMENSION(:) ::   zrcpdt         ! h_su*rho_su*cp_su/dt(h_su being the thick. of surf. layer)
       REAL(wp), POINTER, DIMENSION(:) ::   zts_old        ! previous surface temperature
-      REAL(wp), POINTER, DIMENSION(:) ::   zidsn , z1midsn , zidsnic ! tempory variables
+      REAL(wp), POINTER, DIMENSION(:) ::   zidsn , z1midsn , zidsnic ! temporary variables
       REAL(wp), POINTER, DIMENSION(:) ::   zfnet          ! net heat flux at the top surface( incl. conductive heat flux)
       REAL(wp), POINTER, DIMENSION(:) ::   zsprecip       ! snow accumulation
@@ -98,62 +99,68 @@
       REAL(wp), POINTER, DIMENSION(:) ::   zep            ! internal temperature of the 2nd layer of the snow/ice system
       REAL(wp), DIMENSION(3) :: & 
-          zplediag  &    ! principle diagonal, subdiag. and supdiag. of the 
+            zplediag  &    ! principle diagonal, subdiag. and supdiag. of the 
           , zsubdiag  &    ! tri-diagonal matrix coming from the computation
           , zsupdiag  &    ! of the temperatures inside the snow-ice system
           , zsmbr          ! second member
-       REAL(wp) :: & 
-          zhsu     &     ! thickness of surface layer
-          , zhe      &     ! effective thickness for compu. of equ. thermal conductivity
-          , zheshth  &     ! = zhe / thth
-          , zghe     &     ! correction factor of the thermal conductivity
-          , zumsb    &     ! parameter for numerical method to solve heat-diffusion eq.
-          , zkhsn    &     ! conductivity at the snow layer
-          , zkhic    &     ! conductivity at the ice layers
-          , zkint    &     ! equivalent conductivity at the snow-ice interface
-          , zkhsnint &     ! = zkint*dt / (hsn*rhosn*cpsn)  
-          , zkhicint &     ! = 2*zkint*dt / (hic*rhoic*cpic)
-          , zpiv1 , zpiv2  &       ! tempory scalars used to solve the tri-diagonal system
-          , zb2 , zd2 , zb3 , zd3 &
+       REAL(wp) ::    &
+            zhsu      &    ! thickness of surface layer
+          , zhe       &    ! effective thickness for compu. of equ. thermal conductivity
+          , zheshth   &    ! = zhe / thth
+          , zghe      &    ! correction factor of the thermal conductivity
+          , zumsb     &    ! parameter for numerical method to solve heat-diffusion eq.
+          , zkhsn     &    ! conductivity at the snow layer
+          , zkhic     &    ! conductivity at the ice layers
+          , zkint     &    ! equivalent conductivity at the snow-ice interface
+          , zkhsnint  &    ! = zkint*dt / (hsn*rhosn*cpsn)  
+          , zkhicint  &    ! = 2*zkint*dt / (hic*rhoic*cpic)
+          , zpiv1, zpiv2 & ! temporary scalars used to solve the tri-diagonal system
+          , zb2, zd2  &    ! temporary scalars used to solve the tri-diagonal system
+          , zb3, zd3  &    ! temporary scalars used to solve the tri-diagonal system
           , ztint          ! equivalent temperature at the snow-ice interface
-       REAL(wp) :: & 
-          zexp      &     ! exponential function of the ice thickness
-          , zfsab     &     ! part of solar radiation stored in brine pockets
-          , zfts      &     ! value of energy balance function when the temp. equal surf. temp.
-          , zdfts     &     ! value of derivative of ztfs when the temp. equal surf. temp.
-          , zdts      &     ! surface temperature increment
-          , zqsnw_mlt &     ! energy needed to melt snow
-          , zdhsmlt   &     ! change in snow thickness due to melt
-          , zhsn      &     ! snow thickness (previous+accumulation-melt)
-          , zqsn_mlt_rem &  ! remaining heat coming from snow melting
-          , zqice_top_mlt & ! energy used to melt ice at top surface
-          , zdhssub      &  ! change in snow thick. due to sublimation or evaporation
-          , zdhisub      &  ! change in ice thick. due to sublimation or evaporation    
-          , zdhsn        &  ! snow ice thickness increment
-          , zdtsn        &  ! snow internal temp. increment
-          , zdtic        &  ! ice internal temp. increment
+       REAL(wp) ::    & 
+            zexp      &    ! exponential function of the ice thickness
+          , zfsab     &    ! part of solar radiation stored in brine pockets
+          , zfts      &    ! value of energy balance function when the temp. equal surf. temp.
+          , zdfts     &    ! value of derivative of ztfs when the temp. equal surf. temp.
+          , zdts      &    ! surface temperature increment
+          , zqsnw_mlt &    ! energy needed to melt snow
+          , zdhsmlt   &    ! change in snow thickness due to melt
+          , zhsn      &    ! snow thickness (previous+accumulation-melt)
+          , zqsn_mlt_rem & ! remaining heat coming from snow melting
+          , zqice_top_mlt &! energy used to melt ice at top surface
+          , zdhssub   &    ! change in snow thick. due to sublimation or evaporation
+          , zdhisub   &    ! change in ice thick. due to sublimation or evaporation    
+          , zdhsn     &    ! snow ice thickness increment
+          , zdtsn     &    ! snow internal temp. increment
+          , zdtic     &    ! ice internal temp. increment
           , zqnes          ! conductive energy due to ice melting in the first ice layer
-       REAL(wp) :: & 
-          ztbot     &      ! temperature at the bottom surface
-          , zfcbot    &      ! conductive heat flux at bottom surface
-          , zqice_bot &      ! energy used for bottom melting/growing
-          , zqice_bot_mlt &  ! energy used for bottom melting
-          , zqstbif_bot  &  ! part of energy stored in brine pockets used for bottom melting
-          , zqstbif_old  &  ! tempory var. for zqstbif_bot
-          , zdhicmlt      &  ! change in ice thickness due to bottom melting
-          , zdhicm        &  ! change in ice thickness var. 
-          , zdhsnm        &  ! change in snow thickness var. 
-          , zhsnfi        &  ! snow thickness var. 
-          , zc1, zpc1, zc2, zpc2, zp1, zp2 & ! tempory variables
-          , ztb2, ztb3
-       REAL(wp) :: & 
-          zdrmh         &   ! change in snow/ice thick. after snow-ice formation
-          , zhicnew       &   ! new ice thickness
-          , zhsnnew       &   ! new snow thickness
-          , zquot , ztneq &   ! tempory temp. variables
-          , zqice, zqicetot & ! total heat inside the snow/ice system
-          , zdfrl         &   ! change in ice concentration
-          , zdvsnvol      &   ! change in snow volume
-          , zdrfrl1, zdrfrl2 &  ! tempory scalars
-          , zihsn, zidhb, zihic, zihe, zihq, ziexp, ziqf, zihnf, zibmlt, ziqr, zihgnew, zind
+       REAL(wp) ::    & 
+            ztbot     &    ! temperature at the bottom surface
+          , zfcbot    &    ! conductive heat flux at bottom surface
+          , zqice_bot &    ! energy used for bottom melting/growing
+          , zqice_bot_mlt &! energy used for bottom melting
+          , zqstbif_bot  & ! part of energy stored in brine pockets used for bottom melting
+          , zqstbif_old  & ! temporary var. for zqstbif_bot
+          , zdhicmlt  &    ! change in ice thickness due to bottom melting
+          , zdhicm    &    ! change in ice thickness var. 
+          , zdhsnm    &    ! change in snow thickness var. 
+          , zhsnfi    &    ! snow thickness var. 
+          , zc1, zpc1 &    ! temporary variables
+          , zc2, zpc2 &    ! temporary variables
+          , zp1, zp2  &    ! temporary variables
+          , ztb2, ztb3     ! temporary variables
+       REAL(wp) ::    & 
+            zdrmh     &    ! change in snow/ice thick. after snow-ice formation
+          , zhicnew   &    ! new ice thickness
+          , zhsnnew   &    ! new snow thickness
+          , zquot     &
+          , ztneq     &    ! temporary temp. variables
+          , zqice     &
+          , zqicetot  &    ! total heat inside the snow/ice system
+          , zdfrl     &    ! change in ice concentration
+          , zdvsnvol  &    ! change in snow volume
+          , zdrfrl1, zdrfrl2, zihsn, zidhb, zihic &  ! temporary scalars
+          , zihe, zihq, ziexp, ziqf, zihnf        &  ! temporary scalars
+          , zibmlt, ziqr, zihgnew, zind, ztmp        ! temporary scalars
        !!----------------------------------------------------------------------
        CALL wrk_alloc( jpij, ztsmlt, ztbif  , zksn    , zkic    , zksndh , zfcsu  , zfcsudt , zi0      , z1mi0   , zqmax    )
@@ -169,16 +176,14 @@
        
        DO ji = kideb , kiut
+          ! do nothing if the snow (ice) thickness falls below its minimum thickness
           zihsn = MAX( zzero , SIGN( zone , hsndif - h_snow_1d(ji) ) )
           zihic = MAX( zzero , SIGN( zone , hicdif - h_ice_1d(ji) ) )
-          !--computation of energy due to surface melting
-          zqcmlts(ji) = ( MAX ( zzero ,  &
-             &                   rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn )
-          !--computation of energy due to bottom melting
-          zqcmltb(ji) = ( MAX( zzero , &
-             &                  rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) &
-             &           + MAX( zzero , &
-             &                  rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) &
-             &           ) * ( 1.0 - zihic  )
-          !--limitation of  snow/ice system internal temperature
+          !--energy required to bring snow to its melting point (rt0_snow)
+          zqcmlts(ji) = ( MAX ( zzero , rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn )
+          !--energy required to bring ice to its melting point (rt0_ice)
+          zqcmltb(ji) = ( MAX( zzero , rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) )  &
+             &          + MAX( zzero , rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) )  &
+             &          ) * ( 1.0 - zihic  )
+          !--limitation of snow/ice system internal temperature
           tbif_1d(ji,1)   = MIN( rt0_snow, tbif_1d(ji,1) )
           tbif_1d(ji,2)   = MIN( rt0_ice , tbif_1d(ji,2) )
@@ -480,5 +485,8 @@
           dvsbq_1d(ji) =  ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnw_old(ji) - zsprecip(ji) )
           dvsbq_1d(ji) =  MIN( zzero , dvsbq_1d(ji) )
-          rdmsnif_1d(ji) =  rhosn * dvsbq_1d(ji)
+          ztmp = rhosn * dvsbq_1d(ji)
+          rdm_snw_1d(ji) =  ztmp
+          !--heat content of the water provided to the ocean (referenced to rt0)
+          rdq_snw_1d(ji) =  cpic * ztmp * ( rt0_snow - rt0 )
           !-- If the snow is completely melted the remaining heat is used to melt ice
           zqsn_mlt_rem  = MAX( zzero , -zhsn ) * xlsn
@@ -623,9 +631,17 @@
           !---updating new ice thickness and computing the newly formed ice mass
           zhicnew   =  zihgnew * zhicnew
-          rdmicif_1d(ji) =  rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) ) * rhoic
+          ztmp    =  ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) ) * rhoic
+          rdm_ice_1d(ji) =  rdm_ice_1d(ji) + ztmp
+          !---heat content of the water provided to the ocean (referenced to rt0)
+          !   use of rt0_ice is OK for melting ice; in the case of freezing, tfu_1d should be used. 
+          !   This is done in 9.5 section (see below)
+          rdq_ice_1d(ji) =  cpic * ztmp * ( rt0_ice - rt0 )
           !---updating new snow thickness and computing the newly formed snow mass
           zhsnfi   = zhsn + zdhsnm
           h_snow_1d(ji) = MAX( zzero , zhsnfi )
-          rdmsnif_1d(ji) =  rdmsnif_1d(ji) + ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsn ) * rhosn
+          ztmp = ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsn ) * rhosn
+          rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp
+          !---updating the heat content of the water provided to the ocean (referenced to rt0)
+          rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( rt0_snow - rt0 )
           !--remaining energy in case of total ablation
           zqocea(ji) = - ( zihsn * xlic * zdhicm + xlsn * ( zhsnfi - h_snow_1d(ji) ) ) * ( 1.0 - frld_1d(ji) )
@@ -659,4 +675,8 @@
           tbif_1d(ji,3) =  zihgnew * ztb3 + ( 1.0 - zihgnew ) * tfu_1d(ji)
           h_ice_1d(ji)  =  zhicnew
+          ! update the ice heat content given to the ocean in freezing case 
+          ! (part due to difference between rt0_ice and tfu_1d)
+          ztmp = ( 1. - zidhb ) * rhoic * dvbbq_1d(ji)
+          rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * ( tfu_1d(ji) - rt0_ice )
        END DO
 
@@ -700,6 +720,11 @@
           dmgwi_1d(ji) = dmgwi_1d(ji) + ( 1.0 -frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnnew ) * rhosn
           !---  volume change of ice and snow (used for ocean-ice freshwater flux computation)
-          rdmicif_1d(ji) = rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhicnew - h_ice_1d (ji) ) * rhoic
-          rdmsnif_1d(ji) = rdmsnif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhsnnew - h_snow_1d(ji) ) * rhosn
+          ztmp = ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d (ji) ) * rhoic
+          rdm_ice_1d(ji) = rdm_ice_1d(ji) + ztmp
+          rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * ( tfu_1d(ji) - rt0 )
+          !!gm BUG ??   snow ==>  only needed for nn_ice_embd == 0  (standard levitating sea-ice)
+          ztmp = ( 1.0 - frld_1d(ji) )   * ( zhsnnew - h_snow_1d(ji) ) * rhosn
+          rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp
+          rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( rt0_snow - rt0 )
 
           !---  Actualize new snow and ice thickness.
@@ -748,8 +773,28 @@
           !--variation of ice volume and ice mass 
           dvlbq_1d(ji)   = zihic * ( zfrl_old(ji) - frld_1d(ji) ) * h_ice_1d(ji)
-          rdmicif_1d(ji) = rdmicif_1d(ji) + dvlbq_1d(ji) * rhoic
+          ztmp = dvlbq_1d(ji) * rhoic
+          rdm_ice_1d(ji) = rdm_ice_1d(ji) + ztmp
+!!gm
+!!gm   This should be split in two parts:
+!!gm         1-  heat required to bring sea-ice to tfu  : this part should be added to the heat flux taken from the ocean
+!!gm                 cpic * ztmp * 0.5 * ( tbif_1d(ji,2) + tbif_1d(ji,3) - 2.* rt0_ice )
+!!gm         2-  heat content of lateral ablation referenced to rt0 : this part only put in rdq_ice_1d
+!!gm                 cpic * ztmp * ( rt0_ice - rt0 )
+!!gm   Currently we put all the heat in rdq_ice_1d
+          rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * 0.5 * ( tbif_1d(ji,2) + tbif_1d(ji,3) - 2.* rt0 )
+          !
           !--variation of snow volume and snow mass 
-          zdvsnvol    = zihsn * ( zfrl_old(ji) - frld_1d(ji) ) * h_snow_1d(ji)
-          rdmsnif_1d(ji) = rdmsnif_1d(ji) + zdvsnvol * rhosn
+          zdvsnvol = zihsn * ( zfrl_old(ji) - frld_1d(ji) ) * h_snow_1d(ji)
+          ztmp     = zdvsnvol * rhosn
+          rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp
+!!gm
+!!gm   This should be split in two parts:
+!!gm         1-  heat required to bring snow to tfu  : this part should be added to the heat flux taken from the ocean
+!!gm                 cpic * ztmp * ( tbif_1d(ji,1) - rt0_snow )
+!!gm         2-  heat content of lateral ablation referenced to rt0 : this part only put in rdq_snw_1d
+!!gm                 cpic * ztmp * ( rt0_snow - rt0 )
+!!gm   Currently we put all the heat in rdq_snw_1d
+          rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( tbif_1d(ji,1) - rt0 )
+
           h_snow_1d(ji)  = ziqf * h_snow_1d(ji)
 
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limwri_2.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limwri_2.F90	(revision 3605)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limwri_2.F90	(revision 3625)
@@ -13,6 +13,6 @@
    !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
-   !!   lim_wri_2      : write of the diagnostics variables in ouput file 
-   !!   lim_wri_init_2 : initialization and namelist read
+   !!   lim_wri_2       : write of the diagnostics variables in ouput file 
+   !!   lim_wri_init_2  : initialization and namelist read
    !!   lim_wri_state_2 : write for initial state or/and abandon:
    !!                     > output.init.nc (if ninist = 1 in namelist)
@@ -26,11 +26,12 @@
    USE ice_2
 
-   USE dianam          ! build name of file (routine)
+   USE dianam           ! build name of file (routine)
    USE lbclnk
    USE in_out_manager
-   USE lib_mpp         ! MPP library
-   USE wrk_nemo        ! work arrays
+   USE lib_mpp          ! MPP library
+   USE wrk_nemo         ! work arrays
    USE iom
    USE ioipsl
+   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 
    IMPLICIT NONE
@@ -173,5 +174,5 @@
             zcmo(ji,jj,13) = qns(ji,jj)
             ! See thersf for the coefficient
-            zcmo(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce    !!gm ???
+            zcmo(ji,jj,14) = - sfx (ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce    !!gm ???
             zcmo(ji,jj,15) = utau_ice(ji,jj)
             zcmo(ji,jj,16) = vtau_ice(ji,jj)
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limwri_dimg_2.h90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limwri_dimg_2.h90	(revision 3605)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limwri_dimg_2.h90	(revision 3625)
@@ -118,5 +118,5 @@
           zcmo(ji,jj,13) = qns(ji,jj)
           ! See thersf for the coefficient
-          zcmo(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce
+          zcmo(ji,jj,14) = - sfx (ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce
           zcmo(ji,jj,15) = utau_ice(ji,jj)
           zcmo(ji,jj,16) = vtau_ice(ji,jj)
@@ -161,5 +161,5 @@
                 rcmoy(ji,jj,13) = qns(ji,jj)
                 ! See thersf for the coefficient
-                rcmoy(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce
+                rcmoy(ji,jj,14) = - sfx (ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce
                 rcmoy(ji,jj,15) = utau_ice(ji,jj)
                 rcmoy(ji,jj,16) = vtau_ice(ji,jj)
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/thd_ice_2.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/thd_ice_2.F90	(revision 3605)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/thd_ice_2.F90	(revision 3625)
@@ -68,6 +68,8 @@
       qstbif_1d   ,     &  !:    "                  "      qstoif
       fbif_1d     ,     &  !:    "                  "      fbif
-      rdmicif_1d  ,     &  !:    "                  "      rdmicif
-      rdmsnif_1d  ,     &  !:    "                  "      rdmsnif
+      rdm_ice_1d  ,     &  !:    "                  "      rdm_ice
+      rdq_ice_1d  ,     &  !:    "                  "      rdq_ice
+      rdm_snw_1d  ,     &  !:    "                  "      rdm_snw
+      rdq_snw_1d  ,     &  !:    "                  "      rdq_snw
       qlbbq_1d    ,     &  !:    "                  "      qlbsbq
       dmgwi_1d    ,     &  !:    "                  "      dmgwi
@@ -108,5 +110,6 @@
          &      qstbif_1d(jpij),  fbif_1d(jpij),  Stat=ierr(2))
          !
-      ALLOCATE( rdmicif_1d(jpij), rdmsnif_1d(jpij), qlbbq_1d(jpij),   &
+      ALLOCATE( rdm_ice_1d(jpij), rdq_ice_1d(jpij)                  , &
+         &      rdm_snw_1d(jpij), rdq_snw_1d(jpij), qlbbq_1d(jpij)  , &
          &      dmgwi_1d(jpij)  , dvsbq_1d(jpij)  , rdvomif_1d(jpij), &
          &      dvbbq_1d(jpij)  , dvlbq_1d(jpij)  , dvnbq_1d(jpij)  , &
