Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 3605)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 3625)
@@ -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
@@ -64,8 +65,9 @@
       INTEGER, INTENT( in ) ::   kn_fwb   ! ocean time-step index
       !
-      INTEGER  ::   inum, ikty, iyear   ! local integers
-      REAL(wp) ::   z_fwf, z_fwf_nsrf, zsum_fwf, zsum_erp   ! local scalars
-      REAL(wp) ::   zsurf_neg, zsurf_pos, zsurf_tospread    !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:) ::   ztmsk_neg, ztmsk_pos, ztmsk_tospread, z_wgt, zerp_cor
+      INTEGER  ::   inum, ikty, iyear     ! local integers
+      REAL(wp) ::   z_fwf, z_fwf_nsrf, zsum_fwf, zsum_erp                ! local scalars
+      REAL(wp) ::   zsurf_neg, zsurf_pos, zsurf_tospread, zcoef          !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztmsk_neg, ztmsk_pos, z_wgt ! 2D workspaces
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztmsk_tospread, zerp_cor    !   -      -
       !!----------------------------------------------------------------------
       !
@@ -87,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
       
@@ -95,12 +103,13 @@
          !
          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
-            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) ) ) / area   ! sum over the global domain
-            emp (:,:) = emp (:,:) - z_fwf 
-            emps(:,:) = emps(:,:) - z_fwf 
+            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) -  snwice_fmass(:,:) ) ) / area   ! sum over the global domain
+            zcoef = z_fwf * rcp
+            emp(:,:) = emp(:,:) - z_fwf 
+            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
+         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. )
@@ -117,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 
@@ -125,10 +135,11 @@
          ENDIF
          ! 
-         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN      ! correct the freshwater fluxes
-            emp (:,:) = emp (:,:) + fwfold
-            emps(:,:) = emps(:,:) + fwfold
-         ENDIF
-         !
-         IF( kt == nitend .AND. lwp ) THEN         ! save fwfold value in a file
+         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN         ! correct the freshwater fluxes
+            zcoef = fwfold * rcp
+            emp(:,:) = emp(:,:) + fwfold
+            qns(:,:) = qns(:,:) - zcoef * sst_m(:,:)  ! account for change to the heat budget due to fw correction
+         ENDIF
+         !
+         IF( kt == nitend .AND. lwp ) THEN            ! save fwfold value in a file
             CALL ctl_opn( inum, 'EMPave.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
             WRITE( inum, "(24X,I8,2ES24.16)" ) nyear, a_fwb_b, a_fwb
@@ -143,8 +154,8 @@
             ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:)
             !
-            zsurf_neg = glob_sum( e1e2t(:,:)*ztmsk_neg(:,:) )   ! Area filled by <0 and >0 erp 
+            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
@@ -160,5 +171,5 @@
             z_fwf_nsrf =  zsum_fwf / ( zsurf_tospread + rsmall )
             !                                                  ! weight to respect erp field 2D structure 
-            zsum_erp = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e1e2t(:,:) )
+            zsum_erp   = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e1e2t(:,:) )
             z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( zsum_erp + rsmall )
             !                                                  ! final correction term to apply
@@ -168,7 +179,7 @@
             CALL lbc_lnk( zerp_cor, 'T', 1. )
             !
-            emp (:,:) = emp (:,:) + zerp_cor(:,:)
-            emps(:,:) = emps(:,:) + zerp_cor(:,:)
-            erp (:,:) = erp (:,:) + zerp_cor(:,:)
+            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
