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_radacc
21 :
22 : use math_lib
23 : use op_def
24 : use const_def, only: dp
25 :
26 : implicit none
27 :
28 : private
29 : public :: rd
30 : public :: abund
31 : public :: ross
32 : public :: mix
33 : public :: interp
34 :
35 : logical, parameter :: dbg = .false.
36 :
37 : contains
38 :
39 0 : subroutine abund(nel, izz, kk, iz1, fa, kzz, flmu, am1, fmu1, nkz)
40 : implicit none
41 : integer, intent(in) :: nel, izz(ipe), kk, iz1(nrad)
42 : real, intent(in) :: fa(ipe)
43 : integer, intent(out) :: kzz(nrad), nkz(ipe)
44 : real, intent(out) :: flmu, am1(nrad), fmu1(nrad)
45 :
46 : integer :: k, k1, k2, m
47 0 : real :: fmu, a1, c1, fmu0, amamu(ipe)
48 :
49 : ! Get k1,get amamu(k)
50 0 : do k = 1, nel
51 0 : nkz(k) = -1
52 0 : do m = 1, ipe
53 0 : if (izz(k) == kz(m)) then
54 0 : amamu(k) = amass(m)
55 0 : nkz(k) = m
56 0 : exit
57 : end if
58 : end do
59 0 : if (nkz(k) == -1) then
60 0 : write(*,*) ' k=', k, ', izz(k)=', izz(k)
61 0 : write(*,*) ' kz(m) not found'
62 : end if
63 : end do
64 :
65 0 : k1 = -1
66 0 : do k2 = 1, kk
67 0 : do k = 1, nel
68 0 : if (izz(k) == iz1(k2)) then
69 0 : k1 = k
70 : end if
71 : end do
72 0 : kzz(k2) = k1
73 0 : am1(k2) = amamu(k1)
74 : end do
75 :
76 : ! Mean atomic weight = fmu
77 0 : fmu = 0.
78 0 : do k = 1, nel
79 0 : fmu = fmu + fa(k)*amamu(k)
80 : end do
81 :
82 0 : do k2 = 1, kk
83 0 : a1 = fa(kzz(k2))
84 0 : c1 = 1./(1.-a1)
85 0 : fmu0 = c1*(fmu - fa(kzz(k2))*amamu(kzz(k2)))
86 0 : fmu1(k2) = a1*(am1(k2) - fmu0)*1.660531d-24 !dmu/dlog xi
87 : end do
88 :
89 0 : fmu = fmu*1.660531d-24 ! Convert to cgs
90 0 : flmu = log10(dble(fmu))
91 0 : end subroutine abund
92 :
93 : ! *********************************************************************
94 0 : subroutine rd(i3, kk, kzz, nel, nkz, izz, ilab, jh, ntot, umesh, semesh, ff, rr, ta, fac)
95 : implicit none
96 : integer, intent(in) :: i3, kk, kzz(nrad), nel, nkz(ipe), izz(ipe), ilab(4), jh(4), ntot
97 : real, intent(in) :: umesh(:), semesh(:) ! (nptot)
98 : real(dp), intent(in) :: fac(nel)
99 : real, intent(out) :: ff(:, :, :, :) ! (nptot, ipe, 4, 4)
100 : real, intent(out) :: ta(:, :, :, :) ! (nptot, nrad, 4, 4)
101 : real, intent(out) :: rr(28, ipe, 4, 4)
102 :
103 : integer :: i, j, k, k2, l, n, m, itt, jnn, izp, ne1, ne2, ne, ib, ia
104 0 : real :: ya, yb, d, se, u, fion(-1:28) !, ff_temp(nptot), ta_temp(nptot)
105 :
106 : ! declare variables in common block, by default: real (a-h, o-z), integer (i-n)
107 : ! integer :: ite1, ite2, ite3, jn1, jn2, jne3, ntotp, nc, nf, int, ne1p, ne2p, np, kp1, kp2, kp3, npp, mx, nx
108 : ! real :: umin, umax, epatom, oplnck, fionp, yy1, yy2, yx
109 : ! common /atomdata/ ite1,ite2,ite3,jn1(91),jn2(91),jne3,umin,umax,ntotp, &
110 : ! nc,nf,int(17),epatom(17,91,25),oplnck(17,91,25),ne1p(17,91,25), &
111 : ! ne2p(17,91,25),fionp(-1:28,28,91,25),np(17,91,25),kp1(17,91,25), &
112 : ! kp2(17,91,25),kp3(17,91,25),npp(17,91,25),mx(33417000), &
113 : ! yy1(33417000),yy2(120000000),nx(19305000),yx(19305000)
114 : ! save /atomdata/
115 :
116 : include 'formats'
117 :
118 : ! i=temperature index
119 : ! j=density index
120 : ! k=frequency index
121 : ! n=element index
122 : ! Get:
123 : ! mono opacity cross-section ff(k,n,i,j)
124 : ! modified cross-section for selected element, ta(k,i,j)
125 : ! zet(i,j) for diffusion coefficient
126 :
127 : ! Initializations
128 : if (dbg) then
129 : write (*, 2) 'size(ff,dim=1)', size(ff, dim=1)
130 : write (*, 2) 'size(ff,dim=2)', size(ff, dim=2)
131 : write (*, 2) 'size(ff,dim=3)', size(ff, dim=3)
132 : write (*, 2) 'size(ff,dim=4)', size(ff, dim=4)
133 : write (*, 2) 'size(ta,dim=1)', size(ta, dim=1)
134 : write (*, 2) 'size(ta,dim=2)', size(ta, dim=2)
135 : write (*, 2) 'size(ta,dim=3)', size(ta, dim=3)
136 : write (*, 2) 'size(ta,dim=4)', size(ta, dim=4)
137 : write (*, 2) 'size(rr,dim=1)', size(rr, dim=1)
138 : write (*, 2) 'size(rr,dim=2)', size(rr, dim=2)
139 : write (*, 2) 'size(rr,dim=3)', size(rr, dim=3)
140 : write (*, 2) 'size(rr,dim=4)', size(rr, dim=4)
141 : write (*, 2) 'ipe', ipe
142 : write (*, *) 'rd: ff = 0'
143 : end if
144 :
145 0 : ff = 0.
146 : if (dbg) write (*, *) 'rd: rr = 0'
147 0 : rr = 0.
148 : if (dbg) write (*, *) 'rd: ta = 0'
149 0 : ta = 0.
150 :
151 : ! Start loop on i (temperature index)
152 0 : do i = 1, 4
153 : if (dbg) write (*, 2) 'rd: i', i
154 0 : itt = (ilab(i) - ite1)/2 + 1
155 : ! Read mono opacities
156 0 : do j = 1, 4
157 : if (dbg) write (*, 2) 'rd: j', j
158 0 : jnn = (jh(j)*i3 - jn1(itt))/2 + 1
159 0 : do n = 1, nel
160 : if (dbg) write (*, 2) 'rd: n', n
161 0 : izp = izz(n)
162 0 : ne1 = ne1p(nkz(n), itt, jnn)
163 0 : ne2 = ne2p(nkz(n), itt, jnn)
164 0 : do ne = ne1, ne2
165 0 : fion(ne) = fionp(ne, nkz(n), itt, jnn)
166 0 : if (ne <= min(ne2, izp - 2)) rr(izp - 1 - ne, n, i, j) = fion(ne)
167 : end do
168 :
169 0 : do k = 1, ntot
170 0 : ff(k, n, i, j) = yy2(k + kp2(nkz(n), itt, jnn))
171 : end do
172 :
173 0 : if (fac(n) /= 1d0) then
174 0 : do k = 1, size(ff, dim=1)
175 0 : ff(k, n, i, j) = fac(n)*ff(k, n, i, j)
176 : end do
177 : end if
178 :
179 0 : do k2 = 1, kk
180 : if (dbg) write (*, 2) 'rd: k2', k2
181 0 : if (kzz(k2) == n) then
182 0 : ib = 1
183 0 : yb = yx(1 + kp3(nkz(kzz(k2)), itt, jnn))
184 0 : u = umesh(ib)
185 0 : se = semesh(ib)
186 0 : ta(1, k2, i, j) = se*ff(1, n, i, j) - yb
187 0 : do m = 2, npp(nkz(kzz(k2)), itt, jnn)
188 0 : ia = ib
189 0 : ya = yb
190 0 : ib = nx(m + kp3(nkz(kzz(k2)), itt, jnn))
191 0 : yb = yx(m + kp3(nkz(kzz(k2)), itt, jnn))
192 0 : d = (yb - ya)/float(ib - ia)
193 0 : do l = ia + 1, ib - 1
194 0 : u = umesh(l)
195 0 : se = semesh(l)
196 0 : ta(l, k2, i, j) = se*ff(l, n, i, j) - (ya + (l - ia)*d)
197 : end do
198 0 : u = umesh(ib)
199 0 : se = semesh(ib)
200 0 : ta(ib, k2, i, j) = se*ff(ib, n, i, j) - yb
201 : end do
202 : exit ! get out of k2-loop
203 : end if
204 : end do !k2
205 : end do !n
206 : end do !j
207 : end do !i
208 0 : end subroutine rd
209 :
210 0 : subroutine mix(kk, kzz, ntot, nel, fa, ff, rr, rs, rion)
211 : implicit none
212 : integer, intent(in) :: kk, kzz(nrad), ntot, nel
213 : real, intent(in) :: fa(ipe), rr(28, 17, 4, 4)
214 : real, intent(in) :: ff(:, :, :, :) ! (nptot, ipe, 4, 4)
215 : real, intent(out) :: rs(:, :, :) ! (nptot, 4, 4)
216 : real, intent(out) :: rion(28, 4, 4)
217 :
218 : integer :: i, j, k, n, m, k2
219 :
220 0 : do i = 1, 4
221 0 : do j = 1, 4
222 0 : do n = 1, ntot
223 0 : rs(n, i, j) = dot_product(ff(n, 1:nel, i, j), fa(1:nel))
224 : end do
225 :
226 0 : do m = 1, 28
227 0 : rion(m, i, j) = dot_product(rr(m, 1:nel, i, j), fa(1:nel))
228 : end do
229 : end do
230 : end do
231 0 : end subroutine mix
232 :
233 0 : subroutine ross(kk, flmu, fmu1, dv, ntot, rs, rossl, gaml, ta)
234 : implicit none
235 : integer, intent(in) :: kk, ntot
236 : real, intent(in) :: rs(:, :, :) ! (nptot, 4, 4)
237 : real, intent(in) :: ta(:, :, :, :) ! (nptot, nrad, 4, 4)
238 : real, intent(in) :: flmu, dv, fmu1(nrad)
239 : real, intent(out) :: rossl(4, 4), gaml(4, 4, nrad)
240 :
241 : integer :: k2, i, j, n
242 0 : real(dp) :: drs, dd, dgm(nrad)
243 0 : real :: fmu, oross, tt, exp10_flmu
244 :
245 : ! oross=cross-section in a.u.
246 : ! rossl=log10(ROSS in cgs)
247 0 : exp10_flmu = real(exp10(dble(flmu)))
248 0 : do i = 1, 4
249 0 : do j = 1, 4
250 0 : drs = 0.d0
251 0 : dgm(:) = 0.d0
252 0 : do n = 1, ntot !10000
253 0 : dd = 1.d0/rs(n, i, j)
254 0 : drs = drs + dd
255 0 : do k2 = 1, kk
256 0 : tt = ta(n, k2, i, j)
257 0 : dgm(k2) = dgm(k2) + tt*dd
258 : end do
259 : end do
260 0 : oross = 1./(drs*dv)
261 0 : rossl(i, j) = log10(dble(oross)) - 16.55280 - flmu
262 0 : do k2 = 1, kk
263 0 : if (dgm(k2) > 0) then
264 0 : dgm(k2) = dgm(k2)*dv
265 0 : gaml(i, j, k2) = log10(dgm(k2))
266 : else
267 0 : gaml(i, j, k2) = -30.
268 : end if
269 : end do
270 : end do !j
271 : end do !i
272 0 : end subroutine ross
273 :
274 0 : subroutine interp(nel, kk, rossl, gaml, xi, eta, g, i3, f)
275 : use op_common, only: fint, fintp
276 : implicit none
277 : integer, intent(in) :: nel, kk, i3
278 : real, intent(in) :: gaml(4, 4, nrad), rossl(4, 4), eta, xi
279 : real, intent(out) :: f(nrad), g
280 : integer :: l, i, j, k2
281 0 : real :: v(4), u(4), vyi(4)
282 :
283 : ! interpolation of g (=rosseland mean opacity)
284 0 : do i = 1, 4
285 0 : do j = 1, 4
286 0 : u(j) = rossl(i, j)
287 : end do
288 0 : v(i) = fint(u, eta)
289 0 : vyi(i) = fintp(u, eta)
290 : end do
291 0 : g = fint(v, xi)
292 :
293 0 : do k2 = 1, kk
294 : ! Interpolation of gam (=dimensionless parameter in g_rad)
295 : ! f=log10(gamma)
296 0 : do i = 1, 4
297 0 : do j = 1, 4
298 : ! HH: This gives irregularities, perhaps it is preferable to assign nonnegative values of
299 : ! neighbouring interpolation points (in subroutine "ross") ?
300 0 : u(j) = gaml(i, j, k2)
301 : end do
302 0 : v(i) = fint(u, eta)
303 0 : vyi(i) = fintp(u, eta)
304 : end do
305 :
306 0 : f(k2) = fint(v, xi)
307 : end do
308 0 : end subroutine interp
309 :
310 : end module op_radacc
|
