nuclear_heating.f90

Go to the documentation of this file.

!> @file nuclear_heating.f90
!!
!! The error file code for this file is ***W31***.
!!
!! @brief Module \ref nuclear_heating for self-consistent temperature updates
!!
 
!> Takes care of the temperature and entropy updates
!! due to the nuclear energy release
#include "macros.h"
module nuclear_heating
  use ls_timmes_eos_module
  use error_msg_class, only: raise_exception, num_to_str, int_to_str
  implicit none
 
  real(r_kind), private :: engen_nuc     !< Total energy, generated by nuclear reactions
  integer     , private :: engen_unit    !< Engen output file
  integer     , private :: debug_temp    !< Debug file for temperature integration
  integer     , private :: debug_heatfrac!< Debug file for the heating fraction
 
  real(r_kind)                         , public, save :: qdot_bet = 0 !< Total energy radiated away by
  real(r_kind)                         , public, save :: qdot_nu  = 0 !< Total energy added to the system by neutrino
  real(r_kind)                         , public, save :: qdot_th  = 0 !< Total energy lost by thermal neutrinos
 
  ! Debug variables
  real(r_kind)                         , private :: debug_ffn_qd=0  !< debug variable for the qdot caused by ffn rates
  real(r_kind)                         , private :: debug_nu_qd=0   !< debug variable for the qdot caused by nu rates
  real(r_kind)                         , private :: debug_bet_qd=0  !< debug variable for the qdot caused by beta decay rates
  !
  ! Public and private fields and methods of the module.
  !
  public:: &
      output_nuclear_heating, nuclear_heating_update, nuclear_heating_init
  private:: &
      nuclear_heating_energy, nuclear_heating_entropy,&
      start_heating
 
contains
 
!>
!! Nuclear heating initialization routine
!!
!! This subroutine creates some files (e.g., debug_engen.dat_) and
!! prepares the header for the _engen.dat_ file
!!
!! \b Edited:
!!    - MR: 20.01.21 - Gave toplist file a header
!!    - MR: 07.05.24 - Removed toplist and moved it to analysis module
!! .
subroutine nuclear_heating_init (T9, rho, Ye, Y, entropy)
  use global_class,    only: net_size,isotope,nreac,rrate,heating_switch
  use parameter_class, only: top_engen_every, engen, heating_mode, &
                             heating_density
  use file_handling_class
  use single_zone_vars, only: t9h, t9h_p
  implicit none
 
  real(r_kind),intent(in)                     :: t9       !< temperature [GK]
  real(r_kind),intent(in)                     :: rho      !< density [g/cm^3]
  real(r_kind),intent(in)                     :: ye       !< electron fraction
  real(r_kind),dimension(net_size),intent(in) :: y        !< abundances
  real(r_kind),intent(out)                    :: entropy  !< [kB/baryon]
 
  info_entry("nuclear_heating_init")
 
  ! Check if heating mode is implemented
  if ((heating_mode .gt. 3) .or. (heating_mode .lt. 0)) then
    call raise_exception("Heating_mode not implemented, got "//int_to_str(heating_mode)//".",&
                         'nuclear_heating_init',310004)
  end if
 
  ! Get the initial entropy
  call start_heating (t9, rho, ye, y, entropy)
 
  engen_nuc= 0d0
  engen = 0d0
 
  if (heating_mode .eq. 1) then
    if (verbose_level .ge. 2) then
        engen_unit= open_outfile("debug_engen.dat")
        write(engen_unit,"(A)") "# 1:iteration 2:time[s] 3:timestep[s] 4:E_gen(tot) [erg/g] 5:E_gen(t) [erg/(g*s)] 6:E_gen(t) [MeV/(baryon*s)]"
    end if
  elseif (heating_mode .eq. 2) then
      t9h   = t9
      t9h_p = t9
  end if
 
 
  if (verbose_level .ge. 2) then
    call output_debug_nufrac(rrate,nreac)
    ! Initialize debug variables
    debug_ffn_qd = 0; debug_nu_qd = 0; debug_bet_qd = 0
    debug_heatfrac= open_outfile("debug_heatfrac_split.dat")
    write(debug_heatfrac,"(A)") "# 1:time[s] 2:qdot 3: qdot_ffn 4: qdot_nu 5: qdot_bet 6: qdot_ffn/total 7: qdot_nu/total 8: qdot_bet/total"
  end if
 
 
  if (verbose_level .ge. 2) then
    debug_temp = open_outfile("debug_heating_temp.dat")
    write(debug_temp,"(A)") "# 1:Timestep[s] 2:NR_cnt 3:T9_new 4:T9_old 5:T9_taken"
  end if
 
  ! Check if heating is switched on from the beginning.
  heating_switch = .false.
  if (heating_density .eq. -1) then
    heating_switch = .true.
  else if (heating_density .gt. rho) then
    heating_switch = .true.
  elseif (heating_density .le. rho) then
    heating_switch = .false.
  end if
 
  info_exit("nuclear_heating_init")
