Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 3402)
@@ -39,4 +39,5 @@
    USE prtctl           ! Print control
    USE cpl_oasis3, ONLY : lk_cpl
+   USE oce, ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass 
 
    IMPLICIT NONE
@@ -92,5 +93,5 @@
       !!              - qns     : sea heat flux    : non solar (including heat content of the mass flux)
       !!              - emp     : freshwater budget: mass flux 
-      !!              - emps    : freshwater budget: salt flux due to Freezing/Melting
+      !!              - sfx     : freshwater budget: salt flux due to Freezing/Melting
       !!              - utau    : sea surface i-stress (ocean referential)
       !!              - vtau    : sea surface j-stress (ocean referential)
@@ -119,11 +120,10 @@
       CALL wrk_alloc( jpi, jpj, 1, zalb, zalbp )
 
-      zswitch = 1                                ! Default standard levitating sea-ice (salt exchanges only)
-!!gm ice embedment
-!     SELECT CASE( nn_ice_embd )                 ! levitating/embedded sea-ice option (not yet activated)
-!       CASE( 0    )   ;   zswitch = 1           ! standard levitating sea-ice : salt exchange only
-!       CASE( 1, 2 )   ;   zswitch = 0           ! other levitating sea-ice or embedded sea-ice : salt and volume fluxes
-!     END SELECT                                 !    
-!!gm end embedment
+      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                                 !    
+
       !------------------------------------------!
       !      heat flux at the ocean surface      !
@@ -147,5 +147,5 @@
 !!$            i1mfr   = 1.0 - AINT(  frld(ji,jj) )                   ! = 0. if ice-free ocean else 1. (after ice thermo)
 !!$
-!!$            IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda      ! = 1. if there was snow and ice before the ice thermo. which has been completely melted (possibly overmelted)
+!!$            IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda      ! = zinda if previous thermodynamic step overmelted the ice???
 !!$            ELSE                             ;   ifvt = 0.         ! 
 !!$            ENDIF
@@ -222,9 +222,19 @@
             ! salt flux only       : add concentration dilution term in salt flux  and no  F-M term in volume flux
             ! salt and mass fluxes : non concentartion dilution term in salt flux  and add F-M term in volume flux
-            emps(ji,jj) = zfsalt +                  zswitch  * zcd   ! salt flux (+ C/D if no ice/ocean mass exchange)
+            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 no 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(:,:) )      
@@ -255,5 +265,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 )
@@ -451,4 +461,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 belew 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_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limwri_2.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limwri_2.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limwri_2.F90	(revision 3402)
@@ -173,5 +173,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_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limwri_dimg_2.h90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limwri_dimg_2.h90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limwri_dimg_2.h90	(revision 3402)
@@ -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_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limdia.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limdia.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limdia.F90	(revision 3402)
@@ -113,5 +113,5 @@
                vinfor(31) = vinfor(31) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + & 
                   v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12 
-               vinfor(53) = vinfor(53) + emps(ji,jj)*aire(ji,jj) * 1.e-12_wp !salt flux
+               vinfor(53) = vinfor(53) + sfx (ji,jj)*aire(ji,jj) * 1.e-12_wp !salt flux
                vinfor(55) = vinfor(55) + fsbri(ji,jj)*aire(ji,jj) * 1.e-12_wp !brine drainage flux
                vinfor(57) = vinfor(57) + fseqv(ji,jj)*aire(ji,jj) * 1.e-12_wp !equivalent salt flux
@@ -264,5 +264,5 @@
                vinfor(32) = vinfor(32) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + & 
                   v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12 !ice vel
-               vinfor(54) = vinfor(54) + at_i(ji,jj)*emps(ji,jj)*aire(ji,jj) * 1.e-12_wp ! Total salt flux
+               vinfor(54) = vinfor(54) + at_i(ji,jj)*sfx (ji,jj)*aire(ji,jj) * 1.e-12_wp ! Total salt flux
                vinfor(56) = vinfor(56) + at_i(ji,jj)*fsbri(ji,jj)*aire(ji,jj) * 1.e-12_wp ! Brine drainage salt flux
                vinfor(58) = vinfor(58) + at_i(ji,jj)*fseqv(ji,jj)*aire(ji,jj) * 1.e-12_wp ! Equivalent salt flux
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 3402)
@@ -1545,5 +1545,5 @@
                !                               rhoic * v_i(ji,jj,jl) / rdt_ice
 
-               !           emps(i,j)      = emps(i,j)      + xtmp
+               !           sfx (i,j)      = sfx (i,j)      + xtmp
                !           fsalt_hist(i,j) = fsalt_hist(i,j) + xtmp
 
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 3402)
@@ -86,5 +86,5 @@
       !!              - qns     : sea heat flux: non solar
       !!              - emp     : freshwater budget: volume flux 
-      !!              - emps    : freshwater budget: concentration/dillution 
+      !!              - sfx     : freshwater budget: concentration/dillution 
       !!              - fr_i    : ice fraction
       !!              - tn_ice  : sea-ice surface temperature
@@ -263,7 +263,7 @@
 
       IF( num_sal == 2 ) THEN      ! variable ice salinity: brine drainage included in the salt flux
