!--------1---------2---------3---------4---------5---------6---------7---------8
subroutine flux_melt_below(flux_in) 
!--------1---------2---------3---------4---------5---------6---------7---------8
  use global_constants, only : rkind, dt, time, year, density, Myr &
    & , W_ALTERATION
  use global_parameters, only : temp, radius, dr, surface_arear, halfdr, surface_area
  use heat_parameters, only : thermal_conductivity, latent_heat_init, specific_heat
  use heat_radiogenic, only : iAl, ftime
  use material_parameters, only : mn, iw, ir &
    & , melt, rock_water, dr_rock, radius_rock &
    & , radius_water, three_water, drock, dr_water, mass_ratio_rock &
    & , mass_ratio_ice, mass_ratio_water, material_select &
    & , mass_ratio_anhydrous_rock, mass_ratio_hydrous_rock &
    & , frez_start, radius_ice_water
  use ice_parameters, only : temperature_melting_I
  use energy, only : lat_energy, flux_rock
  use aqueous, only : aqueous_alteration_below, flux_alt
  implicit none
  real(rkind),intent(in) :: flux_in
  real(rkind) :: coef_latent, dlat_energy
  if(mn(iw+1) == 7 .OR. mn(iw+1) == 8) then
    mn(iw+1) = 8
  elseif(mn(iw+1) == 11) then
    mn(iw+1) = 11
  elseif(mn(iw+1) == 9) then
    frez_start=.FALSE.
    call material_select(iw+1)
    coef_latent = latent_heat_init&
      & *surface_arear(radius_water)*density(mn(iw+1))
    lat_energy(iw+1) = lat_energy(iw+1) + flux_in + flux_alt
    dr_water = (flux_in+flux_alt)*dt/(coef_latent)
    radius_water = radius_water + dr_water
    if(radius_water >= radius_ice_water) mn(iw+1) = 11
    return
  elseif(mn(iw+1) == 3 .OR. mn(iw+1)==10) then
    frez_start=.FALSE.
    call material_select(iw+1)
    coef_latent = latent_heat_init&
      & *surface_arear(radius_water)*density(mn(iw+1))
    lat_energy(iw+1) = lat_energy(iw+1) + flux_in + flux_alt
    dr_water = (flux_in+flux_alt)*dt/(coef_latent)
    radius_water = radius_water + dr_water
    return
  else
    write(999,*) 'error occurs at flux_check below!', (time/year+ftime(iAl))/Myr, mn(iw+1), iw+1
    call log_error(iw+1,'flux_check_below mn(iw+1))')
    return
  endif
  if(mn(iw+1) == 9) then
    call log_error(iw+1,'flux_check_below slip mn9')
    write(6,'(a)') radius_water, radius_ice_water
    stop
  endif
  call material_select(iw+1)
  coef_latent = latent_heat_init*mass_ratio_ice(iw+1)&
    & *surface_arear(radius_water)*density(mn(iw+1))
  lat_energy(iw+1) = lat_energy(iw+1) + flux_in + flux_alt
  dr_water = (flux_in+flux_alt)*dt/(coef_latent)
  radius_water = radius_water + dr_water
  if(W_ALTERATION) then
    call aqueous_alteration_below(iw+1,dr_water)
  elseif(.NOT.W_ALTERATION) then
    drock = drock + dr_water*(mass_ratio_anhydrous_rock(iw+1)/density(4)+&
      & mass_ratio_hydrous_rock(iw+1)/density(4))*density(6)*radius_water*radius_water
  endif
  if(radius_water >= radius(iw+1)+halfdr) then
    mass_ratio_ice(iw+1) = 0._rkind
    mass_ratio_water(iw+1) = 1._rkind
    dlat_energy = coef_latent*(radius_water-(radius(iw+1)+halfdr))/dt
    lat_energy(iw+1) = lat_energy(iw+1) - dlat_energy
    mn(iw+1) = 1
    temp(iw+1) = temperature_melting_I
    melt(iw+2) = .TRUE.
    lat_energy(iw+2) = lat_energy(iw+2) + dlat_energy
    if(mn(iw+2) == 7) then
      mn(iw+2) = 8
    elseif(mn(iw+2) == 6) then
      mn(iw+2) = 11
    else
      write(999,*) 'error occurs at flux_check below mn(iw+2)!', log10(time/year), mn(iw+2), iw+2, temp(iw+2)
      call log_error(iw+2,'flux_check_below mn(iw+2)')
      return
    endif
  endif
end subroutine flux_melt_below