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