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 3419)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 3488)
@@ -221,5 +221,5 @@
             !
             ! 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
+            ! 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 no ice/ocean mass exchange)
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 3419)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 3488)
@@ -34,6 +34,6 @@
    USE prtctl           ! Print control
    USE cpl_oasis3, ONLY : lk_cpl
-   USE oce, ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass
-   USE dom_ice, ONLY : tms
+   USE oce,        ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass, sshu_b, sshv_b, sshu_n, sshv_n, sshf_n
+   USE dom_ice,    ONLY : tms
 
    IMPLICIT NONE
@@ -395,4 +395,6 @@
       !!-------------------------------------------------------------------
       !
+      INTEGER  ::   ji, jj                          ! dummy loop indices
+      REAL(wp) ::   zcoefu, zcoefv, zcoeff          ! local scalar
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'lim_sbc_init : LIM-3 sea-ice - surface boundary condition'
@@ -422,8 +424,39 @@
          snwice_mass_b(:,:) = 0.0_wp         ! 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
+      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
+         !
+         ! Note: Changed the initial values of sshb and sshn=>  need to recompute ssh[u,v,f]_[b,n] 
+         !       which were previously set in domvvl
+         IF ( lk_vvl ) THEN            ! Is this necessary? embd 2 should be restricted to vvl only???
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
+                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
+                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
+                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
+                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
+                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
+                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
+                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
+                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
+                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
+                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
+                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
+               END DO
+            END DO
+            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
+            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1      ! NO Vector Opt.
+                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
+                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
+                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
+                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
+               END DO
+            END DO
+            CALL lbc_lnk( sshf_n, 'F', 1. )
+          ENDIF
       ENDIF
       !
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 3419)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 3488)
@@ -372,18 +372,19 @@
          CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-!!$#if defined key_lim3 || defined key_lim2 
-!!$         ! sowaflup = sowaflep + sorunoff + sowafldp + a term associated to
-!!$         !    internal damping to Levitus that can be diagnosed from others
-!!$         ! sowaflcd = sowaflep + sorunoff + sowafldp + iowaflup
-!!$         CALL histdef( nid_T, "iowaflup", "Ice=>ocean net freshwater"          , "kg/m2/s",   &  ! fsalt
-!!$            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-!!$         CALL histdef( nid_T, "sowaflep", "atmos=>ocean net freshwater"        , "kg/m2/s",   &  ! fmass
-!!$            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-!!$#endif
+>>>>>>>>>>>
+      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
+      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux 
+                                                                                  ! (includes virtual salt flux beneath ice 
+                                                                                  ! in linear free surface case)
+#if ! defined key_vvl
+      zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem) * tmask(:,:,1)
+      CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )             ! c/d term on sst
+      zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
+      CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )             ! c/d term on sss
+#endif
+<<<<<<<<<<
          CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-!!$         CALL histdef( nid_T, "sorunoff", "Runoffs"                            , "Kg/m2/s",   &  ! runoffs
-!!$            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "sowaflcd", "concentration/dilution water flux"  , "kg/m2/s",   &  ! (emps-rnf)
+         CALL histdef( nid_T, "sofldow", "downward salt flux"                  , "PSU/m2/s",  &  ! sfx
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
          CALL histdef( nid_T, "sosalflx", "Surface Salt Flux"                  , "Kg/m2/s",   &  ! (emps-rnf) * sn
@@ -540,13 +541,14 @@
       CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal), ndim_hT, ndex_hT )   ! sea surface salinity
       CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
-!!$#if  defined key_lim3 || defined key_lim2 
-!!$      CALL histwrite( nid_T, "iowaflup", it, fsalt(:,:)    , ndim_hT, ndex_hT )   ! ice=>ocean water flux
-!!$      CALL histwrite( nid_T, "sowaflep", it, fmass(:,:)    , ndim_hT, ndex_hT )   ! atmos=>ocean water flux
-!!$#endif
       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, ( 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, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux 
+                                                                                  ! (includes virtual salt flux beneath ice 
+                                                                                  ! in linear free surface case)
+#if ! defined key_vvl
+      zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
+      CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )             ! c/d term on sst
+      zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
+      CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )             ! c/d term on sss
+#endif
       CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
       CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
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 3419)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90	(revision 3488)
@@ -178,5 +178,5 @@
       !
       INTEGER                     ::   ji, jj, jc, jn   ! dummy loop indices
