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 :
1096 : subroutine relax_Tsurf_factor(s)
1097 : type (star_info), pointer :: s
1098 : real(dp) :: next
1099 : include 'formats'
1100 : write(*,*) 'relax_to_this_Tsurf_factor < s% Tsurf_factor', &
1101 : s% job% relax_to_this_Tsurf_factor < s% Tsurf_factor
1102 : write(*,1) 'relax_to_this_Tsurf_factor', s% job% relax_to_this_Tsurf_factor
1103 : write(*,1) 's% Tsurf_factor', s% Tsurf_factor
1104 : if (s% job% relax_to_this_Tsurf_factor < s% Tsurf_factor) then
1105 : next = exp10(safe_log10(s% Tsurf_factor) - s% job% dlogTsurf_factor)
1106 : if (next < s% job% relax_to_this_Tsurf_factor) &
1107 : next = s% job% relax_to_this_Tsurf_factor
1108 : else
1109 : next = exp10(safe_log10(s% Tsurf_factor) + s% job% dlogTsurf_factor)
1110 : if (next > s% job% relax_to_this_Tsurf_factor) &
1111 : next = s% job% relax_to_this_Tsurf_factor
1112 : end if
1113 : s% Tsurf_factor = next
1114 : write(*,1) 'relax_Tsurf_factor', next, s% job% relax_to_this_Tsurf_factor
1115 : end subroutine relax_Tsurf_factor
1116 :
1117 :
1118 1 : subroutine check_if_want_to_stop_warnings(s)
1119 : use utils_lib
1120 : type (star_info), pointer :: s
1121 : character (len=200) :: fname
1122 : integer :: iounit, ierr
1123 1 : ierr = 0
1124 1 : if (s% warn_when_large_rel_run_E_err < 1d2) then
1125 1 : fname = trim(mesa_dir) // '/stop_warnings_for_rel_E_err'
1126 : open(newunit=iounit, file=trim(fname), &
1127 1 : status='old', action='read', iostat=ierr)
1128 1 : if (ierr == 0) then
1129 0 : close(iounit)
1130 0 : s% warn_when_large_rel_run_E_err = 1d99
1131 0 : write(*,*) 'turn off warnings for rel_run_E_err'
1132 : end if
1133 : end if
1134 : ierr = 0
1135 1 : end subroutine check_if_want_to_stop_warnings
1136 :
1137 :
1138 11 : logical function stop_is_requested(s)
1139 : type (star_info), pointer :: s
1140 : logical :: file_exists
1141 11 : stop_is_requested = .false.
1142 11 : if (mod(s% model_number,100) /= 0) return
1143 1 : if (len_trim(s% job% stop_if_this_file_exists) == 0) return
1144 1 : inquire(file=trim(s% job% stop_if_this_file_exists), exist=file_exists)
1145 1 : if (.not. file_exists) return
1146 : write(*,*) 'stopping because found file ' // &
1147 0 : trim(s% job% stop_if_this_file_exists)
1148 0 : stop_is_requested = .true.
1149 0 : end function stop_is_requested
1150 :
1151 :
1152 118 : logical function failed(str,ierr)
1153 : character (len=*), intent(in) :: str
1154 : integer, intent(in) :: ierr
1155 118 : failed = (ierr /= 0)
1156 118 : if (failed) write(*, *) trim(str) // ' ierr', ierr
1157 118 : end function failed
1158 :
1159 :
1160 0 : subroutine show_times(id, s)
1161 : use utils_lib, only: utils_OMP_GET_MAX_THREADS
1162 : use num_lib, only: qsort
1163 :
1164 : integer, intent(in) :: id
1165 : type (star_info), pointer :: s
1166 :
1167 : integer, parameter :: max_num_items = 50
1168 : character(len=60) :: item_names(max_num_items)
1169 : real(dp) :: item_values(max_num_items)
1170 : integer, target :: index_arry(max_num_items)
1171 : integer, pointer :: index(:)
1172 : integer :: ierr, omp_num_threads, item_num, num_items, i, j
1173 : real(dp) :: total, tmp
1174 : include 'formats'
1175 0 : ierr = 0
1176 0 : omp_num_threads = utils_OMP_GET_MAX_THREADS()
1177 : s% time_total = s% job% check_before_step_timing + &
1178 0 : s% job% check_step_loop_timing + s% job% check_after_step_timing
1179 :
1180 0 : write(*,'(A)')
1181 0 : write(*,'(a50,i18)') 'nz', s% nz
1182 0 : write(*,'(a50,i18)') 'nvar_total', s% nvar_total
1183 0 : write(*,'(a50,i18)') trim(s% net_name) // ' species', s% species
1184 0 : write(*,'(a50,i18)') 'total_num_solver_iterations', &
1185 0 : s% total_num_solver_iterations
1186 0 : write(*,'(a50,i18)') 'timing_num_get_eos_calls', &
1187 0 : s% timing_num_get_eos_calls
1188 0 : write(*,'(a50,i18)') 'timing_num_solve_eos_calls', &
1189 0 : s% timing_num_solve_eos_calls
1190 0 : write(*,'(a50,i18)') 'timing_num_get_kap_calls', &
1191 0 : s% timing_num_get_kap_calls
1192 0 : write(*,'(A)')
1193 0 : write(*,'(a50,i18)') 'threads', omp_num_threads
1194 0 : total = 0
1195 0 : item_num = 0
1196 0 : call save1('remesh', s% time_remesh, total)
1197 0 : call save1('adjust_mass', s% time_adjust_mass, total)
1198 0 : call save1('conv_premix', s% time_conv_premix, total)
1199 0 : call save1('element_diffusion', s% time_element_diffusion, total)
1200 0 : call save1('burn', s% time_solve_burn, total)
1201 0 : call save1('mix', s% time_solve_mix, total)
1202 0 : call save1('solve', s% time_struct_burn_mix, total)
1203 0 : call save1('matrix', s% time_solver_matrix, total)
1204 0 : call save1('omega_mix', s% time_solve_omega_mix, total)
1205 0 : call save1('eos', s% time_eos, total)
1206 0 : call save1('neu_and_kap', s% time_neu_kap, total)
1207 0 : call save1('net', s% time_nonburn_net, total)
1208 0 : call save1('mlt', s% time_mlt, total)
1209 0 : call save1('hydro_vars', s% time_set_hydro_vars, total)
1210 0 : call save1('mixing_info', s% time_set_mixing_info, total)
1211 0 : call save1('evolve_step', s% time_evolve_step, total)
1212 0 : call save1('run1_star', s% job% elapsed_time - total, total)
1213 0 : tmp = 0
1214 0 : call save1('total', total, tmp)
1215 :
1216 0 : num_items = item_num
1217 0 : index(1:num_items) => index_arry(1:num_items)
1218 0 : call qsort(index, num_items, item_values)
1219 :
1220 0 : write(*,'(A)')
1221 0 : write(*,'(A)')
1222 0 : do i=1,num_items
1223 0 : j = index(num_items+1-i)
1224 0 : if (item_values(j) == 0d0) cycle
1225 0 : write(*,'(a50,2f9.3)') trim(item_names(j)), &
1226 0 : item_values(j), item_values(j)/total
1227 0 : if (j == num_items) write(*,*)
1228 : end do
1229 :
1230 0 : if (s% job% step_loop_timing/s% job% elapsed_time < 0.9d0) then
1231 0 : write(*,'(A)')
1232 0 : write(*,'(A)')
1233 0 : write(*,1) 'before_step', s% job% before_step_timing/s% job% elapsed_time
1234 0 : write(*,1) 'step_loop', s% job% step_loop_timing/s% job% elapsed_time
1235 0 : write(*,1) 'after_step', s% job% after_step_timing/s% job% elapsed_time
1236 0 : write(*,'(A)')
1237 : end if
1238 0 : write(*,'(A)')
1239 0 : write(*,'(A)')
1240 :
1241 :
1242 : contains
1243 :
1244 :
1245 0 : subroutine save1(name, value, total)
1246 0 : use utils_lib, only: is_bad_num
1247 : character (len=*), intent(in) :: name
1248 : real(dp), intent(in) :: value
1249 : real(dp), intent(inout) :: total
1250 : include 'formats'
1251 0 : item_num = item_num + 1
1252 0 : item_names(item_num) = name
1253 0 : item_values(item_num) = value
1254 0 : total = total + value
1255 0 : end subroutine save1
1256 :
1257 :
1258 : end subroutine show_times
1259 :
1260 :
1261 11 : subroutine do_saves(id, ierr)
1262 : integer, intent(in) :: id
1263 : type (star_info), pointer :: s
1264 :
1265 : integer :: ierr
1266 : ierr = 0
1267 :
1268 11 : call star_ptr(id, s, ierr)
1269 11 : if (ierr/=0) return
1270 :
1271 11 : if (s% model_number == s% job% save_model_number) then
1272 0 : call star_write_model(id, s% job% save_model_filename, ierr)
1273 0 : if (failed('star_write_model',ierr)) return
1274 0 : write(*, *) 'model saved to ' // trim(s% job% save_model_filename)
1275 : end if
1276 :
1277 11 : if (s% model_number == s% job% save_photo_number) then
1278 1 : call star_write_photo(id, s% job% save_photo_filename, ierr)
1279 1 : if (failed('star_write_photo',ierr)) return
1280 1 : if (len_trim(s% job% save_photo_filename) > 0) &
1281 0 : write(*, *) 'photo saved to ' // trim(s% job% save_photo_filename)
1282 : end if
1283 :
1284 11 : if (s% model_number == s% job% save_pulse_data_for_model_number) then
1285 : call star_export_pulse_data(id, s%pulse_data_format, s%job%save_pulse_data_filename, &
1286 : s%add_center_point_to_pulse_data, s%keep_surface_point_for_pulse_data, &
1287 0 : s%add_atmosphere_to_pulse_data, ierr)
1288 0 : if (failed('star_export_pulse_data',ierr)) return
1289 : write(*, *) 'pulsation data saved to ' // &
1290 0 : trim(s% job% save_pulse_data_filename)
1291 : end if
1292 :
1293 11 : if (s% model_number == s% job% profile_model_number) then
1294 0 : write(*, '(a, i7)') 'save profile for model number', s% model_number
1295 0 : call save_profile(id, 3, ierr)
1296 0 : if (failed('save_profile',ierr)) return
1297 : end if
1298 :
1299 : end subroutine do_saves
1300 :
1301 :
1302 : subroutine write_colors_info(id, s, ierr)
1303 : use colors_lib
1304 : use colors_def
1305 : integer, intent(in) :: id
1306 : type (star_info), pointer :: s
1307 : integer, intent(out) :: ierr
1308 :
1309 : !TODO: implement me
1310 :
1311 : end subroutine write_colors_info
1312 :
1313 :
1314 : subroutine read_masses(filename, masses, nmasses, ierr)
1315 : character (len=*), intent(in) :: filename
1316 : real(dp), pointer, intent(inout) :: masses(:)
1317 : integer, intent(out) :: nmasses, ierr
1318 : call read_items(filename, masses, nmasses, 'masses', ierr)
1319 : end subroutine read_masses
1320 :
1321 :
1322 : subroutine read_items(filename, items, nitems, name, ierr)
1323 : use utils_lib
1324 : use utils_def
1325 : character (len=*), intent(in) :: filename, name
1326 : real(dp), pointer, intent(inout) :: items(:)
1327 : integer, intent(out) :: nitems, ierr
1328 :
1329 : integer :: iounit, n, i, t, capacity
1330 : character (len=strlen) :: buffer, string
1331 :
1332 : nitems = 0
1333 : if (.not. associated(items)) then
1334 : capacity = 10
1335 : allocate(items(capacity))
1336 : else
1337 : capacity = size(items,dim=1)
1338 : end if
1339 :
1340 : ierr = 0
1341 :
1342 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
1343 : if (ierr /= 0) then
1344 : write(*,*) 'failed to open file ' // trim(filename)
1345 : return
1346 : end if
1347 :
1348 : n = 0
1349 : i = 0
1350 :
1351 : do
1352 : t = token(iounit, n, i, buffer, string)
1353 : select case(t)
1354 : case(name_token)
1355 : if (string == name) then
1356 : call do_read_items(ierr)
1357 : if (ierr /= 0) then
1358 : return
1359 : end if
1360 : exit ! for now, nothing else to be read
1361 : end if
1362 : call error; return
1363 : case(eof_token)
1364 : exit
1365 : case default
1366 : call error; return
1367 : end select
1368 :
1369 : end do
1370 :
1371 : close(iounit)
1372 :
1373 : contains
1374 :
1375 :
1376 : subroutine error
1377 : ierr = -1
1378 : write(*,*) 'error in reading file' // trim(filename)
1379 : close(iounit)
1380 : end subroutine error
1381 :
1382 :
1383 : subroutine do_read_items(ierr)
1384 : integer, intent(out) :: ierr
1385 : real(dp) :: mass
1386 : ierr = 0
1387 : t = token(iounit, n, i, buffer, string)
1388 : if (t /= left_paren_token) then
1389 : call error; return
1390 : end if
1391 : mass_loop: do
1392 : t = token(iounit, n, i, buffer, string)
1393 : if (t /= name_token) then
1394 : call error; return
1395 : end if
1396 : read(string,fmt=*,iostat=ierr) mass
1397 : if (ierr /= 0) then
1398 : call error; return
1399 : end if
1400 : nitems = nitems+1
1401 : if (nitems > capacity) then
1402 : capacity = capacity + 10
1403 : call realloc_double(items,capacity,ierr)
1404 : if (ierr /= 0) then
1405 : call error; return
1406 : end if
1407 : end if
1408 : items(nitems) = mass
1409 : t = token(iounit, n, i, buffer, string)
1410 : if (t == right_paren_token) exit mass_loop
1411 : if (t /= comma_token) then
1412 : call error; return
1413 : end if
1414 : end do mass_loop
1415 : end subroutine do_read_items
1416 :
1417 :
1418 : end subroutine read_items
1419 :
1420 :
1421 0 : subroutine do_report_mass_not_fe56(s)
1422 : use const_def, only: dp
1423 : type (star_info), pointer :: s
1424 : integer :: k, fe56
1425 : real(dp) :: sumdq
1426 : include 'formats'
1427 0 : fe56 = s% net_iso(ife56)
1428 0 : if (fe56 == 0) return
1429 0 : sumdq = 0
1430 0 : do k = 1, s% nz
1431 0 : sumdq = sumdq + s% dq(k)*(1-s% xa(fe56,k))
1432 : end do
1433 0 : write(*,1) 'R', s% r(1)
1434 0 : write(*,1) 'g', s% cgrav(1)*s% mstar/(s% r(1)*s% r(1))
1435 0 : write(*,1) 'mass non fe56', s% xmstar*sumdq, sumdq
1436 0 : write(*,1) 'M_center (Msun)', s% M_center/Msun
1437 0 : write(*,1) 'xmstar (g)', s% xmstar
1438 0 : do k=1,s% nz
1439 0 : if (fe56 == maxloc(s% xa(:,k),dim=1)) then
1440 0 : write(*,2) 'mass exterior to fe56 (g)', k, (1d0 - s% q(k))*s% xmstar
1441 0 : write(*,2) 'mass coord top of fe56 (g)', k, s% q(k)*s% xmstar
1442 0 : return
1443 : end if
1444 : end do
1445 : end subroutine do_report_mass_not_fe56
1446 :
1447 :
1448 0 : subroutine do_report_cell_for_xm(s)
1449 : use const_def, only: dp
1450 : type (star_info), pointer :: s
1451 : integer :: k
1452 : real(dp) :: sumdq, dq
1453 : include 'formats'
1454 0 : dq = s% job% report_cell_for_xm/s% xmstar
1455 0 : if (dq > 1) then
1456 0 : write(*,2) 'report_cell_for_xm > xmstar', s% nz
1457 0 : return
1458 : end if
1459 0 : sumdq = 0
1460 0 : do k = 1, s% nz
1461 0 : sumdq = sumdq + s% dq(k)
1462 0 : if (sumdq >= dq) then
1463 0 : write(*,'(A)')
1464 0 : write(*,2) 'total mass in cells from 1 to k', k, sumdq*s% xmstar
1465 0 : write(*,2) 'logT(k)', k, s% lnT(k)/ln10
1466 0 : write(*,2) 'logRho(k)', k, s% lnd(k)/ln10
1467 0 : write(*,2) 'entropy(k)', k, exp(s% lnS(k))*amu/kerg
1468 0 : write(*,2) 'xmstar*q(k)', k, s% xmstar*s% q(k)
1469 0 : write(*,2) 'q(k)', k, s% q(k)
1470 0 : write(*,'(A)')
1471 0 : return
1472 : end if
1473 : end do
1474 0 : write(*,2) 'total mass in cells from 1 to nz', s% nz, s% xmstar
1475 : end subroutine do_report_cell_for_xm
1476 :
1477 2 : subroutine set_rate_factors(id, ierr)
1478 : use net_lib, only: get_net_reaction_table_ptr
1479 : use rates_lib, only: rates_reaction_id
1480 : integer, intent(in) :: id
1481 : integer, intent(out) :: ierr
1482 : type (star_info), pointer :: s
1483 : integer :: j, i, ir
1484 2 : integer, pointer :: net_reaction_ptr(:)
1485 : logical :: error
1486 :
1487 : include 'formats'
1488 :
1489 : ierr = 0
1490 2 : error = .false.
1491 2 : call star_ptr(id, s, ierr)
1492 2 : if (ierr /= 0) return
1493 :
1494 636 : s% rate_factors(:) = 1
1495 2 : if (s% job% num_special_rate_factors <= 0) return
1496 :
1497 : ! Dont error if we are changing net
1498 0 : if ((s% job% change_initial_net .or. s% job% change_net) .and. &
1499 : trim(s% job% new_net_name)/=trim(s% net_name)) then
1500 : !write(*,*) "Not changing special rates until net change"
1501 : return
1502 : end if
1503 :
1504 :
1505 0 : call get_net_reaction_table_ptr(s% net_handle, net_reaction_ptr, ierr)
1506 0 : if (ierr /= 0) return
1507 :
1508 0 : do i=1,s% job% num_special_rate_factors
1509 0 : if (len_trim(s% job% reaction_for_special_factor(i)) == 0) cycle
1510 0 : ir = rates_reaction_id(s% job% reaction_for_special_factor(i))
1511 0 : j = 0
1512 0 : if (ir > 0) j = net_reaction_ptr(ir)
1513 0 : if (j <= 0) then
1514 : write(*,*) 'Failed to find reaction_for_special_factor ' // &
1515 0 : trim(s% job% reaction_for_special_factor(i)), &
1516 0 : j, s% job% special_rate_factor(i)
1517 0 : error = .true.
1518 0 : cycle
1519 : end if
1520 0 : s% rate_factors(j) = s% job% special_rate_factor(i)
1521 : write(*,*) 'set special rate factor for ' // &
1522 0 : trim(s% job% reaction_for_special_factor(i)), &
1523 0 : j, s% job% special_rate_factor(i)
1524 : end do
1525 :
1526 0 : if(error) call mesa_error(__FILE__,__LINE__)
1527 :
1528 4 : end subroutine set_rate_factors
1529 :
1530 :
1531 1 : subroutine do_star_job_controls_before(id, s, restart, ierr)
1532 :
1533 2 : use rates_lib, only: rates_warning_init
1534 : use atm_support, only: get_atm_tau_base
1535 :
1536 : integer, intent(in) :: id
1537 : type (star_info), pointer :: s
1538 : logical, intent(in) :: restart
1539 : integer, intent(out) :: ierr
1540 : logical, parameter :: kap_use_cache = .true.
1541 : include 'formats'
1542 :
1543 : ierr = 0
1544 :
1545 1 : s% set_rate_factors => set_rate_factors ! will be called after net is defined
1546 :
1547 1 : call get_atm_tau_base(s, s% tau_base, ierr)
1548 1 : if (failed('atm_tau_base',ierr)) return
1549 :
1550 : call rates_warning_init( &
1551 1 : s% warn_rates_for_high_temp, s% max_safe_logT_for_rates)
1552 :
1553 1 : end subroutine do_star_job_controls_before
1554 :
1555 :
1556 4 : subroutine do_read_star_job_and_return_id(filename, id, ierr)
1557 : character(*), intent(in) :: filename
1558 : integer, intent(out) :: id
1559 : integer, intent(out) :: ierr
1560 : type (star_info), pointer :: s
1561 : character(len=strlen) :: inlist_fname
1562 :
1563 : include 'formats'
1564 1 : ierr = 0
1565 :
1566 1 : if (id_from_read_star_job /= 0) then
1567 0 : write(*,2) 'id_from_read_star_job', id_from_read_star_job
1568 0 : ierr = -1
1569 0 : return
1570 : end if
1571 :
1572 1 : call alloc_star(id, ierr)
1573 1 : if (ierr /= 0) then
1574 0 : write(*,*) 'do_read_star_job failed in alloc_star'
1575 0 : return
1576 : end if
1577 :
1578 1 : call star_ptr(id, s, ierr)
1579 1 : if (ierr /= 0) then
1580 0 : write(*,*) 'do_read_star_job failed in star_ptr'
1581 0 : return
1582 : end if
1583 :
1584 1 : call resolve_inlist_fname(inlist_fname,filename)
1585 1 : call read_star_job(s, inlist_fname, ierr)
1586 1 : if (ierr /= 0) then
1587 0 : write(*,*) 'ierr from read_star_job ' // trim(inlist_fname)
1588 0 : return
1589 : end if
1590 :
1591 1 : id_from_read_star_job = id
1592 :
1593 1 : if (s% job% save_star_job_namelist) then
1594 0 : call write_star_job(s, s% job% star_job_namelist_name, ierr)
1595 0 : if (ierr /= 0) then
1596 : write(*,*) 'ierr from write_star_job ' // &
1597 0 : trim(s% job% star_job_namelist_name)
1598 0 : return
1599 : end if
1600 : end if
1601 :
1602 1 : end subroutine do_read_star_job_and_return_id
1603 :
1604 : ! in a perfect world, we'd pass s as an arg to this routine.
1605 : ! but for backward compatibility for a large number of users
1606 : ! we do it this strange way instead.
1607 1 : subroutine do_read_star_job(filename, ierr)
1608 : character(*), intent(in) :: filename
1609 : integer, intent(out) :: ierr
1610 : integer :: id
1611 1 : call do_read_star_job_and_return_id(filename, id, ierr)
1612 1 : end subroutine do_read_star_job
1613 :
1614 :
1615 1 : subroutine do_load1_star(id, s, restart, restart_filename, ierr)
1616 : integer, intent(in) :: id
1617 : type (star_info), pointer :: s
1618 : logical, intent(in) :: restart
1619 : character (len=*), intent(in) :: restart_filename
1620 : integer, intent(out) :: ierr
1621 :
1622 1 : if (restart) then
1623 0 : call star_load_restart_photo(id, restart_filename, ierr)
1624 0 : if (failed('star_load_restart_photo',ierr)) return
1625 1 : else if (s% job% load_saved_photo) then
1626 0 : write(*,'(a)') 'load saved photo ' // trim(s% job% saved_photo_name)
1627 0 : write(*,'(A)')
1628 0 : call star_load_restart_photo(id, s% job% saved_photo_name, ierr)
1629 0 : if (failed('star_load_restart_photo',ierr)) return
1630 1 : else if (s% job% load_saved_model) then
1631 0 : if (s% job% create_merger_model) then
1632 0 : write(*,*) 'you have both load_saved_model and create_merger_model set true'
1633 0 : write(*,*) 'please pick one and try again'
1634 0 : call mesa_error(__FILE__,__LINE__)
1635 : end if
1636 0 : if (s% job% create_pre_main_sequence_model) then
1637 : write(*,*) 'you have both load_saved_model and ' // &
1638 0 : 'create_pre_main_sequence_model set true'
1639 0 : write(*,*) 'please pick one and try again'
1640 0 : call mesa_error(__FILE__,__LINE__)
1641 : end if
1642 0 : if (s% job% create_initial_model) then
1643 0 : write(*,*) 'you have both load_saved_model and create_initial_model set true'
1644 0 : write(*,*) 'please pick one and try again'
1645 0 : call mesa_error(__FILE__,__LINE__)
1646 : end if
1647 0 : write(*,'(a)') 'load saved model ' // trim(s% job% load_model_filename)
1648 0 : write(*,'(A)')
1649 0 : call star_read_model(id, s% job% load_model_filename, ierr)
1650 0 : if (failed('star_read_model',ierr)) return
1651 1 : else if (s% job% create_merger_model) then
1652 0 : call create_merger_model(s, ierr)
1653 0 : if (failed('create_merger_model',ierr)) return
1654 1 : else if (s% job% create_pre_main_sequence_model) then
1655 0 : if (.not. restart) write(*, *) 'create pre-main-sequence model'
1656 0 : if (s% job% create_initial_model) then
1657 : write(*,*) 'you have both create_pre_main_sequence_model ' // &
1658 0 : 'and create_initial_model set true'
1659 0 : write(*,*) 'please pick one and try again'
1660 0 : call mesa_error(__FILE__,__LINE__)
1661 : end if
1662 : call star_create_pre_ms_model( &
1663 : id, s% job% pre_ms_T_c, s% job% pre_ms_guess_rho_c, &
1664 : s% job% pre_ms_d_log10_P, s% job% pre_ms_logT_surf_limit, &
1665 : s% job% pre_ms_logP_surf_limit, s% job% initial_zfracs, &
1666 : s% job% dump_missing_metals_into_heaviest, &
1667 : (s% job% change_net .or. (s% job% change_initial_net .and. .not. restart)), &
1668 0 : s% job% new_net_name, s% job% pre_ms_relax_num_steps, ierr)
1669 0 : if (failed('star_create_pre_ms_model',ierr)) return
1670 1 : else if (s% job% create_RSP_model) then
1671 0 : if (.not. restart) write(*, *) 'create initial RSP model'
1672 0 : call star_create_RSP_model(id, ierr)
1673 0 : if (failed('star_create_RSP_model',ierr)) return
1674 1 : else if (s% job% create_RSP2_model) then
1675 0 : if (.not. restart) write(*, *) 'create initial RSP2 model'
1676 0 : call star_create_RSP2_model(id, ierr)
1677 0 : if (failed('star_create_RSP_model',ierr)) return
1678 1 : else if (s% job% create_initial_model) then
1679 0 : if (.not. restart) write(*, *) 'create initial model'
1680 0 : if (s% job% create_pre_main_sequence_model) then
1681 : write(*,*) 'you have both create_initial_model and ' // &
1682 0 : 'create_pre_main_sequence_model set true'
1683 0 : write(*,*) 'please pick one and try again'
1684 0 : call mesa_error(__FILE__,__LINE__)
1685 : end if
1686 : call star_create_initial_model(id, &
1687 : s% job% radius_in_cm_for_create_initial_model, &
1688 : s% job% mass_in_gm_for_create_initial_model, &
1689 : s% job% center_logP_1st_try_for_create_initial_model, &
1690 : s% job% entropy_1st_try_for_create_initial_model, &
1691 : s% job% max_tries_for_create_initial_model, &
1692 : s% job% abs_e01_tolerance_for_create_initial_model, &
1693 : s% job% abs_e02_tolerance_for_create_initial_model, &
1694 : s% job% initial_zfracs, &
1695 : s% job% dump_missing_metals_into_heaviest, &
1696 : (s% job% change_net .or. (s% job% change_initial_net .and. .not. restart)), &
1697 : s% job% new_net_name, s% job% initial_model_relax_num_steps, &
1698 : s% job% initial_model_eps, &
1699 0 : ierr)
1700 0 : if (failed('star_create_initial_model',ierr)) return
1701 : else
1702 1 : call star_load_zams(id, ierr)
1703 1 : if (failed('star_load_zams',ierr)) return
1704 : end if
1705 :
1706 : end subroutine do_load1_star
1707 :
1708 :
1709 0 : subroutine create_merger_model(s, ierr)
1710 : use ctrls_io, only : store_controls
1711 : type (star_info), pointer :: s
1712 : integer, intent(out) :: ierr
1713 :
1714 : integer :: id, id_aux, i, j, k
1715 : type (star_info), pointer :: s_aux
1716 0 : real(dp), pointer :: xq(:), xa(:,:)
1717 : real(dp) :: total_mass, partial_mass
1718 :
1719 : include 'formats'
1720 0 : ierr = 0
1721 0 : id = s% id
1722 :
1723 0 : if (s% job% create_pre_main_sequence_model) then
1724 : write(*,*) 'you have both load_saved_model and ' // &
1725 0 : 'create_pre_main_sequence_model set true'
1726 0 : write(*,*) 'please pick one and try again'
1727 0 : call mesa_error(__FILE__,__LINE__)
1728 : end if
1729 0 : if (s% job% create_initial_model) then
1730 0 : write(*,*) 'you have both load_saved_model and create_initial_model set true'
1731 0 : write(*,*) 'please pick one and try again'
1732 0 : call mesa_error(__FILE__,__LINE__)
1733 : end if
1734 : !load first star
1735 0 : call star_read_model(id, s% job% saved_model_for_merger_1, ierr)
1736 0 : if (failed('star_read_model',ierr)) return
1737 :
1738 : !load second star
1739 0 : call alloc_star(id_aux, ierr)
1740 0 : if (failed('alloc_star',ierr)) return
1741 0 : call star_ptr(id_aux, s_aux, ierr)
1742 0 : if (failed('star_ptr',ierr)) return
1743 0 : call init_starting_star_data(s_aux, ierr)
1744 0 : if (failed('init_starting_star_data',ierr)) return
1745 0 : call star_set_kap_and_eos_handles(id_aux, ierr)
1746 0 : if (failed('set_star_kap_and_eos_handles',ierr)) return
1747 0 : call star_set_colors_handles(id_aux, ierr)
1748 0 : if (failed('star_set_colors_handles',ierr)) return
1749 0 : call store_controls(s_aux, ierr)
1750 0 : if (failed('store_controls',ierr)) return
1751 0 : call do_star_job_controls_before(id_aux, s_aux, .false., ierr)
1752 0 : if (ierr /= 0) return
1753 0 : call star_read_model(id_aux, s% job% saved_model_for_merger_2, ierr)
1754 0 : if (failed('star_read_model',ierr)) return
1755 :
1756 : ! create composition and q array through an entropy sorting
1757 0 : total_mass = s% mstar + s_aux% mstar
1758 0 : partial_mass = 0
1759 0 : i = 1
1760 0 : j = 1
1761 0 : allocate(xq(s% nz + s_aux% nz), xa(s% species, s% nz + s_aux% nz))
1762 0 : do while (i <= s% nz .or. j <= s_aux% nz)
1763 0 : if (j > s_aux% nz .or. (i <= s% nz .and. &
1764 : s% entropy(i) >= s_aux% entropy(j))) then
1765 0 : partial_mass = partial_mass + s% dm(i)
1766 0 : do k=1, s% species
1767 0 : xa(k, i+j-1) = s% xa(k, i)
1768 : end do
1769 0 : i = i + 1
1770 0 : else if (i > s% nz .or. (j <= s_aux% nz .and. &
1771 : s_aux% entropy(j) > s% entropy(i))) then
1772 0 : partial_mass = partial_mass + s_aux% dm(j)
1773 0 : do k=1, s% species
1774 0 : xa(k, i+j-1) = s_aux% xa(k, j)
1775 : end do
1776 0 : j = j + 1
1777 : end if
1778 0 : xq(i+j-2) = partial_mass / total_mass
1779 : !write(*,*) "check", i+j-2, xq(i+j-2), xa(1, i+j-2), xa(2, i+j-2), xa(3, i+j-2)
1780 : end do
1781 : ! Relax composition first, then composition mass
1782 : ! Turn off rotation for relaxation
1783 0 : call star_set_rotation_flag(id, .false., ierr)
1784 0 : if (failed('star_set_rotation_flag',ierr)) then
1785 0 : deallocate(xq,xa)
1786 0 : return
1787 : end if
1788 0 : write(*,*) "Relaxing composition to merger composition"
1789 : call star_relax_composition( &
1790 0 : id, s% job% num_steps_to_relax_composition, s% nz + s_aux% nz, s% species, xa, xq, ierr)
1791 0 : if (failed('star_relax_composition',ierr)) then
1792 0 : deallocate(xq,xa)
1793 0 : return
1794 : end if
1795 0 : write(*,*) "Relaxing star mass to total merger mass"
1796 : call star_relax_mass_scale( &
1797 : id, total_mass/Msun, s% job% dlgm_per_step, &
1798 0 : s% job% change_mass_years_for_dt, ierr)
1799 0 : deallocate(xq,xa)
1800 0 : if (failed('star_relax_mass_scale',ierr)) return
1801 :
1802 0 : end subroutine create_merger_model
1803 :
1804 :
1805 11 : subroutine extend_net(s, ierr)
1806 0 : use net_def
1807 : use chem_def
1808 : type (star_info), pointer :: s
1809 : integer, intent(out) :: ierr
1810 : real(dp), parameter :: tiny = 1d-10, small = 1d-2
1811 :
1812 : real(dp) :: cntr_h, cntr_he
1813 :
1814 : include 'formats'
1815 :
1816 11 : ierr = 0
1817 :
1818 : !write(*,2) 'extend_net: current net ' // trim(s% net_name), s% model_number
1819 :
1820 11 : if (s% net_name == s% job% adv_net) return
1821 :
1822 11 : if (s% net_name == s% job% co_net) then
1823 0 : if (s% log_max_temperature > 9d0 .or. s% log_center_density > 9d0) then
1824 0 : call change_net(s% job% adv_net)
1825 0 : if (len_trim(s% job% profile_columns_file) > 0) &
1826 0 : write(*,*) 'read ' // trim(s% job% profile_columns_file)
1827 : call star_set_profile_columns( &
1828 0 : s% id, s% job% profile_columns_file, .true., ierr)
1829 : end if
1830 0 : return
1831 : end if
1832 :
1833 22 : if (s% net_name == s% job% h_he_net) then
1834 11 : cntr_h = current_abundance_at_point(s% id, ih1, s% nz, ierr)
1835 : !write(*,2) 'cntr_h', s% model_number, cntr_h, tiny
1836 11 : if (ierr /= 0) return
1837 11 : if (cntr_h > tiny) return
1838 0 : cntr_he = current_abundance_at_point(s% id, ihe4, s% nz, ierr)
1839 : !write(*,2) 'cntr_he', s% model_number, cntr_he, small
1840 0 : if (ierr /= 0) return
1841 0 : if (cntr_he > small) return
1842 0 : if (s% log_max_temperature > 8.3d0 .or. s% log_center_density > 8.5d0) then
1843 0 : call change_net(s% job% co_net)
1844 0 : if (len_trim(s% job% profile_columns_file) > 0) &
1845 0 : write(*,*) 'read ' // trim(s% job% profile_columns_file)
1846 : call star_set_profile_columns( &
1847 0 : s% id, s% job% profile_columns_file, .true., ierr)
1848 : end if
1849 : end if
1850 :
1851 :
1852 : contains
1853 :
1854 :
1855 0 : subroutine change_net(net_name)
1856 11 : use const_def, only: dp
1857 : character (len=*), intent(in) :: net_name
1858 :
1859 : include 'formats'
1860 :
1861 : call star_change_to_new_net( &
1862 0 : s% id, s% job% adjust_abundances_for_new_isos, net_name, ierr)
1863 0 : if (ierr /= 0) then
1864 0 : write(*,*) 'failed in star_change_to_new_net ' // trim(net_name)
1865 0 : call mesa_error(__FILE__,__LINE__,'change_net')
1866 : return
1867 : end if
1868 :
1869 0 : if (net_name /= s% net_name) then
1870 0 : write(*,*) ' new net_name ', trim(net_name)
1871 0 : write(*,*) 'old s% net_name ', trim(s% net_name)
1872 0 : write(*,*) 'failed to change'
1873 0 : call mesa_error(__FILE__,__LINE__,'change_net')
1874 : end if
1875 :
1876 0 : write(*,'(a)') ' new net = ' // trim(s% net_name)
1877 : !do j=1,s% species
1878 : ! write(*,fmt='(a,x)',advance='no') trim(chem_isos% name(s% chem_id(j)))
1879 : !end do
1880 : !write(*,*)
1881 0 : s% dt_next = s% dt_next/5
1882 : !write(*,1) 'reduce timestep', log10(s% dt_next/secyer)
1883 0 : write(*,'(A)')
1884 : end subroutine change_net
1885 :
1886 :
1887 : end subroutine extend_net
1888 :
1889 :
1890 1 : subroutine before_evolve(id, ierr)
1891 : integer, intent(in) :: id
1892 : integer, intent(out) :: ierr
1893 1 : ierr = 0
1894 : end subroutine before_evolve
1895 :
1896 :
1897 1 : subroutine do_star_job_controls_after(id, s, restart, pgstar_ok, ierr)
1898 : use const_def, only: dp
1899 : use rates_def
1900 : use rates_lib
1901 : use utils_lib, only: utils_OMP_GET_MAX_THREADS
1902 :
1903 : integer, intent(in) :: id
1904 : type (star_info), pointer :: s
1905 : logical, intent(in) :: restart, pgstar_ok
1906 : integer, intent(out) :: ierr
1907 :
1908 : real(dp) :: log_m, log_lifetime, max_dt, max_timestep
1909 : integer :: i, j, nzlo, nzhi, chem_id, chem_id1, chem_id2
1910 : logical :: change_v, change_u
1911 : include 'formats'
1912 :
1913 1 : if (s% job% change_net .or. (s% job% change_initial_net .and. .not. restart)) then
1914 : call star_change_to_new_net( &
1915 0 : id, s% job% adjust_abundances_for_new_isos, s% job% new_net_name, ierr)
1916 0 : if (failed('star_change_to_new_net',ierr)) return
1917 : end if
1918 :
1919 1 : if (s% job% change_small_net .or. &
1920 : (s% job% change_initial_small_net .and. .not. restart)) then
1921 0 : write(*,*) 'change small net to ' // trim(s% job% new_small_net_name)
1922 : call star_change_to_new_small_net( &
1923 0 : id, s% job% adjust_abundances_for_new_isos, s% job% new_small_net_name, ierr)
1924 0 : if (failed('star_change_to_new_small_net',ierr)) return
1925 0 : write(*,*) 'number of species', s% species
1926 : end if
1927 :
1928 :
1929 1 : if (len_trim(s% job% history_columns_file) > 0) &
1930 0 : write(*,*) 'read ' // trim(s% job% history_columns_file)
1931 1 : call star_set_history_columns(id, s% job% history_columns_file, .true., ierr)
1932 1 : if (failed('star_set_history_columns',ierr)) return
1933 :
1934 1 : if (len_trim(s% job% profile_columns_file) > 0) &
1935 0 : write(*,*) 'read ' // trim(s% job% profile_columns_file)
1936 1 : call star_set_profile_columns(id, s% job% profile_columns_file, .true., ierr)
1937 1 : if (failed('star_set_profile_columns',ierr)) return
1938 :
1939 1 : if (pgstar_ok) then
1940 1 : if (s% job% clear_pgstar_history .or. &
1941 : (s% job% clear_initial_pgstar_history .and. .not. restart)) then
1942 1 : call start_new_run_for_pgstar(s, ierr)
1943 1 : if (failed('start_new_run_for_pgstar',ierr)) return
1944 : else
1945 0 : call restart_run_for_pgstar(s, ierr)
1946 0 : if (failed('restart_run_for_pgstar',ierr)) return
1947 : end if
1948 : end if
1949 :
1950 1 : if (s% job% set_tau_factor .or. &
1951 : (s% job% set_initial_tau_factor .and. .not. restart)) then
1952 0 : write(*,1) 'set_tau_factor', s% job% set_to_this_tau_factor
1953 0 : s% tau_factor = s% job% set_to_this_tau_factor
1954 : end if
1955 :
1956 1 : if (s% job% set_Tsurf_factor .or. &
1957 : (s% job% set_initial_Tsurf_factor .and. .not. restart)) then
1958 0 : write(*,1) 'set_Tsurf_factor', s% job% set_to_this_Tsurf_factor
1959 0 : s% Tsurf_factor = s% job% set_to_this_Tsurf_factor
1960 : end if
1961 :
1962 1 : if (s% job% set_initial_age .and. .not. restart) then
1963 0 : write(*,1) 'set_initial_age', s% job% initial_age ! in years
1964 0 : call star_set_age(id, s% job% initial_age, ierr)
1965 0 : if (failed('star_set_age',ierr)) return
1966 : end if
1967 :
1968 1 : if (s% job% set_initial_dt .and. .not. restart) then
1969 0 : if (s% job% years_for_initial_dt > 0d0) then
1970 0 : write(*,1) 'set_initial_dt (years)', s% job% years_for_initial_dt
1971 0 : s% dt_next = s% job% years_for_initial_dt*secyer
1972 0 : else if (s% job% seconds_for_initial_dt > 0d0) then
1973 0 : write(*,1) 'set_initial_dt (seconds)', s% job% seconds_for_initial_dt
1974 0 : s% dt_next = s% job% seconds_for_initial_dt
1975 : end if
1976 : end if
1977 :
1978 1 : if (s% job% limit_initial_dt .and. .not. restart) then
1979 0 : if (s% job% years_for_initial_dt > 0d0) then
1980 0 : write(*,1) 'limit_initial_dt (years)', s% job% years_for_initial_dt
1981 0 : s% dt_next = min(s% dt_next, s% job% years_for_initial_dt*secyer)
1982 0 : else if (s% job% seconds_for_initial_dt > 0d0) then
1983 0 : write(*,1) 'limit_initial_dt (seconds)', s% job% seconds_for_initial_dt
1984 0 : s% dt_next = min(s% dt_next, s% job% seconds_for_initial_dt)
1985 : end if
1986 :
1987 : end if
1988 :
1989 : ! enforce max_timestep on first step
1990 :
1991 1 : if (s% max_years_for_timestep > 0) then
1992 0 : max_timestep = secyer*s% max_years_for_timestep
1993 0 : if (s% max_timestep > 0 .and. s% max_timestep < max_timestep) &
1994 0 : max_timestep = s% max_timestep
1995 : else
1996 1 : max_timestep = s% max_timestep
1997 : end if
1998 :
1999 1 : if (max_timestep > 0 .and. max_timestep < s% dt_next) then
2000 0 : write(*,1) 'max_timestep (seconds)', max_timestep
2001 0 : s% dt_next = max_timestep
2002 : end if
2003 :
2004 1 : if (s% job% set_initial_model_number .and. .not. restart) then
2005 0 : write(*,2) 'set_initial_model_number', s% job% initial_model_number
2006 0 : s% model_number = s% job% initial_model_number
2007 0 : s% init_model_number = s% model_number
2008 : end if
2009 :
2010 1 : if (s% job% set_initial_number_retries .and. .not. restart) then
2011 1 : write(*,2) 'set_initial_number_retries', s% job% initial_number_retries
2012 1 : s% num_retries = s% job% initial_number_retries
2013 : end if
2014 :
2015 1 : if (s% job% steps_to_take_before_terminate >= 0) then
2016 0 : s% max_model_number = s% model_number + s% job% steps_to_take_before_terminate
2017 0 : write(*,2) 'steps_to_take_before_terminate', &
2018 0 : s% job% steps_to_take_before_terminate
2019 0 : write(*,2) 'max_model_number', s% max_model_number
2020 : end if
2021 :
2022 1 : if (s% job% steps_before_start_timing > 0) then
2023 0 : s% job% first_model_for_timing = s% model_number + s% job% steps_before_start_timing
2024 0 : write(*,2) 'steps_before_start_timing', &
2025 0 : s% job% steps_before_start_timing
2026 : end if
2027 :
2028 1 : if (abs(s% job% T9_weaklib_full_off - T9_weaklib_full_off) > 1d-6) then
2029 0 : write(*,1) 'set T9_weaklib_full_off', s% job% T9_weaklib_full_off
2030 0 : T9_weaklib_full_off = s% job% T9_weaklib_full_off
2031 : end if
2032 :
2033 1 : if (abs(s% job% T9_weaklib_full_on - T9_weaklib_full_on) > 1d-6) then
2034 0 : write(*,1) 'set T9_weaklib_full_on', s% job% T9_weaklib_full_on
2035 0 : T9_weaklib_full_on = s% job% T9_weaklib_full_on
2036 : end if
2037 :
2038 1 : if (s% job% weaklib_blend_hi_Z /= weaklib_blend_hi_Z) then
2039 0 : write(*,1) 'set weaklib_blend_hi_Z', s% job% weaklib_blend_hi_Z
2040 0 : weaklib_blend_hi_Z = s% job% weaklib_blend_hi_Z
2041 : end if
2042 :
2043 1 : if (abs(s% job% T9_weaklib_full_off_hi_Z - T9_weaklib_full_off_hi_Z) > 1d-6) then
2044 0 : write(*,1) 'set T9_weaklib_full_off_hi_Z', s% job% T9_weaklib_full_off_hi_Z
2045 0 : T9_weaklib_full_off_hi_Z = s% job% T9_weaklib_full_off_hi_Z
2046 : end if
2047 :
2048 1 : if (abs(s% job% T9_weaklib_full_on_hi_Z - T9_weaklib_full_on_hi_Z) > 1d-6) then
2049 0 : write(*,1) 'set T9_weaklib_full_on_hi_Z', s% job% T9_weaklib_full_on_hi_Z
2050 0 : T9_weaklib_full_on_hi_Z = s% job% T9_weaklib_full_on_hi_Z
2051 : end if
2052 :
2053 : ! set up coulomb corrections for the special weak rates
2054 1 : which_mui_coulomb = get_mui_value(s% job% ion_coulomb_corrections)
2055 1 : which_vs_coulomb = get_vs_value(s% job% electron_coulomb_corrections)
2056 :
2057 : change_v = s% job% change_v_flag .or. &
2058 1 : (s% job% change_initial_v_flag .and. .not. restart)
2059 : change_u = s% job% change_u_flag .or. &
2060 1 : (s% job% change_initial_u_flag .and. .not. restart)
2061 1 : if (change_v .or. change_u) then
2062 : ! do add new before remove old so can set initial values
2063 0 : if (change_v .and. s% job% new_v_flag) then
2064 0 : write(*,*) 'new_v_flag', s% job% new_v_flag
2065 0 : call star_set_v_flag(id, s% job% new_v_flag, ierr)
2066 0 : if (failed('star_set_v_flag',ierr)) return
2067 : end if
2068 0 : if (change_u .and. s% job% new_u_flag) then
2069 0 : write(*,*) 'new_u_flag', s% job% new_u_flag
2070 0 : call star_set_u_flag(id, s% job% new_u_flag, ierr)
2071 0 : if (failed('star_set_u_flag',ierr)) return
2072 : end if
2073 0 : if (change_v .and. .not. s% job% new_v_flag) then
2074 0 : write(*,*) 'new_v_flag', s% job% new_v_flag
2075 0 : call star_set_v_flag(id, s% job% new_v_flag, ierr)
2076 0 : if (failed('star_set_v_flag',ierr)) return
2077 : end if
2078 0 : if (change_u .and. .not. s% job% new_u_flag) then
2079 0 : write(*,*) 'new_u_flag', s% job% new_u_flag
2080 0 : call star_set_u_flag(id, s% job% new_u_flag, ierr)
2081 0 : if (failed('star_set_u_flag',ierr)) return
2082 : end if
2083 : end if
2084 :
2085 1 : if (s% job% change_RTI_flag .or. &
2086 : (s% job% change_initial_RTI_flag .and. .not. restart)) then
2087 0 : write(*,*) 'new_RTI_flag', s% job% new_RTI_flag
2088 0 : call star_set_RTI_flag(id, s% job% new_RTI_flag, ierr)
2089 0 : if (failed('star_set_RTI_flag',ierr)) return
2090 : end if
2091 :
2092 1 : if (s% job% change_RSP2_flag .or. &
2093 : (s% job% change_initial_RSP2_flag .and. .not. restart)) then
2094 0 : write(*,*) 'new_RSP2_flag', s% job% new_RSP2_flag
2095 0 : call star_set_RSP2_flag(id, s% job% new_RSP2_flag, ierr)
2096 0 : if (failed('star_set_RSP2_flag',ierr)) return
2097 : end if
2098 :
2099 1 : if (s% job% change_RSP_flag .or. &
2100 : (s% job% change_initial_RSP_flag .and. .not. restart)) then
2101 0 : write(*,*) 'new_RSP_flag', s% job% new_RSP_flag
2102 0 : call star_set_RSP_flag(id, s% job% new_RSP_flag, ierr)
2103 0 : if (failed('star_set_RSP_flag',ierr)) return
2104 : end if
2105 :
2106 1 : if (s% job% change_w_div_wc_flag .or. &
2107 : (s% job% change_initial_w_div_wc_flag .and. .not. restart)) then
2108 0 : write(*,*) 'new_w_div_wc_flag', s% job% new_w_div_wc_flag
2109 0 : call star_set_w_div_wc_flag(id, s% job% new_w_div_wc_flag, ierr)
2110 0 : if (failed('star_set_w_div_wc_flag',ierr)) return
2111 : end if
2112 :
2113 1 : if (s% job% change_j_rot_flag .or. &
2114 : (s% job% change_initial_j_rot_flag .and. .not. restart)) then
2115 0 : write(*,*) 'new_j_rot_flag', s% job% new_j_rot_flag
2116 0 : call star_set_j_rot_flag(id, s% job% new_j_rot_flag, ierr)
2117 0 : if (failed('star_set_j_rot_flag',ierr)) return
2118 : end if
2119 :
2120 1 : if (s% job% change_D_omega_flag .or. &
2121 : (s% job% change_initial_D_omega_flag .and. .not. restart)) then
2122 0 : call star_set_D_omega_flag(id, s% job% new_D_omega_flag, ierr)
2123 0 : if (failed('star_set_D_omega_flag',ierr)) return
2124 : end if
2125 :
2126 1 : if (s% job% change_am_nu_rot_flag .or. &
2127 : (s% job% change_initial_am_nu_rot_flag .and. .not. restart)) then
2128 0 : call star_set_am_nu_rot_flag(id, s% job% new_am_nu_rot_flag, ierr)
2129 0 : if (failed('star_set_am_nu_rot_flag',ierr)) return
2130 : end if
2131 :
2132 1 : if (s% job% change_rotation_flag .or. &
2133 : (s% job% change_initial_rotation_flag .and. .not. restart)) then
2134 0 : write(*,*) 'new_rotation_flag', s% job% new_rotation_flag
2135 0 : call star_set_rotation_flag(id, s% job% new_rotation_flag, ierr)
2136 0 : if (failed('star_set_rotation_flag',ierr)) return
2137 : end if
2138 :
2139 1 : if (s% rotation_flag .and. s% job% set_omega) then
2140 0 : write(*,1) 'new_omega', s% job% new_omega
2141 0 : call star_set_uniform_omega(id, s% job% new_omega, ierr)
2142 0 : if (failed('star_set_uniform_omega',ierr)) return
2143 : end if
2144 :
2145 1 : if (s% rotation_flag .and. s% job% set_initial_omega .and. .not. restart) then
2146 0 : write(*,1) 'new_omega', s% job% new_omega
2147 0 : call star_set_uniform_omega(id, s% job% new_omega, ierr)
2148 0 : if (failed('star_set_uniform_omega',ierr)) return
2149 : end if
2150 :
2151 1 : if (s% rotation_flag .and. s% job% set_surface_rotation_v) then
2152 0 : s% job% new_omega = s% job% new_surface_rotation_v*1d5/s% r(1)
2153 0 : write(*,1) 'new_surface_rotation_v', &
2154 0 : s% job% new_surface_rotation_v, s% job% new_omega
2155 0 : call star_set_uniform_omega(id, s% job% new_omega, ierr)
2156 0 : if (failed('star_set_uniform_omega',ierr)) return
2157 : end if
2158 :
2159 : if (s% rotation_flag .and. &
2160 1 : s% job% set_initial_surface_rotation_v .and. .not. restart) then
2161 0 : s% job% new_omega = s% job% new_surface_rotation_v*1d5/s% r(1)
2162 0 : write(*,2) 'new_surface_rotation_v', &
2163 0 : s% model_number, s% job% new_surface_rotation_v, s% job% new_omega
2164 0 : call star_set_uniform_omega(id, s% job% new_omega, ierr)
2165 0 : if (failed('star_set_uniform_omega',ierr)) return
2166 : end if
2167 :
2168 1 : if (s% rotation_flag .and. s% job% set_omega_div_omega_crit) then
2169 : s% job% new_omega = &
2170 0 : s% job% new_omega_div_omega_crit*star_surface_omega_crit(id, ierr)
2171 0 : if (failed('star_surface_omega_crit',ierr)) return
2172 0 : write(*,2) 'new_omega_div_omega_crit', &
2173 0 : s% model_number, s% job% new_omega_div_omega_crit, s% job% new_omega
2174 0 : call star_set_uniform_omega(id, s% job% new_omega, ierr)
2175 0 : if (failed('star_set_uniform_omega',ierr)) return
2176 : end if
2177 :
2178 : if (s% rotation_flag .and. &
2179 1 : s% job% set_initial_omega_div_omega_crit .and. .not. restart) then
2180 : s% job% new_omega = &
2181 0 : s% job% new_omega_div_omega_crit*star_surface_omega_crit(id, ierr)
2182 0 : if (failed('star_surface_omega_crit',ierr)) return
2183 0 : write(*,2) 'new_omega_div_omega_crit', &
2184 0 : s% model_number, s% job% new_omega_div_omega_crit, s% job% new_omega
2185 0 : call star_set_uniform_omega(id, s% job% new_omega, ierr)
2186 0 : if (failed('star_set_uniform_omega',ierr)) return
2187 : end if
2188 :
2189 1 : if (s% job% set_to_xa_for_accretion .or. &
2190 : (s% job% set_initial_to_xa_for_accretion .and. .not. restart)) then
2191 0 : write(*,*) 'set_to_xa_for_accretion'
2192 0 : call change_to_xa_for_accretion(id, s% job% set_nzlo, s% job% set_nzhi, ierr)
2193 0 : if (failed('set_to_xa_for_accretion',ierr)) return
2194 : end if
2195 :
2196 1 : if (s% job% first_model_for_timing > 0) &
2197 0 : write(*,2) 'first_model_for_timing', s% job% first_model_for_timing
2198 :
2199 1 : if (s% job% set_uniform_initial_composition .and. .not. restart) then
2200 0 : write(*,'(A)')
2201 0 : write(*,1) 'set_uniform_initial_composition'
2202 0 : write(*,1) 'initial_h1', s% job% initial_h1
2203 0 : write(*,1) 'initial_h2', s% job% initial_h2
2204 0 : write(*,1) 'initial_he3', s% job% initial_he3
2205 0 : write(*,1) 'initial_he4', s% job% initial_he4
2206 0 : select case(s% job% initial_zfracs)
2207 : case (AG89_zfracs)
2208 0 : write(*,1) 'metals AG89'
2209 : case (GN93_zfracs)
2210 0 : write(*,1) 'metals GN93'
2211 : case (GS98_zfracs)
2212 0 : write(*,1) 'metals GS98'
2213 : case (L03_zfracs)
2214 0 : write(*,1) 'metals L03'
2215 : case (AGS05_zfracs)
2216 0 : write(*,1) 'metals AGS05'
2217 : case (AGSS09_zfracs)
2218 0 : write(*,1) 'metals AGSS09'
2219 : case (L09_zfracs)
2220 0 : write(*,1) 'metals L09'
2221 : case (A09_Prz_zfracs)
2222 0 : write(*,1) 'metals A09_Prz'
2223 : case (MB22_photospheric_zfracs)
2224 0 : write(*,1) 'metals MB22_photospheric'
2225 : case (AAG21_photospheric_zfracs)
2226 0 : write(*,1) 'metals AAG21_photospheric'
2227 : case default
2228 0 : write(*,2) 'unknown value for initial_zfracs', s% job% initial_zfracs
2229 : end select
2230 : call star_set_standard_composition( &
2231 : id, s% job% initial_h1, s% job% initial_h2, &
2232 : s% job% initial_he3, s% job% initial_he4, s% job% initial_zfracs, &
2233 0 : s% job% dump_missing_metals_into_heaviest, ierr)
2234 0 : if (failed('set_uniform_initial_composition',ierr)) return
2235 : end if
2236 :
2237 1 : if (s% job% relax_initial_composition .and. .not. restart) then
2238 0 : call do_relax_initial_composition(ierr)
2239 0 : if (failed('do_relax_initial_composition',ierr)) return
2240 : end if
2241 :
2242 1 : if (s% job% relax_initial_to_xaccrete .and. .not. restart) then
2243 0 : call star_relax_to_xaccrete(id, s% job% num_steps_to_relax_composition, ierr)
2244 0 : if (failed('star_relax_to_xaccrete',ierr)) return
2245 : end if
2246 :
2247 1 : if (s% job% set_uniform_xa_from_file) then
2248 0 : call star_uniform_xa_from_file(id, s% job% file_for_uniform_xa, ierr)
2249 0 : if (failed('star_uniform_xa_from_file',ierr)) return
2250 : end if
2251 :
2252 1 : if (s% job% relax_initial_angular_momentum .and. .not. restart) then
2253 0 : call do_relax_initial_angular_momentum(ierr)
2254 0 : if (failed('do_relax_initial_angular_momentum',ierr)) return
2255 : end if
2256 :
2257 1 : if (s% job% relax_initial_entropy .and. .not. restart) then
2258 0 : call do_relax_initial_entropy(ierr)
2259 0 : if (failed('do_relax_initial_entropy',ierr)) return
2260 : end if
2261 :
2262 1 : if (s% job% mix_section .or. &
2263 : (s% job% mix_initial_section .and. .not. restart)) then
2264 0 : write(*,*) 'mix_section'
2265 : call uniform_mix_section( &
2266 0 : id, s% job% mix_section_nzlo, s% job% mix_section_nzhi, ierr)
2267 0 : if (failed('uniform_mix_section',ierr)) return
2268 : end if
2269 :
2270 1 : if (s% job% mix_initial_envelope_down_to_T > 0d0 .and. .not. restart) then
2271 0 : call uniform_mix_envelope_down_to_T(id, s% job% mix_initial_envelope_down_to_T, ierr)
2272 0 : if (failed('uniform_mix_envelope_down_to_T',ierr)) return
2273 : end if
2274 :
2275 1 : if (s% job% mix_envelope_down_to_T > 0d0) then
2276 0 : call uniform_mix_envelope_down_to_T(id, s% job% mix_envelope_down_to_T, ierr)
2277 0 : if (failed('uniform_mix_envelope_down_to_T',ierr)) return
2278 : end if
2279 :
2280 1 : if (s% job% mix_initial_envelope_down_to_T > 0d0) then
2281 0 : call uniform_mix_envelope_down_to_T(id, s% job% mix_initial_envelope_down_to_T, ierr)
2282 0 : if (failed('uniform_mix_envelope_down_to_T',ierr)) return
2283 : end if
2284 :
2285 1 : if (s% job% set_uniform_initial_xa_from_file .and. .not. restart) then
2286 0 : call star_uniform_xa_from_file(id, s% job% file_for_uniform_xa, ierr)
2287 0 : if (failed('star_uniform_xa_from_file',ierr)) return
2288 : end if
2289 :
2290 : ! do change Z before change Y since changing Z can change Y
2291 1 : if (s% job% change_Z) then
2292 0 : call star_set_z(id, s% job% new_Z, ierr)
2293 0 : if (failed('star_set_z',ierr)) return
2294 : end if
2295 :
2296 1 : if (s% job% change_initial_Z .and. .not. restart) then
2297 0 : call star_set_z(id, s% job% new_Z, ierr)
2298 0 : if (failed('star_set_z',ierr)) return
2299 : end if
2300 :
2301 1 : if (s% job% change_Y) then
2302 0 : call star_set_y(id, s% job% new_Y, ierr)
2303 0 : if (failed('change_Y',ierr)) return
2304 : end if
2305 :
2306 1 : if (s% job% change_initial_Y .and. .not. restart) then
2307 0 : call star_set_y(id, s% job% new_Y, ierr)
2308 0 : if (failed('change_initial_Y',ierr)) return
2309 : end if
2310 :
2311 1 : if (s% job% zero_alpha_RTI .or. &
2312 : (s% job% zero_initial_alpha_RTI .and. .not. restart)) then
2313 0 : call star_zero_alpha_RTI(id, ierr)
2314 0 : if (failed('star_zero_alpha_RTI',ierr)) return
2315 : end if
2316 :
2317 1 : if (s% job% set_abundance .or. &
2318 : (s% job% set_initial_abundance .and. .not. restart)) then
2319 0 : nzlo = s% job% set_abundance_nzlo
2320 0 : nzhi = s% job% set_abundance_nzhi
2321 0 : if (nzhi <= 0) nzhi = s% nz
2322 0 : if (nzlo <= 0) nzlo = 1
2323 0 : write(*, *) 'set_abundance of ', &
2324 0 : trim(s% job% chem_name), s% job% new_frac, nzlo, nzhi
2325 0 : chem_id = get_nuclide_index(s% job% chem_name)
2326 0 : if (chem_id <= 0) then
2327 0 : write(*,*) 'failed to find ' // trim(s% job% chem_name)
2328 0 : write(*,*) 'check valid chem_isos% names in chem/public/chem_def.f'
2329 : end if
2330 0 : call set_abundance_in_section(id, chem_id, s% job% new_frac, nzlo, nzhi, ierr)
2331 0 : if (failed('set_abundance_in_section',ierr)) return
2332 : end if
2333 :
2334 1 : if (s% job% replace_element .or. &
2335 : (s% job% replace_initial_element .and. .not. restart)) then
2336 0 : write(*, *) 'replace_element ', &
2337 0 : trim(s% job% chem_name1), ' by ', trim(s% job% chem_name2)
2338 0 : chem_id1 = get_nuclide_index(s% job% chem_name1)
2339 0 : chem_id2 = get_nuclide_index(s% job% chem_name2)
2340 0 : if (chem_id1 <= 0) then
2341 0 : write(*,*) 'failed to find ' // trim(s% job% chem_name1)
2342 0 : write(*,*) 'check valid chem_isos% names in chem/public/chem_def.f'
2343 : end if
2344 0 : if (chem_id2 <= 0) then
2345 0 : write(*,*) 'failed to find ' // trim(s% job% chem_name2)
2346 0 : write(*,*) 'check valid chem_isos% names in chem/public/chem_def.f'
2347 : end if
2348 0 : nzhi = s% job% replace_element_nzhi
2349 0 : nzlo = s% job% replace_element_nzlo
2350 0 : if (nzhi <= 0) nzhi = s% nz
2351 0 : if (nzlo <= 0) nzlo = 1
2352 0 : write(*, *) 'in section', nzlo, nzhi
2353 : call replace_element_in_section( &
2354 0 : id, chem_id1, chem_id2, nzlo, nzhi, ierr)
2355 0 : if (failed('replace_element_in_section',ierr)) return
2356 : end if
2357 :
2358 1 : if (s% job% set_irradiation .or. &
2359 : (s% job% set_initial_irradiation .and. .not. restart)) then
2360 0 : write(*,2) 'set_irradiation'
2361 0 : s% irradiation_flux = s% job% set_to_this_irrad_flux
2362 0 : s% column_depth_for_irradiation = s% job% irrad_col_depth
2363 : end if
2364 :
2365 1 : if (s% job% do_special_test) then
2366 0 : write(*, *) 'do_special_test'
2367 0 : call star_special_test(id, ierr)
2368 0 : if (failed('star_special_test',ierr)) return
2369 : end if
2370 :
2371 1 : if (s% job% set_v_center .or. &
2372 : (s% job% set_initial_v_center .and. .not. restart)) then
2373 0 : write(*, 1) 'set_v_center', s% job% new_v_center
2374 0 : s% v_center = s% job% new_v_center
2375 : end if
2376 :
2377 1 : if (s% job% set_L_center .or. &
2378 : (s% job% set_initial_L_center .and. .not. restart)) then
2379 0 : write(*, 1) 'set_L_center', s% job% new_L_center
2380 0 : s% L_center = s% job% new_L_center*Lsun
2381 : end if
2382 :
2383 : ! do "set" before "relax"
2384 :
2385 : ! must do relax Z before relax Y since relax Z can change Y
2386 : ! (Warrick Ball pointed out this requirement)
2387 1 : if (s% job% relax_initial_Z .and. .not. restart) then
2388 0 : write(*,1) 'relax_initial_Z', s% job% new_Z
2389 : call star_relax_Z(id, s% job% new_Z, s% relax_dlnZ, &
2390 0 : s% job% relax_Z_minq, s% job% relax_Z_maxq, ierr)
2391 0 : if (failed('star_relax_Z',ierr)) return
2392 0 : write(*, 1) 'new z', get_current_z(id, ierr)
2393 0 : if (failed('get_current_z',ierr)) return
2394 : end if
2395 :
2396 1 : if (s% job% relax_Z) then
2397 0 : write(*,1) 'relax_Z', s% job% new_Z
2398 : call star_relax_Z(id, s% job% new_Z, s% relax_dlnZ, &
2399 0 : s% job% relax_Z_minq, s% job% relax_Z_maxq, ierr)
2400 0 : if (failed('star_relax_Z',ierr)) return
2401 0 : write(*, 1) 'new z', get_current_z(id, ierr)
2402 0 : if (failed('get_current_z',ierr)) return
2403 : end if
2404 :
2405 1 : if (s% job% relax_initial_Y .and. .not. restart) then
2406 0 : write(*,1) 'relax_initial_Y', s% job% new_Y
2407 : call star_relax_Y(id, s% job% new_Y, s% relax_dY, &
2408 0 : s% job% relax_Y_minq, s% job% relax_Y_maxq, ierr)
2409 0 : if (failed('star_relax_Y',ierr)) return
2410 0 : write(*, 1) 'new y', get_current_y(id, ierr)
2411 0 : if (failed('get_current_y',ierr)) return
2412 : end if
2413 :
2414 1 : if (s% job% relax_Y) then
2415 0 : write(*,1) 'relax_Y', s% job% new_Y
2416 : call star_relax_Y(id, s% job% new_Y, s% relax_dY, &
2417 0 : s% job% relax_Y_minq, s% job% relax_Y_maxq, ierr)
2418 0 : if (failed('star_relax_Y',ierr)) return
2419 0 : write(*, 1) 'new y', get_current_y(id, ierr)
2420 0 : if (failed('get_current_y',ierr)) return
2421 : end if
2422 :
2423 1 : if (s% job% relax_mass) then
2424 0 : write(*, 1) 'relax_mass', s% job% new_mass
2425 0 : call star_relax_mass(id, s% job% new_mass, s% job% lg_max_abs_mdot, ierr)
2426 0 : if (failed('star_relax_mass',ierr)) return
2427 : end if
2428 :
2429 1 : if (s% job% relax_mass_to_remove_H_env) then
2430 0 : write(*, 1) 'relax_mass_to_remove_H_env_mass'
2431 : call star_relax_mass_to_remove_H_env( &
2432 0 : id, s% job% extra_mass_retained_by_remove_H_env, s% job% lg_max_abs_mdot, ierr)
2433 0 : if (failed('star_relax_mass_to_remove_H_env',ierr)) return
2434 : end if
2435 :
2436 1 : if (s% job% relax_dxdt_nuc_factor .or. &
2437 : (s% job% relax_initial_dxdt_nuc_factor .and. .not. restart)) then
2438 0 : write(*, 1) 'relax_dxdt_nuc_factor', s% job% new_dxdt_nuc_factor
2439 : call star_relax_dxdt_nuc_factor( &
2440 0 : id, s% job% new_dxdt_nuc_factor, s% job% dxdt_nuc_factor_multiplier, ierr)
2441 0 : if (failed('star_relax_dxdt_nuc_factor',ierr)) return
2442 : end if
2443 :
2444 1 : if (s% job% relax_eps_nuc_factor .or. &
2445 : (s% job% relax_initial_eps_nuc_factor .and. .not. restart)) then
2446 0 : write(*, 1) 'relax_eps_nuc_factor', s% job% new_eps_nuc_factor
2447 : call star_relax_eps_nuc_factor( &
2448 0 : id, s% job% new_eps_nuc_factor, s% job% eps_nuc_factor_multiplier, ierr)
2449 0 : if (failed('star_relax_eps_nuc_factor',ierr)) return
2450 : end if
2451 :
2452 1 : if (s% job% relax_opacity_max .or. &
2453 : (s% job% relax_initial_opacity_max .and. .not. restart)) then
2454 0 : write(*, 1) 'relax_opacity_max', s% job% new_opacity_max
2455 : call star_relax_opacity_max( &
2456 0 : id, s% job% new_opacity_max, s% job% opacity_max_multiplier, ierr)
2457 0 : if (failed('star_relax_opacity_max',ierr)) return
2458 : end if
2459 :
2460 1 : if (s% job% relax_max_surf_dq .or. &
2461 : (s% job% relax_initial_max_surf_dq .and. .not. restart)) then
2462 0 : write(*, 1) 'relax_max_surf_dq', s% job% new_max_surf_dq
2463 : call star_relax_max_surf_dq( &
2464 0 : id, s% job% new_max_surf_dq, s% job% max_surf_dq_multiplier, ierr)
2465 0 : if (failed('star_relax_max_surf_dq',ierr)) return
2466 : end if
2467 :
2468 1 : if (s% job% relax_initial_mass .and. .not. restart) then
2469 0 : write(*, 1) 'relax_initial_mass to new_mass', s% job% new_mass
2470 0 : call star_relax_mass(id, s% job% new_mass, s% job% lg_max_abs_mdot, ierr)
2471 0 : if (failed('relax_initial_mass',ierr)) return
2472 : end if
2473 :
2474 1 : if (s% job% relax_initial_mass_to_remove_H_env .and. .not. restart) then
2475 0 : write(*, 1) 'relax_initial_mass_to_remove_H_env'
2476 : call star_relax_mass_to_remove_H_env( &
2477 0 : id, s% job% extra_mass_retained_by_remove_H_env, s% job% lg_max_abs_mdot, ierr)
2478 0 : if (failed('relax_initial_mass_to_remove_H_env',ierr)) return
2479 : end if
2480 :
2481 1 : if (s% job% relax_mass_scale .or. &
2482 : (s% job% relax_initial_mass_scale .and. .not. restart)) then
2483 0 : write(*, 1) 'relax_mass_scale', s% job% new_mass
2484 : call star_relax_mass_scale( &
2485 : id, s% job% new_mass, s% job% dlgm_per_step, &
2486 0 : s% job% change_mass_years_for_dt, ierr)
2487 0 : if (failed('star_relax_mass_scale',ierr)) return
2488 : end if
2489 :
2490 1 : if (s% job% relax_core .or. &
2491 : (s% job% relax_initial_core .and. .not. restart)) then
2492 0 : write(*, 1) 'relax_core', s% job% new_core_mass
2493 : call star_relax_core( &
2494 : id, s% job% new_core_mass, s% job% dlg_core_mass_per_step, &
2495 : s% job% relax_core_years_for_dt, &
2496 0 : s% job% core_avg_rho, s% job% core_avg_eps, ierr)
2497 0 : if (failed('star_relax_core',ierr)) return
2498 : end if
2499 :
2500 1 : call do_remove_center(id, s, restart, ierr)
2501 1 : if (ierr /= 0) return
2502 :
2503 1 : if (s% job% relax_M_center .or. &
2504 : (s% job% relax_initial_M_center .and. .not. restart)) then
2505 0 : write(*, 1) 'relax_M_center', s% job% new_mass
2506 : call star_relax_M_center( &
2507 0 : id, s% job% new_mass, s% job% dlgm_per_step, s% job% relax_M_center_dt, ierr)
2508 0 : if (failed('star_relax_M_center',ierr)) return
2509 : end if
2510 :
2511 1 : if (s% job% relax_R_center .or. &
2512 : (s% job% relax_initial_R_center .and. .not. restart)) then
2513 0 : write(*, 1) 'relax_R_center', s% job% new_R_center
2514 : call star_relax_R_center( &
2515 0 : id, s% job% new_R_center, s% job% dlgR_per_step, s% job% relax_R_center_dt, ierr)
2516 0 : if (failed('star_relax_R_center',ierr)) return
2517 : end if
2518 :
2519 1 : if (s% job% relax_v_center .or. &
2520 : (s% job% relax_initial_v_center .and. .not. restart)) then
2521 0 : write(*, 1) 'relax_v_center', s% job% new_v_center
2522 : call star_relax_v_center( &
2523 0 : id, s% job% new_v_center, s% job% dv_per_step, s% job% relax_v_center_dt, ierr)
2524 0 : if (failed('star_relax_v_center',ierr)) return
2525 : end if
2526 :
2527 1 : if (s% job% relax_L_center .or. &
2528 : (s% job% relax_initial_L_center .and. .not. restart)) then
2529 0 : write(*, 1) 'relax_L_center', s% job% new_L_center
2530 : call star_relax_L_center( &
2531 0 : id, s% job% new_L_center, s% job% dlgL_per_step, s% job% relax_L_center_dt, ierr)
2532 0 : if (failed('star_relax_L_center',ierr)) return
2533 : end if
2534 :
2535 1 : if (s% job% relax_Tsurf_factor .or. &
2536 : (s% job% relax_initial_Tsurf_factor .and. .not. restart)) then
2537 0 : write(*,1) 'relax_Tsurf_factor', s% job% relax_to_this_Tsurf_factor
2538 : call star_relax_Tsurf_factor( &
2539 0 : id, s% job% relax_to_this_Tsurf_factor, s% job% dlogTsurf_factor, ierr)
2540 0 : if (failed('star_relax_Tsurf_factor',ierr)) return
2541 : end if
2542 :
2543 1 : if (s% job% relax_tau_factor .or. &
2544 : (s% job% relax_initial_tau_factor .and. .not. restart)) then
2545 0 : write(*,1) 'relax_tau_factor', s% job% relax_to_this_tau_factor
2546 : call star_relax_tau_factor( &
2547 0 : id, s% job% relax_to_this_tau_factor, s% job% dlogtau_factor, ierr)
2548 0 : if (failed('star_relax_tau_factor',ierr)) return
2549 : end if
2550 :
2551 1 : if (s% job% relax_opacity_factor .or. &
2552 : (s% job% relax_initial_opacity_factor .and. .not. restart)) then
2553 0 : write(*,1) 'relax_opacity_factor', s% job% relax_to_this_opacity_factor
2554 : call star_relax_opacity_factor( &
2555 0 : id, s% job% relax_to_this_opacity_factor, s% job% d_opacity_factor, ierr)
2556 0 : if (failed('star_relax_opacity_factor',ierr)) return
2557 : end if
2558 :
2559 1 : if (s% job% relax_irradiation .or. &
2560 : (s% job% relax_initial_irradiation .and. .not. restart)) then
2561 0 : write(*,2) 'relax_irradiation -- min steps', s% job% relax_irradiation_min_steps
2562 0 : write(*,1) 'relax_irradiation -- max yrs dt', s% job% relax_irradiation_max_yrs_dt
2563 : call star_relax_irradiation(id, &
2564 : s% job% relax_irradiation_min_steps, &
2565 : s% job% relax_to_this_irrad_flux, s% job% irrad_col_depth, &
2566 0 : s% job% relax_irradiation_max_yrs_dt, ierr)
2567 0 : if (failed('star_relax_irradiation',ierr)) return
2568 : end if
2569 :
2570 1 : if (s% job% relax_mass_change .or. &
2571 : (s% job% relax_initial_mass_change .and. .not. restart)) then
2572 0 : write(*,2) 'relax_mass_change -- min steps', &
2573 0 : s% job% relax_mass_change_min_steps
2574 0 : write(*,1) 'relax_mass_change -- max yrs dt', &
2575 0 : s% job% relax_mass_change_max_yrs_dt
2576 0 : write(*,1) 'relax_mass_change -- initial_mass_change', &
2577 0 : s% job% relax_mass_change_init_mdot
2578 0 : write(*,1) 'relax_mass_change -- final_mass_change', &
2579 0 : s% job% relax_mass_change_final_mdot
2580 : call star_relax_mass_change(id, &
2581 : s% job% relax_mass_change_min_steps, &
2582 : s% job% relax_mass_change_init_mdot, &
2583 : s% job% relax_mass_change_final_mdot, &
2584 0 : s% job% relax_mass_change_max_yrs_dt, ierr)
2585 0 : if (failed('star_relax_mass_change',ierr)) return
2586 : end if
2587 :
2588 1 : call do_remove_initial_surface(id, s, restart, ierr)
2589 1 : if (ierr /= 0) return
2590 :
2591 1 : call do_remove_surface(id, s, ierr)
2592 1 : if (ierr /= 0) return
2593 :
2594 1 : if (s% rotation_flag .and. s% job% relax_omega) then
2595 0 : write(*,1) 'new_omega', s% job% new_omega
2596 : call star_relax_uniform_omega( &
2597 : id, relax_to_new_omega, &
2598 : s% job% new_omega, 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 : end if
2602 :
2603 1 : if (s% rotation_flag .and. s% job% relax_initial_omega .and. .not. restart) then
2604 : call star_relax_uniform_omega( &
2605 : id, relax_to_new_omega, &
2606 : s% job% new_omega, s% job% num_steps_to_relax_rotation,&
2607 0 : s% job% relax_omega_max_yrs_dt, ierr)
2608 0 : if (failed('star_relax_uniform_omega',ierr)) return
2609 0 : write(*,1) 'new_omega', s% job% new_omega
2610 : end if
2611 :
2612 1 : if (s% rotation_flag .and. s% job% relax_omega_div_omega_crit) then
2613 0 : if (failed('star_surface_omega_crit',ierr)) return
2614 : call star_relax_uniform_omega( &
2615 : id, relax_to_new_omega_div_omega_crit, &
2616 : s% job% new_omega_div_omega_crit, &
2617 : s% job% num_steps_to_relax_rotation,&
2618 0 : s% job% relax_omega_max_yrs_dt, ierr)
2619 0 : if (failed('star_relax_uniform_omega',ierr)) return
2620 0 : write(*,2) 'new_omega_div_omega_crit', &
2621 0 : s% model_number, s% job% new_omega_div_omega_crit
2622 : end if
2623 :
2624 : if (s% rotation_flag .and. &
2625 1 : s% job% relax_initial_omega_div_omega_crit .and. .not. restart) then
2626 0 : if (failed('star_surface_omega_crit',ierr)) return
2627 : call star_relax_uniform_omega( &
2628 : id, relax_to_new_omega_div_omega_crit, &
2629 : s% job% new_omega_div_omega_crit, &
2630 : s% job% num_steps_to_relax_rotation,&
2631 0 : s% job% relax_omega_max_yrs_dt, ierr)
2632 0 : if (failed('star_relax_uniform_omega',ierr)) return
2633 0 : write(*,2) 'new_omega_div_omega_crit', &
2634 0 : s% model_number, s% job% new_omega_div_omega_crit
2635 : end if
2636 :
2637 1 : if (s% rotation_flag .and. s% job% relax_surface_rotation_v) then
2638 : call star_relax_uniform_omega( &
2639 : id, relax_to_new_surface_rotation_v, &
2640 : s% job% new_surface_rotation_v, s% job% num_steps_to_relax_rotation,&
2641 0 : s% job% relax_omega_max_yrs_dt, ierr)
2642 0 : if (failed('star_relax_uniform_omega',ierr)) return
2643 0 : s% job% new_omega = s% job% new_surface_rotation_v*1d5/s% r(1)
2644 0 : write(*,1) 'new_surface_rotation_v', &
2645 0 : s% job% new_surface_rotation_v, s% job% new_omega
2646 : end if
2647 :
2648 : if (s% rotation_flag .and. &
2649 1 : s% job% relax_initial_surface_rotation_v .and. .not. restart) then
2650 0 : write(*,1) 'new_omega', s% job% new_omega
2651 0 : write(*,*) 'call star_relax_uniform_omega'
2652 : call star_relax_uniform_omega( &
2653 : id, relax_to_new_surface_rotation_v, &
2654 : s% job% new_surface_rotation_v, s% job% num_steps_to_relax_rotation,&
2655 0 : s% job% relax_omega_max_yrs_dt, ierr)
2656 0 : if (failed('star_relax_uniform_omega',ierr)) return
2657 0 : write(*,2) 'new_surface_rotation_v', &
2658 0 : s% model_number, s% job% new_surface_rotation_v
2659 : end if
2660 :
2661 1 : if (s% job% set_max_dt_to_frac_lifetime) then
2662 0 : log_m = log10(s% star_mass) ! in Msun units
2663 0 : log_lifetime = 9.921d0 - (3.6648d0 + (1.9697d0 - 0.9369d0*log_m)*log_m)*log_m
2664 : ! Iben & Laughlin (1989) as quoted in H&K (eqn 2.3)
2665 0 : max_dt = s% job% max_frac_of_lifetime_per_step*secyer*exp10(log_lifetime)
2666 0 : if (max_dt < s% max_timestep) then
2667 0 : s% max_timestep = max_dt
2668 0 : write(*, *) 'set_max_dt_to_frac_lifetime: lg(maxdt/secyer)', &
2669 0 : log10(s% max_timestep/secyer)
2670 : end if
2671 : end if
2672 :
2673 : ! print out info about selected non-standard parameter settings
2674 :
2675 1 : write(*,*) 'net name ' // trim(s% net_name)
2676 :
2677 1 : if (s% do_element_diffusion) &
2678 0 : write(*,*) 'do_element_diffusion', s% do_element_diffusion
2679 :
2680 1 : if (s% RSP_flag) &
2681 0 : write(*,*) 'RSP_flag', s% RSP_flag
2682 :
2683 1 : if (s% v_flag) &
2684 0 : write(*,*) 'v_flag', s% v_flag
2685 :
2686 1 : if (s% u_flag) &
2687 0 : write(*,*) 'u_flag', s% u_flag
2688 :
2689 1 : if (s% rotation_flag) &
2690 0 : write(*,*) 'rotation_flag', s% rotation_flag
2691 :
2692 1 : if (s% w_div_wc_flag) &
2693 0 : write(*,*) 'w_div_wc_flag', s% w_div_wc_flag
2694 :
2695 1 : if (s% j_rot_flag) &
2696 0 : write(*,*) 'j_rot_flag', s% j_rot_flag
2697 :
2698 1 : if (s% mix_factor /= 1d0) &
2699 0 : write(*,1) 'mix_factor', s% mix_factor
2700 :
2701 1 : if (abs(s% tau_base - 2d0/3d0) > 1d-4) &
2702 0 : write(*,1) 'tau_base', s% tau_base
2703 :
2704 1 : if (abs(s% tau_factor - 1) > 1d-4) &
2705 0 : write(*,1) 'tau_factor', s% tau_factor
2706 :
2707 1 : if (s% eps_grav_factor /= 1) &
2708 0 : write(*,1) 'eps_grav_factor', s% eps_grav_factor
2709 :
2710 1 : if (s% eps_mdot_factor /= 1) &
2711 0 : write(*,1) 'eps_mdot_factor', s% eps_mdot_factor
2712 :
2713 1 : if (s% dxdt_nuc_factor /= 1) &
2714 0 : write(*,1) 'dxdt_nuc_factor', s% dxdt_nuc_factor
2715 :
2716 1 : if (.NOT. ( &
2717 : s% atm_option == 'T_tau' .AND. &
2718 : s% atm_T_tau_relation == 'Eddington' .AND. &
2719 : s% atm_T_tau_opacity == 'fixed')) &
2720 0 : write(*,1) 'atm_option: ' // trim(s% atm_option)
2721 :
2722 1 : if (s% M_center /= 0) then
2723 0 : write(*,1) 'xmstar/mstar', s% xmstar/s% mstar
2724 0 : write(*,1) 'xmstar (g)', s% xmstar
2725 0 : write(*,1) 'M_center (g)', s% M_center
2726 0 : write(*,1) 'xmstar/Msun', s% xmstar/Msun
2727 0 : write(*,1) 'M_center/Msun', s% M_center/Msun
2728 : end if
2729 :
2730 1 : if (s% v_flag .or. s% u_flag) then
2731 0 : if (s% v_center /= 0) &
2732 0 : write(*,1) 'v_center (cm/s)', s% v_center
2733 : end if
2734 :
2735 1 : if (s% R_center /= 0) then
2736 0 : write(*,1) 'R_center (cm)', s% R_center
2737 0 : write(*,1) 'R_center/Rsun', s% R_center/Rsun
2738 0 : write(*,1) 'core density', &
2739 0 : s% M_center/(4*pi/3*s% R_center*s% R_center*s% R_center)
2740 : end if
2741 :
2742 1 : if (s% L_center /= 0) &
2743 0 : write(*,1) 'L_center/Lsun', s% L_center/Lsun
2744 :
2745 1 : if (s% opacity_max > 0) &
2746 0 : write(*,1) 'opacity_max', s% opacity_max
2747 :
2748 1 : if (s% opacity_min > 0) &
2749 0 : write(*,1) 'opacity_min', s% opacity_min
2750 :
2751 1 : if (s% job% show_net_reactions_info) then
2752 0 : write(*,'(a)') ' net reactions '
2753 0 : call show_net_reactions_and_info(s% net_handle, 6, ierr)
2754 0 : if (failed('show_net_reactions_and_info',ierr)) return
2755 : end if
2756 :
2757 1 : if (s% job% list_net_reactions) then
2758 0 : write(*,'(a)') ' net reactions '
2759 0 : call show_net_reactions(s% net_handle, 6, ierr)
2760 0 : if (failed('show_net_reactions',ierr)) return
2761 : end if
2762 :
2763 : if (s% job% set_cumulative_energy_error .or. &
2764 1 : (s% job% set_initial_cumulative_energy_error .and. .not. restart) .or. &
2765 : (s% model_number == s% job% set_cumulative_energy_error_at_step)) then
2766 0 : write(*,1) 'set_cumulative_energy_error', s% job% new_cumulative_energy_error
2767 0 : s% cumulative_energy_error = s% job% new_cumulative_energy_error
2768 : end if
2769 :
2770 1 : if (s% job% show_net_species_info) then
2771 0 : write(*,'(a)') ' species'
2772 0 : do j=1,s% species
2773 0 : write(*,'(i6,3x,a)') j, chem_isos% name(s% chem_id(j))
2774 : end do
2775 0 : write(*,'(A)')
2776 : end if
2777 :
2778 1 : if (s% job% show_eqns_and_vars_names) then
2779 0 : do i=1,s% nvar_total
2780 0 : write(*,*) i, s% nameofvar(i), s% nameofequ(i)
2781 : end do
2782 0 : write(*,'(A)')
2783 : end if
2784 :
2785 1 : write(*,*) 'kap_option ' // trim(kap_option_str(s% kap_rq% kap_option))
2786 1 : write(*,*) 'kap_CO_option ' // trim(kap_CO_option_str(s% kap_rq% kap_CO_option))
2787 1 : write(*,*) 'kap_lowT_option ' // trim(kap_lowT_option_str(s% kap_rq% kap_lowT_option))
2788 1 : write(*,2) 'OMP_NUM_THREADS', utils_omp_get_max_threads()
2789 :
2790 1 : call check_if_want_to_stop_warnings(s)
2791 :
2792 : contains
2793 :
2794 0 : subroutine do_relax_initial_composition(ierr)
2795 1 : use utils_lib
2796 : integer, intent(out) :: ierr
2797 0 : real(dp), pointer :: xq(:), xa(:,:)
2798 : integer :: num_pts, num_species, i, iounit
2799 : include 'formats'
2800 :
2801 0 : write(*,'(A)')
2802 0 : write(*,1) 'relax_initial_composition'
2803 :
2804 : open(newunit=iounit, file=trim(s% job% relax_composition_filename), &
2805 0 : status='old', action='read', iostat=ierr)
2806 0 : if (ierr /= 0) then
2807 0 : write(*,*) 'open failed', ierr, iounit
2808 0 : write(*, '(a)') 'failed to open ' // trim(s% job% relax_composition_filename)
2809 0 : return
2810 : end if
2811 0 : read(iounit, *, iostat=ierr) num_pts, num_species
2812 0 : if (ierr /= 0) then
2813 0 : close(iounit)
2814 : write(*, '(a)') 'failed while trying to read 1st line of ' // &
2815 0 : trim(s% job% relax_composition_filename)
2816 0 : return
2817 : end if
2818 0 : if(num_species /= s% species) then
2819 0 : write(*,*) 'Error in ',trim(s% job% relax_composition_filename)
2820 0 : write(*,'(a,I4,a)') 'got ',num_species,' species'
2821 0 : write(*,'(a,I4,a)') 'expected ', s% species,' species'
2822 0 : write(*,'(A)')
2823 0 : ierr=-1
2824 0 : return
2825 : end if
2826 0 : allocate(xq(num_pts), xa(num_species,num_pts))
2827 0 : do i = 1, num_pts
2828 0 : read(iounit,*,iostat=ierr) xq(i), xa(1:num_species,i)
2829 0 : if (ierr /= 0) then
2830 0 : close(iounit)
2831 : write(*, '(a)') &
2832 0 : 'failed while trying to read ' // trim(s% job% relax_composition_filename)
2833 0 : write(*,*) 'line', i+1
2834 0 : write(*,*) 'perhaps wrong info in 1st line?'
2835 0 : write(*,*) '1st line must have num_pts and num_species in that order'
2836 0 : deallocate(xq,xa)
2837 0 : return
2838 : end if
2839 : end do
2840 0 : close(iounit)
2841 :
2842 : call star_relax_composition( &
2843 0 : id, s% job% num_steps_to_relax_composition, num_pts, num_species, xa, xq, ierr)
2844 0 : deallocate(xq,xa)
2845 :
2846 0 : end subroutine do_relax_initial_composition
2847 :
2848 0 : subroutine do_relax_initial_angular_momentum(ierr)
2849 : use utils_lib
2850 : integer, intent(out) :: ierr
2851 0 : real(dp), pointer :: xq(:), angular_momentum(:)
2852 : integer :: num_pts, i, iounit
2853 : include 'formats'
2854 :
2855 0 : write(*,'(A)')
2856 0 : write(*,1) 'relax_initial_angular_momentum'
2857 :
2858 : open(newunit=iounit, file=trim(s% job% relax_angular_momentum_filename), &
2859 0 : status='old', action='read', iostat=ierr)
2860 0 : if (ierr /= 0) then
2861 0 : write(*,*) 'open failed', ierr, iounit
2862 0 : write(*, '(a)') 'failed to open "' // trim(s% job% relax_angular_momentum_filename)//'"'
2863 0 : return
2864 : end if
2865 0 : read(iounit, *, iostat=ierr) num_pts
2866 0 : if (ierr /= 0) then
2867 0 : close(iounit)
2868 : write(*, '(a)') 'failed while trying to read 1st line of ' // &
2869 0 : trim(s% job% relax_angular_momentum_filename)
2870 0 : return
2871 : end if
2872 0 : allocate(xq(num_pts), angular_momentum(num_pts))
2873 0 : do i = 1, num_pts
2874 0 : read(iounit,*,iostat=ierr) xq(i), angular_momentum(i)
2875 0 : if (ierr /= 0) then
2876 0 : close(iounit)
2877 : write(*, '(a)') &
2878 0 : 'failed while trying to read ' // trim(s% job% relax_angular_momentum_filename)
2879 0 : write(*,*) 'line', i+1
2880 0 : write(*,*) 'perhaps wrong info in 1st line?'
2881 0 : write(*,*) '1st line must have num_pts'
2882 0 : deallocate(xq,angular_momentum)
2883 0 : return
2884 : end if
2885 : end do
2886 0 : close(iounit)
2887 : call star_relax_angular_momentum(id, s% job% max_steps_to_relax_angular_momentum, &
2888 0 : num_pts, angular_momentum, xq, ierr)
2889 0 : deallocate(xq,angular_momentum)
2890 0 : end subroutine do_relax_initial_angular_momentum
2891 :
2892 0 : subroutine do_relax_initial_entropy(ierr)
2893 : use utils_lib
2894 : use eos_def
2895 : integer, intent(out) :: ierr
2896 : ! arrays into which data from the input file is read.
2897 : ! in case any of the eos* options is used, input from the
2898 : ! file is read into var1 and var2, and the chosen eos function
2899 : ! is used to extract the entropy from that pair.
2900 0 : real(dp), pointer :: xq(:), entropy(:)
2901 : real(dp) :: var1, var2
2902 : integer :: num_pts, i, k, iounit
2903 : ! these are needed to call eosPT_get
2904 : real(dp) :: Rho, log10Rho, dlnRho_dlnPgas_const_T, dlnRho_dlnT_const_Pgas
2905 : real(dp) :: log10T
2906 : ! these are needed to call eosDT_get_T
2907 : real(dp) :: T_guess_gas, T_guess_rad, logT_guess
2908 : integer :: eos_calls
2909 : ! these are used for all eos calls
2910 : real(dp), dimension(num_eos_basic_results) :: res, d_dlnd, d_dlnT
2911 0 : real(dp), dimension(num_eos_d_dxa_results, s% species) :: d_dxa
2912 : real(dp), parameter :: logT_tol = 1d-8, logE_tol = 1d-8
2913 : integer, parameter :: MAX_ITERS = 20
2914 : include 'formats'
2915 :
2916 0 : write(*,'(A)')
2917 0 : write(*,1) 'relax_initial_entropy'
2918 :
2919 : open(newunit=iounit, file=trim(s% job% relax_entropy_filename), &
2920 0 : status='old', action='read', iostat=ierr)
2921 0 : if (ierr /= 0) then
2922 0 : write(*,*) 'open failed', ierr, iounit
2923 0 : write(*, '(a)') 'failed to open "' // trim(s% job% relax_entropy_filename)//'"'
2924 0 : return
2925 : end if
2926 0 : read(iounit, *, iostat=ierr) num_pts
2927 0 : if (ierr /= 0) then
2928 0 : close(iounit)
2929 : write(*, '(a)') 'failed while trying to read 1st line of ' // &
2930 0 : trim(s% job% relax_entropy_filename)
2931 0 : return
2932 : end if
2933 0 : if (.not. (s% job% get_entropy_for_relax_from_eos == '' .or. &
2934 : s% job% get_entropy_for_relax_from_eos == 'eosDT' .or. &
2935 : s% job% get_entropy_for_relax_from_eos == 'eosPT' .or. &
2936 : s% job% get_entropy_for_relax_from_eos == 'eosDE')) then
2937 0 : ierr = 1
2938 0 : write(*,*) 'invalid value for get_entropy_for_relax_from_eos =', &
2939 0 : s% job% get_entropy_for_relax_from_eos
2940 : end if
2941 0 : allocate(xq(num_pts), entropy(num_pts))
2942 0 : do i = 1, num_pts
2943 0 : if (s% job% get_entropy_for_relax_from_eos == '') then
2944 0 : read(iounit,*,iostat=ierr) xq(i), entropy(i)
2945 : else
2946 0 : read(iounit,*,iostat=ierr) xq(i), var1, var2
2947 : ! get nearest value matching xq for the composition TODO: interpolate
2948 0 : do k=1, s% nz-1
2949 0 : if(1-s% q(k) <= xq(i) .and. 1-s% q(k+1) >= xq(i)) then
2950 : exit
2951 : end if
2952 : end do
2953 : ! get entropy
2954 0 : if (s% job% get_entropy_for_relax_from_eos == 'eosDT') then
2955 : call eosDT_get( &
2956 : s% eos_handle, &
2957 : s% species, s% chem_id, s% net_iso, s% xa(:,k), &
2958 : var1, log10(var1), var2, log10(var2), &
2959 0 : res, d_dlnd, d_dlnT, d_dxa, ierr)
2960 0 : if (ierr /= 0) then
2961 0 : write(*,*) "failed in eosDT_get"
2962 0 : return
2963 : end if
2964 0 : entropy(i) = exp(res(i_lnS))
2965 0 : else if (s% job% get_entropy_for_relax_from_eos == 'eosPT') then
2966 : call eosPT_get( &
2967 : s% eos_handle, &
2968 : s% species, s% chem_id, s% net_iso, s% xa(:,k), &
2969 : var1, log10(var1), var2, log10(var2), &
2970 : Rho, log10Rho, dlnRho_dlnPgas_const_T, dlnRho_dlnT_const_Pgas, &
2971 0 : res, d_dlnd, d_dlnT, d_dxa, ierr)
2972 0 : if (ierr /= 0) then
2973 0 : write(*,*) "failed in eosPT_get"
2974 0 : return
2975 : end if
2976 0 : entropy(i) = exp(res(i_lnS))
2977 : else
2978 0 : T_guess_gas = 2*var2*s% abar(k)*mp/(3*kerg*(1+s% zbar(k))) ! ideal gas (var2=energy)
2979 0 : T_guess_rad = pow(var2/crad,0.25d0)
2980 0 : logT_guess = log10(min(T_guess_gas,T_guess_rad))
2981 : call eosDT_get_T( &
2982 : s% eos_handle, &
2983 : s% species, s% chem_id, s% net_iso, s% xa(:,k), &
2984 : log10(var1), i_lnE, log10(var2)*ln10, &
2985 : logT_tol, logE_tol*ln10, MAX_ITERS, logT_guess, &
2986 : arg_not_provided, arg_not_provided, arg_not_provided, arg_not_provided, &
2987 : log10T, res, d_dlnd, d_dlnT, d_dxa, &
2988 0 : eos_calls, ierr)
2989 0 : if (ierr /= 0) then
2990 0 : write(*,*) "failed in eosDT_get_T (as eosDE)"
2991 0 : return
2992 : end if
2993 0 : entropy(i) = exp(res(i_lnS))
2994 : end if
2995 : end if
2996 0 : if (ierr /= 0) then
2997 0 : close(iounit)
2998 : write(*, '(a)') &
2999 0 : 'failed while trying to read ' // trim(s% job% relax_entropy_filename)
3000 0 : write(*,*) 'line', i+1
3001 0 : write(*,*) 'perhaps wrong info in 1st line?'
3002 0 : write(*,*) '1st line must have num_pts'
3003 0 : deallocate(xq,entropy)
3004 0 : return
3005 : end if
3006 : end do
3007 0 : close(iounit)
3008 0 : call star_relax_entropy(id, s% job% max_steps_to_relax_entropy, num_pts, entropy, xq, ierr)
3009 0 : deallocate(xq,entropy)
3010 0 : end subroutine do_relax_initial_entropy
3011 :
3012 : end subroutine do_star_job_controls_after
3013 :
3014 :
3015 1 : subroutine do_remove_center(id, s, restart, ierr)
3016 : integer, intent(in) :: id
3017 : type (star_info), pointer :: s
3018 : logical, intent(in) :: restart
3019 : integer, intent(out) :: ierr
3020 : include 'formats'
3021 :
3022 1 : if (s% job% remove_center_by_temperature > 0) then
3023 0 : write(*, 1) 'remove_center_by_temperature', s% job% remove_center_by_temperature
3024 : call star_remove_center_by_temperature( &
3025 0 : id, s% job% remove_center_by_temperature, ierr)
3026 0 : if (failed('star_remove_center_by_temperature',ierr)) return
3027 : end if
3028 :
3029 1 : if (s% job% remove_initial_center_by_temperature > 0 .and. .not. restart) then
3030 0 : write(*, 1) 'remove_initial_center_by_temperature', &
3031 0 : s% job% remove_initial_center_by_temperature
3032 : call star_remove_center_by_temperature( &
3033 0 : id, s% job% remove_initial_center_by_temperature, ierr)
3034 0 : if (failed('star_remove_center_by_temperature',ierr)) return
3035 : end if
3036 :
3037 1 : if (s% job% remove_center_by_radius_cm > s% R_center .and. &
3038 : s% job% remove_center_by_radius_cm < s% r(1)) then
3039 0 : write(*, 1) 'remove_center_by_radius_cm', &
3040 0 : s% job% remove_center_by_radius_cm
3041 : call star_remove_center_by_radius_cm( &
3042 0 : id, s% job% remove_center_by_radius_cm, ierr)
3043 0 : if (failed('star_remove_center_by_radius_cm',ierr)) return
3044 : end if
3045 :
3046 : if (s% job% remove_initial_center_by_radius_cm > s% R_center .and. &
3047 1 : s% job% remove_initial_center_by_radius_cm < s% r(1) .and. .not. restart) then
3048 0 : write(*, 1) 'remove_initial_center_by_radius_cm', &
3049 0 : s% job% remove_initial_center_by_radius_cm
3050 : call star_remove_center_by_radius_cm( &
3051 0 : id, s% job% remove_initial_center_by_radius_cm, ierr)
3052 0 : if (failed('star_remove_center_by_radius_cm',ierr)) return
3053 : end if
3054 :
3055 : if (s% job% remove_initial_center_by_he4 > 0d0 .and. &
3056 : s% job% remove_initial_center_by_he4 < 1d0 &
3057 1 : .and. .not. restart) then
3058 0 : write(*, 1) 'remove_initial_center_by_he4', &
3059 0 : s% job% remove_initial_center_by_he4
3060 : call star_remove_center_by_he4( &
3061 0 : id, s% job% remove_initial_center_by_he4, ierr)
3062 0 : if (failed('star_remove_initial_center_by_he4',ierr)) return
3063 : end if
3064 :
3065 1 : if (s% job% remove_center_by_he4 > 0d0 .and. &
3066 : s% job% remove_center_by_he4 < 1d0) then
3067 0 : write(*, 1) 'remove_center_by_he4', &
3068 0 : s% job% remove_center_by_he4
3069 : call star_remove_center_by_he4( &
3070 0 : id, s% job% remove_center_by_he4, ierr)
3071 0 : if (failed('star_remove_center_by_he4',ierr)) return
3072 : end if
3073 :
3074 : if (s% job% remove_initial_center_by_c12_o16 > 0d0 .and. &
3075 : s% job% remove_initial_center_by_c12_o16 < 1d0 &
3076 1 : .and. .not. restart) then
3077 0 : write(*, 1) 'remove_initial_center_by_c12_o16', &
3078 0 : s% job% remove_initial_center_by_c12_o16
3079 : call star_remove_center_by_c12_o16( &
3080 0 : id, s% job% remove_initial_center_by_c12_o16, ierr)
3081 0 : if (failed('star_remove_initial_center_by_c12_o16',ierr)) return
3082 : end if
3083 :
3084 1 : if (s% job% remove_center_by_c12_o16 > 0d0 .and. &
3085 : s% job% remove_center_by_c12_o16 < 1d0) then
3086 0 : write(*, 1) 'remove_center_by_c12_o16', &
3087 0 : s% job% remove_center_by_c12_o16
3088 : call star_remove_center_by_c12_o16( &
3089 0 : id, s% job% remove_center_by_c12_o16, ierr)
3090 0 : if (failed('star_remove_center_by_c12_o16',ierr)) return
3091 : end if
3092 :
3093 : if (s% job% remove_initial_center_by_si28 > 0d0 .and. &
3094 : s% job% remove_initial_center_by_si28 < 1d0 &
3095 1 : .and. .not. restart) then
3096 0 : write(*, 1) 'remove_initial_center_by_si28', &
3097 0 : s% job% remove_initial_center_by_si28
3098 : call star_remove_center_by_si28( &
3099 0 : id, s% job% remove_initial_center_by_si28, ierr)
3100 0 : if (failed('star_remove_initial_center_by_si28',ierr)) return
3101 : end if
3102 :
3103 1 : if (s% job% remove_center_by_si28 > 0d0 .and. &
3104 : s% job% remove_center_by_si28 < 1d0) then
3105 0 : write(*, 1) 'remove_center_by_si28', &
3106 0 : s% job% remove_center_by_si28
3107 : call star_remove_center_by_si28( &
3108 0 : id, s% job% remove_center_by_si28, ierr)
3109 0 : if (failed('star_remove_center_by_si28',ierr)) return
3110 : end if
3111 :
3112 : if (s% job% remove_initial_center_to_reduce_co56_ni56 > 0d0 &
3113 1 : .and. .not. restart) then
3114 0 : write(*, 1) 'remove_initial_center_to_reduce_co56_ni56', &
3115 0 : s% job% remove_initial_center_to_reduce_co56_ni56
3116 : call star_remove_center_to_reduce_co56_ni56( &
3117 0 : id, s% job% remove_initial_center_to_reduce_co56_ni56, ierr)
3118 0 : if (failed('star_remove_initial_center_to_reduce_co56_ni56',ierr)) return
3119 : end if
3120 :
3121 1 : if (s% job% remove_center_to_reduce_co56_ni56 > 0d0) then
3122 0 : write(*, 1) 'remove_center_to_reduce_co56_ni56', &
3123 0 : s% job% remove_center_to_reduce_co56_ni56
3124 : call star_remove_center_to_reduce_co56_ni56( &
3125 0 : id, s% job% remove_center_to_reduce_co56_ni56, ierr)
3126 0 : if (failed('star_remove_center_to_reduce_co56_ni56',ierr)) return
3127 : end if
3128 :
3129 : if (s% job% remove_initial_center_by_ye > 0d0 .and. &
3130 : s% job% remove_initial_center_by_ye < 1d0 &
3131 1 : .and. .not. restart) then
3132 0 : write(*, 1) 'remove_initial_center_by_ye', &
3133 0 : s% job% remove_initial_center_by_ye
3134 : call star_remove_center_by_ye( &
3135 0 : id, s% job% remove_initial_center_by_ye, ierr)
3136 0 : if (failed('star_remove_initial_center_by_ye',ierr)) return
3137 : end if
3138 :
3139 1 : if (s% job% remove_center_by_ye > 0d0 .and. &
3140 : s% job% remove_center_by_ye < 1d0) then
3141 0 : write(*, 1) 'remove_center_by_ye', &
3142 0 : s% job% remove_center_by_ye
3143 : call star_remove_center_by_ye( &
3144 0 : id, s% job% remove_center_by_ye, ierr)
3145 0 : if (failed('star_remove_center_by_ye',ierr)) return
3146 : end if
3147 :
3148 : if (s% job% remove_initial_center_by_entropy > 0d0 &
3149 1 : .and. .not. restart) then
3150 0 : write(*, 1) 'remove_initial_center_by_entropy', &
3151 0 : s% job% remove_initial_center_by_entropy
3152 : call star_remove_center_by_entropy( &
3153 0 : id, s% job% remove_initial_center_by_entropy, ierr)
3154 0 : if (failed('star_remove_initial_center_by_entropy',ierr)) return
3155 : end if
3156 :
3157 1 : if (s% job% remove_center_by_entropy > 0d0) then
3158 0 : write(*, 1) 'remove_center_by_entropy', &
3159 0 : s% job% remove_center_by_entropy
3160 : call star_remove_center_by_entropy( &
3161 0 : id, s% job% remove_center_by_entropy, ierr)
3162 0 : if (failed('star_remove_center_by_entropy',ierr)) return
3163 : end if
3164 :
3165 : if (s% job% remove_initial_center_by_infall_kms /= 0d0 &
3166 1 : .and. .not. restart) then
3167 0 : write(*, 1) 'remove_initial_center_by_infall_kms', &
3168 0 : s% job% remove_initial_center_by_infall_kms
3169 : call star_remove_center_by_infall_kms( &
3170 0 : id, s% job% remove_initial_center_by_infall_kms, ierr)
3171 0 : if (failed('star_remove_initial_center_by_infall_kms',ierr)) return
3172 : end if
3173 :
3174 1 : if (s% job% remove_center_by_infall_kms /= 0d0) then
3175 0 : write(*, 1) 'remove_center_by_infall_kms', &
3176 0 : s% job% remove_center_by_infall_kms
3177 : call star_remove_center_by_infall_kms( &
3178 0 : id, s% job% remove_center_by_infall_kms, ierr)
3179 0 : if (failed('star_remove_center_by_infall_kms',ierr)) return
3180 : end if
3181 :
3182 : if (s% job% remove_initial_center_at_inner_max_abs_v &
3183 1 : .and. .not. restart) then
3184 0 : write(*, 1) 'remove_initial_center_at_inner_max_abs_v'
3185 0 : call star_remove_center_at_inner_max_abs_v(id, ierr)
3186 0 : if (failed('remove_center_at_inner_max_abs_v',ierr)) return
3187 : end if
3188 :
3189 1 : if (s% job% remove_center_at_inner_max_abs_v) then
3190 0 : write(*, 1) 'remove_initial_center_at_inner_max_abs_v'
3191 0 : call star_remove_center_at_inner_max_abs_v(id, ierr)
3192 0 : if (failed('remove_center_at_inner_max_abs_v',ierr)) return
3193 : end if
3194 :
3195 1 : if (s% job% remove_initial_fe_core .and. .not. restart) then
3196 0 : write(*, 1) 'remove_initial_fe_core'
3197 0 : call star_remove_fe_core(id, ierr)
3198 0 : if (failed('remove_fe_core',ierr)) return
3199 : end if
3200 :
3201 1 : if (s% job% remove_fe_core) then
3202 0 : write(*, 1) 'remove_initial_fe_core'
3203 0 : call star_remove_fe_core(id, ierr)
3204 0 : if (failed('remove_fe_core',ierr)) return
3205 : end if
3206 :
3207 : if (s% job% remove_initial_center_by_mass_fraction_q > 0d0 .and. &
3208 : s% job% remove_initial_center_by_mass_fraction_q < 1d0 &
3209 1 : .and. .not. restart) then
3210 0 : write(*, 1) 'remove_initial_center_by_mass_fraction_q', &
3211 0 : s% job% remove_initial_center_by_mass_fraction_q
3212 : call star_remove_center_by_mass_fraction_q( &
3213 0 : id, s% job% remove_initial_center_by_mass_fraction_q, ierr)
3214 0 : if (failed('star_remove_initial_center_by_mass_fraction_q',ierr)) return
3215 : end if
3216 :
3217 1 : if (s% job% remove_center_by_mass_fraction_q > 0d0 .and. &
3218 : s% job% remove_center_by_mass_fraction_q < 1d0) then
3219 0 : write(*, 1) 'remove_center_by_mass_fraction_q', &
3220 0 : s% job% remove_center_by_mass_fraction_q
3221 : call star_remove_center_by_mass_fraction_q( &
3222 0 : id, s% job% remove_center_by_mass_fraction_q, ierr)
3223 0 : if (failed('star_remove_center_by_mass_fraction_q',ierr)) return
3224 : end if
3225 :
3226 1 : if (s% job% remove_center_by_delta_mass_gm > 0) then
3227 0 : write(*, 1) 'remove_center_by_delta_mass_gm', &
3228 0 : s% job% remove_center_by_delta_mass_gm
3229 : call star_remove_center_by_mass_gm(id, &
3230 0 : s% M_center + s% job% remove_center_by_delta_mass_gm, ierr)
3231 0 : if (failed('star_remove_center_by_mass',ierr)) return
3232 : end if
3233 :
3234 1 : if (s% job% remove_initial_center_by_delta_mass_gm > 0 .and. &
3235 : .not. restart) then
3236 0 : write(*, 1) 'remove_initial_center_by_delta_mass_gm', &
3237 0 : s% job% remove_initial_center_by_delta_mass_gm
3238 : call star_remove_center_by_mass_gm(id, &
3239 0 : s% M_center + s% job% remove_initial_center_by_delta_mass_gm, ierr)
3240 0 : if (failed('star_remove_center_by_mass',ierr)) return
3241 : end if
3242 :
3243 1 : if (s% job% remove_center_by_delta_mass_Msun > 0) then
3244 0 : write(*, 1) 'remove_center_by_delta_mass_Msun', &
3245 0 : s% job% remove_center_by_delta_mass_Msun
3246 : call star_remove_center_by_mass_gm(id, &
3247 0 : s% M_center + s% job% remove_center_by_delta_mass_Msun*Msun, ierr)
3248 0 : if (failed('star_remove_center_by_mass',ierr)) return
3249 : end if
3250 :
3251 1 : if (s% job% remove_initial_center_by_delta_mass_Msun > 0 .and. &
3252 : .not. restart) then
3253 0 : write(*, 1) 'remove_initial_center_by_delta_mass_Msun', &
3254 0 : s% job% remove_initial_center_by_delta_mass_Msun
3255 : call star_remove_center_by_mass_gm(id, &
3256 0 : s% M_center + s% job% remove_initial_center_by_delta_mass_Msun*Msun, ierr)
3257 0 : if (failed('star_remove_center_by_mass',ierr)) return
3258 : end if
3259 :
3260 1 : if (s% job% remove_center_by_mass_gm > s% M_center .and. &
3261 : s% job% remove_center_by_mass_gm < s% m(1)) then
3262 0 : write(*, 1) 'remove_center_by_mass_gm', &
3263 0 : s% job% remove_center_by_mass_gm
3264 : call star_remove_center_by_mass_gm( &
3265 0 : id, s% job% remove_center_by_mass_gm, ierr)
3266 0 : if (failed('star_remove_center_by_mass_gm',ierr)) return
3267 : end if
3268 :
3269 : if (s% job% remove_initial_center_by_mass_gm > s% M_center .and. &
3270 1 : s% job% remove_initial_center_by_mass_gm < s% m(1) .and. .not. restart) then
3271 0 : write(*, 1) 'remove_initial_center_by_mass_gm', &
3272 0 : s% job% remove_initial_center_by_mass_gm
3273 : call star_remove_center_by_mass_gm( &
3274 0 : id, s% job% remove_initial_center_by_mass_gm, ierr)
3275 0 : if (failed('star_remove_center_by_mass_gm',ierr)) return
3276 : end if
3277 :
3278 1 : if (s% job% remove_center_by_mass_Msun > s% M_center/Msun .and. &
3279 : s% job% remove_center_by_mass_Msun < s% m(1)/Msun) then
3280 0 : write(*, 1) 'remove_center_by_mass_Msun', &
3281 0 : s% job% remove_center_by_mass_Msun
3282 : call star_remove_center_by_mass_gm( &
3283 0 : id, s% job% remove_center_by_mass_Msun*Msun, ierr)
3284 0 : if (failed('star_remove_center_by_mass_Msun',ierr)) return
3285 : end if
3286 :
3287 : if (s% job% remove_initial_center_by_mass_Msun > s% M_center/Msun .and. &
3288 1 : s% job% remove_initial_center_by_mass_Msun < s% m(1)/Msun .and. &
3289 : .not. restart) then
3290 0 : write(*, 1) 'remove_initial_center_by_mass_Msun', &
3291 0 : s% job% remove_initial_center_by_mass_Msun
3292 : call star_remove_center_by_mass_gm( &
3293 0 : id, s% job% remove_initial_center_by_mass_Msun*Msun, ierr)
3294 0 : if (failed('star_remove_center_by_mass_Msun',ierr)) return
3295 : end if
3296 :
3297 1 : if (s% job% remove_center_by_radius_Rsun > s% R_center/Rsun .and. &
3298 : s% job% remove_center_by_radius_Rsun < s% r(1)/Rsun) then
3299 0 : write(*, 1) 'remove_center_by_radius_Rsun', &
3300 0 : s% job% remove_center_by_radius_Rsun
3301 : call star_remove_center_by_radius_cm( &
3302 0 : id, s% job% remove_center_by_radius_Rsun*Rsun, ierr)
3303 0 : if (failed('star_remove_center_by_radius_Rsun',ierr)) return
3304 : end if
3305 :
3306 : if (s% job% remove_initial_center_by_radius_Rsun > s% R_center/Rsun .and. &
3307 1 : s% job% remove_initial_center_by_radius_Rsun < s% r(1)/Rsun .and. &
3308 : .not. restart) then
3309 0 : write(*, 1) 'remove_initial_center_by_radius_Rsun', &
3310 0 : s% job% remove_initial_center_by_radius_Rsun
3311 : call star_remove_center_by_radius_cm( &
3312 0 : id, s% job% remove_initial_center_by_radius_Rsun*Rsun, ierr)
3313 0 : if (failed('star_remove_center_by_radius_Rsun',ierr)) return
3314 : end if
3315 :
3316 1 : if (s% job% remove_initial_center_at_cell_k > 0 .and. .not. restart .and. &
3317 : s% job% remove_initial_center_at_cell_k <= s% nz) then
3318 0 : write(*, 2) 'remove_initial_center_at_cell_k', s% job% remove_initial_center_at_cell_k
3319 : call star_remove_center_at_cell_k( &
3320 0 : id, s% job% remove_initial_center_at_cell_k, ierr)
3321 0 : if (failed('star_remove_center_at_cell_k',ierr)) return
3322 : end if
3323 :
3324 1 : if (s% job% remove_center_at_cell_k > 0 .and. &
3325 : s% job% remove_center_at_cell_k <= s% nz) then
3326 0 : write(*, 2) 'remove_center_at_cell_k', s% job% remove_center_at_cell_k
3327 0 : call star_remove_center_at_cell_k(id, s% job% remove_center_at_cell_k, ierr)
3328 0 : if (failed('star_remove_center_at_cell_k',ierr)) return
3329 : end if
3330 :
3331 : end subroutine do_remove_center
3332 :
3333 :
3334 1 : subroutine do_remove_initial_surface(id,s,restart,ierr)
3335 : integer, intent(in) :: id
3336 : type (star_info), pointer :: s
3337 : logical, intent(in) :: restart
3338 : integer, intent(out) :: ierr
3339 :
3340 : include 'formats'
3341 :
3342 1 : ierr = 0
3343 :
3344 1 : if (s% job% remove_initial_surface_at_he_core_boundary > 0 .and. .not. restart) then
3345 0 : write(*, 1) 'remove_initial_surface_at_he_core_boundary', &
3346 0 : s% job% remove_initial_surface_at_he_core_boundary
3347 : call star_remove_surface_at_he_core_boundary( &
3348 0 : id, s% job% remove_initial_surface_at_he_core_boundary, ierr)
3349 0 : if (failed('star_remove_surface_at_he_core_boundary',ierr)) return
3350 : end if
3351 :
3352 1 : if (s% job% remove_initial_surface_by_optical_depth > 0 .and. .not. restart) then
3353 0 : write(*, 1) 'remove_initial_surface_by_optical_depth', &
3354 0 : s% job% remove_initial_surface_by_optical_depth
3355 : call star_remove_surface_by_optical_depth( &
3356 0 : id, s% job% remove_initial_surface_by_optical_depth, ierr)
3357 0 : if (failed('star_remove_surface_by_optical_depth',ierr)) return
3358 : end if
3359 :
3360 1 : if (s% job% remove_initial_surface_by_density > 0 .and. .not. restart) then
3361 0 : write(*, 1) 'remove_initial_surface_by_density', &
3362 0 : s% job% remove_initial_surface_by_density
3363 : call star_remove_surface_by_density( &
3364 0 : id, s% job% remove_initial_surface_by_density, ierr)
3365 0 : if (failed('star_remove_surface_by_density',ierr)) return
3366 : end if
3367 :
3368 1 : if (s% job% remove_initial_surface_by_pressure > 0 .and. .not. restart) then
3369 0 : write(*, 1) 'remove_initial_surface_by_pressure', &
3370 0 : s% job% remove_initial_surface_by_pressure
3371 : call star_remove_surface_by_pressure( &
3372 0 : id, s% job% remove_initial_surface_by_pressure, ierr)
3373 0 : if (failed('star_remove_surface_by_pressure',ierr)) return
3374 : end if
3375 :
3376 : if (s% job% remove_initial_surface_by_radius_cm > s% R_center .and. &
3377 1 : s% job% remove_initial_surface_by_radius_cm < s% r(1) .and. .not. restart) then
3378 0 : write(*, 1) 'remove_initial_surface_by_radius_cm', &
3379 0 : s% job% remove_initial_surface_by_radius_cm
3380 : call star_remove_surface_by_radius_cm( &
3381 0 : id, s% job% remove_initial_surface_by_radius_cm, ierr)
3382 0 : if (failed('star_remove_surface_by_radius_cm',ierr)) return
3383 : end if
3384 :
3385 : if (s% job% remove_initial_surface_by_mass_fraction_q > 0d0 .and. &
3386 : s% job% remove_initial_surface_by_mass_fraction_q < 1d0 &
3387 1 : .and. .not. restart) then
3388 0 : write(*, 1) 'remove_initial_surface_by_mass_fraction_q', &
3389 0 : s% job% remove_initial_surface_by_mass_fraction_q
3390 : call star_remove_surface_by_mass_fraction_q( &
3391 0 : id, s% job% remove_initial_surface_by_mass_fraction_q, ierr)
3392 0 : if (failed('star_remove_initial_surface_by_mass_fraction_q',ierr)) return
3393 : end if
3394 :
3395 : if (s% job% remove_initial_surface_by_mass_gm > s% M_center .and. &
3396 1 : s% job% remove_initial_surface_by_mass_gm < s% m(1) .and. .not. restart) then
3397 0 : write(*, 1) 'remove_initial_surface_by_mass_gm', &
3398 0 : s% job% remove_initial_surface_by_mass_gm
3399 : call star_remove_surface_by_mass_gm( &
3400 0 : id, s% job% remove_initial_surface_by_mass_gm, ierr)
3401 0 : if (failed('star_remove_surface_by_mass_gm',ierr)) return
3402 : end if
3403 :
3404 : if (s% job% remove_initial_surface_by_radius_Rsun > s% R_center/Rsun .and. &
3405 1 : s% job% remove_initial_surface_by_radius_Rsun < s% r(1)/Rsun .and. &
3406 : .not. restart) then
3407 0 : write(*, 1) 'remove_initial_surface_by_radius_Rsun', &
3408 0 : s% job% remove_initial_surface_by_radius_Rsun
3409 : call star_remove_surface_by_radius_cm( &
3410 0 : id, s% job% remove_initial_surface_by_radius_Rsun*Rsun, ierr)
3411 0 : if (failed('star_remove_surface_by_radius_Rsun',ierr)) return
3412 : end if
3413 :
3414 : if (s% job% remove_initial_surface_by_mass_Msun > s% M_center/Msun .and. &
3415 1 : s% job% remove_initial_surface_by_mass_Msun < s% m(1)/Msun .and. &
3416 : .not. restart) then
3417 0 : write(*, 1) 'remove_initial_surface_by_mass_Msun', &
3418 0 : s% job% remove_initial_surface_by_mass_Msun
3419 : call star_remove_surface_by_mass_gm( &
3420 0 : id, s% job% remove_initial_surface_by_mass_Msun*Msun, ierr)
3421 0 : if (failed('star_remove_surface_by_mass_Msun',ierr)) return
3422 : end if
3423 :
3424 1 : if (s% job% remove_initial_surface_by_v_surf_km_s > 0 .and. .not. restart) then
3425 0 : write(*, 2) 'remove_initial_surface_by_v_surf_km_s', &
3426 0 : s% job% remove_initial_surface_by_v_surf_km_s
3427 : call star_remove_surface_by_v_surf_km_s( &
3428 0 : id, s% job% remove_initial_surface_by_v_surf_km_s, ierr)
3429 0 : if (failed('star_remove_surface_by_v_surf_km_s',ierr)) return
3430 : end if
3431 :
3432 1 : if (s% job% remove_initial_surface_by_v_surf_div_cs > 0 .and. .not. restart) then
3433 0 : write(*, 2) 'remove_initial_surface_by_v_surf_div_cs', &
3434 0 : s% job% remove_initial_surface_by_v_surf_div_cs
3435 : call star_remove_surface_by_v_surf_div_cs( &
3436 0 : id, s% job% remove_initial_surface_by_v_surf_div_cs, ierr)
3437 0 : if (failed('star_remove_surface_by_v_surf_div_cs',ierr)) return
3438 : end if
3439 :
3440 1 : if (s% job% remove_initial_surface_by_v_surf_div_v_escape > 0 .and. .not. restart) then
3441 0 : write(*, 2) 'remove_initial_surface_by_v_surf_div_v_escape', &
3442 0 : s% job% remove_initial_surface_by_v_surf_div_v_escape
3443 : call star_remove_surface_by_v_surf_div_v_escape( &
3444 0 : id, s% job% remove_initial_surface_by_v_surf_div_v_escape, ierr)
3445 0 : if (failed('star_remove_surface_by_v_surf_div_v_escape',ierr)) return
3446 : end if
3447 :
3448 1 : if (s% job% remove_initial_surface_at_cell_k > 0 .and. .not. restart .and. &
3449 : s% job% remove_initial_surface_at_cell_k <= s% nz) then
3450 0 : write(*, 2) 'remove_initial_surface_at_cell_k', s% job% remove_initial_surface_at_cell_k
3451 : call star_remove_surface_at_cell_k( &
3452 0 : id, s% job% remove_initial_surface_at_cell_k, ierr)
3453 0 : if (failed('star_remove_surface_at_cell_k',ierr)) return
3454 : end if
3455 :
3456 : end subroutine do_remove_initial_surface
3457 :
3458 :
3459 12 : subroutine do_remove_surface(id,s,ierr)
3460 : integer, intent(in) :: id
3461 : type (star_info), pointer :: s
3462 : integer, intent(out) :: ierr
3463 :
3464 : include 'formats'
3465 :
3466 12 : ierr = 0
3467 :
3468 12 : if (s% job% remove_surface_at_he_core_boundary > 0) then
3469 : !write(*, 1) 'remove_surface_at_he_core_boundary', s% job% remove_surface_at_he_core_boundary
3470 : call star_remove_surface_at_he_core_boundary( &
3471 0 : id, s% job% remove_surface_at_he_core_boundary, ierr)
3472 0 : if (failed('star_remove_surface_at_he_core_boundary',ierr)) return
3473 : end if
3474 :
3475 12 : if (s% job% remove_surface_by_optical_depth > 0) then
3476 : !write(*, 1) 'remove_surface_by_optical_depth', s% job% remove_surface_by_optical_depth
3477 : call star_remove_surface_by_optical_depth( &
3478 0 : id, s% job% remove_surface_by_optical_depth, ierr)
3479 0 : if (failed('star_remove_surface_by_optical_depth',ierr)) return
3480 : end if
3481 :
3482 12 : if (s% job% remove_surface_by_density > 0) then
3483 : !write(*, 1) 'remove_surface_by_density', s% job% remove_surface_by_density
3484 : call star_remove_surface_by_density( &
3485 0 : id, s% job% remove_surface_by_density, ierr)
3486 0 : if (failed('star_remove_surface_by_density',ierr)) return
3487 : end if
3488 :
3489 12 : if (s% job% remove_surface_by_pressure > 0) then
3490 : !write(*, 1) 'remove_surface_by_pressure', s% job% remove_surface_by_pressure
3491 : call star_remove_surface_by_pressure( &
3492 0 : id, s% job% remove_surface_by_pressure, ierr)
3493 0 : if (failed('star_remove_surface_by_pressure',ierr)) return
3494 : end if
3495 :
3496 12 : if (s% job% remove_surface_by_radius_cm > s% R_center .and. &
3497 : s% job% remove_surface_by_radius_cm < s% r(1)) then
3498 : !write(*, 1) 'remove_surface_by_radius_cm', s% job% remove_surface_by_radius_cm
3499 : call star_remove_surface_by_radius_cm( &
3500 0 : id, s% job% remove_surface_by_radius_cm, ierr)
3501 0 : if (failed('star_remove_surface_by_radius_cm',ierr)) return
3502 : end if
3503 :
3504 12 : if (s% job% remove_surface_by_mass_fraction_q > 0d0 .and. &
3505 : s% job% remove_surface_by_mass_fraction_q < 1d0) then
3506 : !write(*, 1) 'remove_surface_by_mass_fraction_q', &
3507 : ! s% job% remove_surface_by_mass_fraction_q
3508 : call star_remove_surface_by_mass_fraction_q( &
3509 0 : id, s% job% remove_surface_by_mass_fraction_q, ierr)
3510 0 : if (failed('star_remove_surface_by_mass_fraction_q',ierr)) return
3511 : end if
3512 :
3513 12 : if (s% job% remove_surface_by_mass_gm > s% M_center .and. &
3514 : s% job% remove_surface_by_mass_gm < s% m(1)) then
3515 : !write(*, 1) 'remove_surface_by_mass_gm', &
3516 : ! s% job% remove_surface_by_mass_gm
3517 : call star_remove_surface_by_mass_gm( &
3518 0 : id, s% job% remove_surface_by_mass_gm, ierr)
3519 0 : if (failed('star_remove_surface_by_mass_gm',ierr)) return
3520 : end if
3521 :
3522 12 : if (s% job% remove_surface_by_radius_Rsun > s% R_center/Rsun .and. &
3523 : s% job% remove_surface_by_radius_Rsun < s% r(1)/Rsun) then
3524 : !write(*, 1) 'remove_surface_by_radius_Rsun', &
3525 : ! s% job% remove_surface_by_radius_Rsun
3526 : call star_remove_surface_by_radius_cm( &
3527 0 : id, s% job% remove_surface_by_radius_Rsun*Rsun, ierr)
3528 0 : if (failed('star_remove_surface_by_radius_Rsun',ierr)) return
3529 : end if
3530 :
3531 12 : if (s% job% remove_surface_by_mass_Msun > s% M_center/Msun .and. &
3532 : s% job% remove_surface_by_mass_Msun < s% m(1)/Msun) then
3533 : !write(*, 1) 'remove_surface_by_mass_Msun', &
3534 : ! s% job% remove_surface_by_mass_Msun
3535 : call star_remove_surface_by_mass_gm( &
3536 0 : id, s% job% remove_surface_by_mass_Msun*Msun, ierr)
3537 0 : if (failed('star_remove_surface_by_mass_Msun',ierr)) return
3538 : end if
3539 :
3540 12 : if (s% job% remove_surface_by_v_surf_km_s > 0) then
3541 : !write(*, 2) 'remove_surface_by_v_surf_km_s', s% job% remove_surface_by_v_surf_km_s
3542 0 : call star_remove_surface_by_v_surf_km_s(id, s% job% remove_surface_by_v_surf_km_s, ierr)
3543 0 : if (failed('star_remove_surface_by_v_surf_km_s',ierr)) return
3544 : end if
3545 :
3546 12 : if (s% job% remove_surface_by_v_surf_div_cs > 0) then
3547 : !write(*, 1) 'remove_surface_by_v_surf_div_cs', s% job% remove_surface_by_v_surf_div_cs
3548 0 : call star_remove_surface_by_v_surf_div_cs(id, s% job% remove_surface_by_v_surf_div_cs, ierr)
3549 0 : if (failed('star_remove_surface_by_v_surf_div_cs',ierr)) return
3550 : end if
3551 :
3552 12 : if (s% job% remove_surface_by_v_surf_div_v_escape > 0) then
3553 : !write(*, 2) 'remove_surface_by_v_surf_div_v_escape', s% job% remove_surface_by_v_surf_div_v_escape
3554 0 : call star_remove_surface_by_v_surf_div_v_escape(id, s% job% remove_surface_by_v_surf_div_v_escape, ierr)
3555 0 : if (failed('star_remove_surface_by_v_surf_div_v_escape',ierr)) return
3556 : end if
3557 :
3558 12 : if (s% job% remove_surface_at_cell_k > 0 .and. &
3559 : s% job% remove_surface_at_cell_k <= s% nz) then
3560 : !write(*, 2) 'remove_surface_at_cell_k', s% job% remove_surface_at_cell_k
3561 0 : call star_remove_surface_at_cell_k(id, s% job% remove_surface_at_cell_k, ierr)
3562 0 : if (failed('star_remove_surface_at_cell_k',ierr)) return
3563 : end if
3564 :
3565 : end subroutine do_remove_surface
3566 :
3567 :
3568 2 : subroutine resolve_inlist_fname(inlist_out,inlist_opt)
3569 :
3570 : use ISO_FORTRAN_ENV
3571 :
3572 : character(len=*),intent(out) :: inlist_out
3573 : character(len=*),optional :: inlist_opt
3574 :
3575 : integer :: status
3576 :
3577 : ! initialize inlist_out as empty
3578 2 : inlist_out = ''
3579 :
3580 2 : if (.not. MESA_INLIST_RESOLVED) then
3581 2 : if (COMMAND_ARGUMENT_COUNT() >= 1) then
3582 :
3583 : ! Get filename from the first command-line argument
3584 :
3585 0 : call GET_COMMAND_ARGUMENT(1, inlist_out, STATUS=status)
3586 0 : if (status /= 0) inlist_out = ''
3587 :
3588 : else
3589 :
3590 : ! Get filename from the MESA_INLIST environment variable
3591 :
3592 2 : call GET_ENVIRONMENT_VARIABLE('MESA_INLIST', inlist_out, STATUS=status)
3593 2 : if (status /= 0) inlist_out = ''
3594 :
3595 : end if
3596 : end if
3597 :
3598 2 : if (inlist_out == '') then
3599 :
3600 2 : if (PRESENT(inlist_opt)) then
3601 1 : inlist_out = inlist_opt
3602 : else
3603 1 : inlist_out = 'inlist'
3604 : end if
3605 :
3606 : end if
3607 :
3608 2 : return
3609 :
3610 2 : end subroutine resolve_inlist_fname
3611 :
3612 :
3613 1 : subroutine add_fpe_checks(id, s, ierr)
3614 : integer, intent(in) :: id
3615 : type (star_info), pointer :: s
3616 : integer, intent(out) :: ierr
3617 :
3618 : character(len=1) :: fpe_check
3619 : integer :: status
3620 :
3621 : include 'formats'
3622 :
3623 1 : ierr = 0
3624 :
3625 1 : call GET_ENVIRONMENT_VARIABLE('MESA_FPE_CHECKS_ON', fpe_check, STATUS=status)
3626 1 : if (status /= 0) return
3627 :
3628 0 : if (fpe_check(1:1)=="1") then
3629 0 : write(*,*) "FPE checking is on"
3630 0 : s% fill_arrays_with_nans = .true.
3631 : end if
3632 :
3633 : end subroutine add_fpe_checks
3634 :
3635 :
3636 1 : subroutine multiply_tolerances(id, s, ierr)
3637 : integer, intent(in) :: id
3638 : type (star_info), pointer :: s
3639 : integer, intent(out) :: ierr
3640 : integer :: status
3641 :
3642 : real(dp), save :: test_suite_res_factor = 1
3643 : character(len=20) :: test_suite_resolution_factor_str
3644 :
3645 : include 'formats'
3646 :
3647 1 : ierr = 0
3648 : call GET_ENVIRONMENT_VARIABLE('MESA_TEST_SUITE_RESOLUTION_FACTOR', &
3649 1 : test_suite_resolution_factor_str, STATUS=status)
3650 1 : if (status /= 0) return
3651 :
3652 0 : if (test_suite_resolution_factor_str /= "") then
3653 0 : read(test_suite_resolution_factor_str, *) test_suite_res_factor
3654 0 : write(*,*) ""
3655 0 : write(*,*) "***"
3656 0 : write(*,*) "MESA_TEST_SUITE_RESOLUTION_FACTOR set to", test_suite_res_factor
3657 0 : write(*,*) "***"
3658 0 : write(*,*) "Warning: This environment variable is for testing purposes"
3659 0 : write(*,*) " and should be set to 1 during normal MESA use."
3660 0 : write(*,*) "***"
3661 0 : write(*,*) "Multiplying mesh_delta_coeff and time_delta_coeff by this factor,"
3662 0 : write(*,*) "and max_model_number by its inverse twice:"
3663 0 : write(*,*) ""
3664 0 : write(*,*) " old mesh_delta_coeff = ", s% mesh_delta_coeff
3665 0 : s% mesh_delta_coeff = test_suite_res_factor * s% mesh_delta_coeff
3666 0 : write(*,*) " new mesh_delta_coeff = ", s% mesh_delta_coeff
3667 0 : write(*,*) ""
3668 0 : write(*,*) " old time_delta_coeff = ", s% time_delta_coeff
3669 0 : s% time_delta_coeff = test_suite_res_factor * s% time_delta_coeff
3670 0 : write(*,*) " new time_delta_coeff = ", s% time_delta_coeff
3671 0 : write(*,*) ""
3672 0 : write(*,*) " old max_model_number = ", s% max_model_number
3673 0 : s% max_model_number = s% max_model_number / test_suite_res_factor / test_suite_res_factor
3674 0 : write(*,*) " new max_model_number = ", s% max_model_number
3675 0 : write(*,*) ""
3676 : end if
3677 :
3678 : end subroutine multiply_tolerances
3679 :
3680 :
3681 1 : subroutine pgstar_env_check(id, s, ierr)
3682 : integer, intent(in) :: id
3683 : type (star_info), pointer :: s
3684 : integer, intent(out) :: ierr
3685 : character(len=5) :: flag
3686 : integer :: status
3687 :
3688 : include 'formats'
3689 :
3690 1 : ierr = 0
3691 :
3692 1 : call get_environment_variable('MESA_FORCE_PGSTAR_FLAG', flag, STATUS=status)
3693 1 : if (status /= 0) return
3694 :
3695 0 : select case (trim(flag))
3696 : case ("TRUE", "true")
3697 0 : write(*,*) "PGSTAR forced on"
3698 0 : s% job% pgstar_flag = .true.
3699 : case ("FALSE", "false")
3700 0 : write(*,*) "PGSTAR forced off"
3701 0 : s% job% pgstar_flag = .false.
3702 : end select
3703 :
3704 : end subroutine pgstar_env_check
3705 :
3706 : end module run_star_support
|
