Write Hydrogen-Impurity reaction rates to a HDF5 file
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer(kind=HID_T), | intent(inout) | :: | id | HDF5 file or group object id |
||
| character(len=*), | intent(in) | :: | namelist_file | Namelist file that contains settings |
||
| integer, | intent(in) | :: | n_max | Number of initial atomic energy states to calculate |
||
| integer, | intent(in) | :: | m_max | Number of final atomic energy states to calculate |
subroutine write_bt_H_Aq(id, namelist_file, n_max, m_max)
!+ Write Hydrogen-Impurity reaction rates to a HDF5 file
integer(HID_T), intent(inout) :: id
!+ HDF5 file or group object id
character(len=*), intent(in) :: namelist_file
!+ Namelist file that contains settings
integer, intent(in) :: n_max
!+ Number of initial atomic energy states to calculate
integer, intent(in) :: m_max
!+ Number of final atomic energy states to calculate
logical :: calculate
integer :: q(10) = 0
real(Float64) :: mass
real(Float64) :: emin
real(Float64) :: emax
integer :: nenergy
real(Float64) :: tmin
real(Float64) :: tmax
integer :: ntemp
real(Float64) :: eb
real(Float64) :: dlogE
real(Float64), dimension(:), allocatable :: ebarr
real(Float64) :: ti
real(Float64) :: dlogT
real(Float64), dimension(:), allocatable :: tarr
real(Float64), dimension(:,:,:,:), allocatable :: ioniz, cx
real(Float64), dimension(:,:,:,:,:), allocatable :: excit
integer(HID_T) :: gid
integer(HSIZE_T), dimension(1) :: dim1
integer(HSIZE_T), dimension(2) :: dim2
integer(HSIZE_T), dimension(3) :: dim3
integer(HSIZE_T), dimension(4) :: dim4
integer(HSIZE_T), dimension(5) :: dim5
integer :: iq, ie, it, ia, n, m, error, cnt
real(Float64) :: rate
integer, parameter :: n_bt_amu = 3
real(Float64), dimension(2,n_bt_amu) :: a
character(len=10) :: aname
character(len=5) :: asym
logical :: exis
NAMELIST /H_Aq_rates/ calculate, q, mass, nenergy, emin, emax, ntemp, tmin, tmax
calculate = .True. ; q(1) = 6 ; mass = C_amu
nenergy = 100; emin = 1.d-3 ; emax = 4.d2
ntemp = 100; tmin = 1.d-3 ; tmax = 2.d1
inquire(file=namelist_file,exist=exis)
if(.not.exis) then
write(*,'(a,a)') 'WRITE_BT_H_Aq: Input file does not exist: ',trim(namelist_file)
write(*,'(a)') 'Continuing with default settings...'
else
open(13,file=namelist_file)
read(13,NML=H_Aq_rates)
close(13)
endif
if(.not.calculate) return
q_loop: do iq=1, size(q)
if(q(iq).eq.0) cycle q_loop
allocate(ebarr(nenergy))
allocate(tarr(ntemp))
allocate(ioniz(n_max,nenergy,ntemp,n_bt_amu))
allocate(cx(n_max,nenergy,ntemp,n_bt_amu))
allocate(excit(n_max,m_max,nenergy,ntemp,n_bt_amu))
select case (q(iq))
case (5)
aname = "Boron"
asym = "H_B5"
case (6)
aname = "Carbon"
asym = "H_C6"
case DEFAULT
write(aname,'("Impurity-",i1)') q(iq)
write(asym,'("H_A",i1)') q(iq)
end select
ebarr = 0.d0
ioniz = 0.d0
cx = 0.d0
excit = 0.d0
a(:,1) = [H1_amu, mass]
a(:,2) = [H2_amu, mass]
a(:,3) = [H3_amu, mass]
dlogE = (log10(emax) - log10(emin))/(nenergy - 1)
do ie=1, nenergy
ebarr(ie) = 10.d0**(log10(emin) + (ie-1)*dlogE)
enddo
dlogT = (log10(tmax) - log10(tmin))/(ntemp - 1)
do it=1, ntemp
tarr(it) = 10.d0**(log10(tmin) + (it-1)*dlogT)
enddo
if(verbose) then
write(*,'(a)') "---- H-"//trim(adjustl(aname))//" reaction rates settings ----"
write(*,'(T2,"q = ", i2)') q(iq)
write(*,'(T2,"mass = ",f7.2, " amu")') mass
write(*,'(T2,"Emin = ",e9.2, " keV")') emin
write(*,'(T2,"Emax = ",e9.2, " keV")') emax
write(*,'(T2,"Nenergy = ", i4)') nenergy
write(*,'(T2,"Tmin = ",e9.2, " keV")') tmin
write(*,'(T2,"Tmax = ",e9.2, " keV")') tmax
write(*,'(T2,"Ntemp = ", i4)') ntemp
write(*,*) ''
endif
cnt = 0
!$OMP PARALLEL DO private(ie, it, ia, n, m, eb, ti, rate)
do ie=istart, nenergy, istep
eb = ebarr(ie)
do it=1, ntemp
ti = tarr(it)
do ia=1, n_bt_amu
do n=1, n_max
do m=1, m_max
if(m.gt.n) then
call bt_maxwellian(Aq_excit_n_m, q(iq), ti, eb, &
a(2,ia), a(1,ia), n, m, rate)
excit(n,m,ie,it,ia) = rate
call bt_maxwellian(Aq_excit_n_m, q(iq), ti, eb, &
a(2,ia), a(1,ia), n, m, &
rate, deexcit=.True.)
excit(m,n,ie,it,ia) = rate
endif
enddo
call bt_maxwellian(Aq_cx_n, q(iq), ti, eb, &
a(2,ia), a(1,ia), n, rate)
cx(n,ie,it,ia) = rate
call bt_maxwellian(Aq_ioniz_n, q(iq), ti, eb, &
a(2,ia), a(1,ia), n, rate)
ioniz(n,ie,it,ia) = rate
enddo
enddo
cnt = cnt + 1
if(verbose) WRITE(*,'(f7.2,"%",a,$)') 100*istep*cnt/real(nenergy*ntemp),char(13)
enddo
enddo
!$OMP END PARALLEL DO
#ifdef _MPI
call parallel_sum(cx)
call parallel_sum(excit)
call parallel_sum(ioniz)
#endif
if(verbose) then
call h5gcreate_f(id, trim(adjustl(asym)), gid, error)
dim1 = [1]
dim4 = [n_max, nenergy, ntemp, n_bt_amu]
dim5 = [n_max, m_max, nenergy, ntemp, n_bt_amu]
call h5ltmake_dataset_int_f(gid, "nenergy", 0, dim1, [nenergy], error)
call h5ltmake_dataset_int_f(gid, "ntemp", 0, dim1, [ntemp], error)
call h5ltmake_dataset_int_f(gid, "n_bt_amu", 0, dim1, [n_bt_amu], error)
call h5ltmake_dataset_int_f(gid, "n_max", 0, dim1, [n_max], error)
call h5ltmake_dataset_int_f(gid, "m_max", 0, dim1, [m_max], error)
call h5ltmake_dataset_double_f(gid, "dlogE", 0, dim1, [dlogE], error)
call h5ltmake_dataset_double_f(gid, "emin", 0, dim1, [emin], error)
call h5ltmake_dataset_double_f(gid, "emax", 0, dim1, [emax], error)
call h5ltmake_dataset_double_f(gid, "dlogT", 0, dim1, [dlogT], error)
call h5ltmake_dataset_double_f(gid, "tmin", 0, dim1, [tmin], error)
call h5ltmake_dataset_double_f(gid, "tmax", 0, dim1, [tmax], error)
dim2 = [2,n_bt_amu]
call h5ltmake_compressed_dataset_double_f(gid, "bt_amu", 2, dim2, a, error)
dim1 = [nenergy]
call h5ltmake_compressed_dataset_double_f(gid, "energy", 1, dim1, ebarr, error)
dim1 = [ntemp]
call h5ltmake_compressed_dataset_double_f(gid, "temperature", 1, dim1, tarr, error)
call h5ltmake_compressed_dataset_double_f(gid, "cx", 4, dim4, cx, error)
call h5ltmake_compressed_dataset_double_f(gid, "ionization", 4, dim4, ioniz, error)
call h5ltmake_compressed_dataset_double_f(gid, "excitation", 5, dim5, excit, error)
call h5ltset_attribute_string_f(id, trim(adjustl(asym)), "description", &
"Beam-Target reaction rates for Hydrogen(beam)-"//trim(adjustl(aname))// &
"(target) interactions", error)
call h5ltset_attribute_string_f(gid, "nenergy", "description", &
"Number of energy values", error)
call h5ltset_attribute_string_f(gid, "ntemp", "description", &
"Number of target temperature values", error)
call h5ltset_attribute_string_f(gid, "n_bt_amu", "description", &
"Number of beam-target amu combinations", error)
call h5ltset_attribute_string_f(gid, "n_max", "description", &
"Number of initial energy levels", error)
call h5ltset_attribute_string_f(gid, "m_max", "description", &
"Number of final energy levels", error)
call h5ltset_attribute_string_f(gid, "bt_amu", "description", &
"Combinations of beam-target amu's e.g. b_amu, t_amu = bt_amu[:,i]", error)
call h5ltset_attribute_string_f(gid, "energy", "description", &
"Energy values", error)
call h5ltset_attribute_string_f(gid, "energy", "units", "keV", error)
call h5ltset_attribute_string_f(gid, "dlogE", "description", &
"Energy spacing in log-10", error)
call h5ltset_attribute_string_f(gid, "dlogE", "units", "log10(keV)", error)
call h5ltset_attribute_string_f(gid, "emin","description", &
"Minimum energy", error)
call h5ltset_attribute_string_f(gid, "emin", "units", "keV", error)
call h5ltset_attribute_string_f(gid, "emax","description", &
"Maximum energy", error)
call h5ltset_attribute_string_f(gid, "emax", "units", "keV", error)
call h5ltset_attribute_string_f(gid, "temperature", "description", &
"Target temperature values", error)
call h5ltset_attribute_string_f(gid, "temperature", "units", "keV", error)
call h5ltset_attribute_string_f(gid, "dlogT", "description", &
"Temperature spacing in log-10", error)
call h5ltset_attribute_string_f(gid, "dlogT", "units", "log10(keV)", error)
call h5ltset_attribute_string_f(gid, "tmin","description", &
"Minimum temperature", error)
call h5ltset_attribute_string_f(gid, "tmin", "units", "keV", error)
call h5ltset_attribute_string_f(gid, "tmax","description", &
"Maximum temperature", error)
call h5ltset_attribute_string_f(gid, "tmax", "units", "keV", error)
call h5ltset_attribute_string_f(gid, "cx", "description", &
"n-resolved charge exchange reaction rates: cx(n,energy,temp,bt_amu)", error)
call h5ltset_attribute_string_f(gid, "cx", "units", "cm^3/s", error)
call h5ltset_attribute_string_f(gid, "cx", "reaction", &
"A(q+) + H(n) -> A((q-1)+) + H(+)", error)
call h5ltset_attribute_string_f(gid, "excitation", "description", &
"n/m resolved (de-)excitation reaction rates: excitation(n,m,energy,temp,bt_amu)", error)
call h5ltset_attribute_string_f(gid, "excitation", "units", "cm^3/s", error)
call h5ltset_attribute_string_f(gid, "excitation", "reaction", &
"A(q+) + H(n) -> A(q+) + H(m); m > n excitation, m < n de-excitation", error)
call h5ltset_attribute_string_f(gid, "ionization", "description", &
"n resolved ionization reaction rates: ionization(n,energy,temp,bt_amu)", error)
call h5ltset_attribute_string_f(gid, "ionization", "units", "cm^3/s", error)
call h5ltset_attribute_string_f(gid, "ionization", "reaction", &
"A(q+) + H(n) -> A(q+) + H(+) + e-", error)
call h5gclose_f(gid, error)
endif
deallocate(ebarr, tarr, excit, ioniz, cx)
enddo q_loop
end subroutine write_bt_H_Aq