Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2012-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 : module kap_eval_fixed
21 :
22 : use kap_eval_support
23 : use const_def, only: dp, ln10
24 : use math_lib
25 : use utils_lib, only: mesa_error
26 :
27 : implicit none
28 :
29 : private
30 : public :: Get1_kap_fixed_metal_Results
31 :
32 : contains
33 :
34 346574 : subroutine Get1_kap_fixed_metal_Results( &
35 : z_tables, num_Zs, rq, Z, X, Rho, logRho, T, logT, &
36 : logKap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
37 : use kap_def
38 : use const_def, only: dp, ln10
39 :
40 : ! INPUT
41 : type (Kap_Z_Table), dimension(:), pointer :: z_tables
42 : integer, intent(in) :: num_Zs
43 : type (Kap_General_Info), pointer :: rq
44 : real(dp), intent(in) :: Z, X
45 : real(dp), intent(inout) :: Rho, logRho ! can be modified to clip to table boundaries
46 : real(dp), intent(inout) :: T, logT
47 :
48 : ! OUTPUT
49 : real(dp), intent(out) :: logKap, dlnkap_dlnRho, dlnkap_dlnT
50 : integer, intent(out) :: ierr ! 0 means AOK.
51 :
52 : integer :: iz, i
53 86312 : real(dp) :: Z0, Z1, alfa, beta, lnZ, lnZ0, lnZ1
54 : real(dp) :: K0, logK0, dlogK0_dlogRho, dlogK0_dlogT
55 : real(dp) :: logK1, dlogK1_dlogRho, dlogK1_dlogT
56 : real(dp) :: res
57 : character (len=256) :: message
58 :
59 : logical :: dbg
60 :
61 : include 'formats'
62 :
63 86312 : dbg = .false.
64 : if (dbg) write(*,1) 'Get1_kap_fixed_metal_Results logT', logT
65 :
66 86312 : ierr = 0
67 86312 : logKap = 0d0
68 86312 : dlnkap_dlnRho = 0d0
69 86312 : dlnkap_dlnT = 0d0
70 :
71 86312 : if (num_Zs > 1) then
72 86312 : if (z_tables(1)% Z >= z_tables(2)% Z) then
73 0 : ierr = -3
74 36358 : return
75 : end if
76 : end if
77 :
78 86312 : if (num_Zs == 1 .or. Z >= z_tables(num_Zs)% Z) then ! use the largest Z
79 : call Get_Kap_for_X( &
80 : z_tables, rq, num_Zs, X, logRho, logT, &
81 0 : logKap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
82 : if (dbg) write(*,1) 'logKap logT', logKap, logT
83 0 : return
84 : end if
85 :
86 86312 : if (Z <= z_tables(1)% Z) then ! use the smallest Z
87 : if (dbg) then
88 : write(*,*) 'use the smallest Z'
89 : write(*,*) 'Z', Z
90 : write(*,*) 'z_tables(1)% Z', z_tables(1)% Z
91 : end if
92 : call Get_Kap_for_X( &
93 : z_tables, rq, 1, X, logRho, logT, &
94 0 : logKap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
95 0 : return
96 : end if
97 :
98 : if (dbg) then
99 : do iz = 1, num_Zs
100 : write(*,*) 'Z(iz)', iz, z_tables(iz)% Z
101 : end do
102 : end if
103 :
104 646380 : do iz = 1, num_Zs-1
105 646380 : if (Z < z_tables(iz+1)% Z) exit
106 : end do
107 :
108 86312 : Z0 = z_tables(iz)% Z
109 86312 : Z1 = z_tables(iz+1)% Z
110 :
111 : if (dbg) then
112 : write(*,*) 'Z0', Z0
113 : write(*,*) 'Z ', Z
114 : write(*,*) 'Z1', Z1
115 : end if
116 :
117 86312 : if (Z1 <= Z0) then
118 0 : ierr = 1
119 0 : return
120 : end if
121 :
122 86312 : if (Z <= Z0) then ! use the Z0 table
123 : if (dbg) then
124 : write(*,*) 'use the Z0 table', iz, Z0, Z, Z-Z0
125 : do i = 1, z_tables(iz)% num_Xs
126 : write(*,*) 'X', i, z_tables(iz)% x_tables(i)% X
127 : end do
128 : end if
129 : call Get_Kap_for_X( &
130 : z_tables, rq, iz, X, logRho, logT, &
131 36358 : logKap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
132 36358 : return
133 : end if
134 :
135 49954 : if (Z >= Z1) then ! use the Z1 table
136 : if (dbg) write(*,*) 'use the Z1 table', Z1
137 : call Get_Kap_for_X( &
138 : z_tables, rq, iz+1, X, logRho, logT, &
139 0 : logKap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
140 0 : return
141 : end if
142 :
143 49954 : if (num_Zs >= 4 .and. rq% cubic_interpolation_in_Z) then
144 : if (dbg) write(*,*) 'call Get_Kap_for_Z_cubic'
145 : call Get_Kap_for_Z_cubic(z_tables, num_Zs, rq, iz, Z, X, logRho, logT, &
146 0 : logKap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
147 0 : if (ierr /= 0) then
148 : if (dbg) write(*,*) 'failed in Get_Kap_for_Z_cubic'
149 : return
150 : end if
151 : else ! linear
152 : if (dbg) write(*,*) 'call Get_Kap_for_Z_linear'
153 : call Get_Kap_for_Z_linear(z_tables, rq, iz, Z, Z0, Z1, X, logRho, logT, &
154 49954 : logKap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
155 49954 : if (ierr /= 0) then
156 : if (dbg) write(*,*) 'failed in Get_Kap_for_Z_linear'
157 : return
158 : end if
159 : end if
160 :
161 : if (dbg) then
162 : write(*,1) 'final logK at X and Z', logKap, logT, logRho, X, Z
163 : write(*,'(A)')
164 : end if
165 :
166 : end subroutine Get1_kap_fixed_metal_Results
167 :
168 :
169 : ! use tables iz-1 to iz+2 to do piecewise monotonic cubic interpolation in Z
170 0 : subroutine Get_Kap_for_Z_cubic( &
171 : z_tables, num_Zs, rq, iz, Z, X, logRho, logT, &
172 : logK, dlnkap_dlnRho, dlnkap_dlnT, ierr)
173 : use kap_def
174 : use interp_1d_def, only: pm_work_size
175 : use interp_1d_lib, only: interpolate_vector_autodiff, interp_pm_autodiff
176 : use auto_diff
177 :
178 : type (Kap_Z_Table), dimension(:), pointer :: z_tables
179 : type (Kap_General_Info), pointer :: rq
180 : integer, intent(in) :: num_Zs, iz
181 : real(dp), intent(in) :: Z, X
182 : real(dp), intent(inout) :: logK, logRho, logT
183 : real(dp), intent(out) :: dlnkap_dlnRho, dlnkap_dlnT
184 : integer, intent(out) :: ierr
185 :
186 : integer, parameter :: n_old = 4, n_new = 1
187 0 : real(dp), dimension(n_old) :: logKs, dlogKs_dlogRho, dlogKs_dlogT
188 : type(auto_diff_real_2var_order1), dimension(n_old) :: logKs_ad
189 : type(auto_diff_real_2var_order1) :: logK_ad
190 : type(auto_diff_real_2var_order1) :: z_old(n_old), z_new(n_new)
191 : type(auto_diff_real_2var_order1), target :: work_ary(n_old*pm_work_size)
192 : type(auto_diff_real_2var_order1), pointer :: work(:)
193 : integer :: i, i1, izz
194 :
195 : logical, parameter :: dbg = .false.
196 :
197 : 11 format(a40,e20.10)
198 :
199 0 : ierr = 0
200 0 : work => work_ary
201 :
202 0 : if (iz+2 > num_Zs) then
203 0 : i1 = num_Zs-2
204 0 : else if (iz == 1) then
205 : i1 = 2
206 : else
207 0 : i1 = iz
208 : end if
209 :
210 : if (dbg) write(*,*) 'iz', iz
211 :
212 0 : do i=1,n_old
213 0 : izz = i1-2+i
214 : if (dbg) write(*,*) 'izz', izz
215 0 : z_old(i) %val = z_tables(izz)% Z
216 0 : z_old(i) % d1val1 = 0d0
217 0 : z_old(i) % d1val2 = 0d0
218 : call Get_Kap_for_X( &
219 : z_tables, rq, izz, X, logRho, logT, &
220 0 : logKs(i), dlogKs_dlogRho(i), dlogKs_dlogT(i), ierr)
221 0 : if (ierr /= 0) then
222 : if (dbg) write(*,11) 'logRho', logRho
223 : if (dbg) write(*,11) 'logT', logT
224 : return
225 : end if
226 : ! now pack into auto_diff form
227 0 : logKs_ad(i) % val = logKs(i)
228 0 : logKs_ad(i) % d1val1 = dlogKs_dlogT(i)
229 0 : logKs_ad(i) % d1val2 = dlogKs_dlogRho(i)
230 : end do
231 0 : z_new(1) % val = Z
232 0 : z_new(1) % d1val1 = 0d0
233 0 : z_new(1) % d1val2 = 0d0
234 :
235 0 : call interp1(logKs_ad, logK_ad, ierr)
236 0 : if (ierr /= 0) then
237 0 : call mesa_error(__FILE__,__LINE__,'failed in interp1 for logK')
238 0 : return
239 : end if
240 :
241 : ! unpack auto_diff pack into output reals
242 0 : logK = logK_ad % val
243 0 : dlnkap_dlnT = logK_ad % d1val1
244 0 : dlnkap_dlnRho = logK_ad % d1val2
245 :
246 0 : if (dbg) then
247 :
248 : do i=1,n_old
249 : write(*,*) 'z_old(i)', z_old(i)
250 : end do
251 : write(*,*) 'z_new(1)', z_new(1)
252 : write(*,'(A)')
253 :
254 : do i=1,n_old
255 : write(*,*) 'logKs(i)', logKs(i)
256 : end do
257 : write(*,*) 'logK', logK
258 : write(*,'(A)')
259 :
260 : do i=1,n_old
261 : write(*,*) 'dlogKs_dlogRho(i)', dlogKs_dlogRho(i)
262 : end do
263 : write(*,*) 'dlnkap_dlnRho', dlnkap_dlnRho
264 : write(*,'(A)')
265 :
266 : do i=1,n_old
267 : write(*,*) 'dlogKs_dlogT(i)', dlogKs_dlogT(i)
268 : end do
269 : write(*,*) 'dlnkap_dlnT', dlnkap_dlnT
270 : write(*,'(A)')
271 :
272 : end if
273 :
274 : contains
275 :
276 0 : subroutine interp1(old, new, ierr)
277 : type(auto_diff_real_2var_order1), intent(in) :: old(n_old)
278 : type(auto_diff_real_2var_order1), intent(out) :: new
279 : integer, intent(out) :: ierr
280 : type(auto_diff_real_2var_order1) :: v_old(n_old), v_new(n_new)
281 : integer :: i
282 0 : do i = 1, n_old
283 0 : v_old(i) = old(i)
284 : end do
285 : call interpolate_vector_autodiff( &
286 : n_old, z_old, n_new, z_new, v_old, v_new, interp_pm_autodiff, pm_work_size, work, &
287 0 : 'Get_Kap_for_Z_cubic', ierr)
288 0 : new = v_new(1)
289 0 : end subroutine interp1
290 :
291 : end subroutine Get_Kap_for_Z_cubic
292 :
293 :
294 49954 : subroutine Get_Kap_for_Z_linear(z_tables, rq, iz, Z, Z0, Z1, X, logRho, logT, &
295 : logKap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
296 : use kap_def
297 : type (Kap_Z_Table), dimension(:), pointer :: z_tables
298 : type (Kap_General_Info), pointer :: rq
299 : integer, intent(in) :: iz
300 : real(dp), intent(in) :: Z, Z0, Z1, X
301 : real(dp), intent(inout) :: logRho, logT
302 : real(dp), intent(out) :: logKap, dlnkap_dlnRho, dlnkap_dlnT
303 : integer, intent(out) :: ierr
304 :
305 : real(dp) :: logK0, dlogK0_dlogRho, dlogK0_dlogT
306 49954 : real(dp) :: logK1, dlogK1_dlogRho, dlogK1_dlogT
307 49954 : real(dp) :: alfa, beta
308 :
309 : logical, parameter :: dbg = .false.
310 :
311 : ierr = 0
312 : call Get_Kap_for_X( &
313 : z_tables, rq, iz, X, logRho, logT, &
314 49954 : logK0, dlogK0_dlogRho, dlogK0_dlogT, ierr)
315 49954 : if (ierr /= 0) return
316 : if (dbg) write(*,*) 'logK0', logK0
317 :
318 : call Get_Kap_for_X( &
319 : z_tables, rq, iz+1, X, logRho, logT, &
320 49954 : logK1, dlogK1_dlogRho, dlogK1_dlogT, ierr)
321 49954 : if (ierr /= 0) return
322 : if (dbg) write(*,*) 'logK1', logK1
323 :
324 : ! Z0 result in logK0, Z1 result in logK1
325 49954 : beta = (Z - Z1) / (Z0 - Z1) ! beta -> 1 as Z -> Z0
326 49954 : alfa = 1d0 - beta
327 :
328 49954 : logKap = beta * logK0 + alfa * logK1
329 49954 : dlnkap_dlnRho = beta * dlogK0_dlogRho + alfa * dlogK1_dlogRho
330 49954 : dlnkap_dlnT = beta * dlogK0_dlogT + alfa * dlogK1_dlogT
331 :
332 : end subroutine Get_Kap_for_Z_linear
333 :
334 :
335 136266 : subroutine Get_Kap_for_X(z_tables, rq, iz, X, logRho, logT, &
336 : logK, dlogK_dlogRho, dlogK_dlogT, ierr)
337 : use kap_def
338 : type (Kap_Z_Table), dimension(:), pointer :: z_tables
339 : type (Kap_General_Info), pointer :: rq
340 : integer, intent(in) :: iz
341 : real(dp), intent(in) :: X
342 : real(dp), intent(inout) :: logRho, logT
343 : real(dp), intent(out) :: logK, dlogK_dlogRho, dlogK_dlogT
344 : integer, intent(out) :: ierr
345 :
346 : type (Kap_X_Table), dimension(:), pointer :: x_tables
347 : real(dp) :: logK0, dlogK0_dlogRho, dlogK0_dlogT
348 : real(dp) :: logK1, dlogK1_dlogRho, dlogK1_dlogT
349 136266 : real(dp) :: X0, X1
350 : integer :: ix, i, num_Xs
351 :
352 : logical, parameter :: dbg = .false.
353 :
354 : include 'formats'
355 :
356 136266 : ierr = 0
357 136266 : x_tables => z_tables(iz)% x_tables
358 136266 : num_Xs = z_tables(iz)% num_Xs
359 :
360 136266 : if (X < 0 .or. X > 1) then
361 0 : ierr = -3
362 12270 : return
363 : end if
364 :
365 136266 : if (num_Xs > 1) then
366 136234 : if (x_tables(1)% X >= x_tables(2)% X) then
367 0 : ierr = -3
368 0 : write(*,*) 'x_tables must have increasing X values for Get_Kap_for_X'
369 0 : write(*,*) 'lowT', z_tables(iz)% lowT_flag
370 0 : write(*,2) 'Z', iz, z_tables(iz)% Z
371 0 : write(*,2) 'num_Xs', num_Xs
372 0 : do i=1,num_Xs
373 0 : write(*,2) 'X', i, x_tables(i)% X
374 : end do
375 0 : stop
376 :
377 : return
378 : end if
379 : end if
380 :
381 136266 : if (num_Xs == 1 .or. X <= x_tables(1)% X) then ! use the first X
382 : call Get_Kap_for_logRho_logT( &
383 : z_tables, rq, iz, x_tables, 1, &
384 32 : logRho, logT, logK, dlogK_dlogRho, dlogK_dlogT, ierr)
385 32 : return
386 : end if
387 :
388 136234 : if (X >= x_tables(num_Xs)% X) then ! use the last X
389 : if (dbg) write(*,*) 'use the last X: call Get_Kap_for_logRho_logT'
390 : call Get_Kap_for_logRho_logT( &
391 : z_tables, rq, iz, x_tables, num_Xs, &
392 0 : logRho, logT, logK, dlogK_dlogRho, dlogK_dlogT, ierr)
393 0 : return
394 : end if
395 :
396 : ! search for the X
397 : !write(*,*) 'search for the X'
398 136234 : ix = num_Xs
399 693408 : do i = 1, num_Xs-1
400 693408 : if (X < x_tables(i+1)% X) then
401 136234 : ix = i; exit
402 : end if
403 : end do
404 :
405 136234 : if (ix == num_Xs) then
406 : if (dbg) write(*,*) 'ix == num_Xs: call Get_Kap_for_logRho_logT'
407 : call Get_Kap_for_logRho_logT( &
408 : z_tables, rq, iz, x_tables, num_Xs, &
409 0 : logRho, logT, logK, dlogK_dlogRho, dlogK_dlogT, ierr)
410 0 : return
411 : end if
412 :
413 136234 : X0 = x_tables(ix)% X
414 136234 : X1 = x_tables(ix+1)% X
415 :
416 136234 : if (X1 <= X0) then
417 0 : ierr = 1
418 0 : return
419 : end if
420 :
421 136234 : if (X0 >= X) then ! use the X0 table
422 : if (dbg) write(*,*) 'use the X0 table'
423 : call Get_Kap_for_logRho_logT( &
424 : z_tables, rq, iz, x_tables, ix, &
425 12238 : logRho, logT, logK, dlogK_dlogRho, dlogK_dlogT, ierr)
426 12238 : return
427 : end if
428 :
429 123996 : if (X1 <= X) then ! use the X1 table
430 : if (dbg) write(*,*) 'use the X1 table'
431 : call Get_Kap_for_logRho_logT( &
432 : z_tables, rq, iz, x_tables, ix+1, &
433 0 : logRho, logT, logK, dlogK_dlogRho, dlogK_dlogT, ierr)
434 0 : return
435 : end if
436 :
437 123996 : if (num_Xs >= 4 .and. rq% cubic_interpolation_in_X) then
438 : !write(*,*) 'call Get_Kap_for_X_cubic'
439 : call Get_Kap_for_X_cubic( &
440 : z_tables, rq, iz, ix, num_Xs, x_tables, X, logRho, logT, &
441 0 : logK, dlogK_dlogRho, dlogK_dlogT, ierr)
442 0 : if (ierr /= 0) then
443 : !write(*,*) 'failed in Get_Kap_for_X_cubic'
444 : return
445 : end if
446 : else ! linear
447 : !write(*,*) 'call Get_Kap_for_X_linear'
448 : call Get_Kap_for_X_linear( &
449 : z_tables, rq, iz, ix, x_tables, X, X0, X1, logRho, logT, &
450 123996 : logK, dlogK_dlogRho, dlogK_dlogT, ierr)
451 123996 : if (ierr /= 0) then
452 : !write(*,*) 'failed in Get_Kap_for_X_linear'
453 : return
454 : end if
455 : end if
456 :
457 : if (.false.) then
458 : write(*,1) 'logK at X for Z', logK, logT, logRho, X, z_tables(iz)% Z
459 : write(*,'(A)')
460 : end if
461 :
462 : end subroutine Get_Kap_for_X
463 :
464 :
465 : ! use tables ix-1 to ix+2 to do piecewise monotonic cubic interpolation in X
466 0 : subroutine Get_Kap_for_X_cubic( &
467 : z_tables, rq, iz, ix, num_Xs, x_tables, X, logRho, logT, &
468 : logK, dlogK_dlogRho, dlogK_dlogT, ierr)
469 : use kap_def
470 : use interp_1d_def, only: pm_work_size
471 : use interp_1d_lib, only: interpolate_vector_autodiff, interp_pm_autodiff
472 : use auto_diff
473 :
474 : type (Kap_Z_Table), dimension(:), pointer :: z_tables
475 : type (Kap_General_Info), pointer :: rq
476 : integer, intent(in) :: iz, ix, num_Xs
477 : type (Kap_X_Table), dimension(:), pointer :: x_tables
478 : real(dp), intent(in) :: X
479 : real(dp), intent(inout) :: logRho, logT
480 : real(dp), intent(out) :: logK, dlogK_dlogRho, dlogK_dlogT
481 : integer, intent(out) :: ierr
482 :
483 : integer, parameter :: n_old = 4, n_new = 1
484 0 : real(dp), dimension(n_old) :: logKs, dlogKs_dlogRho, dlogKs_dlogT
485 : type(auto_diff_real_2var_order1), dimension(n_old) :: logKs_ad
486 : type(auto_diff_real_2var_order1) :: logK_ad
487 : type(auto_diff_real_2var_order1) :: x_old(n_old), x_new(n_new)
488 : type(auto_diff_real_2var_order1), target :: work_ary(n_old*pm_work_size)
489 : type(auto_diff_real_2var_order1), pointer :: work(:)
490 : integer :: i, i1, ixx
491 :
492 : logical, parameter :: dbg = .false.
493 :
494 : 11 format(a40,e20.10)
495 :
496 0 : ierr = 0
497 0 : work => work_ary
498 :
499 0 : if (ix+2 > num_Xs) then
500 0 : i1 = num_Xs-2
501 0 : else if (ix == 1) then
502 : i1 = 2
503 : else
504 0 : i1 = ix
505 : end if
506 :
507 : if (dbg) write(*,*) 'ix', ix
508 :
509 0 : do i=1,n_old
510 0 : ixx = i1-2+i
511 : if (dbg) write(*,*) 'ixx', ixx
512 0 : x_old(i) % val = x_tables(ixx)% X
513 0 : x_old(i) % d1val1 = 0d0
514 0 : x_old(i) % d1val2 = 0d0
515 : call Get_Kap_for_logRho_logT( &
516 : z_tables, rq, iz, x_tables, ixx, &
517 0 : logRho, logT, logKs(i), dlogKs_dlogRho(i), dlogKs_dlogT(i), ierr)
518 0 : if (ierr /= 0) then
519 : if (dbg) write(*,11) 'logRho', logRho
520 : if (dbg) write(*,11) 'logT', logT
521 : return
522 : end if
523 : ! now pack into auto_diff form
524 0 : logKs_ad(i) % val = logKs(i)
525 0 : logKs_ad(i) % d1val1 = dlogKs_dlogT(i)
526 0 : logKs_ad(i) % d1val2 = dlogKs_dlogRho(i)
527 : end do
528 0 : x_new(1) % val = X
529 0 : x_new(1) % d1val1 = 0d0
530 0 : x_new(1) % d1val2 = 0d0
531 :
532 0 : call interp1(logKs_ad, logK_ad, ierr)
533 0 : if (ierr /= 0) then
534 0 : call mesa_error(__FILE__,__LINE__,'failed in interp1 for logK')
535 0 : return
536 : end if
537 :
538 : ! unpack auto_diff pack into output reals
539 0 : logK = logK_ad % val
540 0 : dlogK_dlogT = logK_ad % d1val1
541 0 : dlogK_dlogRho = logK_ad % d1val2
542 :
543 0 : if (dbg) then
544 :
545 : do i=1,n_old
546 : write(*,*) 'x_old(i)', x_old(i)
547 : end do
548 : write(*,*) 'x_new(1)', x_new(1)
549 : write(*,'(A)')
550 :
551 : do i=1,n_old
552 : write(*,*) 'logKs(i)', logKs(i)
553 : end do
554 : write(*,*) 'logK', logK
555 : write(*,'(A)')
556 :
557 : do i=1,n_old
558 : write(*,*) 'dlogKs_dlogRho(i)', dlogKs_dlogRho(i)
559 : end do
560 : write(*,*) 'dlogK_dlogRho', dlogK_dlogRho
561 : write(*,'(A)')
562 :
563 : do i=1,n_old
564 : write(*,*) 'dlogKs_dlogT(i)', dlogKs_dlogT(i)
565 : end do
566 : write(*,*) 'dlogK_dlogT', dlogK_dlogT
567 : write(*,'(A)')
568 :
569 : end if
570 :
571 : contains
572 :
573 0 : subroutine interp1(old, new, ierr)
574 : type(auto_diff_real_2var_order1), intent(in) :: old(n_old)
575 : type(auto_diff_real_2var_order1), intent(out) :: new
576 : integer, intent(out) :: ierr
577 : type(auto_diff_real_2var_order1) :: v_old(n_old), v_new(n_new)
578 : integer :: i
579 0 : do i = 1, n_old
580 0 : v_old(i) = old(i)
581 : end do
582 : call interpolate_vector_autodiff( &
583 : n_old, x_old, n_new, x_new, v_old, v_new, interp_pm_autodiff, pm_work_size, work, &
584 0 : 'Get_Kap_for_X_cubic', ierr)
585 0 : new = v_new(1)
586 0 : end subroutine interp1
587 :
588 : end subroutine Get_Kap_for_X_cubic
589 :
590 :
591 : ! use tables ix and ix+1 to do linear interpolation in X
592 123996 : subroutine Get_Kap_for_X_linear( &
593 : z_tables, rq, iz, ix, x_tables, X, X0, X1, logRho, logT, &
594 : logK, dlogK_dlogRho, dlogK_dlogT, ierr)
595 : use kap_def
596 : type (Kap_Z_Table), dimension(:), pointer :: z_tables
597 : type (Kap_General_Info), pointer :: rq
598 : integer, intent(in) :: iz, ix
599 : type (Kap_X_Table), dimension(:), pointer :: x_tables
600 : real(dp), intent(in) :: X, X0, X1
601 : real(dp), intent(inout) :: logRho, logT
602 : real(dp), intent(out) :: logK, dlogK_dlogRho, dlogK_dlogT
603 : integer, intent(out) :: ierr
604 :
605 : real(dp) :: logK0, dlogK0_dlogRho, dlogK0_dlogT
606 123996 : real(dp) :: logK1, dlogK1_dlogRho, dlogK1_dlogT
607 123996 : real(dp) :: alfa, beta
608 :
609 : ierr = 0
610 :
611 : call Get_Kap_for_logRho_logT( &
612 : z_tables, rq, iz, x_tables, ix, &
613 123996 : logRho, logT, logK0, dlogK0_dlogRho, dlogK0_dlogT, ierr)
614 123996 : if (ierr /= 0) return
615 :
616 : call Get_Kap_for_logRho_logT( &
617 : z_tables, rq, iz, x_tables, ix+1, &
618 123996 : logRho, logT, logK1, dlogK1_dlogRho, dlogK1_dlogT, ierr)
619 123996 : if (ierr /= 0) return
620 :
621 : ! X0 result in logK0, X1 result in logK1
622 123996 : beta = (X - X1) / (X0 - X1) ! beta -> 1 as X -> X0
623 123996 : alfa = 1d0 - beta
624 :
625 123996 : logK = beta * logK0 + alfa * logK1
626 123996 : dlogK_dlogRho = beta * dlogK0_dlogRho + alfa * dlogK1_dlogRho
627 123996 : dlogK_dlogT = beta * dlogK0_dlogT + alfa * dlogK1_dlogT
628 :
629 : end subroutine Get_Kap_for_X_linear
630 :
631 :
632 260262 : subroutine Get_Kap_for_logRho_logT( &
633 : z_tables, rq, iz, x_tables, ix, &
634 : logRho, logT, logK, dlogK_dlogRho, dlogK_dlogT, ierr)
635 : use load_kap, only: load_one
636 : use kap_def
637 : type (Kap_Z_Table), dimension(:), pointer :: z_tables
638 : type (Kap_General_Info), pointer :: rq
639 : type (Kap_X_Table), dimension(:), pointer :: x_tables
640 : integer :: iz, ix
641 : real(dp), intent(inout) :: logRho, logT
642 : real(dp), intent(out) :: logK, dlogK_dlogRho, dlogK_dlogT
643 : integer, intent(out) :: ierr
644 :
645 260262 : real(dp) :: logR0, logR1, logT0, logT1, logR, logR_in
646 260262 : real(dp) :: df_dx, df_dy
647 : integer :: iR, jtemp, i, num_logRs, num_logTs
648 : logical :: clipped_logR
649 : logical, parameter :: read_later = .false.
650 :
651 : logical, parameter :: dbg = .false.
652 :
653 : include 'formats'
654 :
655 260262 : ierr = 0
656 260262 : if (x_tables(ix)% not_loaded_yet) then ! avoid doing critical section if possible
657 30 : !$omp critical (load_kap_x_table)
658 30 : if (x_tables(ix)% not_loaded_yet) then
659 : call load_one(rq, &
660 : z_tables, iz, ix, x_tables(ix)% X, x_tables(ix)% Z, &
661 18 : read_later, ierr)
662 : end if
663 : !$omp end critical (load_kap_x_table)
664 : end if
665 260262 : if (ierr /= 0) then
666 : !call mesa_error(__FILE__,__LINE__,'load_one failed in Get_Kap_for_logRho_logT')
667 : return
668 : end if
669 :
670 : ! logR from inputs
671 260262 : logR_in = logRho - 3d0*logT + 18d0
672 :
673 : ! blends at higher levels MUST prevent
674 : ! these tables from being called off their
675 : ! high/low T and low R edges
676 :
677 260262 : if (logT > x_tables(ix)% logT_max) then
678 0 : ierr = -1
679 0 : return
680 : end if
681 :
682 260262 : if (logT < x_tables(ix)% logT_min) then
683 0 : ierr = -1
684 0 : return
685 : end if
686 :
687 260262 : if (logR_in < x_tables(ix)% logR_min) then
688 0 : ierr = -1
689 0 : return
690 : end if
691 :
692 :
693 : ! off the high R edge, we use the input temperature
694 : ! but clip logR to the table edge value
695 :
696 260262 : if (logR_in > x_tables(ix)% logR_max) then
697 0 : logR = x_tables(ix)% logR_max
698 0 : clipped_logR = .true.
699 : else
700 260262 : logR = logR_in
701 260262 : clipped_logR = .false.
702 : end if
703 :
704 :
705 260262 : num_logRs = x_tables(ix)% num_logRs
706 260262 : num_logTs = x_tables(ix)% num_logTs
707 :
708 260262 : if (num_logRs <= 0) then
709 0 : write(*,*) 'num_logRs', num_logRs
710 0 : write(*,*) 'ix', ix
711 0 : call mesa_error(__FILE__,__LINE__,'Get_Kap_for_logRho_logT')
712 : end if
713 :
714 260262 : if (num_logTs <= 0) then
715 0 : write(*,*) 'num_logTs', num_logRs
716 0 : write(*,*) 'ix', ix
717 0 : call mesa_error(__FILE__,__LINE__,'Get_Kap_for_logRho_logT')
718 : end if
719 :
720 : call Locate_logR( &
721 : rq, num_logRs, x_tables(ix)% logR_min, x_tables(ix)% logR_max, &
722 260262 : x_tables(ix)% ili_logRs, x_tables(ix)% logRs, logR, iR, logR0, logR1, ierr)
723 260262 : if (ierr /= 0) then
724 0 : write(*,1) 'x_tables(ix)% logR_min', x_tables(ix)% logR_min
725 0 : write(*,1) 'x_tables(ix)% logR_max', x_tables(ix)% logR_max
726 0 : write(*,2) 'num_logRs', num_logRs
727 0 : write(*,2) 'iR', iR
728 0 : write(*,1) 'logR', logR
729 0 : write(*,1) 'logR0', logR0
730 0 : write(*,1) 'logR1', logR1
731 0 : do i=1,num_logRs
732 0 : write(*,2) 'logR', i, x_tables(ix)% logRs(i)
733 : end do
734 0 : write(*,*) 'clip_to_kap_table_boundaries', clip_to_kap_table_boundaries
735 0 : call mesa_error(__FILE__,__LINE__,'failed in Locate_logR')
736 0 : return
737 : end if
738 :
739 : call Locate_logT( &
740 : rq, num_logTs, x_tables(ix)% logT_min, x_tables(ix)% logT_max, &
741 260262 : x_tables(ix)% ili_logTs, x_tables(ix)% logTs, logT, jtemp, logT0, logT1, ierr)
742 260262 : if (ierr /= 0) return
743 :
744 : call Do_Kap_Interpolations( &
745 : x_tables(ix)% kap1, num_logRs, num_logTs, &
746 : iR, jtemp, logR0, logR, logR1, logT0, logT, logT1, &
747 260262 : logK, df_dx, df_dy)
748 260262 : if (clipped_logR) df_dx = 0
749 :
750 : if (dbg) write(*,1) 'Do_Kap_Interpolations: logK', logK
751 :
752 : ! convert df_dx and df_dy to dlogK_dlogRho, dlogK_dlogT
753 260262 : dlogK_dlogRho = df_dx
754 260262 : dlogK_dlogT = df_dy - 3d0*df_dx
755 :
756 : if (dbg) then
757 : write(*,1) 'logK', logK, logT, logRho, x_tables(ix)% X, z_tables(iz)% Z
758 : end if
759 :
760 260262 : end subroutine Get_Kap_for_logRho_logT
761 :
762 : end module kap_eval_fixed
|
