Index: branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 9867)
+++ branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 9987)
@@ -94,5 +94,5 @@
       !!              - fr_i    : ice fraction
       !!              - tn_ice  : sea-ice surface temperature
-      !!              - alb_ice : sea-ice albedo (only useful in coupled mode)
+      !!              - alb_ice : sea-ice albedo (recomputed only for coupled mode)
       !!
       !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
@@ -106,8 +106,10 @@
       REAL(wp) ::   zqsr                                           ! New solar flux received by the ocean
       !
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb_cs, zalb_os     ! 2D/3D workspace
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb_cs, zalb_os     ! 3D workspace
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zalb                 ! 2D workspace
       !!---------------------------------------------------------------------
 
       ! make calls for heat fluxes before it is modified
+      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
       IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) )                                   !     solar flux at ocean surface
       IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) )                   ! non-solar flux at ocean surface
@@ -118,8 +120,24 @@
       IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   &
          &                                                      * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) )
-      IF( iom_use('qemp_oce' ) )   CALL iom_put( "qemp_oce"  , qemp_oce(:,:) )  
-      IF( iom_use('qemp_ice' ) )   CALL iom_put( "qemp_ice"  , qemp_ice(:,:) )  
-
-      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
+      IF( iom_use('qemp_oce') )  CALL iom_put( "qemp_oce" , qemp_oce(:,:) )  
+      IF( iom_use('qemp_ice') )  CALL iom_put( "qemp_ice" , qemp_ice(:,:) )  
+      IF( iom_use('emp_oce' ) )  CALL iom_put( "emp_oce"  , emp_oce(:,:) )   ! emp over ocean (taking into account the snow blown away from the ice)
+      IF( iom_use('emp_ice' ) )  CALL iom_put( "emp_ice"  , emp_ice(:,:) )   ! emp over ice   (taking into account the snow blown away from the ice)
+
+      ! albedo output
+      CALL wrk_alloc( jpi,jpj, zalb )    
+
+      zalb(:,:) = 0._wp
+      WHERE     ( SUM( a_i_b, dim=3 ) <= epsi06 )  ;  zalb(:,:) = 0.066_wp
+      ELSEWHERE                                    ;  zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 )
+      END WHERE
+      IF( iom_use('alb_ice' ) )  CALL iom_put( "alb_ice"  , zalb(:,:) )           ! ice albedo output
+
+      zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + 0.066_wp * ( 1._wp - SUM( a_i_b, dim=3 ) )      
+      IF( iom_use('albedo'  ) )  CALL iom_put( "albedo"  , zalb(:,:) )           ! ice albedo output
+
+      CALL wrk_dealloc( jpi,jpj, zalb )    
+      !
+      
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -140,10 +158,11 @@
             hfx_out(ji,jj) = hfx_out(ji,jj) + zqmass + zqsr
 
-            ! Add the residual from heat diffusion equation (W.m-2)
-            !-------------------------------------------------------
-            hfx_out(ji,jj) = hfx_out(ji,jj) + hfx_err_dif(ji,jj)
+            ! Add the residual from heat diffusion equation and sublimation (W.m-2)
+            !----------------------------------------------------------------------
+            hfx_out(ji,jj) = hfx_out(ji,jj) + hfx_err_dif(ji,jj) +   &
+               &           ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) )
 
             ! New qsr and qns used to compute the oceanic heat flux at the next time step
-            !---------------------------------------------------
+            !----------------------------------------------------------------------------
             qsr(ji,jj) = zqsr                                      
             qns(ji,jj) = hfx_out(ji,jj) - zqsr              
@@ -165,7 +184,6 @@
 
             ! mass flux at the ocean/ice interface
-            fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) ) * r1_rdtice  ! F/M mass flux save at least for biogeochemical model
-            emp(ji,jj)    = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange)
-            
+            fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_err_sub(ji,jj) )              ! F/M mass flux save at least for biogeochemical model
+            emp(ji,jj)    = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange)
          END DO
       END DO
@@ -175,5 +193,5 @@
       !------------------------------------------!
       sfx(:,:) = sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:)   &
-         &     + sfx_res(:,:) + sfx_dyn(:,:) + sfx_bri(:,:)
+         &     + sfx_res(:,:) + sfx_dyn(:,:) + sfx_bri(:,:) + sfx_sub(:,:)
 
       !-------------------------------------------------------------!
