analysis.f90

Go to the documentation of this file.

!> @file analysis.f90
!!
!! The error file code for this file is ***W10***.
!! @brief Module \ref analysis with various analysis subroutines
!!
 
!> Contains various routines for analysis and diagnostic output
#include "macros.h"
module analysis
 
  use parameter_class, only: nrdiag_every,timescales_every
  use global_class,    only: net_size, ihe4, ineu, ipro
  use file_handling_class
 
  implicit none
 
  real(r_kind),dimension(:),allocatable :: sn                                    !< array of n-separation energies
  real(r_kind)                          :: sn_ave                                !< average n-separation energy
  integer                               :: mainout_unit                          !< file descriptor for the main analysis output file
  integer                               :: nrdiag_unit                           !< diagnostic output inside the Newton-Rhapson loop
  integer                               :: track_unit                            !< file id for tracked nuclei
  integer                               :: toplist_unit                          !< file id for top nuclear energy contributors
  integer                               :: nu_loss_gain_unit                     !< file id for neutrino loss/gain
  integer                               :: cl_start, cl_end, cl_rate, cl_cmax    !< system clock variables
  integer,private                       :: tsfile                                !< timescale file unit
  integer,private                       :: nuc_heat_file                         !< nuclear heating file unit
  real(r_kind)                          :: tsnp, t9snp, rosnp                    !< for snapshot printing triggers
  integer                               :: snapcnt                               !< snapshot count
  integer                               :: flowcnt                               !< flow printing counter
  real,dimension(:),allocatable         :: snapshot_time                         !< point in time for custom snapshot
  integer                               :: snapshot_amount                       !< Amount of custom snapshots
  integer,parameter                     :: toplist_amount=10                     !< Amount of reactions in the toplist
 
  !
  ! Public and private fields and methods of the module
  !
  public:: &
     analysis_init, output_initial_step, output_final_step, output_iteration, &
     analysis_finalize, output_nr_diagnostic, output_final_stats, finab
 
  private:: &
     output_track_nuclei, output_mainout, output_snapshot, printsnap, &
     calc_nseparation_energy, calc_nuclear_heating, calc_av_timescales, &
     finab_sort, output_nu_loss
 
contains
 
!>
!! Open files, write headers, allocate space etc.
!!
subroutine analysis_init()
  use parameter_class
  use global_class,         only: net_names,track_nuclei_nr,track_nuclei_indices,nreac
  use error_msg_class,      only: int_to_str
  use benam_class,          only: findaz
  use ls_timmes_eos_module
#ifdef USE_HDF5
  use hdf5_module,          only: init_hdf5_module
#endif
  implicit none
 
  ! neutron separation energies
  real(r_kind)  :: val                        !< aux variable for reading Sn file
  integer       :: ai, ni, zi, snfile, fstat  !< temporary variables needed for reading the Sn file
  character     :: nl, buf(30)                !< newline character, buffer
  integer       :: i,k
  character*10,dimension(track_nuclei_nr) :: tmp_nucleinames !< Name of the nuclei that will be tracked
  character*1500 :: head_toplist
  character*5    :: help_char
  info_entry("analysis_init")
 
#ifdef USE_HDF5
   ! Only initialize if it is really possible to execute it
   ! (i.e., correct compiler h5fc)
   call init_hdf5_module()
#endif
 
  cl_cmax= 100000
  snapcnt= 1
  flowcnt= 1
  nl = new_line('A')
 
  ! main output file: create, write header
  if (mainout_every .gt. 0) then
     mainout_unit= open_outfile("mainout.dat")
     write(mainout_unit,'(3A)') &
            '# 1:iteration 2:time[s], 3:T[GK], 4:rho[g/cm3], 5:Ye '//nl &
          , '# 6:R[km], 7:Y_n, 8:Y_p, 9:Y_alpha, 10:Y_lights '//nl  &
          , '# 11:Y_heavies 12:<Z> 13:<A> 14:entropy [kB/baryon] (15:Sn [MeV])'
  endif
 
  ! read neutron separation energies
  if (calc_nsep_energy) then
     snfile= open_outfile(nsep_energies_file)
 
     allocate(sn(net_size))
     sn(1:net_size) = 0d0
     read(snfile, *) buf ! skip the first line
     do
        read(snfile, *, iostat=fstat) zi, ni, ai, val
        if (fstat /= 0) exit
        k = findaz(ai, zi)
        if (k>0) then
           sn(k) = val
        endif
     enddo
     call close_io_file(snfile,nsep_energies_file)
  endif
 
  ! Open the timescale file and write a header
  if (timescales_every .gt. 0) then
    tsfile = open_outfile('timescales.dat')
    write(tsfile, '((A1,(A22,3x),*(A23,3x)))') &
       "#","time[s]","Temperature [GK]","tau_ga [s]","tau_ag [s]","tau_ng [s]","tau_gn [s]",   &
       "tau_pg [s]", "tau_gp [s]", "tau_np [s]", "tau_pn [s]", &
       "tau_an [s]", "tau_na [s]", "tau_ap [s]", "tau_pa [s]", "tau_beta [s]", "tau_alpha [s]", &
       "tau_nfiss [s]", "tau_sfiss [s]", "tau_bfiss [s]"
  end if
 
  ! Open the track_nuclei file and write a header
  if (track_nuclei_every .gt. 0) then
     ! Convert indices to names
    do i=1, track_nuclei_nr
      tmp_nucleinames(i) = "Y( "//net_names(track_nuclei_indices(i))//" )"
    end do
    ! Write the header
    track_unit = open_outfile('tracked_nuclei.dat')
    write(track_unit, '((A,20X,*(A,18X)))') "# time[s]",tmp_nucleinames(:)
  end if
 
  ! Write the header for the toplist
  if (top_engen_every .gt. 0) then
    toplist_unit= open_outfile("toplist.dat")
 
    ! Give some hint on what is what
    write (toplist_unit,"('#',31X,A,434X,A,413X,'|',5X,A,185X,A)") "Top nuclear energy producers for each reaction given in erg/g/s", &
                                                               "Top nuclear energy consumers for each reaction given in erg/g/s", &
                                                               "Top nuclear energy producers for each isotope given in erg/g/s", &
                                                               "Top nuclear energy consumers for each isotope given in erg/g/s"
 
    ! Construct the rather complicated header of the file
    head_toplist ="#    Time[s]   Etot[erg/g/s]    Reaction_1"
    do i=2,toplist_amount
        head_toplist = trim(adjustl(head_toplist))//"                       Reaction_"//trim(adjustl(int_to_str(i)))
    end do
    head_toplist = trim(adjustl(head_toplist))//"                           Ereac_1"
    do i=2,toplist_amount
        head_toplist = trim(adjustl(head_toplist))//"        Ereac_"//trim(adjustl(int_to_str(i)))
    end do
    head_toplist = trim(adjustl(head_toplist))//"         Reaction_1"
    do i=2,toplist_amount
        head_toplist = trim(adjustl(head_toplist))//"                       Reaction_"//trim(adjustl(int_to_str(i)))
    end do
    head_toplist = trim(adjustl(head_toplist))//"                          Ereac_1"
    do i=2,toplist_amount
        head_toplist = trim(adjustl(head_toplist))//"        Ereac_"//trim(adjustl(int_to_str(i)))
    end do
    head_toplist = trim(adjustl(head_toplist))//"         |      Iso1"
    do i=2,toplist_amount
        help_char = "Iso"//trim(adjustl(int_to_str(i)))
        head_toplist = trim(adjustl(head_toplist))//"   "//adjustr(help_char)
    end do
    head_toplist = trim(adjustl(head_toplist))//"        Eiso_1"
    do i=2,toplist_amount
        head_toplist = trim(adjustl(head_toplist))//"         Eiso_"//trim(adjustl(int_to_str(i)))
    end do
    head_toplist = trim(adjustl(head_toplist))//"           Iso1"
    do i=2,toplist_amount
        help_char = "Iso"//trim(adjustl(int_to_str(i)))
        head_toplist = trim(adjustl(head_toplist))//"   "//adjustr(help_char)
    end do
    head_toplist = trim(adjustl(head_toplist))//"        Eiso_1"
    do i=2,toplist_amount
        head_toplist = trim(adjustl(head_toplist))//"         Eiso_"//trim(adjustl(int_to_str(i)))
    end do
    ! Finally write it out
    write (toplist_unit,"(A)") trim(adjustl(head_toplist))
  end if
 
 
  if (nrdiag_every.gt.0) then
     nrdiag_unit= open_outfile("nrdiag.dat")
     write(nrdiag_unit,'(11A)') &
            '# Newton-Raphson diagnostic output'//nl               &
          , '# 1 : global iteration count (cnt)'//nl               &
          , '# 2 : global time [s]'//nl                            &
          , '# 3 : temperature T9 [GK]'//nl                        &
          , '# 4 : adapt stepsize loop counter (k)'//nl            &
          , '# 5 : NR loop counter (nr_count)'//nl                 &
          , '# 6 : global timestep [s]'//nl                        &
          , '# 7 : mass conservation error'//nl                    &
          , '# 8 : maximal abundance change (eps)'//nl             &
          , '# 9 : most rapidly evolving isotope (epsl)'//nl &
          , '# 10: abundance of the isotope epsl, y(epsl)'
  end if
 
  ! Open the energy generation file and write a header
  if (engen_every .gt. 0) then
    nuc_heat_file = open_outfile('generated_energy.dat')
    ! Header
    write(nuc_heat_file, '((A,5X))') "# Generated Energy for each nuclear reaction type given in erg/g/s"
    write(nuc_heat_file, '(17(A23,3X))') &
       "#               time[s]", "Engen(Total)", "S_src","Engen(weak)","Engen(n,g)",&
       "Engen(p,g)","Engen(a,g)","Engen(n,p)","Engen(n,a)", &
       "Engen(p,a)", "Engen(fiss)"
  end if
 
 
 
  if ((nu_loss_every .gt. 0)) then
    nu_loss_gain_unit = open_outfile('nu_loss_gain.dat')
     write(nu_loss_gain_unit, '(A)') "# File containing the neutrino loss/gain for the nuclear heating."
     write(nu_loss_gain_unit, '(A)') "# Negative values mean that energy is added to the system."
     write(nu_loss_gain_unit, '(A)') "# Neutrino energies are given in Erg/g/s."
     write(nu_loss_gain_unit, '((A1,(A18,3x),*(A19,3x)))') &
        "#","time[s]","T[GK]", "rho[g/cm3]", "R[km]", "nu_total", "nu_beta",&
        "nu_nuheat", "nu_thermal"
  end if
 
 
  info_exit("analysis_init")
 