end subroutine nuclear_heating_init
 
 
!> Initialize the entrop and temperature for the heating
subroutine start_heating (T9, rho, Ye, Y, entropy)
  use global_class, only: net_size, isotope
  real(r_kind),intent(in)                     :: t9       !< temperature [GK]
  real(r_kind),intent(in)                     :: rho      !< density [g/cm^3]
  real(r_kind),intent(in)                     :: ye       !< electron fraction
  real(r_kind),dimension(net_size),intent(in) :: y        !< abundances
  real(r_kind),intent(out)                    :: entropy  !< [kB/baryon]
  !
  type(timmes_eos_state)   :: state    ! Timmes EOS state
  integer :: eos_status
  real(r_kind) :: ysum, yasum
 
  ! Get the initial entropy
  ysum= sum(y(1:net_size))
  yasum= sum(y(1:net_size)*isotope(1:net_size)%mass)
  state%abar= yasum / ysum
  call timmes_eos(ink,t9*1.d9,rho,ye,state,eos_status)
 
  if (eos_status .ne. 0 ) then
    call raise_exception("The EOS did not converge, T="//num_to_str(t9)//"GK, "//&
                         "rho="//num_to_str(rho)//"g/ccm, Ye="//num_to_str(ye),&
                         "nuclear_heating_init",310003)
  end if
 
  entropy = state%s
 
end subroutine start_heating
 
 
 
!> Reset the qdot variable.
!!
!! This is called by the Jacobian init whenever recalculating the
!! jacobian.
!!
!! @author M. Reichert
subroutine reset_qdot(time)
    use parameter_class, only: heating_mode, nu_loss_every, h_nu_loss_every
    implicit none
    real(r_kind),intent(in)  :: time  !< Current time
    real(r_kind)             :: denom !< helper variable
    real(r_kind)             :: qdot  !< helper variable
 
    if ((heating_mode .gt. 0) .or. (nu_loss_every .gt. 0) .or. &
        (h_nu_loss_every .gt. 0)) then
        ! Output status for debugging, qdot splitted in contributions
        if (verbose_level .gt. 2) then
            qdot = qdot_bet + qdot_nu
            ! Calculate denominator, ensure not to divide by zero
            if (qdot .eq. 0) then
                denom = 1
            else
                denom = abs(qdot)
            end if
 
            ! Output
            write(debug_heatfrac, '(es16.5,1x,es16.5,1x,es16.5,1x,es16.5,1x,es16.5,1x,f5.2,1x,f5.2,1x,f5.2)') &
                                time, qdot, debug_ffn_qd, debug_nu_qd, debug_bet_qd, &
                                debug_ffn_qd/denom, debug_nu_qd/denom, debug_bet_qd/denom
 
            ! Reset variables
            debug_ffn_qd = 0; debug_nu_qd = 0; debug_bet_qd = 0
        end if
 
        qdot_bet = 0
        qdot_nu  = 0
    end if
end subroutine reset_qdot
 
 
 