-      REAL(wp)                    ::   zze2
+      REAL(wp)                    ::   zze2, zcoef, zcoef1
       REAL(wp), DIMENSION (jpncs) ::   zfwf 
       !!----------------------------------------------------------------------
@@ -235,5 +235,8 @@
          IF( ncstt(jc) == 0 ) THEN 
             ! water/evap excess is shared by all open ocean
-            emp (:,:) = emp (:,:) + zfwf(jc) / surf(jpncs+1)
+            zcoef  = zfwf(jc) / surf(jpncs+1)
+            zcoef1 = rcp * zcoef
+            emp(:,:) = emp(:,:) + zcoef
+            qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
          ELSEIF( ncstt(jc) == 1 ) THEN 
             ! Excess water in open sea, at outflow location, excess evap shared
@@ -244,10 +247,15 @@
                   IF (      ji > 1 .AND. ji < jpi   &
                       .AND. jj > 1 .AND. jj < jpj ) THEN 
-                      emp (ji,jj) = emp (ji,jj) + zfwf(jc) /   &
-                         (FLOAT(ncsnr(jc)) * e1t(ji,jj) * e2t(ji,jj))
+                      zcoef  = zfwf(jc) / ( REAL(ncsnr(jc), wp) * e1t(ji,jj) * e2t(ji,jj) )
+                      zcoef1 = rcp * zcoef
+                      emp(ji,jj) = emp(ji,jj) + zcoef
+                      qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj)
                   END IF 
                 END DO 
             ELSE 
-                emp (:,:) = emp (:,:) + zfwf(jc) / surf(jpncs+1)
+                zcoef  = zfwf(jc) / surf(jpncs+1)
+                zcoef1 = rcp * zcoef
+                emp(:,:) = emp(:,:) + zcoef
+                qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
             ENDIF
          ELSEIF( ncstt(jc) == 2 ) THEN 
@@ -258,6 +266,8 @@
                   ji = mi0(ncsir(jc,jn))
                   jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean
-                  emp (ji,jj) = emp (ji,jj) + zfwf(jc)   &
-                      / (FLOAT(ncsnr(jc)) *  e1t(ji,jj) * e2t(ji,jj) )
+                  zcoef  = zfwf(jc) / ( REAL(ncsnr(jc), wp) * e1t(ji,jj) * e2t(ji,jj) )
+                  zcoef1 = rcp * zcoef
+                  emp(ji,jj) = emp(ji,jj) + zcoef
+                  qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj)
                 END DO 
             ENDIF 
@@ -266,5 +276,8 @@
          DO jj = ncsj1(jc), ncsj2(jc)
             DO ji = ncsi1(jc), ncsi2(jc)
-               emp (ji,jj) = emp (ji,jj) - zfwf(jc) / surf(jc)
+               zcoef  = zfwf(jc) / surf(jc)
+               zcoef1 = rcp * zcoef
+               emp(ji,jj) = emp(ji,jj) - zcoef
+               qns(ji,jj) = qns(ji,jj) + zcoef1 * sst_m(ji,jj)
             END DO  
          END DO 
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 3419)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 3488)
@@ -39,6 +39,9 @@
    LOGICAL , PUBLIC ::   ln_ssr      = .FALSE.   !: Sea Surface restoring on SST and/or SSS      
    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_ice      = 0         !: flag for ice in the surface boundary condition (=0/1/2/3)
