Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90	(revision 9273)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90	(revision 9274)
@@ -31,5 +31,5 @@
 
    PUBLIC   ice_thd_dh        ! called by ice_thd
-   PUBLIC   ice_thd_snwblow   ! called in sbcblk/sbcclio/sbccpl and here
+   PUBLIC   ice_thd_snwblow   ! called in sbcblk/sbccpl and here
 
    INTERFACE ice_thd_snwblow
@@ -48,18 +48,17 @@
       !!                ***  ROUTINE ice_thd_dh  ***
       !!
-      !! ** Purpose :   compute ice and snow thickness changes due to growing/melting
+      !! ** Purpose :   compute ice and snow thickness changes due to growth/melting
       !!
       !! ** Method  :   Ice/Snow surface melting arises from imbalance in surface fluxes
-      !!              Bottom accretion/ablation arises from flux budget
-      !!              Snow thickness can increase by precipitation and decrease by sublimation
-      !!              If snow load excesses Archmiede limit, snow-ice is formed by
-      !!              the flooding of sea-water in the snow
+      !!                Bottom accretion/ablation arises from flux budget
+      !!                Snow thickness can increase by precipitation and decrease by sublimation
+      !!                If snow load excesses Archmiede limit, snow-ice is formed by
+      !!                the flooding of sea-water in the snow
       !!
-      !!                 1) Compute available flux of heat for surface ablation
-      !!                 2) Compute snow and sea ice enthalpies
-      !!                 3) Surface ablation and sublimation
-      !!                 4) Bottom accretion/ablation
-      !!                 5) Case of Total ablation
-      !!                 6) Snow ice formation
+      !!                - Compute available flux of heat for surface ablation
+      !!                - Compute snow and sea ice enthalpies
+      !!                - Surface ablation and sublimation
+      !!                - Bottom accretion/ablation
+      !!                - Snow ice formation
       !!
       !! References : Bitz and Lipscomb, 1999, J. Geophys. Res.
@@ -126,30 +125,30 @@
       END DO
       !
-      !------------------------------------------------------------------------------!
-      !  1) Calculate available heat for surface and bottom ablation                 !
-      !------------------------------------------------------------------------------!
-      !
-      SELECT CASE ( nice_jules )
-      !
-      CASE ( np_jules_ACTIVE )
+      ! ---------------------------------------------- !
+      ! Available heat for surface and bottom ablation !
+      ! ---------------------------------------------- !
+      !
+      SELECT CASE( nice_jules )
+      !
+      CASE( np_jules_ACTIVE )
          !
          DO ji = 1, npti
-            zq_su(ji)        =   MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
+            zq_su(ji)      = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
          END DO
          !
-      CASE ( np_jules_EMULE )
+      CASE( np_jules_EMULE )
          !
          DO ji = 1, npti
-            zdum             =   ( qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji) ) 
-            qml_ice_1d(ji)   =   zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
-            zq_su(ji)        =   MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
+            zdum           = qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji)
+            qml_ice_1d(ji) = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
+            zq_su(ji)      = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
          END DO
          !
-      CASE ( np_jules_OFF ) 
+      CASE( np_jules_OFF ) 
          !
          DO ji = 1, npti
-            zdum             =   ( qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji) ) 
-            zdum             =   zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
-            zq_su(ji)        =   MAX( 0._wp, zdum             * rdt_ice )
+            zdum           = qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji) 
+            zdum           = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
+            zq_su(ji)      = MAX( 0._wp, zdum * rdt_ice )
          END DO
          !
@@ -157,11 +156,10 @@
       !
       DO ji = 1, npti
-         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 
-         zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
-      END DO
-
-      !--------------------------------------------!
-      !  2) Computing layer thicknesses            !
-      !--------------------------------------------!
+         zf_tt(ji)         = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 
+         zq_bo(ji)         = MAX( 0._wp, zf_tt(ji) * rdt_ice )
+      END DO
+
+      ! Ice and snow layer thicknesses
+      !-------------------------------
       DO jk = 1, nlay_i
          DO ji = 1, npti