end subroutine analysis_init
 
 
!>
!! Output first step
!!
!! This routine calls output_iteration at the first step
!!
!! @see output_final_step
subroutine output_initial_step(ctime,T9,rho_b,entropy,Rkm,Y,pf)
  implicit none
  real(r_kind),intent(in) :: ctime    !< actual time [s]
  real(r_kind),intent(in) :: t9       !< initial temperature [GK]
  real(r_kind),intent(in) :: rho_b    !< initial density [gcc]
  real(r_kind),intent(in) :: entropy  !< entropy [kB/nucleon]
  real(r_kind),intent(in) :: rkm      !< Radius [km]
  real(r_kind),dimension(net_size),intent(in)   :: y !< array of abundances (Y_i)
  real(r_kind),dimension(0:net_size),intent(inout) :: pf !< partition functions
 
  info_entry('output_initial_step')
  call output_iteration(0,ctime,t9,rho_b,entropy,rkm,y,pf)
  info_exit('output_initial_step')
 
end subroutine output_initial_step
 
 
!>
!! Output final step
!!
!! This routine calls output_iteration at the last step
!!
!! @see output_final_step
!!
!! \b Edited: OK 26.08.2017
subroutine output_final_step(cnt,ctime,T9,rho_b,entropy,Rkm,Y,pf)
  use parameter_class, only: out_every,snapshot_every,flow_every,mainout_every, &
                             track_nuclei_every,termination_criterion, &
                             timescales_every,nu_loss_every,top_engen_every, &
                             engen_every, h_snapshot_every, h_mainout_every, &
                             h_track_nuclei_every, h_timescales_every, h_flow_every, &
                             h_nu_loss_every, h_engen_every
  implicit none
  integer,intent(in)      :: cnt      !< current iteration counter
  real(r_kind),intent(in) :: ctime    !< actual time [s]
  real(r_kind),intent(in) :: t9       !< initial temperature [GK]
  real(r_kind),intent(in) :: rho_b    !< initial density [gcc]
  real(r_kind),intent(in) :: entropy  !< entropy [kB/nucleon]
  real(r_kind),intent(in) :: rkm      !< Radius [km]
  real(r_kind),dimension(net_size),intent(in)   :: y  !< array of abundances (Y_i)
  real(r_kind),dimension(0:net_size),intent(inout) :: pf !< partition functions
 
  info_entry('output_final_step')
 
  ! hack output frequency parameters to output the final step out-of-order
  if (out_every.gt.0)          out_every= 1
  if (mainout_every.gt.0)      mainout_every= 1
  if (snapshot_every.gt.0)     snapshot_every= 1
  if (flow_every.gt.0)         flow_every= 1
  if (timescales_every.gt.0)   timescales_every= 1
  if (track_nuclei_every.gt.0) track_nuclei_every= 1
  if (nu_loss_every.gt.0)      nu_loss_every= 1
  if (top_engen_every.gt.0)    top_engen_every= 1
  if (engen_every.gt.0)        engen_every= 1
  ! Same for hdf5 output
  if (h_snapshot_every.gt.0)     h_snapshot_every= 1
  if (h_mainout_every.gt.0)      h_mainout_every= 1
  if (h_track_nuclei_every.gt.0) h_track_nuclei_every= 1
  if (h_timescales_every.gt.0)   h_timescales_every= 1
  if (h_flow_every.gt.0)         h_flow_every= 1
  if (h_nu_loss_every.gt.0)      h_nu_loss_every= 1
  if (h_engen_every.gt.0)        h_engen_every= 1
 
  print '("---------------")'
  call output_iteration(cnt,ctime,t9,rho_b,entropy,rkm,y,pf)
  print '("============================")'
  select case(termination_criterion)
  case(0)
     print'(A)',"End of trajectory file reached."
  case(1)
     print'(A)',"Final time reached."
  case(2)
     print'(A)',"Final temperature reached."
  case(3)
     print'(A)',"Final density reached."
  case(4)
     print'(A)',"Freeze-out reached."
  endselect
 
  info_exit('output_final_step')
 
end subroutine output_final_step
 
 
!>
!! Periodic output of analysis data for current iteration
!!
!! Prints various outputs depending on the input parameter.
!!
!! \b Edited:
!!       - MR : 19.01.21
!! .
subroutine output_iteration(cnt,ctime,T9,rho_b,entropy,Rkm,Y,pf)
  use parameter_class,  only: out_every,snapshot_every,mainout_every, &
                              track_nuclei_every,heating_mode, &
                              custom_snapshots,flow_every,engen_every,&
                              top_engen_every,h_snapshot_every,h_mainout_every,&
                              h_custom_snapshots,h_track_nuclei_every,&
                              h_timescales_every,h_flow_every,h_engen_every,&
                              nu_loss_every, h_nu_loss_every, unit, use_thermal_nu_loss
  use nucstuff_class,   only: inter_partf
  use nuclear_heating,  only: output_nuclear_heating,&
                              qdot_nu, qdot_bet, qdot_th
  use thermal_neutrino_module, only: thermal_neutrinos
  use single_zone_vars, only: evolution_mode, stepsize
  use global_class,     only: isotope
  use flow_module,      only: flowcalc, flows
#ifdef USE_HDF5
  use hdf5_module
#endif
  implicit none
 
  integer,intent(in)      :: cnt      !< current iteration counter
  real(r_kind),intent(in) :: ctime    !< actual time [s]
  real(r_kind),intent(in) :: t9       !< temperature [GK]
  real(r_kind),intent(in) :: rho_b    !< density [gcc]
  real(r_kind),intent(in) :: entropy  !< entropy [kB/nucleon]
  real(r_kind),intent(in) :: rkm      !< length scale [km]
  real(r_kind),dimension(net_size),intent(in)      :: y  !< array of abundances (Y_i)
  real(r_kind),dimension(0:net_size),intent(inout) :: pf !< partition functions
  integer                 :: i              !< Loop variable
  real(r_kind)            :: ysum,yasum,abar,ye
 
  info_entry("output_iteration")
 
  ! screen output
  if (out_every.gt.0) then
     if (mod(cnt,50*out_every)==0) then
        print *
        print '(A)',"#-- 1:it 2:time[s] 3:T9[GK] 4:rho[g/cm3] &
                     5:entropy[kB/baryon]"
     endif
     if (mod(cnt,out_every) .eq. 0) then
        print '(I8,X,ES14.7,X,F10.7,X,ES14.7,X,F12.7)', &
              cnt,ctime,t9,rho_b,entropy
     endif
 
     if(mod(cnt,out_every).eq.0 .and. (heating_mode .eq. 1)) then
        if (verbose_level.ge.2) then
            call output_nuclear_heating(cnt,ctime)
        end if
     endif
  endif
 
  ! This is only done if the correct compiler (h5fc) was used
#ifdef USE_HDF5
  ! Store snapshots in hdf5
  if (h_snapshot_every.gt.0) then
     if (mod(cnt,h_snapshot_every).eq. 0) then
        ! Ye = sum(isotope(1:net_size)%p_nr*Y(1:net_size))
        call extend_snaps(ctime,y)
     end if
  end if
 
  ! Store mainout in hdf5
  if (h_mainout_every.gt.0) then
     if (mod(cnt,h_mainout_every).eq. 0) then
        call extend_mainout(cnt,ctime,t9,rho_b,entropy,rkm,y)
     end if
  end if
 
  ! Snapshot output for specific times
  if (h_custom_snapshots) then
     ! Test if a snapshot should be done
     h_make_snap: do i=1,snapshot_amount
        ! Time is within certain tolerance? Do the snapshot!
        if (abs(ctime-snapshot_time(i)) .le. 1e-20) then
           call extend_snaps(ctime,y)
           exit h_make_snap
        end if
     enddo h_make_snap
  end if
 
  ! Output the tracked nuclei
  if (h_track_nuclei_every .gt. 0) then
      if (mod(cnt,h_track_nuclei_every).eq.0) then
         call extend_track_nuclei(ctime,y)
      end if
  end if
 
  ! Output average timescales
  if (h_timescales_every .gt. 0) then
     if ((mod(cnt,h_timescales_every).eq.0) .and. (evolution_mode .ne. em_nse)) then
        ! calculate averaged timescales of various types of reactions
        ! calculate partition function first
        call inter_partf (t9, pf)
        call calc_av_timescales(ctime,t9,y,.true.)
     end if
  end if
 
  ! Output the flow
  if ((h_flow_every .gt. 0) .and. (cnt .ne. 0)) then
     if ((mod(cnt,h_flow_every).eq.0)) then
        call flowcalc(y)
        call extend_flow(ctime,t9,rho_b,flows,size(flows),y)
     end if
  end if
 
  ! Calculate and output generated energy
  if (h_engen_every .gt. 0) then
     if (mod(cnt,h_engen_every).eq.0) then
        call calc_nuclear_heating(.true., .true., .false.)
     end if
  end if
 
  ! Calculate and output generated energy
  if ((h_nu_loss_every .gt. 0)) then
     if (mod(cnt,h_nu_loss_every).eq.0) then
 
        if ((heating_mode .eq. 0) .and. (use_thermal_nu_loss)) then
            ysum= sum(y(1:net_size))
            yasum= sum(y(1:net_size)*isotope(1:net_size)%mass)
            ye = sum(isotope(1:net_size)%p_nr*y(1:net_size))
            abar= yasum / ysum
            call thermal_neutrinos(abar, ye, t9, rho_b, stepsize, qdot_th)
            qdot_th = qdot_th * unit%hix ! erg/g -> MeV/baryon
        else
            qdot_th = 0d0
        end if
 
        call extend_nu_loss_gain(ctime,t9,rho_b,rkm,qdot_th/unit%hix/stepsize,&
                                                    qdot_nu/unit%hix/stepsize,&
                                                    qdot_bet/unit%hix/stepsize)
     end if
  end if
 
#endif
 
  ! main data log
  if (mainout_every .gt.0) then
     if (mod(cnt,mainout_every).eq.0) then
        call output_mainout(cnt,ctime,t9,rho_b,entropy,rkm,y)
     end if
  end if
 
  ! Iterative Snapshot and flow output
  if(snapshot_every.gt.0) then
     if(mod(cnt,snapshot_every).eq.0) then
        call output_snapshot(ctime,t9,rho_b,y)
     end if
  end if
 
  if((flow_every.gt.0) .and. (cnt .ne. 0)) then
     if(mod(cnt,flow_every).eq.0) then
        call output_flow(ctime,t9,rho_b,y)
     end if
  end if
 
  if ((nu_loss_every .gt. 0)) then
     if (mod(cnt,nu_loss_every).eq.0) then
 
        if ((heating_mode .eq. 0) .and. (use_thermal_nu_loss)) then
            ysum= sum(y(1:net_size))
            yasum= sum(y(1:net_size)*isotope(1:net_size)%mass)
            ye = sum(isotope(1:net_size)%p_nr*y(1:net_size))
            abar= yasum / ysum
            call thermal_neutrinos(abar, ye, t9, rho_b, stepsize, qdot_th)
            qdot_th = qdot_th * unit%hix ! erg/g -> MeV/baryon
        else
            qdot_th = 0d0
        end if
 
        call output_nu_loss(ctime,t9,rho_b,rkm)
     end if
  end if
 
 
  ! Snapshot output for specific times
  if (custom_snapshots) then
     ! Test if a snapshot should be done
     make_snap: do i=1,snapshot_amount
        ! Time is within certain tolerance? Do the snapshot!
        if (abs(ctime-snapshot_time(i)) .le. 1e-20) then
           call output_snapshot(ctime,t9,rho_b,y)
           exit make_snap
        end if
     enddo make_snap
  end if
 
  if (timescales_every .gt. 0) then
     if ((mod(cnt,timescales_every).eq.0) .and. (evolution_mode .ne. em_nse)) then
        ! calculate averaged timescales of various types of reactions
        ! calculate partition function first
        call inter_partf (t9, pf)
        call calc_av_timescales(ctime,t9,y)
     end if
  end if
 
 ! Output the tracked nuclei
  if (track_nuclei_every .gt. 0) then
     if (mod(cnt,track_nuclei_every).eq.0) then
        call output_track_nuclei(ctime,y)
     end if
  end if
 
  ! Calculate generated energy
  if (engen_every .gt. 0) then
     if (mod(cnt,engen_every).eq.0) then
        call calc_nuclear_heating(.false., .true., .false.)
     end if
  end if
 
    ! Calculate generated energy
  if (top_engen_every .gt. 0) then
    if (mod(cnt,top_engen_every).eq.0) then
       call calc_nuclear_heating(.false., .false., .true.)
    end if
  end if
 
 
  info_exit("output_iteration")
