Index: trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 4688)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 4689)
@@ -109,4 +109,9 @@
                                                               !: = 1 read it in a NetCDF file
 #endif
+#if defined key_lim3
+   REAL,    DIMENSION(jp_bdy) ::   rn_ice_tem               !: choice of the temperature of incoming sea ice
+   REAL,    DIMENSION(jp_bdy) ::   rn_ice_sal               !: choice of the salinity    of incoming sea ice
+   REAL,    DIMENSION(jp_bdy) ::   rn_ice_age               !: choice of the age         of incoming sea ice
+#endif
    !
    
Index: trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 4688)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 4689)
@@ -676,5 +676,5 @@
                CALL iom_close ( inum )
                !CALL fld_clopn ( bn_a_i, nyear, nmonth, nday, ldstop=.TRUE. )
-               !CALL iom_open ( bn_a_i %clname, inum )
+               !CALL iom_open ( bn_a_i%clname, inum )
                !id1 = iom_varid ( bn_a_i%num, bn_a_i%clvar, kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE. )
       	       IF ( zndims == 4 ) THEN
@@ -907,4 +907,2 @@
    !!==============================================================================
 END MODULE bdydta
-
-
Index: trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 4688)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 4689)
@@ -30,5 +30,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  !
-   USE domvvl
+   USE domvvl          ! variable volume
 
    IMPLICIT NONE
Index: trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90	(revision 4688)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90	(revision 4689)
@@ -99,7 +99,4 @@
       REAL(wp) ::   zinda, ztmelts, zdh
 
-      REAL(wp), PARAMETER  ::   zsal = 6.3    ! arbitrary salinity    for incoming ice
-      REAL(wp), PARAMETER  ::   ztem = 270.0  ! arbitrary temperature for incoming ice
-      REAL(wp), PARAMETER  ::   zage = 30.0   ! arbitrary age         for incoming ice
       !!------------------------------------------------------------------------------
       !
@@ -233,13 +230,13 @@
 
                ! Ice salinity, age, temperature
-               sm_i(ji,jj,jl)   = zinda * zsal  + ( 1.0 - zinda ) * s_i_min
-               o_i(ji,jj,jl)    = zinda * zage  + ( 1.0 - zinda )
-               t_su(ji,jj,jl)   = zinda * ztem  + ( 1.0 - zinda ) * ztem
+               sm_i(ji,jj,jl)   = zinda * rn_ice_sal(ib_bdy)  + ( 1.0 - zinda ) * s_i_min
+               o_i(ji,jj,jl)    = zinda * rn_ice_age(ib_bdy)  + ( 1.0 - zinda )
+               t_su(ji,jj,jl)   = zinda * rn_ice_tem(ib_bdy)  + ( 1.0 - zinda ) * rn_ice_tem(ib_bdy)
                DO jk = 1, nlay_s
-                  t_s(ji,jj,jk,jl) = zinda * ztem + ( 1.0 - zinda ) * rtt
+                  t_s(ji,jj,jk,jl) = zinda * rn_ice_tem(ib_bdy) + ( 1.0 - zinda ) * rtt
                END DO
                DO jk = 1, nlay_i
-                  t_i(ji,jj,jk,jl) = zinda * ztem + ( 1.0 - zinda ) * rtt 
-                  s_i(ji,jj,jk,jl) = zinda * zsal + ( 1.0 - zinda ) * s_i_min
+                  t_i(ji,jj,jk,jl) = zinda * rn_ice_tem(ib_bdy) + ( 1.0 - zinda ) * rtt 
+                  s_i(ji,jj,jk,jl) = zinda * rn_ice_sal(ib_bdy) + ( 1.0 - zinda ) * s_i_min
                END DO
                
@@ -259,4 +256,10 @@
 
             END SELECT
+
+            ! if salinity is constant, then overwrite rn_ice_sal
+            IF( num_sal == 1 ) THEN
+               sm_i(ji,jj,jl)   = bulk_sal
+               s_i (ji,jj,:,jl) = bulk_sal
+            ENDIF
 
             ! contents
@@ -338,5 +341,5 @@
       DO ib_bdy=1, nb_bdy
          !
-         SELECT CASE( nn_ice_lim(ib_bdy) )
+         SELECT CASE( cn_ice_lim(ib_bdy) )
 
          CASE('none')
@@ -355,5 +358,5 @@
                   ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd)
                   jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd)
-                  zflag = idx_bdy(ib_bdy)%flagu(jb)
+                  zflag = idx_bdy(ib_bdy)%flagu(jb,jgrd)
                   
                   IF ( ABS( zflag ) == 1. ) THEN  ! eastern and western boundaries
