bug in snow-ice management

If the snow load exceeds a certain threshold (when the ice-snow interface is depleted below sea surface), then snow-ice formation occurs. However, in SI3 this formation is interrupted and snow is put in the ocean directly when SST becomes positive but it should not. This bug has been introduced long time ago and hence must be corrected in 4.0, 4.2 and main branches. The following lines in icethd_dh.F90 must be erased:

     ! Snow load on ice
     ! -----------------
     ! When snow load exceeds Archimede's limit and sst is positive,
     ! snow-ice formation (next bloc) can lead to negative ice enthalpy.
     ! Therefore we consider here that this excess of snow falls into the ocean
     zdeltah = h_s_1d(ji) + h_i_1d(ji) * (rhoi-rho0) * r1_rhos
     DO jk = 0, nlay_s
        IF( zdeltah > 0._wp .AND. sst_1d(ji) > 0._wp ) THEN
           ! snow layer thickness that falls into the ocean
           zdum = MIN( zdeltah , zh_s(jk) )
           ! mass & energy loss to the ocean
           hfx_res_1d(ji) = hfx_res_1d(ji) - ze_s(jk) * zdum * a_i_1d(ji) * r1_Dt_ice  ! heat flux to the ocean [W.m-2], < 0
           wfx_res_1d(ji) = wfx_res_1d(ji) + rhos     * zdum * a_i_1d(ji) * r1_Dt_ice  ! mass flux
           ! update thickness and energy
           h_s_1d(ji) = MAX( 0._wp, h_s_1d(ji) - zdum )
           zh_s  (jk) = MAX( 0._wp, zh_s  (jk) - zdum )
           ! update snow thickness that still has to fall
           zdeltah   = MAX( 0._wp, zdeltah - zdum )
        ENDIF
     END DO