end subroutine output_iteration
 
 
!> Output the abundances of specific nuclei
!!
!! First column is always the time, followed by the abundances of
!! nuclei given in "track_nuclei_file". An example of the file looks like:
!! \l_file{
!! # time[s]                    Y(     n )                  Y(   o24 )                  Y(   f24 )
!!  0.0000000000000000E+00      0.0000000000000000E+00      2.0833333333333332E-02      0.0000000000000000E+00
!!  1.9999999999999999E-20      8.5150206007805894E-19      2.0833333333333332E-02      1.8661570495808601E-21
!!  5.9999999999999994E-20      2.5545061802341771E-18      2.0833333333333332E-02      5.5984711487425806E-21
!!  1.3999999999999998E-19      5.9605144205464133E-18      2.0833333333333332E-02      1.3063099347066021E-20
!! ...}
!!
!! @author: Moritz Reichert
!!
!! \b Edited:
!!      - 09.06.17
!! .
subroutine output_track_nuclei(ctime,Y)
   use global_class,       only: track_nuclei_indices,track_nuclei_nr
   implicit none
   real(r_kind),intent(in)                          :: ctime    !< Actual time
   real(r_kind),dimension(net_size),intent(in)      :: y        !< array of abundances (Y_i)
   real(r_kind),dimension(track_nuclei_nr)          :: tmp_y    !< Stores the abundance of the tracked nuclei
   integer                                          :: i        !< Loop variable
 
   info_entry("output_track_nuclei")
   ! Store the abundance temporary
   do i=1, track_nuclei_nr
      tmp_y(i) = y(track_nuclei_indices(i))
      if (tmp_y(i) .lt. 1d-99) tmp_y(i) = 0.
   end do
   ! Output
   write(track_unit,'(*(es23.16,5x))') ctime,tmp_y(:)
   info_exit("output_track_nuclei")
 
end subroutine output_track_nuclei
 
 
!> Output the energy that enters and leaves the system
!!
!! Outputs the energy that enters and leaves the system when
!! nuclear_heating is turned on.
!! This output is controlled by the parameter \ref nu_loss_every .
!! The output is written to the file "nu_loss_gain_file".
!! An example of the file could look like:
!! \file{
!! \# File containing the neutrino loss/gain for the nuclear heating.
!! \# Negative values mean that energy is added to the system.
!! \# Neutrino energies are given in MeV/baryon/s.
!! \#           time[s]                 T[GK]            rho[g/cm3]                 R[km]              nu_total  ...
!!         1.10385E+00           1.00000E+01           4.89334E+07           3.55626E+02           0.00000E+00  ...
!!         1.10385E+00           1.00000E+01           4.89334E+07           3.55626E+02          -5.05524E+00  ...
!!         1.10385E+00           1.00000E+01           4.89334E+07           3.55626E+02          -5.05524E+00  ...
!!         1.10385E+00           1.00000E+01           4.89334E+07           3.55626E+02          -5.05524E+00  ...
!! ...}
!! @author Moritz Reichert
!! @date   12.04.2023
subroutine output_nu_loss(ctime,temp,rho,Rkm)
    use single_zone_vars,only: stepsize
    use nuclear_heating, only: qdot_bet, qdot_nu, qdot_th
    use parameter_class, only: unit
    implicit none
    real(r_kind),intent(in) :: ctime !< actual time (s)
    real(r_kind),intent(in) :: temp  !< temperature (GK)
    real(r_kind),intent(in) :: rho   !< density (g/cm3)
    real(r_kind),intent(in) :: rkm   !< radius (km)
    ! local variables
    real(r_kind) :: qdot_tot
 
    qdot_tot = qdot_bet + qdot_nu + qdot_th
 
    write(nu_loss_gain_unit, '(*(es19.5,3x))') ctime, temp, rho, rkm, &
                                               qdot_tot/unit%hix/stepsize, &
                                               qdot_bet/unit%hix/stepsize, &
                                               qdot_nu/unit%hix/stepsize,  &
                                               qdot_th/unit%hix/stepsize
 
end subroutine output_nu_loss
 
 
!>
!! Output Mainout file
!!
!! This file contains all basic quantities such as neutron abundance Yn,
!! time, temperature, density, Abar, and much more. An example file
!! looks like:
!! \l_file{
!! # 1:iteration 2:time[s], 3:T[GK], 4:rho[g/cm3], 5:Ye
!! # 6:R[km], 7:Y_n, 8:Y_p, 9:Y_alpha, 10:Y_lights
!! # 11:Y_heavies 12:<Z> 13:<A> 14:entropy [kB/baryon] (15:Sn [MeV])
!!    0  0.0000000000000000E+00   1.0964999999999998E+01   8.7095999999999795E+12   3.4880000000001680E-02   4.9272999999999975E+01   8.7430471311642388E-01   1.8154340816658842E-15   1.8313910002116333E-13   8.4293473893924100E-09   1.6154109228777289E-03   3.9820982195819934E-02   1.1416565996507526E+00   9.4638587929828290E-03
!!   10  4.0940000000000001E-09   1.0965038071994945E+01   8.7093701186135752E+12   3.4880006806212963E-02   4.9273176232580205E+01   8.7430472419630612E-01   1.8157922458135375E-15   1.8317495847108614E-13   8.4302150496944692E-09   1.6154086944791908E-03   3.9820989563731306E-02   1.1416565881127740E+00   9.4639764072160116E-03
!! ...}
!!
!!
subroutine output_mainout(cnt,ctime,T9,rho_b,entropy,Rkm,Y)
  use parameter_class, only: calc_nsep_energy
  use global_class,    only: isotope
  implicit none
 
  integer,intent(in)      :: cnt                         !< current iteration counter
  real(r_kind),intent(in) :: ctime                       !< current time
  real(r_kind),intent(in) :: t9                          !< temperature [GK]
  real(r_kind),intent(in) :: rho_b                       !< density [gcc]
  real(r_kind),intent(in) :: entropy                     !< entropy [kB/nucleon]
  real(r_kind),intent(in) :: rkm                         !< length scale [km]
  real(r_kind),dimension(net_size),intent(in)      :: y  !< array of abundances (Y_i)
  !
  real(r_kind) :: ye,yneu,ypro,yhe4,ylight,yheavies,ysum,yasum,yzsum,abar,zbar
 
  ! calculate neutron separation energy
  if (calc_nsep_energy) then
     call calc_nseparation_energy(y)
  endif
 
  ! output diagnostics
  ye = sum(isotope(1:net_size)%p_nr*y(1:net_size))
  ylight=   max(1.0d-99,sum(y(ipro+1:ihe4-1)))
  yheavies= 0.0
  if (ihe4.ge.1) yheavies= max(1.0d-99,sum(y(ihe4+1:net_size)))
 
  ysum  = sum(y(1:net_size))
  yasum = sum(y(1:net_size)*isotope(1:net_size)%mass)
  yzsum = sum(y(1:net_size)*isotope(1:net_size)%p_nr)
  abar  = yasum/ysum
  zbar  = yzsum/ysum
  !zbar = yzsum/sum(Y(2:net_size)) !as neutrons have Z=0
 
  yneu= 0.0
  if (ineu.ge.1) yneu= max(1.0d-99,y(ineu))
 
  ypro= 0.0
  if (ipro.ge.1) ypro= max(1.0d-99,y(ipro))
 
  yhe4= 0.0
  if (ihe4.ge.1) yhe4= max(1.0d-99,y(ihe4))
 
  if (calc_nsep_energy) then
    write (mainout_unit,'(i8,X,100(es23.16,2X))') cnt,ctime,t9,rho_b,ye,rkm, &
       yneu,ypro,yhe4,ylight,yheavies,zbar,abar,entropy,sn_ave
  else
    write (mainout_unit,'(i8,X,100(es23.16,2X))') cnt,ctime,t9,rho_b,ye,rkm, &
       yneu,ypro,yhe4,ylight,yheavies,zbar,abar,entropy
  endif
 
end subroutine output_mainout
 
 
!>
!! Output Snapshot files
!!
!! Snapshots store the abundances and mass fractions of all nuclei
!! at certain times. An example file looks like:
!! \file{
!! time    temp    dens
!! 9.77198412010000E-02   9.99927332947922E+00   3.88380373937259E+06
!!  nin     zin       y       x
!! 1    0     5.72194135572689E-01   5.72194135572689E-01
!! 0    1     3.92194484636783E-01   3.92194484636783E-01
!! 1    1     7.83033243870819E-05   1.56606648774164E-04
!! ...}
!!
!! \b Edited:
!!           - 03.04.18, M. Jacobi
!! .
subroutine output_snapshot(ctime,T9,rho_b,Y)
  implicit none
 
  real(r_kind),intent(in) :: ctime                 !< current time
  real(r_kind),intent(in) :: t9                    !< temperature [GK]
  real(r_kind),intent(in) :: rho_b                 !< density [gcc]
  real(r_kind),dimension(net_size),intent(in):: y  !< array of abundances Y_i
 
 
  ! Output Snapshots
  call printsnap (ctime,t9,rho_b,y)
 
end subroutine output_snapshot
 
 
!>
!! Output flow files
!!
!! \b Edited:
!!           - 19.02.20, M. Jacobi
subroutine output_flow(ctime,T9,rho_b,Y)
  use flow_module,     only: flowprint, flowcalc
  implicit none
 
  real(r_kind),intent(in) :: ctime                 !< current time
  real(r_kind),intent(in) :: T9                    !< temperature [GK]
  real(r_kind),intent(in) :: rho_b                 !< density [gcc]
  real(r_kind),dimension(net_size),intent(in):: Y  !< array of abundances Y_i
 
  call flowcalc(y)
  call flowprint(ctime,t9,rho_b,y,flowcnt)
  flowcnt = flowcnt + 1
