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 history
21 :
22 : use star_private_def
23 : use star_history_def
24 : use const_def, only: dp, pi, pi4, ln10, one_third, msun, rsun, lsun, clight, secday, secyer, &
25 : no_mixing, semiconvective_mixing, thermohaline_mixing
26 : use chem_def
27 : use history_specs
28 : use star_utils
29 :
30 : implicit none
31 :
32 : private
33 : public :: get_history_specs, get_history_values, get1_hist_value, &
34 : do_get_data_for_history_columns, write_history_info
35 :
36 : logical, parameter :: open_close_log = .true.
37 :
38 : contains
39 :
40 0 : subroutine do_get_data_for_history_columns(s, ierr)
41 : type (star_info), pointer :: s
42 : integer, intent(out) :: ierr
43 : logical, parameter :: write_flag = .false.
44 0 : call do_history_info(s, write_flag, ierr)
45 0 : end subroutine do_get_data_for_history_columns
46 :
47 :
48 4 : subroutine write_history_info(s, ierr)
49 : type (star_info), pointer :: s
50 : integer, intent(out) :: ierr
51 : logical, parameter :: write_flag = .true.
52 4 : call do_history_info(s, write_flag, ierr)
53 4 : end subroutine write_history_info
54 :
55 :
56 8 : subroutine do_history_info(s, write_flag, ierr)
57 : use utils_lib, only : integer_dict_create_hash, integer_dict_free, mkdir, folder_exists
58 : use chem_def, only : category_name
59 : use math_lib, only : math_backend
60 : use net_def, only : Net_General_Info, get_net_ptr
61 : use colors_def, only: Colors_General_Info, get_colors_ptr
62 : use colors_lib, only: how_many_colors_history_columns, data_for_colors_history_columns
63 : type (star_info), pointer :: s
64 :
65 : logical, intent(in) :: write_flag
66 : integer, intent(out) :: ierr
67 :
68 : character (len = strlen) :: fname, dbl_fmt, int_fmt, txt_fmt
69 : integer :: numcols, io, i, nz, col, j, i0, &
70 : num_extra_cols, num_binary_cols, num_extra_binary_cols, num_colors_cols,num_extra_header_items, n
71 : integer, parameter :: num_epsnuc_out = 12
72 : real(dp) :: &
73 : epsnuc_out(num_epsnuc_out), csound_surf, v_surf, envelope_fraction_left, m_div_h, &
74 : min_m_div_h, max_m_div_h
75 : integer :: mixing_regions, mix_relr_regions, burning_regions, burn_relr_regions, rti_regions
76 4 : integer, pointer :: mixing_type(:), mix_relr_type(:), burning_type(:), burn_relr_type(:), rti_type(:)
77 : character (len = maxlen_history_column_name), pointer, dimension(:) :: &
78 4 : extra_col_names, binary_col_names, extra_binary_col_names, colors_col_names, extra_header_item_names
79 : real(dp), pointer, dimension(:) :: &
80 4 : extra_col_vals, binary_col_vals, extra_binary_col_vals, colors_col_vals, extra_header_item_vals
81 :
82 : character (len = maxlen_history_column_name), pointer :: &
83 4 : names(:) ! (num_history_columns)
84 4 : real(dp), pointer :: vals(:) ! (num_history_columns)
85 4 : logical, pointer :: is_int(:) ! (num_history_columns)
86 : type(net_general_info), pointer :: g => null()
87 : type(colors_general_info), pointer :: colors_settings => null()
88 :
89 : logical :: history_file_exists
90 :
91 : include 'formats'
92 :
93 4 : dbl_fmt = s% star_history_dbl_format
94 4 : int_fmt = s% star_history_int_format
95 4 : txt_fmt = s% star_history_txt_format
96 :
97 4 : ierr = 0
98 :
99 4 : call get_net_ptr(s% net_handle, g, ierr)
100 4 : if(ierr/=0) return
101 :
102 4 : call get_colors_ptr(s% colors_handle, colors_settings, ierr)
103 4 : if(ierr/=0) return
104 :
105 4 : if (.not. associated(s% history_column_spec)) then
106 0 : numcols = 0
107 : else
108 4 : numcols = size(s% history_column_spec, dim = 1)
109 : end if
110 4 : num_extra_cols = s% how_many_extra_history_columns(s% id)
111 4 : if (s% include_binary_history_in_log_file) then
112 0 : num_binary_cols = s% how_many_binary_history_columns(s% binary_id)
113 0 : num_extra_binary_cols = s% how_many_extra_binary_history_columns(s% binary_id)
114 : else
115 4 : num_binary_cols = 0
116 4 : num_extra_binary_cols = 0
117 : end if
118 4 : if (colors_settings% use_colors) then
119 0 : num_colors_cols = how_many_colors_history_columns(s% colors_handle)
120 : else
121 4 : num_colors_cols = 0
122 : end if
123 4 : n = numcols + num_extra_cols + num_binary_cols + num_extra_binary_cols + num_colors_cols
124 4 : if (n == 0) then
125 0 : write(*, *) 'WARNING: do not have any output specified for logs.'
126 0 : ierr = -1
127 0 : return
128 : end if
129 :
130 4 : if (s% number_of_history_columns < 0) then
131 1 : s% number_of_history_columns = n
132 3 : else if (s% number_of_history_columns /= n) then
133 0 : if (associated(s% history_values)) then
134 0 : deallocate(s% history_values)
135 0 : nullify(s% history_values)
136 : end if
137 0 : if (associated(s% history_names)) then
138 0 : deallocate(s% history_names)
139 0 : nullify(s% history_names)
140 : end if
141 0 : if (associated(s% history_value_is_integer)) then
142 0 : deallocate(s% history_value_is_integer)
143 0 : nullify(s% history_value_is_integer)
144 : end if
145 0 : if (associated(s% history_names_dict)) then
146 0 : call integer_dict_free(s% history_names_dict)
147 0 : nullify(s% history_names_dict)
148 : end if
149 0 : s% need_to_set_history_names_etc = .true.
150 0 : s% number_of_history_columns = n
151 : end if
152 :
153 4 : if (.not. associated(s% history_values)) then
154 3 : allocate(s% history_values(n))
155 3 : else if (size(s% history_values, dim = 1) /= n) then
156 0 : ierr = -1
157 0 : write(*, 3) 'bad size s% history_values', &
158 0 : size(s% history_values, dim = 1), n
159 : end if
160 4 : vals => s% history_values
161 :
162 4 : if (.not. associated(s% history_names)) then
163 3 : allocate(s% history_names(n))
164 3 : else if (size(s% history_names, dim = 1) /= n) then
165 0 : ierr = -1
166 0 : write(*, 3) 'bad size s% history_names', &
167 0 : size(s% history_names, dim = 1), n
168 : end if
169 :
170 4 : if (s% need_to_set_history_names_etc) then
171 1 : names => s% history_names
172 : else
173 3 : nullify(names)
174 : end if
175 :
176 4 : if (.not. associated(s% history_value_is_integer)) then
177 3 : allocate(s% history_value_is_integer(n))
178 3 : else if (size(s% history_value_is_integer, dim = 1) /= n) then
179 0 : ierr = -1
180 0 : write(*, 2) 'bad size s% history_value_is_integer', &
181 0 : size(s% history_value_is_integer, dim = 1), n
182 : end if
183 4 : if (s% need_to_set_history_names_etc) then
184 1 : is_int => s% history_value_is_integer
185 : else
186 3 : nullify(is_int)
187 : end if
188 :
189 4 : nz = s% nz
190 4 : nullify(mixing_type)
191 4 : nullify(mix_relr_type)
192 4 : nullify(burning_type)
193 4 : nullify(burn_relr_type)
194 4 : nullify(rti_type)
195 4 : nullify(extra_col_names)
196 4 : nullify(extra_col_vals)
197 4 : nullify(binary_col_names)
198 4 : nullify(binary_col_vals)
199 4 : nullify(extra_binary_col_names)
200 4 : nullify(extra_binary_col_vals)
201 4 : nullify(colors_col_names)
202 4 : nullify(colors_col_vals)
203 :
204 4 : if (num_extra_cols > 0) then
205 : allocate(&
206 0 : extra_col_names(num_extra_cols), extra_col_vals(num_extra_cols), stat = ierr)
207 0 : if (ierr /= 0) then
208 0 : call dealloc
209 0 : return
210 : end if
211 0 : extra_col_names(1:num_extra_cols) = 'unknown'
212 0 : extra_col_vals(1:num_extra_cols) = -1d99
213 : call s% data_for_extra_history_columns(&
214 0 : s% id, num_extra_cols, extra_col_names, extra_col_vals, ierr)
215 0 : if (ierr /= 0) then
216 0 : call dealloc
217 0 : return
218 : end if
219 0 : do i = 1, num_extra_cols
220 0 : if(trim(extra_col_names(i))=='unknown') then
221 0 : write(*, *) "Warning empty history name for extra_history_column ", i
222 : end if
223 : end do
224 : end if
225 :
226 4 : if (num_binary_cols > 0) then
227 : allocate(&
228 0 : binary_col_names(num_binary_cols), binary_col_vals(num_binary_cols), stat = ierr)
229 0 : if (ierr /= 0) then
230 0 : call dealloc
231 0 : return
232 : end if
233 : call s% data_for_binary_history_columns(&
234 0 : s% binary_id, num_binary_cols, binary_col_names, binary_col_vals, ierr)
235 0 : if (ierr /= 0) then
236 0 : call dealloc
237 0 : return
238 : end if
239 : end if
240 :
241 4 : if (num_extra_binary_cols > 0) then
242 : allocate(&
243 0 : extra_binary_col_names(num_extra_binary_cols), extra_binary_col_vals(num_extra_binary_cols), stat = ierr)
244 0 : if (ierr /= 0) then
245 0 : call dealloc
246 0 : return
247 : end if
248 0 : extra_binary_col_names(1:num_extra_binary_cols) = 'unknown'
249 0 : extra_binary_col_vals(1:num_extra_binary_cols) = -1d99
250 : call s% data_for_extra_binary_history_columns(&
251 0 : s% binary_id, num_extra_binary_cols, extra_binary_col_names, extra_binary_col_vals, ierr)
252 0 : if (ierr /= 0) then
253 0 : call dealloc
254 0 : return
255 : end if
256 0 : do i = 1, num_extra_binary_cols
257 0 : if(trim(extra_binary_col_names(i))=='unknown') then
258 0 : write(*, *) "Warning empty history name for extra_binary_history_column ", i
259 : end if
260 : end do
261 : end if
262 :
263 4 : if (num_colors_cols > 0) then
264 : allocate(&
265 0 : colors_col_names(num_colors_cols), colors_col_vals(num_colors_cols), stat = ierr)
266 0 : if (ierr /= 0) then
267 0 : call dealloc
268 0 : return
269 : end if
270 0 : colors_col_names(1:num_colors_cols) = 'unknown'
271 0 : colors_col_vals(1:num_colors_cols) = -1d99
272 :
273 : ! Here we compute [Fe/H], and thencall colors to compute history columns.
274 0 : if (colors_settings% z_over_x_ref <= 0d0) then
275 0 : write(*, *) 'colors error: z_over_x_ref must be positive'
276 0 : ierr = -1
277 0 : call dealloc
278 0 : return
279 : end if
280 :
281 0 : min_m_div_h = minval(colors_settings% lu_meta)
282 0 : max_m_div_h = maxval(colors_settings% lu_meta)
283 :
284 : ! Map the current photospheric Z/X onto the atmosphere table metallicity axis.
285 0 : if (s% X(s% photosphere_cell_k) > 0d0 .and. s% Z(s% photosphere_cell_k) > 0d0) then
286 : m_div_h = log10((s% Z(s% photosphere_cell_k)/s% X(s% photosphere_cell_k)) / &
287 0 : colors_settings% z_over_x_ref)
288 0 : m_div_h = max(min_m_div_h, min(max_m_div_h, m_div_h))
289 : else
290 0 : m_div_h = min_m_div_h
291 : end if
292 :
293 : call data_for_colors_history_columns(s% Teff, s% photosphere_logg, s% photosphere_r * Rsun, &
294 : m_div_h, s% model_number, s% colors_handle, &
295 0 : num_colors_cols, colors_col_names, colors_col_vals, ierr)
296 0 : if (ierr /= 0) then
297 0 : call dealloc
298 0 : return
299 : end if
300 0 : do i = 1, num_colors_cols
301 0 : if(trim(colors_col_names(i))=='unknown') then
302 0 : write(*, *) "Warning empty history name for color_history_column ", i
303 : end if
304 : end do
305 : end if
306 :
307 :
308 4 : i0 = 1
309 4 : if (write_flag .and. (open_close_log .or. s% model_number == -100)) then
310 4 : if(.not. folder_exists(trim(s% log_directory))) call mkdir(trim(s% log_directory))
311 :
312 4 : fname = trim(s% log_directory) // '/' // trim(s% star_history_name)
313 4 : inquire(file = trim(fname), exist = history_file_exists)
314 : if ((.not. history_file_exists) .or. &
315 4 : s% doing_first_model_of_run .or. s% need_to_set_history_names_etc) then
316 1 : ierr = 0
317 1 : if (len_trim(s% star_history_header_name) > 0) then
318 0 : fname = trim(s% log_directory) // '/' // trim(s% star_history_header_name)
319 : end if
320 1 : open(newunit = io, file = trim(fname), action = 'write', iostat = ierr)
321 : else
322 3 : i0 = 3
323 3 : open(newunit = io, file = trim(fname), action = 'write', position = 'append', iostat = ierr)
324 : end if
325 4 : if (ierr /= 0) then
326 0 : write(*, *) 'failed to open ' // trim(fname)
327 0 : call dealloc
328 0 : return
329 : end if
330 : end if
331 :
332 4 : csound_surf = eval_csound(s, 1, ierr)
333 4 : if (ierr /= 0) then
334 0 : call dealloc
335 0 : return
336 : end if
337 :
338 4 : if (s% u_flag) then
339 0 : v_surf = s% u(1)
340 4 : else if (s% v_flag) then
341 0 : v_surf = s% v(1)
342 : else
343 4 : v_surf = 0d0
344 : end if
345 :
346 4 : if (s% initial_mass > s% he_core_mass) then
347 : envelope_fraction_left = &
348 4 : (s% star_mass - s% he_core_mass) / (s% initial_mass - s% he_core_mass)
349 : else
350 0 : envelope_fraction_left = 1
351 : end if
352 :
353 4 : if (write_flag .and. i0 == 1) then ! write values at start of log
354 :
355 1 : num_extra_header_items = s% how_many_extra_history_header_items(s% id)
356 1 : if (num_extra_header_items > 0) then
357 : allocate(&
358 : extra_header_item_names(num_extra_header_items), &
359 0 : extra_header_item_vals(num_extra_header_items), stat = ierr)
360 0 : if (ierr /= 0) then
361 0 : call dealloc
362 0 : return
363 : end if
364 0 : extra_header_item_names(1:num_extra_header_items) = 'unknown'
365 0 : extra_header_item_vals(1:num_extra_header_items) = -1d99
366 : call s% data_for_extra_history_header_items(&
367 : s% id, num_extra_header_items, &
368 0 : extra_header_item_names, extra_header_item_vals, ierr)
369 0 : if (ierr /= 0) then
370 0 : call dealloc
371 0 : return
372 : end if
373 0 : do i = 1, num_extra_header_items
374 0 : if(trim(extra_header_item_names(i))=='unknown') then
375 0 : write(*, *) "Warning empty history name for extra_history_header ", i
376 : end if
377 : end do
378 : end if
379 :
380 4 : do i = 1, 3
381 3 : col = 0
382 3 : call write_string(io, col, i, 'version_number', version_number)
383 3 : call write_string(io, col, i, 'compiler', compiler_name)
384 3 : call write_string(io, col, i, 'build', compiler_version_name)
385 3 : call write_string(io, col, i, 'MESA_SDK_version', mesasdk_version_name)
386 3 : call write_string(io, col, i, 'math_backend', math_backend)
387 3 : call write_string(io, col, i, 'date', date)
388 : !call write_val(io, col, i, 'initial_mass', s% initial_mass)
389 : !call write_val(io, col, i, 'initial_z', s% initial_z)
390 3 : call write_val(io, col, i, 'burn_min1', s% burn_min1)
391 3 : call write_val(io, col, i, 'burn_min2', s% burn_min2)
392 :
393 3 : call write_val(io, col, i, 'msun', msun)
394 3 : call write_val(io, col, i, 'rsun', rsun)
395 3 : call write_val(io, col, i, 'lsun', lsun)
396 :
397 3 : do j = 1, num_extra_header_items
398 : call write_val(io, col, i, &
399 3 : extra_header_item_names(j), extra_header_item_vals(j))
400 : end do
401 :
402 4 : write(io, '(A)')
403 :
404 : end do
405 :
406 1 : write(io, '(A)')
407 :
408 1 : if (num_extra_header_items > 0) &
409 0 : deallocate(extra_header_item_names, extra_header_item_vals)
410 :
411 : end if
412 :
413 : ! count the number of declared log regions in pass2. (= param in history_columns.list)
414 : ! write values for the current regions in pass3.
415 4 : burning_regions = 0
416 4 : mixing_regions = 0
417 4 : mix_relr_regions = 0
418 4 : burn_relr_regions = 0
419 4 : rti_regions = 0
420 :
421 10 : do i = i0, 3 ! add a row to the log
422 :
423 6 : col = 0
424 6 : if (i==3) then
425 :
426 4 : if (write_flag .and. i0 == 1 .and. len_trim(s% star_history_header_name) > 0) then
427 0 : close(io)
428 0 : fname = trim(s% log_directory) // '/' // trim(s% star_history_name)
429 0 : open(newunit = io, file = trim(fname), action = 'write', status = 'replace', iostat = ierr)
430 0 : if (ierr /= 0) then
431 0 : call dealloc
432 0 : return
433 : end if
434 : end if
435 :
436 20 : epsnuc_out(1:4) = s% burn_zone_mass(1:4, 1)
437 20 : epsnuc_out(5:8) = s% burn_zone_mass(1:4, 2)
438 20 : epsnuc_out(9:12) = s% burn_zone_mass(1:4, 3)
439 :
440 4 : rti_regions = count_output_mix_regions(rti_offset)
441 4 : if (rti_regions > 0) then
442 0 : allocate(rti_type(nz), stat = ierr)
443 0 : if (ierr /= 0) exit
444 0 : rti_type = 0
445 0 : if (s% RTI_flag) then
446 0 : where (s% alpha_RTI(1:nz) > 0d0) rti_type = 1
447 : end if
448 : end if
449 :
450 4 : mixing_regions = count_output_mix_regions(mixing_offset)
451 4 : if (mixing_regions > 0) then
452 0 : allocate(mixing_type(nz), stat = ierr)
453 0 : if (ierr /= 0) exit
454 0 : call set_mix_types(mixing_type, mixing_regions)
455 0 : call prune_weak_mixing_regions(mixing_type, mixing_regions)
456 : end if
457 :
458 4 : mix_relr_regions = count_output_mix_regions(mix_relr_offset)
459 4 : if (mix_relr_regions > 0) then
460 0 : allocate(mix_relr_type(nz), stat = ierr)
461 0 : if (ierr /= 0) exit
462 0 : call set_mix_types(mix_relr_type, mix_relr_regions)
463 0 : call prune_weak_mixing_regions(mix_relr_type, mix_relr_regions)
464 : end if
465 :
466 4 : burning_regions = count_output_burn_regions(burning_offset)
467 4 : if (burning_regions > 0) then
468 0 : allocate(burning_type(nz), stat = ierr)
469 0 : if (ierr /= 0) exit
470 0 : call set_burn_types(burning_type, burning_offset)
471 : end if
472 :
473 4 : burn_relr_regions = count_output_burn_regions(burn_relr_offset)
474 4 : if (burn_relr_regions > 0) then
475 0 : allocate(burn_relr_type(nz), stat = ierr)
476 0 : if (ierr /= 0) exit
477 0 : call set_burn_types(burn_relr_type, burn_relr_offset)
478 : end if
479 :
480 : end if
481 :
482 366 : do j = 1, numcols
483 366 : call do_col(i, j)
484 : end do
485 6 : do j = 1, num_extra_cols
486 6 : call do_extra_col(i, j, numcols)
487 : end do
488 6 : do j = 1, num_binary_cols
489 6 : call do_binary_col(i, j, numcols + num_extra_cols)
490 : end do
491 6 : do j = 1, num_extra_binary_cols
492 6 : call do_extra_binary_col(i, j, numcols + num_extra_cols + num_binary_cols)
493 : end do
494 6 : do j = 1, num_colors_cols
495 6 : call do_color_col(i, j, numcols + num_extra_cols + num_binary_cols + num_extra_binary_cols)
496 : end do
497 :
498 10 : if (write_flag) write(io, *)
499 :
500 : end do
501 :
502 4 : if (open_close_log .and. write_flag) close(io)
503 :
504 4 : call dealloc
505 :
506 4 : s% model_number_of_history_values = s% model_number
507 :
508 4 : if (s% need_to_set_history_names_etc) then
509 1 : call integer_dict_create_hash(s% history_names_dict, ierr)
510 : end if
511 :
512 8 : s% need_to_set_history_names_etc = .false.
513 :
514 :
515 : contains
516 :
517 :
518 4 : subroutine dealloc
519 4 : if (associated(mixing_type)) deallocate(mixing_type)
520 4 : if (associated(mix_relr_type)) deallocate(mix_relr_type)
521 4 : if (associated(burning_type)) deallocate(burning_type)
522 4 : if (associated(burn_relr_type)) deallocate(burn_relr_type)
523 4 : if (associated(rti_type)) deallocate(rti_type)
524 4 : if (associated(extra_col_names)) deallocate(extra_col_names)
525 4 : if (associated(extra_col_vals)) deallocate(extra_col_vals)
526 4 : if (associated(binary_col_names)) deallocate(binary_col_names)
527 4 : if (associated(binary_col_vals)) deallocate(binary_col_vals)
528 4 : if (associated(extra_binary_col_names)) deallocate(extra_binary_col_names)
529 4 : if (associated(extra_binary_col_vals)) deallocate(extra_binary_col_vals)
530 4 : if (associated(colors_col_names)) deallocate(colors_col_names)
531 4 : if (associated(colors_col_vals)) deallocate(colors_col_vals)
532 :
533 4 : nullify(mixing_type)
534 4 : nullify(mix_relr_type)
535 4 : nullify(burning_type)
536 4 : nullify(burn_relr_type)
537 4 : nullify(rti_type)
538 4 : nullify(extra_col_names)
539 4 : nullify(extra_col_vals)
540 4 : nullify(binary_col_names)
541 4 : nullify(binary_col_vals)
542 4 : nullify(extra_binary_col_names)
543 4 : nullify(extra_binary_col_vals)
544 4 : nullify(colors_col_names)
545 4 : nullify(colors_col_vals)
546 :
547 4 : end subroutine dealloc
548 :
549 :
550 0 : subroutine do_extra_col(pass, j, col_offset)
551 : integer, intent(in) :: pass, j, col_offset
552 : include 'formats'
553 0 : if (pass == 1) then
554 0 : if (write_flag) write(io, fmt = int_fmt, advance = 'no') j + col_offset
555 0 : else if (pass == 2) then
556 0 : call do_name(j + col_offset, extra_col_names(j))
557 0 : else if (pass == 3) then
558 0 : call do_val(j + col_offset, extra_col_vals(j))
559 : end if
560 0 : end subroutine do_extra_col
561 :
562 :
563 0 : subroutine do_binary_col(pass, j, col_offset)
564 : integer, intent(in) :: pass, j, col_offset
565 0 : if (pass == 1) then
566 0 : if (write_flag) write(io, fmt = int_fmt, advance = 'no') j + col_offset
567 0 : else if (pass == 2) then
568 0 : call do_name(j + col_offset, binary_col_names(j))
569 0 : else if (pass == 3) then
570 0 : call do_val(j + col_offset, binary_col_vals(j))
571 : end if
572 0 : end subroutine do_binary_col
573 :
574 :
575 0 : subroutine do_extra_binary_col(pass, j, col_offset)
576 : integer, intent(in) :: pass, j, col_offset
577 0 : if (pass == 1) then
578 0 : if (write_flag) write(io, fmt = int_fmt, advance = 'no') j + col_offset
579 0 : else if (pass == 2) then
580 0 : call do_name(j + col_offset, extra_binary_col_names(j))
581 0 : else if (pass == 3) then
582 0 : call do_val(j + col_offset, extra_binary_col_vals(j))
583 : end if
584 0 : end subroutine do_extra_binary_col
585 :
586 :
587 0 : subroutine do_color_col(pass, j, col_offset)
588 : integer, intent(in) :: pass, j, col_offset
589 0 : if (pass == 1) then
590 0 : if (write_flag) write(io, fmt = int_fmt, advance = 'no') j + col_offset
591 0 : else if (pass == 2) then
592 0 : call do_name(j + col_offset, colors_col_names(j))
593 0 : else if (pass == 3) then
594 0 : call do_val(j + col_offset, colors_col_vals(j))
595 : end if
596 0 : end subroutine do_color_col
597 :
598 :
599 60 : subroutine do_name(j, col_name)
600 : use utils_lib, only : integer_dict_define
601 : integer, intent(in) :: j
602 : integer :: ierr
603 : character (len = *), intent(in) :: col_name
604 60 : if (write_flag) write(io, fmt = txt_fmt, advance = 'no') trim(col_name)
605 60 : if (associated(names)) names(j) = trim(col_name)
606 60 : if (s% need_to_set_history_names_etc) then
607 60 : call integer_dict_define(s% history_names_dict, col_name, j, ierr)
608 60 : if (ierr /= 0) write(*, *) 'failed in integer_dict_define ' // trim(col_name)
609 : end if
610 60 : end subroutine do_name
611 :
612 :
613 8 : integer function count_output_burn_regions(b_regions)
614 : integer, intent(in) :: b_regions
615 : integer :: j, cnt, c, i, ii
616 8 : cnt = 0
617 488 : do j = 1, numcols
618 480 : c = s% history_column_spec(j)
619 488 : if (c > b_regions .and. c <= b_regions + idel) then
620 0 : i = c - b_regions
621 0 : ii = (i + 1) / 2
622 480 : if (ii > cnt) cnt = ii
623 : end if
624 : end do
625 8 : count_output_burn_regions = cnt
626 8 : end function count_output_burn_regions
627 :
628 :
629 0 : subroutine set_burn_types(b_type, b_regions)
630 : integer :: b_type(:), b_regions
631 : integer :: cnt, min_ktop, min_kbot, imax, i, prev_cnt, k
632 : real(dp) :: val
633 : logical, parameter :: dbg = .false.
634 0 : do k = 1, nz
635 0 : val = s% eps_nuc(k) - s% non_nuc_neu(k)
636 0 : b_type(k) = int(sign(1d0, val) * max(0d0, 1d0 + safe_log10(abs(val))))
637 : end do
638 : ! remove smallest regions until <= b_regions remain
639 : imax = nz
640 0 : do i = 1, imax
641 0 : call count_regions(b_type, cnt, min_ktop, min_kbot)
642 : if (dbg) write(*, *) 'count_regions', cnt
643 0 : if (cnt <= b_regions) return
644 0 : if (i > 1 .and. cnt >= prev_cnt) then
645 0 : write(*, *) 'bug in set_burn_types: cnt, prev_cnt', cnt, prev_cnt
646 : if (dbg) call mesa_error(__FILE__, __LINE__, 'debug: set_burn_types')
647 0 : return
648 : end if
649 0 : prev_cnt = cnt
650 : if (dbg) write(*, *) 'remove_region', min_ktop, min_kbot, cnt
651 0 : call remove_region(b_type, min_ktop, min_kbot)
652 : end do
653 : if (dbg) call mesa_error(__FILE__, __LINE__, 'debug: set_burn_types')
654 : end subroutine set_burn_types
655 :
656 :
657 12 : integer function count_output_mix_regions(mx_offset)
658 : integer, intent(in) :: mx_offset
659 : integer :: j, cnt, c, i, ii
660 12 : cnt = 0
661 732 : do j = 1, numcols
662 720 : c = s% history_column_spec(j)
663 732 : if (c > mx_offset .and. c < mx_offset + idel) then
664 0 : i = c - mx_offset
665 0 : ii = (i + 1) / 2
666 720 : if (ii > cnt) cnt = ii
667 : end if
668 : end do
669 12 : count_output_mix_regions = cnt
670 12 : end function count_output_mix_regions
671 :
672 :
673 0 : subroutine prune_weak_mixing_regions(mx_type, mx_regions)
674 : integer :: mx_type(:), mx_regions
675 : real(dp) :: D_max_in_region, D_cutoff
676 : integer :: min_ktop, min_kbot
677 : integer :: i
678 : logical, parameter :: dbg = .false.
679 : include 'formats'
680 0 : D_cutoff = s% mixing_D_limit_for_log
681 0 : do i = 1, 1000
682 : call find_weakest_mixing_region(&
683 0 : mx_type, mx_regions, D_max_in_region, min_ktop, min_kbot)
684 0 : if (D_max_in_region > D_cutoff) then
685 : if (dbg) write(*, 3) 'done', min_ktop, min_kbot, D_max_in_region, D_cutoff
686 0 : return
687 : end if
688 0 : if (mx_type(min_ktop) == no_mixing) then
689 : if (dbg) write(*, 3) 'no mixing regions left'
690 : return
691 : end if
692 : if (dbg) write(*, 3) 'prune_weak_mixing_regions', min_ktop, min_kbot, D_max_in_region
693 : if (dbg) write(*, 3) 'i model mass', &
694 : i, s% model_number, s% m(min_kbot) / Msun, s% m(min_ktop) / Msun
695 : if (dbg) write(*, *)
696 0 : mx_type(min_ktop:min_kbot) = no_mixing
697 : end do
698 : end subroutine prune_weak_mixing_regions
699 :
700 :
701 0 : subroutine find_weakest_mixing_region(&
702 0 : mx_type, mx_regions, D_max_in_region, min_ktop, min_kbot)
703 : integer :: mx_type(:), mx_regions
704 : real(dp), intent(out) :: D_max_in_region
705 : integer, intent(out) :: min_ktop, min_kbot
706 :
707 : integer :: k, kbot, cur_type
708 : real(dp) :: D_max
709 : logical, parameter :: dbg = .false.
710 : include 'formats'
711 0 : min_ktop = 1
712 0 : min_kbot = nz
713 0 : D_max_in_region = 1d99
714 0 : kbot = nz
715 0 : cur_type = mx_type(nz)
716 0 : do k = nz - 1, 1, -1
717 0 : if (cur_type == mx_type(k)) cycle
718 : ! k is bottom of new region; k+1 is top of region cnt
719 : if (cur_type /= no_mixing &
720 : .and. cur_type /= thermohaline_mixing &
721 0 : .and. cur_type /= semiconvective_mixing) then
722 0 : D_max = maxval(s% D_mix_non_rotation(k + 1:kbot))
723 0 : if (D_max < D_max_in_region) then
724 0 : D_max_in_region = D_max; min_ktop = k + 1; min_kbot = kbot
725 : end if
726 : end if
727 : kbot = k
728 0 : cur_type = mx_type(k)
729 : end do
730 0 : if (cur_type == no_mixing) then
731 0 : D_max = maxval(s% D_mix_non_rotation(1:kbot))
732 0 : if (D_max < D_max_in_region) then
733 0 : D_max_in_region = D_max; min_ktop = 1; min_kbot = kbot
734 : end if
735 : end if
736 :
737 0 : end subroutine find_weakest_mixing_region
738 :
739 :
740 0 : subroutine set_mix_types(mx_type, mx_regions)
741 : integer :: mx_type(:), mx_regions
742 :
743 : integer :: cnt, min_ktop, min_kbot, imax, i, prev_cnt, k
744 : logical, parameter :: dbg = .false.
745 :
746 0 : do k = 1, nz
747 0 : mx_type(k) = s% mixing_type(k)
748 : end do
749 : ! remove smallest regions until <= mixing_regions remain
750 : imax = nz
751 0 : do i = 1, imax
752 0 : call count_regions(mx_type, cnt, min_ktop, min_kbot)
753 : if (dbg) write(*, *) 'count_regions', cnt
754 0 : if (cnt <= mx_regions) exit
755 0 : if (i > 1 .and. cnt >= prev_cnt) then
756 0 : write(*, *) 'bug in set_mix_types: cnt, prev_cnt', cnt, prev_cnt
757 : if (dbg) call mesa_error(__FILE__, __LINE__, 'set_mix_types')
758 0 : return
759 : end if
760 0 : prev_cnt = cnt
761 : if (dbg) write(*, *) 'remove_region', min_ktop, min_kbot, cnt
762 0 : call remove_region(mx_type, min_ktop, min_kbot)
763 : end do
764 : if (dbg) call show_regions(mx_type)
765 : end subroutine set_mix_types
766 :
767 :
768 0 : subroutine count_regions(mx_type, cnt, min_ktop, min_kbot)
769 : integer, intent(in) :: mx_type(:)
770 : integer, intent(out) :: cnt, min_ktop, min_kbot
771 : integer :: k, kbot, cur_type
772 : real(dp) :: prev_qtop, qtop, dq, min_dq
773 : logical, parameter :: dbg = .false.
774 : include 'formats'
775 0 : cnt = 1
776 0 : min_ktop = 1
777 0 : min_kbot = nz
778 0 : prev_qtop = 0
779 0 : min_dq = 1
780 0 : kbot = nz
781 0 : cur_type = mx_type(nz)
782 0 : do k = nz - 1, 1, -1
783 0 : if (cur_type == mx_type(k)) cycle
784 : ! k is bottom of new region; k+1 is top of region cnt
785 0 : qtop = s% q(k + 1)
786 0 : dq = qtop - prev_qtop
787 0 : if (dq < min_dq) then
788 0 : min_dq = dq; min_ktop = k + 1; min_kbot = kbot
789 : if (dbg) write(*, 3) &
790 : 'loop min_ktop min_kbot min_dq', min_ktop, min_kbot, min_dq
791 : end if
792 0 : kbot = k
793 0 : cnt = cnt + 1
794 0 : cur_type = mx_type(k)
795 0 : prev_qtop = qtop
796 : end do
797 0 : qtop = 1
798 0 : dq = qtop - prev_qtop
799 0 : if (dq < min_dq) then
800 0 : min_dq = dq; min_ktop = 1; min_kbot = kbot
801 : if (dbg) write(*, 3) &
802 : 'final min_ktop min_kbot min_dq', min_ktop, min_kbot, min_dq
803 : end if
804 0 : end subroutine count_regions
805 :
806 :
807 : subroutine show_regions(mx_type)
808 : integer, intent(in) :: mx_type(:)
809 : integer :: k, kbot, cur_type
810 : include 'formats'
811 : kbot = nz
812 : cur_type = mx_type(nz)
813 : do k = nz - 1, 1, -1
814 : if (cur_type == mx_type(k)) cycle
815 : ! k is bottom of new region; k+1 is top of region cnt
816 : write(*, 4) 'mix region', cur_type, k + 1, kbot, s% m(k + 1) / Msun, s% m(kbot) / Msun
817 : kbot = k
818 : cur_type = mx_type(k)
819 : end do
820 : write(*, 4) 'mix region', cur_type, 1, kbot, s% m(1) / Msun, s% m(kbot) / Msun
821 : end subroutine show_regions
822 :
823 :
824 0 : subroutine remove_region(mx_type, min_ktop, min_kbot)
825 : integer, intent(inout) :: mx_type(:)
826 : integer, intent(in) :: min_ktop, min_kbot
827 : integer :: new_type
828 : logical, parameter :: dbg = .false.
829 : include 'formats'
830 : if (dbg) then
831 : write(*, 2) 'q top', min_ktop, s% q(min_ktop)
832 : write(*, 2) 'q bot', min_kbot, s% q(min_kbot)
833 : write(*, 2) 'dq', min_kbot, s% q(min_ktop) - s% q(min_kbot)
834 : end if
835 0 : if (min_ktop > 1) then ! merge with above
836 0 : new_type = mx_type(min_ktop - 1)
837 0 : else if (min_kbot < nz) then ! merge with below
838 0 : new_type = mx_type(min_kbot + 1)
839 : else
840 0 : write(*, *) 'confusion in args for remove_region', min_ktop, min_kbot
841 : return
842 : end if
843 0 : mx_type(min_ktop:min_kbot) = new_type
844 : if (dbg) write(*, 2) 'new type', min_kbot, new_type
845 : end subroutine remove_region
846 :
847 :
848 0 : integer function region_top(mx_type, j)
849 : integer, intent(in) :: mx_type(:), j
850 : integer :: k, cnt, cur_type
851 0 : cnt = 1
852 0 : cur_type = mx_type(nz)
853 0 : do k = nz - 1, 1, -1
854 0 : if (cur_type == mx_type(k)) cycle
855 : ! k is start of new region; k+1 is top of region cnt
856 0 : if (cnt == j) then
857 0 : region_top = k + 1
858 0 : return
859 : end if
860 0 : cnt = cnt + 1
861 0 : cur_type = mx_type(k)
862 : end do
863 0 : if (cnt == j) then
864 : region_top = 1
865 : else
866 0 : region_top = -1
867 : end if
868 : end function region_top
869 :
870 :
871 0 : real(dp) function interpolate_burn_bdy_q(k) result(val)
872 : use num_lib, only : find0
873 : integer, intent(in) :: k
874 : integer :: bv, bv0, bv1
875 : real(dp) :: eps, eps0, eps1, d0, d1, q0, q1
876 : include 'formats'
877 0 : if (k <= 0) then
878 0 : val = -1; return
879 : end if
880 0 : val = s% q(k)
881 0 : if (k == 1) return
882 0 : eps1 = s% eps_nuc(k) - s% non_nuc_neu(k)
883 0 : bv1 = sign(1d0, eps1) * log10(max(1d0, abs(eps1)))
884 0 : eps0 = s% eps_nuc(k - 1) - s% non_nuc_neu(k - 1)
885 0 : bv0 = sign(1d0, eps0) * log10(max(1d0, abs(eps0)))
886 0 : bv = max(bv0, bv1)
887 0 : eps = pow(10d0, bv)
888 0 : d0 = eps0 - eps
889 0 : d1 = eps1 - eps
890 0 : if (d0 * d1 >= 0) return
891 0 : q1 = s% q(k) - s% dq(k) * 0.5d0
892 0 : q0 = s% q(k - 1) - s% dq(k - 1) * 0.5d0
893 0 : val = find0(q1, d1, q0, d0)
894 0 : end function interpolate_burn_bdy_q
895 :
896 :
897 0 : real(dp) function interpolate_burn_bdy_r(k) result(val)
898 : use num_lib, only : find0
899 : integer, intent(in) :: k
900 : integer :: bv, bv0, bv1
901 : real(dp) :: eps, eps0, eps1, d0, d1
902 : include 'formats'
903 0 : if (k <= 0) then
904 0 : val = -1; return
905 : end if
906 0 : val = s% r(k) / s%r(1)
907 0 : if (k == 1) return
908 : ! Do not subtract s% eps_nuc_neu_total(k) eps_nuc already contains it
909 0 : eps1 = s% eps_nuc(k) - s% non_nuc_neu(k)
910 0 : bv1 = sign(1d0, eps1) * log10(max(1d0, abs(eps1)))
911 : ! Do not subtract s% eps_nuc_neu_total(k) eps_nuc already contains it
912 0 : eps0 = s% eps_nuc(k - 1) - s% non_nuc_neu(k - 1)
913 0 : bv0 = sign(1d0, eps0) * log10(max(1d0, abs(eps0)))
914 0 : bv = max(bv0, bv1)
915 0 : eps = pow(10d0, bv)
916 0 : d0 = eps0 - eps
917 0 : d1 = eps1 - eps
918 0 : if (d0 * d1 >= 0) return
919 0 : val = find0(s%rmid(k), d1, s%rmid(k - 1), d0) / s%r(1)
920 0 : end function interpolate_burn_bdy_r
921 :
922 360 : subroutine do_col(pass, j)
923 : integer, intent(in) :: pass, j
924 360 : if (pass == 1) then
925 60 : call do_col_pass1
926 300 : else if (pass == 2) then
927 60 : call do_col_pass2(j)
928 240 : else if (pass == 3) then
929 240 : call do_col_pass3(s% history_column_spec(j))
930 : end if
931 360 : end subroutine do_col
932 :
933 :
934 60 : subroutine do_col_pass1 ! write the column number
935 60 : col = col + 1
936 60 : if (write_flag) write(io, fmt = int_fmt, advance = 'no') col
937 60 : end subroutine do_col_pass1
938 :
939 :
940 : ! The order of if statements matter, they should be in reverse order
941 : ! to the order in history_specs
942 60 : subroutine do_col_pass2(j) ! get the column name
943 : use colors_lib, only : get_bc_name_by_id
944 : use rates_def, only : reaction_name
945 : integer, intent(in) :: j
946 : character (len = 100) :: col_name
947 : character (len = 10) :: str
948 : integer :: c, i, ii, ir
949 60 : c = s% history_column_spec(j)
950 60 : if (c > rti_offset) then
951 0 : i = c - rti_offset
952 0 : ii = (i + 1) / 2
953 0 : write(str, '(i0)') ii
954 0 : if (mod(i, 2)==1) then
955 0 : col_name = 'rti_type_' // trim(str)
956 : else
957 0 : col_name = 'rti_qtop_' // trim(str)
958 : end if
959 60 : else if (c > burn_relr_offset) then
960 0 : i = c - burn_relr_offset
961 0 : ii = (i + 1) / 2
962 0 : if (ii > burn_relr_regions) burn_relr_regions = ii ! count the regions in pass2
963 0 : if (ii < 10) then
964 0 : write(str, '(i1)') ii
965 0 : else if (ii < 100) then
966 0 : write(str, '(i2)') ii
967 : else
968 0 : write(str, '(i3)') ii
969 : end if
970 0 : if (mod(i, 2)==1) then ! burning type
971 0 : col_name = 'burn_relr_type_' // trim(str)
972 : else ! location of top
973 0 : col_name = 'burn_relr_top_' // trim(str)
974 : end if
975 60 : else if (c > burning_offset) then
976 0 : i = c - burning_offset
977 0 : ii = (i + 1) / 2
978 0 : if (ii > burning_regions) burning_regions = ii ! count the regions in pass2
979 0 : if (ii < 10) then
980 0 : write(str, '(i1)') ii
981 0 : else if (ii < 100) then
982 0 : write(str, '(i2)') ii
983 : else
984 0 : write(str, '(i3)') ii
985 : end if
986 0 : if (mod(i, 2)==1) then ! burning type
987 0 : col_name = 'burn_type_' // trim(str)
988 : else ! location of top
989 0 : col_name = 'burn_qtop_' // trim(str)
990 : end if
991 60 : else if (c > mix_relr_offset) then
992 0 : i = c - mix_relr_offset
993 0 : ii = (i + 1) / 2
994 0 : if (ii > mix_relr_regions) mix_relr_regions = ii ! count the regions in pass2
995 0 : if (ii < 10) then
996 0 : write(str, '(i1)') ii
997 0 : else if (ii < 100) then
998 0 : write(str, '(i2)') ii
999 : else
1000 0 : write(str, '(i3)') ii
1001 : end if
1002 0 : if (mod(i, 2)==1) then
1003 0 : col_name = 'mix_relr_type_' // trim(str)
1004 : else ! location of top
1005 0 : col_name = 'mix_relr_top_' // trim(str)
1006 : end if
1007 60 : else if (c > mixing_offset) then
1008 0 : i = c - mixing_offset
1009 0 : ii = (i + 1) / 2
1010 0 : if (ii > mixing_regions) mixing_regions = ii ! count the regions in pass2
1011 0 : if (ii < 10) then
1012 0 : write(str, '(i1)') ii
1013 0 : else if (ii < 100) then
1014 0 : write(str, '(i2)') ii
1015 : else
1016 0 : write(str, '(i3)') ii
1017 : end if
1018 0 : if (mod(i, 2)==1) then ! mixing type
1019 0 : col_name = 'mix_type_' // trim(str)
1020 : else ! location of top
1021 0 : col_name = 'mix_qtop_' // trim(str)
1022 : end if
1023 60 : else if (c > eps_neu_rate_offset) then
1024 0 : i = c - eps_neu_rate_offset
1025 0 : ir = g % reaction_id(i)
1026 0 : col_name = 'eps_neu_rate_' // trim(reaction_name(ir))
1027 60 : else if (c > eps_nuc_rate_offset) then
1028 0 : i = c - eps_nuc_rate_offset
1029 0 : ir = g % reaction_id(i)
1030 0 : col_name = 'eps_nuc_rate_' // trim(reaction_name(ir))
1031 60 : else if (c > screened_rate_offset) then
1032 0 : i = c - screened_rate_offset
1033 0 : ir = g % reaction_id(i)
1034 0 : col_name = 'screened_rate_' // trim(reaction_name(ir))
1035 60 : else if (c > raw_rate_offset) then
1036 0 : i = c - raw_rate_offset
1037 0 : ir = g % reaction_id(i)
1038 0 : col_name = 'raw_rate_' // trim(reaction_name(ir))
1039 60 : else if (c > log_lum_band_offset) then
1040 0 : i = c - log_lum_band_offset
1041 0 : col_name = 'log_lum_band_' // trim(get_bc_name_by_id(i, ierr))
1042 60 : else if (c > lum_band_offset) then
1043 0 : i = c - lum_band_offset
1044 0 : col_name = 'lum_band_' // trim(get_bc_name_by_id(i, ierr))
1045 60 : else if (c > abs_mag_offset) then
1046 0 : i = c - abs_mag_offset
1047 0 : col_name = 'abs_mag_' // trim(get_bc_name_by_id(i, ierr))
1048 60 : else if (c > bc_offset) then
1049 0 : i = c - bc_offset
1050 0 : col_name = 'bc_' // trim(get_bc_name_by_id(i, ierr))
1051 60 : else if (c > c_log_eps_burn_offset) then
1052 0 : i = c - c_log_eps_burn_offset
1053 0 : col_name = 'c_log_eps_burn_' // trim(category_name(i))
1054 60 : else if (c > max_eps_nuc_offset) then
1055 0 : i = c - max_eps_nuc_offset
1056 0 : col_name = 'max_eps_nuc_log_' // trim(chem_isos% name(i))
1057 60 : else if (c > cz_top_max_offset) then
1058 0 : i = c - cz_top_max_offset
1059 0 : col_name = 'cz_top_log_' // trim(chem_isos% name(i))
1060 60 : else if (c > cz_max_offset) then
1061 0 : i = c - cz_max_offset
1062 0 : col_name = 'cz_log_' // trim(chem_isos% name(i))
1063 60 : else if (c > log_surface_xa_offset) then
1064 0 : i = c - log_surface_xa_offset
1065 0 : col_name = 'log_surface_' // trim(chem_isos% name(i))
1066 60 : else if (c > log_center_xa_offset) then
1067 0 : i = c - log_center_xa_offset
1068 0 : col_name = 'log_center_' // trim(chem_isos% name(i))
1069 60 : else if (c > log_average_xa_offset) then
1070 0 : i = c - log_average_xa_offset
1071 0 : col_name = 'log_average_' // trim(chem_isos% name(i))
1072 60 : else if (c > log_total_mass_offset) then
1073 0 : i = c - log_total_mass_offset
1074 0 : col_name = 'log_total_mass_' // trim(chem_isos% name(i))
1075 60 : else if (c > total_mass_offset) then
1076 2 : i = c - total_mass_offset
1077 2 : col_name = 'total_mass_' // trim(chem_isos% name(i))
1078 58 : else if (c > category_offset) then
1079 3 : i = c - category_offset
1080 3 : col_name = category_name(i)
1081 55 : else if (c > average_xa_offset) then
1082 0 : i = c - average_xa_offset
1083 0 : col_name = 'average_' // trim(chem_isos% name(i))
1084 55 : else if (c > surface_xa_offset) then
1085 2 : i = c - surface_xa_offset
1086 2 : col_name = 'surface_' // trim(chem_isos% name(i))
1087 53 : else if (c > center_xa_offset) then
1088 4 : i = c - center_xa_offset
1089 4 : col_name = 'center_' // trim(chem_isos% name(i))
1090 : else
1091 49 : col_name = trim(history_column_name(c))
1092 : end if
1093 60 : call do_name(j, col_name)
1094 60 : end subroutine do_col_pass2
1095 :
1096 :
1097 240 : subroutine do_col_pass3(c) ! get the column value
1098 : use rates_def
1099 : integer, intent(in) :: c
1100 : integer :: i, ii, k, int_val
1101 : logical :: is_int_val
1102 : real(dp) :: val, frac
1103 240 : int_val = 0; val = 0; is_int_val = .false.
1104 :
1105 240 : if (c > rti_offset) then
1106 0 : i = c - rti_offset
1107 0 : ii = (i + 1) / 2
1108 0 : k = region_top(rti_type, ii)
1109 0 : if (mod(i, 2)==1) then
1110 0 : is_int_val = .true.
1111 0 : if (k > 0) then
1112 0 : int_val = rti_type(k)
1113 : else
1114 0 : int_val = -1
1115 : end if
1116 : else
1117 0 : if (k > 0) then
1118 0 : val = s% q(k) ! outer face of the outermost cell in this RTI region
1119 : else
1120 0 : val = 1d0
1121 : end if
1122 : end if
1123 240 : else if (c > burn_relr_offset) then
1124 0 : i = c - burn_relr_offset
1125 0 : ii = (i + 1) / 2
1126 0 : k = region_top(burn_relr_type, ii)
1127 0 : if (mod(i, 2)==1) then ! burning type
1128 0 : is_int_val = .true.
1129 0 : if (k > 0) then
1130 0 : int_val = burn_relr_type(k)
1131 : else
1132 0 : int_val = -9999
1133 : end if
1134 : else ! location of top
1135 0 : val = interpolate_burn_bdy_r(k)
1136 : end if
1137 240 : else if (c > burning_offset) then
1138 0 : i = c - burning_offset
1139 0 : ii = (i + 1) / 2
1140 0 : k = region_top(burning_type, ii)
1141 0 : if (mod(i, 2)==1) then ! burning type
1142 0 : is_int_val = .true.
1143 0 : if (k > 0) then
1144 0 : int_val = burning_type(k)
1145 : else
1146 0 : int_val = -9999
1147 : end if
1148 : else ! location of top
1149 0 : val = interpolate_burn_bdy_q(k)
1150 : end if
1151 240 : else if (c > mix_relr_offset) then
1152 0 : i = c - mix_relr_offset
1153 0 : ii = (i + 1) / 2
1154 0 : k = region_top(mix_relr_type, ii)
1155 0 : if (mod(i, 2)==1) then ! mixing type
1156 0 : is_int_val = .true.
1157 0 : if (k > 0) then
1158 0 : int_val = mix_relr_type(k)
1159 : else
1160 0 : int_val = -1
1161 : end if
1162 : else ! r/rstar location of boundary
1163 0 : if (k <= 1) then
1164 0 : val = 1d0
1165 : else
1166 0 : frac = s% cz_bdy_dq(k - 1) / s% dq(k - 1)
1167 0 : val = (1d0 - frac) * pow3(s% r(k - 1)) + frac * pow3(s% r(k))
1168 0 : val = pow(val, one_third) / s% r(1)
1169 : end if
1170 : end if
1171 240 : else if (c > mixing_offset) then
1172 0 : i = c - mixing_offset
1173 0 : ii = (i + 1) / 2
1174 0 : k = region_top(mixing_type, ii)
1175 0 : if (mod(i, 2)==1) then ! mixing type
1176 0 : is_int_val = .true.
1177 0 : if (k > 0) then
1178 0 : int_val = mixing_type(k)
1179 : else
1180 0 : int_val = -1
1181 : end if
1182 : else ! q location of boundary
1183 0 : if (k <= 1) then
1184 0 : val = 1d0
1185 : else
1186 0 : val = s% q(k - 1) - s% cz_bdy_dq(k - 1)
1187 : end if
1188 : end if
1189 : else
1190 : call history_getval(&
1191 : s, c, val, int_val, is_int_val, &
1192 240 : nz, v_surf, csound_surf, envelope_fraction_left, epsnuc_out, ierr)
1193 240 : if (ierr /= 0) then
1194 0 : write(*, *) 'missing log info for ' // trim(history_column_name(c)), j, k
1195 0 : int_val = -99999999
1196 0 : is_int_val = .true.
1197 0 : ierr = -1
1198 : end if
1199 : end if
1200 240 : if (is_int_val) then
1201 16 : call do_int_val(j, int_val)
1202 : else
1203 224 : call do_val(j, val)
1204 : end if
1205 240 : end subroutine do_col_pass3
1206 :
1207 :
1208 224 : subroutine do_val(j, val)
1209 : integer, intent(in) :: j
1210 : real(dp), intent(in) :: val
1211 224 : if (write_flag) then
1212 224 : if (is_bad_num(val)) then
1213 0 : write(io, fmt = dbl_fmt, advance = 'no') -1d99
1214 : else
1215 224 : write(io, fmt = dbl_fmt, advance = 'no') val
1216 : end if
1217 : end if
1218 224 : if (associated(vals)) vals(j) = val
1219 224 : if (associated(is_int)) is_int(j) = .false.
1220 224 : end subroutine do_val
1221 :
1222 :
1223 16 : subroutine do_int_val(j, val)
1224 : integer, intent(in) :: j
1225 : integer, intent(in) :: val
1226 16 : if (write_flag) write(io, fmt = int_fmt, advance = 'no') val
1227 16 : if (associated(vals)) vals(j) = dble(val)
1228 16 : if (associated(is_int)) is_int(j) = .true.
1229 16 : end subroutine do_int_val
1230 :
1231 18 : subroutine write_string(io, col, pass, name, val) !for header items only
1232 : integer, intent(in) :: io, pass
1233 : integer, intent(inout) :: col
1234 : character(len = *), intent(in) :: name, val
1235 : character(len = strlen) :: my_val
1236 :
1237 18 : my_val = '"' // trim(val) // '"'
1238 18 : if (pass == 1) then
1239 6 : col = col + 1
1240 6 : write(io, fmt = int_fmt, advance = 'no') col
1241 12 : else if (pass == 2) then
1242 6 : write(io, fmt = txt_fmt, advance = 'no') trim(name)
1243 6 : else if (pass == 3) then
1244 6 : write(io, fmt = txt_fmt, advance = 'no') trim(my_val)
1245 : end if
1246 18 : end subroutine write_string
1247 :
1248 :
1249 : subroutine write_integer(io, col, pass, name, val) ! for header items only
1250 : integer, intent(in) :: io, pass
1251 : integer, intent(inout) :: col
1252 : character (len = *), intent(in) :: name
1253 : integer, intent(in) :: val
1254 : if (pass == 1) then
1255 : col = col + 1
1256 : write(io, fmt = int_fmt, advance = 'no') col
1257 : else if (pass == 2) then
1258 : write(io, fmt = txt_fmt, advance = 'no') trim(name)
1259 : else if (pass == 3) then
1260 : write(io, fmt = int_fmt, advance = 'no') val
1261 : end if
1262 : end subroutine write_integer
1263 :
1264 :
1265 15 : subroutine write_val(io, col, pass, name, val) ! for header items only
1266 : integer, intent(in) :: io, pass
1267 : integer, intent(inout) :: col
1268 : character (len = *), intent(in) :: name
1269 : real(dp), intent(in) :: val
1270 15 : if (pass == 1) then
1271 5 : col = col + 1
1272 5 : write(io, fmt = int_fmt, advance = 'no') col
1273 10 : else if (pass == 2) then
1274 5 : write(io, fmt = txt_fmt, advance = 'no') trim(name)
1275 5 : else if (pass == 3) then
1276 5 : write(io, fmt = dbl_fmt, advance = 'no') val
1277 : end if
1278 15 : end subroutine write_val
1279 :
1280 :
1281 : end subroutine do_history_info
1282 :
1283 240 : subroutine history_getval(&
1284 : s, c, val, int_val, is_int_val, &
1285 240 : nz, v_surf, csound_surf, envelope_fraction_left, epsnuc_out, ierr)
1286 : use colors_lib, only : get_abs_mag_by_id, get_bc_by_id, get_lum_band_by_id
1287 : use chem_lib, only : chem_M_div_h
1288 : use rsp_def, only : rsp_phase_time0
1289 : use gravity_darkening
1290 :
1291 : type (star_info), pointer :: s
1292 : integer, intent(in) :: c, nz
1293 : real(dp), intent(in) :: &
1294 : v_surf, csound_surf, envelope_fraction_left, epsnuc_out(:)
1295 : real(dp), intent(out) :: val
1296 : integer, intent(out) :: int_val
1297 : logical, intent(out) :: is_int_val
1298 : integer, intent(out) :: ierr
1299 :
1300 : integer :: k, i, min_k, k2
1301 : real(dp) :: Ledd, phi_Joss, power_photo, tmp, r, m_div_h, w_div_w_Kep, deltam, csound_face
1302 240 : real(dp), pointer :: v(:)
1303 : logical :: v_flag
1304 :
1305 : include 'formats'
1306 :
1307 240 : ierr = 0
1308 240 : is_int_val = .false.
1309 240 : int_val = 0
1310 240 : val = 0
1311 :
1312 240 : v_flag = .true.
1313 240 : v(1:nz) => s% v(1:nz) ! need this outside of conditional to keep compiler happy
1314 240 : if (s% u_flag) then
1315 0 : v(1:nz) => s% u(1:nz)
1316 240 : else if (.not. s% v_flag) then
1317 240 : v_flag = .false.
1318 : end if
1319 :
1320 480 : if (c > eps_neu_rate_offset) then
1321 0 : i = c - eps_neu_rate_offset
1322 : val = 0
1323 0 : do k = 1, s% nz
1324 0 : val = val + s% eps_neu_rate(i, k) * s% dm(k)
1325 : end do
1326 240 : else if (c > eps_nuc_rate_offset) then
1327 0 : i = c - eps_nuc_rate_offset
1328 : val = 0
1329 0 : do k = 1, s% nz
1330 0 : val = val + s% eps_nuc_rate(i, k) * s% dm(k)
1331 : end do
1332 240 : else if (c > screened_rate_offset) then
1333 0 : i = c - screened_rate_offset
1334 : val = 0
1335 0 : do k = 1, s% nz
1336 0 : val = val + s% screened_rate(i, k) * s% dm(k)
1337 : end do
1338 240 : else if (c > raw_rate_offset) then
1339 0 : i = c - raw_rate_offset
1340 : ! i is the reaction id
1341 : val = 0
1342 0 : do k = 1, s% nz
1343 0 : val = val + s% raw_rate(i, k) * s% dm(k)
1344 : end do
1345 240 : else if (c > log_lum_band_offset) then
1346 : ! We want log Teff, Log g, M/H, Lum/lsun at the photosphere
1347 0 : k = s% photosphere_cell_k
1348 0 : m_div_h = chem_M_div_h(s% X(k), s% Z(k), s% job% initial_zfracs)
1349 0 : if (k > 0) then
1350 0 : i = c - log_lum_band_offset
1351 : val = get_lum_band_by_id(i, safe_log10(s% Teff), &
1352 0 : s% photosphere_logg, m_div_h, s% photosphere_L, ierr)
1353 0 : if (ierr /= 0) return
1354 : end if
1355 0 : val = safe_log10(val * lsun)
1356 240 : else if (c > lum_band_offset) then
1357 : ! We want log Teff, Log g, M/H, Lum/lsun at the photosphere
1358 0 : k = s% photosphere_cell_k
1359 0 : m_div_h = chem_M_div_h(s% X(k), s% Z(k), s% job% initial_zfracs)
1360 0 : if (k > 0) then
1361 0 : i = c - lum_band_offset
1362 : !val = get_lum_band_by_id(i,safe_log10(s% T(k)),safe_log10(s% grav(k)),m_div_h,s% L(k)/lsun, ierr)
1363 : val = get_lum_band_by_id(i, safe_log10(s% Teff), &
1364 0 : s% photosphere_logg, m_div_h, s% photosphere_L, ierr)
1365 0 : val = val * lsun
1366 0 : if (ierr /= 0) return
1367 : end if
1368 240 : else if (c > abs_mag_offset) then
1369 : ! We want log Teff, Log g, M/H, Lum/lsun at the photosphere
1370 0 : k = s% photosphere_cell_k
1371 0 : m_div_h = chem_M_div_h(s% X(k), s% Z(k), s% job% initial_zfracs)
1372 0 : if (k > 0) then
1373 0 : i = c - abs_mag_offset
1374 : !val = get_abs_mag_by_id(i,safe_log10(s% T(k)),safe_log10(s% grav(k)),m_div_h,s% L(k)/lsun, ierr)
1375 : val = get_abs_mag_by_id(i, safe_log10(s% Teff), &
1376 0 : s% photosphere_logg, m_div_h, s% photosphere_L, ierr)
1377 0 : if (ierr /= 0) return
1378 : end if
1379 240 : else if (c > bc_offset) then
1380 : ! We want log Teff, Log g, M/H at the photosphere
1381 0 : k = s% photosphere_cell_k
1382 0 : m_div_h = chem_M_div_h(s% X(k), s% Z(k), s% job% initial_zfracs)
1383 0 : if (k > 0) then
1384 0 : i = c - bc_offset
1385 : !val = get_bc_by_id(i,safe_log10(s% T(k)),safe_log10(s% grav(k)),m_div_h, ierr)
1386 : val = get_bc_by_id(i, safe_log10(s% Teff), &
1387 0 : s% photosphere_logg, m_div_h, ierr)
1388 0 : if (ierr /= 0) return
1389 : end if
1390 240 : else if (c > c_log_eps_burn_offset) then
1391 0 : i = c - c_log_eps_burn_offset
1392 0 : val = safe_log10(abs(s% center_eps_burn(i))) ! abs is for photo
1393 240 : else if (c > max_eps_nuc_offset) then
1394 0 : i = c - max_eps_nuc_offset
1395 0 : val = safe_log10(max_eps_nuc_log_x(s% net_iso(i)))
1396 240 : else if (c > cz_top_max_offset) then
1397 0 : i = c - cz_top_max_offset
1398 0 : val = safe_log10(cz_top_max_log_x(s% net_iso(i)))
1399 240 : else if (c > cz_max_offset) then
1400 0 : i = c - cz_max_offset
1401 0 : val = safe_log10(cz_max_log_x(s% net_iso(i)))
1402 240 : else if (c > log_surface_xa_offset) then
1403 0 : i = c - log_surface_xa_offset
1404 0 : val = safe_log10(surface_avg_x(s, s% net_iso(i)))
1405 240 : else if (c > log_center_xa_offset) then
1406 0 : i = c - log_center_xa_offset
1407 0 : val = safe_log10(center_avg_x(s, s% net_iso(i)))
1408 240 : else if (c > log_average_xa_offset) then
1409 0 : i = c - log_average_xa_offset
1410 0 : val = safe_log10(star_avg_x(s, s% net_iso(i)))
1411 240 : else if (c > log_total_mass_offset) then
1412 0 : i = c - log_total_mass_offset
1413 0 : val = safe_log10(star_avg_x(s, s% net_iso(i)) * s% xmstar / Msun)
1414 240 : else if (c > total_mass_offset) then
1415 8 : i = c - total_mass_offset
1416 8 : val = star_avg_x(s, s% net_iso(i)) * s% xmstar / Msun
1417 232 : else if (c > category_offset) then
1418 12 : i = c - category_offset
1419 12 : val = category_L(i)
1420 220 : else if (c > average_xa_offset) then
1421 0 : i = c - average_xa_offset
1422 0 : val = star_avg_x(s, s% net_iso(i))
1423 220 : elseif (c > surface_xa_offset) then
1424 8 : i = c - surface_xa_offset
1425 8 : val = surface_avg_x(s, s% net_iso(i))
1426 212 : else if (c > center_xa_offset) then
1427 16 : i = c - center_xa_offset
1428 16 : val = center_avg_x(s, s% net_iso(i))
1429 : else
1430 :
1431 4 : select case(c)
1432 :
1433 : case(h_model_number)
1434 4 : is_int_val = .true.
1435 4 : int_val = s% model_number
1436 :
1437 : case(h_log_star_age)
1438 0 : val = safe_log10(s% star_age)
1439 : case(h_star_age)
1440 4 : val = s% star_age
1441 : case(h_log_star_age_sec)
1442 0 : val = safe_log10(s% star_age * secyer)
1443 : case(h_star_age_sec)
1444 0 : val = s% star_age * secyer
1445 : case(h_star_age_min)
1446 0 : val = s% star_age * secyer / 60
1447 : case(h_star_age_hr)
1448 0 : val = s% star_age * secyer / 60 / 60
1449 : case(h_star_age_day)
1450 0 : val = s% star_age * secyer / 60 / 60 / 24
1451 : case(h_star_age_yr)
1452 0 : val = s% star_age
1453 : case(h_day)
1454 0 : val = s% star_age * secyer / 60 / 60 / 24
1455 :
1456 : case(h_time_step)
1457 0 : val = s% time_step
1458 : case(h_log_dt)
1459 4 : val = safe_log10(s% time_step)
1460 : case(h_time_step_sec)
1461 0 : val = s% dt
1462 : case(h_log_dt_sec)
1463 0 : val = safe_log10(s% time_step * secyer)
1464 : case(h_time_step_days)
1465 0 : val = s% time_step * secyer / 60 / 60 / 24
1466 : case(h_log_dt_days)
1467 0 : val = safe_log10(s% time_step * secyer / 60 / 60 / 24)
1468 :
1469 : case(h_log_star_mass)
1470 0 : val = safe_log10(s% star_mass)
1471 : case(h_star_mass)
1472 4 : val = s% star_mass
1473 : case(h_log_xmstar)
1474 4 : val = safe_log10(s% xmstar)
1475 : case(h_delta_mass)
1476 0 : val = s% star_mass - s% initial_mass
1477 : case(h_star_mdot)
1478 0 : val = s% star_mdot
1479 : case(h_log_abs_mdot)
1480 4 : val = safe_log10(abs(s% star_mdot))
1481 :
1482 : case(h_m_center)
1483 0 : val = s% M_center / Msun
1484 : case(h_r_center)
1485 0 : val = s% R_center / Rsun
1486 : case(h_m_center_gm)
1487 0 : val = s% M_center
1488 : case(h_r_center_cm)
1489 0 : val = s% R_center
1490 : case(h_r_center_km)
1491 0 : val = s% R_center * 1d-5
1492 : case(h_L_center)
1493 0 : val = s% L_center / Lsun
1494 : case(h_log_L_center_ergs_s)
1495 0 : val = safe_log10(s% L_center)
1496 : case(h_log_L_center)
1497 0 : val = safe_log10(s% L_center / Lsun)
1498 : case(h_v_center)
1499 0 : val = s% v_center
1500 : case(h_v_center_kms)
1501 0 : val = s% v_center * 1d-5
1502 : case(h_infall_div_cs)
1503 0 : if (s% v_center < 0d0) val = -s% v_center / s% csound(s% nz)
1504 :
1505 : case(h_mdot_timescale)
1506 0 : val = s% star_mass / max(1d-99, abs(s% star_mdot))
1507 :
1508 : case(h_kh_div_mdot_timescales)
1509 : val = s% kh_timescale / &
1510 0 : (s% star_mass / max(1d-99, abs(s% star_mdot)))
1511 : case(h_dlnR_dlnM)
1512 0 : if (abs(s% star_mdot) > 1d-99) &
1513 : val = (s% lnR(1) - s% lnR_start(1)) / &
1514 0 : ((s% mstar - s% mstar_old) / (0.5d0 * (s% mstar + s% mstar_old)))
1515 : case(h_star_gravitational_mass)
1516 0 : val = s% m_grav(1) / Msun
1517 : case(h_star_mass_grav_div_mass)
1518 0 : val = s% m_grav(1) / s% m(1)
1519 :
1520 : case(h_e_thermal)
1521 0 : val = sum(s% dm(1:nz) * s% T(1:nz) * s% cp(1:nz))
1522 : case(h_total_angular_momentum)
1523 0 : val = s% total_angular_momentum
1524 : case(h_log_total_angular_momentum)
1525 0 : val = safe_log10(s% total_angular_momentum)
1526 : case(h_species)
1527 0 : int_val = s% species
1528 0 : is_int_val = .true.
1529 : case(h_tau_factor)
1530 0 : val = s% tau_factor
1531 : case(h_tau_surface)
1532 0 : val = s% tau_factor * s% tau_base
1533 : case(h_log_tau_center)
1534 0 : val = safe_log10(s% tau(s% nz))
1535 :
1536 : case(h_logT_max)
1537 0 : val = s% log_max_temperature
1538 : case(h_gamma1_min)
1539 0 : val = s% min_gamma1
1540 :
1541 : case(h_logQ_max)
1542 0 : val = maxval(s% lnd(1:nz) / ln10 - 2 * s% lnT(1:nz) / ln10 + 12)
1543 : case(h_logQ_min)
1544 0 : val = minval(s% lnd(1:nz) / ln10 - 2 * s% lnT(1:nz) / ln10 + 12)
1545 :
1546 : case(h_num_zones)
1547 4 : int_val = nz
1548 4 : is_int_val = .true.
1549 : case(h_num_retries)
1550 4 : int_val = s% num_retries
1551 4 : is_int_val = .true.
1552 :
1553 : case(h_avg_skipped_setvars_per_step)
1554 0 : val = dble(s% num_skipped_setvars) / max(1, s% model_number)
1555 : case(h_avg_setvars_per_step)
1556 0 : val = dble(s% num_setvars) / max(1, s% model_number)
1557 : case(h_avg_solver_setvars_per_step)
1558 0 : val = dble(s% num_solver_setvars) / max(1, s% model_number)
1559 :
1560 : case(h_total_num_solver_iterations)
1561 0 : int_val = s% total_num_solver_iterations
1562 0 : is_int_val = .true.
1563 : case(h_total_num_solver_calls_made)
1564 0 : int_val = s% total_num_solver_calls_made
1565 0 : is_int_val = .true.
1566 : case(h_total_num_solver_calls_converged)
1567 0 : int_val = s% total_num_solver_calls_converged
1568 0 : is_int_val = .true.
1569 : case(h_total_num_solver_calls_failed)
1570 : int_val = s% total_num_solver_calls_made - &
1571 0 : s% total_num_solver_calls_converged
1572 0 : is_int_val = .true.
1573 :
1574 : case(h_total_num_solver_relax_iterations)
1575 0 : int_val = s% total_num_solver_relax_iterations
1576 0 : is_int_val = .true.
1577 : case(h_total_num_solver_relax_calls_made)
1578 0 : int_val = s% total_num_solver_relax_calls_made
1579 0 : is_int_val = .true.
1580 : case(h_total_num_solver_relax_calls_converged)
1581 0 : int_val = s% total_num_solver_relax_calls_converged
1582 0 : is_int_val = .true.
1583 : case(h_total_num_solver_relax_calls_failed)
1584 : int_val = s% total_num_solver_relax_calls_made - &
1585 0 : s% total_num_solver_relax_calls_converged
1586 0 : is_int_val = .true.
1587 :
1588 : case(h_total_step_attempts)
1589 0 : int_val = s% total_step_attempts
1590 0 : is_int_val = .true.
1591 : case(h_total_step_retries)
1592 0 : int_val = s% total_step_retries
1593 0 : is_int_val = .true.
1594 : case(h_total_step_redos)
1595 0 : int_val = s% total_step_redos
1596 0 : is_int_val = .true.
1597 : case(h_total_steps_taken)
1598 : int_val = s% total_step_attempts - &
1599 0 : s% total_step_retries - s% total_step_redos
1600 0 : is_int_val = .true.
1601 : case(h_total_steps_finished)
1602 0 : int_val = s% total_steps_finished
1603 0 : is_int_val = .true.
1604 :
1605 : case(h_total_relax_step_attempts)
1606 0 : int_val = s% total_relax_step_attempts
1607 0 : is_int_val = .true.
1608 : case(h_total_relax_step_retries)
1609 0 : int_val = s% total_relax_step_retries
1610 0 : is_int_val = .true.
1611 : case(h_total_relax_step_redos)
1612 0 : int_val = s% total_relax_step_redos
1613 0 : is_int_val = .true.
1614 : case(h_total_relax_steps_taken)
1615 : int_val = s% total_relax_step_attempts - &
1616 0 : s% total_relax_step_retries - s% total_relax_step_redos
1617 0 : is_int_val = .true.
1618 : case(h_total_relax_steps_finished)
1619 0 : int_val = s% total_relax_steps_finished
1620 0 : is_int_val = .true.
1621 :
1622 : case(h_avg_num_solver_iters)
1623 : val = dble(s% total_num_solver_iterations) / &
1624 0 : dble(s% total_num_solver_calls_made)
1625 :
1626 : case(h_num_solver_iterations)
1627 0 : int_val = s% num_solver_iterations
1628 0 : is_int_val = .true.
1629 : case(h_num_iters)
1630 4 : int_val = s% num_solver_iterations
1631 4 : is_int_val = .true.
1632 :
1633 : case(h_h1_czb_mass)
1634 0 : val = s% h1_czb_mass
1635 : case(h_surf_c12_minus_o16)
1636 0 : val = s% surface_c12 - s% surface_o16
1637 : case(h_surf_num_c12_div_num_o16)
1638 0 : val = (16d0 / 12d0) * s% surface_c12 / max(1d-99, s% surface_o16)
1639 : case(h_conv_mx1_top)
1640 4 : val = s% conv_mx1_top
1641 : case(h_conv_mx1_bot)
1642 4 : val = s% conv_mx1_bot
1643 : case(h_conv_mx2_top)
1644 4 : val = s% conv_mx2_top
1645 : case(h_conv_mx2_bot)
1646 4 : val = s% conv_mx2_bot
1647 : case(h_mx1_top)
1648 4 : val = s% mx1_top
1649 : case(h_mx1_bot)
1650 4 : val = s% mx1_bot
1651 : case(h_mx2_top)
1652 4 : val = s% mx2_top
1653 : case(h_mx2_bot)
1654 4 : val = s% mx2_bot
1655 : case(h_conv_mx1_top_r)
1656 0 : val = s% conv_mx1_top_r
1657 : case(h_conv_mx1_bot_r)
1658 0 : val = s% conv_mx1_bot_r
1659 : case(h_conv_mx2_top_r)
1660 0 : val = s% conv_mx2_top_r
1661 : case(h_conv_mx2_bot_r)
1662 0 : val = s% conv_mx2_bot_r
1663 : case(h_mx1_top_r)
1664 0 : val = s% mx1_top_r
1665 : case(h_mx1_bot_r)
1666 0 : val = s% mx1_bot_r
1667 : case(h_mx2_top_r)
1668 0 : val = s% mx2_top_r
1669 : case(h_mx2_bot_r)
1670 0 : val = s% mx2_bot_r
1671 : case(h_epsnuc_M_1)
1672 4 : val = epsnuc_out(1)
1673 : case(h_epsnuc_M_2)
1674 4 : val = epsnuc_out(2)
1675 : case(h_epsnuc_M_3)
1676 4 : val = epsnuc_out(3)
1677 : case(h_epsnuc_M_4)
1678 4 : val = epsnuc_out(4)
1679 : case(h_epsnuc_M_5)
1680 4 : val = epsnuc_out(5)
1681 : case(h_epsnuc_M_6)
1682 4 : val = epsnuc_out(6)
1683 : case(h_epsnuc_M_7)
1684 4 : val = epsnuc_out(7)
1685 : case(h_epsnuc_M_8)
1686 4 : val = epsnuc_out(8)
1687 :
1688 : case(h_power_h_burn)
1689 0 : val = s% power_h_burn
1690 : case(h_log_LH)
1691 4 : val = safe_log10(s% power_h_burn)
1692 :
1693 : case(h_power_he_burn)
1694 0 : val = s% power_he_burn
1695 : case(h_log_LHe)
1696 4 : val = safe_log10(s% power_he_burn)
1697 :
1698 : case(h_power_photo)
1699 0 : val = s% power_photo
1700 : case(h_Lnuc_photo)
1701 0 : val = safe_log10(abs(s% power_photo))
1702 :
1703 : case(h_power_z_burn)
1704 0 : val = s% power_z_burn
1705 : case(h_log_LZ)
1706 4 : val = safe_log10(s% power_z_burn)
1707 :
1708 : case(h_log_Lneu)
1709 0 : val = safe_log10(s% power_neutrinos)
1710 : case(h_log_Lneu_nuc)
1711 0 : val = safe_log10(s% power_nuc_neutrinos)
1712 : case(h_log_Lneu_nonnuc)
1713 0 : val = safe_log10(s% power_nonnuc_neutrinos)
1714 :
1715 : case(h_Lsurf_m)
1716 0 : val = s% m(1) / Msun
1717 : case(h_luminosity)
1718 0 : val = s% L_surf
1719 : case(h_log_L)
1720 4 : val = safe_log10(s% L_surf)
1721 : case(h_luminosity_ergs_s)
1722 0 : val = s% L_surf * Lsun
1723 : case(h_log_L_ergs_s)
1724 0 : val = safe_log10(s% L_surf * Lsun)
1725 :
1726 : case(h_photosphere_cell_log_density)
1727 0 : val = s% lnd(s% photosphere_cell_k) / ln10
1728 : case(h_photosphere_cell_density)
1729 0 : val = s% rho(s% photosphere_cell_k)
1730 : case(h_photosphere_cell_log_opacity)
1731 0 : val = safe_log10(s% opacity(s% photosphere_cell_k))
1732 : case(h_photosphere_cell_opacity)
1733 0 : val = s% opacity(s% photosphere_cell_k)
1734 : case(h_photosphere_cell_log_free_e)
1735 0 : val = s% lnfree_e(s% photosphere_cell_k) / ln10
1736 : case(h_photosphere_cell_free_e)
1737 0 : val = exp(s% lnfree_e(s% photosphere_cell_k))
1738 :
1739 : case(h_log_Teff)
1740 4 : val = safe_log10(s% Teff)
1741 : case(h_Teff)
1742 0 : val = s% Teff
1743 : case(h_effective_T)
1744 0 : val = s% Teff
1745 :
1746 : case(h_photosphere_cell_k)
1747 0 : int_val = s% photosphere_cell_k
1748 0 : is_int_val = .true.
1749 : case(h_photosphere_cell_log_T)
1750 0 : val = s% lnT(s% photosphere_cell_k) / ln10
1751 : case(h_photosphere_cell_T)
1752 0 : val = s% T(s% photosphere_cell_k)
1753 : case(h_photosphere_T)
1754 0 : val = s% photosphere_T
1755 : case(h_photosphere_black_body_T)
1756 0 : val = s% photosphere_black_body_T
1757 : case(h_photosphere_logg)
1758 0 : val = s% photosphere_logg
1759 : case(h_photosphere_m)
1760 0 : val = s% photosphere_m
1761 : case(h_photosphere_xm)
1762 0 : val = s% star_mass - s% photosphere_m
1763 : case(h_photosphere_L)
1764 0 : val = s% photosphere_L
1765 : case(h_photosphere_r)
1766 0 : val = s% photosphere_r
1767 : case(h_photosphere_log_L)
1768 0 : val = safe_log10(s% photosphere_L)
1769 : case(h_photosphere_log_r)
1770 0 : val = safe_log10(s% photosphere_r)
1771 : case(h_photosphere_csound)
1772 0 : val = s% photosphere_csound
1773 : case(h_photosphere_opacity)
1774 0 : val = s% photosphere_opacity
1775 : case(h_photosphere_column_density)
1776 0 : val = s% photosphere_column_density
1777 : case(h_photosphere_log_column_density)
1778 0 : val = safe_log10(s% photosphere_column_density)
1779 : case(h_photosphere_v_km_s)
1780 0 : val = s% photosphere_v / 1d5
1781 : case(h_v_phot_km_s)
1782 0 : val = s% photosphere_v / 1d5
1783 : case(h_photosphere_v_div_cs)
1784 0 : val = s% photosphere_v / s% photosphere_csound
1785 :
1786 : case(h_one_div_yphot)
1787 0 : val = 1d0 / s% photosphere_column_density
1788 : case(h_log_one_div_yphot)
1789 0 : val = safe_log10(1d0 / s% photosphere_column_density)
1790 : case(h_min_opacity)
1791 0 : val = minval(s% opacity(1:s% nz))
1792 : case(h_log_min_opacity)
1793 0 : val = safe_log10(minval(s% opacity(1:s% nz)))
1794 :
1795 : case(h_radius_cm)
1796 0 : val = s% R(1)
1797 : case(h_log_R_cm)
1798 0 : val = safe_log10(s% R(1))
1799 : case(h_radius)
1800 0 : val = s% R(1) / Rsun
1801 : case(h_log_R)
1802 4 : val = safe_log10(s% R(1) / Rsun)
1803 :
1804 : case(h_gravity)
1805 0 : val = s% grav(1)
1806 : case(h_log_g)
1807 4 : val = safe_log10(s% grav(1))
1808 :
1809 : case(h_log_cntr_dr_cm)
1810 0 : val = safe_log10(s% r(s% nz) - s% R_center)
1811 : case(h_log_max_T)
1812 0 : val = s% log_max_temperature
1813 : case(h_log_cntr_T)
1814 4 : val = s% log_center_temperature
1815 : case(h_log_cntr_Rho)
1816 4 : val = s% log_center_density
1817 : case(h_log_cntr_P)
1818 4 : val = s% log_center_pressure
1819 : case(h_log_center_T)
1820 0 : val = s% log_center_temperature
1821 : case(h_log_center_Rho)
1822 0 : val = s% log_center_density
1823 : case(h_log_center_P)
1824 0 : val = s% log_center_pressure
1825 :
1826 : case(h_max_T)
1827 0 : val = exp10(s% log_max_temperature)
1828 : case(h_center_T)
1829 0 : val = exp10(s% log_center_temperature)
1830 : case(h_center_Rho)
1831 0 : val = exp10(s% log_center_density)
1832 : case(h_center_P)
1833 0 : val = exp10(s% log_center_pressure)
1834 :
1835 : case(h_log_mesh_adjust_IE_conservation)
1836 0 : val = safe_log10(s% mesh_adjust_IE_conservation)
1837 : case(h_log_mesh_adjust_PE_conservation)
1838 0 : val = safe_log10(s% mesh_adjust_PE_conservation)
1839 : case(h_log_mesh_adjust_KE_conservation)
1840 0 : val = safe_log10(s% mesh_adjust_KE_conservation)
1841 :
1842 : case(h_total_IE_div_IE_plus_KE)
1843 : val = s% total_internal_energy_end / &
1844 0 : (s% total_internal_energy_end + s% total_radial_kinetic_energy_end)
1845 : case(h_total_entropy)
1846 0 : val = dot_product(s% dm(1:nz), s% entropy(1:nz))
1847 :
1848 : case(h_total_internal_energy_after_adjust_mass)
1849 0 : val = s% total_internal_energy_after_adjust_mass
1850 : case(h_total_gravitational_energy_after_adjust_mass)
1851 0 : val = s% total_gravitational_energy_after_adjust_mass
1852 : case(h_total_radial_kinetic_energy_after_adjust_mass)
1853 0 : val = s% total_radial_kinetic_energy_after_adjust_mass
1854 : case(h_total_rotational_kinetic_energy_after_adjust_mass)
1855 0 : val = s% total_rotational_kinetic_energy_after_adjust_mass
1856 : case(h_total_turbulent_energy_after_adjust_mass)
1857 0 : val = s% total_turbulent_energy_after_adjust_mass
1858 : case(h_total_energy_after_adjust_mass)
1859 0 : val = s% total_energy_after_adjust_mass
1860 :
1861 : case(h_total_internal_energy)
1862 0 : val = s% total_internal_energy_end
1863 : case(h_total_gravitational_energy)
1864 0 : val = s% total_gravitational_energy_end
1865 : case(h_total_radial_kinetic_energy)
1866 0 : val = s% total_radial_kinetic_energy_end
1867 : case(h_total_rotational_kinetic_energy)
1868 0 : val = s% total_rotational_kinetic_energy_end
1869 : case(h_total_turbulent_energy)
1870 0 : val = s% total_turbulent_energy_end
1871 : case(h_total_energy)
1872 0 : val = s% total_energy_end
1873 : case(h_total_energy_foe)
1874 0 : val = s% total_energy_end * 1d-51
1875 :
1876 : case(h_log_total_internal_energy)
1877 0 : val = safe_log10(s% total_internal_energy_end)
1878 : case(h_log_total_gravitational_energy)
1879 0 : val = safe_log10(abs(s% total_gravitational_energy_end))
1880 : case(h_log_total_radial_kinetic_energy)
1881 0 : val = safe_log10(s% total_radial_kinetic_energy_end)
1882 : case(h_log_total_rotational_kinetic_energy)
1883 0 : val = safe_log10(s% total_rotational_kinetic_energy_end)
1884 : case(h_log_total_turbulent_energy)
1885 0 : val = safe_log10(s% total_turbulent_energy_end)
1886 : case(h_log_total_energy)
1887 0 : val = safe_log10(abs(s% total_energy_end))
1888 :
1889 : case(h_avg_abs_v_div_cs)
1890 0 : if (v_flag) &
1891 0 : val = sum(abs(v(1:nz)) / s% csound(1:nz)) / nz
1892 : case(h_log_avg_abs_v_div_cs)
1893 0 : if (v_flag) &
1894 0 : val = safe_log10(sum(abs(v(1:nz)) / s% csound(1:nz)) / nz)
1895 : case(h_max_abs_v_div_cs)
1896 0 : if (v_flag) &
1897 0 : val = maxval(abs(v(1:nz)) / s% csound(1:nz))
1898 : case(h_log_max_abs_v_div_cs)
1899 0 : if (v_flag) &
1900 0 : val = safe_log10(maxval(abs(v(1:nz)) / s% csound(1:nz)))
1901 :
1902 : case(h_avg_abs_v)
1903 0 : if (v_flag) &
1904 0 : val = sum(abs(v(1:nz))) / nz
1905 : case(h_log_avg_abs_v)
1906 0 : if (v_flag) &
1907 0 : val = safe_log10(sum(abs(v(1:nz))) / nz)
1908 : case(h_max_abs_v)
1909 0 : if (v_flag) &
1910 0 : val = maxval(abs(v(1:nz)))
1911 : case(h_log_max_abs_v)
1912 0 : if (v_flag) &
1913 0 : val = safe_log10(maxval(abs(v(1:nz))))
1914 :
1915 : case(h_virial_thm_P_avg)
1916 0 : val = s% virial_thm_P_avg
1917 : case(h_virial_thm_rel_err)
1918 0 : val = s% virial_thm_P_avg
1919 0 : val = (val + s% total_gravitational_energy_end) / val
1920 :
1921 : case(h_total_eps_grav)
1922 0 : val = s% total_eps_grav
1923 : case(h_work_outward_at_surface)
1924 0 : val = s% work_outward_at_surface
1925 : case(h_work_inward_at_center)
1926 0 : val = s% work_inward_at_center
1927 : case(h_total_nuclear_heating)
1928 0 : val = s% total_nuclear_heating
1929 : case(h_total_non_nuc_neu_cooling)
1930 0 : val = s% total_non_nuc_neu_cooling
1931 : case(h_total_irradiation_heating)
1932 0 : val = s% total_irradiation_heating
1933 : case(h_total_WD_sedimentation_heating)
1934 0 : if (s% do_element_diffusion) val = s% total_WD_sedimentation_heating
1935 : case(h_total_extra_heating)
1936 0 : val = s% total_extra_heating
1937 :
1938 : case(h_total_energy_sources_and_sinks)
1939 0 : val = s% total_energy_sources_and_sinks
1940 :
1941 : case(h_error_in_energy_conservation)
1942 0 : val = s% error_in_energy_conservation
1943 : case(h_rel_error_in_energy_conservation)
1944 0 : val = s% error_in_energy_conservation / abs(s% total_energy_end)
1945 : case(h_log_rel_error_in_energy_conservation)
1946 0 : val = safe_log10(abs(s% error_in_energy_conservation / s% total_energy_end))
1947 :
1948 : case(h_tot_E_equ_err)
1949 0 : val = sum(s% ergs_error(1:nz))
1950 : case(h_tot_E_err)
1951 0 : val = s% error_in_energy_conservation
1952 : case(h_rel_E_err)
1953 0 : if (s% total_energy_end /= 0d0) &
1954 0 : val = s% error_in_energy_conservation / abs(s% total_energy_end)
1955 : case(h_abs_rel_E_err)
1956 0 : if (s% total_energy_end /= 0d0) &
1957 0 : val = abs(s% error_in_energy_conservation / s% total_energy_end)
1958 : case(h_log_rel_E_err)
1959 0 : if (s% total_energy_end /= 0d0) &
1960 0 : val = safe_log10(abs(s% error_in_energy_conservation / s% total_energy_end))
1961 :
1962 : case(h_cumulative_energy_error)
1963 0 : val = s% cumulative_energy_error
1964 : case(h_rel_cumulative_energy_error)
1965 0 : if (s% total_energy_end /= 0d0) &
1966 0 : val = s% cumulative_energy_error / abs(s% total_energy_end)
1967 : case(h_log_rel_cumulative_energy_error)
1968 0 : if (s% total_energy_end /= 0d0) &
1969 0 : val = safe_log10(abs(s% cumulative_energy_error / s% total_energy_end))
1970 : case(h_rel_run_E_err)
1971 0 : if (s% total_energy_end /= 0d0) &
1972 0 : val = s% cumulative_energy_error / s% total_energy_end
1973 : case(h_log_rel_run_E_err)
1974 0 : if (s% total_energy_end /= 0d0) &
1975 0 : val = safe_log10(abs(s% cumulative_energy_error / s% total_energy_end))
1976 :
1977 : case(h_u_surf_km_s)
1978 0 : if (s% u_flag) val = s% u_face_ad(1)%val * 1d-5
1979 : case(h_u_surf)
1980 0 : if (s% u_flag) val = s% u_face_ad(1)%val
1981 : case(h_u_div_csound_max)
1982 0 : if (s% u_flag) val = maxval(abs(s% u(1:nz)) / s% csound(1:nz))
1983 : case(h_u_div_csound_surf)
1984 0 : if (s% u_flag) val = s% u_face_ad(1)%val / s% csound_face(1)
1985 :
1986 : case(h_surf_escape_v)
1987 0 : val = sqrt(2 * s% cgrav(1) * s% m(1) / (s% r(1)))
1988 : case(h_v_surf_div_escape_v)
1989 0 : val = v_surf / sqrt(2 * s% cgrav(1) * s% m(1) / (s% r(1)))
1990 : case(h_v_div_vesc)
1991 0 : val = v_surf / sqrt(2 * s% cgrav(1) * s% m(1) / (s% r(1)))
1992 : case(h_v_surf_km_s)
1993 0 : val = v_surf * 1d-5
1994 : case(h_v_surf)
1995 0 : val = v_surf
1996 : case(h_v_surf_div_v_kh)
1997 0 : val = v_surf / (s% photosphere_r / s% kh_timescale)
1998 : case(h_v_div_csound_max)
1999 0 : if (v_flag) val = maxval(abs(v(1:nz)) / s% csound_face(1:nz))
2000 : case(h_v_div_csound_surf)
2001 4 : val = v_surf / csound_surf
2002 : case(h_v_div_cs)
2003 0 : if (s% u_flag) then
2004 0 : val = s% u(1) / s% csound(1)
2005 0 : else if (s% v_flag) then
2006 0 : val = s% v(1) / s% csound(1)
2007 : else
2008 : val = 0d0
2009 : end if
2010 : case(h_remnant_M)
2011 0 : val = get_remnant_mass(s) / Msun
2012 : case(h_ejecta_M)
2013 0 : val = get_ejecta_mass(s) / Msun
2014 :
2015 : case(h_log_L_div_Ledd)
2016 0 : Ledd = eval_Ledd(s, ierr)
2017 0 : if (ierr /= 0 .or. Ledd == 0d0) then
2018 0 : ierr = 0
2019 : else
2020 0 : val = safe_log10(s% L_surf * Lsun / Ledd)
2021 : end if
2022 :
2023 : case(h_lum_div_Ledd)
2024 0 : Ledd = eval_Ledd(s, ierr)
2025 0 : if (ierr /= 0 .or. Ledd == 0d0) then
2026 0 : ierr = 0
2027 : else
2028 0 : val = s% L_surf * Lsun / Ledd
2029 : end if
2030 :
2031 : case(h_gradT_excess_alpha)
2032 0 : val = s% gradT_excess_alpha
2033 : case(h_gradT_excess_min_beta)
2034 0 : val = s% gradT_excess_min_beta
2035 : case(h_gradT_excess_max_lambda)
2036 0 : val = s% gradT_excess_max_lambda
2037 :
2038 : case(h_max_L_rad_div_Ledd)
2039 0 : do k = 1, nz
2040 0 : tmp = get_Lrad_div_Ledd(s, k)
2041 0 : if (tmp > val) val = tmp
2042 : end do
2043 :
2044 : case(h_max_L_rad_div_Ledd_div_phi_Joss)
2045 0 : do k = 1, nz
2046 0 : tmp = get_Lrad_div_Ledd(s, k)
2047 0 : phi_Joss = get_phi_Joss(s, k)
2048 0 : if (tmp / phi_Joss > val) val = tmp / phi_Joss
2049 : end do
2050 :
2051 : case(h_i_rot_total)
2052 0 : if(s% rotation_flag) then
2053 : val = 0d0
2054 0 : do k = 1, s% nz
2055 0 : val = val + s% dm_bar(k) * s%i_rot(k)% val
2056 : end do
2057 : end if
2058 : case(h_surf_avg_j_rot)
2059 0 : val = if_rot(s% j_rot_avg_surf)
2060 : case(h_surf_avg_omega)
2061 0 : val = if_rot(s% omega_avg_surf)
2062 : case(h_surf_avg_omega_crit)
2063 0 : val = if_rot(s% omega_crit_avg_surf)
2064 : case(h_surf_avg_omega_div_omega_crit)
2065 0 : val = if_rot(s% w_div_w_crit_avg_surf)
2066 :
2067 : case(h_surf_avg_v_rot)
2068 0 : val = if_rot(s% v_rot_avg_surf) * 1d-5 ! km/sec
2069 : case(h_surf_avg_v_crit)
2070 0 : val = if_rot(s% v_crit_avg_surf) * 1d-5 ! km/sec
2071 : case(h_surf_avg_v_div_v_crit)
2072 0 : val = if_rot(s% v_div_v_crit_avg_surf)
2073 :
2074 : case(h_surf_avg_Lrad_div_Ledd)
2075 0 : val = s% Lrad_div_Ledd_avg_surf
2076 : case(h_surf_avg_opacity)
2077 0 : val = s% opacity_avg_surf
2078 : case(h_surf_avg_logT)
2079 0 : val = s% logT_avg_surf
2080 : case(h_surf_avg_logRho)
2081 0 : val = s% logRho_avg_surf
2082 :
2083 : case(h_v_wind_Km_per_s)
2084 : val = 1d-5 * s% opacity(1) * max(0d0, -s% mstar_dot) / &
2085 0 : (pi4 * s% photosphere_r * Rsun * s% tau_base)
2086 :
2087 : case (h_kh_mdot_limit)
2088 0 : if(s% rotation_flag) then
2089 0 : val = s% rotational_mdot_kh_fac * s% star_mass / s% kh_timescale
2090 : else
2091 : val = 0d0
2092 : end if
2093 : case (h_log_rotational_mdot_boost)
2094 0 : val = safe_log10(if_rot(s% rotational_mdot_boost))
2095 : case (h_rotational_mdot_boost)
2096 0 : val = if_rot(s% rotational_mdot_boost)
2097 :
2098 : case(h_min_Pgas_div_P)
2099 0 : val = minval(s% Pgas(1:nz) / s% Peos(1:nz))
2100 :
2101 : case(h_center_degeneracy)
2102 0 : val = s% center_degeneracy
2103 : case(h_log_center_eps_nuc)
2104 0 : val = safe_log10(s% center_eps_nuc)
2105 : case(h_center_eps_nuc)
2106 0 : val = s% center_eps_nuc
2107 : case(h_d_center_eps_nuc_dlnT)
2108 0 : val = s% d_center_eps_nuc_dlnT
2109 : case(h_d_center_eps_nuc_dlnd)
2110 0 : val = s% d_center_eps_nuc_dlnd
2111 : case(h_center_eps_grav)
2112 0 : val = center_value(s, s% eps_grav_ad(1:nz)% val)
2113 : case(h_center_non_nuc_neu)
2114 0 : val = s% center_non_nuc_neu
2115 : case(h_center_gamma)
2116 0 : val = center_value(s, s% gam)
2117 : case(h_center_zbar)
2118 0 : val = s% center_zbar
2119 : case(h_center_abar)
2120 4 : val = s% center_abar
2121 : case(h_center_mu)
2122 4 : val = s% center_mu
2123 : case(h_center_ye)
2124 4 : val = s% center_ye
2125 : case(h_center_entropy)
2126 0 : val = s% center_entropy
2127 : case(h_max_entropy)
2128 0 : val = s% max_entropy
2129 : case(h_compactness)
2130 0 : if (s% m(1) > 2.5d0 * Msun) then
2131 0 : do k = nz - 1, 1, -1
2132 0 : if (s% m(k) > 2.5d0 * Msun) exit
2133 : end do
2134 0 : r = s% r(k + 1) + (s% r(k) - s% r(k + 1)) * (2.5d0 * Msun - s% m(k + 1)) / s% dm(k)
2135 0 : val = 2.5d0 / (r / 1d8)
2136 : end if
2137 : case(h_compactness_parameter)
2138 0 : if (s% m(1) > 2.5d0 * Msun) then
2139 0 : do k = nz - 1, 1, -1
2140 0 : if (s% m(k) > 2.5d0 * Msun) exit
2141 : end do
2142 0 : r = s% r(k + 1) + (s% r(k) - s% r(k + 1)) * (2.5d0 * Msun - s% m(k + 1)) / s% dm(k)
2143 0 : val = 2.5d0 / (r / 1d8)
2144 : end if
2145 : case(h_mu4)
2146 0 : deltam = 0.3d0 * msun ! Ertl et al 2016
2147 0 : if (s% entropy(1) > 4.0d0) then
2148 0 : do k = nz - 1, 1, -1
2149 0 : if (s% entropy(k) > 4.d0) exit
2150 : end do
2151 0 : do k2 = nz - 1, 1, -1
2152 0 : if (s% m(k2) > s%m(k) + deltam) exit
2153 : end do
2154 :
2155 0 : val = (deltam / msun) / ((s% r(k2) - s% r(k)) / 1d8)
2156 : end if
2157 : case(h_m4)
2158 0 : if (s% entropy(1) > 4.0d0) then
2159 0 : do k = nz - 1, 1, -1
2160 0 : if (s% entropy(k) > 4.d0) exit
2161 : end do
2162 0 : val = s%m(k) / msun
2163 : end if
2164 : case(h_max_infall_speed)
2165 0 : if (s% u_flag) then
2166 0 : val = -minval(s% u(1:s% nz)) * 1d-5 ! convert to km/sec
2167 0 : else if (s% v_flag) then
2168 0 : val = -minval(s% v(1:s% nz)) * 1d-5 ! convert to km/sec
2169 : end if
2170 : case(h_fe_core_infall)
2171 0 : val = s% fe_core_infall * 1d-5 ! convert to km/sec
2172 : case(h_non_fe_core_infall)
2173 0 : val = s% non_fe_core_infall * 1d-5 ! convert to km/sec
2174 : case(h_non_fe_core_rebound)
2175 0 : val = s% non_fe_core_rebound * 1d-5 ! convert to km/sec
2176 : case(h_center_omega)
2177 0 : val = if_rot(s% center_omega)
2178 : case(h_center_omega_div_omega_crit)
2179 0 : val = if_rot(s% center_omega_div_omega_crit)
2180 :
2181 : case(h_surf_r_equatorial_div_r_polar)
2182 0 : if(s%rotation_flag) then
2183 0 : val = s% r_equatorial(1) / s% r_polar(1)
2184 : else
2185 0 : val = 1.0d0
2186 : end if
2187 : case(h_surf_r_equatorial_div_r)
2188 0 : if(s%rotation_flag) then
2189 0 : val = s% r_equatorial(1) / s% r(1)
2190 : else
2191 0 : val = 1.0d0
2192 : end if
2193 : case(h_surf_r_polar_div_r)
2194 0 : if(s%rotation_flag) then
2195 0 : val = s% r_polar(1) / s% r(1)
2196 : else
2197 0 : val = 1.0d0
2198 : end if
2199 : case(h_h_rich_layer_mass)
2200 0 : val = s% star_mass - s% he_core_mass
2201 : case(h_he_rich_layer_mass)
2202 0 : val = max(0d0, s% he_core_mass - s% co_core_mass)
2203 : case(h_co_rich_layer_mass)
2204 0 : val = max(0d0, s% co_core_mass - s% he_core_mass)
2205 :
2206 : case(h_he_core_mass)
2207 4 : val = s% he_core_mass
2208 : case(h_he_core_radius)
2209 0 : val = s% he_core_radius
2210 : case(h_he_core_lgT)
2211 0 : val = s% he_core_lgT
2212 : case(h_he_core_lgRho)
2213 0 : val = s% he_core_lgRho
2214 : case(h_he_core_L)
2215 0 : val = s% he_core_L
2216 : case(h_he_core_v)
2217 0 : val = s% he_core_v
2218 : case(h_he_core_omega)
2219 0 : val = if_rot(s% he_core_omega)
2220 : case(h_he_core_omega_div_omega_crit)
2221 0 : val = if_rot(s% he_core_omega_div_omega_crit)
2222 : case(h_he_core_k)
2223 0 : int_val = s% he_core_k
2224 0 : is_int_val = .true.
2225 :
2226 : case(h_co_core_mass)
2227 4 : val = s% co_core_mass
2228 : case(h_co_core_radius)
2229 0 : val = s% co_core_radius
2230 : case(h_co_core_lgT)
2231 0 : val = s% co_core_lgT
2232 : case(h_co_core_lgRho)
2233 0 : val = s% co_core_lgRho
2234 : case(h_co_core_L)
2235 0 : val = s% co_core_L
2236 : case(h_co_core_v)
2237 0 : val = s% co_core_v
2238 : case(h_co_core_omega)
2239 0 : val = if_rot(s% co_core_omega)
2240 : case(h_co_core_omega_div_omega_crit)
2241 0 : val = if_rot(s% co_core_omega_div_omega_crit)
2242 : case(h_co_core_k)
2243 0 : int_val = s% co_core_k
2244 0 : is_int_val = .true.
2245 :
2246 : case(h_one_core_mass)
2247 4 : val = s% one_core_mass
2248 : case(h_one_core_radius)
2249 0 : val = s% one_core_radius
2250 : case(h_one_core_lgT)
2251 0 : val = s% one_core_lgT
2252 : case(h_one_core_lgRho)
2253 0 : val = s% one_core_lgRho
2254 : case(h_one_core_L)
2255 0 : val = s% one_core_L
2256 : case(h_one_core_v)
2257 0 : val = s% one_core_v
2258 : case(h_one_core_omega)
2259 0 : val = if_rot(s% one_core_omega)
2260 : case(h_one_core_omega_div_omega_crit)
2261 0 : val = if_rot(s% one_core_omega_div_omega_crit)
2262 : case(h_one_core_k)
2263 0 : int_val = s% one_core_k
2264 0 : is_int_val = .true.
2265 :
2266 : case(h_fe_core_mass)
2267 4 : val = s% fe_core_mass
2268 : case(h_fe_core_radius)
2269 0 : val = s% fe_core_radius
2270 : case(h_fe_core_lgT)
2271 0 : val = s% fe_core_lgT
2272 : case(h_fe_core_lgRho)
2273 0 : val = s% fe_core_lgRho
2274 : case(h_fe_core_L)
2275 0 : val = s% fe_core_L
2276 : case(h_fe_core_v)
2277 0 : val = s% fe_core_v
2278 : case(h_fe_core_omega)
2279 0 : val = if_rot(s% fe_core_omega)
2280 : case(h_fe_core_omega_div_omega_crit)
2281 0 : val = if_rot(s% fe_core_omega_div_omega_crit)
2282 : case(h_fe_core_k)
2283 0 : int_val = s% fe_core_k
2284 0 : is_int_val = .true.
2285 :
2286 : case(h_neutron_rich_core_mass)
2287 4 : val = s% neutron_rich_core_mass
2288 : case(h_neutron_rich_core_radius)
2289 0 : val = s% neutron_rich_core_radius
2290 : case(h_neutron_rich_core_lgT)
2291 0 : val = s% neutron_rich_core_lgT
2292 : case(h_neutron_rich_core_lgRho)
2293 0 : val = s% neutron_rich_core_lgRho
2294 : case(h_neutron_rich_core_L)
2295 0 : val = s% neutron_rich_core_L
2296 : case(h_neutron_rich_core_v)
2297 0 : val = s% neutron_rich_core_v
2298 : case(h_neutron_rich_core_omega)
2299 0 : val = if_rot(s% neutron_rich_core_omega)
2300 : case(h_neutron_rich_core_omega_div_omega_crit)
2301 0 : val = if_rot(s% neutron_rich_core_omega_div_omega_crit)
2302 : case(h_neutron_rich_core_k)
2303 0 : int_val = s% neutron_rich_core_k
2304 0 : is_int_val = .true.
2305 :
2306 : case(h_envelope_mass)
2307 0 : val = s% star_mass - s% he_core_mass
2308 : case(h_envelope_fraction_left)
2309 0 : val = envelope_fraction_left
2310 : case(h_dynamic_timescale)
2311 0 : val = s% dynamic_timescale
2312 : case(h_kh_timescale)
2313 0 : val = s% kh_timescale
2314 : case(h_nuc_timescale)
2315 0 : val = s% nuc_timescale
2316 : case(h_dt_div_max_tau_conv)
2317 0 : if(s% max_conv_time_scale > 0d0) val = s% dt / s% max_conv_time_scale
2318 : case(h_dt_div_min_tau_conv)
2319 0 : if(s% min_conv_time_scale > 0d0) val = s% dt / s% min_conv_time_scale
2320 : case(h_max_tau_conv)
2321 0 : val = s% max_conv_time_scale
2322 : case(h_min_tau_conv)
2323 0 : val = s% min_conv_time_scale
2324 : case(h_log_max_tau_conv)
2325 0 : val = safe_log10(s% max_conv_time_scale)
2326 : case(h_log_min_tau_conv)
2327 0 : val = safe_log10(s% min_conv_time_scale)
2328 : case(h_tau_QHSE_yrs)
2329 0 : if (s% v_flag .or. s% u_flag) then
2330 0 : do k = 1, nz
2331 0 : if (s% q(k) > s% max_q_for_QHSE_timescale) cycle
2332 0 : if (s% q(k) < s% min_q_for_QHSE_timescale) exit
2333 0 : val = max(val, QHSE_time_scale(s,k))
2334 : end do
2335 : end if
2336 0 : val = val / secyer
2337 : case(h_eps_grav_integral)
2338 0 : val = dot_product(s% dm(1:nz), s% eps_grav_ad(1:nz)% val) / Lsun
2339 : case(h_extra_L)
2340 0 : val = dot_product(s% dm(1:nz), s% extra_heat(1:nz)%val) / Lsun
2341 : case(h_log_extra_L)
2342 0 : val = safe_log10(dot_product(s% dm(1:nz), s% extra_heat(1:nz)%val) / Lsun)
2343 : case(h_log_abs_Lgrav)
2344 0 : val = safe_log10(abs(dot_product(s% dm(1:nz), s% eps_grav_ad(1:nz)%val) / Lsun))
2345 : case(h_log_Lnuc)
2346 5444 : power_photo = dot_product(s% dm(1:nz), s% eps_nuc_categories(iphoto, 1:nz)) / Lsun
2347 4 : val = safe_log10(s% power_nuc_burn - power_photo)
2348 : case(h_Lnuc)
2349 5444 : power_photo = dot_product(s% dm(1:nz), s% eps_nuc_categories(iphoto, 1:nz)) / Lsun
2350 4 : val = s% power_nuc_burn - power_photo
2351 : case(h_log_Lnuc_ergs_s)
2352 0 : power_photo = dot_product(s% dm(1:nz), s% eps_nuc_categories(iphoto, 1:nz))
2353 0 : val = safe_log10(s% power_nuc_burn * Lsun - power_photo)
2354 : case(h_log_power_nuc_burn)
2355 4 : val = safe_log10(s% power_nuc_burn)
2356 : case(h_power_nuc_burn)
2357 4 : val = s% power_nuc_burn
2358 : case(h_log_Lnuc_sub_log_L)
2359 0 : power_photo = dot_product(s% dm(1:nz), s% eps_nuc_categories(iphoto, 1:nz)) / Lsun
2360 0 : val = safe_log10(s% power_nuc_burn - power_photo)
2361 0 : val = val - safe_log10(s% L_surf * Lsun / Lsun)
2362 : case(h_mass_loc_of_max_eps_nuc)
2363 0 : k = maxloc(s% eps_nuc(1:nz), dim = 1)
2364 0 : val = (s% m(k) - s% dm(k) / 2) / Msun
2365 : case(h_mass_ext_to_max_eps_nuc)
2366 0 : k = maxloc(s% eps_nuc(1:nz), dim = 1)
2367 0 : val = (1d0 - s% q(k) + 0.5d0 * s% dq(k)) * s% xmstar / Msun
2368 :
2369 : case(h_diffusion_time_H_He_bdy)
2370 0 : if (s% he_core_k > 0) then
2371 : val = (s% tau(s% he_core_k) - s% tau_factor * s% tau_base) * &
2372 0 : s% r(s% he_core_k) / clight
2373 : end if
2374 : case(h_temperature_H_He_bdy)
2375 0 : if (s% he_core_k > 0) val = s% T(s% he_core_k)
2376 :
2377 : case(h_total_ni_co_56)
2378 0 : if (s% net_iso(ico56) > 0 .and. s% net_iso(ini56) > 0) &
2379 : val = dot_product(s% dm(1:nz), &
2380 : s% xa(s% net_iso(ico56), 1:nz) + &
2381 0 : s% xa(s% net_iso(ini56), 1:nz)) / Msun
2382 :
2383 : case(h_shock_velocity)
2384 0 : if (s% shock_k > 0) val = s% shock_velocity
2385 : case(h_shock_csound)
2386 0 : if (s% shock_k > 0) val = s% shock_csound
2387 : case(h_shock_v_div_cs)
2388 0 : if (s% shock_csound > 0) &
2389 0 : val = s% shock_velocity / s% shock_csound
2390 : case(h_shock_lgT)
2391 0 : if (s% shock_k > 0) val = s% shock_lgT
2392 : case(h_shock_lgRho)
2393 0 : if (s% shock_k > 0) val = s% shock_lgRho
2394 : case(h_shock_lgP)
2395 0 : if (s% shock_k > 0) val = s% shock_lgP
2396 : case(h_shock_q)
2397 0 : if (s% shock_k > 0) val = s% shock_q
2398 : case(h_shock_tau)
2399 0 : if (s% shock_k > 0) val = s% shock_tau
2400 : case(h_shock_mass)
2401 0 : if (s% shock_k > 0) val = s% shock_mass
2402 : case(h_shock_mass_gm)
2403 0 : if (s% shock_k > 0) val = s% shock_mass * Msun
2404 : case(h_shock_radius)
2405 0 : if (s% shock_k > 0) val = s% shock_radius
2406 : case(h_shock_radius_cm)
2407 0 : if (s% shock_k > 0) val = s% shock_radius * Rsun
2408 : case(h_shock_gamma1)
2409 0 : if (s% shock_k > 0) val = s% shock_gamma1
2410 : case(h_shock_entropy)
2411 0 : if (s% shock_k > 0) val = s% shock_entropy
2412 : case(h_shock_pre_lgRho)
2413 0 : if (s% shock_k > 0) val = s% shock_pre_lgRho
2414 : case(h_shock_k)
2415 0 : if (s% shock_k > 0) int_val = s% shock_k
2416 0 : is_int_val = .true.
2417 :
2418 : case(h_surface_optical_depth)
2419 0 : val = s% tau_base * s% tau_factor
2420 : case(h_log_surf_optical_depth)
2421 0 : val = safe_log10(s% tau_base * s% tau_factor)
2422 :
2423 : case(h_log_surf_cell_opacity)
2424 0 : val = safe_log10(s% opacity(1))
2425 : case(h_log_surf_cell_density)
2426 0 : val = s% lnd(1) / ln10
2427 : case(h_surface_cell_temperature)
2428 0 : val = s% T(1)
2429 : case(h_log_surf_cell_temperature)
2430 0 : val = s% lnT(1) / ln10
2431 : case(h_surface_cell_entropy)
2432 0 : val = s% entropy(1)
2433 : case(h_log_surf_cell_P)
2434 0 : val = s% lnPeos(1) / ln10
2435 : case(h_log_surf_cell_pressure)
2436 0 : val = s% lnPeos(1) / ln10
2437 : case(h_log_surf_cell_z)
2438 : val = 0
2439 0 : if (s% net_iso(ih1) /= 0) val = val + s% xa(s% net_iso(ih1), 1)
2440 0 : if (s% net_iso(ih2) /= 0) val = val + s% xa(s% net_iso(ih2), 1)
2441 0 : if (s% net_iso(ihe3) /= 0) val = val + s% xa(s% net_iso(ihe3), 1)
2442 0 : if (s% net_iso(ihe4) /= 0) val = val + s% xa(s% net_iso(ihe4), 1)
2443 0 : val = safe_log10(1d0 - val)
2444 :
2445 : case(h_log_Ledd)
2446 0 : Ledd = eval_Ledd(s, ierr)
2447 0 : if (ierr /= 0 .or. Ledd <= 0d0) then
2448 0 : ierr = 0
2449 : else
2450 0 : val = safe_log10(Ledd / Lsun)
2451 : end if
2452 :
2453 : case(h_dt_div_dt_cell_collapse)
2454 0 : val = s% dt / eval_min_cell_collapse_time(s, 2, nz, min_k, ierr)
2455 : case(h_dt_cell_collapse)
2456 0 : val = eval_min_cell_collapse_time(s, 2, nz, min_k, ierr)
2457 :
2458 : case(h_min_dr_div_cs_k)
2459 0 : val = min_dr_div_cs(s, int_val)
2460 0 : val = dble(int_val)
2461 0 : is_int_val = .true.
2462 : case(h_min_dr_div_cs)
2463 0 : val = min_dr_div_cs(s, min_k)
2464 : case(h_log_min_dr_div_cs)
2465 0 : val = safe_log10(min_dr_div_cs(s, min_k))
2466 : case(h_min_dr_div_cs_yr)
2467 0 : val = min_dr_div_cs(s, min_k) / secyer
2468 : case(h_log_min_dr_div_cs_yr)
2469 0 : val = safe_log10(min_dr_div_cs(s, min_k) / secyer)
2470 : case(h_dt_div_min_dr_div_cs)
2471 0 : val = min_dr_div_cs(s, min_k)
2472 0 : if (min_k <= 0) then
2473 0 : val = 1d99
2474 : else
2475 0 : val = s% dt / val
2476 : end if
2477 : case(h_log_dt_div_min_dr_div_cs)
2478 0 : val = min_dr_div_cs(s, min_k)
2479 0 : if (min_k <= 0) then
2480 0 : val = 1d99
2481 : else
2482 0 : val = safe_log10(s% dt / val)
2483 : end if
2484 :
2485 : case(h_cz_bot_mass)
2486 0 : if (s% largest_conv_mixing_region /= 0) then
2487 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2488 0 : if (k == nz) then
2489 0 : val = s% M_center / Msun
2490 : else
2491 0 : val = s% m(k) / Msun
2492 : end if
2493 : end if
2494 : case(h_cz_mass)
2495 0 : if (s% largest_conv_mixing_region /= 0) then
2496 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2497 0 : val = s% m(k) / Msun
2498 : end if
2499 : case(h_cz_log_xmass)
2500 0 : if (s% largest_conv_mixing_region == 0) then
2501 0 : val = -99
2502 : else
2503 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2504 0 : val = safe_log10(s% xmstar * sum(s% dq(1:k - 1)))
2505 : end if
2506 : case(h_cz_log_xmsun)
2507 0 : if (s% largest_conv_mixing_region == 0) then
2508 0 : val = -99
2509 : else
2510 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2511 0 : val = safe_log10(s% xmstar * sum(s% dq(1:k - 1)) / Msun)
2512 : end if
2513 : case(h_cz_xm)
2514 0 : if (s% largest_conv_mixing_region /= 0) then
2515 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2516 0 : val = s% xmstar * sum(s% dq(1:k - 1)) / Msun
2517 : end if
2518 : case(h_cz_logT)
2519 0 : if (s% largest_conv_mixing_region /= 0) then
2520 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2521 0 : val = s% lnT(k) / ln10
2522 : end if
2523 : case(h_cz_logRho)
2524 0 : if (s% largest_conv_mixing_region /= 0) then
2525 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2526 0 : val = s% lnd(k) / ln10
2527 : end if
2528 : case(h_cz_logP)
2529 0 : if (s% largest_conv_mixing_region /= 0) then
2530 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2531 0 : val = s% lnPeos(k) / ln10
2532 : end if
2533 : case(h_cz_log_column_depth)
2534 0 : if (s% largest_conv_mixing_region /= 0) then
2535 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2536 0 : val = safe_log10(s% xmstar * sum(s% dq(1:k - 1)) / (pi4 * s% r(k) * s% r(k)))
2537 : end if
2538 : case(h_cz_log_radial_depth)
2539 0 : if (s% largest_conv_mixing_region /= 0) then
2540 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2541 0 : val = safe_log10(s% r(1) - s% r(k))
2542 : end if
2543 : case(h_cz_luminosity)
2544 0 : if (s% largest_conv_mixing_region /= 0) then
2545 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2546 0 : val = s% L(k) / Lsun
2547 : end if
2548 : case(h_cz_log_tau)
2549 0 : if (s% largest_conv_mixing_region /= 0) then
2550 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2551 0 : val = safe_log10(s% tau(k))
2552 : end if
2553 : case(h_cz_opacity)
2554 0 : if (s% largest_conv_mixing_region /= 0) then
2555 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2556 0 : val = s% opacity(k)
2557 : end if
2558 : case(h_cz_log_eps_nuc)
2559 0 : if (s% largest_conv_mixing_region /= 0) then
2560 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2561 0 : val = safe_log10(s% eps_nuc(k))
2562 : end if
2563 : case(h_cz_t_heat)
2564 0 : if (s% largest_conv_mixing_region /= 0) then
2565 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2566 0 : if (s% eps_nuc(k) <= 0) then
2567 0 : val = 1d99
2568 : else
2569 0 : val = s% Cp(k) * s% T(k) / s% eps_nuc(k)
2570 : end if
2571 : end if
2572 : case(h_cz_eta)
2573 0 : if (s% largest_conv_mixing_region /= 0) then
2574 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2575 0 : val = s% eta(k)
2576 : end if
2577 : case(h_cz_csound)
2578 0 : if (s% largest_conv_mixing_region /= 0) then
2579 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2580 0 : val = s% csound(k)
2581 : end if
2582 : case(h_cz_scale_height)
2583 0 : if (s% largest_conv_mixing_region /= 0) then
2584 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2585 0 : val = s% scale_height(k)
2586 : end if
2587 : case(h_cz_grav)
2588 0 : if (s% largest_conv_mixing_region /= 0) then
2589 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2590 0 : val = s% grav(k)
2591 : end if
2592 : case(h_cz_bot_radius)
2593 0 : if (s% largest_conv_mixing_region /= 0) then
2594 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2595 0 : if (k == nz) then
2596 0 : val = s% R_center / Rsun
2597 : else
2598 0 : val = s% R(k) / Rsun
2599 : end if
2600 : end if
2601 : case(h_cz_zone)
2602 0 : if (s% largest_conv_mixing_region == 0) then
2603 0 : k = 0
2604 : else
2605 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2606 : end if
2607 0 : int_val = k
2608 0 : is_int_val = .true.
2609 : case(h_cz_omega)
2610 0 : if (s% largest_conv_mixing_region /= 0 .and. s% rotation_flag) then
2611 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2612 0 : val = s% omega(k)
2613 : end if
2614 : case(h_cz_omega_div_omega_crit)
2615 0 : if (s% largest_conv_mixing_region /= 0 .and. s% rotation_flag) then
2616 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
2617 0 : val = s% omega(k) / omega_crit(s, k)
2618 : end if
2619 :
2620 : case(h_cz_top_mass)
2621 0 : if (s% largest_conv_mixing_region /= 0) then
2622 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2623 0 : val = s% m(k) / Msun
2624 : end if
2625 : case(h_cz_top_log_xmass)
2626 0 : if (s% largest_conv_mixing_region == 0) then
2627 0 : val = -99
2628 : else
2629 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2630 0 : val = safe_log10(s% xmstar * sum(s% dq(1:k - 1)))
2631 : end if
2632 : case(h_cz_top_log_xmsun)
2633 0 : if (s% largest_conv_mixing_region == 0) then
2634 0 : val = -99
2635 : else
2636 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2637 0 : val = safe_log10(s% xmstar * sum(s% dq(1:k - 1)) / Msun)
2638 : end if
2639 : case(h_cz_top_xm)
2640 0 : if (s% largest_conv_mixing_region /= 0) then
2641 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2642 0 : val = s% xmstar * sum(s% dq(1:k - 1)) / Msun
2643 : end if
2644 : case(h_cz_top_logT)
2645 0 : if (s% largest_conv_mixing_region /= 0) then
2646 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2647 0 : val = s% lnT(k) / ln10
2648 : end if
2649 : case(h_cz_top_logRho)
2650 0 : if (s% largest_conv_mixing_region /= 0) then
2651 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2652 0 : val = s% lnd(k) / ln10
2653 : end if
2654 : case(h_cz_top_logP)
2655 0 : if (s% largest_conv_mixing_region /= 0) then
2656 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2657 0 : val = s% lnPeos(k) / ln10
2658 : end if
2659 : case(h_cz_top_log_column_depth)
2660 0 : if (s% largest_conv_mixing_region /= 0) then
2661 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2662 0 : val = safe_log10(s% xmstar * sum(s% dq(1:k - 1)) / (pi4 * s% r(k) * s% r(k)))
2663 : end if
2664 : case(h_cz_top_log_radial_depth)
2665 0 : if (s% largest_conv_mixing_region /= 0) then
2666 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2667 0 : val = safe_log10(s% r(1) - s% r(k))
2668 : end if
2669 : case(h_cz_top_luminosity)
2670 0 : if (s% largest_conv_mixing_region /= 0) then
2671 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2672 0 : val = s% L(k) / Lsun
2673 : end if
2674 : case(h_cz_top_log_tau)
2675 0 : if (s% largest_conv_mixing_region /= 0) then
2676 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2677 0 : val = safe_log10(s% tau(k))
2678 : end if
2679 : case(h_cz_top_opacity)
2680 0 : if (s% largest_conv_mixing_region /= 0) then
2681 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2682 0 : val = s% opacity(k)
2683 : end if
2684 : case(h_cz_top_log_eps_nuc)
2685 0 : if (s% largest_conv_mixing_region /= 0) then
2686 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2687 0 : val = safe_log10(s% eps_nuc(k))
2688 : end if
2689 : case(h_cz_top_t_heat)
2690 0 : if (s% largest_conv_mixing_region /= 0) then
2691 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2692 0 : if (s% eps_nuc(k) <= 0) then
2693 0 : val = 1d99
2694 : else
2695 0 : val = s% Cp(k) * s% T(k) / s% eps_nuc(k)
2696 : end if
2697 : end if
2698 : case(h_cz_top_eta)
2699 0 : if (s% largest_conv_mixing_region /= 0) then
2700 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2701 0 : val = s% eta(k)
2702 : end if
2703 : case(h_cz_top_csound)
2704 0 : if (s% largest_conv_mixing_region /= 0) then
2705 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2706 0 : val = s% csound(k)
2707 : end if
2708 : case(h_cz_top_scale_height)
2709 0 : if (s% largest_conv_mixing_region /= 0) then
2710 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2711 0 : val = s% scale_height(k)
2712 : end if
2713 : case(h_cz_top_grav)
2714 0 : if (s% largest_conv_mixing_region /= 0) then
2715 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2716 0 : val = s% grav(k)
2717 : end if
2718 : case(h_cz_top_radius)
2719 0 : if (s% largest_conv_mixing_region /= 0) then
2720 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2721 0 : val = s% R(k) / Rsun
2722 : end if
2723 :
2724 : case(h_mass_conv_core)
2725 4 : val = s% mass_conv_core
2726 : case(h_mass_semiconv_core)
2727 0 : val = s% mass_semiconv_core
2728 :
2729 : case(h_cz_top_zone_logdq)
2730 0 : if (s% largest_conv_mixing_region /= 0) then
2731 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2732 0 : val = safe_log10(s% dq(k))
2733 : end if
2734 : case(h_cz_top_zone)
2735 0 : if (s% largest_conv_mixing_region == 0) then
2736 0 : k = 0
2737 : else
2738 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2739 : end if
2740 0 : int_val = k
2741 0 : is_int_val = .true.
2742 : case(h_cz_top_omega)
2743 0 : if (s% largest_conv_mixing_region /= 0 .and. s% rotation_flag) then
2744 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2745 0 : val = s% omega(k)
2746 : end if
2747 : case(h_cz_top_omega_div_omega_crit)
2748 0 : if (s% largest_conv_mixing_region /= 0 .and. s% rotation_flag) then
2749 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
2750 0 : val = s% omega(k) / omega_crit(s, k)
2751 : end if
2752 :
2753 : case(h_max_gradT_div_grada)
2754 : val = 0
2755 0 : do k = 2, nz
2756 0 : if (s% grada_face(k) == 0) cycle
2757 0 : if (s% gradT(k) / s% grada_face(k) > val) &
2758 0 : val = s% gradT(k) / s% grada_face(k)
2759 : end do
2760 : case(h_max_gradT_sub_grada)
2761 : val = 0
2762 0 : do k = 2, nz
2763 0 : if (s% gradT(k) - s% grada_face(k) > val) &
2764 0 : val = s% gradT(k) - s% grada_face(k)
2765 : end do
2766 : case(h_min_log_mlt_Gamma)
2767 0 : val = 1d99
2768 0 : do k = 2, nz
2769 0 : if (s% mlt_Gamma(k) > 0 .and. s% mlt_Gamma(k) < val) val = s% mlt_Gamma(k)
2770 : end do
2771 0 : val = safe_log10(val)
2772 :
2773 : case(h_max_conv_vel_div_csound, h_max_eff_mix_speed_div_csound, h_max_effective_mixing_speed_div_csound)
2774 : val = 0
2775 0 : do k = 2, nz
2776 0 : if (s% q(k) > s% max_conv_vel_div_csound_maxq) cycle
2777 0 : if (s% use_face_reconstruction) then
2778 0 : if (.not. s% reconstructed_face_state_valid(k)) cycle
2779 0 : csound_face = s% reconstructed_csound_face(k)
2780 : else
2781 0 : csound_face = s% csound_face(k)
2782 : end if
2783 0 : if (csound_face > 0d0) then
2784 0 : if (c == h_max_conv_vel_div_csound) then
2785 0 : if (.not. s% RSP_flag) val = max(val, s% mlt_vc(k)/csound_face)
2786 : else
2787 0 : val = max(val, get_effective_mixing_speed(s,k)/csound_face)
2788 : end if
2789 : end if
2790 : end do
2791 :
2792 : case(h_min_t_eddy)
2793 : val = 0d0
2794 0 : if (.not. s% RSP_flag) then
2795 0 : val = 1d99
2796 0 : do k = 2, nz
2797 0 : if (s% conv_vel(k) <= 0) cycle
2798 0 : if (s% scale_height(k) / s% conv_vel(k) < val) &
2799 0 : val = s% scale_height(k) / s% conv_vel(k)
2800 : end do
2801 : end if
2802 :
2803 : case(h_elapsed_time)
2804 0 : val = s% total_elapsed_time
2805 :
2806 : case(h_delta_nu)
2807 0 : if (.not. s% get_delta_nu_from_scaled_solar) then
2808 0 : val = 1d6 / (2 * s% photosphere_acoustic_r) ! microHz
2809 : else
2810 : val = &
2811 : s% delta_nu_sun * sqrt(s% star_mass) * pow3(s% Teff / s% astero_Teff_sun) / &
2812 0 : pow(s% L_phot, 0.75d0)
2813 : end if
2814 : case(h_delta_Pg)
2815 0 : if (s% calculate_Brunt_N2) val = s% delta_Pg
2816 : case(h_log_delta_Pg)
2817 0 : if (s% calculate_Brunt_N2) val = safe_log10(s% delta_Pg)
2818 : case(h_nu_max)
2819 0 : val = s% nu_max
2820 : case(h_nu_max_3_4th_div_delta_nu)
2821 0 : val = pow(s% nu_max, 0.75d0) / (1d6 / (2 * s% photosphere_acoustic_r))
2822 : case(h_acoustic_cutoff)
2823 0 : val = s% acoustic_cutoff
2824 : case(h_acoustic_radius)
2825 0 : val = s% photosphere_acoustic_r
2826 : case(h_gs_per_delta_nu)
2827 0 : if (s% calculate_Brunt_N2 .and. s% nu_max > 0 .and. s% delta_Pg >= 0) then
2828 0 : val = 1d6 / (2 * s% photosphere_acoustic_r) ! delta_nu
2829 0 : val = 1d6 * val / (s% nu_max * s% nu_max * s% delta_Pg)
2830 : end if
2831 : case(h_ng_for_nu_max)
2832 0 : if (s% calculate_Brunt_N2 .and. s% nu_max > 0 .and. s% delta_Pg >= 0) then
2833 0 : val = 1d6 / (s% nu_max * s% delta_Pg)
2834 : end if
2835 :
2836 : case(h_int_k_r_dr_nu_max_Sl1)
2837 0 : if (s% calculate_Brunt_N2) val = get_int_k_r_dr(s, 1, 1.0d0)
2838 : case(h_int_k_r_dr_2pt0_nu_max_Sl1)
2839 0 : if (s% calculate_Brunt_N2) val = get_int_k_r_dr(s, 1, 2d0)
2840 : case(h_int_k_r_dr_0pt5_nu_max_Sl1)
2841 0 : if (s% calculate_Brunt_N2) val = get_int_k_r_dr(s, 1, 0.5d0)
2842 : case(h_int_k_r_dr_nu_max_Sl2)
2843 0 : if (s% calculate_Brunt_N2) val = get_int_k_r_dr(s, 2, 1.0d0)
2844 : case(h_int_k_r_dr_2pt0_nu_max_Sl2)
2845 0 : if (s% calculate_Brunt_N2) val = get_int_k_r_dr(s, 2, 2d0)
2846 : case(h_int_k_r_dr_0pt5_nu_max_Sl2)
2847 0 : if (s% calculate_Brunt_N2) val = get_int_k_r_dr(s, 2, 0.5d0)
2848 : case(h_int_k_r_dr_nu_max_Sl3)
2849 0 : if (s% calculate_Brunt_N2) val = get_int_k_r_dr(s, 3, 1.0d0)
2850 : case(h_int_k_r_dr_2pt0_nu_max_Sl3)
2851 0 : if (s% calculate_Brunt_N2) val = get_int_k_r_dr(s, 3, 2d0)
2852 : case(h_int_k_r_dr_0pt5_nu_max_Sl3)
2853 0 : if (s% calculate_Brunt_N2) val = get_int_k_r_dr(s, 3, 0.5d0)
2854 :
2855 : case (h_k_below_const_q)
2856 0 : int_val = s% k_below_const_q
2857 0 : is_int_val = .true.
2858 : case (h_q_below_const_q)
2859 0 : if(s% k_below_const_q>0) val = s% q(s% k_below_const_q)
2860 : case (h_logxq_below_const_q)
2861 0 : if(s% k_below_const_q>0) val = safe_log10(sum(s% dq(1:s% k_below_const_q - 1)))
2862 :
2863 : case (h_k_const_mass)
2864 0 : int_val = s% k_const_mass
2865 0 : is_int_val = .true.
2866 : case (h_q_const_mass)
2867 0 : val = s% q(s% k_const_mass)
2868 : case (h_logxq_const_mass)
2869 0 : val = safe_log10(sum(s% dq(1:s% k_const_mass - 1)))
2870 :
2871 : case (h_k_below_just_added)
2872 0 : int_val = s% k_below_just_added
2873 0 : is_int_val = .true.
2874 : case (h_q_below_just_added)
2875 0 : val = s% q(s% k_below_just_added)
2876 : case (h_logxq_below_just_added)
2877 0 : val = safe_log10(sum(s% dq(1:s% k_below_just_added - 1)))
2878 :
2879 : case (h_k_for_test_CpT_absMdot_div_L)
2880 0 : int_val = s% k_for_test_CpT_absMdot_div_L
2881 0 : is_int_val = .true.
2882 : case (h_q_for_test_CpT_absMdot_div_L)
2883 0 : if (s% k_for_test_CpT_absMdot_div_L == nz) then
2884 : val = 0d0
2885 : else
2886 0 : val = s% q(s% k_for_test_CpT_absMdot_div_L)
2887 : end if
2888 : case (h_logxq_for_test_CpT_absMdot_div_L)
2889 0 : if (s% k_for_test_CpT_absMdot_div_L == nz) then
2890 : val = 0d0
2891 : else
2892 0 : val = safe_log10(sum(s% dq(1:s% k_for_test_CpT_absMdot_div_L - 1)))
2893 : end if
2894 :
2895 : case (h_rotation_solver_steps)
2896 0 : int_val = s% num_rotation_solver_steps
2897 0 : is_int_val = .true.
2898 :
2899 : case (h_burn_solver_maxsteps)
2900 0 : if (s% op_split_burn) &
2901 0 : int_val = maxval(s% burn_num_iters(1:s% nz))
2902 0 : is_int_val = .true.
2903 :
2904 : case (h_diffusion_solver_steps)
2905 0 : int_val = s% num_diffusion_solver_steps
2906 0 : is_int_val = .true.
2907 :
2908 : case (h_diffusion_solver_iters)
2909 0 : int_val = s% num_diffusion_solver_iters
2910 0 : is_int_val = .true.
2911 :
2912 : case(h_tot_IE_div_IE_plus_KE)
2913 : val = s% total_internal_energy_end / &
2914 0 : (s% total_internal_energy_end + s% total_radial_kinetic_energy_end)
2915 :
2916 : case(h_tot_E)
2917 0 : val = s% total_energy_end
2918 : case(h_log_tot_E)
2919 0 : val = safe_log10(abs(s% total_energy_end))
2920 :
2921 : case(h_tot_KE)
2922 0 : val = s% total_radial_kinetic_energy_end
2923 : case(h_log_tot_KE)
2924 0 : val = safe_log10(s% total_radial_kinetic_energy_end)
2925 :
2926 : case(h_tot_PE)
2927 0 : val = s% total_gravitational_energy_end
2928 : case(h_log_tot_PE)
2929 0 : val = safe_log10(abs(s% total_gravitational_energy_end))
2930 :
2931 : case(h_tot_IE)
2932 0 : val = s% total_internal_energy_end
2933 : case(h_log_tot_IE)
2934 0 : val = safe_log10(s% total_internal_energy_end)
2935 :
2936 : case(h_tot_Et)
2937 0 : val = s% total_turbulent_energy_end
2938 : case(h_log_tot_Et)
2939 0 : val = safe_log10(s% total_turbulent_energy_end)
2940 :
2941 : case(h_num_hydro_merges)
2942 0 : int_val = s% num_hydro_merges
2943 0 : is_int_val = .true.
2944 : case(h_num_hydro_splits)
2945 0 : int_val = s% num_hydro_splits
2946 0 : is_int_val = .true.
2947 :
2948 : case(h_RSP_DeltaR)
2949 0 : if (s% RSP_flag) val = s% rsp_DeltaR
2950 : case(h_RSP_DeltaMag)
2951 0 : if (s% RSP_flag) val = s% rsp_DeltaMag
2952 : case(h_RSP_GREKM)
2953 0 : if (s% RSP_flag) val = s% rsp_GREKM
2954 :
2955 : case(h_RSP_phase)
2956 0 : if (s% RSP_flag) val = (s% time - rsp_phase_time0()) / s% RSP_period
2957 : case(h_RSP_period_in_days)
2958 0 : if (s% RSP_flag) val = s% RSP_period / secday ! days
2959 : case(h_RSP_num_periods)
2960 0 : if (s% RSP_flag) int_val = s% RSP_num_periods
2961 0 : is_int_val = .true.
2962 :
2963 : case(h_grav_dark_L_polar) ! pole is at inclination = 0
2964 0 : if(s% rotation_flag) then
2965 0 : w_div_w_Kep = if_rot(s% omega(1) * sqrt(pow3(s% r_equatorial(1)) / (s% cgrav(1) * s% m(1))))
2966 0 : val = gravity_darkening_L_coeff(w_div_w_Kep, 0.0d0) * s% L_surf
2967 : else
2968 : val = 0d0
2969 : end if
2970 : case(h_grav_dark_Teff_polar)
2971 0 : if(s% rotation_flag) then
2972 0 : w_div_w_Kep = if_rot(s% omega(1) * sqrt(pow3(s% r_equatorial(1)) / (s% cgrav(1) * s% m(1))))
2973 0 : val = gravity_darkening_Teff_coeff(w_div_w_Kep, 0.0d0) * s% Teff
2974 : else
2975 : val = 0d0
2976 : end if
2977 : case(h_grav_dark_L_equatorial) ! equator is at inclination = pi/2
2978 0 : if(s% rotation_flag) then
2979 0 : w_div_w_Kep = if_rot(s% omega(1) * sqrt(pow3(s% r_equatorial(1)) / (s% cgrav(1) * s% m(1))))
2980 0 : val = gravity_darkening_L_coeff(w_div_w_Kep, 0.5d0 * pi) * s% L_surf
2981 : else
2982 : val = 0d0
2983 : end if
2984 : case(h_grav_dark_Teff_equatorial)
2985 0 : if(s% rotation_flag) then
2986 0 : w_div_w_Kep = if_rot(s% omega(1) * sqrt(pow3(s% r_equatorial(1)) / (s% cgrav(1) * s% m(1))))
2987 0 : val = gravity_darkening_Teff_coeff(w_div_w_Kep, 0.5d0 * pi) * s% Teff
2988 : else
2989 : val = 0d0
2990 : end if
2991 :
2992 : case(h_apsidal_constant_k2)
2993 0 : val = apsidal_constant(s, 2)
2994 :
2995 : ! following items correspond to names on terminal output lines
2996 :
2997 : case(h_lg_Lnuc_tot) ! _tot indicates the inclusion of photodisintegations
2998 0 : val = safe_log10(s% power_nuc_burn)
2999 :
3000 : case(h_H_rich)
3001 0 : val = s% star_mass - max(s% he_core_mass, s% co_core_mass)
3002 :
3003 : case(h_N_cntr)
3004 0 : val = s% center_n14
3005 :
3006 : case(h_lg_Lneu)
3007 0 : val = safe_log10(abs(s% power_neutrinos))
3008 :
3009 : case(h_He_core)
3010 0 : val = s% he_core_mass
3011 :
3012 : case(h_O_cntr)
3013 0 : val = s% center_o16
3014 :
3015 : case(h_lg_Lphoto)
3016 0 : val = safe_log10(abs(s% power_photo))
3017 :
3018 : case(h_CO_core)
3019 0 : val = s% co_core_mass
3020 :
3021 : case(h_Fe_core)
3022 0 : val = s% fe_core_mass
3023 :
3024 : case(h_Ne_cntr)
3025 0 : val = s% center_ne20
3026 :
3027 : case(h_Mass)
3028 0 : val = s% star_mass
3029 :
3030 : case(h_H_cntr)
3031 0 : val = s% center_h1
3032 :
3033 : case(h_Si_cntr)
3034 0 : val = s% center_si28
3035 :
3036 : case(h_lg_Mdot)
3037 0 : val = safe_log10(abs(s% star_mdot))
3038 :
3039 : case(h_He_cntr)
3040 0 : val = s% center_he3 + s% center_he4
3041 :
3042 : case(h_eta_cntr)
3043 0 : val = s% eta(nz)
3044 :
3045 : case(h_gam_cntr)
3046 0 : val = s% gam(nz)
3047 :
3048 : case(h_lg_Dsurf)
3049 0 : val = s% lnd(1) / ln10
3050 :
3051 : case(h_C_cntr)
3052 0 : val = s% center_c12
3053 :
3054 : case(h_phase_of_evolution)
3055 0 : int_val = s% phase_of_evolution
3056 0 : is_int_val = .true.
3057 :
3058 : case(h_zones)
3059 0 : int_val = nz
3060 0 : is_int_val = .true.
3061 :
3062 : case(h_retries)
3063 0 : int_val = s% num_retries
3064 0 : is_int_val = .true.
3065 :
3066 : case (h_TDC_num_cells)
3067 : int_val = 0
3068 0 : do k = 1, nz
3069 0 : if (s% tdc_num_iters(k) > 0) then
3070 0 : int_val = int_val + 1
3071 : end if
3072 : end do
3073 0 : is_int_val = .true.
3074 :
3075 : case default
3076 196 : ierr = -1
3077 :
3078 : end select
3079 :
3080 : end if
3081 :
3082 :
3083 : contains
3084 :
3085 :
3086 0 : real(dp) function max_eps_nuc_log_x(j)
3087 : integer, intent(in) :: j
3088 : integer :: k
3089 0 : max_eps_nuc_log_x = 0
3090 0 : if (j == 0) return
3091 0 : k = maxloc(s% eps_nuc(1:nz), dim = 1)
3092 0 : if (k < 1 .or. k > nz) return
3093 0 : max_eps_nuc_log_x = s% xa(j, k)
3094 0 : end function max_eps_nuc_log_x
3095 :
3096 :
3097 0 : real(dp) function cz_top_max_log_x(j)
3098 : integer, intent(in) :: j
3099 : integer :: k
3100 0 : cz_top_max_log_x = 0
3101 0 : if (s% largest_conv_mixing_region == 0) return
3102 0 : k = s% mixing_region_top(s% largest_conv_mixing_region)
3103 0 : if (j == 0 .or. k <= 1) return
3104 0 : cz_top_max_log_x = s% xa(j, k - 1)
3105 0 : end function cz_top_max_log_x
3106 :
3107 :
3108 0 : real(dp) function cz_max_log_x(j)
3109 : integer, intent(in) :: j
3110 : integer :: k
3111 0 : cz_max_log_x = 0
3112 0 : if (s% largest_conv_mixing_region == 0) return
3113 0 : k = s% mixing_region_bottom(s% largest_conv_mixing_region)
3114 0 : if (j == 0 .or. k <= 1) return
3115 0 : cz_max_log_x = s% xa(j, k - 1)
3116 0 : end function cz_max_log_x
3117 :
3118 :
3119 12 : real(dp) function category_L(i)
3120 : integer, intent(in) :: i
3121 12 : if (i == 0) then
3122 12 : category_L = 0
3123 : return
3124 : end if
3125 : category_L = &
3126 16332 : safe_log10(dot_product(s% eps_nuc_categories(i, 1:nz), s% dm(1:nz)) / Lsun)
3127 12 : end function category_L
3128 :
3129 :
3130 0 : real(dp) function if_rot(v, alt)
3131 : real(dp), intent(in) :: v
3132 : real(dp), optional, intent(in) :: alt
3133 0 : if (s% rotation_flag) then
3134 0 : if_rot = v
3135 : else
3136 : if (present(alt)) then
3137 : if_rot = alt
3138 : else
3139 : if_rot = 0
3140 : end if
3141 : end if
3142 : end function if_rot
3143 :
3144 :
3145 : end subroutine history_getval
3146 :
3147 :
3148 0 : real(dp) function get_int_k_r_dr(s, el, nu_factor)
3149 : use utils_lib, only : is_bad_num
3150 : use chem_def, only : ih1
3151 : type (star_info), pointer :: s
3152 : integer, intent(in) :: el
3153 : real(dp), intent(in) :: nu_factor
3154 :
3155 : real(dp) :: integral, cs2, r2, n2, sl2, omega2, &
3156 : L2, kr2, dr, r0_outer, r0_inner, xh1
3157 : integer :: k, k1, k_inner, k_outer, h1
3158 :
3159 : logical :: dbg
3160 :
3161 : include 'formats'
3162 :
3163 0 : dbg = .false. !(el == 1 .and. nu_factor == 1d0)
3164 :
3165 : get_int_k_r_dr = 0
3166 0 : L2 = el * (el + 1)
3167 0 : omega2 = pow2(1d-6 * 2 * pi * s% nu_max * nu_factor)
3168 :
3169 : ! k_inner and k_outer are bounds of evanescent region
3170 :
3171 : ! k_outer is outermost k where Sl2 <= omega2 at k-1 and Sl2 > omega2 at k
3172 : ! 1st find outermost where Sl2 <= omega2
3173 :
3174 0 : h1 = s% net_iso(ih1)
3175 :
3176 0 : k1 = 0
3177 0 : do k = 1, s% nz
3178 0 : r2 = s% r(k) * s% r(k)
3179 0 : cs2 = s% csound_face(k) * s% csound_face(k)
3180 0 : sl2 = L2 * cs2 / r2
3181 0 : if (sl2 <= omega2) then
3182 : k1 = k; exit
3183 : end if
3184 : end do
3185 0 : if (k1 == 0) return
3186 : ! then find next k where Sl2 >= omega2
3187 0 : k_outer = 0
3188 0 : do k = k1 + 1, s% nz
3189 0 : r2 = s% r(k) * s% r(k)
3190 0 : cs2 = s% csound_face(k) * s% csound_face(k)
3191 0 : sl2 = L2 * cs2 / r2
3192 0 : if (sl2 > omega2) then
3193 0 : k_outer = k; exit
3194 : end if
3195 : end do
3196 0 : if (k_outer == 0) return
3197 :
3198 : ! k_inner is next k where N2 >= omega2 at k+1 and N2 < omega2 at k
3199 0 : k_inner = 0
3200 0 : do k = k_outer + 1, s% nz
3201 0 : if ((s% brunt_N2(k) - s% brunt_N2_composition_term(k)) >= omega2) then
3202 : ! Use the thermal component of the Brunt as starting point
3203 0 : k_inner = k; exit
3204 : end if
3205 : end do
3206 0 : if (k_inner == 0) return
3207 :
3208 0 : integral = 0
3209 0 : do k = k_inner - 1, k_outer, -1
3210 0 : r2 = s% r(k) * s% r(k)
3211 0 : cs2 = s% csound_face(k) * s% csound_face(k)
3212 0 : n2 = s% brunt_N2(k)
3213 0 : sl2 = L2 * cs2 / r2
3214 0 : xh1 = s% xa(h1, k)
3215 0 : dr = s% rmid(k - 1) - s% rmid(k)
3216 0 : kr2 = (1 - n2 / omega2) * (1 - Sl2 / omega2) * omega2 / cs2
3217 0 : if (kr2 < 0 .and. n2 < omega2 .and. omega2 < Sl2) &
3218 0 : integral = integral + sqrt(-kr2) * dr
3219 : end do
3220 :
3221 0 : if (integral == 0) return
3222 :
3223 0 : get_int_k_r_dr = integral
3224 :
3225 : if (dbg) write(*, 3) 'r0 inner outer', &
3226 : k_inner, k_outer, r0_inner / Rsun, r0_outer / Rsun, get_int_k_r_dr
3227 :
3228 0 : if (.not. is_bad_num(get_int_k_r_dr)) return
3229 :
3230 0 : write(*, 2) 'el', el
3231 0 : write(*, 1) 'nu_factor', nu_factor
3232 0 : write(*, 1) 's% nu_max*nu_factor', s% nu_max * nu_factor
3233 0 : write(*, 1) 'log10 nu_max*nu_factor', log10(s% nu_max * nu_factor)
3234 0 : write(*, 1) 'Radius at k_inner', s% r(k_inner)
3235 0 : write(*, 1) 'Radius at k_outer', s% r(k_outer)
3236 :
3237 0 : write(*, 1) 'get_int_k_r_dr', get_int_k_r_dr
3238 0 : write(*, 1) 'integral', integral
3239 0 : write(*, 2) 'k_inner', k_inner
3240 0 : write(*, 2) 'k_outer', k_outer
3241 :
3242 0 : call mesa_error(__FILE__, __LINE__, 'get_int_k_r_dr')
3243 :
3244 : end function get_int_k_r_dr
3245 :
3246 :
3247 0 : real(dp) function apsidal_constant(s, j)
3248 : type (star_info), pointer :: s
3249 : integer, intent(in) :: j
3250 : real(dp) :: y, dy, rho_bar, dr_div_r, fprmid3
3251 : integer :: k
3252 :
3253 0 : y = j - 2 ! value at r = 0
3254 0 : do k = s% nz, 1, -1
3255 0 : fprmid3 = pi4 * pow(s% rmid(k), 3)
3256 0 : rho_bar = 3d0 * s% m(k) / fprmid3
3257 0 : if (k == s% nz) then
3258 0 : dr_div_r = s% r(k) / s% rmid(k) ! r(nz+1) would be zero
3259 : else
3260 0 : dr_div_r = (s% r(k) - s% r(k + 1)) / s% rmid(k)
3261 : end if
3262 0 : dy = dr_div_r * (j * (j + 1) - (6d0 * s% rho(k) / rho_bar) * (y + 1d0) - y * (y - 1d0))
3263 0 : y = y + dy
3264 : end do
3265 :
3266 0 : apsidal_constant = (j + 1d0 - y) / (2d0 * (j + y))
3267 :
3268 0 : end function apsidal_constant
3269 :
3270 :
3271 8 : real(dp) function star_avg_x(s, j)
3272 : type (star_info), pointer :: s
3273 : integer, intent(in) :: j
3274 8 : if (j == 0) then
3275 8 : star_avg_x = 0
3276 : return
3277 : end if
3278 21768 : star_avg_x = dot_product(s% xa(j, 1:s% nz), s% dq(1:s% nz)) / sum(s% dq(1:s% nz))
3279 8 : end function star_avg_x
3280 :
3281 :
3282 0 : subroutine get_history_specs(s, num, names, specs, report)
3283 :
3284 : use utils_lib
3285 : use utils_def
3286 :
3287 : type (star_info), pointer :: s
3288 : integer, intent(in) :: num
3289 : character (len = *), intent(in) :: names(:)
3290 : integer, intent(out) :: specs(:)
3291 : logical, intent(in) :: report
3292 :
3293 : integer :: i, ierr, n, j, iounit, t
3294 : logical :: special_case
3295 : character (len = strlen) :: buffer, string
3296 :
3297 : include 'formats'
3298 0 : ierr = 0
3299 0 : if (num <= 0) return
3300 0 : iounit = -1
3301 0 : specs(1:num) = 0
3302 0 : do i = 1, num
3303 0 : buffer = names(i)
3304 0 : n = len_trim(buffer) + 1
3305 0 : buffer(n:n) = ' '
3306 0 : j = 0
3307 0 : t = token(iounit, n, j, buffer, string)
3308 0 : if (t /= name_token) then
3309 0 : if (len_trim(names(i)) > 0 .and. report) &
3310 0 : write(*, *) 'bad value for name of history item ' // trim(names(i))
3311 0 : specs(i) = -1
3312 0 : ierr = 0
3313 0 : cycle
3314 : end if
3315 : special_case = .false.
3316 : specs(i) = do1_history_spec(&
3317 0 : s, iounit, t, n, j, string, buffer, special_case, report, ierr)
3318 0 : if (ierr /= 0 .or. special_case) then
3319 0 : if (report) write(*, *) 'get_history_specs failed for ' // trim(names(i))
3320 0 : specs(i) = -1
3321 0 : ierr = 0
3322 : end if
3323 : end do
3324 :
3325 : end subroutine get_history_specs
3326 :
3327 :
3328 0 : logical function get1_hist_value(s, name, val)
3329 : ! includes other_history_columns from run_star_extras
3330 : use utils_lib, only : integer_dict_lookup
3331 : type (star_info), pointer :: s
3332 : character (len = *) :: name
3333 : real(dp), intent(out) :: val
3334 : integer :: i, ierr, num_extra_cols, num_binary_cols
3335 : character (len = 80), pointer, dimension(:) :: &
3336 0 : extra_col_names, binary_col_names
3337 : real(dp), pointer, dimension(:) :: &
3338 0 : extra_col_vals, binary_col_vals
3339 : include 'formats'
3340 :
3341 0 : get1_hist_value = .false.
3342 :
3343 0 : call integer_dict_lookup(s% history_names_dict, name, i, ierr)
3344 0 : if (ierr /= 0 .or. i <= 0) return ! didn't find it
3345 0 : if (associated(s% pg% pgstar_hist)) then
3346 0 : if (associated(s% pg% pgstar_hist% vals)) then
3347 0 : if (size(s% pg% pgstar_hist% vals, dim = 1) >= i) then
3348 0 : val = s% pg% pgstar_hist% vals(i)
3349 0 : get1_hist_value = .true.
3350 0 : return
3351 : end if
3352 : end if
3353 : end if
3354 :
3355 : ! try extras 1st
3356 0 : if (associated(s% how_many_extra_history_columns) .and. &
3357 : associated(s% data_for_extra_history_columns)) then
3358 0 : num_extra_cols = s% how_many_extra_history_columns(s% id)
3359 0 : if (num_extra_cols > 0) then
3360 : allocate(&
3361 : extra_col_names(num_extra_cols), &
3362 0 : extra_col_vals(num_extra_cols), stat = ierr)
3363 : call s% data_for_extra_history_columns(&
3364 0 : s% id, num_extra_cols, extra_col_names, extra_col_vals, ierr)
3365 0 : do i = 1, num_extra_cols
3366 0 : if (extra_col_names(i) == name) then
3367 0 : val = extra_col_vals(i)
3368 0 : get1_hist_value = .true.
3369 0 : exit
3370 : end if
3371 : end do
3372 0 : deallocate(extra_col_names, extra_col_vals)
3373 0 : if (get1_hist_value) return
3374 : end if
3375 : end if
3376 :
3377 : ! try binary history
3378 0 : num_binary_cols = s% how_many_binary_history_columns(s% binary_id)
3379 0 : if (num_binary_cols > 0) then
3380 : allocate(&
3381 : binary_col_names(num_binary_cols), &
3382 0 : binary_col_vals(num_binary_cols))
3383 : call s% data_for_binary_history_columns(&
3384 0 : s% binary_id, num_binary_cols, binary_col_names, binary_col_vals, ierr)
3385 0 : if (ierr == 0) then
3386 0 : do i = 1, num_binary_cols
3387 0 : if (binary_col_names(i) == name) then
3388 0 : val = binary_col_vals(i)
3389 0 : get1_hist_value = .true.
3390 0 : exit
3391 : end if
3392 : end do
3393 : end if
3394 0 : deallocate(binary_col_names, binary_col_vals)
3395 0 : if (get1_hist_value) return
3396 : end if
3397 :
3398 0 : end function get1_hist_value
3399 :
3400 :
3401 0 : subroutine get_history_values(s, num, specs, &
3402 0 : is_int_value, int_values, values, failed_to_find_value)
3403 : ! note: this doesn't handle user-defined extra columns
3404 :
3405 : use utils_lib
3406 : use utils_def
3407 :
3408 : type (star_info), pointer :: s
3409 : integer, intent(in) :: num
3410 : integer, intent(in) :: specs(:)
3411 : logical, intent(out) :: is_int_value(:)
3412 : integer, intent(out) :: int_values(:)
3413 : real(dp), intent(inout) :: values(:)
3414 : logical, intent(out) :: failed_to_find_value(:)
3415 :
3416 : integer :: i, c, ierr
3417 : real(dp) :: epsnuc_out(12), v_surf, csound_surf, envelope_fraction_left
3418 :
3419 : include 'formats'
3420 : ierr = 0
3421 0 : if (num <= 0) return
3422 :
3423 0 : epsnuc_out(1:4) = s% burn_zone_mass(1:4, 1)
3424 0 : epsnuc_out(5:8) = s% burn_zone_mass(1:4, 2)
3425 0 : epsnuc_out(9:12) = s% burn_zone_mass(1:4, 3)
3426 0 : csound_surf = eval_csound(s, 1, ierr)
3427 :
3428 0 : if (s% u_flag) then
3429 0 : v_surf = s% u(1)
3430 0 : else if (s% v_flag) then
3431 0 : v_surf = s% v(1)
3432 : else
3433 0 : v_surf = 0d0
3434 : end if
3435 :
3436 0 : if (s% initial_mass > s% he_core_mass) then
3437 : envelope_fraction_left = &
3438 0 : (s% star_mass - s% he_core_mass) / (s% initial_mass - s% he_core_mass)
3439 : else
3440 0 : envelope_fraction_left = 1
3441 : end if
3442 :
3443 0 : do i = 1, num
3444 0 : failed_to_find_value(i) = .false.
3445 0 : c = specs(i)
3446 0 : if (c <= 0) then
3447 0 : failed_to_find_value(i) = .true.
3448 : else
3449 : call history_getval(&
3450 : s, c, values(i), int_values(i), is_int_value(i), &
3451 0 : s% nz, v_surf, csound_surf, envelope_fraction_left, epsnuc_out, ierr)
3452 0 : if (ierr /= 0) then
3453 0 : failed_to_find_value(i) = .true.
3454 : ierr = 0
3455 : end if
3456 : end if
3457 : end do
3458 :
3459 : end subroutine get_history_values
3460 :
3461 :
3462 : subroutine get_iso_val(s, str, val, ierr)
3463 : use chem_lib, only : chem_get_iso_id
3464 : type (star_info), pointer :: s
3465 : character (len = *), intent(in) :: str
3466 : real(dp), intent(out) :: val
3467 : integer, intent(out) :: ierr
3468 : integer :: n, split, id, i
3469 : ierr = 0
3470 : val = 0
3471 : n = len_trim(str)
3472 : split = 0
3473 : do i = 1, n
3474 : if (str(i:i) == ' ') then
3475 : split = i
3476 : exit
3477 : end if
3478 : end do
3479 : if (split <= 1 .or. split >= n) then ! no interior space to split str
3480 : ierr = -1
3481 : return
3482 : end if
3483 : id = chem_get_iso_id(str(split + 1:n))
3484 : if (id <= 0) then ! not a valid iso name
3485 : ierr = -1
3486 : return
3487 : end if
3488 : i = s% net_iso(id)
3489 : select case (str(1:split - 1))
3490 : case('center')
3491 : val = center_avg_x(s, i)
3492 : case('surface')
3493 : val = surface_avg_x(s, i)
3494 : case('average')
3495 : val = star_avg_x(s, i)
3496 : case('total_mass')
3497 : val = star_avg_x(s, i) * s% xmstar / Msun
3498 : case('log_total_mass')
3499 : val = safe_log10(star_avg_x(s, i) * s% xmstar / Msun)
3500 : case('log_average')
3501 : val = safe_log10(star_avg_x(s, i))
3502 : case('log_center')
3503 : val = safe_log10(center_avg_x(s, i))
3504 : case('log_surface')
3505 : val = safe_log10(surface_avg_x(s, i))
3506 : case default
3507 : ierr = -1
3508 : end select
3509 : end subroutine get_iso_val
3510 :
3511 : end module history
|
