Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2018-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 rsp_eval_eos_and_kap
21 :
22 : use eos_def
23 : use eos_lib
24 : use chem_def
25 : use chem_lib, only: chem_Xsol, basic_composition_info
26 : use kap_lib
27 : use kap_def
28 : use const_def, only: dp, ln10, arg_not_provided
29 : use utils_lib
30 : use star_utils, only: &
31 : store_rho_in_xh, store_lnd_in_xh, get_rho_and_lnd_from_xh, &
32 : store_T_in_xh, store_lnT_in_xh, get_T_and_lnT_from_xh
33 : use star_private_def
34 : use rsp_def, only: xa, X, Z, Y, &
35 : abar, zbar, z53bar, XC, XN, XO, Xne
36 :
37 : implicit none
38 :
39 : integer :: species
40 : integer, pointer, dimension(:) :: net_iso, chem_id
41 : integer :: eos_handle, kap_handle
42 :
43 : contains
44 :
45 0 : subroutine restart_rsp_eos_and_kap(s)
46 : type (star_info), pointer :: s
47 0 : eos_handle = s% eos_handle
48 0 : kap_handle = s% kap_handle
49 0 : net_iso => s% net_iso
50 0 : chem_id => s% chem_id
51 0 : species = s% species
52 0 : end subroutine restart_rsp_eos_and_kap
53 :
54 0 : subroutine init_for_rsp_eos_and_kap(s)
55 : use adjust_xyz, only: get_xa_for_standard_metals
56 : type (star_info), pointer :: s
57 :
58 : integer :: i, iz, ierr
59 : real(dp) :: initial_z, initial_y, initial_x, &
60 0 : initial_h1, initial_h2, initial_he3, initial_he4, &
61 0 : xsol_he3, xsol_he4, z2bar, ye, mass_correction, sumx
62 : include 'formats'
63 0 : ierr = 0
64 0 : eos_handle = s% eos_handle
65 0 : kap_handle = s% kap_handle
66 0 : net_iso => s% net_iso
67 0 : chem_id => s% chem_id
68 0 : species = s% species
69 :
70 0 : initial_x = max(0d0, min(1d0, s% RSP_X))
71 0 : initial_z = max(0d0, min(1d0, s% RSP_Z))
72 0 : initial_y = max(0d0,1d0 - (initial_x + initial_z))
73 0 : initial_h1 = initial_x
74 0 : initial_h2 = 0
75 0 : if (s% initial_he3 < 0d0) then
76 0 : xsol_he3 = chem_Xsol('he3')
77 0 : xsol_he4 = chem_Xsol('he4')
78 0 : initial_he3 = initial_y*xsol_he3/(xsol_he3 + xsol_he4)
79 0 : initial_he4 = initial_y*xsol_he4/(xsol_he3 + xsol_he4)
80 0 : else if (s% initial_he3 < initial_y) then
81 0 : initial_he3 = s% initial_he3
82 0 : initial_he4 = initial_y - s% initial_he3
83 : else
84 0 : write(*,*) 'ERROR: initial_he3 is larger than initial_y'
85 : ierr = -1
86 0 : return
87 : end if
88 : call get_xa_for_standard_metals(s, &
89 : species, chem_id, net_iso, &
90 : initial_h1, initial_h2, initial_he3, initial_he4, &
91 : s% job% initial_zfracs, s% initial_dump_missing_heaviest, &
92 0 : xa, ierr)
93 0 : if (ierr /= 0) then
94 0 : write(*,*) 'failed in get_xa_for_standard_metals'
95 0 : return
96 : end if
97 0 : if (abs(1-sum(xa(:))) > 1d-8) then
98 0 : write(*,1) 'initial_h1', initial_h1
99 0 : write(*,1) 'initial_h2', initial_h2
100 0 : write(*,1) 'initial_he3', initial_he3
101 0 : write(*,1) 'initial_he4', initial_he4
102 0 : write(*,1) 'initial_y', initial_y
103 0 : write(*,1) 'initial_z', initial_z
104 0 : write(*,1) 'initial_h1+h2+he3+he4+z', &
105 0 : initial_h1 + initial_h2 + initial_he3 + initial_he4 + initial_z
106 0 : write(*,1) 'sum(xa(:))', sum(xa(:))
107 0 : write(*,*) 'init_for_rsp_eos_and_kap'
108 : ierr = -1
109 0 : return
110 : end if
111 : call basic_composition_info( &
112 : species, s% chem_id, xa(:), X, Y, Z, abar, zbar, z2bar, z53bar, ye, &
113 0 : mass_correction, sumx)
114 0 : if (is_bad(zbar)) then
115 0 : write(*,1) 'basic_composition_info initial_x', initial_x
116 0 : write(*,1) 'basic_composition_info initial_z', initial_z
117 0 : write(*,1) 'basic_composition_info initial_y', initial_y
118 0 : write(*,1) 'basic_composition_info zbar', zbar
119 0 : do i=1,species
120 0 : write(*,2) 'xa', i, xa(i)
121 : end do
122 0 : stop
123 : end if
124 0 : xc = 0; xn = 0; xo = 0; xne = 0
125 0 : do i=1, species
126 0 : iz = chem_isos% Z(chem_id(i))
127 0 : select case(iz)
128 : case (6)
129 0 : xc = xc + xa(i)
130 : case (7)
131 0 : xn = xn + xa(i)
132 : case (8)
133 0 : xo = xo + xa(i)
134 : case (10)
135 0 : xne = xne + xa(i)
136 : end select
137 : end do
138 :
139 : return
140 :
141 : write(*,1) 'init_for_rsp_eos_and_kap X', X
142 : write(*,1) 'Y', Y
143 : write(*,1) 'Z', Z
144 :
145 :
146 : write(*,1) 'abar', abar
147 : write(*,1) 'zbar', zbar
148 : write(*,1) 'XC', XC
149 : write(*,1) 'XN', XN
150 : write(*,1) 'XO', XO
151 : write(*,1) 'Xne', Xne
152 : do i=1,species
153 : write(*,2) 'xa', i, xa(i)
154 : end do
155 : write(*,*) trim(s% net_name)
156 : call mesa_error(__FILE__,__LINE__,'init_for_rsp_eos_and_kap')
157 0 : end subroutine init_for_rsp_eos_and_kap
158 :
159 :
160 0 : subroutine eval_mesa_eos_and_kap(&
161 : s,k,T_in,V, &
162 : Pg,d_Pg_dV,d_Pg_dT,Pr,d_Pr_dT, &
163 : egas,d_egas_dV,d_egas_dT,erad,d_erad_dV,d_erad_dT, &
164 : CSND,CP,CPV,CPT,Q,QV,QT,OP,OPV,OPT,ierr)
165 : type (star_info), pointer :: s
166 : integer, intent(in) :: k
167 : real(dp), intent(in) :: T_in,V
168 : real(dp), intent(out) :: &
169 : Pg,d_Pg_dV,d_Pg_dT,Pr,d_Pr_dT, &
170 : egas,d_egas_dV,d_egas_dT,erad,d_erad_dV,d_erad_dT, &
171 : CSND,CP,CPV,CPT,Q,QV,QT,OP,OPV,OPT
172 : integer, intent(out) :: ierr
173 : call eval1_mesa_eos_and_kap(&
174 : s,k,.false.,T_in,V, &
175 : Pg,d_Pg_dV,d_Pg_dT,Pr,d_Pr_dT, &
176 : egas,d_egas_dV,d_egas_dT,erad,d_erad_dV,d_erad_dT, &
177 0 : CSND,CP,CPV,CPT,Q,QV,QT,OP,OPV,OPT,ierr)
178 0 : end subroutine eval_mesa_eos_and_kap
179 :
180 :
181 0 : subroutine eval1_mesa_eos_and_kap( &
182 : s,k,skip_kap,T_in,V, &
183 : Pgas,d_Pg_dV,d_Pg_dT,Prad,d_Pr_dT, &
184 : egas,d_egas_dV,d_egas_dT,erad,d_erad_dV,d_erad_dT, &
185 : CSND,CP,CPV,CPT,Q,QV,QT,OP,OPV,OPT,ierr)
186 : use star_utils, only: write_eos_call_info
187 : use eos_support, only: get_eos
188 : use micro, only: store_eos_for_cell
189 : type (star_info), pointer :: s
190 : integer, intent(in) :: k
191 : logical, intent(in) :: skip_kap
192 : real(dp), intent(in) :: T_in,V
193 : real(dp), intent(out) :: &
194 : Pgas,d_Pg_dV,d_Pg_dT,Prad,d_Pr_dT, &
195 : egas,d_egas_dV,d_egas_dT,erad,d_erad_dV,d_erad_dT, &
196 : CSND,CP,CPV,CPT,Q,QV,QT,OP,OPV,OPT
197 : integer, intent(out) :: ierr
198 :
199 : integer :: j
200 0 : real(dp) :: logT, logRho, T, Rho, E, dE_dV, dE_dT, &
201 0 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
202 0 : dlnd_dV, chiRho, chiT, dchiRho_dlnd, dchiRho_dlnT, &
203 0 : dchiT_dlnd, dchiT_dlnT, dQ_dlnd, dQ_dlnT
204 0 : real(dp), dimension(num_eos_basic_results) :: res, d_dlnd, d_dlnT
205 0 : real(dp) :: d_dxa(num_eos_d_dxa_results,s% species)
206 :
207 : include 'formats'
208 :
209 0 : rho = 1d0/V
210 0 : T = T_in
211 0 : logRho = log10(rho)
212 0 : dlnd_dV = -rho
213 0 : logT = log10(T)
214 :
215 0 : if (k > 0 .and. k <= s% nz) then
216 0 : call store_rho_in_xh(s, k, rho)
217 0 : call get_rho_and_lnd_from_xh(s, k, s% rho(k), s% lnd(k))
218 0 : call store_T_in_xh(s, k, T)
219 0 : call get_T_and_lnT_from_xh(s, k, s% T(k), s% lnT(k))
220 0 : s% abar(k) = abar
221 0 : s% zbar(k) = zbar
222 0 : s% z53bar(k) = z53bar
223 : call get_eos( &
224 : s, k, xa, &
225 : Rho, logRho, T, logT, &
226 : res, d_dlnd, d_dlnT, &
227 0 : d_dxa, ierr)
228 0 : if (ierr == 0) then
229 0 : s% lnPgas(k) = res(i_lnPgas)
230 0 : s% Pgas(k) = exp(s% lnPgas(k))
231 0 : do j=1,num_eos_basic_results
232 0 : s% d_eos_dlnd(j,k) = d_dlnd(j)
233 0 : s% d_eos_dlnT(j,k) = d_dlnT(j)
234 : end do
235 0 : call store_eos_for_cell(s, k, res, d_dlnd, d_dlnT, d_dxa, ierr)
236 0 : CSND = s% csound(k)
237 : end if
238 : else ! k <= 0 or k > nz
239 : call get_eos( &
240 : s, 0, xa, &
241 : Rho, logRho, T, logT, &
242 : res, d_dlnd, d_dlnT, &
243 0 : d_dxa, ierr)
244 0 : if (ierr == 0) then
245 0 : Prad = crad*T**4/3d0
246 0 : Pgas = exp(res(i_lnPgas))
247 0 : CSND = sqrt(max(0d0, res(i_gamma1)*(Prad+Pgas)/Rho))
248 : end if
249 : end if
250 :
251 0 : if (ierr /= 0) then
252 0 : !$omp critical (rsp_eval_eos_and_kap_1)
253 0 : if (k > 0 .and. k < s% nz) call write_eos_call_info(s,k)
254 0 : write(*,2) 'X', k, X
255 0 : write(*,2) 'Z', k, Z
256 0 : write(*,2) 'zbar', k, zbar
257 0 : write(*,2) 'abar', k, abar
258 0 : write(*,2) 'V', k, V
259 0 : write(*,2) 'rho', k, rho
260 0 : write(*,2) 'T', k, T
261 0 : write(*,2) 'logRho', k, logRho
262 0 : write(*,2) 'logT', k, logT
263 0 : if (s% stop_for_bad_nums .and. is_bad(logRho+logT)) call mesa_error(__FILE__,__LINE__,'do_eos_for_cell')
264 : !$omp end critical (rsp_eval_eos_and_kap_1)
265 : !return
266 0 : call mesa_error(__FILE__,__LINE__,'RSP failed in get_eos')
267 : end if
268 :
269 0 : if (skip_kap) then
270 0 : OP=0; OPV=0; OPT=0
271 : else
272 : call eval1_kap(s, k, skip_kap, &
273 : V, logRho, dlnd_dV, T, logT, species, chem_id, net_iso, xa, &
274 0 : res, d_dlnd, d_dlnT, OP, OPV, OPT, ierr)
275 : end if
276 :
277 0 : lnfree_e = res(i_lnfree_e)
278 0 : d_lnfree_e_dlnRho = d_dlnd(i_lnfree_e)
279 0 : d_lnfree_e_dlnT = d_dlnT(i_lnfree_e)
280 :
281 0 : Prad = crad*T**4/3d0 ! erg/cm^2
282 0 : d_Pr_dT = 4d0*Prad/T
283 :
284 0 : erad = 3d0*Prad/rho ! 3*Prad*V erg/gm
285 0 : d_erad_dT = 3d0*d_Pr_dT/rho
286 0 : d_erad_dV = 3d0*Prad
287 :
288 0 : E = exp(res(i_lnE))
289 0 : dE_dV = E*d_dlnd(i_lnE)*dlnd_dV
290 0 : dE_dT = E*d_dlnT(i_lnE)/T
291 0 : egas = E - erad
292 0 : d_egas_dV = dE_dV - d_erad_dV
293 0 : d_egas_dT = dE_dT - d_erad_dT
294 :
295 0 : Pgas = exp(res(i_lnPgas))
296 0 : d_Pg_dV = Pgas*d_dlnd(i_lnPgas)*dlnd_dV
297 0 : d_Pg_dT = Pgas*d_dlnT(i_lnPgas)/T
298 :
299 0 : CP = res(i_Cp)
300 0 : CPV = d_dlnd(i_Cp)*dlnd_dV
301 0 : CPT = d_dlnT(i_Cp)/T
302 :
303 0 : chiT = res(i_chiT)
304 0 : dchiT_dlnd = d_dlnd(i_chiT)
305 0 : dchiT_dlnT = d_dlnT(i_chiT)
306 :
307 0 : chiRho = res(i_chiRho)
308 0 : dchiRho_dlnd = d_dlnd(i_chiRho)
309 0 : dchiRho_dlnT = d_dlnT(i_chiRho)
310 :
311 0 : Q = chiT/(rho*T*chiRho) ! thermal expansion coefficient
312 0 : dQ_dlnd = Q*(dchiT_dlnd/chiT - dchiRho_dlnd/chiRho - 1d0)
313 0 : dQ_dlnT = Q*(dchiT_dlnT/chiT - dchiRho_dlnT/chiRho - 1d0)
314 0 : QV = dQ_dlnd*dlnd_dV
315 0 : QT = dQ_dlnT/T
316 :
317 0 : if (is_bad(egas) .or. egas <= 0d0) then
318 0 : ierr = -1
319 : return
320 : end if
321 :
322 0 : end subroutine eval1_mesa_eos_and_kap
323 :
324 :
325 0 : subroutine eval1_mesa_eosDEgas_and_kap( &
326 : s, k, skip_kap, egas, V, T, Pgas, CSND, CP, Q, OP, ierr)
327 0 : use star_utils, only: write_eos_call_info
328 : use eos_support, only: solve_eos_given_DEgas
329 : use micro, only: store_eos_for_cell
330 : type (star_info), pointer :: s
331 : integer, intent(in) :: k
332 : logical, intent(in) :: skip_kap
333 : real(dp), intent(in) :: egas, V
334 : real(dp), intent(out) :: T, Pgas, CSND, CP, Q, OP
335 : integer, intent(out) :: ierr
336 :
337 : integer :: j, eos_calls
338 0 : real(dp) :: rho, logRho, dlnd_dV, egas_tol, logT, &
339 0 : logT_guess, logT_tol, new_erad, new_egas, OPV, OPT
340 : real(dp), dimension(num_eos_basic_results) :: &
341 0 : res, d_dlnd, d_dlnT
342 0 : real(dp) :: d_dxa(num_eos_d_dxa_results, species)
343 :
344 : include 'formats'
345 0 : ierr = 0
346 :
347 0 : rho = 1d0/V
348 0 : logRho = log10(rho)
349 0 : dlnd_dV = -rho
350 :
351 0 : if (egas <= 0d0 .or. is_bad(egas)) then
352 0 : !$OMP critical (RSP_eosDEgas)
353 0 : write(*,2) 'egas', k, egas
354 0 : write(*,*) 'called eval1_mesa_eosDEgas_and_kap with bad value for egas'
355 0 : call mesa_error(__FILE__,__LINE__,'eval1_mesa_eosDEgas_and_kap')
356 : !$OMP end critical (RSP_eosDEgas)
357 0 : ierr = -1
358 0 : return
359 : end if
360 :
361 0 : if (k > 0 .and. k <= s% nz) then
362 0 : egas_tol = egas*1d-11
363 0 : s% rho(k) = rho
364 0 : s% lnd(k) = logRho*ln10
365 0 : s% xh(s% i_lnd,k) = s% lnd(k)
366 0 : s% abar(k) = abar
367 0 : s% zbar(k) = zbar
368 0 : logT_guess = s% lnT(k)/ln10
369 0 : logT_tol = 1d-11
370 : call solve_eos_given_DEgas( &
371 : s, k, xa, &
372 : logRho, egas, logT_guess, logT_tol, egas_tol, &
373 : logT, res, d_dlnd, d_dlnT, d_dxa, &
374 0 : ierr)
375 0 : if (ierr == 0) then
376 0 : call store_lnT_in_xh(s, k, logT*ln10)
377 0 : call get_T_and_lnT_from_xh(s, k, s% T(k), s% lnT(k))
378 0 : T = s% T(k)
379 0 : new_erad = crad*T**4/rho
380 0 : new_egas = exp(res(i_lnE)) - new_erad
381 0 : if (is_bad(new_egas) .or. new_egas <= 0d0 .or. &
382 : abs(new_egas - egas) > 1d3*egas_tol) then
383 0 : !$OMP critical (RSP_eosDEgas)
384 0 : write(*,1) 'logRho', s% lnd(k)/ln10
385 0 : write(*,1) 'logT_guess', logT_guess
386 0 : write(*,1) 'egas', egas
387 0 : write(*,1) 'Z', Z
388 0 : write(*,1) 'X', X
389 0 : write(*,1) 'abar', abar
390 0 : write(*,1) 'zbar', zbar
391 0 : write(*,1) 'logT_tol', logT_tol
392 0 : write(*,1) 'egas_tol', egas_tol
393 0 : write(*,'(A)')
394 0 : write(*,1) 'guess logT', logT_guess
395 0 : write(*,1) 'found logT', logT
396 0 : write(*,1) 'wanted egas', egas
397 0 : write(*,1) 'got egas', new_egas
398 0 : write(*,1) '(want - got)/got', (egas - new_egas)/new_egas
399 0 : write(*,'(A)')
400 0 : write(*,2) 'eos_calls', eos_calls
401 0 : write(*,'(A)')
402 0 : write(*,2) 'failed eval1_mesa_eosDEgas_and_kap', k
403 0 : write(*,'(A)')
404 0 : write(*,*) 'is_bad(new_egas)', is_bad(new_egas)
405 0 : write(*,*) 'new_egas <= 0d0', new_egas <= 0d0
406 0 : write(*,*) 'abs(new_egas - egas) > egas_tol', &
407 0 : abs(new_egas - egas) > egas_tol, &
408 0 : abs(new_egas - egas) - egas_tol, &
409 0 : abs(new_egas - egas), egas_tol
410 0 : call mesa_error(__FILE__,__LINE__,'eval1_mesa_eosDEgas_and_kap')
411 : !$OMP end critical (RSP_eosDEgas)
412 : end if
413 0 : s% lnPgas(k) = res(i_lnPgas)
414 0 : s% Pgas(k) = exp(s% lnPgas(k))
415 0 : Pgas = s% Pgas(k)
416 0 : do j=1,num_eos_basic_results
417 0 : s% d_eos_dlnd(j,k) = d_dlnd(j)
418 0 : s% d_eos_dlnT(j,k) = d_dlnT(j)
419 : end do
420 0 : call store_eos_for_cell(s, k, res, d_dlnd, d_dlnT, d_dxa, ierr)
421 0 : CSND = s% csound(k)
422 : end if
423 : else ! k <= 0 or k > nz
424 0 : write(*,*) 'cannot call eval1_mesa_eosDEgas_and_kap with k <= 0 or k > nz'
425 0 : ierr = -1
426 0 : return
427 : end if
428 :
429 0 : if (ierr /= 0) then
430 0 : !$OMP critical (RSP_eosDEgas)
431 0 : if (k > 0 .and. k < s% nz) call write_eos_call_info(s,k)
432 0 : write(*,2) 'X', k, X
433 0 : write(*,2) 'Z', k, Z
434 0 : write(*,2) 'zbar', k, zbar
435 0 : write(*,2) 'abar', k, abar
436 0 : write(*,2) 'V', k, V
437 0 : write(*,2) 'rho', k, rho
438 0 : write(*,2) 'egas', k, egas
439 0 : write(*,2) 'logRho', k, logRho
440 0 : write(*,2) 'T', k, T
441 0 : write(*,2) 'logT', k, logT
442 : if (s% stop_for_bad_nums .and. egas <= 0d0) call mesa_error(__FILE__,__LINE__,'do_eos_for_cell')
443 : !$OMP end critical (RSP_eosDEgas)
444 0 : return
445 : call mesa_error(__FILE__,__LINE__,'RSP failed in eval1_mesa_eosDEgas_and_kap')
446 : end if
447 :
448 0 : if (skip_kap) then
449 0 : OP=0
450 : else
451 : call eval1_kap(s, k, skip_kap, &
452 : V, logRho, dlnd_dV, T, logT, species, chem_id, net_iso, xa, &
453 0 : res, d_dlnd, d_dlnT, OP, OPV, OPT, ierr)
454 : end if
455 :
456 0 : Pgas = exp(res(i_lnPgas))
457 0 : CP = res(i_Cp)
458 0 : Q = res(i_chiT)/(rho*T*res(i_chiRho)) ! thermal expansion coefficient
459 :
460 0 : end subroutine eval1_mesa_eosDEgas_and_kap
461 :
462 :
463 0 : subroutine eval1_mesa_eosDE_and_kap( & ! for eos, energy = egas + erad
464 : s, k, skip_kap, energy, V, T, Pgas, CSND, CP, Q, OP, ierr)
465 0 : use star_utils, only: write_eos_call_info
466 : use eos_support, only: solve_eos_given_DE
467 : use micro, only: store_eos_for_cell
468 : type (star_info), pointer :: s
469 : integer, intent(in) :: k
470 : logical, intent(in) :: skip_kap
471 : real(dp), intent(in) :: energy, V
472 : real(dp), intent(out) :: T, Pgas, CSND, CP, Q, OP
473 : integer, intent(out) :: ierr
474 :
475 : integer :: j, eos_calls
476 0 : real(dp) :: rho, logRho, dlnd_dV, logE, logE_want, logE_tol, logT, &
477 0 : logT_guess, logT_tol, new_erad, new_egas, OPV, OPT
478 : real(dp), dimension(num_eos_basic_results) :: &
479 0 : res, d_dlnd, d_dlnT
480 0 : real(dp) :: d_dxa(num_eos_d_dxa_results, species)
481 :
482 : include 'formats'
483 0 : ierr = 0
484 :
485 0 : rho = 1d0/V
486 0 : logRho = log10(rho)
487 0 : logE_want = log10(energy)
488 0 : dlnd_dV = -rho
489 :
490 0 : if (energy <= 0d0 .or. is_bad(energy)) then
491 0 : !$OMP critical (RSP_eosDE)
492 0 : write(*,2) 'energy', k, energy
493 0 : write(*,*) 'called eval1_mesa_eosDE_and_kap with bad value for energy'
494 0 : call mesa_error(__FILE__,__LINE__,'eval1_mesa_eosDE_and_kap')
495 : !$OMP end critical (RSP_eosDE)
496 0 : ierr = -1
497 0 : return
498 : end if
499 :
500 0 : if (k > 0 .and. k <= s% nz) then
501 0 : s% rho(k) = rho
502 0 : s% lnd(k) = logRho*ln10
503 0 : s% xh(s% i_lnd,k) = s% lnd(k)
504 0 : s% abar(k) = abar
505 0 : s% zbar(k) = zbar
506 0 : logT_guess = s% lnT(k)/ln10
507 0 : logT_tol = 1d-11
508 0 : logE_tol = 1d-11
509 :
510 : call solve_eos_given_DE( &
511 : s, k, xa, &
512 : logRho, logE_want, logT_guess, logT_tol, logE_tol, &
513 : logT, res, d_dlnd, d_dlnT, d_dxa, &
514 0 : ierr)
515 0 : if (ierr == 0) then
516 0 : call store_lnT_in_xh(s, k, logT*ln10)
517 0 : call get_T_and_lnT_from_xh(s, k, s% T(k), s% lnT(k))
518 0 : T = s% T(k)
519 0 : new_erad = crad*T**4/rho
520 0 : new_egas = exp(res(i_lnE)) - new_erad
521 0 : logE = res(i_lnE)/ln10
522 0 : if (is_bad(logE) .or. logE <= -20d0) then
523 0 : !$OMP critical (RSP_eosDE)
524 0 : write(*,1) 'logRho', s% lnd(k)/ln10
525 0 : write(*,1) 'Z', Z
526 0 : write(*,1) 'X', X
527 0 : write(*,1) 'abar', abar
528 0 : write(*,1) 'zbar', zbar
529 0 : write(*,1) 'logT_tol', logT_tol
530 0 : write(*,1) 'logE_tol', logE_tol
531 0 : write(*,'(A)')
532 0 : write(*,1) 'guess logT', logT_guess
533 0 : write(*,1) 'found logT', logT
534 0 : write(*,'(A)')
535 0 : write(*,1) 'wanted logE', logE_want
536 0 : write(*,1) 'got logE', logE
537 0 : write(*,'(A)')
538 0 : write(*,2) 'eos_calls', eos_calls
539 0 : write(*,'(A)')
540 0 : write(*,2) 'failed eval1_mesa_eosDE_and_kap', k
541 0 : write(*,'(A)')
542 0 : call mesa_error(__FILE__,__LINE__,'eval1_mesa_eosDE_and_kap')
543 : !$OMP end critical (RSP_eosDE)
544 : end if
545 0 : s% lnPgas(k) = res(i_lnPgas)
546 0 : s% Pgas(k) = exp(s% lnPgas(k))
547 0 : Pgas = s% Pgas(k)
548 0 : do j=1,num_eos_basic_results
549 0 : s% d_eos_dlnd(j,k) = d_dlnd(j)
550 0 : s% d_eos_dlnT(j,k) = d_dlnT(j)
551 : end do
552 0 : call store_eos_for_cell(s, k, res, d_dlnd, d_dlnT, d_dxa, ierr)
553 0 : CSND = s% csound(k)
554 : end if
555 : else ! k <= 0 or k > nz
556 0 : write(*,*) 'cannot call eval1_mesa_eosDE_and_kap with k <= 0 or k > nz'
557 0 : ierr = -1
558 0 : return
559 : end if
560 :
561 0 : if (ierr /= 0) then
562 0 : !$OMP critical (RSP_eosDE)
563 0 : if (k > 0 .and. k < s% nz) call write_eos_call_info(s,k)
564 0 : write(*,2) 'X', k, X
565 0 : write(*,2) 'Z', k, Z
566 0 : write(*,2) 'zbar', k, zbar
567 0 : write(*,2) 'abar', k, abar
568 0 : write(*,2) 'V', k, V
569 0 : write(*,2) 'rho', k, rho
570 0 : write(*,2) 'logE', k, logE
571 0 : write(*,2) 'logRho', k, logRho
572 0 : write(*,2) 'T', k, T
573 0 : write(*,2) 'logT', k, logT
574 : !$OMP end critical (RSP_eosDE)
575 0 : return
576 : call mesa_error(__FILE__,__LINE__,'RSP failed in eval1_mesa_eosDE_and_kap')
577 : end if
578 :
579 0 : if (skip_kap) then
580 0 : OP=0
581 : else
582 : call eval1_kap(s, k, skip_kap, &
583 : V, logRho, dlnd_dV, T, logT, species, chem_id, net_iso, xa, &
584 0 : res, d_dlnd, d_dlnT, OP, OPV, OPT, ierr)
585 : end if
586 :
587 0 : Pgas = exp(res(i_lnPgas))
588 0 : CP = res(i_Cp)
589 0 : Q = res(i_chiT)/(rho*T*res(i_chiRho)) ! thermal expansion coefficient
590 :
591 0 : end subroutine eval1_mesa_eosDE_and_kap
592 :
593 :
594 0 : subroutine eval1_kap(s, k, skip_kap, &
595 0 : V, logRho, dlnd_dV, T, logT, species, chem_id, net_iso, xa, &
596 : res, d_dlnd, d_dlnT, OP, OPV, OPT, ierr)
597 : type (star_info), pointer :: s
598 : integer, intent(in) :: k
599 : logical, intent(in) :: skip_kap
600 : integer, intent(in) :: species
601 : integer, pointer :: chem_id(:), net_iso(:)
602 : real(dp), intent(in) :: xa(:)
603 : real(dp), intent(in) :: V, logRho, T, logT, dlnd_dV
604 : real(dp), intent(in), dimension(num_eos_basic_results) :: res, d_dlnd, d_dlnT
605 : real(dp), intent(out) :: OP, OPV, OPT
606 : integer, intent(out) :: ierr
607 :
608 0 : real(dp) :: lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
609 0 : eta, d_eta_dlnRho, d_eta_dlnT, &
610 0 : frac_Type2, opacity, dlnkap_dlnd, dlnkap_dlnT, opacity_factor
611 0 : real(dp) :: kap_fracs(num_kap_fracs), dlnkap_dxa(s% species)
612 :
613 : include 'formats'
614 0 : ierr = 0
615 :
616 0 : if (s% RSP_kap_density_factor > 0d0) then
617 0 : OP = s% RSP_kap_density_factor*V
618 0 : OPV = s% RSP_kap_density_factor
619 0 : OPT = 0d0
620 : return
621 : end if
622 :
623 0 : lnfree_e = res(i_lnfree_e)
624 0 : d_lnfree_e_dlnRho = d_dlnd(i_lnfree_e)
625 0 : d_lnfree_e_dlnT = d_dlnT(i_lnfree_e)
626 :
627 0 : eta = res(i_eta)
628 0 : d_eta_dlnRho = d_dlnd(i_eta)
629 0 : d_eta_dlnT = d_dlnT(i_eta)
630 :
631 0 : if (s% use_other_kap) then
632 : call s% other_kap_get( &
633 : s% id, k, kap_handle, species, chem_id, net_iso, xa, &
634 : logRho, logT, &
635 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
636 : eta, d_eta_dlnRho, d_eta_dlnT, &
637 0 : kap_fracs, opacity, dlnkap_dlnd, dlnkap_dlnT, dlnkap_dxa, ierr)
638 : else
639 : call kap_get( &
640 : kap_handle, species, chem_id, net_iso, xa, &
641 : logRho, logT, &
642 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
643 : eta, d_eta_dlnRho, d_eta_dlnT, &
644 0 : kap_fracs, opacity, dlnkap_dlnd, dlnkap_dlnT, dlnkap_dxa, ierr)
645 : end if
646 0 : frac_Type2 = kap_fracs(i_frac_Type2)
647 0 : if (ierr /= 0) then
648 0 : !$omp critical (rsp_eval_eos_and_kap_2)
649 0 : write(*,*) 'failed in eval1_mesa_eos_and_kap get kap'
650 0 : write(*,2) 'logRho', k, logRho
651 0 : write(*,2) 'logT', k, logT
652 0 : write(*,2) 'lnfree_e', k, lnfree_e
653 0 : write(*,2) 'zbar', k, zbar
654 0 : write(*,2) 'X', k, X
655 0 : write(*,2) 'Z', k, Z
656 0 : call mesa_error(__FILE__,__LINE__,'eval1_kap')
657 : !$omp end critical (rsp_eval_eos_and_kap_2)
658 : !return
659 0 : call mesa_error(__FILE__,__LINE__)
660 : end if
661 :
662 0 : if (k > 0 .and. k <= s% nz .and. s% use_other_opacity_factor) then
663 0 : opacity_factor = s% extra_opacity_factor(k)
664 : else
665 0 : opacity_factor = s% opacity_factor
666 : end if
667 :
668 :
669 0 : if (opacity_factor /= 1d0) then
670 0 : if (s% min_logT_for_opacity_factor_off > 0) then
671 0 : if (logT >= s% max_logT_for_opacity_factor_off .or. &
672 : logT <= s% min_logT_for_opacity_factor_off) then
673 : opacity_factor = 1
674 0 : else if (logT > s% max_logT_for_opacity_factor_on) then
675 : opacity_factor = 1 + (opacity_factor-1)* &
676 : (logT - s% max_logT_for_opacity_factor_off)/ &
677 0 : (s% max_logT_for_opacity_factor_on - s% max_logT_for_opacity_factor_off)
678 0 : else if (logT < s% min_logT_for_opacity_factor_on) then
679 : opacity_factor = 1 + (opacity_factor-1)* &
680 : (logT - s% min_logT_for_opacity_factor_off)/ &
681 0 : (s% min_logT_for_opacity_factor_on - s% min_logT_for_opacity_factor_off)
682 : end if
683 : end if
684 : end if
685 :
686 0 : OP = opacity*opacity_factor
687 0 : OPV = dlnkap_dlnd*opacity*dlnd_dV
688 0 : OPT = dlnkap_dlnT*opacity/T
689 0 : if (k > 0 .and. k <= s% nz) then
690 0 : s% opacity(k) = opacity
691 : end if
692 :
693 0 : end subroutine eval1_kap
694 :
695 :
696 0 : subroutine eval1_mesa_T_given_DP(s, k, Vol, P, T_guess, T, ierr)
697 : use eos_support, only: solve_eos_given_DP
698 : type (star_info), pointer :: s
699 : integer, intent(in) :: k
700 : real(dp), intent(in) :: Vol, P, T_guess
701 : real(dp), intent(out) :: T
702 : integer, intent(out) :: ierr
703 0 : real(dp) :: rho, logRho, logP, logT_guess, &
704 0 : logT_tol, logP_tol, logT
705 : real(dp), dimension(num_eos_basic_results) :: &
706 0 : res, d_dlnd, d_dlnT
707 0 : real(dp) :: d_dxa(num_eos_d_dxa_results, species)
708 : include 'formats'
709 0 : ierr = 0
710 0 : rho = 1d0/Vol
711 0 : logRho = log10(rho)
712 0 : logP = log10(P)
713 0 : logT_guess = log10(T_guess)
714 0 : logT_tol = 1d-11
715 0 : logP_tol = 1d-11
716 : call solve_eos_given_DP( &
717 : s, k, xa, &
718 : logRho, logP, logT_guess, logT_tol, logP_tol, &
719 : logT, res, d_dlnd, d_dlnT, d_dxa, &
720 0 : ierr)
721 0 : T = exp10(logT)
722 0 : end subroutine eval1_mesa_T_given_DP
723 :
724 :
725 0 : subroutine eval1_mesa_Rho_given_PT(s, k, P, T, rho_guess, rho, ierr)
726 0 : use eos_support, only: solve_eos_given_PT
727 : type (star_info), pointer :: s
728 : integer, intent(in) :: k
729 : real(dp), intent(in) :: P, T, rho_guess
730 : real(dp), intent(out) :: rho
731 : integer, intent(out) :: ierr
732 :
733 0 : real(dp) :: logT, logP, logRho_guess, logRho, &
734 0 : logRho_tol, logP_tol
735 : real(dp), dimension(num_eos_basic_results) :: &
736 0 : res, d_dlnd, d_dlnT
737 0 : real(dp) :: d_dxa(num_eos_d_dxa_results, species)
738 : integer :: iter
739 :
740 : include 'formats'
741 :
742 0 : ierr = 0
743 :
744 0 : logT = log10(T)
745 0 : logP = log10(P)
746 0 : logRho_guess = log10(rho_guess)
747 :
748 0 : logRho_tol = 5d-10
749 0 : logP_tol = 5d-13
750 :
751 : ! in some parts of (logRho,logT) cannot meet tight tolerances.
752 : ! so be prepared to relax them.
753 0 : do iter = 1, 9
754 : call solve_eos_given_PT( &
755 : s, k, xa, &
756 : logT, logP, logRho_guess, logRho_tol, logP_tol, &
757 : logRho, res, d_dlnd, d_dlnT, d_dxa, &
758 0 : ierr)
759 0 : if (ierr == 0) exit
760 0 : logRho_tol = logRho_tol*3d0
761 0 : logP_tol = logP_tol*3d0
762 : end do
763 :
764 0 : rho = exp10(logRho)
765 :
766 0 : if (ierr /= 0) then
767 : return
768 : write(*,2) 'Z', k, Z
769 : write(*,2) 'X', k, X
770 : write(*,2) 'abar', k, abar
771 : write(*,2) 'zbar', k, zbar
772 : write(*,2) 'logT', k, logT
773 : write(*,2) 'logP', k, logP
774 : write(*,2) 'logRho_guess', k, logRho_guess
775 : end if
776 :
777 0 : end subroutine eval1_mesa_Rho_given_PT
778 :
779 :
780 0 : real(dp) function eval1_gamma_PT_getRho(s, k, P, T, ierr)
781 0 : use eos_lib, only: eos_gamma_PT_get
782 : type (star_info), pointer :: s
783 : integer, intent(in) :: k
784 : real(dp), intent(in) :: P, T
785 : integer, intent(out) :: ierr
786 : real(dp), dimension(num_eos_basic_results) :: &
787 0 : res, d_dlnd, d_dlnT
788 : real(dp) :: logP, logT, logRho, rho, gamma
789 : include 'formats'
790 0 : logP = log10(P)
791 0 : logT = log10(T)
792 0 : gamma = 5d0/3d0
793 : call eos_gamma_PT_get( &
794 : eos_handle, abar, P, logP, T, logT, gamma, &
795 : rho, logRho, res, d_dlnd, d_dlnT, &
796 0 : ierr)
797 0 : eval1_gamma_PT_getRho = rho
798 0 : end function eval1_gamma_PT_getRho
799 :
800 :
801 0 : subroutine get_surf_P_T_kap(s, &
802 : M, R, L, tau, kap_guess, &
803 : T, P, kap, Teff, ierr)
804 :
805 0 : use atm_support, only: get_atm_PT_legacy_grey_and_kap
806 :
807 : type (star_info), pointer :: s
808 : real(dp), intent(in) :: M, R, L, tau, kap_guess
809 : real(dp), intent(out) :: T, P, kap, Teff
810 : integer, intent(out) :: ierr
811 :
812 : real(dp) :: &
813 : lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
814 : lnP, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap
815 :
816 : ierr = 0
817 :
818 : call get_atm_PT_legacy_grey_and_kap( &
819 : s, tau, L, R, M, s% cgrav(1), .TRUE., &
820 : Teff, kap, &
821 : lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
822 : lnP, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
823 0 : ierr)
824 :
825 0 : T = exp(lnT)
826 0 : P = exp(lnP)
827 :
828 0 : end subroutine get_surf_P_T_kap
829 :
830 :
831 0 : subroutine update_eos_and_kap(s,kk,ierr)
832 : type (star_info), pointer :: s
833 : integer, intent(in) :: kk
834 : integer, intent(out) :: ierr
835 : real(dp) :: &
836 : T,V,Pgas,d_Pg_dV,d_Pg_dT,Prad,d_Pr_dT, &
837 : egas,d_egas_dV,d_egas_dT,erad,d_erad_dV,d_erad_dT, &
838 : CSND,CP,CPV,CPT,Q,QV,QT,OP,OPV,OPT
839 0 : T = s% T(kk)
840 0 : V = 1d0/s% rho(kk)
841 : call eval_mesa_eos_and_kap( &
842 : s,kk,T,V, &
843 : Pgas,d_Pg_dV,d_Pg_dT,Prad,d_Pr_dT, &
844 : egas,d_egas_dV,d_egas_dT,erad,d_erad_dV,d_erad_dT, &
845 0 : CSND,CP,CPV,CPT,Q,QV,QT,OP,OPV,OPT,ierr)
846 0 : end subroutine update_eos_and_kap
847 :
848 :
849 0 : subroutine set_Rho_for_new_Pgas(s,kk,ierr)
850 : type (star_info), pointer :: s
851 : integer, intent(in) :: kk
852 : integer, intent(out) :: ierr
853 : real(dp) :: &
854 0 : other_value, other_tol, logRho_bnd1, other_at_bnd1, &
855 0 : logRho_bnd2, other_at_bnd2, logRho_guess, logRho_result, logRho_tol
856 : real(dp), dimension(num_eos_basic_results) :: &
857 0 : res, d_dlnd, d_dlnT
858 0 : real(dp), dimension(num_eos_d_dxa_results, species) :: d_dxa
859 : integer :: max_iter, which_other, eos_calls
860 :
861 : include 'formats'
862 0 : ierr = 0
863 0 : max_iter = 100
864 0 : which_other = i_lnPgas
865 0 : other_value = log(s% Pgas(kk))
866 0 : other_tol = 1d-11
867 0 : logRho_tol = 1d-11
868 0 : logRho_bnd1 = arg_not_provided
869 0 : other_at_bnd1 = arg_not_provided
870 0 : logRho_bnd2 = arg_not_provided
871 0 : other_at_bnd2 = arg_not_provided
872 0 : logRho_guess = log10(s% rho(kk))
873 0 : call store_T_in_xh(s, kk, s% T(kk))
874 0 : call get_T_and_lnT_from_xh(s, kk, s% T(kk), s% lnT(kk))
875 :
876 : call eosDT_get_Rho( &
877 : eos_handle, &
878 : species, chem_id, net_iso, xa, &
879 : s% lnT(kk)/ln10, which_other, other_value, &
880 : logRho_tol, other_tol, max_iter, logRho_guess, &
881 : logRho_bnd1, logRho_bnd2, other_at_bnd1, other_at_bnd2, &
882 : logRho_result, res, d_dlnd, d_dlnT, &
883 0 : d_dxa, eos_calls, ierr)
884 0 : if (ierr /= 0) return
885 :
886 0 : s% lnd(kk) = logRho_result*ln10
887 0 : s% xh(s% i_lnd,kk) = s% lnd(kk)
888 0 : s% rho(kk) = exp(s% lnd(kk))
889 0 : s% Vol(kk) = 1d0/s% rho(kk)
890 :
891 : end subroutine set_Rho_for_new_Pgas
892 :
893 :
894 0 : subroutine set_T_for_new_egas(s,kk,ierr) ! uses s% T(kk), s% egas(kk) and s% lnd(kk)
895 : type (star_info), pointer :: s
896 : integer, intent(in) :: kk
897 : integer, intent(out) :: ierr
898 :
899 : real(dp) :: &
900 0 : egas_tol, logT_bnd1, egas_at_bnd1, new_egas, egas_want, &
901 0 : logT_bnd2, egas_at_bnd2, logT_guess, logT_result, logT_tol
902 : real(dp), dimension(num_eos_basic_results) :: &
903 0 : res, d_dlnd, d_dlnT
904 0 : real(dp) :: d_dxa(num_eos_d_dxa_results, species)
905 : integer :: max_iter, eos_calls
906 :
907 : include 'formats'
908 0 : ierr = 0
909 0 : max_iter = 100
910 0 : egas_want = s% egas(kk)
911 0 : egas_tol = egas_want*1d-12
912 0 : logT_tol = 1d-11
913 0 : logT_bnd1 = arg_not_provided
914 0 : egas_at_bnd1 = arg_not_provided
915 0 : logT_bnd2 = arg_not_provided
916 0 : egas_at_bnd2 = arg_not_provided
917 0 : logT_guess = log10(s% T(kk))
918 :
919 : call eosDT_get_T( &
920 : eos_handle, &
921 : species, chem_id, net_iso, xa, &
922 : s% lnd(kk)/ln10, i_egas, egas_want, &
923 : logT_tol, egas_tol, max_iter, logT_guess, &
924 : logT_bnd1, logT_bnd2, egas_at_bnd1, egas_at_bnd2, &
925 : logT_result, res, d_dlnd, d_dlnT, &
926 0 : d_dxa, eos_calls, ierr)
927 0 : if (ierr /= 0) return
928 :
929 0 : call store_lnT_in_xh(s, kk, logT_result*ln10)
930 0 : call get_T_and_lnT_from_xh(s, kk, s% T(kk), s% lnT(kk))
931 :
932 0 : new_egas = exp(res(i_lnE)) - crad*s% T(kk)**4/s% rho(kk)
933 0 : if (is_bad(new_egas) .or. new_egas <= 0d0 .or. &
934 : abs(new_egas - egas_want) > 1d3*egas_tol) then
935 0 : write(*,1) 'logRho', s% lnd(kk)/ln10
936 0 : write(*,1) 'logT_guess', logT_guess
937 0 : write(*,1) 'egas_want', egas_want
938 0 : write(*,1) 'Z', Z
939 0 : write(*,1) 'X', X
940 0 : write(*,1) 'abar', abar
941 0 : write(*,1) 'zbar', zbar
942 0 : write(*,1) 'logT_tol', logT_tol
943 0 : write(*,1) 'egas_tol', egas_tol
944 0 : write(*,'(A)')
945 0 : write(*,1) 'guess logT', logT_guess
946 0 : write(*,1) 'found logT', logT_result
947 0 : write(*,1) 'wanted egas', egas_want
948 0 : write(*,1) 'got egas', new_egas
949 0 : write(*,1) '(want - got)/got', (egas_want - new_egas)/new_egas
950 0 : write(*,'(A)')
951 0 : write(*,*) 'eos_calls', eos_calls
952 0 : write(*,'(A)')
953 0 : write(*,2) 'failed set_T_for_new_egas', kk
954 0 : write(*,'(A)')
955 0 : write(*,*) 'is_bad(new_egas)', is_bad(new_egas)
956 0 : write(*,*) 'new_egas <= 0d0', new_egas <= 0d0
957 0 : write(*,*) 'abs(new_egas - egas_want) > egas_tol', &
958 0 : abs(new_egas - egas_want) > egas_tol, &
959 0 : abs(new_egas - egas_want) - egas_tol, &
960 0 : abs(new_egas - egas_want), egas_tol
961 0 : call mesa_error(__FILE__,__LINE__,'set_T_for_new_egas')
962 : end if
963 :
964 : end subroutine set_T_for_new_egas
965 :
966 :
967 0 : subroutine set_T_for_new_Pgas(s,kk,ierr)
968 : type (star_info), pointer :: s
969 : integer, intent(in) :: kk
970 : integer, intent(out) :: ierr
971 : real(dp) :: &
972 0 : other_value, other_tol, logT_bnd1, other_at_bnd1, &
973 0 : logT_bnd2, other_at_bnd2, logT_guess, logT_result, logT_tol
974 : real(dp), dimension(num_eos_basic_results) :: &
975 0 : res, d_dlnd, d_dlnT
976 0 : real(dp), dimension(num_eos_d_dxa_results, species) :: d_dxa
977 : integer :: max_iter, which_other, eos_calls
978 :
979 : include 'formats'
980 0 : ierr = 0
981 0 : max_iter = 100
982 0 : which_other = i_lnPgas
983 0 : other_value = log(s% Pgas(kk))
984 0 : other_tol = 1d-11
985 0 : logT_tol = 1d-11
986 0 : logT_bnd1 = arg_not_provided
987 0 : other_at_bnd1 = arg_not_provided
988 0 : logT_bnd2 = arg_not_provided
989 0 : other_at_bnd2 = arg_not_provided
990 0 : logT_guess = log10(s% T(kk))
991 :
992 : call eosDT_get_T( &
993 : eos_handle, &
994 : species, chem_id, net_iso, xa, &
995 : s% lnd(kk)/ln10, which_other, other_value, &
996 : logT_tol, other_tol, max_iter, logT_guess, &
997 : logT_bnd1, logT_bnd2, other_at_bnd1, other_at_bnd2, &
998 : logT_result, res, d_dlnd, d_dlnT, &
999 0 : d_dxa, eos_calls, ierr)
1000 0 : if (ierr /= 0) return
1001 :
1002 0 : call store_lnT_in_xh(s, kk, logT_result*ln10)
1003 0 : call get_T_and_lnT_from_xh(s, kk, s% T(kk), s% lnT(kk))
1004 :
1005 : end subroutine set_T_for_new_Pgas
1006 :
1007 :
1008 0 : subroutine set_T_for_new_energy(s,kk,logT_tol,other_tol,ierr)
1009 : use eos_lib, only: eosDT_get_T
1010 : type (star_info), pointer :: s
1011 : integer, intent(in) :: kk
1012 : real(dp), intent(in) :: logT_tol, other_tol
1013 : integer, intent(out) :: ierr
1014 : real(dp) :: &
1015 0 : other_value, logT_bnd1, other_at_bnd1, &
1016 0 : logT_bnd2, other_at_bnd2, logT_guess, logT_result
1017 : real(dp), dimension(num_eos_basic_results) :: &
1018 0 : res, d_dlnd, d_dlnT
1019 0 : real(dp), dimension(num_eos_d_dxa_results, species) :: d_dxa
1020 : integer :: max_iter, which_other, eos_calls
1021 0 : ierr = 0
1022 0 : max_iter = 100
1023 0 : which_other = i_lnE
1024 0 : other_value = log(s% energy(kk))
1025 : !other_tol = 1d-12
1026 : !logT_tol = 1d-12
1027 0 : logT_bnd1 = arg_not_provided
1028 0 : other_at_bnd1 = arg_not_provided
1029 0 : logT_bnd2 = arg_not_provided
1030 0 : other_at_bnd2 = arg_not_provided
1031 0 : logT_guess = log10(s% T(kk))
1032 :
1033 : call eosDT_get_T( &
1034 : eos_handle, &
1035 : species, chem_id, net_iso, xa, &
1036 : s% lnd(kk)/ln10, which_other, other_value, &
1037 : logT_tol, other_tol, max_iter, logT_guess, &
1038 : logT_bnd1, logT_bnd2, other_at_bnd1, other_at_bnd2, &
1039 : logT_result, res, d_dlnd, d_dlnT, &
1040 0 : d_dxa, eos_calls, ierr)
1041 0 : if (ierr /= 0) return
1042 :
1043 0 : call store_lnT_in_xh(s, kk, logT_result*ln10)
1044 0 : call get_T_and_lnT_from_xh(s, kk, s% T(kk), s% lnT(kk))
1045 0 : end subroutine set_T_for_new_energy
1046 :
1047 : end module rsp_eval_eos_and_kap
|
