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