Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2018-2020 Aaron Dotter, Josiah Schwab & 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_aesopus
21 :
22 : use forum_m, only: hdf5io_t, OPEN_FILE_RO
23 : use math_lib
24 : use kap_def
25 :
26 : implicit none
27 :
28 : private
29 : public :: kap_aesopus_init
30 : public :: AESOPUS_get
31 :
32 : logical, parameter :: debug = .false.
33 :
34 : ! for 2-D interpolation
35 : integer :: ibcx=0, ibcy=0, ilinx=1, iliny=1
36 : integer :: ict(6) = [ 1, 1, 1, 0, 0, 0 ]
37 : real(dp), parameter :: bc(100) = 0d0
38 :
39 : contains
40 :
41 2 : subroutine kap_aesopus_init(rq, ierr)
42 : use kap_def, only: kap_aesopus_is_initialized
43 : type (Kap_General_Info), pointer :: rq
44 : integer, intent(out) :: ierr
45 :
46 : ! real(dp) :: kap, dlnkap_dlnRho, dlnkap_dlnT
47 :
48 2 : call read_kap_aesopus_tables(rq, ierr)
49 :
50 : ! call AESOPUS_get(0.02d0, 0.7d0, 0d0, 0d0, 0d0, -10d0, 4d0, &
51 : ! kap, dlnkap_dlnRho, dlnkap_dlnT,ierr)
52 :
53 : ! write(*,*) kap, dlnkap_dlnRho, dlnkap_dlnT
54 :
55 2 : if (ierr == 0) kap_aesopus_is_initialized = .true.
56 :
57 2 : end subroutine kap_aesopus_init
58 :
59 :
60 : subroutine kap_aesopus_shutdown
61 : use kap_def, only: kap_aesopus_is_initialized
62 : kap_aesopus_is_initialized = .false.
63 : end subroutine kap_aesopus_shutdown
64 :
65 :
66 2 : subroutine read_kap_aesopus_tables(rq, ierr)
67 :
68 : use interp_2d_lib_db, only: interp_mkbicub_db
69 :
70 : type (Kap_General_Info), pointer :: rq
71 : integer, intent(out) :: ierr
72 :
73 : character(len=256) :: filename
74 : integer :: n
75 : integer :: iX, iCO, iC, iN
76 :
77 2 : type(hdf5io_t) :: hi, hi_ts
78 :
79 2 : real(dp), allocatable :: table_data(:,:,:,:,:,:)
80 : integer :: table_size
81 :
82 : character(len=80) :: group_name ! output buffer
83 :
84 : character(len=30), parameter :: efmt = '(A14, 99ES12.3)'
85 : character(len=30), parameter :: ffmt = '(A14, 99F8.3)'
86 : character(len=30), parameter :: ifmt = '(A14, I4)'
87 :
88 : logical :: file_exists
89 :
90 : ! get the filename
91 2 : filename = trim(aesopus_filename)
92 2 : if (filename == '') then
93 0 : write(*,*) 'failed to specify AESOPUS_filename'
94 0 : ierr = -1
95 0 : return
96 : end if
97 :
98 : ! first try local directory
99 2 : inquire(file=trim(filename), exist=file_exists)
100 2 : if (.not. file_exists) then
101 : ! then try MESA directory
102 2 : filename = trim(kap_dir) // '/' // filename
103 2 : inquire(file=trim(filename), exist=file_exists)
104 2 : if (.not. file_exists) then
105 0 : write(*,*) 'failed to open AESOPUS file ' // trim(filename)
106 0 : ierr = -1
107 0 : return
108 : end if
109 : end if
110 :
111 2 : if (rq% show_info) write(*,*) 'read filename <' // trim(filename) // '>'
112 :
113 2 : ierr = 0
114 :
115 : ! open file (read-only)
116 2 : hi = hdf5io_t(filename, OPEN_FILE_RO)
117 :
118 2 : if (rq% show_info) write(*,*) 'AESOPUS composition parameters'
119 :
120 : ! read composition parameters
121 :
122 2 : call hi% read_dset('Zsun', kA% Zsun)
123 2 : if (rq% show_info) write(*,efmt) 'Zsun =', kA % Zsun
124 :
125 2 : call hi% read_dset('fCO_ref', kA% fCO_ref)
126 2 : if (rq% show_info) write(*,ffmt) 'fCO_ref =', kA % fCO_ref
127 :
128 2 : call hi% read_dset('fC_ref', kA% fC_ref)
129 2 : if (rq% show_info) write(*,ffmt) 'fC_ref =', kA % fC_ref
130 :
131 2 : call hi% read_dset('fN_ref', kA% fN_ref)
132 2 : if (rq% show_info) write(*,ffmt) 'fN_ref =', kA % fN_ref
133 :
134 2 : if (rq% show_info) then
135 0 : write(*,'(A)')
136 0 : write(*,*) 'AESOPUS logT and logR range (logR = logRho - 3 * logT + 18)'
137 : end if
138 :
139 : ! read logT
140 :
141 2 : call hi% alloc_read_dset('logTs', kA% logTs)
142 2 : kA% num_logTs = SIZE(kA% logTs)
143 2 : if (rq% show_info) write(*,ifmt) "num logTs =", kA% num_logTs
144 :
145 88 : kA% min_logT = minval(kA% logTs)
146 88 : kA% max_logT = maxval(kA% logTs)
147 :
148 2 : if (rq% show_info) then
149 0 : write(*,ffmt) 'logTs =', kA% logTs
150 : end if
151 :
152 : ! read logR
153 :
154 2 : call hi% alloc_read_dset('logRs', kA% logRs)
155 2 : kA% num_logRs = SIZE(kA% logRs)
156 2 : if (rq% show_info) write(*,ifmt) "num logRs =", kA% num_logRs
157 :
158 42 : kA% min_logR = minval(kA% logRs)
159 42 : kA% max_logR = maxval(kA% logRs)
160 :
161 2 : if (rq% show_info) then
162 0 : write(*,ffmt) 'logRs =', kA% logRs
163 : end if
164 :
165 2 : if (rq% show_info) then
166 0 : write(*,'(A)')
167 0 : write(*,*) 'AESOPUS metallicities'
168 0 : write(*,*) '(These Z values are the reference metallicities)'
169 : end if
170 :
171 : ! read Zs
172 :
173 2 : call hi% alloc_read_dset('Zs', kA% Zs)
174 2 : kA% num_Zs = SIZE(kA% Zs)
175 2 : if (rq% show_info) write(*,ifmt) "num Zs =", kA% num_Zs
176 :
177 : if (debug) write(*,*) 'Zs', kA% Zs
178 :
179 2 : if (rq% show_info) then
180 0 : write(*,efmt) "Zs =", kA% Zs
181 : end if
182 :
183 : ! pre-compute logZs
184 :
185 2 : allocate(kA% logZs(kA % num_Zs))
186 28 : do n = 1, kA % num_Zs
187 28 : kA% logZs(n) = safe_log10(kA% Zs(n))
188 : end do
189 :
190 : ! now, there is one group for each Z
191 : ! walk through these and construct a AESOPUS_TableSet for each one
192 :
193 28 : allocate(kA% ts(kA% num_Zs))
194 :
195 28 : do n = 1, kA% num_Zs
196 :
197 2 : associate(ts => kA% ts(n))
198 :
199 : ! get group name and open group
200 :
201 26 : write(group_name, 100) kA% Zs(n)
202 : 100 format(F8.6)
203 :
204 26 : if (rq% show_info) then
205 0 : write(*,'(A)')
206 0 : write(*,'(3A, ES9.3)') "Table ", trim(group_name), ": Z = ", kA% Zs(n)
207 : end if
208 :
209 26 : hi_ts = hdf5io_t(hi, group_name)
210 :
211 : ! read Xs
212 :
213 26 : call hi_ts% alloc_read_dset('Xs', ts% Xs)
214 : if (debug) write(*,*) "Xs", ts% Xs
215 26 : ts% num_Xs = SIZE(ts% Xs)
216 26 : if (rq% show_info) write(*,ifmt) "num Xs =", ts% num_Xs
217 :
218 26 : if (rq% show_info) then
219 0 : write(*,ffmt) "Xs =", ts% Xs
220 : end if
221 :
222 : ! read fCOs
223 :
224 26 : call hi_ts% alloc_read_dset('fCOs', ts% fCOs)
225 : if (debug) write(*,*) "fCOs", ts% fCOs
226 26 : ts% num_fCOs = SIZE(ts% fCOs)
227 26 : if (rq% show_info) write(*,ifmt) "num fCOs =", ts% num_fCOs
228 :
229 26 : if (rq% show_info) then
230 0 : write(*,ffmt) "fCOs =", ts% fCOs
231 : end if
232 :
233 : ! read fCs
234 :
235 26 : call hi_ts% alloc_read_dset('fCs', ts% fCs)
236 : if (debug) write(*,*) "fCs", ts% fCs
237 26 : ts% num_fCs = SIZE(ts% fCs)
238 26 : if (rq% show_info) write(*,ifmt) "num fCs =", ts% num_fCs
239 :
240 26 : if (rq% show_info) then
241 0 : write(*,ffmt) "fCs =", ts% fCs
242 : end if
243 :
244 : ! read fNs
245 :
246 26 : call hi_ts% alloc_read_dset('fNs', ts% fNs)
247 : if (debug) write(*,*) "fNs", ts% fNs
248 26 : ts% num_fNs = SIZE(ts% fNs)
249 26 : if (rq% show_info) write(*,ifmt) "num fNs =", ts% num_fNs
250 :
251 26 : if (rq% show_info) then
252 0 : write(*,ffmt) "fNs =", ts% fNs
253 : end if
254 :
255 : ! read opacities
256 :
257 26 : call hi_ts% alloc_read_dset('kap', table_data)
258 :
259 35984 : allocate(ts% t(ts % num_Xs, ts % num_fCOs, ts % num_fCs, ts % num_fNs))
260 :
261 156 : do iX = 1, ts % num_Xs
262 1326 : do iCO = 1, ts % num_fCOs
263 7150 : do iC = 1, ts % num_fCs
264 36270 : do iN = 1, ts % num_fNs
265 :
266 5850 : associate(t=> ts% t(iX, iCO, iC, iN))
267 :
268 : !allocate(t% kap(4, kA% num_logRs, kA% num_logTs))
269 :
270 29250 : table_size = kA% num_logRs*kA% num_logTs
271 29250 : allocate(t% kap(4*table_size))
272 :
273 : ! insert data
274 23400000 : t% kap(1:4*table_size:4) = reshape(table_data(iN, iC, iCO, iX, :, :), [table_size])
275 29250 : t% X = ts% Xs(iX)
276 29250 : t% Z = kA% Zs(n)
277 29250 : t% fCO = ts% fCOs(iCO)
278 29250 : t% fC = ts% fCs(iC)
279 29250 : t% fN = ts% fNs(iN)
280 :
281 : ! create interpolant
282 : call interp_mkbicub_db(kA% logRs, kA% num_logRs, &
283 : ka% logTs, kA% num_logTs, &
284 : t% kap, ka% num_logRs, &
285 : ibcx, bc, ibcx, bc, &
286 : ibcy, bc, ibcy, bc, &
287 58500 : ilinx, iliny, ierr)
288 :
289 : ! if (debug) call interp_evbicub_db(0d0, 3.85d0, &
290 : ! kA% logRs, kA% num_logRs, &
291 : ! ka% logTs, kA% num_logTs, &
292 : ! ilinx, iliny, &
293 : ! t% kap, ka% num_logRs, &
294 : ! ict, res, ierr)
295 : ! if (debug) write(*,'(10F10.4)') kA% logRs(17), kA% logTs(54), &
296 : ! t% X, t% fCO, t% fC, t% fN, &
297 : ! table_data(iN, iC, iCO, iX, 17, 54), res(1)
298 :
299 : end associate
300 :
301 : end do
302 : end do
303 : end do
304 : end do
305 :
306 26 : deallocate(table_data)
307 :
308 52 : call hi_ts%final()
309 :
310 : end associate
311 :
312 : end do
313 :
314 : ! close file
315 :
316 2 : call hi%final()
317 :
318 2 : if (rq% show_info) then
319 0 : write(*,'(A)')
320 0 : write(*,*) 'Finished reading AESOPUS tables'
321 0 : write(*,'(A)')
322 : end if
323 :
324 :
325 2 : end subroutine read_kap_aesopus_tables
326 :
327 :
328 2 : subroutine AESOPUS_interp(Zref, X, XC, XN, XO, logR, logT, res, ierr)
329 : !result=(logKappa, dlogKappa/dlogR, dlogKappa/dlogT)
330 : use interp_1d_def, only: pm_work_size
331 : use interp_1d_lib, only: interpolate_vector, interp_pm
332 : use num_lib, only: binary_search
333 : real(dp), intent(in) :: Zref, X, XC, XN, XO, logR, logT
334 : real(dp), intent(out) :: res(3)
335 : integer, intent(out) :: ierr
336 2 : real(dp), pointer :: work1(:)
337 : integer, parameter :: nZ = 3
338 : integer :: i,iZ
339 30 : real(dp) :: my_Z, raw_res(3, nZ), x_new(1), v_new(1)
340 : character(len=32) :: junk
341 : logical :: clipped
342 :
343 : ! result=0d0; iZ=0; ierr=0
344 :
345 2 : if (outside_R_and_T_bounds(logR,logT)) then
346 0 : write(*,*) 'AESOPUS_interp: logR, logT outside of table bounds'
347 0 : ierr = -1
348 0 : return
349 : end if
350 :
351 : ! restrict to range
352 2 : clipped = .false.
353 2 : if (Zref <= kA% Zs(1)) then
354 0 : my_Z = kA% Zs(1)
355 0 : iZ = 1
356 : clipped = .true.
357 2 : else if (Zref >= kA% Zs(kA% num_Zs)) then
358 0 : my_Z = kA% Zs(kA% num_Zs)
359 0 : iZ = kA% num_Zs
360 : clipped = .true.
361 : end if
362 :
363 : ! if clipped in Z, then just need one call
364 0 : if (clipped) then
365 0 : call AESOPUS_interp_fixedZref(iZ, X, XC, XN, XO, logR, logT, res, ierr)
366 0 : return
367 : end if
368 :
369 2 : my_Z = Zref
370 :
371 : ! it might be easier just to do linear interpolation
372 : ! but for now, we use an adapted version of what kapCN does
373 :
374 : ! require at least 3 Zs for interpolation
375 2 : if (nZ > kA% num_Zs) then
376 0 : write(*,*) 'AESOPUS_interp: insufficient number of Z values for interpolation'
377 0 : write(*,'(I2, A, I2, A)') nZ, ' values are required; ', kA% num_Zs, ' were provided'
378 0 : ierr = -1
379 0 : return
380 : end if
381 :
382 : ! binary_search returns iZ between 1 and num_Zs-1
383 : ! such that Zs(iZ) <= Z < Zs(iZ+1)
384 2 : iZ = binary_search(kA% num_Zs, kA% Zs, 0, my_Z)
385 :
386 : ! make sure this is in an acceptable range
387 : ! unless that would go off the ends
388 : ! this check is hard-coded assuming nZ = 3
389 2 : iZ = min(kA% num_Zs-1, max(iZ, 2))
390 :
391 : ! want to call at [iZ-1, iZ, iZ+1]
392 8 : do i = 1, nZ
393 8 : call AESOPUS_interp_fixedZref(iZ+i-2, X, XC, XN, XO, logR, logT, raw_res(:,i), ierr)
394 : end do
395 :
396 : nullify(work1)
397 2 : allocate(work1(nZ*pm_work_size))
398 2 : x_new(1)=safe_log10(my_Z)
399 :
400 : ! do the Z interpolation
401 : ! loop does over kap and its derivatives
402 8 : do i = 1, 3
403 : call interpolate_vector(nZ, kA% logZs(iZ-1:iZ+1), 1, &
404 : x_new, raw_res(i,:), v_new, &
405 : interp_pm, pm_work_size, &
406 6 : work1, junk, ierr)
407 8 : res(i) = v_new(1)
408 : end do
409 :
410 2 : deallocate(work1)
411 :
412 4 : end subroutine AESOPUS_interp
413 :
414 :
415 2 : logical function outside_R_and_T_bounds(logR,logT)
416 : real(dp), intent(in) :: logR, logT
417 : outside_R_and_T_bounds = &
418 : logR < kA% min_logR .or. logR > kA% max_logR .or. &
419 2 : logT < kA% min_logT .or. logT > kA% max_logT
420 2 : end function outside_R_and_T_bounds
421 :
422 :
423 6 : subroutine AESOPUS_interp_fixedZref(iZ, X, XC, XN, XO, logR, logT, res, ierr)
424 : ! simple interpolation in each of X, fCO, fC, fN
425 : ! at present, interpolation is linear (order = 2)
426 : integer, intent(in) :: iZ
427 : real(dp), intent(in) :: X, XC, XN, XO, logR, logT
428 : real(dp), intent(out) :: res(3)
429 : integer, intent(out) :: ierr
430 :
431 : integer, parameter :: npts = 2
432 54 : real(dp), dimension(npts) :: wX, wfCO, wfC, wfN
433 24 : real(dp) :: raw_res(3)
434 : integer :: iX, ifCO, ifC, ifN
435 : integer :: i, j, k, l
436 :
437 6 : real(dp) :: Zref, fCO, fC, fN
438 : logical :: clipped_X, clipped_fCO, clipped_fC, clipped_fN
439 :
440 6 : ierr = 0
441 6 : res = 0d0
442 : iX=0; ifCO=0; ifC=0; ifN=0
443 :
444 : ! AESOPUS defines these quantities as follows
445 :
446 : ! fCO=log10(XC/XO)-log10(XC/XO)ref
447 : ! fC=log10(XC)-log10(XC)ref
448 : ! fN=log10(XN)-log10(XN)ref
449 :
450 : ! Note that the reference values are different for different solar abundance patterns
451 :
452 6 : Zref = kA% Zs(iZ)
453 6 : fCO = safe_log10(XC/XO) - kA% fCO_ref
454 6 : fC = safe_log10(XC/Zref) - kA% fC_ref
455 6 : fN = safe_log10(XN/Zref) - kA% fN_ref
456 :
457 : if (debug) then
458 : write(*,*) 'call to AESOPUS_interp_RT'
459 : write(*,*) 'logR = ', logR
460 : write(*,*) 'logT = ', logT
461 : write(*,*) 'Zref = ', Zref
462 : write(*,*) 'X = ', X
463 : write(*,*) 'fCO = ', fCO
464 : write(*,*) 'fC = ', fC
465 : write(*,*) 'fN = ', fN
466 : end if
467 :
468 24 : associate(ts => kA% ts(iZ))
469 :
470 : ! get weights for a clipped linear interpolation in each parameter
471 6 : call clipped_linear_weights(X, ts% num_Xs, ts% Xs, iX, wX, clipped_X)
472 6 : call clipped_linear_weights(fCO, ts% num_fCOs, ts% fCOs, ifCO, wfCO, clipped_fCO)
473 6 : call clipped_linear_weights(fC, ts% num_fCs, ts% fCs, ifC, wfC, clipped_fC)
474 6 : call clipped_linear_weights(fN, ts% num_fNs, ts% fNs, ifN, wfN, clipped_fN)
475 :
476 : ! cycles prevent wastefully calling interp_RT with zero weights
477 24 : do i = 1, npts
478 12 : if (wX(i) == 0) cycle
479 42 : do j = 1, npts
480 24 : if (wfCO(j) == 0) cycle
481 84 : do k = 1, npts
482 48 : if (wfC(k) == 0) cycle
483 168 : do l = 1, npts
484 96 : if (wfN(l) == 0) cycle
485 :
486 : if (debug) then
487 : write(*,*) 'call to AESOPUS_interp_RT'
488 : write(*,*) iX+i-1, ifCO+j-1, ifC+k-1,ifN+l-1
489 : write(*,*) 'X = ', X, ts% Xs(iX+i-1), wX(i)
490 : write(*,*) 'fCO = ', fCO, ts% fCOs(ifCO+j-1), wfCO(j)
491 : write(*,*) 'fC = ', fC, ts% fCs(ifC+k-1), wfC(k)
492 : write(*,*) 'fN = ', fN, ts% fNs(ifN+l-1), wfN(l)
493 : end if
494 :
495 : ! now do the call and collect the results
496 :
497 96 : call AESOPUS_interp_RT(ts% t(iX+i-1,ifCO+j-1,ifC+k-1,ifN+l-1), logR, logT, raw_res, ierr)
498 96 : if (ierr /= 0) return
499 :
500 528 : res = res + wX(i)*wfCO(j)*wfC(k)*wfN(l) * raw_res
501 :
502 : end do
503 : end do
504 : end do
505 : end do
506 :
507 : end associate
508 :
509 : contains
510 :
511 24 : subroutine clipped_linear_weights(val, len, vec, loc, weights, clipped)
512 :
513 : ! calculate the weights for a linear interpolation
514 : ! clip to table edges
515 :
516 : use num_lib, only: binary_search
517 :
518 : real(dp), intent(in) :: val ! value
519 : integer, intent(in) :: len ! number of tabulated values
520 : real(dp), dimension(:), intent(in) :: vec ! tabulated values
521 : integer, intent(out) :: loc ! vec(loc) <= val <= vec(loc+1)
522 : real(dp), dimension(2), intent(out) :: weights ! for linear interpolation
523 : logical, intent(out) :: clipped ! did we clip? if so, only locs(1)/weights(1) matter
524 :
525 24 : weights = 0
526 :
527 : ! clip to range, if needed
528 24 : clipped = .false.
529 24 : if (val <= vec(1)) then
530 0 : loc = 1
531 0 : weights(1) = 1d0
532 : weights(2) = 0d0
533 0 : clipped = .true.
534 24 : else if (val >= vec(len)) then
535 0 : loc = len
536 0 : weights(1) = 1d0
537 : weights(2) = 0d0
538 0 : clipped = .true.
539 : end if
540 :
541 : ! find location and calculate linear weights
542 24 : if (.not. clipped) then
543 : ! binary_search returns k between 1 and n-1 such that vec(k) <= val < vec(k+1)
544 24 : loc = binary_search(len, vec, len/2, val)
545 24 : weights(2) = (val - vec(loc)) / (vec(loc+1) - vec(loc))
546 24 : weights(1) = 1d0 - weights(2)
547 : end if
548 :
549 24 : end subroutine clipped_linear_weights
550 :
551 : end subroutine AESOPUS_interp_fixedZref
552 :
553 :
554 96 : subroutine AESOPUS_interp_RT(t, logR, logT, res, ierr)
555 : use interp_2d_lib_db, only: interp_evbicub_db
556 : type(AESOPUS_Table) :: t
557 : real(dp), intent(in) :: logR, logT
558 : real(dp), intent(out) :: res(3)
559 : integer, intent(out) :: ierr
560 672 : real(dp) :: raw_res(6)
561 :
562 : if (debug) then
563 : write(*,*) 'inside call to AESOPUS_interp_RT'
564 : write(*,*) 'X = ', t% X
565 : write(*,*) 'Z = ', t% Z
566 : write(*,*) 'fCO = ', t% fCO
567 : write(*,*) 'fC = ', t% fC
568 : write(*,*) 'fN = ', t% fN
569 : end if
570 :
571 : call interp_evbicub_db(logR, logT, &
572 : kA% logRs, kA% num_logRs, &
573 : ka% logTs, kA% num_logTs, &
574 : ilinx, iliny, &
575 : t% kap, ka% num_logRs, &
576 96 : ict, raw_res, ierr)
577 :
578 384 : res = raw_res(1:3)
579 :
580 96 : end subroutine AESOPUS_interp_RT
581 :
582 :
583 2 : subroutine AESOPUS_get(Zref, X, XC, XN, XO, logRho, logT, kap, &
584 : dlnkap_dlnRho, dlnkap_dlnT,ierr)
585 : real(dp), intent(in) :: Zref, X, XC, XN, XO
586 : real(dp), intent(in) :: logRho, logT
587 : real(dp), intent(out) :: kap, dlnkap_dlnRho, dlnkap_dlnT
588 : integer, intent(out) :: ierr
589 8 : real(dp) :: logR, res(3)
590 :
591 : ierr = 0
592 2 : logR = logRho - 3d0*logT + 18d0
593 :
594 : if (debug) write(*,*) Zref, X, XC, XN, XO, logR, logT
595 2 : call AESOPUS_interp(Zref, X, XC, XN, XO, logR, logT, res, ierr)
596 2 : if (ierr == 0) then
597 2 : kap = exp10(res(1))
598 2 : dlnkap_dlnRho = res(2)
599 2 : dlnkap_dlnT = res(3) - 3d0*res(2)
600 : else
601 0 : kap = -1d0
602 0 : dlnkap_dlnRho = 0d0
603 0 : dlnkap_dlnT = 0d0
604 : end if
605 :
606 2 : end subroutine AESOPUS_get
607 :
608 : end module kap_aesopus
|
