!--------1---------2---------3---------4---------5---------6---------7---------8
module energy 
  use global_constants, only : rkind
  use global_parameters, only : irmax
  implicit none
  private
  public :: energy_check, lat_energy, rad_energy, set_energy, energy_init &
    & , energy_rock, flux_rock, time_variation, temp_old &
    & , alt_energy, alt_lat_energy, dehyd_energy &
    & , sum_cooling_rate, isum_cooling_rate
  real(rkind) :: temp_old(0:irmax), lat_energy(0:irmax), rad_energy(0:irmax) &
    & , energy_rock(0:1), flux_rock, time_variation(0:irmax), temp_max(0:irmax) &
    & , alt_energy(0:irmax), alt_lat_energy(0:irmax), dehyd_energy(0:irmax) &
    & , sum_cooling_rate(0:irmax)
  integer :: isum_cooling_rate(0:irmax)
contains
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine energy_init 
    temp_old(0:irmax) = 0._rkind
    sum_cooling_rate(0:irmax) = 0._rkind
    isum_cooling_rate(0:irmax) = 0
    temp_max(0:irmax) = 0._rkind
    lat_energy(0:irmax) = 0._rkind
    rad_energy(0:irmax) = 0._rkind
    energy_rock(0:1) = 0._rkind 
    time_variation(0:irmax) = 0._rkind
    alt_energy(0:irmax) = 0._rkind
    alt_lat_energy(0:irmax) = 0._rkind
    dehyd_energy(0:irmax) = 0._rkind
  end  subroutine energy_init
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine set_energy 
  use global_constants, only : volume_ratio_rock, time, dt
  use global_parameters, only : dr, surface_area, halfdr
  use material_parameters, only : mn, ice_water_rock, rock_water, ir_old, iw &
    & , material_select, ice_reach_rock
  integer :: i
  real(rkind) :: heat_decay
  rad_energy(0) = heat_decay(time+dt,0)
  lat_energy(0) = 0._rkind
  alt_energy(0) = 0._rkind
  alt_lat_energy(0) = 0._rkind
  dehyd_energy(0) = 0._rkind
  do i = 1, irmax-1
    lat_energy(i) = 0._rkind
    alt_energy(i) = 0._rkind
    alt_lat_energy(i) = 0._rkind
    dehyd_energy(i) = 0._rkind
    if(.NOT.ice_reach_rock) then 
      if(iw /= 0 .AND. (i==ir_old.OR.i==ir_old+1)) cycle
    endif
    call material_select(i)
    rad_energy(i) = heat_decay(time,i)*surface_area(i)*dr 
  enddo
  alt_energy(irmax) = 0._rkind
  lat_energy(irmax) = 0._rkind
  alt_lat_energy(irmax) = 0._rkind
  dehyd_energy(irmax) = 0._rkind
  rad_energy(irmax) = heat_decay(time,irmax)&
    & *surface_area(irmax)*halfdr 
  end subroutine set_energy
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine energy_check 
  use global_constants, only : rkind, density, time, year, dt, pi, imr, Myr, Rocky_Core
  use global_parameters, only : irmax, temp, dr, surface_area, surface_arear &
    & , fname_E, fname_Eall, radius
  use heat_parameters, only : thermal_conductivity, thermometric_diffusivity
  use heat_radiogenic, only : iAl, ftime
  use material_parameters, only : mn, iw, ir, material_select&
    & , ir_old, iw_old, ii, radius_rock, radius_ice, ice_reach_rock
  use ice_parameters, only : temperature_melting_I
  implicit none
  integer :: i
  real(rkind) :: energy_sum,  volume, time_variation_sum, heat_sum, use_sum &
    & , rad_heat_sum, lat_heat_sum,idt, surface_heat, alt_sum, alt_lat_sum &
    & , dehyd_sum, surface_heat_rock &
    & , Tthermal_conductivity, thermometric_conductivity, Tspecific_heat, Tdensity
  open(12, file = fname_E)
  open(13, file = fname_Eall)
  idt = 1._rkind/dt
  volume = 4._rkind/3._rkind*pi*dr*dr*dr 
  call material_select(0)
  time_variation_sum = (temp(0)-temp_old(0))*volume &
      & *Tspecific_heat(0)*Tdensity(0)*idt
  rad_heat_sum = rad_energy(0)*volume
  lat_heat_sum = lat_energy(0) 
  alt_lat_sum = alt_lat_energy(0)
  alt_sum = alt_energy(0)
  dehyd_sum = dehyd_energy(0)
  surface_heat = -Tthermal_conductivity(irmax)*surface_area(irmax)*0.5_rkind &
    & *(3._rkind*temp(irmax)-4._rkind*temp(irmax-1)+temp(irmax-2))/dr
  if(Rocky_Core) then
  if(.NOT.ice_reach_rock) then 
    flux_rock = 0._rkind
    surface_heat_rock = flux_rock*surface_arear(radius_rock)
    surface_heat = surface_heat + surface_heat_rock
  endif
  endif
  do i = 1, irmax-1
    if(.NOT.ice_reach_rock) then 
      if(iw /= 0 .AND. (i==ir_old.OR.i==ir_old+1)) cycle
    endif
    call material_select(i)
    volume = surface_area(i)*dr
    time_variation(i) = (temp(i)-temp_old(i))*volume &
      & *Tspecific_heat(i)*Tdensity(i)*idt
  enddo
  time_variation(irmax) = 0._rkind
  do i = 1, irmax
    rad_heat_sum = rad_heat_sum + rad_energy(i)
    lat_heat_sum = lat_heat_sum + lat_energy(i) 
    alt_lat_sum = alt_lat_sum + alt_lat_energy(i)
    alt_sum = alt_sum + alt_energy(i)
    dehyd_sum = dehyd_sum + dehyd_energy(i)
    time_variation_sum = time_variation_sum + time_variation(i)
    write(13,'(f11.7,4i4,6e11.3)') (time/year+ftime(iAl))/Myr, i, ir, iw, mn(i) &
      & , rad_energy(i), time_variation(i), lat_energy(i), temp(i)-temp_old(i) &
      & , alt_energy(i), dehyd_energy(i)
  enddo
  volume = 4._rkind/3._rkind*pi*radius(irmax)**3._rkind
  energy_sum = rad_heat_sum - time_variation_sum - lat_heat_sum - surface_heat&
    & + alt_sum - alt_lat_sum - dehyd_sum
  heat_sum = rad_heat_sum-surface_heat
  use_sum = time_variation_sum+lat_heat_sum-alt_sum+alt_lat_sum + dehyd_sum
  write(12,'(f11.7,9e12.4)') (time/year+ftime(iAl))/Myr, energy_sum/volume &
    & , heat_sum/volume, use_sum/volume, rad_heat_sum/volume, -time_variation_sum/volume &
    & , -lat_heat_sum/volume-alt_lat_sum, -surface_heat/volume, alt_sum, -dehyd_sum
  call flush(12)
  call flush(13)
  end subroutine energy_check
end module energy