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 star_private_def
21 :
22 : use star_def
23 : use math_lib
24 : use utils_lib, only: is_bad_num, is_nan, mesa_error
25 :
26 : implicit none
27 :
28 : real(dp), parameter :: del_cntr_rho = 1d0
29 : real(dp), parameter :: min_cntr_rho = 3d0
30 : real(dp), parameter :: no_he_ignition_limit = 0.75d0
31 : real(dp), parameter :: no_cntr_T_drops_limit = 6.5d0
32 :
33 : real(dp), parameter :: center_h_gone = 1d-3
34 : real(dp), parameter :: center_h_going = 1d0/3d0
35 : real(dp), parameter :: center_he_going = 5d-2
36 :
37 : contains
38 :
39 2 : subroutine star_private_def_init
40 : use num_def
41 : integer :: i, im1
42 : logical :: okay
43 :
44 : include 'formats'
45 :
46 1 : okay = .true.
47 :
48 34 : auto_diff_star_d1_names(1:auto_diff_star_num_vars) = ''
49 1 : auto_diff_star_d1_names(i_lnd_m1) = 'i_lnd_m1'
50 1 : auto_diff_star_d1_names(i_lnd_00) = 'i_lnd_00'
51 1 : auto_diff_star_d1_names(i_lnd_p1) = 'i_lnd_p1'
52 1 : auto_diff_star_d1_names(i_lnT_m1) = 'i_lnT_m1'
53 1 : auto_diff_star_d1_names(i_lnT_00) = 'i_lnT_00'
54 1 : auto_diff_star_d1_names(i_lnT_p1) = 'i_lnT_p1'
55 1 : auto_diff_star_d1_names(i_w_m1) = 'i_w_m1'
56 1 : auto_diff_star_d1_names(i_w_00) = 'i_w_00'
57 1 : auto_diff_star_d1_names(i_w_p1) = 'i_w_p1'
58 1 : auto_diff_star_d1_names(i_lnR_m1) = 'i_lnR_m1'
59 1 : auto_diff_star_d1_names(i_lnR_00) = 'i_lnR_00'
60 1 : auto_diff_star_d1_names(i_lnR_p1) = 'i_lnR_p1'
61 1 : auto_diff_star_d1_names(i_v_m1) = 'i_v_m1'
62 1 : auto_diff_star_d1_names(i_v_00) = 'i_v_00'
63 1 : auto_diff_star_d1_names(i_v_p1) = 'i_v_p1'
64 1 : auto_diff_star_d1_names(i_L_m1) = 'i_L_m1'
65 1 : auto_diff_star_d1_names(i_L_00) = 'i_L_00'
66 1 : auto_diff_star_d1_names(i_L_p1) = 'i_L_p1'
67 1 : auto_diff_star_d1_names(i_Hp_m1) = 'i_Hp_m1'
68 1 : auto_diff_star_d1_names(i_Hp_00) = 'i_Hp_00'
69 1 : auto_diff_star_d1_names(i_Hp_p1) = 'i_Hp_p1'
70 1 : auto_diff_star_d1_names(i_w_div_wc_m1) = 'i_w_div_wc_m1'
71 1 : auto_diff_star_d1_names(i_w_div_wc_00) = 'i_w_div_wc_00'
72 1 : auto_diff_star_d1_names(i_w_div_wc_p1) = 'i_w_div_wc_p1'
73 1 : auto_diff_star_d1_names(i_jrot_m1) = 'i_jrot_m1'
74 1 : auto_diff_star_d1_names(i_jrot_00) = 'i_jrot_00'
75 1 : auto_diff_star_d1_names(i_jrot_p1) = 'i_jrot_p1'
76 1 : auto_diff_star_d1_names(i_xtra1_m1) = 'i_xtra1_m1'
77 1 : auto_diff_star_d1_names(i_xtra1_00) = 'i_xtra1_00'
78 1 : auto_diff_star_d1_names(i_xtra1_p1) = 'i_xtra1_p1'
79 1 : auto_diff_star_d1_names(i_xtra2_m1) = 'i_xtra2_m1'
80 1 : auto_diff_star_d1_names(i_xtra2_00) = 'i_xtra2_00'
81 1 : auto_diff_star_d1_names(i_xtra2_p1) = 'i_xtra2_p1'
82 :
83 134 : termination_code_str(1:num_termination_codes) = ''
84 :
85 1 : termination_code_str(t_max_age) = 'max_age'
86 1 : termination_code_str(t_max_omega_div_omega_crit) = 'max_omega_div_omega_crit'
87 1 : termination_code_str(t_peak_burn_vconv_div_cs_limit) = 'peak_burn_vconv_div_cs_limit'
88 1 : termination_code_str(t_max_model_number) = 'max_model_number'
89 1 : termination_code_str(t_eta_center_limit) = 'eta_center_limit'
90 1 : termination_code_str(t_log_center_temp_upper_limit) = 'log_center_temp_upper_limit'
91 1 : termination_code_str(t_log_center_temp_lower_limit) = 'log_center_temp_lower_limit'
92 1 : termination_code_str(t_center_entropy_upper_limit) = 'center_entropy_upper_limit'
93 1 : termination_code_str(t_center_entropy_lower_limit) = 'center_entropy_lower_limit'
94 1 : termination_code_str(t_max_entropy_upper_limit) = 'max_entropy_upper_limit'
95 1 : termination_code_str(t_max_entropy_lower_limit) = 'max_entropy_lower_limit'
96 1 : termination_code_str(t_log_center_density_upper_limit) = 'log_center_density_upper_limit'
97 1 : termination_code_str(t_log_center_density_lower_limit) = 'log_center_density_lower_limit'
98 1 : termination_code_str(t_gamma_center_limit) = 'gamma_center_limit'
99 1 : termination_code_str(t_log_max_temp_upper_limit) = 'log_max_temp_upper_limit'
100 1 : termination_code_str(t_log_max_temp_lower_limit) = 'log_max_temp_lower_limit'
101 1 : termination_code_str(t_HB_limit) = 'HB_limit'
102 1 : termination_code_str(t_star_mass_min_limit) = 'star_mass_min_limit'
103 1 : termination_code_str(t_star_mass_max_limit) = 'star_mass_max_limit'
104 1 : termination_code_str(t_remnant_mass_min_limit) = 'remnant_mass_min_limit'
105 1 : termination_code_str(t_ejecta_mass_max_limit) = 'ejecta_mass_max_limit'
106 :
107 1 : termination_code_str(t_star_species_mass_min_limit) = 'star_species_mass_min_limit'
108 1 : termination_code_str(t_star_species_mass_max_limit) = 'star_species_mass_max_limit'
109 :
110 1 : termination_code_str(t_xmstar_min_limit) = 'xmstar_min_limit'
111 1 : termination_code_str(t_xmstar_max_limit) = 'xmstar_max_limit'
112 1 : termination_code_str(t_envelope_mass_limit) = 'envelope_mass_limit'
113 1 : termination_code_str(t_envelope_fraction_left_limit) = 'envelope_fraction_left_limit'
114 :
115 1 : termination_code_str(t_he_core_mass_limit) = 'he_core_mass_limit'
116 1 : termination_code_str(t_co_core_mass_limit) = 'co_core_mass_limit'
117 1 : termination_code_str(t_one_core_mass_limit) = 'one_core_mass_limit'
118 1 : termination_code_str(t_fe_core_mass_limit) = 'fe_core_mass_limit'
119 1 : termination_code_str(t_neutron_rich_core_mass_limit) = 'neutron_rich_core_mass_limit'
120 :
121 1 : termination_code_str(t_he_layer_mass_lower_limit) = 'he_layer_mass_lower_limit'
122 1 : termination_code_str(t_abs_diff_lg_LH_lg_Ls_limit) = 'abs_diff_lg_LH_lg_Ls_limit'
123 1 : termination_code_str(t_Teff_lower_limit) = 'Teff_lower_limit'
124 1 : termination_code_str(t_Teff_upper_limit) = 'Teff_upper_limit'
125 1 : termination_code_str(t_delta_nu_lower_limit) = 'delta_nu_lower_limit'
126 1 : termination_code_str(t_delta_nu_upper_limit) = 'delta_nu_upper_limit'
127 1 : termination_code_str(t_delta_Pg_lower_limit) = 'delta_Pg_lower_limit'
128 1 : termination_code_str(t_delta_Pg_upper_limit) = 'delta_Pg_upper_limit'
129 1 : termination_code_str(t_shock_mass_upper_limit) = 'shock_mass_upper_limit'
130 1 : termination_code_str(t_photosphere_m_sub_M_center_limit) = 'photosphere_m_sub_M_center_limit'
131 1 : termination_code_str(t_photosphere_m_lower_limit) = 'photosphere_m_lower_limit'
132 1 : termination_code_str(t_photosphere_m_upper_limit) = 'photosphere_m_upper_limit'
133 1 : termination_code_str(t_photosphere_r_lower_limit) = 'photosphere_r_lower_limit'
134 1 : termination_code_str(t_photosphere_r_upper_limit) = 'photosphere_r_upper_limit'
135 1 : termination_code_str(t_log_Teff_lower_limit) = 'log_Teff_lower_limit'
136 1 : termination_code_str(t_log_Teff_upper_limit) = 'log_Teff_upper_limit'
137 1 : termination_code_str(t_log_Tsurf_lower_limit) = 'log_Tsurf_lower_limit'
138 1 : termination_code_str(t_log_Tsurf_upper_limit) = 'log_Tsurf_upper_limit'
139 1 : termination_code_str(t_log_Rsurf_lower_limit) = 'log_Rsurf_lower_limit'
140 1 : termination_code_str(t_log_Rsurf_upper_limit) = 'log_Rsurf_upper_limit'
141 1 : termination_code_str(t_log_Psurf_lower_limit) = 'log_Psurf_lower_limit'
142 1 : termination_code_str(t_log_Psurf_upper_limit) = 'log_Psurf_upper_limit'
143 1 : termination_code_str(t_log_Dsurf_lower_limit) = 'log_Dsurf_lower_limit'
144 1 : termination_code_str(t_log_Dsurf_upper_limit) = 'log_Dsurf_upper_limit'
145 1 : termination_code_str(t_log_L_lower_limit) = 'log_L_lower_limit'
146 1 : termination_code_str(t_log_L_upper_limit) = 'log_L_upper_limit'
147 1 : termination_code_str(t_log_g_lower_limit) = 'log_g_lower_limit'
148 1 : termination_code_str(t_log_g_upper_limit) = 'log_g_upper_limit'
149 1 : termination_code_str(t_power_nuc_burn_upper_limit) = 'power_nuc_burn_upper_limit'
150 1 : termination_code_str(t_power_h_burn_upper_limit) = 'power_h_burn_upper_limit'
151 1 : termination_code_str(t_power_he_burn_upper_limit) = 'power_he_burn_upper_limit'
152 1 : termination_code_str(t_power_z_burn_upper_limit) = 'power_z_burn_upper_limit'
153 1 : termination_code_str(t_power_nuc_burn_lower_limit) = 'power_nuc_burn_lower_limit'
154 1 : termination_code_str(t_power_h_burn_lower_limit) = 'power_h_burn_lower_limit'
155 1 : termination_code_str(t_power_he_burn_lower_limit) = 'power_he_burn_lower_limit'
156 1 : termination_code_str(t_power_z_burn_lower_limit) = 'power_z_burn_lower_limit'
157 1 : termination_code_str(t_center_R_lower_limit) = 'center_R_lower_limit'
158 1 : termination_code_str(t_center_Ye_lower_limit) = 'center_Ye_lower_limit'
159 1 : termination_code_str(t_fe_core_infall_limit) = 'fe_core_infall_limit'
160 1 : termination_code_str(t_non_fe_core_infall_limit) = 'non_fe_core_infall_limit'
161 1 : termination_code_str(t_non_fe_core_rebound_limit) = 'non_fe_core_rebound_limit'
162 1 : termination_code_str(t_v_div_csound_max_limit) = 'v_div_csound_max_limit'
163 1 : termination_code_str(t_v_div_csound_surf_limit) = 'v_div_csound_surf_limit'
164 1 : termination_code_str(t_gamma1_limit) = 'gamma1_limit'
165 1 : termination_code_str(t_Pgas_div_P_limit) = 'Pgas_div_P_limit'
166 1 : termination_code_str(t_Lnuc_div_L_lower_limit) = 'Lnuc_div_L_lower_limit'
167 1 : termination_code_str(t_Lnuc_div_L_upper_limit) = 'Lnuc_div_L_upper_limit'
168 1 : termination_code_str(t_v_surf_div_v_kh_lower_limit) = 'v_surf_div_v_kh_lower_limit'
169 1 : termination_code_str(t_v_surf_div_v_kh_upper_limit) = 'v_surf_div_v_kh_upper_limit'
170 1 : termination_code_str(t_v_surf_div_v_esc_limit) = 'v_surf_div_v_esc_limit'
171 1 : termination_code_str(t_v_surf_kms_limit) = 'v_surf_kms_limit'
172 1 : termination_code_str(t_Lnuc_div_L_zams_limit) = 'Lnuc_div_L_zams_limit'
173 1 : termination_code_str(t_phase_PreMS) = 'phase_PreMS'
174 1 : termination_code_str(t_phase_ZAMS) = 'phase_ZAMS'
175 1 : termination_code_str(t_phase_IAMS) = 'phase_IAMS'
176 1 : termination_code_str(t_phase_TAMS) = 'phase_TAMS'
177 1 : termination_code_str(t_phase_He_Burn) = 'phase_He_Burn'
178 1 : termination_code_str(t_phase_ZACHeB) = 'phase_ZACHeB'
179 1 : termination_code_str(t_phase_TACHeB) = 'phase_TACHeB'
180 1 : termination_code_str(t_phase_TP_AGB) = 'phase_TP_AGB'
181 1 : termination_code_str(t_phase_C_Burn) = 'phase_C_Burn'
182 1 : termination_code_str(t_phase_Ne_Burn) = 'phase_Ne_Burn'
183 1 : termination_code_str(t_phase_O_Burn) = 'phase_O_Burn'
184 1 : termination_code_str(t_phase_Si_Burn) = 'phase_Si_Burn'
185 1 : termination_code_str(t_phase_WDCS) = 'phase_WDCS'
186 1 : termination_code_str(t_xa_central_lower_limit) = 'xa_central_lower_limit'
187 1 : termination_code_str(t_xa_central_upper_limit) = 'xa_central_upper_limit'
188 1 : termination_code_str(t_xa_surface_lower_limit) = 'xa_surface_lower_limit'
189 1 : termination_code_str(t_xa_surface_upper_limit) = 'xa_surface_upper_limit'
190 1 : termination_code_str(t_xa_average_lower_limit) = 'xa_average_lower_limit'
191 1 : termination_code_str(t_xa_average_upper_limit) = 'xa_average_upper_limit'
192 1 : termination_code_str(t_surface_accel_div_grav_limit) = 'surface_accel_div_grav_limit'
193 1 : termination_code_str(t_adjust_mesh_failed) = 'adjust_mesh_failed'
194 1 : termination_code_str(t_dt_is_zero) = 'dt_is_zero'
195 1 : termination_code_str(t_min_timestep_limit) = 'min_timestep_limit'
196 1 : termination_code_str(t_failed_prepare_for_new_try) = 'failed_prepare_for_new_try'
197 1 : termination_code_str(t_negative_total_angular_momentum) = 'negative_total_angular_momentum'
198 1 : termination_code_str(t_max_number_retries) = 'max_number_retries'
199 1 : termination_code_str(t_redo_limit) = 'redo_limit'
200 1 : termination_code_str(t_solve_burn) = 'solve_burn'
201 1 : termination_code_str(t_solve_hydro) = 'solve_hydro'
202 1 : termination_code_str(t_solve_mix) = 'solve_mix'
203 1 : termination_code_str(t_solve_omega_mix) = 'solve_omega_mix'
204 1 : termination_code_str(t_timestep_controller) = 'timestep_controller'
205 1 : termination_code_str(t_relax_finished_okay) = 'relax_finished_okay'
206 1 : termination_code_str(t_delta_total_energy) = 'delta total energy'
207 1 : termination_code_str(t_cumulative_extra_heating_limit) = 'cumulative extra heating'
208 1 : termination_code_str(t_max_explicit_hydro_nsteps) = 'reached max explicit hydro nsteps'
209 1 : termination_code_str(t_max_period_number) = 'reached max number of periods'
210 1 : termination_code_str(t_max_abs_rel_run_E_err) = 'exceeded max abs rel_run_E_err'
211 :
212 1 : termination_code_str(t_extras_check_model) = 'extras_check_model'
213 1 : termination_code_str(t_extras_finish_step) = 'extras_finish_step'
214 :
215 1 : termination_code_str(t_xtra1) = 'customize by setting termination_code_str(t_xtra1)'
216 1 : termination_code_str(t_xtra2) = 'customize by setting termination_code_str(t_xtra2)'
217 1 : termination_code_str(t_xtra3) = 'customize by setting termination_code_str(t_xtra3)'
218 1 : termination_code_str(t_xtra4) = 'customize by setting termination_code_str(t_xtra4)'
219 1 : termination_code_str(t_xtra5) = 'customize by setting termination_code_str(t_xtra5)'
220 1 : termination_code_str(t_xtra6) = 'customize by setting termination_code_str(t_xtra6)'
221 1 : termination_code_str(t_xtra7) = 'customize by setting termination_code_str(t_xtra7)'
222 1 : termination_code_str(t_xtra8) = 'customize by setting termination_code_str(t_xtra8)'
223 1 : termination_code_str(t_xtra9) = 'customize by setting termination_code_str(t_xtra9)'
224 :
225 134 : do i=1,num_termination_codes
226 134 : if (len_trim(termination_code_str(i)) == 0) then
227 0 : if (i > 1) then
228 0 : im1 = i-1
229 : write(*,2) 'missing termination_code_str following ' // &
230 0 : trim(termination_code_str(im1)), i
231 : else
232 0 : write(*,2) 'missing termination_code_str 1'
233 : end if
234 : okay = .false.
235 : end if
236 : end do
237 :
238 59 : dt_why_str(1:numTlim) = ''
239 :
240 1 : dt_why_str(Tlim_struc) = 'varcontrol'
241 1 : dt_why_str(Tlim_max_timestep_factor) = 'max increase'
242 1 : dt_why_str(Tlim_min_timestep_factor) = 'max decrease'
243 1 : dt_why_str(Tlim_solver) = 'solver'
244 1 : dt_why_str(Tlim_num_burn_steps) = 'burn steps'
245 1 : dt_why_str(Tlim_num_diff_solver_steps) = 'diff steps'
246 1 : dt_why_str(Tlim_num_diff_solver_iters) = 'diff iters'
247 1 : dt_why_str(Tlim_dX) = 'dX'
248 1 : dt_why_str(Tlim_dX_div_X) = 'dX/X'
249 1 : dt_why_str(Tlim_dL_div_L) = 'dL/L'
250 1 : dt_why_str(Tlim_dlgP) = 'lgP'
251 1 : dt_why_str(Tlim_dlgRho) = 'lgRho'
252 1 : dt_why_str(Tlim_dlgT) = 'lgT'
253 1 : dt_why_str(Tlim_dlgE) = 'lgE'
254 1 : dt_why_str(Tlim_dlgR) = 'lgR'
255 1 : dt_why_str(Tlim_dlgL_nuc_cat) = 'Lnuc_cat'
256 1 : dt_why_str(Tlim_dlgL_H) = 'Lnuc_H'
257 1 : dt_why_str(Tlim_dlgL_He) = 'Lnuc_He'
258 1 : dt_why_str(Tlim_dlgL_z) = 'Lnuc_z'
259 1 : dt_why_str(Tlim_dlgL_nuc) = 'Lnuc'
260 1 : dt_why_str(Tlim_dlgTeff) = 'lgTeff'
261 1 : dt_why_str(Tlim_dlgRho_cntr) = 'lgRho_cntr'
262 1 : dt_why_str(Tlim_dlgT_max) = 'lgT_max'
263 1 : dt_why_str(Tlim_dlgT_max_at_high_T) = 'lgT_max_hi_T'
264 1 : dt_why_str(Tlim_dlgT_cntr) = 'lgT_cntr'
265 1 : dt_why_str(Tlim_dlgP_cntr) = 'lgP_cntr'
266 1 : dt_why_str(Tlim_dlog_eps_nuc) = 'log_eps_nuc'
267 1 : dt_why_str(Tlim_lg_XH_cntr) = 'lg_XH_cntr'
268 1 : dt_why_str(Tlim_dmstar) = 'delta_mstar'
269 1 : dt_why_str(Tlim_dt_div_dt_cell_collapse) = 'dt_collapse'
270 1 : dt_why_str(Tlim_dt_div_min_dr_div_cs) = 'min_dr_div_cs'
271 1 : dt_why_str(Tlim_lgL) = 'lgL'
272 1 : dt_why_str(Tlim_force_timestep) = 'force_dt'
273 1 : dt_why_str(Tlim_max_timestep) = 'max_dt'
274 1 : dt_why_str(Tlim_timestep_hold) = 'hold'
275 1 : dt_why_str(Tlim_dX_nuc_drop) = 'dX_nuc'
276 1 : dt_why_str(Tlim_dX_div_X_cntr) = 'dX_div_X_cntr'
277 1 : dt_why_str(Tlim_lg_XHe_cntr) = 'lg_XHe_cntr'
278 1 : dt_why_str(Tlim_lg_XC_cntr) = 'lg_XC_cntr'
279 : dt_why_str(Tlim_lg_XNe_cntr) = 'lg_XNe_cntr'
280 : dt_why_str(Tlim_lg_XO_cntr) = 'lg_XO_cntr'
281 1 : dt_why_str(Tlim_lg_XSi_cntr) = 'lg_XSi_cntr'
282 1 : dt_why_str(Tlim_XH_cntr) = 'XH_cntr'
283 1 : dt_why_str(Tlim_XHe_cntr) = 'XHe_cntr'
284 1 : dt_why_str(Tlim_XC_cntr) = 'XC_cntr'
285 1 : dt_why_str(Tlim_XNe_cntr) = 'XNe_cntr'
286 1 : dt_why_str(Tlim_XO_cntr) = 'XO_cntr'
287 1 : dt_why_str(Tlim_XSi_cntr) = 'XSi_cntr'
288 1 : dt_why_str(Tlim_neg_X) = 'neg_mass_frac'
289 1 : dt_why_str(Tlim_bad_Xsum) = 'bad_X_sum'
290 1 : dt_why_str(Tlim_dlgL_power_photo) = 'lgL_power_phot'
291 1 : dt_why_str(Tlim_delta_HR) = 'delta_HR'
292 1 : dt_why_str(Tlim_del_mdot) = 'delta mdot'
293 1 : dt_why_str(Tlim_adjust_J_q) = 'adjust_J_q'
294 1 : dt_why_str(Tlim_delta_Ye_highT) = 'highT del Ye'
295 1 : dt_why_str(Tlim_error_in_energy_conservation) = 'rel_E_err'
296 1 : dt_why_str(Tlim_retry) = 'retry'
297 1 : dt_why_str(Tlim_binary) = 'binary'
298 1 : dt_why_str(Tlim_error_other) = 'error_other'
299 1 : dt_why_str(Tlim_other_timestep_limit) = 'other_limit'
300 :
301 59 : do i=1,numTlim
302 59 : if (len_trim(dt_why_str(i)) == 0) then
303 0 : if (i > 1) then
304 0 : im1 = i-1
305 0 : write(*,2) 'missing dt_why_str following ' // trim(dt_why_str(im1)), i
306 : else
307 0 : write(*,2) 'missing dt_why_str 1'
308 : end if
309 : okay = .false.
310 : end if
311 : end do
312 :
313 1 : if (.not. okay) call mesa_error(__FILE__,__LINE__,'star_private_def_init')
314 :
315 1 : end subroutine star_private_def_init
316 :
317 :
318 2 : subroutine alloc_star(id, ierr)
319 : integer, intent(out) :: id, ierr
320 : integer :: i
321 : type (star_info), pointer :: s
322 :
323 2 : ierr = 0
324 2 : id = -1
325 4 : !$omp critical (star_handle)
326 2 : call init_star_handles()
327 2 : do i = 1, max_star_handles
328 2 : if (.not. star_handles(i)% in_use) then
329 2 : star_handles(i)% in_use = .true.
330 2 : id = i
331 2 : exit
332 : end if
333 : end do
334 : !$omp end critical (star_handle)
335 2 : if (id == -1) then
336 0 : ierr = -1
337 0 : return
338 : end if
339 2 : if (star_handles(id)% id /= id) then
340 0 : ierr = -1
341 0 : return
342 : end if
343 2 : s => star_handles(id)
344 :
345 : end subroutine alloc_star
346 :
347 :
348 2 : subroutine init_star_handles()
349 : integer :: i
350 :
351 2 : if (.not. have_initialized_star_handles) then
352 22 : do i = 1, max_star_handles
353 20 : star_handles(i)% id = i
354 22 : star_handles(i)% in_use = .false.
355 : end do
356 2 : have_initialized_star_handles = .true.
357 : end if
358 :
359 2 : end subroutine init_star_handles
360 :
361 :
362 0 : integer function find_next_star_id()
363 : integer :: id
364 0 : id = 0
365 0 : !$omp critical (star_handle_next)
366 0 : if (have_initialized_star_handles) then
367 0 : do id = 1, max_star_handles
368 0 : if (star_handles(id)% in_use .eqv. .false.) exit
369 : end do
370 : end if
371 : !$omp end critical (star_handle_next)
372 0 : find_next_star_id = id
373 0 : end function find_next_star_id
374 :
375 :
376 12 : integer function how_many_allocated_star_ids()
377 : integer :: id
378 12 : how_many_allocated_star_ids = 0
379 12 : if (have_initialized_star_handles) then
380 132 : do id = 1, max_star_handles
381 120 : if (star_handles(id)% in_use .eqv. .true.) &
382 24 : how_many_allocated_star_ids = how_many_allocated_star_ids+1
383 : end do
384 : end if
385 12 : end function how_many_allocated_star_ids
386 :
387 :
388 1 : subroutine free_star(s)
389 : type (star_info), pointer :: s
390 1 : star_handles(s% id)% in_use = .false.
391 1 : end subroutine free_star
392 :
393 :
394 1 : subroutine stardata_init( &
395 : my_mesa_dir, chem_isotopes_filename, &
396 : net_reaction_filename, jina_reaclib_filename, &
397 : use_suzuki_weak_rates, &
398 : use_special_weak_rates, special_weak_states_file, special_weak_transitions_file, &
399 : reaclib_min_T9_in, &
400 : rate_tables_dir, rates_cache_suffix, &
401 : ionization_file_prefix, ionization_Z1_suffix, &
402 : eosDT_cache_dir, &
403 : ionization_cache_dir, kap_cache_dir, rates_cache_dir, &
404 1 : color_num_files,color_file_names,color_num_colors,&
405 : ierr)
406 : use colors_lib, only : colors_init
407 : use kap_lib, only: kap_init
408 : use eos_lib, only: eos_init
409 : use rates_lib, only: rates_init
410 : use rates_def, only: reaclib_min_T9
411 : use net_lib, only: net_init
412 : use ionization_lib, only: ionization_init
413 : use atm_lib
414 : use chem_lib
415 : use const_def, only: dp, mesa_data_dir
416 : use const_lib, only: const_init
417 : use utils_lib
418 : use star_history_def, only: history_column_names_init
419 : use star_profile_def, only: profile_column_names_init
420 : character (len=*), intent(in) :: &
421 : my_mesa_dir, chem_isotopes_filename, net_reaction_filename, &
422 : jina_reaclib_filename, rate_tables_dir, &
423 : special_weak_states_file, special_weak_transitions_file, &
424 : rates_cache_suffix, &
425 : ionization_file_prefix, ionization_Z1_suffix, &
426 : eosDT_cache_dir, &
427 : ionization_cache_dir, kap_cache_dir, rates_cache_dir
428 : integer, intent(in) :: color_num_files
429 : character (len=*), intent(in) :: color_file_names(:)
430 : integer, intent(in) :: color_num_colors(:)
431 : real(dp), intent(in) :: reaclib_min_T9_in
432 : logical, intent(in) :: use_suzuki_weak_rates, use_special_weak_rates
433 : integer, intent(out) :: ierr
434 :
435 : logical, parameter :: use_cache = .false.
436 : character (len=strlen) :: fname
437 : integer :: iounit
438 :
439 : logical, parameter :: dbg = .false.
440 :
441 : include 'formats'
442 :
443 1 : ierr = 0
444 :
445 1 : rate_tables_dir_for_star = rate_tables_dir
446 1 : rates_cache_suffix_for_star = rates_cache_suffix
447 :
448 : if (dbg) write(*,*) 'call const_init'
449 1 : call const_init(my_mesa_dir,ierr)
450 1 : if (ierr /= 0) return
451 :
452 : if (dbg) write(*,*) 'call math_init'
453 1 : call math_init()
454 :
455 1 : call star_private_def_init
456 1 : call result_reason_init
457 :
458 1 : call history_column_names_init(ierr)
459 1 : if (ierr /= 0) return
460 :
461 1 : call profile_column_names_init(ierr)
462 1 : if (ierr /= 0) return
463 :
464 : if (dbg) write(*,*) 'call chem_init'
465 1 : call chem_init(chem_isotopes_filename, ierr)
466 1 : if (ierr /= 0) return
467 :
468 : if (dbg) write(*,*) 'call colors_init'
469 1 : call colors_init(color_num_files,color_file_names,color_num_colors,ierr)
470 1 : if (ierr /= 0) return
471 :
472 : if (dbg) write(*,*) 'call eos_init'
473 : !write(*,*) 'eos_file_prefix "' // trim(eos_file_prefix) // '"'
474 : !write(*,*) 'eosDT_cache_dir "' // trim(eosDT_cache_dir) // '"'
475 : call eos_init( &
476 1 : eosDT_cache_dir, use_cache, ierr)
477 1 : if (ierr /= 0) return
478 :
479 : if (dbg) write(*,*) 'call kap_init'
480 : !write(*,*) 'kap_cache_dir "' // trim(kap_cache_dir) // '"'
481 1 : call kap_init(use_cache, kap_cache_dir, ierr)
482 1 : if (ierr /= 0) return
483 :
484 : if (dbg) write(*,*) 'call rates_init'
485 : call rates_init( &
486 : net_reaction_filename, jina_reaclib_filename, &
487 : rate_tables_dir_for_star, &
488 : use_suzuki_weak_rates, &
489 : use_special_weak_rates, &
490 : special_weak_states_file, &
491 : special_weak_transitions_file, &
492 : rates_cache_dir, &
493 1 : ierr)
494 :
495 1 : if (ierr /= 0) return
496 :
497 1 : if (reaclib_min_T9_in > 0 .and. reaclib_min_T9_in /= reaclib_min_T9) then
498 0 : reaclib_min_T9 = reaclib_min_T9_in
499 0 : write(*,'(A)')
500 0 : write(*,'(A)')
501 0 : write(*,'(A)')
502 0 : write(*,'(A)')
503 0 : write(*,1) 'change reaclib_min_T9', reaclib_min_T9
504 0 : write(*,1) 'must clear data/rates_data/cache of old reaclib rates'
505 0 : write(*,'(A)')
506 0 : write(*,'(A)')
507 0 : write(*,'(A)')
508 0 : write(*,'(A)')
509 0 : write(*,'(A)')
510 : end if
511 :
512 : if (dbg) write(*,*) 'call net_init'
513 1 : call net_init(ierr)
514 1 : if (ierr /= 0) return
515 :
516 : if (dbg) write(*,*) 'call atm_init'
517 1 : call atm_init(use_cache, ierr)
518 1 : if (ierr /= 0) return
519 :
520 : if (dbg) write(*,*) 'call ionization_init'
521 : call ionization_init( &
522 : ionization_file_prefix, ionization_Z1_suffix, &
523 1 : ionization_cache_dir, use_cache, ierr)
524 1 : if (ierr /= 0) return
525 :
526 1 : version_number = ''
527 1 : fname = trim(mesa_data_dir) // '/version_number'
528 1 : open(newunit=iounit, file=trim(fname), action='read', status='old', iostat=ierr)
529 1 : if (ierr /= 0) then
530 0 : write(*, *) 'failed to open ' // trim(fname)
531 0 : return
532 : end if
533 :
534 1 : read(iounit, '(A)', iostat=ierr) version_number
535 :
536 1 : if (ierr /= 0) then
537 0 : close(iounit)
538 0 : return
539 : end if
540 :
541 1 : close(iounit)
542 :
543 1 : write(*,'(1x,a,1x,a)') 'version_number', trim(version_number)
544 :
545 : ! here we store useful information about the compiler and SDK
546 1 : call get_compiler_version(compiler_name,compiler_version_name)
547 1 : call get_mesasdk_version(mesasdk_version_name,ierr)
548 1 : call date_and_time(date=date)
549 :
550 2 : end subroutine stardata_init
551 :
552 : end module star_private_def
|