@@ -384,5 +387,5 @@
                   ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd)
                   jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd)
-                  zflag = idx_bdy(ib_bdy)%flagv(jb)
+                  zflag = idx_bdy(ib_bdy)%flagv(jb,jgrd)
                   
                   IF ( ABS( zflag ) == 1. ) THEN  ! northern and southern boundaries
Index: trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 4688)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 4689)
@@ -102,4 +102,7 @@
 #if ( defined key_lim2 || defined key_lim3 )
          &             cn_ice_lim, nn_ice_lim_dta,                           &
+#endif
+#if defined key_lim3
+         &             rn_ice_tem, rn_ice_sal, rn_ice_age,                 &
 #endif
          &             ln_vol, nn_volctl, nn_rimwidth
@@ -359,4 +362,7 @@
         ENDIF
         IF(lwp) WRITE(numout,*)
+        IF(lwp) WRITE(numout,*) '      tem of bdy sea-ice = ', rn_ice_tem(ib_bdy)         
+        IF(lwp) WRITE(numout,*) '      sal of bdy sea-ice = ', rn_ice_sal(ib_bdy)         
+        IF(lwp) WRITE(numout,*) '      age of bdy sea-ice = ', rn_ice_age(ib_bdy)         
 #endif
 
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 4688)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 4689)
@@ -54,6 +54,6 @@
    REAL(wp), PUBLIC ::   r1_rau0                     !: = 1. / rau0                   [m3/kg]
    REAL(wp), PUBLIC ::   rauw     = 1000._wp         !: volumic mass of pure water    [m3/kg]
-   REAL(wp), PUBLIC ::   rcp      =    4.e3_wp       !: ocean specific heat           [J/Kelvin]
-   REAL(wp), PUBLIC ::   r1_rcp                      !: = 1. / rcp                    [Kelvin/J]
+   REAL(wp), PUBLIC ::   rcp      =    4.e3_wp       !: ocean specific heat           [J/kg/K]
+   REAL(wp), PUBLIC ::   r1_rcp                      !: = 1. / rcp                    [kg.K/J]
    REAL(wp), PUBLIC ::   r1_rau0_rcp                 !: = 1. / ( rau0 * rcp )
 
@@ -69,10 +69,10 @@
 #if defined key_lim3 || defined key_cice
    REAL(wp), PUBLIC ::   rhoic    =  917._wp         !: volumic mass of sea ice                               [kg/m3]
-   REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: thermal conductivity of fresh ice
-   REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow
-   REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat for ice 
+   REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: thermal conductivity of fresh ice                     [W/m/K]
+   REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow                          [W/m/K] 
+   REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat for ice                                 [J/kg/K]
    REAL(wp), PUBLIC ::   lsub     =    2.834e+6_wp   !: pure ice latent heat of sublimation                   [J/kg]
    REAL(wp), PUBLIC ::   lfus     =    0.334e+6_wp   !: latent heat of fusion of fresh ice                    [J/kg]
-   REAL(wp), PUBLIC ::   tmut     =    0.054_wp      !: decrease of seawater meltpoint with salinity
+   REAL(wp), PUBLIC ::   tmut     =    0.054_wp      !: decrease of seawater meltpoint with salinity          [degC/ppt]
    REAL(wp), PUBLIC ::   xlsn                        !: = lfus*rhosn (volumetric latent heat fusion of snow)  [J/m3]
 #else
Index: trunk/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 4688)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 4689)
@@ -32,4 +32,5 @@
    USE trc_oce, ONLY :   nn_dttrc        !  !: frequency of step on passive tracers
    USE icb_oce, ONLY :   nclasses, class_num       !  !: iceberg classes
+   USE par_ice
    USE domngb          ! ocean space and time domain
    USE phycst          ! physical constants
@@ -49,5 +50,5 @@
 #endif
    PUBLIC iom_init, iom_swap, iom_open, iom_close, iom_setkt, iom_varid, iom_get, iom_gettime, iom_rstput, iom_put
-   PUBLIC iom_getatt, iom_context_finalize
+   PUBLIC iom_getatt, iom_use, iom_context_finalize
 
    PRIVATE iom_rp0d, iom_rp1d, iom_rp2d, iom_rp3d
@@ -143,4 +144,7 @@
       CALL iom_set_axis_attr( "nfloat", (/ (REAL(ji,wp), ji=1,nfloat) /) )
 # endif
