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