!< Calculate the energy radiated away by neutrinos or added to the system
!!
!! This subroutine calculates the energy radiated away by neutrinos or added to the system
!! and saves it in the variable qdot.
!! The fraction of energy is either set by a user defined input file (\ref neutrino_loss_file),
!! by the beta decay file (\ref beta_decay_file),
!! by the ffn_type rates (\ref weak_rates_file),
!! or by a user defined variable (\ref heating_frac).
!!
!! @author M. Reichert
!! @date 02.04.2023
subroutine calculate_qdot(rrate,Y,h)
    use global_class,    only: net_size, reactionrate_type
    use parameter_class, only: heating_mode, nu_loss_every, h_nu_loss_every
    implicit none
    type(reactionrate_type),intent(in)          :: rrate !< Reaction rate
    real(r_kind),dimension(net_size),intent(in) :: Y     !< Abundances
    real(r_kind),intent(in)                     :: h     !< timestep
    !-- Internal variables
    integer      :: group        !< Reaclib chapter
    real(r_kind) :: qdot_tmp     !< Temporary energy released by beta-decay
 
    if ((rrate%is_weak) .and. (heating_mode .gt. 0) .or. (nu_loss_every .gt. 0) .or. &
        (h_nu_loss_every .gt. 0)) then
        ! Get reaclib group
        group = rrate%group
        select case(group)
            case(1:3,11)
                qdot_tmp = rrate%cached * rrate%q_value * rrate%nu_frac * &
                           y(rrate%parts(1)) * h
            case(4:7)
                qdot_tmp = rrate%cached * rrate%q_value * rrate%nu_frac * &
                           y(rrate%parts(1))*y(rrate%parts(2)) * h
            case(8:9)
                qdot_tmp = rrate%cached * rrate%q_value * rrate%nu_frac * &
                           y(rrate%parts(1))*y(rrate%parts(2))*y(rrate%parts(3)) * h
            case(10)
                qdot_tmp = rrate%cached * rrate%q_value * rrate%nu_frac * &
                           y(rrate%parts(1))*y(rrate%parts(2))*y(rrate%parts(3))*y(rrate%parts(4)) * h
        end select
 
        if (verbose_level .ge. 2) then
            ! Split the energy dependent on its reaction type
            if (rrate%reac_src .eq. rrs_twr) then
                debug_ffn_qd = debug_ffn_qd + qdot_tmp
            elseif (rrate%reac_src .eq. rrs_nu) then
                debug_nu_qd = debug_nu_qd + qdot_tmp
            else
                debug_bet_qd = debug_bet_qd + qdot_tmp
            end if
        end if
 
        if (rrate%reac_src .eq. rrs_nu) then
            qdot_nu = qdot_nu + qdot_tmp
        else
            qdot_bet = qdot_bet + qdot_tmp
        end if
    end if
end subroutine calculate_qdot
 
 
!> Output the heating fractions to a debug file
!!
!! This subroutine outputs the heating fractions to a debug file. This is
!! useful for debugging purposes. The file is called debug_heating_nufrac.dat
!! and is located in the output directory. An example could look like
!! \file{
!! Group     Isotope      Source     Nu_frac     Q_value
!!     1           n        wc12     0.40000     0.78230
!!     1           t        wc12     0.40000     0.01860
!!     1         he3          ec     0.40000    -0.01900
!!     1         he6        wc12     0.55394     3.50510
!!     1         li9        wc12     0.45403    13.60600
!!...}
!!
!! @author M. Reichert
!! @date   30.03.2023
subroutine output_debug_nufrac(rrate_array,length)
    use global_class,        only: net_size, reactionrate_type, net_names
    use file_handling_class, only: open_outfile
    implicit none
    integer,intent(in)                                     :: length      !< Amount of reactions
    type(reactionrate_type),dimension(length),intent(in)   :: rrate_array !< Reaction rate array
    !-- Internal variables
    integer :: i
    integer :: file_id
 
    ! Open the debug file
    file_id = open_outfile("debug_heating_nufrac.dat")
 
    ! Write the header of the file, align the columns in the correct format
    write(file_id,"(A5,2X,A10,2X,A10,2X,A10,2X,A10)") "Group","Isotope","Source","Nu_frac","Q_value"
 
    ! Write the data, note that the ffn rates are variable and the correct value will not be
    ! reflected in the debug file. Same for the neutrino rates.
    do i=1,length
        if (rrate_array(i)%is_weak) then
            write(file_id,"(2X,I3,2X,A10,2X,A10,2X,F10.5,2X,F10.5)") rrate_array(i)%group,&
                  net_names(rrate_array(i)%parts(1)),rrate_array(i)%source,&
                  rrate_array(i)%nu_frac,rrate_array(i)%q_value
        end if
    end do
    close(file_id)
 
