snow_melting Subroutine

  subroutine snow_melting(j, jea, ida)
    integer, intent(in) :: ida
    integer, intent(in) :: j, jea
    real(dp) :: ph, smli, sml0, ulmax, tsn, sfreez, svz, svl, snoa0
    real(dp) :: precipe0, smr

    sfreez = 0.
    tsn = 0.
    ph = 0.
    svz = 0.
    svl = 0.
    smli = 0.
    sml0 = 0.
    hsn = 0.
    vsn = 0.
    !    gla_melt_total = 0.
    !#### snow pack depth, hsn, hsn0 (last step), cm
    ! snow, mm
    ! hsn0, snow depth, cm
    ! hsn, cm
    hsn = hsn0 + snowVal * .1 / rnew ! actual snow depth, cm
    snoa0 = hsn0 * rsn(j, jea) * 10. ! previous day snow water content, mm
    precipe0 = precipe ! current day precipitation

    !#### different snowmeltrate for forest and opensite (cm/degree)
    smr = smrate(j) ! cm / °C

    !#### snow temperature
    if (tmit < 0.) tsn = tmit
    !#### calculate ulmax (holding capacity of snowpack)
    ulmax = ulmax0 ! 0 or 1
    if (rsn(j, jea) /= 0. .AND. ulmax0 == 1.) ulmax = - 0.11 * rsn(j, jea) + 0.11
    if (rsn(j, jea) == 0. .AND. ulmax0 == 1.) ulmax = 0.
    ulmax = max(0.0_dp, ulmax)

    !#### calculate intensity of snow melt (cm/day)
    if (tmit > tmelt(j) ) then
      smli = smr * (tmit - tmelt(j)) ! cm (cm / °C * °C )
    end if

    !#### calculate intensity of freezing of liquid water in snowpack (cm/day)
    if (tmit < 0.) sfreez = smr * sqrt(- tmit) ! cm / day
    if (sfreez >= (sul(j, jea) * hsn0) ) sfreez = sul(j, jea) * hsn0 ! is actually g / cm2, similar to cm / day?

    !#### calculate the compression sml0 (cm/day)
    !#### the snowmelt water percolates into pack, hits the ice, and freezes until
    !#### the temperature of the pack is 0.
    if (suz(j, jea) > 0.) then
      if (hsn0 > 0.) then
          ph = .15 * rsn(j, jea) * hsn0 ** 2. / exp(- 0.08 * tsn + 21. * rsn(j, jea))
      end if
      sml0 = smli / (suz(j, jea) * 0.917) + ph
      if (sml0 > hsn) sml0 = hsn
      hsn = hsn - sml0
    else
      smli = 0.
      hsn = snowVal * 0.1 / rnew
    end if

    if (hsn .ne. 0.) then ! SL: equals or lower equal?
      !#### recalculate suz (volumetric content of ice in snow g/cm3)
      svz = 0.917 * suz(j, jea) * hsn0 + snowVal * 0.1 + sfreez
      if (smli .gt. svz) smli = svz
      suz(j, jea) = (svz - smli) / (hsn * 0.917)

      !#### recalculate sul (volumetric content of liquid water in snow g/cm3, svl cm)
      svl = sul(j, jea) * hsn0 + precipe * 0.1 + smli
      if (sfreez > svl) sfreez = svl
      sul(j, jea) = (svl - sfreez) / hsn

      !#### water outflow from snowpack
      if (sul(j, jea) > ulmax ) then
        vsn = hsn * 10. * (sul(j, jea) - ulmax)
        if (vsn > precipe0) smle = vsn - precipe0
        sul(j, jea) = ulmax
      end if
    else
      sul(j, jea) = 0.
      suz(j, jea) = 0.
      vsn = rsn(j, jea) * hsn0 * 10. + precipe0
      smle = rsn(j, jea) * hsn0 * 10.
    end if

    !#### recalculate rsn
    rsn(j, jea) = sul(j, jea) + 0.917 * suz(j, jea)
    if (hsn == 0.) rsn(j, jea) = 0.

    !#### recalculate the total precipitation from snowpack
    precipe = vsn
    snoa(j, jea) = hsn * 10. * rsn(j, jea)
    balanc = snowVal + precipe0 + (snoa0 - snoa(j, jea)) - vsn ! SL: balanc is never used hereafter

    !#### conversion from snow to ice
    !   if(ida.eq.273) then
    if (ida == gla_day_out) then
      ! SL: Maybe one should not convert all the snow into ice at day gla_day_out
      !     but every day a fraction only (snow2ice) !
      !if ( snoa(j, jea) > 2000. ) then
      gla(j, jea) = gla(j, jea) + snoa(j, jea) ! * snow2ice
      snoa(j, jea) = 0. !snoa(j, jea) * (1. - snow2ice)
      hsn = 0. !snoa(j, jea) * .1 / rsn(j, jea)
      rsn(j, jea) = 0.
    end if
  !end if

  end subroutine snow_melting