Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010 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 net_lib
21 : ! library for calculating nuclear reaction rates and energy production
22 : ! the data interface for the library is defined in net_def
23 :
24 : use chem_def
25 : use const_def, only: dp, i8
26 :
27 : implicit none
28 :
29 : contains ! the procedure interface for the library
30 : ! client programs should only call these routines.
31 :
32 :
33 : ! call this routine to initialize the net module.
34 : ! only needs to be done once at start of run.
35 :
36 19 : subroutine net_init(ierr)
37 : use net_def, only : do_net_def_init
38 : use net_initialize, only : init_special_case_reaction_info
39 : integer, intent(out) :: ierr ! 0 means AOK.
40 19 : ierr = 0
41 19 : call do_net_def_init
42 19 : call init_special_case_reaction_info
43 :
44 19 : end subroutine net_init
45 :
46 :
47 3 : subroutine net_shutdown
48 19 : end subroutine net_shutdown
49 :
50 :
51 : ! after net_init has finished, you can allocate a "handle".
52 :
53 19 : integer function alloc_net_handle(ierr)
54 : use net_def, only: do_alloc_net, init_net_handle_data
55 : integer, intent(out) :: ierr
56 19 : alloc_net_handle = do_alloc_net(ierr)
57 : if (ierr /= 0) return
58 : end function alloc_net_handle
59 :
60 15 : subroutine free_net_handle(handle)
61 : ! frees the handle and all associated data
62 : use net_def, only: do_free_net
63 : integer, intent(in) :: handle
64 15 : call do_free_net(handle)
65 15 : end subroutine free_net_handle
66 :
67 : ! if you want to access the Net_General_Info record directly,
68 : ! you'll need a pointer to it.
69 19 : subroutine net_ptr(handle, g, ierr)
70 : use net_def, only: Net_General_Info, get_net_ptr
71 : integer, intent(in) :: handle ! from alloc_net_handle
72 : type (Net_General_Info), pointer :: g
73 : integer, intent(out):: ierr
74 19 : call get_net_ptr(handle, g, ierr)
75 19 : end subroutine net_ptr
76 :
77 :
78 : ! routines for defining the net isos and reactions
79 :
80 :
81 : ! call this before starting to define
82 : ! the set of isotopes and reactions for the net.
83 19 : subroutine net_start_def(handle, ierr)
84 : use net_initialize, only:start_net_def
85 : integer, intent(in) :: handle
86 : integer, intent(out) :: ierr
87 19 : call start_net_def(handle, ierr)
88 19 : end subroutine net_start_def
89 :
90 :
91 : ! call this after you've finished defining
92 : ! the set of isotopes and reactions for the net.
93 19 : subroutine net_finish_def(handle, ierr)
94 19 : use net_initialize, only: finish_net_def
95 : integer, intent(in) :: handle
96 : integer, intent(out) :: ierr
97 19 : call finish_net_def(handle, ierr)
98 19 : end subroutine net_finish_def
99 :
100 : ! note: after net_finish_def returns,
101 : ! you have the option of reordering the isotopes
102 : ! before you set up the full set of tables for the net.
103 : ! use get_chem_id_table_ptr and get_net_iso_table_ptr
104 : ! and change both tables to permute the set of isotopes.
105 : ! the default isotope ordering is by increasing chem_id number.
106 :
107 :
108 : ! read_net_file first tries opening the filename in the current directory.
109 : ! if doesn't find that file, then tries the data_dir from the call on net_init.
110 : ! i.e., looks for <data_dir>/net_data/nets/<filename>
111 : ! check net_data/nets/README for info about net files.
112 19 : subroutine read_net_file(filename, handle, ierr)
113 19 : use net_initialize, only: do_read_net_file
114 : character (len=*), intent(in) :: filename
115 : integer, intent(in) :: handle
116 : integer, intent(out) :: ierr
117 19 : call do_read_net_file(filename, handle, ierr)
118 19 : end subroutine read_net_file
119 :
120 :
121 0 : subroutine net_add_iso(handle, iso_id, ierr)
122 19 : use net_initialize, only:add_net_iso
123 : integer, intent(in) :: handle
124 : integer, intent(in) :: iso_id
125 : integer, intent(out) :: ierr
126 0 : call add_net_iso(handle, iso_id, ierr)
127 0 : end subroutine net_add_iso
128 :
129 :
130 0 : subroutine net_add_isos(handle, num_isos, iso_ids, ierr)
131 0 : use net_initialize, only:add_net_isos
132 : integer, intent(in) :: handle
133 : integer, intent(in) :: num_isos, iso_ids(num_isos)
134 : integer, intent(out) :: ierr
135 0 : call add_net_isos(handle, num_isos, iso_ids, ierr)
136 0 : end subroutine net_add_isos
137 :
138 :
139 0 : subroutine net_remove_iso(handle, iso_id, ierr)
140 0 : use net_initialize, only:remove_net_iso
141 : integer, intent(in) :: handle
142 : integer, intent(in) :: iso_id
143 : integer, intent(out) :: ierr
144 0 : call remove_net_iso(handle, iso_id, ierr)
145 0 : end subroutine net_remove_iso
146 :
147 :
148 0 : subroutine net_remove_isos(handle, num_isos, iso_ids, ierr)
149 0 : use net_initialize, only:remove_net_isos
150 : integer, intent(in) :: handle
151 : integer, intent(in) :: num_isos, iso_ids(num_isos)
152 : integer, intent(out) :: ierr
153 0 : call remove_net_isos(handle, num_isos, iso_ids, ierr)
154 0 : end subroutine net_remove_isos
155 :
156 :
157 0 : subroutine net_add_reaction(handle, reaction_id, ierr)
158 0 : use net_initialize, only:add_net_reaction
159 : integer, intent(in) :: handle
160 : integer, intent(in) :: reaction_id
161 : integer, intent(out) :: ierr
162 0 : call add_net_reaction(handle, reaction_id, ierr)
163 0 : end subroutine net_add_reaction
164 :
165 :
166 0 : subroutine net_add_reactions(handle, num_reactions, reaction_ids, ierr)
167 0 : use net_initialize, only:add_net_reactions
168 : integer, intent(in) :: handle
169 : integer, intent(in) :: num_reactions, reaction_ids(num_reactions)
170 : integer, intent(out) :: ierr
171 0 : call add_net_reactions(handle, num_reactions, reaction_ids, ierr)
172 0 : end subroutine net_add_reactions
173 :
174 :
175 0 : subroutine net_remove_reaction(handle, reaction_id, ierr)
176 0 : use net_initialize, only:remove_net_reaction
177 : integer, intent(in) :: handle
178 : integer, intent(in) :: reaction_id
179 : integer, intent(out) :: ierr
180 0 : call remove_net_reaction(handle, reaction_id, ierr)
181 0 : end subroutine net_remove_reaction
182 :
183 :
184 0 : subroutine net_remove_reactions(handle, num_reactions, reaction_ids, ierr)
185 0 : use net_initialize, only:remove_net_reactions
186 : integer, intent(in) :: handle
187 : integer, intent(in) :: num_reactions, reaction_ids(num_reactions)
188 : integer, intent(out) :: ierr
189 0 : call remove_net_reactions(handle, num_reactions, reaction_ids, ierr)
190 0 : end subroutine net_remove_reactions
191 :
192 :
193 0 : subroutine show_net_reactions(handle, iounit, ierr)
194 0 : use net_def
195 : use rates_def, only: reaction_Name
196 : integer, intent(in) :: handle
197 : integer, intent(in) :: iounit
198 : integer, intent(out) :: ierr
199 : integer :: i, id
200 : type (Net_General_Info), pointer :: g
201 : ierr = 0
202 0 : call get_net_ptr(handle, g, ierr)
203 0 : if (ierr /= 0) return
204 0 : do i = 1, g% num_reactions
205 0 : id = g% reaction_id(i)
206 0 : if (id > 0) write(iounit,'(a)') trim(reaction_Name(id))
207 : end do
208 0 : end subroutine show_net_reactions
209 :
210 :
211 0 : subroutine show_net_reactions_and_info(handle, iounit, ierr)
212 0 : use net_def
213 : use rates_def, only: &
214 : reaction_Name, reaction_Info, maxlen_reaction_Info, &
215 : std_reaction_Qs, std_reaction_neuQs, reaction_categories
216 : integer, intent(in) :: handle
217 : integer, intent(in) :: iounit
218 : integer, intent(out) :: ierr
219 : integer :: i, id, weak_id, icat
220 0 : real(dp) :: Q, Qneu
221 : logical :: weaklib_reaction
222 : character (len=maxlen_reaction_Info) :: info
223 : type (Net_General_Info), pointer :: g
224 : ierr = 0
225 0 : call get_net_ptr(handle, g, ierr)
226 0 : if (ierr /= 0) return
227 0 : write(iounit,'(A)')
228 0 : write(iounit,'(4x,a30,3a16,4x,a4,4x,a20)') 'name', 'Qtotal', 'Qneu', 'category', 'info', 'source'
229 0 : do i = 1, g% num_reactions
230 0 : id = g% reaction_id(i)
231 0 : weaklib_reaction = .false.
232 0 : weak_id = g% weak_reaction_index(i)
233 0 : if (weak_id > 0) then
234 0 : if (g% weaklib_ids(weak_id) > 0) weaklib_reaction = .true.
235 : end if
236 0 : if (id > 0) then
237 0 : info = reaction_Info(id)
238 0 : if (reaction_Name(id) == info) info = ''
239 0 : Q = std_reaction_Qs(id)
240 0 : Qneu = std_reaction_neuQs(id)
241 0 : icat = reaction_categories(id)
242 0 : if (weaklib_reaction) then
243 0 : write(iounit,'(i4,a30,16x,a7,9x,4x,a10,4x,a66)') i, &
244 0 : trim(reaction_Name(id)), 'weaklib', trim(category_name(icat)), info
245 0 : else if (Qneu /= 0) then
246 0 : write(iounit,'(i4,a30,2f16.6,4x,a10,4x,a66)') i, trim(reaction_Name(id)), &
247 0 : Q, Qneu, trim(category_name(icat)), info
248 0 : else if (Q /= 0) then
249 0 : write(iounit,'(i4,a30,f16.6,16x,4x,a10,4x,a66)') i, trim(reaction_Name(id)), &
250 0 : Q, trim(category_name(icat)), info
251 : else
252 0 : write(iounit,'(i4,a30,16x,16x,4x,a10,4x,a66)') i, trim(reaction_Name(id)), &
253 0 : trim(category_name(icat)), info
254 : end if
255 : end if
256 : end do
257 0 : write(iounit,'(A)')
258 0 : end subroutine show_net_reactions_and_info
259 :
260 :
261 0 : subroutine show_net_species(handle, iounit, ierr)
262 0 : use net_def
263 : use chem_def
264 : integer, intent(in) :: handle
265 : integer, intent(in) :: iounit
266 : integer, intent(out) :: ierr
267 : integer :: i, id
268 : type (Net_General_Info), pointer :: g
269 : ierr = 0
270 0 : call get_net_ptr(handle, g, ierr)
271 0 : if (ierr /= 0) return
272 0 : do i = 1, g% num_isos
273 0 : id = g% chem_id(i)
274 0 : if (id > 0) write(iounit,'(i4,a10)') i, trim(chem_isos% name(id))
275 : end do
276 0 : end subroutine show_net_species
277 :
278 :
279 0 : subroutine show_net_params(handle, iounit, ierr)
280 0 : use net_def
281 : integer, intent(in) :: handle
282 : integer, intent(in) :: iounit
283 : integer, intent(out) :: ierr
284 : integer :: i, id
285 : type (Net_General_Info), pointer :: g
286 : include 'formats'
287 : ierr = 0
288 0 : call get_net_ptr(handle, g, ierr)
289 0 : if (ierr /= 0) return
290 0 : write(iounit,2) 'logTcut_lo =', g% logTcut_lo
291 0 : write(iounit,2) 'logTcut_lim =', g% logTcut_lim
292 : end subroutine show_net_params
293 :
294 1 : subroutine net_set_fe56ec_fake_factor( &
295 : handle, fe56ec_fake_factor, min_T_for_fe56ec_fake_factor, ierr)
296 : use net_def, only: do_net_set_fe56ec_fake_factor
297 : integer, intent(in) :: handle
298 : real(dp), intent(in) :: fe56ec_fake_factor, min_T_for_fe56ec_fake_factor
299 : integer, intent(out) :: ierr
300 : ierr = 0
301 : call do_net_set_fe56ec_fake_factor( &
302 1 : handle, fe56ec_fake_factor, min_T_for_fe56ec_fake_factor, ierr)
303 : if (ierr /= 0) return
304 : end subroutine net_set_fe56ec_fake_factor
305 :
306 :
307 1 : subroutine net_set_logTcut(handle, logTcut_lo, logTcut_lim, ierr)
308 : use net_def, only: do_net_set_logTcut
309 : integer, intent(in) :: handle
310 : real(dp), intent(in) :: logTcut_lo
311 : real(dp), intent(in) :: logTcut_lim
312 : integer, intent(out) :: ierr
313 : ierr = 0
314 1 : call do_net_set_logTcut(handle, logTcut_lo, logTcut_lim, ierr)
315 : if (ierr /= 0) return
316 : end subroutine net_set_logTcut
317 :
318 :
319 0 : subroutine net_set_eps_nuc_cancel(handle, logT_lo_eps_nuc_cancel, logT_hi_eps_nuc_cancel, ierr)
320 : use net_def, only: do_net_set_eps_nuc_cancel
321 : integer, intent(in) :: handle
322 : real(dp), intent(in) :: logT_lo_eps_nuc_cancel, logT_hi_eps_nuc_cancel
323 : integer, intent(out) :: ierr
324 : ierr = 0
325 0 : call do_net_set_eps_nuc_cancel(handle, logT_lo_eps_nuc_cancel, logT_hi_eps_nuc_cancel, ierr)
326 : if (ierr /= 0) return
327 : end subroutine net_set_eps_nuc_cancel
328 :
329 :
330 : ! call this after you have finished defining the net
331 19 : subroutine net_setup_tables(handle, cache_suffix, ierr)
332 : ! This routine fills in a Net_General_Info structure.
333 : use net_initialize, only : alloc_net_general_info
334 : use rates_lib, only: read_raw_rates_records
335 : use net_def, only: Net_General_Info, get_net_ptr
336 : integer, intent(in) :: handle
337 : character (len=*), intent(in) :: cache_suffix
338 : integer, intent(out) :: ierr
339 : type (Net_General_Info), pointer :: g
340 : ierr = 0
341 19 : call read_raw_rates_records(ierr)
342 19 : if (ierr /= 0) then
343 0 : write(*,*) 'read_raw_rates_records failed in net_setup_tables'
344 0 : return
345 : end if
346 19 : call alloc_net_general_info(handle, cache_suffix, ierr)
347 19 : if (ierr /= 0) then
348 0 : write(*,*) 'alloc_net_general_info failed in net_setup_tables'
349 0 : return
350 : end if
351 19 : end subroutine net_setup_tables
352 :
353 :
354 : ! general info about the net
355 :
356 :
357 1 : integer function net_num_isos(handle, ierr) ! total number in current net
358 19 : use net_def, only: Net_General_Info, get_net_ptr
359 : integer, intent(in) :: handle
360 : integer, intent(out) :: ierr
361 : type (Net_General_Info), pointer :: g
362 1 : net_num_isos = 0
363 1 : call get_net_ptr(handle, g, ierr)
364 1 : if (ierr /= 0) then
365 0 : write(*,*) 'invalid handle for net_num_isos -- did you call alloc_net_handle?'
366 0 : return
367 : end if
368 1 : net_num_isos = g% num_isos
369 1 : end function net_num_isos
370 :
371 1 : integer function net_num_reactions(handle, ierr) ! total number of rates for net
372 : use net_def, only: Net_General_Info, get_net_ptr
373 : integer, intent(in) :: handle
374 : integer, intent(out) :: ierr
375 : type (Net_General_Info), pointer :: g
376 1 : net_num_reactions = 0
377 1 : call get_net_ptr(handle, g, ierr)
378 1 : if (ierr /= 0) then
379 0 : write(*,*) 'invalid handle for net_num_reactions -- did you call alloc_net_handle?'
380 0 : return
381 : end if
382 1 : net_num_reactions = g% num_reactions
383 1 : end function net_num_reactions
384 :
385 32 : subroutine get_chem_id_table(handle, num_isos, chem_id, ierr)
386 : use net_def, only: Net_General_Info, get_net_ptr
387 : integer, intent(in) :: handle, num_isos ! num_isos must be number of isos in current net
388 : integer, intent(out) :: chem_id(num_isos) ! maps net iso number to chem id
389 : ! index from 1 to num_isos in current net
390 : ! value is between 1 and num_chem_isos
391 : integer, intent(out) :: ierr
392 : type (Net_General_Info), pointer :: g
393 : ierr = 0
394 32 : call get_net_ptr(handle, g, ierr)
395 32 : if (ierr /= 0) then
396 0 : write(*,*) 'invalid handle for get_chem_id_table -- did you call alloc_net_handle?'
397 0 : return
398 : end if
399 32 : if (num_isos /= g% num_isos) then
400 0 : ierr = -1
401 0 : return
402 : end if
403 580 : chem_id(1:num_isos) = g% chem_id(1:num_isos)
404 32 : ierr = 0
405 : end subroutine get_chem_id_table
406 :
407 1 : subroutine get_chem_id_table_ptr(handle, chem_id_ptr, ierr)
408 : use net_def, only: Net_General_Info, get_net_ptr
409 : integer, intent(in) :: handle
410 : integer, pointer :: chem_id_ptr(:) ! maps net iso number to chem id
411 : ! index from 1 to num_isos in current net
412 : ! value is between 1 and num_chem_isos
413 : integer, intent(out) :: ierr
414 : type (Net_General_Info), pointer :: g
415 : ierr = 0
416 1 : call get_net_ptr(handle, g, ierr)
417 1 : if (ierr /= 0) then
418 0 : write(*,*) 'invalid handle for get_chem_id_table_ptr -- did you call alloc_net_handle?'
419 0 : return
420 : end if
421 1 : chem_id_ptr => g% chem_id
422 : end subroutine get_chem_id_table_ptr
423 :
424 18 : subroutine get_net_iso_table(handle, net_iso_table, ierr)
425 : use net_def, only: Net_General_Info, get_net_ptr
426 : integer, intent(in) :: handle
427 : integer, intent(out) :: net_iso_table(num_chem_isos) ! maps chem id to net iso number
428 : ! index from 1 to num_chem_isos
429 : ! value is 0 if the iso is not in the current net
430 : ! else is value between 1 and num_isos in current net
431 : integer, intent(out) :: ierr
432 : type (Net_General_Info), pointer :: g
433 : ierr = 0
434 18 : call get_net_ptr(handle, g, ierr)
435 18 : if (ierr /= 0) then
436 0 : write(*,*) 'invalid handle for get_net_iso_table -- did you call alloc_net_handle?'
437 0 : return
438 : end if
439 141426 : net_iso_table(1:num_chem_isos) = g% net_iso(1:num_chem_isos)
440 18 : ierr = 0
441 : end subroutine get_net_iso_table
442 :
443 1 : subroutine get_net_iso_table_ptr(handle, net_iso_ptr, ierr)
444 : use net_def, only: Net_General_Info, get_net_ptr
445 : integer, intent(in) :: handle
446 : integer, pointer :: net_iso_ptr(:) ! maps chem id to net iso number
447 : ! index from 1 to num_chem_isos
448 : ! value is 0 if the iso is not in the current net
449 : ! else is value between 1 and num_isos in current net
450 : integer, intent(out) :: ierr
451 : type (Net_General_Info), pointer :: g
452 : ierr = 0
453 1 : call get_net_ptr(handle, g, ierr)
454 1 : if (ierr /= 0) then
455 0 : write(*,*) 'invalid handle for get_net_iso_table_ptr -- did you call alloc_net_handle?'
456 0 : return
457 : end if
458 1 : net_iso_ptr => g% net_iso
459 1 : ierr = 0
460 : end subroutine get_net_iso_table_ptr
461 :
462 18 : subroutine get_reaction_id_table(handle, num_reactions, reaction_id, ierr)
463 : use net_def, only: Net_General_Info, get_net_ptr
464 : integer, intent(in) :: handle, num_reactions
465 : ! num_reactions must be number of reactions in current net
466 : integer, intent(out) :: reaction_id(num_reactions)
467 : ! maps net reaction number to reaction id
468 : ! index from 1 to num_reactions in current net
469 : ! value is between 1 and num_reactions
470 : integer, intent(out) :: ierr
471 : type (Net_General_Info), pointer :: g
472 : ierr = 0
473 18 : call get_net_ptr(handle, g, ierr)
474 18 : if (ierr /= 0) then
475 0 : write(*,*) 'invalid handle for get_reaction_id_table -- did you call alloc_net_handle?'
476 0 : return
477 : end if
478 18 : if (num_reactions /= g% num_reactions) then
479 0 : ierr = -1
480 0 : return
481 : end if
482 :
483 1298 : reaction_id(1:num_reactions) = g% reaction_id(1:num_reactions)
484 18 : ierr = 0
485 : end subroutine get_reaction_id_table
486 :
487 0 : subroutine get_reaction_id_table_ptr(handle, reaction_id_ptr, ierr)
488 : use net_def, only: Net_General_Info, get_net_ptr
489 : integer, intent(in) :: handle
490 : integer, pointer :: reaction_id_ptr(:) ! maps net reaction number to reaction id
491 : ! index from 1 to num_reactions in current net
492 : ! value is between 1 and rates_reaction_id_max
493 : integer, intent(out) :: ierr
494 : type (Net_General_Info), pointer :: g
495 : ierr = 0
496 0 : call get_net_ptr(handle, g, ierr)
497 0 : if (ierr /= 0) then
498 0 : write(*,*) 'invalid handle for get_reaction_id_table_ptr -- did you call alloc_net_handle?'
499 0 : return
500 : end if
501 0 : reaction_id_ptr => g% reaction_id
502 : end subroutine get_reaction_id_table_ptr
503 :
504 28 : subroutine get_net_reaction_table(handle, net_reaction_table, ierr)
505 : use net_def, only: Net_General_Info, get_net_ptr
506 : use rates_def, only: rates_reaction_id_max
507 : integer, intent(in) :: handle
508 : integer, intent(out) :: net_reaction_table(rates_reaction_id_max)
509 : ! maps reaction id to net reaction number
510 : ! index from 1 to rates_reaction_id_max
511 : ! value is 0 if the reaction is not in the current net
512 : ! else is value between 1 and rates_reaction_id_max in current net
513 : integer, intent(out) :: ierr
514 : type (Net_General_Info), pointer :: g
515 : ierr = 0
516 14 : call get_net_ptr(handle, g, ierr)
517 14 : if (ierr /= 0) then
518 : write(*,*) &
519 0 : 'invalid handle for get_net_reaction_table -- did you call alloc_net_handle?'
520 : return
521 : end if
522 4524 : net_reaction_table(1:rates_reaction_id_max) = g% net_reaction(1:rates_reaction_id_max)
523 14 : ierr = 0
524 14 : end subroutine get_net_reaction_table
525 :
526 4 : subroutine get_net_reaction_table_ptr(handle, net_reaction_ptr, ierr)
527 14 : use net_def, only: Net_General_Info, get_net_ptr
528 : integer, intent(in) :: handle
529 : integer, pointer :: net_reaction_ptr(:)
530 : ! maps reaction id to net reaction number
531 : ! index from 1 to num_reactions
532 : ! value is 0 if the reaction is not in the current net
533 : ! else is value between 1 and num_reactions in current net
534 : integer, intent(out) :: ierr
535 : type (Net_General_Info), pointer :: g
536 : ierr = 0
537 4 : call get_net_ptr(handle, g, ierr)
538 4 : if (ierr /= 0) then
539 : write(*,*) &
540 0 : 'invalid handle for get_net_reaction_table_ptr -- did you call alloc_net_handle?'
541 0 : return
542 : end if
543 4 : net_reaction_ptr => g% net_reaction
544 4 : ierr = 0
545 : end subroutine get_net_reaction_table_ptr
546 :
547 :
548 : ! net evaluation routines
549 :
550 80008 : subroutine net_get( &
551 : handle, just_dxdt, n, num_isos, num_reactions, &
552 80008 : x, temp, log10temp, rho, log10rho, &
553 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
554 : rate_factors, weak_rate_factor, &
555 : reaction_Qs, reaction_neuQs, &
556 80008 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
557 80008 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
558 : screening_mode, &
559 80008 : eps_nuc_categories, eps_neu_total, &
560 : ierr)
561 : use chem_def, only: num_categories
562 : use net_eval, only: eval_net
563 : use net_def, only: Net_General_Info, Net_Info, get_net_ptr
564 :
565 : use rates_def, only: num_rvs
566 :
567 : ! provide T or logT or both (the code needs both, so pass 'em if you've got 'em!)
568 : ! same for Rho and logRho
569 :
570 : integer, intent(in) :: handle
571 : logical, intent(in) :: just_dxdt
572 : type (Net_Info) :: n
573 : integer, intent(in) :: num_isos
574 : integer, intent(in) :: num_reactions
575 : real(dp), intent(in) :: x(:) ! (num_isos)
576 : real(dp), intent(in) :: temp, log10temp ! log10 of temp
577 : real(dp), intent(in) :: rho, log10rho ! log10 of rho
578 : real(dp), intent(in) :: abar ! mean number of nucleons per nucleus
579 : real(dp), intent(in) :: zbar ! mean charge per nucleus
580 : real(dp), intent(in) :: z2bar ! mean charge squared per nucleus
581 : real(dp), intent(in) :: ye
582 : ! mean number free electrons per nucleon, assuming complete ionization
583 : ! d_dxdt_dx(i, j) is d_dxdt(i)_dx(j),
584 : ! i.e., partial derivative of rate for i'th isotope wrt j'th isotope abundance
585 : real(dp), intent(in) :: eta, d_eta_dlnT, d_eta_dlnRho ! electron degeneracy from eos.
586 : ! this arg is only used for prot(e-nu)neut and neut(e+nu)prot.
587 : ! if your net doesn't include those, you can safely ignore this arg.
588 : real(dp), intent(in), pointer :: rate_factors(:) ! (num_reactions)
589 : ! when rates are calculated, they are multiplied by the
590 : ! corresponding values in this array.
591 : ! rate_factors array is indexed by reaction number.
592 : ! use net_reaction_table to map reaction id to reaction number.
593 : real(dp), intent(in) :: weak_rate_factor
594 : real(dp), pointer, intent(in) :: reaction_Qs(:) ! (rates_reaction_id_max)
595 : real(dp), pointer, intent(in) :: reaction_neuQs(:) ! (rates_reaction_id_max)
596 :
597 : real(dp), intent(out) :: eps_nuc ! ergs/g/s from burning after including losses from reaction neutrinos
598 : real(dp), intent(out) :: d_eps_nuc_dT
599 : real(dp), intent(out) :: d_eps_nuc_dRho
600 : real(dp), intent(inout) :: d_eps_nuc_dx(:) ! (num_isos)
601 : ! partial derivatives wrt mass fractions
602 :
603 : real(dp), intent(inout) :: dxdt(:) ! (num_isos)
604 : ! rate of change of mass fractions caused by nuclear reactions
605 : real(dp), intent(inout) :: d_dxdt_dRho(:) ! (num_isos)
606 : real(dp), intent(inout) :: d_dxdt_dT(:) ! (num_isos)
607 : real(dp), intent(inout) :: d_dxdt_dx(:,:) ! (num_isos, num_isos)
608 : ! partial derivatives of rates wrt mass fractions
609 :
610 : real(dp), intent(inout) :: eps_nuc_categories(:) ! (num_categories)
611 : ! eps_nuc subtotals for each reaction category
612 :
613 : real(dp), intent(out) :: eps_neu_total ! ergs/g/s neutrinos from weak reactions
614 :
615 : integer, intent(in) :: screening_mode
616 :
617 : integer, intent(out) :: ierr ! 0 means okay
618 :
619 : integer(i8) :: time0, time1
620 : type (Net_General_Info), pointer :: g
621 80008 : real(dp), pointer, dimension(:) :: actual_Qs, actual_neuQs
622 80008 : logical, pointer :: from_weaklib(:) ! ignore if null
623 : logical, parameter :: symbolic = .false.
624 : logical, parameter :: rates_only = .false.
625 80008 : actual_Qs => null()
626 80008 : actual_neuQs => null()
627 80008 : from_weaklib => null()
628 :
629 : ierr = 0
630 80008 : call get_net_ptr(handle, g, ierr)
631 80008 : if (ierr /= 0) then
632 0 : write(*,*) 'invalid handle for net_get -- did you call alloc_net_handle?'
633 : return
634 : end if
635 :
636 80008 : if (g% doing_timing) then
637 0 : call system_clock(time0)
638 : else
639 : time0 = 0
640 : end if
641 :
642 : call eval_net( &
643 : n, g, rates_only, just_dxdt, num_isos, num_reactions, g% num_wk_reactions, &
644 : x, temp, log10temp, rho, log10rho, &
645 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
646 : rate_factors, weak_rate_factor, &
647 : reaction_Qs, reaction_neuQs, &
648 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
649 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
650 : screening_mode, &
651 : eps_nuc_categories, eps_neu_total, &
652 : actual_Qs, actual_neuQs, from_weaklib, symbolic, &
653 80008 : ierr)
654 80008 : if (g% doing_timing) then
655 0 : call system_clock(time1)
656 0 : g% clock_net_get = g% clock_net_get + (time1 - time0)
657 : end if
658 :
659 160016 : end subroutine net_get
660 :
661 0 : subroutine net_get_rates_only( &
662 : handle, n, num_isos, num_reactions, &
663 0 : x, temp, log10temp, rho, log10rho, &
664 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
665 : rate_factors, weak_rate_factor, &
666 : reaction_Qs, reaction_neuQs, &
667 : screening_mode, &
668 : ierr)
669 80008 : use chem_def, only: num_categories
670 : use net_eval, only: eval_net
671 : use net_def, only: Net_General_Info, Net_Info, get_net_ptr
672 : use rates_def, only: num_rvs
673 :
674 : ! provide T or logT or both (the code needs both, so pass 'em if you've got 'em!)
675 : ! same for Rho and logRho
676 :
677 : integer, intent(in) :: handle
678 : type (Net_Info) :: n
679 : integer, intent(in) :: num_isos
680 : integer, intent(in) :: num_reactions
681 : real(dp), intent(in) :: x(:) ! (num_isos)
682 : real(dp), intent(in) :: temp, log10temp ! log10 of temp
683 : real(dp), intent(in) :: rho, log10rho ! log10 of rho
684 : real(dp), intent(in) :: abar ! mean number of nucleons per nucleus
685 : real(dp), intent(in) :: zbar ! mean charge per nucleus
686 : real(dp), intent(in) :: z2bar ! mean charge squared per nucleus
687 : real(dp), intent(in) :: ye
688 : ! mean number free electrons per nucleon, assuming complete ionization
689 : ! d_dxdt_dx(i, j) is d_dxdt(i)_dx(j),
690 : ! i.e., partial derivative of rate for i'th isotope wrt j'th isotope abundance
691 : real(dp), intent(in) :: eta, d_eta_dlnT, d_eta_dlnRho ! electron degeneracy from eos.
692 : ! this arg is only used for prot(e-nu)neut and neut(e+nu)prot.
693 : ! if your net doesn't include those, you can safely ignore this arg.
694 : real(dp), intent(in), pointer :: rate_factors(:) ! (num_reactions)
695 : ! when rates are calculated, they are multiplied by the
696 : ! corresponding values in this array.
697 : ! rate_factors array is indexed by reaction number.
698 : ! use net_reaction_table to map reaction id to reaction number.
699 : real(dp), intent(in) :: weak_rate_factor
700 : real(dp), pointer, intent(in) :: reaction_Qs(:) ! (rates_reaction_id_max)
701 : real(dp), pointer, intent(in) :: reaction_neuQs(:) ! (rates_reaction_id_max)
702 :
703 : ! rate_raw and rate_screened are described in the declaration of the Net_Info derived type
704 :
705 : integer, intent(in) :: screening_mode
706 :
707 : integer, intent(out) :: ierr ! 0 means okay
708 :
709 : integer(i8) :: time0, time1
710 : type (Net_General_Info), pointer :: g
711 0 : real(dp), pointer, dimension(:) :: actual_Qs, actual_neuQs
712 0 : logical, pointer :: from_weaklib(:) ! ignore if null
713 : logical, parameter :: symbolic = .false.
714 :
715 0 : real(dp) :: eps_nuc, d_eps_nuc_dT, d_eps_nuc_dRho, eps_neu_total
716 :
717 : real(dp), target :: empty_array1(0), empty_array2(0,0)
718 : real(dp), pointer, dimension(:) :: &
719 : d_eps_nuc_dx, dxdt, d_dxdt_dRho, d_dxdt_dT, eps_nuc_categories
720 : real(dp), pointer, dimension(:,:) :: d_dxdt_dx
721 : logical, parameter :: rates_only = .true.
722 : logical, parameter :: just_dxdt = .true.
723 :
724 0 : d_eps_nuc_dx => empty_array1
725 0 : dxdt => empty_array1
726 0 : d_dxdt_dRho => empty_array1
727 0 : d_dxdt_dT => empty_array1
728 :
729 0 : d_dxdt_dx => empty_array2
730 0 : eps_nuc_categories => empty_array1
731 :
732 0 : actual_Qs => null()
733 0 : actual_neuQs => null()
734 0 : from_weaklib => null()
735 :
736 : ierr = 0
737 0 : call get_net_ptr(handle, g, ierr)
738 0 : if (ierr /= 0) then
739 0 : write(*,*) 'invalid handle for net_get -- did you call alloc_net_handle?'
740 : return
741 : end if
742 :
743 0 : if (g% doing_timing) then
744 0 : call system_clock(time0)
745 : else
746 : time0 = 0
747 : end if
748 :
749 : call eval_net( &
750 : n, g, rates_only, just_dxdt, num_isos, num_reactions, g% num_wk_reactions, &
751 : x, temp, log10temp, rho, log10rho, &
752 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
753 : rate_factors, weak_rate_factor, &
754 : reaction_Qs, reaction_neuQs, &
755 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
756 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
757 : screening_mode, &
758 : eps_nuc_categories, eps_neu_total, &
759 : actual_Qs, actual_neuQs, from_weaklib, symbolic, &
760 0 : ierr)
761 0 : if (g% doing_timing) then
762 0 : call system_clock(time1)
763 0 : g% clock_net_get = g% clock_net_get + (time1 - time0)
764 : end if
765 :
766 0 : end subroutine net_get_rates_only
767 :
768 :
769 : ! this sets d_dxdt_dx to 1 in locations where can have a nonzero partial
770 : ! it doesn't set other things such as eps_nuc or rates.
771 : ! takes the same set of args as net_get even though doesn't use them all.
772 0 : subroutine net_get_symbolic_d_dxdt_dx( &
773 : handle, n, num_isos, num_reactions, &
774 0 : x, temp, log10temp, rho, log10rho, &
775 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
776 : rate_factors, weak_rate_factor, &
777 : reaction_Qs, reaction_neuQs, &
778 0 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
779 0 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
780 : screening_mode, &
781 0 : eps_nuc_categories, eps_neu_total, &
782 : ierr)
783 0 : use chem_def, only: num_categories
784 : use net_eval, only: eval_net
785 : use net_def, only: Net_General_Info, Net_Info, get_net_ptr
786 : use rates_def, only: num_rvs
787 :
788 : integer, intent(in) :: handle
789 : type (Net_Info) :: n
790 : integer, intent(in) :: num_isos
791 : integer, intent(in) :: num_reactions
792 : real(dp), intent(in) :: x(:) ! (num_isos)
793 : real(dp), intent(in) :: temp, log10temp ! log10 of temp
794 : real(dp), intent(in) :: rho, log10rho ! log10 of rho
795 : real(dp), intent(in) :: abar ! mean number of nucleons per nucleus
796 : real(dp), intent(in) :: zbar ! mean charge per nucleus
797 : real(dp), intent(in) :: z2bar ! mean charge squared per nucleus
798 : real(dp), intent(in) :: ye
799 : ! mean number free electrons per nucleon, assuming complete ionization
800 : ! d_dxdt_dx(i, j) is d_dxdt(i)_dx(j),
801 : ! i.e., partial derivative of rate for i'th isotope wrt j'th isotope abundance
802 : real(dp), intent(in) :: eta, d_eta_dlnT, d_eta_dlnRho ! electron degeneracy from eos.
803 : ! this arg is only used for prot(e-nu)neut and neut(e+nu)prot.
804 : ! if your net doesn't include those, you can safely ignore this arg.
805 : real(dp), intent(in), pointer :: rate_factors(:) ! (num_reactions)
806 : ! when rates are calculated, they are multiplied by the
807 : ! corresponding values in this array.
808 : ! rate_factors array is indexed by reaction number.
809 : ! use net_reaction_table to map reaction id to reaction number.
810 : real(dp), intent(in) :: weak_rate_factor
811 : real(dp), pointer, intent(in) :: reaction_Qs(:) ! (rates_reaction_id_max)
812 : real(dp), pointer, intent(in) :: reaction_neuQs(:) ! (rates_reaction_id_max)
813 :
814 : real(dp), intent(out) :: eps_nuc ! ergs/g/s from burning after including reaction neutrinos
815 : real(dp), intent(out) :: d_eps_nuc_dT
816 : real(dp), intent(out) :: d_eps_nuc_dRho
817 : real(dp), intent(inout) :: d_eps_nuc_dx(:) ! (num_isos)
818 : ! partial derivatives wrt mass fractions
819 :
820 : real(dp), intent(inout) :: dxdt(:) ! (num_isos)
821 : ! rate of change of mass fractions caused by nuclear reactions
822 : real(dp), intent(inout) :: d_dxdt_dRho(:) ! (num_isos)
823 : real(dp), intent(inout) :: d_dxdt_dT(:) ! (num_isos)
824 : real(dp), intent(inout) :: d_dxdt_dx(:,:) ! (num_isos, num_isos)
825 : ! partial derivatives of rates wrt mass fractions
826 :
827 : real(dp), intent(inout) :: eps_nuc_categories(:) ! (num_categories)
828 : ! eps_nuc subtotals for each reaction category
829 :
830 : real(dp), intent(out) :: eps_neu_total ! ergs/g/s neutrinos from weak reactions
831 :
832 : ! rate_raw and rate_screened are described in the declaration of the Net_Info derived type
833 :
834 : integer, intent(in) :: screening_mode ! Selects which screening mode to use, see rates_def for definition
835 :
836 : integer, intent(out) :: ierr
837 : ! ierr = 0 means AOK
838 : ! ierr = -1 means mass fractions don't add to something very close to 1.0
839 : ! ierr = -2 means neither T nor logT were provided
840 : ! ierr = -3 means neither Rho nor logRho were provided
841 :
842 : type (Net_General_Info), pointer :: g
843 0 : real(dp), pointer, dimension(:) :: actual_Qs, actual_neuQs
844 0 : logical, pointer :: from_weaklib(:) ! ignore if null
845 : logical, parameter :: symbolic = .true.
846 : integer :: num_rates_reduced
847 : real(dp) :: max_old_rate_div_new_rate
848 : logical, parameter :: rates_only = .false.
849 : logical, parameter :: just_dxdt = .false.
850 :
851 0 : actual_Qs => null()
852 0 : actual_neuQs => null()
853 0 : from_weaklib => null()
854 :
855 : ierr = 0
856 0 : call get_net_ptr(handle, g, ierr)
857 0 : if (ierr /= 0) then
858 : write(*,'(a)') &
859 0 : 'invalid handle for net_get_symbolic_d_dxdt_dx -- did you call alloc_net_handle?'
860 : return
861 : end if
862 :
863 : call eval_net( &
864 : n, g, rates_only, just_dxdt, num_isos, num_reactions, g% num_wk_reactions, &
865 : x, temp, log10temp, rho, log10rho, &
866 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
867 : rate_factors, weak_rate_factor, &
868 : reaction_Qs, reaction_neuQs, &
869 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
870 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
871 : screening_mode, &
872 : eps_nuc_categories, eps_neu_total, &
873 : actual_Qs, actual_neuQs, from_weaklib, symbolic, &
874 0 : ierr)
875 :
876 0 : end subroutine net_get_symbolic_d_dxdt_dx
877 :
878 896 : subroutine net_get_with_Qs( &
879 : handle, just_dxdt, n, num_isos, num_reactions, &
880 896 : x, temp, log10temp, rho, log10rho, &
881 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
882 : rate_factors, weak_rate_factor, &
883 : reaction_Qs, reaction_neuQs, &
884 896 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
885 896 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
886 : screening_mode, &
887 896 : eps_nuc_categories, eps_neu_total, &
888 : actual_Qs, actual_neuQs, from_weaklib, &
889 : ierr)
890 0 : use chem_def, only: num_categories
891 : use net_eval, only: eval_net
892 : use net_def, only: Net_General_Info, Net_Info, get_net_ptr
893 : use rates_def, only: num_rvs
894 : integer, intent(in) :: handle
895 : logical, intent(in) :: just_dxdt
896 : type (Net_Info) :: n
897 : integer, intent(in) :: num_isos
898 : integer, intent(in) :: num_reactions
899 : real(dp), intent(in) :: x(:) ! (num_isos)
900 : real(dp), intent(in) :: temp, log10temp ! log10 of temp
901 : real(dp), intent(in) :: rho, log10rho ! log10 of rho
902 : real(dp), intent(in) :: abar ! mean number of nucleons per nucleus
903 : real(dp), intent(in) :: zbar ! mean charge per nucleus
904 : real(dp), intent(in) :: z2bar ! mean charge squared per nucleus
905 : real(dp), intent(in) :: ye
906 : real(dp), intent(in) :: eta, d_eta_dlnT, d_eta_dlnRho ! electron degeneracy from eos.
907 : real(dp), intent(in), pointer :: rate_factors(:) ! (num_reactions)
908 : real(dp), intent(in) :: weak_rate_factor
909 : real(dp), pointer, intent(in) :: reaction_Qs(:) ! (rates_reaction_id_max)
910 : real(dp), pointer, intent(in) :: reaction_neuQs(:) ! (rates_reaction_id_max)
911 : real(dp), intent(out) :: eps_nuc ! ergs/g/s from burning after including reaction neutrinos
912 : real(dp), intent(out) :: d_eps_nuc_dT
913 : real(dp), intent(out) :: d_eps_nuc_dRho
914 : real(dp), intent(inout) :: d_eps_nuc_dx(:) ! (num_isos)
915 : real(dp), intent(inout) :: dxdt(:) ! (num_isos)
916 : real(dp), intent(inout) :: d_dxdt_dRho(:) ! (num_isos)
917 : real(dp), intent(inout) :: d_dxdt_dT(:) ! (num_isos)
918 : real(dp), intent(inout) :: d_dxdt_dx(:,:) ! (num_isos, num_isos)
919 : real(dp), intent(inout) :: eps_nuc_categories(:) ! (num_categories)
920 : real(dp), intent(out) :: eps_neu_total ! ergs/g/s neutrinos from weak reactions
921 : integer, intent(in) :: screening_mode
922 : real(dp), pointer, dimension(:) :: actual_Qs, actual_neuQs ! ignore if null (num_reactions)
923 : logical, pointer :: from_weaklib(:) ! ignore if null
924 : integer, intent(out) :: ierr
925 :
926 : logical, parameter :: rates_only = .false.
927 : logical, parameter :: symbolic = .false.
928 : integer(i8) :: time0, time1
929 : type (Net_General_Info), pointer :: g
930 :
931 : ierr = 0
932 896 : call get_net_ptr(handle, g, ierr)
933 896 : if (ierr /= 0) then
934 0 : write(*,*) 'invalid handle for net_get_with_Qs -- did you call alloc_net_handle?'
935 : return
936 : end if
937 :
938 896 : if (g% doing_timing) then
939 0 : call system_clock(time0)
940 : else
941 : time0 = 0
942 : end if
943 :
944 : call eval_net( &
945 : n, g, rates_only, just_dxdt, num_isos, num_reactions, g% num_wk_reactions, &
946 : x, temp, log10temp, rho, log10rho, &
947 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
948 : rate_factors, weak_rate_factor, &
949 : reaction_Qs, reaction_neuQs, &
950 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
951 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
952 : screening_mode, &
953 : eps_nuc_categories, eps_neu_total, &
954 : actual_Qs, actual_neuQs, from_weaklib, symbolic, &
955 896 : ierr)
956 896 : if (g% doing_timing) then
957 0 : call system_clock(time1)
958 0 : g% clock_net_get = g% clock_net_get + (time1 - time0)
959 : end if
960 :
961 896 : end subroutine net_get_with_Qs
962 :
963 : ! a 1-zone integrator for nets -- for given temperature and density as functions of time
964 2 : subroutine net_1_zone_burn( &
965 2 : net_handle, eos_handle, num_isos, num_reactions, t_start, t_end, starting_x, &
966 : num_times_for_interpolation, times, log10Ts_f1, log10Rhos_f1, etas_f1, &
967 : dxdt_source_term, rate_factors, &
968 : weak_rate_factor, reaction_Qs, reaction_neuQs, &
969 : screening_mode, &
970 : stptry, max_steps, eps, odescal, &
971 : use_pivoting, trace, dbg, burner_finish_substep, &
972 : ! results
973 2 : ending_x, eps_nuc_categories, avg_eps_nuc, eps_neu_total, &
974 : nfcn, njac, nstep, naccpt, nrejct, ierr)
975 896 : use net_burn, only: burn_1_zone
976 : use net_def
977 : use chem_def, only: num_categories
978 :
979 : integer, intent(in) :: net_handle, eos_handle
980 : integer, intent(in) :: num_isos
981 : integer, intent(in) :: num_reactions
982 : real(dp), intent(in) :: t_start, t_end, starting_x(:) ! (num_isos)
983 :
984 : integer, intent(in) :: num_times_for_interpolation
985 : ! ending time is times(num_times); starting time is 0
986 : real(dp), pointer, intent(in) :: times(:) ! (num_times)
987 : real(dp), pointer, intent(in) :: log10Ts_f1(:) ! =(4,numtimes) interpolant for log10T(time)
988 : real(dp), pointer, intent(in) :: log10Rhos_f1(:) ! =(4,numtimes) interpolant for log10Rho(time)
989 : real(dp), pointer, intent(in) :: etas_f1(:) ! =(4,numtimes) interpolant for eta(time)
990 : real(dp), pointer, intent(in) :: dxdt_source_term(:) ! (num_isos) or null if no source term.
991 : real(dp), intent(in), pointer :: rate_factors(:) ! (num_reactions)
992 : real(dp), intent(in) :: weak_rate_factor
993 : real(dp), pointer, intent(in) :: reaction_Qs(:) ! (rates_reaction_id_max)
994 : real(dp), pointer, intent(in) :: reaction_neuQs(:) ! (rates_reaction_id_max)
995 : integer, intent(in) :: screening_mode ! see screen_def
996 : real(dp), intent(in) :: stptry ! try this for 1st step. 0 means try in 1 step.
997 : integer, intent(in) :: max_steps ! maximal number of allowed steps.
998 : real(dp), intent(in) :: eps, odescal ! tolerances. e.g., set both to 1d-6
999 : logical, intent(in) :: use_pivoting ! for matrix solves
1000 : logical, intent(in) :: trace, dbg
1001 : interface
1002 : include 'burner_finish_substep.inc'
1003 : end interface
1004 : real(dp), intent(inout) :: ending_x(:) ! (num_isos)
1005 : real(dp), intent(inout) :: eps_nuc_categories(:) ! (num_categories)
1006 : real(dp), intent(out) :: avg_eps_nuc, eps_neu_total
1007 : integer, intent(out) :: nfcn ! number of function evaluations
1008 : integer, intent(out) :: njac ! number of jacobian evaluations
1009 : integer, intent(out) :: nstep ! number of computed steps
1010 : integer, intent(out) :: naccpt ! number of accepted steps
1011 : integer, intent(out) :: nrejct ! number of rejected steps
1012 : integer, intent(out) :: ierr
1013 :
1014 : call burn_1_zone( &
1015 : net_handle, eos_handle, num_isos, num_reactions, t_start, t_end, starting_x, &
1016 : num_times_for_interpolation, times, log10Ts_f1, log10Rhos_f1, etas_f1, &
1017 : dxdt_source_term, rate_factors, weak_rate_factor, &
1018 : reaction_Qs, reaction_neuQs, screening_mode, &
1019 : stptry, max_steps, eps, odescal, &
1020 : use_pivoting, trace, dbg, burner_finish_substep, &
1021 : ending_x, eps_nuc_categories, avg_eps_nuc, eps_neu_total, &
1022 2 : nfcn, njac, nstep, naccpt, nrejct, ierr)
1023 :
1024 2 : end subroutine net_1_zone_burn
1025 :
1026 :
1027 : ! a 1-zone integrator for nets -- for given density
1028 : ! evolve lnT according to dlnT/dt = eps_nuc/(Cv*T)
1029 0 : subroutine net_1_zone_burn_const_density( &
1030 : net_handle, eos_handle, num_isos, num_reactions, t_start, t_end, &
1031 0 : starting_x, starting_log10T, log10Rho, &
1032 : get_eos_info_for_burn_at_const_density, &
1033 : rate_factors, weak_rate_factor, reaction_Qs, reaction_neuQs, &
1034 : screening_mode, &
1035 : stptry, max_steps, eps, odescal, &
1036 : use_pivoting, trace, dbg, burner_finish_substep, &
1037 : ! results
1038 0 : ending_x, eps_nuc_categories, ending_log10T, &
1039 : avg_eps_nuc, ending_eps_neu_total, &
1040 : nfcn, njac, nstep, naccpt, nrejct, ierr)
1041 2 : use net_burn_const_density, only: burn_const_density_1_zone
1042 : use net_def
1043 : use chem_def, only: num_categories
1044 :
1045 : integer, intent(in) :: net_handle, eos_handle, num_isos, num_reactions
1046 : real(dp), intent(in) :: t_start, t_end, starting_x(:) ! (num_isos)
1047 : real(dp), intent(in) :: starting_log10T, log10Rho
1048 : interface
1049 : include 'burner_const_density_get_eos_info.inc'
1050 : end interface
1051 : real(dp), intent(in), pointer :: rate_factors(:) ! (num_reactions)
1052 : real(dp), intent(in) :: weak_rate_factor
1053 : real(dp), pointer, intent(in) :: reaction_Qs(:) ! (rates_reaction_id_max)
1054 : real(dp), pointer, intent(in) :: reaction_neuQs(:) ! (rates_reaction_id_max)
1055 : integer, intent(in) :: screening_mode ! see screen_def
1056 : real(dp), intent(in) :: stptry ! try this for 1st step. 0 means try in 1 step.
1057 : integer, intent(in) :: max_steps ! maximal number of allowed steps.
1058 : real(dp), intent(in) :: eps, odescal ! tolerances. e.g., set both to 1d-6
1059 : logical, intent(in) :: use_pivoting ! for matrix solves
1060 : logical, intent(in) :: trace, dbg
1061 : interface
1062 : include 'burner_finish_substep.inc'
1063 : end interface
1064 : real(dp), intent(inout) :: ending_x(:) ! (num_isos)
1065 : real(dp), intent(inout) :: eps_nuc_categories(:) ! (num_categories)
1066 : real(dp), intent(out) :: ending_log10T, avg_eps_nuc, ending_eps_neu_total
1067 : integer, intent(out) :: nfcn ! number of function evaluations
1068 : integer, intent(out) :: njac ! number of jacobian evaluations
1069 : integer, intent(out) :: nstep ! number of computed steps
1070 : integer, intent(out) :: naccpt ! number of accepted steps
1071 : integer, intent(out) :: nrejct ! number of rejected steps
1072 : integer, intent(out) :: ierr
1073 :
1074 : call burn_const_density_1_zone( &
1075 : net_handle, eos_handle, num_isos, num_isos+1, num_reactions, t_start, t_end, &
1076 : starting_x, starting_log10T, log10Rho, &
1077 : get_eos_info_for_burn_at_const_density, &
1078 : rate_factors, weak_rate_factor, reaction_Qs, reaction_neuQs, &
1079 : screening_mode, &
1080 : stptry, max_steps, eps, odescal, &
1081 : use_pivoting, trace, dbg, burner_finish_substep, &
1082 : ending_x, eps_nuc_categories, ending_log10T, avg_eps_nuc, ending_eps_neu_total, &
1083 0 : nfcn, njac, nstep, naccpt, nrejct, ierr)
1084 :
1085 0 : end subroutine net_1_zone_burn_const_density
1086 :
1087 : ! evolve T according to dT/dt = eps_nuc/Cp while using given P.
1088 : ! then find new Rho and Cp to match P and new T.
1089 2 : subroutine net_1_zone_burn_const_P( &
1090 : net_handle, eos_handle, num_isos, num_reactions, &
1091 : which_solver, starting_temp, starting_x, clip, &
1092 : ! for interpolating log10P wrt time &
1093 : num_times_for_interpolation, times, log10Ps_f1, &
1094 : ! other args for net_get &
1095 : rate_factors, weak_rate_factor, &
1096 : reaction_Qs, reaction_neuQs, screening_mode, &
1097 : ! args to control the solver &
1098 : h, max_step_size, max_steps, rtol, atol, itol, x_min, x_max, which_decsol, &
1099 : ! results &
1100 : caller_id, solout, iout, &
1101 : ending_x, ending_temp, ending_rho, ending_lnS, initial_rho, initial_lnS, &
1102 : nfcn, njac, nstep, naccpt, nrejct, time_doing_net, time_doing_eos, ierr)
1103 0 : use net_burn_const_P, only: burn_1_zone_const_P
1104 : use chem_def, only: num_categories
1105 : use rates_def, only: num_rvs
1106 :
1107 : integer, intent(in) :: net_handle, eos_handle
1108 : integer, intent(in) :: num_isos
1109 : integer, intent(in) :: num_reactions
1110 : real(dp), pointer, intent(in) :: starting_x(:) ! (num_isos)
1111 : real(dp), intent(in) :: starting_temp
1112 : logical, intent(in) :: clip ! if true, set negative x's to zero during burn.
1113 :
1114 : integer, intent(in) :: which_solver ! as defined in num_def.f
1115 : integer, intent(in) :: num_times_for_interpolation ! ending time is times(num_times); starting time is 0
1116 : real(dp), pointer, intent(in) :: times(:) ! (num_times)
1117 : real(dp), pointer, intent(in) :: log10Ps_f1(:) ! =(4,numtimes) interpolant for log10P(time)
1118 :
1119 : real(dp), intent(in), pointer :: rate_factors(:) ! (num_reactions)
1120 : real(dp), intent(in) :: weak_rate_factor
1121 : real(dp), pointer, intent(in) :: reaction_Qs(:) ! (rates_reaction_id_max)
1122 : real(dp), pointer, intent(in) :: reaction_neuQs(:) ! (rates_reaction_id_max)
1123 : integer, intent(in) :: screening_mode
1124 :
1125 : ! args to control the solver -- see num/public/num_isolve.dek
1126 : real(dp), intent(inout) :: h
1127 : real(dp), intent(in) :: max_step_size ! maximal step size.
1128 : integer, intent(in) :: max_steps ! maximal number of allowed steps.
1129 : ! absolute and relative error tolerances
1130 : real(dp), pointer :: rtol(:) ! relative error tolerance (num_isos)
1131 : real(dp), pointer :: atol(:) ! absolute error tolerance (num_isos)
1132 : integer, intent(in) :: itol ! switch for rtol and atol
1133 : real(dp), intent(in) :: x_min, x_max ! bounds on allowed values
1134 : integer, intent(in) :: which_decsol ! from mtx_def
1135 : integer, intent(in) :: caller_id
1136 : interface ! subroutine called after each successful step
1137 : include "num_solout.dek"
1138 : end interface
1139 : integer, intent(in) :: iout
1140 : real(dp), intent(inout), pointer :: ending_x(:)
1141 : real(dp), intent(out) :: ending_temp, ending_rho, ending_lnS, initial_rho, initial_lnS
1142 : integer, intent(out) :: nfcn ! number of function evaluations
1143 : integer, intent(out) :: njac ! number of jacobian evaluations
1144 : integer, intent(out) :: nstep ! number of computed steps
1145 : integer, intent(out) :: naccpt ! number of accepted steps
1146 : integer, intent(out) :: nrejct ! number of rejected steps
1147 : real(dp), intent(inout) :: time_doing_net
1148 : ! if < 0, then ignore
1149 : ! else on return has input value plus time spent doing eval_net
1150 : real(dp), intent(inout) :: time_doing_eos
1151 : ! if < 0, then ignore
1152 : ! else on return has input value plus time spent doing eos
1153 : integer, intent(out) :: ierr
1154 :
1155 : call burn_1_zone_const_P( &
1156 : net_handle, eos_handle, num_isos, num_reactions, &
1157 : which_solver, starting_temp, starting_x, clip, &
1158 : num_times_for_interpolation, times, log10Ps_f1, &
1159 : rate_factors, weak_rate_factor, &
1160 : reaction_Qs, reaction_neuQs, screening_mode, &
1161 : h, max_step_size, max_steps, rtol, atol, itol, x_min, x_max, which_decsol, &
1162 : caller_id, solout, iout, &
1163 : ending_x, ending_temp, ending_rho, ending_lnS, initial_rho, initial_lnS, &
1164 2 : nfcn, njac, nstep, naccpt, nrejct, time_doing_net, time_doing_eos, ierr)
1165 :
1166 2 : end subroutine net_1_zone_burn_const_P
1167 :
1168 : ! approximate beta decay neutrino energies (in MeV)
1169 : ! Fowler, Caughlan, Zimmerman, Annual Review Astro. Astrophys., 1975.12:69-112. eqn (1).
1170 0 : real(dp) function eval_neutrino_Q(i1, i2)
1171 2 : use net_initialize, only:neutrino_Q
1172 : integer, intent(in) :: i1, i2 ! i1 decays to i2. e.g., i1=in13 and i2=ic13
1173 0 : eval_neutrino_Q = neutrino_Q(i1, i2)
1174 0 : end function eval_neutrino_Q
1175 :
1176 :
1177 : ! for calculating reaction Q
1178 0 : real(dp) function isoB(ci)
1179 0 : use chem_def, only: del_Mp, del_Mn
1180 : integer, intent(in) :: ci
1181 0 : isoB = chem_isos% binding_energy(ci) - chem_isos% Z(ci)*del_Mp - chem_isos% N(ci)*del_Mn
1182 0 : end function isoB
1183 :
1184 :
1185 10 : subroutine clean_up_fractions(nzlo, nzhi, species, nz, xa, max_sum_abs, xsum_tol, ierr)
1186 : ! make sure all fractions are okay and sum to 1.0
1187 0 : use net_eval, only: do_clean_up_fractions
1188 : integer, intent(in) :: nzlo, nzhi, species, nz
1189 : real(dp), intent(inout) :: xa(:,:) ! (species, nz) ! mass fractions
1190 : real(dp), intent(in) :: max_sum_abs
1191 : ! if any k has sum(abs(xa(:,k))) greater than this, set ierr = -1 and return.
1192 : ! else clip each element abundance to the range 0 to 1 and continue.
1193 : real(dp), intent(in) :: xsum_tol
1194 : ! if any sum of abundances is now different from 1 by more than this, set ierr = -1
1195 : ! otherwise, rescale the abundances so that they sum to 1.
1196 : integer, intent(out) :: ierr
1197 10 : call do_clean_up_fractions(nzlo, nzhi, species, nz, xa, max_sum_abs, xsum_tol, ierr)
1198 10 : end subroutine clean_up_fractions
1199 :
1200 :
1201 0 : subroutine clean1(species, xa, max_sum_abs, xsum_tol, ierr)
1202 10 : use net_eval, only: do_clean1
1203 : integer, intent(in) :: species
1204 : real(dp), intent(inout) :: xa(:) ! (species)
1205 : real(dp), intent(in) :: max_sum_abs, xsum_tol
1206 : integer, intent(out) :: ierr
1207 0 : call do_clean1(species, xa, 1, max_sum_abs, xsum_tol, ierr)
1208 0 : end subroutine clean1
1209 :
1210 : end module net_lib
|