end subroutine output_debug_nufrac
 
 
 
 
!>
!! Modifies the temperature and entropy to account for nuclear heating
!!
!! The temperature/entropy is solved (explicitely) in the same Newton-Raphson
!! scheme as used for the abundances in e.g.,
!! \ref timestep_module::advance_implicit_euler and
!! \ref timestep_module::advance_gear.
!!
!! There are 3 different heating modes:
!! <table>
!! <caption id="multi_row">Heating modes</caption>
!! <tr><th> Value     <th> Evolved variable       <th> Additional source term
!! <tr><td> 0         <td> No heating             <td> -
!! <tr><td> 1         <td> Entropy                <td> adiabatic
!! <tr><td> 2         <td> Temperature            <td> trajectory
!! <tr><td> 3         <td> Temperature            <td> adiabatic
!! </table>
!!
!! @author M. Reichert
!! @date 10.02.23
!! .
subroutine nuclear_heating_update(nr_count,rho,Ye,pf,Y_p,Y,entropy_p,entropy,T9_p,T9,T9_tr_p,T9_tr,timestep,T9conv)
  use global_class,    only: net_size,isotope
  use nucstuff_class,  only: inter_partf
  use parameter_class, only: engen, top_engen_every, unit, heating_frac, &
                             heating_mode, freeze_rate_temp, use_thermal_nu_loss
  use thermal_neutrino_module, only: thermal_neutrinos
  implicit none
 
  integer,intent(in)                               :: nr_count     !< NR-iteration count
  real(r_kind),intent(in)                          :: rho          !< density
  real(r_kind),intent(in)                          :: ye           !< electron fraction
  real(r_kind),dimension(0:net_size),intent(out)   :: pf           !< partition functions
  real(r_kind),dimension(net_size),intent(in)      :: y_p          !< abundances
  real(r_kind),dimension(net_size),intent(in)      :: y            !< abundances at the next step
  real(r_kind),intent(in)                          :: entropy_p    !< entropy
  real(r_kind),intent(out)                         :: entropy      !< entropy for the next step
  real(r_kind),intent(in)                          :: t9_tr_p      !< temperature of the trajectory/analytic (prev.)
  real(r_kind),intent(in)                          :: t9_tr        !< temperature of the trajectory/analytic
  real(r_kind),intent(in)                          :: t9_p         !< temperature at the previous step
  real(r_kind),intent(inout)                       :: t9           !< temperature at the next step
  real(r_kind),intent(in)                          :: timestep     !< timestep (s)
  real(r_kind),intent(out)                         :: t9conv       !< Convergence of the temperature
 
  real(r_kind) :: mic2, mui, kbt     !< auxiliary variables and constants
  type(timmes_eos_state) :: state    !< Timmes EOS state
  real(r_kind) :: gi                 !< statistical weight
  real(r_kind) :: twopi_h2c2         !< 2\pi/(hc)^2 [MeV^-2 cm^-2]
  real(r_kind) :: entropy_inc        !< entropy increment during single timestep due to nuclear heating
  real(r_kind) :: engen_i            !< generated energy (per baryon) due to abundance changes of species i
  real(r_kind) :: t9_old             !< last temperature difference
  real(r_kind) :: new_tdiff          !< Temperature difference
  real(r_kind) :: abar               !< average mass number
 
  real(r_kind) :: ysum,yasum, temp_next
  real(r_kind) :: t9_calc
  integer :: i ! Loop variables
  integer :: eos_status
 
  info_entry("nuclear_heating_update")
 
  ! Save the old temperature
  t9_old = t9
 
  ! Only heat for T9 > freeze_rate_temp
  if (t9 .gt. freeze_rate_temp*1.000001)then
 
    if (use_thermal_nu_loss) then
        ysum= sum(y(1:net_size))
        yasum= sum(y(1:net_size)*isotope(1:net_size)%mass)
        abar= yasum / ysum
        call thermal_neutrinos(abar, ye, t9, rho, timestep, qdot_th)
        qdot_th = qdot_th * unit%hix ! erg/g -> MeV/baryon
    else
        qdot_th = 0d0
    end if
 
    ! Different heating treatments
    if (heating_mode .eq. 1) then
        call nuclear_heating_entropy(rho,ye,pf,y_p,y,entropy_p,entropy,t9,eos_status)
    elseif ((heating_mode .eq. 2) .or. (heating_mode .eq. 3)) then
        call nuclear_heating_energy(rho,ye,pf,y_p,y,entropy_p,entropy,t9_p,t9,t9_tr_p,t9_tr,eos_status)
    else
        call raise_exception("Heating_mode not implemented, got "//int_to_str(heating_mode)//".",&
                             'nuclear_heating_update',310004)
    end if
 
    t9_calc = t9
    ! Avoid a crash here and do not care about the temperature
    ! if the EOS fails
    if (eos_status .ne. 0) then
        t9 = t9_p
    else
        ! Check if there is a convergence problem
        ! Sometimes there are alternating temperatures
        ! obtained.
        ! Therefore, if not converged after the first 3 iterations
        ! limit the shooting!
        new_tdiff = t9-t9_old
 
        if (nr_count .ge. 3) then
            t9 = t9_old + new_tdiff*max((1d0/(float(nr_count)-1d0)),0.1d0)
        end if
 
    end if
 
    ! Don't allow T9 to go below freeze_rate_temp and
    ! avoid crazy unreasonable numbers
    t9 = max(t9,freeze_rate_temp,t9_p/2d0)
 
    ! Update the entropy
    ysum= sum(y(1:net_size))
    yasum= sum(y(1:net_size)*isotope(1:net_size)%mass)
    state%abar= yasum / ysum
    call timmes_eos(ink,t9*1d9,rho,ye,state,eos_status)
 
    entropy = state%S
 
  else
    t9 = t9
    t9_calc = t9
  end if
 
  if (verbose_level .ge. 2) then
    write(debug_temp,"(es16.5,1x,I5,1x,es16.5,1x,es16.5,1x,es16.5)") &
                       timestep, nr_count, t9_calc, t9_old, t9
  end if
 
  t9conv = dabs(max(t9_calc,t9_old)/min(t9_calc,t9_old) - 1d0)
 
  info_exit("nuclear_heating_update")
 