@@ -176,79 +174,62 @@
       END DO
 
-      !------------------------------------------------------------------------------!
-      ! If snow temperature is above freezing point, then snow melts 
-      ! (should not happen but sometimes it does)
-      !------------------------------------------------------------------------------!
+      !                       ! ============ !
+      !                       !     Snow     !
+      !                       ! ============ !
+      !
+      ! Internal melting
+      ! ----------------
+      ! IF snow temperature is above freezing point, THEN snow melts (should not happen but sometimes it does)
       DO jk = 1, nlay_s
          DO ji = 1, npti
-            IF( t_s_1d(ji,jk) > rt0 ) THEN !!! Internal melting
-               ! Contribution to heat flux to the ocean [W.m-2], < 0  
-               hfx_res_1d(ji) = hfx_res_1d(ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice
-               ! Contribution to mass flux
-               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice
-               dh_s_mlt(ji)       = dh_s_mlt(ji) - zh_s(ji,jk)
+            IF( t_s_1d(ji,jk) > rt0 ) THEN
+               hfx_res_1d    (ji) = hfx_res_1d    (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice   ! heat flux to the ocean [W.m-2], < 0
+               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn         * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice   ! mass flux
                ! updates
-               h_s_1d(ji)    = h_s_1d(ji) - zh_s(ji,jk)
-               zh_s  (ji,jk) = 0._wp
-               e_s_1d(ji,jk) = 0._wp
-               t_s_1d(ji,jk) = rt0
+               dh_s_mlt(ji)    = dh_s_mlt(ji) - zh_s(ji,jk)
+               h_s_1d  (ji)    = h_s_1d(ji) - zh_s(ji,jk)
+               zh_s    (ji,jk) = 0._wp
+               e_s_1d  (ji,jk) = 0._wp
+               t_s_1d  (ji,jk) = rt0
             END IF
          END DO
-      END DO
-         
-
-      !------------------------------------------------------------------------------|
-      !  3) Surface ablation and sublimation                                         |
-      !------------------------------------------------------------------------------|
-      !
-      !-------------------------
-      ! 3.1 Snow precips / melt
-      !-------------------------
-      ! Snow accumulation in one thermodynamic time step
-      ! snowfall is partitionned between leads and ice
-      ! if snow fall was uniform, a fraction (1-at_i) would fall into leads
-      ! but because of the winds, more snow falls on leads than on sea ice
-      ! and a greater fraction (1-at_i)^beta of the total mass of snow 
-      ! (beta < 1) falls in leads.
-      ! In reality, beta depends on wind speed, 
-      ! and should decrease with increasing wind speed but here, it is 
-      ! considered as a constant. an average value is 0.66
-      ! Martin Vancoppenolle, December 2006
-
-      CALL ice_thd_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing
+      END DO         
+
+      ! Snow precipitation
+      !-------------------
+      CALL ice_thd_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) )   ! snow distribution over ice after wind blowing
 
       zdeltah(1:npti,:) = 0._wp
       DO ji = 1, npti
-         !-----------
-         ! Snow fall
-         !-----------
-         ! thickness change
-         zdh_s_pre(ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji)
-         ! enthalpy of the precip (>0, J.m-3)
-         zqprec   (ji) = - qprec_ice_1d(ji)   
-         IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp
-         ! heat flux from snow precip (>0, W.m-2)
-         hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
-         ! mass flux, <0
-         wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice
-
-         !---------------------
-         ! Melt of falling snow
-         !---------------------
-         ! thickness change
-         rswitch        = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )
-         zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 )
-         zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting 
-         ! heat used to melt snow (W.m-2, >0)
-         hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
-         ! snow melting only = water into the ocean (then without snow precip), >0
-         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
-         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1)
-         ! updates available heat + precipitations after melting
-         zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )      
-         zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1)
-
-         ! update thickness
-         h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) )
+         IF( sprecip_1d(ji) > 0._wp ) THEN
+            !
+            ! --- precipitation ---
+            zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji)   ! thickness change
+            zqprec    (ji) = - qprec_ice_1d(ji)                                             ! enthalpy of the precip (>0, J.m-3)
+            !
+            hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji)    * r1_rdtice   ! heat flux from snow precip (>0, W.m-2)
+            wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn         * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice   ! mass flux, <0
+            
+            ! --- melt of falling snow ---
+            rswitch              = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )
+            zdeltah       (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 )   ! thickness change
+            zdeltah       (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) )                ! bound melting 
+            hfx_snw_1d    (ji)   = hfx_snw_1d    (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji)    * r1_rdtice   ! heat used to melt snow (W.m-2, >0)
+            wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhosn         * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip), >0
+            
+            ! updates available heat + precipitations after melting
+            dh_s_mlt (ji) = dh_s_mlt(ji) + zdeltah(ji,1)
+            zq_su    (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )      
+            zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
+            
+            ! update thickness
+            h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) )
+            !
+         ELSE
+            !
+            zdh_s_pre(ji) = 0._wp
+            zqprec   (ji) = 0._wp
+            !
+         ENDIF
       END DO
 