+   INTEGER , PUBLIC ::   nn_ice_embd = 0         !: flag for levitating/embedding sea-ice in the ocean
+   !                                             !: =0 levitating ice (no mass exchange, concentration/dilution effect)
+   !                                             !: =1 levitating ice with mass and salt exchange but no presure effect
+   !                                             !: =2 embedded sea-ice (full salt and mass exchanges and pressure)
    INTEGER , PUBLIC ::   nn_fwb      = 0         !: FreshWater Budget: 
    !                                             !:  = 0 unchecked 
@@ -62,5 +65,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(:,:) ::   sfx    , emps_b   !: freshwater budget: salt flux                 [Kg/m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx    , sfx_b    !: salt flux                                    [PSU/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]  
@@ -106,8 +109,8 @@
          &      vtau(jpi,jpj) , vtau_b(jpi,jpj) , wndm(jpi,jpj) , STAT=ierr(1) ) 
          !
-      ALLOCATE( qns_tot(jpi,jpj) , qns   (jpi,jpj) , qns_b(jpi,jpj),        &
-         &      qsr_tot(jpi,jpj) , qsr   (jpi,jpj) ,                        &
-         &      emp    (jpi,jpj) , emp_b (jpi,jpj) ,                        &
-         &      sfx    (jpi,jpj) , emps_b(jpi,jpj) , emp_tot(jpi,jpj) , STAT=ierr(2) )
+      ALLOCATE( qns_tot(jpi,jpj) , qns  (jpi,jpj) , qns_b(jpi,jpj),        &
+         &      qsr_tot(jpi,jpj) , qsr  (jpi,jpj) ,                        &
+         &      emp    (jpi,jpj) , emp_b(jpi,jpj) ,                        &
+         &      sfx    (jpi,jpj) , sfx_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/sbccpl.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 3419)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 3488)
@@ -1145,6 +1145,5 @@
 
       zicefr(:,:) = 1.- p_frld(:,:)
-      !! zcptn(:,:) = rcp * sst_m(:,:)                       !! safer to use instantaneous sst value in case coupling frequency does not match nn_fsbc frequency ??
-      zcptn(:,:) = rcp * tsn(:,:,1,jp_tem)
+      zcptn(:,:) = rcp * sst_m(:,:)
       !
       !                                                      ! ========================= !
@@ -1241,6 +1240,6 @@
       qns_tot(:,:) = qns_tot(:,:)                         &            ! qns_tot update over free ocean with:
          &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting
-         &          + (  emp_tot(:,:)                     &            ! remove the heat content of mass flux (assumed to be at SST)
-         &             - emp_ice(:,:) * p_frld(:,:,1)  ) * zcptn(:,:) 
+         &          - (  emp_tot(:,:)                     &            ! remove the heat content of mass flux (assumed to be at SST)
+         &             - emp_ice(:,:) * zicefr(:,:)  ) * zcptn(:,:) 
       IF( lk_diaar5 )   CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
 !!gm
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 3419)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 3488)
@@ -106,11 +106,11 @@
             zcoef = z_fwf * rcp
             emp(:,:) = emp(:,:) - z_fwf 
-            qns(:,:) = qns(:,:) + zcoef * sst_m(:,:)                          ! ensure fw correction does not change the heat budget
+            qns(:,:) = qns(:,:) + zcoef * sst_m(:,:)  ! account for change to the heat budget due to fw correction
          ENDIF
          !
       CASE ( 2 )                             !==  fwf budget adjusted from the previous year  ==!
          !
