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