@@ -260,4 +241,6 @@
       END DO
 
+      ! Snow melting
+      ! ------------
       ! If heat still available (zq_su > 0), then melt more snow
       zdeltah(1:npti,:) = 0._wp
@@ -265,44 +248,53 @@
       DO jk = 1, nlay_s
          DO ji = 1, npti
-            ! thickness change
-            rswitch          = 1._wp - MAX( 0._wp, SIGN( 1._wp, - zh_s(ji,jk) ) ) 
-            rswitch          = rswitch * ( MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) ) ) 
-            zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( e_s_1d(ji,jk), epsi20 )
-            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) ) ! bound melting
-            zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)    
-            ! heat used to melt snow(W.m-2, >0)
-            hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d(ji,jk) * r1_rdtice 
-            ! snow melting only = water into the ocean (then without snow precip)
-            wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
-            dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,jk)
-            ! updates available heat + thickness
-            zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
-            h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) )
-            zh_s(ji,jk)= MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) )
-         END DO
-      END DO
-
-      !------------------------------
-      ! 3.2 Sublimation (part1: snow) 
-      !------------------------------
+            IF( zh_s(ji,jk) > 0._wp .AND. zq_su(ji) > 0._wp ) THEN
+               !
+               rswitch          = MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) )
+               zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( e_s_1d(ji,jk), epsi20 )   ! thickness change
+               zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) )                  ! bound melting
+               zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)
+               
+               hfx_snw_1d(ji)     = hfx_snw_1d(ji)     - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_rdtice   ! heat used to melt snow(W.m-2, >0)
+               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn          * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip)
+               
+               ! updates available heat + thickness
+               dh_s_mlt(ji)    = dh_s_mlt(ji) + zdeltah(ji,jk)
+               zq_su   (ji)    = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
+               h_s_1d  (ji)    = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) )
+               zh_s    (ji,jk) = MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) )
+               !
+            ENDIF
+         END DO
+      END DO
+
+      ! Snow sublimation 
+      !-----------------
       ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
       !    comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean)
       zdeltah(1:npti,:) = 0._wp
       DO ji = 1, npti
-         zdh_s_sub(ji)  = MAX( - h_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) ) * e_s_1d(ji,1)  &
-            &                              ) * a_i_1d(ji) * r1_rdtice
-         ! Mass flux by sublimation
-         wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice
-
-         ! new snow thickness
-         h_s_1d(ji)    =  MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) )
-         ! update precipitations after sublimation and correct sublimation
-         zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
-         zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1)
+         IF( evap_ice_1d(ji) > 0._wp ) THEN
+            !
+            zdh_s_sub (ji)   = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice )
+            zevap_rema(ji)   = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhosn   ! remaining evap in kg.m-2 (used for ice melting later on)
+            zdeltah   (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
+            
+            hfx_sub_1d    (ji) = hfx_sub_1d(ji) + &   ! Heat flux by sublimation [W.m-2], < 0 (sublimate snow that had fallen, then pre-existing snow)
+               &                 ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) )  &
+               &                 * a_i_1d(ji) * r1_rdtice
+            wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice   ! Mass flux by sublimation
+            
+            ! new snow thickness
+            h_s_1d(ji)    =  MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) )
+            ! update precipitations after sublimation and correct sublimation
+            zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
+            zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1)
+            !
+         ELSE
+            !
+            zdh_s_sub (ji) = 0._wp
+            zevap_rema(ji) = 0._wp
+            !
+         ENDIF
       END DO
       
