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