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 : ! FORTRAN 90 module for calculation of radiative accelerations,
21 : ! based on the Opacity Project (OP) code "OPserver".
22 : ! See CHANGES_HU for changes made to the original code.
23 :
24 : ! Haili Hu 2010
25 :
26 : module op_eval
27 :
28 : use op_load
29 : use const_def, only: dp
30 : use math_lib
31 : use kap_def, only: kap_test_partials, kap_test_partials_val, kap_test_partials_dval_dx
32 :
33 : implicit none
34 :
35 : private
36 : public :: eval_op_radacc
37 : public :: eval_op_ev
38 : public :: eval_alt_op
39 :
40 : logical, parameter :: dbg = .false.
41 :
42 : contains
43 :
44 : ! HH: Based on "op_ax.f"
45 : ! Input: kk = number of elements to calculate g_rad for
46 : ! iz1(kk) = charge of element to calculate g_rad for
47 : ! nel = number of elements in mixture
48 : ! izzp(nel) = charge of elements
49 : ! fap(nel) = number fractions of elements
50 : ! fac(nel) = scale factors for element opacity
51 : ! flux = local radiative flux (Lrad/4*pi*r^2)
52 : ! fltp = log T
53 : ! flrhop = log rho
54 : ! screening if true, use screening corrections
55 : ! Output: g1 = log kappa
56 : ! gx1 = d(log kappa)/d(log T)
57 : ! gy1 = d(log kappa)/d(log rho)
58 : ! gp1(kk) = d(log kappa)/d(log xi)
59 : ! grl1(kk) = log grad
60 : ! fx1(kk) = d(log grad)/d(log T)
61 : ! fy1(kk) = d(log grad)/d(log rho)
62 : ! grlp1(kk) = d(log grad)/d(log xi)
63 : ! meanZ(nel) = average ionic charge of elements
64 : ! zetx1(nel) = d(meanZ)/d(log T)
65 : ! zety1(nel) = d(meanZ)/d(log rho)
66 : ! ierr = 0 for correct use
67 0 : subroutine eval_op_radacc(kk, izk, nel, izzp, fap, fac, flux, fltp, flrhop, screening, g1, grl1, umesh, semesh, ff, ta, rs, ierr)
68 : use op_radacc
69 : use op_load, only: msh
70 : use op_common
71 : integer, intent(in) :: kk, nel
72 : integer, intent(in) :: izk(kk), izzp(nel)
73 : real(dp), intent(in) :: fap(nel), fac(nel)
74 : real(dp), intent(in) :: flux, fltp, flrhop
75 : logical, intent(in) :: screening
76 : real(dp), intent(out) :: g1
77 : real(dp), intent(inout) :: grl1(kk)
78 : real, pointer :: umesh(:), semesh(:), ff(:, :, :, :), ta(:, :, :, :), rs(:, :, :)
79 : ! umesh(nptot)
80 : ! semesh(nptot)
81 : ! ff(nptot, ipe, 4, 4)
82 : ! ta(nptot, nrad, 4, 4),
83 : ! rs(nptot, 4, 4)
84 : integer, intent(out) :: ierr
85 : ! local variables
86 : integer :: n, i, k2, i3, ntot, jhmin, jhmax
87 : integer :: ih(4), jh(4), ilab(4), kzz(nrad), nkz(ipe), izz(ipe), iz1(nrad)
88 0 : real :: const, gx, gy, flt, flrho, flmu, dscat, dv, xi, flne, epa, eta, ux, uy, g
89 0 : real :: uf(0:100), rion(28, 4, 4), rossl(4, 4), flr(4, 4), rr(28, ipe, 4, 4), &
90 0 : fa(ipe), gaml(4, 4, nrad), f(nrad), am1(nrad), fmu1(nrad)
91 :
92 : include 'formats'
93 :
94 0 : ierr = 0
95 0 : if (nel <= 0 .or. nel > ipe) then
96 0 : write (6, *) 'OP - NUMBER OF ELEMENTS OUT OF RANGE:', nel
97 0 : ierr = 1
98 0 : return
99 : end if
100 :
101 : if (dbg) write (*, *) 'start eval_op_radacc'
102 :
103 : ! Get i3 for mesh type q='m'
104 0 : i3 = 2
105 :
106 : ! HH: k2 loops over elements for which to calculate grad.
107 0 : do k2 = 1, kk
108 0 : do n = 1, ipe
109 0 : if (izk(k2) == kz(n)) then
110 0 : iz1(k2) = izk(k2)
111 : exit
112 : end if
113 0 : if (n == ipe) then
114 0 : write (6, *) 'OP - SELECTED ELEMENT CANNOT BE TREATED: Z = ', izk(k2)
115 0 : write (6, *) 'OP supports C, N, O, Ne, Na, Mg, Al, Si, S, Ar, Ca, Cr, Mn, Fe, and Ni'
116 0 : ierr = 5
117 0 : return
118 : end if
119 : end do
120 : end do
121 :
122 0 : outer: do i = 1, nel
123 0 : inner: do n = 1, ipe
124 0 : if (izzp(i) == kz(n)) then
125 0 : izz(i) = izzp(i)
126 0 : fa(i) = fap(i)
127 0 : if (fa(i) < 0.0) then
128 0 : write (6, *) 'OP - NEGATIVE FRACTIONAL ABUNDANCE:', fa(i)
129 0 : ierr = 7
130 0 : return
131 : end if
132 : cycle outer
133 : end if
134 : end do inner
135 0 : write (6, *) 'OP - CHEM. ELEMENT CANNOT BE INCLUDED: Z = ', izzp(i)
136 0 : ierr = 8
137 0 : return
138 : end do outer
139 :
140 : ! Calculate mean atomic weight (flmu) and
141 : ! array kzz indicating elements for which to calculate g_rad
142 : if (dbg) write (*, *) 'call abund'
143 0 : call abund(nel, izz, kk, iz1, fa, kzz, flmu, am1, fmu1, nkz)
144 :
145 : ! Other initializations
146 : ! dv = interval in frequency variable v
147 : ! ntot=number of frequency points
148 : ! umesh, values of u=(h*nu/k*T) on mesh points
149 : ! semesh, values of 1-exp(-u) on mesh points
150 : ! uf, dscat used in scattering correction
151 : if (dbg) write (*, *) 'call msh'
152 0 : call msh(dv, ntot, umesh, semesh, uf, dscat) !output variables
153 :
154 : ! Start loop on temperature-density points
155 : ! flt=log10(T, K)
156 : ! flrho=log10(rho, cgs)
157 :
158 0 : flt = fltp
159 0 : flrho = flrhop
160 : ! Get temperature indices
161 : ! Let ite(i) be temperature index used in mono files
162 : ! Put ite(i)=2*ih(i)
163 : ! Use ih(i), i=1 to 4
164 : ! xi=interpolation variable
165 : ! log10(T)=flt=0.025*(ite(1)+xi+3)
166 : ! ilab(i) is temperature label
167 : if (dbg) write (*, *) 'call xindex'
168 0 : call xindex(flt, ilab, xi, ih, i3, ierr)
169 0 : if (ierr /= 0) then
170 0 : write (*, *) "xindex errored in radacc"
171 0 : return
172 : end if
173 :
174 : ! Get density indices
175 : ! Let jne(j) be density index used in mono files
176 : ! Put jne(j)=2*jh(j)
177 : ! Use jh(j), j=1 to 4
178 : ! Get extreme range for jh
179 : if (dbg) write (*, *) 'call jrange'
180 0 : call jrange(ih, jhmin, jhmax, i3)
181 :
182 : ! Get electron density flne=log10(Ne) for specified mass density flrho
183 : ! Also: UY=0.25*[d log10(rho)]/[d log10(Ne)]
184 : ! epa=electrons per atom
185 : if (dbg) write (*, *) 'call findne'
186 0 : call findne(ilab, fa, nel, nkz, jhmin, jhmax, ih, flrho, flt, xi, flne, flmu, flr, epa, uy, i3, ierr)
187 0 : if (ierr /= 0) then
188 0 : write (*, *) "findne encountered error in radacc"
189 0 : return
190 : end if
191 :
192 : ! Get density indices jh(j), j=1 to 4,
193 : ! Interpolation variable eta
194 : ! log10(Ne)=flne=0.25*(jne(1)+eta+3)
195 : if (dbg) write (*, *) 'call yindex'
196 0 : call yindex(jhmin, jhmax, flne, jh, i3, eta)
197 :
198 : ! Get ux=0.025*[d log10(rho)]/[d log10(T)]
199 : if (dbg) write (*, *) 'call findux'
200 0 : call findux(flr, xi, eta, ux)
201 :
202 : ! rossl(i,j)=log10(Rosseland mean) on mesh points (i,j)
203 : ! Get new mono opacities, ff(n,k,i,j)
204 : if (dbg) write (*, *) 'call rd'
205 0 : call rd(i3, kk, kzz, nel, nkz, izz, ilab, jh, ntot, umesh, semesh, ff, rr, ta, fac)
206 :
207 : ! Get rs = weighted sum of monochromatic opacity cross sections
208 : if (dbg) write (*, *) 'call mix'
209 0 : call mix(kk, kzz, ntot, nel, fa, ff, rr, rs, rion)
210 :
211 : ! Screening corrections
212 0 : if (screening) then
213 : ! Get Boercker scattering correction
214 : if (dbg) write (*, *) 'call scatt'
215 0 : call scatt(ih, jh, rion, uf, rs, umesh, semesh, dscat, ntot, epa, ierr)
216 0 : if (ierr /= 0) return
217 : ! Get correction for Debye screening
218 : if (dbg) write (*, *) 'call screen1'
219 0 : call screen1(ih, jh, rion, umesh, ntot, epa, rs)
220 : end if
221 :
222 : ! Get rossl, array of log10(Rosseland mean in cgs)
223 : if (dbg) write (*, *) 'call ross'
224 0 : call ross(kk, flmu, fmu1, dv, ntot, rs, rossl, gaml, ta)
225 :
226 : ! Interpolate to required flt, flrho
227 : ! g=log10(ross, cgs)
228 : if (dbg) write (*, *) 'call interp'
229 0 : call interp(nel, kk, rossl, gaml, xi, eta, g, i3, f)
230 : if (dbg) write (*, *) 'done interp'
231 :
232 : ! Write grad in terms of local radiative flux instead of (Teff, r/R*):
233 0 : const = 13.30295d0 + log10(flux) ! = -log10(c) - log10(amu) + log10(flux)
234 0 : do k2 = 1, kk
235 0 : grl1(k2) = const + flmu - log10(dble(am1(k2))) + f(k2) + g ! log g_rad
236 : end do
237 :
238 0 : g1 = g ! log kappa
239 :
240 : if (dbg) write (*, *) 'done eval_op_radacc'
241 :
242 0 : end subroutine eval_op_radacc
243 :
244 : ! ***********************************************************************
245 : ! HH: Based on "op_mx.f", opacity calculations to be used for stellar evolution calculations
246 : ! Input: nel = number of elements in mixture
247 : ! izzp(nel) = charge of elements
248 : ! fap(nel) = number fractions of elements
249 : ! fac(nel) = scale factors for element opacity
250 : ! fltp = log (temperature)
251 : ! flrhop = log (mass density)
252 : ! screening if true, use screening corrections
253 : ! Output: g1 = log kappa
254 : ! gx1 = d(log kappa)/d(log T)
255 : ! gy1 = d(log kappa)/d(log rho)
256 : ! ierr = 0 for correct use
257 0 : subroutine eval_op_ev(nel, izzp, fap, fac, fltp, flrhop, screening, g1, gx1, gy1, umesh, semesh, ff, rs, ierr)
258 0 : use op_ev
259 : use op_load, only: msh
260 : use op_common
261 : integer, intent(in) :: nel
262 : integer, intent(in) :: izzp(nel)
263 : real(dp), intent(in) :: fap(nel), fac(nel)
264 : real(dp), intent(in) :: fltp, flrhop
265 : logical, intent(in) :: screening
266 : real(dp), intent(inout) :: g1, gx1, gy1
267 : real, pointer :: umesh(:), semesh(:), ff(:, :, :, :), rs(:, :, :)
268 : ! umesh(nptot)
269 : ! semesh(nptot)
270 : ! ff(nptot, ipe, 4, 4)
271 : ! rs(nptot, 4, 4)
272 : ! s(nptot, nrad, 4, 4)
273 : integer, intent(out) :: ierr
274 : ! local variables
275 : integer :: n, i, i3, jhmin, jhmax, ntot
276 : integer :: ih(4), jh(4), ilab(4), izz(ipe), nkz(ipe)
277 0 : real :: flt, flrho, flmu, flne, dv, dscat, const, gx, gy, g, eta, epa, xi, ux, uy
278 0 : real :: uf(0:100), rion(28, 1:4, 1:4), rossl(4, 4), flr(4, 4), fa(ipe), rr(28, ipe, 4, 4), fmu1(nrad)
279 :
280 : ! Initializations
281 0 : ierr = 0
282 0 : if (nel <= 0 .or. nel > ipe) then
283 0 : write (6, *) 'OP - NUMBER OF ELEMENTS OUT OF RANGE:', nel
284 0 : ierr = 1
285 0 : return
286 : end if
287 :
288 : ! Get i3 for mesh type q='m'
289 0 : i3 = 2
290 :
291 0 : outer: do i = 1, nel
292 0 : inner: do n = 1, ipe
293 0 : if (izzp(i) == kz(n)) then
294 0 : izz(i) = izzp(i)
295 0 : fa(i) = fap(i)
296 0 : if (fa(i) < 0.d0) then
297 0 : write (6, *) 'OP - NEGATIVE FRACTIONAL ABUNDANCE:', fa(i)
298 0 : ierr = 7
299 0 : return
300 : end if
301 : cycle outer
302 : end if
303 : end do inner
304 0 : write (6, *) 'OP - CHEM. ELEMENT CANNOT BE INCLUDED: Z = ', izzp(i)
305 0 : ierr = 8
306 0 : return
307 : end do outer
308 :
309 : ! Calculate mean atomic weight (flmu)
310 0 : call abund(nel, izz, fa, flmu, fmu1, nkz)
311 :
312 : ! Other initializations
313 : ! dv = interval in frequency variable v
314 : ! ntot=number of frequency points
315 : ! umesh, values of u=(h*nu/k*T) on mesh points
316 : ! semesh, values of 1-exp(-u) on mesh points
317 : ! uf, dscat used in scattering correction
318 0 : call msh(dv, ntot, umesh, semesh, uf, dscat)
319 :
320 : ! Start loop on temperature-density points
321 : ! flt=log10(T, K)
322 : ! flrho=log10(rho, cgs)
323 0 : flt = fltp
324 0 : flrho = flrhop
325 : ! Get temperature indices
326 : ! Let ite(i) be temperature index used in mono files
327 : ! Put ite(i)=2*ih(i)
328 : ! Use ih(i), i=1 to 4
329 : ! xi=interpolation variable
330 : ! log10(T)=flt=0.025*(ite(1)+xi+3)
331 : ! ilab(i) is temperature label
332 0 : call xindex(flt, ilab, xi, ih, i3, ierr)
333 0 : if (ierr /= 0) then
334 0 : write (*, *) "eval_op_ev failed in xindex"
335 0 : return
336 : end if
337 :
338 : ! Get density indices
339 : ! Let jne(j) be density index used in mono files
340 : ! Put jne(j)=2*jh(j)
341 : ! Use jh(j), j=1 to 4
342 : ! Get extreme range for jh
343 0 : call jrange(ih, jhmin, jhmax, i3)
344 :
345 : ! Get electron density flne=log10(Ne) for specified mass density flrho
346 : ! Also: UY=0.25*[d log10(rho)]/[d log10(Ne)]
347 : ! epa=electrons per atom
348 0 : call findne(ilab, fa, nel, nkz, jhmin, jhmax, ih, flrho, flt, xi, flne, flmu, flr, epa, uy, i3, ierr)
349 0 : if (ierr /= 0) then
350 0 : write (*, *) "eval_op_ev failed in findne"
351 0 : return
352 : end if
353 :
354 : ! Get density indices jh(j), j=1 to 4,
355 : ! Interpolation variable eta
356 : ! log10(Ne)=flne=0.25*(jne(1)+eta+3)
357 0 : call yindex(jhmin, jhmax, flne, jh, i3, eta)
358 :
359 : ! Get ux=0.025*[d log10(rho)]/[d log10(T)]
360 0 : call findux(flr, xi, eta, ux)
361 :
362 : ! rossl(i,j)=log10(Rosseland mean) on mesh points (i,j)
363 : ! Get new mono opacities, ff(n,k,i,j)
364 0 : call rd(nel, nkz, izz, ilab, jh, ntot, ff, rr, i3, umesh, fac)
365 :
366 : ! Up-date mixture
367 0 : call mix(ntot, nel, fa, ff, rs, rr, rion)
368 :
369 0 : if (screening) then
370 : ! Get Boercker scattering correction
371 0 : call scatt(ih, jh, rion, uf, rs, umesh, semesh, dscat, ntot, epa, ierr)
372 0 : if (ierr /= 0) then
373 0 : write (*, *) "scattering correction failed in op_eval_ev"
374 0 : return
375 : end if
376 : ! Get correction for Debye screening
377 0 : call screen1(ih, jh, rion, umesh, ntot, epa, rs)
378 : end if
379 :
380 : ! Get rossl, array of log10(Rosseland mean in cgs)
381 0 : call ross(flmu, fmu1, dv, ntot, rs, rossl)
382 :
383 : ! Interpolate to required flt, flrho
384 : ! g=log10(ross, cgs)
385 0 : call interp(nel, rossl, xi, eta, g, i3, ux, uy, gx, gy)
386 :
387 0 : g1 = g ! log kappa
388 0 : gx1 = gx ! dlogkappa/dt
389 0 : gy1 = gy ! dlogkappa/drho
390 :
391 0 : end subroutine eval_op_ev
392 :
393 : ! ***********************************************************************
394 :
395 : ! HH: Based on "op_mx.f", opacity calculations to be used for non-adiabatic pulsation calculations
396 : ! Special care is taken to ensure smoothness of opacity derivatives
397 : ! Input: nel = number of elements in mixture
398 : ! izzp(nel) = charge of elements
399 : ! fap(nel) = number fractions of elements
400 : ! fac(nel) = scale factors for element opacity
401 : ! fltp = log (temperature)
402 : ! flrhop = log (mass density)
403 : ! screening if true, use screening corrections
404 : ! Output: g1 = log kappa
405 : ! gx1 = d(log kappa)/d(log T)
406 : ! gy1 = d(log kappa)/d(log rho)
407 : ! ierr = 0 for correct use
408 0 : subroutine eval_alt_op(nel, izzp, fap, fac, fltp, flrhop, screening, g1, gx1, gy1, umesh, semesh, ff, rs, ierr)
409 0 : use op_osc
410 : use op_load, only: msh
411 : integer, intent(in) :: nel
412 : integer, intent(in) :: izzp(nel)
413 : real(dp), intent(in) :: fap(nel), fac(nel)
414 : real(dp), intent(in) :: fltp, flrhop
415 : logical, intent(in) :: screening
416 : real(dp), intent(out) :: g1, gx1, gy1
417 : ! real(dp), intent(out) :: meanZ(nel)
418 : real, pointer :: umesh(:), semesh(:), ff(:, :, :, :), rs(:, :, :)
419 : ! umesh(nptot)
420 : ! semesh(nptot)
421 : ! ff(nptot, ipe, 0:5, 0:5)
422 : ! rs(nptot, 0:5, 0:5)
423 : integer, intent(out) :: ierr
424 : ! local variables
425 : integer :: n, i, i3, jhmin, jhmax, ntot
426 : integer :: ih(0:5), jh(0:5), ilab(0:5), izz(ipe), nkz(ipe)
427 0 : real :: flt, flrho, flmu, flne, dv, dscat, const, gx, gy, g, eta, epa, xi, ux, uy
428 0 : real :: uf(0:100), rion(28, 0:5, 0:5), rossl(0:5, 0:5), flr(4, 4), fa(ipe), rr(28, ipe, 0:5, 0:5)
429 :
430 : ! Initializations
431 0 : ierr = 0
432 0 : if (nel <= 0 .or. nel > ipe) then
433 0 : write (6, *) 'OP - NUMBER OF ELEMENTS OUT OF RANGE:', nel
434 0 : ierr = 1
435 0 : return
436 : end if
437 :
438 : ! Get i3 for mesh type q='m'
439 0 : i3 = 2
440 :
441 0 : outer: do i = 1, nel
442 0 : inner: do n = 1, ipe
443 0 : if (izzp(i) == kz(n)) then
444 0 : izz(i) = izzp(i)
445 0 : fa(i) = fap(i)
446 0 : if (fa(i) < 0.0) then
447 0 : write (6, *) 'OP - NEGATIVE FRACTIONAL ABUNDANCE:', fa(i)
448 0 : ierr = 7
449 0 : return
450 : end if
451 : cycle outer
452 : end if
453 : end do inner
454 0 : write (6, *) 'OP - CHEM. ELEMENT CANNOT BE INCLUDED: Z = ', izzp(i)
455 0 : ierr = 8
456 0 : return
457 : end do outer
458 :
459 : ! Calculate mean atomic weight (flmu)
460 0 : call abund(nel, izz, fa, flmu, nkz)
461 :
462 : ! Other initializations
463 : ! dv = interval in frequency variable v
464 : ! ntot=number of frequency points
465 : ! umesh, values of u=(h*nu/k*T) on mesh points
466 : ! semesh, values of 1-exp(-u) on mesh points
467 : ! uf, dscat used in scattering correction
468 0 : call msh(dv, ntot, umesh, semesh, uf, dscat)
469 :
470 : ! Start loop on temperature-density points
471 : ! flt=log10(T, K)
472 : ! flrho=log10(rho, cgs)
473 0 : flt = fltp
474 0 : flrho = flrhop
475 : ! Get temperature indices
476 : ! Let ite(i) be temperature index used in mono files
477 : ! Put ite(i)=2*ih(i)
478 : ! Use ih(i), i=1 to 4
479 : ! xi=interpolation variable
480 : ! log10(T)=flt=0.025*(ite(1)+xi+3)
481 : ! ilab(i) is temperature label
482 0 : call xindex(flt, ilab, xi, ih, i3, ierr)
483 0 : if (ierr /= 0) return
484 :
485 : ! Get density indices
486 : ! Let jne(j) be density index used in mono files
487 : ! Put jne(j)=2*jh(j)
488 : ! Use jh(j), j=1 to 4
489 : ! Get extreme range for jh
490 0 : call jrange(ih, jhmin, jhmax, i3)
491 :
492 : ! Get electron density flne=log10(Ne) for specified mass density flrho
493 : ! Also: UY=0.25*[d log10(rho)]/[d log10(Ne)]
494 : ! epa=electrons per atom
495 0 : call findne(ilab, fa, nel, nkz, jhmin, jhmax, ih, flrho, flt, xi, flne, flmu, flr, epa, uy, i3, ierr)
496 0 : if (ierr /= 0) return
497 :
498 : ! Get density indices jh(j), j=1 to 4,
499 : ! Interpolation variable eta
500 : ! log10(Ne)=flne=0.25*(jne(1)+eta+3)
501 0 : call yindex(jhmin, jhmax, flne, jh, i3, eta)
502 :
503 : ! Get ux=0.025*[d log10(rho)]/[d log10(T)]
504 0 : call findux(flr, xi, eta, ux)
505 :
506 : ! rossl(i,j)=log10(Rosseland mean) on mesh points (i,j)
507 : ! Get new mono opacities, ff(n,k,i,j)
508 0 : call rd(nel, nkz, izz, ilab, jh, ntot, ff, rr, i3, umesh, fac)
509 :
510 : ! Up-date mixture
511 0 : call mix(ntot, nel, fa, ff, rs, rr, rion)
512 :
513 0 : if (screening) then
514 : ! Get Boercker scattering correction
515 0 : call scatt(ih, jh, rion, uf, rs, umesh, semesh, dscat, ntot, epa, ierr)
516 0 : if (ierr /= 0) return
517 : ! Get correction for Debye screening
518 0 : call screen1(ih, jh, rion, umesh, ntot, epa, rs)
519 : end if
520 :
521 : ! Get rossl, array of log10(Rosseland mean in cgs)
522 0 : call ross(flmu, dv, ntot, rs, rossl)
523 :
524 : ! Interpolate to required flt, flrho
525 : ! g=log10(ross, cgs)
526 0 : call interp(nel, rossl, xi, eta, g, i3, ux, uy, gx, gy)
527 :
528 0 : g1 = g ! log kappa
529 0 : gx1 = gx ! dlogkappa/dt
530 0 : gy1 = gy ! dlogkappa/drho
531 :
532 0 : end subroutine eval_alt_op
533 :
534 : end module op_eval
|