+#if defined key_lim3
+      CALL iom_set_axis_attr( "ncatice", (/ (REAL(ji,wp), ji=1,jpl) /) )
+#endif
       CALL iom_set_axis_attr( "icbcla", class_num )
       
@@ -1015,6 +1019,9 @@
       CHARACTER(LEN=*), INTENT(in) ::   cdname
       REAL(wp)        , INTENT(in) ::   pfield0d
+      REAL(wp)        , DIMENSION(jpi,jpj) ::   zz     ! masson
 #if defined key_iomput
-      CALL xios_send_field(cdname, (/pfield0d/))
+      zz(:,:)=pfield0d
+      CALL xios_send_field(cdname, zz)
+      !CALL xios_send_field(cdname, (/pfield0d/)) 
 #else
       IF( .FALSE. )   WRITE(numout,*) cdname, pfield0d   ! useless test to avoid compilation warnings
@@ -1207,9 +1214,10 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(1,1) ::   zz = 1.
+      REAL(wp), DIMENSION(1) ::   zz = 1.
       !!----------------------------------------------------------------------
       CALL iom_set_domain_attr('scalarpoint', ni_glo=jpnij, nj_glo=1, ibegin=narea, jbegin=1, ni=1, nj=1)
-      CALL iom_set_domain_attr('scalarpoint', data_dim=1)
-      CALL iom_set_domain_attr('scalarpoint', lonvalue=(/ zz /), latvalue=(/ zz /))
+      CALL iom_set_domain_attr('scalarpoint', data_dim=2, data_ibegin = 1, data_ni = 1, data_jbegin = 1, data_nj = 1)
+      zz=REAL(narea,wp)
+      CALL iom_set_domain_attr('scalarpoint', lonvalue=zz, latvalue=zz)
 
    END SUBROUTINE set_scalar
@@ -1499,4 +1507,13 @@
 
 #endif
+
+   LOGICAL FUNCTION iom_use( cdname )
+      CHARACTER(LEN=*), INTENT(in) ::   cdname
+#if defined key_iomput
+      iom_use = xios_field_is_active( cdname )
+#else
+      iom_use = .FALSE.
+#endif
+   END FUNCTION iom_use
    
    !!======================================================================
Index: trunk/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90	(revision 4688)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90	(revision 4689)
@@ -217,5 +217,5 @@
          CALL iom_nf90_check(NF90_Inquire_Variable(if90id, ivarid, dimids = idimid(1:i_nvd)), clinfo)   ! dimensions ids
          iom_file(kiomid)%luld(kiv) = .FALSE.   ! default value
-         iom_file(kiomid)%dimsz(:,kiv) = 0   ! reset dimsz in case previously used
+         iom_file(kiomid)%dimsz(:,kiv) = 0      ! reset dimsz in case previously used
          DO ji = 1, i_nvd                       ! dimensions size
             CALL iom_nf90_check(NF90_Inquire_Dimension(if90id, idimid(ji), len = iom_file(kiomid)%dimsz(ji,kiv)), clinfo)   
Index: trunk/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 4688)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 4689)
@@ -120,4 +120,5 @@
                      CALL iom_rstput( kt, nitrst, numrow, 'hdivb'  , hdivb     )
                      CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb      )
+      IF( lk_vvl )   CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
                      !
                      CALL iom_rstput( kt, nitrst, numrow, 'un'     , un        )     ! now fields
@@ -210,4 +211,5 @@
          CALL iom_get( numror, jpdom_autoglo, 'hdivb'  , hdivb   )
          CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb    )
+         IF( lk_vvl )   CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
       ELSE
          neuler = 0
@@ -245,7 +247,12 @@
          hdivb(:,:,:)   = hdivn(:,:,:)
          sshb (:,:)     = sshn (:,:)
-      ENDIF
-      !
-      IF( lk_lim3 ) THEN 
+         IF( lk_vvl ) THEN
+            DO jk = 1, jpk
+               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
+            END DO
+         ENDIF
+      ENDIF
+      !
+      IF( lk_lim3 ) THEN
          CALL iom_get( numror, jpdom_autoglo, 'iatte' , iatte ) ! clem modif
          CALL iom_get( numror, jpdom_autoglo, 'oatte' , oatte ) ! clem modif
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 4688)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 4689)
@@ -563,4 +563,6 @@
       zcoef_dqsb   = rhoa * cpa * Cice
       zcoef_frca   = 1.0  - 0.3
+      ! MV 2014 the proper cloud fraction (mean summer months from the CLIO climato, NH+SH) is 0.19
+      zcoef_frca   = 1.0  - 0.19
 
 !!gm brutal....
