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