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 net_screen
21 :
22 : use const_def,only: dp, pi, ln10
23 : use math_lib
24 : use chem_def, only: chem_isos, ih1, num_chem_isos
25 : use net_def, only: Net_General_Info, Net_Info
26 : use rates_def
27 : use utils_lib, only: mesa_error
28 :
29 : implicit none
30 :
31 :
32 : contains
33 :
34 :
35 19 : subroutine make_screening_tables(n, ierr)
36 : type (Net_Info) :: n
37 : integer, intent(out) :: ierr
38 149283 : real(dp) :: y(num_chem_isos)
39 149283 : y = 0
40 : call screen_net( &
41 : n% g, num_chem_isos, y, 1d0, 1d0, 0d0, 0d0, .true., &
42 : n% rate_raw, n% rate_raw_dT, n% rate_raw_dRho, &
43 : n% rate_screened, n% rate_screened_dT, n% rate_screened_dRho, &
44 : n% screening_mode, &
45 19 : 0d0, 0d0, 0d0, 1d0, ierr)
46 19 : end subroutine make_screening_tables
47 :
48 :
49 89127 : subroutine screen_net( &
50 89127 : g, num_isos, y, temp, den, logT, logRho, init, &
51 534762 : rate_raw, rate_raw_dT, rate_raw_dRho, &
52 89127 : rate_screened, rate_screened_dT, rate_screened_dRho, &
53 : screening_mode, &
54 : zbar, abar, z2bar, ye, ierr)
55 :
56 : use rates_def, only: Screen_Info, reaction_name
57 : use rates_lib, only: screen_set_context
58 :
59 : type (Net_General_Info), pointer :: g
60 : integer, intent(in) :: num_isos, screening_mode
61 : real(dp), intent(in) :: y(:), temp, den, logT, logRho, &
62 : zbar, abar, z2bar, ye
63 : real(dp), intent(inout), dimension(:) :: &
64 : rate_raw, rate_raw_dT, rate_raw_dRho, &
65 : rate_screened, rate_screened_dT, rate_screened_dRho
66 : logical, intent(in) :: init
67 : integer, intent(out) :: ierr
68 :
69 : type (Screen_Info) :: sc
70 : integer :: num_reactions, i, ir, j, op_err
71 89127 : real(dp) :: Tfactor, dTfactordt
72 : logical :: all_okay
73 :
74 : include 'formats'
75 :
76 89127 : ierr = 0
77 :
78 89127 : if (.not. init) then
79 : call screen_set_context( &
80 : sc, temp, den, logT, logRho, zbar, abar, z2bar, &
81 89108 : screening_mode, num_isos, y, g% z158)
82 : end if
83 :
84 89127 : num_reactions = g% num_reactions
85 :
86 3417460 : do i = 1, num_reactions
87 3328333 : ir = g% reaction_id(i)
88 3328333 : if (ir == 0) then
89 0 : write(*,*) 'g% reaction_id(i) == 0', i, num_reactions
90 0 : call mesa_error(__FILE__,__LINE__,'screen_net')
91 : end if
92 3417460 : if (reaction_screening_info(3,ir) > 0) then
93 : call eval_screen_triple( &
94 : init, i, &
95 : reaction_screening_info(1,ir), &
96 : reaction_screening_info(2,ir), &
97 : reaction_screening_info(3,ir), &
98 89127 : i, sc, ir, ierr)
99 89127 : if (ierr /= 0) then
100 : write(*,*) 'screen_net failed in eval_screen_triple ' // &
101 0 : trim(reaction_name(ir))
102 0 : return
103 : end if
104 3239206 : else if (reaction_screening_info(2,ir) > 0) then
105 : call eval_screen_pair( &
106 : init, i, &
107 : reaction_screening_info(1,ir), &
108 : reaction_screening_info(2,ir), &
109 2403532 : i, sc, ir, ierr)
110 2403532 : if (ierr /= 0) then
111 : write(*,*) 'screen_net failed in eval_screen_pair ' // &
112 0 : trim(reaction_name(ir))
113 0 : return
114 : end if
115 : else
116 835674 : rate_screened(i) = rate_raw(i)
117 835674 : rate_screened_dT(i) = rate_raw_dT(i)
118 835674 : rate_screened_dRho(i) = rate_raw_dRho(i)
119 : end if
120 : end do
121 89127 : if (ierr /= 0) return
122 :
123 89127 : if(init) return
124 :
125 89108 : call set_combo_screen_rates(num_isos, y, sc, ierr)
126 89108 : if (ierr /= 0) then
127 0 : write(*,*) 'screen_net failed in set_combo_screen_rates'
128 0 : return
129 : end if
130 :
131 178235 : if (nrattab > 1 .and. (logT < g% logTcut_lim .or. logT <= g% logTcut_lo)) then
132 : ! strong rates cutoff smoothly for logT < logTcut_lim
133 2 : if (logT <= g% logTcut_lo) then
134 190 : do i = 1, num_reactions
135 188 : if (g% weak_reaction_index(i) > 0) cycle
136 180 : rate_screened(i) = 0
137 180 : rate_screened_dT(i) = 0
138 190 : rate_screened_dRho(i) = 0
139 : end do
140 : else
141 0 : Tfactor = (logT - g% logTcut_lo)/(g% logTcut_lim - g% logTcut_lo)
142 0 : Tfactor = 0.5d0*(1 - cospi(Tfactor*Tfactor))
143 : dTfactordt = 0.5d0 * pi * sinpi(Tfactor*Tfactor) * &
144 0 : 2.d0/((g% logTcut_lim - g% logTcut_lo) * temp * ln10)
145 0 : do i = 1, num_reactions
146 0 : if (g% weak_reaction_index(i) > 0) cycle
147 0 : rate_screened_dT(i) = Tfactor * rate_screened_dT(i) + dTfactordt * rate_screened(i)
148 0 : rate_screened_dRho(i) = Tfactor * rate_screened_dRho(i)
149 89106 : rate_screened(i) = Tfactor * rate_screened(i)
150 : end do
151 : end if
152 : end if
153 :
154 :
155 : contains
156 :
157 2581786 : subroutine screening_pair( &
158 : init, ir, jscr, sc, cid1, a1, z1, cid2, a2, z2, scor, scordt, scordd, ierr)
159 89127 : use rates_lib, only: screen_init_AZ_info, screen_pair
160 : use rates_def, only: Screen_Info
161 : use chem_def, only: ih1, ih2, ihe4, ico55, ico57
162 : logical, intent(in) :: init
163 : integer, intent(in) :: ir, jscr
164 : type (Screen_Info) :: sc
165 : integer, intent(in) :: cid1, cid2
166 : real(dp), intent(in) :: a1, z1, a2, z2
167 : real(dp), intent(out) :: scor, scordt, scordd
168 : integer, intent(out) :: ierr
169 :
170 : integer :: i1, i2
171 : include 'formats'
172 2581786 : ierr = 0
173 2581786 : if (init) then
174 : call screen_init_AZ_info( &
175 : a1, z1, a2, z2, &
176 : g% zs13(jscr), g% zhat(jscr), g% zhat2(jscr), g% lzav(jscr), &
177 : g% aznut(jscr), g% zs13inv(jscr), &
178 844 : ierr)
179 844 : if (ierr /= 0) write(*,*) 'screen_init_AZ_info failed in screening_pair ' // &
180 0 : trim(reaction_name(ir))
181 : else
182 2580942 : if (cid1 > 0 .and. cid2 > 0) then
183 2491834 : i1 = g% net_iso(cid1)
184 2491834 : i2 = g% net_iso(cid2)
185 2491834 : if (i1 == 0 .or. i2 == 0) then ! not in current net
186 431228 : if (g% doing_approx21 .and. &
187 : .not. (cid1 == ico55 .or. cid1 == ico57 .or. &
188 : cid2 == ico55 .or. cid2 == ico57 .or. &
189 : cid2 == ih2 .or. cid2 == ih2)) then
190 : ! this skips screening things like al27 + p for approx21
191 154802 : scor = 1d0
192 154802 : scordt = 0d0
193 154802 : scordd = 0d0
194 : return
195 : end if
196 : end if
197 : end if
198 : call screen_pair( &
199 : sc, a1, z1, a2, z2, screening_mode, &
200 : g% zs13(jscr), g% zhat(jscr), g% zhat2(jscr), g% lzav(jscr), &
201 : g% aznut(jscr), g% zs13inv(jscr), g% logTcut_lo, &
202 2426140 : scor, scordt, scordd, ierr)
203 2426140 : if (ierr /= 0) write(*,*) 'screen_pair failed in screening_pair ' // &
204 0 : trim(reaction_name(ir))
205 : end if
206 2581786 : end subroutine screening_pair
207 :
208 2491834 : subroutine set_rate_screening(i, sc1a, sc1adt, sc1add)
209 : integer, intent(in) :: i
210 : real(dp), intent(in) :: sc1a, sc1adt, sc1add
211 : include 'formats'
212 2491834 : if (i == 0) return
213 2491834 : rate_screened(i) = rate_raw(i)*sc1a
214 2491834 : rate_screened_dT(i) = rate_raw_dT(i)*sc1a + rate_raw(i)*sc1adt
215 2491834 : rate_screened_dRho(i) = rate_raw_dRho(i)*sc1a + rate_raw(i)*sc1add
216 2581786 : end subroutine set_rate_screening
217 :
218 2403532 : subroutine eval_screen_pair(init, jscr, i1, i2, i, sc, ir, ierr)
219 : use rates_def, only: Screen_Info
220 : logical, intent(in) :: init
221 : integer, intent(in) :: jscr
222 : type (Screen_Info) :: sc
223 : integer, intent(in) :: i1, i2 ! chem id's for the isotopes
224 : integer, intent(in) :: i ! rate number
225 : integer, intent(in) :: ir
226 : integer, intent(out) :: ierr
227 : real(dp) :: sc1a, sc1adt, sc1add, a1, z1, a2, z2
228 : include 'formats'
229 : ierr = 0
230 2403532 : a1 = chem_isos% Z_plus_N(i1)
231 2403532 : z1 = dble(chem_isos% Z(i1))
232 2403532 : a2 = chem_isos% Z_plus_N(i2)
233 2403532 : z2 = dble(chem_isos% Z(i2))
234 : !if (z1 == 0d0 .or. z2 == 0d0) return
235 : ! with this, get bad burn result, but okay for restart
236 : ! without it, reversed. get good burn, bad restart.
237 : ! bad restart max diff ~ 5e-15 in abundances of n17, n18, o14
238 : ! perhaps tiny difference in some screening factor?
239 : call screening_pair( &
240 2403532 : init, ir, jscr, sc, i1, a1, z1, i2, a2, z2, sc1a, sc1adt, sc1add, ierr)
241 2403532 : if (ierr /= 0) return
242 2403532 : if (init) return
243 2402726 : call set_rate_screening(i, sc1a, sc1adt, sc1add)
244 2403532 : end subroutine eval_screen_pair
245 :
246 89127 : subroutine eval_screen_triple(init, jscr, i1_in, i2_in, i3_in, i, sc, ir, ierr)
247 2403532 : use rates_def, only: Screen_Info
248 : logical, intent(in) :: init
249 : integer, intent(in) :: jscr
250 : type (Screen_Info) :: sc
251 : integer, intent(in) :: i1_in, i2_in, i3_in ! chem id's for the isotopes
252 : integer, intent(in) :: i ! rate number
253 : integer, intent(in) :: ir
254 : integer, intent(out) :: ierr
255 : integer :: i1, i2, i3, ii
256 89127 : real(dp) :: sc1, sc1dt, sc1dd
257 89127 : real(dp) :: sc2, sc2dt, sc2dd
258 89127 : real(dp) :: scor, scordt, scordd
259 : real(dp) :: a1, z1, a2, z2, a3, z3
260 : include 'formats'
261 89127 : ierr = 0
262 89127 : i1 = i1_in; i2 = i2_in; i3 = i3_in
263 89127 : a1 = chem_isos% Z_plus_N(i1)
264 89127 : z1 = dble(chem_isos% Z(i1))
265 89127 : a2 = chem_isos% Z_plus_N(i2)
266 89127 : z2 = dble(chem_isos% Z(i2))
267 89127 : a3 = chem_isos% Z_plus_N(i3)
268 89127 : z3 = dble(chem_isos% Z(i3))
269 89127 : if (z2 == 0) then
270 0 : if (z1 == 0) return ! n + n + A
271 : ! have A + n + B
272 : ! swap 1 and 2 so have n + A + B
273 0 : ii = i2; i2 = i1; i1 = ii
274 0 : a1 = chem_isos% Z_plus_N(i1)
275 0 : z1 = dble(chem_isos% Z(i1))
276 0 : a2 = chem_isos% Z_plus_N(i2)
277 0 : z2 = dble(chem_isos% Z(i2))
278 : end if
279 89127 : if (z3 == 0) then ! have A + B + n
280 : ! swap 1 and 3 so have n + A + B
281 0 : ii = i1; i1 = i3; i3 = ii
282 0 : a1 = chem_isos% Z_plus_N(i1)
283 0 : z1 = dble(chem_isos% Z(i1))
284 0 : a3 = chem_isos% Z_plus_N(i3)
285 0 : z3 = dble(chem_isos% Z(i3))
286 : end if
287 : call screening_pair( &
288 89127 : init, ir, jscr, sc, i2, a2, z2, i3, a3, z3, sc2, sc2dt, sc2dd, ierr)
289 89127 : if (ierr /= 0) return
290 89127 : if (z1 == 0) then
291 0 : if (init) return
292 0 : call set_rate_screening(i, sc2, sc2dt, sc2dd)
293 0 : return ! n + (A + B)
294 : end if
295 89127 : i2 = 0 ! 0 for cid to disable caching
296 89127 : a2 = a2 + a3
297 89127 : z2 = z2 + z3
298 : call screening_pair( &
299 89127 : init, ir, jscr, sc, i1, a1, z1, i2, a2, z2, sc1, sc1dt, sc1dd, ierr)
300 89127 : if (init) return
301 89108 : scor = sc1*sc2
302 89108 : scordt = sc1*sc2dt + sc1dt*sc2
303 89108 : scordd = sc1*sc2dd + sc1dd*sc2
304 89108 : call set_rate_screening(i, scor, scordt, scordd)
305 :
306 : if (.false.) write(*,2) 'scr 3 ' // trim(reaction_Name(ir)) &
307 : // ' ' // trim(chem_isos% name(i1)) &
308 : // ' ' // trim(chem_isos% name(i2)) &
309 : // ' ' // trim(chem_isos% name(i3)), &
310 : ir, scor
311 :
312 89127 : end subroutine eval_screen_triple
313 :
314 89108 : subroutine set_combo_screen_rates(num_isos, y, sc, ierr)
315 89127 : use rates_def, only: Screen_Info
316 : integer, intent(in) :: num_isos
317 : real(dp), intent(in) :: y(:)
318 : type (Screen_Info) :: sc
319 : integer, intent(out) :: ierr
320 :
321 89108 : integer, pointer :: rtab(:)
322 89108 : real(dp) :: rateII, rateIII, rsum, fII, fIII
323 :
324 : include 'formats'
325 :
326 89108 : rtab => g% net_reaction
327 89108 : ierr = 0
328 :
329 80000 : if (rtab(ir34_pp2) /= 0 .and. rtab(ir34_pp3) /= 0) then
330 80000 : if (rate_screened(rtab(ir34_pp2)) /= &
331 : rate_screened(rtab(ir34_pp3))) then
332 0 : ierr = -1
333 : return
334 : end if
335 80000 : if (rtab(ir_be7_wk_li7) /= 0) then
336 0 : rateII = rate_screened(rtab(ir_be7_wk_li7))
337 80000 : else if (rtab(irbe7ec_li7_aux) /= 0) then
338 80000 : rateII = rate_screened(rtab(irbe7ec_li7_aux))
339 : else
340 0 : write(*,*) 'need either r_be7_wk_li7 or rbe7ec_li7_aux'
341 0 : call mesa_error(__FILE__,__LINE__,'set_combo_screen_rates')
342 : end if
343 80000 : if (rtab(ir_be7_pg_b8) /= 0) then
344 0 : rateIII = y(g% net_iso(ih1)) * rate_screened(rtab(ir_be7_pg_b8))
345 80000 : else if (rtab(irbe7pg_b8_aux) /= 0) then
346 80000 : rateIII = y(g% net_iso(ih1)) * rate_screened(rtab(irbe7pg_b8_aux))
347 : else
348 0 : write(*,*) 'need either r_be7_pg_b8 or rbe7pg_b8_aux'
349 0 : call mesa_error(__FILE__,__LINE__,'set_combo_screen_rates')
350 : end if
351 80000 : rsum = rateII + rateIII
352 80000 : if (rsum < 1d-50) then
353 : fII = 0.5d0
354 : else
355 52200 : fII = rateII / rsum
356 : end if
357 80000 : fIII = 1d0 - fII
358 :
359 80000 : rate_screened(rtab(ir34_pp2)) = fII*rate_screened(rtab(ir34_pp2))
360 80000 : rate_screened_dT(rtab(ir34_pp2)) = fII*rate_screened_dT(rtab(ir34_pp2))
361 80000 : rate_screened_dRho(rtab(ir34_pp2)) = fII*rate_screened_dRho(rtab(ir34_pp2))
362 :
363 80000 : rate_screened(rtab(ir34_pp3)) = fIII*rate_screened(rtab(ir34_pp3))
364 80000 : rate_screened_dT(rtab(ir34_pp3)) = fIII*rate_screened_dT(rtab(ir34_pp3))
365 80000 : rate_screened_dRho(rtab(ir34_pp3)) = fIII*rate_screened_dRho(rtab(ir34_pp3))
366 :
367 : end if
368 :
369 89108 : if (rtab(irn14_to_c12) /= 0) &
370 : call rate_for_pg_pa_branches( &
371 : rtab(irn14pg_aux), rtab(irn15pg_aux), rtab(irn15pa_aux), &
372 80000 : 0, rtab(irn14_to_c12))
373 :
374 89108 : if (rtab(irn14_to_o16) /= 0) &
375 : call rate_for_pg_pa_branches( &
376 : rtab(irn14pg_aux), rtab(irn15pg_aux), rtab(irn15pa_aux), &
377 80000 : rtab(irn14_to_o16), 0)
378 :
379 89108 : if (rtab(ir1616ppa) /= 0) &
380 : call rate_for_pg_pa_branches( &
381 : rtab(ir1616p_aux), rtab(irp31pg_aux), rtab(irp31pa_aux), &
382 0 : 0, rtab(ir1616ppa))
383 :
384 89108 : if (rtab(ir1616ppg) /= 0) &
385 : call rate_for_pg_pa_branches( &
386 : rtab(ir1616p_aux), rtab(irp31pg_aux), rtab(irp31pa_aux), &
387 0 : rtab(ir1616ppg), 0)
388 :
389 : call rate_for_alpha_ap( &
390 : irc12ap_aux, irn15pg_aux, irn15pa_aux, &
391 89108 : irc12ap_to_o16)
392 :
393 : call rate_for_alpha_gp( &
394 : iro16gp_aux, irn15pg_aux, irn15pa_aux, &
395 89108 : iro16gp_to_c12)
396 :
397 : call rate_for_alpha_ap( &
398 : iro16ap_aux, irf19pg_aux, irf19pa_aux, &
399 89108 : iro16ap_to_ne20)
400 :
401 : call rate_for_alpha_gp( &
402 : irne20gp_aux, irf19pg_aux, irf19pa_aux, &
403 89108 : irne20gp_to_o16)
404 :
405 : call rate_for_alpha_ap( &
406 : irne20ap_aux, irna23pg_aux, irna23pa_aux, &
407 89108 : irne20ap_to_mg24)
408 :
409 : call rate_for_alpha_gp( &
410 : irmg24gp_aux, irna23pg_aux, irna23pa_aux, &
411 89108 : irmg24gp_to_ne20)
412 :
413 : call rate_for_alpha_ap( &
414 : irmg24ap_aux, iral27pg_aux, iral27pa_aux, &
415 89108 : irmg24ap_to_si28)
416 :
417 : call rate_for_alpha_gp( &
418 : irsi28gp_aux, iral27pg_aux, iral27pa_aux, &
419 89108 : irsi28gp_to_mg24)
420 :
421 : call rate_for_alpha_ap( &
422 : irsi28ap_aux, irp31pg_aux, irp31pa_aux, &
423 89108 : irsi28ap_to_s32)
424 :
425 : call rate_for_alpha_gp( &
426 : irs32gp_aux, irp31pg_aux, irp31pa_aux, &
427 89108 : irs32gp_to_si28)
428 :
429 : call rate_for_alpha_ap( &
430 : irs32ap_aux, ircl35pg_aux, ircl35pa_aux, &
431 89108 : irs32ap_to_ar36)
432 :
433 : call rate_for_alpha_gp( &
434 : irar36gp_aux, ircl35pg_aux, ircl35pa_aux, &
435 89108 : irar36gp_to_s32)
436 :
437 : call rate_for_alpha_ap( &
438 : irar36ap_aux, irk39pg_aux, irk39pa_aux, &
439 89108 : irar36ap_to_ca40)
440 :
441 : call rate_for_alpha_gp( &
442 : irca40gp_aux, irk39pg_aux, irk39pa_aux, &
443 89108 : irca40gp_to_ar36)
444 :
445 : call rate_for_alpha_ap( &
446 : irca40ap_aux, irsc43pg_aux, irsc43pa_aux, &
447 89108 : irca40ap_to_ti44)
448 :
449 : call rate_for_alpha_gp( &
450 : irti44gp_aux, irsc43pg_aux, irsc43pa_aux, &
451 89108 : irti44gp_to_ca40)
452 :
453 : call rate_for_alpha_ap( &
454 : irti44ap_aux, irv47pg_aux, irv47pa_aux, &
455 89108 : irti44ap_to_cr48)
456 :
457 : call rate_for_alpha_gp( &
458 : ircr48gp_aux, irv47pg_aux, irv47pa_aux, &
459 89108 : ircr48gp_to_ti44)
460 :
461 : call rate_for_alpha_ap( &
462 : ircr48ap_aux, irmn51pg_aux, irmn51pa_aux, &
463 89108 : ircr48ap_to_fe52)
464 :
465 : call rate_for_alpha_gp( &
466 : irfe52gp_aux, irmn51pg_aux, irmn51pa_aux, &
467 89108 : irfe52gp_to_cr48)
468 :
469 :
470 178216 : end subroutine set_combo_screen_rates
471 :
472 891080 : subroutine rate_for_alpha_ap(ir_start, irpg, irpa, ir_with_pg)
473 : integer, intent(in) :: ir_start, irpg, irpa, ir_with_pg
474 891080 : integer, pointer :: rtab(:)
475 : include 'formats'
476 651078 : if (ir_start == 0) return
477 891080 : rtab => g% net_reaction
478 891080 : if (rtab(ir_with_pg) == 0) return
479 : call rate_for_pg_pa_branches( &
480 240002 : rtab(ir_start), rtab(irpg), rtab(irpa), rtab(ir_with_pg), 0)
481 980188 : end subroutine rate_for_alpha_ap
482 :
483 891080 : subroutine rate_for_alpha_gp(ir_start, irpg, irpa, ir_with_pa)
484 : integer, intent(in) :: ir_start, irpg, irpa, ir_with_pa
485 891080 : integer, pointer :: rtab(:)
486 891080 : if (ir_start == 0) return
487 891080 : rtab => g% net_reaction
488 891080 : if (rtab(ir_with_pa) == 0) return
489 : call rate_for_pg_pa_branches( &
490 0 : rtab(ir_start), rtab(irpg), rtab(irpa), 0, rtab(ir_with_pa))
491 891080 : end subroutine rate_for_alpha_gp
492 :
493 400002 : subroutine rate_for_pg_pa_branches(ir_start, irpg, irpa, ir_with_pg, ir_with_pa)
494 : integer, intent(in) :: ir_start, irpg, irpa, ir_with_pg, ir_with_pa
495 :
496 400002 : real(dp) :: pg_raw_rate, pa_raw_rate, pg_frac, pa_frac
497 400002 : real(dp) :: d_pg_frac_dT, d_pg_frac_dRho, d_pa_frac_dT, d_pa_frac_dRho
498 400002 : real(dp) :: r, drdT, drdd, x
499 :
500 400002 : if (ir_start == 0) then
501 0 : write(*,*) 'ir_start', ir_start
502 0 : if (irpg /= 0) write(*,*) trim(reaction_Name(g% reaction_id(irpg))) // ' irpg'
503 0 : if (irpa /= 0) write(*,*) trim(reaction_Name(g% reaction_id(irpa))) // ' irpa'
504 0 : if (ir_with_pg /= 0) write(*,*) trim(reaction_Name(g% reaction_id(ir_with_pg))) // ' ir_with_pg'
505 0 : if (ir_with_pa /= 0) write(*,*) trim(reaction_Name(g% reaction_id(ir_with_pa))) // ' ir_with_pa'
506 0 : call mesa_error(__FILE__,__LINE__,'rate_for_pg_pa_branches')
507 : end if
508 :
509 400002 : if (irpg == 0) then
510 0 : write(*,*) 'irpg', irpg
511 0 : if (ir_with_pg /= 0) write(*,*) trim(reaction_Name(g% reaction_id(ir_with_pg))) // ' ir_with_pg'
512 0 : if (ir_with_pa /= 0) write(*,*) trim(reaction_Name(g% reaction_id(ir_with_pa))) // ' ir_with_pa'
513 0 : call mesa_error(__FILE__,__LINE__,'rate_for_pg_pa_branches')
514 : end if
515 :
516 400002 : if (irpa == 0) then
517 0 : write(*,*) 'irpg', irpg
518 0 : if (ir_with_pg /= 0) write(*,*) trim(reaction_Name(g% reaction_id(ir_with_pg))) // ' ir_with_pg'
519 0 : if (ir_with_pa /= 0) write(*,*) trim(reaction_Name(g% reaction_id(ir_with_pa))) // ' ir_with_pa'
520 0 : call mesa_error(__FILE__,__LINE__,'rate_for_pg_pa_branches')
521 : end if
522 :
523 400002 : pg_raw_rate = rate_raw(irpg)
524 400002 : pa_raw_rate = rate_raw(irpa)
525 :
526 400002 : if (pg_raw_rate + pa_raw_rate < 1d-99) then ! avoid divide by 0
527 245796 : pg_raw_rate = 1; pa_raw_rate = 1
528 : end if
529 :
530 400002 : pg_frac = pg_raw_rate / (pg_raw_rate + pa_raw_rate)
531 400002 : pa_frac = 1 - pg_frac
532 :
533 400002 : x = pg_raw_rate + pa_raw_rate
534 : d_pg_frac_dT = &
535 400002 : (pa_raw_rate*rate_raw_dT(irpg) - pg_raw_rate*rate_raw_dT(irpa)) / (x*x)
536 400002 : d_pa_frac_dT = -d_pg_frac_dT
537 :
538 : d_pg_frac_dRho = &
539 400002 : (pa_raw_rate*rate_raw_dRho(irpg) - pg_raw_rate*rate_raw_dRho(irpa)) / (x*x)
540 400002 : d_pa_frac_dRho = -d_pg_frac_dRho
541 :
542 400002 : r = rate_screened(ir_start)
543 400002 : drdT = rate_screened_dT(ir_start)
544 400002 : drdd = rate_screened_dRho(ir_start)
545 :
546 400002 : if (ir_with_pg /= 0) then
547 320002 : rate_screened(ir_with_pg) = r*pg_frac
548 320002 : rate_screened_dT(ir_with_pg) = r*d_pg_frac_dT + drdT*pg_frac
549 320002 : rate_screened_dRho(ir_with_pg) = r*d_pg_frac_dRho + drdd*pg_frac
550 : end if
551 :
552 400002 : if (ir_with_pa /= 0) then
553 80000 : rate_screened(ir_with_pa) = r*pa_frac
554 80000 : rate_screened_dT(ir_with_pa) = r*d_pa_frac_dT + drdT*pa_frac
555 80000 : rate_screened_dRho(ir_with_pa) = r*d_pa_frac_dRho + drdd*pa_frac
556 : end if
557 :
558 400002 : end subroutine rate_for_pg_pa_branches
559 :
560 :
561 : end subroutine screen_net
562 :
563 :
564 :
565 : end module net_screen
566 :
|