-         emps(:,:) = fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)
+         sfx (:,:) = fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)
       ELSE                         ! constant ice salinity:
-         emps(:,:) =              fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)
+         sfx (:,:) =              fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)
       ENDIF
 
@@ -285,5 +285,5 @@
       IF(ln_ctl) THEN
          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=fr_i  , clinfo1=' lim_sbc: fr_i   : ' )
          CALL prt_ctl( tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl )
@@ -402,4 +402,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 belew sea-ice area
+         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
+         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
+      ENDIF
       !
    END SUBROUTINE lim_sbc_init
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limwri_dimg.h90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limwri_dimg.h90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limwri_dimg.h90	(revision 3402)
@@ -111,5 +111,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)
@@ -154,5 +154,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_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/ASM/asmtrj.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/ASM/asmtrj.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/ASM/asmtrj.F90	(revision 3402)
@@ -213,5 +213,5 @@
          !                                   ! Output trajectory fields
          CALL iom_rstput( it, it, inum, 'emp'   , emp    )
-         CALL iom_rstput( it, it, inum, 'emps'  , emps   )
+         CALL iom_rstput( it, it, inum, 'sfx'   , sfx    )
          CALL iom_rstput( it, it, inum, 'un'    , un     )
          CALL iom_rstput( it, it, inum, 'vn'    , vn     )
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 3402)
@@ -546,6 +546,6 @@
       CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
 !!$      CALL histwrite( nid_T, "sorunoff", it, runoff        , ndim_hT, ndex_hT )   ! runoff
-      CALL histwrite( nid_T, "sowaflcd", it, ( emps-rnf )  , ndim_hT, ndex_hT )   ! c/d water flux
-      zw2d(:,:) = ( emps(:,:) - rnf(:,:) ) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
+      CALL histwrite( nid_T, "sowaflcd", it, ( sfx -rnf )  , ndim_hT, ndex_hT )   ! c/d water flux
+      zw2d(:,:) = ( sfx (:,:) - rnf(:,:) ) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
       CALL histwrite( nid_T, "sosalflx", it, zw2d          , ndim_hT, ndex_hT )   ! c/d salt flux
       CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri_dimg.h90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri_dimg.h90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri_dimg.h90	(revision 3402)
@@ -54,5 +54,5 @@
     !!  level 14: qct(:,:)                 equivalent flux due to treshold SST
     !!  level 15: fbt(:,:)                 feedback term .
-    !!  level 16: ( emps(:,:) - rnf(:,:) ) concentration/dilution water flux
+    !!  level 16: ( sfx (:,:) - rnf(:,:) ) concentration/dilution water flux
     !!  level 17: fsalt(:,:)               Ice=>ocean net freshwater
     !!  level 18: gps(:,:)                 the surface pressure (m).
@@ -62,5 +62,5 @@
     !! History:  OPA  ! 1997-02 ( Clipper Group ) dimg files
     !!            -   ! 2003-12 ( J.M. Molines) f90, mpp output for OPA9.0
-    !!   NEMO    1.0  ! 2005-05  (S. Theetten) add emps fsalt move gps spgu spgv 2 lines below
+    !!   NEMO    1.0  ! 2005-05  (S. Theetten) add sfx  fsalt move gps spgu spgv 2 lines below
     !!            -   ! 2005-11  (V. Garnier) Surface pressure gradient organization
     !!----------------------------------------------------------------------
@@ -172,5 +172,5 @@
        !        fsel(:,:,14) = fsel(:,:,14) + qct(:,:)
        !        fsel(:,:,15) = fsel(:,:,15) + fbt(:,:)
-       fsel(:,:,16) = fsel(:,:,16) + ( emps(:,:)-rnf(:,:) ) 
+       fsel(:,:,16) = fsel(:,:,16) + ( sfx (:,:)-rnf(:,:) ) 
        !
        ! Output of dynamics and tracer fields and selected fields
@@ -240,5 +240,5 @@
           !         fsel(:,:,14) =  qct(:,:)
           !         fsel(:,:,15) =  fbt(:,:)
-          fsel(:,:,16) = ( emps(:,:)-rnf(:,:) ) * tmask(:,:,1) 
+          fsel(:,:,16) = ( sfx (:,:)-rnf(:,:) ) * tmask(:,:,1) 
           !
           !         qct(:,:) = 0._wp
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90	(revision 3402)
@@ -173,5 +173,5 @@
       !!      put as run-off in open ocean.
       !!
-      !! ** Action  :   emp, emps   updated surface freshwater fluxes at kt
+      !! ** Action  :   emp   updated surface freshwater flux at kt
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean model time step
@@ -214,5 +214,5 @@
       ENDIF
       !                                                   !--------------------!
-      !                                                   !  update emp, emps  !
+      !                                                   !  update emp        !
       zfwf = 0.e0                                         !--------------------!
       DO jc = 1, jpncs
@@ -236,5 +236,4 @@
             ! water/evap excess is shared by all open ocean
             emp (:,:) = emp (:,:) + zfwf(jc) / surf(jpncs+1)
