Index: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8156)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8158)
@@ -252,19 +252,21 @@
    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   [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_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     !: sublimation of snow/ice    [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sub !: snow sublimation           [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice_sub !: ice sublimation            [kg.m-2.s-1]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_dyn !: dynamical component of wfx_snw    [kg.m-2.s-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sum !: surface melt component of wfx_snw [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(:,:) ::   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]
@@ -467,11 +469,11 @@
 
       ii = ii + 1
-      ALLOCATE( sist   (jpi,jpj) , t_bo   (jpi,jpj) ,                                           &
-         &      frld   (jpi,jpj) , pfrld  (jpi,jpj) , phicif (jpi,jpj) ,                        &
-         &      wfx_snw(jpi,jpj) , wfx_snw_dyn(jpi,jpj), wfx_snw_sum(jpi,jpj) ,                 &
-         &      wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) ,                                           &
-         &      wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,     &
-         &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,     &
-         &      afx_tot(jpi,jpj) , afx_thd(jpi,jpj),  afx_dyn(jpi,jpj) ,                        &
+      ALLOCATE( sist   (jpi,jpj) , t_bo   (jpi,jpj) ,                                               &
+         &      frld   (jpi,jpj) , pfrld  (jpi,jpj) , phicif (jpi,jpj) ,                            &
+         &      wfx_snw(jpi,jpj) , wfx_snw_dyn(jpi,jpj), wfx_snw_sum(jpi,jpj) ,                     &
+         &      wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) , wfx_snw_sub(jpi,jpj) , wfx_ice_sub(jpi,jpj) , &
+         &      wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,         &
+         &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,         &
+         &      afx_tot(jpi,jpj) , afx_thd(jpi,jpj),  afx_dyn(jpi,jpj) ,                            &
          &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl), pahu3D(jpi,jpj,jpl+1), pahv3D(jpi,jpj,jpl+1),            &
          &      qlead  (jpi,jpj) , rn_amax_2d(jpi,jpj),                                         &
Index: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8156)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8158)
@@ -152,4 +152,25 @@
             zfric_u      = MAX( SQRT( ust2s(ji,jj) ), zfric_umin ) 
             fhtur(ji,jj) = MAX( 0._wp, rswitch * rau0 * rcp * zch  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2
+
+            !-> MV
+            !-> the following line makes hfxtur often crazily negative
+            !-> if sst is warmer than freezing, evt is fine 
+            !-> but if sst is super-cooled, fhtur = 0 and zqfr > 0
+            !-> so fhtur becomes - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 )
+            !-> if at_i is small, fhtur can be big
+            !-> an we will impose the ice to give energy to the ocean to absorb all super cooling!
+            !-> makes very high values for small ice concentrations, and grows ice vertically very rapidly
+            !-> eg 5th category in hudson bay
+            !-> should check after 10 years, but within a year fhtur is too negative
+            !-> that may help to build up towers
+            ! could we test to put this energy in zqld
+            ! zqfr_pos kept in fhtur
+            ! zqfr_neg added to zqld
+            ! zqfr_pos = MAX( tmask(ji,jj,1) * rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ), 0. )
+            ! zqfr_neg = MIN( tmask(ji,jj,1) * rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ), 0. )
+            ! fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - zqfr_pos * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ) 
+            ! how do we deal with qlead ?
+
+
             fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
             ! upper bound for fhtur: the heat retrieved from the ocean must be smaller than the heat necessary to reach 