@@ -648,10 +650,5 @@
                p_qsr(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
                ! Long  Wave (lw)
-               ! iovino
-               IF( ff(ji,jj) .GT. 0._wp ) THEN
-                  z_qlw(ji,jj,jl) = ( 0.95 * sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
-               ELSE
-                  z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
-               ENDIF
+               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
                ! lw sensitivity
                z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3                                               
@@ -668,5 +665,11 @@
                   &                         * (  11637800. * EXP( -5897.8 / pst(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) )
                ! Latent heat sensitivity for ice (Dqla/Dt)
-               p_dqla(ji,jj,jl) = rn_efac * zcoef_dqla * z_wnds_t(ji,jj) / ( zst2 ) * EXP( -5897.8 / pst(ji,jj,jl) )
+               ! MV we also have to cap the sensitivity if the flux is zero
+               IF ( p_qla(ji,jj,jl) .GT. 0.0 ) THEN
+                  p_dqla(ji,jj,jl) = rn_efac * zcoef_dqla * z_wnds_t(ji,jj) / ( zst2 ) * EXP( -5897.8 / pst(ji,jj,jl) )
+               ELSE
+                  p_dqla(ji,jj,jl) = 0.0
+               ENDIF
+                             
                ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
                z_dqsb(ji,jj,jl) = zcoef_dqsb * z_wnds_t(ji,jj)
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 4688)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 4689)
@@ -57,4 +57,6 @@
       !!                =1 global mean of emp set to zero at each nn_fsbc time step
       !!                =2 annual global mean corrected from previous year
+      !!                =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 
       !!----------------------------------------------------------------------
@@ -81,5 +83,8 @@
             IF( kn_fwb == 1 )   WRITE(numout,*) '          instantaneously set to zero'
             IF( kn_fwb == 2 )   WRITE(numout,*) '          adjusted from previous year budget'
-         ENDIF
+            IF( kn_fwb == 3 )   WRITE(numout,*) '          fwf set to zero and spread out over erp area'
+         ENDIF
+         !
+         IF( kn_fwb == 3 .AND. nn_sssr /= 2 )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 requires nn_sssr = 2, we stop ' )
          !
          area = glob_sum( e1e2t(:,:) )           ! interior global domain surface
@@ -142,6 +147,58 @@
          ENDIF
          !
+      CASE ( 3 )                             !==  global fwf set to zero and spread out over erp area  ==!
+         !
+         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
+            ztmsk_pos(:,:) = tmask_i(:,:)                      ! Select <0 and >0 area of erp
+            WHERE( erp < 0._wp )   ztmsk_pos = 0._wp
+            ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:)
+            !
+            zsurf_neg = glob_sum( e1e2t(:,:)*ztmsk_neg(:,:) )  ! Area filled by <0 and >0 erp 
+            zsurf_pos = glob_sum( e1e2t(:,:)*ztmsk_pos(:,:) )
+            !                                                  ! 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
+                zsurf_tospread      = zsurf_pos
+                ztmsk_tospread(:,:) = ztmsk_pos(:,:)
+            ELSE                             ! spread out over <0 erp area to increase precipitation
+                zsurf_tospread      = zsurf_neg
+                ztmsk_tospread(:,:) = ztmsk_neg(:,:)
+            ENDIF
+            !
+            zsum_fwf   = glob_sum( e1e2t(:,:) * z_fwf )         ! fwf global mean over <0 or >0 erp area
+!!gm :  zsum_fwf   = z_fwf * area   ???  it is right?  I think so....
+            z_fwf_nsrf =  zsum_fwf / ( zsurf_tospread + rsmall )
+            !                                                  ! weight to respect erp field 2D structure 
+            zsum_erp   = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e1e2t(:,:) )
+            z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( zsum_erp + rsmall )
+            !                                                  ! final correction term to apply
+            zerp_cor(:,:) = -1. * z_fwf_nsrf * zsurf_tospread * z_wgt(:,:)
+            !
+!!gm   ===>>>>  lbc_lnk should be useless as all the computation is done over the whole domain !
+            CALL lbc_lnk( zerp_cor, 'T', 1. )
+            !
+            emp(:,:) = emp(:,:) + zerp_cor(:,:)
+            qns(:,:) = qns(:,:) - zerp_cor(:,:) * rcp * sst_m(:,:)  ! account for change to the heat budget due to fw correction
+            erp(:,:) = erp(:,:) + zerp_cor(:,:)
+            !
+            IF( nprint == 1 .AND. lwp ) THEN                   ! control print
+               IF( z_fwf < 0._wp ) THEN
+                  WRITE(numout,*)'   z_fwf < 0'
+                  WRITE(numout,*)'   SUM(erp+)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
+               ELSE
+                  WRITE(numout,*)'   z_fwf >= 0'
+                  WRITE(numout,*)'   SUM(erp-)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
+               ENDIF
+               WRITE(numout,*)'   SUM(empG)     = ', SUM( z_fwf*e1e2t(:,:) )*1.e-9,' Sv'
+               WRITE(numout,*)'   z_fwf         = ', z_fwf      ,' Kg/m2/s'
+               WRITE(numout,*)'   z_fwf_nsrf    = ', z_fwf_nsrf ,' Kg/m2/s'
+               WRITE(numout,*)'   MIN(zerp_cor) = ', MINVAL(zerp_cor) 
+               WRITE(numout,*)'   MAX(zerp_cor) = ', MAXVAL(zerp_cor) 
+            ENDIF
+         ENDIF
+         !
       CASE DEFAULT                           !==  you should never be there  ==!
