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