-            emps(:,:) = emps(:,:) + zfwf(jc) / surf(jpncs+1)
          ELSEIF( ncstt(jc) == 1 ) THEN 
             ! Excess water in open sea, at outflow location, excess evap shared
@@ -247,11 +246,8 @@
                       emp (ji,jj) = emp (ji,jj) + zfwf(jc) /   &
                          (FLOAT(ncsnr(jc)) * e1t(ji,jj) * e2t(ji,jj))
-                      emps(ji,jj) = emps(ji,jj) + zfwf(jc) /   &
-                          (FLOAT(ncsnr(jc)) * e1t(ji,jj) * e2t(ji,jj))
                   END IF 
                 END DO 
             ELSE 
                 emp (:,:) = emp (:,:) + zfwf(jc) / surf(jpncs+1)
-                emps(:,:) = emps(:,:) + zfwf(jc) / surf(jpncs+1)
             ENDIF
          ELSEIF( ncstt(jc) == 2 ) THEN 
@@ -264,6 +260,4 @@
                   emp (ji,jj) = emp (ji,jj) + zfwf(jc)   &
                       / (FLOAT(ncsnr(jc)) *  e1t(ji,jj) * e2t(ji,jj) )
-                  emps(ji,jj) = emps(ji,jj) + zfwf(jc)   &
-                      / (FLOAT(ncsnr(jc)) *  e1t(ji,jj) * e2t(ji,jj) )
                 END DO 
             ENDIF 
@@ -273,5 +267,4 @@
             DO ji = ncsi1(jc), ncsi2(jc)
                emp (ji,jj) = emp (ji,jj) - zfwf(jc) / surf(jc)
-               emps(ji,jj) = emps(ji,jj) - zfwf(jc) / surf(jc)
             END DO  
          END DO 
@@ -280,5 +273,4 @@
       !
       CALL lbc_lnk( emp , 'T', 1. )
-      CALL lbc_lnk( emps, 'T', 1. )
       !
    END SUBROUTINE sbc_clo
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 3402)
@@ -81,6 +81,7 @@
       !
       INTEGER  ::   ji, jj, jk                             ! dummy loop indices
-      REAL(wp) ::   z2dt, zg_2                             ! temporary scalar
+      REAL(wp) ::   z2dt, zg_2, zintp, zgrau0r             ! temporary scalar
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zpice
       !!----------------------------------------------------------------------
       !
@@ -117,4 +118,28 @@
             END DO
          END DO
+      ENDIF
+
+      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: Pressure gradient due to snow-ice mass ==!
+         CALL wrk_alloc( jpi, jpj, zpice )
+         !                                            
+         zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc )
+         zgrau0r     = - grav * r1_rau0
+         zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrau0r
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               spgu(ji,jj) = ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj)
+               spgv(ji,jj) = ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj)
+            END DO
+         END DO
+         DO jk = 1, jpkm1                             ! Add the surface pressure trend to the general trend
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
+                  va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
+               END DO
+            END DO
+         END DO
+         !
+         CALL wrk_dealloc( jpi, jpj, zpice )
       ENDIF
 
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 3402)
@@ -94,6 +94,6 @@
    ! finally, arrays corresponding to different ice categories
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i                !: category ice fraction
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   topmelt           !: category topmelt
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   botmelt           !: category botmelt
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   topmelt            !: category topmelt
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   botmelt            !: category botmelt
 #endif
 
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 3402)
@@ -40,4 +40,5 @@
    LOGICAL , PUBLIC ::   ln_apr_dyn  = .FALSE.   !: Atmospheric pressure forcing used on dynamics (ocean & ice)
    INTEGER , PUBLIC ::   nn_ice      = 0         !: flag on ice in the surface boundary condition (=0/1/2/3)
+   INTEGER , PUBLIC ::   nn_ice_embd = 0         !: flag on ice embedding in the ocean (fully/partially/levitating) (=0/1/2)
    INTEGER , PUBLIC ::   nn_fwb      = 0         !: FreshWater Budget: 
    !                                             !:  = 0 unchecked 
