!--------1---------2---------3---------4---------5---------6---------7---------8
subroutine flux_check 
!--------1---------2---------3---------4---------5---------6---------7---------8
  use global_constants, only : rkind, time, year, volume_ratio_rock, Myr
  use global_parameters, only : temp, dr, surface_arear, halfdr, surface_area, radius
  use heat_radiogenic, only : iAl, ftime
  use material_parameters, only : iw, ir, ii, radius_rock, radius_ice, mn &
    & , radius_water, drock, material_select
  use heat_parameters, only : thermal_conductivity
  use ice_parameters, only : temperature_melting_I
  use energy, only : flux_rock
  implicit none
  real(rkind) :: flux_in, flux_out, drr, K
  flux_in = flux_rock*surface_arear(radius_rock)
  if(ii-iw >= 2 ) then
    call flux_melt_below(flux_in)
  elseif(ii-iw <= 1) then  
    if(radius_water /= radius(iw+1)) then
      drr = (dr+radius(iw+2)-radius_water)/(radius(iw+2)-radius_water)
    elseif(radius_water == radius(iw+1)) then
      drr = 2._rkind
    endif
    call material_select(iw+2)
    K = thermal_conductivity(mn(iw+2))
    if(mn(iw+2)==2) then
      call material_select(iw+3)
      K = (K+thermal_conductivity(mn(iw+3)))*0.5_rkind
    endif
    flux_out = -K*surface_arear(radius_water) &
      & *(-(drr*drr-1._rkind)*temperature_melting_I &
      & +drr*drr*temp(iw+2)-temp(iw+3))/(drr*dr)
    if((flux_in < flux_out) .AND. (abs(flux_in-flux_out)/flux_in > 5.e-4_rkind ) ) then
      call flux_frez(flux_out-flux_in)
   else
      call flux_melt_below(flux_in-flux_out)
    endif
    if(radius_water < radius(iw+1)) then
       call quadric_temp_water
    elseif((ir/=iw) .AND. (radius_water >= radius(iw+1))) then
      temp(iw+1) = temperature_melting_I
    endif
  endif
  call flush(3)
  return
end subroutine flux_check
!--------1---------2---------3---------4---------5---------6---------7---------8
subroutine quadric_temp_water
  use global_constants, only : rkind, dt, time, year
  use global_parameters, only : temp, radius, dr
  use material_parameters, only : iw, radius_water
  use ice_parameters, only : temperature_melting_I
  implicit none
  real(rkind) :: a, b, c, dr1, dr2, drb
  c = temperature_melting_I
  drb = radius(iw+1) - radius_water
  dr1 = radius(iw+2) - radius_water
  dr2 = radius(iw+3) - radius_water
  a = (-c*(dr2-dr1) + dr2*temp(iw+2) - dr1*temp(iw+3))/(dr1*dr2*(dr1-dr2))
  b = (temp(iw+2)-c-a*dr1*dr1)/dr1
  temp(iw+1) = a*drb*drb + b*drb + c
end subroutine quadric_temp_water