@@ -531,4 +552,6 @@
          CALL tab_2d_1d( nbpb, wfx_snw_sum_1d(1:nbpb), wfx_snw_sum  , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub         , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d( nbpb, wfx_snw_sub_1d(1:nbpb), wfx_snw_sub  , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d( nbpb, wfx_ice_sub_1d(1:nbpb), wfx_ice_sub  , jpi, jpj, npb(1:nbpb) )
          
          CALL tab_2d_1d( nbpb, wfx_bog_1d (1:nbpb), wfx_bog         , jpi, jpj, npb(1:nbpb) )
@@ -587,4 +610,6 @@
          CALL tab_1d_2d( nbpb, wfx_snw_sum   , npb, wfx_snw_sum_1d(1:nbpb),jpi, jpj )
          CALL tab_1d_2d( nbpb, wfx_sub       , npb, wfx_sub_1d(1:nbpb)   , jpi, jpj )
+         CALL tab_1d_2d( nbpb, wfx_snw_sub   , npb, wfx_snw_sub_1d(1:nbpb), jpi, jpj )
+         CALL tab_1d_2d( nbpb, wfx_ice_sub   , npb, wfx_ice_sub_1d(1:nbpb), jpi, jpj )
          
          CALL tab_1d_2d( nbpb, wfx_bog       , npb, wfx_bog_1d(1:nbpb)   , jpi, jpj )
Index: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 8156)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 8158)
@@ -282,5 +282,6 @@
             &                              ) * a_i_1d(ji) * r1_rdtice
          ! Mass flux by sublimation
-         wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice
+         wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice
+
          ! new snow thickness
          ht_s_1d(ji)    =  MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) )
@@ -383,5 +384,6 @@
             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
+            wfx_ice_sub_1d(ji) =  wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice
+
             ! update remaining mass flux
             zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoic
@@ -621,7 +623,8 @@
       END DO
 
-      ! add snow melt to mass snow mass flux to the ocean
-      DO ji = kideb, kiut
-         wfx_snw_1d(ji) = wfx_snw_1d(ji) + wfx_snw_sum_1d(ji)
+      ! Water fluxes
+      DO ji = kideb, kiut
+         wfx_snw_1d(ji) = wfx_snw_1d(ji)     + wfx_snw_sum_1d(ji) ! add snow melt to snow mass flux to the ocean
+         wfx_sub_1d(ji) = wfx_snw_sub_1d(ji) + wfx_ice_sub_1d(ji) ! sum ice and snow sublimation contributions
       END DO
 
Index: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 8156)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 8158)
@@ -276,5 +276,9 @@
       IF  ( iom_use( "dmisum"   ) ) CALL iom_put( "dmisum"      , - wfx_sum                  )          ! Sea-ice mass change through surface melting
       IF  ( iom_use( "dmibom"   ) ) CALL iom_put( "dmibom"      , - wfx_bom                  )          ! Sea-ice mass change through bottom melting
+
       IF  ( iom_use( "dmtsub"   ) ) CALL iom_put( "dmtsub"      , - wfx_sub                  )          ! Sea-ice mass change through evaporation and sublimation
+      IF  ( iom_use( "dmssub"   ) ) CALL iom_put( "dmssub"      , - wfx_snw_sub              )          ! Snow mass change through sublimation
+      IF  ( iom_use( "dmisub"   ) ) CALL iom_put( "dmisub"      , - wfx_ice_sub              )          ! Sea-ice mass change through sublimation
+
       IF  ( iom_use( "dmsspr"   ) ) CALL iom_put( "dmsspr"      , - wfx_spr                  )          ! Snow mass change through snow fall
       IF  ( iom_use( "dmsssi"   ) ) CALL iom_put( "dmsssi"      ,   wfx_sni*rhosn/rhoic      )          ! Snow mass change through snow-to-ice conversion
Index: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 8156)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 8158)
@@ -68,4 +68,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_snw_sum_1d 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_sub_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_snw_sub_1d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_ice_sub_1d 
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_bog_1d    
@@ -164,4 +166,5 @@
          &      fhld_1d    (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d (jpij) , wfx_bom_1d(jpij) ,  &
          &      wfx_sum_1d(jpij)  , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) , wfx_res_1d(jpij) ,  &
+         &      wfx_snw_sub_1d(jpij), wfx_ice_sub_1d(jpij) ,                                    &
          &      dqns_ice_1d(jpij) , evap_ice_1d (jpij),                                         &
          &      qprec_ice_1d(jpij), qevap_ice_1d(jpij), i0         (jpij) ,                     &  
Index: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 8156)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 8158)
@@ -597,4 +597,5 @@
 
       wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp
+      wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp
       
       hfx_thd(:,:) = 0._wp   ;   