-         IF( kt == nit000 ) THEN                   ! initialisation
-            !                                      ! Read the corrective factor on precipitations (fwfold)
+         IF( kt == nit000 ) THEN                      ! initialisation
+            !                                         ! Read the corrective factor on precipitations (fwfold)
             CALL ctl_opn( inum, 'EMPave_old.dat', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
             READ ( inum, "(24X,I8,2ES24.16)" ) iyear, a_fwb_b, a_fwb
@@ -138,5 +138,5 @@
             zcoef = fwfold * rcp
             emp(:,:) = emp(:,:) + fwfold
-            qns(:,:) = qns(:,:) - zcoef * sst_m(:,:)  ! ensure fw correction does not change the heat budget
+            qns(:,:) = qns(:,:) - zcoef * sst_m(:,:)  ! account for change to the heat budget due to fw correction
          ENDIF
          !
@@ -180,5 +180,5 @@
             !
             emp(:,:) = emp(:,:) + zerp_cor(:,:)
-            qns(:,:) = qns(:,:) - zerp_cor(:,:) * rcp * sst_m(:,:)  
+            qns(:,:) = qns(:,:) - zerp_cor(:,:) * rcp * sst_m(:,:)  ! account for change to the heat budget due to fw correction
             erp(:,:) = erp(:,:) + zerp_cor(:,:)
             !
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 3419)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 3488)
@@ -137,4 +137,5 @@
 
    END SUBROUTINE sbc_ice_cice
+
 
    SUBROUTINE cice_sbc_init (nsbc)
@@ -477,4 +478,28 @@
       ENDIF
 
+!! ACC this block needs attention. sfx has replaced emps but its meaning is now
+!! different. Need the equivalent of this block:
+!     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
+!
+!     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)
+!
+! Here zfsalt is salt flux between ice and ocean due to freezing and melting (PSU/m2/s)
+!      zcd    is the virtual salt flux that appears in the standard levitating case only (0 otherwise)
+!             it generates a change in SSS equivalent to that which would occur if a true mass exchange happened
+!      zemp   is: 
+!              in coupled mode:  net mass flux over the grid cell (ice+ocean area) minus the mass flux intercepted by sea-ice
+!              in forced  mode:  mass flux budget (emp) over open ocean fraction minus liquid precip. over the ice 
+!                                (assumed to instantaneous drain into the ocean).
+!    zemp_snw is the snow melt that enters the ocean as pure water (no associated salt flux) 
+!    zfmm     is freezing minus melting
+!
+! with zswitch = 1 ( nn_ice_embd = 0 ) the results should be equivalent to the original CICE code. 
+
+!!
 ! Subtract fluxes from CICE to get freshwater equivalent flux used in 
 ! salinity calculation
@@ -498,4 +523,5 @@
 ! This should not be done in the variable volume case
 
+
       IF (.NOT. lk_vvl) THEN
 
@@ -510,8 +536,14 @@
 
       ENDIF
+!! ACC end of questionable code
 
       CALL lbc_lnk( emp , 'T', 1. )
       CALL lbc_lnk( sfx  , 'T', 1. )
 
+!! ACC Now the latent heat for snow melting is already accounted for in the bulk formulea and coupled interfaces.
+!!     For the non-solar heat flux, in LIM2, code changes were needed to account for the heat content of the mass exchanged 
+!      between ice and ocean. It was not necessary to make changes for LIM3 since all mass exchanges are referenced to 
+!      zero degrees; this is most likely to be the case in CICE too??
+!! 
 ! Solar penetrative radiation and non solar surface heat flux
 
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 3419)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 3488)
@@ -12,4 +12,5 @@
    !!             -   ! 2010-10  (J. Chanut, C. Bricaud, G. Madec)  add the surface pressure forcing
    !!            3.4  ! 2011-11  (C. Harris) CICE added as an option
+   !!            3.5  ! 2012-11  (A. Coward, G. Madec) Rethink of heat, mass and salt surface fluxes
    !!----------------------------------------------------------------------
 
@@ -162,4 +163,6 @@
       IF( nn_ice == 4 .AND. ( .NOT. ( cp_cfg == 'orca' ) .OR. lk_agrif ) )   &
          &   CALL ctl_stop( 'CICE sea-ice model currently only available in a global ORCA configuration without AGRIF' )
+      IF( nn_ice == 3 .AND. nn_ice_embd == 0 )   &
+         &   CALL ctl_stop( 'LIM3 sea-ice model requires nn_ice_embd = 2 or 3' )
       
       IF( ln_dm2dc )   nday_qsr = -1   ! initialisation flag
@@ -222,6 +225,6 @@
       !! ** Action  : - set the ocean surface boundary condition at before and now 
       !!                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  , sfx   , qrp  , erp
