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 profile_getval
21 :
22 : use star_private_def
23 : use star_profile_def
24 : use const_def, only: dp, qe, kerg, avo, amu, boltz_sigma, secday, secyer, standard_cgrav, &
25 : clight, four_thirds_pi, ln10, lsun, msun, pi, pi4, rsun, sqrt_2_div_3, one_third, &
26 : convective_mixing, &
27 : overshoot_mixing, &
28 : semiconvective_mixing, &
29 : thermohaline_mixing, &
30 : minimum_mixing, &
31 : anonymous_mixing, &
32 : leftover_convective_mixing
33 : use star_utils
34 : use utils_lib
35 : use auto_diff_support, only: get_w, get_etrb
36 :
37 : implicit none
38 :
39 : integer, parameter :: idel = 10000
40 : integer, parameter :: add_abundances = idel
41 : integer, parameter :: add_log_abundances = add_abundances + 1
42 : integer, parameter :: category_offset = add_log_abundances + 1
43 : integer, parameter :: abundance_offset = category_offset + idel
44 : integer, parameter :: log_abundance_offset = abundance_offset + idel
45 : integer, parameter :: xadot_offset = log_abundance_offset + idel
46 : integer, parameter :: xaprev_offset = xadot_offset + idel
47 : integer, parameter :: ionization_offset = xaprev_offset + idel
48 : integer, parameter :: typical_charge_offset = ionization_offset + idel
49 : integer, parameter :: edv_offset = typical_charge_offset + idel
50 : integer, parameter :: extra_diffusion_factor_offset = edv_offset + idel
51 : integer, parameter :: v_rad_offset = extra_diffusion_factor_offset + idel
52 : integer, parameter :: log_g_rad_offset = v_rad_offset + idel
53 : integer, parameter :: log_concentration_offset = log_g_rad_offset + idel
54 : integer, parameter :: diffusion_dX_offset = log_concentration_offset + idel
55 : integer, parameter :: diffusion_D_offset = diffusion_dX_offset + idel
56 : integer, parameter :: raw_rate_offset = diffusion_D_offset + idel
57 : integer, parameter :: screened_rate_offset = raw_rate_offset + idel
58 : integer, parameter :: eps_nuc_rate_offset = screened_rate_offset + idel
59 : integer, parameter :: eps_neu_rate_offset = eps_nuc_rate_offset + idel
60 : integer, parameter :: extra_offset = eps_neu_rate_offset + idel
61 : integer, parameter :: max_profile_offset = extra_offset + idel
62 :
63 : contains
64 :
65 12 : integer function do1_profile_spec( &
66 : s, iounit, n, i, string, buffer, report, ierr) result(spec)
67 :
68 : use utils_lib
69 : use utils_def
70 : use chem_def
71 : use chem_lib
72 : use net_def
73 :
74 : type(star_info), pointer :: s
75 : integer :: iounit, n, i, num, t
76 :
77 : character (len=*) :: string, buffer
78 : logical, intent(in) :: report
79 : integer, intent(out) :: ierr
80 :
81 : integer :: id
82 : type(Net_General_Info), pointer :: g
83 :
84 : ierr = 0
85 12 : spec = -1
86 :
87 12 : call get_net_ptr(s% net_handle, g, ierr)
88 12 : if(ierr/=0) return
89 :
90 12 : id = do_get_profile_id(string)
91 12 : if (id > 0) then
92 : spec = id
93 : return
94 : end if
95 :
96 0 : select case(string)
97 :
98 : case ('xadot')
99 0 : call do1_nuclide(xadot_offset)
100 :
101 : case ('xaprev')
102 0 : call do1_nuclide(xaprev_offset)
103 :
104 : case ('ionization')
105 0 : call do1_nuclide(ionization_offset)
106 :
107 : case ('typical_charge')
108 0 : call do1_nuclide(typical_charge_offset)
109 :
110 : case ('edv')
111 0 : call do1_nuclide(edv_offset)
112 :
113 : case ('extra_diffusion_factor')
114 0 : call do1_nuclide(extra_diffusion_factor_offset)
115 :
116 : case ('v_rad')
117 0 : call do1_nuclide(v_rad_offset)
118 :
119 : case ('log_g_rad')
120 0 : call do1_nuclide(log_g_rad_offset)
121 :
122 : case ('log_concentration')
123 0 : call do1_nuclide(log_concentration_offset)
124 :
125 : case ('diffusion_dX')
126 0 : call do1_nuclide(diffusion_dX_offset)
127 :
128 : case ('diffusion_D')
129 0 : call do1_nuclide(diffusion_D_offset)
130 :
131 : case ('log') ! add log of abundance
132 0 : call do1_nuclide(log_abundance_offset)
133 :
134 : case ('eps_neu_rate')
135 0 : call do1_rate(eps_neu_rate_offset)
136 :
137 : case ('eps_nuc_rate')
138 0 : call do1_rate(eps_nuc_rate_offset)
139 :
140 : case ('screened_rate')
141 0 : call do1_rate(screened_rate_offset)
142 :
143 : case ('raw_rate')
144 0 : call do1_rate(raw_rate_offset)
145 :
146 : case ('extra')
147 :
148 0 : t = token(iounit, n, i, buffer, string)
149 0 : if (t /= name_token) then
150 0 : ierr = -1; return
151 : end if
152 0 : read(string,fmt=*,iostat=ierr) num
153 0 : if (ierr /= 0 .or. num <= 0 .or. num > max_num_profile_extras) then
154 0 : write(*,*) 'failed to find valid integer for extra: ' // trim(string)
155 0 : ierr = -1
156 : end if
157 0 : spec = extra_offset + num
158 :
159 : case default
160 :
161 3 : id = chem_get_iso_id(string)
162 3 : if (id > 0) then
163 0 : spec = abundance_offset + id
164 0 : return
165 : end if
166 3 : id = rates_category_id(string)
167 3 : if (id > 0) then
168 3 : spec = category_offset + id
169 3 : return
170 : end if
171 0 : if (report) &
172 0 : write(*,*) 'failed to recognize item for profile columns: ' // trim(string)
173 3 : ierr = -1
174 :
175 : end select
176 :
177 :
178 : contains
179 :
180 :
181 0 : subroutine do1_nuclide(offset)
182 : integer, intent(in) :: offset
183 : integer :: t, id
184 0 : t = token(iounit, n, i, buffer, string)
185 0 : if (t /= name_token) then
186 0 : ierr = -1; return
187 : end if
188 0 : id = chem_get_iso_id(string)
189 0 : if (id > 0) then
190 0 : spec = offset + id
191 0 : return
192 : end if
193 0 : write(*,*) 'bad iso name: ' // trim(string)
194 0 : ierr = -1
195 : end subroutine do1_nuclide
196 :
197 0 : subroutine do1_rate(offset) ! raw_rate, screened_rate, eps_nuc_rate, eps_neu_rate
198 : use rates_lib, only: rates_reaction_id
199 : integer, intent(in) :: offset
200 : integer :: t, id
201 0 : t = token(iounit, n, i, buffer, string)
202 0 : if (t /= name_token) then
203 0 : ierr = -1; return
204 : end if
205 0 : id = rates_reaction_id(string)
206 0 : id = g% net_reaction(id) ! Convert to net id not the global rate id
207 0 : if (id > 0) then
208 0 : spec = offset + id
209 0 : return
210 : end if
211 0 : write(*,*) 'bad rate name: ' // trim(string)
212 0 : ierr = -1
213 : end subroutine do1_rate
214 :
215 :
216 : end function do1_profile_spec
217 :
218 :
219 0 : integer function get_profile_id(s, name) result(spec)
220 : use utils_lib, only: token
221 : use utils_def, only: name_token
222 : type (star_info), pointer :: s
223 : character (len=*), intent(in) :: name
224 : character (len=strlen) :: buffer, string
225 : integer :: i, n, iounit, ierr, t
226 0 : iounit = -1
227 0 : ierr = 0
228 0 : buffer = name
229 0 : n = len_trim(buffer) + 1
230 0 : buffer(n:n) = ' '
231 0 : i = 0
232 0 : t = token(iounit, n, i, buffer, string)
233 0 : if (t /= name_token) then
234 0 : spec = -1; return
235 : end if
236 0 : spec = do1_profile_spec(s, iounit, n, i, string, buffer, .false., ierr)
237 0 : if (ierr == 0) return
238 : ! check to see if it is one of the extra profile columns
239 0 : do i=1,s% num_extra_profile_cols
240 0 : if (name == s% extra_profile_col_names(i)) then
241 0 : spec = i + max_profile_offset
242 0 : return
243 : end if
244 : end do
245 0 : spec = -1
246 0 : end function get_profile_id
247 :
248 :
249 0 : real(dp) function get_profile_val(s, id, k)
250 : type (star_info), pointer :: s
251 : integer, intent(in) :: id, k
252 : integer :: int_val
253 : logical :: int_flag
254 0 : if (id > max_profile_offset) then ! get from extras
255 0 : get_profile_val = s% extra_profile_col_vals(k, id - max_profile_offset)
256 0 : return
257 : end if
258 0 : call getval_for_profile(s, id, k, get_profile_val, int_flag, int_val)
259 0 : if (int_flag) get_profile_val = dble(int_val)
260 0 : end function get_profile_val
261 :
262 :
263 41472 : subroutine getval_for_profile(s, c, k, val, int_flag, int_val)
264 : use chem_def
265 : use rates_def
266 : use ionization_def
267 : use mod_typical_charge, only: eval_typical_charge
268 : use rsp_def, only: rsp_WORK, rsp_WORKQ, rsp_WORKT, rsp_WORKC
269 :
270 : type (star_info), pointer :: s
271 : integer, intent(in) :: c, k
272 : real(dp), intent(out) :: val
273 : integer, intent(out) :: int_val
274 : logical, intent(inout) :: int_flag
275 :
276 : real(dp) :: cno, z, x, L_rad, L_edd, &
277 : r00, rp1, v00, vp1, Ap1, &
278 : r00_start, rp1_start, dr3, dr3_start, &
279 : d_dlnR00, d_dlnRp1, d_dv00, d_dvp1
280 : integer :: j, nz, ionization_k, klo, khi, i, ii, ierr
281 : real(dp) :: f, lgT, full_on, full_off, am_nu_factor
282 : logical :: rsp_or_w, reconstructed_face_state_active
283 : include 'formats'
284 :
285 41472 : if (s% rotation_flag) then
286 0 : full_on = s% D_mix_rotation_max_logT_full_on
287 0 : full_off = s% D_mix_rotation_min_logT_full_off
288 0 : lgT = s% lnT(k)/ln10
289 0 : if (lgT <= full_on) then
290 : f = 1d0
291 0 : else if (lgT >= full_off) then
292 : f = 0d0
293 : else ! lgT > full_on and < full_off
294 0 : f = (lgT - full_on) / (full_off - full_on)
295 : end if
296 0 : am_nu_factor = f*s% am_nu_factor
297 : else
298 : am_nu_factor = 1d0
299 : end if
300 :
301 41472 : val = 0; int_val = 0; int_flag = .false.
302 41472 : nz = s% nz
303 41472 : ionization_k = 0
304 :
305 : int_flag = .false.
306 41472 : rsp_or_w = s% RSP_flag .or. s% RSP2_flag
307 41472 : reconstructed_face_state_active = s% use_face_reconstruction
308 :
309 41472 : if (c > extra_offset) then
310 0 : i = c - extra_offset
311 0 : val = s% profile_extra(k,i)
312 : ! TODO: implement eps_neu_rate, eps_nuc_rate, screened_rate
313 41472 : else if (c > eps_neu_rate_offset) then
314 0 : i = c - eps_neu_rate_offset
315 0 : val = s% eps_neu_rate(i,k) * s% dm(k)
316 41472 : else if (c > eps_nuc_rate_offset) then
317 0 : i = c - eps_nuc_rate_offset
318 0 : val = s% eps_nuc_rate(i,k) * s% dm(k)
319 41472 : else if (c > screened_rate_offset) then
320 0 : i = c - screened_rate_offset
321 0 : val = s% screened_rate(i,k) * s% dm(k)
322 41472 : else if (c > raw_rate_offset) then
323 0 : i = c - raw_rate_offset
324 0 : val = s% raw_rate(i,k) * s% dm(k)
325 41472 : else if (c > diffusion_D_offset) then
326 0 : i = c - diffusion_D_offset
327 0 : ii = s% net_iso(i)
328 0 : if (ii > 0 .and. s% do_element_diffusion) val = s% diffusion_D_self(ii,k)
329 41472 : else if (c > diffusion_dX_offset) then
330 0 : i = c - diffusion_dX_offset
331 0 : ii = s% net_iso(i)
332 0 : if (ii > 0 .and. s% do_element_diffusion) val = s% diffusion_dX(ii,k)
333 41472 : else if (c > log_concentration_offset) then
334 0 : i = c - log_concentration_offset
335 0 : ii = s% net_iso(i)
336 0 : if (ii > 0) val = get_log_concentration(s,ii,k)
337 41472 : else if (c > log_g_rad_offset) then
338 0 : i = c - log_g_rad_offset
339 0 : ii = s% net_iso(i)
340 0 : if (ii > 0 .and. s% do_element_diffusion) val = safe_log10(s% g_rad(ii,k))
341 41472 : else if (c > v_rad_offset) then
342 0 : i = c - v_rad_offset
343 0 : ii = s% net_iso(i)
344 0 : if (ii > 0 .and. s% do_element_diffusion) val = s% v_rad(ii,k)
345 41472 : else if (c > extra_diffusion_factor_offset) then
346 0 : i = c - extra_diffusion_factor_offset
347 0 : ii = s% net_iso(i)
348 0 : if (ii > 0 .and. s% do_element_diffusion) val = s% extra_diffusion_factor(ii,k)
349 41472 : else if (c > edv_offset) then
350 0 : i = c - edv_offset
351 0 : ii = s% net_iso(i)
352 0 : if (ii > 0 .and. s% do_element_diffusion) val = s% edv(ii,k)
353 41472 : else if (c > typical_charge_offset) then
354 0 : i = c - typical_charge_offset
355 0 : ii = s% net_iso(i)
356 0 : if (ii > 0 .and. s% do_element_diffusion) val = s% typical_charge(ii,k)
357 41472 : else if (c > ionization_offset) then
358 0 : i = c - ionization_offset
359 0 : ii = s% net_iso(i)
360 : val = eval_typical_charge( &
361 : i, s% abar(k), exp(s% lnfree_e(k)), &
362 0 : s% T(k), s% lnT(k)/ln10, s% rho(k), s% lnd(k)/ln10)
363 41472 : else if (c > xaprev_offset) then
364 0 : i = c - xaprev_offset
365 0 : ii = s% net_iso(i)
366 0 : if (ii > 0) val = s% xa_start(ii,k)
367 41472 : else if (c > xadot_offset) then
368 0 : i = c - xadot_offset
369 0 : ii = s% net_iso(i)
370 0 : if (ii > 0) val = s% xa(ii,k) - s% xa_start(ii,k)
371 41472 : else if (c > log_abundance_offset) then
372 0 : i = c - log_abundance_offset
373 0 : ii = s% net_iso(i)
374 0 : if (ii > 0) then
375 0 : val = safe_log10(s% xa(ii,k))
376 : else
377 0 : val = -99d0
378 : end if
379 41472 : else if (c > abundance_offset) then
380 0 : i = c - abundance_offset
381 0 : ii = s% net_iso(i)
382 0 : if (ii > 0) val = s% xa(ii,k)
383 41472 : else if (c > category_offset) then
384 10368 : i = c - category_offset
385 10368 : val = s% eps_nuc_categories(i,k)
386 : else
387 :
388 3456 : select case(c)
389 : case (p_zone)
390 3456 : val = dble(k)
391 3456 : int_val = k
392 3456 : int_flag = .true.
393 : case (p_k)
394 0 : val = dble(k)
395 0 : int_val = k
396 0 : int_flag = .true.
397 : case (p_conv_L_div_L)
398 0 : if (s% L(k) > 0d0) val = get_Lconv(s,k)/s% L(k)
399 : case (p_log_conv_L_div_L)
400 0 : if (s% L(k) > 0d0) val = safe_log10(get_Lconv(s,k)/s% L(k))
401 : case (p_lum_erg_s)
402 0 : val = s% L(k)
403 : case (p_L)
404 0 : val = s% L(k)/Lsun
405 : case (p_luminosity)
406 0 : val = s% L(k)/Lsun
407 : case (p_log_abs_lum_erg_s)
408 0 : val = safe_log10(abs(s% L(k)))
409 : case (p_lum_adv)
410 0 : val = get_Ladv(s,k)
411 : case (p_lum_plus_lum_adv)
412 0 : val = s% L(k) + get_Ladv(s,k)
413 : case (p_lum_rad)
414 0 : L_rad = get_Lrad(s,k)
415 0 : val = L_rad/Lsun
416 : case (p_lum_conv)
417 0 : val = get_Lconv(s,k)/Lsun
418 : case (p_lum_conv_MLT)
419 0 : val = s% L_conv(k)/Lsun
420 :
421 : case(p_Frad_div_cUrad)
422 : val = ((s% L(k) - s% L_conv(k)) / (4._dp*pi*pow2(s%r(k)))) &
423 0 : /(clight * s% Prad(k) *3._dp)
424 : case (p_lum_rad_div_L_Edd_sub_fourPrad_div_PchiT)
425 0 : val = get_Lrad_div_Ledd(s,k) - 4*s% Prad(k)/(s% Peos(k)*s% chiT(k))
426 : case (p_lum_rad_div_L_Edd)
427 0 : val = get_Lrad_div_Ledd(s,k)
428 : case (p_lum_conv_div_lum_Edd)
429 0 : L_rad = get_Lrad(s,k)
430 0 : L_edd = get_Ledd(s,k)
431 0 : val = (s% L(k) - L_rad)/L_edd
432 :
433 : case (p_lum_conv_div_lum_rad)
434 0 : L_rad = get_Lrad(s,k)
435 0 : val = (s% L(k) - L_rad)/L_rad
436 :
437 : case (p_lum_rad_div_L)
438 0 : L_rad = get_Lrad(s,k)
439 0 : val = L_rad/max(1d0,s% L(k))
440 : case (p_lum_conv_div_L)
441 0 : L_rad = get_Lrad(s,k)
442 0 : val = (s% L(k) - L_rad)/max(1d0,s% L(k))
443 :
444 : case (p_log_Lrad)
445 0 : L_rad = get_Lrad(s,k)
446 0 : val = safe_log10(L_rad/Lsun)
447 : case (p_log_Lconv)
448 0 : L_rad = get_Lrad(s,k)
449 0 : val = safe_log10((s% L(k) - L_rad)/Lsun)
450 : case (p_log_Lconv_div_L)
451 0 : L_rad = get_Lrad(s,k)
452 0 : val = safe_log10((s% L(k) - L_rad)/s% L(k))
453 :
454 : case (p_log_Lrad_div_L)
455 0 : L_rad = get_Lrad(s,k)
456 0 : val = safe_log10(L_rad/s% L(k))
457 : case (p_log_Lrad_div_Ledd)
458 0 : val = safe_log10(get_Lrad_div_Ledd(s,k))
459 :
460 : case (p_log_g)
461 0 : val = safe_log10(s% grav(k))
462 : case (p_grav)
463 0 : val = s% grav(k)
464 : case (p_g_div_r)
465 0 : val = s% grav(k)/s% r(k)
466 : case (p_r_div_g)
467 0 : val = s% r(k)/s% grav(k)
468 : case (p_signed_log_eps_grav)
469 0 : val = s% eps_grav_ad(k)% val
470 0 : val = sign(1d0,val)*log10(max(1d0,abs(val)))
471 : case (p_net_nuclear_energy)
472 : ! Do not subtract s% eps_nuc_neu_total(k) eps_nuc already contains it
473 0 : val = s% eps_nuc(k) - s% non_nuc_neu(k)
474 0 : val = sign(1d0,val)*log10(max(1d0,abs(val)))
475 : case (p_eps_nuc_plus_nuc_neu)
476 : ! eps_nuc subtracts eps_nuc_neu so this is just the total eenrgy from nuclear burning without neutrinos
477 0 : val = s% eps_nuc(k) + s% eps_nuc_neu_total(k)
478 : case (p_eps_nuc_minus_non_nuc_neu)
479 0 : val = s% eps_nuc(k) - s% non_nuc_neu(k)
480 : case (p_net_energy)
481 0 : val = s% eps_nuc(k) - s% non_nuc_neu(k) + s% eps_grav_ad(k)% val
482 0 : val = sign(1d0,val)*log10(max(1d0,abs(val)))
483 : case (p_signed_log_power)
484 0 : val = s% L(k)
485 0 : val = sign(1d0,val)*log10(max(1d0,abs(val)))
486 : case (p_logL)
487 0 : val = safe_log10(max(1d-12,s% L(k)/Lsun))
488 : case (p_log_Ledd)
489 0 : val = safe_log10(get_Ledd(s,k)/Lsun)
490 : case (p_lum_div_Ledd)
491 0 : val = s% L(k)/get_Ledd(s,k)
492 : case (p_log_L_div_Ledd)
493 0 : val = safe_log10(max(1d-12,s% L(k)/get_Ledd(s,k)))
494 : case (p_log_abs_v)
495 0 : if (s% u_flag) then
496 0 : val = safe_log10(abs(s% u(k)))
497 0 : else if (s% v_flag) then
498 0 : val = safe_log10(abs(s% v(k)))
499 : end if
500 :
501 : case (p_superad_reduction_factor)
502 0 : val = s% superad_reduction_factor(k)
503 : case (p_gradT_excess_effect)
504 0 : val = s% gradT_excess_effect(k)
505 : case (p_diff_grads)
506 0 : val = s% gradr(k) - s% gradL(k) ! convective if this is > 0
507 : case (p_log_diff_grads)
508 0 : val = safe_log10(abs(s% gradr(k) - s% gradL(k)))
509 : case (p_v)
510 0 : if (s% u_flag) then
511 0 : val = s% u(k)
512 0 : else if (s% v_flag) then
513 0 : val = s% v(k)
514 : end if
515 : case (p_velocity)
516 0 : if (s% u_flag) then
517 0 : val = s% u(k)
518 0 : else if (s% v_flag) then
519 0 : val = s% v(k)
520 : end if
521 : case (p_v_kms)
522 0 : if (s% u_flag) then
523 0 : val = s% u(k)*1d-5
524 0 : else if (s% v_flag) then
525 0 : val = s% v(k)*1d-5
526 : end if
527 : case (p_vel_km_per_s)
528 0 : if (s% u_flag) then
529 0 : val = s% u(k)*1d-5
530 0 : else if (s% v_flag) then
531 0 : val = s% v(k)*1d-5
532 : end if
533 : case (p_v_div_r)
534 0 : if (s% u_flag) then
535 0 : val = s% u_face_ad(k)%val/s% r(k)
536 0 : else if (s% v_flag) then
537 0 : val = s% v(k)/s% r(k)
538 : end if
539 :
540 : case (p_v_times_t_div_r)
541 0 : if (s% u_flag) then
542 0 : val = s% u_face_ad(k)%val*s% time/s% r(k)
543 0 : else if (s% v_flag) then
544 0 : val = s% v(k)*s% time/s% r(k)
545 : end if
546 : case (p_radius)
547 0 : val = s% r(k)/Rsun
548 : case (p_radius_cm)
549 0 : val = s% r(k)
550 : case (p_radius_km)
551 0 : val = s% r(k)*1d-5
552 : case (p_rmid)
553 0 : val = s% rmid(k)/Rsun
554 : case (p_logR_cm)
555 0 : val = safe_log10(s% r(k))
556 : case (p_logR)
557 3456 : val = safe_log10(s% r(k)/Rsun)
558 : case (p_psi_roche)
559 0 : if (.not. associated(s% binary_get_roche_potential)) then
560 0 : val = -99d0
561 : else
562 0 : call s% binary_get_roche_potential(s% id, s% r(k), val, ierr)
563 : end if
564 :
565 : case (p_q)
566 0 : val = s% q(k)
567 : case (p_log_q)
568 0 : val = safe_log10(s% q(k))
569 : case (p_dq)
570 0 : val = s% dq(k)
571 : case (p_log_dq)
572 0 : val = safe_log10(s% dq(k))
573 : case (p_mass)
574 3456 : val = s% m(k)/Msun
575 : case (p_log_mass)
576 0 : val = safe_log10(s% m(k)/Msun)
577 : case (p_mass_grams)
578 0 : val = s% m(k)
579 : case (p_mmid)
580 0 : val = (s% M_center + s% xmstar*(s% q(k) - s% dq(k)/2))/Msun
581 :
582 : case (p_dm)
583 0 : val = s% dm(k)
584 : case (p_dm_bar)
585 0 : val = s% dm_bar(k)
586 :
587 : case (p_m_div_r)
588 0 : val = s% m(k)/s% r(k)
589 : case (p_dmbar_m_div_r)
590 0 : val = s% dm_bar(k)*s% m(k)/s% r(k)
591 : case (p_log_dmbar_m_div_r)
592 0 : val = safe_log10(s% dm_bar(k)*s% m(k)/s% r(k))
593 :
594 : case (p_m_grav)
595 0 : val = s% m_grav(k)/Msun
596 : case (p_m_grav_div_m_baryonic)
597 0 : val = s% m_grav(k)/s% m(k)
598 : case (p_mass_correction_factor)
599 0 : val = s% mass_correction(k)
600 :
601 : case (p_xr)
602 0 : val = (s% r(1) - s% r(k))/Rsun
603 : case (p_xr_cm)
604 0 : val = s% r(1) - s% r(k)
605 : case (p_xr_div_R)
606 0 : val = (s% r(1) - s% r(k))/s% r(1)
607 : case (p_log_xr)
608 0 : val = safe_log10((s% r(1) - s% r(k))/Rsun)
609 : case (p_log_xr_cm)
610 0 : val = safe_log10(s% r(1) - s% r(k))
611 : case (p_log_xr_div_R)
612 0 : val = safe_log10((s% r(1) - s% r(k))/s% r(1))
613 :
614 : case (p_x)
615 0 : val = s% X(k)
616 : case (p_log_x)
617 0 : val = safe_log10(s% X(k))
618 : case (p_y)
619 0 : val = s% Y(k)
620 : case (p_log_y)
621 0 : val = safe_log10(s% Y(k))
622 : case (p_z)
623 0 : val = s% Z(k)
624 : case (p_log_z)
625 0 : val = safe_log10(s% Z(k))
626 : case (p_xm)
627 0 : val = sum(s% dm(1:k-1))/Msun
628 : case (p_logxm)
629 0 : val = safe_log10(sum(s% dm(1:k-1))/Msun)
630 : case (p_xq)
631 0 : val = sum(s% dq(1:k-1))
632 : case (p_logxq)
633 0 : val = safe_log10(sum(s% dq(1:k-1)))
634 : case (p_logdq)
635 0 : val = safe_log10(s% dq(k))
636 : case (p_log_column_depth)
637 0 : val = safe_log10(s% xmstar*sum(s% dq(1:k-1))/(pi4*s% r(k)*s% r(k)))
638 : case (p_log_radial_depth)
639 0 : val = safe_log10(s% r(1) - s% r(k))
640 :
641 : case (p_r_div_R)
642 0 : val = s% r(k)/s% r(1)
643 : case (p_log_dr)
644 0 : if (k == s% nz) then
645 0 : val = s% r(k) - s% R_center
646 : else
647 0 : val = s% r(k) - s% r(k+1)
648 : end if
649 0 : val = safe_log10(val)
650 : case (p_dlogR)
651 0 : if (k == s% nz) then
652 0 : val = s% lnR(k) - log(max(1d0,s% R_center))
653 : else
654 0 : val = s% lnR(k) - s% lnR(k+1)
655 : end if
656 0 : val = val/ln10
657 : case (p_dr_div_rmid)
658 0 : if (k == s% nz) then
659 0 : val = s% r(k) - s% R_center
660 : else
661 0 : val = s% r(k) - s% r(k+1)
662 : end if
663 0 : val = val/s% rmid(k)
664 : case (p_log_dr_div_rmid)
665 0 : if (k == s% nz) then
666 0 : val = s% r(k) - s% R_center
667 : else
668 0 : val = s% r(k) - s% r(k+1)
669 : end if
670 0 : val = safe_log10(val/s% rmid(k))
671 : case (p_log_acoustic_radius)
672 0 : val = safe_log10(sum(s% dr_div_csound(k:nz)))
673 : case (p_acoustic_radius)
674 0 : val = sum(s% dr_div_csound(k:nz))
675 : case (p_log_acoustic_depth)
676 0 : if (k > 1) &
677 0 : val = sum(s% dr_div_csound(1:k-1))
678 0 : val = safe_log10(val)
679 : case (p_acoustic_depth)
680 0 : if (k > 1) &
681 0 : val = sum(s% dr_div_csound(1:k-1))
682 : case (p_acoustic_r_div_R_phot)
683 0 : val = sum(s% dr_div_csound(k:nz))/s% photosphere_acoustic_r
684 :
685 : case (p_ergs_error)
686 0 : val = s% ergs_error(k)
687 : case (p_log_rel_E_err)
688 0 : val = safe_log10(abs(s% ergs_error(k)/s% total_energy_start))
689 : case (p_ergs_error_integral)
690 0 : val = sum(s% ergs_error(1:k))
691 : case (p_ergs_rel_error_integral)
692 0 : if (s% total_energy_end /= 0d0) &
693 0 : val = sum(s% ergs_error(1:k))/s% total_energy_end
694 :
695 : case (p_cell_internal_energy_fraction)
696 0 : val = s% energy(k)*s% dm(k)/s% total_internal_energy_end
697 : case (p_cell_internal_energy_fraction_start)
698 0 : val = s% energy_start(k)*s% dm(k)/s% total_internal_energy_start
699 :
700 : case (p_dr_div_R)
701 0 : if (k < s% nz) then
702 0 : val = (s% r(k) - s% r(k+1))/s% r(1)
703 : else
704 0 : val = (s% r(k) - s% r_center)/s% r(1)
705 : end if
706 : case (p_dRstar_div_dr)
707 0 : if (k < s% nz) then
708 0 : val = (s% r(1) - s% R_center)/(s% r(k) - s% r(k+1))
709 : else
710 0 : val = (s% r(1) - s% R_center)/(s% r(k) - s% r_center)
711 : end if
712 : case (p_log_dr_div_R)
713 0 : if (k < s% nz) then
714 0 : val = (s% r(k) - s% r(k+1))/s% r(1)
715 : else
716 0 : val = (s% r(k) - s% r_center)/s% r(1)
717 : end if
718 0 : val = safe_log10(val)
719 :
720 : case(p_t_rad)
721 0 : val = 1d0/(clight*s% opacity(k)*s% rho(k))
722 : case(p_log_t_rad)
723 0 : val = log10(1d0/(clight*s% opacity(k)*s% rho(k)))
724 : case (p_dt_cs_div_dr)
725 0 : if (k < s% nz) then
726 0 : val = s% r(k) - s% r(k+1)
727 : else
728 0 : val = s% r(k) - s% r_center
729 : end if
730 0 : val = s% dt*s% csound(k)/val
731 : case (p_log_dt_cs_div_dr)
732 0 : if (k < s% nz) then
733 0 : val = s% r(k) - s% r(k+1)
734 : else
735 0 : val = s% r(k) - s% r_center
736 : end if
737 0 : val = s% dt*s% csound(k)/val
738 0 : val = safe_log10(val)
739 : case (p_dr_div_cs)
740 0 : if (k == s% nz) then
741 0 : val = s% r(k) - s% R_center
742 : else
743 0 : val = s% r(k) - s% r(k+1)
744 : end if
745 0 : val = val/s% csound(k)
746 : case (p_log_dr_div_cs)
747 0 : if (k == s% nz) then
748 0 : val = s% r(k) - s% R_center
749 : else
750 0 : val = s% r(k) - s% r(k+1)
751 : end if
752 0 : val = safe_log10(val/s% csound(k))
753 :
754 : case (p_dr_div_cs_yr)
755 0 : if (k == s% nz) then
756 0 : val = s% r(k) - s% R_center
757 : else
758 0 : val = s% r(k) - s% r(k+1)
759 : end if
760 0 : val = val/s% csound(k)/secyer
761 : case (p_log_dr_div_cs_yr)
762 0 : if (k == s% nz) then
763 0 : val = s% r(k) - s% R_center
764 : else
765 0 : val = s% r(k) - s% r(k+1)
766 : end if
767 0 : val = safe_log10(val/s% csound(k)/secyer)
768 :
769 : case (p_pgas_div_ptotal)
770 0 : val = s% Pgas(k)/s% Peos(k)
771 : case (p_prad_div_pgas)
772 0 : val = s% Prad(k)/s% Pgas(k)
773 : case(p_prad_div_pgas_div_L_div_Ledd)
774 0 : val = (s% Prad(k)/s% Pgas(k))/max(1d-12,s% L(k)/get_Ledd(s,k))
775 : case (p_pgas_div_p)
776 0 : val = s% Pgas(k)/s% Peos(k)
777 : case (p_flux_limit_R)
778 0 : if (s% use_dPrad_dm_form_of_T_gradient_eqn .and. k > 1) &
779 0 : val = s% flux_limit_R(k)
780 : case (p_flux_limit_lambda)
781 0 : if (s% use_dPrad_dm_form_of_T_gradient_eqn .and. k > 1) then
782 0 : val = s% flux_limit_lambda(k)
783 : else
784 0 : val = 1d0
785 : end if
786 : case (p_cell_collapse_time)
787 0 : if (s% v_flag) then
788 0 : if (k == s% nz) then
789 0 : rp1 = s% R_center
790 0 : vp1 = s% v_center
791 : else
792 0 : rp1 = s% r(k+1)
793 0 : vp1 = s% v(k+1)
794 : end if
795 0 : r00 = s% r(k)
796 0 : v00 = s% v(k)
797 0 : if (vp1 > v00) val = (r00 - rp1)/(vp1 - v00)
798 : end if
799 :
800 : case (p_log_cell_collapse_time)
801 0 : if (s% v_flag) then
802 0 : if (k == s% nz) then
803 0 : rp1 = s% R_center
804 0 : vp1 = s% v_center
805 : else
806 0 : rp1 = s% r(k+1)
807 0 : vp1 = s% v(k+1)
808 : end if
809 0 : r00 = s% r(k)
810 0 : v00 = s% v(k)
811 0 : if (vp1 > v00) val = (r00 - rp1)/(vp1 - v00)
812 : end if
813 0 : val = safe_log10(val)
814 :
815 : case (p_dq_ratio)
816 0 : if (k == 1 .or. k == s% nz) then
817 0 : val = 1
818 : else
819 0 : val = s% dq(k-1)/s% dq(k)
820 : end if
821 :
822 : case (p_compression_gradient)
823 0 : if (k == 1) then
824 : val = s% rho_start(1)*&
825 0 : ((1/s% rho(1) - 1/s% rho_start(1)) - (1/s% rho(2) - 1/s% rho_start(2)))
826 : else
827 : val = s% rho_start(k-1)*&
828 0 : ((1/s% rho(k-1) - 1/s% rho_start(k-1)) - (1/s% rho(k) - 1/s% rho_start(k)))
829 : end if
830 :
831 : case (p_tau)
832 0 : val = s% tau(k)
833 : case (p_logtau)
834 0 : val = safe_log(s% tau(k))/ln10
835 : case (p_xtau)
836 0 : val = s% tau(nz) - s% tau(k)
837 : case (p_xlogtau)
838 0 : val = safe_log10(s% tau(nz) - s% tau(k))
839 : case (p_logtau_sub_xlogtau)
840 0 : val = safe_log10(s% tau(k)) - safe_log10(s% tau(nz) - s% tau(k))
841 :
842 : case (p_tau_eff)
843 0 : val = tau_eff(s,k)
844 : case (p_tau_eff_div_tau)
845 0 : val = tau_eff(s,k)/s% tau(k)
846 :
847 : case (p_kap_frac_lowT)
848 0 : val = s% kap_frac_lowT(k)
849 : case (p_kap_frac_highT)
850 0 : val = s% kap_frac_highT(k)
851 : case (p_kap_frac_Type2)
852 0 : val = s% kap_frac_Type2(k)
853 : case (p_kap_frac_Compton)
854 0 : val = s% kap_frac_Compton(k)
855 : case (p_kap_frac_op_mono)
856 0 : val = s% kap_frac_op_mono(k)
857 : case (p_log_kap)
858 0 : val = safe_log10(s% opacity(k))
859 : case (p_log_opacity)
860 0 : val = safe_log10(s% opacity(k))
861 : case (p_extra_opacity_factor)
862 0 : val = s% extra_opacity_factor(k)
863 : case (p_log_kap_times_factor)
864 0 : val = safe_log10(s% opacity(k)*s% extra_opacity_factor(k))
865 : case (p_energy)
866 0 : val = s% energy(k)
867 : case (p_logM)
868 0 : val = safe_log10(s% m(k)/Msun)
869 : case (p_temperature)
870 0 : val = s% T(k)
871 : case (p_logT)
872 3456 : val = s% lnT(k)/ln10
873 :
874 : case (p_logT_face)
875 0 : if (k == 1) then
876 0 : val = safe_log10(s% T_surf)
877 : else
878 : val = (s% dq(k-1)*s% lnT(k) + &
879 0 : s% dq(k)*s% lnT(k-1))/(s% dq(k-1) + s% dq(k))/ln10
880 : end if
881 : case (p_logT_bb)
882 : val = safe_log10( &
883 0 : pow(s% L(k)/(pi4*s% r(k)*s% r(k)*boltz_sigma), 0.25d0))
884 : case (p_logT_face_div_logT_bb)
885 0 : if (k == 1) then
886 0 : val = safe_log10(s% Teff)
887 : else
888 : val = (s% dq(k-1)*s% lnT(k) + &
889 0 : s% dq(k)*s% lnT(k-1))/(s% dq(k-1) + s% dq(k))/ln10
890 : end if
891 : val = val / safe_log10( &
892 0 : pow(s% L(k)/(pi4*s% r(k)*s% r(k)*boltz_sigma), 0.25d0))
893 :
894 : case (p_density)
895 0 : val = s% rho(k)
896 : case (p_rho)
897 0 : val = s% rho(k)
898 : case (p_logRho)
899 3456 : val = s% lnd(k)/ln10
900 : case (p_pgas)
901 0 : val = s% Pgas(k)
902 : case (p_logPgas)
903 0 : val = s% lnPgas(k)/ln10
904 : case (p_prad)
905 0 : val = s% Prad(k)
906 : case (p_pressure)
907 0 : val = s% Peos(k)
908 : case (p_logP)
909 3456 : val = s% lnPeos(k)/ln10
910 : case (p_logE)
911 0 : val = s% lnE(k)/ln10
912 : case (p_grada)
913 0 : val = s% grada(k)
914 : case (p_dE_dRho)
915 0 : val = s% dE_dRho(k)
916 : case (p_cv)
917 0 : val = s% Cv(k)
918 : case (p_cp)
919 0 : val = s% cp(k)
920 :
921 : case (p_thermal_time_to_surface)
922 0 : if (s% L(1) > 0) &
923 0 : val = sum(s% dm(1:k)*s% cp(1:k)*s% T(1:k))/s% L(1)
924 : case (p_log_thermal_time_to_surface)
925 0 : if (s% L(1) > 0) then
926 0 : val = sum(s% dm(1:k)*s% cp(1:k)*s% T(1:k))/s% L(1)
927 0 : val = safe_log10(val)
928 : end if
929 :
930 : case (p_log_CpT)
931 0 : val = safe_log10(s% cp(k)*s% T(k))
932 : case (p_log_CpT_absMdot_div_L)
933 0 : val = safe_log10(s% cp(k)*s% T(k)*abs(s% mstar_dot)/max(1d-99,s% L(k)))
934 : case (p_logS)
935 0 : val = s% lnS(k)/ln10
936 : case (p_logS_per_baryon)
937 0 : val = s% lnS(k)/ln10 + log10(amu)
938 : case (p_gamma1)
939 0 : val = s% gamma1(k)
940 : case (p_gamma3)
941 0 : val = s% gamma3(k)
942 : case (p_eta)
943 0 : val = s% eta(k)
944 : case (p_gam)
945 0 : val = s% gam(k)
946 : case (p_mu)
947 0 : val = s% mu(k)
948 :
949 : case (p_eos_frac_OPAL_SCVH)
950 0 : val = s% eos_frac_OPAL_SCVH(k)
951 : case (p_eos_frac_HELM)
952 0 : val = s% eos_frac_HELM(k)
953 : case (p_eos_frac_Skye)
954 0 : val = s% eos_frac_Skye(k)
955 : case (p_eos_frac_PC)
956 0 : val = s% eos_frac_PC(k)
957 : case (p_eos_frac_FreeEOS)
958 0 : val = s% eos_frac_FreeEOS(k)
959 : case (p_eos_frac_CMS)
960 0 : val = s% eos_frac_CMS(k)
961 : case (p_eos_frac_ideal)
962 0 : val = s% eos_frac_ideal(k)
963 :
964 : case (p_log_c_div_tau)
965 0 : val = safe_log10(clight/s% tau(k))
966 : case (p_log_v_escape)
967 0 : val = safe_log10(sqrt(2*s% cgrav(k)*s% m(k)/(s% r(k))))
968 : case (p_v_div_vesc)
969 0 : if (s% u_flag) then
970 0 : val = s% u(k)
971 0 : else if (s% v_flag) then
972 0 : val = s% v(k)
973 : end if
974 0 : val = val/sqrt(2*s% cgrav(k)*s% m(k)/(s% r(k)))
975 : case (p_v_div_v_escape)
976 0 : if (s% u_flag) then
977 0 : val = s% u_face_ad(k)%val
978 0 : else if (s% v_flag) then
979 0 : val = s% v(k)
980 : end if
981 0 : val = val/sqrt(2d0*s% cgrav(k)*s% m(k)/(s% r(k)))
982 : case (p_v_div_cs)
983 0 : val = s% v_div_csound(k)
984 : case (p_v_div_csound)
985 0 : val = s% v_div_csound(k)
986 : case (p_log_csound)
987 0 : val = safe_log10(s% csound(k))
988 : case (p_csound)
989 0 : val = s% csound(k)
990 : case (p_csound_face)
991 0 : val = s% csound_face(k)
992 : case (p_scale_height)
993 0 : val = s% scale_height(k)/Rsun
994 : case (p_entropy)
995 0 : val = s% entropy(k)
996 : case (p_free_e)
997 0 : val = exp(s% lnfree_e(k))
998 : case (p_logfree_e)
999 0 : val = s% lnfree_e(k)/ln10
1000 : case (p_chiRho)
1001 0 : val = s% chiRho(k)
1002 : case (p_chiT)
1003 0 : val = s% chiT(k)
1004 : case (p_QQ)
1005 0 : val = s% QQ(k)
1006 :
1007 : case (p_eos_phase)
1008 0 : val = s% phase(k)
1009 : case (p_latent_ddlnT)
1010 0 : val = s% latent_ddlnT(k)
1011 : case (p_latent_ddlnRho)
1012 0 : val = s% latent_ddlnRho(k)
1013 :
1014 : case (p_chiRho_for_partials)
1015 0 : val = s% chiRho_for_partials(k)
1016 : case (p_chiT_for_partials)
1017 0 : val = s% chiT_for_partials(k)
1018 : case (p_rel_diff_chiRho_for_partials)
1019 0 : val = (s% chiRho_for_partials(k) - s% chiRho(k))/s% chiRho(k)
1020 : case (p_rel_diff_chiT_for_partials)
1021 0 : val = (s% chiT_for_partials(k) - s% chiT(k))/s% chiT(k)
1022 :
1023 : case (p_x_mass_fraction_H)
1024 3456 : val = s% X(k)
1025 : case (p_y_mass_fraction_He)
1026 3456 : val = s% Y(k)
1027 : case (p_z_mass_fraction_metals)
1028 3456 : val = s% Z(k)
1029 :
1030 : case (p_abar)
1031 0 : val = s% abar(k)
1032 : case (p_zbar)
1033 0 : val = s% zbar(k)
1034 : case (p_z2bar)
1035 0 : val = s% z2bar(k)
1036 : case (p_ye)
1037 0 : val = s% ye(k)
1038 : case (p_opacity)
1039 0 : val = s% opacity(k)
1040 : case (p_dkap_dlnrho_face)
1041 0 : val = interp_val_to_pt(s% d_opacity_dlnd,k,nz,s% dq,'p_dkap_dlnrho_face')
1042 : case (p_dkap_dlnT_face)
1043 0 : val = interp_val_to_pt(s% d_opacity_dlnT,k,nz,s% dq,'p_dkap_dlnT_face')
1044 :
1045 : case (p_eps_nuc)
1046 0 : val = s% eps_nuc(k)
1047 : case (p_signed_log_eps_nuc)
1048 0 : val = s% eps_nuc(k)
1049 0 : val = sign(1d0,val)*log10(max(1d0,abs(val)))
1050 : case (p_log_abs_eps_nuc)
1051 0 : val = safe_log10(abs(s% eps_nuc(k)))
1052 : case (p_d_epsnuc_dlnd)
1053 0 : val = s% d_epsnuc_dlnd(k)
1054 : case (p_d_lnepsnuc_dlnd)
1055 0 : val = s% d_epsnuc_dlnd(k)/max(1d0,abs(s% eps_nuc(k)))
1056 : case (p_d_epsnuc_dlnT)
1057 0 : val = s% d_epsnuc_dlnT(k)
1058 : case (p_d_lnepsnuc_dlnT)
1059 0 : val = s% d_epsnuc_dlnT(k)/max(1d0,abs(s% eps_nuc(k)))
1060 :
1061 : case (p_deps_dlnd_face)
1062 0 : val = interp_val_to_pt(s% d_epsnuc_dlnd,k,nz,s% dq,'p_deps_dlnd_face')
1063 : case (p_deps_dlnT_face)
1064 0 : val = interp_val_to_pt(s% d_epsnuc_dlnT,k,nz,s% dq,'p_deps_dlnT_face')
1065 : case (p_eps_nuc_neu_total)
1066 0 : val = s% eps_nuc_neu_total(k)
1067 : case (p_non_nuc_neu)
1068 0 : val = s% non_nuc_neu(k)
1069 : case (p_nonnucneu_plas)
1070 0 : val = s% nonnucneu_plas(k)
1071 : case (p_nonnucneu_brem)
1072 0 : val = s% nonnucneu_brem(k)
1073 : case (p_nonnucneu_phot)
1074 0 : val = s% nonnucneu_phot(k)
1075 : case (p_nonnucneu_pair)
1076 0 : val = s% nonnucneu_pair(k)
1077 : case (p_nonnucneu_reco)
1078 0 : val = s% nonnucneu_reco(k)
1079 :
1080 : case (p_log_irradiation_heat)
1081 0 : val = safe_log10(s% irradiation_heat(k))
1082 : case (p_cgrav_factor)
1083 0 : val = s% cgrav(k)/standard_cgrav
1084 : case (p_alpha_mlt)
1085 0 : val = s% alpha_mlt(k)
1086 :
1087 : case (p_extra_jdot)
1088 0 : val = s% extra_jdot(k)
1089 : case (p_extra_omegadot)
1090 0 : val = s% extra_omegadot(k)
1091 : case (p_extra_heat)
1092 0 : val = s% extra_heat(k)%val
1093 : case (p_extra_grav)
1094 0 : val = s% extra_grav(k)%val
1095 : case (p_extra_L)
1096 0 : val = dot_product(s% dm(k:s% nz),s% extra_heat(k:s% nz)%val)/Lsun
1097 : case (p_log_extra_L)
1098 : val = safe_log10( &
1099 0 : dot_product(s% dm(k:s% nz),s% extra_heat(k:s% nz)%val)/Lsun)
1100 :
1101 : case (p_log_abs_eps_grav_dm_div_L)
1102 : val = safe_log10( &
1103 0 : abs(s% eps_grav_ad(k)% val)*s% dm(k)/max(1d0,abs(s% L(k))))
1104 :
1105 : case (p_eps_grav_composition_term)
1106 0 : if (s% include_composition_in_eps_grav) &
1107 0 : val = s% eps_grav_composition_term(k)
1108 :
1109 : case (p_eps_grav_plus_eps_mdot)
1110 0 : val = s% eps_grav_ad(k)% val + s% eps_mdot(k)
1111 : case (p_ergs_eps_grav_plus_eps_mdot)
1112 0 : val = (s% eps_grav_ad(k)% val + s% eps_mdot(k))*s% dm(k)*s% dt
1113 :
1114 : case (p_eps_mdot)
1115 0 : val = s% eps_mdot(k)
1116 : case (p_ergs_mdot)
1117 0 : val = s% eps_mdot(k)*s% dm(k)*s% dt
1118 :
1119 : case (p_div_v)
1120 0 : if (s% v_flag) then
1121 0 : if (k == s% nz) then
1122 0 : vp1 = s% V_center
1123 0 : Ap1 = pi4*s% R_center*s% R_center
1124 : else
1125 0 : vp1 = s% v(k+1)
1126 0 : Ap1 = pi4*s% r(k+1)*s% r(k+1)
1127 : end if
1128 0 : val = (pi4*s% r(k)*s% r(k)*s% v(k) - Ap1*vp1)*s% rho(k)/s% dm(k)
1129 : end if
1130 :
1131 : case (p_d_v_div_r_dm)
1132 0 : if (s% v_flag) then
1133 0 : if (k == s% nz) then
1134 0 : vp1 = s% V_center
1135 0 : rp1 = s% R_center
1136 : else
1137 0 : vp1 = s% v(k+1)
1138 0 : rp1 = s% r(k+1)
1139 : end if
1140 0 : v00 = s% v(k)
1141 0 : r00 = s% r(k)
1142 0 : if (rp1 > 0) then
1143 0 : val = (v00/r00 - vp1/rp1)/s% dm(k)
1144 : end if
1145 : end if
1146 :
1147 : case (p_d_v_div_r_dr)
1148 0 : if (s% v_flag) then
1149 0 : if (k == s% nz) then
1150 0 : vp1 = s% V_center
1151 0 : rp1 = s% R_center
1152 : else
1153 0 : vp1 = s% v(k+1)
1154 0 : rp1 = s% r(k+1)
1155 : end if
1156 0 : v00 = s% v(k)
1157 0 : r00 = s% r(k)
1158 0 : if (rp1 > 0) then
1159 : val = pi4*s% rmid(k)*s% rmid(k)*s% rho(k)* &
1160 0 : (v00/r00 - vp1/rp1)/s% dm(k)
1161 : end if
1162 : end if
1163 :
1164 : case (p_rho_times_r3)
1165 0 : val = s% rho_face(k)*s% r(k)*s% r(k)*s% r(k)
1166 : case (p_log_rho_times_r3)
1167 0 : val = safe_log10(s% rho_face(k)*s% r(k)*s% r(k)*s% r(k))
1168 :
1169 : case(p_du)
1170 0 : if (s% u_flag) then
1171 0 : if (k == s% nz) then
1172 0 : val = s% u(k)
1173 : else
1174 0 : val = s% u(k) - s% u(k+1)
1175 : end if
1176 : end if
1177 :
1178 : case(p_P_face)
1179 0 : if (s% u_flag) val = s% P_face_ad(k)%val
1180 : case(p_log_P_face)
1181 0 : if (s% u_flag) val = safe_log10(s% P_face_ad(k)%val)
1182 :
1183 : case (p_dPdr_div_grav)
1184 0 : if (k > 1 .and. k < nz .and. s% cgrav(k) > 0d0 .and. s% RTI_flag) then
1185 0 : val = s% dPdr_info(k)/s% rho_face(k)
1186 : end if
1187 :
1188 : case (p_gradP_div_rho)
1189 0 : if (k > 1) val = pi4*s% r(k)*s% r(k)*(s% Peos(k-1) - s% Peos(k))/s% dm_bar(k)
1190 : case (p_dlnP_dlnR)
1191 0 : if (k > 1) val = log(s% P_face_ad(k-1)%val/s% P_face_ad(k)%val) / (s% lnR(k-1) - s% lnR(k))
1192 : case (p_dlnRho_dlnR)
1193 0 : if (k > 1) val = log(s% rho_face(k-1)/s% rho_face(k)) / (s% lnR(k-1) - s% lnR(k))
1194 :
1195 : case (p_dvdt_grav)
1196 0 : val = -s% cgrav(k)*s% m(k)/(s% r(k)*s% r(k))
1197 : case (p_grav_eff)
1198 0 : int_val = if_rot_ad(s% fp_rot,k, alt=1.0d0)
1199 0 : val = s% dxh_v(k)/s%dt / (int_val * s% cgrav(k) * s% m(k) /(s% r(k)*s% r(k)))
1200 : case (p_dvdt_dPdm)
1201 0 : if (k > 1) val = -pi4*s% r(k)*s% r(k)*(s% Peos(k-1) - s% Peos(k))/s% dm_bar(k)
1202 :
1203 : case (p_dm_eps_grav)
1204 0 : val = s% eps_grav_ad(k)% val*s% dm(k)
1205 : case (p_eps_grav)
1206 0 : val = s% eps_grav_ad(k)% val
1207 :
1208 : case (p_log_xm_div_delta_m)
1209 0 : if(abs(s% dt*s% mstar_dot) > 0) val = safe_log10((s% m(1) - s% m(k))/abs(s% dt*s% mstar_dot))
1210 : case (p_xm_div_delta_m)
1211 0 : if(abs(s% dt*s% mstar_dot) > 0) val = (s% m(1) - s% m(k))/abs(s% dt*s% mstar_dot)
1212 :
1213 : case (p_env_eps_grav)
1214 : val = -s% gradT_sub_grada(k)*s% grav(k)*s% mstar_dot*s% Cp(k)*s% T(k) / &
1215 0 : (pi4*s% r(k)*s% r(k)*s% Peos(k))
1216 :
1217 : case (p_mlt_mixing_type)
1218 0 : int_val = s% mlt_mixing_type(k)
1219 0 : val = dble(int_val)
1220 0 : int_flag = .true.
1221 : case (p_mlt_mixing_length)
1222 0 : val = s% mlt_mixing_length(k)
1223 : case (p_mlt_Gamma)
1224 0 : val = s% mlt_Gamma(k)
1225 : case (p_mlt_Zeta)
1226 0 : if (abs(s% gradr(k) - s% grada_face(k)) > 1d-20) &
1227 0 : val = (s% gradr(k) - s% gradT(k))/(s% gradr(k) - s% grada_face(k))
1228 : case (p_mlt_Pturb)
1229 0 : if (s% mlt_Pturb_factor > 0d0 .and. s% okay_to_set_mlt_vc) then
1230 0 : if (s% mlt_vc_old(k) > 0d0) &
1231 0 : val = s% mlt_Pturb_factor*pow2(s% mlt_vc(k))*get_rho_face_val(s,k)/3d0
1232 : end if
1233 : case (p_grad_density)
1234 0 : val = s% grad_density(k)
1235 : case (p_grad_temperature)
1236 0 : val = s% grad_temperature(k)
1237 :
1238 : case (p_gradL_sub_gradr)
1239 0 : val = s% gradL(k) - s% gradr(k)
1240 : case (p_grada_sub_gradr)
1241 0 : val = s% grada_face(k) - s% gradr(k)
1242 :
1243 : case (p_gradL)
1244 0 : val = s% gradL(k)
1245 : case (p_sch_stable)
1246 0 : if (s% grada(k) > s% gradr(k)) val = 1
1247 : case (p_ledoux_stable)
1248 0 : if (s% gradL(k) > s% gradr(k)) val = 1
1249 :
1250 : case (p_eps_nuc_start)
1251 0 : val = s% eps_nuc_start(k)
1252 :
1253 : case (p_dominant_isoA_for_thermohaline)
1254 0 : int_val = chem_isos% Z_plus_N(s% dominant_iso_for_thermohaline(k))
1255 0 : int_flag = .true.
1256 : case (p_dominant_isoZ_for_thermohaline)
1257 0 : int_val = chem_isos% Z(s% dominant_iso_for_thermohaline(k))
1258 0 : int_flag = .true.
1259 : case (p_gradL_composition_term)
1260 0 : val = s% gradL_composition_term(k)
1261 :
1262 : case (p_log_D_conv)
1263 0 : if (s% mixing_type(k) == convective_mixing) then
1264 0 : val = safe_log10(s% D_mix_non_rotation(k))
1265 : else
1266 0 : val = -99
1267 : end if
1268 : case (p_log_D_leftover)
1269 0 : if (s% mixing_type(k) == leftover_convective_mixing) then
1270 0 : val = safe_log10(s% D_mix_non_rotation(k))
1271 : else
1272 0 : val = -99
1273 : end if
1274 : case (p_log_D_semi)
1275 0 : if (s% mixing_type(k) == semiconvective_mixing) then
1276 0 : val = safe_log10(s% D_mix_non_rotation(k))
1277 : else
1278 0 : val = -99
1279 : end if
1280 : case (p_log_D_ovr)
1281 0 : if (s% mixing_type(k) == overshoot_mixing) then
1282 0 : val = safe_log10(s% D_mix_non_rotation(k))
1283 : else
1284 0 : val = -99
1285 : end if
1286 : case (p_log_D_rayleigh_taylor)
1287 0 : if(s% RTI_flag) then
1288 0 : val = safe_log10(s% eta_RTI(k))
1289 : else
1290 0 : val =-99
1291 : end if
1292 : case (p_log_D_anon)
1293 0 : if (s% mixing_type(k) == anonymous_mixing) then
1294 0 : val = safe_log10(s% D_mix_non_rotation(k))
1295 : else
1296 0 : val = -99
1297 : end if
1298 : case (p_log_D_thrm)
1299 0 : if (s% mixing_type(k) == thermohaline_mixing) then
1300 0 : val = safe_log10(s% D_mix_non_rotation(k))
1301 : else
1302 0 : val = -99
1303 : end if
1304 :
1305 : case (p_log_D_minimum)
1306 0 : if (s% mixing_type(k) == minimum_mixing) then
1307 0 : val = safe_log10(s% D_mix(k))
1308 : else
1309 0 : val = -99
1310 : end if
1311 :
1312 : case (p_log_lambda_RTI_div_Hrho)
1313 0 : if (s% RTI_flag) val = safe_log10( &
1314 0 : sqrt(s% alpha_RTI(k))*s% r(k)/s% rho(k)*abs(s% dRhodr_info(k)))
1315 : case (p_lambda_RTI)
1316 0 : if (s% RTI_flag) val = sqrt(s% alpha_RTI(k))*s% r(k)
1317 : case (p_dPdr_info)
1318 0 : if (s% RTI_flag) val = s% dPdr_info(k)
1319 : case (p_dRhodr_info)
1320 0 : if (s% RTI_flag) val = s% dRhodr_info(k)
1321 :
1322 : case (p_source_plus_alpha_RTI)
1323 0 : if (s% RTI_flag) val = s% source_plus_alpha_RTI(k)
1324 : case (p_log_source_plus_alpha_RTI)
1325 0 : if (s% RTI_flag) val = safe_log10(s% source_plus_alpha_RTI(k))
1326 : case (p_log_source_RTI)
1327 0 : if (s% RTI_flag) val = safe_log10(s% source_plus_alpha_RTI(k))
1328 : case (p_source_minus_alpha_RTI)
1329 0 : if (s% RTI_flag) val = s% source_minus_alpha_RTI(k)
1330 : case (p_log_source_minus_alpha_RTI)
1331 0 : if (s% RTI_flag) val = safe_log10(abs(s% source_minus_alpha_RTI(k)))
1332 :
1333 : case (p_dudt_RTI)
1334 0 : if (s% RTI_flag) val = s% dudt_RTI(k)
1335 : case (p_dedt_RTI)
1336 0 : if (s% RTI_flag) val = s% dedt_RTI(k)
1337 :
1338 : case (p_eta_RTI)
1339 0 : if (s% RTI_flag) val = s% eta_RTI(k)
1340 : case (p_log_eta_RTI)
1341 0 : if (s% RTI_flag) val = safe_log10(abs(s% eta_RTI(k)))
1342 : case (p_boost_for_eta_RTI)
1343 0 : if (s% RTI_flag) val = s% boost_for_eta_RTI(k)
1344 : case (p_log_boost_for_eta_RTI)
1345 0 : if (s% RTI_flag) val = safe_log10(abs(s% boost_for_eta_RTI(k)))
1346 :
1347 : case (p_alpha_RTI)
1348 0 : if (s% RTI_flag) val = s% alpha_RTI(k)
1349 : case (p_log_alpha_RTI)
1350 0 : if (s% RTI_flag) val = safe_log10(s% alpha_RTI(k))
1351 : case (p_log_etamid_RTI)
1352 0 : if (s% RTI_flag) val = safe_log10(s% etamid_RTI(k))
1353 :
1354 : case (p_log_sig_RTI)
1355 0 : if (s% RTI_flag) val = safe_log10(s% sig_RTI(k))
1356 : case (p_log_sigmid_RTI)
1357 0 : if (s% RTI_flag) val = safe_log10(s% sigmid_RTI(k))
1358 :
1359 : case (p_log_D_omega)
1360 0 : if (s% rotation_flag) val = safe_log10(s% D_omega(k))
1361 :
1362 : case (p_log_D_mix_non_rotation)
1363 0 : val = safe_log10(s% D_mix_non_rotation(k))
1364 : case (p_log_D_mix_rotation)
1365 0 : val = safe_log10(s% D_mix(k) - s% D_mix_non_rotation(k))
1366 : case (p_log_D_mix)
1367 0 : val = safe_log10(s% D_mix(k))
1368 : case (p_log_sig_mix)
1369 0 : val = safe_log10(s% sig(k))
1370 : case (p_log_sig_raw_mix)
1371 0 : val = safe_log10(s% sig_raw(k))
1372 :
1373 : case (p_burn_avg_epsnuc)
1374 0 : if (s% op_split_burn) val = s% burn_avg_epsnuc(k)
1375 : case (p_log_burn_avg_epsnuc)
1376 0 : if (s% op_split_burn) &
1377 0 : val = safe_log10(abs(s% burn_avg_epsnuc(k)))
1378 : case (p_burn_num_iters)
1379 0 : if (s% op_split_burn) then
1380 0 : int_val = s% burn_num_iters(k); val = dble(int_val)
1381 : else
1382 : int_val = 0; val = 0
1383 : end if
1384 0 : int_flag = .true.
1385 :
1386 : case (p_conv_vel_div_mlt_vc)
1387 0 : if (s% mlt_vc(k) > 0d0) val = s% conv_vel(k)/s% mlt_vc(k)
1388 :
1389 : case (p_conv_vel)
1390 0 : val = s% conv_vel(k)
1391 : case (p_dt_times_conv_vel_div_mixing_length)
1392 0 : val = s% dt*s% conv_vel(k)/s% mlt_mixing_length(k)
1393 : case (p_log_dt_times_conv_vel_div_mixing_length)
1394 0 : val = safe_log10(s% dt*s% conv_vel(k)/s% mlt_mixing_length(k))
1395 : case (p_log_conv_vel)
1396 0 : val = safe_log10(s% conv_vel(k))
1397 : case (p_conv_vel_div_L_vel)
1398 0 : val = s% conv_vel(k)/max(1d0,get_L_vel(k))
1399 : case (p_conv_vel_div_csound)
1400 0 : val = s% conv_vel(k)/s% csound_face(k)
1401 : case (p_dvc_dt_TDC_div_g)
1402 0 : val = s%dvc_dt_TDC(k) / s%grav(k)
1403 : case (p_mix_type)
1404 0 : val = dble(s% mixing_type(k))
1405 0 : int_val = s% mixing_type(k)
1406 0 : int_flag = .true.
1407 : case (p_mixing_type)
1408 0 : val = dble(s% mixing_type(k))
1409 0 : int_val = s% mixing_type(k)
1410 0 : int_flag = .true.
1411 : case (p_log_mlt_D_mix)
1412 0 : val = safe_log10(s% mlt_D(k))
1413 : case (p_log_t_thermal)
1414 0 : val = safe_log10(s% Cp(k)*s% T(k)*(s% m(1) - s% m(k))/s% L(k))
1415 : case (p_log_cp_T_div_t_sound)
1416 : val = safe_log10( &
1417 0 : s% Cp(k)*s% T(k)/(s% Peos(k)/(s% rho(k)*s% grav(k))/s% csound(k)))
1418 : case (p_log_t_sound)
1419 0 : val = safe_log10(s% Peos(k)/(s% rho(k)*s% grav(k))/s% csound(k))
1420 : case (p_pressure_scale_height)
1421 0 : val = s% Peos(k)/(s% rho(k)*s% grav(k))/Rsun
1422 : case (p_pressure_scale_height_cm)
1423 0 : val = s% Peos(k)/(s% rho(k)*s% grav(k))
1424 : case (p_gradT)
1425 0 : val = s% gradT(k)
1426 : case (p_gradr)
1427 0 : val = s% gradr(k)
1428 : case (p_grada_sub_gradT)
1429 0 : val = s% grada_face(k) - s% gradT(k)
1430 :
1431 : case (p_omega)
1432 0 : val = if_rot(s% omega,k)
1433 :
1434 : case (p_log_omega)
1435 0 : val = safe_log10(if_rot(s% omega,k))
1436 : case (p_log_j_rot)
1437 0 : val = safe_log10(if_rot(s% j_rot,k))
1438 : case (p_log_J_inside)
1439 0 : if (s% rotation_flag) then
1440 0 : val = safe_log10(dot_product(s% j_rot(k:s% nz), s% dm(k:s% nz)))
1441 : else
1442 0 : val = -99.0d0
1443 : end if
1444 : case (p_log_J_div_M53)
1445 0 : if (s% rotation_flag) then
1446 : val = safe_log10(&
1447 : dot_product(s% j_rot(k:s% nz), s% dm(k:s% nz)) * &
1448 0 : 1d-50/pow(s% m(k)/Msun,5d0/3d0))
1449 : else
1450 0 : val = -99.0d0
1451 : end if
1452 :
1453 : case (p_shear)
1454 0 : val = if_rot(s% omega_shear,k)
1455 : case (p_log_abs_shear)
1456 0 : if (s% rotation_flag) then
1457 0 : val = safe_log10(s% omega_shear(k))
1458 0 : if (is_bad(val)) then
1459 0 : write(*,2) 'val', k, val
1460 0 : write(*,2) 's% omega_shear(k)', k, s% omega_shear(k)
1461 0 : call mesa_error(__FILE__,__LINE__,'profile')
1462 : end if
1463 : else
1464 0 : val = -99
1465 : end if
1466 : case (p_log_abs_dlnR_domega)
1467 0 : if (s% rotation_flag) then
1468 0 : val = -safe_log10(s% omega_shear(k))
1469 0 : if (is_bad(val)) then
1470 0 : write(*,2) 'val', k, val
1471 0 : write(*,2) 's% omega_shear(k)', k, s% omega_shear(k)
1472 0 : call mesa_error(__FILE__,__LINE__,'profile')
1473 : end if
1474 : else
1475 0 : val = -99
1476 : end if
1477 : case (p_i_rot)
1478 0 : val = if_rot_ad(s% i_rot,k)
1479 : case (p_j_rot)
1480 0 : val = if_rot(s% j_rot,k)
1481 : case (p_v_rot)
1482 0 : val = if_rot(s% omega,k)*if_rot(s% r_equatorial,k)*1d-5 ! km/sec
1483 : case (p_fp_rot)
1484 0 : val = if_rot_ad(s% fp_rot,k, alt=1.0d0)
1485 : case (p_ft_rot)
1486 0 : val = if_rot_ad(s% ft_rot,k, alt=1.0d0)
1487 : case (p_ft_rot_div_fp_rot)
1488 0 : if(s% rotation_flag) then
1489 0 : val = s% ft_rot(k)% val/s% fp_rot(k)% val
1490 : else
1491 0 : val = 1.0d0
1492 : end if
1493 : case (p_w_div_w_crit_roche)
1494 0 : val = if_rot(s% w_div_w_crit_roche,k)
1495 : case (p_w_div_w_crit_roche2)
1496 0 : val = if_rot(s% xh(s% i_w_div_wc,:),k)
1497 : case (p_log_am_nu_non_rot)
1498 0 : val = safe_log10(if_rot(s% am_nu_non_rot,k))
1499 : case (p_log_am_nu_rot)
1500 0 : val = safe_log10(if_rot(s% am_nu_rot,k))
1501 : case (p_log_am_nu)
1502 0 : val = safe_log10(if_rot(s% am_nu_rot,k) + if_rot(s% am_nu_non_rot,k))
1503 :
1504 : case (p_r_polar)
1505 0 : val = if_rot(s% r_polar,k, alt=s% r(k))/Rsun
1506 : case (p_log_r_polar)
1507 0 : val = safe_log10(if_rot(s% r_polar,k, alt=s% r(k))/Rsun)
1508 : case (p_r_equatorial)
1509 0 : val = if_rot(s% r_equatorial,k, alt=s% r(k))/Rsun
1510 : case (p_log_r_equatorial)
1511 0 : val = safe_log10(if_rot(s% r_equatorial,k, alt=s% r(k))/Rsun)
1512 : case (p_r_e_div_r_p)
1513 0 : if (s% rotation_flag) then
1514 0 : if(s% r_polar(k) > 1) val = s% r_equatorial(k)/s% r_polar(k)
1515 : end if
1516 : case (p_omega_crit)
1517 0 : val = omega_crit(s,k)
1518 : case (p_omega_div_omega_crit)
1519 0 : if (s% rotation_flag) then
1520 0 : val = omega_crit(s,k)
1521 0 : if (val < 1d-50) then
1522 0 : val = 0
1523 : else
1524 0 : val = s% omega(k)/val
1525 : end if
1526 : end if
1527 :
1528 : case (p_eps_phase_separation)
1529 0 : if (s% do_phase_separation .and. s% do_phase_separation_heating) val = s% eps_phase_separation(k)
1530 :
1531 : case (p_eps_WD_sedimentation)
1532 0 : if (s% do_element_diffusion) val = s% eps_WD_sedimentation(k)
1533 : case (p_log_eps_WD_sedimentation)
1534 0 : if (s% do_element_diffusion) val = safe_log10(s% eps_WD_sedimentation(k))
1535 :
1536 : case (p_eps_diffusion)
1537 0 : if (s% do_element_diffusion) val = s% eps_diffusion(k)
1538 : case (p_log_eps_diffusion)
1539 0 : if (s% do_element_diffusion) val = safe_log10(s% eps_diffusion(k))
1540 :
1541 : case (p_e_field)
1542 0 : if (s% do_element_diffusion) val = s% E_field(k)
1543 : case (p_log_e_field)
1544 0 : if (s% do_element_diffusion) val = safe_log10(s% E_field(k))
1545 :
1546 : case (p_g_field_element_diffusion)
1547 0 : if (s% do_element_diffusion) val = s% g_field_element_diffusion(k)
1548 : case (p_log_g_field_element_diffusion)
1549 0 : if (s% do_element_diffusion) &
1550 0 : val = safe_log10(s% g_field_element_diffusion(k))
1551 :
1552 : case (p_eE_div_mg_element_diffusion)
1553 0 : if (s% do_element_diffusion) then
1554 0 : if ( s% g_field_element_diffusion(k) /= 0d0) then
1555 0 : val = qe * s% E_field(k)/(amu * s% g_field_element_diffusion(k))
1556 : else
1557 : val = 0d0
1558 : end if
1559 : end if
1560 : case (p_log_eE_div_mg_element_diffusion)
1561 0 : if (s% do_element_diffusion) &
1562 0 : val = safe_log10(qe * s% E_field(k)/(amu * s% g_field_element_diffusion(k)))
1563 :
1564 : case (p_richardson_number)
1565 0 : val = if_rot(s% richardson_number,k)
1566 : case (p_am_domega_dlnR)
1567 0 : val = if_rot(s% domega_dlnR,k)
1568 :
1569 : case (p_am_log_sig) ! == am_log_sig_omega
1570 0 : val = safe_log10(if_rot(s% am_sig_omega,k))
1571 : case (p_am_log_sig_omega)
1572 0 : val = safe_log10(if_rot(s% am_sig_omega,k))
1573 : case (p_am_log_sig_j)
1574 0 : val = safe_log10(if_rot(s% am_sig_j,k))
1575 :
1576 : case (p_am_log_nu_omega)
1577 0 : val = safe_log10(if_rot(s% am_nu_omega,k))
1578 : case (p_am_log_nu_j)
1579 0 : val = safe_log10(if_rot(s% am_nu_j,k))
1580 :
1581 : case (p_am_log_nu_rot)
1582 0 : val = safe_log10(if_rot(s% am_nu_rot,k))
1583 : case (p_am_log_nu_non_rot)
1584 0 : val = safe_log10(if_rot(s% am_nu_non_rot,k))
1585 :
1586 : case (p_am_log_D_visc)
1587 0 : if (s% am_nu_visc_factor >= 0) then
1588 : f = s% am_nu_visc_factor
1589 : else
1590 0 : f = s% D_visc_factor
1591 : end if
1592 0 : val = safe_log10(am_nu_factor*f*if_rot(s% D_visc,k))
1593 : case (p_am_log_D_DSI)
1594 0 : if (s% am_nu_DSI_factor >= 0) then
1595 : f = s% am_nu_DSI_factor
1596 : else
1597 0 : f = s% D_DSI_factor
1598 : end if
1599 0 : val = safe_log10(am_nu_factor*f*if_rot(s% D_DSI,k))
1600 : case (p_am_log_D_SH)
1601 0 : if (s% am_nu_SH_factor >= 0) then
1602 : f = s% am_nu_SH_factor
1603 : else
1604 0 : f = s% D_SH_factor
1605 : end if
1606 0 : val = safe_log10(am_nu_factor*f*if_rot(s% D_SH,k))
1607 : case (p_am_log_D_SSI)
1608 0 : if (s% am_nu_SSI_factor >= 0) then
1609 : f = s% am_nu_SSI_factor
1610 : else
1611 0 : f = s% D_SSI_factor
1612 : end if
1613 0 : val = safe_log10(am_nu_factor*f*if_rot(s% D_SSI,k))
1614 :
1615 : case (p_am_log_D_ES)
1616 0 : if (s% am_nu_ES_factor >= 0) then
1617 : f = s% am_nu_ES_factor
1618 : else
1619 0 : f = s% D_ES_factor
1620 : end if
1621 0 : val = safe_log10(am_nu_factor*f*if_rot(s% D_ES,k))
1622 : case (p_am_log_D_GSF)
1623 0 : if (s% am_nu_GSF_factor >= 0) then
1624 : f = s% am_nu_GSF_factor
1625 : else
1626 0 : f = s% D_GSF_factor
1627 : end if
1628 0 : val = safe_log10(am_nu_factor*f*if_rot(s% D_GSF,k))
1629 : case (p_am_log_D_ST)
1630 0 : if (s% am_nu_ST_factor >= 0) then
1631 : f = s% am_nu_ST_factor
1632 : else
1633 0 : f = s% D_ST_factor
1634 : end if
1635 0 : val = safe_log10(am_nu_factor*f*if_rot(s% D_ST,k))
1636 : case (p_am_log_nu_ST)
1637 0 : if (s% am_nu_ST_factor >= 0) then
1638 : f = s% am_nu_ST_factor
1639 : else
1640 0 : f = s% D_ST_factor
1641 : end if
1642 0 : val = safe_log10(am_nu_factor*f*if_rot(s% nu_ST,k))
1643 :
1644 : case (p_dynamo_log_B_r)
1645 0 : val = safe_log10(if_rot(s% dynamo_B_r,k))
1646 : case (p_dynamo_log_B_phi)
1647 0 : val = safe_log10(if_rot(s% dynamo_B_phi,k))
1648 :
1649 : case (p_grada_face)
1650 0 : val = s% grada_face(k)
1651 : case (p_gradr_div_grada)
1652 0 : val = s% gradr(k)/s% grada_face(k)
1653 : case (p_gradr_sub_grada)
1654 0 : val = s% gradr(k) - s% grada_face(k)
1655 : case (p_gradT_sub_a)
1656 0 : val = s% gradT(k) - s% grada_face(k)
1657 : case (p_gradT_sub_grada)
1658 0 : val = s% gradT(k) - s% grada_face(k)
1659 : case (p_gradT_div_grada)
1660 0 : val = s% gradT(k) / s% grada_face(k)
1661 : case (p_gradr_sub_gradT)
1662 0 : val = s% gradr(k) - s% gradT(k)
1663 : case (p_gradT_sub_gradr)
1664 0 : val = s% gradT(k) - s% gradr(k)
1665 :
1666 : case (p_gradT_rel_err)
1667 0 : if (k > 1) then
1668 0 : val = (s% lnT(k-1) - s% lnT(k))/(s% lnPeos(k-1) - s% lnPeos(k))
1669 0 : val = (s% gradT(k) - val)/s% gradT(k)
1670 : end if
1671 :
1672 : case (p_gradT_div_gradr)
1673 0 : if (abs(s% gradr(k)) < 1d-99) then
1674 0 : val = 1d0
1675 : else
1676 0 : val = s% gradT(k) / s% gradr(k)
1677 : end if
1678 : case (p_log_gradT_div_gradr)
1679 0 : if (abs(s% gradr(k)) < 1d-99) then
1680 : val = 0d0
1681 : else
1682 0 : val = safe_log10(s% gradT(k) / s% gradr(k))
1683 : end if
1684 :
1685 : case (p_log_mlt_Gamma)
1686 0 : val = safe_log10(s% mlt_Gamma(k))
1687 : case (p_log_mlt_vc)
1688 0 : val = safe_log10(s% mlt_vc(k))
1689 : case (p_mlt_vc)
1690 0 : val = s% mlt_vc(k)
1691 : case (p_mlt_D)
1692 0 : val = s% mlt_D(k)
1693 : case (p_mlt_gradT)
1694 0 : val = s% mlt_gradT(k)
1695 : case (p_mlt_Y_face)
1696 0 : val = s% Y_face(k)
1697 : case (p_mlt_log_abs_Y)
1698 0 : val = safe_log10(abs(s% Y_face(k)))
1699 : case (p_tdc_num_iters)
1700 0 : int_val = s% tdc_num_iters(k); val = dble(int_val)
1701 0 : int_flag = .true.
1702 : case(p_COUPL)
1703 0 : val = s% COUPL(k)
1704 : case(p_SOURCE)
1705 0 : val = s% SOURCE(k)
1706 : case(p_DAMP)
1707 0 : val = s% DAMP(k)
1708 : case(p_DAMPR)
1709 0 : val = s% DAMPR(k)
1710 :
1711 : case (p_delta_r)
1712 0 : val = s% r(k) - s% r_start(k)
1713 : case (p_delta_L)
1714 0 : val = s% L(k) - s% L_start(k)
1715 : case (p_delta_cell_vol)
1716 0 : if (k == s% nz) then
1717 0 : rp1 = s% R_center
1718 0 : rp1_start = s% R_center_old
1719 : else
1720 0 : rp1 = s% r(k+1)
1721 0 : rp1_start = s% r_start(k+1)
1722 : end if
1723 0 : r00 = s% r(k)
1724 0 : r00_start = s% r_start(k)
1725 0 : dr3 = r00*r00*r00 - rp1*rp1*rp1
1726 0 : dr3_start = r00_start*r00_start*r00_start - rp1_start*rp1_start*rp1_start
1727 0 : val = four_thirds_pi*(dr3 - dr3_start)
1728 : case (p_delta_entropy)
1729 0 : val = s% entropy(k) - exp(s% lnS_start(k))/(avo*kerg)
1730 : case (p_delta_T)
1731 0 : val = s% T(k) - s% T_start(k)
1732 : case (p_delta_rho)
1733 0 : val = s% rho(k) - exp(s% lnd_start(k))
1734 : case (p_delta_eps_nuc)
1735 0 : val = s% eps_nuc(k) - s% eps_nuc_start(k)
1736 : case (p_delta_mu)
1737 0 : val = s% mu(k) - s% mu_start(k)
1738 :
1739 : case (p_cno_div_z)
1740 : cno = s% xa(s% net_iso(ic12),k) + &
1741 0 : s% xa(s% net_iso(in14),k) + s% xa(s% net_iso(io16),k)
1742 0 : z = 1 - (s% xa(s% net_iso(ih1),k) + s% xa(s% net_iso(ihe4),k))
1743 0 : if (z > 1d-50) then
1744 0 : val = cno/z
1745 : else
1746 : val = 0
1747 : end if
1748 : case (p_dE)
1749 0 : val = s% energy(k) - s% energy_start(k)
1750 : case (p_dr)
1751 0 : if (k < s% nz) then
1752 0 : val = s% r(k) - s% r(k+1)
1753 : else
1754 0 : val = s% r(k) - s% R_center
1755 : end if
1756 : case (p_dr_ratio)
1757 0 : if (k == 1 .or. k == s% nz) then
1758 0 : val = 1
1759 : else
1760 0 : val = (s% r(k-1) - s% r(k))/(s% r(k) - s% r(k+1))
1761 : end if
1762 : case (p_dv)
1763 0 : if (.not. s% v_flag) then
1764 : val = 0
1765 0 : else if (k < s% nz) then
1766 0 : val = s% v(k+1) - s% v(k)
1767 : else
1768 0 : val = -s% v(k)
1769 : end if
1770 : case (p_dt_dv_div_dr)
1771 0 : if (.not. s% v_flag) then
1772 : val = 0
1773 0 : else if (k < s% nz) then
1774 0 : val = s% dt*(s% v(k+1) - s% v(k))/(s% r(k) - s% r(k+1))
1775 : else
1776 0 : val = -s% dt*s% v(k)/s% r(k)
1777 : end if
1778 :
1779 : case (p_dlog_h1_dlogP)
1780 0 : val = get_dlogX_dlogP(ih1, k)
1781 : case (p_dlog_he3_dlogP)
1782 0 : val = get_dlogX_dlogP(ihe3, k)
1783 : case (p_dlog_he4_dlogP)
1784 0 : val = get_dlogX_dlogP(ihe4, k)
1785 : case (p_dlog_c12_dlogP)
1786 0 : val = get_dlogX_dlogP(ic12, k)
1787 : case (p_dlog_c13_dlogP)
1788 0 : val = get_dlogX_dlogP(ic13, k)
1789 : case (p_dlog_n14_dlogP)
1790 0 : val = get_dlogX_dlogP(in14, k)
1791 : case (p_dlog_o16_dlogP)
1792 0 : val = get_dlogX_dlogP(io16, k)
1793 : case (p_dlog_ne20_dlogP)
1794 0 : val = get_dlogX_dlogP(ine20, k)
1795 : case (p_dlog_mg24_dlogP)
1796 0 : val = get_dlogX_dlogP(img24, k)
1797 : case (p_dlog_si28_dlogP)
1798 0 : val = get_dlogX_dlogP(isi28, k)
1799 :
1800 : case (p_dlog_pp_dlogP)
1801 0 : val = get_dlog_eps_dlogP(ipp, k)
1802 : case (p_dlog_cno_dlogP)
1803 0 : val = get_dlog_eps_dlogP(icno, k)
1804 : case (p_dlog_3alf_dlogP)
1805 0 : val = get_dlog_eps_dlogP(i3alf, k)
1806 :
1807 : case (p_dlog_burn_c_dlogP)
1808 0 : val = get_dlog_eps_dlogP(i_burn_c, k)
1809 : case (p_dlog_burn_n_dlogP)
1810 0 : val = get_dlog_eps_dlogP(i_burn_n, k)
1811 : case (p_dlog_burn_o_dlogP)
1812 0 : val = get_dlog_eps_dlogP(i_burn_o, k)
1813 :
1814 : case (p_dlog_burn_ne_dlogP)
1815 0 : val = get_dlog_eps_dlogP(i_burn_ne, k)
1816 : case (p_dlog_burn_na_dlogP)
1817 0 : val = get_dlog_eps_dlogP(i_burn_na, k)
1818 : case (p_dlog_burn_mg_dlogP)
1819 0 : val = get_dlog_eps_dlogP(i_burn_mg, k)
1820 :
1821 : case (p_dlog_cc_dlogP)
1822 0 : val = get_dlog_eps_dlogP(icc, k)
1823 : case (p_dlog_co_dlogP)
1824 0 : val = get_dlog_eps_dlogP(ico, k)
1825 : case (p_dlog_oo_dlogP)
1826 0 : val = get_dlog_eps_dlogP(ioo, k)
1827 :
1828 : case (p_dlog_burn_si_dlogP)
1829 0 : val = get_dlog_eps_dlogP(i_burn_si, k)
1830 : case (p_dlog_burn_s_dlogP)
1831 0 : val = get_dlog_eps_dlogP(i_burn_s, k)
1832 : case (p_dlog_burn_ar_dlogP)
1833 0 : val = get_dlog_eps_dlogP(i_burn_ar, k)
1834 : case (p_dlog_burn_ca_dlogP)
1835 0 : val = get_dlog_eps_dlogP(i_burn_ca, k)
1836 : case (p_dlog_burn_ti_dlogP)
1837 0 : val = get_dlog_eps_dlogP(i_burn_ti, k)
1838 : case (p_dlog_burn_cr_dlogP)
1839 0 : val = get_dlog_eps_dlogP(i_burn_cr, k)
1840 : case (p_dlog_burn_fe_dlogP)
1841 0 : val = get_dlog_eps_dlogP(i_burn_fe, k)
1842 : case (p_dlog_pnhe4_dlogP)
1843 0 : val = get_dlog_eps_dlogP(ipnhe4, k)
1844 : case (p_dlog_photo_dlogP)
1845 0 : val = get_dlog_eps_dlogP(iphoto, k)
1846 : case (p_dlog_other_dlogP)
1847 0 : val = get_dlog_eps_dlogP(iother, k)
1848 :
1849 : case(p_d_u_div_rmid)
1850 0 : if (s% u_flag .and. k > 1) &
1851 0 : val = s% u(k-1)/s% rmid(k-1) - s% u(k)/s% rmid(k)
1852 : case(p_d_u_div_rmid_start)
1853 0 : if (s% u_flag .and. k > 1) &
1854 0 : val = s% u(k-1)/s% rmid_start(k-1) - s% u(k)/s% rmid_start(k)
1855 :
1856 : case(p_Ptrb)
1857 0 : if (s% RSP2_flag) then
1858 0 : val = get_etrb(s,k)*s% rho(k)
1859 0 : else if (s% RSP_flag) then
1860 0 : val = s% RSP_Et(k)*s% rho(k)
1861 : end if
1862 : case(p_log_Ptrb)
1863 0 : if (s% RSP2_flag) then
1864 0 : val = safe_log10(get_etrb(s,k)*s% rho(k))
1865 0 : else if (s% RSP_flag) then
1866 0 : val = safe_log10(s% RSP_Et(k)*s% rho(k))
1867 : end if
1868 : case(p_w)
1869 0 : if (s% RSP2_flag) then
1870 0 : val = get_w(s,k)
1871 0 : else if (s% RSP_flag) then
1872 0 : val = s% RSP_w(k)
1873 : else
1874 0 : val = s% mlt_vc(k)/sqrt_2_div_3
1875 : end if
1876 : case(p_log_w)
1877 0 : if (s% RSP2_flag) then
1878 0 : val = get_w(s,k)
1879 0 : else if (s% RSP_flag) then
1880 0 : val = s% RSP_w(k)
1881 : else
1882 0 : val = s% mlt_vc(k)/sqrt_2_div_3
1883 : end if
1884 0 : val = safe_log10(val)
1885 : case(p_etrb)
1886 0 : if (s% RSP2_flag) then
1887 0 : val = get_etrb(s,k)
1888 0 : else if (s% RSP_flag) then
1889 0 : val = s% RSP_Et(k)
1890 : end if
1891 : case(p_log_etrb)
1892 0 : if (s% RSP2_flag) then
1893 0 : val = safe_log10(get_etrb(s,k))
1894 0 : else if (s% RSP_flag) then
1895 0 : val = safe_log10(s% RSP_Et(k))
1896 : end if
1897 : case(p_Pvsc)
1898 0 : if (s% use_Pvsc_art_visc .or. s% RSP_flag) val = s% Pvsc(k)
1899 : case(p_Hp_face)
1900 0 : if (rsp_or_w) val = s% Hp_face(k)
1901 : case(p_Y_face)
1902 0 : if (rsp_or_w) val = s% Y_face(k)
1903 : case(p_PII_face)
1904 0 : if (rsp_or_w) val = s% PII(k)
1905 : case(p_Chi)
1906 0 : val = s% Chi(k)
1907 : case(p_Eq)
1908 0 : val = s% Eq(k)
1909 : case(p_Uq)
1910 0 : val = s% Uq(k)
1911 : case(p_Lr)
1912 0 : val = get_Lrad(s,k)
1913 : case(p_Lr_div_L)
1914 0 : val = get_Lrad(s,k)/s% L(k)
1915 : case(p_Lc)
1916 0 : val = get_Lconv(s,k)
1917 : case(p_Lc_div_L)
1918 0 : val = get_Lconv(s,k)/s% L(k)
1919 : case(p_Lt)
1920 0 : if (rsp_or_w) val = s% Lt(k)
1921 : case(p_Lt_div_L)
1922 0 : if (rsp_or_w) val = s% Lt(k)/s% L(k)
1923 : case(p_reconstructed_T_face)
1924 0 : if (reconstructed_face_state_active .and. s% reconstructed_face_state_valid(k)) val = s% reconstructed_T_face_ad(k)% val
1925 : case(p_reconstructed_rho_face)
1926 0 : if (reconstructed_face_state_active .and. s% reconstructed_face_state_valid(k)) val = s% reconstructed_rho_face_ad(k)% val
1927 : case(p_reconstructed_P_face)
1928 0 : if (reconstructed_face_state_active .and. s% reconstructed_face_state_valid(k)) val = s% reconstructed_P_face_ad(k)% val
1929 : case(p_reconstructed_Cp_face)
1930 0 : if (reconstructed_face_state_active .and. s% reconstructed_face_state_valid(k)) val = s% reconstructed_Cp_face_ad(k)% val
1931 : case(p_reconstructed_ChiRho_face)
1932 0 : if (reconstructed_face_state_active .and. s% reconstructed_face_state_valid(k)) val = s% reconstructed_ChiRho_face_ad(k)% val
1933 : case(p_reconstructed_ChiT_face)
1934 0 : if (reconstructed_face_state_active .and. s% reconstructed_face_state_valid(k)) val = s% reconstructed_ChiT_face_ad(k)% val
1935 : case(p_reconstructed_grada_face)
1936 0 : if (reconstructed_face_state_active .and. s% reconstructed_face_state_valid(k)) val = s% reconstructed_grada_face_ad(k)% val
1937 : case(p_reconstructed_opacity_face)
1938 0 : if (reconstructed_face_state_active .and. s% reconstructed_face_state_valid(k)) val = s% reconstructed_opacity_face_ad(k)% val
1939 : case(p_reconstructed_scale_height_face)
1940 0 : if (reconstructed_face_state_active .and. s% reconstructed_face_state_valid(k)) val = s% reconstructed_scale_height_face_ad(k)% val
1941 : case(p_reconstructed_gradr_face)
1942 0 : if (reconstructed_face_state_active .and. s% reconstructed_face_state_valid(k)) val = s% reconstructed_gradr_face_ad(k)% val
1943 :
1944 :
1945 : case(p_rsp_Et)
1946 0 : if (s% rsp_flag) val = s% RSP_Et(k)
1947 : case(p_rsp_logEt)
1948 0 : if (s% rsp_flag) &
1949 0 : val = safe_log10(s% RSP_Et(k))
1950 : case(p_rsp_Pt)
1951 0 : if (s% rsp_flag) val = s% Ptrb(k)
1952 : case(p_rsp_Eq)
1953 0 : if (s% rsp_flag) val = s% Eq(k)
1954 : case(p_rsp_src_snk)
1955 0 : if (s% rsp_flag) val = s% COUPL(k)
1956 : case(p_rsp_src)
1957 0 : if (s% rsp_flag) val = s% SOURCE(k)
1958 : case(p_rsp_sink)
1959 0 : if (s% rsp_flag) val = s% DAMP(k) + s% DAMPR(k)
1960 : case(p_rsp_damp)
1961 0 : if (s% rsp_flag) val = s% DAMP(k)
1962 : case(p_rsp_dampR)
1963 0 : if (s% rsp_flag) val = s% DAMPR(k)
1964 : case(p_rsp_Hp_face)
1965 0 : if (s% rsp_flag) val = s% Hp_face(k)
1966 : case(p_rsp_Chi)
1967 0 : if (s% rsp_flag) val = s% Chi(k)
1968 : case(p_rsp_Pvsc)
1969 0 : if (s% rsp_flag) val = s% Pvsc(k)
1970 : case(p_rsp_erad)
1971 0 : if (s% rsp_flag) val = s% erad(k)
1972 : case(p_rsp_log_erad)
1973 0 : if (s% rsp_flag) val = safe_log10(s% erad(k))
1974 : case(p_rsp_log_dt_div_heat_exchange_timescale)
1975 0 : if (s% rsp_flag) val = safe_log10(s% dt*clight*s% opacity(k)*s% rho(k))
1976 : case(p_rsp_heat_exchange_timescale)
1977 0 : if (s% rsp_flag) val = 1d0/(clight*s% opacity(k)*s% rho(k))
1978 : case(p_rsp_log_heat_exchange_timescale)
1979 0 : if (s% rsp_flag) &
1980 0 : val = safe_log10(1d0/(clight*s% opacity(k)*s% rho(k)))
1981 : case(p_rsp_Y_face)
1982 0 : if (s% rsp_flag) then
1983 0 : if (k > 1) then
1984 0 : val = s% Y_face(k)
1985 : else ! for plotting, use value at k=2
1986 0 : val = s% Y_face(2)
1987 : end if
1988 : end if
1989 : case(p_rsp_gradT)
1990 0 : if (s% rsp_flag) then
1991 0 : if (k > 1) then ! Y is superadiabatic gradient
1992 0 : val = s% Y_face(k) + 0.5d0*(s% grada(k-1) + s% grada(k))
1993 : else ! for plotting, use value at k=2
1994 0 : val = s% Y_face(2) + 0.5d0*(s% grada(1) + s% grada(2))
1995 : end if
1996 : end if
1997 : case(p_rsp_Uq)
1998 0 : if (s% rsp_flag) then
1999 0 : if (k > 1) then
2000 0 : val = s% Uq(k)
2001 : else ! for plotting, use value at k=2
2002 0 : val = s% Uq(2)
2003 : end if
2004 : end if
2005 : case(p_rsp_Lr)
2006 0 : if (s% rsp_flag) val = s% Fr(k)*pi4*s% r(k)*s% r(k)
2007 : case(p_rsp_Lr_div_L)
2008 0 : if (s% rsp_flag) val = s% Fr(k)*pi4*s% r(k)*s% r(k)/s% L(k)
2009 : case(p_rsp_Lc)
2010 0 : if (s% rsp_flag) then
2011 0 : val = s% Lc(k)
2012 0 : if (k > 1) then
2013 : val = s% Lc(k)
2014 : else ! for plotting, use value at k=2
2015 0 : val = s% Lc(2)
2016 : end if
2017 : end if
2018 : case(p_rsp_Lc_div_L)
2019 0 : if (s% rsp_flag) then
2020 0 : if (k > 1) then
2021 0 : val = s% Lc(k)/s% L(k)
2022 : else ! for plotting, use value at k=2
2023 0 : val = s% Lc(2)/s% L(2)
2024 : end if
2025 : end if
2026 : case(p_rsp_Lt)
2027 0 : if (s% rsp_flag) then
2028 0 : if (k > 1) then
2029 0 : val = s% Lt(k)
2030 : else ! for plotting, use value at k=2
2031 0 : val = s% Lt(2)
2032 : end if
2033 : end if
2034 : case(p_rsp_Lt_div_L)
2035 0 : if (s% rsp_flag) then
2036 0 : if (k > 1) then
2037 0 : val = s% Lt(k)/s% L(k)
2038 : else ! for plotting, use value at k=2
2039 0 : val = s% Lt(2)/s% L(2)
2040 : end if
2041 : end if
2042 :
2043 : case (p_total_energy) ! specific total energy at k
2044 0 : val = eval_cell_section_total_energy(s,k,k)/s% dm(k)
2045 : case (p_total_energy_sign) ! specific total energy at k
2046 0 : val = eval_cell_section_total_energy(s,k,k)
2047 0 : if (val > 0d0) then
2048 0 : int_val = 1
2049 0 : else if (val < 0d0) then
2050 0 : int_val = -1
2051 : else
2052 0 : int_val = 0
2053 : end if
2054 0 : val = dble(int_val)
2055 0 : int_flag = .true.
2056 : case (p_dwork_dm)
2057 : ! differential work per unit mass per unit mass*time dW/dm
2058 : ! W = dwork_dm*dm*dt
2059 0 : val = s% dwork_dm(k) ! returns (dw/dt)/dm
2060 :
2061 : case (p_cell_specific_IE)
2062 0 : val = s% energy(k)
2063 : case (p_cell_ie_div_star_ie)
2064 0 : val = s% energy(k)*s% dm(k)/s% total_internal_energy_end
2065 : case (p_log_cell_specific_IE)
2066 0 : val = safe_log10(s% energy(k))
2067 : case (p_log_cell_ie_div_star_ie)
2068 0 : val = safe_log10(s% energy(k)*s% dm(k)/s% total_internal_energy_end)
2069 :
2070 : case (p_cell_specific_PE)
2071 0 : val = cell_specific_PE(s,k,d_dlnR00,d_dlnRp1)
2072 :
2073 : case (p_cell_specific_KE)
2074 0 : val = cell_specific_KE(s,k,d_dv00,d_dvp1)
2075 :
2076 : case (p_cell_IE_div_IE_plus_KE)
2077 0 : val = s% energy(k)/(s% energy(k) + cell_specific_KE(s,k,d_dv00,d_dvp1))
2078 :
2079 : case (p_cell_KE_div_IE_plus_KE)
2080 0 : f = cell_specific_KE(s,k,d_dv00,d_dvp1)
2081 0 : val = f/(s% energy(k) + f)
2082 :
2083 : case (p_dlnX_dr)
2084 0 : klo = max(1,k-1)
2085 0 : khi = min(nz,k+1)
2086 : val = log(max(1d-99,max(1d-99,s% X(klo))/max(1d-99,s% X(khi)))) &
2087 0 : / (s% rmid(klo) - s% rmid(khi))
2088 : case (p_dlnY_dr)
2089 0 : klo = max(1,k-1)
2090 0 : khi = min(nz,k+1)
2091 : val = log(max(1d-99,max(1d-99,s% Y(klo))/max(1d-99,s% Y(khi)))) &
2092 0 : / (s% rmid(klo) - s% rmid(khi))
2093 : case (p_dlnRho_dr)
2094 0 : klo = max(1,k-1)
2095 0 : khi = min(nz,k+1)
2096 0 : val = (s% lnd(klo) - s% lnd(khi))/(s% rmid(klo) - s% rmid(khi))
2097 :
2098 : case (p_brunt_B)
2099 0 : if (s% calculate_Brunt_N2) val = s% brunt_B(k)
2100 : case (p_brunt_nonB)
2101 0 : if (s% calculate_Brunt_N2) val = -s% gradT_sub_grada(k)
2102 : case (p_log_brunt_B)
2103 0 : val = log10(max(1d-99,s% brunt_B(k)))
2104 : case (p_log_brunt_nonB)
2105 0 : if (s% calculate_Brunt_N2) val = log10(max(1d-99,-s% gradT_sub_grada(k)))
2106 :
2107 : case (p_brunt_N2)
2108 0 : if (s% calculate_Brunt_N2) val = s% brunt_N2(k)
2109 : case (p_brunt_N2_composition_term)
2110 0 : if (s% calculate_Brunt_N2) val = s% brunt_N2_composition_term(k)
2111 : case (p_brunt_N2_structure_term)
2112 0 : if (s% calculate_Brunt_N2) val = s% brunt_N2(k) - s% brunt_N2_composition_term(k)
2113 : case (p_log_brunt_N2_composition_term)
2114 0 : if (s% calculate_Brunt_N2) val = &
2115 0 : safe_log10(s% brunt_N2_composition_term(k))
2116 : case (p_log_brunt_N2_structure_term)
2117 0 : if (s% calculate_Brunt_N2) val = &
2118 0 : safe_log10(s% brunt_N2(k) - s% brunt_N2_composition_term(k))
2119 :
2120 : case (p_brunt_A)
2121 0 : if (s% calculate_Brunt_N2) val = s% brunt_N2(k)*s% r(k)/s% grav(k)
2122 : case (p_brunt_A_div_x2)
2123 0 : x = s% r(k)/s% r(1)
2124 0 : if (s% calculate_Brunt_N2) val = s% brunt_N2(k)*s% r(k)/s% grav(k)/x/x
2125 : case (p_log_brunt_N2_dimensionless)
2126 0 : if (s% calculate_Brunt_N2) val = &
2127 0 : safe_log10(s% brunt_N2(k)/(3*s% cgrav(1)*s% m_grav(1)/pow3(s% r(1))))
2128 : case (p_brunt_N2_dimensionless)
2129 0 : if (s% calculate_Brunt_N2) val = &
2130 0 : s% brunt_N2(k)/(3*s% cgrav(1)*s% m_grav(1)/pow3(s% r(1)))
2131 : case (p_brunt_N_dimensionless)
2132 0 : if (s% calculate_Brunt_N2) val = &
2133 0 : sqrt(max(0d0,s% brunt_N2(k))/(3*s% cgrav(1)*s% m_grav(1)/pow3(s% r(1))))
2134 : case (p_brunt_N)
2135 0 : if (s% calculate_Brunt_N2) val = sqrt(max(0d0,s% brunt_N2(k)))
2136 : case (p_brunt_frequency) ! cycles per day
2137 0 : if (s% calculate_Brunt_N2) val = &
2138 0 : (secday/(2*pi))*sqrt(max(0d0,s% brunt_N2(k)))
2139 : case (p_log_brunt_N)
2140 0 : if (s% calculate_Brunt_N2) val = safe_log10(sqrt(max(0d0,s% brunt_N2(k))))
2141 : case (p_log_brunt_N2)
2142 0 : if (s% calculate_Brunt_N2) val = safe_log10(s% brunt_N2(k))
2143 :
2144 : case (p_brunt_nu) ! micro Hz
2145 0 : if (s% calculate_Brunt_N2) val = s% brunt_N2(k)
2146 0 : val = (1d6/(2*pi))*sqrt(max(0d0,val))
2147 : case (p_log_brunt_nu) ! micro Hz
2148 0 : if (s% calculate_Brunt_N2) &
2149 0 : val = safe_log10((1d6/(2*pi))*sqrt(max(0d0,s% brunt_N2(k))))
2150 :
2151 : case (p_lamb_S)
2152 0 : val = sqrt(2d0)*s% csound_face(k)/s% r(k) ! for l=1
2153 : case (p_lamb_S2)
2154 0 : val = 2d0*pow2(s% csound_face(k)/s% r(k)) ! for l=1
2155 :
2156 : case (p_lamb_Sl1)
2157 0 : val = (1d6/(2*pi))*sqrt(2d0)*s% csound_face(k)/s% r(k) ! microHz
2158 : case (p_lamb_Sl2)
2159 0 : val = (1d6/(2*pi))*sqrt(6d0)*s% csound_face(k)/s% r(k) ! microHz
2160 : case (p_lamb_Sl3)
2161 0 : val = (1d6/(2*pi))*sqrt(12d0)*s% csound_face(k)/s% r(k) ! microHz
2162 : case (p_lamb_Sl10)
2163 0 : val = (1d6/(2*pi))*sqrt(110d0)*s% csound_face(k)/s% r(k) ! microHz
2164 :
2165 : case (p_log_lamb_Sl1)
2166 0 : val = safe_log10((1d6/(2*pi))*sqrt(2d0)*s% csound_face(k)/s% r(k)) ! microHz
2167 : case (p_log_lamb_Sl2)
2168 0 : val = safe_log10((1d6/(2*pi))*sqrt(6d0)*s% csound_face(k)/s% r(k)) ! microHz
2169 : case (p_log_lamb_Sl3)
2170 0 : val = safe_log10((1d6/(2*pi))*sqrt(12d0)*s% csound_face(k)/s% r(k)) ! microHz
2171 : case (p_log_lamb_Sl10)
2172 0 : val = safe_log10((1d6/(2*pi))*sqrt(110d0)*s% csound_face(k)/s% r(k)) ! microHz
2173 :
2174 : case (p_brunt_N_div_r_integral)
2175 0 : if (s% calculate_Brunt_N2) val = get_brunt_N_div_r_integral(k)
2176 : case (p_sign_brunt_N2)
2177 0 : if (s% calculate_Brunt_N2) val = sign(1d0,s% brunt_N2(k))
2178 :
2179 : case (p_k_r_integral)
2180 0 : if (s% calculate_Brunt_N2) val = get_k_r_integral(k,1,1d0)
2181 :
2182 : case (p_brunt_N2_sub_omega2)
2183 0 : if (s% calculate_Brunt_N2) then
2184 0 : val = s% brunt_N2(k) - pow2(2*pi*s% nu_max/1d6)
2185 0 : if (val > 0d0) then
2186 0 : val = 1
2187 : else
2188 0 : val = 0
2189 : end if
2190 : end if
2191 : case (p_sl2_sub_omega2)
2192 0 : if (s% calculate_Brunt_N2) then
2193 0 : val = 2*pow2(s% csound_face(k)/s% r(k)) - pow2(2*pi*s% nu_max/1d6)
2194 0 : if (val >= 0d0) then
2195 0 : val = 1
2196 : else
2197 0 : val = 0
2198 : end if
2199 : end if
2200 :
2201 : case (p_cs_at_cell_bdy)
2202 0 : val = s% csound_face(k)
2203 : case (p_log_mdot_cs) ! log10(4 Pi r^2 csound rho / (Msun/year))
2204 0 : val = safe_log10(pi4*s% r(k)*s% r(k)*s% csound(k)*s% rho(k)/(Msun/secyer))
2205 : case (p_log_mdot_v) ! log10(4 Pi r^2 v rho / (Msun/year))
2206 0 : if (s% u_flag) then
2207 0 : val = safe_log10(4*pi*s% r(k)*s% r(k)*s% u_face_ad(k)%val*s% rho(k)/(Msun/secyer))
2208 0 : else if (s% v_flag) then
2209 0 : val = safe_log10(pi4*s% r(k)*s% r(k)*s% v(k)*s% rho(k)/(Msun/secyer))
2210 : end if
2211 : case (p_log_L_div_CpTMdot)
2212 0 : if (s% star_mdot == 0) then
2213 : val = 0
2214 : else
2215 0 : val = safe_log10(s% L(k)/(s% cp(k)*s% T(k)*abs(s% star_mdot)*(Msun/secyer)))
2216 : end if
2217 : case (p_logR_kap)
2218 0 : val = s% lnd(k)/ln10 - 3d0*s% lnT(k)/ln10 + 18d0
2219 : case (p_logW)
2220 0 : val = s% lnPgas(k)/ln10 - 4d0*s% lnT(k)/ln10
2221 : case (p_logQ)
2222 0 : val = s% lnd(k)/ln10 - 2d0*s% lnT(k)/ln10 + 12d0
2223 : case (p_logV)
2224 0 : val = s% lnd(k)/ln10 - 0.7d0*s% lnE(k)/ln10 + 20d0
2225 :
2226 : case (p_log_zFe)
2227 : val = 0d0
2228 0 : do j=1,s% species
2229 0 : if (chem_isos% Z(s% chem_id(j)) >= 24) val = val + s% xa(j,k)
2230 : end do
2231 0 : val = safe_log10(val)
2232 : case (p_zFe)
2233 : val = 0d0
2234 0 : do j=1,s% species
2235 0 : if (chem_isos% Z(s% chem_id(j)) >= 24) val = val + s% xa(j,k)
2236 : end do
2237 : case(p_u)
2238 0 : if (s% u_flag) val = s% u(k)
2239 : case(p_u_face)
2240 0 : if (s% u_flag) val = s% u_face_ad(k)%val
2241 : case (p_dPdr_dRhodr_info)
2242 0 : if (s% RTI_flag) val = s% dPdr_dRhodr_info(k)
2243 : case(p_RTI_du_diffusion_kick)
2244 0 : if (s% u_flag) val = s% RTI_du_diffusion_kick(k)
2245 : case(p_log_du_kick_div_du)
2246 0 : if (s% u_flag .and. k > 1) then
2247 0 : if (abs(s% u_face_ad(k)%val) > 1d0) &
2248 0 : val = safe_log10(abs(s% RTI_du_diffusion_kick(k)/s% u_face_ad(k)%val))
2249 : end if
2250 :
2251 : case(p_log_dt_div_tau_conv)
2252 0 : val = safe_log10(s% dt/max(1d-20,conv_time_scale(s,k)))
2253 : case(p_dt_div_tau_conv)
2254 0 : val = s% dt/max(1d-20,conv_time_scale(s,k))
2255 : case(p_tau_conv)
2256 0 : val = conv_time_scale(s,k)
2257 : case(p_tau_qhse)
2258 0 : val = QHSE_time_scale(s,k)
2259 : case(p_tau_epsnuc)
2260 0 : val = eps_nuc_time_scale(s,k)
2261 : case(p_tau_cool)
2262 0 : val = cooling_time_scale(s,k)
2263 :
2264 : case(p_max_abs_xa_corr)
2265 0 : val = s% max_abs_xa_corr(k)
2266 :
2267 : case default
2268 0 : write(*,*) 'FATAL ERROR in profile_getval', c, k
2269 : write(*,*) 'between ' // trim(profile_column_name(c-1)) // ' and ' // &
2270 0 : trim(profile_column_name(c+1)), c-1, c+1
2271 0 : val = 0
2272 31104 : call mesa_error(__FILE__,__LINE__,'profile_getval')
2273 :
2274 : end select
2275 :
2276 : end if
2277 :
2278 :
2279 : contains
2280 :
2281 :
2282 0 : real(dp) function get_L_vel(k) result(v) ! velocity if L carried by convection
2283 : integer, intent(in) :: k
2284 : real(dp) :: rho_face
2285 : integer :: j
2286 0 : if (k == 1) then
2287 0 : j = 2
2288 : else
2289 0 : j = k
2290 : end if
2291 0 : rho_face = interp_val_to_pt(s% rho,j,nz,s% dq,'profile get_L_vel')
2292 0 : v = pow(max(1d0,s% L(k))/(pi4*s% r(k)*s% r(k)*rho_face),one_third)
2293 0 : end function get_L_vel
2294 :
2295 :
2296 0 : real(dp) function get_k_r_integral(k_in, el, nu_factor)
2297 : integer, intent(in) :: k_in
2298 : integer, intent(in) :: el
2299 : real(dp), intent(in) :: nu_factor
2300 : real(dp) :: integral, integral_for_k, &
2301 : cs2, r2, n2, sl2, omega2, L2, kr2, dr
2302 : integer :: k, k1, k_inner, k_outer
2303 : include 'formats'
2304 :
2305 0 : if (k_in == 1) then
2306 : get_k_r_integral = 1
2307 : return
2308 : end if
2309 :
2310 : get_k_r_integral = 0
2311 0 : L2 = el*(el+1)
2312 0 : omega2 = pow2(1d-6*2*pi*s% nu_max*nu_factor)
2313 :
2314 : ! k_inner and k_outer are bounds of evanescent region
2315 :
2316 : ! k_outer is outermost k where Sl2 <= omega2 at k-1 and Sl2 > omega2 at k
2317 : ! 1st find outermost where Sl2 <= omega2
2318 0 : k1 = 0
2319 0 : do k = 2, s% nz
2320 0 : r2 = s% r(k)*s% r(k)
2321 0 : cs2 = s% csound_face(k)*s% csound_face(k)
2322 0 : sl2 = L2*cs2/r2
2323 0 : if (sl2 <= omega2) then
2324 : k1 = k; exit
2325 : end if
2326 : end do
2327 0 : if (k1 == 0) return
2328 : ! then find next k where Sl2 >= omega2
2329 0 : k_outer = 0
2330 0 : do k = k1+1, s% nz
2331 0 : r2 = s% r(k)*s% r(k)
2332 0 : cs2 = s% csound_face(k)*s% csound_face(k)
2333 0 : sl2 = L2*cs2/r2
2334 0 : if (sl2 > omega2) then
2335 : k_outer = k; exit
2336 : end if
2337 : end do
2338 0 : if (k_outer == 0) return
2339 0 : if (k_in <= k_outer) then
2340 : get_k_r_integral = 1
2341 : return
2342 : end if
2343 :
2344 : ! k_inner is next k where N2 >= omega2 at k+1 and N2 < omega2 at k
2345 0 : k_inner = 0
2346 0 : do k = k_outer+1, s% nz
2347 0 : if (s% brunt_N2(k) >= omega2) then
2348 : k_inner= k; exit
2349 : end if
2350 : end do
2351 0 : if (k_inner == 0) return
2352 0 : if (k_in > k_inner) then
2353 : get_k_r_integral = 1
2354 : return
2355 : end if
2356 :
2357 0 : integral = 0; integral_for_k = 0
2358 : get_k_r_integral = 0
2359 0 : do k = k_inner, k_outer, -1
2360 0 : r2 = s% r(k)*s% r(k)
2361 0 : cs2 = s% csound_face(k)*s% csound_face(k)
2362 0 : n2 = s% brunt_N2(k)
2363 0 : sl2 = L2*cs2/r2
2364 0 : kr2 = (1 - n2/omega2)*(1 - Sl2/omega2)/cs2
2365 0 : dr = s% rmid(k-1) - s% rmid(k)
2366 0 : if (kr2 < 0 .and. omega2 < Sl2) integral = integral + sqrt(-kr2)*dr
2367 0 : if (k == k_in) integral_for_k = integral
2368 : end do
2369 0 : if (integral < 1d-99) return
2370 0 : get_k_r_integral = integral_for_k/integral
2371 :
2372 0 : if (is_bad(get_k_r_integral)) then
2373 0 : write(*,2) 'get_k_r_integral', k_in, integral_for_k, integral
2374 0 : call mesa_error(__FILE__,__LINE__,'get_k_r_integral')
2375 : end if
2376 :
2377 : end function get_k_r_integral
2378 :
2379 :
2380 0 : real(dp) function get_brunt_N_div_r_integral(k_in)
2381 : integer, intent(in) :: k_in
2382 : real(dp) :: integral, integral_for_k, dr
2383 : integer :: k
2384 0 : integral = 0
2385 0 : integral_for_k = 0
2386 0 : get_brunt_N_div_r_integral = 1
2387 0 : if (k_in == 1) return
2388 0 : get_brunt_N_div_r_integral = 0
2389 0 : do k = s% nz, 2, -1
2390 0 : dr = s% rmid(k-1) - s% rmid(k)
2391 0 : if (s% brunt_N2(k) > 0) &
2392 0 : integral = integral + sqrt(s% brunt_N2(k))*dr/s% r(k)
2393 0 : if (k == k_in) integral_for_k = integral
2394 : end do
2395 0 : if (integral < 1d-99) return
2396 0 : get_brunt_N_div_r_integral = integral_for_k/integral
2397 0 : end function get_brunt_N_div_r_integral
2398 :
2399 :
2400 0 : real(dp) function get_dlogX_dlogP(j, k)
2401 : integer, intent(in) :: j, k
2402 : integer :: ii, i
2403 : real(dp) :: x00, xm1, dlogP, dlogX
2404 : include 'formats'
2405 0 : get_dlogx_dlogp = 0
2406 0 : if (k > 1) then
2407 : ii = k
2408 : else
2409 : ii = 2
2410 : end if
2411 0 : i = s% net_iso(j)
2412 0 : if (i == 0) return
2413 0 : x00 = s% xa(i,ii)
2414 0 : xm1 = s% xa(i,ii-1)
2415 0 : if (x00 < 1d-20 .or. xm1 < 1d-20) return
2416 0 : dlogP = (s% lnPeos(ii) - s% lnPeos(ii-1))/ln10
2417 0 : if (dlogP <= 0d0) return
2418 0 : dlogX = log10(x00/xm1)
2419 0 : get_dlogX_dlogP = dlogX/dlogP
2420 0 : end function get_dlogX_dlogP
2421 :
2422 :
2423 0 : real(dp) function get_dlog_eps_dlogP(cat, k)
2424 : integer, intent(in) :: cat, k
2425 : integer :: ii
2426 : real(dp) :: eps, epsm1, dlogP, dlog_eps
2427 0 : get_dlog_eps_dlogP = 0
2428 0 : if (k > 1) then
2429 : ii = k
2430 : else
2431 : ii = 2
2432 : end if
2433 0 : eps = s% eps_nuc_categories(cat,ii)
2434 0 : epsm1 = s% eps_nuc_categories(cat,ii-1)
2435 0 : if (eps < 1d-3 .or. epsm1 < 1d-3) return
2436 0 : dlogP = (s% lnPeos(ii) - s% lnPeos(ii-1))/ln10
2437 0 : if (dlogP <= 0d0) return
2438 0 : dlog_eps = log10(eps/epsm1)
2439 0 : get_dlog_eps_dlogP = dlog_eps/dlogP
2440 0 : end function get_dlog_eps_dlogP
2441 :
2442 :
2443 : real(dp) function pt(v,k)
2444 : integer, intent(in) :: k
2445 : real(dp), pointer :: v(:)
2446 : if (k == 1) then
2447 : pt = v(k)
2448 : else
2449 : pt = (v(k)*s% dq(k-1) + v(k-1)*s% dq(k))/(s% dq(k-1) + s% dq(k))
2450 : end if
2451 : end function pt
2452 :
2453 :
2454 0 : real(dp) function if_rot(v,k, alt)
2455 : real(dp),dimension(:), intent(in) :: v
2456 : integer, intent(in) :: k
2457 : real(dp), optional, intent(in) :: alt
2458 0 : if (s% rotation_flag) then
2459 0 : if_rot = v(k)
2460 : else
2461 0 : if (present(alt)) then
2462 0 : if_rot = alt
2463 : else
2464 : if_rot = 0
2465 : end if
2466 : end if
2467 0 : end function if_rot
2468 :
2469 :
2470 0 : real(dp) function if_rot_ad(v,k, alt)
2471 : type(auto_diff_real_star_order1), dimension(:), pointer :: v
2472 : integer, intent(in) :: k
2473 : real(dp), optional, intent(in) :: alt
2474 0 : if (s% rotation_flag) then
2475 0 : if_rot_ad = v(k)% val
2476 : else
2477 0 : if (present(alt)) then
2478 0 : if_rot_ad = alt
2479 : else
2480 : if_rot_ad = 0
2481 : end if
2482 : end if
2483 0 : end function if_rot_ad
2484 :
2485 : end subroutine getval_for_profile
2486 :
2487 : end module profile_getval
|