@@ -61,5 +62,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qns_tot           !: total non solar heat flux (over sea and ice) [W/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   emp    , emp_b    !: freshwater budget: volume flux               [Kg/m2/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   emps   , emps_b   !: freshwater budget: concentration/dillution   [Kg/m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx    , emps_b   !: freshwater budget: salt flux                 [Kg/m2/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   emp_tot           !: total E-P over ocean and ice                 [Kg/m2/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rnf    , rnf_b    !: river runoff   [Kg/m2/s]  
@@ -108,5 +109,5 @@
          &      qsr_tot(jpi,jpj) , qsr   (jpi,jpj) ,                        &
          &      emp    (jpi,jpj) , emp_b (jpi,jpj) ,                        &
-         &      emps   (jpi,jpj) , emps_b(jpi,jpj) , emp_tot(jpi,jpj) , STAT=ierr(2) )
+         &      sfx    (jpi,jpj) , emps_b(jpi,jpj) , emp_tot(jpi,jpj) , STAT=ierr(2) )
          !
       ALLOCATE( rnf  (jpi,jpj) , sbc_tsc  (jpi,jpj,jpts) , qsr_hc  (jpi,jpj,jpk) ,     &
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90	(revision 3402)
@@ -90,4 +90,5 @@
          !
          emp (:,:) = rn_emp0
+         sfx (:,:) = 0.0_wp
          qns (:,:) = rn_qns0 - emp(:,:) * sst_m(:,:) * rcp      ! including heat content associated with mass flux at SST
          qsr (:,:) = rn_qsr0
@@ -129,5 +130,5 @@
       !!
       !! ** Action  : - set the ocean surface boundary condition, i.e.   
-      !!                   utau, vtau, taum, wndm, qns, qsr, emp
+      !!                   utau, vtau, taum, wndm, qns, qsr, emp, sfx
       !!
       !! Reference : Hazeleger, W., and S. Drijfhout, JPO, 30, 677-695, 2000.
@@ -224,4 +225,5 @@
       ! freshwater (mass flux) and update of qns with heat content of emp
       emp (:,:) = emp(:,:) - zsumemp * tmask(:,:,1)        ! freshwater flux (=0 in domain average)
+      sfx (:,:) = 0.0_wp                                   ! no salt flux
       qns (:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp   ! evap and precip are at SST
 
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 3402)
@@ -118,4 +118,6 @@
       !!              - qsr         solar heat flux
       !!              - emp         upward mass flux (evap. - precip)
+      !!              - sfx         salt flux; set to zero at nit000 but possibly non-zero
+      !!                            if ice is present (computed in limsbc(_2).F90)
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in  ) ::   kt   ! ocean time step
@@ -172,4 +174,6 @@
          ALLOCATE( sbudyko(jpi,jpj) , stauc(jpi,jpj), STAT=ierr3 )
          IF( ierr3 > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_clio: unable to allocate arrays' )
+         !
+         sfx(:,:) = 0.0_wp                      ! salt flux; zero unless ice is present (computed in limsbc(_2).F90)
          !
       ENDIF
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 3402)
@@ -111,6 +111,8 @@
       !!              - utau, vtau  i- and j-component of the wind stress
       !!              - taum, wndm  wind stress and 10m wind modules at T-point
-      !!              - qns, qsr    non-solar and solar heat flux
+      !!              - qns, qsr    non-solar and solar heat fluxes
       !!              - emp         upward mass flux (evapo. - precip.)
+      !!              - sfx         salt flux due to freezing/melting (non-zero only if ice is present)
+      !!                            (set in limsbc(_2).F90)
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
@@ -177,9 +179,11 @@
          CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' )
          !
-      ENDIF
-
-      CALL fld_read( kt, nn_fsbc, sf )        ! input fields provided at the current time-step
-
-      !                                                        ! surface ocean fluxes computed with CLIO bulk formulea
+         sfx(:,:) = 0.0_wp                         ! salt flux; zero unless ice is present (computed in limsbc(_2).F90)
+         !
+      ENDIF
+
+      CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
+
+      !                                            ! compute the surface ocean fluxes using CORE bulk formulea
       IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce_core( sf, sst_m, ssu_m, ssv_m )
 
@@ -248,5 +252,5 @@
       zcoef_qsatw = 0.98 * 640380. / rhoa
       
-      zst(:,:) = pst(:,:) + rt0      ! converte Celcius to Kelvin (and set minimum value far above 0 K)
+      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
 
       ! ----------------------------------------------------------------------------- !
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90	(revision 3402)
@@ -84,5 +84,5 @@
       !!              - wndm        10m wind module at T-point
       !!              - qns, qsr    non-slor and solar heat flux
-      !!              - emp, emps   evaporation minus precipitation
+      !!              - emp         evaporation minus precipitation
       !!----------------------------------------------------------------------
       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  sh_now   ! specific humidity at T-point 
@@ -258,5 +258,4 @@
            emp (:,:) = evap(:,:) - sf(jp_prec)%fnow(:,:,1) * tmask(:,:,1)
 !CDIR COLLAPSE
-           emps(:,:) = emp(:,:)
 
          CALL iom_put( "qlw_oce",   qbw  )                 ! output downward longwave heat over the ocean
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 3402)
@@ -823,5 +823,5 @@
          !                                                   ! ========================= !
          !
-         !                                                       ! total freshwater fluxes over the ocean (emp, emps)
+         !                                                       ! total freshwater fluxes over the ocean (emp, sfx )
          SELECT CASE( TRIM( sn_rcv_emp%cldes ) )                                    ! evaporation - precipitation
          CASE( 'conservative' )
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 3402)
@@ -70,4 +70,6 @@
       !!              - qsr         solar heat flux
       !!              - emp         upward mass flux (evap. - precip.)
+      !!              - sfx         salt flux; set to zero at nit000 but possibly non-zero
+      !!                            if ice is present (computed in limsbc(_2).F90)
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
@@ -120,4 +122,6 @@
          !                                         ! fill sf with slf_i and control print
          CALL fld_fill( sf, slf_i, cn_dir, 'sbc_flx', 'flux formulation for ocean surface boundary condition', 'namsbc_flx' )
+         !
+         sfx(:,:) = 0.0_wp                         ! salt flux due to freezing/melting (non-zero only if ice is present; set in limsbc(_2).F90)
          !
       ENDIF
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 3402)
@@ -59,4 +59,5 @@
       !!                =3 global mean of emp set to zero at each nn_fsbc time step
       !!                   & spread out over erp area depending its sign
