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 photo_in
21 :
22 : use star_private_def
23 :
24 : implicit none
25 :
26 : private
27 : public :: read_star_photo
28 :
29 : contains
30 :
31 0 : subroutine read_star_photo(s, fname, ierr)
32 : use utils_lib, only: integer_dict_define, integer_dict_create_hash
33 : use rates_def, only: num_rvs
34 : use alloc, only: set_var_info, allocate_star_info_arrays, &
35 : alloc_star_info_old_arrays
36 : use rsp_def, only: rsp_photo_in
37 : type (star_info), pointer :: s
38 : character (len=*), intent(in) :: fname
39 : integer, intent(out) :: ierr
40 :
41 : integer :: iounit, i, j, k, version, part_number, &
42 : len_history_col_spec, nz
43 : logical, parameter :: dbg = .false.
44 :
45 : include 'formats'
46 :
47 0 : ierr = 0
48 0 : part_number = 0 ! part_numbers are just a consistency check
49 :
50 0 : write(*, *) 'read ', trim(fname)
51 : open(newunit=iounit, file=trim(fname), action='read', &
52 0 : status='old', iostat=ierr, form='unformatted')
53 0 : if (ierr /= 0) then
54 0 : if (s% report_ierr) write(*, *) 'Failed to open ', trim(fname)
55 0 : return
56 : end if
57 :
58 0 : read(iounit, iostat=ierr) version
59 0 : if (failed('version')) return
60 0 : if (version /= star_def_version) then
61 : write(*,'(/,a,/)') ' FAILURE: the restart data' // &
62 0 : ' is from a previous version of the code and is no longer usable.'
63 0 : ierr = -1
64 0 : return
65 : end if
66 :
67 0 : call read_part_number(iounit)
68 0 : if (failed('initial_z')) return
69 :
70 : read(iounit, iostat=ierr) &
71 0 : s% initial_z, & ! need this since read_model can change what is in the inlist
72 0 : s% total_num_solver_iterations, &
73 0 : s% nz, s% nvar_hydro, s% nvar_chem, s% nvar_total, &
74 0 : s% v_flag, s% u_flag, s% rotation_flag, s% RSP2_flag, s% RSP_flag, &
75 0 : s% RTI_flag, s% w_div_wc_flag, s% j_rot_flag, s% D_omega_flag, s% am_nu_rot_flag, &
76 0 : s% have_mlt_vc, s% species, s% num_reactions, &
77 0 : s% model_number, s% star_mass, &
78 0 : s% mstar, s% xmstar, s% M_center, s% v_center, s% R_center, s% L_center, &
79 0 : s% time, s% dt, s% have_previous_conv_vel, &
80 0 : s% was_in_implicit_wind_limit, &
81 0 : s% using_revised_net_name, &
82 0 : s% revised_net_name, &
83 0 : s% using_revised_max_yr_dt, &
84 0 : s% revised_max_yr_dt, &
85 0 : s% astero_using_revised_max_yr_dt, &
86 0 : s% astero_revised_max_yr_dt, &
87 0 : s% cumulative_energy_error, s% cumulative_extra_heating, &
88 0 : s% have_initial_energy_integrals, s% total_energy_initial, &
89 0 : s% force_tau_factor, s% force_Tsurf_factor, s% force_opacity_factor, &
90 0 : s% crystal_core_boundary_mass
91 :
92 0 : if (failed('initial_y')) return
93 0 : s% nz_old = s% nz ! needed by alloc
94 :
95 0 : if (s% force_tau_factor > 0 .and. s% tau_factor /= s% force_tau_factor .and. &
96 : s% tau_factor /= s% job% set_to_this_tau_factor) then
97 0 : s% tau_factor = s% force_tau_factor
98 0 : write(*,1) 'set tau_factor to photo value', s% tau_factor
99 : end if
100 :
101 0 : if (s% force_Tsurf_factor > 0 .and. s% Tsurf_factor /= s% force_Tsurf_factor .and. &
102 : s% Tsurf_factor /= s% job% set_to_this_Tsurf_factor) then
103 0 : s% Tsurf_factor = s% force_Tsurf_factor
104 0 : write(*,1) 'set Tsurf_factor to photo value', s% Tsurf_factor
105 : end if
106 :
107 0 : if (s% force_opacity_factor > 0 .and. s% opacity_factor /= s% force_opacity_factor .and. &
108 : s% opacity_factor /= s% job% relax_to_this_opacity_factor) then
109 0 : s% opacity_factor = s% force_opacity_factor
110 0 : write(*,1) 'set opacity_factor to photo value', s% opacity_factor
111 : end if
112 :
113 0 : if (s% using_revised_net_name) &
114 0 : s% net_name = s% revised_net_name
115 :
116 0 : if (s% using_revised_max_yr_dt) &
117 0 : s% max_years_for_timestep = s% revised_max_yr_dt
118 :
119 0 : if (s% astero_using_revised_max_yr_dt) &
120 0 : s% max_years_for_timestep = s% astero_revised_max_yr_dt
121 :
122 0 : nz = s% nz
123 :
124 0 : read(iounit, iostat=ierr) s% net_name
125 0 : if (failed('net_name')) return
126 :
127 0 : call set_var_info(s, ierr)
128 0 : if (failed('set_var_info')) return
129 :
130 : ! allocate everything
131 0 : call allocate_star_info_arrays(s, ierr)
132 0 : if (failed('allocate_star_info_arrays')) return
133 :
134 0 : call alloc_star_info_old_arrays(s, ierr)
135 0 : if (failed('alloc_star_info_old_arrays')) return
136 :
137 0 : call read_part_number(iounit)
138 0 : if (failed('dq')) return
139 :
140 : read(iounit, iostat=ierr) &
141 0 : s% dq(1:nz), s% xa(:,1:nz), s% xh(:,1:nz), &
142 0 : s% omega(1:nz), s% j_rot(1:nz), s% mlt_vc(1:nz), s% conv_vel(1:nz), &
143 0 : s% D_ST_start(1:nz), s% nu_ST_start(1:nz), & ! needed for ST time smoothing
144 0 : s% have_ST_start_info
145 :
146 0 : call read_part_number(iounit)
147 0 : if (failed('rsp_num_periods')) return
148 :
149 : read(iounit, iostat=ierr) &
150 0 : s% rsp_num_periods, s% rsp_dt, s% rsp_period, s% RSP_have_set_velocities, &
151 0 : s% dt_limit_ratio, s% tau_base
152 0 : if (failed('cz_top_mass_old')) return
153 :
154 : read(iounit, iostat=ierr) &
155 0 : s% i_lnd, s% i_lnT, s% i_lnR, s% i_lum, s% i_Et_RSP, s% i_erad_RSP, s% i_Fr_RSP, &
156 0 : s% i_v, s% i_u, s% i_alpha_RTI, s% i_w, s% i_Hp, s% i_w_div_wc, s% i_j_rot, &
157 0 : s% i_dv_dt, s% i_equL, s% i_dlnd_dt, s% i_dlnE_dt, &
158 0 : s% i_dEt_RSP_dt, s% i_derad_RSP_dt, s% i_dFr_RSP_dt, s% i_du_dt, s% i_dlnR_dt, &
159 0 : s% i_dalpha_RTI_dt, s% i_detrb_dt, s% i_equ_Hp
160 0 : if (failed('i_dalpha_RTI_dt')) return
161 :
162 : read(iounit, iostat=ierr) &
163 0 : s% model_controls_filename, s% model_data_filename, &
164 0 : s% most_recent_profile_filename, s% most_recent_controls_filename, &
165 0 : s% most_recent_model_data_filename
166 0 : if (failed('most_recent_model_data_filename')) return
167 :
168 0 : call read_part_number(iounit)
169 0 : if (failed('recent_log_header')) return
170 :
171 : read(iounit, iostat=ierr) &
172 0 : s% recent_log_header, s% phase_of_evolution, s% dt_next, s% dt_next_unclipped
173 0 : if (failed('eps_nuc_neu_total')) return
174 0 : if (s% dt_next <= 0d0) s% dt_next = s% dt_next_unclipped
175 :
176 0 : call read_part_number(iounit)
177 0 : if (failed('read_part_number')) return
178 :
179 : read(iounit, iostat=ierr) &
180 0 : s% num_skipped_setvars, s% num_retries, s% num_setvars, &
181 0 : s% total_num_solver_iterations, s% total_num_solver_relax_iterations, &
182 0 : s% total_num_solver_calls_made, s% total_num_solver_relax_calls_made, &
183 0 : s% total_num_solver_calls_converged, s% total_num_solver_relax_calls_converged, &
184 0 : s% total_step_attempts, s% total_relax_step_attempts, &
185 0 : s% total_step_retries, s% total_relax_step_retries, &
186 0 : s% total_step_redos, s% total_relax_step_redos, &
187 0 : s% total_steps_finished, s% total_relax_steps_finished, &
188 0 : s% num_hydro_merges, s% num_hydro_splits, s% num_solver_setvars, &
189 0 : s% mesh_call_number, s% solver_call_number, s% diffusion_call_number, &
190 0 : s% gradT_excess_alpha, s% gradT_excess_alpha_old, s% Teff, s% mstar_dot, &
191 0 : s% power_nuc_burn, s% power_h_burn, s% power_he_burn, s% power_z_burn, s% power_photo, &
192 0 : s% why_Tlim, s% dt_why_count(1:numTlim), s% dt_why_retry_count(1:numTlim), &
193 0 : s% timestep_hold, s% model_number_for_last_retry, s% model_number_for_last_retry_old, &
194 0 : s% init_model_number, s% most_recent_photo_name, &
195 0 : s% rand_i97, s% rand_j97, s% rand_u(1:rand_u_len), s% rand_c, s% rand_cd, s% rand_cm
196 0 : if (failed('most_recent_photo_name')) return
197 :
198 0 : call read_part_number(iounit)
199 0 : if (failed('len_extra_iwork')) return
200 :
201 0 : read(iounit, iostat=ierr) s% len_extra_iwork, s% len_extra_work
202 0 : if (failed('len_extra_work')) return
203 :
204 0 : if (s% len_extra_iwork > 0) then
205 : allocate( &
206 : s% extra_iwork(s% len_extra_iwork), &
207 : s% extra_iwork_old(s% len_extra_iwork), &
208 0 : stat=ierr)
209 0 : if (failed('allocate extra_iwork')) return
210 0 : read(iounit, iostat=ierr) s% extra_iwork(1:s% len_extra_iwork)
211 0 : if (failed('read extra_iwork')) return
212 : !read(iounit, iostat=ierr) s% extra_iwork_old(1:s% len_extra_iwork)
213 : !if (failed('allocate extra_iwork_old')) return
214 : else
215 0 : nullify(s% extra_iwork, s% extra_iwork_old)
216 : end if
217 :
218 0 : if (s% len_extra_work > 0) then
219 : allocate( &
220 : s% extra_work(s% len_extra_work), &
221 : s% extra_work_old(s% len_extra_work), &
222 0 : stat=ierr)
223 0 : if (failed('allocate extra_work')) return
224 0 : read(iounit, iostat=ierr) s% extra_work(1:s% len_extra_work)
225 0 : if (failed('read extra_work')) return
226 : !read(iounit, iostat=ierr) s% extra_work_old(1:s% len_extra_work)
227 : !if (failed('read extra_work_old')) return
228 : else
229 0 : nullify(s% extra_work, s% extra_work_old)
230 : end if
231 :
232 0 : read(iounit, iostat=ierr) s% ixtra
233 0 : if (failed('ixtra')) return
234 0 : read(iounit, iostat=ierr) s% xtra
235 0 : if (failed('xtra')) return
236 0 : read(iounit, iostat=ierr) s% lxtra
237 0 : if (failed('lxtra')) return
238 :
239 0 : read(iounit, iostat=ierr) len_history_col_spec
240 0 : if (failed('len_history_col_spec')) return
241 0 : if (len_history_col_spec > 0) then
242 0 : allocate(s% history_column_spec(len_history_col_spec), stat=ierr)
243 0 : if (failed('alloc history_column_spec')) return
244 0 : read(iounit, iostat=ierr) s% history_column_spec(1:len_history_col_spec)
245 0 : if (failed('read history_column_spec')) return
246 : end if
247 :
248 : read(iounit, iostat=ierr) &
249 0 : s% number_of_history_columns, s% model_number_of_history_values, &
250 0 : s% need_to_set_history_names_etc
251 0 : if (failed('number_of_history_columns')) return
252 :
253 0 : if (s% number_of_history_columns > 0) then
254 :
255 0 : allocate(s% history_value_is_integer(s% number_of_history_columns), stat=ierr)
256 0 : if (failed('alloc history_value_is_integer')) return
257 0 : read(iounit, iostat=ierr) s% history_value_is_integer(1:s% number_of_history_columns)
258 0 : if (failed('read history_value_is_integer')) return
259 :
260 0 : allocate(s% history_names(s% number_of_history_columns), stat=ierr)
261 0 : if (failed('alloc history_names')) return
262 0 : do k=1,s% number_of_history_columns
263 0 : read(iounit, iostat=ierr) s% history_names(k)
264 0 : if (failed('read history_names')) return
265 : end do
266 :
267 : ! rebuild the history_names_dict
268 0 : do j = 1, s% number_of_history_columns
269 0 : call integer_dict_define(s% history_names_dict, s% history_names(j), j, ierr)
270 0 : if (failed('integer_dict_define history_names_dict')) return
271 : end do
272 0 : call integer_dict_create_hash(s% history_names_dict, ierr)
273 0 : if (failed('integer_dict_create_hash history_names_dict')) return
274 :
275 : end if
276 :
277 0 : if (s% rsp_flag) then
278 0 : call rsp_photo_in(s, iounit, ierr)
279 0 : if (failed('after rsp_photo_in')) return
280 : end if
281 :
282 0 : call read_part_number(iounit)
283 0 : if (failed('before other_photo_read')) return
284 :
285 0 : call s% other_photo_read(s% id, iounit, ierr)
286 0 : if (failed('after other_photo_read')) return
287 :
288 0 : call read_part_number(iounit)
289 0 : if (failed('final read_part_number')) return
290 :
291 0 : s% need_to_setvars = .true. ! set this after photo out or photo in
292 :
293 0 : close(iounit)
294 :
295 : contains
296 :
297 0 : subroutine read_part_number(iounit)
298 : integer, intent(in) :: iounit
299 : integer :: i
300 0 : part_number = part_number + 1
301 0 : read(iounit, iostat=ierr) i
302 0 : if (ierr /= 0) return
303 0 : if (i /= part_number) ierr = -1
304 : !write(*,*) 'part_number', part_number
305 0 : end subroutine read_part_number
306 :
307 0 : logical function failed(str)
308 : character (len=*), intent(in) :: str
309 0 : i = i+1
310 0 : if (ierr /= 0) then
311 0 : write(*, *) 'read_star_photo failed for ' // trim(str)
312 0 : failed = .true.
313 0 : return
314 : end if
315 0 : failed = .false.
316 : end function failed
317 :
318 : end subroutine read_star_photo
319 :
320 : end module photo_in
|
