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 5 : 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 5 : call extract_rates_from_reaclib(reaclib,nreaclib,nuclides,rates,set,use_weaklib,ierr)
44 5 : if (failed('extract_rates_from_reaclib')) return
45 :
46 : if (dbg) write(*,*) 'call set_up_network_information'
47 5 : call set_up_network_information(rates)
48 :
49 : if (dbg) write(*,*) 'call assign_weights'
50 5 : call assign_weights(rates)
51 :
52 : if (dbg) write(*,*) 'call compute_rev_ratio'
53 5 : call compute_rev_ratio(rates,nuclides)
54 :
55 5 : if (dbg) write(*,*) 'return from do_extract_rates'
56 :
57 : contains
58 :
59 5 : logical function failed(msg)
60 : character (len=*), intent(in) :: msg
61 5 : 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 5 : end function failed
68 :
69 : end subroutine do_extract_rates
70 :
71 :
72 5 : 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 5 : ierr = 0
91 :
92 5 : cache_filename = trim(rates_cache_dir) // '/jina_reaclib.bin'
93 :
94 5 : filename = trim(reaclib_dir) // '/' //trim(reaclib_filename)
95 5 : open(newunit=reaclib_unitno, file=filename, iostat=ierr, status="old", action="read")
96 5 : if ( ierr /= 0 ) then
97 0 : write(*,*) 'Error opening ' // filename
98 0 : stop
99 : end if
100 :
101 : ! allocate the library
102 5 : call allocate_reaclib_data(reaclib,max_nreaclib,ierr)
103 5 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'unable to allocate storage for reaclib')
104 :
105 : if (use_cache) then
106 5 : ios = 0
107 : open(newunit=cache_io_unit,file=trim(cache_filename),action='read', &
108 5 : status='old',iostat=ios,form='unformatted')
109 5 : if (ios == 0) then ! opened it okay
110 4 : call read_reaclib_cache(cache_io_unit,ios)
111 4 : close(cache_io_unit)
112 4 : if (ios == 0) then ! read it okay
113 4 : close(reaclib_unitno)
114 4 : 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 : 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 5 : 100 call mesa_error(__FILE__,__LINE__,'error in do_read_reaclib')
181 :
182 :
183 : contains
184 :
185 :
186 4 : subroutine read_reaclib_cache(io,ios)
187 : integer, intent(in) :: io
188 : integer, intent(out) :: ios
189 : integer :: n
190 4 : ios = 0
191 4 : read(io,iostat=ios) nreaclib
192 4 : if (ios /= 0) return
193 4 : n = nreaclib
194 : read(io,iostat=ios) &
195 4 : reaclib% chapter(1:n), &
196 4 : reaclib% species(1:max_species_per_reaction,1:n), &
197 4 : reaclib% label(1:n), &
198 4 : reaclib% reaction_flag(1:n), &
199 4 : reaclib% reverse_flag(1:n), &
200 4 : reaclib% Qvalue(1:n), &
201 8 : reaclib% coefficients(1:ncoefficients,1:n)
202 : if (ios /= 0) return
203 5 : 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 : write(io) &
212 1 : reaclib% chapter(1:n), &
213 1 : reaclib% species(1:max_species_per_reaction,1:n), &
214 1 : reaclib% label(1:n), &
215 1 : reaclib% reaction_flag(1:n), &
216 1 : reaclib% reverse_flag(1:n), &
217 1 : reaclib% Qvalue(1:n), &
218 2 : reaclib% coefficients(1:ncoefficients,1:n)
219 1 : end subroutine write_reaclib_cache
220 :
221 :
222 : end subroutine do_read_reaclib
223 :
224 :
225 5 : 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 5 : 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 5 : call allocate_reaction_data(r,max_nreaclib,max_weaklib_rates,ierr)
261 5 : if (ierr /= 0) then
262 0 : print *,'unable to allocate temporary storage for rates'
263 0 : return
264 : end if
265 :
266 5 : count = 0
267 5 : if (use_weaklib) then ! add weaklib rates first
268 : if (dbg) write(*,*) 'add weaklib rates'
269 39285 : do i = 1, nuclides% nnuclides
270 308622965 : do j = 1, nuclides% nnuclides
271 308583680 : if (nuclides% Z(i)+nuclides% N(i) /= nuclides% Z(j)+nuclides% N(j)) cycle
272 1030650 : if (abs(nuclides% Z(i) - nuclides% Z(j)) /= 1) cycle
273 : ! check if this one is in weaklib
274 74710 : call get_weaklib_name(i, name_i)
275 74710 : call get_weaklib_name(j, name_j)
276 74710 : indx = do_get_weak_rate_id(name_i, name_j)
277 74710 : if (indx == 0) cycle
278 2323 : count = count + 1
279 2323 : r% weaklib_ids(count) = indx
280 2323 : r% also_in_reaclib(count) = .false.
281 2323 : r% chapter(count) = 1
282 2323 : r% pspecies(1,count) = get_nuclide_index_in_set(name_i,set)
283 2323 : r% pspecies(2,count) = get_nuclide_index_in_set(name_j,set)
284 11615 : r% pspecies(3:max_species_per_reaction,count) = 0
285 18584 : r% coefficients(:,count) = 0
286 2323 : r% coefficients(1,count) = -99
287 2323 : r% reaction_flag(count) = ''
288 : call get_reaction_handle( &
289 : 1, 1, r% pspecies(:,count), nuclides, &
290 2323 : r% reaction_flag(count), r% reaction_handle(count))
291 2323 : r% reverse_handle(count) = ''
292 2323 : r% Q(count) = 0
293 2323 : r% Qneu(count) = 0
294 2323 : r% reaction_flag(count) = 'w'
295 : if (.false.) write(*,*) 'found weaklib rate for ' // name_i &
296 308620637 : // name_j // ' ' // trim(r% reaction_handle(count)), count
297 : end do
298 : end do
299 : end if
300 :
301 5 : weaklib_count = count
302 5 : num_skip_for_weaklib = 0
303 5 : num_from_reaclib = 0
304 5 : loc_count = 0
305 5 : rate_ipp = 0; rate_ipep = 0
306 :
307 : if (dbg) write(*,*) 'loop_over_rates'
308 410870 : loop_over_rates: do i = 1,nreaclib
309 410865 : include_this_rate = .true.
310 410865 : pspecies(:) = 0
311 410865 : chapter = reaclib% chapter(i)
312 410865 : num_in = Nin(chapter)
313 410865 : num_out = Nout(chapter)
314 410865 : nt = num_in+num_out
315 1782060 : loop_over_nuclides: do j = 1, nt
316 1376380 : l = get_nuclide_index_in_set(reaclib% species(j,i),set)
317 1782060 : if (l == nuclide_not_found) then
318 : include_this_rate = .false.
319 : exit loop_over_nuclides
320 : else
321 1371195 : pspecies(j) = l
322 : end if
323 :
324 : end do loop_over_nuclides
325 :
326 410865 : 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 410865 : if(include_this_rate .and. .not. is_weak) then
332 1665520 : if(sum(nuclides% binding_energy(pspecies(1:num_in))) - &
333 : sum(nuclides% binding_energy(pspecies(num_in+1:Nt))) > 0 ) then
334 : cycle loop_over_rates
335 : end if
336 : end if
337 246600 : if (include_this_rate) then
338 241410 : if (use_weaklib .and. num_in == 1 .and. num_out == 1) then
339 35775 : already_included_from_weaklib = .false.
340 35775 : name1 = adjustl(reaclib% species(1,i))
341 35775 : name2 = adjustl(reaclib% species(2,i))
342 35775 : indx = do_get_weak_rate_id(name1, name2)
343 35775 : if (indx /= 0) then
344 : already_included_from_weaklib = .true.
345 : else
346 34632 : indx = do_get_weak_rate_id(name2, name1)
347 34632 : 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 1143 : found_it = .false.
353 263891 : do j=1,weaklib_count
354 263891 : if (r% weaklib_ids(j) == indx) then
355 9144 : r% coefficients(:,j) = reaclib% coefficients(:,i)
356 1143 : 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 240267 : count = count + 1
370 240267 : num_from_reaclib = num_from_reaclib + 1
371 240267 : r% chapter(count) = chapter
372 1681869 : r% pspecies(:,count) = pspecies(:)
373 240267 : r% reaction_flag(count) = reaclib% reaction_flag(i)
374 1922136 : r% coefficients(:,count) = reaclib% coefficients(:,i)
375 : ! mark the ec rates so they'll get an extra factor of Ye*Rho
376 240267 : label = adjustl(reaclib% label(i))
377 240267 : 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 240267 : r% reaction_flag(count), handle)
381 240267 : r% reaction_handle(count) = handle
382 240267 : if (r% reaction_flag(count) == 'e' .or. r% reaction_flag(count) == 'w') then
383 73117 : r% reverse_handle(count) = ''
384 : else
385 : call get_reverse_reaction_handle( &
386 167150 : num_in, num_out, r% pspecies(:,count), nuclides, r% reverse_handle(count))
387 : end if
388 :
389 240267 : r% Q(count) = 0
390 1036396 : do j=1,num_in+num_out
391 796129 : l = pspecies(j)
392 796129 : Q = get_Q(nuclides,l)
393 1036396 : if (j <= num_in) then
394 395232 : r% Q(count) = r% Q(count) - Q
395 : else
396 400897 : r% Q(count) = r% Q(count) + Q
397 : end if
398 : end do
399 :
400 240267 : r% Qneu(count) = 0
401 :
402 240267 : if (handle == 'r_b8_to_he4_he4') then
403 0 : r% Qneu(count) = 0.6735D+01
404 :
405 240267 : else if (handle == 'r_h1_he3_to_he4') then
406 0 : r% Qneu(count) = 9.628D0
407 :
408 240267 : else if (handle == 'r_h1_h1_ec_h2') then
409 5 : rate_ipep = count
410 5 : r% Qneu(count) = 1.445D0
411 :
412 240262 : else if (handle == 'r_h1_h1_wk_h2') then
413 5 : rate_ipp = count
414 5 : r% Qneu(count) = 0.2668D0
415 :
416 240257 : else if (handle == 'r_he3_ec_h3') then
417 5 : r% Qneu(count) = 10D0 ! who knows? who cares?
418 :
419 240252 : else if (adjustl(reaclib% reaction_flag(i)) == 'w') then ! check weak_info list
420 73102 : name1 = reaclib% species(1,i)
421 73102 : if (num_out == 1) then
422 34632 : name2 = reaclib% species(2,i)
423 38470 : else if (num_out == 2) then
424 14945 : name2 = reaclib% species(3,i)
425 : else
426 23525 : name2 = ''
427 : end if
428 73102 : j = do_get_weak_info_list_id(name1, name2)
429 73102 : 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 : 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 5 : call allocate_reaction_data(rates,count,weaklib_count,ierr)
456 5 : if (ierr /= 0) then
457 0 : print *,'unable to allocate storage for rates'
458 0 : return
459 : end if
460 :
461 5 : rates% nreactions = count
462 5 : rates% nuclides => nuclides
463 :
464 5 : rates% num_from_weaklib = weaklib_count
465 2328 : do i=1,weaklib_count
466 2323 : rates% weaklib_ids(i) = r% weaklib_ids(i)
467 2328 : rates% also_in_reaclib(i) = r% also_in_reaclib(i)
468 : end do
469 :
470 242595 : do i=1,count
471 242590 : rates% reaction_handle(i) = r% reaction_handle(i)
472 242590 : rates% reverse_handle(i) = r% reverse_handle(i)
473 242590 : rates% chapter(i) = r% chapter(i)
474 242590 : rates% reaction_flag(i) = r% reaction_flag(i)
475 242590 : rates% Q(i) = r% Q(i)
476 242590 : rates% Qneu(i) = r% Qneu(i)
477 1698130 : do j=1,max_species_per_reaction
478 1698130 : rates% pspecies(j,i) = r% pspecies(j,i)
479 : end do
480 1940725 : do j=1,ncoefficients
481 1940720 : rates% coefficients(j,i) = r% coefficients(j,i)
482 : end do
483 : end do
484 :
485 5 : nullify(rates% reaction_dict)
486 5 : nullify(rates% reverse_dict)
487 242595 : do i=1,count
488 242590 : chapter = rates% chapter(i)
489 242590 : num_in = Nin(chapter)
490 242590 : num_out = Nout(chapter)
491 242590 : max_lhs_Z = -1
492 640145 : do j=1,num_in
493 640145 : max_lhs_Z = max(max_lhs_Z, nuclides% Z(rates% pspecies(j,i)))
494 : end do
495 242590 : min_Z = 999999
496 1043365 : do j=1,num_in+num_out
497 1043365 : min_Z = min(min_Z, nuclides% Z(rates% pspecies(j,i)))
498 : end do
499 242590 : cat = -1
500 : ! NOTE: reaction categories that are used by net are set in rates_initialize
501 242590 : if (chapter == r_one_one) then
502 36955 : if (min_Z > 0) then
503 36935 : if (max_lhs_Z <= 5) then
504 : cat = ipp
505 36840 : else if (max_lhs_Z <= 10) then
506 : cat = icno
507 : end if
508 : end if
509 205635 : else if (chapter == r_two_one .or. chapter == r_two_two) then
510 154720 : if (rates% reaction_handle(i) == 'r_he4_c12_to_o16') then
511 : cat = i_burn_c
512 154720 : else if (rates% reaction_handle(i) == 'r_c12_ag_o16') then
513 : cat = i_burn_c
514 154710 : else if (rates% reaction_handle(i) == 'r_o16_ag_ne20') then
515 : cat = i_burn_o
516 154695 : else if (rates% reaction_handle(i) == 'r_he4_o16_to_ne20') then
517 : cat = i_burn_o
518 154695 : else if (rates% reaction_handle(i) == 'r_c12_c12_to_mg24') then
519 : cat = icc
520 154695 : else if (rates% reaction_handle(i) == 'r_c12_c12_to_he4_ne20') then
521 : cat = icc
522 154690 : else if (rates% reaction_handle(i) == 'r_c12_c12_to_p_na23') then
523 : cat = icc
524 154690 : else if (rates% reaction_handle(i) == 'r_c12_o16_to_si28') then
525 : cat = ico
526 154690 : else if (rates% reaction_handle(i) == 'r_c12_o16_to_he4_mg24') then
527 : cat = ico
528 154685 : else if (rates% reaction_handle(i) == 'r_c12_o16_to_p_al27') then
529 : cat = ico
530 154685 : else if (rates% reaction_handle(i) == 'r_o16_o16_to_s32') then
531 : cat = ioo
532 154685 : else if (rates% reaction_handle(i) == 'r_o16_o16_to_he4_si28') then
533 : cat = ioo
534 154680 : else if (rates% reaction_handle(i) == 'r_o16_o16_to_p_p31') then
535 : cat = ioo
536 154680 : else if (min_Z > 0) then
537 69735 : if (max_lhs_Z <= 5) then
538 : cat = ipp
539 69405 : 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 115 : 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 35 : 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 5 : 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 50730 : 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 190050 : if (max_lhs_Z <= 5) then
564 400 : if (min_Z > 0) then
565 : cat = ipp
566 : else
567 162455 : 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 189255 : else if (max_lhs_Z <= 10) then
576 : cat = i_burn_ne
577 188200 : else if (max_lhs_Z <= 11) then
578 : cat = i_burn_na
579 187200 : else if (max_lhs_Z <= 12) then
580 : cat = i_burn_mg
581 186150 : else if (max_lhs_Z <= 14) then
582 : cat = i_burn_si
583 183795 : else if (max_lhs_Z <= 16) then
584 : cat = i_burn_s
585 181365 : else if (max_lhs_Z <= 18) then
586 : cat = i_burn_ar
587 178650 : else if (max_lhs_Z <= 20) then
588 : cat = i_burn_ca
589 175685 : else if (max_lhs_Z <= 22) then
590 : cat = i_burn_ti
591 172535 : else if (max_lhs_Z <= 24) then
592 : cat = i_burn_cr
593 169150 : else if (max_lhs_Z <= 28) then
594 : cat = i_burn_fe
595 : else
596 162455 : cat = iother
597 : end if
598 : end if
599 242590 : rates% category(i) = cat
600 :
601 242590 : call integer_dict_define(rates% reaction_dict, rates% reaction_handle(i), i, ierr)
602 242590 : 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 485185 : if (len_trim(rates% reverse_handle(i)) > 0) then
607 167150 : call integer_dict_define(rates% reverse_dict, rates% reverse_handle(i), i, ierr)
608 167150 : 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 5 : call integer_dict_create_hash(rates% reaction_dict, ierr)
617 5 : 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 5 : call integer_dict_create_hash(rates% reverse_dict, ierr)
624 5 : 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 5 : call free_reaction_data(r)
631 :
632 : ! don't allow blanks in reaction flag
633 242595 : where (rates% reaction_flag == ' ') rates% reaction_flag = '-'
634 :
635 5 : if (dbg) write(*,*) 'done extract_rates_from_reaclib'
636 :
637 :
638 : contains
639 :
640 :
641 149420 : subroutine get_weaklib_name(i,name)
642 : integer, intent(in) :: i
643 : character (len=iso_name_length), intent(out) :: name
644 149420 : if (nuclides% Z(i) == 0) then
645 20 : name = 'neut'
646 149400 : else if (nuclides% Z(i) == 1 .and. nuclides% N(i) == 0) then
647 20 : name = 'h1'
648 : else
649 149380 : name = nuclides% name(i)
650 : end if
651 5 : 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 0 : 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
|