-         CALL ctl_stop( 'sbc_fwb : wrong nn_fwb value for the FreshWater Budget correction, choose either 1 or 2' )
+         CALL ctl_stop( 'sbc_fwb : wrong nn_fwb value for the FreshWater Budget correction, choose either 1, 2 or 3' )
          !
       END SELECT
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 4688)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 4689)
@@ -59,4 +59,5 @@
    USE prtctl          ! Print control
    USE lib_fortran     ! 
+   USE cpl_oasis3, ONLY : lk_cpl
 
 #if defined key_bdy 
@@ -68,4 +69,5 @@
 
    PUBLIC sbc_ice_lim  ! routine called by sbcmod.F90
+   PUBLIC lim_prt_state
    
    !! * Substitutions
@@ -133,5 +135,5 @@
       INTEGER, INTENT(in) ::   kblk    ! type of bulk (=3 CLIO, =4 CORE)
       !!
-      INTEGER  ::   jl      ! dummy loop index
+      INTEGER  ::   ji, jj, jl, jk      ! dummy loop index
       REAL(wp) ::   zcoef   ! local scalar
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice_os, zalb_ice_cs  ! albedo of the ice under overcast/clear sky
@@ -146,4 +148,6 @@
       REAL(wp), POINTER, DIMENSION(:,:) :: z_dqns_ice_all  ! Mean d(qns)/dT over all categories
       REAL(wp), POINTER, DIMENSION(:,:) :: z_dqla_ice_all  ! Mean d(qla)/dT over all categories
+      REAL(wp) ::   ztmelts           ! clem 2014: for HC diags
+      REAL(wp) ::   epsi20 = 1.e-20   !
       !!----------------------------------------------------------------------
 
@@ -152,11 +156,10 @@
       IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim')
 
-      CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs )
-
-#if defined key_coupled
-      IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_alloc( jpi,jpj,jpl, zalb_ice)
-      IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) &
-         &   CALL wrk_alloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all)
-#endif
+      CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs, zalb_ice )
+
+      IF( lk_cpl ) THEN
+         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) &
+            &   CALL wrk_alloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all)
+      ENDIF
 
       IF( kt == nit000 ) THEN
@@ -168,5 +171,5 @@
          !
          IF( ln_nicep ) THEN      ! control print at a given point
-            jiindx = 177   ;   jjindx = 112
+            jiindx = 15    ;   jjindx =  44
             IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx
          ENDIF
@@ -176,12 +179,13 @@
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only  !
          !                                     !----------------------!
-         !                                           !  Bulk Formulea !
+         !                                           !  Bulk Formulae !
          !                                           !----------------!
          !
          u_oce(:,:) = ssu_m(:,:)                     ! mean surface ocean current at ice velocity point
          v_oce(:,:) = ssv_m(:,:)                     ! (C-grid dynamics :  U- & V-points as the ocean)
-         !
-         t_bo(:,:) = tfreez( sss_m ) +  rt0          ! masked sea surface freezing temperature [Kelvin]
-         !                                           ! (set to rt0 over land)
+
+         ! masked sea surface freezing temperature [Kelvin]
+         t_bo(:,:) = ( tfreez( sss_m ) +  rt0 ) * tmask(:,:,1) + rt0 * ( 1. - tmask(:,:,1) )
+
          CALL albedo_ice( t_su, ht_i, ht_s, zalb_ice_cs, zalb_ice_os )  ! ... ice albedo
 
