Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010 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 kap_def
21 : use const_def, only: dp, strlen
22 :
23 : implicit none
24 :
25 : ! interfaces for procedure pointers
26 : abstract interface
27 :
28 : subroutine other_elect_cond_opacity_interface( &
29 : handle, &
30 : zbar, logRho, logT, &
31 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
32 : use const_def, only: dp
33 : implicit none
34 : integer, intent(in) :: handle ! kap handle; from star, pass s% kap_handle
35 : real(dp), intent(in) :: zbar ! average ionic charge (for electron conduction)
36 : real(dp), intent(in) :: logRho ! the density
37 : real(dp), intent(in) :: logT ! the temperature
38 : real(dp), intent(out) :: kap ! electron conduction opacity
39 : real(dp), intent(out) :: dlnkap_dlnRho ! partial derivative at constant T
40 : real(dp), intent(out) :: dlnkap_dlnT ! partial derivative at constant Rho
41 : integer, intent(out) :: ierr ! 0 means AOK.
42 : end subroutine other_elect_cond_opacity_interface
43 :
44 : subroutine other_compton_opacity_interface( &
45 : handle, &
46 : Rho, T, lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
47 : eta, d_eta_dlnRho, d_eta_dlnT, &
48 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
49 : use const_def, only: dp
50 : implicit none
51 : integer, intent(in) :: handle ! kap handle; from star, pass s% kap_handle
52 : real(dp), intent(in) :: Rho, T
53 : real(dp), intent(in) :: lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT
54 : ! free_e := total combined number per nucleon of free electrons and positrons
55 : real(dp), intent(in) :: eta, d_eta_dlnRho, d_eta_dlnT
56 : ! eta := electron degeneracy parameter from eos
57 : real(dp), intent(out) :: kap ! electron conduction opacity
58 : real(dp), intent(out) :: dlnkap_dlnRho, dlnkap_dlnT
59 : integer, intent(out) :: ierr ! 0 means AOK.
60 : end subroutine other_compton_opacity_interface
61 :
62 : subroutine other_radiative_opacity_interface( &
63 : handle, &
64 : X, Z, XC, XN, XO, XNe, logRho, logT, &
65 : frac_lowT, frac_highT, frac_Type2, kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
66 : use const_def, only: dp
67 : implicit none
68 : integer, intent(in) :: handle ! kap handle; from star, pass s% kap_handle
69 : real(dp), intent(in) :: X, Z, XC, XN, XO, XNe ! composition
70 : real(dp), intent(in) :: logRho ! density
71 : real(dp), intent(in) :: logT ! temperature
72 : real(dp), intent(out) :: frac_lowT, frac_highT, frac_Type2
73 : real(dp), intent(out) :: kap ! opacity
74 : real(dp), intent(out) :: dlnkap_dlnRho ! partial derivative at constant T
75 : real(dp), intent(out) :: dlnkap_dlnT ! partial derivative at constant Rho
76 : integer, intent(out) :: ierr ! 0 means AOK.
77 : end subroutine other_radiative_opacity_interface
78 :
79 : end interface
80 :
81 :
82 : logical, parameter :: show_allocations = .false. ! for debugging memory usage
83 :
84 : ! for kap output
85 : integer, parameter :: num_kap_fracs = 4
86 : integer, parameter :: i_frac_lowT = 1
87 : integer, parameter :: i_frac_highT = i_frac_lowT + 1
88 : integer, parameter :: i_frac_Type2 = i_frac_highT + 1
89 : integer, parameter :: i_frac_Compton = i_frac_Type2 + 1
90 :
91 : ! info about op_mono elements
92 : integer, parameter :: num_op_mono_elements = 17
93 : integer :: op_mono_element_Z(num_op_mono_elements)
94 : character(len=2) :: op_mono_element_name(num_op_mono_elements)
95 : real(dp) :: op_mono_element_mass(num_op_mono_elements)
96 :
97 :
98 : integer, parameter :: kap_table_fixed_metal_form = 1
99 : integer, parameter :: kap_table_co_enhanced_form = 2
100 :
101 :
102 : ! for fixed metal tables (no enhancements)
103 : type Kap_Z_Table ! holds pointers to all the X tables for a particular Z
104 : logical :: lowT_flag
105 : real(dp) :: Z
106 : integer :: num_Xs ! number of X's for this Z
107 : type (Kap_X_Table), dimension(:), pointer :: x_tables => null() ! in order of increasing X
108 : end type Kap_Z_Table
109 :
110 :
111 : ! note: logR = logRho - 3*logT + 18
112 :
113 : type Kap_X_Table
114 : logical :: not_loaded_yet
115 : real(dp) :: X
116 : real(dp) :: Z
117 : real(dp) :: logR_min
118 : real(dp) :: logR_max
119 : integer :: num_logRs
120 : integer :: ili_logRs ! =1 if logRs are evenly spaced
121 : real(dp), pointer :: logRs(:) => null() ! indexed from 1 to num_logRs
122 : real(dp) :: logT_min
123 : real(dp) :: logT_max
124 : integer :: num_logTs
125 : integer :: ili_logTs ! =1 if logTs are evenly spaced
126 : real(dp), pointer :: logTs(:) => null() ! indexed from 1 to num_logTs
127 : real(dp), pointer :: kap1(:) => null()
128 : end type Kap_X_Table
129 :
130 :
131 : ! for C/O enhanced tables
132 : type Kap_CO_Z_Table
133 : real(dp) :: Zbase, Zfrac_C, Zfrac_N, Zfrac_O, Zfrac_Ne
134 : real(dp) :: log10_Zbase ! log10(Zbase)
135 : type (Kap_CO_X_Table), dimension(:), pointer :: x_tables => null() ! stored in order of increasing X
136 : ! the X tables need not be equally spaced
137 : end type Kap_CO_Z_Table
138 :
139 : integer, parameter :: num_kap_CO_dXs = 8
140 : real(dp), parameter :: kap_CO_dXs(num_kap_CO_dXs) = &
141 : [ 0.00d0, 0.01d0, 0.03d0, 0.10d0, 0.20d0, 0.40d0, 0.60d0, 1.0d0 ]
142 :
143 : type Kap_CO_Table
144 : integer :: table_num ! the table number from the data file
145 : real(dp) :: X
146 : real(dp) :: Z
147 : real(dp) :: dXC
148 : real(dp) :: dXO
149 : real(dp) :: dXC_lookup
150 : real(dp) :: dXO_lookup
151 : real(dp), dimension(:), pointer :: kap1 => null()
152 : end type Kap_CO_Table
153 :
154 :
155 : integer, parameter :: max_num_CO_tables = 70
156 :
157 : ! standard number of CO tables for each X+Z combo
158 : ! X
159 : ! Z 0.0 0.1 0.35 0.7
160 : ! 0.0 58 58 51 32
161 : ! 0.001 58 58 51 30
162 : ! 0.004 58 58 51 30
163 : ! 0.010 58 58 51 30
164 : ! 0.020 58 58 51 30
165 : ! 0.030 58 58 49 30
166 : ! 0.100 58 53 43 26
167 :
168 : ! 1362 tables total
169 :
170 :
171 :
172 : type Kap_CO_X_Table
173 :
174 : logical :: not_loaded_yet
175 : real(dp) :: X
176 : real(dp) :: Z
177 : real(dp) :: logR_min
178 : real(dp) :: logR_max
179 : integer :: num_logRs
180 : integer :: ili_logRs ! =1 if logRs are evenly spaced
181 : real(dp), dimension(:), pointer :: logRs => null() ! indexed from 1 to num_logRs
182 : real(dp) :: logT_min
183 : real(dp) :: logT_max
184 : integer :: num_logTs
185 : integer :: ili_logTs ! =1 if logTs are evenly spaced
186 : real(dp), dimension(:), pointer :: logTs => null() ! indexed from 1 to num_logTs
187 :
188 : integer :: num_CO_tables
189 : ! the tables are in 3 groups
190 : ! 1) tables with dXC > dXO, ordered by increasing dXO, and by increasing dXC within same dXO.
191 : ! 2) tables with dXC = dXO, ordered by increasing value.
192 : ! 3) tables with dXC < dXO, ordered by increasing dXC, and by increasing dXO within same dXC.
193 : ! the spacing of dXC's is the same as dXO's, so there are as many tables in 3) as in 1).
194 : integer :: num_dXC_gt_dXO ! the number of tables with dXC > dXO
195 : integer :: CO_table_numbers(num_kap_CO_dXs,num_kap_CO_dXs)
196 : ! entry (i,j) is the co_index for table with dXC=Xs(i) and dXO=Xs(j), or -1 if no such table.
197 : integer :: next_dXO_table(max_num_CO_tables)
198 : ! entry (i) is the co_index for the table with same dXC and next larger dXO, or -1 if none such.
199 : integer :: next_dXC_table(max_num_CO_tables)
200 : ! entry (i) is the co_index for the table with same dXO and next larger dXC, or -1 if none such.
201 : type (Kap_CO_Table), dimension(:), pointer :: co_tables => null()
202 :
203 : end type Kap_CO_X_Table
204 :
205 : type Kap_General_Info
206 :
207 : real(dp) :: Zbase
208 :
209 : integer :: kap_option, kap_CO_option, kap_lowT_option
210 :
211 : ! blending in T is done between the following limits
212 : real(dp) :: kap_blend_logT_upper_bdy ! = 3.88d0 ! old value was 4.1d0
213 : real(dp) :: kap_blend_logT_lower_bdy ! = 3.80d0 ! old value was 4.0d0
214 : ! last time I looked, the table bottom for the higher T tables was logT = 3.75
215 : ! while max logT for the lower T Freeman tables was 4.5
216 : ! so for those, you need to keep kap_blend_logT_upper_bdy < 4.5
217 : ! and kap_blend_logT_lower_bdy > 3.75
218 : ! it is also probably a good idea to keep the blend away from H ionization
219 : ! logT upper of about 3.9 or a bit less will do that.
220 :
221 : logical :: cubic_interpolation_in_X
222 : logical :: cubic_interpolation_in_Z
223 :
224 : ! conductive opacities
225 : logical :: include_electron_conduction
226 : logical :: use_blouin_conductive_opacities
227 :
228 : logical :: use_Zbase_for_Type1
229 : logical :: use_Type2_opacities
230 :
231 : ! switch to Type1 if X too large
232 : real(dp) :: kap_Type2_full_off_X ! Type2 full off for X >= this
233 : real(dp) :: kap_Type2_full_on_X ! Type2 full on for X <= this
234 :
235 : ! switch to Type1 if dZ too small (dZ = Z - Zbase)
236 : real(dp) :: kap_Type2_full_off_dZ ! Type2 is full off for dZ <= this
237 : real(dp) :: kap_Type2_full_on_dZ ! Type2 can be full on for dZ >= this
238 :
239 : real(dp) :: logT_Compton_blend_hi, logR_Compton_blend_lo
240 :
241 : logical :: show_info
242 :
243 : ! debugging info
244 : logical :: dbg
245 : real(dp) :: logT_lo, logT_hi
246 : real(dp) :: logRho_lo, logRho_hi
247 : real(dp) :: X_lo, X_hi
248 : real(dp) :: Z_lo, Z_hi
249 :
250 : ! bookkeeping
251 : integer :: handle
252 : logical :: in_use
253 :
254 : ! User supplied inputs
255 : real(dp) :: kap_ctrl(10)
256 : integer :: kap_integer_ctrl(10)
257 : logical :: kap_logical_ctrl(10)
258 : character(len=strlen) :: kap_character_ctrl(10)
259 :
260 : ! other hooks
261 :
262 : logical :: use_other_elect_cond_opacity
263 : procedure (other_elect_cond_opacity_interface), pointer, nopass :: &
264 : other_elect_cond_opacity => null()
265 :
266 : logical :: use_other_compton_opacity
267 : procedure (other_compton_opacity_interface), pointer, nopass :: &
268 : other_compton_opacity => null()
269 :
270 : logical :: use_other_radiative_opacity
271 : procedure (other_radiative_opacity_interface), pointer, nopass :: &
272 : other_radiative_opacity => null()
273 :
274 : end type Kap_General_Info
275 :
276 : ! kap_options
277 : integer, parameter :: &
278 : kap_gn93 = 1, &
279 : kap_gs98 = 2, &
280 : kap_a09 = 3, &
281 : kap_OP_gs98 = 4, &
282 : kap_OP_a09 = 5, &
283 : kap_oplib_gs98 = 6, &
284 : kap_oplib_agss09 = 7, &
285 : kap_oplib_aag21 = 8, &
286 : kap_oplib_mb22 = 9, &
287 : kap_user = 10, &
288 : kap_test = 11, &
289 : kap_options_max = 11
290 :
291 :
292 : integer, parameter :: kap_max_dim = 50 !change this to make even larger grids in X and/or Z
293 :
294 : integer, dimension(kap_options_max) :: num_kap_Xs = 0
295 : real(dp), dimension(kap_max_dim, kap_options_max) :: kap_Xs = -1d0
296 :
297 : integer, dimension(kap_options_max) :: num_kap_Zs = 0
298 : real(dp), dimension(kap_max_dim, kap_options_max) :: kap_Zs= -1d0
299 :
300 : integer, dimension(kap_max_dim, kap_options_max) :: num_kap_Xs_for_this_Z = 0
301 :
302 :
303 : ! kap_CO_options
304 : integer, parameter :: &
305 : kap_CO_gn93 = 1, &
306 : kap_CO_gs98 = 2, &
307 : kap_CO_a09 = 3, &
308 : kap_CO_user = 4, &
309 : kap_CO_test = 5, &
310 : kap_CO_options_max = 5
311 :
312 :
313 : integer, dimension(kap_CO_options_max) :: num_kap_CO_Xs = 0
314 : real(dp), dimension(kap_max_dim, kap_CO_options_max) :: kap_CO_Xs = -1d0
315 :
316 : integer, dimension(kap_CO_options_max) :: num_kap_CO_Zs = 0
317 : real(dp), dimension(kap_max_dim, kap_CO_options_max) :: kap_CO_Zs= -1d0
318 :
319 : integer, dimension(kap_max_dim, kap_CO_options_max) :: num_kap_CO_Xs_for_this_Z = 0
320 :
321 :
322 : ! kap_lowT_options
323 : integer, parameter :: &
324 : kap_lowT_fa05_mb22= 1, &
325 : kap_lowT_fa05_aag21 = 2, &
326 : kap_lowT_Freedman11 = 3, &
327 : kap_lowT_fa05_gs98 = 4, &
328 : kap_lowT_fa05_gn93 = 5, &
329 : kap_lowT_fa05_a09p = 6, &
330 : kap_lowT_af94_gn93 = 7, &
331 : kap_lowT_rt14_ag89 = 8, &
332 : kap_lowT_kapCN = 9, &
333 : kap_lowT_AESOPUS = 10, &
334 : kap_lowT_user = 11, &
335 : kap_lowT_test = 12, &
336 : kap_lowT_options_max = 12
337 :
338 :
339 : integer, dimension(kap_lowT_options_max) :: num_kap_lowT_Xs = 0
340 : real(dp), dimension(kap_max_dim, kap_lowT_options_max) :: kap_lowT_Xs = -1d0
341 :
342 : integer, dimension(kap_lowT_options_max) :: num_kap_lowT_Zs = 0
343 : real(dp), dimension(kap_max_dim, kap_lowT_options_max) :: kap_lowT_Zs= -1d0
344 :
345 : integer, dimension(kap_max_dim, kap_lowT_options_max) :: num_kap_lowT_Xs_for_this_Z = 0
346 :
347 :
348 : character (len=256) :: &
349 : kap_option_str(kap_options_max), &
350 : kap_CO_option_str(kap_CO_options_max), &
351 : kap_lowT_option_str(kap_lowT_options_max)
352 :
353 : type Kap_Z_Table_Array ! in order of increasing Z
354 : type (Kap_Z_Table), dimension(:), pointer :: ar
355 : end type Kap_Z_Table_Array
356 :
357 : type Kap_CO_Z_Table_Array ! in order of increasing Z
358 : type (Kap_CO_Z_Table), dimension(:), pointer :: ar
359 : end type Kap_CO_Z_Table_Array
360 :
361 : type (Kap_Z_Table_Array), dimension(kap_options_max) :: kap_z_tables
362 : type (Kap_Z_Table_Array), dimension(kap_lowT_options_max) :: kap_lowT_z_tables
363 : type (Kap_CO_Z_Table_Array), dimension(kap_CO_options_max) :: kap_co_z_tables
364 :
365 :
366 : ! NOTE: in the following, "log" means base 10, "ln" means natural log, and units are cgs.
367 :
368 : integer, parameter :: sz_per_kap_point = 4
369 : !
370 : ! function f(x,y) with samples f(i,j) has bicubic spline fit s(x,y).
371 : ! compact representation of spline fit uses 4 entries as follows:
372 : !
373 : ! d(1,i,j) = s(i,j)
374 : ! d(2,i,j) = d2s_dx2(i,j)
375 : ! d(3,i,j) = d2s_dy2(i,j)
376 : ! d(4,i,j) = d4s_dx2_dy2(i,j)
377 : !
378 : ! given f(i,j), the spline fitting code can compute the other entries
379 : !
380 : ! given d(1:4,i,j), spline interpolation code can compute s(x,y)
381 : ! and also the partials ds_dx(x,y) and ds_dy(x,y)
382 : !
383 :
384 :
385 : logical :: kap_is_initialized = .false.
386 :
387 : character (len=1000) :: kap_dir, kap_cache_dir, kap_temp_cache_dir
388 : logical :: kap_use_cache = .true.
389 : logical :: kap_read_after_write_cache = .true.
390 : logical :: clip_to_kap_table_boundaries = .true. ! typically, this should be set true.
391 : ! if this is set true, then temperature and density args are
392 : ! clipped to the boundaries of the table.
393 : real(dp), parameter :: kap_min_logRho = -40d0
394 : ! below this, clip logRho and set partials wrt logRho to zero
395 :
396 :
397 : integer, parameter :: max_kap_handles = 10
398 : type (Kap_General_Info), target :: kap_handles(max_kap_handles)
399 :
400 : !for kapCN
401 : logical, parameter :: kapCN_use_cache = .true.
402 : integer, parameter :: num_kapCN_Xs = 3
403 : integer, parameter :: num_kapCN_Zs = 14
404 : integer, parameter :: num_kapCN_fCs = 7
405 : integer, parameter :: num_kapCN_fNs = 3
406 : integer, parameter :: kapCN_num_logT = 18
407 : integer, parameter :: kapCN_num_logR = 17
408 : integer, parameter :: kapCN_tbl_size = kapCN_num_logR*kapCN_num_logT ! 306
409 : integer, parameter :: kapCN_num_tbl = num_kapCN_Xs*num_kapCN_fCs*num_kapCN_fNs ! 63
410 :
411 : real(dp), target :: kapCN_Z(num_kapCN_Zs)
412 : real(dp), target :: kapCN_fN(num_kapCN_fNs,num_kapCN_Zs)
413 : real(dp), target :: kapCN_fC(num_kapCN_fCs,num_kapCN_Zs)
414 : real(dp), target :: kapCN_X(num_kapCN_Xs) = [ 0.5d0, 0.7d0, 0.8d0 ]
415 : real(dp), target :: kapCN_logT(kapCN_num_logT), kapCN_logR(kapCN_num_logR)
416 : real(dp) :: kapCN_min_logR, kapCN_max_logR, kapCN_min_logT, kapCN_max_logT
417 : logical :: kapCN_is_initialized = .false.
418 :
419 : real(dp), parameter :: kapCN_ZC=0.1644d0, kapCN_ZN=0.0532d0
420 :
421 : !stores a table for one {Z,X,fC,fN}
422 : type KapCN_Table
423 : real(dp) :: X
424 : real(dp) :: Z, Zbase
425 : real(dp) :: fC
426 : real(dp) :: fN
427 : real(dp) :: falpha
428 : real(dp), pointer :: kap(:) => null()
429 : end type KapCN_Table
430 :
431 : !stores a set of tables for one Z
432 : type KapCN_Set
433 : character(len=13) :: filename
434 : logical :: not_loaded_yet
435 : real(dp) :: Zbase
436 : real(dp) :: fC(num_kapCN_fCs)
437 : real(dp) :: fN(num_kapCN_fNs)
438 : type(KapCN_Table) :: t(kapCN_num_tbl)
439 : end type KapCN_Set
440 : type(KapCN_Set), target :: kCN(num_kapCN_Zs)
441 :
442 :
443 : logical :: kap_aesopus_is_initialized = .false.
444 : character(len=256) :: aesopus_filename
445 :
446 : ! stores a table for one {Z,X,fCO,fC,fN}
447 : type AESOPUS_Table
448 : real(dp) :: X
449 : real(dp) :: Z, Zbase
450 : real(dp) :: fCO
451 : real(dp) :: fC
452 : real(dp) :: fN
453 :
454 : ! this has to be a pointer because of the interp2d routines
455 : real(dp), dimension(:), pointer :: kap => null()
456 :
457 : end type AESOPUS_Table
458 :
459 : ! stores a set of tables for one Z
460 : type AESOPUS_TableSet
461 :
462 : integer :: num_Xs
463 : integer :: num_fCOs
464 : integer :: num_fCs
465 : integer :: num_fNs
466 :
467 : real(dp), dimension(:), allocatable :: Xs, fCOs, fCs, fNs
468 :
469 : type(AESOPUS_Table), dimension(:,:,:,:), allocatable :: t
470 :
471 : end type AESOPUS_TableSet
472 :
473 :
474 : type kapAESOPUS
475 :
476 : integer :: num_logTs
477 : integer :: num_logRs
478 :
479 : real(dp), dimension(:), allocatable :: logTs(:), logRs(:)
480 :
481 : real(dp) :: min_logR, max_logR
482 : real(dp) :: min_logT, max_logT
483 :
484 : integer :: num_Zs
485 : real(dp), dimension(:), allocatable :: Zs, logZs
486 :
487 : real(dp) :: Zsun
488 : real(dp) :: fCO_ref, fC_ref, fN_ref
489 :
490 : type(AESOPUS_TableSet), dimension(:), allocatable :: ts
491 :
492 : end type kapAESOPUS
493 :
494 : type(kapAESOPUS) :: kA
495 :
496 : logical :: kap_test_partials
497 : real(dp) :: kap_test_partials_val, kap_test_partials_dval_dx
498 :
499 : contains
500 :
501 :
502 5 : subroutine kap_def_init(kap_cache_dir_in)
503 : use chem_def
504 : use utils_lib, only : mkdir
505 : use const_def, only: mesa_data_dir, mesa_caches_dir, mesa_temp_caches_dir, use_mesa_temp_cache
506 : character (*), intent(in) :: kap_cache_dir_in
507 : integer :: i
508 :
509 5 : kap_test_partials = .false.
510 :
511 5 : if (len_trim(kap_cache_dir_in) > 0) then
512 0 : kap_cache_dir = kap_cache_dir_in
513 5 : else if (len_trim(mesa_caches_dir) > 0) then
514 0 : kap_cache_dir = trim(mesa_caches_dir) // '/kap_cache'
515 : else
516 5 : kap_cache_dir = trim(mesa_data_dir) // '/kap_data/cache'
517 : end if
518 5 : call mkdir(kap_cache_dir)
519 :
520 5 : clip_to_kap_table_boundaries = .true.
521 :
522 55 : do i=1,max_kap_handles
523 50 : kap_handles(i)% handle = i
524 55 : kap_handles(i)% in_use = .false.
525 : end do
526 :
527 : op_mono_element_Z(1:num_op_mono_elements) = [ &
528 5 : 1, 2, 6, 7, 8, 10, 11, 12, 13, 14, 16, 18, 20, 24, 25, 26, 28 ]
529 : op_mono_element_name(1:num_op_mono_elements) = [ &
530 : 'h ', 'he', 'c ', 'n ', 'o ', 'ne', 'na', 'mg', 'al', &
531 90 : 'si', 's ', 'ar', 'ca', 'cr', 'mn', 'fe', 'ni' ]
532 : op_mono_element_mass(1:num_op_mono_elements) = [ &
533 : 1.0080d0, 4.0026d0, 12.0111d0, 14.0067d0, 15.9994d0, 20.179d0, &
534 : 22.9898d0, 24.305d0, 26.9815d0, 28.086d0, 32.06d0, 39.948d0, &
535 5 : 40.08d0, 51.996d0, 54.9380d0, 55.847d0, 58.71d0 ]
536 :
537 5 : kap_temp_cache_dir=trim(mesa_temp_caches_dir)//'/kap_cache'
538 5 : if(use_mesa_temp_cache) call mkdir(kap_temp_cache_dir)
539 :
540 5 : kap_option_str(kap_gn93) = 'gn93'
541 5 : kap_option_str(kap_gs98) = 'gs98'
542 5 : kap_option_str(kap_a09) = 'a09'
543 5 : kap_option_str(kap_OP_gs98) = 'OP_gs98'
544 5 : kap_option_str(kap_OP_a09) = 'OP_a09_nans_removed_by_hand'
545 5 : kap_option_str(kap_oplib_gs98) = 'oplib_gs98'
546 5 : kap_option_str(kap_oplib_agss09) = 'oplib_agss09'
547 5 : kap_option_str(kap_oplib_aag21) = 'oplib_aag21'
548 5 : kap_option_str(kap_oplib_mb22) = 'oplib_mb22'
549 5 : kap_option_str(kap_test) = 'test'
550 :
551 5 : kap_CO_option_str(kap_CO_gn93) = 'gn93_co'
552 5 : kap_CO_option_str(kap_CO_gs98) = 'gs98_co'
553 5 : kap_CO_option_str(kap_CO_a09) = 'a09_co'
554 5 : kap_CO_option_str(kap_CO_test) = 'test_co'
555 :
556 5 : kap_lowT_option_str(kap_lowT_fa05_mb22) = 'lowT_fa05_mb22'
557 5 : kap_lowT_option_str(kap_lowT_fa05_aag21) = 'lowT_fa05_aag21'
558 5 : kap_lowT_option_str(kap_lowT_Freedman11) = 'lowT_Freedman11'
559 5 : kap_lowT_option_str(kap_lowT_fa05_gs98) = 'lowT_fa05_gs98'
560 5 : kap_lowT_option_str(kap_lowT_fa05_gn93) = 'lowT_fa05_gn93'
561 5 : kap_lowT_option_str(kap_lowT_fa05_a09p) = 'lowT_fa05_a09p'
562 5 : kap_lowT_option_str(kap_lowT_af94_gn93) = 'lowT_af94_gn93'
563 5 : kap_lowT_option_str(kap_lowT_rt14_ag89) = 'lowT_rt14_ag89'
564 5 : kap_lowT_option_str(kap_lowT_kapCN) = 'kapCN'
565 5 : kap_lowT_option_str(kap_lowT_AESOPUS) = 'AESOPUS'
566 5 : kap_lowT_option_str(kap_lowT_test) = 'lowT_test'
567 :
568 60 : do i=1,kap_options_max
569 60 : nullify(kap_z_tables(i)% ar)
570 : end do
571 65 : do i=1,kap_lowT_options_max
572 65 : nullify(kap_lowT_z_tables(i)% ar)
573 : end do
574 30 : do i=1,kap_CO_options_max
575 30 : nullify(kap_co_z_tables(i)% ar)
576 : end do
577 :
578 60 : do i=1, kap_options_max
579 5 : select case (i)
580 : case (6:9)
581 20 : num_kap_Xs(i) = 30
582 : kap_Xs(1:num_kap_Xs(i), i) = [0.0d0, 0.000001d0, 0.00001d0, 0.000d0, 0.001d0, 0.01d0, &
583 : 0.05d0, 0.1d0, 0.15d0, 0.2d0, 0.25d0, 0.3d0, 0.35d0, 0.4d0, 0.45d0, 0.5d0, 0.55d0, &
584 : 0.6d0, 0.65d0, 0.7d0, 0.75d0, 0.8d0, 0.85d0, 0.9d0, 0.91d0, 0.92d0, 0.93d0, 0.94d0, &
585 620 : 0.95d0,1.0d0]
586 20 : num_kap_Zs(i) = 41
587 : kap_Zs(1:num_kap_Zs(i), i) = [0.0d0, 0.000001d0, 0.00001d0, 0.00003d0, 0.00007d0, 0.0001d0, &
588 : 0.0003d0 ,0.0007d0, 0.001d0, 0.002d0, 0.003d0, 0.004d0, 0.006d0, 0.008d0, 0.01d0, 0.012d0, &
589 : 0.014d0, 0.015d0, 0.016d0, 0.017d0, 0.018d0, 0.019d0, 0.02d0, 0.021d0, 0.022d0, 0.023d0, &
590 : 0.024d0, 0.025d0, 0.026d0, 0.028d0, 0.03d0, 0.035d0, 0.04d0, 0.05d0, 0.06d0, 0.07d0, 0.08d0, &
591 840 : 0.09d0, 0.1d0, 0.15d0, 0.2d0]
592 : num_kap_Xs_for_this_Z(1:num_kap_Zs(i), i) = [ num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), &
593 : num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), &
594 : num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), &
595 : num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), &
596 : num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), &
597 : num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), &
598 : num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), &
599 : num_kap_Xs(i), num_kap_Xs(i)-1, num_kap_Xs(i)-2, num_kap_Xs(i)-3, num_kap_Xs(i)-4, &
600 980 : num_kap_Xs(i)-5, num_kap_Xs(i)-6, num_kap_Xs(i)-7, num_kap_Xs(i)-8]
601 : case DEFAULT
602 35 : num_kap_Xs(i) = 10
603 : kap_Xs(1:num_kap_Xs(i), i) = [0.00d0, 0.10d0, 0.20d0, &
604 385 : 0.35d0, 0.50d0, 0.70d0, 0.80d0, 0.90d0, 0.95d0, 1d0]
605 35 : num_kap_Zs(i) = 13
606 : kap_Zs(1:num_kap_Zs(i), i) = [0.000d0, 0.0001d0, 0.0003d0, &
607 : 0.001d0, 0.002d0, 0.004d0, 0.01d0, 0.02d0, 0.03d0, &
608 490 : 0.04d0, 0.06d0, 0.08d0, 0.100d0 ]
609 : num_kap_Xs_for_this_Z(1:num_kap_Zs(i), i) = [ num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), &
610 : num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), num_kap_Xs(i), &
611 545 : num_kap_Xs(i), num_kap_Xs(i)-1, num_kap_Xs(i)-1, num_kap_Xs(i)-2 ]
612 : end select
613 : end do
614 :
615 65 : do i=1, kap_lowT_options_max
616 5 : select case (i)
617 : case(kap_lowT_Freedman11)
618 5 : num_kap_lowT_Xs(i) = 1
619 10 : kap_lowT_Xs(1:num_kap_lowT_Xs(i), i) = [ 0.00d0 ]
620 5 : num_kap_lowT_Zs(i) = 7
621 : kap_lowT_Zs(1:num_kap_lowT_Zs(i), i) = &
622 40 : [ 0.01d0, 0.02d0, 0.04d0, 0.100d0, 0.200d0, 0.63d0, 1.00d0 ]
623 40 : num_kap_lowT_Xs_for_this_Z(1:num_kap_lowT_Zs(i), i) = num_kap_lowT_Xs(i)
624 : case DEFAULT
625 55 : num_kap_lowT_Xs(i) = 10
626 : kap_lowT_Xs(1:num_kap_lowT_Xs(i), i) = [0.00d0, 0.10d0, 0.20d0, &
627 605 : 0.35d0, 0.50d0, 0.70d0, 0.80d0, 0.90d0, 0.95d0, 1d0]
628 55 : num_kap_lowT_Zs(i) = 13
629 : kap_lowT_Zs(1:num_kap_lowT_Zs(i), i) = [0.000d0, 0.0001d0, 0.0003d0, &
630 : 0.001d0, 0.002d0, 0.004d0, 0.01d0, 0.02d0, 0.03d0, &
631 770 : 0.04d0, 0.06d0, 0.08d0, 0.100d0 ]
632 : num_kap_lowT_Xs_for_this_Z(1:num_kap_lowT_Zs(i), i) = [ num_kap_lowT_Xs(i), num_kap_lowT_Xs(i), num_kap_lowT_Xs(i), &
633 : num_kap_lowT_Xs(i), num_kap_lowT_Xs(i), num_kap_lowT_Xs(i), num_kap_lowT_Xs(i), num_kap_lowT_Xs(i), &
634 830 : num_kap_lowT_Xs(i), num_kap_lowT_Xs(i), num_kap_lowT_Xs(i)-1, num_kap_lowT_Xs(i)-1, num_kap_lowT_Xs(i)-2 ]
635 : end select
636 : end do
637 :
638 30 : do i=1, kap_CO_options_max
639 5 : select case (i)
640 : case DEFAULT
641 25 : num_kap_CO_Xs(i) = 5
642 : kap_CO_Xs(1:num_kap_CO_Xs(i), i) = &
643 150 : [ 0.00d0, 0.03d0, 0.10d0, 0.35d0, 0.70d0 ]
644 25 : num_kap_CO_Zs(i) = 8
645 : kap_CO_Zs(1:num_kap_CO_Zs(i),i) = &
646 225 : [ 0.000d0, 0.001d0, 0.004d0, 0.010d0, 0.020d0, 0.030d0, 0.050d0, 0.100d0 ]
647 225 : num_kap_CO_Xs_for_this_Z(1:num_kap_CO_Zs(i), i) = num_kap_CO_Xs(i)
648 : end select
649 : end do
650 :
651 :
652 10 : end subroutine kap_def_init
653 :
654 :
655 :
656 9 : integer function do_alloc_kap(ierr)
657 : integer, intent(out) :: ierr
658 : integer :: i
659 9 : ierr = 0
660 9 : do_alloc_kap = -1
661 18 : !$omp critical (kap_handle)
662 15 : do i = 1, max_kap_handles
663 15 : if (.not. kap_handles(i)% in_use) then
664 9 : kap_handles(i)% in_use = .true.
665 9 : do_alloc_kap = i
666 9 : exit
667 : end if
668 : end do
669 : !$omp end critical (kap_handle)
670 9 : if (do_alloc_kap == -1) then
671 0 : ierr = -1
672 0 : return
673 : end if
674 9 : if (kap_handles(do_alloc_kap)% handle /= do_alloc_kap) then
675 0 : ierr = -1
676 0 : return
677 : end if
678 5 : end function do_alloc_kap
679 :
680 :
681 1 : subroutine do_free_kap(handle)
682 : integer, intent(in) :: handle
683 1 : if (handle >= 1 .and. handle <= max_kap_handles) &
684 1 : kap_handles(handle)% in_use = .false.
685 1 : end subroutine do_free_kap
686 :
687 :
688 86270 : subroutine get_kap_ptr(handle,rq,ierr)
689 : integer, intent(in) :: handle
690 : type (Kap_General_Info), pointer :: rq
691 : integer, intent(out):: ierr
692 86270 : if (handle < 1 .or. handle > max_kap_handles) then
693 0 : ierr = -1
694 0 : return
695 : end if
696 86270 : rq => kap_handles(handle)
697 86270 : ierr = 0
698 : end subroutine get_kap_ptr
699 :
700 :
701 5470 : subroutine get_output_string(x,xstr,ierr) !works with X and Z
702 : real(dp), intent(in) :: x
703 : character(len=*), intent(out) :: xstr
704 : integer, intent(out) :: ierr
705 : character(len=1), parameter :: c(9)=['1','2','3','4','5','6','7','8','9']
706 : character(len=8) :: str
707 : integer :: i, j, k
708 :
709 5470 : if(X < 0d0.or.X>1d0) then
710 0 : xstr='bad'
711 0 : ierr=-1
712 0 : return
713 : end if
714 5470 : ierr=0
715 5470 : write(str,'(f8.6)') X
716 5470 : k=0
717 54700 : do i=1,9
718 49230 : j=index(str,c(i),back=.true.)
719 54700 : k=max(k,j)
720 : end do
721 5470 : xstr=str(1:max(k,3))
722 5470 : end subroutine get_output_string
723 :
724 :
725 3 : subroutine do_Free_Kap_Tables
726 : integer :: i
727 :
728 36 : do i=1,kap_options_max
729 33 : if (associated(kap_z_tables(i)% ar)) &
730 8 : call free_z_tables(kap_z_tables(i)% ar)
731 : end do
732 39 : do i=1,kap_lowT_options_max
733 36 : if (associated(kap_lowT_z_tables(i)% ar)) &
734 6 : call free_z_tables(kap_lowT_z_tables(i)% ar)
735 : end do
736 18 : do i=1,kap_CO_options_max
737 15 : if (associated(kap_co_z_tables(i)% ar)) &
738 10 : call free_co_z_tables(kap_co_z_tables(i)% ar)
739 : end do
740 :
741 : contains
742 :
743 8 : subroutine free_z_tables(z_tables)
744 : type (Kap_Z_Table), dimension(:), pointer :: z_tables
745 : integer :: num_Zs
746 : integer :: iz
747 8 : num_Zs = size(z_tables,dim=1)
748 112 : do iz = 1, num_Zs
749 104 : if (associated(z_tables(iz)% x_tables)) &
750 112 : call free_x_tables(z_tables(iz)% x_tables, z_tables(iz)% num_Xs)
751 : end do
752 8 : deallocate(z_tables)
753 : nullify(z_tables)
754 8 : end subroutine free_z_tables
755 :
756 104 : subroutine free_x_tables(x_tables, num_Xs)
757 : type (Kap_X_Table), dimension(:), pointer :: x_tables
758 : integer, intent(in) :: num_Xs
759 : integer :: ix
760 1112 : do ix = 1, num_Xs
761 1008 : if (associated(x_tables(ix)% logRs)) then
762 4 : deallocate(x_tables(ix)% logRs)
763 4 : nullify(x_tables(ix)% logRs)
764 : end if
765 1008 : if (associated(x_tables(ix)% logTs)) then
766 4 : deallocate(x_tables(ix)% logTs)
767 4 : nullify(x_tables(ix)% logTs)
768 : end if
769 1112 : if (associated(x_tables(ix)% kap1)) then
770 4 : deallocate(x_tables(ix)% kap1)
771 4 : nullify(x_tables(ix)% kap1)
772 : end if
773 : end do
774 104 : if (associated(x_tables)) then
775 104 : deallocate(x_tables)
776 : nullify(x_tables)
777 : end if
778 104 : end subroutine free_x_tables
779 :
780 7 : subroutine free_co_z_tables(co_z_tables)
781 : type (Kap_CO_Z_Table), dimension(:), pointer :: co_z_tables
782 : integer :: num_Zs
783 : integer :: iz
784 7 : num_Zs = size(co_z_tables,dim=1)
785 63 : do iz = 1, num_Zs
786 63 : if (associated(co_z_tables(iz)% x_tables)) then
787 56 : call free_co_x_tables(co_z_tables(iz)% x_tables)
788 : end if
789 : end do
790 7 : deallocate(co_z_tables)
791 : nullify(co_z_tables)
792 7 : end subroutine free_co_z_tables
793 :
794 56 : subroutine free_co_x_tables(x_tables)
795 : type (Kap_CO_X_Table), dimension(:), pointer :: x_tables
796 : ! stored in order of increasing X
797 : integer :: num_Xs
798 : integer :: ix
799 56 : num_Xs = size(x_tables,dim=1)
800 336 : do ix = 1, num_Xs
801 280 : call free_co_table(x_tables(ix)% co_tables, x_tables(ix)% num_CO_tables)
802 280 : if (associated(x_tables(ix)% logRs)) then
803 4 : deallocate(x_tables(ix)% logRs)
804 4 : nullify(x_tables(ix)% logRs)
805 : end if
806 336 : if (associated(x_tables(ix)% logTs)) then
807 4 : deallocate(x_tables(ix)% logTs)
808 4 : nullify(x_tables(ix)% logTs)
809 : end if
810 : end do
811 56 : if (associated(x_tables))then
812 56 : deallocate(x_tables)
813 : nullify(x_tables)
814 : end if
815 56 : end subroutine free_co_x_tables
816 :
817 280 : subroutine free_co_table(co_tables, num_COs)
818 : type (Kap_CO_Table), dimension(:), pointer :: co_tables
819 : integer, intent(in) :: num_COs
820 : integer :: ico
821 14336 : do ico = 1, num_COs
822 14336 : if (associated(co_tables(ico)% kap1)) then
823 232 : deallocate(co_tables(ico)% kap1)
824 232 : nullify(co_tables(ico)% kap1)
825 : end if
826 : end do
827 280 : if (associated(co_tables)) then
828 280 : deallocate(co_tables)
829 : nullify(co_tables)
830 : end if
831 280 : end subroutine free_co_table
832 :
833 : end subroutine do_Free_Kap_Tables
834 :
835 0 : end module kap_def
|