+      !! Note: if sea ice is embedded it is taken into account when computing the budget 
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
@@ -88,4 +89,10 @@
          !
          area = glob_sum( e1e2t(:,:) )           ! interior global domain surface
+         !
+#if ! defined key_lim2 &&  ! defined key_lim3 && ! defined key_cice 
+         snwice_mass_b(:,:) = 0.e0               ! no sea-ice model is being used : no snow+ice mass
+         snwice_mass  (:,:) = 0.e0
+#endif
+         !
       ENDIF
       
@@ -96,5 +103,5 @@
          !
          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
-            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) ) ) / area   ! sum over the global domain
+            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) -  snwice_fmass(:,:) ) ) / area   ! sum over the global domain
             zcoef = z_fwf * rcp
             emp(:,:) = emp(:,:) - z_fwf 
@@ -119,6 +126,7 @@
          ikty = 365 * 86400 / rdttra(1)    !!bug  use of 365 days leap year or 360d year !!!!!!!
          IF( MOD( kt, ikty ) == 0 ) THEN
-            a_fwb_b = a_fwb
-            a_fwb   = glob_sum( e1e2t(:,:) * sshn(:,:) )   ! sum over the global domain
+            a_fwb_b = a_fwb                           ! mean sea level taking into account the ice+snow
+                                                      ! sum over the global domain
+            a_fwb   = glob_sum( e1e2t(:,:) * ( sshn(:,:) + snwice_mass(:,:) * r1_rau0 ) )
             a_fwb   = a_fwb * 1.e+3 / ( area * 86400. * 365. )     ! convert in Kg/m3/s = mm/s
 !!gm        !                                                      !!bug 365d year 
@@ -148,6 +156,6 @@
             zsurf_neg = glob_sum( e1e2t(:,:)*ztmsk_neg(:,:) )  ! Area filled by <0 and >0 erp 
             zsurf_pos = glob_sum( e1e2t(:,:)*ztmsk_pos(:,:) )
-            !                                                  ! fwf global mean 
-            z_fwf     = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) ) ) / area
+            !                                                  ! fwf global mean (excluding ocean to ice/snow exchanges) 
+            z_fwf     = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) - snwice_fmass(:,:) ) ) / area
             !            
             IF( z_fwf < 0._wp ) THEN         ! spread out over >0 erp area to increase evaporation
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 3402)
@@ -107,5 +107,5 @@
       !! ** Action  : - time evolution of the CICE sea-ice model
       !!              - update all sbc variables below sea-ice:
-      !!                utau, vtau, qns , qsr, emp , emps
+      !!                utau, vtau, qns , qsr, emp , sfx 
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt      ! ocean time step
@@ -480,5 +480,5 @@
 ! salinity calculation
       CALL cice2nemo(fresh_gbm,ztmp,'T', 1. )
-      emps(:,:)=emp(:,:)-ztmp(:,:)
+      sfx (:,:)=emp(:,:)-ztmp(:,:)
 ! Note the 1000.0 is to convert from kg salt to g salt (needed for PSU)
       CALL cice2nemo(fsalt_gbm,ztmp,'T', 1. )
@@ -486,5 +486,5 @@
          DO ji=1,jpi
             IF (sss_m(ji,jj).gt.0.0) THEN
-               emps(ji,jj)=emps(ji,jj)+ztmp(ji,jj)*1000.0/sss_m(ji,jj)
+               sfx (ji,jj)=sfx (ji,jj)+ztmp(ji,jj)*1000.0/sss_m(ji,jj)
             ENDIF
          ENDDO
@@ -494,5 +494,5 @@
 ! weight of the precip will affect the free surface even if it falls on the ice
 ! (same to the argument that freezing / melting of ice doesn't change the free surface) 
-! Sublimation from the ice is treated in a similar way (included in emp but not emps)  
+! Sublimation from the ice is treated in a similar way (included in emp but not sfx )  
 !
 ! This should not be done in the variable volume case
@@ -512,5 +512,5 @@
 
       CALL lbc_lnk( emp , 'T', 1. )
-      CALL lbc_lnk( emps , 'T', 1. )
+      CALL lbc_lnk( sfx  , 'T', 1. )
 
 ! Solar penetrative radiation and non solar surface heat flux
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90	(revision 3402)
@@ -51,5 +51,5 @@
       !!                taum, wndm : remain unchanged
       !!                qns, qsr   : update heat flux below sea-ice
-      !!                emp, emps  : update freshwater flux below sea-ice
+      !!                emp, sfx   : update freshwater flux below sea-ice
       !!                fr_i       : update the ice fraction
       !!---------------------------------------------------------------------
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 3402)
@@ -88,5 +88,5 @@
       !! ** Action  : - time evolution of the LIM sea-ice model
       !!              - update all sbc variables below sea-ice:
