!--------1---------2---------3---------4---------5---------6---------7---------8
subroutine check_cool 
!--------1---------2---------3---------4---------5---------6---------7---------8
  use global_constants, only : rkind, year, kilometer, time, Myr, dt &
    & , cool_temperature, cool_temperature_TOP, cool_temperature_BOTTOM
  use global_parameters, only : irmax, radius, temp &
    & , cooling_rate, cooled_time, temp_max
  use heat_radiogenic, only : iAl, ftime
  use energy, only : temp_old, sum_cooling_rate, isum_cooling_rate
  implicit none
  real(rkind) :: time_Myr
  integer :: i
  time_Myr = (time/year+ftime(iAl))/Myr
  do i = 0, irmax
    if(temp_old(i) < temp(i)) cycle
    if(temp_max(i) < cool_temperature) cycle
    if(temp_old(i) < cool_temperature_TOP .AND. temp(i) > cool_temperature_bottom) then
      cooling_rate(i) = (temp_old(i) - temp(i))/dt*Myr*year
      sum_cooling_rate(i) = sum_cooling_rate(i) + cooling_rate(i)
      isum_cooling_rate(i) = isum_cooling_rate(i) + 1
    endif
    if(temp_old(i) > cool_temperature .AND. temp(i) < cool_temperature) then
      cooled_time(i) = (time/year+ftime(iAl))/Myr
    endif
  enddo
  return
end subroutine check_cool
!--------1---------2---------3---------4---------5---------6---------7---------8
subroutine graph_radius 
!--------1---------2---------3---------4---------5---------6---------7---------8
  use global_constants, only : rkind, year, kilometer, time, Myr, temperature_Kelvin 
  use global_parameters, only : irmax, radius, temp &
    & , temp_max, cooling_rate, cooled_time, fname_radius
  use heat_radiogenic, only : iAl, ftime
  use energy, only : sum_cooling_rate, isum_cooling_rate
  implicit none
  real(rkind) :: time_Myr, average_cooling_rate
  integer :: i
  open(10, file = fname_radius)
  time_Myr = (time/year+ftime(iAl))/Myr
  do i = 0, irmax
    average_cooling_rate = sum_cooling_rate(i)/real(isum_cooling_rate(i))
    if(isum_cooling_rate(i) == 0) average_cooling_rate = 0._rkind
    write(10,'(f9.3,3f10.4,f9.3,e12.3,i10)') radius(i)/kilometer, temp_max(i) - temperature_Kelvin &
      & , average_cooling_rate, cooled_time(i), (radius(irmax)-radius(i))/kilometer 
  enddo
end subroutine graph_radius