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 : module kap_eval_co
21 : use utils_lib,only: is_bad, mesa_error
22 : use kap_eval_support
23 : use const_def, only: dp
24 : use math_lib
25 :
26 : implicit none
27 :
28 : private
29 : public :: Get1_kap_CO_Results
30 :
31 : contains
32 :
33 2 : subroutine Get1_kap_CO_Results( &
34 : rq, Zbase, X, dXC, dXO, Rho, logRho, T, logT, &
35 : logKap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
36 : use kap_def
37 : use const_def, only: dp
38 :
39 : ! INPUT
40 : type (Kap_General_Info), pointer :: rq
41 : real(dp), intent(in) :: Zbase, X, dXC, dXO
42 : real(dp), intent(inout) :: Rho, logRho ! can be modified to clip to table boundaries
43 : real(dp), intent(inout) :: T, logT
44 :
45 : ! OUTPUT
46 : real(dp), intent(out) :: logKap, dlnkap_dlnRho, dlnkap_dlnT
47 : integer, intent(out) :: ierr ! 0 means AOK.
48 :
49 : integer :: iz, use_iz, num_Zs, CO_option
50 2 : real(dp) :: Z0, Z1,log10_Zbase, log10_Z0, log10_Z1
51 : real(dp) :: alfa, beta
52 : real(dp) :: logK0, dlogK0_dlogRho, dlogK0_dlogT, logK1, dlogK1_dlogRho, dlogK1_dlogT
53 : character (len=256) :: message
54 :
55 : logical, parameter :: use_closest_Z = .false.
56 :
57 : logical, parameter :: dbg = .false.
58 :
59 : include 'formats'
60 :
61 2 : CO_option = rq% kap_CO_option
62 2 : num_Zs = num_kap_CO_Zs(CO_option)
63 :
64 2 : if (num_Zs > 1) then
65 2 : if (kap_co_z_tables(CO_option)% ar(1)% Zbase >= kap_co_z_tables(CO_option)% ar(2)% Zbase) then
66 0 : ierr = -3
67 0 : return
68 : end if
69 : end if
70 :
71 2 : if (num_Zs == 1 .or. &
72 : Zbase >= kap_co_z_tables(CO_option)% ar(num_Zs)% Zbase) then ! use the largest Zbase
73 : if (dbg) write(*,*) 'use the largest Zbase', &
74 : num_kap_CO_Zs(CO_option), kap_co_z_tables(CO_option)% ar(num_Zs)% Zbase
75 : call Get_Kap_for_CO_X( &
76 : rq, dXC, dXO, num_Zs, X, logRho, logT, &
77 0 : logKap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
78 0 : return
79 : end if
80 :
81 2 : if (Zbase <= kap_co_z_tables(CO_option)% ar(1)% Zbase) then ! use the smallest Zbase
82 : if (dbg) then
83 : write(*,*) 'use the smallest Zbase'
84 : write(*,*) 'Zbase', Zbase
85 : write(*,*) 'kap_co_z_tables(1)% Zbase', kap_co_z_tables(CO_option)% ar(1)% Zbase
86 : end if
87 : call Get_Kap_for_CO_X( &
88 : rq, dXC, dXO, 1, X, logRho, logT, &
89 0 : logKap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
90 0 : return
91 : end if
92 :
93 8 : do iz = 1, num_Zs-1
94 8 : if (Zbase < kap_co_z_tables(CO_option)% ar(iz+1)% Zbase) exit
95 : end do
96 :
97 2 : Z0 = kap_co_z_tables(CO_option)% ar(iz)% Zbase
98 2 : Z1 = kap_co_z_tables(CO_option)% ar(iz+1)% Zbase
99 :
100 2 : if (Zbase <= Z0) then ! use the Z0 table
101 : if (dbg) write(*,*) 'use the Z0 table', Z0
102 : call Get_Kap_for_CO_X( &
103 : rq, dXC, dXO, iz, X, logRho, logT, &
104 0 : logKap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
105 0 : return
106 : end if
107 :
108 2 : if (Zbase >= Z1) then ! use the Z1 table
109 : if (dbg) write(*,*) 'use the Z1 table', Z1
110 : call Get_Kap_for_CO_X( &
111 : rq, dXC, dXO, iz+1, X, logRho, logT, &
112 0 : logKap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
113 0 : return
114 : end if
115 :
116 : if (use_closest_Z) then
117 : log10_Z0 = kap_co_z_tables(CO_option)% ar(iz)% log10_Zbase
118 : log10_Z1 = kap_co_z_tables(CO_option)% ar(iz+1)% log10_Zbase
119 : log10_Zbase = log10(dble(Zbase))
120 : if (log10_Z1 - log10_Zbase > log10_Zbase - log10_Z0) then ! use the Z0 table
121 : use_iz = iz
122 : else
123 : use_iz = iz+1
124 : end if
125 : if (dbg) write(*,*) 'use the Z0 table', Z0
126 : call Get_Kap_for_CO_X( &
127 : rq, dXC, dXO, use_iz, X, logRho, logT, &
128 : logKap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
129 : return
130 : end if
131 :
132 : if (dbg) then
133 : write(*,*) 'iz', iz
134 : write(*,*) ' Z0', Z0
135 : write(*,*) 'Zbase', Zbase
136 : write(*,*) ' Z1', Z1
137 : write(*,'(A)')
138 : end if
139 :
140 2 : if (num_Zs >= 4 .and. rq% cubic_interpolation_in_Z) then
141 : if (dbg) write(*,*) 'call Get_Kap_for_CO_Z_cubic'
142 : call Get_Kap_for_CO_Z_cubic(rq, iz, Zbase, X, dXC, dXO, logRho, logT, &
143 0 : logKap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
144 0 : if (ierr /= 0) then
145 : if (dbg) write(*,*) 'failed in Get_Kap_for_CO_Z_cubic'
146 : return
147 : end if
148 : else ! linear
149 : if (dbg) write(*,*) 'call Get_Kap_for_CO_Z_linear'
150 : call Get_Kap_for_CO_Z_linear(rq, iz, Zbase, Z0, Z1, X, dXC, dXO, logRho, logT, &
151 2 : logKap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
152 2 : if (ierr /= 0) then
153 : if (dbg) write(*,*) 'failed in Get_Kap_for_CO_Z_linear'
154 : return
155 : end if
156 : end if
157 :
158 : end subroutine Get1_kap_CO_Results
159 :
160 :
161 : ! use tables iz-1 to iz+2 to do piecewise monotonic cubic interpolation in Z
162 0 : subroutine Get_Kap_for_CO_Z_cubic(rq, iz, Z, X, dXC, dXO, logRho, logT, &
163 : logK, dlnkap_dlnRho, dlnkap_dlnT, ierr)
164 : use kap_def
165 : use interp_1d_def, only: pm_work_size
166 : use interp_1d_lib, only: interpolate_vector_autodiff, interp_pm_autodiff
167 : use auto_diff
168 :
169 : type (Kap_General_Info), pointer :: rq
170 : integer, intent(in) :: iz
171 : real(dp), intent(in) :: Z, X, dXC, dXO
172 : real(dp), intent(inout) :: logRho, logT
173 : real(dp), intent(out) :: logK, dlnkap_dlnRho, dlnkap_dlnT
174 : integer, intent(out) :: ierr
175 :
176 : integer, parameter :: n_old = 4, n_new = 1
177 0 : real(dp), dimension(n_old) :: logKs, dlogKs_dlogRho, dlogKs_dlogT
178 : type(auto_diff_real_2var_order1), dimension(n_old) :: logKs_ad
179 : type(auto_diff_real_2var_order1) :: logK_ad
180 : type(auto_diff_real_2var_order1) :: z_old(n_old), z_new(n_new)
181 : type(auto_diff_real_2var_order1), target :: work_ary(n_old*pm_work_size)
182 : type(auto_diff_real_2var_order1), pointer :: work(:)
183 : integer :: i, i1, izz, num_Zs, CO_option
184 :
185 : logical, parameter :: dbg = .false.
186 :
187 : 11 format(a40,e20.10)
188 :
189 0 : ierr = 0
190 0 : work => work_ary
191 0 : CO_option = rq% kap_CO_option
192 0 : num_Zs = num_kap_CO_Zs(CO_option)
193 :
194 0 : if (iz+2 > num_Zs) then
195 0 : i1 = num_Zs-2
196 0 : else if (iz == 1) then
197 : i1 = 2
198 : else
199 0 : i1 = iz
200 : end if
201 :
202 : if (dbg) then
203 : write(*,*) 'n_old', n_old
204 : write(*,*) 'i1', i1
205 : write(*,*) 'iz', iz
206 : write(*,*) 'Z', Z
207 : write(*,'(A)')
208 : end if
209 :
210 0 : do i=1,n_old
211 0 : izz = i1-2+i
212 0 : z_old(i) %val= kap_co_z_tables(CO_option)% ar(izz)% Zbase
213 0 : z_old(i) % d1val1 = 0d0
214 0 : z_old(i) % d1val2 = 0d0
215 : if (dbg) then
216 : write(*,*) 'izz', izz
217 : write(*,*) 'z_old', i, z_old(i)
218 : end if
219 : call Get_Kap_for_CO_X( &
220 : rq, dXC, dXO, izz, X, logRho, logT, &
221 0 : logKs(i), dlogKs_dlogRho(i), dlogKs_dlogT(i), ierr)
222 : if (dbg) write(*,11) 'logK', logKs(i)
223 0 : if (ierr /= 0) then
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(*,'(A)')
252 : write(*,*) 'z_new(1)', z_new(1)
253 : write(*,'(A)')
254 :
255 : do i=1,n_old
256 : write(*,*) 'logK', i, logKs(i)
257 : end do
258 : write(*,'(A)')
259 : write(*,*) 'logK', logK
260 : write(*,'(A)')
261 :
262 : do i=1,n_old
263 : write(*,*) 'dlogKs_dlogRho', i, dlogKs_dlogRho(i)
264 : end do
265 : write(*,'(A)')
266 : write(*,*) 'dlnkap_dlnRho', dlnkap_dlnRho
267 : write(*,'(A)')
268 :
269 : do i=1,n_old
270 : write(*,*) 'dlogKs_dlogT', i, dlogKs_dlogT(i)
271 : end do
272 : write(*,'(A)')
273 : write(*,*) 'dlnkap_dlnT', dlnkap_dlnT
274 : write(*,'(A)')
275 :
276 : end if
277 :
278 : contains
279 :
280 0 : subroutine interp1(old, new, ierr)
281 : type(auto_diff_real_2var_order1), intent(in) :: old(n_old)
282 : type(auto_diff_real_2var_order1), intent(out) :: new
283 : integer, intent(out) :: ierr
284 : type(auto_diff_real_2var_order1) :: v_old(n_old), v_new(n_new)
285 : integer :: i
286 0 : do i = 1, n_old
287 0 : v_old(i) = old(i)
288 : end do
289 : call interpolate_vector_autodiff( &
290 : n_old, z_old, n_new, z_new, v_old, v_new, interp_pm_autodiff, pm_work_size, work, &
291 0 : 'Get_Kap_for_CO_Z_cubic', ierr)
292 0 : new = v_new(1)
293 0 : end subroutine interp1
294 :
295 : end subroutine Get_Kap_for_CO_Z_cubic
296 :
297 :
298 2 : subroutine Get_Kap_for_CO_Z_linear( &
299 : rq, iz, Z, Z0, Z1, X, dXC, dXO, logRho, logT, &
300 : logKap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
301 : use kap_def
302 : type (Kap_General_Info), pointer :: rq
303 : integer, intent(in) :: iz
304 : real(dp), intent(in) :: Z, Z0, Z1, X, dXC, dXO
305 : real(dp), intent(inout) :: logRho, logT
306 : real(dp), intent(out) :: logKap, dlnkap_dlnRho, dlnkap_dlnT
307 : integer, intent(out) :: ierr
308 :
309 2 : real(dp) :: logK0, dlogK0_dlogRho, dlogK0_dlogT, logK1, dlogK1_dlogRho, dlogK1_dlogT
310 2 : real(dp) :: alfa, beta
311 :
312 : logical, parameter :: dbg = .false.
313 :
314 : ierr = 0
315 : if (dbg) write(*,*) 'call Get_Kap_for_CO_X'
316 : call Get_Kap_for_CO_X( &
317 : rq, dXC, dXO, iz, X, logRho, logT, &
318 2 : logK0, dlogK0_dlogRho, dlogK0_dlogT, ierr)
319 2 : if (ierr /= 0) return
320 : if (dbg) write(*,*) 'logK0', logK0
321 :
322 : call Get_Kap_for_CO_X( &
323 : rq, dXC, dXO, iz+1, X, logRho, logT, &
324 2 : logK1, dlogK1_dlogRho, dlogK1_dlogT, ierr)
325 2 : if (ierr /= 0) return
326 : if (dbg) write(*,*) 'logK1', logK1
327 :
328 : ! Z0 result in logK0, Z1 result in logK1
329 2 : beta = (Z - Z1) / (Z0 - Z1) ! beta -> 1 as Z -> Z0
330 2 : alfa = 1d0 - beta
331 2 : logKap = beta*logK0 + alfa*logK1
332 2 : dlnkap_dlnRho = beta*dlogK0_dlogRho + alfa*dlogK1_dlogRho
333 2 : dlnkap_dlnT = beta*dlogK0_dlogT + alfa*dlogK1_dlogT
334 :
335 : end subroutine Get_Kap_for_CO_Z_linear
336 :
337 :
338 4 : subroutine Get_Kap_for_CO_X(rq, dXC, dXO, iz, X, logRho, logT, &
339 : logK, dlogK_dlogRho, dlogK_dlogT, ierr)
340 : use kap_def
341 : ! return opacity from Z table number iz for the given X, dXC, dXO
342 : type (Kap_General_Info), pointer :: rq
343 : integer, intent(in) :: iz
344 : real(dp), intent(in) :: X, dXC, dXO
345 : real(dp), intent(inout) :: logRho, logT
346 : real(dp), intent(out) :: logK, dlogK_dlogRho, dlogK_dlogT
347 : integer, intent(out) :: ierr
348 :
349 : type (Kap_CO_X_Table), dimension(:), pointer :: x_tables
350 : real(dp) :: logK0, dlogK0_dlogRho, dlogK0_dlogT, logK1, dlogK1_dlogRho, dlogK1_dlogT
351 4 : real(dp) :: X0, X1
352 : real(dp) :: alfa, beta
353 : integer :: ix, i, num_Xs, CO_option
354 :
355 : logical, parameter :: dbg = .false.
356 :
357 4 : CO_option = rq% kap_CO_option
358 4 : num_Xs = num_kap_CO_Xs(CO_option)
359 4 : x_tables => kap_co_z_tables(CO_option)% ar(iz)% x_tables
360 :
361 4 : if (X < 0 .or. X > 1) then
362 0 : ierr = -3
363 4 : return
364 : end if
365 :
366 4 : if (num_Xs > 1) then
367 4 : if (x_tables(1)% X >= x_tables(2)% X) then
368 0 : ierr = -3
369 0 : return
370 : end if
371 : end if
372 :
373 4 : if (X >= x_tables(num_Xs)% X) then ! use the last X
374 : if (dbg) write(*,*) 'use the last X'
375 : call Get_Kap_for_dXCO( &
376 : rq, iz, x_tables, dXC, dXO, num_Xs, &
377 0 : logRho, logT, logK, dlogK_dlogRho, dlogK_dlogT, ierr)
378 0 : return
379 : end if
380 :
381 4 : if (X <= x_tables(1)% X) then ! use the first X
382 : if (dbg) write(*,*) 'use the first X'
383 : call Get_Kap_for_dXCO( &
384 : rq, iz, x_tables, dXC, dXO, 1, &
385 4 : logRho, logT, logK, dlogK_dlogRho, dlogK_dlogT, ierr)
386 4 : return
387 : end if
388 :
389 : ! search for the X
390 : if (dbg) write(*,*) 'search for the X'
391 0 : ix = num_Xs
392 0 : do i = 1, num_Xs-1
393 0 : if (X < x_tables(i+1)% X) then
394 0 : ix = i; exit
395 : end if
396 : end do
397 :
398 0 : if (ix == num_Xs) then
399 : call Get_Kap_for_dXCO( &
400 : rq, iz, x_tables, dXC, dXO, num_Xs, &
401 0 : logRho, logT, logK, dlogK_dlogRho, dlogK_dlogT, ierr)
402 0 : return
403 : end if
404 :
405 0 : X0 = x_tables(ix)% X
406 0 : X1 = x_tables(ix+1)% X
407 :
408 0 : if (X1 <= X0) then
409 0 : ierr = 1
410 0 : return
411 : end if
412 :
413 0 : if (X0 >= X) then ! use the X0 table
414 : call Get_Kap_for_dXCO( &
415 : rq, iz, x_tables, dXC, dXO, ix, &
416 0 : logRho, logT, logK, dlogK_dlogRho, dlogK_dlogT, ierr)
417 0 : return
418 : end if
419 :
420 0 : if (X1 <= X) then ! use the X1 table
421 : call Get_Kap_for_dXCO( &
422 : rq, iz, x_tables, dXC, dXO, ix+1, &
423 0 : logRho, logT, logK, dlogK_dlogRho, dlogK_dlogT, ierr)
424 0 : return
425 : end if
426 :
427 0 : if (num_Xs >= 4 .and. rq% cubic_interpolation_in_X) then
428 : call Get_Kap_for_CO_X_cubic( &
429 : rq, iz, ix, dXC, dXO, x_tables, X, logRho, logT, &
430 0 : logK, dlogK_dlogRho, dlogK_dlogT, ierr)
431 : else ! linear
432 : call Get_Kap_for_CO_X_linear( &
433 : rq, iz, ix, dXC, dXO, x_tables, X, X0, X1, logRho, logT, &
434 0 : logK, dlogK_dlogRho, dlogK_dlogT, ierr)
435 : end if
436 :
437 : end subroutine Get_Kap_for_CO_X
438 :
439 :
440 : ! use tables ix-1 to ix+2 to do piecewise monotonic cubic interpolation in X
441 0 : subroutine Get_Kap_for_CO_X_cubic( &
442 : rq, iz, ix, dXC, dXO, x_tables, X, logRho, logT, &
443 : logK, dlogK_dlogRho, dlogK_dlogT, ierr)
444 : use kap_def
445 : use interp_1d_def, only: pm_work_size
446 : use interp_1d_lib, only: interpolate_vector_autodiff, interp_pm_autodiff
447 : use auto_diff
448 :
449 : ! return opacity from Z table number iz for the given X, XC, XO
450 : type (Kap_General_Info), pointer :: rq
451 : integer, intent(in) :: iz, ix
452 : type (Kap_CO_X_Table), dimension(:), pointer :: x_tables
453 : real(dp), intent(in) :: X, dXC, dXO
454 : real(dp), intent(inout) :: logRho, logT
455 : real(dp), intent(out) :: logK, dlogK_dlogRho, dlogK_dlogT
456 : integer, intent(out) :: ierr
457 :
458 : integer, parameter :: n_old = 4, n_new = 1
459 0 : real(dp), dimension(n_old) :: logKs, dlogKs_dlogRho, dlogKs_dlogT
460 : type(auto_diff_real_2var_order1), dimension(n_old) :: logKs_ad
461 : type(auto_diff_real_2var_order1) :: logK_ad
462 : type(auto_diff_real_2var_order1) :: x_old(n_old), x_new(n_new)
463 : type(auto_diff_real_2var_order1), target :: work_ary(n_old*pm_work_size)
464 : type(auto_diff_real_2var_order1), pointer :: work(:)
465 : integer :: i, i1, ixx, num_Xs
466 :
467 : logical, parameter :: dbg = .false.
468 :
469 : 11 format(a40,e20.10)
470 :
471 0 : ierr = 0
472 0 : work => work_ary
473 :
474 0 : num_Xs = num_kap_CO_Xs(rq% kap_CO_option)
475 :
476 0 : if (ix+2 > num_Xs) then
477 0 : i1 = num_Xs-2
478 0 : else if (ix == 1) then
479 : i1 = 2
480 : else
481 0 : i1 = ix
482 : end if
483 :
484 : if (dbg) write(*,*) 'ix', ix
485 :
486 0 : do i=1,n_old
487 0 : ixx = i1-2+i
488 : if (dbg) write(*,*) 'ixx', ixx
489 0 : x_old(i) % val = x_tables(ixx)% X
490 0 : x_old(i) % d1val1 = 0d0
491 0 : x_old(i) % d1val2 = 0d0
492 : call Get_Kap_for_dXCO( &
493 : rq, iz, x_tables, dXC, dXO, ixx, &
494 0 : logRho, logT, logKs(i), dlogKs_dlogRho(i), dlogKs_dlogT(i), ierr)
495 0 : if (ierr /= 0) then
496 : if (dbg) write(*,11) 'logRho', logRho
497 : if (dbg) write(*,11) 'logT', logT
498 : return
499 : end if
500 : ! now pack into auto_diff form
501 0 : logKs_ad(i) % val = logKs(i)
502 0 : logKs_ad(i) % d1val1 = dlogKs_dlogT(i)
503 0 : logKs_ad(i) % d1val2 = dlogKs_dlogRho(i)
504 : end do
505 0 : x_new(1) % val = X
506 0 : x_new(1) % d1val1 = 0d0
507 0 : x_new(1) % d1val2 = 0d0
508 :
509 0 : call interp1(logKs_ad, logK_ad, ierr)
510 0 : if (ierr /= 0) then
511 0 : call mesa_error(__FILE__,__LINE__,'failed in interp1 for logK')
512 0 : return
513 : end if
514 :
515 : ! unpack auto_diff pack into output reals
516 0 : logK = logK_ad % val
517 0 : dlogK_dlogT = logK_ad % d1val1
518 0 : dlogK_dlogRho = logK_ad % d1val2
519 :
520 0 : if (dbg) then
521 :
522 : do i=1,n_old
523 : write(*,*) 'x_old(i)', x_old(i)
524 : end do
525 : write(*,*) 'x_new(1)', x_new(1)
526 : write(*,'(A)')
527 :
528 : do i=1,n_old
529 : write(*,*) 'logKs(i)', logKs(i)
530 : end do
531 : write(*,*) 'logK', logK
532 : write(*,'(A)')
533 :
534 : do i=1,n_old
535 : write(*,*) 'dlogKs_dlogRho(i)', dlogKs_dlogRho(i)
536 : end do
537 : write(*,*) 'dlogK_dlogRho', dlogK_dlogRho
538 : write(*,'(A)')
539 :
540 : do i=1,n_old
541 : write(*,*) 'dlogKs_dlogT(i)', dlogKs_dlogT(i)
542 : end do
543 : write(*,*) 'dlogK_dlogT', dlogK_dlogT
544 : write(*,'(A)')
545 :
546 : end if
547 :
548 : contains
549 :
550 0 : subroutine interp1(old, new, ierr)
551 : type(auto_diff_real_2var_order1), intent(in) :: old(n_old)
552 : type(auto_diff_real_2var_order1), intent(out) :: new
553 : integer, intent(out) :: ierr
554 : type(auto_diff_real_2var_order1) :: v_old(n_old), v_new(n_new)
555 : integer :: i
556 0 : do i = 1, n_old
557 0 : v_old(i) = old(i)
558 : end do
559 : call interpolate_vector_autodiff( &
560 : n_old, x_old, n_new, x_new, v_old, v_new, interp_pm_autodiff, pm_work_size, work, &
561 0 : 'Get_Kap_for_CO_X_cubic', ierr)
562 0 : new = v_new(1)
563 0 : end subroutine interp1
564 :
565 : end subroutine Get_Kap_for_CO_X_cubic
566 :
567 :
568 0 : subroutine Get_Kap_for_CO_X_linear( &
569 : rq, iz, ix, dXC, dXO, x_tables, X, X0, X1, logRho, logT, &
570 : logK, dlogK_dlogRho, dlogK_dlogT, ierr)
571 : use kap_def
572 : ! return opacity from Z table number iz for the given X, dXC, dXO
573 : type (Kap_General_Info), pointer :: rq
574 : integer, intent(in) :: iz, ix
575 : type (Kap_CO_X_Table), dimension(:), pointer :: x_tables
576 : real(dp), intent(in) :: dXC, dXO, X, X0, X1
577 : real(dp), intent(inout) :: logRho, logT
578 : real(dp), intent(out) :: logK, dlogK_dlogRho, dlogK_dlogT
579 : integer, intent(out) :: ierr
580 :
581 0 : real(dp) :: logK0, dlogK0_dlogRho, dlogK0_dlogT, logK1, dlogK1_dlogRho, dlogK1_dlogT
582 0 : real(dp) :: alfa, beta
583 : integer :: i
584 :
585 : logical, parameter :: dbg = .false.
586 :
587 : ierr = 0
588 : call Get_Kap_for_dXCO( &
589 : rq, iz, x_tables, dXC, dXO, ix, &
590 0 : logRho, logT, logK0, dlogK0_dlogRho, dlogK0_dlogT, ierr)
591 0 : if (ierr /= 0) return
592 :
593 : call Get_Kap_for_dXCO( &
594 : rq, iz, x_tables, dXC, dXO, ix+1, &
595 0 : logRho, logT, logK1, dlogK1_dlogRho, dlogK1_dlogT, ierr)
596 0 : if (ierr /= 0) return
597 :
598 : ! X0 result in logK0, X1 result in logK1
599 0 : beta = (X - X1) / (X0 - X1) ! beta -> 1 as X -> X0
600 0 : alfa = 1d0 - beta
601 :
602 0 : logK = beta*logK0 + alfa*logK1
603 0 : dlogK_dlogRho = beta*dlogK0_dlogRho + alfa*dlogK1_dlogRho
604 0 : dlogK_dlogT = beta*dlogK0_dlogT + alfa*dlogK1_dlogT
605 :
606 : end subroutine Get_Kap_for_CO_X_linear
607 :
608 :
609 4 : subroutine Get_Kap_for_dXCO( &
610 : rq, iz, x_tables, dXC_in, dXO_in, ix, &
611 : logRho, logT, logK, dlogK_dlogRho, dlogK_dlogT, ierr)
612 : ! return value from xtable number ix for the given dXC and dXO
613 : use kap_def
614 : use load_CO_kap
615 : type (Kap_General_Info), pointer :: rq
616 : type (Kap_CO_X_Table), dimension(:), pointer :: x_tables
617 : integer :: iz, ix
618 : real(dp), intent(in) :: dXC_in, dXO_in
619 : real(dp), intent(inout) :: logRho, logT
620 : real(dp), intent(out) :: logK, dlogK_dlogRho, dlogK_dlogT
621 : integer, intent(out) :: ierr
622 :
623 4 : type (Kap_CO_Table), dimension(:), pointer :: co_tables ! stored by table number
624 4 : real(dp) :: dXC, dXO, fac, dXCO_max, Z, dXC_lookup, dXO_lookup
625 : integer :: num_CO_tables, num_dXC_gt_dXO, i1, i2, i3, i4
626 4 : real(dp) :: alfa, beta
627 4 : real(dp) :: logK1, dlogK1_dlogRho, dlogK1_dlogT, logK2, dlogK2_dlogRho, dlogK2_dlogT
628 4 : real(dp) :: logK3, dlogK3_dlogRho, dlogK3_dlogT, logK4, dlogK4_dlogRho, dlogK4_dlogT
629 4 : real(dp) :: dXC1_lookup, dXO1_lookup, dXC2_lookup, dXO2_lookup, dXC_2_4_lookup, dXO_2_4_lookup
630 4 : real(dp) :: dXC3_lookup, dXO3_lookup, dXC4_lookup, dXO4_lookup, dXC_1_3_lookup, dXO_1_3_lookup
631 4 : real(dp) :: logK_2_4, dlogK_2_4_dlogRho, dlogK_2_4_dlogT
632 4 : real(dp) :: logK_1_3, dlogK_1_3_dlogRho, dlogK_1_3_dlogT
633 : logical, parameter :: read_later = .false., dbg = .false.
634 :
635 : include 'formats'
636 :
637 4 : ierr = 0
638 :
639 : if (dbg) write(*,1) 'enter Get_Kap_for_dXCO dXC_in dXO_in', dXC_in, dXO_in
640 :
641 4 : dXC = max(0.0_dp, dXC_in)
642 4 : dXO = max(0.0_dp, dXO_in)
643 :
644 4 : if (x_tables(ix)% not_loaded_yet) then ! avoid doing critical section if possible
645 4 : !$omp critical (load_co_table)
646 4 : if (x_tables(ix)% not_loaded_yet) then
647 4 : call load_one_CO(rq,kap_co_z_tables(rq% kap_CO_option)% ar,iz,ix,read_later,ierr)
648 : end if
649 : !$omp end critical (load_co_table)
650 : end if
651 4 : if (ierr /= 0) return
652 :
653 4 : co_tables => x_tables(ix)% co_tables
654 4 : num_CO_tables = x_tables(ix)% num_CO_tables
655 : if (dbg) write(*,2) 'num_CO_tables', num_CO_tables
656 :
657 4 : if (num_CO_tables < 1) then
658 0 : ierr = -1
659 0 : write(*,2) 'num_CO_tables', num_CO_tables
660 0 : return
661 : end if
662 :
663 4 : num_dXC_gt_dXO = x_tables(ix)% num_dXC_gt_dXO
664 4 : Z = x_tables(ix)% Z
665 :
666 4 : dXCO_max = 1 - ((x_tables(ix)% X) + (x_tables(ix)% Z))
667 4 : if (dXC + dXO > dXCO_max) then
668 0 : fac = dXCO_max / (dXC + dXO)
669 0 : dXC = fac*dXC
670 0 : dXO = fac*dXO
671 : end if
672 :
673 4 : dXC_lookup = get_dX_lookup(dXC, Z)
674 4 : dXO_lookup = get_dX_lookup(dXO, Z)
675 :
676 : if (dbg) write(*,2) 'call Find_CO_Tables', ix
677 : call Find_CO_Tables(rq, x_tables, ix, x_tables(ix)% CO_table_numbers, &
678 : x_tables(ix)% next_dXO_table, x_tables(ix)% next_dXC_table, &
679 : co_tables, num_CO_tables, num_dXC_gt_dXO, &
680 4 : dXCO_max, dXC, dXO, dXC_lookup, dXO_lookup, i1, i2, i3, i4, ierr)
681 4 : if (ierr /= 0) then
682 0 : write(*,*) 'kap failed in Find_CO_Tables'
683 0 : return
684 : end if
685 :
686 4 : if (i1 > 0 .and. i2 <= 0 .and. i3 <= 0 .and. i4 <= 0) then
687 : call Get_CO_Kap_for_logRho_logT(rq, x_tables, ix, co_tables, i1, logRho, logT, &
688 0 : logK, dlogK_dlogRho, dlogK_dlogT, ierr)
689 0 : return
690 : end if
691 :
692 4 : if (i1 <= 0 .or. i2 <= 0 .or. i3 <= 0) call mesa_error(__FILE__,__LINE__,'error in result from Find_CO_Tables')
693 :
694 4 : if (matches_table(i2)) then
695 : call Get_CO_Kap_for_logRho_logT(rq, x_tables, ix, co_tables, i2, logRho, logT, &
696 0 : logK, dlogK_dlogRho, dlogK_dlogT, ierr)
697 0 : return
698 : end if
699 :
700 4 : if (matches_table(i3)) then
701 : call Get_CO_Kap_for_logRho_logT(rq, x_tables, ix, co_tables, i3, logRho, logT, &
702 0 : logK, dlogK_dlogRho, dlogK_dlogT, ierr)
703 0 : return
704 : end if
705 :
706 4 : if (i4 > 0) then
707 4 : if (matches_table(i4)) then
708 : call Get_CO_Kap_for_logRho_logT(rq, x_tables, ix, co_tables, i4, logRho, logT, &
709 0 : logK, dlogK_dlogRho, dlogK_dlogT, ierr)
710 0 : return
711 : end if
712 : end if
713 :
714 : call Get_CO_Kap_for_logRho_logT(rq, x_tables, ix, co_tables, i1, logRho, logT, &
715 4 : logK1, dlogK1_dlogRho, dlogK1_dlogT, ierr)
716 4 : if (ierr /= 0) return
717 4 : dXC1_lookup = co_tables(i1)% dXC_lookup
718 4 : dXO1_lookup = co_tables(i1)% dXO_lookup
719 :
720 : call Get_CO_Kap_for_logRho_logT(rq, x_tables, ix, co_tables, i2, logRho, logT, &
721 4 : logK2, dlogK2_dlogRho, dlogK2_dlogT, ierr)
722 4 : if (ierr /= 0) return
723 4 : dXC2_lookup = co_tables(i2)% dXC_lookup
724 4 : dXO2_lookup = co_tables(i2)% dXO_lookup
725 :
726 : call Get_CO_Kap_for_logRho_logT(rq, x_tables, ix, co_tables, i3, logRho, logT, &
727 4 : logK3, dlogK3_dlogRho, dlogK3_dlogT, ierr)
728 4 : if (ierr /= 0) return
729 4 : dXC3_lookup = co_tables(i3)% dXC_lookup
730 4 : dXO3_lookup = co_tables(i3)% dXO_lookup
731 :
732 4 : if (i4 > 0) then
733 : call Get_CO_Kap_for_logRho_logT(rq, x_tables, ix, co_tables, i4, logRho, logT, &
734 4 : logK4, dlogK4_dlogRho, dlogK4_dlogT, ierr)
735 4 : if (ierr /= 0) return
736 4 : dXC4_lookup = co_tables(i4)% dXC_lookup
737 4 : dXO4_lookup = co_tables(i4)% dXO_lookup
738 : else ! copy i3 results
739 0 : logK4 = logK3
740 0 : dlogK4_dlogRho = dlogK3_dlogRho
741 0 : dlogK4_dlogT = dlogK3_dlogT
742 0 : dXC4_lookup = dXC3_lookup
743 0 : dXO4_lookup = dXO3_lookup
744 : end if
745 :
746 4 : if (dXC >= dXO) then ! use values on lines i1-i3 and i2-i4 at dXO
747 :
748 : call Get_Kap_at_dXO(dXO_lookup, &
749 : dXC2_lookup, dXO2_lookup, logK2, dlogK2_dlogRho, dlogK2_dlogT, &
750 : dXC4_lookup, dXO4_lookup, logK4, dlogK4_dlogRho, dlogK4_dlogT, &
751 4 : logK_2_4, dlogK_2_4_dlogRho, dlogK_2_4_dlogT, dXC_2_4_lookup)
752 :
753 : call Get_Kap_at_dXO(dXO_lookup, &
754 : dXC1_lookup, dXO1_lookup, logK1, dlogK1_dlogRho, dlogK1_dlogT, &
755 : dXC3_lookup, dXO3_lookup, logK3, dlogK3_dlogRho, dlogK3_dlogT, &
756 4 : logK_1_3, dlogK_1_3_dlogRho, dlogK_1_3_dlogT, dXC_1_3_lookup)
757 4 : if (dXC_1_3_lookup == dXC_2_4_lookup) then
758 : alfa = 0d0
759 : else
760 4 : alfa = (dXC_lookup - dXC_2_4_lookup) / (dXC_1_3_lookup - dXC_2_4_lookup)
761 : end if
762 :
763 : else ! use values on lines i1-i3 and i2-i4 at dXC
764 :
765 : call Get_Kap_at_dXC(dXC_lookup, &
766 : dXC2_lookup, dXO2_lookup, logK2, dlogK2_dlogRho, dlogK2_dlogT, &
767 : dXC4_lookup, dXO4_lookup, logK4, dlogK4_dlogRho, dlogK4_dlogT, &
768 0 : logK_2_4, dlogK_2_4_dlogRho, dlogK_2_4_dlogT, dXO_2_4_lookup)
769 :
770 : call Get_Kap_at_dXC(dXC_lookup, &
771 : dXC1_lookup, dXO1_lookup, logK1, dlogK1_dlogRho, dlogK1_dlogT, &
772 : dXC3_lookup, dXO3_lookup, logK3, dlogK3_dlogRho, dlogK3_dlogT, &
773 0 : logK_1_3, dlogK_1_3_dlogRho, dlogK_1_3_dlogT, dXO_1_3_lookup)
774 0 : if (dXO_1_3_lookup == dXO_2_4_lookup) then
775 : alfa = 0d0
776 : else
777 0 : alfa = (dXO_lookup - dXO_2_4_lookup) / (dXO_1_3_lookup - dXO_2_4_lookup)
778 : end if
779 :
780 : end if
781 :
782 4 : beta = 1d0 - alfa
783 :
784 4 : logK = alfa*logK_1_3 + beta*logK_2_4
785 4 : dlogK_dlogRho = alfa*dlogK_1_3_dlogRho + beta*dlogK_2_4_dlogRho
786 4 : dlogK_dlogT = alfa*dlogK_1_3_dlogT + beta*dlogK_2_4_dlogT
787 :
788 4 : if (is_bad(logK)) then
789 0 : ierr = -1
790 0 : return
791 : write(*,1) 'logK', logK
792 : write(*,1) 'logK_1_3', logK_1_3
793 : write(*,1) 'logK_2_4', logK_2_4
794 : write(*,1) 'alfa', alfa
795 : write(*,1) 'beta', beta
796 : write(*,'(A)')
797 : write(*,2) 'dXC1_lookup', i1, dXC1_lookup
798 : write(*,2) 'dXO1_lookup', i1, dXO1_lookup
799 : write(*,2) 'dXC2_lookup', i2, dXC2_lookup
800 : write(*,2) 'dXO2_lookup', i2, dXO2_lookup
801 : write(*,2) 'dXC3_lookup', i3, dXC3_lookup
802 : write(*,2) 'dXO3_lookup', i3, dXO3_lookup
803 : write(*,2) 'dXC4_lookup', i4, dXC4_lookup
804 : write(*,2) 'dXO4_lookup', i4, dXO4_lookup
805 : write(*,1) 'dXC', dXC
806 : write(*,1) 'dXO', dXO
807 : call mesa_error(__FILE__,__LINE__,'Get_Kap_for_dXCO')
808 : end if
809 :
810 : contains
811 :
812 :
813 12 : logical function matches_table(i)
814 : integer :: i
815 12 : if (i < 1 .or. i > num_CO_tables) then
816 0 : write(*,*) 'logRho', logRho
817 0 : write(*,*) 'logT', logT
818 0 : call mesa_error(__FILE__,__LINE__,'bug in kap_eval_co matches_table')
819 0 : matches_table = .false.
820 12 : else if (abs(dXC_lookup - co_tables(i)% dXC_lookup) == 0 .and. &
821 : abs(dXO_lookup - co_tables(i)% dXO_lookup) == 0) then
822 : matches_table = .true.
823 : else
824 12 : matches_table = .false.
825 : end if
826 4 : end function matches_table
827 :
828 :
829 8 : subroutine Get_Kap_at_dXO(dXO_lookup, &
830 : dXC_a_lookup, dXO_a_lookup, logK_a, dlogK_a_dlogRho, dlogK_a_dlogT, &
831 : dXC_b_lookup, dXO_b_lookup, logK_b, dlogK_b_dlogRho, dlogK_b_dlogT, &
832 : logK_a_b, dlogK_a_b_dlogRho, dlogK_a_b_dlogT, dXC_a_b_lookup)
833 : real(dp), intent(in) :: dXO_lookup
834 : real(dp), intent(in) :: dXC_a_lookup, dXO_a_lookup
835 : real(dp), intent(in) :: logK_a, dlogK_a_dlogRho, dlogK_a_dlogT
836 : real(dp), intent(in) :: dXC_b_lookup, dXO_b_lookup
837 : real(dp), intent(in) :: logK_b, dlogK_b_dlogRho, dlogK_b_dlogT
838 : real(dp), intent(out) :: logK_a_b, dlogK_a_b_dlogRho, dlogK_a_b_dlogT
839 : real(dp), intent(out) :: dXC_a_b_lookup
840 :
841 8 : real(dp) :: alfa, beta
842 :
843 8 : if (dXO_a_lookup == dXO_b_lookup) then
844 : alfa = 0d0
845 : else
846 8 : alfa = (dXO_lookup - dXO_b_lookup) / (dXO_a_lookup - dXO_b_lookup)
847 : end if
848 :
849 8 : dXC_a_b_lookup = dXC_b_lookup + (dXC_a_lookup - dXC_b_lookup)*alfa
850 8 : beta = 1d0 - alfa
851 8 : logK_a_b = alfa*logK_a + beta*logK_b
852 8 : dlogK_a_b_dlogRho = alfa*dlogK_a_dlogRho + beta*dlogK_b_dlogRho
853 8 : dlogK_a_b_dlogT = alfa*dlogK_a_dlogT + beta*dlogK_b_dlogT
854 :
855 8 : end subroutine Get_Kap_at_dXO
856 :
857 :
858 0 : subroutine Get_Kap_at_dXC(dXC_lookup, &
859 : dXC_a_lookup, dXO_a_lookup, logK_a, dlogK_a_dlogRho, dlogK_a_dlogT, &
860 : dXC_b_lookup, dXO_b_lookup, logK_b, dlogK_b_dlogRho, dlogK_b_dlogT, &
861 : logK_a_b, dlogK_a_b_dlogRho, dlogK_a_b_dlogT, dXO_a_b_lookup)
862 : real(dp), intent(in) :: dXC_lookup
863 : real(dp), intent(in) :: dXC_a_lookup, dXO_a_lookup
864 : real(dp), intent(in) :: logK_a, dlogK_a_dlogRho, dlogK_a_dlogT
865 : real(dp), intent(in) :: dXC_b_lookup, dXO_b_lookup
866 : real(dp), intent(in) :: logK_b, dlogK_b_dlogRho, dlogK_b_dlogT
867 : real(dp), intent(out) :: logK_a_b, dlogK_a_b_dlogRho, dlogK_a_b_dlogT
868 : real(dp), intent(out) :: dXO_a_b_lookup
869 :
870 0 : real(dp) :: alfa, beta
871 :
872 0 : if (dXC_a_lookup == dXC_b_lookup) then
873 : alfa = 0d0
874 : else
875 0 : alfa = (dXC_lookup - dXC_b_lookup) / (dXC_a_lookup - dXC_b_lookup)
876 : end if
877 :
878 0 : dXO_a_b_lookup = dXO_b_lookup + (dXO_a_lookup - dXO_b_lookup)*alfa
879 0 : beta = 1d0 - alfa
880 :
881 0 : logK_a_b = alfa*logK_a + beta*logK_b
882 0 : dlogK_a_b_dlogRho = alfa*dlogK_a_dlogRho + beta*dlogK_b_dlogRho
883 0 : dlogK_a_b_dlogT = alfa*dlogK_a_dlogT + beta*dlogK_b_dlogT
884 :
885 0 : end subroutine Get_Kap_at_dXC
886 :
887 :
888 : end subroutine Get_Kap_for_dXCO
889 :
890 :
891 4 : subroutine Find_CO_Tables( &
892 : rq, x_tables, ix, CO_table_numbers, next_dXO_table, next_dXC_table, &
893 : co_tables, num_CO_tables, num_dXC_gt_dXO, &
894 : dXCO_max, dXC, dXO, dXC_lookup, dXO_lookup, i1, i2, i3, i4, ierr)
895 :
896 : ! for linear interpolation to be smooth,
897 : ! must use the smallest convex hull around the given point
898 : use kap_def
899 : use load_CO_kap
900 :
901 : type (Kap_General_Info), pointer :: rq
902 : type (Kap_CO_X_Table), dimension(:), pointer :: x_tables
903 : integer :: ix
904 : integer, intent(in) :: CO_table_numbers(num_kap_CO_dXs,num_kap_CO_dXs)
905 : integer, intent(in) :: next_dXO_table(max_num_CO_tables)
906 : integer, intent(in) :: next_dXC_table(max_num_CO_tables)
907 : type (Kap_CO_Table), dimension(:), pointer :: co_tables
908 : integer, intent(in) :: num_CO_tables, num_dXC_gt_dXO
909 : real(dp), intent(in) :: dXCO_max, dXC, dXO, dXC_lookup, dXO_lookup
910 : integer, intent(out) :: i1, i2, i3, i4
911 : integer, intent(out) :: ierr
912 :
913 4 : real(dp) :: dXC2_lookup, dXO2_lookup, dXC4_lookup, dXO4_lookup
914 : integer :: idXC, idXO
915 : real(dp), parameter :: tiny = 1d-7
916 :
917 : logical, parameter :: dbg = .false.
918 :
919 : include 'formats'
920 :
921 : if (dbg) write(*,*) 'enter Find_CO_Tables'
922 : if (dbg) write(*,*) 'associated(co_tables)', associated(co_tables)
923 : if (dbg) write(*,*) 'size(co_tables,dim=1)', size(co_tables,dim=1)
924 : if (dbg) write(*,2) 'num_kap_CO_dXs', num_kap_CO_dXs
925 : if (dbg) write(*,2) 'num_CO_tables', num_CO_tables
926 :
927 4 : ierr = 0
928 :
929 : ! find idXC s.t. kap_CO_dXs(idXC-1) < dXC <= kap_CO_dXs(idXC)
930 8 : do idXC = 2, num_kap_CO_dXs
931 8 : if (kap_CO_dXs(idXC) >= dXC) exit
932 : end do
933 : if (dbg) write(*,2) 'idXC', idXC
934 :
935 : ! find idXO s.t. kap_CO_dXs(idXO-1) < dXO <= kap_CO_dXs(idXO)
936 8 : do idXO = 2, num_kap_CO_dXs
937 8 : if (kap_CO_dXs(idXO) >= dXO) exit
938 : end do
939 : if (dbg) write(*,2) 'idXO', idXO
940 :
941 4 : i1 = CO_table_numbers(idXC-1,idXO-1)
942 : if (dbg) write(*,2) 'i1', i1
943 4 : if (matches_table(i1)) then
944 0 : i2 = -1; i3 = -1; i4 = -1
945 : if (dbg) write(*,2) 'matches_table(i1)', i1
946 0 : return
947 : end if
948 :
949 : if (dbg) write(*,*) '(dXC >= dXO)', (dXC >= dXO)
950 4 : if (dXC >= dXO) then
951 4 : i2 = CO_table_numbers(idXC,idXO-1)
952 : if (dbg) write(*,2) 'i2', i2
953 4 : i3 = CO_table_numbers(idXC-1,idXO)
954 : if (dbg) write(*,2) 'i3', i3
955 : else
956 0 : i2 = CO_table_numbers(idXC-1,idXO)
957 : if (dbg) write(*,2) 'i2', i2
958 0 : i3 = CO_table_numbers(idXC,idXO-1)
959 : if (dbg) write(*,2) 'i3', i3
960 : end if
961 4 : i4 = CO_table_numbers(idXC,idXO)
962 : if (dbg) write(*,2) 'i4', i4
963 :
964 4 : if (i4 > 0) then
965 4 : if (i1 <= 0 .or. i2 <= 0 .or. i3 <= 0) then
966 :
967 0 : write(*,2) 'i1', i1
968 0 : write(*,2) 'i2', i2
969 0 : write(*,2) 'i3', i3
970 0 : write(*,2) 'i4', i4
971 0 : write(*,2) 'idXC', idXC
972 0 : write(*,2) 'idXO', idXO
973 :
974 0 : write(*,1) 'dXCO_max', dble(dXCO_max)
975 0 : write(*,1) 'dXC', dble(dXC)
976 0 : write(*,1) 'dXO', dble(dXO)
977 0 : write(*,1) 'dXC_lookup', dble(dXC_lookup)
978 0 : write(*,1) 'dXO_lookup', dble(dXO_lookup)
979 :
980 0 : call mesa_error(__FILE__,__LINE__,'logical failure1 in looking for CO tables')
981 : end if
982 4 : if (matches_table(i2)) then
983 0 : i1 = i2; i2 = -1; i3 = -1; i4 = -1; return
984 : end if
985 4 : if (matches_table(i3)) then
986 0 : i1 = i3; i2 = -1; i3 = -1; i4 = -1; return
987 : end if
988 4 : if (matches_table(i4)) then
989 0 : i1 = i4; i2 = -1; i3 = -1; i4 = -1; return
990 : end if
991 : return
992 : end if
993 :
994 0 : if (on_midline(i1)) then ! middle triangle
995 0 : if (dXC >= dXO) then
996 0 : i2 = next_dXC_table(i1)
997 0 : i3 = next_dXO_table(i1)
998 : else
999 0 : i2 = next_dXO_table(i1)
1000 0 : i3 = next_dXC_table(i1)
1001 : end if
1002 0 : return
1003 : end if
1004 :
1005 : ! trapezoid or triangle near the diagonal boundary
1006 :
1007 4 : if (dXC >= dXO) then
1008 :
1009 0 : if (i3 <= 0) then
1010 0 : if (on_diagonal(i1)) then
1011 0 : i3 = num_CO_tables
1012 0 : if (i2 > 0) then ! bail -- just use i1
1013 0 : i2 = -1; i3 = -1; i4 = -1; return
1014 : end if
1015 : else
1016 0 : i3 = num_dXC_gt_dXO
1017 : end if
1018 : end if
1019 :
1020 0 : if (.not. on_diagonal(i3)) then
1021 0 : i4 = next_dXC_table(i3)
1022 0 : if (.not. on_diagonal(i4)) then ! bail -- just use i1
1023 0 : i2 = -1; i3 = -1; i4 = -1; return
1024 : end if
1025 : end if
1026 0 : if (i2 <= 0) i2 = next_dXC_table(i1)
1027 0 : if (on_diagonal(i2)) return
1028 :
1029 :
1030 0 : dXC4_lookup = co_tables(i4)% dXC_lookup
1031 0 : if (dXC_lookup <= dXC4_lookup) return
1032 : ! check if on smaller dXC_lookup side of the line from i2 to i4
1033 0 : dXC2_lookup = co_tables(i2)% dXC_lookup
1034 0 : dXO2_lookup = co_tables(i2)% dXO_lookup
1035 0 : dXO4_lookup = co_tables(i4)% dXO_lookup
1036 0 : if (dXC_lookup <= dXC4_lookup+(dXC2_lookup-dXC4_lookup)* &
1037 : (dXO_lookup-dXO4_lookup)/(dXO2_lookup-dXO4_lookup)) return
1038 : ! else we're in the dXC > dXO triangle
1039 0 : i1 = i2; i2 = next_dXC_table(i1)
1040 0 : if (.not. on_diagonal(i2)) then ! bail -- just use i1
1041 0 : i2 = -1; i3 = -1; i4 = -1; return
1042 : end if
1043 0 : i3 = i4; i4 = -1
1044 :
1045 : else ! dXC < dXO
1046 :
1047 : ! reverse roles of dXC and dXO
1048 :
1049 0 : if (i3 <= 0) then ! must be in one of the triangles
1050 0 : if (on_diagonal(i1)) then
1051 0 : i3 = num_dXC_gt_dXO
1052 : else
1053 0 : i3 = num_CO_tables
1054 0 : if (i2 > 0) then ! bail -- just use i1
1055 0 : i2 = -1; i3 = -1; i4 = -1; return
1056 : end if
1057 : end if
1058 : end if
1059 0 : if (.not. on_diagonal(i3)) then
1060 0 : i4 = next_dXO_table(i3)
1061 0 : if (.not. on_diagonal(i4)) then ! bail -- just use i1
1062 0 : i2 = -1; i3 = -1; i4 = -1; return
1063 : end if
1064 : end if
1065 0 : if (i2 <= 0) i2 = next_dXO_table(i1)
1066 0 : if (on_diagonal(i2)) return
1067 0 : dXO4_lookup = co_tables(i4)% dXO_lookup
1068 0 : if (dXO_lookup <= dXO4_lookup) return
1069 : ! check if on smaller dXO_lookup side of the line from i2 to i4
1070 0 : dXC2_lookup = co_tables(i2)% dXC_lookup
1071 0 : dXO2_lookup = co_tables(i2)% dXO_lookup
1072 0 : dXC4_lookup = co_tables(i4)% dXC_lookup
1073 0 : if (dXO_lookup <= dXO4_lookup+(dXO2_lookup-dXO4_lookup)* &
1074 : (dXC_lookup-dXC4_lookup)/(dXC2_lookup-dXC4_lookup)) return
1075 : ! else we're in the dXC < dXO triangle
1076 0 : i1 = i2; i2 = next_dXO_table(i1)
1077 0 : if (.not. on_diagonal(i2)) then ! bail -- just use i1
1078 0 : i2 = -1; i3 = -1; i4 = -1; return
1079 : end if
1080 0 : i3 = num_CO_tables; i4 = -1
1081 :
1082 : end if
1083 :
1084 :
1085 : contains
1086 :
1087 :
1088 16 : logical function matches_table(i)
1089 : integer :: i
1090 : include 'formats'
1091 : if (dbg) write(*,2) 'matches_table i', i
1092 16 : if (abs(dXC_lookup - co_tables(i)% dXC_lookup) < tiny .and. &
1093 : abs(dXO_lookup - co_tables(i)% dXO_lookup) < tiny) then
1094 : matches_table = .true.
1095 : else
1096 16 : matches_table = .false.
1097 : end if
1098 : if (dbg) write(*,*) 'matches_table', matches_table
1099 4 : end function matches_table
1100 :
1101 :
1102 0 : logical function on_midline(i)
1103 : integer :: i
1104 0 : if (abs(co_tables(i)% dXC - co_tables(i)% dXO) < tiny) then
1105 : on_midline = .true.
1106 : else
1107 0 : on_midline = .false.
1108 : end if
1109 0 : end function on_midline
1110 :
1111 :
1112 0 : logical function on_diagonal(i)
1113 : integer :: i
1114 0 : if (abs((co_tables(i)% dXC + co_tables(i)% dXO) - dXCO_max) < tiny) then
1115 : on_diagonal = .true.
1116 : else
1117 0 : on_diagonal = .false.
1118 : end if
1119 0 : end function on_diagonal
1120 :
1121 :
1122 : end subroutine Find_CO_Tables
1123 :
1124 :
1125 64 : subroutine Get_CO_Kap_for_logRho_logT( &
1126 : rq, x_tables, ix, co_tables, ico, &
1127 : logRho, logT, logK, dlogK_dlogRho, dlogK_dlogT, ierr)
1128 : use kap_def
1129 : type (Kap_General_Info), pointer :: rq
1130 : type (Kap_CO_X_Table), dimension(:), pointer :: x_tables
1131 : integer :: ix
1132 : type (Kap_CO_Table), dimension(:), pointer :: co_tables
1133 : integer, intent(in) :: ico
1134 : real(dp), intent(inout) :: logRho, logT
1135 : real(dp), intent(out) :: logK, dlogK_dlogRho, dlogK_dlogT
1136 : integer, intent(out) :: ierr
1137 :
1138 16 : real(dp) :: logR0, logR1, logT0, logT1, logR, logR_in
1139 16 : real(dp) :: df_dx, df_dy
1140 : integer :: iR, jtemp, i, num_logRs, num_logTs
1141 : logical :: clipped_logR
1142 :
1143 : logical, parameter :: dbg = .false.
1144 :
1145 : include 'formats'
1146 :
1147 : if (dbg) write(*,1) 'enter Get_CO_Kap_for_logRho_logT', logRho, logT
1148 :
1149 16 : ierr = 0
1150 :
1151 : ! logR from inputs
1152 16 : logR_in = logRho - 3d0*logT + 18d0
1153 :
1154 : ! blends at higher levels MUST prevent
1155 : ! these tables from being called off their
1156 : ! high/low T and low R edges
1157 :
1158 16 : if (logT > x_tables(ix)% logT_max) then
1159 0 : ierr = -1
1160 0 : return
1161 : end if
1162 :
1163 16 : if (logT < x_tables(ix)% logT_min) then
1164 0 : ierr = -1
1165 0 : return
1166 : end if
1167 :
1168 16 : if (logR_in < x_tables(ix)% logR_min) then
1169 0 : ierr = -1
1170 0 : return
1171 : end if
1172 :
1173 :
1174 : ! off the high R edge, we use the input temperature
1175 : ! but clip logR to the table edge value
1176 :
1177 16 : if (logR_in > x_tables(ix)% logR_max) then
1178 0 : logR = x_tables(ix)% logR_max
1179 0 : clipped_logR = .true.
1180 : else
1181 16 : logR = logR_in
1182 16 : clipped_logR = .false.
1183 : end if
1184 :
1185 :
1186 16 : num_logRs = x_tables(ix)% num_logRs
1187 16 : num_logTs = x_tables(ix)% num_logTs
1188 :
1189 16 : if (num_logRs <= 0) then
1190 0 : write(*,*) 'num_logRs', num_logRs
1191 0 : write(*,*) 'ix', ix
1192 0 : call mesa_error(__FILE__,__LINE__,'Get_Kap_for_logRho_logT')
1193 : end if
1194 :
1195 16 : if (num_logTs <= 0) then
1196 0 : write(*,*) 'num_logTs', num_logRs
1197 0 : write(*,*) 'ix', ix
1198 0 : call mesa_error(__FILE__,__LINE__,'Get_Kap_for_logRho_logT')
1199 : end if
1200 :
1201 : call Locate_logR( &
1202 : rq, num_logRs, x_tables(ix)% logR_min, x_tables(ix)% logR_max, &
1203 16 : x_tables(ix)% ili_logRs, x_tables(ix)% logRs, logR, iR, logR0, logR1, ierr)
1204 16 : if (ierr /= 0) then
1205 0 : write(*,1) 'x_tables(ix)% logR_min', x_tables(ix)% logR_min
1206 0 : write(*,1) 'x_tables(ix)% logR_max', x_tables(ix)% logR_max
1207 0 : write(*,2) 'num_logRs', num_logRs
1208 0 : write(*,2) 'iR', iR
1209 0 : write(*,1) 'logR', logR
1210 0 : write(*,1) 'logR0', logR0
1211 0 : write(*,1) 'logR1', logR1
1212 0 : do i=1,num_logRs
1213 0 : write(*,2) 'logR', i, x_tables(ix)% logRs(i)
1214 : end do
1215 0 : write(*,*) 'clip_to_kap_table_boundaries', clip_to_kap_table_boundaries
1216 0 : call mesa_error(__FILE__,__LINE__,'failed in Locate_logR')
1217 0 : return
1218 : end if
1219 :
1220 : call Locate_logT( &
1221 : rq, num_logTs, x_tables(ix)% logT_min, x_tables(ix)% logT_max, &
1222 16 : x_tables(ix)% ili_logTs, x_tables(ix)% logTs, logT, jtemp, logT0, logT1, ierr)
1223 16 : if (ierr /= 0) return
1224 :
1225 : call Do_Kap_Interpolations( &
1226 : co_tables(ico)% kap1, num_logRs, num_logTs, &
1227 : iR, jtemp, logR0, logR, logR1, logT0, logT, logT1, &
1228 16 : logK, df_dx, df_dy)
1229 16 : if (clipped_logR) df_dx = 0
1230 :
1231 : if (dbg) write(*,1) 'Do_Kap_Interpolations: logK', logK
1232 :
1233 : ! convert df_dx and df_dy to dlogK_dlogRho, dlogK_dlogT
1234 16 : dlogK_dlogRho = df_dx
1235 16 : dlogK_dlogT = df_dy - 3d0*df_dx
1236 :
1237 : if (dbg) write(*,*) 'done Get_CO_Kap_for_logRho_logT'
1238 :
1239 : end subroutine Get_CO_Kap_for_logRho_logT
1240 :
1241 : end module kap_eval_co
|