-      !!                utau, vtau, taum, wndm, qns , qsr, emp , emps
+      !!                utau, vtau, taum, wndm, qns , qsr, emp , sfx 
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt      ! ocean time step
@@ -392,5 +392,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF( ABS( emps(ji,jj) ) .GT. 1.0e-2 ) THEN
+            IF( ABS( sfx (ji,jj) ) .GT. 1.0e-2 ) THEN
                CALL lim_prt_state( ji, jj, 3, ' ALERTE 5 :   High salt flux ' )
                DO jl = 1, jpl
@@ -644,5 +644,5 @@
         WRITE(numout,*)
         WRITE(numout,*) ' - Salt fluxes at bottom interface ***'
-        WRITE(numout,*) ' emps       : ', emps(ki,kj)
+        WRITE(numout,*) ' sfx        : ', sfx (ki,kj)
         WRITE(numout,*) ' emp        : ', emp(ki,kj)
         WRITE(numout,*) ' fsbri      : ', fsbri(ki,kj)
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90	(revision 3402)
@@ -82,5 +82,5 @@
       !! ** Action  : - time evolution of the LIM sea-ice model
       !!              - update all sbc variables below sea-ice:
-      !!                utau, vtau, taum, wndm, qns , qsr, emp , emps
+      !!                utau, vtau, taum, wndm, qns , qsr, emp , sfx 
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt      ! ocean time step
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 3402)
@@ -82,6 +82,7 @@
       INTEGER ::   icpt   ! local integer
       !!
-      NAMELIST/namsbc/ nn_fsbc   , ln_ana , ln_flx  , ln_blk_clio, ln_blk_core, ln_cpl,   &
-         &             ln_blk_mfs, ln_apr_dyn, nn_ice , ln_dm2dc, ln_rnf, ln_ssr     , nn_fwb, ln_cdgw
+      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx,  ln_blk_clio, ln_blk_core, ln_cpl,   &
+         &             ln_blk_mfs, ln_apr_dyn, nn_ice,  nn_ice_embd, ln_dm2dc   , ln_rnf,   &
+         &             ln_ssr    , nn_fwb    , ln_cdgw
       !!----------------------------------------------------------------------
 
@@ -119,4 +120,5 @@
          WRITE(numout,*) '              Patm gradient added in ocean & ice Eqs.    ln_apr_dyn  = ', ln_apr_dyn
          WRITE(numout,*) '              ice management in the sbc (=0/1/2/3)       nn_ice      = ', nn_ice 
+         WRITE(numout,*) '              ice-ocean embedded/levitating (=0/1/2)     nn_ice_embd = ', nn_ice_embd
          WRITE(numout,*) '              daily mean to diurnal cycle qsr            ln_dm2dc    = ', ln_dm2dc 
          WRITE(numout,*) '              runoff / runoff mouths                     ln_rnf      = ', ln_rnf
@@ -134,12 +136,12 @@
          IF( sbc_rnf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_rnf arrays' )
          nkrnf         = 0
-         rnf     (:,:) = 0.e0
-         rnfmsk  (:,:) = 0.e0
-         rnfmsk_z(:)   = 0.e0
+         rnf     (:,:) = 0.0_wp
+         rnfmsk  (:,:) = 0.0_wp
+         rnfmsk_z(:)   = 0.0_wp
       ENDIF
       IF( nn_ice == 0  )   fr_i(:,:) = 0.e0        ! no ice in the domain, ice fraction is always zero
 
-      emps(:,:) = 0.e0                             ! the salt flux will be computed (i.e. will be non-zero) only if 
-      !                                            ! sea-ice is present, or lk_vvl=F, or surface salt restoring is used.
+      sfx(:,:) = 0.0_wp                            ! the salt flux due to freezing/melting will be computed (i.e. will be non-zero) 
+                                                   ! only if sea-ice is present
 
       !                                            ! restartability   
@@ -221,5 +223,5 @@
       !!                time step, i.e.  
       !!                utau_b, vtau_b, qns_b, qsr_b, emp_n, emps_b, qrp_b, erp_b
-      !!                utau  , vtau  , qns  , qsr  , emp  , emps  , qrp  , erp
+      !!                utau  , vtau  , qns  , qsr  , emp  , sfx   , qrp  , erp
       !!              - updte the ice fraction : fr_i
       !!----------------------------------------------------------------------
@@ -238,5 +240,5 @@
          ! qsr_b (:,:) = qsr (:,:)
          emp_b (:,:) = emp (:,:)
-         emps_b(:,:) = emps(:,:)
+         emps_b(:,:) = sfx(:,:)
       ENDIF
       !                                            ! ---------------------------------------- !
@@ -256,5 +258,5 @@
                                                             
       SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition
-      !                                                  ! (i.e. utau,vtau, qns, qsr, emp, emps)
+      !                                                  ! (i.e. utau,vtau, qns, qsr, emp, sfx)
       CASE(  0 )   ;   CALL sbc_gyre    ( kt )                    ! analytical formulation : GYRE configuration
       CASE(  1 )   ;   CALL sbc_ana     ( kt )                    ! analytical formulation : uniform sbc
@@ -316,5 +318,5 @@
             ! qsr_b (:,:) = qsr (:,:)
             emp_b (:,:) = emp (:,:)
-            emps_b(:,:) = emps(:,:)
+            emps_b(:,:) = sfx(:,:)
          ENDIF
       ENDIF
@@ -332,5 +334,5 @@
          ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  )
          CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp  )
