!--------1---------2---------3---------4---------5---------6---------7---------8
module energy_dump 
  use global_constants, only : rkind, time, year, dt, pi
  use global_parameters, only : irmax, temp, dr, surface_area, fw
  use energy, only : lat_energy, rad_energy, time_variation, alt_energy &
    & , temp_old, alt_lat_energy, dehyd_energy
  implicit none
  private
  public :: energy_w_dump, energy_r_dump, it_start_flux_alt, it_start_water_alt
  integer :: it_start_flux_alt, it_start_water_alt
contains
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine energy_r_dump(it) 
    use global_constants, only : itmax, dt, Myr
    use global_parameters, only : fname_Enew, fname_Ecum
    use heat_radiogenic, only : iAl, ftime
    integer,intent(in) :: it
    integer :: i, its
    real(rkind) :: sum, heat_sum, use_sum, tmp, alt_sums, tv_sums, rh_sums &
      & , deh_sums, lh_sums, sh_sums, sum_sums
    real(rkind),allocatable :: al_time(:), al_alt_sum(:) &
      & , al_time_variation_sum(:), al_rad_heat_sum(:), al_dehyd_sum(:) &
      & , al_lat_heat_sum(:), al_surface_heat(:), al_alt_lat_sum(:)
    open(fw, status='unknown', form='unformatted')
    open(14, file = fname_Enew)
    open(15, file = fname_Ecum)
    allocate (al_alt_sum(it),al_time(it), al_time_variation_sum(it) &
      & , al_rad_heat_sum(it), al_lat_heat_sum(it), al_surface_heat(it) &
      & , al_alt_lat_sum(it), al_dehyd_sum(it))
    do i = 1, it
      read(fw) al_time(i), al_rad_heat_sum(i), al_time_variation_sum(i) &
        & , al_lat_heat_sum(i), al_surface_heat(i), al_alt_sum(i) &
        & , al_alt_lat_sum(i), al_dehyd_sum(i)
    enddo
    its = it_start_water_alt
    sum = 0._rkind
    do i = it_start_water_alt, it_start_flux_alt-1
      if(al_alt_sum(i) == 0._rkind) then
        cycle
      elseif(al_alt_sum(i) >= al_alt_sum(i-1)) then
        sum = sum + al_alt_sum(i)
      elseif(al_alt_sum(i) <= al_alt_sum(i-1)) then
        if(i-1==its) cycle
        sum = sum/real(i-1-its)
        al_alt_sum(its:i-1) = sum
        al_alt_lat_sum(its:i-1) = sum 
        sum = al_alt_sum(i)
        its = i
      endif
    enddo
    sum = sum/real(i-1-its)
    al_alt_sum(its:i-1) = sum 
    al_alt_lat_sum(its:i-1) = sum
    alt_sums = 0._rkind
    tv_sums=0._rkind
    rh_sums=0._rkind
    deh_sums=0._rkind
    lh_sums=0._rkind
    sh_sums=0._rkind
    do i = 1, it
      heat_sum = al_rad_heat_sum(i) - al_surface_heat(i)
      use_sum = al_time_variation_sum(i) + al_lat_heat_sum(i) &
        & - al_alt_sum(i) + al_alt_lat_sum(i)  + al_dehyd_sum(i)
      alt_sums = alt_sums + al_alt_sum(i)*dt
      tv_sums = tv_sums + al_time_variation_sum(i)*dt
      rh_sums = rh_sums + al_rad_heat_sum(i)*dt
      deh_sums = deh_sums + al_dehyd_sum(i)*dt
      lh_sums = lh_sums + al_lat_heat_sum(i)*dt + al_alt_lat_sum(i)*dt
      sh_sums = sh_sums + al_surface_heat(i)*dt
      sum_sums = rh_sums - (sh_sums+tv_sums+lh_sums+deh_sums) + alt_sums
      if(mod(i,itmax/1000)==0) then
      write(14,'(f11.7,9e12.4)') (al_time(i)/year+ftime(iAl))/Myr, heat_sum-use_sum &
        & , heat_sum, use_sum, al_rad_heat_sum(i), -al_time_variation_sum(i) &
        & , -al_lat_heat_sum(i) - al_alt_lat_sum(i), -al_surface_heat(i) &
        & , al_alt_sum(i), -al_dehyd_sum(i)
      write(15,'(f11.7,7e13.4)') (al_time(i)/year+ftime(iAl))/Myr, sum_sums, rh_sums, sh_sums &
        & , tv_sums, lh_sums, alt_sums, deh_sums
      endif
      if(i>=it_start_water_alt.AND.i<its) then
      write(14,'(f11.7,9e12.4)') (al_time(i)/year+ftime(iAl))/Myr, heat_sum-use_sum &
        & , heat_sum, use_sum, al_rad_heat_sum(i), -al_time_variation_sum(i) &
        & , -al_lat_heat_sum(i) - al_alt_lat_sum(i), -al_surface_heat(i) &
        & , al_alt_sum(i), -al_dehyd_sum(i)
      endif
    enddo
    deallocate (al_alt_sum,al_time, al_time_variation_sum, al_dehyd_sum &
      & , al_rad_heat_sum, al_lat_heat_sum, al_surface_heat, al_alt_lat_sum)
    return
  end subroutine energy_r_dump
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine energy_w_dump 
    use global_constants, only : dt, pi
    use global_parameters, only : dr, surface_area, surface_arear
    use heat_parameters, only : thermal_conductivity, thermometric_diffusivity
    use material_parameters, only : mn, iw, ir, material_select&
      & , ir_old, iw_old, ii, ice_reach_rock,radius_rock
    use energy, only : flux_rock
    implicit none
    integer :: i
    real(rkind) :: idt, volume
    real(rkind) :: time_variation_sum, rad_heat_sum, lat_heat_sum &
      & , surface_heat, alt_sum, alt_lat_sum, dehyd_sum, surface_heat_rock
    open(fw, status='unknown', form='unformatted')
    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 &
        & *thermal_conductivity(mn(0))/thermometric_diffusivity(mn(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 = -thermal_conductivity(6)*surface_area(irmax)*0.5_rkind &
      & *(3._rkind*temp(irmax)-4._rkind*temp(irmax-1)+temp(irmax-2))/dr
    if(.NOT.ice_reach_rock) then 
      surface_heat_rock = flux_rock*surface_arear(radius_rock)
      surface_heat = surface_heat + surface_heat_rock
    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 &
        & *thermal_conductivity(mn(i))/thermometric_diffusivity(mn(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)
    enddo
    write(fw) time, rad_heat_sum, time_variation_sum &
      & , lat_heat_sum, surface_heat, alt_sum, alt_lat_sum, dehyd_sum
  return
  end subroutine energy_w_dump
end module energy_dump