@@ -312,29 +304,31 @@
       END DO
 
-      !-------------------------------------------
-      ! 3.3 Update temperature, energy
-      !-------------------------------------------
+      ! Update temperature, energy
+      !---------------------------
       ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
       DO jk = 1, nlay_s
          DO ji = 1,npti
             rswitch       = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) )
-            e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) *           &
-              &            ( ( zdh_s_pre(ji)               ) * zqprec(ji) +  &
-              &              ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
-         END DO
-      END DO
-
-      !--------------------------
-      ! 3.4 Surface ice ablation 
-      !--------------------------
+            e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) *            &
+              &             ( ( zdh_s_pre(ji)              ) * zqprec(ji) +  &
+              &               ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
+         END DO
+      END DO
+      
+      !                       ! ============ !
+      !                       !     Ice      !
+      !                       ! ============ !
+
+      ! Surface ice melting 
+      !--------------------
       zdeltah(1:npti,:) = 0._wp ! important
       DO jk = 1, nlay_i
          DO ji = 1, npti
-            ztmelts           = - tmut * sz_i_1d(ji,jk)          ! Melting point of layer k [C]
-            
-            IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !!! Internal melting
+            ztmelts = - tmut * sz_i_1d(ji,jk)   ! Melting point of layer k [C]
+            
+            IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting
 
                zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]       
-               zdE            = 0._wp                                 ! Specific enthalpy difference   (J/kg, <0)
+               zdE            = 0._wp                                 ! Specific enthalpy difference (J/kg, <0)
                                                                       ! set up at 0 since no energy is needed to melt water...(it is already melted)
                zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing     
@@ -346,14 +340,11 @@
                zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0]
 
-               ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) 
-               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
-               
-               ! Contribution to salt flux (using s_i_1d and not sz_i_1d(jk) is ok)
-               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice
-               
-               ! Contribution to mass flux
-               wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
-
-            ELSE                                !!! Surface melting
+               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice                           ! Heat flux to the ocean [W.m-2], <0
+               !                                                                                                  ice enthalpy zEi is "sent" to the ocean
+               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice   ! Salt flux
+               !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok
+               wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                ! Mass flux
+
+            ELSE                                        !-- Surface melting
                
                zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]
@@ -375,30 +366,26 @@
                zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0]
                
-               ! Contribution to salt flux >0 (using s_i_1d and not sz_i_1d(jk) is ok)
-               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice
-               
-               ! Contribution to heat flux [W.m-2], < 0
-               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
-               
-               ! Total heat flux used in this process [W.m-2], > 0  
-               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
-               
-               ! Contribution to mass flux
-               wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
+               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice   ! Salt flux >0
+               !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok)
+               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux [W.m-2], < 0
+               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat flux used in this process [W.m-2], > 0  
+               ! 
+               wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                ! Mass flux
                
             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 (clem: 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 * s_i_1d(ji) * r1_rdtice
-            ! Heat flux [W.m-2], < 0
-            hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice
-            ! Mass flux > 0
-            wfx_ice_sub_1d(ji) =  wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice
+            
+            ! Ice sublimation
+            ! ---------------
+            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
+            
+            sfx_sub_1d(ji)     = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_rdtice ! Salt flux >0
+            !                                                                                          clem: flux is sent to the ocean for simplicity
+            !                                                                                                but salt should remain in the ice except
+            !                                                                                                if all ice is melted. => must be corrected
+            hfx_sub_1d(ji)     = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice      ! Heat flux [W.m-2], < 0
+
+            wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice          ! Mass flux > 0
 
             ! update remaining mass flux
@@ -417,4 +404,5 @@
          END DO
       END DO
+      
       ! update ice thickness
       DO ji = 1, npti
@@ -428,10 +416,5 @@
 
 
-      !------------------------------------------------------------------------------!
-      ! 4) Basal growth / melt                                                       !
-      !------------------------------------------------------------------------------!
-      !
-      !------------------
-      ! 4.1 Basal growth 
+      ! Ice Basal growth 
       !------------------
       ! Basal growth is driven by heat imbalance at the ice-ocean interface,
