Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 Ed Brown & The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module chem_lib
21 :
22 : use chem_def, only: chem_has_been_initialized
23 : use const_def, only: dp
24 : use math_lib
25 :
26 : implicit none
27 :
28 : contains
29 :
30 18 : subroutine chem_init(isotopes_filename, ierr)
31 : ! uses mesa_data_dir from const_def
32 : use chem_def
33 : use chem_isos_io, only: do_read_chem_isos
34 : use lodders_mod, only : read_lodders03_data
35 : character (len=*), intent(in) :: isotopes_filename
36 : integer, intent(out) :: ierr
37 :
38 18 : ierr = 0
39 18 : if (chem_has_been_initialized) return
40 16 : call do_read_chem_isos(isotopes_filename, ierr)
41 16 : if (ierr /= 0) then
42 0 : write(*,*) 'failed in do_read_chem_isos'
43 0 : return
44 : end if
45 16 : call init_chem_tables
46 16 : call read_lodders03_data('lodders03.data',ierr)
47 16 : if (ierr /= 0) then
48 0 : write(*,*) 'failed in read_lodders03_data'
49 0 : return
50 : end if
51 16 : chem_has_been_initialized = .true.
52 16 : call set_some_isos
53 :
54 18 : end subroutine chem_init
55 :
56 :
57 1 : subroutine chem_shutdown ()
58 :
59 18 : use chem_def
60 : use utils_lib, only: integer_dict_free
61 :
62 1 : if (.NOT. chem_has_been_initialized) return
63 :
64 1 : call free_lodders03_table()
65 1 : call free_nuclide_data(chem_isos)
66 :
67 : ! Free dictionaries
68 :
69 1 : if (ASSOCIATED(Xsol_names_dict)) call integer_dict_free(Xsol_names_dict)
70 1 : if (ASSOCIATED(chem_element_names_dict)) call integer_dict_free(chem_element_names_dict)
71 1 : if (ASSOCIATED(chem_isos_dict)) call integer_dict_free(chem_isos_dict)
72 1 : if (ASSOCIATED(category_names_dict)) call integer_dict_free(category_names_dict)
73 :
74 1 : chem_has_been_initialized = .FALSE.
75 :
76 1 : end subroutine chem_shutdown
77 :
78 :
79 572 : function lodders03_element_atom_percent(nuclei) result(percent)
80 : ! Katharina Lodders,
81 : ! "Solar System Abundances and Condensation Temperatures of the Elements",
82 : ! ApJ, 591, 1220-1247 (2003).
83 : ! Table 6: Abundances of the Isotopes in the Solar System
84 :
85 : ! These are element atom percentages (i.e., by number, not by mass)
86 :
87 : ! NOTE: The data here stops at ge -- the table in the paper goes to 92U
88 :
89 : ! TO DO: add the rest of the info from the paper
90 1 : use chem_def
91 : use lodders_mod
92 : character(len=*), intent(in) :: nuclei
93 : real(dp) :: percent
94 : integer :: ierr
95 :
96 286 : if (.not. chem_has_been_initialized) then
97 0 : write(*,*) 'must call chem_init before calling any other routine in chem_lib'
98 0 : percent = -1.0d0
99 0 : return
100 : end if
101 286 : percent = get_lodders03_isotopic_abundance(nuclei, ierr)
102 286 : if (ierr /= 0) percent = 0.0d0
103 286 : end function lodders03_element_atom_percent
104 :
105 :
106 : ! returns the index of a particular nuclide in a particular set
107 : ! returns nuclide_not_found if name not found
108 4143078 : function get_nuclide_index_in_set(nuclei, set) result(indx)
109 286 : use chem_def
110 : use nuclide_set_mod
111 : character(len=*), intent(in) :: nuclei
112 : type(nuclide_set), dimension(:), intent(in) :: set
113 : integer :: indx
114 : character(len=iso_name_length) :: name
115 1381026 : name = nuclei
116 1381026 : indx = rank_in_set(name, set)
117 2762052 : end function get_nuclide_index_in_set
118 :
119 :
120 5 : subroutine generate_nuclide_set(names, set)
121 1381026 : use chem_def
122 : use nuclide_set_mod, only: sort_nuclide_set
123 : character(len=iso_name_length), dimension(:), intent(in) :: names
124 : type(nuclide_set), dimension(size(names)), intent(out) :: set
125 : integer :: i
126 78570 : set = [(nuclide_set(names(i), i), i=1, size(names))]
127 5 : call sort_nuclide_set(set)
128 5 : end subroutine generate_nuclide_set
129 :
130 :
131 394361 : subroutine basic_composition_info( &
132 394361 : num_isos, chem_id, x, xh, xhe, z, &
133 : abar, zbar, z2bar, z53bar, ye, mass_correction, sumx)
134 : integer, intent(in) :: num_isos
135 : integer, intent(in) :: chem_id(:) ! (num_isos) ! the nuclide indices for the entries in x
136 : real(dp), intent(in) :: x(:) ! (num_isos) ! baryon fractions. should sum to 1.0
137 : real(dp), intent(out) :: &
138 : xh, xhe, z, abar, zbar, z2bar, z53bar, ye, mass_correction, sumx
139 : real(dp), dimension(0) :: dabar_dx, dzbar_dx, dmc_dx
140 : call get_composition_info( &
141 : num_isos, chem_id, x, xh, xhe, z, &
142 : abar, zbar, z2bar, z53bar, ye, mass_correction, &
143 394361 : sumx, .true., dabar_dx, dzbar_dx, dmc_dx)
144 5 : end subroutine basic_composition_info
145 :
146 :
147 1019 : subroutine composition_info( &
148 1019 : num_isos, chem_id, x, xh, xhe, z, &
149 : abar, zbar, z2bar, z53bar, ye, mass_correction, &
150 1019 : sumx, dabar_dx, dzbar_dx, dmc_dx)
151 : integer, intent(in) :: num_isos
152 : integer, intent(in) :: chem_id(:) ! (num_isos) ! the nuclide indices for the entries in x
153 : real(dp), intent(in) :: x(:) ! (num_isos) ! baryon fractions. should sum to 1.0
154 : real(dp), intent(out) :: &
155 : xh, xhe, z, abar, zbar, z2bar, z53bar, ye, mass_correction, sumx
156 : real(dp), dimension(:) :: dabar_dx, dzbar_dx, dmc_dx
157 : call get_composition_info( &
158 : num_isos, chem_id, x, xh, xhe, z, &
159 : abar, zbar, z2bar, z53bar, ye, mass_correction, &
160 1019 : sumx, .false., dabar_dx, dzbar_dx, dmc_dx)
161 1019 : end subroutine composition_info
162 :
163 :
164 395380 : subroutine get_composition_info( &
165 395380 : num_isos, chem_id, x, xh, xhe, xz, &
166 : abar, zbar, z2bar, z53bar, ye, mass_correction, &
167 395380 : sumx, skip_partials, dabar_dx, dzbar_dx, dmc_dx)
168 :
169 : ! here's a reminder of definitions:
170 : ! X(i) ion baryon fraction (g/g)
171 : ! A(i) ion atomic mass number
172 : ! W(i) ion atomic weight (g/mole)
173 : ! Z(i) ion charge (number of protons)
174 : ! Y(i) = X(i)/A(i), ion abundance
175 : ! n(i) = rho*avo*Y(i), ion number density (g/cm^3)*(#/mole)*(mole/g) -> (#/cm^3)
176 :
177 : ! abar = sum(n(i)*A(i))/sum(n(i)) -- average atomic mass number
178 : ! zbar = sum(n(i)*Z(i))/sum(n(i)) -- average charge number
179 : ! z2bar = sum(n(i)*Z(i)^2)/sum(n(i)) -- average charge^2
180 : ! n = rho*avo/abar -- total number density (#/cm^3)
181 : ! ye = zbar/abar -- average charge per baryon = proton fraction
182 : ! xh = mass fraction hydrogen
183 : ! xhe = mass fraction helium
184 : ! xz = mass fraction metals
185 : ! mass_correction = sum(n(i)*W(i))/sum(n(i)*A(i)) --
186 : ! (mass density) = (baryon density) * m_u * mass_correction
187 :
188 : use chem_def
189 : integer, intent(in) :: num_isos
190 : integer, intent(in) :: chem_id(:) ! (num_isos) ! the nuclide indices for the entries in x
191 : real(dp), intent(in) :: x(:) ! (num_isos) ! baryon fractions. should sum to 1.0
192 : real(dp), intent(out) :: &
193 : xh, xhe, xz, abar, zbar, z2bar, z53bar, ye, mass_correction, sumx
194 : logical, intent(in) :: skip_partials
195 : real(dp), dimension(:) :: dabar_dx, dzbar_dx, dmc_dx
196 :
197 13943148 : real(dp), dimension(num_isos) :: y, z, w, a
198 : integer :: i, cid, iz
199 395380 : if (.not. chem_has_been_initialized) then
200 0 : write(*,*) 'must call chem_init before calling any other routine in chem_lib'
201 0 : xh=0; xhe=0; abar=0; zbar=0; z2bar=0; z53bar=0; ye=0; mass_correction = 0; sumx=0
202 0 : return
203 : end if
204 395380 : xh = 0d0; xhe = 0d0
205 3782322 : do i=1,num_isos
206 3386942 : cid = chem_id(i)
207 3386942 : iz = chem_isos% Z(cid)
208 3386942 : z(i) = iz
209 3386942 : y(i) = x(i)/dble(chem_isos% Z_plus_N(cid))
210 3386942 : w(i) = chem_isos% W(cid)
211 3386942 : a(i) = chem_isos% Z_plus_N(cid)
212 395380 : select case(iz)
213 : case (1)
214 413539 : xh = xh + x(i)
215 : case (2)
216 3386942 : xhe = xhe + x(i)
217 : end select
218 : end do
219 :
220 395380 : xz = max(0d0, min(1d0, 1d0 - (xh + xhe)))
221 3782322 : sumx = sum(x(1:num_isos)) ! this should be one, always, since we define x as a baryon fraction
222 3782322 : abar = sumx / sum(y(1:num_isos))
223 3782322 : zbar = sum(y(1:num_isos)*z(1:num_isos)) * abar
224 3782322 : z2bar = sum(y(1:num_isos)*z(1:num_isos)*z(1:num_isos)) * abar
225 : z53bar = sum(y(1:num_isos)*z(1:num_isos)*z(1:num_isos)) * abar
226 395380 : ye = zbar/abar ! assume complete ionization
227 3782322 : mass_correction = sum(y(1:num_isos)*w(1:num_isos))
228 :
229 395380 : z53bar = 0d0
230 3782322 : do i=1,num_isos
231 3782322 : z53bar = z53bar + y(i) * chem_isos% Z53(chem_id(i))
232 : end do
233 395380 : z53bar = z53bar * abar
234 :
235 395380 : if (skip_partials) return
236 :
237 20577 : do i=1,num_isos
238 19558 : dabar_dx(i) = abar*(a(i)-abar)/a(i)/sumx
239 19558 : dzbar_dx(i) = abar*(z(i)-zbar)/a(i)/sumx
240 20577 : dmc_dx(i) = w(i)/a(i) - mass_correction
241 : end do
242 :
243 395380 : end subroutine get_composition_info
244 :
245 :
246 : ! Q: Is positron annihilation actually included in Qtotal?
247 : ! A: (from Frank Timmes)
248 : !
249 : ! the formula used in the code is the atomic mass excess - not the nuclear mass excess -
250 : ! in terms of the binding energy and Deltas. as bill notes, this formulation
251 : ! implicitly takes care of the electron masses. i can see why some confusion may
252 : ! arise at first blush because of the notation used in the code - its says “Mp” but
253 : ! really means “Mh” where Mh = m_p + m_e.
254 : !
255 : ! p + p -> d + e^+ + nu
256 : !
257 : ! let’s do it by hand and then by code, neglecting binding energy terms of
258 : ! the atom (13.6 eV and such) as they are a million times smaller than the nuclear terms.
259 : !
260 : ! by hand:
261 : !
262 : ! mass excess left-hand-side = twice hydrogen mass excess =
263 : ! 2 * (m_h - 1 amu) = 2 (m_p + m_e - 1 amu)
264 : !
265 : ! mass excess right-hand-side =
266 : ! m_h + m_n - B(d) - 2 amu = m_p + m_e + m_n - B(d) - 2 amu
267 : !
268 : ! Q = left - right = m_p + m_e - m_n + B(d)
269 : !
270 : ! which may be written (adding and subtracting m_p)
271 : !
272 : ! Q = 2 m_p + m_e - (m_p + m_n - B(d))
273 : ! = (twice nuclear mass of proton) + (electron mass) - (nuclear mass of deuterium)
274 : !
275 : !
276 : ! by code (with interpretation):
277 : !
278 : ! first loop
279 : ! Q = -del_Mp = -(m_h - 1 amu) = -(m_p + m_e - 1 amu)
280 : ! Qtotal = del_Mp = m_p + m_e - 1 amu
281 : !
282 : ! second loop
283 : ! Q = -del_Mp = -(m_h - 1 amu) = -(m_p + m_e - 1 amu)
284 : ! Qtotal = 2 del_Mp = 2 (m_p + m_e - 1 amu) ! note this is the same as the left-hand side term above
285 : !
286 : !
287 : ! third loop
288 : ! Q = B(d) - del_Mp - del_Mn
289 : ! = B(d) - (m_h - 1 amu) - (m_n - 1 amu)
290 : ! = B(d) - (m_p + m_e - 1 amu) - (m_n - 1 amu)
291 : ! = B(d) - (mp + m_e + m_n - 2 amu) ! note this is the negative of the right-hand-side term above
292 : !
293 : ! Qtotal = 2 del_Mp + [B(d) - del_Mp - del_Mn]
294 : ! = m_p + m_e - m_n + B(d) ! same as above
295 : !
296 : !
297 : ! i think the confusion lay between the code nomenclature
298 : ! and nuclear vs atomic mass excesses.
299 :
300 1609 : real(dp) function reaction_Qtotal(num_in,num_out,reactants,nuclides)
301 395380 : use chem_def
302 : integer, intent(in) :: num_in,num_out,reactants(:)
303 : type(nuclide_data), intent(in) :: nuclides
304 : integer :: j, l
305 1609 : real(dp) :: Q
306 1609 : reaction_Qtotal = 0d0
307 6189 : do j=1,num_in+num_out
308 4580 : l = reactants(j)
309 4580 : Q = get_Q(nuclides,l)
310 6189 : if (j <= num_in) then
311 2522 : reaction_Qtotal = reaction_Qtotal - Q
312 : else
313 2058 : reaction_Qtotal = reaction_Qtotal + Q
314 : end if
315 :
316 : end do
317 1609 : end function reaction_Qtotal
318 :
319 :
320 51 : integer function chem_get_element_id(cname)
321 : ! NOTE: this is for elements like 'h', not for isotopes like 'h1'
322 : ! use chem_get_iso_id for looking up isotope names
323 1609 : use chem_def
324 : use utils_lib
325 : character (len=*), intent(in) :: cname
326 : ! name of the element (e.g. 'h', 'he', 'ne')
327 : ! same names as in chem_element_Name
328 : ! returns id for the element if there is a matching name
329 : ! returns 0 otherwise.
330 : integer :: ierr, value
331 51 : if (.not. chem_has_been_initialized) then
332 0 : write(*,*) 'must call chem_init before calling any other routine in chem_lib'
333 0 : chem_get_element_id = -1
334 0 : return
335 : end if
336 51 : call integer_dict_lookup(chem_element_names_dict, cname, value, ierr)
337 51 : if (ierr /= 0) value = -1
338 51 : chem_get_element_id = value
339 51 : end function chem_get_element_id
340 :
341 :
342 306 : real(dp) function chem_Xsol(nam)
343 : character (len=*), intent(in) :: nam
344 : ! name of the isotope (e.g. 'h1', 'he4', 'ne20')
345 153 : real(dp) :: z, a, xelem
346 : integer :: ierr
347 153 : if (.not. chem_has_been_initialized) then
348 0 : write(*,*) 'must call chem_init before calling any other routine in chem_lib'
349 0 : chem_Xsol = -1
350 0 : return
351 : end if
352 153 : call chem_get_solar(nam, z, a, xelem, ierr)
353 153 : if (ierr /= 0) then
354 : chem_Xsol = 0d0
355 : else
356 153 : chem_Xsol = xelem
357 : end if
358 51 : end function chem_Xsol
359 :
360 :
361 153 : subroutine chem_get_solar(nam, z, a, xelem, ierr)
362 : use chem_def
363 : use utils_lib
364 : ! returns data from Anders and Grevesse, 1989
365 : character (len=*), intent(in) :: nam
366 : ! name of the isotope (e.g. 'h1', 'he4', 'ne20')
367 : ! note that these names match those in the nuclear net library iso_Names array
368 : ! but some net isotopes are not here (ex/ be7, n13, o14, o15, f17, f18, ... fe52, ni56 )
369 : ! and some isotopes are here that are not in the nets (ex/ hf176)
370 : real(dp), intent(out) :: z ! charge
371 : real(dp), intent(out) :: a ! number of nucleons (protons and neutrons)
372 : real(dp), intent(out) :: xelem ! elemental mass fraction associated with this isotope
373 : integer, intent(out) :: ierr ! == 0 means AOK; == -1 means did not find the requested name
374 : integer :: i
375 153 : ierr = 0
376 153 : if (.not. chem_has_been_initialized) then
377 0 : write(*,*) 'must call chem_init before calling any other routine in chem_lib'
378 0 : ierr = -1; z = 0; a = 0; xelem = 0; return
379 : end if
380 153 : call integer_dict_lookup(Xsol_names_dict, nam, i, ierr)
381 153 : if (ierr /= 0) then
382 0 : z = 0; a = 0; xelem = 0; return
383 : end if
384 153 : z = dble(izsol(i))
385 153 : a = dble(iasol(i))
386 153 : xelem = solx(i)
387 153 : end subroutine chem_get_solar
388 :
389 :
390 :
391 : ! given an array of Z, A, returns an array of names in chem_isos format
392 0 : subroutine generate_nuclide_names(Z, A, names)
393 153 : use chem_def
394 : integer, dimension(:), intent(in) :: Z, A
395 : character(len=iso_name_length), dimension(size(Z)), intent(out) :: names
396 : integer :: i, ierr, count_isomer
397 : logical :: use_al26_isomers
398 :
399 0 : count_isomer = 0
400 0 : use_al26_isomers = .false.
401 0 : do i = 1, size(Z)
402 0 : if (A(i) == 26 .and. Z(i) == 13) count_isomer = count_isomer + 1
403 : end do
404 0 : if (count_isomer > 1) use_al26_isomers = .true.
405 :
406 0 : count_isomer = 1
407 0 : do i = 1, size(Z)
408 0 : select case(Z(i))
409 : case (0)
410 0 : names(i) = el_name(Z(i))
411 : case (1:max_el_z)
412 0 : write(names(i), '(a, i0)') trim(el_name(Z(i))), A(i)
413 0 : if (A(i) == 26 .and. Z(i) == 13 .and. use_al26_isomers) then
414 0 : count_isomer = count_isomer + 1
415 : !if (count_isomer > 1) names(i) = al_isomers(count_isomer)
416 0 : names(i) = al_isomers(count_isomer)
417 : end if
418 : case default
419 0 : write(*, '(a, i0, a, i0)') 'warning: ', Z(i), ' greater than Zmax = ', max_el_z
420 0 : names(i) = '*****'
421 0 : ierr = nuclide_not_found
422 : end select
423 0 : names(i) = adjustr(names(i))
424 : end do
425 0 : end subroutine generate_nuclide_names
426 :
427 :
428 0 : subroutine generate_long_nuclide_names(Z, A, long_names)
429 0 : use chem_def
430 : integer, dimension(:), intent(in) :: Z, A
431 : character(len=long_name_length), dimension(size(Z)), intent(out) :: long_names
432 : integer :: i, ierr, count_isomer
433 : logical :: use_al26_isomers
434 :
435 0 : count_isomer = 0
436 0 : use_al26_isomers = .false.
437 0 : do i = 1, size(Z)
438 0 : if (A(i) == 26 .and. Z(i) == 13) count_isomer = count_isomer + 1
439 : end do
440 0 : if (count_isomer > 1) use_al26_isomers = .true.
441 :
442 0 : count_isomer = 1
443 0 : do i = 1, size(Z)
444 0 : select case(Z(i))
445 : case (0) ! neutrons are special?
446 0 : long_names(i) = el_long_name(Z(i))
447 : case (1:max_el_z)
448 0 : write(long_names(i), '(a, "-", i0)') trim(el_long_name(Z(i))), A(i)
449 0 : if (A(i) == 26 .and. Z(i) == 13 .and. use_al26_isomers) then
450 0 : count_isomer = count_isomer + 1
451 : !if (count_isomer > 1) long_names(i) = long_al_isomers(count_isomer)
452 0 : long_names(i) = long_al_isomers(count_isomer)
453 : end if
454 : case default
455 0 : write(*, '(a, i0, a, i0)') 'warning: ', Z(i), ' greater than Zmax = ', max_el_z
456 0 : long_names(i) = '********'
457 0 : ierr = nuclide_not_found
458 : end select
459 : end do
460 0 : end subroutine generate_long_nuclide_names
461 :
462 :
463 : ! nuclide information comes from the chem_isos tables
464 : ! the storage container for the data is called 'chem_isos'
465 : ! it has name, A, Z, N, spin, and B for each nuclide
466 : ! use the function chem_get_iso_id to find the index given the name.
467 41257 : integer function chem_get_iso_id(cname)
468 0 : use chem_def, only: get_nuclide_index
469 : character(len=*), intent(in) :: cname
470 41257 : chem_get_iso_id = get_nuclide_index(cname)
471 41257 : end function chem_get_iso_id
472 :
473 :
474 0 : integer function lookup_ZN(Z,N)
475 : integer, intent(in) :: Z, N
476 0 : lookup_ZN = lookup_ZN_isomeric_state(Z,N,0)
477 41257 : end function lookup_ZN
478 :
479 :
480 0 : integer function lookup_ZN_isomeric_state(Z,N,isomeric_state)
481 : use chem_def, only: chem_isos, num_chem_isos
482 : integer, intent(in) :: Z, N, isomeric_state
483 : integer :: cid, i
484 0 : iso_loop: do cid = 1, num_chem_isos
485 0 : if (chem_isos% Z(cid) == Z .and. chem_isos% N(cid) == N) then
486 0 : if (chem_isos% isomeric_state(cid) == isomeric_state) then
487 0 : lookup_ZN_isomeric_state = cid
488 0 : return
489 : end if
490 0 : do i = cid+1, num_chem_isos
491 0 : if (chem_isos% Z(i) /= Z .or. chem_isos% N(i) /= N) exit iso_loop
492 0 : if (chem_isos% isomeric_state(i) == isomeric_state) then
493 0 : lookup_ZN_isomeric_state = i
494 0 : return
495 : end if
496 : end do
497 : end if
498 : end do iso_loop
499 0 : lookup_ZN_isomeric_state = 0 ! indicating failure
500 0 : end function lookup_ZN_isomeric_state
501 :
502 :
503 14 : integer function rates_category_id(cname)
504 0 : use chem_def, only: category_names_dict
505 : use utils_lib
506 : character (len=*), intent(in) :: cname
507 : ! returns id for the category if there is a matching name
508 : ! returns 0 otherwise.
509 : integer :: ierr, value
510 14 : call integer_dict_lookup(category_names_dict, cname, value, ierr)
511 14 : if (ierr /= 0) value = 0
512 14 : rates_category_id = value
513 14 : end function rates_category_id
514 :
515 :
516 0 : function binding_energy(nuclides, Y) result (B)
517 14 : use chem_def
518 : type(nuclide_data), intent(in) :: nuclides
519 : real(dp), dimension(size(nuclides% binding_energy)), intent(in) :: Y
520 : real(dp) :: B
521 0 : B = dot_product(nuclides% binding_energy, Y)
522 0 : end function binding_energy
523 :
524 : ! returns mass excess in MeV
525 1677414 : function get_mass_excess(nuclides,chem_id) result (mass_excess)
526 0 : use chem_def
527 : type(nuclide_data), intent(in) :: nuclides
528 : integer, intent(in) :: chem_id
529 : real(dp) :: mass_excess
530 : logical :: use_nuclides_mass_excess
531 :
532 838707 : use_nuclides_mass_excess = .false.
533 :
534 : ! These should be identical but can have slight ~ulp difference
535 : ! due to floating point maths
536 : if(use_nuclides_mass_excess)then
537 : mass_excess = nuclides% mass_excess(chem_id)
538 : else
539 : mass_excess = nuclides% Z(chem_id)*del_Mp + nuclides% N(chem_id)*del_Mn -&
540 838707 : nuclides% binding_energy(chem_id)
541 : end if
542 :
543 838707 : end function get_mass_excess
544 :
545 1676766 : function get_Q(nuclides,chem_id) result (q)
546 838707 : use chem_def
547 : type(nuclide_data), intent(in) :: nuclides
548 : integer, intent(in) :: chem_id
549 : real(dp) :: q
550 :
551 : !Minus the mass excess
552 838383 : q=-get_mass_excess(nuclides,chem_id)
553 :
554 838383 : end function get_Q
555 :
556 : ! returns the indx corresponding to Tpart just less than T9
557 : ! T9 is the temperature in units of GK
558 : ! returns a value of 0 or npart if value is less than the minimum or maximum of the partition function
559 : ! temperature array Tpart
560 852157 : function get_partition_fcn_indx(T9) result(indx)
561 838383 : use chem_def, only: Tpart, npart
562 : real(dp), intent(in) :: T9
563 : integer :: indx
564 : integer, parameter :: max_iterations = 8
565 : integer :: low, high, mid, i
566 :
567 852157 : low = 1
568 852157 : high = npart
569 :
570 852157 : if (T9 < Tpart(low)) then
571 852157 : indx = low-1
572 : return
573 : end if
574 : if (T9 > Tpart(high)) then
575 3238428 : indx = high + 1
576 : end if
577 :
578 3238428 : do i = 1, max_iterations
579 3238428 : if (high-low <= 1) then
580 577505 : indx = low
581 577505 : return
582 : end if
583 2660923 : mid = (high+low)/2
584 2660923 : if (T9 < Tpart(mid)) then
585 : high = mid
586 : else
587 1493538 : low = mid
588 : end if
589 : end do
590 : ! should never get here
591 0 : indx = low-1
592 :
593 852157 : end function get_partition_fcn_indx
594 :
595 : ! Given a the chem_id's and abundances for a set of isotopes
596 : ! return the abundances of the stable isotopes where the unstable ones
597 : ! have decayed to the stable versions.
598 : ! Code from Frank Timmes "decay.zip"
599 : ! Note this makes some assumptions, firstly that isotopes can only decay to one
600 : ! output (thus there are no branches), also we assume an infinite timescale
601 : ! for decay. If you need a high precision output I suggest you use a one zone
602 : ! burn model rather than this.
603 0 : subroutine get_stable_mass_frac(chem_id,num_species,abun_in,abun_out)
604 852157 : use chem_def
605 : integer,intent(in),dimension(:) :: chem_id
606 : integer,intent(in) :: num_species
607 : real(dp),dimension(:),intent(in) :: abun_in
608 : real(dp),dimension(:),intent(out) :: abun_out
609 : integer :: i,j,a,z
610 : logical :: found
611 :
612 0 : abun_out(1:solsiz)=0.d0
613 :
614 0 : outer: do i=1,num_species
615 0 : Z=chem_isos%Z(chem_id(i))
616 0 : a=z+chem_isos%n(chem_id(i))
617 0 : found =.false.
618 :
619 : ! special case of be8 to he4
620 0 : if (Z == 4 .and. A==8) then ! be8
621 0 : j = 2 ! he4
622 0 : abun_out(j) = abun_out(j) + 2*abun_in(i)
623 0 : found = .true.
624 : end if
625 :
626 0 : inner: do j=1,solsiz
627 0 : if (a /= iasol(j)) then
628 : cycle inner
629 : end if
630 : if ((jcode(j) == 0) .or. &
631 : (z>=izsol(j) .and. jcode(j)==1) .or. &
632 : (z<=izsol(j) .and. jcode(j)==2) .or. &
633 : (z==izsol(j) .and. jcode(j)==3) .or. &
634 : (Z==17 .and. A==36 .and. izsol(j) == 18 .and. iasol(j) == 36) .or. & ! cl36 -> ar36 special case
635 : (Z==21 .and. A==46 .and. izsol(j) == 22 .and. iasol(j) == 46) .or. & ! sc46 -> ti46 special case
636 : (Z==21 .and. A==48 .and. izsol(j) == 22 .and. iasol(j) == 48) .or. & ! sc48 -> ti48 special case
637 : (Z==25 .and. A==54 .and. izsol(j) == 24 .and. iasol(j) == 54) .or. & ! mn54 -> cr54 special case
638 : (Z==27 .and. A==58 .and. izsol(j) == 26 .and. iasol(j) == 58) .or. & ! co58-> fe58 special case
639 : (Z==29 .and. A==64 .and. izsol(j) == 28 .and. iasol(j) == 64) .or. & ! cu64 -> ni64 special case
640 : (Z==31 .and. A==70 .and. izsol(j) == 32 .and. iasol(j) == 70) .or. & ! ga70 -> ge70 special case
641 : (Z==33 .and. A==74 .and. izsol(j) == 32 .and. iasol(j) == 74) .or. & ! as74 -> ge74 special case
642 : (Z==33 .and. A==76 .and. izsol(j) == 34 .and. iasol(j) == 76) .or. & ! as76 -> se76 special case
643 : (Z==35 .and. A==78 .and. izsol(j) == 34 .and. iasol(j) == 78) .or. & ! br78 -> se78 special case
644 : (Z==35 .and. A==80 .and. izsol(j) == 36 .and. iasol(j) == 80) .or. & ! br80 -> kr80 special case
645 0 : (Z==35 .and. A==82 .and. izsol(j) == 36 .and. iasol(j) == 82) .or. & ! br82 -> kr82 special case
646 : (Z==37 .and. A==84 .and. izsol(j) == 36 .and. iasol(j) == 84) & ! rb84 -> kr84 special case
647 0 : ) then
648 0 : abun_out(j) = abun_out(j) + abun_in(i)
649 0 : found = .true.
650 : end if
651 : end do inner
652 :
653 0 : if(found .eqv. .false.) then
654 0 : write(*,*) "Failed to match isotope ",trim(chem_isos%name(chem_id(i)))
655 : end if
656 :
657 : end do outer
658 :
659 : !Normalise results
660 0 : abun_out(1:solsiz)=abun_out(1:solsiz)/sum(abun_out(1:solsiz))
661 :
662 0 : end subroutine get_stable_mass_frac
663 :
664 :
665 0 : real(dp) function chem_M_div_h(x,z,zfrac_choice) ! Returns [M/H]
666 0 : use chem_def
667 : use utils_lib, only: mesa_error
668 : real(dp), intent(in) :: x ! Hydrogen fraction
669 : real(dp), intent(in) :: z ! metal fraction
670 : integer, intent(in) :: zfrac_choice ! See chem_def, *_zfracs options
671 :
672 : real(dp) :: zsolar,ysolar
673 :
674 0 : zsolar = 0d0
675 0 : ysolar = 0d0
676 0 : select case(zfrac_choice)
677 : case(AG89_zfracs)
678 : zsolar = AG89_zsol
679 0 : ysolar = AG89_ysol
680 : case(GN93_zfracs)
681 0 : zsolar = GN93_zsol
682 0 : ysolar = GN93_ysol
683 : case(GS98_zfracs)
684 0 : zsolar = GS98_zsol
685 0 : ysolar = GS98_ysol
686 : case(L03_zfracs)
687 0 : zsolar = L03_zsol
688 0 : ysolar = L03_ysol
689 : case(AGS05_zfracs)
690 0 : zsolar = AGS05_zsol
691 0 : ysolar = AGS05_ysol
692 : case(AGSS09_zfracs)
693 0 : zsolar = AGSS09_zsol
694 0 : ysolar = AGSS09_ysol
695 : case(L09_zfracs)
696 0 : zsolar = L09_zsol
697 0 : ysolar = L09_ysol
698 : case(A09_Prz_zfracs)
699 0 : zsolar = A09_Prz_zsol
700 0 : ysolar = A09_Prz_ysol
701 : case(MB22_photospheric_zfracs)
702 0 : zsolar = MB22_photospheric_zsol
703 0 : ysolar = MB22_photospheric_ysol
704 : case(AAG21_photospheric_zfracs)
705 0 : zsolar = AAG21_photospheric_zsol
706 0 : ysolar = AAG21_photospheric_ysol
707 : case(Custom_zfracs)
708 0 : call mesa_error(__FILE__,__LINE__,"[M/H] not supported with custom zfracs.")
709 : case default
710 0 : call mesa_error(__FILE__,__LINE__,"Bad zfrac_choice")
711 : end select
712 :
713 0 : chem_M_div_h = log10(z/x)-log10(zsolar/(1.d0-zsolar-ysolar))
714 :
715 0 : end function chem_M_div_h
716 :
717 : end module chem_lib
|
