Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 6398)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 6399)
@@ -253,20 +253,20 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhld        !: heat flux from the lead used for bottom melting
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw    !: snow-ocean mass exchange over 1 time step [kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_spr    !: snow precipitation on ice over 1 time step [kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sub    !: snow sublimation over 1 time step [kg/m2]
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice    !: ice-ocean mass exchange over 1 time step [kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sni    !: snow ice growth component of wfx_ice [kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_opw    !: lateral ice growth component of wfx_ice [kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bog    !: bottom ice growth component of wfx_ice [kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_dyn    !: dynamical ice growth component of wfx_ice [kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bom    !: bottom melt component of wfx_ice [kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sum    !: surface melt component of wfx_ice [kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_res    !: residual component of wfx_ice [kg/m2]
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_tot     !: ice concentration tendency (total) [s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw    !: snow-ocean mass exchange   [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_spr    !: snow precipitation on ice  [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sub    !: snow/ice sublimation       [kg.m-2.s-1]
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice    !: ice-ocean mass exchange                   [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sni    !: snow ice growth component of wfx_ice      [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_opw    !: lateral ice growth component of wfx_ice   [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bog    !: bottom ice growth component of wfx_ice    [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_dyn    !: dynamical ice growth component of wfx_ice [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bom    !: bottom melt component of wfx_ice          [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sum    !: surface melt component of wfx_ice         [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_res    !: residual component of wfx_ice             [kg.m-2.s-1]
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_tot     !: ice concentration tendency (total)          [s-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_thd     !: ice concentration tendency (thermodynamics) [s-1]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_dyn     !: ice concentration tendency (dynamics) [s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_dyn     !: ice concentration tendency (dynamics)       [s-1]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bog     !: salt flux due to ice growth/melt                      [PSU/m2/s]
@@ -279,24 +279,27 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_res     !: residual salt flux due to correction of ice thickness [PSU/m2/s]
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bog     !: total heat flux causing bottom ice growth 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bom     !: total heat flux causing bottom ice melt 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sum     !: total heat flux causing surface ice melt 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_opw     !: total heat flux causing open water ice formation
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dif     !: total heat flux causing Temp change in the ice 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_snw     !: heat flux for snow melt 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err     !: heat flux error after heat diffusion 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_dif !: heat flux remaining due to change in non-solar flux
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_rem !: heat flux error after heat remapping 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_in      !: heat flux available for thermo transformations 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_out     !: heat flux remaining at the end of thermo transformations 
-
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_sub     !: salt flux due to ice sublimation
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bog     !: total heat flux causing bottom ice growth        [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bom     !: total heat flux causing bottom ice melt          [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sum     !: total heat flux causing surface ice melt         [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_opw     !: total heat flux causing open water ice formation [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dif     !: total heat flux causing Temp change in the ice   [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_snw     !: heat flux for snow melt                          [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err     !: heat flux error after heat diffusion             [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_dif !: heat flux remaining due to change in non-solar flux [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_rem !: heat flux error after heat remapping             [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_in      !: heat flux available for thermo transformations   [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_out     !: heat flux remaining at the end of thermo transformations  [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_err_sub !: mass flux error after sublimation [kg.m-2.s-1]
+   
    ! heat flux associated with ice-atmosphere mass exchange
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sub     !: heat flux for sublimation 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_spr     !: heat flux of the snow precipitation 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sub     !: heat flux for sublimation  [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_spr     !: heat flux of the snow precipitation  [W.m-2]
 
    ! heat flux associated with ice-ocean mass exchange
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_thd     !: ice-ocean heat flux from thermo processes (limthd_dh) 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dyn     !: ice-ocean heat flux from mecanical processes (limitd_me) 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_res     !: residual heat flux due to correction of ice thickness
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_thd     !: ice-ocean heat flux from thermo processes (limthd_dh)  [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dyn     !: ice-ocean heat flux from mecanical processes (limitd_me)  [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_res     !: residual heat flux due to correction of ice thickness [W.m-2]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice   !: transmitted solar radiation under ice
@@ -441,8 +444,8 @@
          &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl), qlead  (jpi,jpj) ,                     &
          &      rn_amax_2d(jpi,jpj),                                                            &
-         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) ,                        &
+         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) ,                       &
          &      sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,    &
          &      hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj) ,     & 
-         &      hfx_err_dif(jpi,jpj) , hfx_err_rem(jpi,jpj) ,                                   &
+         &      hfx_err_dif(jpi,jpj) , hfx_err_rem(jpi,jpj) , wfx_err_sub(jpi,jpj) ,       &
          &      hfx_in (jpi,jpj) , hfx_out(jpi,jpj) , fhld(jpi,jpj) ,                           &
          &      hfx_sum(jpi,jpj) , hfx_bom(jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) , hfx_opw(jpi,jpj) ,    &
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 6398)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 6399)
@@ -24,5 +24,6 @@
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
    USE sbc_oce , ONLY : sfx  ! Surface boundary condition: ocean fields
-
+   USE sbc_ice , ONLY : qevap_ice
+   
    IMPLICIT NONE
    PRIVATE
@@ -184,5 +185,5 @@
          ! salt flux
          zfs_b  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
-            &                  sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  &
+            &                  sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:)                   &
             &                ) *  e12t(:,:) * tmask(:,:,1) * zconv )
 
@@ -209,5 +210,5 @@
          ! salt flux
          zfs  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
-            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  & 
+            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:)                   & 
             &              ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zfs_b
 
@@ -287,5 +288,6 @@
 #if ! defined key_bdy
       ! heat flux
-      zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - hfx_sub ) * e12t * tmask(:,:,1) * zconv ) 
+      zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - SUM( qevap_ice * a_i_b, dim=3 ) )  &
+         &              * e12t * tmask(:,:,1) * zconv ) 
       ! salt flux
       zsfx  = glob_sum( ( sfx + diag_smvi ) * e12t * tmask(:,:,1) * zconv ) * rday
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 6398)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 6399)
@@ -125,5 +125,5 @@
       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)
 
-      ! clem 2016: albedo output
+      ! albedo output
       CALL wrk_alloc( jpi,jpj, zalb )    
 
@@ -158,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              
@@ -183,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
@@ -193,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(:,:)
 
       !-------------------------------------------------------------!
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 6398)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 6399)
@@ -461,5 +461,5 @@
 
       DO ji = kideb, kiut
-         zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) )
+         zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) + dh_i_sub(ji) )
          IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp )  THEN
             zvi          = a_i_1d(ji) * ht_i_1d(ji)
@@ -470,5 +470,5 @@
             zda_mel     = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) )
             a_i_1d(ji)  = MAX( epsi20, a_i_1d(ji) + zda_mel ) 
-             ! adjust thickness
+            ! adjust thickness
             ht_i_1d(ji) = zvi / a_i_1d(ji)            
             ht_s_1d(ji) = zvs / a_i_1d(ji)            
@@ -514,4 +514,5 @@
          
          CALL tab_2d_1d( nbpb, qprec_ice_1d(1:nbpb), qprec_ice(:,:) , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d( nbpb, qevap_ice_1d(1:nbpb), qevap_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, fr1_i0_1d  (1:nbpb), fr1_i0          , jpi, jpj, npb(1:nbpb) )
@@ -543,5 +544,6 @@
          CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri         , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res         , jpi, jpj, npb(1:nbpb) )
-         
+         CALL tab_2d_1d( nbpb, sfx_sub_1d (1:nbpb), sfx_sub         , jpi, jpj,npb(1:nbpb) )
+ 
          CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd         , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr         , jpi, jpj, npb(1:nbpb) )
@@ -593,5 +595,6 @@
          CALL tab_1d_2d( nbpb, sfx_res       , npb, sfx_res_1d(1:nbpb)   , jpi, jpj )
          CALL tab_1d_2d( nbpb, sfx_bri       , npb, sfx_bri_1d(1:nbpb)   , jpi, jpj )
-         
+         CALL tab_1d_2d( nbpb, sfx_sub       , npb, sfx_sub_1d(1:nbpb)   , jpi, jpj )        
+ 
          CALL tab_1d_2d( nbpb, hfx_thd       , npb, hfx_thd_1d(1:nbpb)   , jpi, jpj )
          CALL tab_1d_2d( nbpb, hfx_spr       , npb, hfx_spr_1d(1:nbpb)   , jpi, jpj )
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 6398)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 6399)
@@ -74,5 +74,5 @@
 
       REAL(wp) ::   ztmelts             ! local scalar
-      REAL(wp) ::   zfdum       
+      REAL(wp) ::   zdum       
       REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
       REAL(wp) ::   zs_snic      ! snow-ice salinity
@@ -95,4 +95,5 @@
       REAL(wp), POINTER, DIMENSION(:) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2)
       REAL(wp), POINTER, DIMENSION(:) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2)
+      REAL(wp), POINTER, DIMENSION(:) ::   zevap_rema  ! remaining mass flux from sublimation        (kg.m-2)
 
       REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel   ! snow melt 
@@ -105,6 +106,4 @@
 
       REAL(wp), POINTER, DIMENSION(:) ::   zqh_i       ! total ice heat content  (J.m-2)
-      REAL(wp), POINTER, DIMENSION(:) ::   zqh_s       ! total snow heat content (J.m-2)
-      REAL(wp), POINTER, DIMENSION(:) ::   zq_s        ! total snow enthalpy     (J.m-3)
       REAL(wp), POINTER, DIMENSION(:) ::   zsnw        ! distribution of snow after wind blowing
 
@@ -122,16 +121,15 @@
       END SELECT
 
-      CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw )
-      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
+      CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema )
+      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i )
       CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i )
       CALL wrk_alloc( jpij, nlay_i, icount )
        
