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 rates_support
21 : use const_def, only: dp, use_mesa_temp_cache, missing_value, ln10
22 : use math_lib
23 : use rates_def
24 : use utils_lib, only: mv, switch_str, mesa_error
25 :
26 : implicit none
27 :
28 : integer, parameter :: cache_version = 4
29 :
30 :
31 : contains
32 :
33 89108 : subroutine do_get_raw_rates( &
34 89108 : num_reactions, reaction_id, rattab, rattab_f1, nT8s, &
35 89108 : ye, logtemp_in, btemp, bden, raw_rate_factor, logttab, &
36 89108 : rate_raw, rate_raw_dT, rate_raw_dRho, ierr)
37 : integer, intent(in) :: num_reactions, reaction_id(:), nT8s
38 : real(dp), intent(in) :: &
39 : ye, logtemp_in, btemp, bden, raw_rate_factor(:), &
40 : rattab(:,:), logttab(:)
41 : real(dp), pointer, intent(in) :: rattab_f1(:)
42 : real(dp), intent(inout), dimension(:) :: rate_raw, rate_raw_dT, rate_raw_dRho
43 : integer, intent(out) :: ierr
44 :
45 : integer :: imax, iat0, iat, ir, i, irho
46 : integer, parameter :: mp = 4
47 : real(dp), allocatable :: dtab(:), ddtab(:)
48 : real(dp), pointer :: rattab_f(:,:,:)
49 89108 : real(dp) :: logtemp, fac
50 :
51 : include 'formats'
52 :
53 89108 : ierr = 0
54 :
55 : nullify(rattab_f)
56 89108 : allocate(dtab(num_reactions), ddtab(num_reactions))
57 :
58 89108 : rattab_f(1:4,1:nT8s,1:num_reactions) => rattab_f1(1:4*nT8s*num_reactions)
59 :
60 3416130 : do i = 1,num_reactions
61 :
62 3327022 : ir = reaction_id(i)
63 : !dtab(i) = ye**reaction_ye_rho_exponents(1,ir)
64 6484926 : select case(reaction_ye_rho_exponents(1,ir))
65 : case (0)
66 3157904 : dtab(i) = 1d0
67 : case (1)
68 169118 : dtab(i) = ye
69 : case (2)
70 3327022 : dtab(i) = ye*ye
71 : end select
72 :
73 : !dtab(i) = dtab(i)*bden**reaction_ye_rho_exponents(2,ir)
74 3327022 : irho = reaction_ye_rho_exponents(2,ir)
75 2813792 : select case(irho)
76 : case (1)
77 2813792 : dtab(i) = dtab(i)*bden
78 : case (2)
79 169110 : dtab(i) = dtab(i)*bden*bden
80 : case (3)
81 0 : dtab(i) = dtab(i)*bden*bden*bden
82 : case (4)
83 3327022 : dtab(i) = dtab(i)*bden*bden*bden*bden
84 : end select
85 :
86 3416130 : ddtab(i) = irho*dtab(i)/bden
87 :
88 : end do
89 :
90 89108 : if(warn_rates_for_high_temp .and. logtemp_in >= max_safe_logT_for_rates) then
91 0 : write(*,'(A,F0.6,A,F0.6,A)') "WARNING: evaluating rates with lgT=",logtemp_in," which is above lgT=",&
92 0 : max_safe_logT_for_rates,", rates have been truncated"
93 : end if
94 :
95 :
96 89108 : if(logtemp_in >= max_safe_logT_for_rates) then
97 0 : logtemp = max_safe_logT_for_rates
98 : else
99 89108 : logtemp = logtemp_in
100 : end if
101 :
102 89108 : if (nrattab > 1) then
103 89108 : imax = nrattab
104 89108 : if (logtemp > rattab_thi) then
105 0 : ierr = -1
106 0 : return
107 : end if
108 89108 : iat0 = int((logtemp - rattab_tlo)/rattab_tstp) + 1
109 89108 : iat = max(1, min(iat0 - mp/2 + 1, imax - mp + 1))
110 89108 : call get_rates_from_table(1, num_reactions)
111 : else ! table only has a single temperature
112 0 : do i=1,num_reactions
113 0 : rate_raw(i) = rattab(i,1)*dtab(i)
114 0 : rate_raw_dT(i) = 0
115 0 : rate_raw_dRho(i) = rate_raw(i)*ddtab(i)/dtab(i)
116 : end do
117 : end if
118 :
119 3416130 : do i=1,num_reactions
120 3416130 : if(raw_rate_factor(i)> max_safe_rate_for_any_temp) then
121 0 : write(*,*) "Rate has exceeded any sensible limit for a reaction rate"
122 0 : write(*,*) trim(reaction_Name(reaction_id(i)))
123 0 : write(*,*) raw_rate_factor(i),max_safe_rate_for_any_temp,raw_rate_factor(i)/max_safe_rate_for_any_temp
124 0 : call mesa_error(__FILE__,__LINE__)
125 : end if
126 : end do
127 :
128 3416130 : do i=1,num_reactions
129 3327022 : fac = raw_rate_factor(i)
130 3327022 : rate_raw(i) = rate_raw(i)*fac
131 3327022 : rate_raw_dT(i) = rate_raw_dT(i)*fac
132 3416130 : rate_raw_dRho(i) = rate_raw_dRho(i)*fac
133 : end do
134 :
135 89108 : if(logtemp >= max_safe_logT_for_rates) then
136 0 : rate_raw_dT(1:num_reactions) = 0d0
137 : end if
138 :
139 89108 : nullify(rattab_f)
140 :
141 : contains
142 :
143 89108 : subroutine get_rates_from_table(r1, r2)
144 : integer, intent(in) :: r1, r2
145 :
146 : integer :: i, k
147 89108 : real(dp) :: dt
148 :
149 : include 'formats'
150 :
151 89108 : k = iat+1 ! starting guess for search
152 107934 : do while (logtemp < logttab(k) .and. k > 1)
153 18826 : k = k-1
154 : end do
155 89108 : do while (logtemp > logttab(k+1) .and. k+1 < nrattab)
156 89108 : k = k+1
157 : end do
158 89108 : dt = logtemp - logttab(k)
159 :
160 3416130 : do i = r1,r2
161 :
162 : rate_raw(i) = &
163 : (rattab_f(1,k,i) + dt*(rattab_f(2,k,i) + &
164 : dt*(rattab_f(3,k,i) + dt*rattab_f(4,k,i))) &
165 3327022 : ) * dtab(i)
166 :
167 3327022 : rate_raw_dRho(i) = rate_raw(i) * ddtab(i) / dtab(i)
168 :
169 : rate_raw_dT(i) = &
170 : (rattab_f(2,k,i) + 2*dt*(rattab_f(3,k,i) + &
171 : 1.5d0*dt*rattab_f(4,k,i)) &
172 3416130 : ) * dtab(i) / (btemp * ln10)
173 :
174 : end do
175 :
176 89108 : end subroutine get_rates_from_table
177 :
178 :
179 : end subroutine do_get_raw_rates
180 :
181 :
182 19 : subroutine do_make_rate_tables( &
183 19 : num_reactions, cache_suffix, net_reaction_id, &
184 19 : rattab, rattab_f1, nT8s, ttab, logttab, ierr)
185 : use interp_1d_lib, only: interp_pm, interp_m3q
186 : use interp_1d_def, only: pm_work_size, mp_work_size
187 : use utils_lib
188 : integer, intent(in) :: nT8s, num_reactions, net_reaction_id(:)
189 : character (len=*), intent(in) :: cache_suffix
190 : real(dp) :: rattab(:,:), ttab(:), logttab(:)
191 : real(dp), pointer :: rattab_f1(:)
192 : integer, intent(out) :: ierr
193 :
194 : integer :: i, j, operr, num_to_add_to_cache,thread_num
195 19 : real(dp) :: logT, btemp
196 : real(dp), pointer :: work1(:)=>null(), f1(:)=>null(), rattab_f(:,:,:)=>null()
197 : integer, pointer :: reaction_id(:) =>null()
198 19 : real(dp), allocatable, target :: work(:,:)
199 :
200 : logical :: all_okay, a_okay, all_in_cache
201 :
202 : include 'formats'
203 :
204 19 : ierr = 0
205 :
206 : rattab_f(1:4,1:nrattab,1:num_reactions) => &
207 19 : rattab_f1(1:4*nrattab*num_reactions)
208 :
209 19 : allocate(reaction_id(num_reactions))
210 1330 : reaction_id(:) = net_reaction_id(:)
211 :
212 19 : num_to_add_to_cache = 0
213 19 : if (nrattab == 1) then
214 : all_in_cache = .false.
215 : else
216 19 : all_in_cache = .true.
217 1330 : do i=1, num_reactions
218 1311 : if (read_reaction_from_cache( net_reaction_id, cache_suffix, i, rattab)) then
219 1133 : reaction_id(i) = 0
220 1133 : cycle
221 : end if
222 178 : all_in_cache = .false.
223 178 : num_to_add_to_cache = num_to_add_to_cache + 1
224 : write(*,'(a)') 'create rate data for ' // &
225 197 : trim(reaction_Name(net_reaction_id(i)))
226 : !stop
227 : end do
228 : end if
229 :
230 19 : if (all_in_cache) then
231 :
232 13 : !$OMP PARALLEL DO PRIVATE(i, logT, btemp)
233 : do i=1, nrattab
234 : logT = rattab_tlo + real(i-1,kind=dp)*rattab_tstp
235 : btemp = exp10(logT)
236 : ttab(i) = btemp
237 : logttab(i) = logT
238 : end do
239 : !$OMP END PARALLEL DO
240 :
241 : else
242 :
243 6 : if (num_to_add_to_cache > 20) then
244 3 : write(*,2) 'number not already in cache:', num_to_add_to_cache
245 3 : if (num_to_add_to_cache > 100) write(*,*) 'this will take some time .....'
246 : end if
247 6 : all_okay = .true.
248 : !x$OMP PARALLEL DO PRIVATE(i, operr, logT, btemp, a_okay, j)
249 : ! Disable parralisation as this can cause bugs in the
250 : ! load tables See github bug #360
251 60012 : do i=1, nrattab
252 60006 : logT = rattab_tlo + real(i-1,kind=dp)*rattab_tstp
253 60006 : btemp = exp10(logT)
254 60006 : ttab(i) = btemp
255 60006 : logttab(i) = logT
256 3660366 : do j=1, num_reactions
257 3600360 : if (reaction_id(j) <= 0) cycle
258 3660366 : rattab(j, i) = missing_value ! so can check
259 : end do
260 : operr = 0
261 : !write(*,2) 'logT', i, logT
262 : call get_net_rates_for_tables( &
263 : reaction_id, logT, btemp, num_reactions, &
264 60006 : rattab(1:num_reactions, i), operr)
265 60006 : if (operr /= 0) then
266 0 : ierr = -1
267 0 : cycle
268 : end if
269 60006 : a_okay = .true.
270 3660366 : do j = 1, num_reactions
271 3600360 : if (reaction_id(j) <= 0) cycle
272 1840184 : if (rattab(j, i) == missing_value) then
273 0 : write(*, '(a,i4,2x,a)') 'missing raw rate for ', &
274 0 : j, trim(reaction_Name(reaction_id(j)))
275 0 : a_okay = .false.
276 : end if
277 : end do
278 120018 : if (.not. a_okay) all_okay = .false.
279 : end do
280 : !x$OMP END PARALLEL DO
281 6 : if (.not. all_okay) call mesa_error(__FILE__,__LINE__,'make_rate_tables')
282 6 : if (ierr /= 0) then
283 0 : write(*,*) 'make_rate_tables failed'
284 0 : deallocate(reaction_id)
285 0 : return
286 : end if
287 : end if
288 :
289 19 : if (nrattab > 1) then ! create interpolants
290 38 : allocate(work(nrattab*mp_work_size,utils_OMP_GET_MAX_THREADS()), stat=ierr)
291 19 : call fill_with_NaNs_2D(work)
292 19 : if (ierr /= 0) return
293 19 : !$OMP PARALLEL DO PRIVATE(i,operr,work1,f1,thread_num)
294 : do i=1,num_reactions
295 : thread_num=utils_OMP_GET_THREAD_NUM()+1
296 : work1(1:nrattab*mp_work_size) => work(1:nrattab*mp_work_size,thread_num)
297 : rattab_f(1,1:nrattab,i) = rattab(i,1:nrattab)
298 : f1(1:4*nT8s) => rattab_f1(1+4*nT8s*(i-1):4*nT8s*i)
299 : call interp_m3q(logttab, nrattab, f1, mp_work_size, work1, &
300 : 'rates do_make_rate_tables', operr)
301 : if (operr /= 0) ierr = -1
302 : nullify(f1,work1)
303 : end do
304 : !$OMP END PARALLEL DO
305 19 : deallocate(work)
306 : end if
307 :
308 19 : if (ierr == 0 .and. nrattab > 1 .and. .not. all_in_cache) then
309 366 : do i=1, num_reactions
310 360 : if (reaction_id(i) <= 0) cycle
311 366 : call write_reaction_to_cache(reaction_id, cache_suffix, i, rattab)
312 : end do
313 : end if
314 :
315 19 : deallocate(reaction_id)
316 :
317 19 : end subroutine do_make_rate_tables
318 :
319 :
320 1489 : subroutine reaction_filename(ir, cache_suffix, which, cache_filename, temp_cache_filename, ierr)
321 : integer, intent(in) :: ir, which
322 : character (len=*), intent(in) :: cache_suffix
323 : character (len=*), intent(out) :: cache_filename, temp_cache_filename
324 : integer, intent(out) :: ierr
325 : character (len=64) :: suffix
326 1489 : ierr = 0
327 0 : if (which == 0 .and. len_trim(cache_suffix) > 0) then
328 0 : suffix = cache_suffix
329 : else
330 1489 : if (which < 0) then
331 0 : ierr = -1
332 0 : suffix = '?'
333 1489 : else if (which >= 100) then
334 0 : write(suffix,'(i3)') which
335 1489 : else if (which >= 10) then
336 0 : write(suffix,'(i2)') which
337 : else
338 1489 : write(suffix,'(i1)') which
339 : end if
340 : end if
341 : write(cache_filename,'(a)') &
342 : trim(rates_cache_dir) // '/' // &
343 1489 : trim(reaction_Name(ir)) // '_' // trim(suffix) // '.bin'
344 :
345 : write(temp_cache_filename,'(a)') &
346 : trim(rates_temp_cache_dir) // '/' // &
347 1489 : trim(reaction_Name(ir)) // '_' // trim(suffix) // '.bin'
348 1508 : end subroutine reaction_filename
349 :
350 :
351 :
352 1311 : logical function read_reaction_from_cache(reaction_id, cache_suffix, i, rattab)
353 : integer, intent(in) :: i, reaction_id(:)
354 : character (len=*), intent(in) :: cache_suffix
355 : real(dp),intent(out) :: rattab(:,:)
356 :
357 : integer :: file_version, file_nrattab, file_which
358 1311 : real(dp) :: file_rattab_thi, file_rattab_tlo, file_rattab_tstp
359 : character (len=256) :: cache_filename, temp_cache_filename
360 : integer :: io_unit, ios, ir, which, j, ierr, rir
361 : real(dp), parameter :: tiny = 1d-6
362 : character (len=maxlen_reaction_Name) :: name
363 :
364 : logical, parameter :: show_read_cache = .false.
365 : logical :: reverse_is_table
366 :
367 1311 : ierr = 0
368 1311 : read_reaction_from_cache = .false.
369 1311 : if (.not. rates_use_cache) return
370 :
371 1311 : ir = reaction_id(i)
372 1311 : which = 1
373 :
374 1311 : reverse_is_table = .false.
375 1311 : rir = reverse_reaction_id(ir)
376 1311 : if(rir>0) reverse_is_table = raw_rates_records(reverse_reaction_id(ir))% use_rate_table
377 :
378 :
379 1311 : if (raw_rates_records(ir)% use_rate_table .or. reverse_is_table) then
380 : which = 0
381 : !Dont read a cached version of a users local rate
382 : return
383 : end if
384 :
385 1311 : call reaction_filename(reaction_id(i), cache_suffix, which, cache_filename, temp_cache_filename, ierr)
386 1311 : if (ierr /= 0) then
387 : if (show_read_cache) write(*,*) 'read cache -- bad reaction_filename ' // trim(cache_filename)
388 : return
389 : end if
390 :
391 1311 : ios = 0
392 : open(newunit=io_unit,file=trim(cache_filename),action='read', &
393 1311 : status='old',iostat=ios,form='unformatted')
394 1311 : if (ios /= 0) then
395 : if (show_read_cache) write(*,*) 'read cache failed for open ' // trim(cache_filename)
396 : return
397 : end if
398 :
399 : read(io_unit, iostat=ios) &
400 1133 : name, file_which, file_version, file_nrattab, &
401 2266 : file_rattab_thi, file_rattab_tlo, file_rattab_tstp
402 1133 : if (ios /= 0) then
403 : if (show_read_cache) write(*,*) 'read cache failed for read header ' // trim(cache_filename)
404 0 : close(io_unit)
405 0 : return
406 : end if
407 :
408 1133 : if (name /= reaction_Name(ir)) then
409 : if (show_read_cache) write(*,*) 'read cache failed for name'
410 0 : close(io_unit)
411 0 : return
412 : end if
413 :
414 1133 : if (which /= file_which) then
415 : if (show_read_cache) write(*,*) 'read cache failed for which reaction'
416 0 : close(io_unit)
417 0 : return
418 : end if
419 :
420 1133 : if (cache_version /= file_version) then
421 : if (show_read_cache) write(*,*) 'read cache failed for version'
422 0 : close(io_unit)
423 0 : return
424 : end if
425 :
426 1133 : if (abs(rattab_thi-file_rattab_thi) > tiny) then
427 : if (show_read_cache) write(*,*) 'read cache failed for rattab_thi'
428 0 : close(io_unit)
429 0 : return
430 : end if
431 :
432 1133 : if (abs(rattab_tlo-file_rattab_tlo) > tiny) then
433 : if (show_read_cache) write(*,*) 'read cache failed for rattab_tlo'
434 0 : close(io_unit)
435 0 : return
436 : end if
437 :
438 1133 : if (abs(rattab_tstp-file_rattab_tstp) > tiny) then
439 : if (show_read_cache) write(*,*) 'read cache failed for rattab_tstp'
440 0 : close(io_unit)
441 0 : return
442 : end if
443 :
444 11332266 : do j = 1, nrattab
445 11331133 : read(io_unit, iostat=ios) rattab(i,j)
446 11332266 : if (ios /= 0) then
447 : if (show_read_cache) write(*,*) 'read cache failed for reaction'
448 0 : close(io_unit)
449 0 : return
450 : end if
451 : end do
452 :
453 1133 : close(io_unit)
454 :
455 1133 : read_reaction_from_cache = .true.
456 :
457 1133 : end function read_reaction_from_cache
458 :
459 :
460 :
461 178 : subroutine write_reaction_to_cache(reaction_id, cache_suffix, i, rattab)
462 : integer, intent(in) :: i
463 : character (len=*), intent(in) :: cache_suffix
464 : integer, intent(in) :: reaction_id(:)
465 : real(dp), intent(in) :: rattab(:,:)
466 :
467 : character (len=256) :: cache_filename, temp_cache_filename
468 : integer :: io_unit, ios, ir, which, ierr, j, rir
469 :
470 : logical, parameter :: show_write_cache = .true.
471 : logical :: reverse_is_table
472 :
473 178 : ierr = 0
474 178 : if (.not. rates_use_cache) return
475 :
476 178 : ir = reaction_id(i)
477 178 : which = 1
478 :
479 178 : reverse_is_table = .false.
480 178 : rir = reverse_reaction_id(ir)
481 178 : if(rir>0) reverse_is_table = raw_rates_records(reverse_reaction_id(ir))% use_rate_table
482 :
483 :
484 178 : if (raw_rates_records(ir)% use_rate_table .or. reverse_is_table) which = 0
485 :
486 :
487 : ! Write cache file to temporary storage that is local to the run,
488 : ! then at the end move the file atomically to the final cache location
489 178 : call reaction_filename(reaction_id(i), cache_suffix, which, cache_filename, temp_cache_filename, ierr)
490 178 : if (ierr /= 0) return
491 :
492 178 : ios = 0
493 : open(newunit=io_unit, file=trim(switch_str(temp_cache_filename, cache_filename, use_mesa_temp_cache)),&
494 178 : iostat=ios, action='write', form='unformatted')
495 178 : if (ios /= 0) then
496 0 : if (show_write_cache) write(*,*) 'write_cache failed to open ', trim(temp_cache_filename)
497 0 : return
498 : end if
499 :
500 178 : if (show_write_cache) write(*,'(a)') 'write ' // trim(cache_filename)
501 :
502 : write(io_unit) &
503 178 : reaction_Name(ir), which, cache_version, nrattab, &
504 356 : rattab_thi, rattab_tlo, rattab_tstp
505 :
506 1780356 : do j = 1, nrattab
507 1780356 : write(io_unit) rattab(i,j)
508 : end do
509 :
510 178 : close(io_unit)
511 178 : if(use_mesa_temp_cache) call mv(temp_cache_filename,cache_filename,.true.)
512 :
513 :
514 : end subroutine write_reaction_to_cache
515 :
516 :
517 0 : subroutine do_show_reaction_from_cache(cache_filename, ierr)
518 : character (len=*) :: cache_filename
519 : integer, intent(out) :: ierr
520 :
521 : integer :: version, nrattab, which
522 0 : real(dp) :: rattab_thi, rattab_tlo, rattab_tstp, rate, T8, logT
523 : integer :: ios, j, io_unit
524 : real(dp), parameter :: tiny = 1d-6
525 : character (len=maxlen_reaction_Name) :: name
526 :
527 0 : ierr = 0
528 0 : ios = 0
529 : open(newunit=io_unit,file=trim(cache_filename),action='read', &
530 0 : status='old',iostat=ios,form='unformatted')
531 0 : if (ios /= 0) then
532 0 : write(*,*) 'read cache failed for open ' // trim(cache_filename)
533 0 : return
534 : end if
535 :
536 : read(io_unit, iostat=ios) &
537 0 : name, which, version, nrattab, &
538 0 : rattab_thi, rattab_tlo, rattab_tstp
539 0 : if (ios /= 0) then
540 0 : write(*,*) 'read cache failed for read header ' // trim(cache_filename)
541 0 : close(io_unit)
542 0 : return
543 : end if
544 :
545 0 : write(*,'(a)') '# T8 rate'
546 0 : write(*,'(A)')
547 0 : write(*,*) nrattab
548 :
549 0 : do j = 1, nrattab
550 0 : read(io_unit, iostat=ios) rate
551 0 : if (ios /= 0) then
552 0 : write(*,*) 'read cache failed for reaction data', j
553 0 : close(io_unit)
554 0 : return
555 : end if
556 0 : logT = rattab_tlo + dble(j-1)*rattab_tstp
557 0 : T8 = exp10(logT - 8d0)
558 0 : write(*,'(1pe26.16,3x,1pe26.16e3)') T8, rate
559 : end do
560 0 : write(*,'(A)')
561 :
562 0 : close(io_unit)
563 :
564 : end subroutine do_show_reaction_from_cache
565 :
566 :
567 60006 : subroutine get_net_rates_for_tables( &
568 60006 : reaction_id, logT, btemp, num_reactions, rates, ierr)
569 : use ratelib, only: tfactors
570 : use raw_rates, only: set_raw_rates
571 : use utils_lib, only: is_bad
572 :
573 : real(dp), intent(in) :: logT, btemp
574 : integer, intent(in) :: num_reactions, reaction_id(:)
575 : real(dp), intent(inout) :: rates(:)
576 : integer, intent(out) :: ierr
577 :
578 : integer :: i, ir
579 : type (T_Factors) :: tf
580 :
581 : include 'formats'
582 :
583 60006 : ierr = 0
584 :
585 60006 : call tfactors(tf, logT, btemp)
586 : call set_raw_rates( &
587 60006 : num_reactions, reaction_id, btemp, tf, rates, ierr)
588 60006 : if (ierr /= 0) return
589 :
590 3660366 : do i = 1, num_reactions
591 3600360 : ir = reaction_id(i)
592 3600360 : if (ir <= 0) cycle
593 1840184 : if (is_bad(rates(i))) then
594 0 : write(*,2) trim(reaction_Name(ir)) // ' rates', i, rates(i)
595 0 : call mesa_error(__FILE__,__LINE__,'get_net_rates_for_tables')
596 : end if
597 : end do
598 :
599 60006 : end subroutine get_net_rates_for_tables
600 :
601 :
602 0 : subroutine do_eval_reaclib_21( &
603 0 : ir, temp, den, rate_raw, reverse_rate_raw, ierr)
604 60006 : use raw_rates, only: get_reaclib_rate_and_dlnT
605 : integer, intent(in) :: ir ! reaction_id
606 : real(dp), intent(in) :: temp, den
607 : real(dp), intent(inout) :: rate_raw(:), reverse_rate_raw(:)
608 : integer, intent(out) :: ierr
609 :
610 : real(dp) :: lambda, dlambda_dlnT, rlambda, drlambda_dlnT
611 :
612 : include 'formats'
613 :
614 : ierr = 0
615 : call get_reaclib_rate_and_dlnT( &
616 0 : ir, temp, lambda, dlambda_dlnT, rlambda, drlambda_dlnT, ierr)
617 0 : if (ierr /= 0) return
618 :
619 0 : if (reaction_ye_rho_exponents(2,ir) /= 1) then
620 0 : ierr = -1
621 0 : return
622 : end if
623 :
624 0 : rate_raw(i_rate) = lambda*den
625 0 : rate_raw(i_rate_dT) = dlambda_dlnT*den/temp
626 0 : rate_raw(i_rate_dRho) = lambda
627 :
628 0 : reverse_rate_raw(i_rate) = rlambda
629 0 : reverse_rate_raw(i_rate_dT) = drlambda_dlnT/temp
630 0 : reverse_rate_raw(i_rate_dRho) = 0d0
631 :
632 0 : return
633 :
634 : !$omp critical (rates_eval_reaclib_21)
635 : write(*,1) 'do_eval_reaclib_21 ' // trim(reaction_Name(ir))
636 : write(*,'(A)')
637 : write(*,1) 'den', den
638 : write(*,1) 'temp', temp
639 : write(*,'(A)')
640 : write(*,1) 'lambda', lambda
641 : write(*,1) 'dlambda_dlnT', dlambda_dlnT
642 : write(*,1) 'rate_raw', rate_raw(1:num_rvs)
643 : write(*,'(A)')
644 : write(*,1) 'rlambda', rlambda
645 : write(*,1) 'drlambda_dlnT', drlambda_dlnT
646 : write(*,1) 'reverse_rate_raw', reverse_rate_raw(1:num_rvs)
647 : write(*,'(A)')
648 : !$omp end critical (rates_eval_reaclib_21)
649 :
650 0 : end subroutine do_eval_reaclib_21
651 :
652 :
653 0 : subroutine do_eval_reaclib_22( &
654 0 : ir, temp, den, rate_raw, reverse_rate_raw, ierr)
655 0 : use raw_rates, only: get_reaclib_rate_and_dlnT
656 : integer, intent(in) :: ir ! reaction_id
657 : real(dp), intent(in) :: temp, den
658 : real(dp), intent(inout) :: rate_raw(:), reverse_rate_raw(:)
659 : integer, intent(out) :: ierr
660 :
661 : real(dp) :: lambda, dlambda_dlnT, rlambda, drlambda_dlnT
662 :
663 : include 'formats'
664 :
665 : ierr = 0
666 : call get_reaclib_rate_and_dlnT( &
667 0 : ir, temp, lambda, dlambda_dlnT, rlambda, drlambda_dlnT, ierr)
668 0 : if (ierr /= 0) return
669 :
670 0 : if (reaction_ye_rho_exponents(2,ir) /= 1) then
671 0 : ierr = -1
672 0 : return
673 : end if
674 :
675 0 : rate_raw(i_rate) = lambda*den
676 0 : rate_raw(i_rate_dT) = dlambda_dlnT*den/temp
677 0 : rate_raw(i_rate_dRho) = lambda
678 :
679 0 : reverse_rate_raw(i_rate) = rlambda*den
680 0 : reverse_rate_raw(i_rate_dT) = drlambda_dlnT*den/temp
681 0 : reverse_rate_raw(i_rate_dRho) = rlambda
682 :
683 0 : return
684 :
685 : !$omp critical (rates_eval_reaclib_22)
686 : write(*,1) 'do_eval_reaclib_22 ' // trim(reaction_Name(ir))
687 : write(*,'(A)')
688 : write(*,1) 'den', den
689 : write(*,1) 'temp', temp
690 : write(*,'(A)')
691 : write(*,1) 'lambda', lambda
692 : write(*,1) 'dlambda_dlnT', dlambda_dlnT
693 : write(*,1) 'rate_raw', rate_raw(1:num_rvs)
694 : write(*,'(A)')
695 : write(*,1) 'rlambda', rlambda
696 : write(*,1) 'drlambda_dlnT', drlambda_dlnT
697 : write(*,1) 'reverse_rate_raw', reverse_rate_raw(1:num_rvs)
698 : write(*,'(A)')
699 : !call mesa_error(__FILE__,__LINE__,'do_eval_reaclib_22')
700 : !$omp end critical (rates_eval_reaclib_22)
701 :
702 0 : end subroutine do_eval_reaclib_22
703 :
704 :
705 : end module rates_support
706 :
|
