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