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 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 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_grav_eff)
1191 0 : int_val = if_rot_ad(s% fp_rot,k, alt=1.0d0)
1192 0 : val = s% dxh_v(k)/s%dt / (int_val * s% cgrav(k) * s% m(k) /(s% r(k)*s% r(k)))
1193 : case (p_dvdt_dPdm)
1194 0 : if (k > 1) val = -pi4*s% r(k)*s% r(k)*(s% Peos(k-1) - s% Peos(k))/s% dm_bar(k)
1195 :
1196 : case (p_dm_eps_grav)
1197 0 : val = s% eps_grav_ad(k)% val*s% dm(k)
1198 : case (p_eps_grav)
1199 0 : val = s% eps_grav_ad(k)% val
1200 :
1201 : case (p_log_xm_div_delta_m)
1202 0 : if(abs(s% dt*s% mstar_dot) > 0) val = safe_log10((s% m(1) - s% m(k))/abs(s% dt*s% mstar_dot))
1203 : case (p_xm_div_delta_m)
1204 0 : if(abs(s% dt*s% mstar_dot) > 0) val = (s% m(1) - s% m(k))/abs(s% dt*s% mstar_dot)
1205 :
1206 : case (p_env_eps_grav)
1207 : val = -s% gradT_sub_grada(k)*s% grav(k)*s% mstar_dot*s% Cp(k)*s% T(k) / &
1208 0 : (pi4*s% r(k)*s% r(k)*s% Peos(k))
1209 :
1210 : case (p_mlt_mixing_type)
1211 0 : int_val = s% mlt_mixing_type(k)
1212 0 : val = dble(int_val)
1213 0 : int_flag = .true.
1214 : case (p_mlt_mixing_length)
1215 0 : val = s% mlt_mixing_length(k)
1216 : case (p_mlt_Gamma)
1217 0 : val = s% mlt_Gamma(k)
1218 : case (p_mlt_Zeta)
1219 0 : if (abs(s% gradr(k) - s% grada_face(k)) > 1d-20) &
1220 0 : val = (s% gradr(k) - s% gradT(k))/(s% gradr(k) - s% grada_face(k))
1221 : case (p_mlt_Pturb)
1222 0 : if (s% mlt_Pturb_factor > 0d0 .and. s% okay_to_set_mlt_vc) then
1223 0 : if (s% mlt_vc_old(k) > 0d0) &
1224 0 : val = s% mlt_Pturb_factor*pow2(s% mlt_vc(k))*get_rho_face_val(s,k)/3d0
1225 : end if
1226 : case (p_grad_density)
1227 0 : val = s% grad_density(k)
1228 : case (p_grad_temperature)
1229 0 : val = s% grad_temperature(k)
1230 :
1231 : case (p_gradL_sub_gradr)
1232 0 : val = s% gradL(k) - s% gradr(k)
1233 : case (p_grada_sub_gradr)
1234 0 : val = s% grada_face(k) - s% gradr(k)
1235 :
1236 : case (p_gradL)
1237 0 : val = s% gradL(k)
1238 : case (p_sch_stable)
1239 0 : if (s% grada(k) > s% gradr(k)) val = 1
1240 : case (p_ledoux_stable)
1241 0 : if (s% gradL(k) > s% gradr(k)) val = 1
1242 :
1243 : case (p_eps_nuc_start)
1244 0 : val = s% eps_nuc_start(k)
1245 :
1246 : case (p_dominant_isoA_for_thermohaline)
1247 0 : int_val = chem_isos% Z_plus_N(s% dominant_iso_for_thermohaline(k))
1248 0 : int_flag = .true.
1249 : case (p_dominant_isoZ_for_thermohaline)
1250 0 : int_val = chem_isos% Z(s% dominant_iso_for_thermohaline(k))
1251 0 : int_flag = .true.
1252 : case (p_gradL_composition_term)
1253 0 : val = s% gradL_composition_term(k)
1254 :
1255 : case (p_log_D_conv)
1256 0 : if (s% mixing_type(k) == convective_mixing) then
1257 0 : val = safe_log10(s% D_mix_non_rotation(k))
1258 : else
1259 0 : val = -99
1260 : end if
1261 : case (p_log_D_leftover)
1262 0 : if (s% mixing_type(k) == leftover_convective_mixing) then
1263 0 : val = safe_log10(s% D_mix_non_rotation(k))
1264 : else
1265 0 : val = -99
1266 : end if
1267 : case (p_log_D_semi)
1268 0 : if (s% mixing_type(k) == semiconvective_mixing) then
1269 0 : val = safe_log10(s% D_mix_non_rotation(k))
1270 : else
1271 0 : val = -99
1272 : end if
1273 : case (p_log_D_ovr)
1274 0 : if (s% mixing_type(k) == overshoot_mixing) then
1275 0 : val = safe_log10(s% D_mix_non_rotation(k))
1276 : else
1277 0 : val = -99
1278 : end if
1279 : case (p_log_D_rayleigh_taylor)
1280 0 : if(s% RTI_flag) then
1281 0 : val = safe_log10(s% eta_RTI(k))
1282 : else
1283 0 : val =-99
1284 : end if
1285 : case (p_log_D_anon)
1286 0 : if (s% mixing_type(k) == anonymous_mixing) then
1287 0 : val = safe_log10(s% D_mix_non_rotation(k))
1288 : else
1289 0 : val = -99
1290 : end if
1291 : case (p_log_D_thrm)
1292 0 : if (s% mixing_type(k) == thermohaline_mixing) then
1293 0 : val = safe_log10(s% D_mix_non_rotation(k))
1294 : else
1295 0 : val = -99
1296 : end if
1297 :
1298 : case (p_log_D_minimum)
1299 0 : if (s% mixing_type(k) == minimum_mixing) then
1300 0 : val = safe_log10(s% D_mix(k))
1301 : else
1302 0 : val = -99
1303 : end if
1304 :
1305 : case (p_log_lambda_RTI_div_Hrho)
1306 0 : if (s% RTI_flag) val = safe_log10( &
1307 0 : sqrt(s% alpha_RTI(k))*s% r(k)/s% rho(k)*abs(s% dRhodr_info(k)))
1308 : case (p_lambda_RTI)
1309 0 : if (s% RTI_flag) val = sqrt(s% alpha_RTI(k))*s% r(k)
1310 : case (p_dPdr_info)
1311 0 : if (s% RTI_flag) val = s% dPdr_info(k)
1312 : case (p_dRhodr_info)
1313 0 : if (s% RTI_flag) val = s% dRhodr_info(k)
1314 :
1315 : case (p_source_plus_alpha_RTI)
1316 0 : if (s% RTI_flag) val = s% source_plus_alpha_RTI(k)
1317 : case (p_log_source_plus_alpha_RTI)
1318 0 : if (s% RTI_flag) val = safe_log10(s% source_plus_alpha_RTI(k))
1319 : case (p_log_source_RTI)
1320 0 : if (s% RTI_flag) val = safe_log10(s% source_plus_alpha_RTI(k))
1321 : case (p_source_minus_alpha_RTI)
1322 0 : if (s% RTI_flag) val = s% source_minus_alpha_RTI(k)
1323 : case (p_log_source_minus_alpha_RTI)
1324 0 : if (s% RTI_flag) val = safe_log10(abs(s% source_minus_alpha_RTI(k)))
1325 :
1326 : case (p_dudt_RTI)
1327 0 : if (s% RTI_flag) val = s% dudt_RTI(k)
1328 : case (p_dedt_RTI)
1329 0 : if (s% RTI_flag) val = s% dedt_RTI(k)
1330 :
1331 : case (p_eta_RTI)
1332 0 : if (s% RTI_flag) val = s% eta_RTI(k)
1333 : case (p_log_eta_RTI)
1334 0 : if (s% RTI_flag) val = safe_log10(abs(s% eta_RTI(k)))
1335 : case (p_boost_for_eta_RTI)
1336 0 : if (s% RTI_flag) val = s% boost_for_eta_RTI(k)
1337 : case (p_log_boost_for_eta_RTI)
1338 0 : if (s% RTI_flag) val = safe_log10(abs(s% boost_for_eta_RTI(k)))
1339 :
1340 : case (p_alpha_RTI)
1341 0 : if (s% RTI_flag) val = s% alpha_RTI(k)
1342 : case (p_log_alpha_RTI)
1343 0 : if (s% RTI_flag) val = safe_log10(s% alpha_RTI(k))
1344 : case (p_log_etamid_RTI)
1345 0 : if (s% RTI_flag) val = safe_log10(s% etamid_RTI(k))
1346 :
1347 : case (p_log_sig_RTI)
1348 0 : if (s% RTI_flag) val = safe_log10(s% sig_RTI(k))
1349 : case (p_log_sigmid_RTI)
1350 0 : if (s% RTI_flag) val = safe_log10(s% sigmid_RTI(k))
1351 :
1352 : case (p_log_D_omega)
1353 0 : if (s% rotation_flag) val = safe_log10(s% D_omega(k))
1354 :
1355 : case (p_log_D_mix_non_rotation)
1356 0 : val = safe_log10(s% D_mix_non_rotation(k))
1357 : case (p_log_D_mix_rotation)
1358 0 : val = safe_log10(s% D_mix(k) - s% D_mix_non_rotation(k))
1359 : case (p_log_D_mix)
1360 0 : val = safe_log10(s% D_mix(k))
1361 : case (p_log_sig_mix)
1362 0 : val = safe_log10(s% sig(k))
1363 : case (p_log_sig_raw_mix)
1364 0 : val = safe_log10(s% sig_raw(k))
1365 :
1366 : case (p_burn_avg_epsnuc)
1367 0 : if (s% op_split_burn) val = s% burn_avg_epsnuc(k)
1368 : case (p_log_burn_avg_epsnuc)
1369 0 : if (s% op_split_burn) &
1370 0 : val = safe_log10(abs(s% burn_avg_epsnuc(k)))
1371 : case (p_burn_num_iters)
1372 0 : if (s% op_split_burn) then
1373 0 : int_val = s% burn_num_iters(k); val = dble(int_val)
1374 : else
1375 : int_val = 0; val = 0
1376 : end if
1377 0 : int_flag = .true.
1378 :
1379 : case (p_conv_vel_div_mlt_vc)
1380 0 : if (s% mlt_vc(k) > 0d0) val = s% conv_vel(k)/s% mlt_vc(k)
1381 :
1382 : case (p_conv_vel)
1383 0 : val = s% conv_vel(k)
1384 : case (p_dt_times_conv_vel_div_mixing_length)
1385 0 : val = s% dt*s% conv_vel(k)/s% mlt_mixing_length(k)
1386 : case (p_log_dt_times_conv_vel_div_mixing_length)
1387 0 : val = safe_log10(s% dt*s% conv_vel(k)/s% mlt_mixing_length(k))
1388 : case (p_log_conv_vel)
1389 0 : val = safe_log10(s% conv_vel(k))
1390 : case (p_conv_vel_div_L_vel)
1391 0 : val = s% conv_vel(k)/max(1d0,get_L_vel(k))
1392 : case (p_conv_vel_div_csound)
1393 0 : val = s% conv_vel(k)/s% csound(k)
1394 : case (p_dvc_dt_TDC_div_g)
1395 0 : val = s%dvc_dt_TDC(k) / s%grav(k)
1396 : case (p_mix_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_mixing_type)
1401 0 : val = dble(s% mixing_type(k))
1402 0 : int_val = s% mixing_type(k)
1403 0 : int_flag = .true.
1404 : case (p_log_mlt_D_mix)
1405 0 : val = safe_log10(s% mlt_D(k))
1406 : case (p_log_t_thermal)
1407 0 : val = safe_log10(s% Cp(k)*s% T(k)*(s% m(1) - s% m(k))/s% L(k))
1408 : case (p_log_cp_T_div_t_sound)
1409 : val = safe_log10( &
1410 0 : s% Cp(k)*s% T(k)/(s% Peos(k)/(s% rho(k)*s% grav(k))/s% csound(k)))
1411 : case (p_log_t_sound)
1412 0 : val = safe_log10(s% Peos(k)/(s% rho(k)*s% grav(k))/s% csound(k))
1413 : case (p_pressure_scale_height)
1414 0 : val = s% Peos(k)/(s% rho(k)*s% grav(k))/Rsun
1415 : case (p_pressure_scale_height_cm)
1416 0 : val = s% Peos(k)/(s% rho(k)*s% grav(k))
1417 : case (p_gradT)
1418 0 : val = s% gradT(k)
1419 : case (p_gradr)
1420 0 : val = s% gradr(k)
1421 : case (p_grada_sub_gradT)
1422 0 : val = s% grada_face(k) - s% gradT(k)
1423 :
1424 : case (p_omega)
1425 0 : val = if_rot(s% omega,k)
1426 :
1427 : case (p_log_omega)
1428 0 : val = safe_log10(if_rot(s% omega,k))
1429 : case (p_log_j_rot)
1430 0 : val = safe_log10(if_rot(s% j_rot,k))
1431 : case (p_log_J_inside)
1432 0 : if (s% rotation_flag) then
1433 0 : val = safe_log10(dot_product(s% j_rot(k:s% nz), s% dm(k:s% nz)))
1434 : else
1435 0 : val = -99.0d0
1436 : end if
1437 : case (p_log_J_div_M53)
1438 0 : if (s% rotation_flag) then
1439 : val = safe_log10(&
1440 : dot_product(s% j_rot(k:s% nz), s% dm(k:s% nz)) * &
1441 0 : 1d-50/pow(s% m(k)/Msun,5d0/3d0))
1442 : else
1443 0 : val = -99.0d0
1444 : end if
1445 :
1446 : case (p_shear)
1447 0 : val = if_rot(s% omega_shear,k)
1448 : case (p_log_abs_shear)
1449 0 : if (s% rotation_flag) then
1450 0 : val = safe_log10(s% omega_shear(k))
1451 0 : if (is_bad(val)) then
1452 0 : write(*,2) 'val', k, val
1453 0 : write(*,2) 's% omega_shear(k)', k, s% omega_shear(k)
1454 0 : call mesa_error(__FILE__,__LINE__,'profile')
1455 : end if
1456 : else
1457 0 : val = -99
1458 : end if
1459 : case (p_log_abs_dlnR_domega)
1460 0 : if (s% rotation_flag) then
1461 0 : val = -safe_log10(s% omega_shear(k))
1462 0 : if (is_bad(val)) then
1463 0 : write(*,2) 'val', k, val
1464 0 : write(*,2) 's% omega_shear(k)', k, s% omega_shear(k)
1465 0 : call mesa_error(__FILE__,__LINE__,'profile')
1466 : end if
1467 : else
1468 0 : val = -99
1469 : end if
1470 : case (p_i_rot)
1471 0 : val = if_rot_ad(s% i_rot,k)
1472 : case (p_j_rot)
1473 0 : val = if_rot(s% j_rot,k)
1474 : case (p_v_rot)
1475 0 : val = if_rot(s% omega,k)*if_rot(s% r_equatorial,k)*1d-5 ! km/sec
1476 : case (p_fp_rot)
1477 0 : val = if_rot_ad(s% fp_rot,k, alt=1.0d0)
1478 : case (p_ft_rot)
1479 0 : val = if_rot_ad(s% ft_rot,k, alt=1.0d0)
1480 : case (p_ft_rot_div_fp_rot)
1481 0 : if(s% rotation_flag) then
1482 0 : val = s% ft_rot(k)% val/s% fp_rot(k)% val
1483 : else
1484 0 : val = 1.0d0
1485 : end if
1486 : case (p_w_div_w_crit_roche)
1487 0 : val = if_rot(s% w_div_w_crit_roche,k)
1488 : case (p_w_div_w_crit_roche2)
1489 0 : val = if_rot(s% xh(s% i_w_div_wc,:),k)
1490 : case (p_log_am_nu_non_rot)
1491 0 : val = safe_log10(if_rot(s% am_nu_non_rot,k))
1492 : case (p_log_am_nu_rot)
1493 0 : val = safe_log10(if_rot(s% am_nu_rot,k))
1494 : case (p_log_am_nu)
1495 0 : val = safe_log10(if_rot(s% am_nu_rot,k) + if_rot(s% am_nu_non_rot,k))
1496 :
1497 : case (p_r_polar)
1498 0 : val = if_rot(s% r_polar,k, alt=s% r(k))/Rsun
1499 : case (p_log_r_polar)
1500 0 : val = safe_log10(if_rot(s% r_polar,k, alt=s% r(k))/Rsun)
1501 : case (p_r_equatorial)
1502 0 : val = if_rot(s% r_equatorial,k, alt=s% r(k))/Rsun
1503 : case (p_log_r_equatorial)
1504 0 : val = safe_log10(if_rot(s% r_equatorial,k, alt=s% r(k))/Rsun)
1505 : case (p_r_e_div_r_p)
1506 0 : if (s% rotation_flag) then
1507 0 : if(s% r_polar(k) > 1) val = s% r_equatorial(k)/s% r_polar(k)
1508 : end if
1509 : case (p_omega_crit)
1510 0 : val = omega_crit(s,k)
1511 : case (p_omega_div_omega_crit)
1512 0 : if (s% rotation_flag) then
1513 0 : val = omega_crit(s,k)
1514 0 : if (val < 1d-50) then
1515 0 : val = 0
1516 : else
1517 0 : val = s% omega(k)/val
1518 : end if
1519 : end if
1520 :
1521 : case (p_eps_phase_separation)
1522 0 : if (s% do_phase_separation .and. s% do_phase_separation_heating) val = s% eps_phase_separation(k)
1523 :
1524 : case (p_eps_WD_sedimentation)
1525 0 : if (s% do_element_diffusion) val = s% eps_WD_sedimentation(k)
1526 : case (p_log_eps_WD_sedimentation)
1527 0 : if (s% do_element_diffusion) val = safe_log10(s% eps_WD_sedimentation(k))
1528 :
1529 : case (p_eps_diffusion)
1530 0 : if (s% do_element_diffusion) val = s% eps_diffusion(k)
1531 : case (p_log_eps_diffusion)
1532 0 : if (s% do_element_diffusion) val = safe_log10(s% eps_diffusion(k))
1533 :
1534 : case (p_e_field)
1535 0 : if (s% do_element_diffusion) val = s% E_field(k)
1536 : case (p_log_e_field)
1537 0 : if (s% do_element_diffusion) val = safe_log10(s% E_field(k))
1538 :
1539 : case (p_g_field_element_diffusion)
1540 0 : if (s% do_element_diffusion) val = s% g_field_element_diffusion(k)
1541 : case (p_log_g_field_element_diffusion)
1542 0 : if (s% do_element_diffusion) &
1543 0 : val = safe_log10(s% g_field_element_diffusion(k))
1544 :
1545 : case (p_eE_div_mg_element_diffusion)
1546 0 : if (s% do_element_diffusion) then
1547 0 : if ( s% g_field_element_diffusion(k) /= 0d0) then
1548 0 : val = qe * s% E_field(k)/(amu * s% g_field_element_diffusion(k))
1549 : else
1550 : val = 0d0
1551 : end if
1552 : end if
1553 : case (p_log_eE_div_mg_element_diffusion)
1554 0 : if (s% do_element_diffusion) &
1555 0 : val = safe_log10(qe * s% E_field(k)/(amu * s% g_field_element_diffusion(k)))
1556 :
1557 : case (p_richardson_number)
1558 0 : val = if_rot(s% richardson_number,k)
1559 : case (p_am_domega_dlnR)
1560 0 : val = if_rot(s% domega_dlnR,k)
1561 :
1562 : case (p_am_log_sig) ! == am_log_sig_omega
1563 0 : val = safe_log10(if_rot(s% am_sig_omega,k))
1564 : case (p_am_log_sig_omega)
1565 0 : val = safe_log10(if_rot(s% am_sig_omega,k))
1566 : case (p_am_log_sig_j)
1567 0 : val = safe_log10(if_rot(s% am_sig_j,k))
1568 :
1569 : case (p_am_log_nu_omega)
1570 0 : val = safe_log10(if_rot(s% am_nu_omega,k))
1571 : case (p_am_log_nu_j)
1572 0 : val = safe_log10(if_rot(s% am_nu_j,k))
1573 :
1574 : case (p_am_log_nu_rot)
1575 0 : val = safe_log10(if_rot(s% am_nu_rot,k))
1576 : case (p_am_log_nu_non_rot)
1577 0 : val = safe_log10(if_rot(s% am_nu_non_rot,k))
1578 :
1579 : case (p_am_log_D_visc)
1580 0 : if (s% am_nu_visc_factor >= 0) then
1581 : f = s% am_nu_visc_factor
1582 : else
1583 0 : f = s% D_visc_factor
1584 : end if
1585 0 : val = safe_log10(am_nu_factor*f*if_rot(s% D_visc,k))
1586 : case (p_am_log_D_DSI)
1587 0 : if (s% am_nu_DSI_factor >= 0) then
1588 : f = s% am_nu_DSI_factor
1589 : else
1590 0 : f = s% D_DSI_factor
1591 : end if
1592 0 : val = safe_log10(am_nu_factor*f*if_rot(s% D_DSI,k))
1593 : case (p_am_log_D_SH)
1594 0 : if (s% am_nu_SH_factor >= 0) then
1595 : f = s% am_nu_SH_factor
1596 : else
1597 0 : f = s% D_SH_factor
1598 : end if
1599 0 : val = safe_log10(am_nu_factor*f*if_rot(s% D_SH,k))
1600 : case (p_am_log_D_SSI)
1601 0 : if (s% am_nu_SSI_factor >= 0) then
1602 : f = s% am_nu_SSI_factor
1603 : else
1604 0 : f = s% D_SSI_factor
1605 : end if
1606 0 : val = safe_log10(am_nu_factor*f*if_rot(s% D_SSI,k))
1607 :
1608 : case (p_am_log_D_ES)
1609 0 : if (s% am_nu_ES_factor >= 0) then
1610 : f = s% am_nu_ES_factor
1611 : else
1612 0 : f = s% D_ES_factor
1613 : end if
1614 0 : val = safe_log10(am_nu_factor*f*if_rot(s% D_ES,k))
1615 : case (p_am_log_D_GSF)
1616 0 : if (s% am_nu_GSF_factor >= 0) then
1617 : f = s% am_nu_GSF_factor
1618 : else
1619 0 : f = s% D_GSF_factor
1620 : end if
1621 0 : val = safe_log10(am_nu_factor*f*if_rot(s% D_GSF,k))
1622 : case (p_am_log_D_ST)
1623 0 : if (s% am_nu_ST_factor >= 0) then
1624 : f = s% am_nu_ST_factor
1625 : else
1626 0 : f = s% D_ST_factor
1627 : end if
1628 0 : val = safe_log10(am_nu_factor*f*if_rot(s% D_ST,k))
1629 : case (p_am_log_nu_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% nu_ST,k))
1636 :
1637 : case (p_dynamo_log_B_r)
1638 0 : val = safe_log10(if_rot(s% dynamo_B_r,k))
1639 : case (p_dynamo_log_B_phi)
1640 0 : val = safe_log10(if_rot(s% dynamo_B_phi,k))
1641 :
1642 : case (p_grada_face)
1643 0 : val = s% grada_face(k)
1644 : case (p_gradr_div_grada)
1645 0 : val = s% gradr(k)/s% grada_face(k)
1646 : case (p_gradr_sub_grada)
1647 0 : val = s% gradr(k) - s% grada_face(k)
1648 : case (p_gradT_sub_a)
1649 0 : val = s% gradT(k) - s% grada_face(k)
1650 : case (p_gradT_sub_grada)
1651 0 : val = s% gradT(k) - s% grada_face(k)
1652 : case (p_gradT_div_grada)
1653 0 : val = s% gradT(k) / s% grada_face(k)
1654 : case (p_gradr_sub_gradT)
1655 0 : val = s% gradr(k) - s% gradT(k)
1656 : case (p_gradT_sub_gradr)
1657 0 : val = s% gradT(k) - s% gradr(k)
1658 :
1659 : case (p_gradT_rel_err)
1660 0 : if (k > 1) then
1661 0 : val = (s% lnT(k-1) - s% lnT(k))/(s% lnPeos(k-1) - s% lnPeos(k))
1662 0 : val = (s% gradT(k) - val)/s% gradT(k)
1663 : end if
1664 :
1665 : case (p_gradT_div_gradr)
1666 0 : if (abs(s% gradr(k)) < 1d-99) then
1667 0 : val = 1d0
1668 : else
1669 0 : val = s% gradT(k) / s% gradr(k)
1670 : end if
1671 : case (p_log_gradT_div_gradr)
1672 0 : if (abs(s% gradr(k)) < 1d-99) then
1673 : val = 0d0
1674 : else
1675 0 : val = safe_log10(s% gradT(k) / s% gradr(k))
1676 : end if
1677 :
1678 : case (p_log_mlt_Gamma)
1679 0 : val = safe_log10(s% mlt_Gamma(k))
1680 : case (p_log_mlt_vc)
1681 0 : val = safe_log10(s% mlt_vc(k))
1682 : case (p_mlt_vc)
1683 0 : val = s% mlt_vc(k)
1684 : case (p_mlt_D)
1685 0 : val = s% mlt_D(k)
1686 : case (p_mlt_gradT)
1687 0 : val = s% mlt_gradT(k)
1688 : case (p_mlt_Y_face)
1689 0 : val = s% Y_face(k)
1690 : case (p_mlt_log_abs_Y)
1691 0 : val = safe_log10(abs(s% Y_face(k)))
1692 : case (p_tdc_num_iters)
1693 0 : int_val = s% tdc_num_iters(k); val = dble(int_val)
1694 0 : int_flag = .true.
1695 : case(p_COUPL)
1696 0 : val = s% COUPL(k)
1697 : case(p_SOURCE)
1698 0 : val = s% SOURCE(k)
1699 : case(p_DAMP)
1700 0 : val = s% DAMP(k)
1701 : case(p_DAMPR)
1702 0 : val = s% DAMPR(k)
1703 :
1704 : case (p_delta_r)
1705 0 : val = s% r(k) - s% r_start(k)
1706 : case (p_delta_L)
1707 0 : val = s% L(k) - s% L_start(k)
1708 : case (p_delta_cell_vol)
1709 0 : if (k == s% nz) then
1710 0 : rp1 = s% R_center
1711 0 : rp1_start = s% R_center_old
1712 : else
1713 0 : rp1 = s% r(k+1)
1714 0 : rp1_start = s% r_start(k+1)
1715 : end if
1716 0 : r00 = s% r(k)
1717 0 : r00_start = s% r_start(k)
1718 0 : dr3 = r00*r00*r00 - rp1*rp1*rp1
1719 0 : dr3_start = r00_start*r00_start*r00_start - rp1_start*rp1_start*rp1_start
1720 0 : val = four_thirds_pi*(dr3 - dr3_start)
1721 : case (p_delta_entropy)
1722 0 : val = s% entropy(k) - exp(s% lnS_start(k))/(avo*kerg)
1723 : case (p_delta_T)
1724 0 : val = s% T(k) - s% T_start(k)
1725 : case (p_delta_rho)
1726 0 : val = s% rho(k) - exp(s% lnd_start(k))
1727 : case (p_delta_eps_nuc)
1728 0 : val = s% eps_nuc(k) - s% eps_nuc_start(k)
1729 : case (p_delta_mu)
1730 0 : val = s% mu(k) - s% mu_start(k)
1731 :
1732 : case (p_cno_div_z)
1733 : cno = s% xa(s% net_iso(ic12),k) + &
1734 0 : s% xa(s% net_iso(in14),k) + s% xa(s% net_iso(io16),k)
1735 0 : z = 1 - (s% xa(s% net_iso(ih1),k) + s% xa(s% net_iso(ihe4),k))
1736 0 : if (z > 1d-50) then
1737 0 : val = cno/z
1738 : else
1739 : val = 0
1740 : end if
1741 : case (p_dE)
1742 0 : val = s% energy(k) - s% energy_start(k)
1743 : case (p_dr)
1744 0 : if (k < s% nz) then
1745 0 : val = s% r(k) - s% r(k+1)
1746 : else
1747 0 : val = s% r(k) - s% R_center
1748 : end if
1749 : case (p_dr_ratio)
1750 0 : if (k == 1 .or. k == s% nz) then
1751 0 : val = 1
1752 : else
1753 0 : val = (s% r(k-1) - s% r(k))/(s% r(k) - s% r(k+1))
1754 : end if
1755 : case (p_dv)
1756 0 : if (.not. s% v_flag) then
1757 : val = 0
1758 0 : else if (k < s% nz) then
1759 0 : val = s% v(k+1) - s% v(k)
1760 : else
1761 0 : val = -s% v(k)
1762 : end if
1763 : case (p_dt_dv_div_dr)
1764 0 : if (.not. s% v_flag) then
1765 : val = 0
1766 0 : else if (k < s% nz) then
1767 0 : val = s% dt*(s% v(k+1) - s% v(k))/(s% r(k) - s% r(k+1))
1768 : else
1769 0 : val = -s% dt*s% v(k)/s% r(k)
1770 : end if
1771 :
1772 : case (p_dlog_h1_dlogP)
1773 0 : val = get_dlogX_dlogP(ih1, k)
1774 : case (p_dlog_he3_dlogP)
1775 0 : val = get_dlogX_dlogP(ihe3, k)
1776 : case (p_dlog_he4_dlogP)
1777 0 : val = get_dlogX_dlogP(ihe4, k)
1778 : case (p_dlog_c12_dlogP)
1779 0 : val = get_dlogX_dlogP(ic12, k)
1780 : case (p_dlog_c13_dlogP)
1781 0 : val = get_dlogX_dlogP(ic13, k)
1782 : case (p_dlog_n14_dlogP)
1783 0 : val = get_dlogX_dlogP(in14, k)
1784 : case (p_dlog_o16_dlogP)
1785 0 : val = get_dlogX_dlogP(io16, k)
1786 : case (p_dlog_ne20_dlogP)
1787 0 : val = get_dlogX_dlogP(ine20, k)
1788 : case (p_dlog_mg24_dlogP)
1789 0 : val = get_dlogX_dlogP(img24, k)
1790 : case (p_dlog_si28_dlogP)
1791 0 : val = get_dlogX_dlogP(isi28, k)
1792 :
1793 : case (p_dlog_pp_dlogP)
1794 0 : val = get_dlog_eps_dlogP(ipp, k)
1795 : case (p_dlog_cno_dlogP)
1796 0 : val = get_dlog_eps_dlogP(icno, k)
1797 : case (p_dlog_3alf_dlogP)
1798 0 : val = get_dlog_eps_dlogP(i3alf, k)
1799 :
1800 : case (p_dlog_burn_c_dlogP)
1801 0 : val = get_dlog_eps_dlogP(i_burn_c, k)
1802 : case (p_dlog_burn_n_dlogP)
1803 0 : val = get_dlog_eps_dlogP(i_burn_n, k)
1804 : case (p_dlog_burn_o_dlogP)
1805 0 : val = get_dlog_eps_dlogP(i_burn_o, k)
1806 :
1807 : case (p_dlog_burn_ne_dlogP)
1808 0 : val = get_dlog_eps_dlogP(i_burn_ne, k)
1809 : case (p_dlog_burn_na_dlogP)
1810 0 : val = get_dlog_eps_dlogP(i_burn_na, k)
1811 : case (p_dlog_burn_mg_dlogP)
1812 0 : val = get_dlog_eps_dlogP(i_burn_mg, k)
1813 :
1814 : case (p_dlog_cc_dlogP)
1815 0 : val = get_dlog_eps_dlogP(icc, k)
1816 : case (p_dlog_co_dlogP)
1817 0 : val = get_dlog_eps_dlogP(ico, k)
1818 : case (p_dlog_oo_dlogP)
1819 0 : val = get_dlog_eps_dlogP(ioo, k)
1820 :
1821 : case (p_dlog_burn_si_dlogP)
1822 0 : val = get_dlog_eps_dlogP(i_burn_si, k)
1823 : case (p_dlog_burn_s_dlogP)
1824 0 : val = get_dlog_eps_dlogP(i_burn_s, k)
1825 : case (p_dlog_burn_ar_dlogP)
1826 0 : val = get_dlog_eps_dlogP(i_burn_ar, k)
1827 : case (p_dlog_burn_ca_dlogP)
1828 0 : val = get_dlog_eps_dlogP(i_burn_ca, k)
1829 : case (p_dlog_burn_ti_dlogP)
1830 0 : val = get_dlog_eps_dlogP(i_burn_ti, k)
1831 : case (p_dlog_burn_cr_dlogP)
1832 0 : val = get_dlog_eps_dlogP(i_burn_cr, k)
1833 : case (p_dlog_burn_fe_dlogP)
1834 0 : val = get_dlog_eps_dlogP(i_burn_fe, k)
1835 : case (p_dlog_pnhe4_dlogP)
1836 0 : val = get_dlog_eps_dlogP(ipnhe4, k)
1837 : case (p_dlog_photo_dlogP)
1838 0 : val = get_dlog_eps_dlogP(iphoto, k)
1839 : case (p_dlog_other_dlogP)
1840 0 : val = get_dlog_eps_dlogP(iother, k)
1841 :
1842 : case(p_d_u_div_rmid)
1843 0 : if (s% u_flag .and. k > 1) &
1844 0 : val = s% u(k-1)/s% rmid(k-1) - s% u(k)/s% rmid(k)
1845 : case(p_d_u_div_rmid_start)
1846 0 : if (s% u_flag .and. k > 1) &
1847 0 : val = s% u(k-1)/s% rmid_start(k-1) - s% u(k)/s% rmid_start(k)
1848 :
1849 : case(p_Ptrb)
1850 0 : if (s% RSP2_flag) then
1851 0 : val = get_etrb(s,k)*s% rho(k)
1852 0 : else if (s% RSP_flag) then
1853 0 : val = s% RSP_Et(k)*s% rho(k)
1854 : end if
1855 : case(p_log_Ptrb)
1856 0 : if (s% RSP2_flag) then
1857 0 : val = safe_log10(get_etrb(s,k)*s% rho(k))
1858 0 : else if (s% RSP_flag) then
1859 0 : val = safe_log10(s% RSP_Et(k)*s% rho(k))
1860 : end if
1861 : case(p_w)
1862 0 : if (s% RSP2_flag) then
1863 0 : val = get_w(s,k)
1864 0 : else if (s% RSP_flag) then
1865 0 : val = s% RSP_w(k)
1866 : else
1867 0 : val = s% mlt_vc(k)/sqrt_2_div_3
1868 : end if
1869 : case(p_log_w)
1870 0 : if (s% RSP2_flag) then
1871 0 : val = get_w(s,k)
1872 0 : else if (s% RSP_flag) then
1873 0 : val = s% RSP_w(k)
1874 : else
1875 0 : val = s% mlt_vc(k)/sqrt_2_div_3
1876 : end if
1877 0 : val = safe_log10(val)
1878 : case(p_etrb)
1879 0 : if (s% RSP2_flag) then
1880 0 : val = get_etrb(s,k)
1881 0 : else if (s% RSP_flag) then
1882 0 : val = s% RSP_Et(k)
1883 : end if
1884 : case(p_log_etrb)
1885 0 : if (s% RSP2_flag) then
1886 0 : val = safe_log10(get_etrb(s,k))
1887 0 : else if (s% RSP_flag) then
1888 0 : val = safe_log10(s% RSP_Et(k))
1889 : end if
1890 : case(p_Pvsc)
1891 0 : if (s% use_Pvsc_art_visc .or. s% RSP_flag) val = s% Pvsc(k)
1892 : case(p_Hp_face)
1893 0 : if (rsp_or_w) val = s% Hp_face(k)
1894 : case(p_Y_face)
1895 0 : if (rsp_or_w) val = s% Y_face(k)
1896 : case(p_PII_face)
1897 0 : if (rsp_or_w) val = s% PII(k)
1898 : case(p_Chi)
1899 0 : val = s% Chi(k)
1900 : case(p_Eq)
1901 0 : val = s% Eq(k)
1902 : case(p_Uq)
1903 0 : val = s% Uq(k)
1904 : case(p_Lr)
1905 0 : val = get_Lrad(s,k)
1906 : case(p_Lr_div_L)
1907 0 : val = get_Lrad(s,k)/s% L(k)
1908 : case(p_Lc)
1909 0 : val = get_Lconv(s,k)
1910 : case(p_Lc_div_L)
1911 0 : val = get_Lconv(s,k)/s% L(k)
1912 : case(p_Lt)
1913 0 : if (rsp_or_w) val = s% Lt(k)
1914 : case(p_Lt_div_L)
1915 0 : if (rsp_or_w) val = s% Lt(k)/s% L(k)
1916 :
1917 :
1918 : case(p_rsp_Et)
1919 0 : if (s% rsp_flag) val = s% RSP_Et(k)
1920 : case(p_rsp_logEt)
1921 0 : if (s% rsp_flag) &
1922 0 : val = safe_log10(s% RSP_Et(k))
1923 : case(p_rsp_Pt)
1924 0 : if (s% rsp_flag) val = s% Ptrb(k)
1925 : case(p_rsp_Eq)
1926 0 : if (s% rsp_flag) val = s% Eq(k)
1927 : case(p_rsp_src_snk)
1928 0 : if (s% rsp_flag) val = s% COUPL(k)
1929 : case(p_rsp_src)
1930 0 : if (s% rsp_flag) val = s% SOURCE(k)
1931 : case(p_rsp_sink)
1932 0 : if (s% rsp_flag) val = s% DAMP(k) + s% DAMPR(k)
1933 : case(p_rsp_damp)
1934 0 : if (s% rsp_flag) val = s% DAMP(k)
1935 : case(p_rsp_dampR)
1936 0 : if (s% rsp_flag) val = s% DAMPR(k)
1937 : case(p_rsp_Hp_face)
1938 0 : if (s% rsp_flag) val = s% Hp_face(k)
1939 : case(p_rsp_Chi)
1940 0 : if (s% rsp_flag) val = s% Chi(k)
1941 : case(p_rsp_Pvsc)
1942 0 : if (s% rsp_flag) val = s% Pvsc(k)
1943 : case(p_rsp_erad)
1944 0 : if (s% rsp_flag) val = s% erad(k)
1945 : case(p_rsp_log_erad)
1946 0 : if (s% rsp_flag) val = safe_log10(s% erad(k))
1947 : case(p_rsp_log_dt_div_heat_exchange_timescale)
1948 0 : if (s% rsp_flag) val = safe_log10(s% dt*clight*s% opacity(k)*s% rho(k))
1949 : case(p_rsp_heat_exchange_timescale)
1950 0 : if (s% rsp_flag) val = 1d0/(clight*s% opacity(k)*s% rho(k))
1951 : case(p_rsp_log_heat_exchange_timescale)
1952 0 : if (s% rsp_flag) &
1953 0 : val = safe_log10(1d0/(clight*s% opacity(k)*s% rho(k)))
1954 : case(p_rsp_Y_face)
1955 0 : if (s% rsp_flag) then
1956 0 : if (k > 1) then
1957 0 : val = s% Y_face(k)
1958 : else ! for plotting, use value at k=2
1959 0 : val = s% Y_face(2)
1960 : end if
1961 : end if
1962 : case(p_rsp_gradT)
1963 0 : if (s% rsp_flag) then
1964 0 : if (k > 1) then ! Y is superadiabatic gradient
1965 0 : val = s% Y_face(k) + 0.5d0*(s% grada(k-1) + s% grada(k))
1966 : else ! for plotting, use value at k=2
1967 0 : val = s% Y_face(2) + 0.5d0*(s% grada(1) + s% grada(2))
1968 : end if
1969 : end if
1970 : case(p_rsp_Uq)
1971 0 : if (s% rsp_flag) then
1972 0 : if (k > 1) then
1973 0 : val = s% Uq(k)
1974 : else ! for plotting, use value at k=2
1975 0 : val = s% Uq(2)
1976 : end if
1977 : end if
1978 : case(p_rsp_Lr)
1979 0 : if (s% rsp_flag) val = s% Fr(k)*pi4*s% r(k)*s% r(k)
1980 : case(p_rsp_Lr_div_L)
1981 0 : if (s% rsp_flag) val = s% Fr(k)*pi4*s% r(k)*s% r(k)/s% L(k)
1982 : case(p_rsp_Lc)
1983 0 : if (s% rsp_flag) then
1984 0 : val = s% Lc(k)
1985 0 : if (k > 1) then
1986 : val = s% Lc(k)
1987 : else ! for plotting, use value at k=2
1988 0 : val = s% Lc(2)
1989 : end if
1990 : end if
1991 : case(p_rsp_Lc_div_L)
1992 0 : if (s% rsp_flag) then
1993 0 : if (k > 1) then
1994 0 : val = s% Lc(k)/s% L(k)
1995 : else ! for plotting, use value at k=2
1996 0 : val = s% Lc(2)/s% L(2)
1997 : end if
1998 : end if
1999 : case(p_rsp_Lt)
2000 0 : if (s% rsp_flag) then
2001 0 : if (k > 1) then
2002 0 : val = s% Lt(k)
2003 : else ! for plotting, use value at k=2
2004 0 : val = s% Lt(2)
2005 : end if
2006 : end if
2007 : case(p_rsp_Lt_div_L)
2008 0 : if (s% rsp_flag) then
2009 0 : if (k > 1) then
2010 0 : val = s% Lt(k)/s% L(k)
2011 : else ! for plotting, use value at k=2
2012 0 : val = s% Lt(2)/s% L(2)
2013 : end if
2014 : end if
2015 :
2016 : case (p_total_energy) ! specific total energy at k
2017 0 : val = eval_cell_section_total_energy(s,k,k)/s% dm(k)
2018 : case (p_total_energy_sign) ! specific total energy at k
2019 0 : val = eval_cell_section_total_energy(s,k,k)
2020 0 : if (val > 0d0) then
2021 0 : int_val = 1
2022 0 : else if (val < 0d0) then
2023 0 : int_val = -1
2024 : else
2025 0 : int_val = 0
2026 : end if
2027 0 : val = dble(int_val)
2028 0 : int_flag = .true.
2029 :
2030 : case (p_cell_specific_IE)
2031 0 : val = s% energy(k)
2032 : case (p_cell_ie_div_star_ie)
2033 0 : val = s% energy(k)*s% dm(k)/s% total_internal_energy_end
2034 : case (p_log_cell_specific_IE)
2035 0 : val = safe_log10(s% energy(k))
2036 : case (p_log_cell_ie_div_star_ie)
2037 0 : val = safe_log10(s% energy(k)*s% dm(k)/s% total_internal_energy_end)
2038 :
2039 : case (p_cell_specific_PE)
2040 0 : val = cell_specific_PE(s,k,d_dlnR00,d_dlnRp1)
2041 :
2042 : case (p_cell_specific_KE)
2043 0 : val = cell_specific_KE(s,k,d_dv00,d_dvp1)
2044 :
2045 : case (p_cell_IE_div_IE_plus_KE)
2046 0 : val = s% energy(k)/(s% energy(k) + cell_specific_KE(s,k,d_dv00,d_dvp1))
2047 :
2048 : case (p_cell_KE_div_IE_plus_KE)
2049 0 : f = cell_specific_KE(s,k,d_dv00,d_dvp1)
2050 0 : val = f/(s% energy(k) + f)
2051 :
2052 : case (p_dlnX_dr)
2053 0 : klo = max(1,k-1)
2054 0 : khi = min(nz,k+1)
2055 : val = log(max(1d-99,max(1d-99,s% X(klo))/max(1d-99,s% X(khi)))) &
2056 0 : / (s% rmid(klo) - s% rmid(khi))
2057 : case (p_dlnY_dr)
2058 0 : klo = max(1,k-1)
2059 0 : khi = min(nz,k+1)
2060 : val = log(max(1d-99,max(1d-99,s% Y(klo))/max(1d-99,s% Y(khi)))) &
2061 0 : / (s% rmid(klo) - s% rmid(khi))
2062 : case (p_dlnRho_dr)
2063 0 : klo = max(1,k-1)
2064 0 : khi = min(nz,k+1)
2065 0 : val = (s% lnd(klo) - s% lnd(khi))/(s% rmid(klo) - s% rmid(khi))
2066 :
2067 : case (p_brunt_B)
2068 0 : if (s% calculate_Brunt_N2) val = s% brunt_B(k)
2069 : case (p_brunt_nonB)
2070 0 : if (s% calculate_Brunt_N2) val = -s% gradT_sub_grada(k)
2071 : case (p_log_brunt_B)
2072 0 : val = log10(max(1d-99,s% brunt_B(k)))
2073 : case (p_log_brunt_nonB)
2074 0 : if (s% calculate_Brunt_N2) val = log10(max(1d-99,-s% gradT_sub_grada(k)))
2075 :
2076 : case (p_brunt_N2)
2077 0 : if (s% calculate_Brunt_N2) val = s% brunt_N2(k)
2078 : case (p_brunt_N2_composition_term)
2079 0 : if (s% calculate_Brunt_N2) val = s% brunt_N2_composition_term(k)
2080 : case (p_brunt_N2_structure_term)
2081 0 : if (s% calculate_Brunt_N2) val = s% brunt_N2(k) - s% brunt_N2_composition_term(k)
2082 : case (p_log_brunt_N2_composition_term)
2083 0 : if (s% calculate_Brunt_N2) val = &
2084 0 : safe_log10(s% brunt_N2_composition_term(k))
2085 : case (p_log_brunt_N2_structure_term)
2086 0 : if (s% calculate_Brunt_N2) val = &
2087 0 : safe_log10(s% brunt_N2(k) - s% brunt_N2_composition_term(k))
2088 :
2089 : case (p_brunt_A)
2090 0 : if (s% calculate_Brunt_N2) val = s% brunt_N2(k)*s% r(k)/s% grav(k)
2091 : case (p_brunt_A_div_x2)
2092 0 : x = s% r(k)/s% r(1)
2093 0 : if (s% calculate_Brunt_N2) val = s% brunt_N2(k)*s% r(k)/s% grav(k)/x/x
2094 : case (p_log_brunt_N2_dimensionless)
2095 0 : if (s% calculate_Brunt_N2) val = &
2096 0 : safe_log10(s% brunt_N2(k)/(3*s% cgrav(1)*s% m_grav(1)/pow3(s% r(1))))
2097 : case (p_brunt_N2_dimensionless)
2098 0 : if (s% calculate_Brunt_N2) val = &
2099 0 : s% brunt_N2(k)/(3*s% cgrav(1)*s% m_grav(1)/pow3(s% r(1)))
2100 : case (p_brunt_N_dimensionless)
2101 0 : if (s% calculate_Brunt_N2) val = &
2102 0 : sqrt(max(0d0,s% brunt_N2(k))/(3*s% cgrav(1)*s% m_grav(1)/pow3(s% r(1))))
2103 : case (p_brunt_N)
2104 0 : if (s% calculate_Brunt_N2) val = sqrt(max(0d0,s% brunt_N2(k)))
2105 : case (p_brunt_frequency) ! cycles per day
2106 0 : if (s% calculate_Brunt_N2) val = &
2107 0 : (secday/(2*pi))*sqrt(max(0d0,s% brunt_N2(k)))
2108 : case (p_log_brunt_N)
2109 0 : if (s% calculate_Brunt_N2) val = safe_log10(sqrt(max(0d0,s% brunt_N2(k))))
2110 : case (p_log_brunt_N2)
2111 0 : if (s% calculate_Brunt_N2) val = safe_log10(s% brunt_N2(k))
2112 :
2113 : case (p_brunt_nu) ! micro Hz
2114 0 : if (s% calculate_Brunt_N2) val = s% brunt_N2(k)
2115 0 : val = (1d6/(2*pi))*sqrt(max(0d0,val))
2116 : case (p_log_brunt_nu) ! micro Hz
2117 0 : if (s% calculate_Brunt_N2) &
2118 0 : val = safe_log10((1d6/(2*pi))*sqrt(max(0d0,s% brunt_N2(k))))
2119 :
2120 : case (p_lamb_S)
2121 0 : val = sqrt(2d0)*s% csound_face(k)/s% r(k) ! for l=1
2122 : case (p_lamb_S2)
2123 0 : val = 2d0*pow2(s% csound_face(k)/s% r(k)) ! for l=1
2124 :
2125 : case (p_lamb_Sl1)
2126 0 : val = (1d6/(2*pi))*sqrt(2d0)*s% csound_face(k)/s% r(k) ! microHz
2127 : case (p_lamb_Sl2)
2128 0 : val = (1d6/(2*pi))*sqrt(6d0)*s% csound_face(k)/s% r(k) ! microHz
2129 : case (p_lamb_Sl3)
2130 0 : val = (1d6/(2*pi))*sqrt(12d0)*s% csound_face(k)/s% r(k) ! microHz
2131 : case (p_lamb_Sl10)
2132 0 : val = (1d6/(2*pi))*sqrt(110d0)*s% csound_face(k)/s% r(k) ! microHz
2133 :
2134 : case (p_log_lamb_Sl1)
2135 0 : val = safe_log10((1d6/(2*pi))*sqrt(2d0)*s% csound_face(k)/s% r(k)) ! microHz
2136 : case (p_log_lamb_Sl2)
2137 0 : val = safe_log10((1d6/(2*pi))*sqrt(6d0)*s% csound_face(k)/s% r(k)) ! microHz
2138 : case (p_log_lamb_Sl3)
2139 0 : val = safe_log10((1d6/(2*pi))*sqrt(12d0)*s% csound_face(k)/s% r(k)) ! microHz
2140 : case (p_log_lamb_Sl10)
2141 0 : val = safe_log10((1d6/(2*pi))*sqrt(110d0)*s% csound_face(k)/s% r(k)) ! microHz
2142 :
2143 : case (p_brunt_N_div_r_integral)
2144 0 : if (s% calculate_Brunt_N2) val = get_brunt_N_div_r_integral(k)
2145 : case (p_sign_brunt_N2)
2146 0 : if (s% calculate_Brunt_N2) val = sign(1d0,s% brunt_N2(k))
2147 :
2148 : case (p_k_r_integral)
2149 0 : if (s% calculate_Brunt_N2) val = get_k_r_integral(k,1,1d0)
2150 :
2151 : case (p_brunt_N2_sub_omega2)
2152 0 : if (s% calculate_Brunt_N2) then
2153 0 : val = s% brunt_N2(k) - pow2(2*pi*s% nu_max/1d6)
2154 0 : if (val > 0d0) then
2155 0 : val = 1
2156 : else
2157 0 : val = 0
2158 : end if
2159 : end if
2160 : case (p_sl2_sub_omega2)
2161 0 : if (s% calculate_Brunt_N2) then
2162 0 : val = 2*pow2(s% csound_face(k)/s% r(k)) - pow2(2*pi*s% nu_max/1d6)
2163 0 : if (val >= 0d0) then
2164 0 : val = 1
2165 : else
2166 0 : val = 0
2167 : end if
2168 : end if
2169 :
2170 : case (p_cs_at_cell_bdy)
2171 0 : val = s% csound_face(k)
2172 : case (p_log_mdot_cs) ! log10(4 Pi r^2 csound rho / (Msun/year))
2173 0 : val = safe_log10(pi4*s% r(k)*s% r(k)*s% csound(k)*s% rho(k)/(Msun/secyer))
2174 : case (p_log_mdot_v) ! log10(4 Pi r^2 v rho / (Msun/year))
2175 0 : if (s% u_flag) then
2176 0 : val = safe_log10(4*pi*s% r(k)*s% r(k)*s% u_face_ad(k)%val*s% rho(k)/(Msun/secyer))
2177 0 : else if (s% v_flag) then
2178 0 : val = safe_log10(pi4*s% r(k)*s% r(k)*s% v(k)*s% rho(k)/(Msun/secyer))
2179 : end if
2180 : case (p_log_L_div_CpTMdot)
2181 0 : if (s% star_mdot == 0) then
2182 : val = 0
2183 : else
2184 0 : val = safe_log10(s% L(k)/(s% cp(k)*s% T(k)*abs(s% star_mdot)*(Msun/secyer)))
2185 : end if
2186 : case (p_logR_kap)
2187 0 : val = s% lnd(k)/ln10 - 3d0*s% lnT(k)/ln10 + 18d0
2188 : case (p_logW)
2189 0 : val = s% lnPgas(k)/ln10 - 4d0*s% lnT(k)/ln10
2190 : case (p_logQ)
2191 0 : val = s% lnd(k)/ln10 - 2d0*s% lnT(k)/ln10 + 12d0
2192 : case (p_logV)
2193 0 : val = s% lnd(k)/ln10 - 0.7d0*s% lnE(k)/ln10 + 20d0
2194 :
2195 : case (p_log_zFe)
2196 : val = 0d0
2197 0 : do j=1,s% species
2198 0 : if (chem_isos% Z(s% chem_id(j)) >= 24) val = val + s% xa(j,k)
2199 : end do
2200 0 : val = safe_log10(val)
2201 : case (p_zFe)
2202 : val = 0d0
2203 0 : do j=1,s% species
2204 0 : if (chem_isos% Z(s% chem_id(j)) >= 24) val = val + s% xa(j,k)
2205 : end do
2206 : case(p_u)
2207 0 : if (s% u_flag) val = s% u(k)
2208 : case(p_u_face)
2209 0 : if (s% u_flag) val = s% u_face_ad(k)%val
2210 : case (p_dPdr_dRhodr_info)
2211 0 : if (s% RTI_flag) val = s% dPdr_dRhodr_info(k)
2212 : case(p_RTI_du_diffusion_kick)
2213 0 : if (s% u_flag) val = s% RTI_du_diffusion_kick(k)
2214 : case(p_log_du_kick_div_du)
2215 0 : if (s% u_flag .and. k > 1) then
2216 0 : if (abs(s% u_face_ad(k)%val) > 1d0) &
2217 0 : val = safe_log10(abs(s% RTI_du_diffusion_kick(k)/s% u_face_ad(k)%val))
2218 : end if
2219 :
2220 : case(p_log_dt_div_tau_conv)
2221 0 : val = safe_log10(s% dt/max(1d-20,conv_time_scale(s,k)))
2222 : case(p_dt_div_tau_conv)
2223 0 : val = s% dt/max(1d-20,conv_time_scale(s,k))
2224 : case(p_tau_conv)
2225 0 : val = conv_time_scale(s,k)
2226 : case(p_tau_qhse)
2227 0 : val = QHSE_time_scale(s,k)
2228 : case(p_tau_epsnuc)
2229 0 : val = eps_nuc_time_scale(s,k)
2230 : case(p_tau_cool)
2231 0 : val = cooling_time_scale(s,k)
2232 :
2233 : case(p_max_abs_xa_corr)
2234 0 : val = s% max_abs_xa_corr(k)
2235 :
2236 : case default
2237 0 : write(*,*) 'FATAL ERROR in profile_getval', c, k
2238 : write(*,*) 'between ' // trim(profile_column_name(c-1)) // ' and ' // &
2239 0 : trim(profile_column_name(c+1)), c-1, c+1
2240 0 : val = 0
2241 31104 : call mesa_error(__FILE__,__LINE__,'profile_getval')
2242 :
2243 : end select
2244 :
2245 : end if
2246 :
2247 :
2248 : contains
2249 :
2250 :
2251 0 : real(dp) function get_L_vel(k) result(v) ! velocity if L carried by convection
2252 : integer, intent(in) :: k
2253 0 : real(dp) :: rho_face
2254 : integer :: j
2255 0 : if (k == 1) then
2256 0 : j = 2
2257 : else
2258 0 : j = k
2259 : end if
2260 0 : rho_face = interp_val_to_pt(s% rho,j,nz,s% dq,'profile get_L_vel')
2261 0 : v = pow(max(1d0,s% L(k))/(pi4*s% r(k)*s% r(k)*rho_face),one_third)
2262 41472 : end function get_L_vel
2263 :
2264 :
2265 0 : real(dp) function get_k_r_integral(k_in, el, nu_factor)
2266 : integer, intent(in) :: k_in
2267 : integer, intent(in) :: el
2268 : real(dp), intent(in) :: nu_factor
2269 0 : real(dp) :: integral, integral_for_k, &
2270 0 : cs2, r2, n2, sl2, omega2, L2, kr2, dr
2271 : integer :: k, k1, k_inner, k_outer
2272 : include 'formats'
2273 :
2274 0 : if (k_in == 1) then
2275 : get_k_r_integral = 1
2276 : return
2277 : end if
2278 :
2279 0 : get_k_r_integral = 0
2280 0 : L2 = el*(el+1)
2281 0 : omega2 = pow2(1d-6*2*pi*s% nu_max*nu_factor)
2282 :
2283 : ! k_inner and k_outer are bounds of evanescent region
2284 :
2285 : ! k_outer is outermost k where Sl2 <= omega2 at k-1 and Sl2 > omega2 at k
2286 : ! 1st find outermost where Sl2 <= omega2
2287 0 : k1 = 0
2288 0 : do k = 2, s% nz
2289 0 : r2 = s% r(k)*s% r(k)
2290 0 : cs2 = s% csound_face(k)*s% csound_face(k)
2291 0 : sl2 = L2*cs2/r2
2292 0 : if (sl2 <= omega2) then
2293 : k1 = k; exit
2294 : end if
2295 : end do
2296 0 : if (k1 == 0) return
2297 : ! then find next k where Sl2 >= omega2
2298 0 : k_outer = 0
2299 0 : do k = k1+1, s% nz
2300 0 : r2 = s% r(k)*s% r(k)
2301 0 : cs2 = s% csound_face(k)*s% csound_face(k)
2302 0 : sl2 = L2*cs2/r2
2303 0 : if (sl2 > omega2) then
2304 : k_outer = k; exit
2305 : end if
2306 : end do
2307 0 : if (k_outer == 0) return
2308 0 : if (k_in <= k_outer) then
2309 : get_k_r_integral = 1
2310 : return
2311 : end if
2312 :
2313 : ! k_inner is next k where N2 >= omega2 at k+1 and N2 < omega2 at k
2314 0 : k_inner = 0
2315 0 : do k = k_outer+1, s% nz
2316 0 : if (s% brunt_N2(k) >= omega2) then
2317 : k_inner= k; exit
2318 : end if
2319 : end do
2320 0 : if (k_inner == 0) return
2321 0 : if (k_in > k_inner) then
2322 : get_k_r_integral = 1
2323 : return
2324 : end if
2325 :
2326 0 : integral = 0; integral_for_k = 0
2327 : get_k_r_integral = 0
2328 0 : do k = k_inner, k_outer, -1
2329 0 : r2 = s% r(k)*s% r(k)
2330 0 : cs2 = s% csound_face(k)*s% csound_face(k)
2331 0 : n2 = s% brunt_N2(k)
2332 0 : sl2 = L2*cs2/r2
2333 0 : kr2 = (1 - n2/omega2)*(1 - Sl2/omega2)/cs2
2334 0 : dr = s% rmid(k-1) - s% rmid(k)
2335 0 : if (kr2 < 0 .and. omega2 < Sl2) integral = integral + sqrt(-kr2)*dr
2336 0 : if (k == k_in) integral_for_k = integral
2337 : end do
2338 0 : if (integral < 1d-99) return
2339 0 : get_k_r_integral = integral_for_k/integral
2340 :
2341 0 : if (is_bad(get_k_r_integral)) then
2342 0 : write(*,2) 'get_k_r_integral', k_in, integral_for_k, integral
2343 0 : call mesa_error(__FILE__,__LINE__,'get_k_r_integral')
2344 : end if
2345 :
2346 : end function get_k_r_integral
2347 :
2348 :
2349 0 : real(dp) function get_brunt_N_div_r_integral(k_in)
2350 : integer, intent(in) :: k_in
2351 0 : real(dp) :: integral, integral_for_k, dr
2352 : integer :: k
2353 0 : integral = 0
2354 0 : integral_for_k = 0
2355 0 : get_brunt_N_div_r_integral = 1
2356 0 : if (k_in == 1) return
2357 0 : get_brunt_N_div_r_integral = 0
2358 0 : do k = s% nz, 2, -1
2359 0 : dr = s% rmid(k-1) - s% rmid(k)
2360 0 : if (s% brunt_N2(k) > 0) &
2361 0 : integral = integral + sqrt(s% brunt_N2(k))*dr/s% r(k)
2362 0 : if (k == k_in) integral_for_k = integral
2363 : end do
2364 0 : if (integral < 1d-99) return
2365 0 : get_brunt_N_div_r_integral = integral_for_k/integral
2366 0 : end function get_brunt_N_div_r_integral
2367 :
2368 :
2369 0 : real(dp) function get_dlogX_dlogP(j, k)
2370 : integer, intent(in) :: j, k
2371 : integer :: ii, i
2372 0 : real(dp) :: x00, xm1, dlogP, dlogX
2373 : include 'formats'
2374 0 : get_dlogx_dlogp = 0
2375 0 : if (k > 1) then
2376 : ii = k
2377 : else
2378 : ii = 2
2379 : end if
2380 0 : i = s% net_iso(j)
2381 0 : if (i == 0) return
2382 0 : x00 = s% xa(i,ii)
2383 0 : xm1 = s% xa(i,ii-1)
2384 0 : if (x00 < 1d-20 .or. xm1 < 1d-20) return
2385 0 : dlogP = (s% lnPeos(ii) - s% lnPeos(ii-1))/ln10
2386 0 : if (dlogP <= 0d0) return
2387 0 : dlogX = log10(x00/xm1)
2388 0 : get_dlogX_dlogP = dlogX/dlogP
2389 0 : end function get_dlogX_dlogP
2390 :
2391 :
2392 0 : real(dp) function get_dlog_eps_dlogP(cat, k)
2393 : integer, intent(in) :: cat, k
2394 : integer :: ii
2395 0 : real(dp) :: eps, epsm1, dlogP, dlog_eps
2396 0 : get_dlog_eps_dlogP = 0
2397 0 : if (k > 1) then
2398 : ii = k
2399 : else
2400 : ii = 2
2401 : end if
2402 0 : eps = s% eps_nuc_categories(cat,ii)
2403 0 : epsm1 = s% eps_nuc_categories(cat,ii-1)
2404 0 : if (eps < 1d-3 .or. epsm1 < 1d-3) return
2405 0 : dlogP = (s% lnPeos(ii) - s% lnPeos(ii-1))/ln10
2406 0 : if (dlogP <= 0d0) return
2407 0 : dlog_eps = log10(eps/epsm1)
2408 0 : get_dlog_eps_dlogP = dlog_eps/dlogP
2409 0 : end function get_dlog_eps_dlogP
2410 :
2411 :
2412 : real(dp) function pt(v,k)
2413 : integer, intent(in) :: k
2414 : real(dp), pointer :: v(:)
2415 : if (k == 1) then
2416 : pt = v(k)
2417 : else
2418 : pt = (v(k)*s% dq(k-1) + v(k-1)*s% dq(k))/(s% dq(k-1) + s% dq(k))
2419 : end if
2420 : end function pt
2421 :
2422 :
2423 0 : real(dp) function if_rot(v,k, alt)
2424 : real(dp),dimension(:), intent(in) :: v
2425 : integer, intent(in) :: k
2426 : real(dp), optional, intent(in) :: alt
2427 0 : if (s% rotation_flag) then
2428 0 : if_rot = v(k)
2429 : else
2430 0 : if (present(alt)) then
2431 0 : if_rot = alt
2432 : else
2433 : if_rot = 0
2434 : end if
2435 : end if
2436 0 : end function if_rot
2437 :
2438 :
2439 0 : real(dp) function if_rot_ad(v,k, alt)
2440 : type(auto_diff_real_star_order1), dimension(:), pointer :: v
2441 : integer, intent(in) :: k
2442 : real(dp), optional, intent(in) :: alt
2443 0 : if (s% rotation_flag) then
2444 0 : if_rot_ad = v(k)% val
2445 : else
2446 0 : if (present(alt)) then
2447 0 : if_rot_ad = alt
2448 : else
2449 : if_rot_ad = 0
2450 : end if
2451 : end if
2452 0 : end function if_rot_ad
2453 :
2454 : end subroutine getval_for_profile
2455 :
2456 : end module profile_getval
|