end subroutine output_flow
 
 
!>
!! Output of the Newton-Raphson loop diagnostics
!!
!! The diagnostic output contains NR iterations, timesteps,
!! the nucleus that restricts the timestep and more.
!! An example file looks like:
!! \l_file{
!! # Newton-Raphson diagnostic output
!! # 1 : global iteration count (cnt)
!! # 2 : global time [s]
!! # 3 : temperature T9 [GK]
!! # 4 : adapt stepsize loop counter (k)
!! # 5 : NR loop counter (nr_count)
!! # 6 : global timestep [s]
!! # 7 : mass conservation error
!! # 8 : maximal abundance change (eps)
!! # 9 : most rapidly evolving isotope (epsl)
!! # 10: abundance of the isotope epsl, y(epsl)
!!       0   5.20740000000E-22   7.42450000000E+00    0    1   5.20740000000E-22   9.32587340685E-15   2.22044604925E-16 cl42   2.43899441583E-08
!!       0   5.20740000000E-22   7.42450000000E+00    0    2   5.20740000000E-22   9.32587340685E-15   2.22044604925E-16  s40   2.73320077046E-08
!!       1   1.56220000000E-21   7.42450000000E+00    0    1   1.04146000000E-21   9.32587340685E-15   2.22044604925E-16  he4   3.13811866442E-03
!!       1   1.56220000000E-21   7.42450000000E+00    0    2   1.04146000000E-21   9.32587340685E-15   2.22044604925E-16 mg30   1.57736344914E-10
!! ...}
subroutine output_nr_diagnostic(cnt,k,nr_c,ctime,T9,stepsize,mtot,Y_p,Y)
  use parameter_class, only: nrdiag_every
  use global_class, only: isotope
  implicit none
  integer,intent(in)      :: cnt               !< global iteration counter
  integer,intent(in)      :: k                 !< counter in adapt_stepsize loop
  integer,intent(in)      :: nr_c              !< counter in the NR loop
  real(r_kind),intent(in) :: ctime             !< global current time
  real(r_kind),intent(in) :: t9                !< current temperature [GK]
  real(r_kind),intent(in) :: stepsize          !< current step size
  real(r_kind),intent(in) :: mtot              !< total mass
  real(r_kind),dimension(net_size) :: y_p,y    !< old/new abundances
  !
  integer      :: i
  real(r_kind) :: eps,epst
  integer      :: epsl
 
  info_entry("output_nr_diagnostic")
 
     ! check periodicity condition
     if (nrdiag_every.gt.0) then
        if (mod(cnt,nrdiag_every) .ne. 0) return
     else
        return
     endif
 
     ! which isotope has changed the most
     eps = 0.d0
     epsl= 1
     do i=1,net_size
        if(y_p(i).lt.1.d-10) cycle
        if(dabs(y_p(i)-y(i)).eq.0.d0) cycle
        epst = dabs(y(i)/y_p(i)-1.d0)
        if (epst.gt.eps) then
           eps = epst
           epsl = i
        end if
     end do
 
     ! output convergence diagnostic
     write(nrdiag_unit, &
        '(i7,2(1x,es19.11),2(1x,i4),3(1x,es19.11),a5,x,es19.11)')           &
        cnt,ctime,t9,k,nr_c,stepsize,dabs(mtot-1d0),eps,isotope(epsl)%name, &
        y(epsl)
 
 
  info_exit("output_nr_diagnostic")
 
end subroutine output_nr_diagnostic
 
 
!>
!! Print final statistics
!!
!! An example looks like:
!! \file{
!! Final time reached.
!! Final time:  2.000000000000E+17
!! Total number of iterations:       5334
!! Elapsed time [s]:  5.41572E+02
!! }
!!
!! \b Edited:
!!  - 2023-02-06, M.R. took care of wrapping of the simulation time
!!                     that could lead to negative simulation times.
!! .
subroutine output_final_stats(cnt,ctime)
  implicit none
  integer,intent(in)      :: cnt
  real(r_kind),intent(in) :: ctime
  !
  real(r_kind) :: simtime
 
  info_entry('output_final_stats')
 
  ! calculate elapsed time, check wrapping at cl_cmax
  if (cl_end .gt. cl_start) then
    simtime= float(cl_end-cl_start)/float(cl_rate)
  else
    simtime= (float(cl_end)+float(cl_cmax)-float(cl_start))/float(cl_rate)
  end if
 
  ! output final statistics
  print '(A,ES19.12)', 'Final time: ', ctime
  print '(A,I10)',     'Total number of iterations: ', cnt
  print '(A,ES12.5)',  'Elapsed time [s]: ', simtime
 
  info_exit('output_final_stats')
 
end subroutine output_final_stats
 
 
!>
!! Print current snapshot
!!
!! Snapshots store the abundances and mass fractions of all nuclei
!! at certain times. An example file looks like:
!! \file{
!! time    temp    dens
!! 9.77198412010000E-02   9.99927332947922E+00   3.88380373937259E+06
!!  nin     zin       y       x
!! 1    0     5.72194135572689E-01   5.72194135572689E-01
!! 0    1     3.92194484636783E-01   3.92194484636783E-01
!! 1    1     7.83033243870819E-05   1.56606648774164E-04
!! ...}
subroutine printsnap (t,T9,rho_b,Y)
  use global_class, only: isotope
  implicit none
 
  real(r_kind),intent(in)                :: t,t9,rho_b
  real(r_kind),dimension(net_size),intent(in) :: y
  character (len=21) :: snapfile
  integer :: snapunit, z
  real(r_kind) :: xz
 
  info_entry("printsnap")
 
  write(snapfile,'("snaps/snapsh_",i4.4,".dat")') snapcnt
  snapunit= open_outfile(snapfile)
 
  write(snapunit,'(3a8)')'time','temp','dens'
  write(snapunit,'(3es23.14)')t,t9,rho_b
  write(snapunit,'(4(a8))')'nin','zin','y','x'
  do z=1,net_size
     xz= y(z)*isotope(z)%mass
     write(snapunit,'(2(i3,2x),2(es23.14))') &
          isotope(z)%n_nr, isotope(z)%p_nr, max(1d-99,y(z)), max(1d-99,xz)
  end do
  call close_io_file(snapunit,snapfile)
  snapcnt= snapcnt + 1
  tsnp= t
  t9snp= t9
  rosnp=rho_b
 
  info_exit("printsnap")
 
end subroutine printsnap
 
 
!>
!! Calculates average neutron separation energy
!!
!! This neutron separation energy is printed in the mainout when
!! parameter_class::calc_nsep_energy was set to "yes" and
!! a valid file for parameter_class::nsep_energies_file.
subroutine calc_nseparation_energy(Y)
  use global_class, only: isotope
  implicit none
 
  real(r_kind),dimension(net_size),intent(in) :: y  !< abundances
  !
  real(r_kind) :: ysum,yasum,yzsum,abar_ions,zbar_ions
  integer :: i
 
  info_entry("calc_nseparation_energy")
 
  abar_ions = 0d0
  zbar_ions = 0d0
  do i=max(1,ihe4),net_size
     abar_ions = abar_ions + isotope(i)%mass * y(i)
     zbar_ions = zbar_ions + isotope(i)%p_nr * y(i)
     !Sn_ave = Sn_ave + Sn(i) * Y(i)
  enddo
  ysum = sum(y(max(1,ihe4):net_size))
  abar_ions = abar_ions / ysum
  zbar_ions = zbar_ions / ysum
  ysum      = 0d0
  sn_ave    = 0d0
  do i=max(1,ihe4),net_size
     if (sn(i) /= 0) then
        if ((isotope(i)%mass.ge.220).and.(isotope(i)%mass.le.260)) then
           sn_ave = sn_ave + sn(i) * y(i)
           ysum = ysum + y(i)
        end if
     endif
  enddo
  sn_ave = sn_ave / ysum
  ysum = sum(y)
  yasum = sum(y*isotope%mass)
  yzsum = sum(y*isotope%p_nr)
 
  info_exit("calc_nseparation_energy")
 
end subroutine calc_nseparation_energy
 
 
 
