Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010 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 run_star_support
21 :
22 : use star_lib
23 : use star_def
24 : use chem_def
25 : use chem_lib
26 : use const_def, only: dp, pi, ln10, amu, mp, secyer, kerg, msun, rsun, lsun, mesa_dir, arg_not_provided
27 : use math_lib
28 : use eos_lib
29 : use kap_def
30 : use net_def
31 : use net_lib
32 : use other_extras
33 :
34 : implicit none
35 :
36 : integer :: id_from_read_star_job = 0
37 :
38 : ! Set MESA_INLIST_RESOLVED to true when you no longer want the routine
39 : ! resolve_inlist_fname to look at the MESA_INLIST environment variable
40 : logical :: MESA_INLIST_RESOLVED = .false.
41 :
42 : private
43 : public :: do_read_star_job, do_read_star_job_and_return_id
44 : public :: run1_star
45 : public :: start_run1_star
46 : public :: do_evolve_one_step
47 : public :: after_evolve_loop
48 : public :: failed
49 : public :: id_from_read_star_job
50 : public :: MESA_INLIST_RESOLVED
51 : public :: do_star_job_controls_before, do_star_job_controls_after
52 :
53 : ! deprecated, but kept around for use by binary
54 : public :: before_evolve_loop, after_step_loop, before_step_loop, do_saves, &
55 : resolve_inlist_fname, terminate_normal_evolve_loop, null_binary_controls
56 :
57 : contains
58 :
59 :
60 3 : subroutine run1_star( &
61 : do_alloc_star, do_free_star, okay_to_restart, &
62 : id, restart, &
63 : extras_controls, &
64 : ierr, &
65 : inlist_fname_arg)
66 :
67 : logical, intent(in) :: do_alloc_star, do_free_star, okay_to_restart
68 : integer, intent(inout) :: id ! input if not do_alloc_star
69 : logical, intent(inout) :: restart ! input if not do_alloc_star
70 : character (len=*) :: inlist_fname_arg
71 : integer, intent(out) :: ierr
72 : optional inlist_fname_arg
73 :
74 : interface
75 :
76 : subroutine extras_controls(id, ierr)
77 : implicit none
78 : integer, intent(in) :: id
79 : integer, intent(out) :: ierr
80 : end subroutine extras_controls
81 :
82 : end interface
83 :
84 : logical :: continue_evolve_loop
85 : type (star_info), pointer :: s
86 : character (len=strlen) :: restart_filename
87 :
88 : logical, parameter :: pgstar_ok = .true.
89 : logical, parameter :: dbg = .false.
90 :
91 : 1 format(a35, 99(1pe26.16))
92 : 2 format(a55, i7, 1pe26.16)
93 : 3 format(a15, 2x, f15.6)
94 : 4 format(a15, 2x, e15.6)
95 :
96 : 11 format(a35, f20.10)
97 :
98 1 : restart_filename = 'restart_photo'
99 : call start_run1_star( &
100 : do_alloc_star, do_free_star, okay_to_restart, &
101 : id, restart, restart_filename, pgstar_ok, dbg, &
102 2 : extras_controls, ierr, inlist_fname_arg)
103 1 : if (failed('do_before_evolve_loop',ierr)) return
104 :
105 1 : call star_ptr(id, s, ierr)
106 1 : if (failed('star_ptr',ierr)) return
107 :
108 : continue_evolve_loop = .true.
109 :
110 : if (dbg) write(*,*) 'start evolve_loop'
111 12 : evolve_loop: do while(continue_evolve_loop) ! evolve one step per loop
112 :
113 11 : continue_evolve_loop = do_evolve_one_step(s, dbg, ierr)
114 12 : if (failed('do_evolve_one_step',ierr)) return
115 :
116 : end do evolve_loop
117 :
118 1 : call after_evolve_loop(s% id, do_free_star, ierr)
119 1 : if (failed('after_evolve_loop',ierr)) return
120 :
121 : end subroutine run1_star
122 :
123 :
124 2 : subroutine start_run1_star( &
125 : do_alloc_star, do_free_star, okay_to_restart, &
126 : id, restart, restart_filename, pgstar_ok, dbg, &
127 : extras_controls, ierr, inlist_fname_arg)
128 :
129 : logical, intent(in) :: do_alloc_star, do_free_star, okay_to_restart
130 : integer, intent(inout) :: id ! input if not do_alloc_star
131 : logical, intent(inout) :: restart ! input if not do_alloc_star
132 : logical, intent(in) :: pgstar_ok, dbg
133 : character (len=*) :: restart_filename, inlist_fname_arg
134 : optional inlist_fname_arg
135 : integer, intent(out) :: ierr
136 :
137 : interface
138 :
139 : subroutine extras_controls(id, ierr)
140 : implicit none
141 : integer, intent(in) :: id
142 : integer, intent(out) :: ierr
143 : end subroutine extras_controls
144 :
145 : end interface
146 :
147 : type (star_info), pointer :: s
148 : character (len=strlen) :: inlist_fname
149 :
150 : include 'formats'
151 :
152 1 : ierr = 0
153 :
154 2 : call resolve_inlist_fname(inlist_fname,inlist_fname_arg)
155 :
156 : ! star is initialized here
157 : call do_before_evolve_loop( &
158 : do_alloc_star, okay_to_restart, restart, pgstar_ok, &
159 : null_binary_controls, extras_controls, &
160 : id_from_read_star_job, inlist_fname, restart_filename, &
161 1 : dbg, 0, id, ierr)
162 1 : if (failed('do_before_evolve_loop',ierr)) return
163 :
164 1 : call star_ptr(id, s, ierr)
165 1 : if (failed('star_ptr',ierr)) return
166 :
167 1 : s% doing_timing = .false.
168 1 : s% job% check_before_step_timing = 0
169 1 : s% job% check_step_loop_timing = 0
170 1 : s% job% check_after_step_timing = 0
171 1 : s% job% time0_initial = 0
172 :
173 : end subroutine start_run1_star
174 :
175 :
176 22 : logical function do_evolve_one_step(s, dbg, ierr) result(continue_evolve_loop)
177 : type (star_info), pointer :: s
178 : logical, intent(in) :: dbg
179 : integer, intent(out) :: ierr
180 :
181 : logical :: first_try
182 : integer :: id
183 : integer :: result, model_number
184 :
185 : include 'formats'
186 :
187 11 : ierr = 0
188 11 : id = s% id
189 11 : continue_evolve_loop = .true.
190 :
191 11 : call before_step_loop(s% id, ierr)
192 11 : if (failed('before_step_loop',ierr)) return
193 :
194 11 : result = s% extras_start_step(id)
195 11 : if (result /= keep_going) then
196 1 : continue_evolve_loop = .false.
197 : return
198 : end if
199 :
200 11 : first_try = .true.
201 :
202 0 : step_loop: do ! may need to repeat this loop
203 :
204 11 : if (stop_is_requested(s)) then
205 0 : continue_evolve_loop = .false.
206 0 : result = terminate
207 0 : exit step_loop
208 : end if
209 :
210 11 : result = star_evolve_step(id, first_try)
211 11 : if (result == keep_going) result = star_check_model(id)
212 11 : if (result == keep_going) result = s% extras_check_model(id)
213 11 : if (result == keep_going) result = star_pick_next_timestep(id)
214 11 : if (result == keep_going) exit step_loop
215 :
216 1 : model_number = get_model_number(id, ierr)
217 1 : if (failed('get_model_number',ierr)) return
218 :
219 1 : if (result == retry .and. s% job% report_retries) then
220 0 : write(*,'(i6,3x,a,/)') model_number, &
221 0 : 'retry reason ' // trim(result_reason_str(s% result_reason))
222 : end if
223 :
224 1 : if (result == redo) then
225 0 : result = star_prepare_to_redo(id)
226 : end if
227 1 : if (result == retry) then
228 0 : result = star_prepare_to_retry(id)
229 : end if
230 1 : if (result == terminate) then
231 : continue_evolve_loop = .false.
232 : exit step_loop
233 : end if
234 10 : first_try = .false.
235 :
236 : end do step_loop
237 :
238 : ! once we get here, the only options are keep_going or terminate.
239 : ! redo or retry must be done inside the step_loop
240 :
241 : call after_step_loop(s% id, s% inlist_fname, &
242 11 : dbg, result, ierr)
243 11 : if (failed('after_step_loop',ierr)) return
244 :
245 11 : if (result /= keep_going) then
246 1 : if (result /= terminate) then
247 0 : write(*,2) 'ERROR in result value in run_star_extras: model', &
248 0 : s% model_number
249 0 : write(*,2) 'extras_finish_step must return keep_going or terminate'
250 0 : write(*,2) 'result', result
251 0 : continue_evolve_loop = .false.
252 0 : return
253 : end if
254 1 : if (s% result_reason == result_reason_normal) then
255 : call terminate_normal_evolve_loop(s% id, &
256 1 : dbg, result, ierr)
257 1 : if (failed('terminate_normal_evolve_loop',ierr)) return
258 : end if
259 1 : continue_evolve_loop = .false.
260 1 : return
261 : end if
262 :
263 10 : call do_saves(id, ierr)
264 10 : if (failed('do_saves',ierr)) return
265 :
266 10 : if (s% doing_timing) then
267 0 : call system_clock(s% job% time1_extra,s% job% clock_rate)
268 : s% job% after_step_timing = s% job% after_step_timing + &
269 0 : dble(s% job% time1_extra - s% job% time0_extra) / s% job% clock_rate
270 0 : s% job% check_time_end = eval_total_times(s% id, ierr)
271 : s% job% check_after_step_timing = s% job% check_after_step_timing + &
272 0 : (s% job% check_time_end - s% job% check_time_start)
273 : end if
274 :
275 : end function do_evolve_one_step
276 :
277 :
278 1 : subroutine null_binary_controls(id, binary_id, ierr)
279 : integer, intent(in) :: id, binary_id
280 : integer, intent(out) :: ierr
281 1 : ierr = 0
282 1 : end subroutine null_binary_controls
283 :
284 :
285 : ! Binary requires to set some controls in here, which is why
286 : ! binary_controls and binary_id are arguments. These do nothing
287 : ! for the case of single star evolution.
288 0 : subroutine before_evolve_loop( &
289 : do_alloc_star, okay_to_restart, restart, &
290 : binary_controls, extras_controls, &
291 : id_from_read_star_job, inlist_fname, restart_filename, &
292 : dbg, binary_id, id, ierr)
293 : logical, intent(in) :: do_alloc_star, okay_to_restart
294 : logical :: restart
295 : interface
296 : subroutine binary_controls(id, binary_id, ierr)
297 : implicit none
298 : integer, intent(in) :: id, binary_id
299 : integer, intent(out) :: ierr
300 : end subroutine binary_controls
301 : subroutine extras_controls(id, ierr)
302 : implicit none
303 : integer, intent(in) :: id
304 : integer, intent(out) :: ierr
305 : end subroutine extras_controls
306 : end interface
307 : integer :: id_from_read_star_job
308 : character (len=*) :: inlist_fname, restart_filename
309 : logical, intent(in) :: dbg
310 : integer, intent(in) :: binary_id
311 : integer, intent(out) :: id, ierr
312 : call do_before_evolve_loop( &
313 : do_alloc_star, okay_to_restart, restart, .true., &
314 : binary_controls, extras_controls, &
315 : id_from_read_star_job, inlist_fname, restart_filename, &
316 0 : dbg, binary_id, id, ierr)
317 0 : end subroutine before_evolve_loop
318 :
319 :
320 14 : subroutine do_before_evolve_loop( &
321 : do_alloc_star, okay_to_restart, restart, pgstar_ok, &
322 : binary_controls, extras_controls, &
323 : id_from_read_star_job, inlist_fname, restart_filename, &
324 : dbg, binary_id, id, ierr)
325 : use utils_lib, only: utils_OMP_SET_NUM_THREADS
326 : logical, intent(in) :: do_alloc_star, okay_to_restart, pgstar_ok
327 : logical :: restart
328 : interface
329 : subroutine binary_controls(id, binary_id, ierr)
330 : implicit none
331 : integer, intent(in) :: id, binary_id
332 : integer, intent(out) :: ierr
333 : end subroutine binary_controls
334 : subroutine extras_controls(id, ierr)
335 : implicit none
336 : integer, intent(in) :: id
337 : integer, intent(out) :: ierr
338 : end subroutine extras_controls
339 : end interface
340 : integer :: id_from_read_star_job
341 : character (len=*) :: inlist_fname, restart_filename
342 : character (len=512) :: temp_fname
343 : logical, intent(in) :: dbg
344 : integer, intent(in) :: binary_id
345 : integer, intent(out) :: id, ierr
346 :
347 : type (star_info), pointer :: s
348 :
349 : include 'formats'
350 :
351 2 : if (do_alloc_star) then
352 1 : if (id_from_read_star_job /= 0) then
353 : ! already allocated by read_star_job
354 1 : id = id_from_read_star_job
355 1 : id_from_read_star_job = 0
356 : else
357 0 : call alloc_star(id, ierr)
358 0 : if (failed('alloc_star',ierr)) return
359 : end if
360 1 : call star_ptr(id, s, ierr)
361 1 : if (failed('star_ptr',ierr)) return
362 : else
363 0 : call star_ptr(id, s, ierr)
364 0 : if (failed('star_ptr',ierr)) return
365 0 : call init_starting_star_data(s, ierr)
366 0 : if (failed('init_starting_star_data',ierr)) return
367 : end if
368 :
369 1 : s% inlist_fname = inlist_fname
370 :
371 1 : if (dbg) write(*,*) 'call starlib_init'
372 1 : call starlib_init(s, ierr) ! okay to do extra calls on this
373 1 : if (failed('star_init',ierr)) return
374 :
375 1 : if (dbg) write(*,*) 'call star_set_kap_and_eos_handles'
376 1 : call star_set_kap_and_eos_handles(id, ierr)
377 1 : if (failed('set_star_kap_and_eos_handles',ierr)) return
378 :
379 1 : if (dbg) write(*,*) 'call star_colors_handles'
380 1 : call star_set_colors_handles(id, ierr)
381 1 : if (failed('set_star_colors_handles',ierr)) return
382 :
383 1 : if (dbg) write(*,*) 'call star_setup'
384 1 : call star_setup(id, inlist_fname, ierr)
385 1 : if (failed('star_setup',ierr)) return
386 :
387 1 : if(dbg) write(*,*) 'call add_fpe_checks'
388 1 : call add_fpe_checks(id, s, ierr)
389 1 : if (failed('add_fpe_checks',ierr)) return
390 :
391 1 : if(dbg) write(*,*) 'call multiply_tolerances'
392 1 : call multiply_tolerances(id, s, ierr)
393 1 : if (failed('multiply_tolerances',ierr)) return
394 :
395 1 : if(dbg) write(*,*) 'call pgstar_env_check'
396 1 : call pgstar_env_check(id, s, ierr)
397 1 : if (failed('pgstar_env_check',ierr)) return
398 :
399 : ! testing module-level (atm/eos/kap/net) partials requires single-threaded execution
400 : if (s% solver_test_atm_partials .or. s% solver_test_eos_partials .or. &
401 1 : s% solver_test_kap_partials .or. s% solver_test_net_partials) then
402 0 : if (s% solver_test_partials_k > 0 .and. s% solver_test_partials_dx_0 > 0) then
403 0 : write(*,*) 'Forcing single-thread mode for testing of module-level partials'
404 0 : call utils_OMP_SET_NUM_THREADS(1)
405 : end if
406 : end if
407 :
408 1 : if (len_trim(s% op_mono_data_path) == 0) &
409 : call get_environment_variable( &
410 1 : "MESA_OP_MONO_DATA_PATH", s% op_mono_data_path)
411 :
412 1 : if (len_trim(s% op_mono_data_cache_filename) == 0) &
413 : call get_environment_variable( &
414 1 : "MESA_OP_MONO_DATA_CACHE_FILENAME", s% op_mono_data_cache_filename)
415 :
416 1 : if (len_trim(s% emesh_data_for_op_mono_path) == 0) &
417 : call get_environment_variable( &
418 1 : "MESA_OP_MONO_MASTER_GRID", s% emesh_data_for_op_mono_path)
419 :
420 1 : s% extras_startup => null_extras_startup
421 1 : s% extras_check_model => null_extras_check_model
422 1 : s% extras_start_step => null_extras_start_step
423 1 : s% extras_finish_step => null_extras_finish_step
424 1 : s% extras_after_evolve => null_extras_after_evolve
425 1 : s% how_many_extra_history_columns => null_how_many_extra_history_columns
426 1 : s% data_for_extra_history_columns => null_data_for_extra_history_columns
427 1 : s% how_many_extra_profile_columns => null_how_many_extra_profile_columns
428 1 : s% data_for_extra_profile_columns => null_data_for_extra_profile_columns
429 :
430 1 : if (dbg) write(*,*) 'call extras_controls'
431 1 : call extras_controls(id, ierr)
432 1 : if (ierr /= 0) return
433 :
434 1 : if (restart_filename /= "restart_photo") then
435 0 : temp_fname = trim(s% photo_directory) // '/' // trim(restart_filename)
436 0 : restart_filename = trim(temp_fname)
437 : end if
438 :
439 1 : if (okay_to_restart) then
440 1 : restart = doing_a_restart(restart_filename)
441 : else
442 0 : restart = .false.
443 : end if
444 :
445 1 : if (s% job% show_log_description_at_start .and. .not. restart) then
446 0 : write(*,'(A)')
447 0 : call show_log_description(id, ierr)
448 0 : if (failed('show_log_description',ierr)) return
449 : end if
450 :
451 1 : if (dbg) write(*,*) 'call binary_controls'
452 1 : call binary_controls(id, binary_id, ierr)
453 1 : if (ierr /= 0) return
454 :
455 1 : if (dbg) write(*,*) 'call do_star_job_controls_before'
456 1 : call do_star_job_controls_before(id, s, restart, ierr)
457 1 : if (ierr /= 0) return
458 :
459 1 : if (dbg) write(*,*) 'call do_load1_star'
460 1 : call do_load1_star(id, s, restart, restart_filename, ierr)
461 1 : if (failed('do_load1_star',ierr)) return
462 :
463 1 : if (dbg) write(*,*) 'call do_star_job_controls_after'
464 1 : call do_star_job_controls_after(id, s, restart, pgstar_ok, ierr)
465 1 : if (failed('do_star_job_controls_after',ierr)) return
466 :
467 1 : write(*,'(A)')
468 1 : write(*,'(A)')
469 :
470 1 : if (.not. restart) then
471 1 : if (dbg) write(*,*) 'call before_evolve'
472 1 : call before_evolve(id, ierr)
473 1 : if (failed('before_evolve',ierr)) return
474 : else
475 0 : call show_terminal_header(id, ierr)
476 0 : if (failed('show_terminal_header',ierr)) return
477 : end if
478 :
479 1 : if (dbg) write(*,*) 'call extras_startup'
480 1 : call s% extras_startup(id, restart, ierr)
481 1 : if (failed('extras_startup',ierr)) return
482 :
483 1 : if (s% job% profile_starting_model .and. .not. restart) then
484 0 : call star_set_vars(id, 0d0, ierr)
485 0 : if (failed('star_set_vars',ierr)) return
486 0 : write(*, '(a, i12)') 'save profile for model number ', s% model_number
487 0 : call save_profile(id,3,ierr)
488 0 : if (failed('save_profile',ierr)) return
489 : end if
490 :
491 1 : if (s% model_number == s% job% save_model_number) then
492 0 : call star_set_vars(id, 0d0, ierr)
493 0 : if (failed('star_set_vars',ierr)) return
494 0 : write(*, '(a, i12)') 'write initial model ', s% model_number
495 0 : call star_write_model(id, 'initial.mod', ierr)
496 0 : if (failed('star_write_model',ierr)) return
497 0 : write(*, *) 'saved to ' // 'initial.mod' ! trim(s% job% save_model_filename)
498 : end if
499 :
500 1 : if (len_trim(s% job% echo_at_start) > 0) then
501 0 : write(*,'(A)')
502 0 : write(*,'(a)') trim(s% job% echo_at_start)
503 0 : write(*,'(A)')
504 : end if
505 :
506 1 : end subroutine do_before_evolve_loop
507 :
508 :
509 11 : subroutine before_step_loop(id, ierr)
510 : integer, intent(in) :: id
511 : type (star_info), pointer :: s
512 : integer, intent(out) :: ierr
513 : integer :: model_number, j
514 : integer :: num_DT, num_FreeEOS
515 :
516 : 1 format(a35, 99(1pe26.16))
517 : 2 format(a35, i7, 1pe26.16)
518 : 3 format(a15, 2x, f15.6)
519 : 4 format(a15, 2x, e15.6)
520 :
521 : 11 format(a35, f20.10)
522 :
523 11 : call star_ptr(id, s, ierr)
524 11 : if (ierr/=0) return
525 :
526 11 : s% result_reason = result_reason_normal
527 :
528 : if (s% job% first_model_for_timing >= 0 .and. &
529 11 : s% model_number >= s% job% first_model_for_timing .and. &
530 : .not. s% doing_timing) then
531 0 : s% doing_timing = .true.
532 0 : write(*,*) 'start timing', s% model_number
533 0 : write(*,'(A)')
534 0 : call system_clock(s% job% time0, s% job% clock_rate)
535 0 : s% job% time0_initial = s% job% time0
536 0 : s% job% step_loop_timing = 0
537 0 : s% job% after_step_timing = 0
538 0 : s% job% before_step_timing = 0
539 : end if
540 :
541 11 : if (s% doing_timing) then
542 0 : call system_clock(s% job% time0_extra,s% job% clock_rate)
543 0 : s% job% check_time_start = eval_total_times(s% id, ierr)
544 : end if
545 :
546 11 : if(s% job% num_steps_for_garbage_collection > 0 .and. s% model_number > 1) then
547 0 : if(mod(s% model_number, s% job% num_steps_for_garbage_collection) == 0)then
548 0 : if (s% job% report_garbage_collection) then
549 : call num_eos_files_loaded( &
550 0 : num_DT, num_FreeEOS)
551 0 : write(*,*) "Start garbage collection model_number", s%model_number,"num eosDT", num_DT, &
552 0 : "num FreeEOS",num_FreeEOS
553 : end if
554 0 : call star_do_garbage_collection(s% id,ierr)
555 0 : if (failed('star_do_garbage_collection',ierr)) return
556 : end if
557 :
558 : ! If reporting, we want to look at the step and the next step (to see the difference)
559 : if(mod(s% model_number-1, s% job% num_steps_for_garbage_collection) == 0 &
560 0 : .and. s% job% report_garbage_collection)then
561 : call num_eos_files_loaded( &
562 0 : num_DT, num_FreeEOS)
563 0 : write(*,*) "End garbage collection model_number ", s%model_number,"num eosDT", num_DT, &
564 0 : "num FreeEOS",num_FreeEOS
565 : end if
566 : end if
567 :
568 11 : if (s% job% enable_adaptive_network) then
569 : call star_adjust_net(s% id, &
570 : s% job% min_x_for_keep, &
571 : s% job% min_x_for_n, &
572 : s% job% min_x_for_add, &
573 : s% job% max_Z_for_add, &
574 : s% job% max_N_for_add, &
575 : s% job% max_A_for_add, &
576 0 : ierr)
577 0 : if (failed('star_adjust_net',ierr)) return
578 : end if
579 :
580 11 : if (s% job% auto_extend_net) then
581 11 : call extend_net(s, ierr)
582 11 : if (failed('extend_net',ierr)) return
583 : end if
584 :
585 11 : if (s% use_other_remove_surface) then
586 0 : call s% other_remove_surface(id, ierr, j)
587 0 : if (failed('other_remove_surface',ierr)) return
588 0 : if (j > 0) then
589 0 : call star_remove_surface_at_cell_k(s% id, j, ierr)
590 : end if
591 0 : call do_remove_surface(id, s, ierr)
592 : else
593 11 : call do_remove_surface(id, s, ierr)
594 11 : if (failed('do_remove_surface',ierr)) return
595 : end if
596 :
597 11 : if (s% job% remove_fallback_at_each_step) then
598 0 : call star_remove_fallback(id,ierr)
599 0 : if (failed('star_remove_fallback',ierr)) return
600 : end if
601 :
602 11 : if (s% job% limit_center_logP_at_each_step > -1d90) then
603 : call star_limit_center_logP( &
604 0 : id, s% job% limit_center_logP_at_each_step, ierr)
605 0 : if (failed('star_limit_center_logP',ierr)) return
606 : end if
607 :
608 11 : if (s% job% remove_center_logRho_limit > -1d90) then
609 : call star_remove_center_by_logRho( &
610 0 : id, s% job% remove_center_logRho_limit, ierr)
611 0 : if (failed('star_remove_center_by_logRho',ierr)) return
612 : end if
613 :
614 : if (s% center_ye <= s% job% center_ye_limit_for_v_flag &
615 11 : .and. (.not. s% v_flag) .and. (.not. s% u_flag)) then
616 0 : write(*,1) 'have reached center ye limit', &
617 0 : s% center_ye, s% job% center_ye_limit_for_v_flag
618 0 : write(*,1) 'set v_flag true'
619 0 : call star_set_v_flag(id, .true., ierr)
620 0 : if (failed('star_set_v_flag',ierr)) return
621 0 : if (ierr /= 0) return
622 : end if
623 :
624 11 : if (s% job% change_RSP2_flag_at_model_number == s% model_number) then
625 0 : write(*,*) 'have reached model number for new_RSP2_flag', &
626 0 : s% model_number, s% job% new_RSP2_flag
627 0 : call star_set_RSP2_flag(id, s% job% new_RSP2_flag, ierr)
628 0 : if (failed('star_set_RSP2_flag',ierr)) return
629 : end if
630 :
631 11 : if (s% job% report_mass_not_fe56) call do_report_mass_not_fe56(s)
632 11 : if (s% job% report_cell_for_xm > 0) call do_report_cell_for_xm(s)
633 :
634 11 : model_number = get_model_number(id, ierr)
635 11 : if (failed('get_model_number',ierr)) return
636 :
637 11 : if (s% star_age < s% job% set_cumulative_energy_error_each_step_if_age_less_than) then
638 0 : if (mod(model_number, s% terminal_interval) == 0) &
639 0 : write(*,1) 'cumulative_energy_error reset to', s% job% new_cumulative_energy_error
640 0 : s% cumulative_energy_error = s% job% new_cumulative_energy_error
641 : end if
642 :
643 11 : if (s% doing_timing) then
644 :
645 0 : call system_clock(s% job% time1_extra, s% job% clock_rate)
646 : s% job% before_step_timing = &
647 : s% job% before_step_timing + &
648 0 : dble(s% job% time1_extra - s% job% time0_extra) / s% job% clock_rate
649 :
650 0 : s% job% check_time_end = eval_total_times(s% id, ierr)
651 : s% job% check_before_step_timing = &
652 : s% job% check_before_step_timing + &
653 0 : (s% job% check_time_end - s% job% check_time_start)
654 :
655 0 : s% job% time0_extra = s% job% time1_extra
656 0 : s% job% check_time_start = s% job% check_time_end
657 :
658 : end if
659 :
660 : end subroutine before_step_loop
661 :
662 :
663 11 : subroutine after_step_loop(id, inlist_fname, dbg, result, ierr)
664 : integer, intent(in) :: id
665 : type (star_info), pointer :: s
666 : character (len=*) :: inlist_fname
667 : logical, intent(in) :: dbg
668 : integer, intent(out) :: ierr
669 : integer, intent(inout) :: result
670 : logical :: will_read_pgstar_inlist
671 :
672 : real(dp) :: tmp
673 :
674 : include 'formats'
675 :
676 11 : call star_ptr(id, s, ierr)
677 11 : if (ierr/=0) return
678 :
679 11 : if (s% doing_timing) then
680 0 : call system_clock(s% job% time1_extra,s% job% clock_rate)
681 : s% job% step_loop_timing = s% job% step_loop_timing + &
682 0 : dble(s% job% time1_extra - s% job% time0_extra) / s% job% clock_rate
683 0 : s% job% check_time_end = eval_total_times(s% id, ierr)
684 : s% job% check_step_loop_timing = s% job% check_step_loop_timing + &
685 0 : (s% job% check_time_end - s% job% check_time_start)
686 0 : s% job% time0_extra = s% job% time1_extra
687 0 : s% job% check_time_start = s% job% check_time_end
688 : end if
689 :
690 11 : if (s% model_number == s% job% set_cumulative_energy_error_at_step) then
691 0 : write(*,1) 'set_cumulative_energy_error', s% job% new_cumulative_energy_error
692 0 : s% cumulative_energy_error = s% job% new_cumulative_energy_error
693 : end if
694 :
695 11 : if (is_bad(s% total_energy_end)) then
696 0 : ierr = 1
697 0 : return
698 : end if
699 :
700 11 : if(s% total_energy_end /= 0d0) then
701 11 : if (abs(s% cumulative_energy_error/s% total_energy_end) > &
702 : s% warn_when_large_rel_run_E_err) then
703 0 : write(*,2) 'WARNING: rel_run_E_err', &
704 0 : s% model_number, abs(s% cumulative_energy_error/s% total_energy_end)
705 : end if
706 : end if
707 :
708 11 : if (.not. (s% rotation_flag .or. s% u_flag .or. s% use_mass_corrections &
709 : .or. s% v_flag .or. s% m_center > 0 .or. s% star_mdot /= 0d0)) then
710 11 : tmp = abs(1d0 + s% total_gravitational_energy_end/s% virial_thm_P_avg)
711 11 : if (tmp > s% warn_when_large_virial_thm_rel_err) then
712 0 : write(*,2) 'WARNING: virial_thm_rel_err', &
713 0 : s% model_number, tmp, s% warn_when_large_virial_thm_rel_err, &
714 0 : abs(s% total_gravitational_energy_end), s% virial_thm_P_avg
715 : end if
716 : end if
717 :
718 11 : if (result == keep_going) then
719 10 : if (s% job% pgstar_flag) then
720 0 : will_read_pgstar_inlist = .false.
721 0 : if (s% pg% pgstar_interval <= 0) then
722 : will_read_pgstar_inlist = .true.
723 0 : else if(mod(s% model_number, s% pg% pgstar_interval) == 0) then
724 : will_read_pgstar_inlist = .true.
725 : end if
726 0 : if(will_read_pgstar_inlist) then
727 0 : call read_pgstar_inlist(s, inlist_fname, ierr)
728 0 : if (failed('read_pgstar_controls',ierr)) return
729 : end if
730 : end if
731 : end if
732 :
733 11 : if (result == keep_going) then
734 10 : result = s% extras_finish_step(id)
735 1 : else if (result == terminate) then
736 : ! call extras_finish_step one last time before terminate
737 1 : result = s% extras_finish_step(id)
738 1 : result = terminate
739 : end if
740 :
741 11 : if (result == keep_going) then
742 10 : if (dbg) write(*,*) 'call star_finish_step'
743 10 : result = star_finish_step(id, ierr)
744 10 : if (failed('star_finish_step',ierr)) return
745 : end if
746 :
747 11 : if (result == keep_going .and. s% job% pgstar_flag) then
748 0 : if (dbg) write(*,*) 'call update_pgstar_plots'
749 0 : call update_pgstar_plots(s, .false., ierr)
750 0 : if (failed('update_pgstar_plots',ierr)) return
751 : end if
752 :
753 11 : if (result == keep_going) then
754 10 : call adjust_tau_factor(s)
755 : if (s% L_nuc_burn_total/s% L_phot >= s% Lnuc_div_L_zams_limit &
756 10 : .and. .not. s% rotation_flag) then
757 10 : call do_rotation_near_zams(s,ierr)
758 10 : if (ierr /= 0) return
759 : end if
760 10 : if (s% rotation_flag) then
761 0 : call do_rotation(s,ierr)
762 0 : if (ierr /= 0) return
763 : end if
764 : end if
765 :
766 : end subroutine after_step_loop
767 :
768 :
769 1 : subroutine terminate_normal_evolve_loop(id, &
770 : dbg, result, ierr)
771 : integer, intent(in) :: id
772 : type (star_info), pointer :: s
773 : logical, intent(in) :: dbg
774 : integer, intent(out) :: result, ierr
775 : integer :: i
776 : include 'formats'
777 :
778 1 : call star_ptr(id, s, ierr)
779 1 : if (ierr/=0) return
780 :
781 1 : if (dbg) write(*,*) 'call star_pick_next_timestep'
782 1 : result = star_pick_next_timestep(id) ! for saved model if any
783 1 : if (dbg) write(*,*) 'call save_profile'
784 1 : call save_profile(id, 3, ierr)
785 1 : s% need_to_save_profiles_now = .false.
786 1 : s% need_to_update_history_now = .true.
787 1 : if (dbg) write(*,*) 'call star_finish_step'
788 1 : result = star_finish_step(id, ierr)
789 1 : if (failed('star_finish_step',ierr)) return
790 1 : if (s% job% save_photo_when_terminate .and. termination_code_string_okay()) &
791 1 : s% job% save_photo_number = s% model_number
792 1 : if (s% job% save_model_when_terminate .and. termination_code_string_okay()) &
793 0 : s% job% save_model_number = s% model_number
794 1 : if (s% job% save_pulse_data_when_terminate) &
795 0 : s% job% save_pulse_data_for_model_number = s% model_number
796 1 : if (s% job% write_profile_when_terminate) then
797 0 : if (len_trim(s% job% filename_for_profile_when_terminate) > 0) then
798 : call star_write_profile_info( &
799 : id, s% job% filename_for_profile_when_terminate, &
800 0 : ierr)
801 0 : if (failed('star_write_profile_info',ierr)) return
802 : else
803 0 : write(*,*) "filename_for_profile_when_terminate must be non empty"
804 0 : ierr = -1
805 0 : return
806 : end if
807 : end if
808 1 : if (s% job% show_retry_counts_when_terminate) then
809 59 : do i=1,numTlim
810 59 : if (s% dt_why_retry_count(i) > 0) then
811 0 : write(*,2) trim(dt_why_str(i)) // ' retries', s% dt_why_retry_count(i)
812 : end if
813 : end do
814 1 : write(*,'(A)')
815 : end if
816 1 : if (s% job% show_timestep_limit_counts_when_terminate) then
817 59 : do i=1,numTlim
818 59 : if (s% dt_why_count(i) > 0) then
819 1 : write(*,2) trim(dt_why_str(i)) // ' dt limit', s% dt_why_count(i)
820 : end if
821 : end do
822 1 : write(*,'(A)')
823 : end if
824 1 : call do_saves(id, ierr)
825 1 : if (failed('do_saves terminate_normal_evolve_loop',ierr)) return
826 :
827 : contains
828 :
829 1 : logical function termination_code_string_okay()
830 : integer :: j, n
831 1 : termination_code_string_okay = .true.
832 1 : if (s% termination_code == 0) return
833 : n = num_termination_code_strings
834 10 : j = maxval(len_trim(s% job% required_termination_code_string(1:n)))
835 1 : if (j == 0) return
836 0 : termination_code_string_okay = .false.
837 0 : do j=1,num_termination_code_strings
838 0 : if (s% job% required_termination_code_string(j) == &
839 0 : termination_code_str(s% termination_code)) then
840 1 : termination_code_string_okay = .true.
841 : return
842 : end if
843 : end do
844 : end function termination_code_string_okay
845 :
846 : end subroutine terminate_normal_evolve_loop
847 :
848 :
849 2 : subroutine after_evolve_loop(id, &
850 : do_free_star, ierr)
851 : integer, intent(in) :: id
852 : type (star_info), pointer :: s
853 : logical, intent(in) :: do_free_star
854 : integer, intent(out) :: ierr
855 :
856 1 : call star_ptr(id, s, ierr)
857 1 : if (ierr/=0) return
858 :
859 1 : if (s% doing_timing) then
860 0 : call system_clock(s% job% time1,s% job% clock_rate)
861 : s% job% elapsed_time = &
862 0 : dble(s% job% time1 - s% job% time0_initial) / s% job% clock_rate
863 0 : call show_times(id,s)
864 : end if
865 :
866 1 : if (s% result_reason /= result_reason_normal) then
867 : write(*, '(a)') 'terminated evolution: ' // &
868 0 : trim(result_reason_str(s% result_reason))
869 : end if
870 :
871 1 : if (s% termination_code > 0 .and. s% termination_code <= num_termination_codes) then
872 : write(*, '(a)') 'termination code: ' // &
873 1 : trim(termination_code_str(s% termination_code))
874 : end if
875 :
876 1 : if (s% job% pause_before_terminate) then
877 0 : write(*,'(a)') 'pause_before_terminate: hit RETURN to continue'
878 0 : read(*,*)
879 : end if
880 :
881 1 : call s% extras_after_evolve(id, ierr)
882 1 : if (failed('after_evolve_extras',ierr)) return
883 :
884 1 : if (s% result_reason == result_reason_normal) then
885 :
886 1 : if (s% job% pgstar_flag) &
887 : call update_pgstar_plots( &
888 : s, s% job% save_pgstar_files_when_terminate, &
889 0 : ierr)
890 1 : if (failed('update_pgstar_plots',ierr)) return
891 :
892 1 : call show_terminal_header(id, ierr)
893 1 : if (failed('show_terminal_header',ierr)) return
894 :
895 1 : call write_terminal_summary(id, ierr)
896 1 : if (failed('write_terminal_summary',ierr)) return
897 :
898 : end if
899 :
900 1 : if (len_trim(s% job% echo_at_end) > 0) then
901 0 : write(*,'(A)')
902 0 : write(*,'(a)') trim(s% job% echo_at_end)
903 0 : write(*,'(A)')
904 : end if
905 :
906 1 : if (do_free_star) then
907 1 : call free_star(id, ierr)
908 1 : if (failed('free_star',ierr)) return
909 : end if
910 :
911 : end subroutine after_evolve_loop
912 :
913 :
914 10 : subroutine adjust_tau_factor(s)
915 : type (star_info), pointer :: s
916 : include 'formats'
917 :
918 10 : if (s% job% adjust_tau_factor_to_surf_density .and. &
919 : s% job% base_for_adjust_tau_factor_to_surf_density > 0d0) then
920 0 : s% tau_factor = s% rho(1)/s% job% base_for_adjust_tau_factor_to_surf_density
921 : !write(*,1) 'adjust_tau_factor_to_surf_density', s% tau_factor
922 0 : s% need_to_setvars = .true.
923 : end if
924 :
925 10 : if (s% job% set_tau_factor_after_core_He_burn > 0 .and. &
926 : abs(s% tau_factor - s% job% set_to_this_tau_factor) > &
927 : 1d-6*max(s% tau_factor, s% job% set_to_this_tau_factor)) then
928 0 : if (check_for_after_He_burn(s, s% job% set_tau_factor_after_core_He_burn)) then
929 0 : s% tau_factor = s% job% set_to_this_tau_factor
930 0 : write(*,1) 'set_tau_factor_after_core_He_burn', s% tau_factor
931 0 : s% need_to_setvars = .true.
932 : end if
933 : end if
934 :
935 10 : if (s% job% set_tau_factor_after_core_C_burn > 0 .and. &
936 : abs(s% tau_factor - s% job% set_to_this_tau_factor) > &
937 : 1d-6*max(s% tau_factor, s% job% set_to_this_tau_factor)) then
938 0 : if (check_for_after_C_burn(s, s% job% set_tau_factor_after_core_C_burn)) then
939 0 : s% tau_factor = s% job% set_to_this_tau_factor
940 0 : write(*,1) 'set_tau_factor_after_core_C_burn', s% tau_factor
941 0 : s% need_to_setvars = .true.
942 : end if
943 : end if
944 :
945 10 : if (s% job% relax_tau_factor_after_core_He_burn > 0 .and. &
946 : abs(s% tau_factor - s% job% relax_to_this_tau_factor) > &
947 : 1d-6*max(s% tau_factor, s% job% relax_to_this_tau_factor)) then
948 0 : if (check_for_after_He_burn(s, s% job% relax_tau_factor_after_core_He_burn)) &
949 0 : call relax_tau_factor(s)
950 : end if
951 :
952 10 : if (s% job% relax_tau_factor_after_core_C_burn > 0 .and. &
953 : abs(s% tau_factor - s% job% relax_to_this_tau_factor) > &
954 : 1d-6*max(s% tau_factor, s% job% relax_to_this_tau_factor)) then
955 0 : if (check_for_after_C_burn(s, s% job% relax_tau_factor_after_core_C_burn)) &
956 0 : call relax_tau_factor(s)
957 : end if
958 :
959 :
960 10 : end subroutine adjust_tau_factor
961 :
962 :
963 0 : subroutine do_rotation(s,ierr)
964 : type (star_info), pointer :: s
965 : integer, intent(out) :: ierr
966 : include 'formats'
967 0 : ierr = 0
968 :
969 0 : if (s% model_number <= s% job% set_surf_rotation_v_step_limit) then
970 0 : s% job% new_omega = s% job% new_surface_rotation_v*1d5/(s% photosphere_r*Rsun)
971 0 : write(*,2) 'surface_rotation_v', s% model_number, s% job% new_surface_rotation_v
972 0 : write(*,2) 'omega', s% model_number, s% job% new_omega
973 0 : call star_set_uniform_omega(s% id, s% job% new_omega, ierr)
974 0 : if (failed('star_set_uniform_omega',ierr)) return
975 :
976 0 : else if (s% model_number <= s% job% set_omega_step_limit) then
977 0 : write(*,2) 'omega', s% model_number, s% job% new_omega
978 0 : if (failed('star_surface_omega_crit',ierr)) return
979 0 : call star_set_uniform_omega(s% id, s% job% new_omega, ierr)
980 0 : if (failed('star_set_uniform_omega',ierr)) return
981 :
982 0 : else if (s% model_number <= s% job% set_omega_div_omega_crit_step_limit) then
983 : s% job% new_omega = &
984 0 : s% job% new_omega_div_omega_crit*star_surface_omega_crit(s% id, ierr)
985 0 : write(*,2) 'omega_div_omega_crit', &
986 0 : s% model_number, s% job% new_omega_div_omega_crit
987 0 : write(*,2) 'omega', s% model_number, s% job% new_omega
988 0 : if (failed('star_surface_omega_crit',ierr)) return
989 0 : call star_set_uniform_omega(s% id, s% job% new_omega, ierr)
990 0 : if (failed('star_set_uniform_omega',ierr)) return
991 : end if
992 : end subroutine do_rotation
993 :
994 :
995 10 : subroutine do_rotation_near_zams(s,ierr)
996 : type (star_info), pointer :: s
997 : integer, intent(out) :: ierr
998 : include 'formats'
999 10 : ierr = 0
1000 :
1001 10 : if (s% job% set_near_zams_surface_rotation_v_steps > 0 .and. &
1002 : s% job% new_surface_rotation_v /= 0d0) then
1003 0 : s% job% new_rotation_flag = .true.
1004 0 : call star_set_rotation_flag(s% id, s% job% new_rotation_flag, ierr)
1005 0 : if (failed('star_set_rotation_flag',ierr)) return
1006 : s% job% set_surf_rotation_v_step_limit = &
1007 0 : s% model_number + s% job% set_near_zams_surface_rotation_v_steps - 1
1008 0 : write(*,2) 'near zams: set_surf_rotation_v_step_limit', &
1009 0 : s% job% set_surf_rotation_v_step_limit
1010 :
1011 10 : else if (s% job% set_near_zams_omega_steps > 0 .and. &
1012 : s% job% new_omega /= 0d0) then
1013 0 : s% job% new_rotation_flag = .true.
1014 0 : call star_set_rotation_flag(s% id, s% job% new_rotation_flag, ierr)
1015 0 : if (failed('star_set_rotation_flag',ierr)) return
1016 : s% job% set_omega_step_limit = &
1017 0 : s% model_number + s% job% set_near_zams_omega_steps - 1
1018 0 : write(*,2) 'near zams: set_omega_step_limit', s% job% set_omega_step_limit
1019 :
1020 10 : else if (s% job% set_near_zams_omega_div_omega_crit_steps > 0 .and. &
1021 : s% job% new_omega_div_omega_crit /= 0d0) then
1022 0 : s% job% new_rotation_flag = .true.
1023 0 : call star_set_rotation_flag(s% id, s% job% new_rotation_flag, ierr)
1024 0 : if (failed('star_set_rotation_flag',ierr)) return
1025 : s% job% set_omega_div_omega_crit_step_limit = &
1026 0 : s% model_number + s% job% set_near_zams_omega_div_omega_crit_steps - 1
1027 0 : write(*,2) 'near zams: set_omega_div_omega_crit_step_limit', &
1028 0 : s% job% set_omega_div_omega_crit_step_limit
1029 :
1030 10 : else if (s% job% near_zams_relax_omega .and. &
1031 : s% job% new_omega /= 0d0) then
1032 0 : s% job% new_rotation_flag = .true.
1033 0 : call star_set_rotation_flag(s% id, s% job% new_rotation_flag, ierr)
1034 0 : if (failed('star_set_rotation_flag',ierr)) return
1035 0 : write(*,2) 'new_omega', s% model_number, s% job% new_omega
1036 : call star_relax_uniform_omega( &
1037 : s% id, relax_to_new_omega, s% job% new_omega, s% job% num_steps_to_relax_rotation,&
1038 0 : s% job% relax_omega_max_yrs_dt, ierr)
1039 0 : if (failed('star_relax_uniform_omega',ierr)) return
1040 :
1041 10 : else if (s% job% near_zams_relax_omega_div_omega_crit .and. &
1042 : s% job% new_omega_div_omega_crit /= 0d0) then
1043 0 : s% job% new_rotation_flag = .true.
1044 0 : call star_set_rotation_flag(s% id, s% job% new_rotation_flag, ierr)
1045 0 : if (failed('star_set_rotation_flag',ierr)) return
1046 0 : write(*,2) 'new_omega_div_omega_crit', &
1047 0 : s% model_number, s% job% new_omega_div_omega_crit
1048 : call star_relax_uniform_omega( &
1049 : s% id, relax_to_new_omega_div_omega_crit, s% job% new_omega_div_omega_crit, &
1050 : s% job% num_steps_to_relax_rotation,&
1051 0 : s% job% relax_omega_max_yrs_dt, ierr)
1052 0 : if (failed('star_relax_uniform_omega',ierr)) return
1053 :
1054 10 : else if (s% job% near_zams_relax_initial_surface_rotation_v .and. &
1055 : s% job% new_surface_rotation_v /= 0d0) then
1056 0 : s% job% new_rotation_flag = .true.
1057 0 : call star_set_rotation_flag(s% id, s% job% new_rotation_flag, ierr)
1058 0 : if (failed('star_set_rotation_flag',ierr)) return
1059 0 : write(*,2) 'new_surface_rotation_v', &
1060 0 : s% model_number, s% job% new_surface_rotation_v
1061 : call star_relax_uniform_omega( &
1062 : s% id, relax_to_new_surface_rotation_v, s% job% new_surface_rotation_v, &
1063 : s% job% num_steps_to_relax_rotation,&
1064 0 : s% job% relax_omega_max_yrs_dt, ierr)
1065 0 : if (failed('star_relax_uniform_omega',ierr)) return
1066 :
1067 : end if
1068 : end subroutine do_rotation_near_zams
1069 :
1070 :
1071 0 : subroutine relax_tau_factor(s)
1072 : type (star_info), pointer :: s
1073 : real(dp) :: next
1074 : include 'formats'
1075 0 : write(*,*) 'relax_to_this_tau_factor < s% tau_factor', &
1076 0 : s% job% relax_to_this_tau_factor < s% tau_factor
1077 0 : write(*,1) 'relax_to_this_tau_factor', s% job% relax_to_this_tau_factor
1078 0 : write(*,1) 's% tau_factor', s% tau_factor
1079 0 : if (s% job% relax_to_this_tau_factor < s% tau_factor) then
1080 0 : next = exp10(safe_log10(s% tau_factor) - s% job% dlogtau_factor)
1081 0 : if (next < s% job% relax_to_this_tau_factor) &
1082 0 : next = s% job% relax_to_this_tau_factor
1083 : else
1084 0 : next = exp10(safe_log10(s% tau_factor) + s% job% dlogtau_factor)
1085 0 : if (next > s% job% relax_to_this_tau_factor) &
1086 0 : next = s% job% relax_to_this_tau_factor
1087 : end if
1088 0 : if (next /= s% tau_factor) then
1089 0 : s% tau_factor = next
1090 0 : write(*,1) 'relax_tau_factor', next, s% job% relax_to_this_tau_factor
1091 0 : s% need_to_setvars = .true.
1092 : end if
1093 0 : end subroutine relax_tau_factor
1094 :
1095 1 : subroutine check_if_want_to_stop_warnings(s)
1096 : use utils_lib
1097 : type (star_info), pointer :: s
1098 : character (len=200) :: fname
1099 : integer :: iounit, ierr
1100 1 : ierr = 0
1101 1 : if (s% warn_when_large_rel_run_E_err < 1d2) then
1102 1 : fname = trim(mesa_dir) // '/stop_warnings_for_rel_E_err'
1103 : open(newunit=iounit, file=trim(fname), &
1104 1 : status='old', action='read', iostat=ierr)
1105 1 : if (ierr == 0) then
1106 0 : close(iounit)
1107 0 : s% warn_when_large_rel_run_E_err = 1d99
1108 0 : write(*,*) 'turn off warnings for rel_run_E_err'
1109 : end if
1110 : end if
1111 : ierr = 0
1112 1 : end subroutine check_if_want_to_stop_warnings
1113 :
1114 :
1115 11 : logical function stop_is_requested(s)
1116 : type (star_info), pointer :: s
1117 : logical :: file_exists
1118 11 : stop_is_requested = .false.
1119 11 : if (mod(s% model_number,100) /= 0) return
1120 1 : if (len_trim(s% job% stop_if_this_file_exists) == 0) return
1121 1 : inquire(file=trim(s% job% stop_if_this_file_exists), exist=file_exists)
1122 1 : if (.not. file_exists) return
1123 : write(*,*) 'stopping because found file ' // &
1124 0 : trim(s% job% stop_if_this_file_exists)
1125 0 : stop_is_requested = .true.
1126 0 : end function stop_is_requested
1127 :
1128 :
1129 118 : logical function failed(str,ierr)
1130 : character (len=*), intent(in) :: str
1131 : integer, intent(in) :: ierr
1132 118 : failed = (ierr /= 0)
1133 118 : if (failed) write(*, *) trim(str) // ' ierr', ierr
1134 118 : end function failed
1135 :
1136 :
1137 0 : subroutine show_times(id, s)
1138 : use utils_lib, only: utils_OMP_GET_MAX_THREADS
1139 : use num_lib, only: qsort
1140 :
1141 : integer, intent(in) :: id
1142 : type (star_info), pointer :: s
1143 :
1144 : integer, parameter :: max_num_items = 50
1145 : character(len=60) :: item_names(max_num_items)
1146 : real(dp) :: item_values(max_num_items)
1147 : integer, target :: index_arry(max_num_items)
1148 : integer, pointer :: index(:)
1149 : integer :: ierr, omp_num_threads, item_num, num_items, i, j
1150 : real(dp) :: total, tmp
1151 : include 'formats'
1152 0 : ierr = 0
1153 0 : omp_num_threads = utils_OMP_GET_MAX_THREADS()
1154 : s% time_total = s% job% check_before_step_timing + &
1155 0 : s% job% check_step_loop_timing + s% job% check_after_step_timing
1156 :
1157 0 : write(*,'(A)')
1158 0 : write(*,'(a50,i18)') 'nz', s% nz
1159 0 : write(*,'(a50,i18)') 'nvar_total', s% nvar_total
1160 0 : write(*,'(a50,i18)') trim(s% net_name) // ' species', s% species
1161 0 : write(*,'(a50,i18)') 'total_num_solver_iterations', &
1162 0 : s% total_num_solver_iterations
1163 0 : write(*,'(a50,i18)') 'timing_num_get_eos_calls', &
1164 0 : s% timing_num_get_eos_calls
1165 0 : write(*,'(a50,i18)') 'timing_num_solve_eos_calls', &
1166 0 : s% timing_num_solve_eos_calls
1167 0 : write(*,'(a50,i18)') 'timing_num_get_kap_calls', &
1168 0 : s% timing_num_get_kap_calls
1169 0 : write(*,'(A)')
1170 0 : write(*,'(a50,i18)') 'threads', omp_num_threads
1171 0 : total = 0
1172 0 : item_num = 0
1173 0 : call save1('remesh', s% time_remesh, total)
1174 0 : call save1('adjust_mass', s% time_adjust_mass, total)
1175 0 : call save1('conv_premix', s% time_conv_premix, total)
1176 0 : call save1('element_diffusion', s% time_element_diffusion, total)
1177 0 : call save1('burn', s% time_solve_burn, total)
1178 0 : call save1('mix', s% time_solve_mix, total)
1179 0 : call save1('solve', s% time_struct_burn_mix, total)
1180 0 : call save1('matrix', s% time_solver_matrix, total)
1181 0 : call save1('omega_mix', s% time_solve_omega_mix, total)
1182 0 : call save1('eos', s% time_eos, total)
1183 0 : call save1('neu_and_kap', s% time_neu_kap, total)
1184 0 : call save1('net', s% time_nonburn_net, total)
1185 0 : call save1('mlt', s% time_mlt, total)
1186 0 : call save1('hydro_vars', s% time_set_hydro_vars, total)
1187 0 : call save1('mixing_info', s% time_set_mixing_info, total)
1188 0 : call save1('evolve_step', s% time_evolve_step, total)
1189 0 : call save1('run1_star', s% job% elapsed_time - total, total)
1190 0 : tmp = 0
1191 0 : call save1('total', total, tmp)
1192 :
1193 0 : num_items = item_num
1194 0 : index(1:num_items) => index_arry(1:num_items)
1195 0 : call qsort(index, num_items, item_values)
1196 :
1197 0 : write(*,'(A)')
1198 0 : write(*,'(A)')
1199 0 : do i=1,num_items
1200 0 : j = index(num_items+1-i)
1201 0 : if (item_values(j) == 0d0) cycle
1202 0 : write(*,'(a50,2f9.3)') trim(item_names(j)), &
1203 0 : item_values(j), item_values(j)/total
1204 0 : if (j == num_items) write(*,*)
1205 : end do
1206 :
1207 0 : if (s% job% step_loop_timing/s% job% elapsed_time < 0.9d0) then
1208 0 : write(*,'(A)')
1209 0 : write(*,'(A)')
1210 0 : write(*,1) 'before_step', s% job% before_step_timing/s% job% elapsed_time
1211 0 : write(*,1) 'step_loop', s% job% step_loop_timing/s% job% elapsed_time
1212 0 : write(*,1) 'after_step', s% job% after_step_timing/s% job% elapsed_time
1213 0 : write(*,'(A)')
1214 : end if
1215 0 : write(*,'(A)')
1216 0 : write(*,'(A)')
1217 :
1218 :
1219 : contains
1220 :
1221 :
1222 0 : subroutine save1(name, value, total)
1223 : use utils_lib, only: is_bad_num
1224 : character (len=*), intent(in) :: name
1225 : real(dp), intent(in) :: value
1226 : real(dp), intent(inout) :: total
1227 : include 'formats'
1228 0 : item_num = item_num + 1
1229 0 : item_names(item_num) = name
1230 0 : item_values(item_num) = value
1231 0 : total = total + value
1232 0 : end subroutine save1
1233 :
1234 :
1235 : end subroutine show_times
1236 :
1237 :
1238 11 : subroutine do_saves(id, ierr)
1239 : integer, intent(in) :: id
1240 : type (star_info), pointer :: s
1241 :
1242 : integer :: ierr
1243 : ierr = 0
1244 :
1245 11 : call star_ptr(id, s, ierr)
1246 11 : if (ierr/=0) return
1247 :
1248 11 : if (s% model_number == s% job% save_model_number) then
1249 0 : call star_write_model(id, s% job% save_model_filename, ierr)
1250 0 : if (failed('star_write_model',ierr)) return
1251 0 : write(*, *) 'model saved to ' // trim(s% job% save_model_filename)
1252 : end if
1253 :
1254 11 : if (s% model_number == s% job% save_photo_number) then
1255 1 : call star_write_photo(id, s% job% save_photo_filename, ierr)
1256 1 : if (failed('star_write_photo',ierr)) return
1257 1 : if (len_trim(s% job% save_photo_filename) > 0) &
1258 0 : write(*, *) 'photo saved to ' // trim(s% job% save_photo_filename)
1259 : end if
1260 :
1261 11 : if (s% model_number == s% job% save_pulse_data_for_model_number) then
1262 : call star_export_pulse_data(id, s%pulse_data_format, s%job%save_pulse_data_filename, &
1263 : s%add_center_point_to_pulse_data, s%keep_surface_point_for_pulse_data, &
1264 0 : s%add_atmosphere_to_pulse_data, ierr)
1265 0 : if (failed('star_export_pulse_data',ierr)) return
1266 : write(*, *) 'pulsation data saved to ' // &
1267 0 : trim(s% job% save_pulse_data_filename)
1268 : end if
1269 :
1270 11 : if (s% model_number == s% job% profile_model_number) then
1271 0 : write(*, '(a, i7)') 'save profile for model number', s% model_number
1272 0 : call save_profile(id, 3, ierr)
1273 0 : if (failed('save_profile',ierr)) return
1274 : end if
1275 :
1276 : end subroutine do_saves
1277 :
1278 :
1279 : subroutine write_colors_info(id, s, ierr)
1280 : use colors_lib
1281 : use colors_def
1282 : integer, intent(in) :: id
1283 : type (star_info), pointer :: s
1284 : integer, intent(out) :: ierr
1285 :
1286 : !TODO: implement me
1287 :
1288 : end subroutine write_colors_info
1289 :
1290 :
1291 : subroutine read_masses(filename, masses, nmasses, ierr)
1292 : character (len=*), intent(in) :: filename
1293 : real(dp), pointer, intent(inout) :: masses(:)
1294 : integer, intent(out) :: nmasses, ierr
1295 : call read_items(filename, masses, nmasses, 'masses', ierr)
1296 : end subroutine read_masses
1297 :
1298 :
1299 : subroutine read_items(filename, items, nitems, name, ierr)
1300 : use utils_lib
1301 : use utils_def
1302 : character (len=*), intent(in) :: filename, name
1303 : real(dp), pointer, intent(inout) :: items(:)
1304 : integer, intent(out) :: nitems, ierr
1305 :
1306 : integer :: iounit, n, i, t, capacity
1307 : character (len=strlen) :: buffer, string
1308 :
1309 : nitems = 0
1310 : if (.not. associated(items)) then
1311 : capacity = 10
1312 : allocate(items(capacity))
1313 : else
1314 : capacity = size(items,dim=1)
1315 : end if
1316 :
1317 : ierr = 0
1318 :
1319 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
1320 : if (ierr /= 0) then
1321 : write(*,*) 'failed to open file ' // trim(filename)
1322 : return
1323 : end if
1324 :
1325 : n = 0
1326 : i = 0
1327 :
1328 : do
1329 : t = token(iounit, n, i, buffer, string)
1330 : select case(t)
1331 : case(name_token)
1332 : if (string == name) then
1333 : call do_read_items(ierr)
1334 : if (ierr /= 0) then
1335 : return
1336 : end if
1337 : exit ! for now, nothing else to be read
1338 : end if
1339 : call error; return
1340 : case(eof_token)
1341 : exit
1342 : case default
1343 : call error; return
1344 : end select
1345 :
1346 : end do
1347 :
1348 : close(iounit)
1349 :
1350 : contains
1351 :
1352 :
1353 : subroutine error
1354 : ierr = -1
1355 : write(*,*) 'error in reading file' // trim(filename)
1356 : close(iounit)
1357 : end subroutine error
1358 :
1359 :
1360 : subroutine do_read_items(ierr)
1361 : integer, intent(out) :: ierr
1362 : real(dp) :: mass
1363 : ierr = 0
1364 : t = token(iounit, n, i, buffer, string)
1365 : if (t /= left_paren_token) then
1366 : call error; return
1367 : end if
1368 : mass_loop: do
1369 : t = token(iounit, n, i, buffer, string)
1370 : if (t /= name_token) then
1371 : call error; return
1372 : end if
1373 : read(string,fmt=*,iostat=ierr) mass
1374 : if (ierr /= 0) then
1375 : call error; return
1376 : end if
1377 : nitems = nitems+1
1378 : if (nitems > capacity) then
1379 : capacity = capacity + 10
1380 : call realloc_double(items,capacity,ierr)
1381 : if (ierr /= 0) then
1382 : call error; return
1383 : end if
1384 : end if
1385 : items(nitems) = mass
1386 : t = token(iounit, n, i, buffer, string)
1387 : if (t == right_paren_token) exit mass_loop
1388 : if (t /= comma_token) then
1389 : call error; return
1390 : end if
1391 : end do mass_loop
1392 : end subroutine do_read_items
1393 :
1394 :
1395 : end subroutine read_items
1396 :
1397 :
1398 0 : subroutine do_report_mass_not_fe56(s)
1399 : use const_def, only: dp
1400 : type (star_info), pointer :: s
1401 : integer :: k, fe56
1402 : real(dp) :: sumdq
1403 : include 'formats'
1404 0 : fe56 = s% net_iso(ife56)
1405 0 : if (fe56 == 0) return
1406 0 : sumdq = 0
1407 0 : do k = 1, s% nz
1408 0 : sumdq = sumdq + s% dq(k)*(1-s% xa(fe56,k))
1409 : end do
1410 0 : write(*,1) 'R', s% r(1)
1411 0 : write(*,1) 'g', s% cgrav(1)*s% mstar/(s% r(1)*s% r(1))
1412 0 : write(*,1) 'mass non fe56', s% xmstar*sumdq, sumdq
1413 0 : write(*,1) 'M_center (Msun)', s% M_center/Msun
1414 0 : write(*,1) 'xmstar (g)', s% xmstar
1415 0 : do k=1,s% nz
1416 0 : if (fe56 == maxloc(s% xa(:,k),dim=1)) then
1417 0 : write(*,2) 'mass exterior to fe56 (g)', k, (1d0 - s% q(k))*s% xmstar
1418 0 : write(*,2) 'mass coord top of fe56 (g)', k, s% q(k)*s% xmstar
1419 0 : return
1420 : end if
1421 : end do
1422 : end subroutine do_report_mass_not_fe56
1423 :
1424 :
1425 0 : subroutine do_report_cell_for_xm(s)
1426 : use const_def, only: dp
1427 : type (star_info), pointer :: s
1428 : integer :: k
1429 : real(dp) :: sumdq, dq
1430 : include 'formats'
1431 0 : dq = s% job% report_cell_for_xm/s% xmstar
1432 0 : if (dq > 1) then
1433 0 : write(*,2) 'report_cell_for_xm > xmstar', s% nz
1434 0 : return
1435 : end if
1436 0 : sumdq = 0
1437 0 : do k = 1, s% nz
1438 0 : sumdq = sumdq + s% dq(k)
1439 0 : if (sumdq >= dq) then
1440 0 : write(*,'(A)')
1441 0 : write(*,2) 'total mass in cells from 1 to k', k, sumdq*s% xmstar
1442 0 : write(*,2) 'logT(k)', k, s% lnT(k)/ln10
1443 0 : write(*,2) 'logRho(k)', k, s% lnd(k)/ln10
1444 0 : write(*,2) 'entropy(k)', k, exp(s% lnS(k))*amu/kerg
1445 0 : write(*,2) 'xmstar*q(k)', k, s% xmstar*s% q(k)
1446 0 : write(*,2) 'q(k)', k, s% q(k)
1447 0 : write(*,'(A)')
1448 0 : return
1449 : end if
1450 : end do
1451 0 : write(*,2) 'total mass in cells from 1 to nz', s% nz, s% xmstar
1452 : end subroutine do_report_cell_for_xm
1453 :
1454 2 : subroutine set_rate_factors(id, ierr)
1455 : use net_lib, only: get_net_reaction_table_ptr
1456 : use rates_lib, only: rates_reaction_id
1457 : integer, intent(in) :: id
1458 : integer, intent(out) :: ierr
1459 : type (star_info), pointer :: s
1460 : integer :: j, i, ir
1461 2 : integer, pointer :: net_reaction_ptr(:)
1462 :
1463 : include 'formats'
1464 :
1465 : ierr = 0
1466 2 : call star_ptr(id, s, ierr)
1467 4 : if (ierr /= 0) return
1468 :
1469 636 : s% rate_factors(:) = 1
1470 2 : if (s% job% num_special_rate_factors <= 0) return
1471 :
1472 : ! Dont error if we are changing net
1473 0 : if ((s% job% change_initial_net .or. s% job% change_net) .and. &
1474 : trim(s% job% new_net_name)/=trim(s% net_name)) then
1475 : !write(*,*) "Not changing special rates until net change"
1476 : return
1477 : end if
1478 :
1479 :
1480 0 : call get_net_reaction_table_ptr(s% net_handle, net_reaction_ptr, ierr)
1481 0 : if (ierr /= 0) return
1482 :
1483 0 : do i=1,s% job% num_special_rate_factors
1484 0 : if (len_trim(s% job% reaction_for_special_factor(i)) == 0) cycle
1485 0 : ir = rates_reaction_id(s% job% reaction_for_special_factor(i))
1486 0 : j = 0
1487 0 : if (ir > 0) j = net_reaction_ptr(ir)
1488 0 : if (j <= 0) then
1489 : write(*,*) 'Failed to find reaction_for_special_factor ' // &
1490 0 : trim(s% job% reaction_for_special_factor(i)), &
1491 0 : j, s% job% special_rate_factor(i)
1492 0 : cycle
1493 : end if
1494 0 : s% rate_factors(j) = s% job% special_rate_factor(i)
1495 : write(*,*) 'set special rate factor for ' // &
1496 0 : trim(s% job% reaction_for_special_factor(i)), &
1497 0 : j, s% job% special_rate_factor(i)
1498 : end do
1499 :
1500 2 : end subroutine set_rate_factors
1501 :
1502 :
1503 1 : subroutine do_star_job_controls_before(id, s, restart, ierr)
1504 :
1505 : use rates_lib, only: rates_warning_init
1506 : use atm_support, only: get_atm_tau_base
1507 :
1508 : integer, intent(in) :: id
1509 : type (star_info), pointer :: s
1510 : logical, intent(in) :: restart
1511 : integer, intent(out) :: ierr
1512 : logical, parameter :: kap_use_cache = .true.
1513 : include 'formats'
1514 :
1515 : ierr = 0
1516 :
1517 1 : s% set_rate_factors => set_rate_factors ! will be called after net is defined
1518 :
1519 1 : call get_atm_tau_base(s, s% tau_base, ierr)
1520 1 : if (failed('atm_tau_base',ierr)) return
1521 :
1522 : call rates_warning_init( &
1523 1 : s% warn_rates_for_high_temp, s% max_safe_logT_for_rates)
1524 :
1525 : end subroutine do_star_job_controls_before
1526 :
1527 :
1528 4 : subroutine do_read_star_job_and_return_id(filename, id, ierr)
1529 : character(*), intent(in) :: filename
1530 : integer, intent(out) :: id
1531 : integer, intent(out) :: ierr
1532 : type (star_info), pointer :: s
1533 : character(len=strlen) :: inlist_fname
1534 :
1535 : include 'formats'
1536 1 : ierr = 0
1537 :
1538 1 : if (id_from_read_star_job /= 0) then
1539 0 : write(*,2) 'id_from_read_star_job', id_from_read_star_job
1540 0 : ierr = -1
1541 0 : return
1542 : end if
1543 :
1544 1 : call alloc_star(id, ierr)
1545 1 : if (ierr /= 0) then
1546 0 : write(*,*) 'do_read_star_job failed in alloc_star'
1547 0 : return
1548 : end if
1549 :
1550 1 : call star_ptr(id, s, ierr)
1551 1 : if (ierr /= 0) then
1552 0 : write(*,*) 'do_read_star_job failed in star_ptr'
1553 0 : return
1554 : end if
1555 :
1556 1 : call resolve_inlist_fname(inlist_fname,filename)
1557 1 : call read_star_job(s, inlist_fname, ierr)
1558 1 : if (ierr /= 0) then
1559 0 : write(*,*) 'ierr from read_star_job ' // trim(inlist_fname)
1560 0 : return
1561 : end if
1562 :
1563 1 : id_from_read_star_job = id
1564 :
1565 1 : if (s% job% save_star_job_namelist) then
1566 0 : call write_star_job(s, s% job% star_job_namelist_name, ierr)
1567 0 : if (ierr /= 0) then
1568 : write(*,*) 'ierr from write_star_job ' // &
1569 0 : trim(s% job% star_job_namelist_name)
1570 0 : return
1571 : end if
1572 : end if
1573 :
1574 : end subroutine do_read_star_job_and_return_id
1575 :
1576 : ! in a perfect world, we'd pass s as an arg to this routine.
1577 : ! but for backward compatibility for a large number of users
1578 : ! we do it this strange way instead.
1579 1 : subroutine do_read_star_job(filename, ierr)
1580 : character(*), intent(in) :: filename
1581 : integer, intent(out) :: ierr
1582 : integer :: id
1583 1 : call do_read_star_job_and_return_id(filename, id, ierr)
1584 1 : end subroutine do_read_star_job
1585 :
1586 :
1587 1 : subroutine do_load1_star(id, s, restart, restart_filename, ierr)
1588 : integer, intent(in) :: id
1589 : type (star_info), pointer :: s
1590 : logical, intent(in) :: restart
1591 : character (len=*), intent(in) :: restart_filename
1592 : integer, intent(out) :: ierr
1593 :
1594 1 : if (restart) then
1595 0 : call star_load_restart_photo(id, restart_filename, ierr)
1596 0 : if (failed('star_load_restart_photo',ierr)) return
1597 1 : else if (s% job% load_saved_photo) then
1598 0 : write(*,'(a)') 'load saved photo ' // trim(s% job% saved_photo_name)
1599 0 : write(*,'(A)')
1600 0 : call star_load_restart_photo(id, s% job% saved_photo_name, ierr)
1601 0 : if (failed('star_load_restart_photo',ierr)) return
1602 1 : else if (s% job% load_saved_model) then
1603 0 : if (s% job% create_merger_model) then
1604 0 : write(*,*) 'you have both load_saved_model and create_merger_model set true'
1605 0 : write(*,*) 'please pick one and try again'
1606 0 : call mesa_error(__FILE__,__LINE__)
1607 : end if
1608 0 : if (s% job% create_pre_main_sequence_model) then
1609 : write(*,*) 'you have both load_saved_model and ' // &
1610 0 : 'create_pre_main_sequence_model set true'
1611 0 : write(*,*) 'please pick one and try again'
1612 0 : call mesa_error(__FILE__,__LINE__)
1613 : end if
1614 0 : if (s% job% create_initial_model) then
1615 0 : write(*,*) 'you have both load_saved_model and create_initial_model set true'
1616 0 : write(*,*) 'please pick one and try again'
1617 0 : call mesa_error(__FILE__,__LINE__)
1618 : end if
1619 0 : write(*,'(a)') 'load saved model ' // trim(s% job% load_model_filename)
1620 0 : write(*,'(A)')
1621 0 : call star_read_model(id, s% job% load_model_filename, ierr)
1622 0 : if (failed('star_read_model',ierr)) return
1623 1 : else if (s% job% create_merger_model) then
1624 0 : call create_merger_model(s, ierr)
1625 0 : if (failed('create_merger_model',ierr)) return
1626 1 : else if (s% job% create_pre_main_sequence_model) then
1627 0 : if (.not. restart) write(*, *) 'create pre-main-sequence model'
1628 0 : if (s% job% create_initial_model) then
1629 : write(*,*) 'you have both create_pre_main_sequence_model ' // &
1630 0 : 'and create_initial_model set true'
1631 0 : write(*,*) 'please pick one and try again'
1632 0 : call mesa_error(__FILE__,__LINE__)
1633 : end if
1634 : call star_create_pre_ms_model( &
1635 : id, s% job% pre_ms_T_c, s% job% pre_ms_guess_rho_c, &
1636 : s% job% pre_ms_d_log10_P, s% job% pre_ms_logT_surf_limit, &
1637 : s% job% pre_ms_logP_surf_limit, s% job% initial_zfracs, &
1638 : s% job% dump_missing_metals_into_heaviest, &
1639 : (s% job% change_net .or. (s% job% change_initial_net .and. .not. restart)), &
1640 0 : s% job% new_net_name, s% job% pre_ms_relax_num_steps, ierr)
1641 0 : if (failed('star_create_pre_ms_model',ierr)) return
1642 1 : else if (s% job% create_RSP_model) then
1643 0 : if (.not. restart) write(*, *) 'create initial RSP model'
1644 0 : call star_create_RSP_model(id, ierr)
1645 0 : if (failed('star_create_RSP_model',ierr)) return
1646 1 : else if (s% job% create_RSP2_model) then
1647 0 : if (.not. restart) write(*, *) 'create initial RSP2 model'
1648 0 : call star_create_RSP2_model(id, ierr)
1649 0 : if (failed('star_create_RSP_model',ierr)) return
1650 1 : else if (s% job% create_initial_model) then
1651 0 : if (.not. restart) write(*, *) 'create initial model'
1652 0 : if (s% job% create_pre_main_sequence_model) then
1653 : write(*,*) 'you have both create_initial_model and ' // &
1654 0 : 'create_pre_main_sequence_model set true'
1655 0 : write(*,*) 'please pick one and try again'
1656 0 : call mesa_error(__FILE__,__LINE__)
1657 : end if
1658 : call star_create_initial_model(id, &
1659 : s% job% radius_in_cm_for_create_initial_model, &
1660 : s% job% mass_in_gm_for_create_initial_model, &
1661 : s% job% center_logP_1st_try_for_create_initial_model, &
1662 : s% job% entropy_1st_try_for_create_initial_model, &
1663 : s% job% max_tries_for_create_initial_model, &
1664 : s% job% abs_e01_tolerance_for_create_initial_model, &
1665 : s% job% abs_e02_tolerance_for_create_initial_model, &
1666 : s% job% initial_zfracs, &
1667 : s% job% dump_missing_metals_into_heaviest, &
1668 : (s% job% change_net .or. (s% job% change_initial_net .and. .not. restart)), &
1669 : s% job% new_net_name, s% job% initial_model_relax_num_steps, &
1670 : s% job% initial_model_eps, &
1671 0 : ierr)
1672 0 : if (failed('star_create_initial_model',ierr)) return
1673 : else
1674 1 : call star_load_zams(id, ierr)
1675 1 : if (failed('star_load_zams',ierr)) return
1676 : end if
1677 :
1678 : end subroutine do_load1_star
1679 :
1680 :
1681 0 : subroutine create_merger_model(s, ierr)
1682 : use ctrls_io, only : store_controls
1683 : type (star_info), pointer :: s
1684 : integer, intent(out) :: ierr
1685 :
1686 : integer :: id, id_aux, i, j, k
1687 : type (star_info), pointer :: s_aux
1688 0 : real(dp), pointer :: xq(:), xa(:,:)
1689 : real(dp) :: total_mass, partial_mass
1690 :
1691 : include 'formats'
1692 0 : ierr = 0
1693 0 : id = s% id
1694 :
1695 0 : if (s% job% create_pre_main_sequence_model) then
1696 : write(*,*) 'you have both load_saved_model and ' // &
1697 0 : 'create_pre_main_sequence_model set true'
1698 0 : write(*,*) 'please pick one and try again'
1699 0 : call mesa_error(__FILE__,__LINE__)
1700 : end if
1701 0 : if (s% job% create_initial_model) then
1702 0 : write(*,*) 'you have both load_saved_model and create_initial_model set true'
1703 0 : write(*,*) 'please pick one and try again'
1704 0 : call mesa_error(__FILE__,__LINE__)
1705 : end if
1706 : !load first star
1707 0 : call star_read_model(id, s% job% saved_model_for_merger_1, ierr)
1708 0 : if (failed('star_read_model',ierr)) return
1709 :
1710 : !load second star
1711 0 : call alloc_star(id_aux, ierr)
1712 0 : if (failed('alloc_star',ierr)) return
1713 0 : call star_ptr(id_aux, s_aux, ierr)
1714 0 : if (failed('star_ptr',ierr)) return
1715 0 : call init_starting_star_data(s_aux, ierr)
1716 0 : if (failed('init_starting_star_data',ierr)) return
1717 0 : call star_set_kap_and_eos_handles(id_aux, ierr)
1718 0 : if (failed('set_star_kap_and_eos_handles',ierr)) return
1719 0 : call star_set_colors_handles(id_aux, ierr)
1720 0 : if (failed('star_set_colors_handles',ierr)) return
1721 0 : call store_controls(s_aux, ierr)
1722 0 : if (failed('store_controls',ierr)) return
1723 0 : call do_star_job_controls_before(id_aux, s_aux, .false., ierr)
1724 0 : if (ierr /= 0) return
1725 0 : call star_read_model(id_aux, s% job% saved_model_for_merger_2, ierr)
1726 0 : if (failed('star_read_model',ierr)) return
1727 :
1728 : ! create composition and q array through an entropy sorting
1729 0 : total_mass = s% mstar + s_aux% mstar
1730 0 : partial_mass = 0
1731 0 : i = 1
1732 0 : j = 1
1733 0 : allocate(xq(s% nz + s_aux% nz), xa(s% species, s% nz + s_aux% nz))
1734 0 : do while (i <= s% nz .or. j <= s_aux% nz)
1735 0 : if (j > s_aux% nz .or. (i <= s% nz .and. &
1736 : s% entropy(i) >= s_aux% entropy(j))) then
1737 0 : partial_mass = partial_mass + s% dm(i)
1738 0 : do k=1, s% species
1739 0 : xa(k, i+j-1) = s% xa(k, i)
1740 : end do
1741 0 : i = i + 1
1742 0 : else if (i > s% nz .or. (j <= s_aux% nz .and. &
1743 : s_aux% entropy(j) > s% entropy(i))) then
1744 0 : partial_mass = partial_mass + s_aux% dm(j)
1745 0 : do k=1, s% species
1746 0 : xa(k, i+j-1) = s_aux% xa(k, j)
1747 : end do
1748 0 : j = j + 1
1749 : end if
1750 0 : xq(i+j-2) = partial_mass / total_mass
1751 : !write(*,*) "check", i+j-2, xq(i+j-2), xa(1, i+j-2), xa(2, i+j-2), xa(3, i+j-2)
1752 : end do
1753 : ! Relax composition first, then composition mass
1754 : ! Turn off rotation for relaxation
1755 0 : call star_set_rotation_flag(id, .false., ierr)
1756 0 : if (failed('star_set_rotation_flag',ierr)) then
1757 0 : deallocate(xq,xa)
1758 0 : return
1759 : end if
1760 0 : write(*,*) "Relaxing composition to merger composition"
1761 : call star_relax_composition( &
1762 0 : id, s% job% num_steps_to_relax_composition, s% nz + s_aux% nz, s% species, xa, xq, ierr)
1763 0 : if (failed('star_relax_composition',ierr)) then
1764 0 : deallocate(xq,xa)
1765 0 : return
1766 : end if
1767 0 : write(*,*) "Relaxing star mass to total merger mass"
1768 : call star_relax_mass_scale( &
1769 : id, total_mass/Msun, s% job% dlgm_per_step, &
1770 0 : s% job% change_mass_years_for_dt, ierr)
1771 0 : deallocate(xq,xa)
1772 0 : if (failed('star_relax_mass_scale',ierr)) return
1773 :
1774 0 : end subroutine create_merger_model
1775 :
1776 :
1777 11 : subroutine extend_net(s, ierr)
1778 : use net_def
1779 : use chem_def
1780 : type (star_info), pointer :: s
1781 : integer, intent(out) :: ierr
1782 : real(dp), parameter :: tiny = 1d-10, small = 1d-2
1783 :
1784 : real(dp) :: cntr_h, cntr_he
1785 :
1786 : include 'formats'
1787 :
1788 11 : ierr = 0
1789 :
1790 : !write(*,2) 'extend_net: current net ' // trim(s% net_name), s% model_number
1791 :
1792 11 : if (s% net_name == s% job% adv_net) return
1793 :
1794 11 : if (s% net_name == s% job% co_net) then
1795 0 : if (s% log_max_temperature > 9d0 .or. s% log_center_density > 9d0) then
1796 0 : call change_net(s% job% adv_net)
1797 0 : if (len_trim(s% job% profile_columns_file) > 0) &
1798 0 : write(*,*) 'read ' // trim(s% job% profile_columns_file)
1799 : call star_set_profile_columns( &
1800 0 : s% id, s% job% profile_columns_file, .true., ierr)
1801 : end if
1802 0 : return
1803 : end if
1804 :
1805 11 : if (s% net_name == s% job% h_he_net) then
1806 11 : cntr_h = current_abundance_at_point(s% id, ih1, s% nz, ierr)
1807 : !write(*,2) 'cntr_h', s% model_number, cntr_h, tiny
1808 11 : if (ierr /= 0) return
1809 11 : if (cntr_h > tiny) return
1810 0 : cntr_he = current_abundance_at_point(s% id, ihe4, s% nz, ierr)
1811 : !write(*,2) 'cntr_he', s% model_number, cntr_he, small
1812 0 : if (ierr /= 0) return
1813 0 : if (cntr_he > small) return
1814 0 : if (s% log_max_temperature > 8.3d0 .or. s% log_center_density > 8.5d0) then
1815 0 : call change_net(s% job% co_net)
1816 0 : if (len_trim(s% job% profile_columns_file) > 0) &
1817 0 : write(*,*) 'read ' // trim(s% job% profile_columns_file)
1818 : call star_set_profile_columns( &
1819 0 : s% id, s% job% profile_columns_file, .true., ierr)
1820 : end if
1821 : end if
1822 :
1823 :
1824 : contains
1825 :
1826 :
1827 0 : subroutine change_net(net_name)
1828 : use const_def, only: dp
1829 : character (len=*), intent(in) :: net_name
1830 :
1831 : include 'formats'
1832 :
1833 : call star_change_to_new_net( &
1834 0 : s% id, s% job% adjust_abundances_for_new_isos, net_name, ierr)
1835 0 : if (ierr /= 0) then
1836 0 : write(*,*) 'failed in star_change_to_new_net ' // trim(net_name)
1837 0 : call mesa_error(__FILE__,__LINE__,'change_net')
1838 : return
1839 : end if
1840 :
1841 0 : if (net_name /= s% net_name) then
1842 0 : write(*,*) ' new net_name ', trim(net_name)
1843 0 : write(*,*) 'old s% net_name ', trim(s% net_name)
1844 0 : write(*,*) 'failed to change'
1845 0 : call mesa_error(__FILE__,__LINE__,'change_net')
1846 : end if
1847 :
1848 0 : write(*,'(a)') ' new net = ' // trim(s% net_name)
1849 : !do j=1,s% species
1850 : ! write(*,fmt='(a,x)',advance='no') trim(chem_isos% name(s% chem_id(j)))
1851 : !end do
1852 : !write(*,*)
1853 0 : s% dt_next = s% dt_next/5
1854 : !write(*,1) 'reduce timestep', log10(s% dt_next/secyer)
1855 0 : write(*,'(A)')
1856 : end subroutine change_net
1857 :
1858 :
1859 : end subroutine extend_net
1860 :
1861 :
1862 1 : subroutine before_evolve(id, ierr)
1863 : integer, intent(in) :: id
1864 : integer, intent(out) :: ierr
1865 1 : ierr = 0
1866 : end subroutine before_evolve
1867 :
1868 :
1869 4 : subroutine do_star_job_controls_after(id, s, restart, pgstar_ok, ierr)
1870 : use const_def, only: dp
1871 : use rates_def
1872 : use rates_lib
1873 : use utils_lib, only: utils_OMP_GET_MAX_THREADS
1874 :
1875 : integer, intent(in) :: id
1876 : type (star_info), pointer :: s
1877 : logical, intent(in) :: restart, pgstar_ok
1878 : integer, intent(out) :: ierr
1879 :
1880 : real(dp) :: log_m, log_lifetime, max_dt, max_timestep
1881 : integer :: i, j, nzlo, nzhi, chem_id, chem_id1, chem_id2
1882 : logical :: change_v, change_u
1883 : include 'formats'
1884 :
1885 1 : if (s% job% change_net .or. (s% job% change_initial_net .and. .not. restart)) then
1886 : call star_change_to_new_net( &
1887 0 : id, s% job% adjust_abundances_for_new_isos, s% job% new_net_name, ierr)
1888 0 : if (failed('star_change_to_new_net',ierr)) return
1889 : end if
1890 :
1891 1 : if (s% job% change_small_net .or. &
1892 : (s% job% change_initial_small_net .and. .not. restart)) then
1893 0 : write(*,*) 'change small net to ' // trim(s% job% new_small_net_name)
1894 : call star_change_to_new_small_net( &
1895 0 : id, s% job% adjust_abundances_for_new_isos, s% job% new_small_net_name, ierr)
1896 0 : if (failed('star_change_to_new_small_net',ierr)) return
1897 0 : write(*,*) 'number of species', s% species
1898 : end if
1899 :
1900 :
1901 1 : if (len_trim(s% job% history_columns_file) > 0) &
1902 0 : write(*,*) 'read ' // trim(s% job% history_columns_file)
1903 1 : call star_set_history_columns(id, s% job% history_columns_file, .true., ierr)
1904 1 : if (failed('star_set_history_columns',ierr)) return
1905 :
1906 1 : if (len_trim(s% job% profile_columns_file) > 0) &
1907 0 : write(*,*) 'read ' // trim(s% job% profile_columns_file)
1908 1 : call star_set_profile_columns(id, s% job% profile_columns_file, .true., ierr)
1909 1 : if (failed('star_set_profile_columns',ierr)) return
1910 :
1911 1 : if (pgstar_ok) then
1912 1 : if (s% job% clear_pgstar_history .or. &
1913 : (s% job% clear_initial_pgstar_history .and. .not. restart)) then
1914 1 : call start_new_run_for_pgstar(s, ierr)
1915 1 : if (failed('start_new_run_for_pgstar',ierr)) return
1916 : else
1917 0 : call restart_run_for_pgstar(s, ierr)
1918 0 : if (failed('restart_run_for_pgstar',ierr)) return
1919 : end if
1920 : end if
1921 :
1922 1 : if (s% job% set_tau_factor .or. &
1923 : (s% job% set_initial_tau_factor .and. .not. restart)) then
1924 0 : write(*,1) 'set_tau_factor', s% job% set_to_this_tau_factor
1925 0 : s% tau_factor = s% job% set_to_this_tau_factor
1926 : end if
1927 :
1928 1 : if (s% job% set_initial_age .and. .not. restart) then
1929 0 : write(*,1) 'set_initial_age', s% job% initial_age ! in years
1930 0 : call star_set_age(id, s% job% initial_age, ierr)
1931 0 : if (failed('star_set_age',ierr)) return
1932 : end if
1933 :
1934 1 : if (s% job% set_initial_dt .and. .not. restart) then
1935 0 : if (s% job% years_for_initial_dt > 0d0) then
1936 0 : write(*,1) 'set_initial_dt (years)', s% job% years_for_initial_dt
1937 0 : s% dt_next = s% job% years_for_initial_dt*secyer
1938 0 : else if (s% job% seconds_for_initial_dt > 0d0) then
1939 0 : write(*,1) 'set_initial_dt (seconds)', s% job% seconds_for_initial_dt
1940 0 : s% dt_next = s% job% seconds_for_initial_dt
1941 : end if
1942 : end if
1943 :
1944 1 : if (s% job% limit_initial_dt .and. .not. restart) then
1945 0 : if (s% job% years_for_initial_dt > 0d0) then
1946 0 : write(*,1) 'limit_initial_dt (years)', s% job% years_for_initial_dt
1947 0 : s% dt_next = min(s% dt_next, s% job% years_for_initial_dt*secyer)
1948 0 : else if (s% job% seconds_for_initial_dt > 0d0) then
1949 0 : write(*,1) 'limit_initial_dt (seconds)', s% job% seconds_for_initial_dt
1950 0 : s% dt_next = min(s% dt_next, s% job% seconds_for_initial_dt)
1951 : end if
1952 :
1953 : end if
1954 :
1955 : ! enforce max_timestep on first step
1956 :
1957 1 : if (s% max_years_for_timestep > 0) then
1958 0 : max_timestep = secyer*s% max_years_for_timestep
1959 0 : if (s% max_timestep > 0 .and. s% max_timestep < max_timestep) &
1960 0 : max_timestep = s% max_timestep
1961 : else
1962 1 : max_timestep = s% max_timestep
1963 : end if
1964 :
1965 1 : if (max_timestep > 0 .and. max_timestep < s% dt_next) then
1966 0 : write(*,1) 'max_timestep (seconds)', max_timestep
1967 0 : s% dt_next = max_timestep
1968 : end if
1969 :
1970 1 : if (s% job% set_initial_model_number .and. .not. restart) then
1971 0 : write(*,2) 'set_initial_model_number', s% job% initial_model_number
1972 0 : s% model_number = s% job% initial_model_number
1973 0 : s% init_model_number = s% model_number
1974 : end if
1975 :
1976 1 : if (s% job% set_initial_number_retries .and. .not. restart) then
1977 1 : write(*,2) 'set_initial_number_retries', s% job% initial_number_retries
1978 1 : s% num_retries = s% job% initial_number_retries
1979 : end if
1980 :
1981 1 : if (s% job% steps_to_take_before_terminate >= 0) then
1982 0 : s% max_model_number = s% model_number + s% job% steps_to_take_before_terminate
1983 0 : write(*,2) 'steps_to_take_before_terminate', &
1984 0 : s% job% steps_to_take_before_terminate
1985 0 : write(*,2) 'max_model_number', s% max_model_number
1986 : end if
1987 :
1988 1 : if (s% job% steps_before_start_timing > 0) then
1989 0 : s% job% first_model_for_timing = s% model_number + s% job% steps_before_start_timing
1990 0 : write(*,2) 'steps_before_start_timing', &
1991 0 : s% job% steps_before_start_timing
1992 : end if
1993 :
1994 1 : if (abs(s% job% T9_weaklib_full_off - T9_weaklib_full_off) > 1d-6) then
1995 0 : write(*,1) 'set T9_weaklib_full_off', s% job% T9_weaklib_full_off
1996 0 : T9_weaklib_full_off = s% job% T9_weaklib_full_off
1997 : end if
1998 :
1999 1 : if (abs(s% job% T9_weaklib_full_on - T9_weaklib_full_on) > 1d-6) then
2000 0 : write(*,1) 'set T9_weaklib_full_on', s% job% T9_weaklib_full_on
2001 0 : T9_weaklib_full_on = s% job% T9_weaklib_full_on
2002 : end if
2003 :
2004 1 : if (s% job% weaklib_blend_hi_Z /= weaklib_blend_hi_Z) then
2005 0 : write(*,1) 'set weaklib_blend_hi_Z', s% job% weaklib_blend_hi_Z
2006 0 : weaklib_blend_hi_Z = s% job% weaklib_blend_hi_Z
2007 : end if
2008 :
2009 1 : if (abs(s% job% T9_weaklib_full_off_hi_Z - T9_weaklib_full_off_hi_Z) > 1d-6) then
2010 0 : write(*,1) 'set T9_weaklib_full_off_hi_Z', s% job% T9_weaklib_full_off_hi_Z
2011 0 : T9_weaklib_full_off_hi_Z = s% job% T9_weaklib_full_off_hi_Z
2012 : end if
2013 :
2014 1 : if (abs(s% job% T9_weaklib_full_on_hi_Z - T9_weaklib_full_on_hi_Z) > 1d-6) then
2015 0 : write(*,1) 'set T9_weaklib_full_on_hi_Z', s% job% T9_weaklib_full_on_hi_Z
2016 0 : T9_weaklib_full_on_hi_Z = s% job% T9_weaklib_full_on_hi_Z
2017 : end if
2018 :
2019 : ! set up coulomb corrections for the special weak rates
2020 1 : which_mui_coulomb = get_mui_value(s% job% ion_coulomb_corrections)
2021 1 : which_vs_coulomb = get_vs_value(s% job% electron_coulomb_corrections)
2022 :
2023 : change_v = s% job% change_v_flag .or. &
2024 1 : (s% job% change_initial_v_flag .and. .not. restart)
2025 : change_u = s% job% change_u_flag .or. &
2026 1 : (s% job% change_initial_u_flag .and. .not. restart)
2027 1 : if (change_v .or. change_u) then
2028 : ! do add new before remove old so can set initial values
2029 0 : if (change_v .and. s% job% new_v_flag) then
2030 0 : write(*,*) 'new_v_flag', s% job% new_v_flag
2031 0 : call star_set_v_flag(id, s% job% new_v_flag, ierr)
2032 0 : if (failed('star_set_v_flag',ierr)) return
2033 : end if
2034 0 : if (change_u .and. s% job% new_u_flag) then
2035 0 : write(*,*) 'new_u_flag', s% job% new_u_flag
2036 0 : call star_set_u_flag(id, s% job% new_u_flag, ierr)
2037 0 : if (failed('star_set_u_flag',ierr)) return
2038 : end if
2039 0 : if (change_v .and. .not. s% job% new_v_flag) then
2040 0 : write(*,*) 'new_v_flag', s% job% new_v_flag
2041 0 : call star_set_v_flag(id, s% job% new_v_flag, ierr)
2042 0 : if (failed('star_set_v_flag',ierr)) return
2043 : end if
2044 0 : if (change_u .and. .not. s% job% new_u_flag) then
2045 0 : write(*,*) 'new_u_flag', s% job% new_u_flag
2046 0 : call star_set_u_flag(id, s% job% new_u_flag, ierr)
2047 0 : if (failed('star_set_u_flag',ierr)) return
2048 : end if
2049 : end if
2050 :
2051 1 : if (s% job% change_RTI_flag .or. &
2052 : (s% job% change_initial_RTI_flag .and. .not. restart)) then
2053 0 : write(*,*) 'new_RTI_flag', s% job% new_RTI_flag
2054 0 : call star_set_RTI_flag(id, s% job% new_RTI_flag, ierr)
2055 0 : if (failed('star_set_RTI_flag',ierr)) return
2056 : end if
2057 :
2058 1 : if (s% job% change_RSP2_flag .or. &
2059 : (s% job% change_initial_RSP2_flag .and. .not. restart)) then
2060 0 : write(*,*) 'new_RSP2_flag', s% job% new_RSP2_flag
2061 0 : call star_set_RSP2_flag(id, s% job% new_RSP2_flag, ierr)
2062 0 : if (failed('star_set_RSP2_flag',ierr)) return
2063 : end if
2064 :
2065 1 : if (s% job% change_RSP_flag .or. &
2066 : (s% job% change_initial_RSP_flag .and. .not. restart)) then
2067 0 : write(*,*) 'new_RSP_flag', s% job% new_RSP_flag
2068 0 : call star_set_RSP_flag(id, s% job% new_RSP_flag, ierr)
2069 0 : if (failed('star_set_RSP_flag',ierr)) return
2070 : end if
2071 :
2072 1 : if (s% job% change_w_div_wc_flag .or. &
2073 : (s% job% change_initial_w_div_wc_flag .and. .not. restart)) then
2074 0 : write(*,*) 'new_w_div_wc_flag', s% job% new_w_div_wc_flag
2075 0 : call star_set_w_div_wc_flag(id, s% job% new_w_div_wc_flag, ierr)
2076 0 : if (failed('star_set_w_div_wc_flag',ierr)) return
2077 : end if
2078 :
2079 1 : if (s% job% change_j_rot_flag .or. &
2080 : (s% job% change_initial_j_rot_flag .and. .not. restart)) then
2081 0 : write(*,*) 'new_j_rot_flag', s% job% new_j_rot_flag
2082 0 : call star_set_j_rot_flag(id, s% job% new_j_rot_flag, ierr)
2083 0 : if (failed('star_set_j_rot_flag',ierr)) return
2084 : end if
2085 :
2086 1 : if (s% job% change_D_omega_flag .or. &
2087 : (s% job% change_initial_D_omega_flag .and. .not. restart)) then
2088 0 : call star_set_D_omega_flag(id, s% job% new_D_omega_flag, ierr)
2089 0 : if (failed('star_set_D_omega_flag',ierr)) return
2090 : end if
2091 :
2092 1 : if (s% job% change_am_nu_rot_flag .or. &
2093 : (s% job% change_initial_am_nu_rot_flag .and. .not. restart)) then
2094 0 : call star_set_am_nu_rot_flag(id, s% job% new_am_nu_rot_flag, ierr)
2095 0 : if (failed('star_set_am_nu_rot_flag',ierr)) return
2096 : end if
2097 :
2098 1 : if (s% job% change_rotation_flag .or. &
2099 : (s% job% change_initial_rotation_flag .and. .not. restart)) then
2100 0 : write(*,*) 'new_rotation_flag', s% job% new_rotation_flag
2101 0 : call star_set_rotation_flag(id, s% job% new_rotation_flag, ierr)
2102 0 : if (failed('star_set_rotation_flag',ierr)) return
2103 : end if
2104 :
2105 1 : if (s% rotation_flag .and. s% job% set_omega) then
2106 0 : write(*,1) 'new_omega', s% job% new_omega
2107 0 : call star_set_uniform_omega(id, s% job% new_omega, ierr)
2108 0 : if (failed('star_set_uniform_omega',ierr)) return
2109 : end if
2110 :
2111 1 : if (s% rotation_flag .and. s% job% set_initial_omega .and. .not. restart) then
2112 0 : write(*,1) 'new_omega', s% job% new_omega
2113 0 : call star_set_uniform_omega(id, s% job% new_omega, ierr)
2114 0 : if (failed('star_set_uniform_omega',ierr)) return
2115 : end if
2116 :
2117 1 : if (s% rotation_flag .and. s% job% set_surface_rotation_v) then
2118 0 : s% job% new_omega = s% job% new_surface_rotation_v*1d5/s% r(1)
2119 0 : write(*,1) 'new_surface_rotation_v', &
2120 0 : s% job% new_surface_rotation_v, s% job% new_omega
2121 0 : call star_set_uniform_omega(id, s% job% new_omega, ierr)
2122 0 : if (failed('star_set_uniform_omega',ierr)) return
2123 : end if
2124 :
2125 : if (s% rotation_flag .and. &
2126 1 : s% job% set_initial_surface_rotation_v .and. .not. restart) then
2127 0 : s% job% new_omega = s% job% new_surface_rotation_v*1d5/s% r(1)
2128 0 : write(*,2) 'new_surface_rotation_v', &
2129 0 : s% model_number, s% job% new_surface_rotation_v, s% job% new_omega
2130 0 : call star_set_uniform_omega(id, s% job% new_omega, ierr)
2131 0 : if (failed('star_set_uniform_omega',ierr)) return
2132 : end if
2133 :
2134 1 : if (s% rotation_flag .and. s% job% set_omega_div_omega_crit) then
2135 : s% job% new_omega = &
2136 0 : s% job% new_omega_div_omega_crit*star_surface_omega_crit(id, ierr)
2137 0 : if (failed('star_surface_omega_crit',ierr)) return
2138 0 : write(*,2) 'new_omega_div_omega_crit', &
2139 0 : s% model_number, s% job% new_omega_div_omega_crit, s% job% new_omega
2140 0 : call star_set_uniform_omega(id, s% job% new_omega, ierr)
2141 0 : if (failed('star_set_uniform_omega',ierr)) return
2142 : end if
2143 :
2144 : if (s% rotation_flag .and. &
2145 1 : s% job% set_initial_omega_div_omega_crit .and. .not. restart) then
2146 : s% job% new_omega = &
2147 0 : s% job% new_omega_div_omega_crit*star_surface_omega_crit(id, ierr)
2148 0 : if (failed('star_surface_omega_crit',ierr)) return
2149 0 : write(*,2) 'new_omega_div_omega_crit', &
2150 0 : s% model_number, s% job% new_omega_div_omega_crit, s% job% new_omega
2151 0 : call star_set_uniform_omega(id, s% job% new_omega, ierr)
2152 0 : if (failed('star_set_uniform_omega',ierr)) return
2153 : end if
2154 :
2155 1 : if (s% job% set_to_xa_for_accretion .or. &
2156 : (s% job% set_initial_to_xa_for_accretion .and. .not. restart)) then
2157 0 : write(*,*) 'set_to_xa_for_accretion'
2158 0 : call change_to_xa_for_accretion(id, s% job% set_nzlo, s% job% set_nzhi, ierr)
2159 0 : if (failed('set_to_xa_for_accretion',ierr)) return
2160 : end if
2161 :
2162 1 : if (s% job% first_model_for_timing > 0) &
2163 0 : write(*,2) 'first_model_for_timing', s% job% first_model_for_timing
2164 :
2165 1 : if (s% job% set_uniform_initial_composition .and. .not. restart) then
2166 0 : write(*,'(A)')
2167 0 : write(*,1) 'set_uniform_initial_composition'
2168 0 : write(*,1) 'initial_h1', s% job% initial_h1
2169 0 : write(*,1) 'initial_h2', s% job% initial_h2
2170 0 : write(*,1) 'initial_he3', s% job% initial_he3
2171 0 : write(*,1) 'initial_he4', s% job% initial_he4
2172 0 : select case(s% job% initial_zfracs)
2173 : case (AG89_zfracs)
2174 0 : write(*,1) 'metals AG89'
2175 : case (GN93_zfracs)
2176 0 : write(*,1) 'metals GN93'
2177 : case (GS98_zfracs)
2178 0 : write(*,1) 'metals GS98'
2179 : case (L03_zfracs)
2180 0 : write(*,1) 'metals L03'
2181 : case (AGS05_zfracs)
2182 0 : write(*,1) 'metals AGS05'
2183 : case (AGSS09_zfracs)
2184 0 : write(*,1) 'metals AGSS09'
2185 : case (L09_zfracs)
2186 0 : write(*,1) 'metals L09'
2187 : case (A09_Prz_zfracs)
2188 0 : write(*,1) 'metals A09_Prz'
2189 : case (MB22_photospheric_zfracs)
2190 0 : write(*,1) 'metals MB22_photospheric'
2191 : case (AAG21_photospheric_zfracs)
2192 0 : write(*,1) 'metals AAG21_photospheric'
2193 : case default
2194 0 : write(*,2) 'unknown value for initial_zfracs', s% job% initial_zfracs
2195 : end select
2196 : call star_set_standard_composition( &
2197 : id, s% job% initial_h1, s% job% initial_h2, &
2198 : s% job% initial_he3, s% job% initial_he4, s% job% initial_zfracs, &
2199 0 : s% job% dump_missing_metals_into_heaviest, ierr)
2200 0 : if (failed('set_uniform_initial_composition',ierr)) return
2201 : end if
2202 :
2203 1 : if (s% job% relax_initial_composition .and. .not. restart) then
2204 0 : call do_relax_initial_composition(ierr)
2205 0 : if (failed('do_relax_initial_composition',ierr)) return
2206 : end if
2207 :
2208 1 : if (s% job% relax_initial_to_xaccrete .and. .not. restart) then
2209 0 : call star_relax_to_xaccrete(id, s% job% num_steps_to_relax_composition, ierr)
2210 0 : if (failed('star_relax_to_xaccrete',ierr)) return
2211 : end if
2212 :
2213 1 : if (s% job% set_uniform_xa_from_file) then
2214 0 : call star_uniform_xa_from_file(id, s% job% file_for_uniform_xa, ierr)
2215 0 : if (failed('star_uniform_xa_from_file',ierr)) return
2216 : end if
2217 :
2218 1 : if (s% job% relax_initial_angular_momentum .and. .not. restart) then
2219 0 : call do_relax_initial_angular_momentum(ierr)
2220 0 : if (failed('do_relax_initial_angular_momentum',ierr)) return
2221 : end if
2222 :
2223 1 : if (s% job% relax_initial_entropy .and. .not. restart) then
2224 0 : call do_relax_initial_entropy(ierr)
2225 0 : if (failed('do_relax_initial_entropy',ierr)) return
2226 : end if
2227 :
2228 1 : if (s% job% mix_section .or. &
2229 : (s% job% mix_initial_section .and. .not. restart)) then
2230 0 : write(*,*) 'mix_section'
2231 : call uniform_mix_section( &
2232 0 : id, s% job% mix_section_nzlo, s% job% mix_section_nzhi, ierr)
2233 0 : if (failed('uniform_mix_section',ierr)) return
2234 : end if
2235 :
2236 1 : if (s% job% mix_initial_envelope_down_to_T > 0d0 .and. .not. restart) then
2237 0 : call uniform_mix_envelope_down_to_T(id, s% job% mix_initial_envelope_down_to_T, ierr)
2238 0 : if (failed('uniform_mix_envelope_down_to_T',ierr)) return
2239 : end if
2240 :
2241 1 : if (s% job% mix_envelope_down_to_T > 0d0) then
2242 0 : call uniform_mix_envelope_down_to_T(id, s% job% mix_envelope_down_to_T, ierr)
2243 0 : if (failed('uniform_mix_envelope_down_to_T',ierr)) return
2244 : end if
2245 :
2246 1 : if (s% job% mix_initial_envelope_down_to_T > 0d0) then
2247 0 : call uniform_mix_envelope_down_to_T(id, s% job% mix_initial_envelope_down_to_T, ierr)
2248 0 : if (failed('uniform_mix_envelope_down_to_T',ierr)) return
2249 : end if
2250 :
2251 1 : if (s% job% set_uniform_initial_xa_from_file .and. .not. restart) then
2252 0 : call star_uniform_xa_from_file(id, s% job% file_for_uniform_xa, ierr)
2253 0 : if (failed('star_uniform_xa_from_file',ierr)) return
2254 : end if
2255 :
2256 : ! do change Z before change Y since changing Z can change Y
2257 1 : if (s% job% change_Z) then
2258 0 : call star_set_z(id, s% job% new_Z, ierr)
2259 0 : if (failed('star_set_z',ierr)) return
2260 : end if
2261 :
2262 1 : if (s% job% change_initial_Z .and. .not. restart) then
2263 0 : call star_set_z(id, s% job% new_Z, ierr)
2264 0 : if (failed('star_set_z',ierr)) return
2265 : end if
2266 :
2267 1 : if (s% job% change_Y) then
2268 0 : call star_set_y(id, s% job% new_Y, ierr)
2269 0 : if (failed('change_Y',ierr)) return
2270 : end if
2271 :
2272 1 : if (s% job% change_initial_Y .and. .not. restart) then
2273 0 : call star_set_y(id, s% job% new_Y, ierr)
2274 0 : if (failed('change_initial_Y',ierr)) return
2275 : end if
2276 :
2277 1 : if (s% job% zero_alpha_RTI .or. &
2278 : (s% job% zero_initial_alpha_RTI .and. .not. restart)) then
2279 0 : call star_zero_alpha_RTI(id, ierr)
2280 0 : if (failed('star_zero_alpha_RTI',ierr)) return
2281 : end if
2282 :
2283 1 : if (s% job% set_abundance .or. &
2284 : (s% job% set_initial_abundance .and. .not. restart)) then
2285 0 : nzlo = s% job% set_abundance_nzlo
2286 0 : nzhi = s% job% set_abundance_nzhi
2287 0 : if (nzhi <= 0) nzhi = s% nz
2288 0 : if (nzlo <= 0) nzlo = 1
2289 0 : write(*, *) 'set_abundance of ', &
2290 0 : trim(s% job% chem_name), s% job% new_frac, nzlo, nzhi
2291 0 : chem_id = get_nuclide_index(s% job% chem_name)
2292 0 : if (chem_id <= 0) then
2293 0 : write(*,*) 'failed to find ' // trim(s% job% chem_name)
2294 0 : write(*,*) 'check valid chem_isos% names in chem/public/chem_def.f'
2295 : end if
2296 0 : call set_abundance_in_section(id, chem_id, s% job% new_frac, nzlo, nzhi, ierr)
2297 0 : if (failed('set_abundance_in_section',ierr)) return
2298 : end if
2299 :
2300 1 : if (s% job% replace_element .or. &
2301 : (s% job% replace_initial_element .and. .not. restart)) then
2302 0 : write(*, *) 'replace_element ', &
2303 0 : trim(s% job% chem_name1), ' by ', trim(s% job% chem_name2)
2304 0 : chem_id1 = get_nuclide_index(s% job% chem_name1)
2305 0 : chem_id2 = get_nuclide_index(s% job% chem_name2)
2306 0 : if (chem_id1 <= 0) then
2307 0 : write(*,*) 'failed to find ' // trim(s% job% chem_name1)
2308 0 : write(*,*) 'check valid chem_isos% names in chem/public/chem_def.f'
2309 : end if
2310 0 : if (chem_id2 <= 0) then
2311 0 : write(*,*) 'failed to find ' // trim(s% job% chem_name2)
2312 0 : write(*,*) 'check valid chem_isos% names in chem/public/chem_def.f'
2313 : end if
2314 0 : nzhi = s% job% replace_element_nzhi
2315 0 : nzlo = s% job% replace_element_nzlo
2316 0 : if (nzhi <= 0) nzhi = s% nz
2317 0 : if (nzlo <= 0) nzlo = 1
2318 0 : write(*, *) 'in section', nzlo, nzhi
2319 : call replace_element_in_section( &
2320 0 : id, chem_id1, chem_id2, nzlo, nzhi, ierr)
2321 0 : if (failed('replace_element_in_section',ierr)) return
2322 : end if
2323 :
2324 1 : if (s% job% set_irradiation .or. &
2325 : (s% job% set_initial_irradiation .and. .not. restart)) then
2326 0 : write(*,2) 'set_irradiation'
2327 0 : s% irradiation_flux = s% job% set_to_this_irrad_flux
2328 0 : s% column_depth_for_irradiation = s% job% irrad_col_depth
2329 : end if
2330 :
2331 1 : if (s% job% do_special_test) then
2332 0 : write(*, *) 'do_special_test'
2333 0 : call star_special_test(id, ierr)
2334 0 : if (failed('star_special_test',ierr)) return
2335 : end if
2336 :
2337 1 : if (s% job% set_v_center .or. &
2338 : (s% job% set_initial_v_center .and. .not. restart)) then
2339 0 : write(*, 1) 'set_v_center', s% job% new_v_center
2340 0 : s% v_center = s% job% new_v_center
2341 : end if
2342 :
2343 1 : if (s% job% set_L_center .or. &
2344 : (s% job% set_initial_L_center .and. .not. restart)) then
2345 0 : write(*, 1) 'set_L_center', s% job% new_L_center
2346 0 : s% L_center = s% job% new_L_center*Lsun
2347 : end if
2348 :
2349 : ! do "set" before "relax"
2350 :
2351 : ! must do relax Z before relax Y since relax Z can change Y
2352 : ! (Warrick Ball pointed out this requirement)
2353 1 : if (s% job% relax_initial_Z .and. .not. restart) then
2354 0 : write(*,1) 'relax_initial_Z', s% job% new_Z
2355 : call star_relax_Z(id, s% job% new_Z, s% relax_dlnZ, &
2356 0 : s% job% relax_Z_minq, s% job% relax_Z_maxq, ierr)
2357 0 : if (failed('star_relax_Z',ierr)) return
2358 0 : write(*, 1) 'new z', get_current_z(id, ierr)
2359 0 : if (failed('get_current_z',ierr)) return
2360 : end if
2361 :
2362 1 : if (s% job% relax_Z) then
2363 0 : write(*,1) 'relax_Z', s% job% new_Z
2364 : call star_relax_Z(id, s% job% new_Z, s% relax_dlnZ, &
2365 0 : s% job% relax_Z_minq, s% job% relax_Z_maxq, ierr)
2366 0 : if (failed('star_relax_Z',ierr)) return
2367 0 : write(*, 1) 'new z', get_current_z(id, ierr)
2368 0 : if (failed('get_current_z',ierr)) return
2369 : end if
2370 :
2371 1 : if (s% job% relax_initial_Y .and. .not. restart) then
2372 0 : write(*,1) 'relax_initial_Y', s% job% new_Y
2373 : call star_relax_Y(id, s% job% new_Y, s% relax_dY, &
2374 0 : s% job% relax_Y_minq, s% job% relax_Y_maxq, ierr)
2375 0 : if (failed('star_relax_Y',ierr)) return
2376 0 : write(*, 1) 'new y', get_current_y(id, ierr)
2377 0 : if (failed('get_current_y',ierr)) return
2378 : end if
2379 :
2380 1 : if (s% job% relax_Y) then
2381 0 : write(*,1) 'relax_Y', s% job% new_Y
2382 : call star_relax_Y(id, s% job% new_Y, s% relax_dY, &
2383 0 : s% job% relax_Y_minq, s% job% relax_Y_maxq, ierr)
2384 0 : if (failed('star_relax_Y',ierr)) return
2385 0 : write(*, 1) 'new y', get_current_y(id, ierr)
2386 0 : if (failed('get_current_y',ierr)) return
2387 : end if
2388 :
2389 1 : if (s% job% relax_mass) then
2390 0 : write(*, 1) 'relax_mass', s% job% new_mass
2391 0 : call star_relax_mass(id, s% job% new_mass, s% job% lg_max_abs_mdot, ierr)
2392 0 : if (failed('star_relax_mass',ierr)) return
2393 : end if
2394 :
2395 1 : if (s% job% relax_mass_to_remove_H_env) then
2396 0 : write(*, 1) 'relax_mass_to_remove_H_env_mass'
2397 : call star_relax_mass_to_remove_H_env( &
2398 0 : id, s% job% extra_mass_retained_by_remove_H_env, s% job% lg_max_abs_mdot, ierr)
2399 0 : if (failed('star_relax_mass_to_remove_H_env',ierr)) return
2400 : end if
2401 :
2402 1 : if (s% job% relax_dxdt_nuc_factor .or. &
2403 : (s% job% relax_initial_dxdt_nuc_factor .and. .not. restart)) then
2404 0 : write(*, 1) 'relax_dxdt_nuc_factor', s% job% new_dxdt_nuc_factor
2405 : call star_relax_dxdt_nuc_factor( &
2406 0 : id, s% job% new_dxdt_nuc_factor, s% job% dxdt_nuc_factor_multiplier, ierr)
2407 0 : if (failed('star_relax_dxdt_nuc_factor',ierr)) return
2408 : end if
2409 :
2410 1 : if (s% job% relax_eps_nuc_factor .or. &
2411 : (s% job% relax_initial_eps_nuc_factor .and. .not. restart)) then
2412 0 : write(*, 1) 'relax_eps_nuc_factor', s% job% new_eps_nuc_factor
2413 : call star_relax_eps_nuc_factor( &
2414 0 : id, s% job% new_eps_nuc_factor, s% job% eps_nuc_factor_multiplier, ierr)
2415 0 : if (failed('star_relax_eps_nuc_factor',ierr)) return
2416 : end if
2417 :
2418 1 : if (s% job% relax_opacity_max .or. &
2419 : (s% job% relax_initial_opacity_max .and. .not. restart)) then
2420 0 : write(*, 1) 'relax_opacity_max', s% job% new_opacity_max
2421 : call star_relax_opacity_max( &
2422 0 : id, s% job% new_opacity_max, s% job% opacity_max_multiplier, ierr)
2423 0 : if (failed('star_relax_opacity_max',ierr)) return
2424 : end if
2425 :
2426 1 : if (s% job% relax_max_surf_dq .or. &
2427 : (s% job% relax_initial_max_surf_dq .and. .not. restart)) then
2428 0 : write(*, 1) 'relax_max_surf_dq', s% job% new_max_surf_dq
2429 : call star_relax_max_surf_dq( &
2430 0 : id, s% job% new_max_surf_dq, s% job% max_surf_dq_multiplier, ierr)
2431 0 : if (failed('star_relax_max_surf_dq',ierr)) return
2432 : end if
2433 :
2434 1 : if (s% job% relax_initial_mass .and. .not. restart) then
2435 0 : write(*, 1) 'relax_initial_mass to new_mass', s% job% new_mass
2436 0 : call star_relax_mass(id, s% job% new_mass, s% job% lg_max_abs_mdot, ierr)
2437 0 : if (failed('relax_initial_mass',ierr)) return
2438 : end if
2439 :
2440 1 : if (s% job% relax_initial_mass_to_remove_H_env .and. .not. restart) then
2441 0 : write(*, 1) 'relax_initial_mass_to_remove_H_env'
2442 : call star_relax_mass_to_remove_H_env( &
2443 0 : id, s% job% extra_mass_retained_by_remove_H_env, s% job% lg_max_abs_mdot, ierr)
2444 0 : if (failed('relax_initial_mass_to_remove_H_env',ierr)) return
2445 : end if
2446 :
2447 1 : if (s% job% relax_mass_scale .or. &
2448 : (s% job% relax_initial_mass_scale .and. .not. restart)) then
2449 0 : write(*, 1) 'relax_mass_scale', s% job% new_mass
2450 : call star_relax_mass_scale( &
2451 : id, s% job% new_mass, s% job% dlgm_per_step, &
2452 0 : s% job% change_mass_years_for_dt, ierr)
2453 0 : if (failed('star_relax_mass_scale',ierr)) return
2454 : end if
2455 :
2456 1 : if (s% job% relax_core .or. &
2457 : (s% job% relax_initial_core .and. .not. restart)) then
2458 0 : write(*, 1) 'relax_core', s% job% new_core_mass
2459 : call star_relax_core( &
2460 : id, s% job% new_core_mass, s% job% dlg_core_mass_per_step, &
2461 : s% job% relax_core_years_for_dt, &
2462 0 : s% job% core_avg_rho, s% job% core_avg_eps, ierr)
2463 0 : if (failed('star_relax_core',ierr)) return
2464 : end if
2465 :
2466 1 : call do_remove_center(id, s, restart, ierr)
2467 1 : if (ierr /= 0) return
2468 :
2469 1 : if (s% job% relax_M_center .or. &
2470 : (s% job% relax_initial_M_center .and. .not. restart)) then
2471 0 : write(*, 1) 'relax_M_center', s% job% new_mass
2472 : call star_relax_M_center( &
2473 0 : id, s% job% new_mass, s% job% dlgm_per_step, s% job% relax_M_center_dt, ierr)
2474 0 : if (failed('star_relax_M_center',ierr)) return
2475 : end if
2476 :
2477 1 : if (s% job% relax_R_center .or. &
2478 : (s% job% relax_initial_R_center .and. .not. restart)) then
2479 0 : write(*, 1) 'relax_R_center', s% job% new_R_center
2480 : call star_relax_R_center( &
2481 0 : id, s% job% new_R_center, s% job% dlgR_per_step, s% job% relax_R_center_dt, ierr)
2482 0 : if (failed('star_relax_R_center',ierr)) return
2483 : end if
2484 :
2485 1 : if (s% job% relax_v_center .or. &
2486 : (s% job% relax_initial_v_center .and. .not. restart)) then
2487 0 : write(*, 1) 'relax_v_center', s% job% new_v_center
2488 : call star_relax_v_center( &
2489 0 : id, s% job% new_v_center, s% job% dv_per_step, s% job% relax_v_center_dt, ierr)
2490 0 : if (failed('star_relax_v_center',ierr)) return
2491 : end if
2492 :
2493 1 : if (s% job% relax_L_center .or. &
2494 : (s% job% relax_initial_L_center .and. .not. restart)) then
2495 0 : write(*, 1) 'relax_L_center', s% job% new_L_center
2496 : call star_relax_L_center( &
2497 0 : id, s% job% new_L_center, s% job% dlgL_per_step, s% job% relax_L_center_dt, ierr)
2498 0 : if (failed('star_relax_L_center',ierr)) return
2499 : end if
2500 :
2501 1 : if (s% job% relax_tau_factor .or. &
2502 : (s% job% relax_initial_tau_factor .and. .not. restart)) then
2503 0 : write(*,1) 'relax_tau_factor', s% job% relax_to_this_tau_factor
2504 : call star_relax_tau_factor( &
2505 0 : id, s% job% relax_to_this_tau_factor, s% job% dlogtau_factor, ierr)
2506 0 : if (failed('star_relax_tau_factor',ierr)) return
2507 : end if
2508 :
2509 1 : if (s% job% relax_opacity_factor .or. &
2510 : (s% job% relax_initial_opacity_factor .and. .not. restart)) then
2511 0 : write(*,1) 'relax_opacity_factor', s% job% relax_to_this_opacity_factor
2512 : call star_relax_opacity_factor( &
2513 0 : id, s% job% relax_to_this_opacity_factor, s% job% d_opacity_factor, ierr)
2514 0 : if (failed('star_relax_opacity_factor',ierr)) return
2515 : end if
2516 :
2517 1 : if (s% job% relax_irradiation .or. &
2518 : (s% job% relax_initial_irradiation .and. .not. restart)) then
2519 0 : write(*,2) 'relax_irradiation -- min steps', s% job% relax_irradiation_min_steps
2520 0 : write(*,1) 'relax_irradiation -- max yrs dt', s% job% relax_irradiation_max_yrs_dt
2521 : call star_relax_irradiation(id, &
2522 : s% job% relax_irradiation_min_steps, &
2523 : s% job% relax_to_this_irrad_flux, s% job% irrad_col_depth, &
2524 0 : s% job% relax_irradiation_max_yrs_dt, ierr)
2525 0 : if (failed('star_relax_irradiation',ierr)) return
2526 : end if
2527 :
2528 1 : if (s% job% relax_mass_change .or. &
2529 : (s% job% relax_initial_mass_change .and. .not. restart)) then
2530 0 : write(*,2) 'relax_mass_change -- min steps', &
2531 0 : s% job% relax_mass_change_min_steps
2532 0 : write(*,1) 'relax_mass_change -- max yrs dt', &
2533 0 : s% job% relax_mass_change_max_yrs_dt
2534 0 : write(*,1) 'relax_mass_change -- initial_mass_change', &
2535 0 : s% job% relax_mass_change_init_mdot
2536 0 : write(*,1) 'relax_mass_change -- final_mass_change', &
2537 0 : s% job% relax_mass_change_final_mdot
2538 : call star_relax_mass_change(id, &
2539 : s% job% relax_mass_change_min_steps, &
2540 : s% job% relax_mass_change_init_mdot, &
2541 : s% job% relax_mass_change_final_mdot, &
2542 0 : s% job% relax_mass_change_max_yrs_dt, ierr)
2543 0 : if (failed('star_relax_mass_change',ierr)) return
2544 : end if
2545 :
2546 1 : call do_remove_initial_surface(id, s, restart, ierr)
2547 1 : if (ierr /= 0) return
2548 :
2549 1 : call do_remove_surface(id, s, ierr)
2550 1 : if (ierr /= 0) return
2551 :
2552 1 : if (s% rotation_flag .and. s% job% relax_omega) then
2553 0 : write(*,1) 'new_omega', s% job% new_omega
2554 : call star_relax_uniform_omega( &
2555 : id, relax_to_new_omega, &
2556 : s% job% new_omega, s% job% num_steps_to_relax_rotation,&
2557 0 : s% job% relax_omega_max_yrs_dt, ierr)
2558 0 : if (failed('star_relax_uniform_omega',ierr)) return
2559 : end if
2560 :
2561 1 : if (s% rotation_flag .and. s% job% relax_initial_omega .and. .not. restart) then
2562 : call star_relax_uniform_omega( &
2563 : id, relax_to_new_omega, &
2564 : s% job% new_omega, s% job% num_steps_to_relax_rotation,&
2565 0 : s% job% relax_omega_max_yrs_dt, ierr)
2566 0 : if (failed('star_relax_uniform_omega',ierr)) return
2567 0 : write(*,1) 'new_omega', s% job% new_omega
2568 : end if
2569 :
2570 1 : if (s% rotation_flag .and. s% job% relax_omega_div_omega_crit) then
2571 0 : if (failed('star_surface_omega_crit',ierr)) return
2572 : call star_relax_uniform_omega( &
2573 : id, relax_to_new_omega_div_omega_crit, &
2574 : s% job% new_omega_div_omega_crit, &
2575 : s% job% num_steps_to_relax_rotation,&
2576 0 : s% job% relax_omega_max_yrs_dt, ierr)
2577 0 : if (failed('star_relax_uniform_omega',ierr)) return
2578 0 : write(*,2) 'new_omega_div_omega_crit', &
2579 0 : s% model_number, s% job% new_omega_div_omega_crit
2580 : end if
2581 :
2582 : if (s% rotation_flag .and. &
2583 1 : s% job% relax_initial_omega_div_omega_crit .and. .not. restart) then
2584 0 : if (failed('star_surface_omega_crit',ierr)) return
2585 : call star_relax_uniform_omega( &
2586 : id, relax_to_new_omega_div_omega_crit, &
2587 : s% job% new_omega_div_omega_crit, &
2588 : s% job% num_steps_to_relax_rotation,&
2589 0 : s% job% relax_omega_max_yrs_dt, ierr)
2590 0 : if (failed('star_relax_uniform_omega',ierr)) return
2591 0 : write(*,2) 'new_omega_div_omega_crit', &
2592 0 : s% model_number, s% job% new_omega_div_omega_crit
2593 : end if
2594 :
2595 1 : if (s% rotation_flag .and. s% job% relax_surface_rotation_v) then
2596 : call star_relax_uniform_omega( &
2597 : id, relax_to_new_surface_rotation_v, &
2598 : s% job% new_surface_rotation_v, s% job% num_steps_to_relax_rotation,&
2599 0 : s% job% relax_omega_max_yrs_dt, ierr)
2600 0 : if (failed('star_relax_uniform_omega',ierr)) return
2601 0 : s% job% new_omega = s% job% new_surface_rotation_v*1d5/s% r(1)
2602 0 : write(*,1) 'new_surface_rotation_v', &
2603 0 : s% job% new_surface_rotation_v, s% job% new_omega
2604 : end if
2605 :
2606 : if (s% rotation_flag .and. &
2607 1 : s% job% relax_initial_surface_rotation_v .and. .not. restart) then
2608 0 : write(*,1) 'new_omega', s% job% new_omega
2609 0 : write(*,*) 'call star_relax_uniform_omega'
2610 : call star_relax_uniform_omega( &
2611 : id, relax_to_new_surface_rotation_v, &
2612 : s% job% new_surface_rotation_v, s% job% num_steps_to_relax_rotation,&
2613 0 : s% job% relax_omega_max_yrs_dt, ierr)
2614 0 : if (failed('star_relax_uniform_omega',ierr)) return
2615 0 : write(*,2) 'new_surface_rotation_v', &
2616 0 : s% model_number, s% job% new_surface_rotation_v
2617 : end if
2618 :
2619 1 : if (s% job% set_max_dt_to_frac_lifetime) then
2620 0 : log_m = log10(s% star_mass) ! in Msun units
2621 0 : log_lifetime = 9.921d0 - (3.6648d0 + (1.9697d0 - 0.9369d0*log_m)*log_m)*log_m
2622 : ! Iben & Laughlin (1989) as quoted in H&K (eqn 2.3)
2623 0 : max_dt = s% job% max_frac_of_lifetime_per_step*secyer*exp10(log_lifetime)
2624 0 : if (max_dt < s% max_timestep) then
2625 0 : s% max_timestep = max_dt
2626 0 : write(*, *) 'set_max_dt_to_frac_lifetime: lg(maxdt/secyer)', &
2627 0 : log10(s% max_timestep/secyer)
2628 : end if
2629 : end if
2630 :
2631 : ! print out info about selected non-standard parameter settings
2632 :
2633 1 : write(*,*) 'net name ' // trim(s% net_name)
2634 :
2635 1 : if (s% do_element_diffusion) &
2636 0 : write(*,*) 'do_element_diffusion', s% do_element_diffusion
2637 :
2638 1 : if (s% RSP_flag) &
2639 0 : write(*,*) 'RSP_flag', s% RSP_flag
2640 :
2641 1 : if (s% v_flag) &
2642 0 : write(*,*) 'v_flag', s% v_flag
2643 :
2644 1 : if (s% u_flag) &
2645 0 : write(*,*) 'u_flag', s% u_flag
2646 :
2647 1 : if (s% rotation_flag) &
2648 0 : write(*,*) 'rotation_flag', s% rotation_flag
2649 :
2650 1 : if (s% w_div_wc_flag) &
2651 0 : write(*,*) 'w_div_wc_flag', s% w_div_wc_flag
2652 :
2653 1 : if (s% j_rot_flag) &
2654 0 : write(*,*) 'j_rot_flag', s% j_rot_flag
2655 :
2656 1 : if (s% mix_factor /= 1d0) &
2657 0 : write(*,1) 'mix_factor', s% mix_factor
2658 :
2659 1 : if (abs(s% tau_base - 2d0/3d0) > 1d-4) &
2660 0 : write(*,1) 'tau_base', s% tau_base
2661 :
2662 1 : if (abs(s% tau_factor - 1) > 1d-4) &
2663 0 : write(*,1) 'tau_factor', s% tau_factor
2664 :
2665 1 : if (s% eps_grav_factor /= 1) &
2666 0 : write(*,1) 'eps_grav_factor', s% eps_grav_factor
2667 :
2668 1 : if (s% eps_mdot_factor /= 1) &
2669 0 : write(*,1) 'eps_mdot_factor', s% eps_mdot_factor
2670 :
2671 1 : if (s% dxdt_nuc_factor /= 1) &
2672 0 : write(*,1) 'dxdt_nuc_factor', s% dxdt_nuc_factor
2673 :
2674 1 : if (.NOT. ( &
2675 : s% atm_option == 'T_tau' .AND. &
2676 : s% atm_T_tau_relation == 'Eddington' .AND. &
2677 : s% atm_T_tau_opacity == 'fixed')) &
2678 0 : write(*,1) 'atm_option: ' // trim(s% atm_option)
2679 :
2680 1 : if (s% M_center /= 0) then
2681 0 : write(*,1) 'xmstar/mstar', s% xmstar/s% mstar
2682 0 : write(*,1) 'xmstar (g)', s% xmstar
2683 0 : write(*,1) 'M_center (g)', s% M_center
2684 0 : write(*,1) 'xmstar/Msun', s% xmstar/Msun
2685 0 : write(*,1) 'M_center/Msun', s% M_center/Msun
2686 : end if
2687 :
2688 1 : if (s% v_flag .or. s% u_flag) then
2689 0 : if (s% v_center /= 0) &
2690 0 : write(*,1) 'v_center (cm/s)', s% v_center
2691 : end if
2692 :
2693 1 : if (s% R_center /= 0) then
2694 0 : write(*,1) 'R_center (cm)', s% R_center
2695 0 : write(*,1) 'R_center/Rsun', s% R_center/Rsun
2696 0 : write(*,1) 'core density', &
2697 0 : s% M_center/(4*pi/3*s% R_center*s% R_center*s% R_center)
2698 : end if
2699 :
2700 1 : if (s% L_center /= 0) &
2701 0 : write(*,1) 'L_center/Lsun', s% L_center/Lsun
2702 :
2703 1 : if (s% opacity_max > 0) &
2704 0 : write(*,1) 'opacity_max', s% opacity_max
2705 :
2706 1 : if (s% opacity_min > 0) &
2707 0 : write(*,1) 'opacity_min', s% opacity_min
2708 :
2709 1 : if (s% job% show_net_reactions_info) then
2710 0 : write(*,'(a)') ' net reactions '
2711 0 : call show_net_reactions_and_info(s% net_handle, 6, ierr)
2712 0 : if (failed('show_net_reactions_and_info',ierr)) return
2713 : end if
2714 :
2715 1 : if (s% job% list_net_reactions) then
2716 0 : write(*,'(a)') ' net reactions '
2717 0 : call show_net_reactions(s% net_handle, 6, ierr)
2718 0 : if (failed('show_net_reactions',ierr)) return
2719 : end if
2720 :
2721 : if (s% job% set_cumulative_energy_error .or. &
2722 1 : (s% job% set_initial_cumulative_energy_error .and. .not. restart) .or. &
2723 : (s% model_number == s% job% set_cumulative_energy_error_at_step)) then
2724 0 : write(*,1) 'set_cumulative_energy_error', s% job% new_cumulative_energy_error
2725 0 : s% cumulative_energy_error = s% job% new_cumulative_energy_error
2726 : end if
2727 :
2728 1 : if (s% job% show_net_species_info) then
2729 0 : write(*,'(a)') ' species'
2730 0 : do j=1,s% species
2731 0 : write(*,'(i6,3x,a)') j, chem_isos% name(s% chem_id(j))
2732 : end do
2733 0 : write(*,'(A)')
2734 : end if
2735 :
2736 1 : if (s% job% show_eqns_and_vars_names) then
2737 0 : do i=1,s% nvar_total
2738 0 : write(*,*) i, s% nameofvar(i), s% nameofequ(i)
2739 : end do
2740 0 : write(*,'(A)')
2741 : end if
2742 :
2743 1 : write(*,*) 'kap_option ' // trim(kap_option_str(s% kap_rq% kap_option))
2744 1 : write(*,*) 'kap_CO_option ' // trim(kap_CO_option_str(s% kap_rq% kap_CO_option))
2745 1 : write(*,*) 'kap_lowT_option ' // trim(kap_lowT_option_str(s% kap_rq% kap_lowT_option))
2746 1 : write(*,2) 'OMP_NUM_THREADS', utils_omp_get_max_threads()
2747 :
2748 1 : call check_if_want_to_stop_warnings(s)
2749 :
2750 : contains
2751 :
2752 0 : subroutine do_relax_initial_composition(ierr)
2753 : use utils_lib
2754 : integer, intent(out) :: ierr
2755 0 : real(dp), pointer :: xq(:), xa(:,:)
2756 : integer :: num_pts, num_species, i, iounit
2757 : include 'formats'
2758 :
2759 0 : write(*,'(A)')
2760 0 : write(*,1) 'relax_initial_composition'
2761 :
2762 : open(newunit=iounit, file=trim(s% job% relax_composition_filename), &
2763 0 : status='old', action='read', iostat=ierr)
2764 0 : if (ierr /= 0) then
2765 0 : write(*,*) 'open failed', ierr, iounit
2766 0 : write(*, '(a)') 'failed to open ' // trim(s% job% relax_composition_filename)
2767 0 : return
2768 : end if
2769 0 : read(iounit, *, iostat=ierr) num_pts, num_species
2770 0 : if (ierr /= 0) then
2771 0 : close(iounit)
2772 : write(*, '(a)') 'failed while trying to read 1st line of ' // &
2773 0 : trim(s% job% relax_composition_filename)
2774 0 : return
2775 : end if
2776 0 : if(num_species /= s% species) then
2777 0 : write(*,*) 'Error in ',trim(s% job% relax_composition_filename)
2778 0 : write(*,'(a,I4,a)') 'got ',num_species,' species'
2779 0 : write(*,'(a,I4,a)') 'expected ', s% species,' species'
2780 0 : write(*,'(A)')
2781 0 : ierr=-1
2782 0 : return
2783 : end if
2784 0 : allocate(xq(num_pts), xa(num_species,num_pts))
2785 0 : do i = 1, num_pts
2786 0 : read(iounit,*,iostat=ierr) xq(i), xa(1:num_species,i)
2787 0 : if (ierr /= 0) then
2788 0 : close(iounit)
2789 : write(*, '(a)') &
2790 0 : 'failed while trying to read ' // trim(s% job% relax_composition_filename)
2791 0 : write(*,*) 'line', i+1
2792 0 : write(*,*) 'perhaps wrong info in 1st line?'
2793 0 : write(*,*) '1st line must have num_pts and num_species in that order'
2794 0 : deallocate(xq,xa)
2795 0 : return
2796 : end if
2797 : end do
2798 0 : close(iounit)
2799 :
2800 : call star_relax_composition( &
2801 0 : id, s% job% num_steps_to_relax_composition, num_pts, num_species, xa, xq, ierr)
2802 0 : deallocate(xq,xa)
2803 :
2804 0 : end subroutine do_relax_initial_composition
2805 :
2806 0 : subroutine do_relax_initial_angular_momentum(ierr)
2807 : use utils_lib
2808 : integer, intent(out) :: ierr
2809 0 : real(dp), pointer :: xq(:), angular_momentum(:)
2810 : integer :: num_pts, i, iounit
2811 : include 'formats'
2812 :
2813 0 : write(*,'(A)')
2814 0 : write(*,1) 'relax_initial_angular_momentum'
2815 :
2816 : open(newunit=iounit, file=trim(s% job% relax_angular_momentum_filename), &
2817 0 : status='old', action='read', iostat=ierr)
2818 0 : if (ierr /= 0) then
2819 0 : write(*,*) 'open failed', ierr, iounit
2820 0 : write(*, '(a)') 'failed to open "' // trim(s% job% relax_angular_momentum_filename)//'"'
2821 0 : return
2822 : end if
2823 0 : read(iounit, *, iostat=ierr) num_pts
2824 0 : if (ierr /= 0) then
2825 0 : close(iounit)
2826 : write(*, '(a)') 'failed while trying to read 1st line of ' // &
2827 0 : trim(s% job% relax_angular_momentum_filename)
2828 0 : return
2829 : end if
2830 0 : allocate(xq(num_pts), angular_momentum(num_pts))
2831 0 : do i = 1, num_pts
2832 0 : read(iounit,*,iostat=ierr) xq(i), angular_momentum(i)
2833 0 : if (ierr /= 0) then
2834 0 : close(iounit)
2835 : write(*, '(a)') &
2836 0 : 'failed while trying to read ' // trim(s% job% relax_angular_momentum_filename)
2837 0 : write(*,*) 'line', i+1
2838 0 : write(*,*) 'perhaps wrong info in 1st line?'
2839 0 : write(*,*) '1st line must have num_pts'
2840 0 : deallocate(xq,angular_momentum)
2841 0 : return
2842 : end if
2843 : end do
2844 0 : close(iounit)
2845 : call star_relax_angular_momentum(id, s% job% max_steps_to_relax_angular_momentum, &
2846 0 : num_pts, angular_momentum, xq, ierr)
2847 0 : deallocate(xq,angular_momentum)
2848 0 : end subroutine do_relax_initial_angular_momentum
2849 :
2850 0 : subroutine do_relax_initial_entropy(ierr)
2851 : use utils_lib
2852 : use eos_def
2853 : integer, intent(out) :: ierr
2854 : ! arrays into which data from the input file is read.
2855 : ! in case any of the eos* options is used, input from the
2856 : ! file is read into var1 and var2, and the chosen eos function
2857 : ! is used to extract the entropy from that pair.
2858 0 : real(dp), pointer :: xq(:), entropy(:)
2859 : real(dp) :: var1, var2
2860 : integer :: num_pts, i, k, iounit
2861 : ! these are needed to call eosPT_get
2862 : real(dp) :: Rho, log10Rho, dlnRho_dlnPgas_const_T, dlnRho_dlnT_const_Pgas
2863 : real(dp) :: log10T
2864 : ! these are needed to call eosDT_get_T
2865 : real(dp) :: T_guess_gas, T_guess_rad, logT_guess
2866 : integer :: eos_calls
2867 : ! these are used for all eos calls
2868 : real(dp), dimension(num_eos_basic_results) :: res, d_dlnd, d_dlnT
2869 0 : real(dp), dimension(num_eos_d_dxa_results, s% species) :: d_dxa
2870 : real(dp), parameter :: logT_tol = 1d-8, logE_tol = 1d-8
2871 : integer, parameter :: MAX_ITERS = 20
2872 : include 'formats'
2873 :
2874 0 : write(*,'(A)')
2875 0 : write(*,1) 'relax_initial_entropy'
2876 :
2877 : open(newunit=iounit, file=trim(s% job% relax_entropy_filename), &
2878 0 : status='old', action='read', iostat=ierr)
2879 0 : if (ierr /= 0) then
2880 0 : write(*,*) 'open failed', ierr, iounit
2881 0 : write(*, '(a)') 'failed to open "' // trim(s% job% relax_entropy_filename)//'"'
2882 0 : return
2883 : end if
2884 0 : read(iounit, *, iostat=ierr) num_pts
2885 0 : if (ierr /= 0) then
2886 0 : close(iounit)
2887 : write(*, '(a)') 'failed while trying to read 1st line of ' // &
2888 0 : trim(s% job% relax_entropy_filename)
2889 0 : return
2890 : end if
2891 0 : if (.not. (s% job% get_entropy_for_relax_from_eos == '' .or. &
2892 : s% job% get_entropy_for_relax_from_eos == 'eosDT' .or. &
2893 : s% job% get_entropy_for_relax_from_eos == 'eosPT' .or. &
2894 : s% job% get_entropy_for_relax_from_eos == 'eosDE')) then
2895 0 : ierr = 1
2896 0 : write(*,*) 'invalid value for get_entropy_for_relax_from_eos =', &
2897 0 : s% job% get_entropy_for_relax_from_eos
2898 : end if
2899 0 : allocate(xq(num_pts), entropy(num_pts))
2900 0 : do i = 1, num_pts
2901 0 : if (s% job% get_entropy_for_relax_from_eos == '') then
2902 0 : read(iounit,*,iostat=ierr) xq(i), entropy(i)
2903 : else
2904 0 : read(iounit,*,iostat=ierr) xq(i), var1, var2
2905 : ! get nearest value matching xq for the composition TODO: interpolate
2906 0 : do k=1, s% nz-1
2907 0 : if(1-s% q(k) <= xq(i) .and. 1-s% q(k+1) >= xq(i)) then
2908 : exit
2909 : end if
2910 : end do
2911 : ! get entropy
2912 0 : if (s% job% get_entropy_for_relax_from_eos == 'eosDT') then
2913 : call eosDT_get( &
2914 : s% eos_handle, &
2915 : s% species, s% chem_id, s% net_iso, s% xa(:,k), &
2916 : var1, log10(var1), var2, log10(var2), &
2917 0 : res, d_dlnd, d_dlnT, d_dxa, ierr)
2918 0 : if (ierr /= 0) then
2919 0 : write(*,*) "failed in eosDT_get"
2920 0 : return
2921 : end if
2922 0 : entropy(i) = exp(res(i_lnS))
2923 0 : else if (s% job% get_entropy_for_relax_from_eos == 'eosPT') then
2924 : call eosPT_get( &
2925 : s% eos_handle, &
2926 : s% species, s% chem_id, s% net_iso, s% xa(:,k), &
2927 : var1, log10(var1), var2, log10(var2), &
2928 : Rho, log10Rho, dlnRho_dlnPgas_const_T, dlnRho_dlnT_const_Pgas, &
2929 0 : res, d_dlnd, d_dlnT, d_dxa, ierr)
2930 0 : if (ierr /= 0) then
2931 0 : write(*,*) "failed in eosPT_get"
2932 0 : return
2933 : end if
2934 0 : entropy(i) = exp(res(i_lnS))
2935 : else
2936 0 : T_guess_gas = 2*var2*s% abar(k)*mp/(3*kerg*(1+s% zbar(k))) ! ideal gas (var2=energy)
2937 0 : T_guess_rad = pow(var2/crad,0.25d0)
2938 0 : logT_guess = log10(min(T_guess_gas,T_guess_rad))
2939 : call eosDT_get_T( &
2940 : s% eos_handle, &
2941 : s% species, s% chem_id, s% net_iso, s% xa(:,k), &
2942 : log10(var1), i_lnE, log10(var2)*ln10, &
2943 : logT_tol, logE_tol*ln10, MAX_ITERS, logT_guess, &
2944 : arg_not_provided, arg_not_provided, arg_not_provided, arg_not_provided, &
2945 : log10T, res, d_dlnd, d_dlnT, d_dxa, &
2946 0 : eos_calls, ierr)
2947 0 : if (ierr /= 0) then
2948 0 : write(*,*) "failed in eosDT_get_T (as eosDE)"
2949 0 : return
2950 : end if
2951 0 : entropy(i) = exp(res(i_lnS))
2952 : end if
2953 : end if
2954 0 : if (ierr /= 0) then
2955 0 : close(iounit)
2956 : write(*, '(a)') &
2957 0 : 'failed while trying to read ' // trim(s% job% relax_entropy_filename)
2958 0 : write(*,*) 'line', i+1
2959 0 : write(*,*) 'perhaps wrong info in 1st line?'
2960 0 : write(*,*) '1st line must have num_pts'
2961 0 : deallocate(xq,entropy)
2962 0 : return
2963 : end if
2964 : end do
2965 0 : close(iounit)
2966 0 : call star_relax_entropy(id, s% job% max_steps_to_relax_entropy, num_pts, entropy, xq, ierr)
2967 0 : deallocate(xq,entropy)
2968 0 : end subroutine do_relax_initial_entropy
2969 :
2970 : end subroutine do_star_job_controls_after
2971 :
2972 :
2973 1 : subroutine do_remove_center(id, s, restart, ierr)
2974 : integer, intent(in) :: id
2975 : type (star_info), pointer :: s
2976 : logical, intent(in) :: restart
2977 : integer, intent(out) :: ierr
2978 : include 'formats'
2979 :
2980 1 : if (s% job% remove_center_by_temperature > 0) then
2981 0 : write(*, 1) 'remove_center_by_temperature', s% job% remove_center_by_temperature
2982 : call star_remove_center_by_temperature( &
2983 0 : id, s% job% remove_center_by_temperature, ierr)
2984 0 : if (failed('star_remove_center_by_temperature',ierr)) return
2985 : end if
2986 :
2987 1 : if (s% job% remove_initial_center_by_temperature > 0 .and. .not. restart) then
2988 0 : write(*, 1) 'remove_initial_center_by_temperature', &
2989 0 : s% job% remove_initial_center_by_temperature
2990 : call star_remove_center_by_temperature( &
2991 0 : id, s% job% remove_initial_center_by_temperature, ierr)
2992 0 : if (failed('star_remove_center_by_temperature',ierr)) return
2993 : end if
2994 :
2995 1 : if (s% job% remove_center_by_radius_cm > s% R_center .and. &
2996 : s% job% remove_center_by_radius_cm < s% r(1)) then
2997 0 : write(*, 1) 'remove_center_by_radius_cm', &
2998 0 : s% job% remove_center_by_radius_cm
2999 : call star_remove_center_by_radius_cm( &
3000 0 : id, s% job% remove_center_by_radius_cm, ierr)
3001 0 : if (failed('star_remove_center_by_radius_cm',ierr)) return
3002 : end if
3003 :
3004 : if (s% job% remove_initial_center_by_radius_cm > s% R_center .and. &
3005 1 : s% job% remove_initial_center_by_radius_cm < s% r(1) .and. .not. restart) then
3006 0 : write(*, 1) 'remove_initial_center_by_radius_cm', &
3007 0 : s% job% remove_initial_center_by_radius_cm
3008 : call star_remove_center_by_radius_cm( &
3009 0 : id, s% job% remove_initial_center_by_radius_cm, ierr)
3010 0 : if (failed('star_remove_center_by_radius_cm',ierr)) return
3011 : end if
3012 :
3013 : if (s% job% remove_initial_center_by_he4 > 0d0 .and. &
3014 : s% job% remove_initial_center_by_he4 < 1d0 &
3015 1 : .and. .not. restart) then
3016 0 : write(*, 1) 'remove_initial_center_by_he4', &
3017 0 : s% job% remove_initial_center_by_he4
3018 : call star_remove_center_by_he4( &
3019 0 : id, s% job% remove_initial_center_by_he4, ierr)
3020 0 : if (failed('star_remove_initial_center_by_he4',ierr)) return
3021 : end if
3022 :
3023 1 : if (s% job% remove_center_by_he4 > 0d0 .and. &
3024 : s% job% remove_center_by_he4 < 1d0) then
3025 0 : write(*, 1) 'remove_center_by_he4', &
3026 0 : s% job% remove_center_by_he4
3027 : call star_remove_center_by_he4( &
3028 0 : id, s% job% remove_center_by_he4, ierr)
3029 0 : if (failed('star_remove_center_by_he4',ierr)) return
3030 : end if
3031 :
3032 : if (s% job% remove_initial_center_by_c12_o16 > 0d0 .and. &
3033 : s% job% remove_initial_center_by_c12_o16 < 1d0 &
3034 1 : .and. .not. restart) then
3035 0 : write(*, 1) 'remove_initial_center_by_c12_o16', &
3036 0 : s% job% remove_initial_center_by_c12_o16
3037 : call star_remove_center_by_c12_o16( &
3038 0 : id, s% job% remove_initial_center_by_c12_o16, ierr)
3039 0 : if (failed('star_remove_initial_center_by_c12_o16',ierr)) return
3040 : end if
3041 :
3042 1 : if (s% job% remove_center_by_c12_o16 > 0d0 .and. &
3043 : s% job% remove_center_by_c12_o16 < 1d0) then
3044 0 : write(*, 1) 'remove_center_by_c12_o16', &
3045 0 : s% job% remove_center_by_c12_o16
3046 : call star_remove_center_by_c12_o16( &
3047 0 : id, s% job% remove_center_by_c12_o16, ierr)
3048 0 : if (failed('star_remove_center_by_c12_o16',ierr)) return
3049 : end if
3050 :
3051 : if (s% job% remove_initial_center_by_si28 > 0d0 .and. &
3052 : s% job% remove_initial_center_by_si28 < 1d0 &
3053 1 : .and. .not. restart) then
3054 0 : write(*, 1) 'remove_initial_center_by_si28', &
3055 0 : s% job% remove_initial_center_by_si28
3056 : call star_remove_center_by_si28( &
3057 0 : id, s% job% remove_initial_center_by_si28, ierr)
3058 0 : if (failed('star_remove_initial_center_by_si28',ierr)) return
3059 : end if
3060 :
3061 1 : if (s% job% remove_center_by_si28 > 0d0 .and. &
3062 : s% job% remove_center_by_si28 < 1d0) then
3063 0 : write(*, 1) 'remove_center_by_si28', &
3064 0 : s% job% remove_center_by_si28
3065 : call star_remove_center_by_si28( &
3066 0 : id, s% job% remove_center_by_si28, ierr)
3067 0 : if (failed('star_remove_center_by_si28',ierr)) return
3068 : end if
3069 :
3070 : if (s% job% remove_initial_center_to_reduce_co56_ni56 > 0d0 &
3071 1 : .and. .not. restart) then
3072 0 : write(*, 1) 'remove_initial_center_to_reduce_co56_ni56', &
3073 0 : s% job% remove_initial_center_to_reduce_co56_ni56
3074 : call star_remove_center_to_reduce_co56_ni56( &
3075 0 : id, s% job% remove_initial_center_to_reduce_co56_ni56, ierr)
3076 0 : if (failed('star_remove_initial_center_to_reduce_co56_ni56',ierr)) return
3077 : end if
3078 :
3079 1 : if (s% job% remove_center_to_reduce_co56_ni56 > 0d0) then
3080 0 : write(*, 1) 'remove_center_to_reduce_co56_ni56', &
3081 0 : s% job% remove_center_to_reduce_co56_ni56
3082 : call star_remove_center_to_reduce_co56_ni56( &
3083 0 : id, s% job% remove_center_to_reduce_co56_ni56, ierr)
3084 0 : if (failed('star_remove_center_to_reduce_co56_ni56',ierr)) return
3085 : end if
3086 :
3087 : if (s% job% remove_initial_center_by_ye > 0d0 .and. &
3088 : s% job% remove_initial_center_by_ye < 1d0 &
3089 1 : .and. .not. restart) then
3090 0 : write(*, 1) 'remove_initial_center_by_ye', &
3091 0 : s% job% remove_initial_center_by_ye
3092 : call star_remove_center_by_ye( &
3093 0 : id, s% job% remove_initial_center_by_ye, ierr)
3094 0 : if (failed('star_remove_initial_center_by_ye',ierr)) return
3095 : end if
3096 :
3097 1 : if (s% job% remove_center_by_ye > 0d0 .and. &
3098 : s% job% remove_center_by_ye < 1d0) then
3099 0 : write(*, 1) 'remove_center_by_ye', &
3100 0 : s% job% remove_center_by_ye
3101 : call star_remove_center_by_ye( &
3102 0 : id, s% job% remove_center_by_ye, ierr)
3103 0 : if (failed('star_remove_center_by_ye',ierr)) return
3104 : end if
3105 :
3106 : if (s% job% remove_initial_center_by_entropy > 0d0 &
3107 1 : .and. .not. restart) then
3108 0 : write(*, 1) 'remove_initial_center_by_entropy', &
3109 0 : s% job% remove_initial_center_by_entropy
3110 : call star_remove_center_by_entropy( &
3111 0 : id, s% job% remove_initial_center_by_entropy, ierr)
3112 0 : if (failed('star_remove_initial_center_by_entropy',ierr)) return
3113 : end if
3114 :
3115 1 : if (s% job% remove_center_by_entropy > 0d0) then
3116 0 : write(*, 1) 'remove_center_by_entropy', &
3117 0 : s% job% remove_center_by_entropy
3118 : call star_remove_center_by_entropy( &
3119 0 : id, s% job% remove_center_by_entropy, ierr)
3120 0 : if (failed('star_remove_center_by_entropy',ierr)) return
3121 : end if
3122 :
3123 : if (s% job% remove_initial_center_by_infall_kms /= 0d0 &
3124 1 : .and. .not. restart) then
3125 0 : write(*, 1) 'remove_initial_center_by_infall_kms', &
3126 0 : s% job% remove_initial_center_by_infall_kms
3127 : call star_remove_center_by_infall_kms( &
3128 0 : id, s% job% remove_initial_center_by_infall_kms, ierr)
3129 0 : if (failed('star_remove_initial_center_by_infall_kms',ierr)) return
3130 : end if
3131 :
3132 1 : if (s% job% remove_center_by_infall_kms /= 0d0) then
3133 0 : write(*, 1) 'remove_center_by_infall_kms', &
3134 0 : s% job% remove_center_by_infall_kms
3135 : call star_remove_center_by_infall_kms( &
3136 0 : id, s% job% remove_center_by_infall_kms, ierr)
3137 0 : if (failed('star_remove_center_by_infall_kms',ierr)) return
3138 : end if
3139 :
3140 : if (s% job% remove_initial_center_at_inner_max_abs_v &
3141 1 : .and. .not. restart) then
3142 0 : write(*, 1) 'remove_initial_center_at_inner_max_abs_v'
3143 0 : call star_remove_center_at_inner_max_abs_v(id, ierr)
3144 0 : if (failed('remove_center_at_inner_max_abs_v',ierr)) return
3145 : end if
3146 :
3147 1 : if (s% job% remove_center_at_inner_max_abs_v) then
3148 0 : write(*, 1) 'remove_initial_center_at_inner_max_abs_v'
3149 0 : call star_remove_center_at_inner_max_abs_v(id, ierr)
3150 0 : if (failed('remove_center_at_inner_max_abs_v',ierr)) return
3151 : end if
3152 :
3153 1 : if (s% job% remove_initial_fe_core .and. .not. restart) then
3154 0 : write(*, 1) 'remove_initial_fe_core'
3155 0 : call star_remove_fe_core(id, ierr)
3156 0 : if (failed('remove_fe_core',ierr)) return
3157 : end if
3158 :
3159 1 : if (s% job% remove_fe_core) then
3160 0 : write(*, 1) 'remove_initial_fe_core'
3161 0 : call star_remove_fe_core(id, ierr)
3162 0 : if (failed('remove_fe_core',ierr)) return
3163 : end if
3164 :
3165 : if (s% job% remove_initial_center_by_mass_fraction_q > 0d0 .and. &
3166 : s% job% remove_initial_center_by_mass_fraction_q < 1d0 &
3167 1 : .and. .not. restart) then
3168 0 : write(*, 1) 'remove_initial_center_by_mass_fraction_q', &
3169 0 : s% job% remove_initial_center_by_mass_fraction_q
3170 : call star_remove_center_by_mass_fraction_q( &
3171 0 : id, s% job% remove_initial_center_by_mass_fraction_q, ierr)
3172 0 : if (failed('star_remove_initial_center_by_mass_fraction_q',ierr)) return
3173 : end if
3174 :
3175 1 : if (s% job% remove_center_by_mass_fraction_q > 0d0 .and. &
3176 : s% job% remove_center_by_mass_fraction_q < 1d0) then
3177 0 : write(*, 1) 'remove_center_by_mass_fraction_q', &
3178 0 : s% job% remove_center_by_mass_fraction_q
3179 : call star_remove_center_by_mass_fraction_q( &
3180 0 : id, s% job% remove_center_by_mass_fraction_q, ierr)
3181 0 : if (failed('star_remove_center_by_mass_fraction_q',ierr)) return
3182 : end if
3183 :
3184 1 : if (s% job% remove_center_by_delta_mass_gm > 0) then
3185 0 : write(*, 1) 'remove_center_by_delta_mass_gm', &
3186 0 : s% job% remove_center_by_delta_mass_gm
3187 : call star_remove_center_by_mass_gm(id, &
3188 0 : s% M_center + s% job% remove_center_by_delta_mass_gm, ierr)
3189 0 : if (failed('star_remove_center_by_mass',ierr)) return
3190 : end if
3191 :
3192 1 : if (s% job% remove_initial_center_by_delta_mass_gm > 0 .and. &
3193 : .not. restart) then
3194 0 : write(*, 1) 'remove_initial_center_by_delta_mass_gm', &
3195 0 : s% job% remove_initial_center_by_delta_mass_gm
3196 : call star_remove_center_by_mass_gm(id, &
3197 0 : s% M_center + s% job% remove_initial_center_by_delta_mass_gm, ierr)
3198 0 : if (failed('star_remove_center_by_mass',ierr)) return
3199 : end if
3200 :
3201 1 : if (s% job% remove_center_by_delta_mass_Msun > 0) then
3202 0 : write(*, 1) 'remove_center_by_delta_mass_Msun', &
3203 0 : s% job% remove_center_by_delta_mass_Msun
3204 : call star_remove_center_by_mass_gm(id, &
3205 0 : s% M_center + s% job% remove_center_by_delta_mass_Msun*Msun, ierr)
3206 0 : if (failed('star_remove_center_by_mass',ierr)) return
3207 : end if
3208 :
3209 1 : if (s% job% remove_initial_center_by_delta_mass_Msun > 0 .and. &
3210 : .not. restart) then
3211 0 : write(*, 1) 'remove_initial_center_by_delta_mass_Msun', &
3212 0 : s% job% remove_initial_center_by_delta_mass_Msun
3213 : call star_remove_center_by_mass_gm(id, &
3214 0 : s% M_center + s% job% remove_initial_center_by_delta_mass_Msun*Msun, ierr)
3215 0 : if (failed('star_remove_center_by_mass',ierr)) return
3216 : end if
3217 :
3218 1 : if (s% job% remove_center_by_mass_gm > s% M_center .and. &
3219 : s% job% remove_center_by_mass_gm < s% m(1)) then
3220 0 : write(*, 1) 'remove_center_by_mass_gm', &
3221 0 : s% job% remove_center_by_mass_gm
3222 : call star_remove_center_by_mass_gm( &
3223 0 : id, s% job% remove_center_by_mass_gm, ierr)
3224 0 : if (failed('star_remove_center_by_mass_gm',ierr)) return
3225 : end if
3226 :
3227 : if (s% job% remove_initial_center_by_mass_gm > s% M_center .and. &
3228 1 : s% job% remove_initial_center_by_mass_gm < s% m(1) .and. .not. restart) then
3229 0 : write(*, 1) 'remove_initial_center_by_mass_gm', &
3230 0 : s% job% remove_initial_center_by_mass_gm
3231 : call star_remove_center_by_mass_gm( &
3232 0 : id, s% job% remove_initial_center_by_mass_gm, ierr)
3233 0 : if (failed('star_remove_center_by_mass_gm',ierr)) return
3234 : end if
3235 :
3236 1 : if (s% job% remove_center_by_mass_Msun > s% M_center/Msun .and. &
3237 : s% job% remove_center_by_mass_Msun < s% m(1)/Msun) then
3238 0 : write(*, 1) 'remove_center_by_mass_Msun', &
3239 0 : s% job% remove_center_by_mass_Msun
3240 : call star_remove_center_by_mass_gm( &
3241 0 : id, s% job% remove_center_by_mass_Msun*Msun, ierr)
3242 0 : if (failed('star_remove_center_by_mass_Msun',ierr)) return
3243 : end if
3244 :
3245 : if (s% job% remove_initial_center_by_mass_Msun > s% M_center/Msun .and. &
3246 1 : s% job% remove_initial_center_by_mass_Msun < s% m(1)/Msun .and. &
3247 : .not. restart) then
3248 0 : write(*, 1) 'remove_initial_center_by_mass_Msun', &
3249 0 : s% job% remove_initial_center_by_mass_Msun
3250 : call star_remove_center_by_mass_gm( &
3251 0 : id, s% job% remove_initial_center_by_mass_Msun*Msun, ierr)
3252 0 : if (failed('star_remove_center_by_mass_Msun',ierr)) return
3253 : end if
3254 :
3255 1 : if (s% job% remove_center_by_radius_Rsun > s% R_center/Rsun .and. &
3256 : s% job% remove_center_by_radius_Rsun < s% r(1)/Rsun) then
3257 0 : write(*, 1) 'remove_center_by_radius_Rsun', &
3258 0 : s% job% remove_center_by_radius_Rsun
3259 : call star_remove_center_by_radius_cm( &
3260 0 : id, s% job% remove_center_by_radius_Rsun*Rsun, ierr)
3261 0 : if (failed('star_remove_center_by_radius_Rsun',ierr)) return
3262 : end if
3263 :
3264 : if (s% job% remove_initial_center_by_radius_Rsun > s% R_center/Rsun .and. &
3265 1 : s% job% remove_initial_center_by_radius_Rsun < s% r(1)/Rsun .and. &
3266 : .not. restart) then
3267 0 : write(*, 1) 'remove_initial_center_by_radius_Rsun', &
3268 0 : s% job% remove_initial_center_by_radius_Rsun
3269 : call star_remove_center_by_radius_cm( &
3270 0 : id, s% job% remove_initial_center_by_radius_Rsun*Rsun, ierr)
3271 0 : if (failed('star_remove_center_by_radius_Rsun',ierr)) return
3272 : end if
3273 :
3274 1 : if (s% job% remove_initial_center_at_cell_k > 0 .and. .not. restart .and. &
3275 : s% job% remove_initial_center_at_cell_k <= s% nz) then
3276 0 : write(*, 2) 'remove_initial_center_at_cell_k', s% job% remove_initial_center_at_cell_k
3277 : call star_remove_center_at_cell_k( &
3278 0 : id, s% job% remove_initial_center_at_cell_k, ierr)
3279 0 : if (failed('star_remove_center_at_cell_k',ierr)) return
3280 : end if
3281 :
3282 1 : if (s% job% remove_center_at_cell_k > 0 .and. &
3283 : s% job% remove_center_at_cell_k <= s% nz) then
3284 0 : write(*, 2) 'remove_center_at_cell_k', s% job% remove_center_at_cell_k
3285 0 : call star_remove_center_at_cell_k(id, s% job% remove_center_at_cell_k, ierr)
3286 0 : if (failed('star_remove_center_at_cell_k',ierr)) return
3287 : end if
3288 :
3289 : end subroutine do_remove_center
3290 :
3291 :
3292 1 : subroutine do_remove_initial_surface(id,s,restart,ierr)
3293 : integer, intent(in) :: id
3294 : type (star_info), pointer :: s
3295 : logical, intent(in) :: restart
3296 : integer, intent(out) :: ierr
3297 :
3298 : include 'formats'
3299 :
3300 1 : ierr = 0
3301 :
3302 1 : if (s% job% remove_initial_surface_at_he_core_boundary > 0 .and. .not. restart) then
3303 0 : write(*, 1) 'remove_initial_surface_at_he_core_boundary', &
3304 0 : s% job% remove_initial_surface_at_he_core_boundary
3305 : call star_remove_surface_at_he_core_boundary( &
3306 0 : id, s% job% remove_initial_surface_at_he_core_boundary, ierr)
3307 0 : if (failed('star_remove_surface_at_he_core_boundary',ierr)) return
3308 : end if
3309 :
3310 1 : if (s% job% remove_initial_surface_by_optical_depth > 0 .and. .not. restart) then
3311 0 : write(*, 1) 'remove_initial_surface_by_optical_depth', &
3312 0 : s% job% remove_initial_surface_by_optical_depth
3313 : call star_remove_surface_by_optical_depth( &
3314 0 : id, s% job% remove_initial_surface_by_optical_depth, ierr)
3315 0 : if (failed('star_remove_surface_by_optical_depth',ierr)) return
3316 : end if
3317 :
3318 1 : if (s% job% remove_initial_surface_by_density > 0 .and. .not. restart) then
3319 0 : write(*, 1) 'remove_initial_surface_by_density', &
3320 0 : s% job% remove_initial_surface_by_density
3321 : call star_remove_surface_by_density( &
3322 0 : id, s% job% remove_initial_surface_by_density, ierr)
3323 0 : if (failed('star_remove_surface_by_density',ierr)) return
3324 : end if
3325 :
3326 1 : if (s% job% remove_initial_surface_by_pressure > 0 .and. .not. restart) then
3327 0 : write(*, 1) 'remove_initial_surface_by_pressure', &
3328 0 : s% job% remove_initial_surface_by_pressure
3329 : call star_remove_surface_by_pressure( &
3330 0 : id, s% job% remove_initial_surface_by_pressure, ierr)
3331 0 : if (failed('star_remove_surface_by_pressure',ierr)) return
3332 : end if
3333 :
3334 : if (s% job% remove_initial_surface_by_radius_cm > s% R_center .and. &
3335 1 : s% job% remove_initial_surface_by_radius_cm < s% r(1) .and. .not. restart) then
3336 0 : write(*, 1) 'remove_initial_surface_by_radius_cm', &
3337 0 : s% job% remove_initial_surface_by_radius_cm
3338 : call star_remove_surface_by_radius_cm( &
3339 0 : id, s% job% remove_initial_surface_by_radius_cm, ierr)
3340 0 : if (failed('star_remove_surface_by_radius_cm',ierr)) return
3341 : end if
3342 :
3343 : if (s% job% remove_initial_surface_by_mass_fraction_q > 0d0 .and. &
3344 : s% job% remove_initial_surface_by_mass_fraction_q < 1d0 &
3345 1 : .and. .not. restart) then
3346 0 : write(*, 1) 'remove_initial_surface_by_mass_fraction_q', &
3347 0 : s% job% remove_initial_surface_by_mass_fraction_q
3348 : call star_remove_surface_by_mass_fraction_q( &
3349 0 : id, s% job% remove_initial_surface_by_mass_fraction_q, ierr)
3350 0 : if (failed('star_remove_initial_surface_by_mass_fraction_q',ierr)) return
3351 : end if
3352 :
3353 : if (s% job% remove_initial_surface_by_mass_gm > s% M_center .and. &
3354 1 : s% job% remove_initial_surface_by_mass_gm < s% m(1) .and. .not. restart) then
3355 0 : write(*, 1) 'remove_initial_surface_by_mass_gm', &
3356 0 : s% job% remove_initial_surface_by_mass_gm
3357 : call star_remove_surface_by_mass_gm( &
3358 0 : id, s% job% remove_initial_surface_by_mass_gm, ierr)
3359 0 : if (failed('star_remove_surface_by_mass_gm',ierr)) return
3360 : end if
3361 :
3362 : if (s% job% remove_initial_surface_by_radius_Rsun > s% R_center/Rsun .and. &
3363 1 : s% job% remove_initial_surface_by_radius_Rsun < s% r(1)/Rsun .and. &
3364 : .not. restart) then
3365 0 : write(*, 1) 'remove_initial_surface_by_radius_Rsun', &
3366 0 : s% job% remove_initial_surface_by_radius_Rsun
3367 : call star_remove_surface_by_radius_cm( &
3368 0 : id, s% job% remove_initial_surface_by_radius_Rsun*Rsun, ierr)
3369 0 : if (failed('star_remove_surface_by_radius_Rsun',ierr)) return
3370 : end if
3371 :
3372 : if (s% job% remove_initial_surface_by_mass_Msun > s% M_center/Msun .and. &
3373 1 : s% job% remove_initial_surface_by_mass_Msun < s% m(1)/Msun .and. &
3374 : .not. restart) then
3375 0 : write(*, 1) 'remove_initial_surface_by_mass_Msun', &
3376 0 : s% job% remove_initial_surface_by_mass_Msun
3377 : call star_remove_surface_by_mass_gm( &
3378 0 : id, s% job% remove_initial_surface_by_mass_Msun*Msun, ierr)
3379 0 : if (failed('star_remove_surface_by_mass_Msun',ierr)) return
3380 : end if
3381 :
3382 1 : if (s% job% remove_initial_surface_by_v_surf_km_s > 0 .and. .not. restart) then
3383 0 : write(*, 2) 'remove_initial_surface_by_v_surf_km_s', &
3384 0 : s% job% remove_initial_surface_by_v_surf_km_s
3385 : call star_remove_surface_by_v_surf_km_s( &
3386 0 : id, s% job% remove_initial_surface_by_v_surf_km_s, ierr)
3387 0 : if (failed('star_remove_surface_by_v_surf_km_s',ierr)) return
3388 : end if
3389 :
3390 1 : if (s% job% remove_initial_surface_by_v_surf_div_cs > 0 .and. .not. restart) then
3391 0 : write(*, 2) 'remove_initial_surface_by_v_surf_div_cs', &
3392 0 : s% job% remove_initial_surface_by_v_surf_div_cs
3393 : call star_remove_surface_by_v_surf_div_cs( &
3394 0 : id, s% job% remove_initial_surface_by_v_surf_div_cs, ierr)
3395 0 : if (failed('star_remove_surface_by_v_surf_div_cs',ierr)) return
3396 : end if
3397 :
3398 1 : if (s% job% remove_initial_surface_by_v_surf_div_v_escape > 0 .and. .not. restart) then
3399 0 : write(*, 2) 'remove_initial_surface_by_v_surf_div_v_escape', &
3400 0 : s% job% remove_initial_surface_by_v_surf_div_v_escape
3401 : call star_remove_surface_by_v_surf_div_v_escape( &
3402 0 : id, s% job% remove_initial_surface_by_v_surf_div_v_escape, ierr)
3403 0 : if (failed('star_remove_surface_by_v_surf_div_v_escape',ierr)) return
3404 : end if
3405 :
3406 1 : if (s% job% remove_initial_surface_at_cell_k > 0 .and. .not. restart .and. &
3407 : s% job% remove_initial_surface_at_cell_k <= s% nz) then
3408 0 : write(*, 2) 'remove_initial_surface_at_cell_k', s% job% remove_initial_surface_at_cell_k
3409 : call star_remove_surface_at_cell_k( &
3410 0 : id, s% job% remove_initial_surface_at_cell_k, ierr)
3411 0 : if (failed('star_remove_surface_at_cell_k',ierr)) return
3412 : end if
3413 :
3414 : end subroutine do_remove_initial_surface
3415 :
3416 :
3417 12 : subroutine do_remove_surface(id,s,ierr)
3418 : integer, intent(in) :: id
3419 : type (star_info), pointer :: s
3420 : integer, intent(out) :: ierr
3421 :
3422 : include 'formats'
3423 :
3424 12 : ierr = 0
3425 :
3426 12 : if (s% job% remove_surface_at_he_core_boundary > 0) then
3427 : !write(*, 1) 'remove_surface_at_he_core_boundary', s% job% remove_surface_at_he_core_boundary
3428 : call star_remove_surface_at_he_core_boundary( &
3429 0 : id, s% job% remove_surface_at_he_core_boundary, ierr)
3430 0 : if (failed('star_remove_surface_at_he_core_boundary',ierr)) return
3431 : end if
3432 :
3433 12 : if (s% job% remove_surface_by_optical_depth > 0) then
3434 : !write(*, 1) 'remove_surface_by_optical_depth', s% job% remove_surface_by_optical_depth
3435 : call star_remove_surface_by_optical_depth( &
3436 0 : id, s% job% remove_surface_by_optical_depth, ierr)
3437 0 : if (failed('star_remove_surface_by_optical_depth',ierr)) return
3438 : end if
3439 :
3440 12 : if (s% job% remove_surface_by_density > 0) then
3441 : !write(*, 1) 'remove_surface_by_density', s% job% remove_surface_by_density
3442 : call star_remove_surface_by_density( &
3443 0 : id, s% job% remove_surface_by_density, ierr)
3444 0 : if (failed('star_remove_surface_by_density',ierr)) return
3445 : end if
3446 :
3447 12 : if (s% job% remove_surface_by_pressure > 0) then
3448 : !write(*, 1) 'remove_surface_by_pressure', s% job% remove_surface_by_pressure
3449 : call star_remove_surface_by_pressure( &
3450 0 : id, s% job% remove_surface_by_pressure, ierr)
3451 0 : if (failed('star_remove_surface_by_pressure',ierr)) return
3452 : end if
3453 :
3454 12 : if (s% job% remove_surface_by_radius_cm > s% R_center .and. &
3455 : s% job% remove_surface_by_radius_cm < s% r(1)) then
3456 : !write(*, 1) 'remove_surface_by_radius_cm', s% job% remove_surface_by_radius_cm
3457 : call star_remove_surface_by_radius_cm( &
3458 0 : id, s% job% remove_surface_by_radius_cm, ierr)
3459 0 : if (failed('star_remove_surface_by_radius_cm',ierr)) return
3460 : end if
3461 :
3462 12 : if (s% job% remove_surface_by_mass_fraction_q > 0d0 .and. &
3463 : s% job% remove_surface_by_mass_fraction_q < 1d0) then
3464 : !write(*, 1) 'remove_surface_by_mass_fraction_q', &
3465 : ! s% job% remove_surface_by_mass_fraction_q
3466 : call star_remove_surface_by_mass_fraction_q( &
3467 0 : id, s% job% remove_surface_by_mass_fraction_q, ierr)
3468 0 : if (failed('star_remove_surface_by_mass_fraction_q',ierr)) return
3469 : end if
3470 :
3471 12 : if (s% job% remove_surface_by_mass_gm > s% M_center .and. &
3472 : s% job% remove_surface_by_mass_gm < s% m(1)) then
3473 : !write(*, 1) 'remove_surface_by_mass_gm', &
3474 : ! s% job% remove_surface_by_mass_gm
3475 : call star_remove_surface_by_mass_gm( &
3476 0 : id, s% job% remove_surface_by_mass_gm, ierr)
3477 0 : if (failed('star_remove_surface_by_mass_gm',ierr)) return
3478 : end if
3479 :
3480 12 : if (s% job% remove_surface_by_radius_Rsun > s% R_center/Rsun .and. &
3481 : s% job% remove_surface_by_radius_Rsun < s% r(1)/Rsun) then
3482 : !write(*, 1) 'remove_surface_by_radius_Rsun', &
3483 : ! s% job% remove_surface_by_radius_Rsun
3484 : call star_remove_surface_by_radius_cm( &
3485 0 : id, s% job% remove_surface_by_radius_Rsun*Rsun, ierr)
3486 0 : if (failed('star_remove_surface_by_radius_Rsun',ierr)) return
3487 : end if
3488 :
3489 12 : if (s% job% remove_surface_by_mass_Msun > s% M_center/Msun .and. &
3490 : s% job% remove_surface_by_mass_Msun < s% m(1)/Msun) then
3491 : !write(*, 1) 'remove_surface_by_mass_Msun', &
3492 : ! s% job% remove_surface_by_mass_Msun
3493 : call star_remove_surface_by_mass_gm( &
3494 0 : id, s% job% remove_surface_by_mass_Msun*Msun, ierr)
3495 0 : if (failed('star_remove_surface_by_mass_Msun',ierr)) return
3496 : end if
3497 :
3498 12 : if (s% job% remove_surface_by_v_surf_km_s > 0) then
3499 : !write(*, 2) 'remove_surface_by_v_surf_km_s', s% job% remove_surface_by_v_surf_km_s
3500 0 : call star_remove_surface_by_v_surf_km_s(id, s% job% remove_surface_by_v_surf_km_s, ierr)
3501 0 : if (failed('star_remove_surface_by_v_surf_km_s',ierr)) return
3502 : end if
3503 :
3504 12 : if (s% job% remove_surface_by_v_surf_div_cs > 0) then
3505 : !write(*, 1) 'remove_surface_by_v_surf_div_cs', s% job% remove_surface_by_v_surf_div_cs
3506 0 : call star_remove_surface_by_v_surf_div_cs(id, s% job% remove_surface_by_v_surf_div_cs, ierr)
3507 0 : if (failed('star_remove_surface_by_v_surf_div_cs',ierr)) return
3508 : end if
3509 :
3510 12 : if (s% job% remove_surface_by_v_surf_div_v_escape > 0) then
3511 : !write(*, 2) 'remove_surface_by_v_surf_div_v_escape', s% job% remove_surface_by_v_surf_div_v_escape
3512 0 : call star_remove_surface_by_v_surf_div_v_escape(id, s% job% remove_surface_by_v_surf_div_v_escape, ierr)
3513 0 : if (failed('star_remove_surface_by_v_surf_div_v_escape',ierr)) return
3514 : end if
3515 :
3516 12 : if (s% job% remove_surface_at_cell_k > 0 .and. &
3517 : s% job% remove_surface_at_cell_k <= s% nz) then
3518 : !write(*, 2) 'remove_surface_at_cell_k', s% job% remove_surface_at_cell_k
3519 0 : call star_remove_surface_at_cell_k(id, s% job% remove_surface_at_cell_k, ierr)
3520 0 : if (failed('star_remove_surface_at_cell_k',ierr)) return
3521 : end if
3522 :
3523 : end subroutine do_remove_surface
3524 :
3525 :
3526 2 : subroutine resolve_inlist_fname(inlist_out,inlist_opt)
3527 :
3528 : use ISO_FORTRAN_ENV
3529 :
3530 : character(len=*),intent(out) :: inlist_out
3531 : character(len=*),optional :: inlist_opt
3532 :
3533 : integer :: status
3534 :
3535 : ! initialize inlist_out as empty
3536 2 : inlist_out = ''
3537 :
3538 2 : if (.not. MESA_INLIST_RESOLVED) then
3539 2 : if (COMMAND_ARGUMENT_COUNT() >= 1) then
3540 :
3541 : ! Get filename from the first command-line argument
3542 :
3543 0 : call GET_COMMAND_ARGUMENT(1, inlist_out, STATUS=status)
3544 0 : if (status /= 0) inlist_out = ''
3545 :
3546 : else
3547 :
3548 : ! Get filename from the MESA_INLIST environment variable
3549 :
3550 2 : call GET_ENVIRONMENT_VARIABLE('MESA_INLIST', inlist_out, STATUS=status)
3551 2 : if (status /= 0) inlist_out = ''
3552 :
3553 : end if
3554 : end if
3555 :
3556 2 : if (inlist_out == '') then
3557 :
3558 2 : if (PRESENT(inlist_opt)) then
3559 1 : inlist_out = inlist_opt
3560 : else
3561 1 : inlist_out = 'inlist'
3562 : end if
3563 :
3564 : end if
3565 :
3566 2 : return
3567 :
3568 2 : end subroutine resolve_inlist_fname
3569 :
3570 :
3571 1 : subroutine add_fpe_checks(id, s, ierr)
3572 : integer, intent(in) :: id
3573 : type (star_info), pointer :: s
3574 : integer, intent(out) :: ierr
3575 :
3576 : character(len=1) :: fpe_check
3577 : integer :: status
3578 :
3579 : include 'formats'
3580 :
3581 1 : ierr = 0
3582 :
3583 1 : call GET_ENVIRONMENT_VARIABLE('MESA_FPE_CHECKS_ON', fpe_check, STATUS=status)
3584 1 : if (status /= 0) return
3585 :
3586 0 : if (fpe_check(1:1)=="1") then
3587 0 : write(*,*) "FPE checking is on"
3588 0 : s% fill_arrays_with_nans = .true.
3589 : end if
3590 :
3591 : end subroutine add_fpe_checks
3592 :
3593 :
3594 1 : subroutine multiply_tolerances(id, s, ierr)
3595 : integer, intent(in) :: id
3596 : type (star_info), pointer :: s
3597 : integer, intent(out) :: ierr
3598 : integer :: status
3599 :
3600 : real(dp), save :: test_suite_res_factor = 1
3601 : character(len=20) :: test_suite_resolution_factor_str
3602 :
3603 : include 'formats'
3604 :
3605 1 : ierr = 0
3606 : call GET_ENVIRONMENT_VARIABLE('MESA_TEST_SUITE_RESOLUTION_FACTOR', &
3607 1 : test_suite_resolution_factor_str, STATUS=status)
3608 1 : if (status /= 0) return
3609 :
3610 0 : if (test_suite_resolution_factor_str /= "") then
3611 0 : read(test_suite_resolution_factor_str, *) test_suite_res_factor
3612 0 : write(*,*) ""
3613 0 : write(*,*) "***"
3614 0 : write(*,*) "MESA_TEST_SUITE_RESOLUTION_FACTOR set to", test_suite_res_factor
3615 0 : write(*,*) "***"
3616 0 : write(*,*) "Warning: This environment variable is for testing purposes"
3617 0 : write(*,*) " and should be set to 1 during normal MESA use."
3618 0 : write(*,*) "***"
3619 0 : write(*,*) "Multiplying mesh_delta_coeff and time_delta_coeff by this factor,"
3620 0 : write(*,*) "and max_model_number by its inverse twice:"
3621 0 : write(*,*) ""
3622 0 : write(*,*) " old mesh_delta_coeff = ", s% mesh_delta_coeff
3623 0 : s% mesh_delta_coeff = test_suite_res_factor * s% mesh_delta_coeff
3624 0 : write(*,*) " new mesh_delta_coeff = ", s% mesh_delta_coeff
3625 0 : write(*,*) ""
3626 0 : write(*,*) " old time_delta_coeff = ", s% time_delta_coeff
3627 0 : s% time_delta_coeff = test_suite_res_factor * s% time_delta_coeff
3628 0 : write(*,*) " new time_delta_coeff = ", s% time_delta_coeff
3629 0 : write(*,*) ""
3630 0 : write(*,*) " old max_model_number = ", s% max_model_number
3631 0 : s% max_model_number = s% max_model_number / test_suite_res_factor / test_suite_res_factor
3632 0 : write(*,*) " new max_model_number = ", s% max_model_number
3633 0 : write(*,*) ""
3634 : end if
3635 :
3636 : end subroutine multiply_tolerances
3637 :
3638 :
3639 1 : subroutine pgstar_env_check(id, s, ierr)
3640 : integer, intent(in) :: id
3641 : type (star_info), pointer :: s
3642 : integer, intent(out) :: ierr
3643 : character(len=5) :: flag
3644 : integer :: status
3645 :
3646 : include 'formats'
3647 :
3648 1 : ierr = 0
3649 :
3650 1 : call get_environment_variable('MESA_FORCE_PGSTAR_FLAG', flag, STATUS=status)
3651 1 : if (status /= 0) return
3652 :
3653 0 : select case (trim(flag))
3654 : case ("TRUE", "true")
3655 0 : write(*,*) "PGSTAR forced on"
3656 0 : s% job% pgstar_flag = .true.
3657 : case ("FALSE", "false")
3658 0 : write(*,*) "PGSTAR forced off"
3659 0 : s% job% pgstar_flag = .false.
3660 : end select
3661 :
3662 : end subroutine pgstar_env_check
3663 :
3664 : end module run_star_support
|