@@ -192,13 +196,13 @@
          IF ( ln_cpl ) zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) )
          
-#if defined key_coupled
-         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
-            !
-            ! Compute mean albedo and temperature
-            zalb_ice_all (:,:) = fice_ice_ave ( zalb_ice (:,:,:) ) 
-            ztem_ice_all (:,:) = fice_ice_ave ( tn_ice   (:,:,:) ) 
-            !
+         IF( lk_cpl ) THEN
+            IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
+               !
+               ! Compute mean albedo and temperature
+               zalb_ice_all (:,:) = fice_ice_ave ( zalb_ice (:,:,:) ) 
+               ztem_ice_all (:,:) = fice_ice_ave ( tn_ice   (:,:,:) ) 
+               !
+            ENDIF
          ENDIF
-#endif
                                                ! Bulk formulea - provides the following fields:
          ! utau_ice, vtau_ice : surface ice stress                     (U- & V-points)   [N/m2]
@@ -218,5 +222,13 @@
             !         
          CASE( 4 )                                       ! CORE bulk formulation
-            CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice_cs,               &
+            ! MV 2014
+            ! We must account for cloud fraction in the computation of the albedo
+            ! The present ref just uses the clear sky value
+            ! The overcast sky value is 0.06 higher, and polar skies are mostly overcast
+            ! CORE has no cloud fraction, hence we must prescribe it
+            ! Mean summer cloud fraction computed from CLIO = 0.81
+            zalb_ice(:,:,:) = 0.19 * zalb_ice_cs(:,:,:) + 0.81 * zalb_ice_os(:,:,:)
+            ! Following line, we replace zalb_ice_cs by simply zalb_ice
+            CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice   ,               &
                &                      utau_ice  , vtau_ice  , qns_ice    , qsr_ice   ,   &
                &                      qla_ice   , dqns_ice  , dqla_ice   ,               &
@@ -239,5 +251,5 @@
 
          ! Average over all categories
-#if defined key_coupled
+         IF( lk_cpl ) THEN
          IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
 
@@ -269,5 +281,5 @@
             END IF
          END IF
-#endif
+         ENDIF
          !                                           !----------------------!
          !                                           ! LIM-3  time-stepping !
@@ -285,8 +297,8 @@
          old_smv_i(:,:,:)   = smv_i(:,:,:)     ! salt content
          old_oa_i (:,:,:)   = oa_i (:,:,:)     ! areal age content
-         !
-         old_u_ice(:,:) = u_ice(:,:)
-         old_v_ice(:,:) = v_ice(:,:)
-         !                                           ! intialisation to zero    !!gm is it truly necessary ???
+         old_u_ice(:,:)     = u_ice(:,:)
+         old_v_ice(:,:)     = v_ice(:,:)
+
+         ! trends    !!gm is it truly necessary ???
          d_a_i_thd  (:,:,:)   = 0._wp   ;   d_a_i_trp  (:,:,:)   = 0._wp
          d_v_i_thd  (:,:,:)   = 0._wp   ;   d_v_i_trp  (:,:,:)   = 0._wp
@@ -296,33 +308,40 @@
          d_smv_i_thd(:,:,:)   = 0._wp   ;   d_smv_i_trp(:,:,:)   = 0._wp
          d_oa_i_thd (:,:,:)   = 0._wp   ;   d_oa_i_trp (:,:,:)   = 0._wp