end subroutine nuclear_heating_update
 
 
!>
!! Modifies the temperature to account for nuclear heating
!!
!! Here, an additional equation is calculated for the temperature.
!! The temperature is calculated by:
!! \f[ \dot{T}_\mathrm{nuc} = \frac{1}{c_v} \sum_i M_i dY_i   \f]
!! with the symbols defined as:
!! <pre style="line-height:.5">
!!  - M_i     : mass excess [MeV] (see \ref global_class::isotope\%mass_exc)
!!  - cv      : Specific heat capacity at constant volume
!! .
!! </pre>
!!
!! For heating_mode 2, the temperature is calculated by:
!! \f[ \dot{T}_\mathrm{nuc} = \dot{T}_\mathrm{nuc} + \dot{T}_\mathrm{tr} \f]
!! where \f$\dot{T}_\mathrm{tr}\f$ is the temperature change of the external conditions,
!! i.e., the trajectory or the analytic equation.
!! For heating_mode 3, the temperature is calculated by:
!! \f[ \dot{T}_\mathrm{nuc} = \dot{T}_\mathrm{nuc} + \dot{T}_\mathrm{ad} \f]
!! where \f$\dot{T}_\mathrm{ad}\f$ is the temperature change assuming adiabatic conditions.
!!
!! @author M. Reichert
!! @date 10.02.23
!!
!! @see [Lippuner & Roberts 2017](https://ui.adsabs.harvard.edu/abs/2017ApJS..233...18L/abstract),
!!      [XNet](https://github.com/starkiller-astro/XNet)
!! .
subroutine nuclear_heating_energy(rho,Ye,pf,Y_p,Y,entropy_p,entropy,T9_p,T9,T9_tr_p,T9_tr,eos_status)
    use global_class,    only: net_size,isotope
    use nucstuff_class,  only: inter_partf
    use parameter_class, only: engen, top_engen_every, unit, &
                               heating_frac, freeze_rate_temp, heating_mode
    use expansion_module,only: expand_count
    implicit none
 
    real(r_kind),intent(in)                          :: rho          !< density
    real(r_kind),intent(in)                          :: ye           !< electron fraction
    real(r_kind),dimension(0:net_size),intent(out)   :: pf           !< partition functions
    real(r_kind),dimension(net_size),intent(in)      :: y_p          !< abundances
    real(r_kind),dimension(net_size),intent(in)      :: y            !< abundances at the next step
    real(r_kind),intent(in)                          :: entropy_p    !< entropy
    real(r_kind),intent(out)                         :: entropy      !< entropy for the next step
    real(r_kind),intent(inout)                       :: t9           !< temperature at the next step
    real(r_kind),intent(in)                          :: t9_p         !< temperature at the previous step
    real(r_kind),intent(in)                          :: t9_tr_p      !< temperature of the trajectory/analytic (prev.)
    real(r_kind),intent(in)                          :: t9_tr        !< temperature of the trajectory/analytic
    integer,intent(out)                              :: eos_status   !< status of the EOS call
 
    real(r_kind) :: mic2, mui, kbt     !< auxiliary variables and constants
    type(timmes_eos_state) :: state    !< Timmes EOS state
    real(r_kind) :: gi                 !< statistical weight
    real(r_kind) :: twopi_h2c2         !< 2\pi/(hc)^2 [MeV^-2 cm^-2]
    real(r_kind) :: entropy_inc        !< entropy increment during single timestep due to nuclear heating
    real(r_kind) :: tdot_i
    real(r_kind) :: tdot
    real(r_kind) :: add_incr
 
    real(r_kind) :: ysum,yasum, temp_next, cv
    integer :: i ! Loop variables
 
 
    ! Call the Timmes EOS to get the specific heat at constant volume
    ysum= sum(y(1:net_size))
    yasum= sum(y(1:net_size)*isotope(1:net_size)%mass)
    state%abar= yasum / ysum
    call timmes_eos(ink,t9*1d9,rho,ye,state,eos_status)
 
    ! If this fails, return
    if (eos_status .ne. 0) then
        return
    end if
 
    ! Calculate the energy increment due to nuclear heating
    tdot = 0
    do i=1,net_size
        tdot_i = -(isotope(i)%mass_exc)*(y(i)-y_p(i)) / (state%cv)
        tdot = tdot + tdot_i
    end do
 
    ! Radiate energy away
    tdot = tdot - (qdot_bet+qdot_nu+qdot_th) / (state%cv)
 
 
    ! Only heat when the temperature is above freeze_rate_temp
    if (t9 .gt. freeze_rate_temp*1.000001)then
 
        if (heating_mode .eq. 2) then
            ! Heating mode 2 follows the trajectory until it ends,
            ! afterwards adiabatic expansion is assumed
            if (expand_count .ge. 1) then
                call timmes_eos(ins,entropy_p,rho,ye,state,eos_status)
 
                ! If this fails, return
                if (eos_status .ne. 0) then
                    return
                end if
 
                add_incr = state%t*1d-9 - t9_p
            else
                add_incr = (t9_tr-t9_tr_p)
            end if
        elseif (heating_mode .eq. 3) then
            ! Heating mode 3 assumes adiabatic expansion
            call timmes_eos(ins,entropy_p,rho,ye,state,eos_status)
 
            ! If this fails, return
            if (eos_status .ne. 0) then
                return
            end if
 
            add_incr = state%t*1d-9 - t9_p
        end if
 
        t9 = t9_p + (tdot*1e-9) + add_incr
 
     else
        eos_status = 0
        t9 = t9
     end if
     ! Todo, maybe add the following as debug file
     !  write(55,'(100(es23.16,2X))')T9,T9_p,T9_tr,T9_tr_p,(heating_frac*engen*1e-9), (T9_tr-T9_tr_p),T9-T9_p
  end subroutine nuclear_heating_energy
 
 
 
!>
!! Modifies the temperature and ent_p to account for nuclear heating
!!
!! Here, an additional equation is calculated for the entropy.
!! The entropy is calculated by:
!! \f[ E_\mathrm{gen} = \sum_i (m_i c^2 + \mu_i + Z_i \mu_e) dY_i  -\Delta q \f]
!! , where \f$ \mu_i \f$ is given by
!! \f[
!! \mu_i = -k_B T \log \left( \frac{g_i G_i}{\rho N_A Y_i}
!!          \left(\frac{2 \pi m_ic^2 k_B T}{h^2 c^2}\right)^{1.5} \right)
!! \f]
!! and the nuclear mass defined by the atomic mass excess:
!! \f[ m_ic^2 = M_i +Am_uc^2 - Z_i m_e c^2 .\f]
!! The symbols are defined as:
!! <pre style="line-height:.5">
!!  - M_i         : mass excess [MeV] (see \ref global_class::isotope\%mass_exc)
!!  - m_ic^2      : nuclear mass (see \ref global_class::isotope\%mass [MeV])
!!  - kB          : Boltzman constant (see \ref parameter_class::unit\%k_mev)
!!  - g_i         : 2 * spin + 1 ( see \ref global_class::isotope\%spin )
!!  - G_i         : normalized partition function of isotope i
!!  - rho         : Density [cgs]
!!  - \f$mu_e\f$  : Electron chemical potential
!!  - Z_i         : Proton number of isotope i
!!  - Navo        : Avogadro constant (see \ref parameter_class::unit\%n_a)
!!  - Y_i         : Abundance [mol/g]
!!  - h           : Planck constant [MeV*s] (see \ref parameter_class::unit\%h_mevs)
!!  - c           : Speed of light [cm/s] (see \ref parameter_class::unit\%clight)
!! .
!! </pre>
!! The entropy is converted to a temperature via the equation of state at the
!! end of this subroutine. Furthermore, \f$ \Delta q \f$ is the energy loss due
!! neutrinos that won't thermalize (\ref parameter_class::heating_frac ,
!! \ref beta_decay_file , \ref neutrino_loss_file).
!!
!! @see
!! [Freiburghaus et al. 1999, Eq. 2](https://ui.adsabs.harvard.edu/abs/1999ApJ...525L.121F/abstract)
!! [Mueller 1986](https://ui.adsabs.harvard.edu/abs/1986A%26A...162..103M/abstract)
!!
!! @author O. Korobkin
!!
!! \b Edited:
!!    - MR: 20.01.21 - Separated top contributor to independent subroutine
!! .
subroutine nuclear_heating_entropy(rho,Ye,pf,Y_p,Y,entropy_p,entropy,T9,eos_status)
    use global_class,    only: net_size,isotope
    use nucstuff_class,  only: inter_partf
    use parameter_class, only: engen, top_engen_every, unit, &
                               heating_frac, freeze_rate_temp
    implicit none
 
    real(r_kind),intent(in)                          :: rho          !< density
    real(r_kind),intent(in)                          :: ye           !< electron fraction
    real(r_kind),dimension(0:net_size),intent(out)   :: pf           !< partition functions
    real(r_kind),dimension(net_size),intent(in)      :: y_p          !< abundances
    real(r_kind),dimension(net_size),intent(in)      :: y            !< abundances at the next step
    real(r_kind),intent(in)                          :: entropy_p    !< entropy
    real(r_kind),intent(out)                         :: entropy      !< entropy for the next step
    real(r_kind),intent(inout)                       :: t9           !< temperature at the next step
    integer,intent(out)                              :: eos_status   !< Status of the EOS call
 
    real(r_kind) :: mic2, mui, kbt     !< auxiliary variables and constants
    type(timmes_eos_state) :: state    !< Timmes EOS state
    real(r_kind) :: gi                 !< statistical weight
    real(r_kind) :: twopi_h2c2         !< 2\pi/(hc)^2 [MeV^-2 cm^-2]
    real(r_kind) :: entropy_inc        !< entropy increment during single timestep due to nuclear heating
    real(r_kind) :: engen_i            !< generated energy (per baryon) due to abundance changes of species i
    real(r_kind) :: etaele,etapos      !< Electron and positron chemical potentials from eos
    real(r_kind) :: mue                !< Electron and positron chemical potentials
 
    real(r_kind) :: ysum,yasum, temp_next
    integer :: i ! Loop variables
 
    info_entry("nuclear_heating_entropy")
 
    ! calculate partition functions for all nuclei
    call inter_partf (t9, pf)
 
    temp_next= t9 * 1.d9
 
 
    ! Get the electron chemical potential
    call chempot(temp_next,rho,ye,etaele,etapos)
    mue = etaele * unit%k_mev * temp_next
    mue = mue + unit%mass_e
    ! calculate <A> based on abundances of the next timestep; this is needed for the EOS later
    ysum= sum(y(1:net_size))
    yasum= sum(y(1:net_size)*isotope(1:net_size)%mass)
    state%abar= yasum / ysum
 
    ! update the temperature, using the entropy increment
    !  entropy  = engen / T
    ! [entropy] = [MeV/baryon] / [K] -> kB [MeV/K] / baryon
    ! calculate energy generated per baryon:
    !
    ! engen = sum_i (m_ic^2 + \mu_i + Z_i \mu_e) dY_i - \Delta q
    !
    ! where the chemical potential is:
    !                   {   g_i G_i       2 \pi m_ic^2 kB T  3/2 }
    ! \mu_i = -kB T log { ------------- ( ---------------- )     }
    !                   { rho Navo Y_i        h^2 c^2            }
    !
    ! M_i     = isotope(i)%mass_exc [MeV]
    ! m_ic^2  = isotope(i)%mass [MeV]
    ! kB      = 8.6173324d-11 [Mev/K]    = unit%k_mev
    ! kB T    = 8.6173324d-11 [Mev/K] * T [K] -> [MeV]
    ! mu_e    = Electron chemical potential [MeV]
    ! g_i     = isotope(i)%spin * 2 + 1
    ! G_i     = normalized partition function of isotope i
    ! rho     = rho [cgs]
    ! Navo    = 6.0221367d23 [1/mol]     = unit%n_a
    ! Y_i     = Y(i) [mol/g] ... me: but also [1/baryon]?
    ! h       = 4.135667516e-21 [MeV*s]  = unit%h_mevs
    ! c       = 2.99792458d10 [cm/s]     = unit%clight
    !
    ! [engen] = MeV/nuc * nuc/baryon = MeV / baryon
 
    engen = 0.d0
    kbt = temp_next * 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) / (rho*y(i)*unit%n_a) * &
               (twopi_h2c2 * isotope(i)%mass*unit%mass_u * kbt)**1.5d0 )
               ! note that Amuc2 is a approximation of mic2. It is used in order to be consistent with what is
               ! used in NSE.
 
          engen_i = -(mic2 + mui + isotope(i)%p_nr*mue)*(y(i)-y_p(i))  ! energy generated by species i [MeV/baryon]
          engen = engen + engen_i
       endif
    end do
 
    ! the total energy from nuc. reactions
    engen_nuc = engen_nuc + engen / unit%hix ! [erg/g]
 
    if (t9 .gt. freeze_rate_temp*1.000001)then
      ! use only a fraction (=heating_frac) of engen: the rest is radiated away, e.g., due to neutrino cooling
      entropy_inc= (engen - (qdot_bet + qdot_nu + qdot_th)) / kbt  ! divide by kB in order to get entropy in units kB per baryon
      entropy= entropy_p + entropy_inc         ! [kB/baryon]
 
      ! Prevent unphysical negative entropy
      ! This mostly happens for regions where the EOS is not
      ! valid in any case, i.e., very high densities.
      if (entropy .lt. 0) then
        ! Complain about it at high verbose levels
        if (verbose_level .gt. 1) then
            write(*,*) "Warning! Restricting entropy as it was negative: ", entropy
        end if
        entropy = 0
      end if
 
      call timmes_eos(ins,entropy,rho,ye,state,eos_status)
 
      if (eos_status .ne. 0 ) then
        continue
        ! Do nothing here as it can only happen for a bad NR iteration
        ! call raise_exception("The EOS did not converge, "//NEW_LINE("A")//&
        !                      "T   ="//num_to_str(T9)//"GK, "//NEW_LINE("A")//&
        !                      "rho ="//num_to_str(rho)//"g/ccm,"//NEW_LINE("A")//&
        !                      "Ye  ="//num_to_str(Ye), "nuclear_heating_update",&
        !                      310003)
      else
        temp_next= state%t     ! new temperature with updated entropy
        t9 = temp_next / 1d9   ! new temperature in [GK]
      end if
      t9 = temp_next / 1d9     ! new temperature in [GK]
    else
      eos_status = 0
      t9= t9
    end if
 
  info_exit("nuclear_heating_entropy")
 
 
end subroutine nuclear_heating_entropy
 
 
 
!>
!! Output data related to nuclear heating
!!
!! For example the energy generation in the file _engen.dat_
!! <div style="width:800px;overflow:scroll;padding:5px;
!! background-color:#e8e8e8;scrollbar-base-color:#DEBB07;border:1px dashed grey;
!! display: inline-block">
!! <pre>
!! # 1:iteration 2:time[s] 3:timestep[s] 4:E_gen(tot) [erg/g] 5:E_gen(t) [erg/(g*s )] 6:E_gen(t) [MeV/(baryon*s)]
!!         0  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
!!         0  2.000000000000000E-12  2.000000000000000E-12  2.100942703338183E+08  5.252356611052194E+19  5.443684205322992E+01
!!         1  6.000000000000000E-12  4.000000000000000E-12  6.302695988823750E+08  5.252191701420075E+19  5.443513288527704E+01
!! ...</pre></div>
!!
!! \b Edited:
!!   - MR: 20.1.21  - Modified format of _toplist.dat_
!!   - MR: 07.5.24  - Removed _toplist.dat_ from here
!! .
subroutine output_nuclear_heating(cnt,ctime)
  use single_zone_vars, only: stepsize
  use parameter_class,  only: engen, unit
  implicit none
  integer,intent(in)      :: cnt   !< current iteration
  real(r_kind),intent(in) :: ctime !< current global time [s]
  !
  real(r_kind)  :: engen_per_s
 
  info_entry("output_nuclear_heating")
  ! output nuclear energy generation
  if ((stepsize .gt. 0.d0)) then  ! engen calculated in nuclear_heating_update
     engen_per_s = engen / stepsize
  else
     engen_per_s = engen                                  ! engen from jacobian_init (already per second) and first entry
  end if
 
  ! output nuclear energy generation
  write(engen_unit,'(i10,5es23.15)') cnt, ctime, stepsize, engen_nuc, engen_per_s / unit%hix, engen_per_s
 
 
  info_exit("output_nuclear_heating")
 
end subroutine output_nuclear_heating
 
 
end module nuclear_heating