!>
!! @brief: Calculate the generated energy for each class of reactions separately.
!!
!! This subroutine calculates the generated energy of the nuclear reactions
!! with the help of the Q-values. This energy is written to an output only.
!! It is not fed back to the temperature and has therefore
!! no impact on the final result. Also, depending on the input the
!! subroutine can calculate the top contributors to the energy generation.
!! These top contributors can be reactions and isotopes.
!!
!! @see \ref nuclear_heating
!!
!! @author: Moritz Reichert
!! @date   29.02.20
!!
!! \b Edited:
!!  - 08.05.24: M.R. Moved the top list calculation from the nuclear heating to this subroutine.
!! .
subroutine calc_nuclear_heating(hdf5_mode, write_engen, write_toplist)
   use global_class,        only: reactionrate_type,rrate,nreac,ihe4,ineu,ipro, net_names
   use fission_rate_module, only: fissionrate_type, fissrate, nfiss
   use nucstuff_class,      only: isotope,inter_partf
   use single_zone_vars,    only: y,y_p,t9,rhob,time,stepsize,evolution_mode,ye
   use parameter_class,     only: unit,fissflag,h_engen_detailed
   use mergesort_module,    only: quicksort
   use benam_class,         only: reaction_string
   use hdf5_module
   implicit none
   ! Declare input variables
   logical, intent(in)     :: hdf5_mode    !< whether to store as hdf5 or ascii
   logical, intent(in)     :: write_engen  !< whether to write the energy generation
   logical, intent(in)     :: write_toplist!< whether to write the energy toplist
   ! Local variables
   type(reactionrate_type) :: rr_tmp     !< auxiliary shorthand for reaction
   type(fissionrate_type)  :: fr_tmp     !< auxiliary shorthand for fission reaction
   character*4  :: src                   !< reaction source
   integer      :: i,j                   !< Loop variable
   integer      :: help_parts            !< Helper variable
   real(r_kind) :: rat                   !< storage of the cached reaction rate
   ! Heating variables
   real(r_kind) :: engen_q_value         !< All reactions
   real(r_kind) :: engen_q_weak          !< weak reactions
   real(r_kind) :: engen_q_ng,engen_q_gn !< gamma-n reactions
   real(r_kind) :: engen_q_pg,engen_q_gp !< gamma-p reactions
   real(r_kind) :: engen_q_ga,engen_q_ag !< gamma-alpha reactions
   real(r_kind) :: engen_q_pa,engen_q_ap !< p-alpha reactions
   real(r_kind) :: engen_q_pn,engen_q_np !< p-n reactions
   real(r_kind) :: engen_q_na,engen_q_an !< n-alpha reactions
   real(r_kind) :: engen_fiss            !< fission reactions
   real(r_kind) :: engen_tot_ncap, engen_tot_pcap, engen_tot_acap, &
                   engen_tot_npcap,engen_tot_nacap, engen_tot_pacap !< net energy
   real(r_kind) :: engen_tot_exc         !< total energy with mass excess formula
   real(r_kind) :: entropy_src_term      !< Source term of entropy * kBT
   real(r_kind) :: mic2,gi,twopi_h2c2,mui!< Helper variables
   real(r_kind) :: kbt                   !< Helper variables
   real(r_kind),dimension(0:net_size)::pf!< partition functions
   real(r_kind) :: etaele,etapos         !< electron and positron chemical potental
   real(r_kind) :: mue                   !< electron chemical potental [MeV]
   ! Helpers
   real(r_kind) :: tmp_energy                         !< Temporary energy per reaction
   real(r_kind),dimension(net_size)  ::det_bet_engen  !< Energy of decays per parent nucleus
   real(r_kind),dimension(net_size)  ::det_ag_engen   !< Energy of (a,g)+(g,a) per parent nucleus
   real(r_kind),dimension(net_size)  ::det_pg_engen   !< Energy of (p,g)+(g,p) per parent nucleus
   real(r_kind),dimension(net_size)  ::det_ng_engen   !< Energy of (n,g)+(g,n) per parent nucleus
   real(r_kind),dimension(net_size)  ::det_pn_engen   !< Energy of (p,n)+(n,p) per parent nucleus
   real(r_kind),dimension(net_size)  ::det_ap_engen   !< Energy of (a,p)+(p,a) per parent nucleus
   real(r_kind),dimension(net_size)  ::det_an_engen   !< Energy of (a,n)+(n,a) per parent nucleus
   real(r_kind),dimension(net_size)  ::det_other_engen!< Energy of other reactions per parent nucleus
   real(r_kind),dimension(net_size)  ::det_fiss_engen !< Energy of fission per parent nucleus
   real(r_kind),dimension(net_size)  ::det_tot_engen  !< Total energy per parent nucleus
   ! Toplist stuff
   integer, dimension(nreac)         ::toplist_idx    !< Toplist indices of reactions
   real(r_kind), dimension(nreac)    ::toplist_energy !< Energy of reactions
   logical, dimension(nreac)         ::toplist_mask   !< Mask which energy has been calculated
   character*30                      ::reac_dummy     !< reaction source
   character*30, dimension(toplist_amount)::tplist_str !< Reduced list of top energies
   real(r_kind), dimension(toplist_amount)::tplist     !< Reduced list of top energies
   character*30, dimension(toplist_amount)::btlist_str !< Reduced list of top energies
   real(r_kind), dimension(toplist_amount)::btlist     !< Reduced list of bottom energies
   real(r_kind), dimension(net_size)      ::toplist_energy_iso  !< Reduced list of bottom energies
   integer, dimension(net_size)           ::toplist_idx_iso     !< Indices of isotopes
   character*5, dimension(toplist_amount) ::tplist_str_iso !< Reduced list of top energies
   real(r_kind), dimension(toplist_amount)::tplist_iso     !< Reduced list of top energies
   character*5, dimension(toplist_amount) ::btlist_str_iso !< Reduced list of top energies
   real(r_kind), dimension(toplist_amount)::btlist_iso     !< Reduced list of bottom energies
 
   info_entry("calc_nuclear_heating")
 
   ! Calculate the heating with help of the q-value
   engen_q_value = 0 ! Total energy [erg/g/s]
   engen_q_weak  = 0 ! Energy generated by weak reactions [erg/g/s]
   engen_q_ng    = 0 ! Energy generated by n-gamma reactions [erg/g/s]
   engen_q_gn    = 0 ! Energy generated by gamma-n reactions [erg/g/s]
   engen_q_pg    = 0 ! Energy generated by p-gamma reactions [erg/g/s]
   engen_q_gp    = 0 ! Energy generated by gamma-p reactions [erg/g/s]
   engen_q_ga    = 0 ! Energy generated by gamma-alpha reactions [erg/g/s]
   engen_q_ag    = 0 ! Energy generated by alpha-gamma reactions [erg/g/s]
   engen_q_np    = 0 ! Energy generated by n-p reactions [erg/g/s]
   engen_q_pn    = 0 ! Energy generated by p-n reactions [erg/g/s]
   engen_q_na    = 0 ! Energy generated by n-alpha reactions [erg/g/s]
   engen_q_an    = 0 ! Energy generated by alpha-n reactions [erg/g/s]
   engen_q_pa    = 0 ! Energy generated by p-alpha reactions [erg/g/s]
   engen_q_ap    = 0 ! Energy generated by alpha-p reactions [erg/g/s]
   engen_fiss    = 0 ! Energy generated by fission [erg/g/s]
 
   ! Initialize decay energy
   if ((h_engen_detailed) .and. (write_engen)) then
       det_bet_engen(:)   = 0.0
       det_ag_engen(:)    = 0.0
       det_pg_engen(:)    = 0.0
       det_ng_engen(:)    = 0.0
       det_pn_engen(:)    = 0.0
       det_ap_engen(:)    = 0.0
       det_an_engen(:)    = 0.0
       det_other_engen(:) = 0.0
       det_fiss_engen(:)  = 0.0
       det_tot_engen(:)   = 0.0
   end if
 
   ! Initialize the energy for the top contributors
   ! This list has one entry per reaction
   if (write_toplist) then
     toplist_energy(:) = 0d0
     toplist_energy_iso(:) = 0d0
     toplist_mask(:)   = .false.
   end if
 
   ! For NSE the reaction rates are not cached and the generated energy would be therefore weird
   ! Set it to zero for this case
   ! Loop through all reactions
   do i=1,nreac
        rr_tmp = rrate(i)
        src = rr_tmp%source
        rat = rr_tmp%cached
 
        if ((evolution_mode .eq. em_nse) .and. (.not. rr_tmp%is_weak)) cycle
 
        ! Do nothing if the rate is zero
        if (rat .eq. 0) cycle
 
        ! dYdt depends on the reaclib chapter and we therefore have to separate them
        select case (rr_tmp%group)
        case(1:3,11)! Reaclib chapter 1-3 and 11
            ! Energy per reaction
            tmp_energy = rat*y(rr_tmp%parts(1)) * rr_tmp%q_value
            ! Total energy
            engen_q_value = engen_q_value + tmp_energy
            ! Seperate the reaction types, weak reactions and photodissociations should be in chapter 1-3
 
            ! weak rates (only contained in chapter 1-3)
            if (rr_tmp%is_weak) then
            engen_q_weak = engen_q_weak + tmp_energy
            if ((h_engen_detailed) .and. (write_engen)) then
                det_bet_engen(rr_tmp%parts(1)) = det_bet_engen(rr_tmp%parts(1))+tmp_energy
            end if
            end if
            ! (gamma,n) reactions
            if(rr_tmp%reac_type.eq.rrt_gn) then
            engen_q_gn = engen_q_gn + tmp_energy
            ! Detailed output
            if ((h_engen_detailed) .and. (write_engen)) then
                if (rr_tmp%parts(2) .eq. ineu) then
                help_parts = rr_tmp%parts(3)
                else
                help_parts = rr_tmp%parts(2)
                end if
                det_ng_engen(help_parts) = det_ng_engen(help_parts)+tmp_energy
            end if
 
            endif
            ! (gamma,p) reactions
            if(rr_tmp%reac_type.eq.rrt_gp) then
            engen_q_gp = engen_q_gp + tmp_energy
            ! Detailed output
            if ((h_engen_detailed) .and. (write_engen)) then
                if (rr_tmp%parts(2) .eq. ipro) then
                help_parts = rr_tmp%parts(3)
                else
                help_parts = rr_tmp%parts(2)
                end if
                det_pg_engen(help_parts) = det_pg_engen(help_parts)+tmp_energy
            end if
            endif
            ! (gamma,alpha) reactions
            if(rr_tmp%reac_type.eq.rrt_ga) then
            engen_q_ga = engen_q_ga + tmp_energy
            ! Detailed output
            if ((h_engen_detailed) .and. (write_engen)) then
                if (rr_tmp%parts(2) .eq. ihe4) then
                help_parts = rr_tmp%parts(3)
                else
                help_parts = rr_tmp%parts(2)
                end if
                det_ag_engen(help_parts) = det_ag_engen(help_parts)+tmp_energy
            end if
            endif
 
 
        if(rr_tmp%reac_type.eq.rrt_o) then
            if ((h_engen_detailed) .and. (write_engen)) then
            det_other_engen(rr_tmp%parts(1)) = det_other_engen(rr_tmp%parts(1))+tmp_energy
            end if
        end if
 
        case(4:7)! Reaclib chapter 4-7
            ! Energy per reaction
            tmp_energy = rat*y(rr_tmp%parts(1)) * y(rr_tmp%parts(2)) &
                        * rr_tmp%q_value
            ! Total energy
            engen_q_value = engen_q_value + tmp_energy
 
            ! Seperate the reaction types
            ! (n-gamma) reactions
            if(rr_tmp%reac_type.eq.rrt_ng) then
            engen_q_ng = engen_q_ng + tmp_energy
            ! Detailed output
            if ((h_engen_detailed) .and. (write_engen)) then
                if (rr_tmp%parts(2) .eq. ineu) then
                help_parts = rr_tmp%parts(1)
                else
                help_parts = rr_tmp%parts(2)
                end if
                det_ng_engen(help_parts) = det_ng_engen(help_parts)+tmp_energy
            end if
            endif
            ! (p-gamma) reactions
            if(rr_tmp%reac_type.eq.rrt_pg) then
            engen_q_pg = engen_q_pg + tmp_energy
            ! Detailed output
            if ((h_engen_detailed) .and. (write_engen)) then
                if (rr_tmp%parts(2) .eq. ipro) then
                help_parts = rr_tmp%parts(1)
                else
                help_parts = rr_tmp%parts(2)
                end if
                det_pg_engen(help_parts) = det_pg_engen(help_parts)+tmp_energy
            end if
            endif
            ! (alpha gamma) reactions
            if(rr_tmp%reac_type.eq.rrt_ag) then
            engen_q_ag = engen_q_ag + tmp_energy
            ! Detailed output
            if ((h_engen_detailed) .and. (write_engen)) then
                if (rr_tmp%parts(2) .eq. ihe4) then
                help_parts = rr_tmp%parts(1)
                else
                help_parts = rr_tmp%parts(2)
                end if
                det_ag_engen(help_parts) = det_ag_engen(help_parts)+tmp_energy
            end if
            endif
            ! (p,alpha) reactions
            if(rr_tmp%reac_type.eq.rrt_pa) then
            engen_q_pa = engen_q_pa + tmp_energy
            ! Detailed output
            if ((h_engen_detailed) .and. (write_engen)) then
                if (rr_tmp%parts(3) .eq. ihe4) then
                help_parts = rr_tmp%parts(4)
                else
                help_parts = rr_tmp%parts(3)
                end if
                det_ap_engen(help_parts) = det_ap_engen(help_parts)+tmp_energy
            end if
            endif
            ! (alpha,p) reactions
            if(rr_tmp%reac_type.eq.rrt_ap) then
            engen_q_ap = engen_q_ap + tmp_energy
            ! Detailed output
            if ((h_engen_detailed) .and. (write_engen)) then
                if (rr_tmp%parts(2) .eq. ihe4) then
                help_parts = rr_tmp%parts(1)
                else
                help_parts = rr_tmp%parts(2)
                end if
                det_ap_engen(help_parts) = det_ap_engen(help_parts)+tmp_energy
            end if
            endif
            ! (n,p) reactions
            if(rr_tmp%reac_type.eq.rrt_np) then
            engen_q_np = engen_q_np + tmp_energy
            ! Detailed output
            if ((h_engen_detailed) .and. (write_engen)) then
                if (rr_tmp%parts(3) .eq. ipro) then
                help_parts = rr_tmp%parts(4)
                else
                help_parts = rr_tmp%parts(3)
                end if
                det_pn_engen(help_parts) = det_pn_engen(help_parts)+tmp_energy
            end if
            endif
            ! (p,n) reactions
            if(rr_tmp%reac_type.eq.rrt_pn) then
            engen_q_pn = engen_q_pn + tmp_energy
            ! Detailed output
            if ((h_engen_detailed) .and. (write_engen)) then
                if (rr_tmp%parts(1) .eq. ipro) then
                help_parts = rr_tmp%parts(2)
                else
                help_parts = rr_tmp%parts(1)
                end if
                det_pn_engen(help_parts) = det_pn_engen(help_parts)+tmp_energy
            end if
            endif
            ! (alpha,n) reactions
            if(rr_tmp%reac_type.eq.rrt_an) then
            engen_q_an = engen_q_an + tmp_energy
            ! Detailed output
            if ((h_engen_detailed) .and. (write_engen)) then
                if (rr_tmp%parts(1) .eq. ihe4) then
                help_parts = rr_tmp%parts(2)
                else
                help_parts = rr_tmp%parts(1)
                end if
                det_an_engen(help_parts) = det_an_engen(help_parts)+tmp_energy
            end if
            endif
            ! (n,alpha) reactions
            if(rr_tmp%reac_type.eq.rrt_na) then
            engen_q_na = engen_q_na + tmp_energy
            ! Detailed output
            if ((h_engen_detailed) .and. (write_engen)) then
                if (rr_tmp%parts(3) .eq. ihe4) then
                help_parts = rr_tmp%parts(4)
                else
                help_parts = rr_tmp%parts(3)
                end if
                det_an_engen(help_parts) = det_an_engen(help_parts)+tmp_energy
            end if
            endif
 
            if(rr_tmp%reac_type.eq.rrt_o) then
            if ((h_engen_detailed) .and. (write_engen)) then
                det_other_engen(rr_tmp%parts(1)) = det_other_engen(rr_tmp%parts(1))+tmp_energy
            end if
            end if
 
        case(8:9)! Reaclib chapter 8-9
            ! Energy per reaction
            tmp_energy = rat*y(rr_tmp%parts(1)) * y(rr_tmp%parts(2)) &
                            *y(rr_tmp%parts(3)) * rr_tmp%q_value
 
            ! Total energy
            engen_q_value = engen_q_value + tmp_energy
 
            if(rr_tmp%reac_type.eq.rrt_o) then
                if ((h_engen_detailed) .and. (write_engen)) then
                    det_other_engen(rr_tmp%parts(1)) = det_other_engen(rr_tmp%parts(1))+tmp_energy
                end if
            end if
 
        case(10)! Reaclib chapter 8
            ! Energy per reaction
            tmp_energy = rat*y(rr_tmp%parts(1)) * y(rr_tmp%parts(2)) &
                            *y(rr_tmp%parts(3)) * y(rr_tmp%parts(4))* rr_tmp%q_value
            ! Total energy
            engen_q_value = engen_q_value + tmp_energy
 
            if(rr_tmp%reac_type.eq.rrt_o) then
                if ((h_engen_detailed) .and. (write_engen)) then
                det_other_engen(rr_tmp%parts(1)) = det_other_engen(rr_tmp%parts(1))+tmp_energy
                end if
            end if
 
 
        case default ! This should not happen
            cycle
        end select
 
        if (write_toplist) then
            ! Store the energy of the top contributors
            ! This list has one entry per reaction and will be later sorted
            if (tmp_energy .ne. 0d0) then
                toplist_energy(i) = tmp_energy
                toplist_mask(i)   = .true.
            end if
        end if
    end do
 
    ! Do the same for the fission reactions
    if (evolution_mode .ne. em_nse) then
        if (fissflag.ne.0) then
            fission:do i=1,nfiss
                fr_tmp = fissrate(i)
                rat = fr_tmp%cached
                if (rat .eq. 0) cycle
                if ((fr_tmp%mode.eq.2).or.(fr_tmp%mode.eq.3)) then  ! spont. and b-delayed fission
                engen_q_value = engen_q_value + rat*y(fissrate(i)%fissparts(1)) * fissrate(i)%averageQ
                engen_fiss    = engen_fiss + rat*y(fissrate(i)%fissparts(1)) * fissrate(i)%averageQ
                tmp_energy    = rat*y(fissrate(i)%fissparts(1)) * fissrate(i)%averageQ
                else if (fr_tmp%mode.eq.1) then ! n-induced fission
                engen_q_value = engen_q_value + rat*y(fissrate(i)%fissparts(1))*y(fissrate(i)%fissparts(2)) * fissrate(i)%averageQ
                engen_fiss    = engen_fiss + rat*y(fissrate(i)%fissparts(1))*y(fissrate(i)%fissparts(2)) * fissrate(i)%averageQ
                tmp_energy    = rat*y(fissrate(i)%fissparts(1))*y(fissrate(i)%fissparts(2)) * fissrate(i)%averageQ
                end if
 
                if ((h_engen_detailed) .and. (write_engen)) then
                det_fiss_engen(fissrate(i)%fissparts(1)) = det_fiss_engen(fissrate(i)%fissparts(1))+tmp_energy
                end if
 
            end do fission
        end if
    end if
 
 
   engen_tot_exc = 0
   do i=1,net_size
     if (y(i).gt.1d-25) then
         tmp_energy = -(isotope(i)%mass_exc)*(y(i)-y_p(i))
         engen_tot_exc = engen_tot_exc + tmp_energy  ! energy generated by species i [MeV/baryon]
 
         if (write_toplist) then
            toplist_energy_iso(i) = tmp_energy
         end if
 
         if ((h_engen_detailed) .and. (write_engen)) then
           det_tot_engen(i) = -(isotope(i)%mass_exc)*(y(i)-y_p(i))
         end if
     endif
   end do
   engen_tot_exc = engen_tot_exc/stepsize / unit%hix
 
   entropy_src_term = 0.d0
   call inter_partf (t9, pf)
 
   ! Take care that one does not go off-grid
   if (rhob .gt. 1d-10) then
    ! Get the electron chemical potential
    call chempot(t9*1d9,rhob,ye,etaele,etapos)
    mue = etaele * unit%k_mev * t9*1d9
   else
    mue = 0.d0
   end if
   mue = mue + unit%mass_e
 
   kbt = t9*1d9 * unit%k_mev ! 8.617332d-11 * T [MeV]
   do i=1,net_size
     if (y(i).gt.1d-25) then
         mic2 = isotope(i)%mass*unit%mass_u - isotope(i)%p_nr * unit%mass_e + isotope(i)%mass_exc    ! [MeV]
         gi = 2d0*isotope(i)%spin + 1d0                               ! statistical weight
         twopi_h2c2 = 2d0*unit%pi/(unit%h_mevs*unit%clight)**2        ! 2\pi/(hc)^2 [MeV^-2 cm^-2]
         mui = -kbt*log( gi*pf(i) / (rhob*y(i)*unit%n_a) * &
             (twopi_h2c2 * mic2 * kbt)**1.5d0 )
         ! entropy generated by species i [MeV/baryon]
         entropy_src_term = entropy_src_term -(mic2 + mui + isotope(i)%p_nr*mue)*(y(i)-y_p(i))
     endif
   end do
 
   entropy_src_term = entropy_src_term/ unit%hix/stepsize
   ! Convert from MeV/Baryon/s - > Erg/g/s
   engen_q_value = engen_q_value / unit%hix
   engen_q_weak  = engen_q_weak  / unit%hix
   engen_q_ng    = engen_q_ng    / unit%hix
   engen_q_gn    = engen_q_gn    / unit%hix
   engen_q_pg    = engen_q_pg    / unit%hix
   engen_q_gp    = engen_q_gp    / unit%hix
   engen_q_ga    = engen_q_ga    / unit%hix
   engen_q_ag    = engen_q_ag    / unit%hix
   engen_q_np    = engen_q_np    / unit%hix
   engen_q_pn    = engen_q_pn    / unit%hix
   engen_q_na    = engen_q_na    / unit%hix
   engen_q_an    = engen_q_an    / unit%hix
   engen_q_pa    = engen_q_pa    / unit%hix
   engen_q_ap    = engen_q_ap    / unit%hix
   engen_fiss    = engen_fiss    / unit%hix
   ! Detailed output
   if ((h_engen_detailed) .and. (write_engen)) then
     det_tot_engen   = det_tot_engen/stepsize / unit%hix
     det_bet_engen   = det_bet_engen   / unit%hix
     det_ag_engen    = det_ag_engen    / unit%hix
     det_pg_engen    = det_pg_engen    / unit%hix
     det_ng_engen    = det_ng_engen    / unit%hix
     det_pn_engen    = det_pn_engen    / unit%hix
     det_ap_engen    = det_ap_engen    / unit%hix
     det_an_engen    = det_an_engen    / unit%hix
     det_other_engen = det_other_engen / unit%hix
     det_fiss_engen  = det_fiss_engen  / unit%hix
   end if
 
   !Calculate net energy
   engen_tot_ncap = engen_q_ng + engen_q_gn
   engen_tot_pcap = engen_q_pg + engen_q_gp
   engen_tot_acap = engen_q_ag + engen_q_ga
   engen_tot_npcap= engen_q_np + engen_q_pn
   engen_tot_nacap= engen_q_na + engen_q_an
   engen_tot_pacap= engen_q_pa + engen_q_ap
 
   if (write_toplist) then
     ! Store the energy of the top contributors
     ! This list has one entry per reaction and will be later sorted
     toplist_energy(:)     = toplist_energy(:) / unit%hix
     toplist_energy_iso(:) = toplist_energy_iso(:) /stepsize/ unit%hix
 
     ! Sort the list
     call quicksort(0,nreac,toplist_energy,toplist_idx)
 
     if (toplist_energy(1) .lt. toplist_energy(2)) then
        call raise_exception("There was a bug in the toplist. Energy is not sorted.", &
                             "calc_nuclear_heating")
     end if
 
     i=0
     j=1
     tplist_str(:) = "---------"
     tplist(:) = 0d0
     do while(i<10)
 
        if (j .gt. nreac) exit
 
        if (toplist_mask(toplist_idx(j))) then
            reac_dummy=trim(adjustl(reaction_string(rrate(toplist_idx(j)))))
            tplist_str(i+1) = reaction_string(rrate(toplist_idx(j)))
            tplist(i+1) = toplist_energy(j)
            i = i+1
        else
            continue
        end if
 
        j = j+1
     end do
 
     if (tplist(1) .lt. tplist(2)) then
        call raise_exception("There was a bug in the toplist. Energy is not sorted.", &
                             "calc_nuclear_heating")
     end if
 
     i=0
     j=1
     btlist_str(:) = "---------"
     btlist(:) = 0d0
     do while(i<10)
 
        if (j .gt. nreac) exit
 
        if (toplist_mask(toplist_idx(nreac+1-j))) then
            reac_dummy=trim(adjustl(reaction_string(rrate(toplist_idx(nreac+1-j)))))
            btlist_str(i+1) = reaction_string(rrate(toplist_idx(nreac+1-j)))
            btlist(i+1) = toplist_energy(nreac+1-j)
            i = i+1
        else
            continue
        end if
 
        j = j+1
     end do
 
 
    ! Do the same for the isotopes
    call quicksort(0,net_size,toplist_energy_iso,toplist_idx_iso)
 
    i=0
    j=1
    tplist_str_iso(:) = "-----"
    tplist_iso(:) = 0d0
    do while(i<10)
 
        if (j .gt. net_size) exit
 
        if (toplist_energy_iso(j) .ne. 0d0) then
            tplist_str_iso(i+1) = net_names(toplist_idx_iso(j))
            tplist_iso(i+1) = toplist_energy_iso(j)
            i = i+1
        else
            continue
        end if
 
        j = j+1
    end do
 
    i=0
    j=1
    btlist_str_iso(:) = "-----"
    btlist_iso(:) = 0d0
    do while(i<10)
 
        if (j .gt. net_size) exit
 
        if (toplist_energy_iso(net_size+1-j) .ne. 0d0) then
            btlist_str_iso(i+1) = net_names(toplist_idx_iso(net_size+1-j))
            btlist_iso(i+1) = toplist_energy_iso(net_size+1-j)
            i = i+1
        else
            continue
        end if
 
        j = j+1
    end do
 
 
 
     write(toplist_unit,"(es12.4, 4X, es12.4, 4X, "//int_to_str(toplist_amount)//"(A,3X), 3X,"//&
                        int_to_str(toplist_amount)//"(es12.4,3X), 4X, "//int_to_str(toplist_amount)//&
                        "(A,3X), 2X,"//int_to_str(toplist_amount)//"(es12.4,3X),'    |     ',"//&
                        int_to_str(toplist_amount)//"(A,3X), 3X,"//int_to_str(toplist_amount)//&
                        "(es12.4,3X), 4X, "//int_to_str(toplist_amount)//&
                        "(A,3X), 3X,"//int_to_str(toplist_amount)//"(es12.4,3X))" )&
                        time,engen_tot_exc, tplist_str, tplist, btlist_str,btlist, tplist_str_iso, tplist_iso, &
                        btlist_str_iso, btlist_iso
 
   end if
 
 
 
   if ((.not. hdf5_mode) .and. write_engen) then
      ! Write the output (detailed)
      !write(nuc_heat_file,'(*(es23.14,3x))') time,engen_q_value,engen_q_weak,&
      !                                       engen_q_ng,engen_q_gn,engen_tot_ncap,&
      !                                       engen_q_pg,engen_q_gp,engen_tot_pcap,&
      !                                       engen_q_ag,engen_q_ga,engen_tot_acap,&
      !                                       engen_q_np,engen_q_pn,engen_tot_npcap,&
      !                                       engen_q_na,engen_q_an,engen_tot_nacap,&
      !                                       engen_q_pa,engen_q_ap,engen_tot_pacap,&
      !                                       engen_fiss
      ! Write the output (more relaxed)
      write(nuc_heat_file,'(*(es23.14E3,3x))') time,engen_tot_exc,entropy_src_term,&
                                               engen_q_weak,&
                                               engen_tot_ncap,&
                                               engen_tot_pcap,&
                                               engen_tot_acap,&
                                               engen_tot_npcap,&
                                               engen_tot_nacap,&
                                               engen_tot_pacap,&
                                               engen_fiss
   end if
 
#ifdef USE_HDF5
   if ((write_engen) .and. (hdf5_mode)) then
      call extend_engen(time,engen_tot_exc, entropy_src_term,engen_tot_ncap, &
                        engen_tot_pcap, engen_tot_acap,engen_tot_npcap, &
                        engen_tot_nacap, engen_tot_pacap, engen_q_weak, &
                        engen_fiss, &
                        det_bet_engen, det_ag_engen, det_pg_engen, &
                        det_ng_engen, det_pn_engen, det_ap_engen, &
                        det_an_engen, det_other_engen, det_fiss_engen, &
                        det_tot_engen, &
                        h_engen_detailed)
   end if
#endif
 
 
   info_exit("calc_nuclear_heating")
 
end subroutine calc_nuclear_heating
 
 
 
 
 
 
!>
!! Calculate average timescales of different classes of reactions
!!
!! These timescales are calculated as in
!! [Arcones et al. 2012](https://ui.adsabs.harvard.edu/abs/2012ApJ...750...18A/abstract)
!! (Eq. 1-6). The (p,gamma) timescale is given by e.g.,
!! \f[
!! \tau_{p\gamma} = \left[ \frac{\rho Y_p}{Y_h}N_A \sum \langle \sigma \nu \rangle_{p,\gamma}(Z,A)Y(Z,A) \right]^{-1}
!! \f]
!!
!! An example file looks like:
!! \l_file{
!! # time[s]     Temperature [GK]     tau_ng [s]     tau_gn [s]     tau_pg [s]     tau_gp [s]     tau_np [s]     tau_pn [s]     tau_an [s]     tau_na [s]     tau_beta [s]     tau_alpha [s]     tau_nfiss [s]     tau_sfiss [s]     tau_bfiss [s]
!!  1.1130000000000000E+00    9.1821750000000009E+00    0.0000000000000000E+00    0.0000000000000000E+00    0.0000000000000000E+00    0.0000000000000000E+00    0.0000000000000000E+00    0.0000000000000000E+00    0.0000000000000000E+00    0.0000000000000000E+00    1.2642767370327558E+01    2.9487624276191804E+00    0.0000000000000000E+00    0.0000000000000000E+00    0.0000000000000000E+00
!!  1.1220000000000001E+00    8.2638489999999951E+00    0.0000000000000000E+00    0.0000000000000000E+00    0.0000000000000000E+00    0.0000000000000000E+00    0.0000000000000000E+00    0.0000000000000000E+00    0.0000000000000000E+00    0.0000000000000000E+00    1.0557999180036711E+02    2.5923655028366777E+01    0.0000000000000000E+00    0.0000000000000000E+00    0.0000000000000000E+00
!!  1.1247618491767180E+00    7.9904974735555152E+00    9.5690472538317728E-13    5.3150487750907188E-11    7.4419112944811684E-10    9.8479480155439862E-08    1.9792689873245490E-12    1.9444248037144314E-12    4.8686627517806918E-10    3.9657685448682009E-12    1.4600540103519526E+01    5.2634156341152767E+01    0.0000000000000000E+00    0.0000000000000000E+00    0.0000000000000000E+00
!! ...}
subroutine calc_av_timescales(ctime,T9,Y,hdf5)
  use global_class,        only: reactionrate_type,rrate,nreac,ihe4,ineu,ipro
  use fission_rate_module, only: fissionrate_type, fissrate, nfiss
  use file_handling_class
#ifdef USE_HDF5
  use hdf5_module,         only: extend_av_timescales
#endif
  implicit none
  real(r_kind),intent(in)       :: ctime !< current time [s]
  real(r_kind),intent(in)       :: t9    !< Temperature [GK]
  logical,optional,intent(in)   :: hdf5  !< Whether the output is written to hdf5 or not
  real(r_kind),dimension(net_size),intent(in)   :: y  !< abundances
 
  ! average timescales
  real(r_kind) :: tau_ng           !< (n,gamma)
  real(r_kind) :: tau_gn           !< (gamma,n)
  real(r_kind) :: tau_pg           !< (p,gamma)
  real(r_kind) :: tau_gp           !< (gamma,p)
  real(r_kind) :: tau_np           !< (n,p)
  real(r_kind) :: tau_pn           !< (p,n)
  real(r_kind) :: tau_an           !< (alpha,n)
  real(r_kind) :: tau_na           !< (n,alpha)
  real(r_kind) :: tau_ag           !< (alpha,gamma)
  real(r_kind) :: tau_ga           !< (gamma,alpha)
  real(r_kind) :: tau_pa           !< (p,alpha)
  real(r_kind) :: tau_ap           !< (alpha,p)
  !real(r_kind) :: tau_npro        !< neutron production timescale
  !real(r_kind) :: tau_fdis        !< photodissociation timescale
  real(r_kind) :: tau_beta         !< beta-decay timescale
  real(r_kind) :: tau_alpha        !< alpha-decay timescale
  real(r_kind) :: tau_nfiss        !< n-induced fission timescale
  real(r_kind) :: tau_sfiss        !< spontaneous fission timescale
  real(r_kind) :: tau_bfiss        !< beta-delayed fission timescale
  integer      :: nuc              !< auxiliary var
  real(r_kind) :: rat              !< Value of reaction rate
  real(r_kind) :: ng_sum,gn_sum,pg_sum,gp_sum,beta_sum
  real(r_kind) :: np_sum,pn_sum,an_sum,na_sum,ag_sum,ga_sum
  real(r_kind) :: ap_sum,pa_sum
  real(r_kind) :: alpha_sum, nfiss_sum, sfiss_sum, bfiss_sum
  real(r_kind) :: ysum               !< Sum of abundances heavier than helium
  integer      :: i                  !< Loop variable
  type(reactionrate_type) :: rr_tmp  !< auxiliary shorthand for reaction
  type(fissionrate_type)  :: fr_tmp  !< auxiliary shorthand for reaction
  character*4  :: src                !< reaction source
  real(r_kind) :: q_sm = 1.e-20      !< small quantity to add to denominator
 
  info_entry("calc_av_timescales")
 
 
  ! to avoid division by zero
  ng_sum    = q_sm
  gn_sum    = q_sm
  pg_sum    = q_sm
  gp_sum    = q_sm
  np_sum    = q_sm
  pn_sum    = q_sm
  an_sum    = q_sm
  na_sum    = q_sm
  ag_sum    = q_sm
  ga_sum    = q_sm
  ap_sum    = q_sm
  pa_sum    = q_sm
  beta_sum  = q_sm
  alpha_sum = q_sm
  nfiss_sum = q_sm
  sfiss_sum = q_sm
  bfiss_sum = q_sm
 
 
  ysum = sum(y(max(1,ihe4):net_size))
 
  do i=1,nreac
     rr_tmp = rrate(i)
     src = rr_tmp%source
     rat = rr_tmp%cached
 
     ! beta-decays (not including beta-delayed fission)
     if ((rr_tmp%group.eq.1).or.(rr_tmp%group.eq.2).or.(rr_tmp%group.eq.3) .or.(rr_tmp%group.eq.11)) then
        if ((rr_tmp%reac_type .eq.rrt_betm).or.(rr_tmp%reac_type .eq.rrt_betp).and.(rr_tmp%source.ne.'ms99')) then
        ! only include significant reactions
           nuc = rr_tmp%parts(1)
           if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20) then
              beta_sum = beta_sum + rat*y(nuc)
           end if
        endif
     endif
 
     ! alpha-decays
     if ((rr_tmp%reac_type.eq.rrt_alpd)) then
        nuc = rr_tmp%parts(1)
        if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20) then
           alpha_sum = alpha_sum + rat*y(nuc)
        end if
     end if
 
     ! (a,p)
     if(rr_tmp%reac_type.eq.rrt_ap) then
        nuc = rr_tmp%parts(2)
        if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20) then
           ap_sum  = ap_sum  + rat*y(nuc)*y(ihe4)
        end if
     endif
 
     ! (p,a)
     if(rr_tmp%reac_type.eq.rrt_pa) then
        nuc = rr_tmp%parts(2)
        if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20) then
           pa_sum  = pa_sum  + rat*y(nuc)*y(ipro)
        end if
     endif
 
     ! (n,gamma)
     if(rr_tmp%reac_type.eq.rrt_ng) then
        nuc = rr_tmp%parts(2)
        if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20) then
           ng_sum  = ng_sum  + rat*y(nuc)*y(ineu)
        end if
     endif
 
     ! (gamma,n)
     if(rr_tmp%reac_type.eq.rrt_gn) then
        nuc = rr_tmp%parts(1)
        if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20) then
           gn_sum  = gn_sum  + rat*y(nuc)
        end if
     endif
 
     ! (p,gamma)
     if(rr_tmp%reac_type.eq.rrt_pg) then
        nuc = rr_tmp%parts(2)
        if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20) then
           pg_sum  = pg_sum  + rat*y(nuc)*y(ipro)
        end if
     endif
 
     ! (gamma,p)
     if(rr_tmp%reac_type.eq.rrt_gp) then
        nuc = rr_tmp%parts(1)
        if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20) then
           gp_sum  = gp_sum  + rat*y(nuc)
        end if
     endif
 
     ! (a,g)
     if(rr_tmp%reac_type.eq.rrt_ag) then
        nuc = rr_tmp%parts(2)
        if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20) then
           ag_sum  = ag_sum  + rat*y(nuc)*y(ihe4)
        end if
     endif
 
     ! (g,a)
     if(rr_tmp%reac_type.eq.rrt_ga) then
        nuc = rr_tmp%parts(1)
        if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20) then
           ga_sum  = ga_sum  + rat*y(nuc)
        end if
     endif
 
     ! (n,p)
     if(rr_tmp%reac_type.eq.rrt_np) then
        nuc = rr_tmp%parts(2)
        if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20) then
           np_sum  = np_sum  + rat*y(nuc)*y(ineu)
        end if
     endif
 
     ! (p,n)
     if(rr_tmp%reac_type.eq.rrt_pn) then
        nuc = rr_tmp%parts(2)
        if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20) then
           pn_sum  = pn_sum  + rat*y(nuc)*y(ipro)
        end if
     endif
 
     ! (a,n)
     if(rr_tmp%reac_type.eq.rrt_an) then
        nuc = rr_tmp%parts(2)
        if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20) then
           an_sum  = an_sum  + rat*y(nuc)*y(ihe4)
        end if
     endif
 
     ! (n,a)
     if(rr_tmp%reac_type.eq.rrt_na) then
        nuc = rr_tmp%parts(2)
        if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20) then
           na_sum  = na_sum  + rat*y(nuc)*y(ineu)
        end if
     endif
 
  enddo
 
  fission:do i=1,nfiss
     fr_tmp = fissrate(i)
     rat = fr_tmp%cached
     ! neutron-induced fission
     if(fr_tmp%reac_type.eq.rrt_nf) then
        nuc = fr_tmp%fissparts(2)
        if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20) then
           nfiss_sum = nfiss_sum + rat*y(nuc)*y(ineu)
        end if
     end if
 
     ! spontaneous fission
     if(fr_tmp%reac_type.eq.rrt_sf) then
        nuc = fr_tmp%fissparts(1)
        if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20) then
           sfiss_sum = sfiss_sum + rat*y(nuc)
        end if
     end if
 
     ! beta-delayed fission
     if(fr_tmp%reac_type.eq.rrt_bf) then
        nuc = fr_tmp%fissparts(1)
        if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20) then
           bfiss_sum = bfiss_sum + rat*y(nuc)
        end if
     end if
 
  end do fission
 
  tau_ng    = ysum / ng_sum
  tau_gn    = ysum / gn_sum
  tau_pg    = ysum / pg_sum
  tau_gp    = ysum / gp_sum
  tau_ag    = ysum / ag_sum
  tau_ga    = ysum / ga_sum
  tau_np    = ysum / np_sum
  tau_pn    = ysum / pn_sum
  tau_an    = ysum / an_sum
  tau_na    = ysum / na_sum
  tau_pa    = ysum / pa_sum
  tau_ap    = ysum / ap_sum
  tau_beta  = ysum / beta_sum
  tau_alpha = ysum / alpha_sum
  tau_nfiss = ysum / nfiss_sum
  tau_sfiss = ysum / sfiss_sum
  tau_bfiss = ysum / bfiss_sum
 
  if (.not. present(hdf5)) then
     ! write it to an ascii file
     write(tsfile, '(*(es23.16,3x))') &
        ctime, t9, tau_ga,tau_ag, tau_ng, tau_gn, tau_pg, tau_gp, tau_np, &
                   tau_pn, tau_an, tau_na,tau_pa,tau_ap, tau_beta, tau_alpha, &
                   tau_nfiss, tau_sfiss, tau_bfiss
  else
