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 :
21 : module kap_lib
22 : ! library for calculating opacities
23 :
24 : ! the data interface for the library is defined in kap_def
25 :
26 : use const_def, only: dp
27 : use math_lib
28 :
29 : implicit none
30 :
31 : integer , pointer, public, save :: izz(:),ite(:),jne(:)
32 : real(dp), pointer, public, save :: sig(:,:,:)
33 : real(dp), pointer, public, save :: epatom(:,:),amamu(:),eumesh(:,:,:)
34 : real(dp), pointer, public, save :: lkap_ross_pcg(:)
35 : integer, public, parameter :: ngp = 2
36 : real(dp), public, save :: lgamm_pcg(ngp,17,1648), lkap_face_pcg(ngp,1648), logT_pcg(ngp,1648), logRho_pcg(ngp,1648)
37 : integer, parameter :: nzm = 3000
38 : real(dp), public, save :: fk_old(nzm,17), fk_old_grad(nzm,17)
39 : real(dp), public, save ::logT_grid_old(nzm,4,4), logRho_grid_old(nzm,4,4), fk_pcg(17), fk_grad_pcg(ngp,17)
40 : real(dp), public, save :: lkap_grid_old(nzm,4,4), logT_cntr_old(nzm), logRho_cntr_old(nzm)
41 : logical :: initialize_fk_old = .true.
42 : contains ! the procedure interface for the library
43 : ! client programs should only call these routines.
44 :
45 :
46 : ! call this routine to initialize the kap module.
47 : ! only needs to be done once at start of run.
48 : ! Reads data from the 'kap' directory in the data_dir.
49 : ! If use_cache is true and there is a 'kap/cache' directory, it will try that first.
50 : ! If it doesn't find what it needs in the cache,
51 : ! it reads the data and writes the cache for next time.
52 5 : subroutine kap_init(use_cache, kap_cache_dir, ierr)
53 : use kap_def, only : kap_def_init, kap_use_cache, kap_is_initialized
54 : logical, intent(in) :: use_cache
55 : character (len=*), intent(in) :: kap_cache_dir ! blank means use default
56 : integer, intent(out) :: ierr ! 0 means AOK.
57 5 : ierr = 0
58 5 : if (kap_is_initialized) return
59 5 : call kap_def_init(kap_cache_dir)
60 5 : kap_use_cache = use_cache
61 5 : kap_is_initialized = .true.
62 : end subroutine kap_init
63 :
64 :
65 3 : subroutine kap_shutdown
66 : use kap_def, only: do_Free_Kap_Tables, kap_is_initialized
67 3 : call do_Free_Kap_Tables()
68 3 : kap_is_initialized = .false.
69 3 : end subroutine kap_shutdown
70 :
71 :
72 : ! after kap_init has finished, you can allocate a "handle".
73 :
74 2 : integer function alloc_kap_handle(ierr) result(handle)
75 : integer, intent(out) :: ierr ! 0 means AOK.
76 : character (len=0) :: inlist
77 2 : handle = alloc_kap_handle_using_inlist(inlist, ierr)
78 2 : end function alloc_kap_handle
79 :
80 9 : integer function alloc_kap_handle_using_inlist(inlist,ierr) result(handle)
81 : use kap_def, only:do_alloc_kap,kap_is_initialized
82 : use kap_ctrls_io, only:read_namelist
83 : character (len=*), intent(in) :: inlist ! empty means just use defaults.
84 : integer, intent(out) :: ierr ! 0 means AOK.
85 : ierr = 0
86 9 : if (.not. kap_is_initialized) then
87 0 : ierr=-1
88 0 : return
89 : end if
90 9 : handle = do_alloc_kap(ierr)
91 9 : if (ierr /= 0) return
92 9 : call read_namelist(handle, inlist, ierr)
93 9 : if (ierr /= 0) return
94 9 : call kap_setup_tables(handle, ierr)
95 9 : call kap_setup_hooks(handle, ierr)
96 18 : end function alloc_kap_handle_using_inlist
97 :
98 1 : subroutine free_kap_handle(handle)
99 : ! frees the handle and all associated data
100 9 : use kap_def,only: Kap_General_Info,do_free_kap
101 : integer, intent(in) :: handle
102 1 : call do_free_kap(handle)
103 1 : end subroutine free_kap_handle
104 :
105 :
106 172482 : subroutine kap_ptr(handle,rq,ierr)
107 : use kap_def,only:Kap_General_Info,get_kap_ptr,kap_is_initialized
108 : integer, intent(in) :: handle ! from alloc_kap_handle
109 : type (Kap_General_Info), pointer :: rq
110 : integer, intent(out):: ierr
111 86241 : if (.not. kap_is_initialized) then
112 0 : ierr=-1
113 0 : return
114 : end if
115 86241 : call get_kap_ptr(handle,rq,ierr)
116 : end subroutine kap_ptr
117 :
118 :
119 11 : subroutine kap_setup_tables(handle, ierr)
120 : use kap_def, only : Kap_General_Info, get_kap_ptr
121 : use load_kap, only : Setup_Kap_Tables
122 : integer, intent(in) :: handle
123 : integer, intent(out):: ierr
124 :
125 : type (Kap_General_Info), pointer :: rq
126 : logical, parameter :: use_cache = .true.
127 : logical, parameter :: load_on_demand = .true.
128 :
129 : ierr = 0
130 11 : call get_kap_ptr(handle,rq,ierr)
131 11 : call Setup_Kap_Tables(rq, use_cache, load_on_demand, ierr)
132 :
133 11 : end subroutine kap_setup_tables
134 :
135 :
136 9 : subroutine kap_setup_hooks(handle, ierr)
137 11 : use kap_def, only : Kap_General_Info, get_kap_ptr
138 : use other_elect_cond_opacity
139 : use other_compton_opacity
140 : integer, intent(in) :: handle
141 : integer, intent(out):: ierr
142 :
143 : type (Kap_General_Info), pointer :: rq
144 :
145 : ierr = 0
146 9 : call get_kap_ptr(handle,rq,ierr)
147 :
148 9 : rq% other_elect_cond_opacity => null_other_elect_cond_opacity
149 9 : rq% other_compton_opacity => null_other_compton_opacity
150 :
151 9 : end subroutine kap_setup_hooks
152 :
153 :
154 : ! kap evaluation
155 : ! you can call these routines after you've setup the tables for the handle.
156 : ! NOTE: the structures referenced via the handle are read-only
157 : ! for the evaluation routines, so you can do multiple evaluations in parallel
158 : ! using the same handle.
159 :
160 :
161 86216 : subroutine kap_get( &
162 86216 : handle, species, chem_id, net_iso, xa, &
163 : logRho, logT, &
164 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
165 : eta, d_eta_dlnRho, d_eta_dlnT , &
166 86216 : kap_fracs, kap, dlnkap_dlnRho, dlnkap_dlnT, dlnkap_dxa, ierr)
167 : use chem_def, only: chem_isos
168 : use chem_lib, only: basic_composition_info
169 : use kap_def, only : kap_is_initialized, Kap_General_Info, num_kap_fracs, i_frac_Type2
170 : use kap_eval, only : Get_kap_Results
171 : ! INPUT
172 : integer, intent(in) :: handle ! from alloc_kap_handle; in star, pass s% kap_handle
173 : integer, intent(in) :: species
174 : integer, pointer :: chem_id(:) ! maps species to chem id
175 : integer, pointer :: net_iso(:) ! maps chem id to species number
176 : real(dp), intent(in) :: xa(:) ! mass fractions
177 : real(dp), intent(in) :: logRho ! density
178 : real(dp), intent(in) :: logT ! temperature
179 : real(dp), intent(in) :: lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT
180 : ! free_e := total combined number per nucleon of free electrons and positrons
181 : real(dp), intent(in) :: eta, d_eta_dlnRho, d_eta_dlnT
182 : ! eta := electron degeneracy parameter from eos
183 :
184 : ! OUTPUT
185 : real(dp), intent(out) :: kap_fracs(num_kap_fracs)
186 : real(dp), intent(out) :: kap ! opacity
187 : real(dp), intent(out) :: dlnkap_dlnRho ! partial derivative at constant T
188 : real(dp), intent(out) :: dlnkap_dlnT ! partial derivative at constant Rho
189 : real(dp), intent(out) :: dlnkap_dxa(:) ! partial derivative w.r.t. species
190 : integer, intent(out) :: ierr ! 0 means AOK.
191 :
192 : type (Kap_General_Info), pointer :: rq
193 :
194 86216 : real(dp) :: X, Y, Z, abar, zbar, z2bar, z53bar, ye, mass_correction, sumx
195 86216 : real(dp) :: XC, XN, XO, XNe
196 : integer :: i, iz
197 :
198 : ierr = 0
199 86216 : if (.not. kap_is_initialized) then
200 0 : ierr=-1
201 0 : return
202 : end if
203 :
204 86216 : call kap_ptr(handle,rq,ierr)
205 86216 : if (ierr /= 0) return
206 :
207 : call basic_composition_info( &
208 : species, chem_id, xa, X, Y, Z, &
209 86216 : abar, zbar, z2bar, z53bar, ye, mass_correction, sumx)
210 :
211 86216 : xc = 0; xn = 0; xo = 0; xne = 0
212 769728 : do i=1, species
213 683512 : iz = chem_isos% Z(chem_id(i))
214 86216 : select case(iz)
215 : case (6)
216 86216 : xc = xc + xa(i)
217 : case (7)
218 86216 : xn = xn + xa(i)
219 : case (8)
220 86216 : xo = xo + xa(i)
221 : case (10)
222 683512 : xne = xne + xa(i)
223 : end select
224 : end do
225 :
226 : call Get_kap_Results( &
227 : rq, zbar, X, Z, XC, XN, XO, XNe, logRho, logT, &
228 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
229 : eta, d_eta_dlnRho, d_eta_dlnT, &
230 86216 : kap_fracs, kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
231 :
232 : ! composition derivatives not implemented
233 769728 : dlnkap_dxa = 0
234 :
235 86216 : end subroutine kap_get
236 :
237 :
238 0 : subroutine kap_get_elect_cond_opacity( &
239 : handle, zbar, logRho, logT, &
240 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
241 86216 : use condint, only: do_electron_conduction
242 : use kap_def, only : kap_is_initialized, Kap_General_Info
243 : integer, intent(in) :: handle ! from alloc_kap_handle; in star, pass s% kap_handle
244 : real(dp), intent(in) :: zbar ! average ionic charge (for electron conduction)
245 : real(dp), intent(in) :: logRho ! the density
246 : real(dp), intent(in) :: logT ! the temperature
247 : real(dp), intent(out) :: kap ! electron conduction opacity
248 : real(dp), intent(out) :: dlnkap_dlnRho ! partial derivative at constant T
249 : real(dp), intent(out) :: dlnkap_dlnT ! partial derivative at constant Rho
250 : integer, intent(out) :: ierr ! 0 means AOK.
251 :
252 : type (Kap_General_Info), pointer :: rq
253 :
254 0 : if (.not. kap_is_initialized) then
255 0 : ierr=-1
256 0 : return
257 : end if
258 : ierr = 0
259 :
260 0 : call kap_ptr(handle,rq,ierr)
261 0 : if (ierr /= 0) return
262 :
263 : call do_electron_conduction( &
264 : rq, zbar, logRho, logT, &
265 0 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
266 :
267 0 : end subroutine kap_get_elect_cond_opacity
268 :
269 :
270 0 : subroutine kap_get_compton_opacity( &
271 : handle, &
272 : Rho, T, lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
273 : eta, d_eta_dlnRho, d_eta_dlnT, &
274 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
275 0 : use kap_eval, only: Compton_Opacity
276 : use kap_def, only : kap_is_initialized, Kap_General_Info
277 : integer, intent(in) :: handle ! kap handle; from star, pass s% kap_handle
278 : real(dp), intent(in) :: Rho, T
279 : real(dp), intent(in) :: lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT
280 : ! free_e := total combined number per nucleon of free electrons and positrons
281 : real(dp), intent(in) :: eta, d_eta_dlnT, d_eta_dlnRho
282 : ! eta := electron degeneracy parameter from eos
283 : real(dp), intent(out) :: kap ! electron conduction opacity
284 : real(dp), intent(out) :: dlnkap_dlnRho, dlnkap_dlnT
285 : integer, intent(out) :: ierr ! 0 means AOK.
286 :
287 : type (Kap_General_Info), pointer :: rq
288 :
289 0 : if (.not. kap_is_initialized) then
290 0 : ierr=-1
291 0 : return
292 : end if
293 : ierr = 0
294 :
295 0 : call kap_ptr(handle,rq,ierr)
296 0 : if (ierr /= 0) return
297 :
298 : call Compton_Opacity(rq, &
299 : Rho, T, lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
300 : eta, d_eta_dlnRho, d_eta_dlnT, &
301 0 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
302 :
303 0 : end subroutine kap_get_compton_opacity
304 :
305 :
306 0 : subroutine kap_get_radiative_opacity( &
307 : handle, &
308 : X, Z, XC, XN, XO, XNe, logRho, logT, &
309 : frac_lowT, frac_highT, frac_Type2, kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
310 :
311 0 : use kap_eval, only: Get_kap_Results_blend_T
312 : use kap_def, only : kap_is_initialized, Kap_General_Info
313 :
314 : ! INPUT
315 : integer, intent(in) :: handle ! kap handle; from star, pass s% kap_handle
316 : real(dp), intent(in) :: X, Z, XC, XN, XO, XNe ! composition
317 : real(dp), intent(in) :: logRho ! density
318 : real(dp), intent(in) :: logT ! temperature
319 :
320 : ! OUTPUT
321 : real(dp), intent(out) :: frac_lowT, frac_highT, frac_Type2
322 : real(dp), intent(out) :: kap ! opacity
323 : real(dp), intent(out) :: dlnkap_dlnRho ! partial derivative at constant T
324 : real(dp), intent(out) :: dlnkap_dlnT ! partial derivative at constant Rho
325 : integer, intent(out) :: ierr ! 0 means AOK.
326 :
327 : type (Kap_General_Info), pointer :: rq
328 :
329 0 : if (.not. kap_is_initialized) then
330 0 : ierr=-1
331 0 : return
332 : end if
333 : ierr = 0
334 :
335 0 : call kap_ptr(handle,rq,ierr)
336 0 : if (ierr /= 0) return
337 :
338 : call Get_kap_Results_blend_T( &
339 : rq, X, Z, XC, XN, XO, XNe, logRho, logT, &
340 0 : frac_lowT, frac_highT, frac_Type2, kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
341 :
342 0 : end subroutine kap_get_radiative_opacity
343 :
344 :
345 0 : subroutine kap_get_op_mono( &
346 : handle, zbar, logRho, logT, &
347 : ! args for op_mono
348 : use_op_mono_alt_get_kap, &
349 0 : nel, izzp, fap, fac, screening, umesh, semesh, ff, rs, &
350 : ! output
351 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
352 0 : use kap_def
353 : use kap_eval, only: combine_rad_with_conduction
354 : integer, intent(in) :: handle ! from alloc_kap_handle; in star, pass s% kap_handle
355 : real(dp), intent(in) :: zbar ! average ionic charge (for electron conduction)
356 : real(dp), intent(in) :: logRho ! the density
357 : real(dp), intent(in) :: logT ! the temperature
358 : ! args for op_mono_get_kap
359 : logical, intent(in) :: use_op_mono_alt_get_kap
360 : integer, intent(in) :: nel
361 : integer, intent(in) :: izzp(:) ! (nel)
362 : real(dp), intent(in) :: fap(:) ! (nel) number fractions of elements
363 : real(dp), intent(in) :: fac(:) ! (nel) scale factors for element opacity
364 : logical, intent(in) :: screening
365 : ! work arrays
366 : real, pointer :: umesh(:), semesh(:), ff(:,:,:,:), rs(:,:,:)
367 : ! umesh(nptot)
368 : ! semesh(nptot)
369 : ! ff(nptot, ipe, 4, 4)
370 : ! rs(nptot, 4, 4)
371 : ! output
372 : real(dp), intent(out) :: kap ! opacity
373 : real(dp), intent(out) :: dlnkap_dlnRho ! partial derivative at constant T
374 : real(dp), intent(out) :: dlnkap_dlnT ! partial derivative at constant Rho
375 : integer, intent(out) :: ierr ! 0 means AOK.
376 :
377 0 : real(dp) :: g1, kap_rad, dlnkap_rad_dlnRho, dlnkap_rad_dlnT
378 0 : real(dp) :: Rho, T
379 : type (Kap_General_Info), pointer :: rq
380 :
381 : ierr = 0
382 0 : if (.not. kap_is_initialized) then
383 0 : ierr=-1
384 0 : return
385 : end if
386 0 : call kap_ptr(handle,rq,ierr)
387 0 : if (ierr /= 0) return
388 :
389 0 : Rho = exp10(logRho)
390 0 : T = exp10(logT)
391 :
392 0 : if (use_op_mono_alt_get_kap) then
393 : call op_mono_alt_get_kap( &
394 : nel, izzp, fap, fac, logT, logRho, screening, &
395 : g1, dlnkap_rad_dlnT, dlnkap_rad_dlnRho, &
396 0 : umesh, semesh, ff, rs, ierr)
397 : else
398 : call op_mono_get_kap( &
399 : nel, izzp, fap, fac, logT, logRho, screening, &
400 : g1, dlnkap_rad_dlnT, dlnkap_rad_dlnRho, &
401 0 : umesh, semesh, ff, rs, ierr)
402 : end if
403 0 : if (ierr /= 0) return
404 :
405 0 : kap_rad = exp10(g1)
406 :
407 : call combine_rad_with_conduction( &
408 : rq, logRho, logT, zbar, &
409 : kap_rad, dlnkap_rad_dlnRho, dlnkap_rad_dlnT, &
410 0 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
411 :
412 0 : end subroutine kap_get_op_mono
413 :
414 :
415 :
416 : ! interface to OP routines as modified by Haili Hu for radiative levitation in diffusion
417 :
418 : ! ref: Hu et al MNRAS 418 (2011)
419 :
420 0 : subroutine load_op_mono_data(op_mono_data_path, op_mono_data_cache_filename, ierr)
421 0 : use kap_def
422 : use op_load, only: op_dload
423 : character (len=*), intent(in) :: op_mono_data_path, op_mono_data_cache_filename
424 : integer, intent(out) :: ierr
425 0 : if (.not. kap_is_initialized) then
426 0 : ierr=-1
427 : return
428 : end if
429 0 : call op_dload(op_mono_data_path, op_mono_data_cache_filename, ierr)
430 0 : end subroutine load_op_mono_data
431 :
432 0 : subroutine call_load_op_master(emesh_data_for_op_mono_path, ierr)
433 0 : use op_load_master, only: load_op_master
434 :
435 : character (len=*), intent(in) :: emesh_data_for_op_mono_path
436 : integer, intent(inout) :: ierr
437 :
438 0 : call load_op_master(emesh_data_for_op_mono_path, izz,ite,jne,epatom,amamu,sig,eumesh,ierr)
439 :
440 0 : end subroutine call_load_op_master
441 :
442 :
443 : ! sizes for work arrays
444 0 : subroutine get_op_mono_params(op_nptot, op_ipe, op_nrad)
445 : use kap_def
446 : use op_def, only: nptot, ipe, nrad
447 : integer, intent(out) :: op_nptot, op_ipe, op_nrad
448 0 : op_nptot = nptot
449 0 : op_ipe = ipe
450 0 : op_nrad = nrad
451 0 : end subroutine get_op_mono_params
452 :
453 :
454 : ! HH: Based on "op_ax.f"
455 : ! Input: kk = number of elements to calculate g_rad for
456 : ! iz1(kk) = charge of element to calculate g_rad for
457 : ! nel = number of elements in mixture
458 : ! izzp(nel) = charge of elements
459 : ! fap(nel) = number fractions of elements
460 : ! fac(nel) = scale factors for element opacity
461 : ! flux = local radiative flux (Lrad/4*pi*r^2)
462 : ! fltp = log10 T
463 : ! flrhop = log10 rho
464 : ! screening if true, use screening corrections
465 : ! Output: g1 = log10 kappa
466 : ! gx1 = d(log kappa)/d(log T)
467 : ! gy1 = d(log kappa)/d(log rho)
468 : ! gp1(kk) = d(log kappa)/d(log xi)
469 : ! grl1(kk) = log10 grad
470 : ! fx1(kk) = d(log grad)/d(log T)
471 : ! fy1(kk) = d(log grad)/d(log rho)
472 : ! grlp1(kk) = d(log grad)/d(log chi),
473 : ! chi is the fraction with which the number fraction is varied, i.e.:
474 : ! chi = nf_new/nf_previous
475 : ! where nf is the number fraction
476 : ! meanZ(nel) = average ionic charge of elements
477 : ! zetx1(nel) = d(meanZ)/d(log T)
478 : ! zety1(nel) = d(meanZ)/d(log rho)
479 : ! ierr = 0 for correct use
480 : ! ierr = 101 for rho out of range for this T
481 : ! ierr = 102 for T out of range
482 0 : subroutine op_mono_get_radacc( &
483 0 : kk, izk, nel, izzp, fap, fac, flux, fltp, flrhop, screening, &
484 0 : g1, grl1, &
485 : umesh, semesh, ff, ta, rs, ierr)
486 : use kap_def
487 : use op_eval, only: eval_op_radacc
488 : integer, intent(in) :: kk, nel
489 : integer, intent(in) :: izk(kk), izzp(nel)
490 : real(dp), intent(in) :: fap(:) ! (nel) number fractions of elements
491 : real(dp), intent(in) :: fac(:) ! (nel) scale factors for element opacity
492 : real(dp), intent(in) :: flux, fltp
493 : real(dp), intent(in) :: flrhop
494 : logical, intent(in) :: screening
495 : real(dp), intent(out) :: g1
496 : real(dp), intent(inout) :: &
497 : grl1(kk)
498 : ! work arrays
499 : real, pointer :: umesh(:), semesh(:), ff(:,:,:,:), ta(:,:,:,:), rs(:,:,:)
500 : ! umesh(nptot)
501 : ! semesh(nptot)
502 : ! ff(nptot, ipe, 4, 4)
503 : ! ta(nptot, nrad, 4, 4),
504 : ! rs(nptot, 4, 4)
505 : integer,intent(out) :: ierr
506 0 : if (.not. kap_is_initialized) then
507 0 : ierr=-1
508 : return
509 : end if
510 : call eval_op_radacc( &
511 : kk, izk, nel, izzp, fap, fac, flux, fltp, flrhop, screening, &
512 : g1, grl1, &
513 0 : umesh, semesh, ff, ta, rs, ierr)
514 0 : end subroutine op_mono_get_radacc
515 :
516 :
517 : ! note: for op mono, elements must come from the set given in op_mono_element_Z in kap_def.
518 :
519 : ! HH: Based on "op_mx.f", opacity calculations to be used for stellar evolution calculations
520 : ! Input: nel = number of elements in mixture
521 : ! izzp(nel) = charge of elements
522 : ! fap(nel) = number fractions of elements
523 : ! fac(nel) = scale factors for element opacity
524 : ! fltp = log10 (temperature)
525 : ! flrhop = log10 (mass density)
526 : ! screening if true, use screening corrections
527 : ! Output: g1 = log10 kappa
528 : ! gx1 = d(log kappa)/d(log T)
529 : ! gy1 = d(log kappa)/d(log rho)
530 : ! ierr = 0 for correct use
531 : ! ierr = 101 for rho out of range for this T
532 : ! ierr = 102 for T out of range
533 0 : subroutine op_mono_get_kap( &
534 0 : nel, izzp, fap, fac, fltp, flrhop, screening, &
535 : g1, gx1, gy1, &
536 : umesh, semesh, ff, rs, ierr)
537 0 : use kap_def
538 : use op_eval, only: eval_op_ev
539 : integer, intent(in) :: nel
540 : integer, intent(in) :: izzp(nel)
541 : real(dp), intent(in) :: fap(:) ! (nel) number fractions of elements
542 : real(dp), intent(in) :: fac(:) ! (nel) scale factors for element opacity
543 : real(dp), intent(in) :: fltp, flrhop
544 : logical, intent(in) :: screening
545 : real(dp), intent(inout) :: g1, gx1, gy1
546 : ! work arrays
547 : real, pointer :: umesh(:), semesh(:), ff(:,:,:,:), rs(:,:,:)
548 : ! umesh(nptot)
549 : ! semesh(nptot)
550 : ! ff(nptot, ipe, 4, 4)
551 : ! rs(nptot, 4, 4)
552 : ! s(nptot, nrad, 4, 4)
553 : integer,intent(out) :: ierr
554 0 : if (.not. kap_is_initialized) then
555 0 : ierr=-1
556 : return
557 : end if
558 : call eval_op_ev( &
559 : nel, izzp, fap, fac, fltp, flrhop, screening, &
560 : g1, gx1, gy1, &
561 0 : umesh, semesh, ff, rs, ierr)
562 0 : end subroutine op_mono_get_kap
563 :
564 :
565 : ! HH: Based on "op_mx.f", opacity calculations to be used for non-adiabatic pulsation calculations
566 : ! Special care is taken to ensure smoothness of opacity derivatives
567 : ! Input: nel = number of elements in mixture
568 : ! izzp(nel) = charge of elements
569 : ! fap(nel) = number fractions of elements
570 : ! fac(nel) = scale factors for element opacity
571 : ! fltp = log10 (temperature)
572 : ! flrhop = log10 (mass density)
573 : ! screening if true, use screening corrections
574 : ! Output: g1 = log10 kappa
575 : ! gx1 = d(log kappa)/d(log T)
576 : ! gy1 = d(log kappa)/d(log rho)
577 : ! ierr = 0 for correct use
578 : ! ierr = 101 for rho out of range for this T
579 : ! ierr = 102 for T out of range
580 0 : subroutine op_mono_alt_get_kap( &
581 0 : nel, izzp, fap, fac, fltp, flrhop, screening, &
582 : g1, gx1, gy1, &
583 : umesh, semesh, ff, rs, ierr)
584 0 : use kap_def
585 : use op_eval, only: eval_alt_op
586 : integer, intent(in) :: nel
587 : integer, intent(in) :: izzp(nel)
588 : real(dp), intent(in) :: fap(:) ! (nel) number fractions of elements
589 : real(dp), intent(in) :: fac(:) ! (nel) scale factors for element opacity
590 : real(dp), intent(in) :: fltp, flrhop
591 : logical, intent(in) :: screening
592 : real(dp), intent(out) :: g1, gx1, gy1
593 : ! work arrays
594 : real, pointer :: umesh(:), semesh(:), ff(:,:,:,:), rs(:,:,:)
595 : ! umesh(nptot)
596 : ! semesh(nptot)
597 : ! ff(nptot, ipe, 0:5, 0:5)
598 : ! rs(nptot, 0:5, 0:5)
599 : integer,intent(out) :: ierr
600 0 : if (.not. kap_is_initialized) then
601 0 : ierr=-1
602 : return
603 : end if
604 : call eval_alt_op( &
605 : nel, izzp, fap, fac, fltp, flrhop, screening, &
606 : g1, gx1, gy1, &
607 0 : umesh, semesh, ff, rs, ierr)
608 0 : end subroutine op_mono_alt_get_kap
609 :
610 :
611 0 : subroutine get_op_mono_args( &
612 0 : species, X, min_X_to_include, chem_id, chem_factors, &
613 0 : nel, izzp, fap, fac, ierr)
614 0 : use chem_def, only: chem_isos
615 : use kap_def
616 : integer, intent(in) :: species, chem_id(:)
617 : real(dp), intent(in) :: X(:) ! mass fractions (assumed baryonic)
618 : real(dp), intent(in) :: min_X_to_include ! skip iso if X < this
619 : real(dp), intent(in) :: chem_factors(:) ! (species)
620 : integer, intent(out) :: nel
621 : integer, intent(inout) :: izzp(:)
622 : real(dp), intent(inout) :: fap(:)
623 : real(dp), intent(inout) :: fac(:)
624 : integer,intent(out) :: ierr
625 :
626 : integer :: i, cid, j, Z, iel
627 0 : real(dp) :: tot
628 :
629 0 : ierr = 0
630 0 : if (.not. kap_is_initialized) then
631 0 : ierr=-1
632 0 : return
633 : end if
634 :
635 0 : nel = 0
636 0 : izzp(:) = 0
637 0 : fap(:) = 0d0
638 :
639 0 : do i=1,species
640 0 : if (X(i) < min_X_to_include) cycle
641 0 : cid = chem_id(i)
642 0 : Z = chem_isos% Z(cid)
643 0 : if (Z == 0) cycle
644 : ! change Z if necessary so that in op set
645 0 : do j=num_op_mono_elements,1,-1
646 0 : if (Z >= op_mono_element_Z(j)) then
647 : Z = op_mono_element_Z(j)
648 : exit
649 : end if
650 : end do
651 0 : iel = 0
652 0 : do j=1,nel
653 0 : if (izzp(j) == Z) then
654 : iel = j
655 : exit
656 : end if
657 : end do
658 0 : if (iel == 0) then
659 0 : nel = nel+1
660 0 : iel = nel
661 0 : izzp(nel) = Z
662 : end if
663 0 : fap(iel) = fap(iel) + X(i)/dble(chem_isos% Z_plus_N(cid))
664 0 : fac(iel) = chem_factors(i)
665 : end do
666 :
667 0 : tot = sum(fap(1:nel))
668 0 : if (tot <= 0d0) then
669 0 : ierr = -1
670 0 : return
671 : end if
672 :
673 0 : do j=1,nel
674 0 : fap(j) = fap(j)/tot ! number fractions
675 : end do
676 :
677 0 : end subroutine get_op_mono_args
678 :
679 :
680 0 : subroutine kap_get_control_namelist(handle, name, val, ierr)
681 0 : use kap_def
682 : use kap_ctrls_io, only: get_kap_controls
683 : integer, intent(in) :: handle ! kap handle; from star, pass s% kap_handle
684 : character(len=*),intent(in) :: name
685 : character(len=*),intent(out) :: val
686 : integer, intent(out) :: ierr
687 : type (kap_General_Info), pointer :: rq
688 : ierr = 0
689 0 : call kap_ptr(handle,rq,ierr)
690 0 : if(ierr/=0) return
691 0 : call get_kap_controls(rq, name, val, ierr)
692 :
693 0 : end subroutine kap_get_control_namelist
694 :
695 0 : subroutine kap_set_control_namelist(handle, name, val, ierr)
696 0 : use kap_def
697 : use kap_ctrls_io, only: set_kap_controls
698 : integer, intent(in) :: handle ! kap handle; from star, pass s% kap_handle
699 : character(len=*),intent(in) :: name
700 : character(len=*),intent(in) :: val
701 : integer, intent(out) :: ierr
702 : type (kap_General_Info), pointer :: rq
703 : ierr = 0
704 0 : call kap_ptr(handle,rq,ierr)
705 0 : if(ierr/=0) return
706 0 : call set_kap_controls(rq, name, val, ierr)
707 :
708 0 : end subroutine kap_set_control_namelist
709 :
710 :
711 :
712 0 : subroutine call_compute_grad_mombarg(k, j, blend, fk, T_face, Rho_face,&
713 : l, r, logKappa, lgrad, ierr)
714 0 : use op_eval_mombarg, only : compute_grad, compute_grad_fast
715 : !use crlibm_lib
716 :
717 : integer, intent(in) :: k, j
718 : real(dp), intent(in) :: fk(:)
719 : real(dp), intent(in) :: blend, T_face, Rho_face, l, r
720 : integer, intent(inout) :: ierr
721 : real(dp), intent(out) :: logKappa, lgrad(17)
722 0 : real(dp) :: logT_face, logRho_face, lgrad1(17), lgrad2(17)
723 :
724 : !write(*,*) 'calling compute_grad'
725 0 : logT_face = log10(T_face)
726 0 : logRho_face = log10(Rho_face)
727 :
728 : !call compute_grad(k, fk, logT_face, logRho_face,&
729 : ! l, r, &
730 : ! logKappa,lgrad, ierr,&
731 : ! izz,ite,jne,epatom,amamu,sig,eumesh)
732 0 : if (j == -1) then
733 : call compute_grad_fast(k, &
734 : fk_grad_pcg(1,:), logT_face, logRho_face, l, r, &
735 : lgrad1, ierr,&
736 0 : ite,jne,epatom,amamu,logT_pcg(1,:),logRho_pcg(1,:),lgamm_pcg(1,:,:),lkap_face_pcg(1,:))
737 : call compute_grad_fast(k, &
738 : fk_grad_pcg(2,:), logT_face, logRho_face, l, r, &
739 : lgrad2, ierr,&
740 0 : ite,jne,epatom,amamu,logT_pcg(2,:),logRho_pcg(2,:),lgamm_pcg(2,:,:),lkap_face_pcg(2,:))
741 0 : lgrad = (lgrad1*blend + lgrad2*(1-blend))
742 :
743 : else
744 : call compute_grad_fast(k, &
745 : fk_grad_pcg(j,:), logT_face, logRho_face, l, r, &
746 : lgrad, ierr,&
747 0 : ite,jne,epatom,amamu,logT_pcg(j,:),logRho_pcg(j,:),lgamm_pcg(j,:,:),lkap_face_pcg(j,:))
748 : end if
749 0 : end subroutine call_compute_grad_mombarg
750 :
751 0 : subroutine call_compute_gamma_grid_mombarg(j, fk, ierr)
752 : use op_eval_mombarg, only : compute_gamma_grid
753 : !use crlibm_lib
754 :
755 : integer, intent(in) :: j
756 : real(dp), intent(in) :: fk(:)
757 : integer, intent(inout) :: ierr
758 :
759 0 : fk_grad_pcg(j,:) = fk
760 :
761 : call compute_gamma_grid(ngp, fk, &
762 : lgamm_pcg(j,:,:), lkap_face_pcg(j,:), logT_pcg(j,:), logRho_pcg(j,:), ierr,&
763 0 : ite,jne,epatom,amamu,sig,eumesh)
764 :
765 :
766 0 : end subroutine call_compute_gamma_grid_mombarg
767 :
768 0 : subroutine call_compute_kappa_mombarg(handle, k,&
769 0 : fk, logT_cntr, logRho_cntr,&
770 : zbar, lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
771 : kap, dlnkap_dlnT, dlnkap_dlnRho, log_kap_rad_cell, ierr)
772 : use kap_eval, only: combine_rad_with_conduction
773 : use kap_def, only : kap_is_initialized, Kap_General_Info
774 : use op_eval_mombarg, only : compute_kappa, compute_kappa_fast, interpolate_kappa
775 : use chem_def, only: chem_isos, ih1, ihe3, ihe4, ic12, in14, io16, ine20, ina23, &
776 : img24, ial27, isi28, is32, iar40, ica40, icr52, imn55, ife56, ini58
777 :
778 : integer, intent(in) :: handle ! from alloc_kap_handle
779 : integer, intent(in) :: k
780 : real(dp), intent(in) :: fk(:)
781 : real(dp), intent(in) :: logT_cntr, logRho_cntr
782 : real(dp), intent(in) :: zbar ! average ionic charge (for electron conduction)
783 : real(dp), intent(in) :: lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT
784 : integer, intent(inout) :: ierr
785 : real(dp), intent(out) :: log_kap_rad_cell, kap, dlnkap_dlnT, dlnkap_dlnRho
786 :
787 : integer :: eid(17), ke
788 0 : real(dp) :: dlnkap_rad_dlnT, dlnkap_rad_dlnRho, kap_rad, delta, delta2 !,fk(17),fk_norm_fac
789 : type (Kap_General_Info), pointer :: rq
790 :
791 : ierr = 0
792 0 : if (.not. kap_is_initialized) then
793 0 : ierr=-1
794 0 : return
795 : end if
796 0 : call kap_ptr(handle,rq,ierr)
797 0 : if (ierr /= 0) return
798 :
799 :
800 : eid = [ ih1, ihe4, ic12, in14, io16, ine20, ina23, &
801 0 : img24, ial27, isi28, is32, iar40, ica40, icr52, imn55, ife56, ini58 ]
802 :
803 0 : if (initialize_fk_old) then
804 0 : fk_old = 0
805 0 : initialize_fk_old = .false.
806 : end if
807 :
808 :
809 0 : delta = MAXVAL(ABS(fk - fk_pcg)/fk_pcg, MASK=fk_pcg>0 )
810 0 : delta2 = MAXVAL(ABS(fk - fk_old(k,:))/fk_old(k,:), MASK=fk_old(k,:)>0 )
811 0 : if (SUM(fk_old(k,:)) == 0) then
812 0 : delta2 = 1d99
813 : end if
814 :
815 :
816 :
817 0 : if (delta > 1d-4) then
818 : if (delta2 > 1d-4 .or. ABS(logT_cntr - logT_cntr_old(k)) > 0.01d0 &
819 0 : .or. ABS(logRho_cntr - logRho_cntr_old(k)) > 0.1d0) then
820 : call compute_kappa(k,&
821 : fk, logT_cntr, logRho_cntr, &
822 : log_kap_rad_cell, dlnkap_rad_dlnT, dlnkap_rad_dlnRho, ierr,&
823 : ite,jne,epatom,amamu,sig,logT_grid_old(k,:,:),logRho_grid_old(k,:,:),&
824 0 : lkap_grid_old(k,:,:))
825 0 : fk_old(k,:) = fk
826 0 : logT_cntr_old(k) = logT_cntr
827 0 : logRho_cntr_old(k) = logRho_cntr
828 : else
829 : call interpolate_kappa(k,&
830 : logT_cntr, logRho_cntr, &
831 : log_kap_rad_cell, dlnkap_rad_dlnT, dlnkap_rad_dlnRho, ierr,&
832 : logT_grid_old(k,:,:),logRho_grid_old(k,:,:),&
833 0 : lkap_grid_old(k,:,:))
834 : end if
835 : else
836 : call compute_kappa_fast(k,&
837 : fk_pcg, logT_cntr, logRho_cntr, &
838 : log_kap_rad_cell, dlnkap_rad_dlnT, dlnkap_rad_dlnRho, ierr,&
839 0 : ite,jne,epatom,amamu,sig,lkap_ross_pcg)
840 :
841 : end if
842 :
843 0 : if (ierr == 1) return
844 :
845 0 : kap_rad = exp10(log_kap_rad_cell)
846 :
847 : call combine_rad_with_conduction( &
848 : rq, logRho_cntr, logT_cntr, zbar, &
849 : kap_rad, dlnkap_rad_dlnRho, dlnkap_rad_dlnT, &
850 0 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
851 :
852 0 : end subroutine call_compute_kappa_mombarg
853 :
854 0 : subroutine call_compute_kappa_grid_mombarg(fk, ierr)
855 0 : use op_eval_mombarg, only : compute_kappa_grid
856 :
857 : real(dp), intent(in) :: fk(:)
858 : integer, intent(inout) :: ierr
859 :
860 0 : fk_pcg = fk
861 :
862 : call compute_kappa_grid(fk, &
863 : lkap_ross_pcg, ierr,&
864 0 : ite,jne,epatom,amamu,sig)
865 :
866 0 : end subroutine call_compute_kappa_grid_mombarg
867 :
868 : end module kap_lib
|
