Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 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 do_one_utils
21 :
22 : use star_private_def
23 : use const_def, only: dp, ln10, secday, dayyer, clight, msun, rsun
24 : use utils_lib, only: is_bad
25 :
26 : implicit none
27 :
28 : private
29 : public :: do_one_check_model
30 : public :: write_terminal_header
31 : public :: do_bare_bones_check_model
32 : public :: do_check_limits
33 : public :: do_show_log_description
34 : public :: do_show_terminal_header
35 : public :: do_terminal_summary
36 :
37 : ! model log priorities
38 : integer, parameter :: delta_priority = 1
39 : integer, parameter :: phase_priority = 2
40 :
41 : contains
42 :
43 11 : logical function model_is_okay(s)
44 : type (star_info), pointer :: s
45 : ! for now, just check for valid number in the final dynamic timescale
46 : model_is_okay = ((s% dynamic_timescale - s% dynamic_timescale) == 0d0) &
47 : .and. ((s% dynamic_timescale + 1d0) > 1d0)
48 : end function model_is_okay
49 :
50 2 : subroutine set_save_profiles_info(s, model_priority)
51 : type (star_info), pointer :: s
52 : integer, intent(in) :: model_priority
53 2 : s% need_to_save_profiles_now = .true.
54 2 : s% save_profiles_model_priority = model_priority
55 2 : end subroutine set_save_profiles_info
56 :
57 :
58 2 : subroutine write_terminal_header(s)
59 : type (star_info), pointer :: s
60 2 : if (s% model_number <= s% recent_log_header) return
61 2 : if (s% just_wrote_terminal_header) return
62 2 : s% recent_log_header = s% model_number
63 2 : call do_show_terminal_header(s)
64 2 : s% just_wrote_terminal_header = .true.
65 : end subroutine write_terminal_header
66 :
67 :
68 0 : subroutine do_show_log_description(id, ierr)
69 : integer, intent(in) :: id
70 : integer, intent(out) :: ierr
71 : type (star_info), pointer :: s
72 : ierr = 0
73 0 : call get_star_ptr(id, s, ierr)
74 0 : if (ierr /= 0) return
75 0 : write(*,'(A)')
76 0 : write(*,'(a)') " The terminal output contains the following information"
77 0 : write(*,'(A)')
78 0 : write(*,'(a)') " 'step' is the number of steps since the start of the run,"
79 0 : write(*,'(a)') " 'lg_dt' is log10 timestep in years,"
80 0 : write(*,'(a)') " 'age_yr' is the simulated years since the start run,"
81 0 : write(*,'(a)') " 'lg_Tcntr' is log10 center temperature (K),"
82 0 : write(*,'(a)') " 'lg_Dcntr' is log10 center density (g/cm^3),"
83 0 : write(*,'(a)') " 'lg_Tmax' is log10 max temperature (K),"
84 0 : write(*,'(a)') " 'Teff' is the surface temperature (K),"
85 0 : write(*,'(a)') " 'lg_R' is log10 surface radius (Rsun),"
86 0 : write(*,'(a)') " 'lg_L' is log10 surface luminosity (Lsun),"
87 0 : write(*,'(a)') " 'lg_LH' is log10 total PP and CNO hydrogen burning power (Lsun),"
88 0 : write(*,'(a)') " 'lg_L3a' is log10 total triple-alpha helium burning power (Lsun),"
89 0 : write(*,'(a)') " 'lg_gsurf' is log10 surface gravity,"
90 0 : write(*,'(a)') " 'lg_Lnuc' is log10 nuclear power (Lsun),"
91 0 : write(*,'(a)') " 'lg_Lneu' is log10 total neutrino power (Lsun),"
92 0 : write(*,'(a)') " 'lg_Lphoto' is log10 total photodisintegration (Lsun),"
93 0 : write(*,'(a)') " 'Mass' is the total stellar baryonic mass (Msun),"
94 0 : write(*,'(a)') " 'lg_Mdot' is log10 magnitude of rate of change of mass (Msun/year),"
95 0 : write(*,'(a)') " 'lg_Dsurf' is log10 surface density (g/cm^3),"
96 0 : write(*,'(a)') " 'H_env' is the amount of mass where H is dominant,"
97 0 : write(*,'(a,e9.2)') " 'He_core' is the largest mass where He is dominant."
98 0 : write(*,'(a,e9.2)') " 'CO_core' is the largest mass where CO is dominant."
99 0 : write(*,'(a)') " 'H_cntr' is the center H1 mass fraction,"
100 0 : write(*,'(a)') " 'He_cntr' is the center He4 mass fraction,"
101 0 : write(*,'(a)') " 'C_cntr' is the center C12 mass fraction,"
102 0 : write(*,'(a)') " 'N_cntr' is the center N14 mass fraction,"
103 0 : write(*,'(a)') " 'O_cntr' is the center O16 mass fraction,"
104 0 : write(*,'(a)') " 'Ne_cntr' is the center Ne20 mass fraction,"
105 0 : write(*,'(a)') " 'gam_cntr' is the center plasma interaction parameter,"
106 0 : write(*,'(a)') " 'eta_cntr' is the center electron degeneracy parameter,"
107 0 : write(*,'(a)') " 'zones' is the number of zones in the current model,"
108 0 : write(*,'(a)') " 'iters' is the number of solver iterations for the current step,"
109 0 : write(*,'(a)') " 'retry' is the number of step retries required during the run,"
110 0 : write(*,'(a)') " 'dt_limit' is an indication of what limited the timestep."
111 0 : write(*,'(A)')
112 0 : write(*,'(a)') " All this and more are saved in the LOGS directory during the run."
113 : end subroutine do_show_log_description
114 :
115 :
116 3 : subroutine do_show_terminal_header(s)
117 : type (star_info), pointer :: s
118 : integer :: ierr, io
119 3 : call output_terminal_header(s,terminal_iounit)
120 3 : if (len_trim(s% extra_terminal_output_file) > 0) then
121 0 : ierr = 0
122 : open(newunit=io, file=trim(s% extra_terminal_output_file), &
123 0 : action='write', position='append', iostat=ierr)
124 0 : if (ierr == 0) then
125 0 : call output_terminal_header(s,io)
126 0 : close(io)
127 : else
128 : write(*,*) 'failed to open extra_terminal_output_file ' // &
129 0 : trim(s% extra_terminal_output_file)
130 : end if
131 : end if
132 3 : end subroutine do_show_terminal_header
133 :
134 :
135 3 : subroutine output_terminal_header(s,io)
136 : use chem_def, only: isi28
137 : type (star_info), pointer :: s
138 : integer, intent(in) :: io
139 : character (len=5) :: iters
140 3 : iters = 'iters'
141 :
142 : include 'formats'
143 :
144 : write(io,'(a)') &
145 : '_______________________________________________________________________' // &
146 3 : '___________________________________________________________________________'
147 3 : write(io,'(A)')
148 : write(io,'(a)') &
149 : ' step lg_Tmax Teff lg_LH lg_Lnuc_tot Mass ' // &
150 3 : 'H_rich H_cntr N_cntr Y_surf eta_cntr zones retry'
151 :
152 : ! note that if the age is in days, then the timestep is automatically in seconds.
153 3 : if (trim(s% terminal_show_timestep_units) == 'seconds' .or. &
154 : trim(s% terminal_show_timestep_units) == 'secs') then
155 0 : if (s% terminal_show_log_dt) then
156 0 : write(io,'(a)',advance='no') ' lg_dt_sec'
157 : else
158 0 : write(io,'(a)',advance='no') ' dt_sec'
159 : end if
160 3 : else if (trim(s% terminal_show_timestep_units) == 'days') then
161 0 : if (s% terminal_show_log_dt) then
162 0 : write(io,'(a)',advance='no') ' lg_dt_days'
163 : else
164 0 : write(io,'(a)',advance='no') ' dt_days'
165 : end if
166 3 : else if (trim(s% terminal_show_timestep_units) == 'years' .or. &
167 : trim(s% terminal_show_timestep_units) == 'yrs') then
168 3 : if (s% terminal_show_log_dt) then
169 3 : write(io,'(a)',advance='no') ' lg_dt_yrs'
170 : else
171 0 : write(io,'(a)',advance='no') ' dt_yrs'
172 : end if
173 : else
174 0 : write(*,*) 'unrecognized option for terminal_show_timestep_units ' // trim(s% terminal_show_timestep_units)
175 0 : return
176 : end if
177 :
178 3 : if (s% initial_z >= 1d-5) then
179 : write(io,'(a)') &
180 : ' lg_Tcntr lg_R lg_L3a lg_Lneu lg_Mdot ' // &
181 3 : 'He_core He_cntr O_cntr Z_surf gam_cntr ' // iters // ' '
182 : else
183 : write(io,'(a)') &
184 : ' lg_Tcntr lg_R lg_L3a lg_Lneu lg_Mdot ' // &
185 0 : 'He_core He_cntr O_cntr lg_Z_surf gam_cntr ' // iters // ' '
186 : end if
187 :
188 3 : if (trim(s% terminal_show_age_units) == 'seconds' .or. &
189 : trim(s% terminal_show_age_units) == 'secs') then
190 0 : if (s% terminal_show_log_age) then
191 0 : write(io,'(a)',advance='no') 'lg_age_secs'
192 : else
193 0 : write(io,'(a)',advance='no') ' age_secs'
194 : end if
195 3 : else if (trim(s% terminal_show_age_units) == 'days') then
196 0 : if (s% terminal_show_log_age) then
197 0 : write(io,'(a)',advance='no') 'lg_age_days'
198 : else
199 0 : write(io,'(a)',advance='no') ' age_days'
200 : end if
201 3 : else if (trim(s% terminal_show_age_units) == 'years' .or. &
202 : trim(s% terminal_show_age_units) == 'yrs') then
203 3 : if (s% terminal_show_log_age) then
204 0 : write(io,'(a)',advance='no') ' lg_age_yrs'
205 : else
206 3 : write(io,'(a)',advance='no') ' age_yrs'
207 : end if
208 : else
209 0 : write(*,*) 'unrecognized option for terminal_show_age_units ' // trim(s% terminal_show_age_units)
210 0 : return
211 : end if
212 :
213 3 : if (s% net_iso(isi28) == 0) then
214 : write(io,'(a)') &
215 : ' lg_Dcntr lg_L lg_LZ lg_Lphoto lg_Dsurf ' // &
216 3 : 'CO_core C_cntr Ne_cntr Z_cntr v_div_cs dt_limit'
217 : else
218 : write(io,'(a)') &
219 : ' lg_Dcntr lg_L lg_LZ lg_Lphoto lg_Dsurf ' // &
220 0 : 'CO_core C_cntr Ne_cntr Si_cntr v_div_cs dt_limit'
221 : end if
222 : write(io,'(a)') &
223 : '_______________________________________________________________________' // &
224 3 : '___________________________________________________________________________'
225 3 : write(io,'(A)')
226 :
227 3 : end subroutine output_terminal_header
228 :
229 :
230 12 : subroutine do_terminal_summary(s)
231 : type (star_info), pointer :: s
232 : integer :: ierr, io
233 12 : call output_terminal_summary(s,terminal_iounit)
234 12 : if (len_trim(s% extra_terminal_output_file) > 0) then
235 0 : ierr = 0
236 : open(newunit=io, file=trim(s% extra_terminal_output_file), &
237 0 : action='write', position='append', iostat=ierr)
238 0 : if (ierr == 0) then
239 0 : call output_terminal_summary(s,io)
240 0 : close(io)
241 : else
242 : write(*,*) 'failed to open extra_terminal_output_file ' // &
243 0 : trim(s% extra_terminal_output_file)
244 : end if
245 : end if
246 3 : end subroutine do_terminal_summary
247 :
248 :
249 12 : subroutine output_terminal_summary(s,io)
250 : use num_def, only:banded_matrix_type
251 : use const_def, only:secyer
252 : use chem_def
253 : use rates_def, only: i_rate
254 : use star_utils, only:eval_current_y, eval_current_z
255 : type (star_info), pointer :: s
256 : integer, intent(in) :: io
257 :
258 12 : real(dp) :: time_step, age, dt, Xmax, v, vsurf_div_csound, tmp, &
259 12 : sum_Lnuc, sum_LH, sum_LHe, sum_Lz, sum_Lphoto
260 : integer :: model, ierr, nz, iters
261 : character (len=3) :: id_str
262 : character (len=32) :: why
263 : character (len=90) :: fmt, fmt1, fmt2, fmt3, fmt4, fmt5
264 :
265 : include 'formats'
266 :
267 12 : age = s% star_age ! in years
268 12 : if (trim(s% terminal_show_age_units) == 'seconds' .or. &
269 : trim(s% terminal_show_age_units) == 'secs') then
270 0 : age = age*secyer
271 12 : else if (trim(s% terminal_show_age_units) == 'days') then
272 0 : age = age*dayyer
273 : end if
274 :
275 12 : time_step = s% time_step ! in years
276 12 : if (trim(s% terminal_show_timestep_units) == 'seconds' .or. &
277 : trim(s% terminal_show_timestep_units) == 'secs') then
278 0 : time_step = time_step*secyer
279 12 : else if (trim(s% terminal_show_timestep_units) == 'days') then
280 0 : time_step = time_step*dayyer
281 : end if
282 :
283 12 : if (s% terminal_show_log_age) age = safe_log10(age)
284 12 : if (s% terminal_show_log_dt) time_step = safe_log10(time_step)
285 :
286 12 : model = s% model_number
287 12 : nz = s% nz
288 :
289 12 : ierr = 0
290 :
291 12 : Xmax = dot_product(s% dq(1:nz), s% xa(s% species,1:nz))
292 :
293 12 : if (s% u_flag) then
294 0 : v = s% u(1)
295 12 : else if (s% v_flag) then
296 0 : v = s% v(1)
297 : else
298 : v = 0d0
299 : end if
300 12 : vsurf_div_csound = v / s% csound(1)
301 :
302 12 : dt = s% time_step*secyer
303 :
304 12 : sum_Lnuc = s% power_nuc_burn
305 12 : sum_LH = s% power_h_burn
306 12 : sum_LHe = s% power_he_burn
307 12 : sum_Lphoto = abs(s% power_photo)
308 12 : sum_Lz = s% power_z_burn
309 :
310 12 : if (how_many_allocated_star_ids() == 1) then
311 12 : id_str = ''
312 : else
313 0 : write(id_str,'(i3)') s% id
314 : end if
315 :
316 12 : fmt1 = '(a3,i8,f11.6,'
317 :
318 12 : if (s% Teff < 1d4) then
319 0 : fmt2 = 'f11.3,'
320 : else
321 12 : fmt2 = '1pe11.3,0p,'
322 : end if
323 :
324 12 : if (s% star_mass >= 1d2) then
325 0 : fmt3 = '2f11.6,2(1pe11.3),0p,'
326 : else
327 12 : fmt3 = '4f11.6,'
328 : end if
329 :
330 12 : if (s% eta(s% nz) >= 1d3) then
331 0 : fmt4 = '3f11.6,e11.3,'
332 : else
333 12 : fmt4 = '3f11.6,f11.6,'
334 : end if
335 12 : fmt5 = '2i7)'
336 :
337 12 : fmt = trim(fmt1) // trim(fmt2) // trim(fmt3) // trim(fmt4) // trim(fmt5)
338 : !write(*,*) 'fmt line1 ' // trim(fmt)
339 : write(io,fmt=fmt) &
340 12 : id_str, model, &
341 12 : s% log_max_temperature, & ! fmt1
342 12 : s% Teff, & ! fmt2
343 12 : safe_log10(sum_LH), & ! fmt3
344 12 : safe_log10(sum_Lnuc), &
345 12 : s% star_mass, &
346 12 : s% star_mass - max(s% he_core_mass, s% co_core_mass), &
347 12 : s% center_h1, & ! fmt4
348 12 : s% center_n14, &
349 12 : s% surface_he3 + s% surface_he4, &
350 12 : s% eta(s% nz), &
351 12 : s% nz, & ! fmt5
352 24 : s% num_retries
353 :
354 12 : tmp = max(0d0, min(1d0, 1 - (s% surface_h1 + s% surface_he3 + s% surface_he4)))
355 12 : if (s% initial_z >= 1d-5) then
356 12 : fmt1 = '(1pe11.4, 0p, 9f11.6, '
357 : else
358 0 : tmp = safe_log10(tmp)
359 0 : fmt1 = '(1pe11.4, 0p, 8f11.6, e11.2, '
360 : end if
361 12 : if (s% gam(s% nz) >= 1d3) then
362 0 : fmt2 = 'e11.3, '
363 : else
364 12 : fmt2 = 'f11.6, '
365 : end if
366 12 : fmt3 = ' i7)'
367 12 : fmt = trim(fmt1) // trim(fmt2) // trim(fmt3)
368 : !write(*,*) 'fmt line2 ' // trim(fmt)
369 12 : iters = s% num_solver_iterations
370 : write(io,fmt=fmt) &
371 12 : time_step, &
372 12 : s% log_center_temperature, &
373 12 : s% log_surface_radius, &
374 12 : safe_log10(sum_LHe), &
375 12 : safe_log10(abs(s% power_neutrinos)), &
376 12 : safe_log10(abs(s% star_mdot)), &
377 12 : s% he_core_mass, &
378 12 : s% center_he3 + s% center_he4, &
379 12 : s% center_o16, &
380 12 : tmp, &
381 12 : s% gam(s% nz), &
382 24 : iters
383 :
384 12 : if (s% why_Tlim < 0) then
385 0 : why = ''
386 12 : else if (s% why_Tlim == 0) then
387 1 : why = 'initial dt'
388 : else
389 11 : why = dt_why_str(min(numTlim,s% why_Tlim))
390 : if (s% why_Tlim == Tlim_dX .and. s% Tlim_dX_species > 0 &
391 11 : .and. s% Tlim_dX_species <= s% species) then
392 : why = trim(dt_why_str(s% why_Tlim)) // ' ' // &
393 0 : trim(chem_isos% name(s% chem_id(s% Tlim_dX_species)))
394 : else if (s% why_Tlim == Tlim_dX_div_X .and. s% Tlim_dX_div_X_species > 0 &
395 11 : .and. s% Tlim_dX_div_X_species <= s% species) then
396 : why = trim(dt_why_str(s% why_Tlim)) // ' ' // &
397 0 : trim(chem_isos% name(s% chem_id(s% Tlim_dX_div_X_species)))
398 : else if (s% why_Tlim == Tlim_dX_nuc_drop .and. s% dX_nuc_drop_max_j > 0 &
399 11 : .and. s% dX_nuc_drop_max_j <= s% species) then
400 : why = trim(dt_why_str(s% why_Tlim)) // ' ' // &
401 0 : trim(chem_isos% name(s% chem_id(s% dX_nuc_drop_max_j)))
402 11 : else if (s% why_Tlim == Tlim_dlgL_nuc_cat) then
403 : if (s% Tlim_dlgL_nuc_category > 0 &
404 0 : .and. s% Tlim_dlgL_nuc_category <= num_categories ) then
405 0 : why = trim(category_name(s% Tlim_dlgL_nuc_category))
406 : else
407 0 : why = '???'
408 : end if
409 : end if
410 : end if
411 :
412 12 : if (s% net_iso(isi28) == 0) then
413 12 : tmp = 1 - (s% center_h1 + s% center_he3 + s% center_he4)
414 12 : fmt = '(1pe11.4, 0p, 5f11.6, 0p4f11.6, 0p, e11.3, a14)'
415 : !write(*,*) 'fmt line3 ' // trim(fmt)
416 : write(io,fmt) &
417 12 : age, &
418 12 : s% log_center_density, &
419 12 : s% log_surface_luminosity, &
420 12 : safe_log10(sum_Lz), &
421 12 : safe_log10(abs(s% power_photo)), &
422 12 : s% lnd(1)/ln10, &
423 12 : s% co_core_mass, &
424 12 : s% center_c12, &
425 12 : s% center_ne20, &
426 12 : tmp, &
427 12 : vsurf_div_csound, &
428 24 : trim(why)
429 : else
430 0 : tmp = s% center_si28
431 0 : fmt = '(1pe11.4, 0p, 5f11.6, 0p4f11.6, 0p, e11.3, a14)'
432 : !write(*,*) 'fmt line3 ' // trim(fmt)
433 : write(io,fmt) &
434 0 : age, &
435 0 : s% log_center_density, &
436 0 : s% log_surface_luminosity, &
437 0 : safe_log10(sum_Lz), &
438 0 : safe_log10(abs(s% power_photo)), &
439 0 : s% lnd(1)/ln10, &
440 0 : s% co_core_mass, &
441 0 : s% center_c12, &
442 0 : s% center_ne20, &
443 0 : tmp, &
444 0 : vsurf_div_csound, &
445 0 : trim(why)
446 : end if
447 :
448 12 : call show_trace_history_values(max(0, s% num_trace_history_values))
449 12 : write(io,'(A)')
450 :
451 12 : s% just_wrote_terminal_header = .false.
452 :
453 :
454 : contains
455 :
456 :
457 12 : subroutine show_trace_history_values(num)
458 : use history, only: do_get_data_for_history_columns, &
459 12 : get_history_specs, get_history_values, get1_hist_value
460 : integer, intent(in) :: num
461 12 : real(dp) :: values(num)
462 24 : integer :: int_values(num), specs(num)
463 24 : logical :: is_int_value(num)
464 12 : logical :: failed_to_find_value(num)
465 12 : real(dp) :: val
466 : integer :: i
467 : include 'formats'
468 12 : if (num == 0) return
469 : call do_get_data_for_history_columns(s, ierr)
470 0 : if (ierr /= 0) return
471 0 : call get_history_specs(s, num, s% trace_history_value_name, specs, .false.)
472 : call get_history_values( &
473 0 : s, num, specs, is_int_value, int_values, values, failed_to_find_value)
474 0 : do i = 1, num
475 0 : if (failed_to_find_value(i)) then
476 0 : if (.not. get1_hist_value(s, s% trace_history_value_name(i), val)) then
477 0 : write(*,*) 'failed to find history value ' // trim(s% trace_history_value_name(i))
478 0 : cycle
479 : end if
480 0 : values(i) = val
481 0 : if (is_bad(values(i))) then
482 0 : call mesa_error(__FILE__,__LINE__,'show_trace_history_values bad from get1_hist_value')
483 : end if
484 0 : else if (is_int_value(i)) then
485 : write(io,'(a40,i14)') &
486 0 : trim(s% trace_history_value_name(i)), int_values(i)
487 0 : cycle
488 : end if
489 0 : if ((values(i) == 0) .or. &
490 0 : (abs(values(i)) > 1d-4 .and. abs(values(i)) < 1d4)) then
491 : write(io,'(a40,99(f26.16))') &
492 0 : trim(s% trace_history_value_name(i)), values(i)
493 : else
494 : write(io,'(a40,99(1pd26.16))') &
495 0 : trim(s% trace_history_value_name(i)), values(i)
496 0 : if (is_bad(values(i))) then
497 0 : call mesa_error(__FILE__,__LINE__,'show_trace_history_values')
498 : end if
499 : end if
500 : end do
501 12 : end subroutine show_trace_history_values
502 :
503 :
504 : end subroutine output_terminal_summary
505 :
506 :
507 11 : logical function get_history_info(s, do_write)
508 : type (star_info), pointer :: s
509 : logical, intent(in) :: do_write
510 : integer :: model
511 : logical :: write_history, write_terminal
512 : include 'formats'
513 11 : model = s% model_number
514 11 : if (s% history_interval > 0) then
515 11 : write_history = (mod(model, s% history_interval) == 0) .or. do_write
516 : else
517 : write_history = .false.
518 : end if
519 11 : if (s% terminal_interval > 0) then
520 11 : write_terminal = (mod(model, s% terminal_interval) == 0) .or. do_write
521 : else
522 : write_terminal = .false.
523 : end if
524 11 : get_history_info = write_history .or. write_terminal
525 11 : if (.not. get_history_info) return
526 11 : if (s% write_header_frequency*s% terminal_interval > 0) then
527 : if ( mod(model, s% write_header_frequency*s% terminal_interval) == 0 &
528 11 : .and. .not. s% doing_first_model_of_run) then
529 1 : write(*,'(A)')
530 1 : call write_terminal_header(s)
531 : end if
532 : end if
533 11 : if (write_terminal) call do_terminal_summary(s)
534 11 : if (write_history) s% need_to_update_history_now = .true.
535 :
536 : end function get_history_info
537 :
538 :
539 0 : integer function do_bare_bones_check_model(id)
540 : integer, intent(in) :: id
541 : integer :: ierr
542 : logical :: logged
543 : type (star_info), pointer :: s
544 0 : call get_star_ptr(id, s, ierr)
545 0 : if (ierr /= 0) then
546 0 : do_bare_bones_check_model = terminate
547 : return
548 : end if
549 0 : logged = get_history_info( s, .false. )
550 0 : do_bare_bones_check_model = keep_going
551 0 : end function do_bare_bones_check_model
552 :
553 :
554 : subroutine save_profile(id, ierr)
555 : integer, intent(in) :: id
556 : integer, intent(out) :: ierr
557 : type (star_info), pointer :: s
558 : ierr = 0
559 : call get_star_ptr(id, s, ierr)
560 : if (ierr /= 0) return
561 : call set_save_profiles_info(s, phase_priority)
562 : end subroutine save_profile
563 :
564 :
565 11 : integer function do_check_limits(id)
566 : use rates_def
567 : use chem_def
568 : use chem_lib, only: chem_get_iso_id
569 : use star_utils
570 : integer, intent(in) :: id
571 : type (star_info), pointer :: s
572 : integer :: ierr, i, j, k, cid, k_omega, nz, max_abs_vel_loc, &
573 : period_number, max_period_number
574 11 : real(dp) :: log_surface_gravity, log_surface_temperature, log_surface_density, &
575 11 : log_surface_pressure, v_div_csound_max, remnant_mass, ejecta_mass, &
576 11 : power_nuc_burn, power_h_burn, power_he_burn, power_z_burn, logQ, max_logQ, min_logQ, &
577 11 : envelope_fraction_left, avg_x, v_surf, csound_surf, delta_nu, v_surf_div_v_esc, &
578 11 : peak_burn_vconv_div_cs, min_pgas_div_p, v_surf_div_v_kh, GREKM_avg_abs, &
579 11 : max_omega_div_omega_crit, omega_div_omega_crit, log_Teff, Lnuc_div_L, max_abs_vel, &
580 11 : species_mass_for_min_limit, species_mass_for_max_limit, center_gamma
581 :
582 : include 'formats'
583 :
584 : ierr = 0
585 11 : call get_star_ptr(id, s, ierr)
586 11 : if (ierr /= 0) then
587 0 : do_check_limits = terminate
588 0 : return
589 : end if
590 :
591 11 : if (s% RSP_flag) then
592 0 : max_period_number = s% RSP_max_num_periods
593 0 : period_number = s% RSP_num_periods
594 0 : GREKM_avg_abs = s% rsp_GREKM_avg_abs
595 : else
596 11 : max_period_number = 0
597 11 : period_number = -1
598 11 : GREKM_avg_abs = 1d99
599 : end if
600 :
601 11 : nz = s% nz
602 11 : do_check_limits = keep_going
603 :
604 11 : species_mass_for_min_limit = get_species_mass(s% star_species_mass_min_limit_iso)
605 11 : species_mass_for_max_limit = get_species_mass(s% star_species_mass_max_limit_iso)
606 :
607 11 : csound_surf = s% csound_face(1)
608 11 : if (s% u_flag) then
609 0 : v_surf = s% u(1)
610 0 : v_div_csound_max = maxval(abs(s% u(1:nz)/s% csound(1:nz)))
611 11 : else if (s% v_flag) then
612 0 : v_surf = s% v(1)
613 0 : v_div_csound_max = maxval(abs(s% v(1:nz)/s% csound_face(1:nz)))
614 : else
615 11 : v_surf = 0d0
616 11 : v_div_csound_max = 0d0
617 : end if
618 :
619 11 : remnant_mass = get_remnant_mass(s)/Msun
620 11 : ejecta_mass = get_ejecta_mass(s)/Msun
621 :
622 11 : if(s%u_flag) then
623 0 : max_abs_vel_loc = maxloc(abs(s%u(1:nz)),dim=1)
624 0 : max_abs_vel = s%u(max_abs_vel_loc)
625 11 : else if (s%v_flag) then
626 0 : max_abs_vel_loc = maxloc(abs(s%v(1:nz)),dim=1)
627 0 : max_abs_vel = s%v(max_abs_vel_loc)
628 : else
629 11 : max_abs_vel_loc = -1
630 11 : max_abs_vel = 0d0
631 : end if
632 :
633 :
634 11 : if (s% photosphere_r > 0d0) then
635 11 : v_surf_div_v_kh = abs(v_surf)*s% kh_timescale/s% photosphere_r
636 11 : if (s% cgrav(1) <= 0d0) then
637 0 : v_surf_div_v_esc = 0d0
638 : else
639 11 : v_surf_div_v_esc = v_surf/sqrt(2*s% cgrav(1)*s% m(1)/(s% photosphere_r*Rsun))
640 : end if
641 : else
642 0 : v_surf_div_v_kh = 0d0
643 0 : v_surf_div_v_esc = 0d0
644 : end if
645 :
646 11 : log_surface_gravity = safe_log10(s%grav(1))
647 11 : log_surface_temperature = s% lnT(1) / ln10
648 11 : log_surface_density = s% lnd(1)/ln10 ! log10(density at surface)
649 11 : log_surface_pressure = s% lnPeos(1)/ln10 ! log10(eos pressure at surface)
650 :
651 11 : center_gamma = center_value(s, s% gam)
652 :
653 11 : power_nuc_burn = s% power_nuc_burn
654 11 : power_h_burn = s% power_h_burn
655 11 : power_he_burn = s% power_he_burn
656 11 : power_z_burn = s% power_z_burn
657 11 : log_Teff = safe_log10(s% Teff)
658 11 : if (s% L_phot > 0d0) then
659 11 : Lnuc_div_L = s% L_nuc_burn_total / s% L_phot
660 11 : if (.not. s% get_delta_nu_from_scaled_solar) then
661 11 : delta_nu = 1d6/(2*s% photosphere_acoustic_r) ! microHz
662 : else
663 : delta_nu = &
664 : s% delta_nu_sun*sqrt(s% star_mass)*pow3(s% Teff/s% astero_Teff_sun) / &
665 0 : pow(s% L_phot,0.75d0)
666 : end if
667 : else
668 0 : Lnuc_div_L = 0d0
669 0 : delta_nu = 0d0
670 : end if
671 :
672 11 : if (s% dxdt_nuc_factor > 0d0) then
673 13109 : k = maxloc(s% eps_nuc(1:nz), dim=1)
674 11 : peak_burn_vconv_div_cs = s% conv_vel(k)/s% csound(k)
675 : else
676 0 : peak_burn_vconv_div_cs = 0d0
677 : end if
678 :
679 11 : if (s% initial_mass > s% he_core_mass) then
680 : envelope_fraction_left = &
681 11 : (s% star_mass - s% he_core_mass)/(s% initial_mass - s% he_core_mass)
682 : else
683 0 : envelope_fraction_left = 1
684 : end if
685 :
686 11 : max_logQ = -99
687 11 : min_logQ = 100
688 13098 : do k = 1, s% nz
689 13087 : logQ = s% lnd(k)/ln10 - 2*s% lnT(k)/ln10 + 12
690 : if (s% lnT(k)/ln10 < 5.5d0) then ! only worry about lower T cases
691 : if (logQ > max_logQ) max_logQ = logQ
692 : end if
693 11 : if (logQ < min_logQ) min_logQ = logQ
694 : end do
695 :
696 11 : min_pgas_div_p = 1d99
697 5617 : do k = s% nz, 1, -1
698 5617 : if (s% q(k) > s% Pgas_div_P_limit_max_q) exit
699 5617 : if (s% pgas(k)/s% Peos(k) < min_pgas_div_p) min_pgas_div_p = s% pgas(k)/s% Peos(k)
700 : end do
701 :
702 11 : max_omega_div_omega_crit = 0; k_omega = 0
703 11 : if (s% rotation_flag .and. s% omega_div_omega_crit_limit > 0) then
704 0 : do k = 1, s% nz
705 0 : omega_div_omega_crit = abs(s% omega(k))/omega_crit(s,k)
706 0 : if (omega_div_omega_crit > max_omega_div_omega_crit) then
707 0 : k_omega = k
708 0 : max_omega_div_omega_crit = omega_div_omega_crit
709 : end if
710 : end do
711 : end if
712 :
713 11 : if(s% max_abs_rel_run_E_err > 0d0) then
714 : if (abs(s% cumulative_energy_error/s% total_energy) > s% max_abs_rel_run_E_err &
715 0 : .and. .not. s% doing_relax) then
716 : write(*, '(/,a,/, 2e20.10)') &
717 0 : 'stop because abs rel_run_E_err exceeds limit (max_abs_rel_run_E_err) ', &
718 0 : abs(s% cumulative_energy_error/s% total_energy), s% max_abs_rel_run_E_err
719 0 : do_check_limits = terminate
720 0 : s% termination_code = t_max_abs_rel_run_E_err
721 0 : s% result_reason = abs_rel_run_E_err
722 0 : return
723 : end if
724 : end if
725 :
726 11 : if (max_abs_vel > clight) then
727 0 : if (s% retry_for_v_above_clight) then
728 : write(*, '(/,a,/, I5,1X,2e20.10)') &
729 0 : 'retry because maximum velocity exceeds speed of light ',max_abs_vel_loc,max_abs_vel,max_abs_vel/clight
730 0 : do_check_limits = retry
731 0 : return
732 : else
733 : write(*, '(/,a,/, I5,1X,2e20.10)') &
734 0 : 'Warning: maximum velocity exceeds speed of light ',max_abs_vel_loc,max_abs_vel,max_abs_vel/clight
735 : end if
736 : end if
737 :
738 11 : if (peak_burn_vconv_div_cs > 0.75d0*s% peak_burn_vconv_div_cs_limit) then
739 0 : write(*,1) 'peak_burn_vconv_div_cs: ', &
740 0 : peak_burn_vconv_div_cs / s% peak_burn_vconv_div_cs_limit, &
741 0 : peak_burn_vconv_div_cs, s% peak_burn_vconv_div_cs_limit
742 0 : k = maxloc(s% eps_nuc(1:nz), dim=1)
743 0 : write(*,2) 'maxloc eps_nuc', k, s% conv_vel(k), s% csound(k), s% eps_nuc(k)
744 0 : call mesa_error(__FILE__,__LINE__,'test do_one_utils')
745 : end if
746 :
747 11 : if (s% fe_core_infall < s% fe_core_infall_limit .and. &
748 : s% fe_core_infall > 0.99d0*s% fe_core_infall_limit) &
749 0 : write(*,1) 'nearing fe_core_infall limit', &
750 0 : s% fe_core_infall, s% fe_core_infall_limit
751 :
752 11 : if (s% non_fe_core_infall < s% non_fe_core_infall_limit .and. &
753 : s% non_fe_core_infall > 0.99d0*s% non_fe_core_infall_limit) &
754 0 : write(*,1) 'nearing non_fe_core_infall limit', &
755 0 : s% non_fe_core_infall, s% non_fe_core_infall_limit
756 :
757 11 : if (s% non_fe_core_rebound > 0.99d0*s% non_fe_core_rebound_limit) &
758 0 : write(*,1) 'nearing non_fe_core_rebound limit', &
759 0 : s% non_fe_core_rebound, s% non_fe_core_rebound_limit
760 :
761 : if (max_omega_div_omega_crit > 0.75d0*s% omega_div_omega_crit_limit .and. &
762 11 : s% omega_div_omega_crit_limit > 0 .and. k_omega > 0) &
763 0 : write(*,2) 'omega_div_omega_crit', k_omega, &
764 0 : max_omega_div_omega_crit, s% omega_div_omega_crit_limit, &
765 0 : s% m(k_omega)/Msun, s% r_equatorial(k_omega)/Rsun, &
766 0 : s% omega(k_omega), &
767 0 : sqrt(s% cgrav(k_omega)*s% m(k_omega)/ pow3(s% r_equatorial(k_omega)))
768 :
769 11 : if (s% star_age >= s% max_age .and. s% max_age > 0) then
770 : call compare_to_target('star_age >= max_age', s% star_age, s% max_age, &
771 1 : t_max_age)
772 :
773 10 : else if (s% time >= s% max_age_in_days*secday .and. s% max_age_in_days > 0) then
774 : call compare_to_target('time >= max_age_in_days', &
775 0 : s% time/secday, s% max_age_in_days, t_max_age)
776 :
777 10 : else if (s% time >= s% max_age_in_seconds .and. s% max_age_in_seconds > 0) then
778 : call compare_to_target('time >= max_age_in_seconds', &
779 0 : s% time, s% max_age_in_seconds, t_max_age)
780 :
781 10 : else if (max_omega_div_omega_crit >= s% omega_div_omega_crit_limit .and. &
782 : s% omega_div_omega_crit_limit > 0) then
783 : write(*, '(/,a,/, 2e20.10)') &
784 0 : 'stop max_omega_div_omega_crit >= omega_div_omega_crit_limit', &
785 0 : max_omega_div_omega_crit, s% omega_div_omega_crit_limit
786 0 : do_check_limits = terminate
787 0 : s% termination_code = t_max_omega_div_omega_crit
788 0 : s% result_reason = result_reason_normal
789 :
790 10 : else if (peak_burn_vconv_div_cs >= s% peak_burn_vconv_div_cs_limit) then
791 : write(*, '(/,a,/, 2e20.10)') &
792 0 : 'stop peak_burn_vconv_div_cs >= peak_burn_vconv_div_cs_limit', &
793 0 : peak_burn_vconv_div_cs, s% peak_burn_vconv_div_cs_limit
794 0 : do_check_limits = terminate
795 0 : s% termination_code = t_peak_burn_vconv_div_cs_limit
796 0 : s% result_reason = result_reason_normal
797 :
798 10 : else if (s% model_number >= s% max_model_number .and. s% max_model_number >= 0) then
799 0 : write(*, '(/,a,/, 2i9)') 'stop because model_number >= max_model_number', &
800 0 : s% model_number, s% max_model_number
801 0 : do_check_limits = terminate
802 0 : s% termination_code = t_max_model_number
803 0 : s% result_reason = result_reason_normal
804 :
805 10 : else if (period_number >= max_period_number .and. max_period_number >= 0) then
806 0 : write(*, '(/,a,/, 2i9)') 'stop because period_number >= max_period_number', &
807 0 : period_number, max_period_number
808 0 : do_check_limits = terminate
809 0 : s% termination_code = t_max_period_number
810 0 : s% result_reason = result_reason_normal
811 :
812 : else if (GREKM_avg_abs < s% RSP_GREKM_avg_abs_limit &
813 : .and. s% RSP_GREKM_avg_abs_limit >= 0 &
814 10 : .and. period_number >= 10) then
815 : write(*, '(/,a,/, 2e20.10)') &
816 0 : 'stop because GREKM_avg_abs < RSP_GREKM_avg_abs_limit', &
817 0 : GREKM_avg_abs, s% RSP_GREKM_avg_abs_limit
818 0 : do_check_limits = terminate
819 0 : s% termination_code = 0
820 0 : s% result_reason = result_reason_normal
821 :
822 10 : else if (s% center_degeneracy >= s% eta_center_limit) then
823 : call compare_to_target('center_degeneracy >= eta_center_limit', &
824 0 : s% center_degeneracy, s% eta_center_limit, t_eta_center_limit)
825 :
826 10 : else if (s% log_center_temperature >= s% log_center_temp_upper_limit) then
827 : call compare_to_target('log_center_temperature >= log_center_temp_upper_limit', &
828 0 : center_value(s, s% lnT)/ln10, s% log_center_temp_upper_limit, t_log_center_temp_upper_limit)
829 :
830 10 : else if (s% log_center_temperature <= s% log_center_temp_lower_limit) then
831 : call compare_to_target('log_center_temperature <= log_center_temp_lower_limit', &
832 : center_value(s, s% lnT)/ln10, s% log_center_temp_lower_limit, &
833 0 : t_log_center_temp_lower_limit)
834 :
835 10 : else if (s% max_entropy >= s% max_entropy_upper_limit) then
836 : call compare_to_target('max_entropy >= max_entropy_upper_limit', &
837 0 : s% max_entropy, s% max_entropy_upper_limit, t_max_entropy_upper_limit)
838 :
839 10 : else if (s% max_entropy <= s% max_entropy_lower_limit) then
840 : call compare_to_target('max_entropy <= max_entropy_lower_limit', &
841 : s% max_entropy, s% max_entropy_lower_limit, &
842 0 : t_max_entropy_lower_limit)
843 :
844 10 : else if (s% center_entropy >= s% center_entropy_upper_limit) then
845 : call compare_to_target('center_entropy >= center_entropy_upper_limit', &
846 0 : s% center_entropy, s% center_entropy_upper_limit, t_center_entropy_upper_limit)
847 :
848 10 : else if (s% center_entropy <= s% center_entropy_lower_limit) then
849 : call compare_to_target('center_entropy <= center_entropy_lower_limit', &
850 : s% center_entropy, s% center_entropy_lower_limit, &
851 0 : t_center_entropy_lower_limit)
852 :
853 10 : else if (s% log_center_density <= s% log_center_density_lower_limit) then
854 : call compare_to_target('log_center_density <= log_center_density_lower_limit', &
855 : center_value(s, s% lnd)/ln10, s% log_center_density_lower_limit, &
856 0 : t_log_center_density_lower_limit)
857 :
858 10 : else if (s% log_center_density >= s% log_center_density_upper_limit) then
859 : call compare_to_target('log_center_density >= log_center_density_upper_limit', &
860 0 : center_value(s, s% lnd)/ln10, s% log_center_density_upper_limit, t_log_center_density_upper_limit)
861 :
862 10 : else if (center_gamma > s% gamma_center_limit) then
863 : call compare_to_target('center_gamma > gamma_center_limit', &
864 0 : center_gamma, s% gamma_center_limit, t_gamma_center_limit)
865 :
866 10 : else if (s% log_max_temperature >= s% log_max_temp_upper_limit) then
867 : call compare_to_target('log_max_temperature >= log_max_temp_upper_limit', &
868 0 : s% log_max_temperature, s% log_max_temp_upper_limit, t_log_max_temp_upper_limit)
869 :
870 10 : else if (s% log_max_temperature <= s% log_max_temp_lower_limit) then
871 : call compare_to_target('log_max_temperature <= log_max_temp_lower_limit', &
872 0 : s% log_max_temperature, s% log_max_temp_lower_limit, t_log_max_temp_lower_limit)
873 :
874 10 : else if (s% center_he4 < s% HB_limit .and. s% center_h1 < 1d-4) then
875 0 : call compare_to_target('center he4 < HB_limit', s% center_he4, s% HB_limit, t_HB_limit)
876 :
877 10 : else if (s% star_mass_min_limit > 0 .and. s% star_mass <= s% star_mass_min_limit) then
878 : call compare_to_target('star_mass <= star_mass_min_limit', &
879 0 : s% star_mass, s% star_mass_min_limit, t_star_mass_min_limit)
880 :
881 10 : else if (s% star_mass_max_limit > 0 .and. s% star_mass >= s% star_mass_max_limit) then
882 : call compare_to_target('star_mass >= star_mass_max_limit', &
883 0 : s% star_mass, s% star_mass_max_limit, t_star_mass_max_limit)
884 :
885 10 : else if (s% remnant_mass_min_limit > 0 .and. remnant_mass <= s% remnant_mass_min_limit) then
886 : call compare_to_target('remnant_mass <= remnant_mass_min_limit', &
887 0 : remnant_mass, s% remnant_mass_min_limit, t_remnant_mass_min_limit)
888 :
889 10 : else if (s% ejecta_mass_max_limit > 0 .and. ejecta_mass >= s% ejecta_mass_max_limit) then
890 : call compare_to_target('ejecta_mass >= ejecta_mass_max_limit', &
891 0 : ejecta_mass, s% ejecta_mass_max_limit, t_ejecta_mass_max_limit)
892 :
893 10 : else if (species_mass_for_min_limit >= 0 .and. &
894 : species_mass_for_min_limit <= s% star_species_mass_min_limit) then
895 : call compare_to_target( &
896 : trim(s% star_species_mass_min_limit_iso) // ' total mass <= star_species_mass_min_limit', &
897 0 : species_mass_for_min_limit, s% star_species_mass_min_limit, t_star_species_mass_min_limit)
898 :
899 10 : else if (species_mass_for_max_limit >= s% star_species_mass_max_limit) then
900 : call compare_to_target( &
901 : trim(s% star_species_mass_max_limit_iso) // ' total mass >= star_species_mass_max_limit', &
902 0 : species_mass_for_max_limit, s% star_species_mass_max_limit, t_star_species_mass_max_limit)
903 :
904 10 : else if (s% xmstar_min_limit > 0 .and. s% xmstar <= s% xmstar_min_limit) then
905 : call compare_to_target('xmstar <= xmstar_min_limit', &
906 0 : s% xmstar, s% xmstar_min_limit, t_xmstar_min_limit)
907 :
908 10 : else if (s% xmstar_max_limit > 0 .and. s% xmstar >= s% xmstar_max_limit) then
909 : call compare_to_target('xmstar >= xmstar_max_limit', &
910 0 : s% xmstar, s% xmstar_max_limit, t_xmstar_max_limit)
911 :
912 10 : else if (s% star_mass - s% he_core_mass < s% envelope_mass_limit) then
913 : call compare_to_target('envelope mass < envelope_mass_limit', &
914 : s% star_mass - s% he_core_mass, s% envelope_mass_limit, &
915 0 : t_envelope_mass_limit)
916 :
917 10 : else if (envelope_fraction_left < s% envelope_fraction_left_limit) then
918 : call compare_to_target('envelope_fraction_left < limit', &
919 : envelope_fraction_left, s% envelope_fraction_left_limit, &
920 0 : t_envelope_fraction_left_limit)
921 :
922 10 : else if (s% he_core_mass >= s% he_core_mass_limit) then
923 : call compare_to_target('he_core_mass >= he_core_mass_limit', &
924 0 : s% he_core_mass, s% he_core_mass_limit, t_he_core_mass_limit)
925 :
926 10 : else if (s% co_core_mass >= s% co_core_mass_limit) then
927 : call compare_to_target('co_core_mass >= co_core_mass_limit', &
928 0 : s% co_core_mass, s% co_core_mass_limit, t_co_core_mass_limit)
929 :
930 10 : else if (s% one_core_mass >= s% one_core_mass_limit) then
931 : call compare_to_target('one_core_mass >= one_core_mass_limit', &
932 0 : s% one_core_mass, s% one_core_mass_limit, t_one_core_mass_limit)
933 :
934 10 : else if (s% fe_core_mass >= s% fe_core_mass_limit) then
935 : call compare_to_target('fe_core_mass >= fe_core_mass_limit', &
936 0 : s% fe_core_mass, s% fe_core_mass_limit, t_fe_core_mass_limit)
937 :
938 10 : else if (s% neutron_rich_core_mass >= s% neutron_rich_core_mass_limit) then
939 : call compare_to_target('neutron_rich_core_mass >= neutron_rich_core_mass_limit', &
940 0 : s% neutron_rich_core_mass, s% neutron_rich_core_mass_limit, t_neutron_rich_core_mass_limit)
941 :
942 : else if ( &
943 : s% he_core_mass >= s% co_core_mass .and. &
944 : s% co_core_mass > 0 .and. &
945 10 : s% center_he4 < 1d-4 .and. &
946 : s% he_core_mass - s% co_core_mass < s% he_layer_mass_lower_limit) then
947 : call compare_to_target('he layer mass < he_layer_mass_lower_limit', &
948 : s% he_core_mass - s% co_core_mass, s% he_layer_mass_lower_limit, &
949 0 : t_he_layer_mass_lower_limit)
950 :
951 : else if (abs(safe_log10(power_h_burn) - s% log_surface_luminosity) <= &
952 : s% abs_diff_lg_LH_lg_Ls_limit &
953 10 : .and. s% abs_diff_lg_LH_lg_Ls_limit > 0) then
954 : call compare_to_target('abs(lg_LH - lg_Ls) <= limit', &
955 : abs(safe_log10(power_h_burn) - s% log_surface_luminosity), &
956 0 : s% abs_diff_lg_LH_lg_Ls_limit, t_abs_diff_lg_LH_lg_Ls_limit)
957 :
958 10 : else if (s% Teff <= s% Teff_lower_limit) then
959 : call compare_to_target('Teff <= Teff_lower_limit', &
960 0 : s% Teff, s% Teff_lower_limit, t_Teff_lower_limit)
961 :
962 10 : else if (s% Teff >= s% Teff_upper_limit) then
963 : call compare_to_target('Teff >= Teff_upper_limit', &
964 0 : s% Teff, s% Teff_upper_limit, t_Teff_upper_limit)
965 :
966 10 : else if (delta_nu <= s% delta_nu_lower_limit .and. s% delta_nu_lower_limit > 0) then
967 : call compare_to_target('delta_nu <= delta_nu_lower_limit', &
968 0 : delta_nu, s% delta_nu_lower_limit, t_delta_nu_lower_limit)
969 :
970 10 : else if (delta_nu >= s% delta_nu_upper_limit .and. s% delta_nu_upper_limit > 0) then
971 : call compare_to_target('delta_nu >= delta_nu_upper_limit', &
972 0 : delta_nu, s% delta_nu_upper_limit, t_delta_nu_upper_limit)
973 :
974 10 : else if (s% delta_Pg <= s% delta_Pg_lower_limit .and. s% delta_Pg_lower_limit > 0) then
975 : call compare_to_target('delta_Pg <= delta_Pg_lower_limit', &
976 0 : s% delta_Pg, s% delta_Pg_lower_limit, t_delta_Pg_lower_limit)
977 :
978 10 : else if (s% delta_Pg >= s% delta_Pg_upper_limit .and. s% delta_Pg_upper_limit > 0) then
979 : call compare_to_target('delta_Pg >= delta_Pg_upper_limit', &
980 0 : s% delta_Pg, s% delta_Pg_upper_limit, t_delta_Pg_upper_limit)
981 :
982 10 : else if (s% photosphere_m - s% M_center/Msun <= s% photosphere_m_sub_M_center_limit) then
983 : call compare_to_target( &
984 : 'photosphere_m - M_center/Msun <= photosphere_m_sub_M_center_limit', &
985 : s% photosphere_m - s% M_center/Msun, &
986 : s% photosphere_m_sub_M_center_limit, &
987 0 : t_photosphere_m_sub_M_center_limit)
988 :
989 10 : else if (s% photosphere_m <= s% photosphere_m_lower_limit) then
990 : call compare_to_target('photosphere_m <= photosphere_m_lower_limit', &
991 0 : s% photosphere_m, s% photosphere_m_lower_limit, t_photosphere_m_lower_limit)
992 :
993 10 : else if (s% photosphere_m >= s% photosphere_m_upper_limit) then
994 : call compare_to_target('photosphere_m >= photosphere_m_upper_limit', &
995 0 : s% photosphere_m, s% photosphere_m_upper_limit, t_photosphere_m_upper_limit)
996 :
997 10 : else if (s% photosphere_r <= s% photosphere_r_lower_limit) then
998 : call compare_to_target('photosphere_r <= photosphere_r_lower_limit', &
999 0 : s% photosphere_r, s% photosphere_r_lower_limit, t_photosphere_r_lower_limit)
1000 :
1001 10 : else if (s% photosphere_r >= s% photosphere_r_upper_limit) then
1002 : call compare_to_target('photosphere_r >= photosphere_r_upper_limit', &
1003 0 : s% photosphere_r, s% photosphere_r_upper_limit, t_photosphere_r_upper_limit)
1004 :
1005 10 : else if (log_Teff <= s% log_Teff_lower_limit) then
1006 : call compare_to_target('log_Teff <= log_Teff_lower_limit', &
1007 0 : log_Teff, s% log_Teff_lower_limit, t_log_Teff_lower_limit)
1008 :
1009 10 : else if (log_Teff >= s% log_Teff_upper_limit) then
1010 : call compare_to_target('log_Teff >= log_Teff_upper_limit', &
1011 0 : log_Teff, s% log_Teff_upper_limit, t_log_Teff_upper_limit)
1012 :
1013 10 : else if (log_surface_temperature <= s% log_Tsurf_lower_limit) then
1014 : call compare_to_target('log_surface_temperature <= log_Tsurf_lower_limit', &
1015 0 : log_surface_temperature, s% log_Tsurf_lower_limit, t_log_Tsurf_lower_limit)
1016 :
1017 10 : else if (log_surface_temperature >= s% log_Tsurf_upper_limit) then
1018 : call compare_to_target('log_surface_temperature >= log_Tsurf_upper_limit', &
1019 0 : log_surface_temperature, s% log_Tsurf_upper_limit, t_log_Tsurf_upper_limit)
1020 :
1021 10 : else if (s% log_surface_radius <= s% log_Rsurf_lower_limit) then
1022 : call compare_to_target('log_surface_radius <= log_Rsurf_lower_limit', &
1023 0 : s% log_surface_radius, s% log_Rsurf_lower_limit, t_log_Rsurf_lower_limit)
1024 :
1025 10 : else if (s% log_surface_radius >= s% log_Rsurf_upper_limit) then
1026 : call compare_to_target('log_surface_radius >= log_Rsurf_upper_limit', &
1027 0 : s% log_surface_radius, s% log_Rsurf_upper_limit, t_log_Rsurf_upper_limit)
1028 :
1029 10 : else if (log_surface_pressure <= s% log_Psurf_lower_limit) then
1030 : call compare_to_target('log_surface_pressure <= log_Psurf_lower_limit', &
1031 0 : log_surface_pressure, s% log_Psurf_lower_limit, t_log_Psurf_lower_limit)
1032 :
1033 10 : else if (log_surface_pressure >= s% log_Psurf_upper_limit) then
1034 : call compare_to_target('log_surface_pressure >= log_Psurf_upper_limit', &
1035 0 : log_surface_pressure, s% log_Psurf_upper_limit, t_log_Psurf_upper_limit)
1036 :
1037 10 : else if (log_surface_density <= s% log_Dsurf_lower_limit) then
1038 : call compare_to_target('log_surface_density <= log_Dsurf_lower_limit', &
1039 0 : log_surface_density, s% log_Dsurf_lower_limit, t_log_Dsurf_lower_limit)
1040 :
1041 10 : else if (log_surface_density >= s% log_Dsurf_upper_limit) then
1042 : call compare_to_target('log_surface_density >= log_Dsurf_upper_limit', &
1043 0 : log_surface_density, s% log_Dsurf_upper_limit, t_log_Dsurf_upper_limit)
1044 :
1045 10 : else if (s% log_surface_luminosity <= s% log_L_lower_limit) then
1046 : call compare_to_target('log_surface_luminosity <= log_L_lower_limit', &
1047 0 : s% log_surface_luminosity, s% log_L_lower_limit, t_log_L_lower_limit)
1048 :
1049 10 : else if (s% log_surface_luminosity >= s% log_L_upper_limit) then
1050 : call compare_to_target('log_surface_luminosity >= log_L_upper_limit', &
1051 0 : s% log_surface_luminosity, s% log_L_upper_limit, t_log_L_upper_limit)
1052 :
1053 10 : else if (log_surface_gravity <= s% log_g_lower_limit) then
1054 : call compare_to_target('log_surface_gravity <= log_g_lower_limit', &
1055 0 : log_surface_gravity, s% log_g_lower_limit, t_log_g_lower_limit)
1056 :
1057 10 : else if (log_surface_gravity >= s% log_g_upper_limit) then
1058 : call compare_to_target('log_surface_gravity >= log_g_upper_limit', &
1059 0 : log_surface_gravity, s% log_g_upper_limit, t_log_g_upper_limit)
1060 :
1061 10 : else if (power_nuc_burn >= s% power_nuc_burn_upper_limit) then
1062 : call compare_to_target('power_nuc_burn >= power_nuc_burn_upper_limit', &
1063 0 : power_nuc_burn, s% power_nuc_burn_upper_limit, t_power_nuc_burn_upper_limit)
1064 :
1065 10 : else if (power_h_burn >= s% power_h_burn_upper_limit) then
1066 : call compare_to_target('power_h_burn >= power_h_burn_upper_limit', &
1067 0 : power_h_burn, s% power_h_burn_upper_limit, t_power_h_burn_upper_limit)
1068 :
1069 10 : else if (power_he_burn >= s% power_he_burn_upper_limit) then
1070 : call compare_to_target('power_he_burn >= power_he_burn_upper_limit', &
1071 0 : power_he_burn, s% power_he_burn_upper_limit, t_power_he_burn_upper_limit)
1072 :
1073 10 : else if (power_z_burn >= s% power_z_burn_upper_limit) then
1074 : call compare_to_target('power_z_burn >= power_z_burn_upper_limit', &
1075 0 : power_z_burn, s% power_z_burn_upper_limit, t_power_z_burn_upper_limit)
1076 :
1077 10 : else if (power_nuc_burn < s% power_nuc_burn_lower_limit) then
1078 : call compare_to_target('power_nuc_burn < power_nuc_burn_lower_limit', &
1079 0 : power_nuc_burn, s% power_nuc_burn_lower_limit, t_power_nuc_burn_lower_limit)
1080 :
1081 10 : else if (power_h_burn < s% power_h_burn_lower_limit) then
1082 : call compare_to_target('power_h_burn < power_h_burn_lower_limit', &
1083 0 : power_h_burn, s% power_h_burn_lower_limit, t_power_h_burn_lower_limit)
1084 :
1085 10 : else if (power_he_burn < s% power_he_burn_lower_limit) then
1086 : call compare_to_target('power_he_burn < power_he_burn_lower_limit', &
1087 0 : power_he_burn, s% power_he_burn_lower_limit, t_power_he_burn_lower_limit)
1088 :
1089 10 : else if (power_z_burn < s% power_z_burn_lower_limit) then
1090 : call compare_to_target('power_z_burn < power_z_burn_lower_limit', &
1091 0 : power_z_burn, s% power_z_burn_lower_limit, t_power_z_burn_lower_limit)
1092 :
1093 10 : else if (s% center_Ye < s% center_Ye_lower_limit) then
1094 : call compare_to_target('center_Ye < center_Ye_lower_limit', &
1095 0 : s% center_Ye, s% center_Ye_lower_limit, t_center_Ye_lower_limit)
1096 :
1097 10 : else if (s% R_center < s% center_R_lower_limit) then
1098 : call compare_to_target('R_center < center_R_lower_limit', &
1099 0 : s% R_center, s% center_R_lower_limit, t_center_R_lower_limit)
1100 :
1101 10 : else if (s% fe_core_infall > s% fe_core_infall_limit) then
1102 0 : if (abs(s% error_in_energy_conservation/s% total_energy_end) < &
1103 : s% hard_limit_for_rel_error_in_energy_conservation) then
1104 0 : do_check_limits = terminate
1105 0 : s% result_reason = result_reason_normal
1106 0 : s% termination_code = t_fe_core_infall_limit
1107 : write(*, '(/,a,/, 99e20.10)') &
1108 0 : 'stop because fe_core_infall > fe_core_infall_limit', &
1109 0 : s% fe_core_infall, s% fe_core_infall_limit
1110 : else
1111 0 : write(*,2) 'rel_E_err too large for fe_core_infall termination', &
1112 0 : s% model_number, s% error_in_energy_conservation/abs(s% total_energy_end)
1113 : end if
1114 :
1115 10 : else if (s% non_fe_core_infall > s% non_fe_core_infall_limit) then
1116 : call compare_to_target('non_fe_core_infall > non_fe_core_infall_limit', &
1117 0 : s% non_fe_core_infall, s% non_fe_core_infall_limit, t_non_fe_core_infall_limit)
1118 :
1119 10 : else if (s% non_fe_core_rebound > s% non_fe_core_rebound_limit) then
1120 : call compare_to_target('non_fe_core_rebound > non_fe_core_rebound_limit', &
1121 0 : s% non_fe_core_rebound, s% non_fe_core_rebound_limit, t_non_fe_core_rebound_limit)
1122 :
1123 10 : else if (v_surf/csound_surf > s% v_div_csound_surf_limit) then
1124 : call compare_to_target('v_surf/csound_surf > v_div_csound_surf_limit', &
1125 0 : v_surf/csound_surf, s% v_div_csound_surf_limit, t_v_div_csound_surf_limit)
1126 :
1127 10 : else if (v_div_csound_max > s% v_div_csound_max_limit) then
1128 : call compare_to_target('v_div_csound_max > v_div_csound_max_limit', &
1129 0 : v_div_csound_max, s% v_div_csound_max_limit, t_v_div_csound_max_limit)
1130 :
1131 10 : else if (s% min_gamma1 < s% gamma1_limit) then
1132 : call compare_to_target('min_gamma1 < gamma1_limit', &
1133 0 : s% min_gamma1, s% gamma1_limit, t_gamma1_limit)
1134 :
1135 10 : else if (min_pgas_div_p < s% Pgas_div_P_limit) then
1136 : call compare_to_target('min_pgas_div_p < Pgas_div_P_limit', &
1137 0 : min_pgas_div_p, s% Pgas_div_P_limit, t_Pgas_div_P_limit)
1138 :
1139 10 : else if (Lnuc_div_L <= s% Lnuc_div_L_lower_limit) then
1140 : call compare_to_target('Lnuc_div_L <= Lnuc_div_L_lower_limit', &
1141 0 : Lnuc_div_L, s% Lnuc_div_L_lower_limit, t_Lnuc_div_L_lower_limit)
1142 :
1143 10 : else if (Lnuc_div_L >= s% Lnuc_div_L_upper_limit) then
1144 : call compare_to_target('Lnuc_div_L >= Lnuc_div_L_upper_limit', &
1145 0 : Lnuc_div_L, s% Lnuc_div_L_upper_limit, t_Lnuc_div_L_upper_limit)
1146 :
1147 10 : else if (v_surf_div_v_kh <= s% v_surf_div_v_kh_lower_limit) then
1148 : call compare_to_target('v_surf_div_v_kh <= v_surf_div_v_kh_lower_limit', &
1149 0 : v_surf_div_v_kh, s% v_surf_div_v_kh_lower_limit, t_v_surf_div_v_kh_lower_limit)
1150 :
1151 10 : else if (v_surf_div_v_kh >= s% v_surf_div_v_kh_upper_limit) then
1152 : call compare_to_target('v_surf_div_v_kh >= v_surf_div_v_kh_upper_limit', &
1153 0 : v_surf_div_v_kh, s% v_surf_div_v_kh_upper_limit, t_v_surf_div_v_kh_upper_limit)
1154 :
1155 10 : else if (v_surf_div_v_esc >= s% v_surf_div_v_esc_limit) then
1156 : call compare_to_target('v_surf_div_v_esc >= v_surf_div_v_esc_limit', &
1157 0 : v_surf_div_v_esc, s% v_surf_div_v_esc_limit, t_v_surf_div_v_esc_limit)
1158 :
1159 10 : else if (v_surf*1d-5 >= s% v_surf_kms_limit) then
1160 : call compare_to_target('v_surf_kms >= v_surf_kms_limit', &
1161 0 : v_surf*1d-5, s% v_surf_kms_limit, t_v_surf_kms_limit)
1162 :
1163 : else if (s% cumulative_extra_heating >= &
1164 10 : s% stop_when_reach_this_cumulative_extra_heating .and.&
1165 : s% stop_when_reach_this_cumulative_extra_heating > 0) then
1166 : call compare_to_target( &
1167 : 'cumulative_extra_heating >= limit', &
1168 : s% cumulative_extra_heating, &
1169 : s% stop_when_reach_this_cumulative_extra_heating, &
1170 0 : t_cumulative_extra_heating_limit)
1171 :
1172 10 : else if (s% stop_near_zams .and. &
1173 : Lnuc_div_L >= s% Lnuc_div_L_zams_limit) then
1174 0 : do_check_limits = terminate
1175 0 : s% termination_code = t_Lnuc_div_L_zams_limit
1176 0 : s% result_reason = result_reason_normal
1177 : write(*, '(/,a,/, 99e20.10)') &
1178 0 : 'stop because Lnuc_div_L >= Lnuc_div_L_zams_limit', Lnuc_div_L, s% Lnuc_div_L_zams_limit
1179 :
1180 10 : else if (s% stop_at_phase_PreMS .and. s% phase_of_evolution == phase_PreMS) then
1181 0 : do_check_limits = terminate
1182 0 : s% termination_code = t_phase_PreMS
1183 0 : s% result_reason = result_reason_normal
1184 0 : write(*, '(/,a,/)') 'stop because phase_of_evolution == phase_PreMS'
1185 :
1186 10 : else if (s% stop_at_phase_ZAMS .and. s% phase_of_evolution == phase_ZAMS) then
1187 0 : do_check_limits = terminate
1188 0 : s% termination_code = t_phase_ZAMS
1189 0 : s% result_reason = result_reason_normal
1190 0 : write(*, '(/,a,/)') 'stop because phase_of_evolution == phase_ZAMS'
1191 :
1192 10 : else if (s% stop_at_phase_IAMS .and. s% phase_of_evolution == phase_IAMS) then
1193 0 : do_check_limits = terminate
1194 0 : s% termination_code = t_phase_IAMS
1195 0 : s% result_reason = result_reason_normal
1196 0 : write(*, '(/,a,/)') 'stop because phase_of_evolution == phase_IAMS'
1197 :
1198 10 : else if (s% stop_at_phase_TAMS .and. s% phase_of_evolution == phase_TAMS) then
1199 0 : do_check_limits = terminate
1200 0 : s% termination_code = t_phase_TAMS
1201 0 : s% result_reason = result_reason_normal
1202 0 : write(*, '(/,a,/)') 'stop because phase_of_evolution == phase_TAMS'
1203 :
1204 10 : else if (s% stop_at_phase_He_Burn .and. s% phase_of_evolution == phase_He_Burn) then
1205 0 : do_check_limits = terminate
1206 0 : s% termination_code = t_phase_He_Burn
1207 0 : s% result_reason = result_reason_normal
1208 0 : write(*, '(/,a,/)') 'stop because phase_of_evolution == phase_He_Burn'
1209 :
1210 10 : else if (s% stop_at_phase_ZACHeB .and. s% phase_of_evolution == phase_ZACHeB) then
1211 0 : do_check_limits = terminate
1212 0 : s% termination_code = t_phase_ZACHeB
1213 0 : s% result_reason = result_reason_normal
1214 0 : write(*, '(/,a,/)') 'stop because phase_of_evolution == phase_ZACHeB'
1215 :
1216 10 : else if (s% stop_at_phase_TACHeB .and. s% phase_of_evolution == phase_TACHeB) then
1217 0 : do_check_limits = terminate
1218 0 : s% termination_code = t_phase_TACHeB
1219 0 : s% result_reason = result_reason_normal
1220 0 : write(*, '(/,a,/)') 'stop because phase_of_evolution == phase_TACHeB'
1221 :
1222 10 : else if (s% stop_at_phase_TP_AGB .and. s% phase_of_evolution == phase_TP_AGB) then
1223 0 : do_check_limits = terminate
1224 0 : s% termination_code = t_phase_TP_AGB
1225 0 : s% result_reason = result_reason_normal
1226 0 : write(*, '(/,a,/)') 'stop because phase_of_evolution == phase_TP_AGB'
1227 :
1228 10 : else if (s% stop_at_phase_C_Burn .and. s% phase_of_evolution == phase_C_Burn) then
1229 0 : do_check_limits = terminate
1230 0 : s% termination_code = t_phase_C_Burn
1231 0 : s% result_reason = result_reason_normal
1232 0 : write(*, '(/,a,/)') 'stop because phase_of_evolution == phase_C_Burn'
1233 :
1234 10 : else if (s% stop_at_phase_Ne_Burn .and. s% phase_of_evolution == phase_Ne_Burn) then
1235 0 : do_check_limits = terminate
1236 0 : s% termination_code = t_phase_Ne_Burn
1237 0 : s% result_reason = result_reason_normal
1238 0 : write(*, '(/,a,/)') 'stop because phase_of_evolution == phase_Ne_Burn'
1239 :
1240 10 : else if (s% stop_at_phase_O_Burn .and. s% phase_of_evolution == phase_O_Burn) then
1241 0 : do_check_limits = terminate
1242 0 : s% termination_code = t_phase_O_Burn
1243 0 : s% result_reason = result_reason_normal
1244 0 : write(*, '(/,a,/)') 'stop because phase_of_evolution == phase_O_Burn'
1245 :
1246 10 : else if (s% stop_at_phase_Si_Burn .and. s% phase_of_evolution == phase_Si_Burn) then
1247 0 : do_check_limits = terminate
1248 0 : s% termination_code = t_phase_Si_Burn
1249 0 : s% result_reason = result_reason_normal
1250 0 : write(*, '(/,a,/)') 'stop because phase_of_evolution == phase_Si_Burn'
1251 :
1252 10 : else if (s% stop_at_phase_WDCS .and. s% phase_of_evolution == phase_WDCS) then
1253 0 : do_check_limits = terminate
1254 0 : s% termination_code = t_phase_WDCS
1255 0 : s% result_reason = result_reason_normal
1256 0 : write(*, '(/,a,/)') 'stop because phase_of_evolution == phase_WDCS'
1257 :
1258 : end if
1259 :
1260 11 : if (s% u_flag .or. s% v_flag) then ! Things that depend on hydro related quantities
1261 0 : if (s% shock_mass >= s% shock_mass_upper_limit .and. s% shock_mass_upper_limit > 0) then
1262 : call compare_to_target('shock_mass >= shock_mass_upper_limit', &
1263 0 : s% shock_mass, s% shock_mass_upper_limit, t_shock_mass_upper_limit)
1264 : end if
1265 : end if
1266 :
1267 11 : if (do_check_limits /= keep_going) return
1268 :
1269 100 : do j=1,num_xa_central_limits
1270 90 : if (s% xa_central_lower_limit(j) <= 0) cycle
1271 0 : if (len_trim(s% xa_central_lower_limit_species(j)) == 0) cycle
1272 0 : cid = chem_get_iso_id(s% xa_central_lower_limit_species(j))
1273 0 : if (cid <= 0) then
1274 0 : write(*,'(A)')
1275 0 : write(*,2) '<' // trim(s% xa_central_lower_limit_species(j)) // '>'
1276 0 : write(*,2) 'is invalid for xa_central_lower_limit_species', j
1277 0 : write(*,'(A)')
1278 0 : do_check_limits = terminate
1279 0 : return
1280 : end if
1281 0 : i = s% net_iso(cid)
1282 0 : if (i == 0) cycle
1283 0 : avg_x = center_avg_x(s,i)
1284 10 : if (avg_x < s% xa_central_lower_limit(j)) then
1285 : call compare_to_target('have dropped below central lower limit for ' // &
1286 : trim(s% xa_central_lower_limit_species(j)), &
1287 0 : avg_x, s% xa_central_lower_limit(j), t_xa_central_lower_limit)
1288 0 : exit
1289 : end if
1290 : end do
1291 :
1292 10 : if (do_check_limits /= keep_going) return
1293 :
1294 100 : do j=1,num_xa_central_limits
1295 90 : if (s% xa_central_upper_limit(j) <= 0) cycle
1296 0 : if (s% xa_central_upper_limit(j) >= 1) cycle
1297 0 : if (len_trim(s% xa_central_upper_limit_species(j)) == 0) cycle
1298 0 : cid = chem_get_iso_id(s% xa_central_upper_limit_species(j))
1299 0 : if (cid <= 0) then
1300 : !cycle
1301 0 : write(*,'(A)')
1302 0 : write(*,2) '<' // trim(s% xa_central_upper_limit_species(j)) // '>'
1303 0 : write(*,2) 'is invalid for xa_central_upper_limit_species', j
1304 0 : write(*,'(A)')
1305 0 : do_check_limits = terminate
1306 0 : return
1307 : end if
1308 0 : i = s% net_iso(cid)
1309 0 : if (i == 0) cycle
1310 0 : avg_x = center_avg_x(s,i)
1311 10 : if (avg_x > s% xa_central_upper_limit(j)) then
1312 : call compare_to_target('have risen above central upper limit for ' // &
1313 : trim(s% xa_central_upper_limit_species(j)), &
1314 0 : avg_x, s% xa_central_upper_limit(j), t_xa_central_upper_limit)
1315 0 : exit
1316 : end if
1317 : end do
1318 :
1319 10 : if (do_check_limits /= keep_going) return
1320 :
1321 100 : do j=1,num_xa_surface_limits
1322 90 : if (s% xa_surface_lower_limit(j) <= 0) cycle
1323 0 : if (len_trim(s% xa_surface_lower_limit_species(j)) == 0) cycle
1324 0 : cid = chem_get_iso_id(s% xa_surface_lower_limit_species(j))
1325 0 : if (cid <= 0) then
1326 0 : write(*,'(A)')
1327 0 : write(*,2) '<' // trim(s% xa_surface_lower_limit_species(j)) // '>'
1328 0 : write(*,2) 'is invalid for xa_surface_lower_limit_species', j
1329 0 : write(*,'(A)')
1330 0 : do_check_limits = terminate
1331 0 : return
1332 : end if
1333 0 : i = s% net_iso(cid)
1334 0 : if (i == 0) cycle
1335 0 : avg_x = surface_avg_x(s,i)
1336 10 : if (avg_x < s% xa_surface_lower_limit(j)) then
1337 : call compare_to_target('have dropped below surface lower limit for ' // &
1338 : trim(s% xa_surface_lower_limit_species(j)), &
1339 0 : avg_x, s% xa_surface_lower_limit(j), t_xa_surface_lower_limit)
1340 0 : exit
1341 : end if
1342 : end do
1343 :
1344 10 : if (do_check_limits /= keep_going) return
1345 :
1346 100 : do j=1,num_xa_surface_limits
1347 90 : if (s% xa_surface_upper_limit(j) <= 0) cycle
1348 0 : if (s% xa_surface_upper_limit(j) >= 1) cycle
1349 0 : if (len_trim(s% xa_surface_upper_limit_species(j)) == 0) cycle
1350 0 : cid = chem_get_iso_id(s% xa_surface_upper_limit_species(j))
1351 0 : if (cid <= 0) then
1352 : !cycle
1353 0 : write(*,'(A)')
1354 0 : write(*,2) '<' // trim(s% xa_surface_upper_limit_species(j)) // '>'
1355 0 : write(*,2) 'is invalid for xa_surface_upper_limit_species', j
1356 0 : write(*,'(A)')
1357 0 : do_check_limits = terminate
1358 0 : return
1359 : end if
1360 0 : i = s% net_iso(cid)
1361 0 : if (i == 0) cycle
1362 0 : avg_x = surface_avg_x(s,i)
1363 10 : if (avg_x > s% xa_surface_upper_limit(j)) then
1364 : call compare_to_target('have risen above surface upper limit for ' // &
1365 : trim(s% xa_surface_upper_limit_species(j)), &
1366 0 : avg_x, s% xa_surface_upper_limit(j), t_xa_surface_upper_limit)
1367 0 : exit
1368 : end if
1369 : end do
1370 :
1371 10 : if (do_check_limits /= keep_going) return
1372 :
1373 100 : do j=1,num_xa_average_limits
1374 90 : if (s% xa_average_lower_limit(j) <= 0) cycle
1375 0 : if (len_trim(s% xa_average_lower_limit_species(j)) == 0) cycle
1376 0 : cid = chem_get_iso_id(s% xa_average_lower_limit_species(j))
1377 0 : if (cid <= 0) then
1378 : !cycle
1379 0 : write(*,'(A)')
1380 0 : write(*,2) '<' // trim(s% xa_average_lower_limit_species(j)) // '>'
1381 0 : write(*,2) 'is invalid for xa_average_lower_limit_species', j
1382 0 : write(*,'(A)')
1383 0 : do_check_limits = terminate
1384 0 : return
1385 : end if
1386 0 : i = s% net_iso(cid)
1387 0 : if (i == 0) cycle
1388 0 : avg_x = dot_product(s% dq(1:nz), s% xa(i,1:nz))
1389 10 : if (avg_x < s% xa_average_lower_limit(j)) then
1390 : call compare_to_target('have dropped below average lower limit for ' // &
1391 : trim(s% xa_average_lower_limit_species(j)), &
1392 0 : avg_x, s% xa_average_lower_limit(j), t_xa_average_lower_limit)
1393 0 : exit
1394 : end if
1395 : end do
1396 :
1397 10 : if (do_check_limits /= keep_going) return
1398 :
1399 111 : do j=1,num_xa_average_limits
1400 90 : if (s% xa_average_upper_limit(j) <= 0) cycle
1401 0 : if (s% xa_average_upper_limit(j) >= 1) cycle
1402 0 : if (len_trim(s% xa_average_upper_limit_species(j)) == 0) cycle
1403 0 : cid = chem_get_iso_id(s% xa_average_upper_limit_species(j))
1404 0 : if (cid <= 0) then
1405 0 : write(*,'(A)')
1406 0 : write(*,2) '<' // trim(s% xa_average_upper_limit_species(j)) // '>'
1407 0 : write(*,2) 'is invalid for xa_average_upper_limit_species', j
1408 0 : write(*,'(A)')
1409 0 : do_check_limits = terminate
1410 0 : return
1411 : end if
1412 0 : i = s% net_iso(cid)
1413 0 : if (i == 0) cycle
1414 0 : avg_x = dot_product(s% dq(1:nz), s% xa(i,1:nz))
1415 10 : if (avg_x > s% xa_average_upper_limit(j)) then
1416 : call compare_to_target('have risen above average upper limit for ' // &
1417 : trim(s% xa_average_upper_limit_species(j)), &
1418 0 : avg_x, s% xa_average_upper_limit(j), t_xa_average_upper_limit)
1419 0 : exit
1420 : end if
1421 : end do
1422 :
1423 :
1424 : contains
1425 :
1426 :
1427 23 : subroutine compare_to_target(str, value, target_value, termination_code)
1428 : character (len=*), intent(in) :: str
1429 : real(dp), intent(in) :: value, target_value
1430 : integer, intent(in) :: termination_code
1431 1 : real(dp) :: err
1432 : include 'formats'
1433 : err = abs(value - target_value)/ &
1434 1 : (s% when_to_stop_atol + s% when_to_stop_rtol*max(abs(value),abs(target_value)))
1435 1 : if (err > 1) then
1436 0 : do_check_limits = redo
1437 0 : s% dt = 0.5d0*s% dt
1438 : write(*,'(/,a,5e20.10)') &
1439 0 : 'redo with smaller timestep to get closer to stopping target ' // trim(str), &
1440 0 : value, target_value
1441 : else
1442 1 : do_check_limits = terminate
1443 1 : s% result_reason = result_reason_normal
1444 1 : s% termination_code = termination_code
1445 1 : write(*, '(/,a,/, 99e20.10)') 'stop because ' // trim(str), value, target_value
1446 : end if
1447 11 : end subroutine compare_to_target
1448 :
1449 :
1450 22 : real(dp) function get_species_mass(str) ! Msun
1451 : use chem_lib, only: chem_get_iso_id
1452 : character(len=*), intent(in) :: str
1453 : integer :: id, j
1454 22 : get_species_mass = -1d0
1455 22 : id = chem_get_iso_id(str)
1456 22 : if (id > 0) then
1457 0 : j = s% net_iso(id)
1458 0 : if (j > 0) then
1459 0 : get_species_mass = dot_product(s% dm(1:s% nz),s% xa(j,1:s% nz))/Msun
1460 : end if
1461 : end if
1462 22 : end function get_species_mass
1463 :
1464 :
1465 : end function do_check_limits
1466 :
1467 :
1468 11 : integer function do_one_check_model(id)
1469 : use rates_def, only: i_rate
1470 : use chem_def, only: i_burn_c
1471 : use star_utils, only: update_time, total_times
1472 : integer, intent(in) :: id
1473 :
1474 : logical :: must_do_profile
1475 : real(dp), parameter :: log_he_temp = 7.8d0
1476 : real(dp), parameter :: d_tau_min = 1d-2, d_tau_max = 1d0
1477 : real(dp), parameter :: little_step_factor = 10d0, little_step_size = 10d0
1478 11 : real(dp) :: v, power_he_burn, power_z_burn, &
1479 11 : power_neutrinos
1480 : integer :: model, profile_priority, ierr
1481 : integer, parameter :: tau_ramp = 50
1482 : type (star_info), pointer :: s
1483 : logical :: logged
1484 : integer :: nz
1485 : logical, parameter :: dbg = .false.
1486 :
1487 : include 'formats'
1488 :
1489 11 : call get_star_ptr(id, s, ierr)
1490 11 : if (ierr /= 0) then
1491 11 : do_one_check_model = terminate
1492 : return
1493 : end if
1494 :
1495 11 : nz = s% nz
1496 11 : must_do_profile = .false.
1497 11 : profile_priority = delta_priority
1498 11 : model = s% model_number
1499 11 : do_one_check_model = keep_going
1500 :
1501 11 : do_one_check_model = do_check_limits(id)
1502 11 : if (do_one_check_model /= keep_going) then
1503 : if (dbg) write(*,*) 'do_check_limits /= keep_going'
1504 1 : write(*,'(A)')
1505 1 : must_do_profile = .true.
1506 : end if
1507 :
1508 11 : if (s% u_flag) then
1509 : v = s% u(1)
1510 : else if (s% v_flag) then
1511 : v = s% v(1)
1512 : else
1513 11 : v = 0d0
1514 : end if
1515 :
1516 11 : power_he_burn = s% power_he_burn
1517 11 : power_z_burn = s% power_z_burn
1518 11 : power_neutrinos = s% power_neutrinos
1519 :
1520 11 : if (must_do_profile) profile_priority = phase_priority
1521 :
1522 11 : logged = get_history_info(s, must_do_profile)
1523 :
1524 11 : if (logged .and. s% write_profiles_flag) then
1525 : if (s% model_number == s% profile_model &
1526 11 : .or. (s% profile_interval > 0 .and. &
1527 : (s% doing_first_model_of_run .or. &
1528 : mod(s% model_number,s% profile_interval) == 0))) then
1529 1 : if (s% write_profiles_flag) must_do_profile = .true.
1530 : if (s% model_number == s% profile_model .or. &
1531 1 : s% doing_first_model_of_run .or. &
1532 : (mod(s% model_number, s% priority_profile_interval) == 0)) then
1533 11 : profile_priority = phase_priority
1534 : end if
1535 : end if
1536 11 : if ( must_do_profile ) then
1537 : if (dbg) write(*,*) 'do_one_check_model: call set_save_profiles_info'
1538 2 : call set_save_profiles_info(s, profile_priority)
1539 : end if
1540 : end if
1541 :
1542 11 : end function do_one_check_model
1543 :
1544 : end module do_one_utils
|