-         !
-         d_u_ice_dyn(:,:) = 0._wp
-         d_v_ice_dyn(:,:) = 0._wp
-         !
-         sfx    (:,:) = 0._wp   ;   sfx_thd  (:,:) = 0._wp
-         sfx_bri(:,:) = 0._wp   ;   sfx_mec  (:,:) = 0._wp   ;   sfx_res  (:,:) = 0._wp
-         fhbri  (:,:) = 0._wp   ;   fheat_mec(:,:) = 0._wp   ;   fheat_res(:,:) = 0._wp
-         fhmec  (:,:) = 0._wp   ;   
-         fmmec  (:,:) = 0._wp
-         fmmflx (:,:) = 0._wp     
-         focea2D(:,:) = 0._wp
-         fsup2D (:,:) = 0._wp
-
-         ! used in limthd.F90
-         rdvosif(:,:) = 0._wp   ! variation of ice volume at surface
-         rdvobif(:,:) = 0._wp   ! variation of ice volume at bottom
-         fdvolif(:,:) = 0._wp   ! total variation of ice volume
-         rdvonif(:,:) = 0._wp   ! lateral variation of ice volume
-         fstric (:,:) = 0._wp   ! part of solar radiation transmitted through the ice
-         ffltbif(:,:) = 0._wp   ! linked with fstric
-         qfvbq  (:,:) = 0._wp   ! linked with fstric
-         rdm_snw(:,:) = 0._wp   ! variation of snow mass per unit area
-         rdm_ice(:,:) = 0._wp   ! variation of ice mass per unit area
-         hicifp (:,:) = 0._wp   ! daily thermodynamic ice production. 
-         !
-         diag_sni_gr(:,:) = 0._wp   ;   diag_lat_gr(:,:) = 0._wp
-         diag_bot_gr(:,:) = 0._wp   ;   diag_dyn_gr(:,:) = 0._wp
-         diag_bot_me(:,:) = 0._wp   ;   diag_sur_me(:,:) = 0._wp
-         diag_res_pr(:,:) = 0._wp   ;   diag_trp_vi(:,:) = 0._wp
+         d_u_ice_dyn(:,:)     = 0._wp   ;   d_v_ice_dyn(:,:)     = 0._wp
+
+         ! salt, heat and mass fluxes
+         sfx    (:,:) = 0._wp   ;
+         sfx_bri(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
+         sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp
+         sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp
+         sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp
+         sfx_res(:,:) = 0._wp
+
+         wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp
+         wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp
+         wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp
+         wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp
+         wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp
+         wfx_spr(:,:) = 0._wp   ;   
+
+         hfx_in (:,:) = 0._wp   ;   hfx_out(:,:) = 0._wp
+         hfx_thd(:,:) = 0._wp   ;   
+         hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp
+         hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp
+         hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp
+         hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp
+         hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp 
+         hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp
+
+         !
+         fhld  (:,:)    = 0._wp 
+         fmmflx(:,:)    = 0._wp     
+         ! part of solar radiation transmitted through the ice
+         ftr_ice(:,:,:) = 0._wp
+
+         ! diags
+         diag_trp_vi  (:,:) = 0._wp  ; diag_trp_vs(:,:) = 0._wp  ;  diag_trp_ei(:,:) = 0._wp  ;  diag_trp_es(:,:) = 0._wp
+         diag_heat_dhc(:,:) = 0._wp  
+
          ! dynamical invariants
          delta_i(:,:) = 0._wp       ;   divu_i(:,:) = 0._wp       ;   shear_i(:,:) = 0._wp
@@ -375,5 +394,4 @@
                           zcoef = rdt_ice /rday           !  Ice natural aging
                           oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef
-                          CALL lim_var_glo2eqv            ! this CALL is maybe not necessary (Martin)
          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print
                           CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion  !
@@ -391,7 +409,7 @@
          !                                           ! Diagnostics and outputs 
          IF (ln_limdiaout) CALL lim_diahsb
-!clem # if ! defined key_iomput
+
                           CALL lim_wri( 1  )              ! Ice outputs 
-!clem # endif
+
          IF( kt == nit000 .AND. ln_rstart )   &
             &             CALL iom_close( numrir )        ! clem: close input ice restart file
@@ -413,12 +431,10 @@
       
 !!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
-      !
-      CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs )
-
-#if defined key_coupled
-      IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice)
-      IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) &
-         &    CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all)
-#endif
+      CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs, zalb_ice )
+
+      IF( lk_cpl ) THEN
+         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) &
+            &    CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all)
+      ENDIF
       !
       IF( nn_timing == 1 )  CALL timing_stop('sbc_ice_lim')
@@ -534,5 +550,4 @@
 !                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj)
 !                 WRITE(numout,*) ' sss                  : ', sss_m(ji,jj)
-!                 WRITE(numout,*) ' s_i_newice           : ', s_i_newice(ji,jj,1:jpl)
 !                 WRITE(numout,*) 
                   inb_alp(ialert_id) = inb_alp(ialert_id) + 1
@@ -591,17 +606,4 @@
                !WRITE(numout,*) ' sst       : ', sst_m(ji,jj)
                !WRITE(numout,*) ' sss       : ', sss_m(ji,jj)
