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 screen
21 : use const_def, only: dp, amu, pi4, one_third, two_thirds, four_thirds, five_thirds
22 : use rates_def
23 : use math_lib
24 :
25 : implicit none
26 :
27 : contains
28 :
29 89116 : subroutine do_screen_set_context( &
30 : sc, temp, den, logT, logRho, zbar, abar, z2bar, &
31 89116 : screening_mode, num_isos, y, iso_z158)
32 : type (Screen_Info) :: sc
33 : integer, intent(in) :: num_isos
34 : real(dp), intent(in) :: &
35 : temp, den, logT, logRho, zbar, abar, z2bar, &
36 : y(:), &
37 : iso_z158(:) ! Z**1.58
38 : integer, intent(in) :: screening_mode
39 :
40 : real(dp), parameter :: x13 = 1.0d0/3.0d0
41 : real(dp), parameter :: x14 = 1.0d0/4.0d0
42 : real(dp), parameter :: x53 = 5.0d0/3.0d0
43 : real(dp), parameter :: x532 = 5.0d0/32.0d0
44 : real(dp), parameter :: x512 = 5.0d0/12.0d0
45 : real(dp), parameter :: fact = 1.25992104989487d0 ! the cube root of 2
46 : real(dp), parameter :: co2 = x13 * 4.248719d3
47 : real(dp) :: qq
48 : integer :: j
49 :
50 : logical, parameter :: debug = .false.
51 : !logical, parameter :: debug = .true.
52 :
53 : include 'formats'
54 :
55 89116 : if (screening_mode == no_screening .or. zbar == 0d0) return
56 :
57 89114 : sc% temp = temp
58 89114 : sc% den = den
59 89114 : sc% logT = logT
60 89114 : sc% logRho = logRho
61 89114 : sc% zbar = zbar
62 89114 : sc% abar = abar
63 89114 : sc% z2bar = z2bar
64 :
65 : ! get the info that depends only on temp, den, and overall composition
66 :
67 89114 : sc% ytot = 1.0d0/abar
68 89114 : sc% rr = den * sc% ytot
69 89114 : sc% tempi = 1.0d0/temp
70 89114 : sc% dtempi = -sc% tempi * sc% tempi
71 89114 : sc% deni = 1.0d0/den
72 89114 : qq = 0d0
73 920446 : do j=1,num_isos
74 831332 : if (iso_z158(j) == 0d0) cycle
75 920446 : qq = qq + iso_z158(j) * y(j)
76 : end do
77 89114 : sc% z1pt58bar = abar * qq
78 89114 : sc% zbar13 = pow(zbar,1d0/3d0)
79 :
80 89114 : sc% pp = sqrt(sc% rr * sc% tempi * (z2bar + zbar))
81 89114 : qq = 0.5d0/(sc% pp) *(z2bar + zbar)
82 89114 : sc% dppdt = qq*sc% rr*sc% dtempi
83 89114 : sc% dppdd = qq*sc% ytot*sc% tempi
84 :
85 89114 : sc% qlam0z = 1.88d8 * sc% tempi * sc% pp
86 89114 : sc% qlam0zdt = 1.88d8 * (sc% dtempi * sc% pp + sc% tempi * sc% dppdt)
87 89114 : sc% qlam0zdd = 1.88d8 * sc% tempi * sc% dppdd
88 :
89 89114 : sc% taufac = co2 * pow(sc% tempi,x13)
90 89114 : sc% taufacdt = -x13*sc% taufac*sc% tempi
91 :
92 89114 : qq = sc% rr*zbar
93 89114 : sc% xni = pow(qq,x13)
94 89114 : sc% dxnidd = x13 * sc% xni * sc% deni
95 :
96 89114 : sc% aa = 2.27493d5 * sc% tempi * sc% xni
97 89114 : sc% daadt = 2.27493d5 * sc% dtempi * sc% xni
98 89114 : sc% daadd = 2.27493d5 * sc% tempi * sc% dxnidd
99 :
100 : ! ion and electron sphere radii (itoh 1979 eq 1-3)
101 89114 : sc% ntot = den / (amu*abar)
102 89114 : sc% a_e = pow((3.d0 /(pi4 * zbar * sc% ntot)),x13)
103 :
104 : end subroutine do_screen_set_context
105 :
106 :
107 : end module screen
108 :
|