@@ -449,8 +432,7 @@
       IF( nn_icesal == 2 )   num_iter_max = 5  ! salinity varying in time
       
-      ! Iterative procedure
       DO ji = 1, npti
          IF(  zf_tt(ji) < 0._wp  ) THEN
-            DO iter = 1, num_iter_max
+            DO iter = 1, num_iter_max   ! iterations
 
                ! New bottom ice salinity (Cox & Weeks, JGR88 )
@@ -458,55 +440,47 @@
                !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7 
                !--- zswi2  if dh/dt > 3.6e-7
-               zgrr               = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) )
-               zswi2              = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
-               zswi12             = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
-               zswi1              = 1. - zswi2 * zswi12
-               zfracs             = MIN ( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
-                  &               + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 )
-
-               s_i_new(ji)        = zswitch_sal * zfracs * sss_1d(ji)  &  ! New ice salinity
-                                  + ( 1. - zswitch_sal ) * s_i_1d(ji) 
-               ! New ice growth
-               ztmelts            = - tmut * s_i_new(ji)                 ! New ice melting point (C)
-
-               zt_i_new           = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
-               
-               zEi                = cpic * ( zt_i_new - (ztmelts+rt0) ) &     ! Specific enthalpy of forming ice (J/kg, <0)      
-                  &               - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) )   &
-                  &               + rcp  * ztmelts          
-
-               zEw                = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
-
-               zdE                = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
-
-               dh_i_bott(ji)      = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
+               zgrr     = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) )
+               zswi2    = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
+               zswi12   = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
+               zswi1    = 1. - zswi2 * zswi12
+               zfracs   = MIN( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
+                  &          + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 )
+
+               s_i_new(ji)   = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji)  ! New ice salinity
+
+               ztmelts       = - tmut * s_i_new(ji)                                                   ! New ice melting point (C)
+
+               zt_i_new      = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
+               
+               zEi           = cpic * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0)
+                  &            - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp  * ztmelts
+
+               zEw           = rcp  * ( t_bo_1d(ji) - rt0 )                                           ! Specific enthalpy of seawater (J/kg, < 0)
+
+               zdE           = zEi - zEw                                                              ! Specific enthalpy difference (J/kg, <0)
+
+               dh_i_bott(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
                
             END DO
             ! Contribution to Energy and Salt Fluxes                                    
-            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0)
-
-            ztmelts        = - tmut * s_i_new(ji)                ! New ice melting point (C)
+            zfmdt          = - rhoic * dh_i_bott(ji)                                                  ! Mass flux x time step (kg/m2, < 0)
+
+            ztmelts        = - tmut * s_i_new(ji)                                                     ! New ice melting point (C)
             
             zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
             
-            zEi            = cpic * ( zt_i_new - (ztmelts+rt0) ) &     ! Specific enthalpy of forming ice (J/kg, <0)      
-               &               - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) )   &
-               &               + rcp  * ztmelts          
-            
-            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
-            
-            zdE            = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
-            
-            ! Contribution to heat flux to the ocean [W.m-2], >0  
-            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
-
-            ! Total heat flux used in this process [W.m-2], <0  
-            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
-            
-            ! Contribution to salt flux, <0
-            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice
-
-            ! Contribution to mass flux, <0
-            wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice
+            zEi            = cpic * ( zt_i_new - (ztmelts+rt0) ) &                                    ! Specific enthalpy of forming ice (J/kg, <0)
+               &           - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp  * ztmelts          
+            
+            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )                                             ! Specific enthalpy of seawater (J/kg, < 0)
+            
+            zdE            = zEi - zEw                                                                ! Specific enthalpy difference (J/kg, <0)
+            
+            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux to the ocean [W.m-2], >0
+            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat flux used in this process [W.m-2], <0
+            
+            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice   ! Salt flux, <0
+
+            wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice                 ! Mass flux, <0
 
             ! update heat content (J.m-2) and layer thickness
