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