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 read_model
21 :
22 : use star_private_def
23 : use const_def, only: dp, msun, secyer
24 :
25 : implicit none
26 :
27 : integer, parameter :: bit_for_zams_file = 0
28 : integer, parameter :: bit_for_lnPgas = 1 ! OBSOLETE: includes lnPgas variables in place of lnd
29 : integer, parameter :: bit_for_2models = 2
30 : integer, parameter :: bit_for_velocity = 3
31 : integer, parameter :: bit_for_rotation = 4
32 : integer, parameter :: bit_for_mlt_vc = 5
33 : integer, parameter :: bit_for_RSP2 = 6
34 : integer, parameter :: bit_for_RTI = 7
35 : !integer, parameter :: = 8
36 : integer, parameter :: bit_for_u = 9
37 : integer, parameter :: bit_for_D_omega = 10
38 : integer, parameter :: bit_for_am_nu_rot = 11
39 : integer, parameter :: bit_for_j_rot = 12
40 : !integer, parameter :: = 13
41 : !integer, parameter :: = 14
42 : integer, parameter :: bit_for_RSP = 15
43 : integer, parameter :: bit_for_no_L_basic_variable = 16
44 :
45 : integer, parameter :: increment_for_rotation_flag = 1
46 : integer, parameter :: increment_for_have_j_rot = 1
47 : integer, parameter :: increment_for_have_mlt_vc = 1
48 : integer, parameter :: increment_for_D_omega_flag = 1
49 : integer, parameter :: increment_for_am_nu_rot_flag = 1
50 : integer, parameter :: increment_for_RTI_flag = 1
51 : integer, parameter :: increment_for_RSP_flag = 3
52 : integer, parameter :: increment_for_RSP2_flag = 1
53 :
54 : integer, parameter :: max_increment = increment_for_rotation_flag &
55 : + increment_for_have_j_rot &
56 : + increment_for_have_mlt_vc &
57 : + increment_for_D_omega_flag &
58 : + increment_for_am_nu_rot_flag &
59 : + increment_for_RTI_flag &
60 : + increment_for_RSP_flag &
61 : + increment_for_RSP2_flag
62 :
63 : integer, parameter :: mesa_zams_file_type = 2**bit_for_zams_file
64 :
65 : character (len=100000) :: buf
66 :
67 : contains
68 :
69 18 : subroutine finish_load_model(s, restart, ierr)
70 : use hydro_vars, only: set_vars
71 : use star_utils, only: set_m_and_dm, set_m_grav_and_grav, set_dm_bar, &
72 : total_angular_momentum, reset_epsnuc_vectors, set_qs
73 : use hydro_rotation, only: use_xh_to_update_i_rot_and_j_rot, &
74 : set_i_rot_from_omega_and_j_rot, use_xh_to_update_i_rot, set_rotation_info
75 : use hydro_RSP2, only: set_RSP2_vars
76 : use RSP, only: RSP_setup_part1, RSP_setup_part2
77 : use report, only: do_report
78 : use alloc, only: fill_ad_with_zeros
79 : use brunt, only: do_brunt_B, do_brunt_N2
80 : type (star_info), pointer :: s
81 : logical, intent(in) :: restart
82 : integer, intent(out) :: ierr
83 : integer :: k, nz
84 : include 'formats'
85 : ierr = 0
86 1 : nz = s% nz
87 2465 : s% brunt_B(1:nz) = 0 ! temporary proxy for brunt_B
88 1 : call set_qs(s, nz, s% q, s% dq, ierr)
89 1 : if (ierr /= 0) then
90 0 : write(*,*) 'set_qs failed in finish_load_model'
91 0 : return
92 : end if
93 1 : call set_m_and_dm(s)
94 1 : call set_m_grav_and_grav(s)
95 1 : call set_dm_bar(s, nz, s% dm, s% dm_bar)
96 :
97 1 : call reset_epsnuc_vectors(s)
98 :
99 1 : s% star_mass = s% mstar/msun
100 :
101 1 : if (s% rotation_flag) then
102 : ! older MESA versions stored only omega in saved models. However, when
103 : ! using rotation dependent moments of inertia one actually needs to store
104 : ! the angular momentum in order to initialize the model. This flag is here
105 : ! to account for the loading of old saved models.
106 0 : if (s% have_j_rot) then
107 : !if (restart) then
108 : ! only need to compute irot, w_div_w_crit_roche is stored in photos
109 : ! call set_i_rot_from_omega_and_j_rot(s)
110 : !else
111 : ! need to set w_div_w_crit_roche as well
112 0 : call use_xh_to_update_i_rot(s)
113 0 : do k=1, s% nz
114 0 : s% omega(k) = s% j_rot(k)/s% i_rot(k)% val
115 : end do
116 : !end if
117 : else
118 : ! need to recompute irot and jrot
119 0 : call use_xh_to_update_i_rot_and_j_rot(s)
120 : end if
121 : ! this ensures fp, ft, r_equatorial and r_polar are set by the end
122 : !call set_rotation_info(s, .true., ierr)
123 : !if (ierr /= 0) then
124 : ! write(*,*) &
125 : ! 'finish_load_model failed in set_rotation_info'
126 : ! return
127 : !end if
128 : end if
129 :
130 : ! clear some just to avoid getting NaNs at start
131 : ! e.g., from profile_starting_model
132 2465 : s% D_mix(1:nz) = 0
133 2465 : s% adjust_mlt_gradT_fraction(1:nz) = -1
134 2465 : s% eps_mdot(1:nz) = 0
135 2465 : s% dvc_dt_TDC(1:nz) = 0
136 1 : call fill_ad_with_zeros(s% eps_grav_ad,1,-1)
137 2465 : s% ergs_error(1:nz) = 0
138 1 : if (.not. restart) s% have_ST_start_info = .false.
139 1 : if (s% do_element_diffusion) s% edv(:,1:nz) = 0
140 1 : if (s% u_flag) then
141 0 : call fill_ad_with_zeros(s% u_face_ad,1,-1)
142 0 : call fill_ad_with_zeros(s% P_face_ad,1,-1)
143 : end if
144 :
145 2465 : s% flux_limit_R(1:nz) = 0
146 2465 : s% flux_limit_lambda(1:nz) = 0
147 :
148 1 : if (s% RSP_flag) then
149 0 : call RSP_setup_part1(s, restart, ierr)
150 0 : if (ierr /= 0) then
151 0 : write(*,*) 'finish_load_model: RSP_setup_part1 returned ierr', ierr
152 0 : return
153 : end if
154 : end if
155 :
156 1 : if (.not. s% have_mlt_vc) then
157 1 : s% okay_to_set_mlt_vc = .true.
158 : end if
159 :
160 1 : s% doing_finish_load_model = .true.
161 1 : call set_vars(s, s% dt, ierr)
162 1 : if (ierr == 0 .and. s% RSP2_flag) call set_RSP2_vars(s,ierr)
163 1 : s% doing_finish_load_model = .false.
164 1 : if (ierr /= 0) then
165 0 : write(*,*) 'finish_load_model: failed in set_vars'
166 0 : return
167 : end if
168 :
169 1 : if (s% rotation_flag) s% total_angular_momentum = total_angular_momentum(s)
170 :
171 1 : if (s% RSP_flag) then
172 0 : call RSP_setup_part2(s, restart, ierr)
173 0 : if (ierr /= 0) then
174 0 : write(*,*) 'finish_load_model: RSP_setup_part2 returned ierr', ierr
175 0 : return
176 : end if
177 : end if
178 :
179 1 : s% doing_finish_load_model = .true.
180 :
181 1 : if(s% calculate_Brunt_B) call do_brunt_B(s, 1, s%nz, ierr)
182 1 : if (ierr /= 0) then
183 0 : write(*,*) 'finish_load_model: failed in do_brunt_b'
184 0 : return
185 : end if
186 :
187 1 : if(s% calculate_Brunt_N2) call do_brunt_N2(s, 1, s%nz, ierr)
188 1 : if (ierr /= 0) then
189 0 : write(*,*) 'finish_load_model: failed in do_brunt_N2'
190 0 : return
191 : end if
192 :
193 1 : call do_report(s, ierr)
194 1 : s% doing_finish_load_model = .false.
195 1 : if (ierr /= 0) then
196 0 : write(*,*) 'finish_load_model: failed in do_report'
197 0 : return
198 : end if
199 :
200 1 : end subroutine finish_load_model
201 :
202 :
203 0 : subroutine do_read_saved_model(s, filename, ierr)
204 1 : use utils_lib
205 : use utils_def
206 : use chem_def
207 : use net, only: set_net
208 : use alloc, only: set_var_info, &
209 : free_star_info_arrays, allocate_star_info_arrays, set_chem_names
210 : use star_utils, only: yrs_for_init_timestep, set_phase_of_evolution
211 : type (star_info), pointer :: s
212 : character (len=*), intent(in) :: filename
213 : integer, intent(out) :: ierr
214 :
215 : integer :: iounit, n, i, t, file_type, &
216 : year_month_day_when_created, nz, species, nvar, count
217 : logical :: do_read_prev, no_L
218 0 : real(dp) :: initial_mass, initial_z, initial_y, &
219 0 : tau_factor, Tsurf_factor, opacity_factor, mixing_length_alpha
220 : character (len=strlen) :: buffer, string
221 : character (len=net_name_len) :: net_name
222 0 : character(len=iso_name_length), pointer :: names(:) ! (species)
223 0 : integer, pointer :: perm(:) ! (species)
224 :
225 : include 'formats'
226 :
227 0 : ierr = 0
228 0 : open(newunit=iounit, file=trim(filename), status='old', action='read', iostat=ierr)
229 0 : if (ierr /= 0) then
230 0 : write(*,*) 'open failed', ierr, iounit
231 0 : write(*, '(a)') 'failed to open ' // trim(filename)
232 0 : return
233 : end if
234 :
235 : ! use token to get file_type so can have comments at start of file
236 0 : n = 0
237 0 : i = 0
238 0 : t = token(iounit, n, i, buffer, string)
239 0 : if (t == eof_token) then
240 0 : write(*, '(a)') 'failed to find file type at start of ' // trim(filename)
241 0 : return
242 : end if
243 0 : if (t /= name_token) then
244 0 : write(*, '(a)') 'failed to find file type at start of ' // trim(filename)
245 0 : return
246 : end if
247 0 : read(string,fmt=*,iostat=ierr) file_type
248 0 : if (ierr /= 0) then
249 0 : write(*, '(a)') 'failed to find file type at start of ' // trim(filename)
250 0 : return
251 : end if
252 :
253 0 : read(iounit, *, iostat=ierr) ! skip the blank line after the file type
254 0 : if (ierr /= 0) then
255 : return
256 : end if
257 :
258 : ! refuse to load old models using lnPgas as a structure variable
259 0 : if (BTEST(file_type, bit_for_lnPgas)) then
260 0 : write(*,'(A)')
261 0 : write(*,*) 'MESA no longer supports models using lnPgas as a structure variable'
262 0 : write(*,'(A)')
263 0 : ierr = -1
264 0 : return
265 : end if
266 :
267 0 : s% model_number = 0
268 0 : s% star_age = 0
269 0 : s% xmstar = -1
270 :
271 0 : tau_factor = s% tau_factor
272 0 : Tsurf_factor = s% Tsurf_factor
273 0 : mixing_length_alpha = s% mixing_length_alpha
274 0 : opacity_factor = s% opacity_factor
275 :
276 : call read_properties(iounit, &
277 : net_name, species, nz, year_month_day_when_created, &
278 : initial_mass, initial_z, initial_y, mixing_length_alpha, &
279 : s% model_number, s% star_age, tau_factor, s% Teff, &
280 : s% power_nuc_burn, s% power_h_burn, s% power_he_burn, s% power_z_burn, s% power_photo, &
281 : Tsurf_factor, opacity_factor, s% crystal_core_boundary_mass, &
282 : s% xmstar, s% R_center, s% L_center, s% v_center, &
283 0 : s% cumulative_energy_error, s% num_retries, ierr)
284 :
285 : if (ierr /= 0 .or. initial_mass < 0 .or. nz < 0 &
286 : .or. initial_z < 0 .or. species < 0 .or. &
287 0 : is_bad(s% xmstar) .or. &
288 : is_bad(initial_mass + initial_z)) then
289 0 : ierr = -1
290 0 : write(*, *) 'do_read_model: missing required properties'
291 0 : write(*,*) 'initial_mass', initial_mass
292 0 : write(*,*) 'xmstar', s% xmstar
293 0 : write(*,*) 'initial_z', initial_z
294 0 : write(*,*) 'nz', nz
295 0 : write(*,*) 'species', species
296 0 : return
297 : end if
298 :
299 0 : s% init_model_number = s% model_number
300 0 : s% time = s% star_age*secyer
301 :
302 0 : if (abs(tau_factor - s% tau_factor) > tau_factor*1d-9 .and. &
303 : s% tau_factor /= s% job% set_to_this_tau_factor) then
304 : ! don't change if just set by inlist
305 0 : write(*,'(A)')
306 0 : write(*,1) 'WARNING: changing to saved tau_factor =', tau_factor
307 0 : write(*,'(A)')
308 0 : s% tau_factor = tau_factor
309 0 : s% force_tau_factor = tau_factor
310 : end if
311 :
312 0 : if (abs(Tsurf_factor - s% Tsurf_factor) > Tsurf_factor*1d-9 .and. &
313 : s% Tsurf_factor /= s% job% set_to_this_Tsurf_factor) then
314 : ! don't change if just set by inlist
315 0 : write(*,'(A)')
316 0 : write(*,1) 'WARNING: changing to saved Tsurf_factor =', Tsurf_factor
317 0 : write(*,'(A)')
318 0 : s% Tsurf_factor = Tsurf_factor
319 0 : s% force_Tsurf_factor = Tsurf_factor
320 : end if
321 :
322 0 : if (abs(opacity_factor - s% opacity_factor) > opacity_factor*1d-9 .and. &
323 : s% opacity_factor /= s% job% relax_to_this_opacity_factor) then
324 : ! don't change if just set by inlist
325 0 : write(*,'(A)')
326 0 : write(*,1) 'WARNING: changing to saved opacity_factor =', opacity_factor
327 0 : write(*,'(A)')
328 0 : s% opacity_factor = opacity_factor
329 0 : s% force_opacity_factor = opacity_factor
330 : end if
331 :
332 0 : if (abs(mixing_length_alpha - s% mixing_length_alpha) > mixing_length_alpha*1d-9) then
333 0 : write(*,'(A)')
334 0 : write(*,1) 'WARNING: model saved with mixing_length_alpha =', mixing_length_alpha
335 0 : write(*,1) 'but current setting for mixing_length_alpha =', s% mixing_length_alpha
336 0 : write(*,'(A)')
337 : end if
338 :
339 0 : s% v_flag = BTEST(file_type, bit_for_velocity)
340 0 : s% u_flag = BTEST(file_type, bit_for_u)
341 0 : s% rotation_flag = BTEST(file_type, bit_for_rotation)
342 0 : s% have_j_rot = BTEST(file_type, bit_for_j_rot)
343 0 : s% have_mlt_vc = BTEST(file_type, bit_for_mlt_vc)
344 0 : s% D_omega_flag = BTEST(file_type, bit_for_D_omega)
345 0 : s% am_nu_rot_flag = BTEST(file_type, bit_for_am_nu_rot)
346 0 : s% RTI_flag = BTEST(file_type, bit_for_RTI)
347 0 : s% RSP_flag = BTEST(file_type, bit_for_RSP)
348 0 : s% RSP2_flag = BTEST(file_type, bit_for_RSP2)
349 0 : no_L = BTEST(file_type, bit_for_no_L_basic_variable)
350 :
351 : if (BTEST(file_type, bit_for_lnPgas)) then
352 : write(*,'(A)')
353 : write(*,*) 'MESA no longer supports models using lnPgas as a structure variable'
354 : write(*,'(A)')
355 : ierr = -1
356 : return
357 : end if
358 :
359 0 : s% net_name = trim(net_name)
360 0 : s% species = species
361 0 : s% initial_z = initial_z
362 :
363 0 : s% mstar = initial_mass*Msun
364 0 : if (s% xmstar < 0) then
365 0 : s% M_center = 0
366 0 : s% xmstar = s% mstar
367 : else
368 0 : s% M_center = s% mstar - s% xmstar
369 : end if
370 0 : if (is_bad(s% M_center)) then
371 0 : write(*,1) 'M_center mstar xmstar initial_mass', &
372 0 : s% M_center, s% mstar, s% xmstar, initial_mass
373 0 : call mesa_error(__FILE__,__LINE__,'do_read_saved_model')
374 : end if
375 :
376 0 : call set_net(s, s% net_name, ierr)
377 0 : if (ierr /= 0) then
378 : write(*,*) &
379 0 : 'do_read_saved_model failed in set_net for net_name = ' // trim(s% net_name)
380 0 : return
381 : end if
382 :
383 0 : call set_var_info(s, ierr)
384 0 : if (ierr /= 0) then
385 0 : write(*,*) 'do_read_saved_model failed in set_var_info'
386 0 : return
387 : end if
388 :
389 : ! fixup chem names now that have nvar_hydro
390 0 : call set_chem_names(s)
391 :
392 0 : s% nz = nz
393 0 : call free_star_info_arrays(s)
394 0 : call allocate_star_info_arrays(s, ierr)
395 0 : if (ierr /= 0) then
396 0 : write(*,*) 'do_read_saved_model failed in allocate_star_info_arrays'
397 0 : return
398 : end if
399 :
400 0 : allocate(names(species), perm(species))
401 0 : call get_chem_col_names(s, iounit, species, names, perm, ierr)
402 0 : if (ierr /= 0) then
403 0 : deallocate(names, perm)
404 0 : write(*,*) 'do_read_saved_model failed in get_chem_col_names'
405 0 : return
406 : end if
407 :
408 0 : count = 0
409 0 : do i=1,species
410 0 : if (perm(i)==0) then
411 0 : count = count+1
412 0 : write(*,*) "Mod file has isotope ",trim(names(i)), " but that is not in the net"
413 : end if
414 : end do
415 0 : if (count/=0) call mesa_error(__FILE__,__LINE__)
416 :
417 0 : nvar = s% nvar_total
418 : call read1_model( &
419 : s, s% species, s% nvar_hydro, nz, iounit, &
420 : s% xh, s% xa, s% q, s% dq, s% omega, s% j_rot, &
421 0 : perm, ierr)
422 0 : deallocate(names, perm)
423 0 : if (ierr /= 0) then
424 0 : write(*,*) 'do_read_saved_model failed in read1_model'
425 0 : return
426 : end if
427 :
428 0 : do_read_prev = BTEST(file_type, bit_for_2models)
429 : if (ierr == 0) then
430 0 : if (do_read_prev) then
431 0 : call read_prev
432 : else
433 0 : s% generations = 1
434 : end if
435 : end if
436 :
437 0 : close(iounit)
438 :
439 :
440 : contains
441 :
442 :
443 0 : subroutine read_prev
444 : integer :: k
445 :
446 0 : do k=1, 3
447 0 : read(iounit, *, iostat=ierr)
448 0 : if (ierr /= 0) return
449 : end do
450 0 : call read_prev_properties
451 0 : if (ierr /= 0) return
452 :
453 : ! we do read_prev_properties to set initial timestep,
454 : ! but we don't use the previous model
455 : ! because we need to have other info about that isn't saved
456 : ! such as conv_vel and mixing_type
457 :
458 0 : s% generations = 1
459 :
460 0 : end subroutine read_prev
461 :
462 :
463 0 : subroutine read_prev_properties
464 : character (len=132) :: line
465 0 : real(dp) :: tmp, skip_val
466 : include 'formats'
467 :
468 0 : ierr = 0
469 0 : s% dt = -1
470 0 : s% mstar_old = -1
471 0 : s% dt_next = -1
472 0 : s% nz_old = -1
473 :
474 : do ! until reach a blank line
475 0 : read(iounit, fmt='(a)', iostat=ierr) line
476 0 : if (ierr /= 0) return
477 :
478 0 : if (len_trim(line) == 0) exit ! blank line
479 :
480 0 : if (match_keyword('previous n_shells', line, tmp)) then
481 0 : s% nz_old = int(tmp)
482 0 : cycle
483 : end if
484 :
485 0 : if (match_keyword('timestep (seconds)', line, s% dt)) then
486 : cycle
487 : end if
488 :
489 0 : if (match_keyword('previous mass (grams)', line, s% mstar_old)) then
490 : cycle
491 : end if
492 :
493 0 : if (match_keyword('dt_next (seconds)', line, s% dt_next)) then
494 : cycle
495 : end if
496 :
497 0 : if (match_keyword('year_month_day_when_created', line, skip_val)) cycle
498 :
499 : end do
500 0 : if (s% dt < 0) then
501 0 : ierr = -1
502 0 : write(*, *) 'missing dt for previous model'
503 : end if
504 0 : if (s% mstar_old < 0) then
505 0 : ierr = -1
506 0 : write(*, *) 'missing mstar_old for previous model'
507 : end if
508 0 : if (s% dt_next < 0) then
509 0 : ierr = -1
510 0 : write(*, *) 'missing dt_next for previous model'
511 : end if
512 :
513 : end subroutine read_prev_properties
514 :
515 :
516 : end subroutine do_read_saved_model
517 :
518 :
519 2 : subroutine read1_model( &
520 : s, species, nvar_hydro, nz, iounit, &
521 2 : xh, xa, q, dq, omega, j_rot, &
522 2 : perm, ierr)
523 : use star_utils, only: set_qs
524 : use chem_def
525 : type (star_info), pointer :: s
526 : integer, intent(in) :: species, nvar_hydro, nz, iounit, perm(:)
527 : real(dp), dimension(:,:), intent(out) :: xh, xa
528 : real(dp), dimension(:), intent(out) :: &
529 : q, dq, omega, j_rot
530 : integer, intent(out) :: ierr
531 :
532 : integer :: j, k, n, i_lnd, i_lnT, i_lnR, i_lum, i_w, i_Hp, &
533 : i_Et_RSP, i_erad_RSP, i_Fr_RSP, i_v, i_u, i_alpha_RTI, ii
534 46 : real(dp), target :: vec_ary(species + nvar_hydro + max_increment)
535 : real(dp), pointer :: vec(:)
536 : integer :: nvec
537 :
538 : include 'formats'
539 :
540 2 : ierr = 0
541 2 : vec => vec_ary
542 :
543 2 : i_lnd = s% i_lnd
544 2 : i_lnT = s% i_lnT
545 2 : i_lnR = s% i_lnR
546 2 : i_lum = s% i_lum
547 2 : i_w = s% i_w
548 2 : i_Hp = s% i_Hp
549 2 : i_v = s% i_v
550 2 : i_u = s% i_u
551 2 : i_alpha_RTI = s% i_alpha_RTI
552 2 : i_Et_RSP = s% i_Et_RSP
553 2 : i_erad_RSP = s% i_erad_RSP
554 2 : i_Fr_RSP = s% i_Fr_RSP
555 :
556 : n = species + nvar_hydro + 1 ! + 1 is for dq
557 : if (s% rotation_flag) n = n+increment_for_rotation_flag ! read omega
558 : if (s% have_j_rot) n = n+increment_for_have_j_rot ! read j_rot
559 : if (s% have_mlt_vc) n = n+increment_for_have_mlt_vc
560 : if (s% D_omega_flag) n = n+increment_for_D_omega_flag ! read D_omega
561 : if (s% am_nu_rot_flag) n = n+increment_for_am_nu_rot_flag ! read am_nu_rot
562 : if (s% RTI_flag) n = n+increment_for_RTI_flag ! read alpha_RTI
563 : if (s% RSP_flag) n = n+increment_for_RSP_flag ! read RSP_et, erad, Fr
564 : if (s% RSP2_flag) n = n+increment_for_RSP2_flag ! read w, Hp
565 :
566 4 : !$omp critical (read1_model_loop)
567 : ! make this a critical section to so don't have to dynamically allocate buf
568 4932 : do k = 1, nz
569 4930 : read(iounit,'(a)',iostat=ierr) buf
570 4930 : if (ierr /= 0) then
571 0 : write(*,3) 'read failed i', k, nz
572 0 : exit
573 : end if
574 4930 : call str_to_vector(buf, vec, nvec, ierr)
575 4930 : if (ierr /= 0) then
576 0 : write(*,*) 'str_to_vector failed'
577 0 : write(*,'(a,i8,1x,a)') 'buf', k, trim(buf)
578 0 : exit
579 : end if
580 4930 : j = int(vec(1))
581 4930 : if (j /= k) then
582 0 : ierr = -1
583 0 : write(*, *) 'error in reading model data j /= k'
584 0 : write(*, *) 'species', species
585 0 : write(*, *) 'j', j
586 0 : write(*, *) 'k', k
587 0 : write(*,'(a,1x,a)') 'buf', trim(buf)
588 0 : exit
589 : end if
590 : j = 1
591 4930 : j=j+1; xh(i_lnd,k) = vec(j)
592 4930 : j=j+1; xh(i_lnT,k) = vec(j)
593 4930 : j=j+1; xh(i_lnR,k) = vec(j)
594 4930 : if (s% RSP_flag) then
595 0 : j=j+1; xh(i_Et_RSP,k) = vec(j)
596 0 : j=j+1; xh(i_erad_RSP,k) = vec(j)
597 0 : j=j+1; xh(i_Fr_RSP,k) = vec(j)
598 4930 : else if (s% RSP2_flag) then
599 0 : j=j+1; xh(i_w,k) = vec(j)
600 0 : j=j+1; xh(i_Hp,k) = vec(j)
601 : end if
602 4930 : if (i_lum /= 0) then
603 4930 : j=j+1; xh(i_lum,k) = vec(j)
604 : else
605 0 : j=j+1; s% L(k) = vec(j)
606 : end if
607 4930 : j=j+1; dq(k) = vec(j)
608 4930 : if (s% v_flag) then
609 0 : j=j+1; xh(i_v,k) = vec(j)
610 : end if
611 4930 : if (s% rotation_flag) then
612 0 : j=j+1; omega(k) = vec(j)
613 : end if
614 4930 : if (s% have_j_rot) then
615 : !NOTE: MESA version 10108 was first to store j_rot in saved files
616 0 : j=j+1; j_rot(k) = vec(j)
617 : end if
618 4930 : if (s% D_omega_flag) then
619 0 : j=j+1; ! skip saving the file data
620 : end if
621 4930 : if (s% am_nu_rot_flag) then
622 0 : j=j+1; ! skip saving the file data
623 : end if
624 4930 : if (s% u_flag) then
625 0 : j=j+1; xh(i_u,k) = vec(j)
626 : end if
627 4930 : if (s% RTI_flag) then
628 0 : j=j+1; xh(i_alpha_RTI,k) = vec(j)
629 : end if
630 4930 : if (s% have_mlt_vc) then
631 0 : j=j+1; s% mlt_vc(k) = vec(j); s% conv_vel(k) = vec(j)
632 : end if
633 4930 : if (j+species > nvec) then
634 0 : ierr = -1
635 0 : write(*, *) 'error in reading model data j+species > nvec'
636 0 : write(*, *) 'j+species', j+species
637 0 : write(*, *) 'nvec', nvec
638 0 : write(*, *) 'j', j
639 0 : write(*, *) 'species', species
640 0 : write(*,'(a,1x,a)') 'buf', trim(buf)
641 0 : exit
642 : end if
643 49302 : do ii=1,species
644 44370 : xa(perm(ii),k) = vec(j+ii)
645 : end do
646 : end do
647 : !$omp end critical (read1_model_loop)
648 2 : if (ierr /= 0) then
649 0 : write(*,*) 'read1_model_loop failed'
650 0 : return
651 : end if
652 :
653 2 : if (s% rotation_flag .and. .not. s% D_omega_flag) &
654 0 : s% D_omega(1:nz) = 0d0
655 :
656 2 : if (s% rotation_flag .and. .not. s% am_nu_rot_flag) &
657 0 : s% am_nu_rot(1:nz) = 0d0
658 :
659 2 : call set_qs(s, nz, q, dq, ierr)
660 2 : if (ierr /= 0) then
661 0 : write(*,*) 'set_qs failed in read1_model sum(dq)', sum(dq(1:nz))
662 0 : return
663 : end if
664 :
665 2 : end subroutine read1_model
666 :
667 :
668 0 : subroutine do_read_saved_model_number(fname, model_number, ierr)
669 : character (len=*), intent(in) :: fname
670 : integer, intent(inout) :: model_number
671 : integer, intent(out) :: ierr
672 : character (len=strlen) :: net_name
673 : integer :: species, n_shells, &
674 : num_retries, year_month_day_when_created
675 0 : real(dp) :: m_div_msun, initial_z, &
676 0 : mixing_length_alpha, star_age, &
677 0 : Teff, tau_factor, Tsurf_factor, opacity_factor, crystal_core_boundary_mass, &
678 0 : power_nuc_burn, power_h_burn, power_he_burn, power_z_burn, power_photo, &
679 0 : xmstar, R_center, L_center, v_center, cumulative_energy_error
680 : call do_read_saved_model_properties(fname, &
681 : net_name, species, n_shells, year_month_day_when_created, &
682 : m_div_msun, initial_z, mixing_length_alpha, &
683 : model_number, star_age, tau_factor, Teff, &
684 : power_nuc_burn, power_h_burn, power_he_burn, power_z_burn, power_photo, &
685 : Tsurf_factor, opacity_factor, crystal_core_boundary_mass, &
686 : xmstar, R_center, L_center, v_center, &
687 0 : cumulative_energy_error, num_retries, ierr)
688 2 : end subroutine do_read_saved_model_number
689 :
690 :
691 0 : subroutine do_read_saved_model_properties(fname, &
692 : net_name, species, n_shells, year_month_day_when_created, &
693 : m_div_msun, initial_z, mixing_length_alpha, &
694 : model_number, star_age, tau_factor, Teff, &
695 : power_nuc_burn, power_h_burn, power_he_burn, power_z_burn, power_photo, &
696 : Tsurf_factor, opacity_factor, crystal_core_boundary_mass, &
697 : xmstar, R_center, L_center, v_center, &
698 : cumulative_energy_error, num_retries, ierr)
699 : use utils_lib
700 : character (len=*), intent(in) :: fname
701 : character (len=*), intent(inout) :: net_name
702 : integer, intent(inout) :: species, n_shells, &
703 : year_month_day_when_created, num_retries, model_number
704 : real(dp), intent(inout) :: m_div_msun, initial_z, &
705 : mixing_length_alpha, star_age, tau_factor, Teff, &
706 : power_nuc_burn, power_h_burn, power_he_burn, power_z_burn, power_photo, &
707 : Tsurf_factor, opacity_factor, crystal_core_boundary_mass, &
708 : xmstar, R_center, L_center, v_center, cumulative_energy_error
709 : integer, intent(out) :: ierr
710 : integer :: iounit
711 0 : real(dp) :: initial_y
712 0 : ierr = 0
713 0 : open(newunit=iounit, file=trim(fname), action='read', status='old', iostat=ierr)
714 0 : if (ierr /= 0) then
715 0 : write(*, *) 'failed to open ' // trim(fname)
716 0 : return
717 : end if
718 0 : read(iounit, *, iostat=ierr)
719 0 : if (ierr /= 0) then
720 0 : close(iounit)
721 0 : return
722 : end if
723 0 : read(iounit, *, iostat=ierr)
724 0 : if (ierr /= 0) then
725 0 : close(iounit)
726 0 : return
727 : end if
728 : call read_properties(iounit, &
729 : net_name, species, n_shells, year_month_day_when_created, &
730 : m_div_msun, initial_z, initial_y, mixing_length_alpha, &
731 : model_number, star_age, tau_factor, Teff, &
732 : power_nuc_burn, power_h_burn, power_he_burn, power_z_burn, power_photo, &
733 : Tsurf_factor, opacity_factor, crystal_core_boundary_mass, &
734 : xmstar, R_center, L_center, v_center, &
735 0 : cumulative_energy_error, num_retries, ierr)
736 0 : close(iounit)
737 : end subroutine do_read_saved_model_properties
738 :
739 :
740 0 : subroutine do_read_net_name(iounit, net_name, ierr)
741 : integer, intent(in) :: iounit
742 : character (len=*), intent(inout) :: net_name
743 : integer, intent(out) :: ierr
744 : integer :: species, n_shells, &
745 : year_month_day_when_created, model_number, num_retries
746 0 : real(dp) :: m_div_msun, initial_z, initial_y, &
747 0 : mixing_length_alpha, star_age, tau_factor, Teff, &
748 0 : power_nuc_burn, power_h_burn, power_he_burn, power_z_burn, power_photo, &
749 0 : Tsurf_factor, opacity_factor, crystal_core_boundary_mass, &
750 0 : xmstar, R_center, L_center, v_center, cumulative_energy_error
751 : call read_properties(iounit, &
752 : net_name, species, n_shells, year_month_day_when_created, &
753 : m_div_msun, initial_z, initial_y, mixing_length_alpha, &
754 : model_number, star_age, tau_factor, Teff, &
755 : power_nuc_burn, power_h_burn, power_he_burn, power_z_burn, power_photo, &
756 : Tsurf_factor, opacity_factor, crystal_core_boundary_mass, &
757 : xmstar, R_center, L_center, v_center, &
758 0 : cumulative_energy_error, num_retries, ierr)
759 0 : end subroutine do_read_net_name
760 :
761 :
762 0 : subroutine do_read_saved_model_age(fname, star_age, ierr)
763 : character (len=*), intent(in) :: fname
764 : real(dp), intent(inout) :: star_age
765 : integer, intent(out) :: ierr
766 : character (len=strlen) :: net_name
767 : integer :: species, n_shells, model_number, &
768 : num_retries, year_month_day_when_created
769 0 : real(dp) :: m_div_msun, initial_z, &
770 0 : mixing_length_alpha, cumulative_energy_error, &
771 0 : Teff, tau_factor, Tsurf_factor, &
772 0 : power_nuc_burn, power_h_burn, power_he_burn, power_z_burn, power_photo, &
773 0 : opacity_factor, crystal_core_boundary_mass, &
774 0 : xmstar, R_center, L_center, v_center
775 : call do_read_saved_model_properties(fname, &
776 : net_name, species, n_shells, year_month_day_when_created, &
777 : m_div_msun, initial_z, mixing_length_alpha, &
778 : model_number, star_age, tau_factor, Teff, &
779 : power_nuc_burn, power_h_burn, power_he_burn, power_z_burn, power_photo, &
780 : Tsurf_factor, opacity_factor, crystal_core_boundary_mass, &
781 : xmstar, R_center, L_center, v_center, &
782 0 : cumulative_energy_error, num_retries, ierr)
783 0 : end subroutine do_read_saved_model_age
784 :
785 :
786 18 : subroutine read_properties(iounit, &
787 : net_name, species, n_shells, year_month_day_when_created, &
788 : m_div_msun, initial_z, initial_y, mixing_length_alpha, &
789 : model_number, star_age, tau_factor, Teff, &
790 : power_nuc_burn, power_h_burn, power_he_burn, power_z_burn, power_photo, &
791 : Tsurf_factor, opacity_factor, crystal_core_boundary_mass, &
792 : xmstar, R_center, L_center, v_center, &
793 : cumulative_energy_error, num_retries, ierr)
794 : integer, intent(in) :: iounit
795 : character (len=*), intent(inout) :: net_name
796 : integer, intent(inout) :: species, n_shells, &
797 : year_month_day_when_created, model_number, num_retries
798 : real(dp), intent(inout) :: m_div_msun, initial_z, initial_y, &
799 : mixing_length_alpha, star_age, tau_factor, Teff, &
800 : power_nuc_burn, power_h_burn, power_he_burn, power_z_burn, power_photo, &
801 : Tsurf_factor, opacity_factor, crystal_core_boundary_mass, &
802 : xmstar, R_center, L_center, v_center, cumulative_energy_error
803 : integer, intent(out) :: ierr
804 : character (len=132) :: line
805 18 : real(dp) :: tmp
806 : ierr = 0
807 : do ! until reach a blank line
808 62 : read(iounit, fmt='(a)', iostat=ierr) line
809 62 : if (ierr /= 0) return
810 62 : if (len_trim(line) == 0) return ! blank line
811 44 : if (match_keyword_for_string('net_name', line, net_name)) then; cycle; end if
812 43 : if (match_keyword('species', line, tmp)) then; species = int(tmp); cycle; end if
813 42 : if (match_keyword('n_shells', line, tmp)) then; n_shells = int(tmp); cycle; end if
814 25 : if (match_keyword('model_number', line, tmp)) then; model_number = int(tmp); cycle; end if
815 25 : if (match_keyword('M/Msun', line, m_div_msun)) cycle
816 8 : if (match_keyword('star_age', line, star_age)) cycle
817 8 : if (match_keyword('initial_z', line, initial_z)) cycle
818 7 : if (match_keyword('initial_y', line, initial_y)) cycle
819 6 : if (match_keyword('mixing_length_alpha', line, mixing_length_alpha)) cycle
820 6 : if (match_keyword('tau_factor', line, tau_factor)) cycle
821 6 : if (match_keyword('Teff', line, Teff)) cycle
822 6 : if (match_keyword('power_nuc_burn', line, power_nuc_burn)) cycle
823 6 : if (match_keyword('power_h_burn', line, power_h_burn)) cycle
824 6 : if (match_keyword('power_he_burn', line, power_he_burn)) cycle
825 6 : if (match_keyword('power_z_burn', line, power_z_burn)) cycle
826 6 : if (match_keyword('power_photo', line, power_photo)) cycle
827 6 : if (match_keyword('Tsurf_factor', line, Tsurf_factor)) cycle
828 6 : if (match_keyword('opacity_factor', line, opacity_factor)) cycle
829 6 : if (match_keyword('crystal_core_boundary_mass', line, crystal_core_boundary_mass)) cycle
830 6 : if (match_keyword('xmstar', line, xmstar)) cycle
831 6 : if (match_keyword('R_center', line, R_center)) cycle
832 6 : if (match_keyword('L_center', line, L_center)) cycle
833 6 : if (match_keyword('v_center', line, v_center)) cycle
834 6 : if (match_keyword('cumulative_energy_error', line, cumulative_energy_error)) cycle
835 6 : if (match_keyword('year_month_day_when_created', line, tmp)) then
836 1 : year_month_day_when_created = int(tmp); cycle; end if
837 5 : if (match_keyword('tau_photosphere', line, tmp)) cycle
838 5 : if (match_keyword('num_retries', line, tmp)) then; num_retries = int(tmp); cycle; end if
839 : end do
840 : end subroutine read_properties
841 :
842 :
843 270 : logical function match_keyword(key, txt, value)
844 : ! returns true if leading non-blank part of txt is same as key.
845 : ! i.e., skips leading blanks in txt before testing equality.
846 : character (len=*), intent(in) :: key, txt
847 : real(dp), intent(inout) :: value
848 : integer :: i, j, k, ierr
849 270 : i = len(key)
850 270 : k = len(txt)
851 270 : j = 1
852 6028 : do while (j <= k .and. txt(j:j) == ' ')
853 6028 : j = j+1
854 : end do
855 270 : match_keyword = (txt(j:j+i-1) == key)
856 270 : ierr = 0
857 270 : if (match_keyword) then
858 38 : read(txt(j+i:k), fmt=*, iostat=ierr) value
859 270 : if (ierr /= 0) match_keyword = .false.
860 : end if
861 270 : end function match_keyword
862 :
863 :
864 44 : logical function match_keyword_for_string(key, txt, value)
865 : ! returns true if leading non-blank part of txt is same as key.
866 : ! i.e., skips leading blanks in txt before testing equality.
867 : character (len=*), intent(in) :: key, txt
868 : character (len=*), intent(inout) :: value
869 : integer :: i, j, k, str_len
870 : logical, parameter :: dbg = .false.
871 44 : i = len(key)
872 44 : k = len(txt)
873 44 : j = 1
874 1099 : do while (j <= k .and. txt(j:j) == ' ')
875 1099 : j = j+1
876 : end do
877 44 : match_keyword_for_string = (txt(j:j+i-1) == key)
878 44 : if (.not. match_keyword_for_string) return
879 : if (dbg) then
880 : write(*,*) 'matching ' // trim(key)
881 : write(*,*) 'txt ' // trim(txt)
882 : end if
883 : j = j+i
884 4 : do while (j <= k .and. txt(j:j) == ' ')
885 4 : j = j+1
886 : end do
887 1 : if (j > k) then
888 44 : match_keyword_for_string = .false.
889 : if (dbg) write(*,*) 'j > k'
890 : return
891 : end if
892 1 : if (txt(j:j) /= '''') then
893 44 : match_keyword_for_string = .false.
894 : if (dbg) write(*,*) 'no leading quote'
895 : return
896 : end if
897 1 : j = j+1
898 1 : i = 1
899 1 : str_len = len(value)
900 10 : do while (j <= k .and. txt(j:j) /= '''')
901 9 : value(i:i) = txt(j:j)
902 9 : i = i+1
903 10 : j = j+1
904 : end do
905 248 : do while (i <= str_len)
906 247 : value(i:i) = ' '
907 247 : i = i+1
908 : end do
909 : if (dbg) write(*,*) 'value <' // trim(value) // ">"
910 : end function match_keyword_for_string
911 :
912 :
913 17 : subroutine get_chem_col_names(s, iounit, species, names, perm, ierr)
914 : use chem_def, only: iso_name_length
915 : use chem_lib, only: chem_get_iso_id
916 : type (star_info), pointer :: s
917 : integer, intent(in) :: iounit, species
918 : character(len=iso_name_length), intent(out) :: names(species)
919 : integer, intent(out) :: perm(species)
920 : integer, intent(out) :: ierr
921 :
922 : character (len=50000) :: buffer
923 : character (len=20) :: string
924 : integer :: n, i, j1, j2, str_len, l, indx, j, num_found
925 :
926 17 : ierr = 0
927 17 : read(iounit,fmt='(a)',iostat=ierr) buffer
928 17 : if (ierr /= 0) return
929 :
930 17 : n = len_trim(buffer)
931 17 : i = 0
932 17 : num_found = 0
933 : token_loop: do ! have non-empty buffer
934 4964 : i = i+1
935 4964 : if (i > n) then
936 0 : write(*,*) 'get_chem_col_names: failed to find all of the names'
937 0 : ierr = -1
938 0 : return
939 : end if
940 4964 : if (buffer(i:i) == char(9)) cycle token_loop ! skip tabs
941 : select case(buffer(i:i))
942 : case (' ')
943 : cycle token_loop ! skip spaces
944 : case default
945 : j1 = i; j2 = i
946 : name_loop: do
947 629 : if (i+1 > n) exit name_loop
948 612 : if (buffer(i+1:i+1) == ' ') exit name_loop
949 408 : if (buffer(i+1:i+1) == '(') exit name_loop
950 : if (buffer(i+1:i+1) == ')') exit name_loop
951 : if (buffer(i+1:i+1) == ',') exit name_loop
952 : i = i+1
953 221 : j2 = i
954 : end do name_loop
955 221 : str_len = len(string)
956 221 : l = j2-j1+1
957 221 : if (l > str_len) then
958 0 : l = str_len
959 0 : j2 = l+j1-1
960 : end if
961 221 : string(1:l) = buffer(j1:j2)
962 4012 : do j = l+1, str_len
963 4012 : string(j:j) = ' '
964 : end do
965 :
966 221 : indx = chem_get_iso_id(string)
967 :
968 5185 : if (indx > 0) then
969 136 : num_found = num_found+1
970 136 : names(num_found) = trim(string)
971 136 : perm(num_found) = s% net_iso(indx)
972 : !write(*,*) trim(string), num_found, perm(num_found)
973 136 : if (num_found == species) return
974 : end if
975 :
976 : end select
977 : end do token_loop
978 :
979 17 : end subroutine get_chem_col_names
980 :
981 : end module read_model
|
