Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 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 reaclib_input
21 : use rates_def
22 : use utils_lib, only: mesa_error
23 :
24 : implicit none
25 :
26 : type(reaclib_data) :: reaclib
27 : integer :: nreaclib
28 :
29 : contains
30 :
31 :
32 135 : subroutine do_extract_rates(set,nuclides,rates,use_weaklib,ierr)
33 : use reaclib_support
34 : use chem_def, only: nuclide_set, nuclide_data
35 : type(nuclide_set), dimension(:), intent(in) :: set
36 : type(nuclide_data), intent(in), target :: nuclides
37 : logical, intent(in) :: use_weaklib
38 : type(reaction_data), intent(out) :: rates
39 : integer, intent(out) :: ierr
40 : logical, parameter :: dbg = .false.
41 :
42 : if (dbg) write(*,*) 'call extract_rates_from_reaclib', nreaclib
43 3 : call extract_rates_from_reaclib(reaclib,nreaclib,nuclides,rates,set,use_weaklib,ierr)
44 3 : if (failed('extract_rates_from_reaclib')) return
45 :
46 : if (dbg) write(*,*) 'call set_up_network_information'
47 3 : call set_up_network_information(rates)
48 :
49 : if (dbg) write(*,*) 'call assign_weights'
50 3 : call assign_weights(rates)
51 :
52 : if (dbg) write(*,*) 'call compute_rev_ratio'
53 3 : call compute_rev_ratio(rates,nuclides)
54 :
55 3 : if (dbg) write(*,*) 'return from do_extract_rates'
56 :
57 : contains
58 :
59 3 : logical function failed(msg)
60 : character (len=*), intent(in) :: msg
61 3 : if (ierr /= 0) then
62 0 : failed = .true.
63 0 : write(*,*) 'do_extract_rates failed in ' // trim(msg)
64 : else
65 : failed = .false.
66 : end if
67 3 : end function failed
68 :
69 : end subroutine do_extract_rates
70 :
71 :
72 3 : subroutine do_read_reaclib(ierr)
73 : use utils_lib, only: integer_dict_define, integer_dict_create_hash
74 : use math_lib, only: str_to_double
75 : integer, intent(out) :: ierr
76 : integer :: i, j, count, reaclib_unitno, cache_io_unit, ios
77 : ! file and table format
78 : character(len=*), parameter :: line0 = '(i2)'
79 : character(len=*), parameter :: line1 = '(5x,6a5,8x,a4,a1,a1,3x,a12)'
80 : !character(len=*), parameter :: line2 = '(4es13.6)'
81 : !character(len=*), parameter :: line3 = '(3es13.6)'
82 : character(len=iso_name_length) :: species
83 : character(len=256) :: filename, cache_filename, buf
84 : character(len=12) :: Qvalue_str
85 :
86 : logical, parameter :: use_cache = .true.
87 :
88 : include 'formats'
89 :
90 3 : ierr = 0
91 :
92 3 : cache_filename = trim(rates_cache_dir) // '/jina_reaclib.bin'
93 :
94 3 : filename = trim(reaclib_dir) // '/' //trim(reaclib_filename)
95 3 : open(newunit=reaclib_unitno, file=filename, iostat=ierr, status="old", action="read")
96 3 : if ( ierr /= 0 ) then
97 0 : write(*,*) 'Error opening ' // filename
98 0 : stop
99 : end if
100 :
101 : ! allocate the library
102 3 : call allocate_reaclib_data(reaclib,max_nreaclib,ierr)
103 3 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'unable to allocate storage for reaclib')
104 :
105 : if (use_cache) then
106 3 : ios = 0
107 : open(newunit=cache_io_unit,file=trim(cache_filename),action='read', &
108 3 : status='old',iostat=ios,form='unformatted')
109 3 : if (ios == 0) then ! opened it okay
110 2 : call read_reaclib_cache(cache_io_unit,ios)
111 2 : close(cache_io_unit)
112 2 : if (ios == 0) then ! read it okay
113 2 : close(reaclib_unitno)
114 2 : return
115 : end if
116 : end if
117 : end if
118 :
119 : count = 0
120 82174 : do i = 1, max_nreaclib
121 82174 : read(unit=reaclib_unitno, fmt=line0, iostat=ierr) reaclib% chapter(i)
122 82174 : if (ierr /= 0 ) then ! assume end of file
123 1 : ierr = 0; exit
124 : end if
125 0 : read(unit=reaclib_unitno,fmt=line1,iostat=ierr,err=100) &
126 82173 : reaclib% species(1,i), reaclib% species(2,i), reaclib% species(3,i), &
127 82173 : reaclib% species(4,i), reaclib% species(5,i), reaclib% species(6,i), &
128 164346 : reaclib% label(i), reaclib% reaction_flag(i), reaclib% reverse_flag(i), Qvalue_str
129 82173 : call str_to_double(Qvalue_str, reaclib% Qvalue(i), ierr)
130 82173 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__)
131 575211 : do j=1,6
132 493038 : species = adjustl(reaclib% species(j,i))
133 575211 : if (species == 'p') then
134 30973 : reaclib% species(j,i) = 'h1'
135 462065 : else if (species == 'n') then
136 50000 : reaclib% species(j,i) = 'neut'
137 412065 : else if (species == 'd') then
138 72 : reaclib% species(j,i) = 'h2'
139 411993 : else if (species== 't') then
140 52 : reaclib% species(j,i) = 'h3'
141 411941 : else if (species== 'al-6') then
142 37 : reaclib% species(j,i) = 'al26-1'
143 411904 : else if (species== 'al*6') then
144 24 : reaclib% species(j,i) = 'al26-2'
145 : end if
146 : end do
147 82173 : read(unit=reaclib_unitno,fmt='(a)',iostat=ierr,err=100) buf
148 82173 : call str_to_double(buf(1:13), reaclib% coefficients(1,i), ierr)
149 82173 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__)
150 82173 : call str_to_double(buf(14:26), reaclib% coefficients(2,i), ierr)
151 82173 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__)
152 82173 : call str_to_double(buf(27:39), reaclib% coefficients(3,i), ierr)
153 82173 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__)
154 82173 : call str_to_double(buf(40:52), reaclib% coefficients(4,i), ierr)
155 82173 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__)
156 82173 : read(unit=reaclib_unitno,fmt='(a)',iostat=ierr,err=100) buf
157 82173 : call str_to_double(buf(1:13), reaclib% coefficients(5,i), ierr)
158 82173 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__)
159 82173 : call str_to_double(buf(14:26), reaclib% coefficients(6,i), ierr)
160 82173 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__)
161 82173 : call str_to_double(buf(27:39), reaclib% coefficients(7,i), ierr)
162 82173 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__)
163 739557 : count = count + 1
164 : end do
165 :
166 1 : nreaclib = count
167 1 : close(reaclib_unitno)
168 :
169 : if (use_cache) then
170 : open(newunit=cache_io_unit, file=trim(cache_filename), iostat=ios, &
171 1 : action='write', form='unformatted')
172 1 : if (ios == 0) then
173 : !write(*,'(a)') 'write ' // trim(cache_filename)
174 1 : call write_reaclib_cache(cache_io_unit)
175 1 : close(cache_io_unit)
176 : end if
177 : end if
178 :
179 : return
180 3 : 100 call mesa_error(__FILE__,__LINE__,'error in do_read_reaclib')
181 :
182 :
183 : contains
184 :
185 :
186 2 : subroutine read_reaclib_cache(io,ios)
187 : integer, intent(in) :: io
188 : integer, intent(out) :: ios
189 : integer :: n
190 2 : ios = 0
191 2 : read(io,iostat=ios) nreaclib
192 2 : if (ios /= 0) return
193 2 : n = nreaclib
194 2 : read(io,iostat=ios) &
195 8 : reaclib% chapter(1:n), &
196 8 : reaclib% species(1:max_species_per_reaction,1:n), &
197 6 : reaclib% label(1:n), &
198 6 : reaclib% reaction_flag(1:n), &
199 6 : reaclib% reverse_flag(1:n), &
200 8 : reaclib% Qvalue(1:n), &
201 8 : reaclib% coefficients(1:ncoefficients,1:n)
202 : if (ios /= 0) return
203 3 : end subroutine read_reaclib_cache
204 :
205 :
206 1 : subroutine write_reaclib_cache(io)
207 : integer, intent(in) :: io
208 : integer :: n
209 1 : write(io) nreaclib
210 1 : n = nreaclib
211 1 : write(io) &
212 4 : reaclib% chapter(1:n), &
213 4 : reaclib% species(1:max_species_per_reaction,1:n), &
214 3 : reaclib% label(1:n), &
215 3 : reaclib% reaction_flag(1:n), &
216 3 : reaclib% reverse_flag(1:n), &
217 4 : reaclib% Qvalue(1:n), &
218 4 : reaclib% coefficients(1:ncoefficients,1:n)
219 1 : end subroutine write_reaclib_cache
220 :
221 :
222 : end subroutine do_read_reaclib
223 :
224 :
225 138 : subroutine extract_rates_from_reaclib(reaclib,nreaclib,nuclides,rates,set,use_weaklib,ierr)
226 : use chem_def, only: &
227 : nuclide_set, nuclide_not_found, &
228 : ipp, icno, i3alf, i_burn_c, i_burn_n, i_burn_o, i_burn_ne, i_burn_na, i_burn_mg, &
229 : i_burn_si, i_burn_s, i_burn_ar, i_burn_ca, i_burn_ti, i_burn_cr, i_burn_fe, icc, &
230 : ico, ioo, iphoto, iother
231 : use chem_lib, only : get_nuclide_index_in_set, get_Q
232 : use utils_lib, only: integer_dict_define, integer_dict_create_hash, integer_dict_lookup
233 : use reaclib_support, only: get_reaction_handle, get_reverse_reaction_handle
234 : type(reaclib_data), intent(in) :: reaclib
235 : integer, intent(in) :: nreaclib
236 : type(nuclide_data), intent(in), target :: nuclides
237 : type(reaction_data), intent(out) :: rates
238 : type(nuclide_set), dimension(:), intent(in) :: set
239 : logical, intent(in) :: use_weaklib
240 : integer, intent(out) :: ierr
241 : type(reaction_data) :: r ! temporary storage
242 : integer :: i,j,l,count,loc_count,nt,indx,cat, &
243 : weaklib_count,chapter,num_in,num_out,num_skip_for_weaklib,num_from_reaclib, &
244 : max_lhs_Z, min_Z
245 : logical :: include_this_rate, already_included_from_weaklib, found_it, is_weak
246 : integer, dimension(max_species_per_reaction) :: pspecies
247 : character(len=max_id_length) :: handle
248 : character(len=iso_name_length) :: name_i, name_j, name1, name2, label
249 : integer, parameter :: max_weaklib_rates = 1000
250 : real(dp) :: Q
251 : integer :: rate_ipp, rate_ipep
252 :
253 : logical, parameter :: dbg = .false.
254 :
255 : include 'formats'
256 :
257 : ierr = 0
258 :
259 : if (dbg) write(*,*) 'call allocate_reaction_data'
260 3 : call allocate_reaction_data(r,max_nreaclib,max_weaklib_rates,ierr)
261 3 : if (ierr /= 0) then
262 0 : print *,'unable to allocate temporary storage for rates'
263 0 : return
264 : end if
265 :
266 3 : count = 0
267 3 : if (use_weaklib) then ! add weaklib rates first
268 : if (dbg) write(*,*) 'add weaklib rates'
269 23571 : do i = 1, nuclides% nnuclides
270 185173779 : do j = 1, nuclides% nnuclides
271 185150208 : if (nuclides% Z(i)+nuclides% N(i) /= nuclides% Z(j)+nuclides% N(j)) cycle
272 618390 : if (abs(nuclides% Z(i) - nuclides% Z(j)) /= 1) cycle
273 : ! check if this one is in weaklib
274 44826 : call get_weaklib_name(i, name_i)
275 44826 : call get_weaklib_name(j, name_j)
276 44826 : indx = do_get_weak_rate_id(name_i, name_j)
277 44826 : if (indx == 0) cycle
278 1394 : count = count + 1
279 1394 : r% weaklib_ids(count) = indx
280 1394 : r% also_in_reaclib(count) = .false.
281 1394 : r% chapter(count) = 1
282 1394 : r% pspecies(1,count) = get_nuclide_index_in_set(name_i,set)
283 1394 : r% pspecies(2,count) = get_nuclide_index_in_set(name_j,set)
284 6970 : r% pspecies(3:max_species_per_reaction,count) = 0
285 11152 : r% coefficients(:,count) = 0
286 1394 : r% coefficients(1,count) = -99
287 1394 : r% reaction_flag(count) = ''
288 : call get_reaction_handle( &
289 : 1, 1, r% pspecies(:,count), nuclides, &
290 1394 : r% reaction_flag(count), r% reaction_handle(count))
291 1394 : r% reverse_handle(count) = ''
292 1394 : r% Q(count) = 0
293 1394 : r% Qneu(count) = 0
294 1394 : r% reaction_flag(count) = 'w'
295 : if (.false.) write(*,*) 'found weaklib rate for ' // name_i &
296 185172382 : // name_j // ' ' // trim(r% reaction_handle(count)), count
297 : end do
298 : end do
299 : end if
300 :
301 3 : weaklib_count = count
302 3 : num_skip_for_weaklib = 0
303 3 : num_from_reaclib = 0
304 3 : loc_count = 0
305 3 : rate_ipp = 0; rate_ipep = 0
306 :
307 : if (dbg) write(*,*) 'loop_over_rates'
308 246522 : loop_over_rates: do i = 1,nreaclib
309 246519 : include_this_rate = .true.
310 246519 : pspecies(:) = 0
311 246519 : chapter = reaclib% chapter(i)
312 246519 : num_in = Nin(chapter)
313 246519 : num_out = Nout(chapter)
314 246519 : nt = num_in+num_out
315 1069236 : loop_over_nuclides: do j = 1, nt
316 825828 : l = get_nuclide_index_in_set(reaclib% species(j,i),set)
317 1069236 : if (l == nuclide_not_found) then
318 : include_this_rate = .false.
319 : exit loop_over_nuclides
320 : else
321 822717 : pspecies(j) = l
322 : end if
323 :
324 : end do loop_over_nuclides
325 :
326 246519 : is_weak = (use_weaklib .and. num_in == 1 .and. num_out == 1)
327 :
328 : ! only include forward rates
329 : ! Define the reverse rate as being the endothermic reaction, always
330 : ! Some photo disintegrations can be exothermic, so dont trust what REACLIB calls a reverse rate
331 246519 : if(include_this_rate .and. .not. is_weak) then
332 1221255 : if(sum(nuclides% binding_energy(pspecies(1:num_in))) - &
333 633414 : sum(nuclides% binding_energy(pspecies(num_in+1:Nt))) > 0 ) then
334 : cycle loop_over_rates
335 : end if
336 : end if
337 147960 : if (include_this_rate) then
338 144846 : if (use_weaklib .and. num_in == 1 .and. num_out == 1) then
339 21465 : already_included_from_weaklib = .false.
340 21465 : name1 = adjustl(reaclib% species(1,i))
341 21465 : name2 = adjustl(reaclib% species(2,i))
342 21465 : indx = do_get_weak_rate_id(name1, name2)
343 21465 : if (indx /= 0) then
344 : already_included_from_weaklib = .true.
345 : else
346 20779 : indx = do_get_weak_rate_id(name2, name1)
347 20779 : if (indx /= 0) then
348 : already_included_from_weaklib = .true.
349 : end if
350 : end if
351 : if (already_included_from_weaklib) then ! find it and store coefficients
352 686 : found_it = .false.
353 158381 : do j=1,weaklib_count
354 158381 : if (r% weaklib_ids(j) == indx) then
355 5488 : r% coefficients(:,j) = reaclib% coefficients(:,i)
356 686 : r% also_in_reaclib(j) = .true.
357 : found_it = .true.
358 : exit
359 : end if
360 : end do
361 : if (.not. found_it) then
362 0 : write(*,*) 'problem in reaclib_input'
363 0 : write(*,*) 'failed to find', indx
364 0 : call mesa_error(__FILE__,__LINE__)
365 : end if
366 : cycle loop_over_rates
367 : end if
368 : end if
369 144160 : count = count + 1
370 144160 : num_from_reaclib = num_from_reaclib + 1
371 144160 : r% chapter(count) = chapter
372 1009120 : r% pspecies(:,count) = pspecies(:)
373 144160 : r% reaction_flag(count) = reaclib% reaction_flag(i)
374 1153280 : r% coefficients(:,count) = reaclib% coefficients(:,i)
375 : ! mark the ec rates so they'll get an extra factor of Ye*Rho
376 144160 : label = adjustl(reaclib% label(i))
377 144160 : if (label(1:2) == 'ec') r% reaction_flag(count) = 'e'
378 : call get_reaction_handle( &
379 : num_in, num_out, r% pspecies(:,count), nuclides, &
380 144160 : r% reaction_flag(count), handle)
381 144160 : r% reaction_handle(count) = handle
382 144160 : if (r% reaction_flag(count) == 'e' .or. r% reaction_flag(count) == 'w') then
383 43870 : r% reverse_handle(count) = ''
384 : else
385 : call get_reverse_reaction_handle( &
386 100290 : num_in, num_out, r% pspecies(:,count), nuclides, r% reverse_handle(count))
387 : end if
388 :
389 144160 : r% Q(count) = 0
390 621837 : do j=1,num_in+num_out
391 477677 : l = pspecies(j)
392 477677 : Q = get_Q(nuclides,l)
393 621837 : if (j <= num_in) then
394 237139 : r% Q(count) = r% Q(count) - Q
395 : else
396 240538 : r% Q(count) = r% Q(count) + Q
397 : end if
398 : end do
399 :
400 144160 : r% Qneu(count) = 0
401 :
402 144160 : if (handle == 'r_b8_to_he4_he4') then
403 0 : r% Qneu(count) = 0.6735D+01
404 :
405 144160 : else if (handle == 'r_h1_he3_to_he4') then
406 0 : r% Qneu(count) = 9.628D0
407 :
408 144160 : else if (handle == 'r_h1_h1_ec_h2') then
409 3 : rate_ipep = count
410 3 : r% Qneu(count) = 1.445D0
411 :
412 144157 : else if (handle == 'r_h1_h1_wk_h2') then
413 3 : rate_ipp = count
414 3 : r% Qneu(count) = 0.2668D0
415 :
416 144154 : else if (handle == 'r_he3_ec_h3') then
417 3 : r% Qneu(count) = 10D0 ! who knows? who cares?
418 :
419 144151 : else if (adjustl(reaclib% reaction_flag(i)) == 'w') then ! check weak_info list
420 43861 : name1 = reaclib% species(1,i)
421 43861 : if (num_out == 1) then
422 20779 : name2 = reaclib% species(2,i)
423 23082 : else if (num_out == 2) then
424 8967 : name2 = reaclib% species(3,i)
425 : else
426 14115 : name2 = ''
427 : end if
428 43861 : j = do_get_weak_info_list_id(name1, name2)
429 43861 : if (j > 0) r% Qneu(count) = weak_info_list_Qneu(j)
430 : if (.false. .and. num_out == 2) &
431 : write(*,2) trim(name1) // ' ' // trim(name2) // &
432 : ' ' // trim(handle) // ' Qneu', count, r% Qneu(count)
433 : end if
434 :
435 144160 : if (reaclib% reaction_flag(i) == 'w' .and. .false.) then
436 : write(*,2) 'reaclib weak ' // trim(handle) // ' Qneu', count, r% Qneu(count)
437 : end if
438 :
439 : end if
440 : end do loop_over_rates
441 :
442 : if (.false.) then
443 : write(*,2) 'num_skip_for_weaklib', num_skip_for_weaklib
444 : write(*,2) 'weaklib_count', weaklib_count
445 : write(*,2) 'num_from_reaclib', num_from_reaclib
446 : write(*,2) 'reaclib+weaklib', num_from_reaclib + weaklib_count
447 : write(*,2) 'total num reactions', count
448 : call mesa_error(__FILE__,__LINE__,'extract_rates_from_reaclib')
449 : end if
450 :
451 :
452 : ! we can now stuff our temporary file into the output and discard the temporary
453 :
454 : if (dbg) write(*,*) 'call allocate_reaction_data'
455 3 : call allocate_reaction_data(rates,count,weaklib_count,ierr)
456 3 : if (ierr /= 0) then
457 0 : print *,'unable to allocate storage for rates'
458 0 : return
459 : end if
460 :
461 3 : rates% nreactions = count
462 3 : rates% nuclides => nuclides
463 :
464 3 : rates% num_from_weaklib = weaklib_count
465 1397 : do i=1,weaklib_count
466 1394 : rates% weaklib_ids(i) = r% weaklib_ids(i)
467 1397 : rates% also_in_reaclib(i) = r% also_in_reaclib(i)
468 : end do
469 :
470 145557 : do i=1,count
471 145554 : rates% reaction_handle(i) = r% reaction_handle(i)
472 145554 : rates% reverse_handle(i) = r% reverse_handle(i)
473 145554 : rates% chapter(i) = r% chapter(i)
474 145554 : rates% reaction_flag(i) = r% reaction_flag(i)
475 145554 : rates% Q(i) = r% Q(i)
476 145554 : rates% Qneu(i) = r% Qneu(i)
477 1018878 : do j=1,max_species_per_reaction
478 1018878 : rates% pspecies(j,i) = r% pspecies(j,i)
479 : end do
480 1164435 : do j=1,ncoefficients
481 1164432 : rates% coefficients(j,i) = r% coefficients(j,i)
482 : end do
483 : end do
484 :
485 3 : nullify(rates% reaction_dict)
486 3 : nullify(rates% reverse_dict)
487 145557 : do i=1,count
488 145554 : chapter = rates% chapter(i)
489 145554 : num_in = Nin(chapter)
490 145554 : num_out = Nout(chapter)
491 145554 : max_lhs_Z = -1
492 384087 : do j=1,num_in
493 384087 : max_lhs_Z = max(max_lhs_Z, nuclides% Z(rates% pspecies(j,i)))
494 : end do
495 145554 : min_Z = 999999
496 626019 : do j=1,num_in+num_out
497 626019 : min_Z = min(min_Z, nuclides% Z(rates% pspecies(j,i)))
498 : end do
499 145554 : cat = -1
500 : ! NOTE: reaction categories that are used by net are set in rates_initialize
501 145554 : if (chapter == r_one_one) then
502 22173 : if (min_Z > 0) then
503 22161 : if (max_lhs_Z <= 5) then
504 : cat = ipp
505 22104 : else if (max_lhs_Z <= 10) then
506 : cat = icno
507 : end if
508 : end if
509 123381 : else if (chapter == r_two_one .or. chapter == r_two_two) then
510 92832 : if (rates% reaction_handle(i) == 'r_he4_c12_to_o16') then
511 : cat = i_burn_c
512 92832 : else if (rates% reaction_handle(i) == 'r_c12_ag_o16') then
513 : cat = i_burn_c
514 92826 : else if (rates% reaction_handle(i) == 'r_o16_ag_ne20') then
515 : cat = i_burn_o
516 92817 : else if (rates% reaction_handle(i) == 'r_he4_o16_to_ne20') then
517 : cat = i_burn_o
518 92817 : else if (rates% reaction_handle(i) == 'r_c12_c12_to_mg24') then
519 : cat = icc
520 92817 : else if (rates% reaction_handle(i) == 'r_c12_c12_to_he4_ne20') then
521 : cat = icc
522 92814 : else if (rates% reaction_handle(i) == 'r_c12_c12_to_p_na23') then
523 : cat = icc
524 92814 : else if (rates% reaction_handle(i) == 'r_c12_o16_to_si28') then
525 : cat = ico
526 92814 : else if (rates% reaction_handle(i) == 'r_c12_o16_to_he4_mg24') then
527 : cat = ico
528 92811 : else if (rates% reaction_handle(i) == 'r_c12_o16_to_p_al27') then
529 : cat = ico
530 92811 : else if (rates% reaction_handle(i) == 'r_o16_o16_to_s32') then
531 : cat = ioo
532 92811 : else if (rates% reaction_handle(i) == 'r_o16_o16_to_he4_si28') then
533 : cat = ioo
534 92808 : else if (rates% reaction_handle(i) == 'r_o16_o16_to_p_p31') then
535 : cat = ioo
536 92808 : else if (min_Z > 0) then
537 41841 : if (max_lhs_Z <= 5) then
538 : cat = ipp
539 41643 : else if (max_lhs_Z <= 9) then
540 : cat = icno
541 : end if
542 : end if
543 : else if (chapter == r_two_three) then
544 69 : if (min_Z > 0 .and. max_lhs_Z <= 5) then
545 : cat = ipp
546 : end if
547 : else if (chapter == r_three_one) then
548 21 : if (rates% reaction_handle(i) == 'r_he4_he4_he4_to_c12') then
549 : cat = i3alf
550 : end if
551 : else if (chapter == r_three_two) then
552 3 : if (rates% reaction_handle(i) == 'r_h1_h1_he4_to_he3_he3') then
553 : cat = ipp
554 : end if
555 : else if (chapter == r_one_two .or. chapter == r_one_three .or. chapter == r_one_four) then
556 30438 : if (rates% reaction_handle(i) == 'r_b8_to_he4_he4') then
557 : cat = ipp
558 : else
559 : cat = iphoto
560 : end if
561 : end if
562 : if (cat == -1) then
563 114030 : if (max_lhs_Z <= 5) then
564 240 : if (min_Z > 0) then
565 : cat = ipp
566 : else
567 97473 : cat = iother
568 : end if
569 : else if (max_lhs_Z <= 6) then
570 : cat = i_burn_c
571 : else if (max_lhs_Z <= 7) then
572 : cat = i_burn_n
573 : else if (max_lhs_Z <= 8) then
574 : cat = i_burn_o
575 113553 : else if (max_lhs_Z <= 10) then
576 : cat = i_burn_ne
577 112920 : else if (max_lhs_Z <= 11) then
578 : cat = i_burn_na
579 112320 : else if (max_lhs_Z <= 12) then
580 : cat = i_burn_mg
581 111690 : else if (max_lhs_Z <= 14) then
582 : cat = i_burn_si
583 110277 : else if (max_lhs_Z <= 16) then
584 : cat = i_burn_s
585 108819 : else if (max_lhs_Z <= 18) then
586 : cat = i_burn_ar
587 107190 : else if (max_lhs_Z <= 20) then
588 : cat = i_burn_ca
589 105411 : else if (max_lhs_Z <= 22) then
590 : cat = i_burn_ti
591 103521 : else if (max_lhs_Z <= 24) then
592 : cat = i_burn_cr
593 101490 : else if (max_lhs_Z <= 28) then
594 : cat = i_burn_fe
595 : else
596 97473 : cat = iother
597 : end if
598 : end if
599 145554 : rates% category(i) = cat
600 :
601 145554 : call integer_dict_define(rates% reaction_dict, rates% reaction_handle(i), i, ierr)
602 145554 : if (ierr /= 0) then
603 0 : write(*,*) 'FATAL ERROR: extract_rates_from_reaclib failed in integer_dict_define'
604 0 : call mesa_error(__FILE__,__LINE__)
605 : end if
606 291111 : if (len_trim(rates% reverse_handle(i)) > 0) then
607 100290 : call integer_dict_define(rates% reverse_dict, rates% reverse_handle(i), i, ierr)
608 100290 : if (ierr /= 0) then
609 0 : write(*,*) 'FATAL ERROR: extract_rates_from_reaclib failed in integer_dict_define'
610 0 : call mesa_error(__FILE__,__LINE__)
611 : end if
612 : end if
613 : end do
614 :
615 : if (dbg) write(*,*) 'call integer_dict_create_hash reaction_dict'
616 3 : call integer_dict_create_hash(rates% reaction_dict, ierr)
617 3 : if (ierr /= 0) then
618 0 : write(*,*) 'FATAL ERROR: extract_rates_from_reaclib failed in integer_dict_create_hash'
619 0 : call mesa_error(__FILE__,__LINE__)
620 : end if
621 :
622 : if (dbg) write(*,*) 'call integer_dict_create_hash reverse_dict'
623 3 : call integer_dict_create_hash(rates% reverse_dict, ierr)
624 3 : if (ierr /= 0) then
625 0 : write(*,*) 'FATAL ERROR: extract_rates_from_reaclib failed in integer_dict_create_hash'
626 0 : call mesa_error(__FILE__,__LINE__)
627 : end if
628 :
629 : if (dbg) write(*,*) 'call free_reaction_data'
630 3 : call free_reaction_data(r)
631 :
632 : ! don't allow blanks in reaction flag
633 145557 : where (rates% reaction_flag == ' ') rates% reaction_flag = '-'
634 :
635 138 : if (dbg) write(*,*) 'done extract_rates_from_reaclib'
636 :
637 :
638 : contains
639 :
640 :
641 89652 : subroutine get_weaklib_name(i,name)
642 : integer, intent(in) :: i
643 : character (len=iso_name_length), intent(out) :: name
644 89652 : if (nuclides% Z(i) == 0) then
645 12 : name = 'neut'
646 89640 : else if (nuclides% Z(i) == 1 .and. nuclides% N(i) == 0) then
647 12 : name = 'h1'
648 : else
649 89628 : name = nuclides% name(i)
650 : end if
651 3 : end subroutine get_weaklib_name
652 :
653 :
654 : end subroutine extract_rates_from_reaclib
655 :
656 :
657 : ! Fowler, Caughlan, Zimmerman, Annual Review Astro. Astrophys., 1975.12:69-112. eqn (1).
658 0 : real(dp) function neutrino_Q(b1, b2)
659 : real(dp), intent(in) :: b1, b2
660 : real(dp) :: sum, sum2
661 0 : sum = b2 - b1 - 0.782d0 - 1.022d0
662 0 : sum = 1.0d0 + sum/0.511d0
663 0 : sum2 = sum*sum
664 : neutrino_Q = abs(0.5d0 * sum * 0.511d0 * (1.0d0 - 1.0d0/sum2) &
665 0 : * (1.0d0 - 1.0d0/(4.0d0*sum) - 1.0d0/(9.0d0*sum2)))
666 0 : end function neutrino_Q
667 :
668 :
669 :
670 : end module reaclib_input
|