-      dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp
+      dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp ; dh_i_sub(:) = 0._wp
       dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp   
 
       zqprec   (:) = 0._wp ; zq_su    (:) = 0._wp ; zq_bo    (:) = 0._wp ; zf_tt(:) = 0._wp
-      zq_rema  (:) = 0._wp ; zsnw     (:) = 0._wp
+      zq_rema  (:) = 0._wp ; zsnw     (:) = 0._wp ; zevap_rema(:) = 0._wp ;
       zdh_s_mel(:) = 0._wp ; zdh_s_pre(:) = 0._wp ; zdh_s_sub(:) = 0._wp ; zqh_i(:) = 0._wp
-      zqh_s    (:) = 0._wp ; zq_s     (:) = 0._wp     
 
       zdeltah(:,:) = 0._wp ; zh_i(:,:) = 0._wp       
@@ -159,8 +157,8 @@
       !
       DO ji = kideb, kiut
-         zfdum      = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
+         zdum       = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
          zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 
 
-         zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
+         zq_su (ji) = MAX( 0._wp, zdum      * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
          zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
       END DO
@@ -187,10 +185,4 @@
       !  2) Computing layer thicknesses and enthalpies.            !
       !------------------------------------------------------------!
-      !
-      DO jk = 1, nlay_s
-         DO ji = kideb, kiut
-            zqh_s(ji) =  zqh_s(ji) + q_s_1d(ji,jk) * ht_s_1d(ji) * r1_nlay_s
-         END DO
-      END DO
       !
       DO jk = 1, nlay_i