#ifdef USE_HDF5
     ! Write it to the hdf5
     call extend_av_timescales(ctime,tau_ga, tau_ag, tau_ng, tau_gn, tau_pg,   &
                                     tau_gp, tau_np, tau_pn, tau_an, tau_na,   &
                                     tau_ap, tau_pa, tau_beta, tau_alpha,  &
                                     tau_nfiss, tau_sfiss, tau_bfiss)
 
#else
     continue
#endif
  endif
 
  info_exit('calc_av_timescales')
 
end subroutine calc_av_timescales
 
 
!>
!! Compute final abundances and write final output
!!
!! There are three files written _finabsum.dat_
!! containing abundances summed over equal A,
!! _finab.dat_ containing all final abundances,
!! _finabelem.dat_ containing all abundances
!! summed over equal Z.
subroutine finab(Y)
  use global_class, only: isotope, net_size
  use file_handling_class
  implicit none
 
  real(r_kind),dimension(net_size),intent(in)  :: y
  integer,dimension(net_size)                  :: list
  real(r_kind),dimension(:),allocatable        :: xa,yz
  integer            :: finfile,finsfile,finzfile
  integer            :: mval,zmval
  integer            :: i,k
  integer            :: mass,protonnum
 
  info_entry('finab')
 
  finfile= open_outfile('finab.dat')
  write(finfile,'(A)') "# 1:A 2:Z 3:N 4:Yi 5:Xi"
 
  call finab_sort(net_size,list)
 
  do k=1,net_size
     i=list(k)
     if (y(i).le.1.d-20) cycle
     write(finfile,'(3(i5,2x),2es23.14)') isotope(i)%mass,   &
          isotope(i)%p_nr,isotope(i)%n_nr,y(i),              &
          y(i)*isotope(i)%mass
  end do
  call close_io_file(finfile,'finab.dat')
 
  mval = maxval(isotope%mass)
  allocate(xa(mval))
  zmval = maxval(isotope%p_nr)
  allocate(yz(zmval))
 
  xa = 0.d0
  yz = 0.d0
 
  do i=1,net_size
     mass=isotope(i)%mass
     protonnum = isotope(i)%p_nr
     if(y(i).gt.1.d-20) then
        xa(mass) = xa(mass) + y(i)*mass
        if (protonnum.gt.0) yz(protonnum) = yz(protonnum) + y(i) ! omit neutrons
     end if
  end do
 
  finsfile= open_outfile('finabsum.dat')
  write(finsfile,'(A)') "# 1:A 2:Y(A) 3:X(A)"
  do i=1,mval
     !i:Mass, XA(i)/i: Abundance, XA(i): Mass fraction
     write(finsfile,'(i5,2x,es23.14,2x,es23.14)')i, xa(i)/i ,xa(i)
  end do
  finzfile= open_outfile('finabelem.dat')
  write(finzfile,'(A)') "# 1:Z 2:Y(Z)"
  do i=1,zmval
     write(finzfile,'(i5,2x,es23.14)')i, yz(i)
  end do
 
  call close_io_file(finsfile,'finabsum.dat')
 
  info_exit('finab')