@@ -518,7 +492,6 @@
       END DO
 
-      !----------------
-      ! 4.2 Basal melt
-      !----------------
+      ! Ice Basal melt
+      !---------------
       zdeltah(1:npti,:) = 0._wp ! important
       DO jk = nlay_i, 1, -1
@@ -528,9 +501,9 @@
                ztmelts = - tmut * sz_i_1d(ji,jk)  ! Melting point of layer jk (C)
 
-               IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !!! Internal melting
+               IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting
 
                   zEi               = - e_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0)
                   zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0)
-                                                                    ! set up at 0 since no energy is needed to melt water...(it is already melted)
+                                                                    !    set up at 0 since no energy is needed to melt water...(it is already melted)
                   zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing     
                                                                     ! this should normally not happen, but sometimes, heat diffusion leads to this
@@ -540,12 +513,9 @@
                   zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0
 
-                  ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) 
-                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
-
-                  ! Contribution to salt flux (using s_i_1d and not sz_i_1d(jk) is ok)
-                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice
-                                    
-                  ! Contribution to mass flux
-                  wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
+                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice                           ! Heat flux to the ocean [W.m-2], <0
+                  !                                                                                                  ice enthalpy zEi is "sent" to the ocean
+                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice   ! Salt flux
+                  !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok
+                  wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                ! Mass flux
 
                   ! update heat content (J.m-2) and layer thickness
@@ -553,35 +523,31 @@
                   h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
 
-               ELSE                               !!! Basal melting
-
-                  zEi             = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)
-                  zEw             = rcp * ztmelts              ! Specific enthalpy of meltwater (J/kg, <0)
-                  zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0)
-
-                  zfmdt           = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0)
-
-                  zdeltah(ji,jk)  = - zfmdt * r1_rhoic         ! Gross thickness change
-
-                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
+               ELSE                                        !-- Basal melting
+
+                  zEi             = - e_i_1d(ji,jk) * r1_rhoic                                ! Specific enthalpy of melting ice (J/kg, <0)
+                  zEw             = rcp * ztmelts                                             ! Specific enthalpy of meltwater (J/kg, <0)
+                  zdE             = zEi - zEw                                                 ! Specific enthalpy difference   (J/kg, <0)
+
+                  zfmdt           = - zq_bo(ji) / zdE                                         ! Mass flux x time step (kg/m2, >0)
+
+                  zdeltah(ji,jk)  = - zfmdt * r1_rhoic                                        ! Gross thickness change
+
+                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) )       ! bound thickness change
                   
-                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors
-
-                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt
-
-                  zfmdt           = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0
-
-                  zQm             = zfmdt * zEw         ! Heat exchanged with ocean
-
-                  ! Contribution to heat flux to the ocean [W.m-2], <0  
-                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
-
-                  ! Contribution to salt flux (using s_i_1d and not sz_i_1d(jk) is ok)
-                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice
+                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE )   ! update available heat. MAX is necessary for roundup errors
+
+                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)                            ! Update basal melt
+
+                  zfmdt           = - zdeltah(ji,jk) * rhoic                                  ! Mass flux x time step > 0
+
+                  zQm             = zfmdt * zEw                                               ! Heat exchanged with ocean
+
+                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux to the ocean [W.m-2], <0  
+                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat used in this process [W.m-2], >0  
+
+                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice  ! Salt flux
+                  !                                                                                                   using s_i_1d and not sz_i_1d(jk) is ok
                   
-                  ! Total heat flux used in this process [W.m-2], >0  
-                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
-                  
-                  ! Contribution to mass flux
-                  wfx_bom_1d(ji)  =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
+                  wfx_bom_1d(ji)  = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                ! Mass flux
 
                   ! update heat content (J.m-2) and layer thickness
