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