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 eosDT_eval
21 :
22 : use eos_def
23 : use const_def, only: avo, crad, ln10, arg_not_provided, mp, kerg, dp, qp
24 : use utils_lib, only: is_bad, mesa_error
25 : use math_lib
26 : use eosDT_load_tables, only: load_single_eosDT_table_by_id
27 : use eos_HELM_eval
28 : use eoscms_eval, only: Get_CMS_alfa, get_CMS_for_eosdt
29 : use skye, only: get_Skye_for_eosdt, get_Skye_alfa, get_Skye_alfa_simple
30 : use ideal, only: get_ideal_for_eosdt
31 :
32 : implicit none
33 :
34 : logical, parameter :: return_ierr_beyond_table_bounds = .true.
35 :
36 : integer, parameter :: use_none = 1
37 : integer, parameter :: use_all = 2
38 : integer, parameter :: blend_in_x = 3
39 : integer, parameter :: blend_in_y = 4
40 : integer, parameter :: blend_corner_out = 5
41 : integer, parameter :: blend_corner_in = 6
42 : integer, parameter :: blend_diagonal = 7
43 : integer, parameter :: blend_in_Z = 8
44 :
45 : abstract interface
46 : subroutine get_values_for_eosdt_interface( &
47 : handle, dbg, Z, X, abar, zbar, &
48 : species, chem_id, net_iso, xa, &
49 : rho, logRho, T, logT, remaining_fraction, &
50 : res, d_dlnd, d_dlnT, d_dxa, &
51 : skip, ierr)
52 : use const_def, only: dp
53 : use eos_def, only: nv
54 : implicit none
55 : integer, intent(in) :: handle
56 : logical, intent(in) :: dbg
57 : real(dp), intent(in) :: &
58 : Z, X, abar, zbar, remaining_fraction
59 : integer, intent(in) :: species
60 : integer, pointer :: chem_id(:), net_iso(:)
61 : real(dp), intent(in) :: xa(:)
62 : real(dp), intent(in) :: rho, logRho, T, logT
63 : real(dp), intent(inout), dimension(nv) :: &
64 : res, d_dlnd, d_dlnT
65 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
66 : logical, intent(out) :: skip
67 : integer, intent(out) :: ierr
68 : end subroutine get_values_for_eosdt_interface
69 : end interface
70 :
71 :
72 : contains
73 :
74 :
75 0 : subroutine Test_one_eosDT_component(rq, which_eos, &
76 : Z, X, abar, zbar, &
77 0 : species, chem_id, net_iso, xa, &
78 : arho, alogrho, atemp, alogtemp, &
79 0 : res, d_dlnd, d_dlnT, d_dxa, ierr)
80 : type (EoS_General_Info), pointer :: rq
81 : real(dp), intent(in) :: Z, X, abar, zbar
82 : integer, intent(in) :: which_eos, species
83 : integer, pointer :: chem_id(:), net_iso(:)
84 : real(dp), intent(in) :: xa(:)
85 : real(dp), intent(in) :: arho, alogrho, atemp, alogtemp
86 : real(dp), intent(inout), dimension(nv) :: res, d_dlnd, d_dlnT
87 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
88 : integer, intent(out) :: ierr
89 :
90 0 : real(dp) :: rho, logRho, T, logT
91 : logical :: skip
92 :
93 : include 'formats'
94 :
95 0 : T = atemp; logT = alogtemp
96 0 : if (atemp == arg_not_provided .and. alogtemp == arg_not_provided) then
97 0 : ierr = -1; return
98 : end if
99 0 : if (alogtemp == arg_not_provided) logT = log10(T)
100 0 : if (atemp == arg_not_provided) T = exp10(logT)
101 :
102 0 : if (T <= 0) then
103 0 : ierr = -1
104 0 : return
105 : end if
106 :
107 0 : Rho = arho; logrho = alogrho
108 0 : if (arho == arg_not_provided .and. alogrho == arg_not_provided) then
109 0 : ierr = -1; return
110 : end if
111 0 : if (alogrho == arg_not_provided) logRho = log10(Rho)
112 0 : if (arho == arg_not_provided) Rho = exp10(logRho)
113 :
114 0 : if (Rho <= 0) then
115 0 : ierr = -1
116 0 : return
117 : end if
118 0 : if (is_bad(Rho) .or. is_bad(T)) then
119 0 : ierr = -1
120 0 : return
121 : end if
122 :
123 0 : select case(which_eos)
124 : case(i_eos_HELM)
125 : call get_helm_for_eosdt( &
126 : rq% handle, dbg, Z, X, abar, zbar, &
127 : species, chem_id, net_iso, xa, &
128 : rho, logRho, T, logT, 1d0, &
129 : res, d_dlnd, d_dlnT, d_dxa, &
130 0 : skip, ierr)
131 : case(i_eos_OPAL_SCVH)
132 : call get_opal_scvh_for_eosdt( &
133 : rq% handle, dbg, Z, X, abar, zbar, &
134 : species, chem_id, net_iso, xa, &
135 : rho, logRho, T, logT, 1d0, &
136 : res, d_dlnd, d_dlnT, d_dxa, &
137 0 : skip, ierr)
138 : case(i_eos_FreeEOS)
139 : call get_FreeEOS_for_eosdt( &
140 : rq% handle, dbg, Z, X, abar, zbar, &
141 : species, chem_id, net_iso, xa, &
142 : rho, logRho, T, logT, 1d0, &
143 : res, d_dlnd, d_dlnT, d_dxa, &
144 0 : skip, ierr)
145 : case(i_eos_PC)
146 : call get_PC_for_eosdt( &
147 : rq% handle, dbg, Z, X, abar, zbar, &
148 : species, chem_id, net_iso, xa, &
149 : rho, logRho, T, logT, 1d0, &
150 : res, d_dlnd, d_dlnT, d_dxa, &
151 0 : skip, ierr)
152 : case(i_eos_Skye)
153 : call get_Skye_for_eosdt( &
154 : rq% handle, dbg, Z, X, abar, zbar, &
155 : species, chem_id, net_iso, xa, &
156 : rho, logRho, T, logT, 1d0, &
157 : res, d_dlnd, d_dlnT, d_dxa, &
158 0 : skip, ierr)
159 : case(i_eos_CMS)
160 : call get_CMS_for_eosdt( &
161 : rq% handle, dbg, Z, X, abar, zbar, &
162 : species, chem_id, net_iso, xa, &
163 : rho, logRho, T, logT, 1d0, &
164 : res, d_dlnd, d_dlnT, d_dxa, &
165 0 : skip, ierr)
166 : case(i_eos_ideal)
167 : call get_ideal_for_eosdt( &
168 : rq% handle, dbg, Z, X, abar, zbar, &
169 : species, chem_id, net_iso, xa, &
170 : rho, logRho, T, logT, 1d0, &
171 : res, d_dlnd, d_dlnT, d_dxa, &
172 0 : skip, ierr)
173 : case default
174 0 : ierr = -1
175 : end select
176 :
177 0 : if (ierr /= 0) then
178 0 : write(*,*) 'failed in Test_one_eosDT_component', which_eos
179 0 : return
180 : end if
181 :
182 0 : if (skip) then
183 0 : write(*,*) 'skipped - no results Test_one_eosDT_component', which_eos
184 0 : return
185 : end if
186 :
187 : end subroutine Test_one_eosDT_component
188 :
189 :
190 241252 : subroutine Get_eosDT_Results(rq, &
191 : Z, X, abar, zbar, &
192 241252 : species, chem_id, net_iso, xa, &
193 : arho, alogrho, atemp, alogtemp, &
194 241252 : res, d_dlnd, d_dlnT, d_dxa, ierr)
195 : type (EoS_General_Info), pointer :: rq
196 : real(dp), intent(in) :: Z, X, abar, zbar
197 : integer, intent(in) :: species
198 : integer, pointer :: chem_id(:), net_iso(:)
199 : real(dp), intent(in) :: xa(:)
200 : real(dp), intent(in) :: arho, alogrho, atemp, alogtemp
201 : real(dp), intent(inout), dimension(nv) :: res, d_dlnd, d_dlnT
202 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
203 : integer, intent(out) :: ierr
204 :
205 241252 : real(dp) :: rho, logRho, T, logT
206 : logical :: skip, dbg
207 :
208 : include 'formats'
209 :
210 :
211 241252 : T = atemp; logT = alogtemp
212 241252 : if (atemp == arg_not_provided .and. alogtemp == arg_not_provided) then
213 0 : ierr = -1; return
214 : end if
215 241252 : if (alogtemp == arg_not_provided) logT = log10(T)
216 241252 : if (atemp == arg_not_provided) T = exp10(logT)
217 :
218 241252 : if (T <= 0) then
219 0 : ierr = -1
220 0 : return
221 : end if
222 :
223 241252 : Rho = arho; logrho = alogrho
224 241252 : if (arho == arg_not_provided .and. alogrho == arg_not_provided) then
225 0 : ierr = -1; return
226 : end if
227 241252 : if (alogrho == arg_not_provided) logRho = log10(Rho)
228 241252 : if (arho == arg_not_provided) Rho = exp10(logRho)
229 :
230 241252 : if (Rho <= 0) then
231 0 : ierr = -1
232 0 : return
233 : end if
234 241252 : if (is_bad(Rho) .or. is_bad(T)) then
235 0 : ierr = -1
236 0 : return
237 : end if
238 :
239 241252 : dbg = rq% dbg
240 241252 : if (dbg) dbg = & ! check limits
241 : logT >= rq% logT_lo .and. logT <= rq% logT_hi .and. &
242 : logRho >= rq% logRho_lo .and. logRho <= rq% logRho_hi .and. &
243 : X >= rq% X_lo .and. X <= rq% X_hi .and. &
244 0 : Z >= rq% Z_lo .and. Z <= rq% Z_hi
245 :
246 : call get_level0_for_eosdt( &
247 : rq% handle, dbg, Z, X, abar, zbar, &
248 : species, chem_id, net_iso, xa, &
249 : rho, logRho, T, logT, 1d0, &
250 : res, d_dlnd, d_dlnT, d_dxa, &
251 241252 : skip, ierr)
252 241252 : if (skip) ierr = -1
253 241252 : if (ierr /= 0) return
254 :
255 : ! opportunity for the user to modify the eos results
256 241252 : if (rq% use_other_eos_results) then
257 : call rq% other_eos_results( &
258 : rq% handle, &
259 : species, chem_id, net_iso, xa, &
260 : Rho, logRho, T, logT, &
261 0 : res, d_dlnd, d_dlnT, d_dxa, ierr)
262 : end if
263 :
264 241252 : if (eos_test_partials) then
265 0 : eos_test_partials_val = abar
266 0 : eos_test_partials_dval_dx = 0
267 0 : write(*,*) 'eos_test_partials'
268 : end if
269 :
270 : end subroutine Get_eosDT_Results
271 :
272 :
273 0 : subroutine get_other_for_eosdt( &
274 : handle, dbg, Z, X, abar, zbar, &
275 0 : species, chem_id, net_iso, xa, &
276 : rho, logRho, T, logT, remaining_fraction, &
277 0 : res, d_dlnd, d_dlnT, d_dxa, &
278 : skip, ierr)
279 : integer, intent(in) :: handle
280 : logical, intent(in) :: dbg
281 : real(dp), intent(in) :: &
282 : Z, X, abar, zbar, remaining_fraction
283 : integer, intent(in) :: species
284 : integer, pointer :: chem_id(:), net_iso(:)
285 : real(dp), intent(in) :: xa(:)
286 : real(dp), intent(in) :: rho, logRho, T, logT
287 : real(dp), intent(inout), dimension(nv) :: &
288 : res, d_dlnd, d_dlnT
289 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
290 : logical, intent(out) :: skip
291 : integer, intent(out) :: ierr
292 :
293 : type (EoS_General_Info), pointer :: rq
294 :
295 0 : rq => eos_handles(handle)
296 :
297 : ierr = 0
298 :
299 : call rq% other_eos_component( &
300 : handle, &
301 : species, chem_id, net_iso, xa, &
302 : rho, logRho, T, logT, &
303 : res, d_dlnd, d_dlnT, d_dxa, &
304 0 : ierr)
305 :
306 : ! zero all frac components
307 0 : res(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
308 0 : d_dlnd(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
309 0 : d_dlnT(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
310 :
311 0 : skip = .false.
312 :
313 0 : end subroutine get_other_for_eosdt
314 :
315 :
316 241252 : subroutine get_level0_for_eosdt( & ! other
317 : handle, dbg, Z, X, abar, zbar, &
318 241252 : species, chem_id, net_iso, xa, &
319 : rho, logRho, T, logT, remaining_fraction, &
320 241252 : res, d_dlnd, d_dlnT, d_dxa, &
321 : skip, ierr)
322 : integer, intent(in) :: handle
323 : logical, intent(in) :: dbg
324 : real(dp), intent(in) :: Z, X, abar, zbar, remaining_fraction
325 : integer, intent(in) :: species
326 : integer, pointer :: chem_id(:), net_iso(:)
327 : real(dp), intent(in) :: xa(:)
328 : real(dp), intent(in) :: rho, logRho, T, logT
329 : real(dp), intent(inout), dimension(nv) :: &
330 : res, d_dlnd, d_dlnT
331 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
332 : logical, intent(out) :: skip
333 : integer, intent(out) :: ierr
334 :
335 241252 : real(dp) :: frac, d_frac_dlogT, d_frac_dlogRho
336 241252 : real(dp) :: alfa, d_alfa_dlogT, d_alfa_dlogRho
337 : type (EoS_General_Info), pointer :: rq
338 : procedure (get_values_for_eosdt_interface), pointer :: get_1st, get_2nd
339 :
340 : include 'formats'
341 :
342 241252 : ierr = 0
343 241252 : rq => eos_handles(handle)
344 :
345 241252 : if (rq% use_other_eos_component) then
346 : call rq% other_eos_frac( &
347 : handle, &
348 : species, chem_id, net_iso, xa, &
349 : rho, logRho, T, logT, &
350 : frac, d_frac_dlogRho, d_frac_dlogT, &
351 0 : ierr)
352 0 : if (ierr /= 0) return
353 0 : alfa = 1d0 - frac
354 0 : d_alfa_dlogT = -d_frac_dlogT
355 0 : d_alfa_dlogRho = -d_frac_dlogRho
356 : else
357 241252 : alfa = 1d0 ! no other
358 241252 : d_alfa_dlogT = 0d0
359 241252 : d_alfa_dlogRho = 0d0
360 : end if
361 :
362 241252 : if (dbg) write(*,1) 'other', (1d0 - alfa)*remaining_fraction
363 :
364 241252 : get_1st => get_other_for_eosdt
365 241252 : get_2nd => get_level1_for_eosdt
366 : call combine_for_eosdt( &
367 : get_1st, get_2nd, alfa*remaining_fraction, &
368 : alfa, d_alfa_dlogT, d_alfa_dlogRho, &
369 : rq, dbg, Z, X, abar, zbar, &
370 : species, chem_id, net_iso, xa, &
371 : rho, logRho, T, logT, &
372 : res, d_dlnd, d_dlnT, d_dxa, &
373 241252 : skip, ierr)
374 241252 : if (ierr /= 0 .and. rq% okay_to_convert_ierr_to_skip) then
375 0 : skip = .true.
376 0 : ierr = 0
377 : end if
378 :
379 : end subroutine get_level0_for_eosdt
380 :
381 :
382 241252 : subroutine get_level1_for_eosdt( & ! CMS
383 : handle, dbg, Z, X, abar, zbar, &
384 241252 : species, chem_id, net_iso, xa, &
385 : rho, logRho, T, logT, remaining_fraction, &
386 241252 : res, d_dlnd, d_dlnT, d_dxa, &
387 : skip, ierr)
388 : use eoscms_eval, only: Get_CMS_alfa, get_CMS_for_eosdt
389 : integer, intent(in) :: handle
390 : logical, intent(in) :: dbg
391 : real(dp), intent(in) :: Z, X, abar, zbar, remaining_fraction
392 : integer, intent(in) :: species
393 : integer, pointer :: chem_id(:), net_iso(:)
394 : real(dp), intent(in) :: xa(:)
395 : real(dp), intent(in) :: rho, logRho, T, logT
396 : real(dp), intent(inout), dimension(nv) :: &
397 : res, d_dlnd, d_dlnT
398 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
399 : logical, intent(out) :: skip
400 : integer, intent(out) :: ierr
401 :
402 241252 : real(dp) :: alfa, d_alfa_dlogT, d_alfa_dlogRho
403 : type (EoS_General_Info), pointer :: rq
404 : procedure (get_values_for_eosdt_interface), pointer :: get_1st, get_2nd
405 :
406 : include 'formats'
407 :
408 241252 : ierr = 0
409 241252 : rq => eos_handles(handle)
410 :
411 241252 : if (rq% use_CMS) then
412 : call Get_CMS_alfa( &
413 : rq, logRho, logT, Z, abar, zbar, &
414 : alfa, d_alfa_dlogT, d_alfa_dlogRho, &
415 0 : ierr)
416 0 : if (ierr /= 0) return
417 : else
418 241252 : alfa = 1d0 ! no CMS
419 241252 : d_alfa_dlogT = 0d0
420 241252 : d_alfa_dlogRho = 0d0
421 : end if
422 :
423 241252 : if (dbg) write(*,1) 'CMS', (1d0 - alfa)*remaining_fraction
424 :
425 241252 : get_1st => get_CMS_for_eosdt
426 241252 : get_2nd => get_level2_for_eosdt
427 : call combine_for_eosdt( &
428 : get_1st, get_2nd, alfa*remaining_fraction, &
429 : alfa, d_alfa_dlogT, d_alfa_dlogRho, &
430 : rq, dbg, Z, X, abar, zbar, &
431 : species, chem_id, net_iso, xa, &
432 : rho, logRho, T, logT, &
433 : res, d_dlnd, d_dlnT, d_dxa, &
434 241252 : skip, ierr)
435 241252 : if (ierr /= 0 .and. rq% okay_to_convert_ierr_to_skip) then
436 0 : skip = .true.
437 0 : ierr = 0
438 : end if
439 :
440 241252 : end subroutine get_level1_for_eosdt
441 :
442 :
443 241252 : subroutine get_level2_for_eosdt( & ! Skye
444 : handle, dbg, Z, X, abar, zbar, &
445 241252 : species, chem_id, net_iso, xa, &
446 : rho, logRho, T, logT, remaining_fraction, &
447 241252 : res, d_dlnd, d_dlnT, d_dxa, &
448 : skip, ierr)
449 : integer, intent(in) :: handle
450 : logical, intent(in) :: dbg
451 : real(dp), intent(in) :: Z, X, abar, zbar, remaining_fraction
452 : integer, intent(in) :: species
453 : integer, pointer :: chem_id(:), net_iso(:)
454 : real(dp), intent(in) :: xa(:)
455 : real(dp), intent(in) :: rho, logRho, T, logT
456 : real(dp), intent(inout), dimension(nv) :: &
457 : res, d_dlnd, d_dlnT
458 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
459 : logical, intent(out) :: skip
460 : integer, intent(out) :: ierr
461 :
462 241252 : real(dp) :: alfa, d_alfa_dlogT, d_alfa_dlogRho
463 : type (EoS_General_Info), pointer :: rq
464 : procedure (get_values_for_eosdt_interface), pointer :: get_1st, get_2nd
465 :
466 : include 'formats'
467 :
468 241252 : ierr = 0
469 241252 : rq => eos_handles(handle)
470 :
471 241252 : if (rq% use_Skye) then
472 241252 : if (rq% use_simple_Skye_blends) then
473 : call Get_Skye_alfa_simple( &
474 : rq, logRho, logT, Z, abar, zbar, &
475 : alfa, d_alfa_dlogT, d_alfa_dlogRho, &
476 0 : ierr)
477 : else
478 : call Get_Skye_alfa( &
479 : rq, logRho, logT, Z, abar, zbar, &
480 : alfa, d_alfa_dlogT, d_alfa_dlogRho, &
481 241252 : ierr)
482 : end if
483 241252 : if (ierr /= 0) return
484 : else
485 0 : alfa = 1d0 ! no Skye
486 0 : d_alfa_dlogT = 0d0
487 0 : d_alfa_dlogRho = 0d0
488 : end if
489 :
490 241252 : if (dbg) write(*,1) 'Skye', (1d0 - alfa)*remaining_fraction
491 241252 : get_1st => get_Skye_for_eosdt
492 :
493 241252 : get_2nd => get_level3_for_eosdt
494 : call combine_for_eosdt( &
495 : get_1st, get_2nd, alfa*remaining_fraction, &
496 : alfa, d_alfa_dlogT, d_alfa_dlogRho, &
497 : rq, dbg, Z, X, abar, zbar, &
498 : species, chem_id, net_iso, xa, &
499 : rho, logRho, T, logT, &
500 : res, d_dlnd, d_dlnT, d_dxa, &
501 241252 : skip, ierr)
502 241252 : if (ierr /= 0 .and. rq% okay_to_convert_ierr_to_skip) then
503 0 : skip = .true.
504 0 : ierr = 0
505 : end if
506 :
507 241252 : end subroutine get_level2_for_eosdt
508 :
509 :
510 210666 : subroutine get_level3_for_eosdt( & ! PC
511 : handle, dbg, Z, X, abar, zbar, &
512 210666 : species, chem_id, net_iso, xa, &
513 : rho, logRho, T, logT, remaining_fraction, &
514 210666 : res, d_dlnd, d_dlnT, d_dxa, &
515 : skip, ierr)
516 : use eospc_eval, only: Get_PC_alfa, get_PC_for_eosdt
517 : integer, intent(in) :: handle
518 : logical, intent(in) :: dbg
519 : real(dp), intent(in) :: Z, X, abar, zbar, remaining_fraction
520 : integer, intent(in) :: species
521 : integer, pointer :: chem_id(:), net_iso(:)
522 : real(dp), intent(in) :: xa(:)
523 : real(dp), intent(in) :: rho, logRho, T, logT
524 : real(dp), intent(inout), dimension(nv) :: &
525 : res, d_dlnd, d_dlnT
526 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
527 : logical, intent(out) :: skip
528 : integer, intent(out) :: ierr
529 :
530 210666 : real(dp) :: alfa, d_alfa_dlogT, d_alfa_dlogRho
531 : type (EoS_General_Info), pointer :: rq
532 : procedure (get_values_for_eosdt_interface), pointer :: get_1st, get_2nd
533 :
534 : include 'formats'
535 :
536 210666 : ierr = 0
537 210666 : rq => eos_handles(handle)
538 :
539 210666 : if (rq% use_PC) then
540 : call Get_PC_alfa( &
541 : rq, logRho, logT, Z, abar, zbar, &
542 : alfa, d_alfa_dlogT, d_alfa_dlogRho, &
543 0 : ierr)
544 0 : if (ierr /= 0) return
545 : else
546 210666 : alfa = 1d0 ! no PC
547 210666 : d_alfa_dlogT = 0d0
548 210666 : d_alfa_dlogRho = 0d0
549 : end if
550 :
551 210666 : if (dbg) write(*,1) 'PC', (1d0 - alfa)*remaining_fraction
552 210666 : get_1st => get_PC_for_eosdt
553 :
554 210666 : get_2nd => get_level4_for_eosdt
555 : call combine_for_eosdt( &
556 : get_1st, get_2nd, alfa*remaining_fraction, &
557 : alfa, d_alfa_dlogT, d_alfa_dlogRho, &
558 : rq, dbg, Z, X, abar, zbar, &
559 : species, chem_id, net_iso, xa, &
560 : rho, logRho, T, logT, &
561 : res, d_dlnd, d_dlnT, d_dxa, &
562 210666 : skip, ierr)
563 210666 : if (ierr /= 0 .and. rq% okay_to_convert_ierr_to_skip) then
564 0 : skip = .true.
565 0 : ierr = 0
566 : end if
567 :
568 210666 : end subroutine get_level3_for_eosdt
569 :
570 :
571 210666 : subroutine get_level4_for_eosdt( & ! FreeEOS
572 : handle, dbg, Z, X, abar, zbar, &
573 210666 : species, chem_id, net_iso, xa, &
574 : rho, logRho, T, logT, remaining_fraction, &
575 210666 : res, d_dlnd, d_dlnT, d_dxa, &
576 : skip, ierr)
577 210666 : use skye, only: get_Skye_for_eosdt, get_Skye_alfa
578 : integer, intent(in) :: handle
579 : logical, intent(in) :: dbg
580 : real(dp), intent(in) :: Z, X, abar, zbar, remaining_fraction
581 : integer, intent(in) :: species
582 : integer, pointer :: chem_id(:), net_iso(:)
583 : real(dp), intent(in) :: xa(:)
584 : real(dp), intent(in) :: rho, logRho, T, logT
585 : real(dp), intent(inout), dimension(nv) :: &
586 : res, d_dlnd, d_dlnT
587 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
588 : logical, intent(out) :: skip
589 : integer, intent(out) :: ierr
590 :
591 210666 : real(dp) :: alfa, d_alfa_dlogT, d_alfa_dlogRho
592 : type (EoS_General_Info), pointer :: rq
593 : procedure (get_values_for_eosdt_interface), pointer :: get_1st, get_2nd
594 :
595 : include 'formats'
596 :
597 210666 : ierr = 0
598 210666 : rq => eos_handles(handle)
599 :
600 210666 : if (rq% use_FreeEOS) then
601 : call Get_FreeEOS_alfa( &
602 : rq, dbg, logRho, logT, Z, abar, zbar, &
603 : alfa, d_alfa_dlogT, d_alfa_dlogRho, &
604 210666 : ierr)
605 210666 : if (ierr /= 0) return
606 : else
607 0 : alfa = 1d0 ! no FreeEOS
608 0 : d_alfa_dlogT = 0d0
609 0 : d_alfa_dlogRho = 0d0
610 : end if
611 :
612 210666 : if (dbg) write(*,1) 'FreeEOS', (1d0 - alfa)*remaining_fraction
613 210666 : get_1st => get_FreeEOS_for_eosdt
614 :
615 210666 : get_2nd => get_level5_for_eosdt
616 : call combine_for_eosdt( &
617 : get_1st, get_2nd, alfa*remaining_fraction, &
618 : alfa, d_alfa_dlogT, d_alfa_dlogRho, &
619 : rq, dbg, Z, X, abar, zbar, &
620 : species, chem_id, net_iso, xa, &
621 : rho, logRho, T, logT, &
622 : res, d_dlnd, d_dlnT, d_dxa, &
623 210666 : skip, ierr)
624 210666 : if (ierr /= 0 .and. rq% okay_to_convert_ierr_to_skip) then
625 0 : skip = .true.
626 0 : ierr = 0
627 : end if
628 :
629 210666 : end subroutine get_level4_for_eosdt
630 :
631 :
632 70 : subroutine get_level5_for_eosdt( & ! OPAL/SCVH
633 : handle, dbg, Z, X, abar, zbar, &
634 70 : species, chem_id, net_iso, xa, &
635 : rho, logRho, T_in, logT_in, remaining_fraction, &
636 70 : res, d_dlnd, d_dlnT, d_dxa, &
637 : skip, ierr)
638 : integer, intent(in) :: handle
639 : logical, intent(in) :: dbg
640 : real(dp), intent(in) :: &
641 : Z, X, abar, zbar, remaining_fraction
642 : integer, intent(in) :: species
643 : integer, pointer :: chem_id(:), net_iso(:)
644 : real(dp), intent(in) :: xa(:)
645 : real(dp), intent(in) :: rho, logRho, T_in, logT_in
646 : real(dp), intent(inout), dimension(nv) :: &
647 : res, d_dlnd, d_dlnT
648 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
649 : logical, intent(out) :: skip
650 : integer, intent(out) :: ierr
651 :
652 70 : real(dp) :: alfa, d_alfa_dlogT, d_alfa_dlogRho, &
653 : logT_HELM, T_HELM, logQ, logQ2, T, logT
654 : type (EoS_General_Info), pointer :: rq
655 : procedure (get_values_for_eosdt_interface), pointer :: get_1st, get_2nd
656 :
657 : include 'formats'
658 :
659 70 : ierr = 0
660 70 : rq => eos_handles(handle)
661 :
662 70 : T = T_in
663 70 : logT = logT_in
664 :
665 70 : if (rq% use_OPAL_SCVH) then
666 : call get_opal_scvh_alfa_and_partials( &
667 : rq, logT, logRho, Z, &
668 70 : alfa, d_alfa_dlogRho, d_alfa_dlogT, ierr)
669 70 : if (ierr /= 0) return
670 : else
671 0 : alfa = 1d0 ! no OPAL_SCVH
672 0 : d_alfa_dlogT = 0d0
673 0 : d_alfa_dlogRho = 0d0
674 : end if
675 :
676 70 : if (dbg) write(*,1) 'OPAL/SCVH', (1d0 - alfa)*remaining_fraction
677 70 : get_1st => get_opal_scvh_for_eosdt
678 :
679 70 : get_2nd => get_level6_for_eosdt
680 : call combine_for_eosdt( &
681 : get_1st, get_2nd, alfa*remaining_fraction, &
682 : alfa, d_alfa_dlogT, d_alfa_dlogRho, &
683 : rq, dbg, Z, X, abar, zbar, &
684 : species, chem_id, net_iso, xa, &
685 : rho, logRho, T, logT, &
686 : res, d_dlnd, d_dlnT, d_dxa, &
687 70 : skip, ierr)
688 70 : if (ierr /= 0 .and. rq% okay_to_convert_ierr_to_skip) then
689 0 : skip = .true.
690 0 : ierr = 0
691 : end if
692 :
693 210666 : end subroutine get_level5_for_eosdt
694 :
695 0 : subroutine get_level6_for_eosdt( & ! HELM/ideal
696 : handle, dbg, Z, X, abar, zbar, &
697 0 : species, chem_id, net_iso, xa, &
698 : rho, logRho, T_in, logT_in, remaining_fraction, &
699 0 : res, d_dlnd, d_dlnT, d_dxa, &
700 : skip, ierr)
701 : integer, intent(in) :: handle
702 : logical, intent(in) :: dbg
703 : real(dp), intent(in) :: &
704 : Z, X, abar, zbar, remaining_fraction
705 : integer, intent(in) :: species
706 : integer, pointer :: chem_id(:), net_iso(:)
707 : real(dp), intent(in) :: xa(:)
708 : real(dp), intent(in) :: rho, logRho, T_in, logT_in
709 : real(dp), intent(inout), dimension(nv) :: &
710 : res, d_dlnd, d_dlnT
711 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
712 : logical, intent(out) :: skip
713 : integer, intent(out) :: ierr
714 :
715 : real(dp) :: alfa, d_alfa_dlogT, d_alfa_dlogRho, &
716 : logT_HELM, T_HELM, logQ, logQ2, T, logT
717 : type (EoS_General_Info), pointer :: rq
718 : procedure (get_values_for_eosdt_interface), pointer :: get_1st, get_2nd
719 :
720 : include 'formats'
721 :
722 : ierr = 0
723 0 : rq => eos_handles(handle)
724 :
725 0 : T = T_in
726 0 : logT = logT_in
727 :
728 0 : call get_HELM_alfa(rq, logRho, logT, alfa, d_alfa_dlogRho, d_alfa_dlogT, ierr)
729 :
730 0 : if (dbg) write(*,1) 'HELM', (1d0 - alfa)*remaining_fraction
731 :
732 0 : get_1st => get_helm_for_eosdt
733 0 : get_2nd => get_ideal_for_eosdt
734 : call combine_for_eosdt( &
735 : get_1st, get_2nd, alfa*remaining_fraction, &
736 : alfa, d_alfa_dlogT, d_alfa_dlogRho, &
737 : rq, dbg, Z, X, abar, zbar, &
738 : species, chem_id, net_iso, xa, &
739 : rho, logRho, T, logT, &
740 : res, d_dlnd, d_dlnT, d_dxa, &
741 0 : skip, ierr)
742 0 : if (ierr /= 0 .and. rq% okay_to_convert_ierr_to_skip) then
743 0 : skip = .true.
744 0 : ierr = 0
745 : end if
746 :
747 0 : end subroutine get_level6_for_eosdt
748 :
749 0 : subroutine get_HELM_alfa( &
750 : rq, logRho, logT, alfa, d_alfa_dlogT, d_alfa_dlogRho, ierr)
751 : use const_def, only: dp
752 : use eos_blend
753 : use auto_diff
754 : type (EoS_General_Info), pointer :: rq
755 : real(dp), intent(in) :: logRho, logT
756 : real(dp), intent(out) :: alfa, d_alfa_dlogT, d_alfa_dlogRho
757 : integer, intent(out) :: ierr
758 :
759 : logical :: contained
760 : type(auto_diff_real_2var_order1) :: p(2), blend, dist
761 :
762 : ! Blend parameters
763 0 : real(dp) :: helm_blend_width
764 : integer, parameter :: num_points = 4
765 0 : real(dp) :: bounds(4,2)
766 : type (Helm_Table), pointer :: ht
767 :
768 0 : ierr = 0
769 0 : ht => eos_ht
770 0 : helm_blend_width = 0.1d0
771 :
772 0 : bounds(1,1) = ht% logdlo
773 0 : bounds(1,2) = ht% logthi
774 :
775 0 : bounds(2,1) = ht% logdlo
776 0 : bounds(2,2) = ht% logtlo
777 :
778 0 : bounds(3,1) = ht% logdhi
779 0 : bounds(3,2) = ht% logtlo
780 :
781 0 : bounds(4,1) = ht% logdhi
782 0 : bounds(4,2) = ht% logthi
783 :
784 : ! Set up auto_diff point
785 0 : p(1) = logRho
786 0 : p(1)%d1val1 = 1d0
787 0 : p(2) = logT
788 0 : p(2)%d1val2 = 1d0
789 :
790 0 : contained = is_contained(num_points, bounds, p)
791 0 : dist = min_distance_to_polygon(num_points, bounds, p)
792 :
793 0 : if (contained) then ! Make distance negative for points inside the polygon
794 0 : dist = -dist
795 : end if
796 :
797 0 : dist = dist / helm_blend_width
798 0 : blend = max(dist, 0d0)
799 0 : blend = min(blend, 1d0)
800 :
801 0 : alfa = blend%val
802 0 : d_alfa_dlogRho = blend%d1val1
803 0 : d_alfa_dlogT = blend%d1val2
804 :
805 0 : end subroutine get_HELM_alfa
806 :
807 :
808 210666 : subroutine Get_FreeEOS_alfa( &
809 : rq, dbg, logRho, logT, Z, abar, zbar, &
810 : alfa, d_alfa_dlogT, d_alfa_dlogRho, &
811 : ierr)
812 0 : use const_def, only: dp
813 : use auto_diff
814 : type (EoS_General_Info), pointer :: rq
815 : logical, intent(in) :: dbg
816 : real(dp), intent(in) :: logRho, logT, Z, abar, zbar
817 : real(dp), intent(out) :: alfa, d_alfa_dlogT, d_alfa_dlogRho
818 : integer, intent(out) :: ierr
819 210666 : real(dp) :: logQ_cut_lo, logQ_cut_hi
820 : type(auto_diff_real_2var_order1) :: logT_auto, logRho_auto, logQ_auto
821 : type(auto_diff_real_2var_order1) :: blend, blend_logT, blend_logRho, blend_logQ
822 : type(auto_diff_real_2var_order1) :: blend_cut, blend_logT_cut, blend_logRho_cut, blend_logQ_cut
823 :
824 : include 'formats'
825 :
826 210666 : ierr = 0
827 :
828 : ! auto diff
829 : ! var1: logRho
830 : ! var2: logT
831 :
832 210666 : logRho_auto% val = logRho
833 210666 : logRho_auto% d1val1 = 1
834 210666 : logRho_auto% d1val2 = 0
835 :
836 210666 : logT_auto% val = logT
837 210666 : logT_auto% d1val1 = 0
838 210666 : logT_auto% d1val2 = 1
839 :
840 210666 : logQ_auto = logRho_auto - 2d0*logT_auto + 12d0
841 :
842 : ! The FreeEOS region usually looks like
843 : !
844 : ! ________C________
845 : ! / |
846 : ! / D
847 : ! / __E__|
848 : ! B /
849 : ! / /
850 : ! / F
851 : ! / /
852 : ! / /
853 : ! / /
854 : ! / |
855 : ! / G
856 : ! /________A_______|
857 : !
858 : ! where blend A and C come from blend_logT,
859 : ! blend B comes from blend_logQ,
860 : ! blend D usually comes from PC/Skye,
861 : ! blend E comes from blend_logT_cut,
862 : ! blend F comes from blend_logQ_cut,
863 : ! blend G comes from blend_logRho_cut.
864 :
865 :
866 : ! logT blend
867 210666 : if (logT_auto < rq% logT_min_FreeEOS_lo) then
868 0 : blend_logT = 0d0
869 210666 : else if (logT_auto < rq% logT_min_FreeEOS_hi) then
870 6 : blend_logT = (logT_auto - rQ% logT_min_FreeEOS_lo) / (rq% logT_min_FreeEOS_hi - rq% logT_min_FreeEOS_lo)
871 210660 : else if (logT_auto < rq% logT_max_FreeEOS_lo) then
872 210660 : blend_logT = 1d0
873 0 : else if (logT_auto < rq% logT_max_FreeEOS_hi) then
874 0 : blend_logT = (logT_auto - rQ% logT_max_FreeEOS_hi) / (rq% logT_max_FreeEOS_lo - rq% logT_max_FreeEOS_hi)
875 : else
876 0 : blend_logT = 0
877 : end if
878 :
879 :
880 : ! logRho blend
881 210666 : if (logRho_auto < rq% logRho_min_FreeEOS_lo) then
882 0 : blend_logRho = 0d0
883 : else if (logRho_auto < rq% logRho_min_FreeEOS_lo) then
884 : blend_logRho = (logRho_auto - rQ% logRho_min_FreeEOS_lo) / (rq% logRho_min_FreeEOS_hi - rq% logRho_min_FreeEOS_lo)
885 210666 : else if (logRho_auto < rq% logRho_max_FreeEOS_lo) then
886 210666 : blend_logRho = 1d0
887 0 : else if (logRho_auto < rq% logRho_max_FreeEOS_hi) then
888 0 : blend_logRho = (logRho_auto - rQ% logRho_max_FreeEOS_hi) / (rq% logRho_max_FreeEOS_lo - rq% logRho_max_FreeEOS_hi)
889 : else
890 0 : blend_logRho = 0
891 : end if
892 :
893 :
894 : ! logQ blend
895 210666 : if (logQ_auto < rq% logQ_min_FreeEOS_lo) then
896 0 : blend_logQ = 0d0
897 210666 : else if (logQ_auto < rq% logQ_min_FreeEOS_hi) then
898 0 : blend_logQ = (logQ_auto - rQ% logQ_min_FreeEOS_lo) / (rq% logQ_min_FreeEOS_hi - rq% logQ_min_FreeEOS_lo)
899 210666 : else if (logQ_auto < rq% logQ_max_FreeEOS_lo) then
900 210666 : blend_logQ = 1d0
901 0 : else if (logQ_auto < rq% logQ_max_FreeEOS_hi) then
902 0 : blend_logQ = (logQ_auto - rQ% logQ_max_FreeEOS_hi) / (rq% logQ_max_FreeEOS_lo - rq% logQ_max_FreeEOS_hi)
903 : else
904 0 : blend_logQ = 0
905 : end if
906 :
907 :
908 : ! cut blend logRho
909 210666 : if (logRho_auto < rq% logRho_cut_FreeEOS_lo) then
910 51021 : blend_logRho_cut = 1d0
911 159645 : else if (logRho_auto < rq% logRho_cut_FreeEOS_hi) then
912 14245 : blend_logRho_cut = (logRho_auto - rQ% logRho_cut_FreeEOS_hi) / (rq% logRho_cut_FreeEOS_lo - rq% logRho_cut_FreeEOS_hi)
913 : else
914 145400 : blend_logRho_cut = 0d0
915 : end if
916 :
917 : ! cut blend logT
918 210666 : if (logT_auto < rq% logT_cut_FreeEOS_lo) then
919 110826 : blend_logT_cut = 0d0
920 99840 : else if (logT_auto < rq% logT_cut_FreeEOS_hi) then
921 3432 : blend_logT_cut = (logT_auto - rQ% logT_cut_FreeEOS_lo) / (rq% logT_cut_FreeEOS_hi - rq% logT_cut_FreeEOS_lo)
922 : else
923 96408 : blend_logT_cut = 1d0
924 : end if
925 :
926 : ! cut blend logQ
927 210666 : if (Z <= rq% logQ_cut_FreeEOS_lo_Z_max) then
928 210666 : logQ_cut_hi = rq% logQ_cut_lo_Z_FreeEOS_hi
929 210666 : logQ_cut_lo = rq% logQ_cut_lo_Z_FreeEOS_lo
930 : else
931 0 : logQ_cut_hi = rq% logQ_cut_hi_Z_FreeEOS_hi
932 0 : logQ_cut_lo = rq% logQ_cut_hi_Z_FreeEOS_lo
933 : end if
934 :
935 210666 : if (logQ_auto < logQ_cut_lo) then
936 210598 : blend_logQ_cut = 1d0
937 68 : else if (logQ_auto < logQ_cut_hi) then
938 10 : blend_logQ_cut = (logQ_auto - logQ_cut_hi) / (logQ_cut_lo - logQ_cut_hi)
939 : else
940 58 : blend_logQ_cut = 0d0
941 : end if
942 :
943 : ! combine cut blends
944 210666 : blend_cut = 0
945 210666 : if (blend_logT_cut == 0) then
946 110826 : if (blend_logQ_cut == 0) then
947 58 : blend_cut = blend_logRho_cut
948 110768 : else if (blend_logQ_cut == 1) then
949 110758 : blend_cut = 1d0
950 : else
951 10 : blend_cut = max(blend_logRho_cut, blend_logQ_cut)
952 : end if
953 99840 : elseif (blend_logT_cut == 1) then
954 96408 : blend_cut = 1d0
955 : else
956 3432 : if (blend_logQ_cut == 0) then
957 0 : blend_cut = blend_logT_cut
958 3432 : else if (blend_logQ_cut == 1) then
959 3432 : blend_cut = 1d0
960 : else
961 0 : blend_cut = max(blend_logT_cut, blend_logQ_cut)
962 : end if
963 : end if
964 :
965 : ! combine all blends
966 210666 : blend = blend_cut * blend_logT * blend_logRho * blend_logQ
967 :
968 : ! unpack auto_diff
969 210666 : alfa = 1d0 - blend% val
970 210666 : d_alfa_dlogRho = -blend% d1val1
971 210666 : d_alfa_dlogT = -blend% d1val2
972 :
973 210666 : end subroutine Get_FreeEOS_alfa
974 :
975 :
976 70 : subroutine get_opal_scvh_for_eosdt( &
977 : handle, dbg, Z, X, abar, zbar, &
978 70 : species, chem_id, net_iso, xa, &
979 : rho, logRho, T, logT, remaining_fraction, &
980 70 : res, d_dlnd, d_dlnT, d_dxa, &
981 : skip, ierr)
982 : integer, intent(in) :: handle
983 : logical, intent(in) :: dbg
984 : real(dp), intent(in) :: &
985 : Z, X, abar, zbar, remaining_fraction
986 : integer, intent(in) :: species
987 : integer, pointer :: chem_id(:), net_iso(:)
988 : real(dp), intent(in) :: xa(:)
989 : real(dp), intent(in) :: rho, logRho, T, logT
990 : real(dp), intent(inout), dimension(nv) :: &
991 : res, d_dlnd, d_dlnT
992 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
993 : logical, intent(out) :: skip
994 : integer, intent(out) :: ierr
995 : call get1_for_eosdt( &
996 : handle, eosdt_OPAL_SCVH, dbg, &
997 : Z, X, abar, zbar, &
998 : species, chem_id, net_iso, xa, &
999 : rho, logRho, T, logT, remaining_fraction, &
1000 : res, d_dlnd, d_dlnT, d_dxa, &
1001 70 : skip, ierr)
1002 :
1003 : ! zero phase information
1004 280 : res(i_phase:i_latent_ddlnRho) = 0d0
1005 280 : d_dlnT(i_phase:i_latent_ddlnRho) = 0d0
1006 280 : d_dlnd(i_phase:i_latent_ddlnRho) = 0d0
1007 :
1008 : ! zero all components
1009 560 : res(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
1010 560 : d_dlnd(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
1011 560 : d_dlnT(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
1012 :
1013 : ! mark this one
1014 70 : res(i_frac_OPAL_SCVH) = 1.0d0
1015 :
1016 210666 : end subroutine get_opal_scvh_for_eosdt
1017 :
1018 :
1019 210602 : subroutine get_FreeEOS_for_eosdt( &
1020 : handle, dbg, Z, X, abar, zbar, &
1021 210602 : species, chem_id, net_iso, xa, &
1022 : rho, logRho, T, logT, remaining_fraction, &
1023 210602 : res, d_dlnd, d_dlnT, d_dxa, &
1024 : skip, ierr)
1025 : integer, intent(in) :: handle
1026 : logical, intent(in) :: dbg
1027 : real(dp), intent(in) :: &
1028 : Z, X, abar, zbar, remaining_fraction
1029 : integer, intent(in) :: species
1030 : integer, pointer :: chem_id(:), net_iso(:)
1031 : real(dp), intent(in) :: xa(:)
1032 : real(dp), intent(in) :: rho, logRho, T, logT
1033 : real(dp), intent(inout), dimension(nv) :: &
1034 : res, d_dlnd, d_dlnT
1035 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
1036 : logical, intent(out) :: skip
1037 : integer, intent(out) :: ierr
1038 : call get1_for_eosdt( &
1039 : handle, eosdt_max_FreeEOS, dbg, Z, X, abar, zbar, &
1040 : species, chem_id, net_iso, xa, &
1041 : rho, logRho, T, logT, remaining_fraction, &
1042 : res, d_dlnd, d_dlnT, d_dxa, &
1043 210602 : skip, ierr)
1044 :
1045 : ! zero phase information
1046 842408 : res(i_phase:i_latent_ddlnRho) = 0d0
1047 842408 : d_dlnT(i_phase:i_latent_ddlnRho) = 0d0
1048 842408 : d_dlnd(i_phase:i_latent_ddlnRho) = 0d0
1049 :
1050 : ! zero all components
1051 1684816 : res(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
1052 1684816 : d_dlnd(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
1053 1684816 : d_dlnT(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
1054 :
1055 : ! mark this one
1056 210602 : res(i_frac_FreeEOS) = 1.0d0
1057 :
1058 210602 : end subroutine get_FreeEOS_for_eosdt
1059 :
1060 :
1061 70 : subroutine get_opal_scvh_alfa_and_partials( &
1062 : rq, logT, logRho, Z, alfa, d_alfa_dlogRho, d_alfa_dlogT, ierr)
1063 : type (EoS_General_Info), pointer :: rq
1064 : real(dp), intent(in) :: logT, logRho, Z
1065 : real(dp), intent(out) :: alfa, d_alfa_dlogRho, d_alfa_dlogT
1066 : integer, intent(out) :: ierr
1067 :
1068 : integer :: iregion
1069 : real(dp) :: logRho1_max, logRho1, logRho2, logRho5, logRho6, logRho7, &
1070 70 : logRho8, logT5, logT6, logT3, logT4
1071 70 : real(dp) :: logQ1, logQ2, logQ3, logQ4, logQmax, Z_all_HELM, Z_no_HELM
1072 70 : real(dp) :: beta, &
1073 70 : logT1, logT2, logT7, logT8, logRho3, logRho4
1074 70 : real(dp) :: logQ, A, B, dA_dlnT, dA_dlnRho, dB_dlnT, dB_dlnRho
1075 70 : real(dp) :: c_dx, c_dy, d_dx_dlogT, d_dx_dlogRho, d_dy_dlogT, d_dy_dlogRho
1076 : real(dp), parameter :: tiny = 1d-20
1077 :
1078 : include 'formats'
1079 :
1080 70 : logRho1_max = 3.71d0
1081 :
1082 70 : logQ1 = 5.69d0 ! SCVH full off for logQ > this
1083 70 : logQ2 = 5.68d0 ! must have logQ2 < logQ1
1084 70 : logQmax = rq% logQ_max_OPAL_SCVH ! 5.3
1085 70 : logQ3 = rq% logQ_min_OPAL_SCVH ! -8.0
1086 70 : logQ4 = rq% logQ_min_OPAL_SCVH ! -8.0
1087 :
1088 70 : logRho5 = rq% logRho_min_OPAL_SCVH_limit ! -14.299
1089 70 : logRho6 = logRho5 - 1d-3 ! -14.3
1090 70 : logRho7 = -14.90d0
1091 70 : logRho8 = -14.99d0
1092 :
1093 70 : logRho1 = rq% logRho1_OPAL_SCVH_limit ! 3.50
1094 70 : if (logRho1 > logRho1_max) then
1095 0 : write(*,*) 'sorry: value for logRho1_OPAL_SCVH_limit is too large. max allowed is ', &
1096 0 : logRho1_max
1097 0 : ierr = -1
1098 0 : return
1099 : end if
1100 70 : logRho2 = rq% logRho2_OPAL_SCVH_limit ! 3.48
1101 :
1102 70 : logT8 = rq% logT_low_all_HELM ! 2.2
1103 70 : logT7 = rq% logT_low_all_SCVH ! 2.3
1104 70 : logT6 = 4.890d0
1105 70 : logT5 = 4.899d0 ! problems with blend here so just jump
1106 70 : logT2 = rq% logT_all_OPAL ! 7.6
1107 70 : logT1 = rq% logT_all_HELM ! 7.7
1108 :
1109 70 : Z_all_HELM = rq% Z_all_HELM
1110 70 : Z_no_HELM = rq% Z_all_OPAL
1111 :
1112 70 : if (logT >= logT1) then ! just use other
1113 :
1114 0 : alfa = 1d0
1115 : beta = 0d0
1116 :
1117 : else
1118 :
1119 70 : logT3 = (logRho1 - logQ1 + 12d0)/2d0
1120 70 : logT4 = (logRho2 - logQ2 + 12d0)/2d0
1121 :
1122 70 : logRho3 = logQ1 + 2*logT7 - 12d0
1123 70 : logRho4 = logQ2 + 2*logT7 - 12d0
1124 :
1125 : if (.false.) then
1126 : write(*,'(A)')
1127 : write(*,1) 'logRho1', logRho1
1128 : write(*,1) 'logRho2', logRho2
1129 : write(*,1) 'logRho3', logRho3
1130 : write(*,1) 'logRho4', logRho4
1131 : write(*,1) 'logRho5', logRho5
1132 : write(*,1) 'logRho6', logRho6
1133 : write(*,1) 'logRho7', logRho7
1134 : write(*,1) 'logRho8', logRho8
1135 : write(*,'(A)')
1136 : write(*,1) 'logT1', logT1
1137 : write(*,1) 'logT2', logT2
1138 : write(*,1) 'logT3', logT3
1139 : write(*,1) 'logT4', logT4
1140 : write(*,1) 'logT5', logT5
1141 : write(*,1) 'logT6', logT6
1142 : write(*,1) 'logT7', logT7
1143 : write(*,1) 'logT8', logT8
1144 : write(*,'(A)')
1145 : write(*,1) 'logQ1', logQ1
1146 : write(*,1) 'logQ2', logQ2
1147 : write(*,'(A)')
1148 : call mesa_error(__FILE__,__LINE__,'eosdt_eval')
1149 : end if
1150 :
1151 : ! check validity of Rho's and T's for region boundaries
1152 : if (logRho1 <= logRho2 .or. logRho2 <= logRho3 .or. &
1153 : logRho3 <= logRho4 .or. logRho4 <= logRho5 .or. &
1154 : logRho5 <= logRho6 .or. logRho6 <= logRho7 .or. &
1155 : logRho7 <= logRho8 .or. &
1156 70 : logT1 <= logT2 .or. logT2 <= logT3 .or. logT3 <= logT4 .or. &
1157 : logT7 <= logT8) then
1158 0 : write(*,'(A)')
1159 0 : write(*,'(A)')
1160 0 : write(*,'(A)')
1161 0 : write(*,'(a)') 'must have strictly decreasing values for eos logT + logRho region boundaries'
1162 0 : write(*,'(A)')
1163 0 : write(*,1) 'logRho1', logRho1
1164 0 : write(*,1) 'logRho2', logRho2
1165 0 : write(*,1) 'logRho3', logRho3
1166 0 : write(*,1) 'logRho4', logRho4
1167 0 : write(*,1) 'logRho5', logRho5
1168 0 : write(*,1) 'logRho6', logRho6
1169 0 : write(*,1) 'logRho7', logRho7
1170 0 : write(*,1) 'logRho8', logRho8
1171 0 : write(*,'(A)')
1172 0 : write(*,1) 'logT1', logT1
1173 0 : write(*,1) 'logT2', logT2
1174 0 : write(*,1) 'logT3', logT3
1175 0 : write(*,1) 'logT4', logT4
1176 0 : write(*,1) 'logT5', logT5
1177 0 : write(*,1) 'logT6', logT6
1178 0 : write(*,1) 'logT7', logT7
1179 0 : write(*,1) 'logT8', logT8
1180 0 : write(*,'(A)')
1181 0 : write(*,1) 'logQ1', logQ1
1182 0 : write(*,1) 'logQ2', logQ2
1183 0 : write(*,'(A)')
1184 0 : write(*,'(A)')
1185 0 : if (logT3 <= logT4) then
1186 0 : write(*,'(a)') 'must have logRho1 > logRho2 + logQ1 - logQ2'
1187 0 : write(*,1) 'must have logRho1 > ', logRho2 + logQ1 - logQ2
1188 0 : write(*,1) 'logRho1', logRho1
1189 0 : write(*,'(a)') 'logRho1_OPAL_SCVH_limit sets logRho1'
1190 0 : write(*,'(a)') 'logRho2_OPAL_SCVH_limit sets logRho2'
1191 0 : write(*,1) 'max allowed logRho1 is', logRho1_max
1192 : end if
1193 0 : write(*,'(A)')
1194 0 : ierr = -1
1195 0 : return
1196 : end if
1197 :
1198 70 : call determine_region_opal_scvh
1199 :
1200 70 : call set_alfa_and_partials
1201 70 : if (ierr /= 0) return
1202 :
1203 : end if
1204 :
1205 :
1206 : contains
1207 :
1208 :
1209 70 : subroutine determine_region_opal_scvh
1210 : logical, parameter :: dbg = .false.
1211 70 : real(dp) :: logRho_hi, logRho_lo, d_logRho_dlogT, &
1212 70 : d_alfa_dlogQ, dlogQ_dlogRho, dlogQ_dlogT, Z_all_HELM
1213 :
1214 : include 'formats'
1215 :
1216 70 : logQ = logRho - 2d0*logT + 12d0
1217 70 : d_dx_dlogRho=0d0
1218 70 : d_dy_dlogT=0d0
1219 :
1220 : ! in high-Z fall back to HELM
1221 70 : if (Z >= rq% Z_all_HELM) then
1222 0 : iregion = use_none
1223 : if (dbg) then
1224 : write(*,1) 'iregion = use_none'
1225 : write(*,1) 'Z Z_all_HELM', Z, rq% Z_all_HELM
1226 : stop
1227 : end if
1228 : return
1229 : end if
1230 :
1231 : ! blends in T/Rho
1232 :
1233 : if (logT >= logT1 .or. logT <= logT8 .or. logRho >= logRho1 .or. &
1234 70 : logQ <= logQ4 .or. logQ >= logQmax .or. &
1235 : (logRho <= logRho6 .and. logT <= logT6)) then
1236 : if (dbg) then
1237 : write(*,*) 'logT >= logT1', logT >= logT1
1238 : write(*,*) 'logT <= logT8', logT <= logT8
1239 : write(*,*) 'logRho >= logRho1', logRho >= logRho1
1240 : write(*,*) 'logQ <= logQ4', logQ <= logQ4, logQ, logQ4
1241 : write(*,*) 'logQ >= logQmax', logQ >= logQmax
1242 : write(*,*) 'logRho <= logRho6 .and. logT <= logT6', logRho <= logRho6 .and. logT <= logT6
1243 : write(*,1) 'iregion = use_none 1 logT logT5 logT6', logT, logT5, logT6
1244 : end if
1245 0 : iregion = use_none
1246 :
1247 70 : else if (logQ <= logQ3 .and. logT >= logT5) then ! blend in Q
1248 0 : d_alfa_dlogQ = 1d0/(logQ4 - logQ3)
1249 0 : alfa = (logQ - logQ3)*d_alfa_dlogQ
1250 0 : dlogQ_dlogRho = 1d0
1251 0 : dlogQ_dlogT = -2d0
1252 0 : d_dx_dlogRho = d_alfa_dlogQ*dlogQ_dlogRho
1253 0 : d_dy_dlogT = d_alfa_dlogQ*dlogQ_dlogT
1254 0 : c_dx = alfa
1255 : if (dbg) then
1256 : write(*,*) 'iregion = blend_diagonal'
1257 : write(*,1) 'logQ3', logQ3
1258 : write(*,1) 'logQ', logQ
1259 : write(*,1) 'logQ4', logQ4
1260 : write(*,1) 'd_alfa_dlogQ', d_alfa_dlogQ
1261 : write(*,1) 'c_dx', c_dx
1262 : write(*,1) 'd_dx_dlogRho', d_dx_dlogRho
1263 : write(*,1) 'd_dy_dlogT', d_dy_dlogT
1264 : end if
1265 0 : iregion = blend_diagonal
1266 :
1267 70 : else if (logT >= logT2) then
1268 : if (dbg) write(*,*) 'logT >= logT2', logT, logT2
1269 0 : if (logT1 - logT2 < 0.01d0) then
1270 : d_dy_dlogT = 0d0
1271 : ! bad blend partials cause problems for 150M_z1m4_pre_ms_to_collapse
1272 : ! have tried to fix, but failed. hence this awful workaround.
1273 : else
1274 0 : d_dy_dlogT = 1/(logT1 - logT2)
1275 : end if
1276 0 : c_dy = (logT - logT2)*d_dy_dlogT
1277 0 : if (logRho > logRho2) then
1278 : if (dbg) write(*,*) 'logRho > logRho2', logRho, logRho2
1279 0 : d_dx_dlogRho = 1/(logRho1 - logRho2)
1280 0 : c_dx = (logRho - logRho2)*d_dx_dlogRho
1281 : if (dbg) write(*,*) 'iregion = blend_corner_out'
1282 0 : iregion = blend_corner_out
1283 : else ! logRho <= logRho2
1284 : if (dbg) write(*,*) 'logRho <= logRho2', logRho, logRho2
1285 : if (dbg) write(*,*) 'iregion = blend_in_y'
1286 0 : iregion = blend_in_y
1287 : end if
1288 :
1289 70 : else if (logT >= logT3) then ! NOTE: this assumes logT3 > logT4
1290 : if (dbg) write(*,*) 'logT >= logT3', logT, logT3
1291 8 : if (logRho > logRho2) then
1292 : if (dbg) write(*,*) 'logRho > logRho2', logRho, logRho2
1293 0 : d_dx_dlogRho = 1/(logRho1 - logRho2)
1294 0 : c_dx = (logRho - logRho2)*d_dx_dlogRho
1295 : if (dbg) write(*,*) 'iregion = blend_in_x'
1296 0 : iregion = blend_in_x
1297 : else
1298 : if (dbg) write(*,*) 'logRho <= logRho2', logRho, logRho2
1299 : if (dbg) write(*,*) 'iregion = use_all'
1300 8 : iregion = use_all
1301 : end if
1302 :
1303 62 : else if (logT >= logT4) then
1304 : if (dbg) write(*,*) 'logT >= logT4', logT, logT4
1305 0 : logRho_hi = logQ1 + 2*logT - 12
1306 0 : if (logRho >= logRho_hi) then
1307 : if (dbg) write(*,*) 'logRho >= logRho_hi', logRho, logRho_hi
1308 : if (dbg) write(*,*) 'iregion = use_none 2'
1309 0 : iregion = use_none
1310 0 : else if (logRho > logRho2) then
1311 : if (dbg) write(*,*) 'logRho > logRho2', logRho, logRho2
1312 0 : d_dx_dlogRho = 1/(logRho_hi - logRho2)
1313 0 : c_dx = (logRho - logRho2)*d_dx_dlogRho
1314 : if (dbg) write(*,*) 'iregion = blend_in_x'
1315 0 : iregion = blend_in_x
1316 : else ! logRho <= logRho2
1317 : if (dbg) write(*,*) 'logRho <= logRho2', logRho, logRho2
1318 : if (dbg) write(*,*) 'iregion = use_all'
1319 0 : iregion = use_all
1320 : end if
1321 :
1322 62 : else if (logRho > logRho4) then
1323 : if (dbg) write(*,*) 'logRho > logRho4', logRho, logRho4
1324 8 : if (logT > logT7) then
1325 8 : A = ((logQ1+2*logT4-12) - logRho3)/(logT4-logT7)
1326 8 : logRho_hi = logRho3 + (logT-logT7)*A
1327 8 : B = (logRho2-logRho4)/(logT4-logT7)
1328 8 : logRho_lo = logRho4 + (logT-logT7)*B
1329 8 : if (logRho >= logRho_hi) then
1330 : if (dbg) write(*,*) 'logRho >= logRho_hi', logRho, logRho_hi
1331 : if (dbg) write(*,*) 'iregion = use_none 3'
1332 0 : iregion = use_none
1333 8 : else if (logRho >= logRho_lo) then
1334 : if (dbg) write(*,*) 'logRho >= logRho_lo', logRho, logRho_lo
1335 0 : c_dx = (logRho - logRho_lo)/(logRho_hi - logRho_lo)
1336 0 : d_dx_dlogRho = 1/(logRho3 - logRho4 + (A - B)*(logT - logT7))
1337 : if (dbg) write(*,*) 'iregion = blend_in_x'
1338 0 : iregion = blend_in_x
1339 : else ! logRho < logRho_lo
1340 : if (dbg) write(*,*) 'logRho < logRho_lo', logRho, logRho_lo
1341 : if (dbg) write(*,*) 'iregion = use_all'
1342 8 : iregion = use_all
1343 : end if
1344 : else ! logT is > logT8
1345 : if (dbg) write(*,*) 'logT > logT8', logT, logT8
1346 0 : if (logRho > logRho3) then
1347 : if (dbg) write(*,*) 'logRho > logRho3', logRho, logRho3
1348 : if (dbg) write(*,*) 'iregion = use_none 4'
1349 0 : iregion = use_none
1350 : else ! logRho is > logRho4
1351 : if (dbg) write(*,*) 'logRho <= logRho3', logRho, logRho3
1352 0 : d_dx_dlogRho = 1/(logRho3 - logRho4)
1353 0 : c_dx = (logRho - logRho4)*d_dx_dlogRho
1354 0 : d_dy_dlogT = 1/(logT8 - logT7)
1355 0 : c_dy = (logT - logT7)*d_dy_dlogT
1356 : if (dbg) write(*,*) 'iregion = blend_corner_out'
1357 0 : iregion = blend_corner_out
1358 : end if
1359 : end if
1360 :
1361 54 : else if (logRho >= logRho5 .or. logT > logT5) then
1362 : if (dbg) write(*,*) 'iregion = use_all'
1363 54 : iregion = use_all
1364 :
1365 0 : else if (logT >= logT6) then
1366 0 : if (logRho <= logRho6) then
1367 0 : d_dy_dlogT = 1/(logT6 - logT5)
1368 0 : c_dy = (logT - logT5)*d_dy_dlogT
1369 : if (dbg) write(*,*) 'iregion = blend_in_y'
1370 0 : iregion = blend_in_y
1371 : else
1372 0 : d_dx_dlogRho = 1/(logRho5 - logRho6)
1373 0 : c_dx = (logRho - logRho6)*d_dx_dlogRho
1374 0 : d_dy_dlogT = 1/(logT5 - logT6)
1375 0 : c_dy = (logT - logT6)*d_dy_dlogT
1376 : if (dbg) write(*,*) 'iregion = blend_corner_in'
1377 0 : iregion = blend_corner_in
1378 : end if
1379 :
1380 : else
1381 : if (dbg) write(*,*) 'logRho > logRho6', logRho, logRho6
1382 0 : d_dx_dlogRho = 1/(logRho6 - logRho5)
1383 0 : c_dx = (logRho - logRho5)*d_dx_dlogRho
1384 : if (dbg) write(*,*) 'iregion = blend_in_x'
1385 0 : iregion = blend_in_x
1386 : end if
1387 :
1388 : if (dbg) call mesa_error(__FILE__,__LINE__,'determine_region')
1389 :
1390 : end subroutine determine_region_opal_scvh
1391 :
1392 :
1393 70 : subroutine set_alfa_and_partials ! alfa = fraction other
1394 : logical, parameter :: dbg = .false.
1395 :
1396 70 : real(dp) :: zfactor
1397 :
1398 : include 'formats'
1399 :
1400 70 : d_alfa_dlogT = 0d0
1401 70 : d_alfa_dlogRho = 0
1402 :
1403 70 : if (iregion == use_none .or. Z >= Z_all_HELM) then
1404 : if (dbg) write(*,*) 'iregion == use_none'
1405 0 : alfa = 1
1406 : else if (iregion == use_all) then
1407 : if (dbg) write(*,*) 'iregion == use_all'
1408 70 : alfa = 0
1409 : else if (iregion == blend_in_y) then
1410 : if (dbg) write(*,*) 'iregion == blend_in_y'
1411 0 : alfa = c_dy
1412 0 : d_alfa_dlogT = d_dy_dlogT
1413 : else if (iregion == blend_in_x) then
1414 : if (dbg) write(*,*) 'iregion == blend_in_x'
1415 0 : alfa = c_dx
1416 0 : d_alfa_dlogRho = d_dx_dlogRho
1417 : else if (iregion == blend_diagonal) then
1418 : if (dbg) write(*,*) 'iregion == blend_diagonal'
1419 0 : alfa = c_dx
1420 0 : d_alfa_dlogRho = d_dx_dlogRho
1421 0 : d_alfa_dlogT = d_dy_dlogT
1422 : else if (iregion == blend_corner_out) then
1423 : if (dbg) write(*,*) 'iregion == blend_corner_out'
1424 0 : alfa = sqrt(c_dx*c_dx + c_dy*c_dy)
1425 0 : if (alfa >= 1d0) then
1426 0 : alfa = 1
1427 0 : else if (alfa < 1d-10) then
1428 0 : alfa = 0
1429 : else
1430 0 : d_alfa_dlogT = c_dy*d_dy_dlogT/alfa
1431 0 : d_alfa_dlogRho = c_dx*d_dx_dlogRho/alfa
1432 : end if
1433 : else if (iregion == blend_corner_in) then
1434 : if (dbg) write(*,*) 'iregion == blend_corner_in'
1435 0 : beta = sqrt(c_dx*c_dx + c_dy*c_dy)
1436 0 : alfa = 1d0 - beta
1437 0 : if (alfa >= 1d0) then
1438 0 : alfa = 1d0
1439 0 : else if (alfa < 1d-10) then
1440 0 : alfa = 0
1441 : else
1442 0 : d_alfa_dlogT = -c_dy*d_dy_dlogT/beta
1443 0 : d_alfa_dlogRho = -c_dx*d_dx_dlogRho/beta
1444 : end if
1445 : else
1446 0 : ierr = -1
1447 0 : return
1448 : end if
1449 :
1450 70 : if (Z > Z_no_HELM .and. Z < Z_all_HELM .and. alfa < 1d0) then
1451 : ! reduce alfa to reduce the HELM fraction
1452 0 : zfactor = (Z - Z_no_HELM)/(Z_all_HELM - Z_no_HELM)
1453 0 : alfa = alfa*zfactor
1454 0 : d_alfa_dlogRho = d_alfa_dlogRho*zfactor
1455 0 : d_alfa_dlogT = d_alfa_dlogT*zfactor
1456 : end if
1457 :
1458 : end subroutine set_alfa_and_partials
1459 :
1460 :
1461 : end subroutine get_opal_scvh_alfa_and_partials
1462 :
1463 :
1464 1145158 : subroutine combine_for_eosdt( &
1465 : get_1st, get_2nd, remaining_fraction, &
1466 : alfa_in, d_alfa_dlogT_in, d_alfa_dlogRho_in, &
1467 : rq, dbg, Z, X, abar, zbar, &
1468 1145158 : species, chem_id, net_iso, xa, &
1469 : rho, logRho, T, logT, &
1470 1145158 : res, d_dlnd, d_dlnT, d_dxa, &
1471 : skip, ierr)
1472 : use eosdt_support, only : Do_Blend
1473 : procedure (get_values_for_eosdt_interface), pointer :: get_1st, get_2nd
1474 : type (EoS_General_Info), pointer :: rq
1475 : logical, intent(in) :: dbg
1476 : real(dp), intent(in) :: Z, X, abar, zbar, &
1477 : alfa_in, d_alfa_dlogT_in, d_alfa_dlogRho_in, remaining_fraction
1478 : integer, intent(in) :: species
1479 : integer, pointer :: chem_id(:), net_iso(:)
1480 : real(dp), intent(in) :: xa(:)
1481 : real(dp), intent(in) :: rho, logRho, T, logT
1482 : real(dp), intent(inout), dimension(nv) :: res, d_dlnd, d_dlnT
1483 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
1484 : logical, intent(out) :: skip
1485 : integer, intent(out) :: ierr
1486 :
1487 : real(dp), dimension(nv) :: &
1488 179789806 : res_1, d_dlnd_1, d_dlnT_1, res_2, d_dlnd_2, d_dlnT_2
1489 1145158 : real(dp), dimension(:,:), allocatable :: d_dxa_1, d_dxa_2
1490 : real(dp) :: alfa, d_alfa_dlogT, d_alfa_dlogRho
1491 : logical :: skip_1st, skip_2nd
1492 : logical, parameter :: linear_blend = .false.
1493 :
1494 : include 'formats'
1495 :
1496 1145158 : ierr = 0
1497 1145158 : skip = .false.
1498 :
1499 1145158 : allocate(d_dxa_1(nv, species), d_dxa_2(nv, species))
1500 :
1501 1145158 : alfa = alfa_in
1502 1145158 : d_alfa_dlogT = d_alfa_dlogT_in
1503 1145158 : d_alfa_dlogRho = d_alfa_dlogRho_in
1504 :
1505 1145158 : if (alfa == 0d0) then ! pure 1st
1506 : call get_1st(rq% handle, dbg, &
1507 : Z, X, abar, zbar, &
1508 : species, chem_id, net_iso, xa, &
1509 : rho, logRho, T, logT, remaining_fraction, &
1510 : res, d_dlnd, d_dlnT, d_dxa, &
1511 241252 : skip_1st, ierr)
1512 241252 : if (ierr /= 0) then
1513 2260436 : if (.not. rq% okay_to_convert_ierr_to_skip) return
1514 0 : if (dbg) write(*,*) 'ierr => skip 1st in combine_for_eosdt'
1515 0 : ierr = 0; skip_1st = .true.
1516 : end if
1517 241252 : if (skip_1st) then ! switch to pure 2nd
1518 0 : alfa = 1d0; d_alfa_dlogT = 0d0; d_alfa_dlogRho = 0d0
1519 : else
1520 : return
1521 : end if
1522 : end if
1523 :
1524 903906 : if (alfa < 1d0) then ! some of 1st
1525 : call get_1st(rq% handle, dbg, &
1526 : Z, X, abar, zbar, &
1527 : species, chem_id, net_iso, xa, &
1528 : rho, logRho, T, logT, remaining_fraction, &
1529 : res_2, d_dlnd_2, d_dlnT_2, d_dxa_2, &
1530 14940 : skip_1st, ierr)
1531 14940 : if (ierr /= 0) then
1532 0 : if (.not. rq% okay_to_convert_ierr_to_skip) return
1533 0 : if (dbg) write(*,*) 'ierr => skip 1st in combine_for_eosdt'
1534 0 : ierr = 0; skip_1st = .true.
1535 : end if
1536 14940 : if (skip_1st) then ! switch to pure 2nd
1537 0 : alfa = 1d0; d_alfa_dlogT = 0d0; d_alfa_dlogRho = 0d0
1538 : end if
1539 : end if
1540 :
1541 903906 : if (alfa == 1d0) then ! no 1st
1542 : call get_2nd(rq% handle, dbg, &
1543 : Z, X, abar, zbar, &
1544 : species, chem_id, net_iso, xa, &
1545 : Rho, logRho, T, logT, remaining_fraction, &
1546 : res, d_dlnd, d_dlnT, d_dxa, &
1547 888966 : skip_2nd, ierr)
1548 888966 : if (ierr /= 0) then
1549 0 : if (.not. rq% okay_to_convert_ierr_to_skip) return
1550 0 : if (dbg) write(*,*) 'ierr => skip 2nd in combine_for_eosdt'
1551 0 : ierr = 0; skip_2nd = .true.
1552 : end if
1553 888966 : if (skip_2nd) skip = .true.
1554 888966 : return
1555 : end if
1556 :
1557 : ! blend 1st and 2nd
1558 :
1559 : call get_2nd( &
1560 : rq% handle, dbg, Z, X, abar, zbar, &
1561 : species, chem_id, net_iso, xa, &
1562 : Rho, logRho, T, logT, remaining_fraction, &
1563 : res_1, d_dlnd_1, d_dlnT_1, d_dxa_1, &
1564 14940 : skip_2nd, ierr)
1565 14940 : if (ierr /= 0) then
1566 0 : if (.not. rq% okay_to_convert_ierr_to_skip) return
1567 0 : if (dbg) write(*,*) 'ierr => skip 2nd in combine_for_eosdt'
1568 0 : ierr = 0; skip_2nd = .true.
1569 : end if
1570 14940 : if (skip_2nd) then
1571 0 : skip = .true.
1572 0 : return
1573 : end if
1574 : call Do_Blend( &
1575 : rq, species, Rho, logRho, T, logT, &
1576 : alfa, d_alfa_dlogT, d_alfa_dlogRho, linear_blend, &
1577 : res_1, d_dlnd_1, d_dlnT_1, d_dxa_1, &
1578 : res_2, d_dlnd_2, d_dlnT_2, d_dxa_2, &
1579 14940 : res, d_dlnd, d_dlnT, d_dxa)
1580 :
1581 3405594 : end subroutine combine_for_eosdt
1582 :
1583 :
1584 210672 : subroutine get1_for_eosdt( &
1585 : handle, which_eosdt, dbg, Z, X, abar, zbar, &
1586 : species, chem_id, net_iso, xa, &
1587 : rho, logRho, T, logT, remaining_fraction, &
1588 210672 : res, d_dlnd, d_dlnT, d_dxa, &
1589 : skip, ierr)
1590 1145158 : use chem_def, only: chem_isos
1591 : integer, intent(in) :: handle
1592 : integer, intent(in) :: which_eosdt
1593 : logical, intent(in) :: dbg
1594 : real(dp), intent(in) :: &
1595 : Z, X, abar, zbar, remaining_fraction
1596 : integer, intent(in) :: species
1597 : integer, pointer :: chem_id(:), net_iso(:)
1598 : real(dp), intent(in) :: xa(:)
1599 : real(dp), intent(in) :: rho, logRho, T, logT
1600 : real(dp), intent(inout), dimension(nv) :: &
1601 : res, d_dlnd, d_dlnT
1602 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
1603 11376288 : real(dp), dimension(nv) :: d_dX, d_dZ
1604 : logical, intent(out) :: skip
1605 : integer, intent(out) :: ierr
1606 : type (EoS_General_Info), pointer :: rq
1607 : type (DT_xz_Info), pointer :: xz
1608 : integer :: i
1609 210672 : rq => eos_handles(handle)
1610 210672 : if (which_eosdt == eosdt_max_FreeEOS) then
1611 210602 : xz => FreeEOS_xz_struct
1612 : else
1613 70 : xz => eosDT_xz_struct
1614 : end if
1615 : call Get1_eosdt_Results( &
1616 : rq, which_eosdt, xz, Z, X, Rho, logRho, T, logT, &
1617 210672 : res, d_dlnd, d_dlnT, d_dX, d_dZ, ierr)
1618 1877236 : do i=1,species
1619 210672 : select case(chem_isos% Z(chem_id(i))) ! charge
1620 : case (1) ! X
1621 5688144 : d_dxa(:,i) = d_dX
1622 : case (2) ! Y
1623 10868202 : d_dxa(:,i) = 0
1624 : case default ! Z
1625 29054080 : d_dxa(:,i) = d_dZ
1626 : end select
1627 : end do
1628 210672 : skip = .false.
1629 210672 : end subroutine get1_for_eosdt
1630 :
1631 :
1632 210672 : subroutine Get1_eosdt_Results( & ! blend in Z
1633 : rq, which_eosdt, xz, Z, X, Rho, logRho, T, logT, &
1634 : res, dlnd, dlnT, dX, dZ, ierr)
1635 210672 : use chem_def
1636 : type (EoS_General_Info), pointer :: rq
1637 : integer, intent(in) :: which_eosdt
1638 : type (DT_xz_Info), pointer :: xz
1639 : real(dp), intent(in) :: Z, X, Rho, logRho, T, logT
1640 : real(dp), intent(inout), dimension(nv) :: res, dlnd, dlnT, dX, dZ
1641 : integer, intent(out) :: ierr
1642 :
1643 45715824 : real(dp), dimension(nv, 2) :: res_zx, dlnd_zx, dlnT_zx, dX_zx
1644 1053360 : real(dp) :: denom, c(2), dcdZ(2), tiny
1645 :
1646 : integer :: iz, j, ci
1647 :
1648 : include 'formats'
1649 :
1650 210672 : ierr = 0
1651 210672 : tiny = rq% tiny_fuzz
1652 :
1653 210672 : if (xz% nZs < 3) then
1654 0 : write(*, *) 'error: Get1_eosdt_Results assumes nZs >= 3'
1655 0 : call mesa_error(__FILE__,__LINE__)
1656 : end if
1657 :
1658 210672 : if (xz% Zs(1) /= 0) then
1659 0 : write(*, *) 'error: Get1_eosdt_Results assumes eos_Zs(1) == 0'
1660 0 : call mesa_error(__FILE__,__LINE__)
1661 : end if
1662 :
1663 210672 : if (abs(xz% Zs(1) - 2*xz% Zs(2) + xz% Zs(3)) > tiny) then
1664 0 : write(*, *) 'error: Get1_eosdt_Results assumes equal spaced Zs(1:3)'
1665 0 : call mesa_error(__FILE__,__LINE__)
1666 : end if
1667 :
1668 210672 : if (Z <= max(1d-20,xz% Zs(1))) then
1669 : call Get1_eosdt_for_X( &
1670 : rq, which_eosdt, xz, 1, X, &
1671 : Rho, logRho, T, logT, &
1672 12 : res, dlnd, dlnT, dX, ierr)
1673 12 : dZ = 0
1674 12 : return
1675 : end if
1676 :
1677 : ! zero these for now
1678 1895940 : res_zx(i_phase:i_latent_ddlnRho,:) = 0d0
1679 3581220 : res_zx(i_frac:i_frac+num_eos_frac_results-1,:) = 0d0
1680 :
1681 1895940 : dlnd_zx(i_phase:i_latent_ddlnRho,:) = 0d0
1682 3581220 : dlnd_zx(i_frac:i_frac+num_eos_frac_results-1,:) = 0d0
1683 :
1684 1895940 : dlnT_zx(i_phase:i_latent_ddlnRho,:) = 0d0
1685 3581220 : dlnT_zx(i_frac:i_frac+num_eos_frac_results-1,:) = 0d0
1686 :
1687 210660 : if (Z >= xz% Zs(xz% nZs)) then
1688 : call Get1_eosdt_for_X( &
1689 : rq, which_eosdt, xz, xz% nZs, X, &
1690 : Rho, logRho, T, logT, &
1691 0 : res, dlnd, dlnT, dX, ierr)
1692 0 : dZ = 0
1693 0 : return
1694 : end if
1695 :
1696 317068 : do iz = 2, xz% nZs
1697 317068 : if (Z < xz% Zs(iz)) then
1698 210660 : call do_interp2(iz-1,iz,ierr)
1699 210660 : if (ierr /= 0) return
1700 : exit
1701 : end if
1702 : end do
1703 :
1704 5898492 : do j=1,nv
1705 :
1706 5477160 : res(j) = c(1)*res_zx(j,1) + c(2)*res_zx(j,2)
1707 :
1708 : dlnd(j) = &
1709 5477160 : c(1)*dlnd_zx(j,1) + c(2)*dlnd_zx(j,2)
1710 :
1711 : dlnT(j) = &
1712 5477160 : c(1)*dlnT_zx(j,1) + c(2)*dlnT_zx(j,2)
1713 :
1714 : dX(j) = &
1715 5477160 : c(1)*dX_zx(j,1) + c(2)*dX_zx(j,2)
1716 :
1717 : dZ(j) = &
1718 5687820 : dcdZ(1)*res_zx(j,1) + dcdZ(2)*res_zx(j,2)
1719 :
1720 : end do
1721 :
1722 : contains
1723 :
1724 210660 : subroutine do_interp2(iz1, iz2, ierr)
1725 : integer, intent(in) :: iz1, iz2
1726 : integer, intent(out) :: ierr
1727 210660 : real(dp) :: Z1, Z2
1728 : include 'formats'
1729 : ierr = 0
1730 210660 : Z1 = xz% Zs(iz1)
1731 210660 : Z2 = xz% Zs(iz2)
1732 210660 : c(2) = (Z - Z1) / (Z2 - Z1)
1733 210660 : c(1) = 1d0 - c(2)
1734 210660 : dcdZ(2) = 1d0/(Z2 - Z1)
1735 210660 : dcdZ(1) = -dcdZ(2)
1736 : call Get1_eosdt_for_X( &
1737 : rq, which_eosdt, xz, iz1, X, Rho, logRho, T, logT, &
1738 : res_zx(:,1), dlnd_zx(:,1), dlnT_zx(:,1), dX_zx(:,1), &
1739 210660 : ierr)
1740 210660 : if (ierr /= 0) return
1741 : call Get1_eosdt_for_X( &
1742 : rq, which_eosdt, xz, iz2, X, Rho, logRho, T, logT, &
1743 : res_zx(:,2), dlnd_zx(:,2), dlnT_zx(:,2), dX_zx(:,2), &
1744 210660 : ierr)
1745 : if (ierr /= 0) return
1746 210672 : end subroutine do_interp2
1747 :
1748 : end subroutine Get1_eosdt_Results
1749 :
1750 :
1751 421332 : subroutine Get1_eosdt_for_X( &
1752 : rq, which_eosdt, xz, iz, X, Rho, logRho, T, logT, &
1753 : res, dlnd, dlnT, d_dX, ierr)
1754 : type (EoS_General_Info), pointer :: rq
1755 : integer, intent(in) :: which_eosdt
1756 : type (DT_xz_Info), pointer :: xz
1757 : integer, intent(in) :: iz ! the index in eos_Zs
1758 : real(dp), intent(in) :: X, Rho, logRho, T, logT
1759 : real(dp), intent(inout), dimension(nv) :: &
1760 : res, dlnd, dlnT, d_dX
1761 : integer, intent(out) :: ierr
1762 :
1763 : real(dp), dimension(nv, 4) :: &
1764 137354232 : res_zx, dlnd_zx, dlnT_zx
1765 2527992 : real(dp) :: dX, dX1, dX2, dX3, c(4), dcdX(4), denom, delX, coef, dcoef_dX, alfa, beta, dalfa_dX, dbeta_dX, tiny
1766 : character (len=256) :: message
1767 : integer :: ix, ix_lo, ix_hi, j, num_Xs
1768 : logical, parameter :: dbg_for_X = dbg ! .or. .true.
1769 : logical :: what_we_use_is_equal_spaced
1770 :
1771 : include 'formats'
1772 :
1773 421332 : ierr = 0
1774 421332 : tiny = rq% tiny_fuzz
1775 :
1776 421332 : num_Xs = xz% nXs_for_Z(iz)
1777 :
1778 421332 : if (xz% Xs_for_Z(1,iz) /= 0d0) then
1779 0 : write(*, *) 'error: Get1_eosdt_for_X assumes xz% nXs_for_Z(1) == 0'
1780 0 : call mesa_error(__FILE__,__LINE__)
1781 : end if
1782 :
1783 421332 : if (X < tiny .or. num_Xs == 1) then
1784 : call Get1_eosdt_XTable_Results( &
1785 : rq, which_eosdt, 1, iz, Rho, logRho, T, logT, &
1786 10 : res, dlnd, dlnT, ierr)
1787 10 : d_dX = 0
1788 16 : return
1789 : end if
1790 :
1791 421322 : if (X >= xz% Xs_for_Z(num_Xs,iz)) then
1792 :
1793 : call Get1_eosdt_XTable_Results( &
1794 : rq, which_eosdt, num_Xs, iz, Rho, logRho, T, logT, &
1795 6 : res, dlnd, dlnT, ierr)
1796 6 : d_dX = 0
1797 :
1798 6 : if (is_bad(res(i_lnS))) then
1799 0 : ierr = -1
1800 0 : if (.not. stop_for_is_bad) return
1801 : write(*,1) 'res(i_lnS), logRho, logT', res(i_lnS), logRho, logT
1802 : call mesa_error(__FILE__,__LINE__,'Get1_eosdt_for_X num_Xs')
1803 : end if
1804 :
1805 : return
1806 : end if
1807 :
1808 421316 : if (rq% eosDT_use_linear_interp_for_X .or. num_Xs == 2) then
1809 0 : call do_linear
1810 0 : return
1811 : end if
1812 :
1813 421316 : ix_hi = -1
1814 421316 : if (X <= xz% Xs_for_Z(2,iz)) then
1815 0 : ix_lo = 1; ix_hi = 3
1816 421316 : else if (X >= xz% Xs_for_Z(num_Xs-1,iz)) then
1817 76 : ix_lo = num_Xs-2; ix_hi = num_Xs
1818 : else
1819 2527300 : do ix = 3, num_Xs-1
1820 2527300 : if (X <= xz% Xs_for_Z(ix,iz)) then
1821 421240 : ix_lo = ix-2; ix_hi = ix+1; exit
1822 : end if
1823 : end do
1824 : end if
1825 :
1826 421316 : if (ix_hi < 0) then
1827 0 : write(*, *) 'X', X
1828 0 : write(*, *) 'ix_lo', ix_lo
1829 0 : write(*, *) 'ix_hi', ix_hi
1830 0 : write(*, *) 'error: Get1_eosdt_for_X logic bug'
1831 0 : call mesa_error(__FILE__,__LINE__)
1832 : end if
1833 :
1834 : if (dbg_for_X) then
1835 : write(*, *) 'X', X
1836 : write(*, *) 'ix_lo', ix_lo
1837 : write(*, *) 'ix_hi', ix_hi
1838 : end if
1839 :
1840 421316 : what_we_use_is_equal_spaced = .true.
1841 421316 : dX1 = xz% Xs_for_Z(ix_lo+1,iz)-xz% Xs_for_Z(ix_lo,iz)
1842 421316 : dX2 = xz% Xs_for_Z(ix_lo+2,iz)-xz% Xs_for_Z(ix_lo+1,iz)
1843 421316 : if (ix_hi-ix_lo==2) then ! check that the 3 table X's are equal spaced
1844 76 : if (abs(dX1 - dX2) > tiny) what_we_use_is_equal_spaced = .false.
1845 : else ! check that the 4 table X's are equal spaced
1846 421240 : dX3 = xz% Xs_for_Z(ix_hi,iz)-xz% Xs_for_Z(ix_lo+2,iz)
1847 421240 : if (abs(dX1 - dX2) > tiny .or. abs(dX2 - dX3) > tiny) &
1848 : what_we_use_is_equal_spaced = .false.
1849 : end if
1850 :
1851 : if (.not. what_we_use_is_equal_spaced) then
1852 0 : call do_linear
1853 0 : if (is_bad(d_dX(1))) then
1854 0 : call mesa_error(__FILE__,__LINE__,'Get1_eosdt_for_X bad d_dX; linear')
1855 : end if
1856 0 : return
1857 : end if
1858 :
1859 2106504 : do ix=ix_lo, ix_hi
1860 1685188 : j = ix-ix_lo+1
1861 : call Get1_eosdt_XTable_Results( &
1862 : rq, which_eosdt, ix, iz, Rho, logRho, T, logT, &
1863 : res_zx(:, j), dlnd_zx(:, j), dlnT_zx(:, j), &
1864 1685188 : ierr)
1865 2106504 : if (ierr /= 0) return
1866 : end do
1867 :
1868 : ! zero these for now
1869 7162372 : res_zx(i_phase:i_latent_ddlnRho,:) = 0d0
1870 13903428 : res_zx(i_frac:i_frac+num_eos_frac_results-1,:) = 0d0
1871 :
1872 7162372 : dlnd_zx(i_phase:i_latent_ddlnRho,:) = 0d0
1873 13903428 : dlnd_zx(i_frac:i_frac+num_eos_frac_results-1,:) = 0d0
1874 :
1875 7162372 : dlnT_zx(i_phase:i_latent_ddlnRho,:) = 0d0
1876 13903428 : dlnT_zx(i_frac:i_frac+num_eos_frac_results-1,:) = 0d0
1877 :
1878 :
1879 421316 : delX = X - xz% Xs_for_Z(ix_lo,iz)
1880 421316 : dX = dX1
1881 :
1882 421316 : if (ix_hi-ix_lo==2) then
1883 :
1884 76 : denom = 2*dX*dX
1885 76 : c(1) = (2*dX*dX - 3*dX*delX + delX*delX)/denom
1886 76 : c(2) = 2*(2*dX-delX)*delX/denom
1887 76 : c(3) = delX*(delX-dX)/denom
1888 2052 : res(:) = c(1)*res_zx(:, 1) + c(2)*res_zx(:, 2) + c(3)*res_zx(:, 3)
1889 :
1890 : dlnd(:) = &
1891 : c(1)*dlnd_zx(:,1) + &
1892 : c(2)*dlnd_zx(:,2) + &
1893 2052 : c(3)*dlnd_zx(:,3)
1894 : dlnT(:) = &
1895 : c(1)*dlnT_zx(:,1) + &
1896 : c(2)*dlnT_zx(:,2) + &
1897 2052 : c(3)*dlnT_zx(:,3)
1898 :
1899 76 : dcdx(1) = (-3*dX + 2*delX)/denom
1900 76 : dcdx(2) = 2*(2*dX-2*delX)/denom
1901 76 : dcdx(3) = (2*delX-dX)/denom
1902 :
1903 : d_dX(:) = &
1904 : dcdX(1)*res_zx(:,1) + &
1905 : dcdX(2)*res_zx(:,2) + &
1906 2052 : dcdX(3)*res_zx(:,3)
1907 :
1908 76 : if (is_bad(d_dX(1))) then
1909 0 : call mesa_error(__FILE__,__LINE__,'Get1_eosdt_for_X bad d_dX; 3')
1910 : end if
1911 :
1912 : else
1913 :
1914 421240 : coef = (X-xz% Xs_for_Z(ix_lo+1,iz))/dX
1915 : ! coef = fractional location of X between 2nd and 3rd X's for fit.
1916 : ! coef is weight for the quadratic based on points 2, 3, 4 of fit.
1917 : ! (1-coef) is weight for quadratic based on points 1, 2, 3 of fit.
1918 421240 : coef = min(1d0,max(0d0,coef))
1919 421240 : c(1) = -coef*(coef-1)*(coef-1)/2
1920 421240 : c(2) = (2 - coef*coef*(5 - 3*coef))/2
1921 421240 : c(3) = coef*(1 + coef*(4 - 3*coef))/2
1922 421240 : c(4) = coef*coef*(coef-1)/2
1923 : res(:) = c(1)*res_zx(:, 1) + &
1924 : (c(2)*res_zx(:, 2) + &
1925 : (c(3)*res_zx(:, 3) + &
1926 11373480 : c(4)*res_zx(:, 4)))
1927 :
1928 : dlnd(:) = &
1929 : c(1)*dlnd_zx(:, 1) + &
1930 : (c(2)*dlnd_zx(:, 2) + &
1931 : (c(3)*dlnd_zx(:, 3) + &
1932 11373480 : c(4)*dlnd_zx(:, 4)))
1933 : dlnT(:) = &
1934 : c(1)*dlnT_zx(:, 1) + &
1935 : (c(2)*dlnT_zx(:, 2) + &
1936 : (c(3)*dlnT_zx(:, 3) + &
1937 11373480 : c(4)*dlnT_zx(:, 4)))
1938 :
1939 421240 : dcoef_dX = 1d0/dX
1940 : dcdX = 0
1941 421240 : dcdX(1) = -(3*coef*coef-4*coef+1)/2*dcoef_dX
1942 421240 : dcdX(2) = (9*coef*coef-10*coef)/2*dcoef_dX
1943 421240 : dcdX(3) = -(9*coef*coef-8*coef-1)/2*dcoef_dX
1944 421240 : dcdX(4) = coef*(3*coef-2)/2*dcoef_dX
1945 :
1946 : d_dX(:) = &
1947 : dcdX(1)*res_zx(:,1) + &
1948 : dcdX(2)*res_zx(:,2) + &
1949 : dcdX(3)*res_zx(:,3) + &
1950 11373480 : dcdX(4)*res_zx(:,4)
1951 :
1952 421240 : if (is_bad(d_dX(1))) then
1953 0 : call mesa_error(__FILE__,__LINE__,'Get1_eosdt_for_X bad d_dX; 4')
1954 : end if
1955 :
1956 : end if
1957 :
1958 : contains
1959 :
1960 0 : subroutine do_linear
1961 :
1962 0 : do ix = 2, num_Xs
1963 0 : if (xz% Xs_for_Z(ix,iz) >= X) exit
1964 : end do
1965 :
1966 0 : j = 1
1967 : call Get1_eosdt_XTable_Results( &
1968 : rq, which_eosdt, ix-1, iz, Rho, logRho, T, logT, &
1969 : res_zx(:,j), dlnd_zx(:,j), dlnT_zx(:,j), &
1970 0 : ierr)
1971 0 : if (ierr /= 0) then
1972 : if (.not. stop_for_is_bad) return
1973 : call mesa_error(__FILE__,__LINE__,'Get1_eosdt_for_X')
1974 : end if
1975 :
1976 0 : j = 2
1977 : call Get1_eosdt_XTable_Results( &
1978 : rq, which_eosdt, ix, iz, Rho, logRho, T, logT, &
1979 : res_zx(:,j), dlnd_zx(:,j), dlnT_zx(:,j), &
1980 0 : ierr)
1981 0 : if (ierr /= 0) then
1982 : if (.not. stop_for_is_bad) return
1983 : call mesa_error(__FILE__,__LINE__,'Get1_eosdt_for_X')
1984 : end if
1985 :
1986 : ! zero these for now
1987 0 : res_zx(i_phase:i_latent_ddlnRho,:) = 0d0
1988 0 : res_zx(i_frac:i_frac+num_eos_frac_results-1,:) = 0d0
1989 :
1990 0 : dlnd_zx(i_phase:i_latent_ddlnRho,:) = 0d0
1991 0 : dlnd_zx(i_frac:i_frac+num_eos_frac_results-1,:) = 0d0
1992 :
1993 0 : dlnT_zx(i_phase:i_latent_ddlnRho,:) = 0d0
1994 0 : dlnT_zx(i_frac:i_frac+num_eos_frac_results-1,:) = 0d0
1995 :
1996 :
1997 0 : alfa = (X - xz% Xs_for_Z(ix,iz))/(xz% Xs_for_Z(ix-1,iz) - xz% Xs_for_Z(ix,iz))
1998 0 : beta = 1d0 - alfa
1999 :
2000 0 : dalfa_dX = 1d0 / (xz% Xs_for_Z(ix-1,iz) - xz% Xs_for_Z(ix,iz))
2001 0 : dbeta_dX = -dalfa_dX
2002 :
2003 0 : do j=1,nv
2004 :
2005 0 : res(j) = alfa*res_zx(j,1) + beta*res_zx(j,2)
2006 :
2007 : dlnd(j) = &
2008 0 : alfa*dlnd_zx(j,1) + beta*dlnd_zx(j,2)
2009 :
2010 : dlnT(j) = &
2011 0 : alfa*dlnT_zx(j,1) + beta*dlnT_zx(j,2)
2012 :
2013 : d_dX(j) = &
2014 0 : dalfa_dX*res_zx(j,1) + dbeta_dX*res_zx(j,2)
2015 :
2016 : end do
2017 :
2018 : end subroutine do_linear
2019 :
2020 : end subroutine Get1_eosdt_for_X
2021 :
2022 :
2023 1685204 : subroutine Locate_logQ(rq, ep, logQ, iQ, logQ0, logQ1, ierr)
2024 : type (EoS_General_Info), pointer :: rq
2025 : type (EosDT_xz_Info), pointer :: ep
2026 : real(dp), intent(inout) :: logQ
2027 : integer, intent(out) :: iQ
2028 : real(dp), intent(out) :: logQ0, logQ1
2029 : integer, intent(out) :: ierr
2030 1685204 : ierr = 0
2031 1685204 : iQ = int((logQ - ep% logQ_min)/ep% del_logQ + 1d-4) + 1
2032 1685204 : if (iQ < 1 .or. iQ >= ep% num_logQs) then
2033 0 : if (iQ < 1) then
2034 0 : iQ = 1
2035 0 : logQ0 = ep% logQ_min
2036 0 : logQ1 = logQ0 + ep% del_logQ
2037 0 : logQ = logQ0
2038 0 : if (return_ierr_beyond_table_bounds) ierr = -1
2039 : else
2040 0 : iQ = ep% num_logQs-1
2041 0 : logQ0 = ep% logQ_min + (iQ-1)*ep% del_logQ
2042 0 : logQ1 = logQ0 + ep% del_logQ
2043 0 : logQ = logQ1
2044 0 : if (return_ierr_beyond_table_bounds) ierr = -1
2045 : end if
2046 : else
2047 1685204 : logQ0 = ep% logQ_min + (iQ-1)*ep% del_logQ
2048 1685204 : logQ1 = logQ0 + ep% del_logQ
2049 : end if
2050 1685204 : end subroutine Locate_logQ
2051 :
2052 :
2053 1685204 : subroutine Locate_logT(rq, ep, logT, iT, logT0, logT1, ierr)
2054 : type (EoS_General_Info), pointer :: rq
2055 : type (EosDT_xz_Info), pointer :: ep
2056 : real(dp), intent(inout) :: logT
2057 : integer, intent(out) :: iT
2058 : real(dp), intent(out) :: logT0, logT1
2059 : integer, intent(out) :: ierr
2060 1685204 : ierr = 0
2061 1685204 : iT = int((logT - ep% logT_min)/ep% del_logT + 1d-4) + 1
2062 1685204 : if (iT < 1 .or. iT >= ep% num_logTs) then
2063 0 : if (iT < 1) then
2064 0 : iT = 1
2065 0 : logT0 = ep% logT_min
2066 0 : logT1 = logT0 + ep% del_logT
2067 0 : logT = logT0
2068 0 : if (return_ierr_beyond_table_bounds) ierr = -1
2069 : else
2070 0 : iT = ep% num_logTs-1
2071 0 : logT0 = ep% logT_min + (iT-1)*ep% del_logT
2072 0 : logT1 = logT0 + ep% del_logT
2073 0 : logT = logT1
2074 0 : if (return_ierr_beyond_table_bounds) ierr = -1
2075 : end if
2076 : else
2077 1685204 : logT0 = ep% logT_min + (iT-1)*ep% del_logT
2078 1685204 : logT1 = logT0 + ep% del_logT
2079 : end if
2080 1685204 : end subroutine Locate_logT
2081 :
2082 :
2083 1685204 : subroutine Get1_eosdt_XTable_Results( &
2084 : rq, which_eosdt, ix, iz, Rho, logRho_in, T, logT_in, &
2085 : res, d_dlnd, d_dlnT, ierr)
2086 : use eosdt_support, only: Do_EoS_Interpolations
2087 : type (EoS_General_Info), pointer :: rq
2088 : integer, intent(in) :: which_eosdt
2089 : integer, intent(in) :: ix, iz
2090 : real(dp), intent(in) :: Rho, logRho_in, T, logT_in
2091 : real(dp), intent(inout), dimension(nv) :: res, d_dlnd, d_dlnT
2092 : integer, intent(out) :: ierr
2093 :
2094 : real(dp), parameter :: ln10sq = ln10*ln10
2095 : real(dp) :: &
2096 133131116 : fval(nv), df_dx(nv), df_dy(nv), &
2097 89315812 : df_dlnd(nv), df_dlnT(nv), &
2098 : energy, entropy, P, Pgas, Prad, x, y, &
2099 : dx_dlnd, dx_dlnT, dy_dlnd, dy_dlnT, &
2100 : chiT, chiRho, Cv, gamma1, numeric, &
2101 : dS_dlnd, dS_dlnT, dE_dlnd, dE_dlnT, &
2102 : dPgas_dlnd, dPgas_dlnT, dP_dlnd, dP_dlnT
2103 1685204 : real(dp) :: logQ0, logQ1, logT0, logT1, logRho0, logRho1
2104 : integer :: iQ, jtemp, k, j, irho
2105 : type (EosDT_xz_Info), pointer :: ep
2106 : logical, parameter :: show = .false.
2107 1685204 : real(dp) :: logRho, logT, logQ
2108 :
2109 : include 'formats'
2110 :
2111 1685204 : logRho = logRho_in
2112 1685204 : logT = logT_in
2113 1685204 : logQ = logRho - 2*logT + 12
2114 :
2115 : ierr = 0
2116 1685204 : call load_single_eosDT_table_by_id(rq, which_eosdt, ep, ix, iz, ierr)
2117 1685204 : if (ierr /= 0) return
2118 :
2119 1685204 : call Locate_logQ(rq, ep, logQ, iQ, logQ0, logQ1, ierr)
2120 1685204 : if (ierr /= 0) then
2121 0 : write(*,1) 'eosDT failed in Locate_logQ', logQ
2122 0 : return
2123 : end if
2124 :
2125 1685204 : call Locate_logT(rq, ep, logT, jtemp, logT0, logT1, ierr)
2126 1685204 : if (ierr /= 0) then
2127 0 : write(*,1) 'eosDT failed in Locate_logT', logT
2128 0 : return
2129 : end if
2130 :
2131 : call Do_EoS_Interpolations( &
2132 : 1, nv, nv, ep% num_logQs, ep% logQs, ep% num_logTs, ep% logTs, &
2133 : ep% tbl1, iQ, jtemp, logQ0, logQ, logQ1, logT0, logT, logT1, &
2134 1685204 : fval, df_dx, df_dy, ierr)
2135 1685204 : if (ierr /= 0) then
2136 0 : write(*,1) 'failed in Do_EoS_Interpolations'
2137 0 : return
2138 : end if
2139 :
2140 1685204 : if (is_bad(fval(i_lnS))) then
2141 0 : ierr = -1
2142 0 : if (.not. stop_for_is_bad) return
2143 : write(*,1) 'fval(i_lnS), logRho, logT', fval(i_lnS), logRho, logT
2144 : call mesa_error(__FILE__,__LINE__,'after Do_Interp_with_2nd_derivs')
2145 : end if
2146 :
2147 1685204 : res(i_lnPgas) = fval(i_lnPgas)
2148 1685204 : res(i_lnE) = fval(i_lnE)
2149 1685204 : res(i_lnS) = fval(i_lnS)
2150 :
2151 1685204 : if (is_bad(res(i_lnS))) then
2152 0 : ierr = -1
2153 0 : if (.not. stop_for_is_bad) return
2154 : write(*,1) 'res(i_lnS), logRho, logT', res(i_lnS), logRho, logT
2155 : call mesa_error(__FILE__,__LINE__,'after interpolation')
2156 : end if
2157 :
2158 1685204 : if (is_bad(res(i_lnS)) .or. res(i_lnS) > ln10*100) then
2159 0 : ierr = -1
2160 0 : if (.not. stop_for_is_bad) return
2161 : write(*,1) 'res(i_lnS), logRho, logT', res(i_lnS), logRho, logT
2162 : call mesa_error(__FILE__,__LINE__,'after interpolation')
2163 : end if
2164 :
2165 1685204 : res(i_grad_ad) = fval(i_grad_ad)
2166 1685204 : res(i_chiRho) = fval(i_chiRho)
2167 1685204 : res(i_chiT) = fval(i_chiT)
2168 :
2169 1685204 : res(i_Cp) = fval(i_Cp)
2170 1685204 : res(i_Cv) = fval(i_Cv)
2171 :
2172 1685204 : res(i_dE_dRho) = fval(i_dE_dRho)
2173 1685204 : res(i_dS_dT) = fval(i_dS_dT)
2174 1685204 : res(i_dS_dRho) = fval(i_dS_dRho)
2175 :
2176 1685204 : res(i_mu) = fval(i_mu)
2177 1685204 : res(i_lnfree_e) = fval(i_lnfree_e)
2178 1685204 : res(i_gamma1) = fval(i_gamma1)
2179 1685204 : res(i_gamma3) = fval(i_gamma3)
2180 1685204 : res(i_eta) = fval(i_eta)
2181 :
2182 : ! convert df_dx and df_dy to df_dlogRho_c_T and df_dlogT_c_Rho
2183 :
2184 : ! df_dx is df_dlogQ at const T
2185 : ! df_dy is df_dlogT_c_Rho at const Q
2186 : ! logQ = logRho - 2*logT + 12
2187 :
2188 : ! f = f(logQ(logRho,logT),logT)
2189 : ! df/dlogRho|T = df/dlogQ|T * dlogQ/dlogRho|T = df_dx
2190 : ! df/dlogT|Rho = df/dlogT|Q + df/dlogQ|T * dlogQ/dlogT|Rho = df_dy - 2*df_dx
2191 :
2192 45500508 : do k=1,nv
2193 43815304 : df_dlnd(k) = df_dx(k)/ln10
2194 45500508 : df_dlnT(k) = df_dy(k)/ln10 - 2d0*df_dlnd(k)
2195 : end do
2196 :
2197 1685204 : d_dlnd(i_lnPgas) = df_dlnd(i_lnPgas)
2198 1685204 : d_dlnd(i_lnE) = df_dlnd(i_lnE)
2199 1685204 : d_dlnd(i_lnS) = df_dlnd(i_lnS)
2200 1685204 : d_dlnd(i_grad_ad) = df_dlnd(i_grad_ad)
2201 1685204 : d_dlnd(i_chiRho) = df_dlnd(i_chiRho)
2202 1685204 : d_dlnd(i_chiT) = df_dlnd(i_chiT)
2203 :
2204 1685204 : d_dlnd(i_Cp) = df_dlnd(i_Cp)
2205 1685204 : d_dlnd(i_Cv) = df_dlnd(i_Cv)
2206 1685204 : d_dlnd(i_dE_dRho) = df_dlnd(i_dE_dRho)
2207 1685204 : d_dlnd(i_dS_dT) = df_dlnd(i_dS_dT)
2208 1685204 : d_dlnd(i_dS_dRho) = df_dlnd(i_dS_dRho)
2209 1685204 : d_dlnd(i_mu) = df_dlnd(i_mu)
2210 1685204 : d_dlnd(i_lnfree_e) = df_dlnd(i_lnfree_e)
2211 1685204 : d_dlnd(i_gamma1) = df_dlnd(i_gamma1)
2212 1685204 : d_dlnd(i_gamma3) = df_dlnd(i_gamma3)
2213 1685204 : d_dlnd(i_eta) = df_dlnd(i_eta)
2214 :
2215 1685204 : d_dlnT(i_lnPgas) = df_dlnT(i_lnPgas)
2216 1685204 : d_dlnT(i_lnE) = df_dlnT(i_lnE)
2217 1685204 : d_dlnT(i_lnS) = df_dlnT(i_lnS)
2218 1685204 : d_dlnT(i_grad_ad) = df_dlnT(i_grad_ad)
2219 1685204 : d_dlnT(i_chiRho) = df_dlnT(i_chiRho)
2220 1685204 : d_dlnT(i_chiT) = df_dlnT(i_chiT)
2221 1685204 : d_dlnT(i_Cp) = df_dlnT(i_Cp)
2222 1685204 : d_dlnT(i_Cv) = df_dlnT(i_Cv)
2223 1685204 : d_dlnT(i_dE_dRho) = df_dlnT(i_dE_dRho)
2224 1685204 : d_dlnT(i_dS_dT) = df_dlnT(i_dS_dT)
2225 1685204 : d_dlnT(i_dS_dRho) = df_dlnT(i_dS_dRho)
2226 1685204 : d_dlnT(i_mu) = df_dlnT(i_mu)
2227 1685204 : d_dlnT(i_lnfree_e) = df_dlnT(i_lnfree_e)
2228 1685204 : d_dlnT(i_gamma1) = df_dlnT(i_gamma1)
2229 1685204 : d_dlnT(i_gamma3) = df_dlnT(i_gamma3)
2230 1685204 : d_dlnT(i_eta) = df_dlnT(i_eta)
2231 :
2232 1685204 : if (is_bad(d_dlnd(i_lnS)) .or. is_bad(d_dlnT(i_lnS))) then
2233 0 : ierr = -1
2234 0 : if (.not. stop_for_is_bad) return
2235 : write(*,1) 'fval(i_lnS)', fval(i_lnS)
2236 : write(*,1) 'd_dlnd(i_lnS)', d_dlnd(i_lnS)
2237 : write(*,1) 'd_dlnT(i_lnS)', d_dlnT(i_lnS)
2238 : call mesa_error(__FILE__,__LINE__,'Get1_eosdt_XTable_Results')
2239 : end if
2240 :
2241 1685204 : end subroutine Get1_eosdt_XTable_Results
2242 :
2243 :
2244 1476 : subroutine get_T( &
2245 : handle, Z, X, abar, zbar, &
2246 1476 : species, chem_id, net_iso, xa, &
2247 : logRho, which_other, other_value, &
2248 : logT_tol, other_tol, max_iter, logT_guess, &
2249 : logT_bnd1, logT_bnd2, other_at_bnd1, other_at_bnd2, &
2250 1476 : logT_result, res, d_dlnRho_c_T, d_dlnT_c_Rho, d_dxa_c_TRho, &
2251 : eos_calls, ierr)
2252 :
2253 : integer, intent(in) :: handle
2254 :
2255 : real(dp), intent(in) :: Z ! the metals mass fraction
2256 : real(dp), intent(in) :: X ! the hydrogen mass fraction
2257 :
2258 : real(dp), intent(in) :: abar, zbar
2259 :
2260 : integer, intent(in) :: species
2261 : integer, pointer :: chem_id(:)
2262 : integer, pointer :: net_iso(:)
2263 : real(dp), intent(in) :: xa(:)
2264 :
2265 : real(dp), intent(in) :: logRho ! log10 of density
2266 : integer, intent(in) :: which_other ! from eos_def. e.g., i_P for pressure
2267 : real(dp), intent(in) :: other_value ! desired value for the other variable
2268 : real(dp), intent(in) :: other_tol
2269 :
2270 : real(dp), intent(in) :: logT_tol
2271 : integer, intent(in) :: max_iter ! max number of iterations
2272 :
2273 : real(dp), intent(in) :: logT_guess
2274 : real(dp), intent(in) :: logT_bnd1, logT_bnd2 ! bounds for logT
2275 : ! set to arg_not_provided if do not know bounds
2276 : real(dp), intent(in) :: other_at_bnd1, other_at_bnd2 ! values at bounds
2277 : ! if don't know these values, just set to arg_not_provided (defined in c_def)
2278 :
2279 : real(dp), intent(out) :: logT_result
2280 : real(dp), intent(inout), dimension(nv) :: res, d_dlnRho_c_T, d_dlnT_c_Rho
2281 : real(dp), intent(inout), dimension(:,:) :: d_dxa_c_TRho
2282 :
2283 : integer, intent(out) :: eos_calls
2284 : integer, intent(out) :: ierr ! 0 means AOK.
2285 :
2286 : logical, parameter :: doing_Rho = .false.
2287 :
2288 : call do_safe_get_Rho_T( &
2289 : handle, Z, X, abar, zbar, &
2290 : species, chem_id, net_iso, xa, &
2291 : logRho, which_other, other_value, doing_Rho, &
2292 : logT_guess, logT_result, logT_bnd1, logT_bnd2, other_at_bnd1, other_at_bnd2, &
2293 : logT_tol, other_tol, max_iter, res, d_dlnRho_c_T, d_dlnT_c_Rho, d_dxa_c_TRho, &
2294 1476 : eos_calls, ierr)
2295 :
2296 1685204 : end subroutine get_T
2297 :
2298 :
2299 12152 : subroutine get_Rho( &
2300 : handle, Z, X, abar, zbar, &
2301 12152 : species, chem_id, net_iso, xa, &
2302 : logT, which_other, other_value, &
2303 : logRho_tol, other_tol, max_iter, logRho_guess, &
2304 : logRho_bnd1, logRho_bnd2, other_at_bnd1, other_at_bnd2, &
2305 12152 : logRho_result, res, d_dlnRho_c_T, d_dlnT_c_Rho, d_dxa_c_TRho, &
2306 : eos_calls, ierr)
2307 :
2308 : use const_def, only: dp
2309 :
2310 : integer, intent(in) :: handle
2311 :
2312 : real(dp), intent(in) :: Z ! the metals mass fraction
2313 : real(dp), intent(in) :: X ! the hydrogen mass fraction
2314 :
2315 : real(dp), intent(in) :: abar, zbar
2316 :
2317 : integer, intent(in) :: species
2318 : integer, pointer :: chem_id(:)
2319 : integer, pointer :: net_iso(:)
2320 : real(dp), intent(in) :: xa(:)
2321 :
2322 : real(dp), intent(in) :: logT ! log10 of temperature
2323 :
2324 : integer, intent(in) :: which_other ! from eos_def.
2325 : real(dp), intent(in) :: other_value ! desired value for the other variable
2326 : real(dp), intent(in) :: other_tol
2327 :
2328 : real(dp), intent(in) :: logRho_tol
2329 :
2330 : integer, intent(in) :: max_iter ! max number of Newton iterations
2331 :
2332 : real(dp), intent(in) :: logRho_guess
2333 : real(dp), intent(in) :: logRho_bnd1, logRho_bnd2 ! bounds for logrho
2334 : ! set to arg_not_provided if do not know bounds
2335 : real(dp), intent(in) :: other_at_bnd1, other_at_bnd2 ! values at bounds
2336 : ! if don't know these values, just set to arg_not_provided (defined in c_def)
2337 :
2338 : real(dp), intent(out) :: logRho_result
2339 : real(dp), intent(inout), dimension(nv) :: res, d_dlnRho_c_T, d_dlnT_c_Rho
2340 : real(dp), intent(inout), dimension(:,:) :: d_dxa_c_TRho
2341 :
2342 : integer, intent(out) :: eos_calls
2343 : integer, intent(out) :: ierr ! 0 means AOK.
2344 :
2345 : logical, parameter :: doing_Rho = .true.
2346 : real(dp) :: Prad
2347 :
2348 : call do_safe_get_Rho_T( &
2349 : handle, Z, X, abar, zbar, &
2350 : species, chem_id, net_iso, xa, &
2351 : logT, which_other, other_value, doing_Rho, &
2352 : logRho_guess, logRho_result, logRho_bnd1, logRho_bnd2, other_at_bnd1, other_at_bnd2, &
2353 : logRho_tol, other_tol, max_iter, res, d_dlnRho_c_T, d_dlnT_c_Rho, d_dxa_c_TRho, &
2354 12152 : eos_calls, ierr)
2355 :
2356 12152 : end subroutine get_Rho
2357 :
2358 :
2359 13628 : subroutine do_safe_get_Rho_T( &
2360 : handle, Z, XH1, abar, zbar, &
2361 13628 : species, chem_id, net_iso, xa, &
2362 : the_other_log, which_other, other_value, doing_Rho, &
2363 : initial_guess, x, xbnd1, xbnd2, other_at_bnd1, other_at_bnd2, &
2364 : xacc, yacc, ntry, &
2365 13628 : res, d_dlnRho_c_T, d_dlnT_c_Rho, d_dxa_c_TRho, &
2366 : eos_calls, ierr)
2367 : use const_def, only: dp
2368 : use chem_def, only: num_chem_isos
2369 : use num_lib, only: safe_root_with_guess
2370 : integer, intent(in) :: handle
2371 : real(dp), intent(in) :: Z, XH1, abar, zbar
2372 : integer, intent(in) :: species
2373 : integer, pointer :: chem_id(:)
2374 : integer, pointer :: net_iso(:)
2375 : real(dp), intent(in) :: xa(:)
2376 : integer, intent(in) :: which_other ! 0 means total P
2377 : real(dp), intent(in) :: other_value
2378 : logical, intent(in) :: doing_Rho
2379 : real(dp), intent(in) :: initial_guess ! for x
2380 : real(dp), intent(out) :: x ! if doing_Rho, then logRho, else logT
2381 : real(dp), intent(in) :: the_other_log
2382 : real(dp), intent(in) :: xbnd1, xbnd2, other_at_bnd1, other_at_bnd2
2383 : real(dp), intent(in) :: xacc, yacc ! tolerances
2384 : integer, intent(in) :: ntry ! max number of iterations
2385 : real(dp), intent(inout), dimension(nv) :: res, d_dlnRho_c_T, d_dlnT_c_Rho
2386 : real(dp), dimension(:,:) :: d_dxa_c_TRho
2387 : integer, intent(out) :: eos_calls, ierr
2388 :
2389 : integer :: i, j, ix, iz
2390 : integer, parameter :: lrpar = 0, lipar = 0, newt_imax = 6
2391 : real(dp), parameter :: dx = 0.1d0
2392 13628 : integer, pointer :: ipar(:)
2393 13628 : real(dp), pointer :: rpar(:)
2394 13628 : real(dp) :: the_other_val, logRho, logT, rho, T, x1, x3, y1, y3
2395 : type (EoS_General_Info), pointer :: rq
2396 :
2397 : include 'formats'
2398 :
2399 : ierr = 0
2400 :
2401 13628 : call get_eos_ptr(handle, rq, ierr)
2402 13628 : if (ierr /= 0) then
2403 0 : write(*, *) 'get_eos_ptr returned ierr', ierr
2404 : return
2405 : end if
2406 :
2407 13628 : x1 = arg_not_provided
2408 13628 : x3 = arg_not_provided
2409 13628 : y1 = arg_not_provided
2410 13628 : y3 = arg_not_provided
2411 :
2412 13628 : eos_calls = 0
2413 13628 : the_other_val = exp10(the_other_log)
2414 13628 : nullify(ipar, rpar)
2415 :
2416 : x = safe_root_with_guess( &
2417 : f, initial_guess, dx, x1, x3, y1, y3, &
2418 : min(ntry,newt_imax), ntry, xacc, yacc, &
2419 27256 : lrpar, rpar, lipar, ipar, ierr)
2420 :
2421 : contains
2422 :
2423 45564 : real(dp) function f(x, dfdx, lrpar, rpar, lipar, ipar, ierr)
2424 : ! returns with ierr = 0 if was able to evaluate f and df/dx at x
2425 : ! if df/dx not available, it is okay to set it to 0
2426 13628 : use const_def, only: dp
2427 : integer, intent(in) :: lrpar, lipar
2428 : real(dp), intent(in) :: x
2429 : real(dp), intent(out) :: dfdx
2430 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2431 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2432 : integer, intent(out) :: ierr
2433 :
2434 45564 : real(dp) :: Pgas, Prad, energy, entropy, dPgas_dlnT, dPrad_dlnT, &
2435 45564 : dPgas_dlnRho, erad, egas, derad_dlnT, degas_dlnT, derad_dlnRho
2436 :
2437 : include 'formats'
2438 45564 : ierr = 0
2439 45564 : eos_calls = eos_calls + 1
2440 45564 : f = 0; dfdx = 0
2441 :
2442 45564 : if (x > 50d0) then
2443 0 : ierr = -1
2444 0 : return
2445 : end if
2446 :
2447 45564 : if (doing_Rho) then
2448 42460 : logRho = x
2449 42460 : rho = exp10(logRho)
2450 42460 : logT = the_other_log
2451 42460 : T = the_other_val
2452 : else
2453 3104 : logT = x
2454 3104 : T = exp10(logT)
2455 3104 : logRho = the_other_log
2456 3104 : rho = the_other_val
2457 : end if
2458 :
2459 : call Get_eosDT_Results(rq, Z, XH1, abar, zbar, &
2460 : species, chem_id, net_iso, xa, &
2461 : rho, logRho, T, logT, &
2462 45564 : res, d_dlnRho_c_T, d_dlnT_c_Rho, d_dxa_c_TRho, ierr)
2463 :
2464 45564 : Pgas = exp(res(i_lnPgas))
2465 45564 : Prad = crad*T*T*T*T/3d0
2466 45564 : energy = exp(res(i_lnE))
2467 45564 : entropy = exp(res(i_lnS))
2468 :
2469 45564 : if (ierr /= 0) then
2470 : if (.false.) then
2471 : write(*,*) 'Get_eosDT_Results returned ierr', ierr
2472 : write(*,1) 'Z', Z
2473 : write(*,1) 'XH1', XH1
2474 : write(*,1) 'abar', abar
2475 : write(*,1) 'zbar', zbar
2476 : write(*,1) 'rho', rho
2477 : write(*,1) 'logRho', logRho
2478 : write(*,1) 'T', T
2479 : write(*,1) 'logT', logT
2480 : write(*,'(A)')
2481 : call mesa_error(__FILE__,__LINE__,'do_safe_get_Rho_T')
2482 : end if
2483 : return
2484 : end if
2485 :
2486 45564 : if (is_bad(res(i_Cv))) then
2487 0 : ierr = -1
2488 0 : if (.not. stop_for_is_bad) return
2489 : write(*,1) 'res(i_Cv)', res(i_Cv)
2490 : write(*,1) 'Z', Z
2491 : write(*,1) 'XH1', XH1
2492 : write(*,1) 'abar', abar
2493 : write(*,1) 'zbar', zbar
2494 : write(*,1) 'rho', rho
2495 : write(*,1) 'logRho', logRho
2496 : write(*,1) 'T', T
2497 : write(*,1) 'logT', logT
2498 : write(*,'(A)')
2499 : call mesa_error(__FILE__,__LINE__,'do_safe_get_Rho_T')
2500 : end if
2501 :
2502 45564 : if (which_other == -1) then ! other_value is egas
2503 0 : erad = crad*pow4(T)/rho
2504 0 : egas = energy - erad
2505 0 : f = egas - other_value
2506 0 : if (doing_Rho) then
2507 0 : derad_dlnRho = -erad
2508 0 : dfdx = energy*d_dlnRho_c_T(i_lnE)*ln10 - derad_dlnRho
2509 : else
2510 0 : derad_dlnT = 4d0*erad
2511 0 : degas_dlnT = energy*d_dlnT_c_Rho(i_lnE) - derad_dlnT
2512 0 : dfdx = degas_dlnT*ln10
2513 : end if
2514 45564 : else if (which_other == 0) then ! other_value is log10P
2515 36 : f = log10(Pgas + Prad) - other_value
2516 36 : if (doing_Rho) then
2517 0 : dPgas_dlnRho = Pgas*d_dlnRho_c_T(i_lnPgas)
2518 0 : dfdx = dPgas_dlnRho/(Pgas + Prad)*ln10
2519 : else
2520 36 : dPgas_dlnT = Pgas*d_dlnT_c_Rho(i_lnPgas)
2521 36 : dPrad_dlnT = 4d0*Prad
2522 36 : dfdx = (dPgas_dlnT + dPrad_dlnT)/(Pgas + Prad)*ln10
2523 : end if
2524 : else
2525 45528 : f = res(which_other) - other_value
2526 45528 : if (doing_Rho) then
2527 42460 : dfdx = d_dlnRho_c_T(which_other)*ln10
2528 : else
2529 3068 : dfdx = d_dlnT_c_Rho(which_other)*ln10
2530 : end if
2531 : end if
2532 :
2533 : end function f
2534 :
2535 : end subroutine do_safe_get_Rho_T
2536 :
2537 : end module eosDT_eval
|