@@ -275,17 +267,15 @@
       END DO
 
-      !----------------------
-      ! 3.2 Snow sublimation 
-      !----------------------
+      !------------------------------
+      ! 3.2 Sublimation (part1: snow) 
+      !------------------------------
       ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
       ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean)
-      ! clem comment: ice should also sublimate
       zdeltah(:,:) = 0._wp
-      ! coupled mode: sublimation is set to 0 (evap_ice = 0) until further notice
-      ! forced  mode: snow thickness change due to sublimation
-      DO ji = kideb, kiut
-         zdh_s_sub(ji)  =  MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice )
-         ! Heat flux by sublimation [W.m-2], < 0
-         !      sublimate first snow that had fallen, then pre-existing snow
+      DO ji = kideb, kiut
+         zdh_s_sub(ji)  = MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice )
+         ! remaining evap in kg.m-2 (used for ice melting later on)
+         zevap_rema(ji)  = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhosn
+         ! Heat flux by sublimation [W.m-2], < 0 (sublimate first snow that had fallen, then pre-existing snow)
          zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
          hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * q_s_1d(ji,1)  &
@@ -309,12 +299,10 @@
       !-------------------------------------------
       ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
-      zq_s(:) = 0._wp 
       DO jk = 1, nlay_s
          DO ji = kideb,kiut
