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 init_model
21 :
22 : use star_private_def
23 : use const_def, only: dp, msun, secyer
24 :
25 : implicit none
26 :
27 : private
28 : public :: get_zams_model
29 :
30 : integer, parameter :: min_when_created = 20081212
31 :
32 : contains
33 :
34 : subroutine get_revised_mass(s, fullname, ierr)
35 : use read_model, only: do_read_saved_model
36 : use relax, only:do_relax_mass
37 : use relax, only: do_relax_Z
38 : type (star_info), pointer :: s
39 : character (len=*), intent(in) :: fullname
40 : integer, intent(out) :: ierr
41 :
42 : real(dp), parameter :: lg_max_abs_mdot = -3.5d0
43 : real(dp) :: init_mass, init_z
44 :
45 : init_mass = s% initial_mass
46 : init_z = s% initial_z
47 :
48 : call do_read_saved_model(s, fullname, ierr)
49 : if (ierr /= 0) return
50 :
51 : if (abs(s% initial_z - init_z) > 1d-3*init_z) then
52 : ierr = -1
53 : return
54 : end if
55 :
56 : call do_relax_mass(s% id, init_mass, lg_max_abs_mdot, ierr)
57 : if (ierr /= 0) return
58 :
59 : s% initial_mass = init_mass
60 : s% dt_next = min(s% dt_next, 1d2*secyer)
61 :
62 : end subroutine get_revised_mass
63 :
64 :
65 1 : subroutine get_zams_model(s, zams_filename, ierr)
66 : use alloc, only: allocate_star_info_arrays
67 : use utils_lib, only: is_bad
68 : use relax, only: do_relax_mass
69 : type (star_info), pointer :: s
70 : character (len=*) :: zams_filename
71 : integer, intent(out) :: ierr
72 :
73 : integer :: nz
74 1 : real(dp), dimension(:,:), pointer :: xh, xa
75 1 : real(dp), dimension(:), pointer :: q, dq, omega
76 1 : real(dp) :: init_mass
77 : logical :: in_range
78 : real(dp), parameter :: lg_max_abs_mdot = -3.5d0
79 :
80 : include 'formats'
81 :
82 1 : ierr = 0
83 1 : if (is_bad(s% initial_mass)) then
84 0 : write(*,1) 's% initial_mass', s% initial_mass
85 0 : call mesa_error(__FILE__,__LINE__,'get_zams_model')
86 : end if
87 :
88 1 : init_mass = s% initial_mass
89 1 : s% mstar = s% initial_mass*Msun
90 1 : s% xmstar = s% mstar
91 1 : s% M_center = 0
92 :
93 : call get1_zams_model( &
94 : s, zams_filename, nz, xh, xa, q, dq, &
95 1 : omega, in_range, ierr)
96 1 : if (ierr /= 0) then
97 0 : write(*,1) 'failed in get1_zams_model'
98 0 : call mesa_error(__FILE__,__LINE__,'get_zams_model')
99 : end if
100 :
101 :
102 1 : s% nz = nz
103 1 : call allocate_star_info_arrays(s, ierr)
104 1 : if (ierr /= 0) then
105 0 : call dealloc
106 0 : return
107 : end if
108 :
109 : ! copy, then deallocate
110 24641 : s% xh(:,1:nz) = xh(:,1:nz)
111 44353 : s% xa(:,1:nz) = xa(:,1:nz)
112 4929 : s% q(1:nz) = q(1:nz)
113 4929 : s% dq(1:nz) = dq(1:nz)
114 4929 : s% omega(1:nz) = omega(1:nz)
115 :
116 1 : call dealloc
117 :
118 1 : if (.not. in_range) then ! have revised s% initial_mass
119 0 : s% mstar = s% initial_mass*Msun
120 0 : s% xmstar = s% mstar
121 0 : s% M_center = 0
122 0 : s% dt_next = 1d2*secyer
123 0 : call do_relax_mass(s% id, init_mass, lg_max_abs_mdot, ierr)
124 0 : if (ierr /= 0) return
125 0 : s% initial_mass = init_mass
126 0 : s% dt_next = min(s% dt_next, 1d2*secyer)
127 : end if
128 :
129 : contains
130 :
131 1 : subroutine dealloc
132 1 : deallocate(xh, xa, q, dq, omega)
133 1 : end subroutine dealloc
134 :
135 : end subroutine get_zams_model
136 :
137 :
138 1 : subroutine get1_zams_model( &
139 : s, zams_filename, nz, xh, xa, q, dq, &
140 : omega, in_range, ierr)
141 : use utils_lib
142 : use const_def, only: mesa_data_dir
143 : use net, only: set_net
144 : type (star_info), pointer :: s
145 : character (len=*), intent(in) :: zams_filename
146 : integer, intent(out) :: nz
147 : real(dp), dimension(:,:), pointer :: xh, xa
148 1 : real(dp), dimension(:), pointer :: q, dq, omega, j_rot
149 : logical, intent(out) :: in_range
150 : integer, intent(out) :: ierr
151 :
152 : integer :: iounit, nz1, nz2, file_type, nvar_hydro, species
153 : character (len=250) :: fname, line
154 1 : real(dp) :: m1, m2, initial_mass
155 : logical :: okay
156 :
157 : include 'formats'
158 :
159 1 : ierr = 0
160 1 : nz = 0
161 :
162 1 : fname = zams_filename
163 1 : open(newunit=iounit, file=trim(fname), action='read', status='old', iostat=ierr)
164 1 : if (ierr /= 0) then
165 1 : ierr = 0
166 1 : fname = trim(mesa_data_dir) // '/star_data/zams_models/' // trim(zams_filename)
167 1 : open(newunit=iounit, file=trim(fname), action='read', status='old', iostat=ierr)
168 1 : if (ierr /= 0) then
169 0 : write(*, *) 'failed to open ' // trim(zams_filename)
170 0 : write(*, *) 'failed to open ' // trim(fname)
171 0 : return
172 : end if
173 : end if
174 :
175 1 : read(iounit, *, iostat=ierr) file_type
176 1 : if (ierr /= 0) then
177 0 : write(*, *) 'ERROR: first line needs to contains file type number: ' // trim(fname)
178 0 : close(iounit)
179 0 : return
180 : end if
181 :
182 1 : initial_mass = s% initial_mass
183 1 : nvar_hydro = s% nvar_hydro
184 :
185 1 : call read_zams_header ! sets net_name
186 1 : if (ierr /= 0) then
187 0 : close(iounit)
188 0 : write(*,*) 'failed in read_zams_header'
189 0 : call mesa_error(__FILE__,__LINE__,'get1_zams_model')
190 0 : return
191 : end if
192 :
193 1 : call set_net(s, s% net_name, ierr)
194 1 : if (ierr /= 0) then
195 0 : write(*,*) 'failed in set_net'
196 0 : return
197 : end if
198 :
199 1 : species = s% species
200 :
201 : ! find mass that is closest to init_m
202 1 : m1 = 0
203 1 : nz1 = 0
204 1 : okay = .false.
205 1 : in_range = .true.
206 :
207 1 : read(iounit, *, iostat=ierr) ! this is the 'M/Msun n_shells' header line
208 1 : if (ierr == 0) then
209 16 : index_loop: do
210 17 : read(iounit, *, iostat=ierr) m2, nz2
211 17 : if (ierr /= 0) exit index_loop
212 17 : if (m2 <= 0) then ! end of list
213 0 : read(iounit, *, iostat=ierr) ! blank line
214 0 : m2 = m1
215 0 : nz2 = nz1
216 0 : in_range = .false.
217 0 : s% initial_mass = m2
218 0 : initial_mass = m2
219 0 : exit index_loop
220 : end if
221 17 : if (m2 >= initial_mass) then
222 1 : if (m1 == 0) then
223 0 : m1 = m2
224 0 : nz1 = nz2
225 0 : in_range = .false.
226 0 : s% initial_mass = m2
227 0 : initial_mass = m2
228 : end if
229 : ! skip to end of index
230 : skip_loop: do
231 20 : read(iounit, fmt='(a)', iostat=ierr) line
232 20 : if (len_trim(line) == 0) then ! blank line indicates end of index
233 1 : okay = .true.
234 : exit index_loop
235 : end if
236 19 : if (ierr /= 0) exit index_loop
237 : end do skip_loop
238 : exit index_loop
239 : end if
240 16 : m1 = m2
241 16 : nz1 = nz2
242 : end do index_loop
243 : end if
244 :
245 1 : nz = nz1
246 :
247 : allocate(xh(nvar_hydro,nz), xa(species,nz), q(nz), dq(nz), &
248 3 : omega(nz), j_rot(nz), stat=ierr)
249 1 : if (ierr /= 0) then
250 0 : close(iounit)
251 0 : return
252 : end if
253 :
254 : call get1_mass( &
255 : s, iounit, m1, nz1, m2, nz2, initial_mass, &
256 : nvar_hydro, species, xh, xa, q, dq, &
257 1 : omega, j_rot, ierr)
258 1 : if (ierr /= 0) then
259 0 : write(*,*) 'failed in get1_mass'
260 0 : call mesa_error(__FILE__,__LINE__,'get_zams_model')
261 : end if
262 1 : close(iounit)
263 4 : deallocate(j_rot)
264 :
265 : contains
266 :
267 1 : subroutine read_zams_header
268 1 : use read_model, only: read_properties
269 : integer :: year_month_day_when_created, iprop
270 1 : real(dp) :: dprop, initial_z, initial_y
271 1 : read(iounit, *, iostat=ierr) ! skip blank line before property list
272 : include 'formats'
273 1 : if (ierr /= 0) return
274 1 : year_month_day_when_created = -1
275 :
276 :
277 : call read_properties(iounit, &
278 : s% net_name, iprop, iprop, year_month_day_when_created, &
279 : dprop, initial_z, initial_y, &
280 : dprop, iprop, dprop, dprop, &
281 : dprop, dprop, dprop, dprop, &
282 : dprop, dprop, dprop, dprop, dprop, dprop, &
283 1 : dprop, dprop, dprop, dprop, iprop, ierr)
284 1 : if (ierr /= 0) then
285 0 : write(*,2) 'year_month_day_when_created', year_month_day_when_created
286 0 : write(*,*) 'net_name' // trim(s% net_name)
287 0 : call mesa_error(__FILE__,__LINE__,'read_zams_header')
288 0 : return
289 : end if
290 :
291 1 : if (year_month_day_when_created < min_when_created) then
292 0 : ierr = -1
293 0 : write(*, *)
294 0 : write(*, *)
295 0 : write(*, *)
296 0 : write(*, *) 'sorry: you need to update the zams data file.'
297 0 : write(*, *) 'found "year_month_day_when_created" =', year_month_day_when_created
298 0 : write(*, *) 'but need at least', min_when_created
299 0 : write(*, *)
300 0 : write(*, *)
301 0 : write(*, *)
302 0 : return
303 : end if
304 1 : if (abs(initial_z - s% initial_z) > 1d-3*s% initial_z) then
305 0 : ierr = -1
306 0 : write(*, *)
307 0 : write(*, *)
308 0 : write(*, *)
309 0 : write(*, *) 'WARNING: requested initial_z does not match zams file initial_z.'
310 0 : write(*, 1) 'zams file initial_z', initial_z
311 0 : write(*, 1) 'requested initial_z', s% initial_z
312 0 : write(*, *)
313 0 : write(*, *)
314 0 : write(*, *)
315 0 : return
316 : end if
317 1 : if (s% initial_y > 0 .and. &
318 : abs(initial_y - s% initial_y) > 1d-3*s% initial_y) then
319 0 : ierr = -1
320 0 : write(*, *)
321 0 : write(*, *)
322 0 : write(*, *)
323 0 : write(*, *) 'WARNING: requested initial_y does not match zams file initial_y.'
324 0 : write(*, 1) 'zams file initial_y', initial_y
325 0 : write(*, 1) 'requested initial_y', s% initial_y
326 0 : write(*, *)
327 0 : write(*, *)
328 0 : write(*, *)
329 0 : return
330 : end if
331 1 : end subroutine read_zams_header
332 :
333 : end subroutine get1_zams_model
334 :
335 :
336 1 : subroutine get1_mass( &
337 : s, iounit, m1, nz1, m2, nz2, initial_mass, &
338 1 : nvar_hydro, species, xh, xa, q, dq, &
339 1 : omega, j_rot, ierr)
340 : use read_model, only: read_properties, read1_model
341 : use chem_def, only: iso_name_length
342 : use read_model, only: get_chem_col_names
343 : use star_utils, only: interp_q
344 : type (star_info), pointer :: s
345 : integer, intent(in) :: iounit, nz1, nz2, nvar_hydro, species
346 : real(dp), intent(in) :: m1, m2, initial_mass
347 : real(dp), intent(inout) :: xh(:,:) ! (nvar_hydro,nz1)
348 : real(dp), intent(inout) :: xa(:,:) ! (species,nz1)
349 : real(dp), intent(inout), dimension(:) :: &
350 : q, dq, omega, j_rot ! (nz1)
351 : integer, intent(out) :: ierr
352 :
353 : integer :: i, k, nz, nz_in, iprop
354 1 : real(dp) :: m_in, m_read, dprop, lnm1, lnm2
355 1 : real(dp), dimension(:, :), pointer :: xh2, xa2
356 : real(dp), dimension(:), pointer :: &
357 1 : q2, dq2, omega2, j_rot2
358 14 : real(dp) :: alfa, struct(nvar_hydro), comp(species)
359 : logical :: okay
360 : character (len=net_name_len) :: net_name
361 1 : character(len=iso_name_length), pointer :: names(:) ! (species)
362 1 : integer, pointer :: perm(:) ! (species)
363 :
364 : include 'formats'
365 :
366 1 : nz = nz1
367 1 : m_read = m1
368 :
369 : allocate( &
370 : xh2(nvar_hydro, nz2), xa2(species, nz2), q2(nz2), dq2(nz2), &
371 : omega2(nz2), j_rot2(nz2), &
372 3 : names(species), perm(species), stat=ierr)
373 1 : if (ierr /= 0) return
374 : okay = .false.
375 17 : mass_loop: do ! loop until find desired mass
376 :
377 17 : m_in = -1; nz_in = -1; net_name = ''
378 : call read_properties(iounit, &
379 : net_name, iprop, nz_in, iprop, m_in, &
380 : dprop, dprop, dprop, iprop, &
381 : dprop, dprop, dprop, dprop, dprop, &
382 : dprop, dprop, dprop, dprop, dprop, dprop, &
383 17 : dprop, dprop, dprop, dprop, dprop, iprop, ierr)
384 17 : if (ierr /= 0 .or. m_in < 0 .or. nz_in < 0) then
385 0 : write(*,*) 'missing required properties'
386 0 : write(*,*) 'ierr', ierr
387 0 : write(*,*) 'm_in', m_in
388 0 : write(*,*) 'nz_in', nz_in
389 0 : ierr = -1
390 0 : exit mass_loop
391 : end if
392 :
393 17 : call get_chem_col_names(s, iounit, species, names, perm, ierr)
394 17 : if (ierr /= 0) exit mass_loop
395 :
396 17 : if (abs(m_in-m_read) > 1d-4) then
397 :
398 50344 : do i = 1, nz_in ! skip this one
399 50329 : read(iounit, *, iostat=ierr)
400 50344 : if (ierr /= 0) exit mass_loop
401 : end do
402 :
403 : else ! store this one
404 :
405 2 : if (nz /= nz_in) then
406 : write(*, '(a, 2i6)') &
407 0 : 'nz /= nz_in: logic error in routine for reading model file', nz, nz_in
408 0 : ierr = -1
409 0 : return
410 : end if
411 :
412 2 : if (m_read == m1) then
413 : call read1_model( &
414 : s, species, nvar_hydro, nz, iounit, &
415 1 : xh, xa, q, dq, omega, j_rot, perm, ierr)
416 1 : if (ierr /= 0) exit mass_loop
417 1 : okay = .true.
418 1 : if (m2 == m1) exit mass_loop
419 1 : m_read = m2
420 1 : nz = nz2
421 : else
422 : call read1_model( &
423 : s, species, nvar_hydro, nz, iounit, &
424 1 : xh2, xa2, q2, dq2, omega2, j_rot2, perm, ierr)
425 1 : if (ierr /= 0) exit mass_loop
426 1 : okay = .true.
427 1 : nz = nz1
428 : exit mass_loop
429 : end if
430 :
431 : end if
432 16 : read(iounit, *, iostat=ierr) ! skip line following the last zone
433 16 : if (ierr /= 0) exit mass_loop
434 :
435 : end do mass_loop
436 :
437 0 : if (.not. okay) then
438 0 : ierr = -1
439 0 : call dealloc
440 0 : return
441 : end if
442 :
443 2 : if (m1 /= m2) then ! interpolate linearly in log(m)
444 1 : lnm1 = log(m1)
445 1 : lnm2 = log(m2)
446 1 : alfa = (log(initial_mass) - lnm2) / (lnm1 - lnm2)
447 2465 : do k=1,nz1
448 : call interp_q( &
449 2464 : nz2, nvar_hydro, species, q(k), xh2, xa2, q2, dq2, struct, comp, ierr)
450 2464 : if (ierr /= 0) then
451 0 : call dealloc
452 0 : return
453 : end if
454 12320 : xh(1:nvar_hydro,k) = alfa*xh(1:nvar_hydro,k) + (1-alfa)*struct(1:nvar_hydro)
455 22177 : xa(1:species,k) = alfa*xa(1:species,k) + (1-alfa)*comp(1:species)
456 : end do
457 1 : call dealloc
458 : end if
459 :
460 : contains
461 :
462 1 : subroutine dealloc
463 1 : deallocate(xh2, xa2, q2, dq2, omega2, j_rot2, names, perm)
464 1 : end subroutine dealloc
465 :
466 : end subroutine get1_mass
467 :
468 : end module init_model
|