+      !!                utau_b, vtau_b, qns_b, qsr_b, emp_n, sfx_b, qrp_b, erp_b
+      !!                utau  , vtau  , qns  , qsr  , emp  , sfx  , qrp  , erp
       !!              - updte the ice fraction : fr_i
       !!----------------------------------------------------------------------
@@ -239,6 +242,6 @@
          ! The 3D heat content due to qsr forcing is treated in traqsr
          ! qsr_b (:,:) = qsr (:,:)
-         emp_b (:,:) = emp (:,:)
-         emps_b(:,:) = sfx(:,:)
+         emp_b(:,:) = emp(:,:)
+         sfx_b(:,:) = sfx(:,:)
       ENDIF
       !                                            ! ---------------------------------------- !
@@ -308,7 +311,12 @@
             CALL iom_get( numror, jpdom_autoglo, 'qns_b' , qns_b  )   ! before non solar heat flux (T-point)
             ! The 3D heat content due to qsr forcing is treated in traqsr
-            ! CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b  )   ! before     solar heat flux (T-point)
-            CALL iom_get( numror, jpdom_autoglo, 'emp_b' , emp_b  )   ! before     freshwater flux (T-point)
-            CALL iom_get( numror, jpdom_autoglo, 'emps_b', emps_b )   ! before C/D freshwater flux (T-point)
+            ! CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b  ) ! before     solar heat flux (T-point)
+            CALL iom_get( numror, jpdom_autoglo, 'emp_b', emp_b  )    ! before     freshwater flux (T-point)
+            ! To ensure restart capability with 3.3x/3.4 restart files    !! to be removed in v3.6
+            IF( iom_varid( numror, 'sfx_b', ldstop = .FALSE. ) > 0 ) THEN
+               CALL iom_get( numror, jpdom_autoglo, 'sfx_b', sfx_b )  ! before salt flux (T-point)
+            ELSE
+               sfx_b (:,:) = sfx(:,:)
+            ENDIF
          ELSE                                                   !* no restart: set from nit000 values
             IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields set to nit000'
@@ -316,7 +324,6 @@
             vtau_b(:,:) = vtau(:,:)
             qns_b (:,:) = qns (:,:)
-            ! qsr_b (:,:) = qsr (:,:)
-            emp_b (:,:) = emp (:,:)
-            emps_b(:,:) = sfx(:,:)
+            emp_b (:,:) = emp(:,:)
+            sfx_b (:,:) = sfx(:,:)
          ENDIF
       ENDIF
@@ -334,5 +341,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' , sfx )
+         CALL iom_rstput( kt, nitrst, numrow, 'sfx_b' , sfx )
       ENDIF
 
@@ -342,5 +349,7 @@
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
          CALL iom_put( "empmr" , emp  - rnf )                   ! upward water flux
-         CALL iom_put( "empsmr", sfx - rnf )                    ! c/d water flux
+         CALL iom_put( "saltflx", sfx  )                        ! downward salt flux  
+                                                                ! (includes virtual salt flux beneath ice 
+                                                                ! in linear free surface case)
          CALL iom_put( "qt"    , qns  + qsr )                   ! total heat flux 
          CALL iom_put( "qns"   , qns        )                   ! solar heat flux
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 3419)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 3488)
@@ -183,4 +183,6 @@
             END DO
          END DO
+         CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) )                          ! c/d term on sst
+         CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) )                          ! c/d term on sss
       ENDIF
       ! Concentration dilution effect on (t,s) due to evapouration, precipitation and qns, but not river runoff  
Index: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 3419)
+++ branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 3488)
@@ -318,7 +318,5 @@
       IF( lk_bdy        )   CALL     tide_init      ! Open boundaries initialisation of tidal harmonic forcing
 
-                            CALL flush(numout)
                             CALL dyn_nept_init  ! simplified form of Neptune effect
-                            CALL flush(numout)
 
                             CALL  istate_init   ! ocean initial state (Dynamics and tracers)