-            rswitch       = MAX(  0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 )  )
-            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *                          &
-              &            ( (   zdh_s_pre(ji)             ) * zqprec(ji) +  &
-              &              (   ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
-            zq_s(ji)     =  zq_s(ji) + q_s_1d(ji,jk)
+            rswitch       = MAX( 0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 ) )
+            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *           &
+              &            ( ( zdh_s_pre(ji)               ) * zqprec(ji) +  &
+              &              ( ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
          END DO
       END DO
@@ -370,5 +358,5 @@
                zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0]
                
-               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
+               ! Contribution to salt flux >0 (clem: using sm_i_1d and not s_i_1d(jk) is ok)
                sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
                
@@ -383,4 +371,20 @@
                
             END IF
+            ! ----------------------
+            ! Sublimation part2: ice
+            ! ----------------------
+            zdum      = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoic )
+            zdeltah(ji,jk) = zdeltah(ji,jk) + zdum
+            dh_i_sub(ji)  = dh_i_sub(ji) + zdum
+            ! Salt flux > 0 (clem2016: flux is sent to the ocean for simplicity but salt should remain in the ice except if all ice is melted.
+            !                          It must be corrected at some point)
+            sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * sm_i_1d(ji) * r1_rdtice
+            ! Heat flux [W.m-2], < 0
+            hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * q_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice
+            ! Mass flux > 0
+            wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice
+            ! update remaining mass flux
+            zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoic
+            
             ! record which layers have disappeared (for bottom melting) 
             !    => icount=0 : no layer has vanished
@@ -389,5 +393,5 @@
             icount(ji,jk) = NINT( rswitch )
             zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
-
+                        
             ! update heat content (J.m-2) and layer thickness
             qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
@@ -397,5 +401,11 @@
       ! update ice thickness
       DO ji = kideb, kiut
-         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) )
+         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) + dh_i_sub(ji) )
+      END DO
+
+      ! remaining "potential" evap is sent to ocean
+      DO ji = kideb, kiut
+         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
+         wfx_err_sub(ii,ij) = wfx_err_sub(ii,ij) - zevap_rema(ji) * a_i_1d(ji) * r1_rdtice  ! <=0 (net evap for the ocean in kg.m-2.s-1)
       END DO
 
@@ -686,6 +696,6 @@
       WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp
       
-      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw )
-      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
+      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema )
+      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i )
       CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i )
       CALL wrk_dealloc( jpij, nlay_i, icount )
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 6398)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 6399)
@@ -45,4 +45,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qns_ice_1d  
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_bo_1d     
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rn_amax_1d
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_sum_1d
@@ -51,5 +52,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_dif_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_opw_1d
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rn_amax_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_snw_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_1d
@@ -84,4 +84,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_res_1d  
 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_sub_1d
+
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sprecip_1d    !: <==> the 2D  sprecip
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   frld_1d       !: <==> the 2D  frld
@@ -92,4 +94,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   evap_ice_1d   !: <==> the 2D  evap_ice
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qprec_ice_1d  !: <==> the 2D  qprec_ice
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qevap_ice_1d  !: <==> the 3D  qevap_ice
    !                                                     ! to reintegrate longwave flux inside the ice thermodynamics
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   i0            !: fraction of radiation transmitted to the ice
@@ -108,4 +111,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_s_tot      !: Snow accretion/ablation        [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_surf     !: Ice surface accretion/ablation [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_sub      !: Ice surface sublimation [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_bott     !: Ice bottom accretion/ablation  [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_snowice    !: Snow ice formation             [m of ice]
@@ -155,7 +159,7 @@
          &      wfx_sum_1d(jpij)  , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) , wfx_res_1d(jpij) ,  &
          &      dqns_ice_1d(jpij) , evap_ice_1d (jpij),                                         &
-         &      qprec_ice_1d(jpij), i0         (jpij) ,                                         &  
+         &      qprec_ice_1d(jpij), qevap_ice_1d(jpij), i0         (jpij) ,                     &  
          &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij),  &
-         &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) ,                     &
+         &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij),  &
          &      dsm_i_fl_1d(jpij) , dsm_i_gd_1d(jpij) , dsm_i_se_1d(jpij) ,                     &     
          &      dsm_i_si_1d(jpij) , hicol_1d    (jpij)                     , STAT=ierr(2) )
@@ -163,6 +167,6 @@
       ALLOCATE( t_su_1d   (jpij) , a_i_1d   (jpij) , ht_i_1d  (jpij) ,    &   
          &      ht_s_1d   (jpij) , fc_su    (jpij) , fc_bo_i  (jpij) ,    &    
-         &      dh_s_tot  (jpij) , dh_i_surf(jpij) , dh_i_bott(jpij) ,    &    
-         &      dh_snowice(jpij) , sm_i_1d  (jpij) , s_i_new  (jpij) ,    &
+         &      dh_s_tot  (jpij) , dh_i_surf(jpij) , dh_i_sub (jpij) ,    &    
+         &      dh_i_bott (jpij) ,dh_snowice(jpij) , sm_i_1d  (jpij) , s_i_new  (jpij) ,    &
          &      t_s_1d(jpij,nlay_s) , t_i_1d(jpij,nlay_i) , s_i_1d(jpij,nlay_i) ,  &            
          &      q_i_1d(jpij,nlay_i+1) , q_s_1d(jpij,nlay_s) ,                        &
