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