-               !WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj)
-               !WRITE(numout,*) ' qldif     : ', qldif(ji,jj)
-               !WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) / rdt_ice
-               !WRITE(numout,*) ' qldif     : ', qldif(ji,jj) / rdt_ice
-               !WRITE(numout,*) ' qfvbq     : ', qfvbq(ji,jj)
-               !WRITE(numout,*) ' qdtcn     : ', qdtcn(ji,jj)
-               !WRITE(numout,*) ' qfvbq / dt: ', qfvbq(ji,jj) / rdt_ice
-               !WRITE(numout,*) ' qdtcn / dt: ', qdtcn(ji,jj) / rdt_ice
-               !WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj) 
-               !WRITE(numout,*) ' fhmec     : ', fhmec(ji,jj) 
-               !WRITE(numout,*) ' fheat_mec : ', fheat_mec(ji,jj) 
-               !WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj) 
-               !WRITE(numout,*) ' fhbri     : ', fhbri(ji,jj) 
                !
                !CALL lim_prt_state( kt, ji, jj, 2, '   ')
@@ -790,14 +792,8 @@
                WRITE(numout,*) ' - Heat / FW fluxes '
                WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ '
-               WRITE(numout,*) ' emp        : ', emp      (ji,jj)
-               WRITE(numout,*) ' sfx        : ', sfx      (ji,jj)
-               WRITE(numout,*) ' sfx_thd    : ', sfx_thd(ji,jj)
-               WRITE(numout,*) ' sfx_bri    : ', sfx_bri  (ji,jj)
-               WRITE(numout,*) ' sfx_mec    : ', sfx_mec  (ji,jj)
-               WRITE(numout,*) ' sfx_res    : ', sfx_res(ji,jj)
-               WRITE(numout,*) ' fmmec      : ', fmmec    (ji,jj)
-               WRITE(numout,*) ' fhmec      : ', fhmec    (ji,jj)
-               WRITE(numout,*) ' fhbri      : ', fhbri    (ji,jj)
-               WRITE(numout,*) ' fheat_mec  : ', fheat_mec(ji,jj)
+               WRITE(numout,*) ' - Heat fluxes in and out the ice ***'
+               WRITE(numout,*) ' qsr_ini       : ', pfrld(ji,jj) * qsr(ji,jj) + SUM( old_a_i(ji,jj,:) * qsr_ice(ji,jj,:) )
+               WRITE(numout,*) ' qns_ini       : ', pfrld(ji,jj) * qns(ji,jj) + SUM( old_a_i(ji,jj,:) * qns_ice(ji,jj,:) )
+               WRITE(numout,*)
                WRITE(numout,*) 
                WRITE(numout,*) ' sst        : ', sst_m(ji,jj)  
@@ -829,21 +825,27 @@
                WRITE(numout,*) ' qsr       : ', qsr(ji,jj)
                WRITE(numout,*) ' qns       : ', qns(ji,jj)
-               WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj)
-               WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) * r1_rdtice
-               WRITE(numout,*) ' qldif     : ', qldif(ji,jj) * r1_rdtice
+               WRITE(numout,*)
+               WRITE(numout,*) ' hfx_mass     : ', hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_snw(ji,jj) + hfx_res(ji,jj)
+               WRITE(numout,*) ' hfx_in       : ', hfx_in(ji,jj)
+               WRITE(numout,*) ' hfx_out      : ', hfx_out(ji,jj)
+               WRITE(numout,*) ' dhc          : ', diag_heat_dhc(ji,jj)              
+               WRITE(numout,*)
+               WRITE(numout,*) ' hfx_dyn      : ', hfx_dyn(ji,jj)
+               WRITE(numout,*) ' hfx_thd      : ', hfx_thd(ji,jj)
+               WRITE(numout,*) ' hfx_res      : ', hfx_res(ji,jj)
+               WRITE(numout,*) ' fhtur        : ', fhtur(ji,jj) 
+               WRITE(numout,*) ' qlead        : ', qlead(ji,jj) * r1_rdtice
                WRITE(numout,*)
                WRITE(numout,*) ' - Salt fluxes at bottom interface ***'
                WRITE(numout,*) ' emp       : ', emp    (ji,jj)
-               WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ji,jj)
                WRITE(numout,*) ' sfx       : ', sfx    (ji,jj)
                WRITE(numout,*) ' sfx_res   : ', sfx_res(ji,jj)
-               WRITE(numout,*) ' sfx_mec   : ', sfx_mec(ji,jj)
-               WRITE(numout,*) ' - Heat fluxes at bottom interface ***'
-               WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj)
+               WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ji,jj)
+               WRITE(numout,*) ' sfx_dyn   : ', sfx_dyn(ji,jj)
                WRITE(numout,*)
                WRITE(numout,*) ' - Momentum fluxes '
                WRITE(numout,*) ' utau      : ', utau(ji,jj) 
                WRITE(numout,*) ' vtau      : ', vtau(ji,jj)
-            ENDIF
+            ENDIF 
             WRITE(numout,*) ' '
             !