end subroutine finab
 
 
!>
!! auxiliary sorting subroutine
!!
subroutine finab_sort(nsize,list)
  use global_class, only:isotope
  use file_handling_class
 
  implicit none
 
  type :: sorttype
     integer :: a
     integer :: z
     logical :: chk
  end type sorttype
 
  integer,intent(in)                        :: nsize
  integer,dimension(nsize),intent(out)      :: list
  type(sorttype),dimension(nsize)           :: tmp
  integer                                   :: i,j
  integer                                   :: amin,zmin
  integer                                   :: nxt
 
  info_entry('finab_sort')
 
  amin = 10000
  zmin = 10000
 
  do i=1,nsize
     tmp(i)%a = isotope(i)%mass
     tmp(i)%z = isotope(i)%p_nr
     tmp(i)%chk = .false.
  end do
 
  do i=1,nsize
     do j=1,nsize
        if(tmp(j)%chk) cycle
        select case (tmp(j)%a-amin)
        case(:-1)
           nxt = j
           amin = tmp(j)%a
           zmin = tmp(j)%z
        case(0)
           if (tmp(j)%z.lt.zmin) then
              nxt = j
              amin = tmp(j)%a
              zmin = tmp(j)%z
           end if
        case default
           cycle
        end select
     end do
     list(i) = nxt
     tmp(nxt)%chk = .true.
     amin=10000
     zmin=10000
  end do
 
  info_exit('finab_sort')
 
  return
 
end subroutine finab_sort
 
 
!>
!! Close files at the end of the calculation
!!
!! @author: Moritz Reichert
!!
!! \b Edited:
!!     - 01.03.17
!!     - 25.01.21
!! .
subroutine analysis_finalize()
   use parameter_class, only: track_nuclei_every, top_engen_every, &
                              timescales_every, nrdiag_every
   implicit none
 
   info_entry('analysis_finalize')
 
   ! Close the timescale file
   if (timescales_every .gt. 0) then
      call close_io_file(tsfile,'timescales.dat')
   end if
 
   ! Close the Newton-Raphson diagnostic file
   if (nrdiag_every.gt.0) then
      call close_io_file(nrdiag_unit,"nrdiag.dat")
   end if
 
   ! Close track nuclei file
   if (track_nuclei_every .gt. 0) then
      call close_io_file(track_unit,"tracked_nuclei.dat")
   end if
 
   ! Cleanup the toplist
   if (top_engen_every .gt. 0) then
    call close_io_file(toplist_unit,"toplist.dat")
  end if
 
   info_exit('analysis_finalize')
 
end subroutine analysis_finalize
 
 
end module analysis