Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2018-2020 Aaron Dotter & 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 : ! kapCN by Aaron Dotter
21 : ! This is a MESA module that reads in and interpolates the low-T
22 : ! opacities by M.T. Lederer & B. Aringer, 2009, A&A, 494, 403
23 : ! see doc/ReadMe for details of the data files; it uses MESA
24 : ! modules for interpolation and a few other things
25 : ! kapCN_ZC=0.1644, kapCN_ZN=0.0532
26 :
27 : module kapcn
28 :
29 : use kap_def
30 : use num_lib, only: binary_search
31 : use const_def, only: dp
32 : use math_lib
33 :
34 : implicit none
35 :
36 : private
37 : public :: kapcn_init
38 : public :: kapCN_get
39 :
40 : !local stuff
41 : logical, parameter :: debug = .false.
42 : character(len=32), parameter :: kapCN_param_file = 'kR_Z_fCN.data'
43 :
44 : !for kap_def
45 :
46 : !for 2-D interpolation
47 : integer :: ibcx=0, ibcy=0, ilinx=1, iliny=1
48 : integer :: ict(6) = [ 1, 1, 1, 0, 0, 0 ]
49 : real(dp), parameter :: bc(kapCN_num_logT) = 0d0
50 :
51 : contains
52 :
53 : !for kapCN
54 0 : subroutine kapCN_init(show_info, ierr)
55 : use const_def, only: mesa_data_dir
56 : use kap_def, only: kapcn_is_initialized
57 : logical, intent(in) :: show_info
58 : integer, intent(out) :: ierr
59 : integer :: iZ
60 : type(kapCN_set), pointer :: k
61 :
62 0 : do iZ=1,num_kapCN_Zs
63 0 : k => kCN(iZ)
64 0 : k% not_loaded_yet = .true.
65 0 : k% Zbase = kapCN_Z(iZ)
66 0 : k% fC = kapCN_fC(:,iZ)
67 0 : k% fN = kapCN_fN(:,iZ)
68 : end do
69 :
70 0 : call read_kapCN_tables(ierr)
71 :
72 0 : if(ierr==0) kapCN_is_initialized = .true.
73 :
74 0 : end subroutine kapCN_init
75 :
76 : subroutine kapCN_shutdown
77 : use kap_def, only: kapcn_is_initialized
78 : kapCN_is_initialized = .false.
79 : end subroutine kapCN_shutdown
80 :
81 0 : subroutine read_kapCN_tables(ierr)
82 : integer, intent(out) :: ierr
83 : character(len=256) :: filename
84 : character(len=4) :: prefix, string_Z(num_kapCN_Zs)
85 : character(len=6) :: tmp_Z, suffix
86 : integer :: i, io
87 :
88 0 : prefix = 'kR_Z'
89 0 : suffix = '.data'
90 :
91 0 : string_Z = ''
92 :
93 0 : filename = trim(kap_dir) // '/' // trim(kapCN_param_file)
94 :
95 : if(debug) write(*,*) ' parameter file = ', trim(filename)
96 :
97 0 : open(newunit=io,file=trim(filename))
98 0 : read(io,*) !skip first line
99 0 : do i=num_kapCN_Zs,1,-1 !read backwards into array so Z is increasing
100 0 : read(io,'(f7.5,11f9.1)') kapCN_Z(i), kapCN_fC(:,i), kapCN_fN(:,i)
101 : end do
102 0 : close(io)
103 :
104 0 : do i=1,num_kapCN_Zs
105 0 : write(tmp_Z,'(1p,1e6.1e1)') kapCN_Z(i)
106 0 : string_Z(i) = tmp_Z(1:1) // tmp_Z(4:6)
107 : end do
108 :
109 0 : kCN(:)% filename = prefix // string_Z // trim(suffix)
110 :
111 0 : do i=1,num_kapCN_Zs
112 0 : call read_one_kapCN_table(i,ierr)
113 : end do
114 :
115 0 : kapCN_min_logR = minval(kapCN_logR)
116 0 : kapCN_max_logR = maxval(kapCN_logR)
117 0 : kapCN_min_logT = minval(kapCN_logT)
118 0 : kapCN_max_logT = maxval(kapCN_logT)
119 :
120 0 : end subroutine read_kapCN_tables
121 :
122 :
123 0 : subroutine read_one_kapCN_table(iZ,ierr)
124 : use interp_2d_lib_db, only: interp_mkbicub_db
125 : integer, intent(in) :: iZ
126 : integer, intent(out) :: ierr
127 : character(len=256) :: ascii_file, cache_file
128 : character(len=10) :: my_Z
129 : integer :: i, io, j
130 0 : real(dp) :: c_div_o, y
131 0 : real(dp) :: table(kapCN_num_logR,kapCN_num_logT)
132 0 : real(dp), pointer :: logR(:), logT(:)
133 : logical :: have_cache
134 : type(kapCN_set), pointer :: k
135 :
136 0 : ierr=0
137 :
138 0 : k => kCN(iZ)
139 0 : logR => kapCN_logR
140 0 : logT => kapCN_logT
141 :
142 : ! data files
143 0 : ascii_file = trim(kap_dir) // '/' // trim(k% filename)
144 0 : cache_file = trim(kap_dir) // '/cache/' // trim(k% filename(1:9)) // 'bin'
145 :
146 : if(debug) write(*,*) ' ascii file = ', trim(ascii_file)
147 : if(debug) write(*,*) ' cache file = ', trim(cache_file)
148 :
149 : !check for cached bin file; if it exists, read it and exit
150 : if(kapCN_use_cache)then
151 0 : inquire(file=trim(cache_file),exist=have_cache)
152 0 : if(have_cache)then
153 : !read cache
154 0 : open(newunit=io,file=trim(cache_file),form='unformatted',iostat=ierr)
155 0 : if(ierr/=0) write(*,*) &
156 0 : ' read_one_kapCN_table: problem opening cache file for read'
157 0 : read(io) logR, logT
158 0 : read(io) k% t(:)% X
159 0 : read(io) k% t(:)% Z
160 0 : read(io) k% t(:)% fC
161 0 : read(io) k% t(:)% fN
162 0 : read(io) k% t(:)% falpha
163 0 : do i=1,kapCN_num_tbl
164 0 : allocate(k% t(i)% kap(4*kapCN_tbl_size))
165 0 : read(io) k% t(i)% kap
166 : end do
167 0 : close(io)
168 : !table is now fully loaded
169 0 : k% not_loaded_yet = .false.
170 0 : return
171 : end if
172 : end if
173 :
174 : !no cache file, so read ascii file
175 0 : open(newunit=io,file=trim(ascii_file),iostat=ierr)
176 0 : if(ierr/=0) write(*,*) &
177 0 : ' read_one_kapCN_table: problem opening ascii file for read'
178 0 : do i=1,4
179 0 : read(io,*) !skip header
180 : end do
181 0 : read(io,'(6x,a10)') my_Z
182 : if(debug) write(*,*) ' has Z = ', my_Z, ' Zbase = ', k% Zbase
183 :
184 0 : do i=1,14
185 0 : read(io,*)
186 : end do
187 :
188 0 : do i=1,kapCN_num_tbl
189 0 : read(io,*) j, k% t(i)% x, y, k% t(i)% z, c_div_o, &
190 0 : k% t(i)% fc, k% t(i)% fn, k% t(i)% falpha
191 : end do
192 :
193 0 : read(io,*) !last line of #'s
194 :
195 0 : do i=1,kapCN_num_tbl
196 0 : do j=1,14
197 0 : read(io,*) !skip individual headers
198 : end do
199 0 : read(io,'(5x,17f7.3)') logR
200 0 : do j=1,kapCN_num_logT
201 0 : read(io,'(f5.3,17f7.3)') logT(j), table(:,j)
202 : end do
203 0 : allocate(k% t(i)% kap(4*kapCN_tbl_size))
204 0 : k% t(i)% kap(1:4*kapCN_tbl_size:4) = reshape(table,[kapCN_tbl_size])
205 : end do
206 :
207 0 : close(io)
208 :
209 : !create interpolants for tables
210 0 : do i=1,kapCN_num_tbl
211 : call interp_mkbicub_db( logR, kapCN_num_logR, &
212 : logT, kapCN_num_logT, &
213 : k% t(i)% kap, kapCN_num_logR, &
214 : ibcx, bc, ibcx, bc, &
215 : ibcy, bc, ibcy, bc, &
216 0 : ilinx, iliny, ierr)
217 : end do
218 :
219 : !table is now fully loaded
220 0 : k% not_loaded_yet = .false.
221 :
222 :
223 : !write cache file
224 0 : open(newunit=io,file=trim(cache_file),form='unformatted',iostat=ierr)
225 0 : if(ierr/=0) write(*,*) &
226 0 : ' read_one_kapCN_table: problem opening cache file for write'
227 0 : write(io) logR, logT
228 0 : write(io) k% t(:)% x
229 0 : write(io) k% t(:)% z
230 0 : write(io) k% t(:)% fC
231 0 : write(io) k% t(:)% fN
232 0 : write(io) k% t(:)% falpha
233 0 : do i=1,kapCN_num_tbl
234 0 : write(io) k% t(i)% kap
235 : end do
236 0 : close(io)
237 :
238 0 : end subroutine read_one_kapCN_table
239 :
240 :
241 : !this one is specifically for use with MESA
242 :
243 :
244 0 : subroutine kapCN_interp(Z,X,fC,fN,logR,logT,result,ierr)
245 : !result=(logKappa, dlogKappa/dlogR, dlogKappa/dlogT)
246 : use interp_1d_def, only: pm_work_size
247 : use interp_1d_lib, only: interpolate_vector, interp_pm
248 : real(dp), intent(in) :: Z, X, fC, fN, logR, logT
249 : real(dp), intent(out) :: result(3)
250 : integer, intent(out) :: ierr
251 0 : real(dp), pointer :: Z_ary(:), work1(:), logZ_ary(:)
252 : integer, parameter :: nZ = 4
253 : integer :: i,iZ
254 0 : real(dp) :: my_Z, res(3,nZ), x_new(1), v_new(1)
255 : character(len=32) :: junk
256 :
257 0 : result=0d0; iZ=0; ierr=0
258 :
259 0 : if(.not.kapCN_is_initialized)then
260 0 : write(*,*) ' kapCN is not initialized; call kapCN_init()'
261 0 : ierr=-99
262 0 : return
263 : end if
264 :
265 0 : if(outside_R_and_T_bounds(logR,logT))then
266 : !write(*,*) 'kapCN_interp: logR, logT outside of table bounds'
267 0 : ierr=-1
268 0 : return
269 : end if
270 :
271 0 : Z_ary => kapCN_Z
272 0 : my_Z = max(min(Z,Z_ary(num_kapCN_Zs)),Z_ary(1))
273 0 : iZ=binary_search(num_kapCN_Zs,Z_ary, iZ, Z)
274 0 : iZ = max(nz/2,min(iZ,num_kapCN_Zs-nz/2))
275 :
276 : !check to see if exact match in Z, then just need one call
277 0 : if(Z==Z_ary(iZ))then
278 0 : call kapCN_interp_fixedZ(iZ,X,fC,fN,logR,logT,result)
279 0 : return
280 : end if
281 :
282 : !else do the full Z interpolation
283 0 : do i=1,nZ
284 0 : call kapCN_interp_fixedZ(iZ+i-2,X,fC,fN,logR,logT,res(:,i))
285 : end do
286 :
287 : nullify(work1)
288 0 : allocate(work1(nZ*pm_work_size))
289 0 : x_new(1)=log10(my_Z)
290 :
291 0 : allocate(logZ_ary(iZ-1:iZ+2))
292 :
293 0 : do i=iZ-1,iZ+2
294 0 : logZ_ary(i) = log10(z_ary(i))
295 : end do
296 :
297 0 : do i=1,3
298 : call interpolate_vector(nZ,logZ_ary, 1, &
299 : x_new, res(i,:), v_new, &
300 : interp_pm, pm_work_size, &
301 0 : work1, junk, ierr )
302 0 : result(i)=v_new(1)
303 : end do
304 :
305 0 : deallocate(work1, logz_ary)
306 :
307 0 : end subroutine kapCN_interp
308 :
309 :
310 0 : logical function outside_R_and_T_bounds(logR,logT)
311 : real(dp), intent(in) :: logR, logT
312 : outside_R_and_T_bounds = &
313 : logR < kapCN_min_logR .or. logR > kapCN_max_logR .or. &
314 0 : logT < kapCN_min_logT .or. logT > kapCN_max_logT
315 0 : end function outside_R_and_T_bounds
316 :
317 :
318 0 : subroutine kapCN_interp_fixedZ(iZ,X,fC,fN,logR,logT,result)
319 : !using simple quadratic interpolation in each of X, fC, fN
320 : integer, intent(in) :: iZ
321 : real(dp), intent(in) :: X, fC, fN, logR, logT
322 : real(dp), intent(out) :: result(3)
323 : integer, parameter :: my_num_kapCN_fCs = 3, f1=num_kapCN_fCs*num_kapCN_Xs, f2=num_kapCN_fCs
324 0 : real(dp) :: my_X, my_fC, my_fN, wX(num_kapCN_Xs), wfN(num_kapCN_fNs), wfC(my_num_kapCN_fCs)
325 0 : real(dp) :: res(3)
326 : integer :: iX, ifC, ifN, i,j,tbl,n,nlo,nhi,ierr
327 : real(dp), pointer :: X_ary(:), fC_ary(:), fN_ary(:)
328 :
329 0 : result = 0d0; iX=0; ifC=0; ifN=0
330 0 : X_ary => kapCN_X; fC_ary => kapCN_fC(:,iZ); fN_ary => kapCN_fN(:,iZ)
331 :
332 : !limit input quantities to lie within tabulated range
333 0 : my_X = max(min(X, X_ary(num_kapCN_Xs)),X_ary(1))
334 0 : my_fC = max(min(fC,fC_ary(num_kapCN_fCs)),fC_ary(1))
335 0 : my_fN = max(min(fN,fN_ary(num_kapCN_fNs)),fN_ary(1))
336 :
337 : !locate inputs in arrays
338 0 : iX = binary_search(num_kapCN_Xs, X_ary, iX, my_X)
339 0 : ifC= binary_search(num_kapCN_fCs, fC_ary, ifC, my_fC)
340 0 : ifN= binary_search(num_kapCN_fNs, fN_ary, ifN, my_fN)
341 :
342 : !make sure 1 < ifC < num_kapCN_fCs
343 0 : ifC = max(2,min(ifC,num_kapCN_fCs-1))
344 :
345 : !interpolation coefficients in X
346 0 : call interp(X_ary,wX,my_X,num_kapCN_Xs)
347 :
348 : !interpolation coefficients in fN
349 0 : call interp(fN_ary(:),wfN,my_fN,num_kapCN_fNs)
350 :
351 : !interpolation coefficients in fC
352 0 : call interp(fC_ary(ifC-1:ifC+1),wfC,my_fC,my_num_kapCN_fCs)
353 :
354 0 : do i=1,num_kapCN_fNs
355 0 : do j=1,num_kapCN_Xs
356 0 : do n=1,my_num_kapCN_fCs
357 0 : res=0d0
358 0 : tbl = f1*(i-1) + f2*(j-1) + (ifC+n-2)
359 0 : nlo = 3*(n-1)+1
360 0 : nhi = nlo+2
361 0 : call kapCN_interp_RT(iZ,tbl,logR,logT,res,ierr)
362 0 : result = result + wfN(i)*wX(j)*wfC(n)*res
363 : end do
364 : end do
365 : end do
366 :
367 0 : if(debug) write(*,'(1p9e12.4)') my_X, my_fC, my_fN, result
368 :
369 : contains
370 :
371 0 : subroutine interp(a,b,x,n)
372 : ! {a} are the tabulated values for use in interpolation
373 : ! {b} are coefficients of the interpolating polynomial
374 : ! x is the abscissa to be interpolated
375 : ! n is the number of points to be used, interpolating polynomial
376 : ! has order n-1
377 : integer, intent(in) :: n
378 : real(dp), intent(in) :: a(n), x
379 : real(dp), intent(out) :: b(n)
380 : integer :: i,j
381 0 : do i=1,n
382 0 : b(i)=1d0
383 0 : do j=1,n
384 0 : if(j/=i) b(i)=b(i)*(x-a(j))/(a(i)-a(j))
385 : end do
386 : end do
387 0 : end subroutine interp
388 :
389 : end subroutine kapCN_interp_fixedZ
390 :
391 :
392 0 : subroutine kapCN_interp_RT(iZ,tbl,logR,logT,result,ierr)
393 : use interp_2d_lib_db, only: interp_evbicub_db
394 : integer, intent(in) :: iZ, tbl
395 : real(dp), intent(in) :: logR, logT
396 : real(dp), intent(out) :: result(3)
397 : integer, intent(out) :: ierr
398 0 : real(dp) :: res(6)
399 : real(dp), pointer :: logR_ary(:), logT_ary(:)
400 : type(kapCN_set), pointer :: k
401 :
402 0 : k => kCN(iZ)
403 0 : logR_ary => kapCN_logR
404 0 : logT_ary => kapCN_logT
405 :
406 : call interp_evbicub_db( logR, logT, &
407 : logR_ary, kapCN_num_logR, &
408 : logT_ary, kapCN_num_logT, &
409 : ilinx, iliny, &
410 : k% t(tbl)% kap, kapCN_num_logR, &
411 0 : ict, res, ierr)
412 :
413 0 : result = res(1:3)
414 :
415 0 : end subroutine kapCN_interp_RT
416 :
417 0 : subroutine kapCN_get(Z,X,fC,fN,logRho,logT,kappa, &
418 : dlnkap_dlnRho,dlnkap_dlnT,ierr)
419 : real(dp), intent(in) :: Z, X, fC, fN, logRho, logT
420 : real(dp), intent(out) :: kappa
421 : real(dp), intent(out) :: dlnkap_dlnRho, dlnkap_dlnT
422 : integer, intent(out) :: ierr
423 0 : real(dp) :: logR, result(3)
424 :
425 : ierr = 0
426 0 : logR = logRho - 3d0*logT + 18d0
427 :
428 0 : call kapCN_interp(Z,X,fC,fN,logR,logT,result,ierr)
429 0 : if(ierr==0)then
430 0 : kappa = exp10(result(1))
431 0 : dlnkap_dlnRho = result(2)
432 0 : dlnkap_dlnT = result(3) - 3*result(2)
433 : else
434 0 : kappa = -1d0
435 0 : dlnkap_dlnRho = 0d0
436 0 : dlnkap_dlnT = 0d0
437 : end if
438 :
439 0 : end subroutine kapCN_get
440 :
441 : end module kapcn
|