@@ -594,35 +560,29 @@
       END DO
 
-      !-------------------------------------------
       ! Update temperature, energy
-      !-------------------------------------------
-      DO ji = 1, npti
-         h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_bott(ji) )
+      ! --------------------------
+      DO ji = 1, npti
+         h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bott(ji) )
       END DO  
 
-      !-------------------------------------------
-      ! 5. What to do with remaining energy
-      !-------------------------------------------
-      ! If heat still available for melting and snow remains, then melt more snow
+      ! If heat still available then melt more snow
       !-------------------------------------------
       zdeltah(1:npti,:) = 0._wp ! important
       DO ji = 1, npti
-         zq_rema(ji)     = zq_su(ji) + zq_bo(ji) 
-         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) )   ! =1 if snow
-         rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) )
-         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 )
-         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting
-         dh_s_tot (ji)   = dh_s_tot(ji)  + zdeltah(ji,1)
-         h_s_1d   (ji)   = h_s_1d(ji)   + zdeltah(ji,1)
+         zq_rema (ji)   = zq_su(ji) + zq_bo(ji) 
+         rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) )   ! =1 if snow
+         rswitch        = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) )
+         zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 )
+         zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting
+         dh_s_tot(ji)   = dh_s_tot(ji) + zdeltah(ji,1)
+         h_s_1d  (ji)   = h_s_1d  (ji) + zdeltah(ji,1)
         
-         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                ! update available heat (J.m-2)
-         ! heat used to melt snow
-         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0)
-         ! Contribution to mass flux
-         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
+         zq_rema(ji)        = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                               ! update available heat (J.m-2)
+         hfx_snw_1d(ji)     = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice   ! Heat used to melt snow, W.m-2 (>0)
+         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice      ! Mass flux
          dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1)
          !    
          ! Remaining heat flux (W.m-2) is sent to the ocean heat budget
-         hfx_out_1d(ji)  = hfx_out_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
+         hfx_out_1d(ji) = hfx_out_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
 
          IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
@@ -630,7 +590,6 @@
 
       !
-      !------------------------------------------------------------------------------|
-      !  6) Snow-Ice formation                                                       |
-      !------------------------------------------------------------------------------|
+      ! Snow-Ice formation
+      ! ------------------
       ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level, 
       ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
@@ -648,9 +607,7 @@
          zQm            = zfmdt * zEw 
          
-         ! Contribution to heat flux
-         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
-
-         ! Contribution to salt flux
-         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice 
+         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux
+
+         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice ! Salt flux
 
          ! Case constant salinity in time: virtual salt flux to keep salinity constant
@@ -660,6 +617,5 @@
          ENDIF
 
-         ! Contribution to mass flux
-         ! All snow is thrown in the ocean, and seawater is taken to replace the volume
+         ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume
          wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice
          wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice
@@ -672,10 +628,9 @@
 
       !
-      !-------------------------------------------
       ! Update temperature, energy
-      !-------------------------------------------
-      DO ji = 1, npti
-         rswitch     =  1.0 - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 
-         t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0
+      ! --------------------------
+      DO ji = 1, npti
+         rswitch     = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 
+         t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1._wp - rswitch ) * rt0
       END DO
 
@@ -698,5 +653,16 @@
    !!--------------------------------------------------------------------------
    !! INTERFACE ice_thd_snwblow
+   !!
    !! ** Purpose :   Compute distribution of precip over the ice
+   !!
+   !!                Snow accumulation in one thermodynamic time step
+   !!                snowfall is partitionned between leads and ice.
+   !!                If snow fall was uniform, a fraction (1-at_i) would fall into leads
+   !!                but because of the winds, more snow falls on leads than on sea ice
+   !!                and a greater fraction (1-at_i)^beta of the total mass of snow 
+   !!                (beta < 1) falls in leads.
+   !!                In reality, beta depends on wind speed, 
+   !!                and should decrease with increasing wind speed but here, it is 
+   !!                considered as a constant. an average value is 0.66
    !!--------------------------------------------------------------------------
 !!gm  I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE....