-         CALL iom_rstput( kt, nitrst, numrow, 'emps_b' , emps )
+         CALL iom_rstput( kt, nitrst, numrow, 'emps_b' , sfx )
       ENDIF
 
@@ -340,5 +342,5 @@
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
          CALL iom_put( "empmr" , emp  - rnf )                   ! upward water flux
-         CALL iom_put( "empsmr", emps - rnf )                   ! c/d water flux
+         CALL iom_put( "empsmr", sfx - rnf )                    ! c/d water flux
          CALL iom_put( "qt"    , qns  + qsr )                   ! total heat flux 
          CALL iom_put( "qns"   , qns        )                   ! solar heat flux
@@ -357,5 +359,5 @@
          CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i     - : ', mask1=tmask, ovlap=1 )
          CALL prt_ctl(tab2d_1=(emp-rnf)        , clinfo1=' emp-rnf  - : ', mask1=tmask, ovlap=1 )
-         CALL prt_ctl(tab2d_1=(emps-rnf)       , clinfo1=' emps-rnf - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=(sfx-rnf)        , clinfo1=' sfx-rnf  - : ', mask1=tmask, ovlap=1 )
          CALL prt_ctl(tab2d_1=qns              , clinfo1=' qns      - : ', mask1=tmask, ovlap=1 )
          CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr      - : ', mask1=tmask, ovlap=1 )
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 3402)
@@ -63,6 +63,6 @@
       !!              - at each nscb time step
       !!                   add a retroaction term on qns    (nn_sstr = 1)
-      !!                   add a damping term on emps       (nn_sssr = 1)
-      !!                   add a damping term on emp & emps (nn_sssr = 2)
+      !!                   add a damping term on sfx        (nn_sssr = 1)
+      !!                   add a damping term on emp        (nn_sssr = 2)
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in   ) ::   kt   ! ocean time step
@@ -168,5 +168,5 @@
             ENDIF
             !
-            IF( nn_sssr == 1 ) THEN                                   !* Salinity damping term (salt flux only (emps))
+            IF( nn_sssr == 1 ) THEN                                   !* Salinity damping term (salt flux only (sfx))
                zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s]
 !CDIR COLLAPSE
@@ -175,6 +175,6 @@
                      zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
                         &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) 
-                     emps(ji,jj) = emps(ji,jj) + zerp                 ! salt flux
-                     erp( ji,jj) = zerp / MAX( sss_m(ji,jj), 1.e-20 ) ! converted into an equivalent volume flux (diagnostic only)
+                     sfx(ji,jj) = sfx(ji,jj) + zerp                 ! salt flux
+                     erp(ji,jj) = zerp / MAX( sss_m(ji,jj), 1.e-20 ) ! converted into an equivalent volume flux (diagnostic only)
                   END DO
                END DO
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 3402)
@@ -61,5 +61,5 @@
       !!      (3) Fwe, tracer carried with the water that is exchanged. 
       !!            - salinity    : salt flux only due to freezing/melting
-      !!            sa = sa +  fsalt / rau0 / e3t  for k=1
+      !!            sa = sa +  sfx / rau0 / e3t  for k=1
       !!
       !!      Fext, flux through the air-sea interface for temperature and salt: 
@@ -89,8 +89,8 @@
       !!         freezing/melting, thus:
       !!            sa = sa + emp * sn / rau0 / e3t   for k=1
-      !!                    + fsalt    / rau0 / e3t
+      !!                    + sfx    / rau0 / e3t
       !!         where emp, the surface freshwater budget (evaporation minus
       !!         precipitation minus runoff) given in kg/m2/s is divided
-      !!         by 1035 kg/m3 (density of ocena water) to obtain m/s.    
+      !!         by rau0 = 1020 kg/m3 (density of sea water) to obtain m/s.    
       !!         Note: even though Fwe does not appear explicitly for 
       !!         temperature in this routine, the heat carried by the water
@@ -169,8 +169,8 @@
             DO ji = 1, jpi 
                sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)                              ! non solar heat flux
-               sbc_tsc(ji,jj,jp_sal) = r1_rau0     * emps(ji,jj)                             ! salt flux (freezing/melting)
+               sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)                              ! salt flux due to freezing/melting
             END DO
          END DO
-      ELSE                                         ! Constant Volume case
+      ELSE                                         ! Constant Volume case ==>> Concentration dilution effect
          DO jj = 2, jpj
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -178,7 +178,7 @@
                sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)                          &   ! non solar heat flux
                   &                  + r1_rau0     * emp(ji,jj)  * tsn(ji,jj,1,jp_tem)       ! concent./dilut. effect
