Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2013-2019 Haili Hu & 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 op_ev
21 :
22 : use math_lib
23 : use op_def
24 : use utils_lib
25 : use const_def, only: dp
26 : use kap_def, only: kap_test_partials, kap_test_partials_val, kap_test_partials_dval_dx
27 :
28 : implicit none
29 :
30 : private
31 : public :: abund
32 : public :: rd
33 : public :: ross
34 : public :: mix
35 : public :: interp
36 :
37 : contains
38 :
39 0 : subroutine abund(nel, izz, fa, flmu, fmu1, nkz)
40 : integer, intent(in) :: nel, izz(ipe)
41 : real, intent(in) :: fa(ipe)
42 : real, intent(out) :: flmu, fmu1(nrad)
43 : integer, intent(out) :: nkz(ipe)
44 : integer :: k, k1, k2, m
45 0 : real :: amamu(ipe), fmu, a1, c1, fmu0
46 : logical :: found
47 :
48 : ! Get k1, get amamu(k)
49 0 : do k = 1, nel
50 0 : found = .false.
51 0 : do m = 1, ipe
52 0 : if (izz(k) == kz(m)) then
53 0 : amamu(k) = amass(m)
54 0 : nkz(k) = m
55 : found = .true.
56 : exit
57 : end if
58 : end do
59 0 : if (.not. found) then
60 0 : write (*, *) ' k=', k, ', izz(k)=', izz(k)
61 0 : call mesa_error(__FILE__, __LINE__, ' kz(m) not found')
62 : end if
63 : end do
64 :
65 : ! Mean atomic weight = fmu
66 0 : fmu = 0.0d0
67 0 : do k = 1, nel
68 0 : fmu = fmu + fa(k)*amamu(k)
69 : end do
70 :
71 0 : do k2 = 1, nel
72 0 : a1 = fa(k2)
73 0 : c1 = 1.0d0/(1.0d0 - a1)
74 0 : fmu0 = c1*(fmu - fa(k2)*amamu(k2))
75 0 : fmu1(k2) = a1*(amamu(k2) - fmu0)*1.660531d-24 ! dmu/dlog xi
76 : end do
77 :
78 0 : fmu = fmu*1.660531d-24 ! Convert to cgs
79 0 : flmu = log10(dble(fmu))
80 :
81 0 : end subroutine abund
82 : ! **********************************************************************
83 :
84 0 : subroutine rd(nel, nkz, izz, ilab, jh, n_tot, ff, rr, i3, umesh, fac)
85 : integer, intent(in) :: nel, nkz(ipe), izz(ipe), ilab(4), jh(4), n_tot, i3
86 : real, intent(in) :: umesh(:) ! (nptot)
87 : real(dp), intent(in) :: fac(nel)
88 : real, intent(out) :: ff(:, :, :, :) ! (nptot, ipe, 4, 4)
89 : real, intent(out) :: rr(28, ipe, 4, 4)
90 : integer :: i, j, k, l, m, n, itt, jnn, izp, ne1, ne2, ne, ib, ia
91 0 : real :: fion(-1:28), yb, ya, d
92 : ! declare variables in common block (instead of by default: real (a-h, o-z), integer (i-n))
93 : ! integer :: ite1, ite2, ite3, jn1, jn2, jne3, ntot, nc, nf, int, ne1p, ne2p, np, kp1, kp2, kp3, npp, mx, nx
94 : ! real :: umin, umax, epatom, oplnck, fionp, yy1, yy2, yx
95 : ! common /atomdata/ ite1,ite2,ite3,jn1(91),jn2(91),jne3,umin,umax,ntot,
96 : ! + nc,nf,int(17),epatom(17,91,25),oplnck(17,91,25),ne1p(17,91,25),
97 : ! + ne2p(17,91,25),fionp(-1:28,28,91,25),np(17,91,25),kp1(17,91,25),
98 : ! + kp2(17,91,25),kp3(17,91,25),npp(17,91,25),mx(33417000),
99 : ! + yy1(33417000),yy2(120000000),nx(19305000),yx(19305000)
100 : ! save /atomdata/
101 :
102 : ! i=temperature index
103 : ! j=density index
104 : ! k=frequency index
105 : ! n=element index
106 : ! Get:
107 : ! mono opacity cross-section ff(k,n,i,j)
108 : ! modified cross-section for selected element, ta(k,i,j)
109 :
110 0 : rr = 0.0d0
111 :
112 : ! Start loop on i (temperature index)
113 0 : do i = 1, 4
114 0 : itt = (ilab(i) - ite1)/2 + 1
115 0 : do j = 1, 4
116 0 : jnn = (jh(j)*i3 - jn1(itt))/2 + 1
117 : ! Read mono opacities
118 0 : do n = 1, nel
119 0 : izp = izz(n)
120 0 : ne1 = ne1p(nkz(n), itt, jnn)
121 0 : ne2 = ne2p(nkz(n), itt, jnn)
122 0 : do ne = ne1, ne2
123 0 : fion(ne) = fionp(ne, nkz(n), itt, jnn)
124 : end do
125 0 : do ne = ne1, min(ne2, izp - 2)
126 0 : rr(izp - 1 - ne, n, i, j) = fion(ne)
127 : end do
128 : ! call zetbarp(izp, ne1, ne2, fion, zet, i3)
129 : ! zetal(n, i, j) = zet
130 0 : do k = 1, n_tot
131 0 : ff(k, n, i, j) = yy2(k + kp2(nkz(n), itt, jnn))
132 : end do
133 :
134 0 : if (fac(n) /= 1d0) then
135 0 : do k = 1, size(ff, dim=1)
136 0 : ff(k, n, i, j) = fac(n)*ff(k, n, i, j)
137 : end do
138 : end if
139 : end do
140 : end do
141 : end do
142 :
143 0 : end subroutine rd
144 :
145 0 : subroutine ross(flmu, fmu1, dv, ntot, rs, rossl)
146 : integer, intent(in) :: ntot
147 : real, intent(in) :: flmu, dv, fmu1(nrad)
148 : real, intent(in) :: rs(:, :, :) ! (nptot, 4, 4)
149 : real, intent(out) :: rossl(4, 4)
150 : integer :: i, j, n, k2
151 0 : real(dp) :: drs, dd, oross
152 0 : real :: fmu, tt, ss, dd2, drsp(nrad), exp10_flmu
153 : real, parameter :: log10_bohr_radius_sqr = -16.55280d0
154 :
155 : ! oross=cross-section in a.u.
156 : ! rossl=log10(ROSS in cgs)
157 0 : exp10_flmu = real(exp10(dble(flmu)))
158 0 : do i = 1, 4
159 0 : do j = 1, 4
160 0 : drs = 0.d0
161 0 : do n = 1, ntot
162 0 : dd = 1.d0/rs(n, i, j)
163 0 : dd2 = dd*dd
164 0 : drs = drs + dd
165 : end do
166 0 : oross = 1.d0/(drs*dv)
167 0 : rossl(i, j) = log10(dble(oross)) + log10_bohr_radius_sqr - flmu !log10(fmu)
168 : end do
169 : end do
170 :
171 0 : end subroutine ross
172 :
173 0 : subroutine mix(ntot, nel, fa, ff, rs, rr, rion)
174 : integer, intent(in) :: ntot, nel
175 : real, intent(in) :: ff(:, :, :, :) ! (nptot, ipe, 4, 4)
176 : real, intent(in) :: fa(ipe), rr(28, 17, 4, 4)
177 : real, intent(out) :: rs(:, :, :) ! (nptot, 4, 4)
178 : real, intent(out) :: rion(28, 4, 4)
179 : integer :: i, j, k, n, m
180 0 : real :: rs_temp, rion_temp
181 :
182 0 : do j = 1, 4
183 0 : do i = 1, 4
184 0 : do n = 1, ntot
185 0 : rs_temp = ff(n, 1, i, j)*fa(1)
186 0 : do k = 2, nel
187 0 : rs_temp = rs_temp + ff(n, k, i, j)*fa(k)
188 : end do
189 0 : rs(n, i, j) = rs_temp
190 : !rs(n, i, j) = dot_product(ff(n,1:nel,i,j),fa(1:nel))
191 : end do
192 0 : do m = 1, 28
193 0 : rion_temp = rr(m, 1, i, j)*fa(1)
194 0 : do k = 2, nel
195 0 : rion_temp = rion_temp + rr(m, k, i, j)*fa(k)
196 : end do
197 0 : rion(m, i, j) = rion_temp
198 : !rion(m,i,j) = dot_product(rr(m,1:nel,i,j),fa(1:nel))
199 : end do
200 : end do
201 : end do
202 :
203 0 : end subroutine mix
204 :
205 0 : subroutine interp(nel, rossl, xi, eta, g, i3, ux, uy, gx, gy)
206 : use op_common, only: fint, fintp
207 : integer, intent(in) :: nel, i3
208 : real, intent(in) :: ux, uy, xi, eta, rossl(4, 4)
209 : real, intent(out) :: gx, gy, g
210 : integer :: i, j, l, k2
211 0 : real :: V(4), U(4), vyi(4)
212 :
213 0 : do i = 1, 4
214 0 : do j = 1, 4
215 0 : U(j) = rossl(i, j)
216 : end do
217 0 : V(i) = fint(U, eta)
218 0 : vyi(i) = fintp(U, eta)
219 : end do
220 0 : g = fint(V, xi)
221 0 : gy = fint(vyi, xi)
222 0 : gx = fintp(V, xi)
223 :
224 0 : gy = gy/uy
225 0 : gx = (80.0d0/real(i3))*(gx - gy*ux)
226 :
227 0 : end subroutine interp
228 :
229 : end module op_ev
|
