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