-               ! salinity    : salt flux + concent./dilut. effect (both in emps)
-               sbc_tsc(ji,jj,jp_sal) = r1_rau0  * (  emps(ji,jj)                         &   ! salt flux (freezing/melting)
-                  &                                + emp (ji,jj) * tsn(ji,jj,1,jp_sal) )     ! concent./dilut. effect
+               ! salinity    : salt flux + concent./dilut. effect (both in sfx)
+               sbc_tsc(ji,jj,jp_sal) = r1_rau0  * (  sfx(ji,jj)                          &   ! salt flux (freezing/melting)
+                  &                                + emp(ji,jj) * tsn(ji,jj,1,jp_sal) )      ! concent./dilut. effect
             END DO
          END DO
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90	(revision 3402)
@@ -426,9 +426,11 @@
             zBosol(ji,jj) = grav * zthermal * qsr(ji,jj)
             ! Non radiative surface buoyancy force
-            zBo   (ji,jj) = grav * zthermal * qns(ji,jj) -  grav * zhalin * ( emps(ji,jj)-rnf(ji,jj) ) 
+            zBo   (ji,jj) = grav * zthermal * qns(ji,jj) -  grav * zhalin * ( emp(ji,jj)-rnf(ji,jj) )  &
+               &                                         -  grav * rbeta * rcs * sfx(ji,jj)
             ! Surface Temperature flux for non-local term
             wt0(ji,jj) = - ( qsr(ji,jj) + qns(ji,jj) )* r1_rau0_rcp * tmask(ji,jj,1)
             ! Surface salinity flux for non-local term
-            ws0(ji,jj) = - ( ( emps(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal) * rcs ) * tmask(ji,jj,1) 
+            ws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal)                          &
+               &             + sfx(ji,jj)                                     ) * rcs * tmask(ji,jj,1) 
          ENDDO
       ENDDO
@@ -1324,5 +1326,5 @@
                DO ji = fs_2, fs_jpim1
                   ! Surface tracer flux for non-local term 
-                  zflx = - ( emps(ji,jj) * tra(ji,jj,1,jn) * rcs ) * tmask(ji,jj,1)
+                  zflx = - ( sfx (ji,jj) * tra(ji,jj,1,jn) * rcs ) * tmask(ji,jj,1)
                   ! compute the trend
                   ztra = - ( ghats(ji,jj,jk  ) * fsavs(ji,jj,jk  )   &
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/oce.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 3401)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 3402)
@@ -47,4 +47,12 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   gru , grv    !: horizontal gradient of rd at bottom u-point
 
+   !! arrays relating to embedding ice in the ocean. These arrays need to be declared 
+   !! even if no ice model is required. In the no ice model or traditional levitating 
+   !! ice cases they contain only zeros
+   !! ---------------------
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   snwice_mass        !: mass of snow and ice at current  ice time step   [Kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   snwice_mass_b      !: mass of snow and ice at previous ice time step   [Kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   snwice_fmass       !: time evolution of mass of snow+ice               [Kg/m2/s]
+
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011)
@@ -58,5 +66,5 @@
       !!                   ***  FUNCTION oce_alloc  ***
       !!----------------------------------------------------------------------
-      INTEGER :: ierr(2)
+      INTEGER :: ierr(3)
       !!----------------------------------------------------------------------
       !
@@ -69,13 +77,16 @@
          &      rn2b (jpi,jpj,jpk)      , rn2  (jpi,jpj,jpk)                              , STAT=ierr(1) )
          !
-      ALLOCATE(rhd (jpi,jpj,jpk) ,                                         &
-         &     rhop(jpi,jpj,jpk) ,                                         &
-         &     sshb  (jpi,jpj)   , sshn  (jpi,jpj) , ssha  (jpi,jpj) ,     &
-         &     sshu_b(jpi,jpj)   , sshu_n(jpi,jpj) , sshu_a(jpi,jpj) ,     &
-         &     sshv_b(jpi,jpj)   , sshv_n(jpi,jpj) , sshv_a(jpi,jpj) ,     &
-         &                         sshf_n(jpi,jpj) ,                       &
-         &     spgu  (jpi,jpj)   , spgv(jpi,jpj)   ,                       &
-         &     gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts),                     &
-         &     gru(jpi,jpj)      , grv(jpi,jpj)                      , STAT=ierr(2) )
+      ALLOCATE( rhd (jpi,jpj,jpk) ,                                         &
+         &      rhop(jpi,jpj,jpk) ,                                         &
+         &      sshb  (jpi,jpj)   , sshn  (jpi,jpj) , ssha  (jpi,jpj) ,     &
+         &      sshu_b(jpi,jpj)   , sshu_n(jpi,jpj) , sshu_a(jpi,jpj) ,     &
+         &      sshv_b(jpi,jpj)   , sshv_n(jpi,jpj) , sshv_a(jpi,jpj) ,     &
+         &                          sshf_n(jpi,jpj) ,                       &
+         &      spgu  (jpi,jpj)   , spgv(jpi,jpj)   ,                       &
+         &      gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts),                     &
+         &      gru(jpi,jpj)      , grv(jpi,jpj)                      , STAT=ierr(2) )
+         !
+      ALLOCATE( snwice_mass(jpi,jpj)  , snwice_mass_b(jpi,jpj),             &
+         &      snwice_fmass(jpi,jpj), STAT= ierr(3) )
          !
       oce_alloc = MAXVAL( ierr )
