Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2018 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 pulse_gyre
21 :
22 : use star_private_def
23 : use const_def, only: dp, pi, four_thirds, rsun
24 : use utils_lib
25 : use atm_def
26 : use atm_support
27 : use eps_grav
28 : use pulse_utils
29 :
30 : implicit none
31 :
32 : private
33 : public :: get_gyre_data
34 : public :: write_gyre_data
35 :
36 : contains
37 :
38 0 : subroutine get_gyre_data (id, &
39 : add_center_point, keep_surface_point, add_atmosphere, global_data, point_data, ierr)
40 :
41 : integer, intent(in) :: id
42 : logical, intent(in) :: add_center_point
43 : logical, intent(in) :: keep_surface_point
44 : logical, intent(in) :: add_atmosphere
45 : real(dp), allocatable, intent(out) :: global_data(:)
46 : real(dp), allocatable, target, intent(out) :: point_data(:,:)
47 : integer, intent(out) :: ierr
48 :
49 : type(star_info), pointer :: s
50 0 : integer, allocatable :: k_a(:)
51 0 : integer, allocatable :: k_b(:)
52 : integer :: n_sg
53 : integer :: n_atm
54 : integer :: nn_atm
55 : integer :: n_env
56 : integer :: nn_env
57 : integer :: nn
58 : real(dp), pointer :: r(:) => NULL()
59 : real(dp), pointer :: m(:) => NULL()
60 : real(dp), pointer :: L(:) => NULL()
61 : real(dp), pointer :: P(:) => NULL()
62 : real(dp), pointer :: T(:) => NULL()
63 : real(dp), pointer :: rho(:) => NULL()
64 : real(dp), pointer :: nabla(:) => NULL()
65 : real(dp), pointer :: N2(:) => NULL()
66 : real(dp), pointer :: Gamma_1(:) => NULL()
67 : real(dp), pointer :: nabla_ad(:) => NULL()
68 : real(dp), pointer :: delta(:) => NULL()
69 : real(dp), pointer :: kap(:) => NULL()
70 : real(dp), pointer :: kap_kap_T(:) => NULL()
71 : real(dp), pointer :: kap_kap_rho(:) => NULL()
72 : real(dp), pointer :: eps_nuc(:) => NULL()
73 : real(dp), pointer :: eps_eps_T(:) => NULL()
74 : real(dp), pointer :: eps_eps_rho(:) => NULL()
75 : real(dp), pointer :: eps_grav(:) => NULL()
76 : real(dp), pointer :: Omega_rot(:) => NULL()
77 0 : real(dp) :: r_outer
78 0 : real(dp) :: m_outer
79 : integer :: j
80 : integer :: k
81 : integer :: sg
82 :
83 : ! Get model data for GYRE output
84 :
85 0 : call get_star_ptr(id, s, ierr)
86 0 : if (ierr /= 0) then
87 0 : write(*,*) 'bad star id for get_gyre_data'
88 0 : return
89 : end if
90 :
91 : ! Set up segment indices
92 :
93 : call set_segment_indices(s, k_a, k_b, add_center_point)
94 :
95 0 : n_sg = SIZE(k_a)
96 :
97 : ! Determine data dimensions
98 :
99 0 : if (add_atmosphere) then
100 0 : call build_atm(s, s%L(1), s%r(1), s%Teff, s%m_grav(1), s%cgrav(1), ierr)
101 0 : if (ierr /= 0) then
102 0 : write(*,*) 'failed in build_atm'
103 0 : return
104 : end if
105 0 : n_atm = s%atm_structure_num_pts
106 0 : nn_atm = n_atm - 1
107 : else
108 : n_atm = 0
109 : nn_atm = 0
110 : end if
111 :
112 0 : n_env = s%nz
113 :
114 0 : if (keep_surface_point) then
115 0 : nn_env = n_env + n_sg - 1
116 : else
117 0 : nn_env = n_env - 1 + n_sg - 1
118 : end if
119 :
120 0 : if (add_center_point) then
121 0 : nn = nn_env + nn_atm + 1
122 : else
123 0 : nn = nn_env + nn_atm
124 : end if
125 :
126 : ! Allocate arrays & set up data pointers
127 :
128 0 : allocate(global_data(3))
129 :
130 : ! Set up data pointers
131 :
132 0 : select case(s%gyre_data_schema)
133 :
134 : case(101,110)
135 :
136 0 : allocate(point_data(18,nn))
137 :
138 0 : r => point_data(1,:)
139 0 : m => point_data(2,:)
140 0 : L => point_data(3,:)
141 0 : P => point_data(4,:)
142 0 : T => point_data(5,:)
143 0 : rho => point_data(6,:)
144 0 : nabla => point_data(7,:)
145 0 : N2 => point_data(8,:)
146 0 : Gamma_1 => point_data(9,:)
147 0 : nabla_ad => point_data(10,:)
148 0 : delta => point_data(11,:)
149 0 : kap => point_data(12,:)
150 0 : kap_kap_T => point_data(13,:)
151 0 : kap_kap_rho => point_data(14,:)
152 0 : eps_nuc => point_data(15,:)
153 0 : eps_eps_T => point_data(16,:)
154 0 : eps_eps_rho => point_data(17,:)
155 0 : Omega_rot => point_data(18,:)
156 :
157 : case(120)
158 :
159 0 : allocate(point_data(19,nn))
160 :
161 0 : r => point_data(1,:)
162 0 : m => point_data(2,:)
163 0 : L => point_data(3,:)
164 0 : P => point_data(4,:)
165 0 : T => point_data(5,:)
166 0 : rho => point_data(6,:)
167 0 : nabla => point_data(7,:)
168 0 : N2 => point_data(8,:)
169 0 : Gamma_1 => point_data(9,:)
170 0 : nabla_ad => point_data(10,:)
171 0 : delta => point_data(11,:)
172 0 : kap => point_data(12,:)
173 0 : kap_kap_T => point_data(13,:)
174 0 : kap_kap_rho => point_data(14,:)
175 0 : eps_nuc => point_data(15,:)
176 0 : eps_eps_T => point_data(16,:)
177 0 : eps_eps_rho => point_data(17,:)
178 0 : eps_grav => point_data(18,:)
179 0 : Omega_rot => point_data(19,:)
180 :
181 : case default
182 :
183 0 : write(*,*) 'invalid gyre_data_schema'
184 0 : ierr = -1
185 0 : return
186 :
187 : end select
188 :
189 : ! If necessary, update the eps_grav data in the model
190 :
191 0 : if (ASSOCIATED(eps_grav)) then
192 :
193 0 : do k = 1, s%nz
194 0 : call eval_eps_grav_and_partials(s, k, ierr)
195 0 : if (ierr /= 0) then
196 0 : write(*,*) 'failed in call to eval_eps_grav_and_partials'
197 0 : return
198 : end if
199 : end do
200 :
201 : end if
202 :
203 : ! Store global data
204 :
205 0 : r_outer = Rsun*s%photosphere_r
206 0 : m_outer = s%m_grav(1)
207 :
208 0 : global_data(1) = m_outer
209 0 : global_data(2) = r_outer
210 0 : global_data(3) = s%L(1)
211 :
212 : ! Store point data
213 :
214 0 : j = 1
215 :
216 : ! Atmosphere (we skip the point at the base of the atm to smooth
217 : ! the transition)
218 :
219 0 : atm_loop : do k = 1, n_atm-1
220 0 : call store_point_data_atm(j, k)
221 0 : j = j + 1
222 : end do atm_loop
223 :
224 : ! Envelope
225 :
226 0 : sg = 1
227 :
228 0 : env_loop : do k = 1, n_env
229 :
230 0 : if (k == 1 .AND. .NOT. keep_surface_point) cycle env_loop
231 :
232 0 : call store_point_data_env(j, k, k_a(sg), k_b(sg))
233 0 : j = j + 1
234 :
235 0 : if (k == k_b(sg)+1) then
236 :
237 0 : sg = sg + 1
238 :
239 0 : call store_point_data_env(j, k, k_a(sg), k_b(sg))
240 0 : j = j + 1
241 :
242 : end if
243 :
244 : end do env_loop
245 :
246 : ! Center
247 :
248 0 : if (add_center_point) then
249 0 : call store_point_data_ctr(j, k_a(n_sg), k_b(n_sg))
250 0 : j = j + 1
251 : end if
252 :
253 : ! Check that all point data has correctly been calculated
254 :
255 0 : if (j /= nn+1) call mesa_error(__FILE__,__LINE__,'Invalid cell index in get_gyre_data')
256 :
257 : ! Reverse the data ordering (GYRE format is center-to-surface)
258 :
259 0 : point_data = point_data(:,nn:1:-1)
260 :
261 : ! Tidy up
262 :
263 0 : if (ASSOCIATED(s%atm_structure)) then
264 0 : deallocate(s%atm_structure)
265 : end if
266 :
267 0 : return
268 :
269 : contains
270 :
271 0 : subroutine store_point_data_atm (j, k)
272 :
273 : integer, intent(in) :: j
274 : integer, intent(in) :: k
275 :
276 0 : real(dp) :: grav
277 :
278 : ! Store data associated with atmosphere point k into the
279 : ! point_data array at position j
280 :
281 0 : r(j) = s%r(1) + s%atm_structure(atm_delta_r,k)
282 0 : m(j) = s%m_grav(1) !+ s%atm_structure(atm_delta_m,k)
283 0 : L(j) = s%L(1)
284 :
285 0 : P(j) = exp(s%atm_structure(atm_lnP,k))
286 0 : rho(j) = exp(s%atm_structure(atm_lnd,k))
287 0 : T(j) = exp(s%atm_structure(atm_lnT,k))
288 :
289 0 : Gamma_1(j) = s%atm_structure(atm_gamma1,k)
290 0 : nabla_ad(j) = s%atm_structure(atm_grada,k)
291 0 : delta(j) = s%atm_structure(atm_chiT,k)/s%atm_structure(atm_chiRho,k)
292 0 : nabla(j) = s%atm_structure(atm_gradT,k)
293 :
294 0 : grav = s%cgrav(1)*m(j)/(r(j)*r(j))
295 0 : N2(j) = grav*grav*(rho(j)/P(j))*delta(j)*(nabla_ad(j) - nabla(j))
296 :
297 0 : kap(j) = s%atm_structure(atm_kap,k)
298 0 : kap_kap_T(j) = kap(j)*s%atm_structure(atm_dlnkap_dlnT,k)
299 0 : kap_kap_rho(j) = kap(j)*s%atm_structure(atm_dlnkap_dlnd,k)
300 :
301 0 : eps_nuc(j) = 0d0
302 0 : eps_eps_T(j) = 0d0
303 0 : eps_eps_rho(j) = 0d0
304 0 : if (ASSOCIATED(eps_grav)) eps_grav(j) = 0d0
305 :
306 0 : if (s%rotation_flag) then
307 0 : Omega_rot(j) = s%omega(1)
308 : else
309 0 : Omega_rot(j) = 0d0
310 : end if
311 :
312 0 : return
313 :
314 : end subroutine store_point_data_atm
315 :
316 :
317 0 : subroutine store_point_data_env (j, k, k_a, k_b)
318 :
319 : integer, intent(in) :: j
320 : integer, intent(in) :: k
321 : integer, intent(in) :: k_a
322 : integer, intent(in) :: k_b
323 :
324 : ! Store data associated with envelope face k into the point_data
325 : ! array at position j
326 :
327 0 : r(j) = s%r(k)
328 0 : m(j) = s%m_grav(k)
329 0 : L(j) = s%L(k)
330 :
331 0 : P(j) = eval_face(s%dq, s%Peos, k, 1, s%nz)
332 0 : if (s%interpolate_rho_for_pulse_data) then
333 0 : rho(j) = eval_face(s%dq, s%rho, k, k_a, k_b)
334 : else
335 0 : rho(j) = eval_face_rho(s, k, k_a, k_b)
336 : end if
337 0 : T(j) = eval_face(s%dq, s%T, k, 1, s%nz)
338 :
339 0 : N2(j) = eval_face_A_ast(s, k, k_a, k_b)*s%grav(k)/s%r(k)
340 0 : Gamma_1(j) = eval_face(s%dq, s%gamma1, k, k_a, k_b)
341 0 : nabla_ad(j) = eval_face(s%dq, s%grada, k, k_a, k_b)
342 0 : delta(j) = eval_face(s%dq, s%chiT, k, k_a, k_b)/eval_face(s%dq, s%chiRho, k, k_a, k_b)
343 0 : nabla(j) = s%gradT(k) ! Not quite right; gradT can be discontinuous
344 :
345 0 : kap(j) = eval_face(s%dq, s%opacity, k, k_a, k_b)
346 0 : kap_kap_T(j) = eval_face(s%dq, s%d_opacity_dlnT, k, k_a, k_b)
347 0 : kap_kap_rho(j) = eval_face(s%dq, s%d_opacity_dlnd, k, k_a, k_b)
348 :
349 0 : eps_nuc(j) = eval_face(s%dq, s%eps_nuc, k, k_a, k_b)
350 0 : eps_eps_T(j) = eval_face(s%dq, s%d_epsnuc_dlnT, k, k_a, k_b)
351 0 : eps_eps_rho(j) = eval_face(s%dq, s%d_epsnuc_dlnd, k, k_a, k_b)
352 0 : if (ASSOCIATED(eps_grav)) eps_grav(j) = eval_face(s%dq, s%eps_grav_ad%val, k, k_a, k_b)
353 :
354 0 : if (s%rotation_flag) then
355 0 : Omega_rot(j) = s%omega(k) ! Not quite right; omega can be discontinuous
356 : else
357 0 : Omega_rot = 0d0
358 : end if
359 :
360 0 : return
361 :
362 : end subroutine store_point_data_env
363 :
364 :
365 0 : subroutine store_point_data_ctr (j, k_a, k_b)
366 :
367 : integer, intent(in) :: j
368 : integer, intent(in) :: k_a
369 : integer, intent(in) :: k_b
370 :
371 0 : real(dp) :: d2P_dr2_c
372 :
373 : ! Store data for the center into the point_data array at position j
374 :
375 0 : r(j) = 0d0
376 0 : m(j) = 0d0
377 0 : L(j) = 0d0
378 :
379 0 : P(j) = eval_center(s%rmid, s%Peos, 1, s%nz)
380 0 : if (s%interpolate_rho_for_pulse_data) then
381 0 : rho(j) = eval_center(s%rmid, s%rho, k_a, k_b)
382 : else
383 0 : rho(j) = eval_center_rho(s, k_b)
384 : end if
385 : ! at the centre d²P/dr² = -4πGρ²/3
386 0 : d2P_dr2_c = -four_thirds*pi*s% cgrav(s% nz)*rho(j)**2
387 0 : P(j) = s%Peos(s% nz) - 0.5d0*d2P_dr2_c*s% rmid(s% nz)**2
388 0 : T(j) = eval_center(s%rmid, s%T, 1, s%nz)
389 :
390 0 : N2(j) = 0d0
391 0 : Gamma_1(j) = eval_center(s%rmid, s%gamma1, k_a, k_b)
392 0 : nabla_ad(j) = eval_center(s%rmid, s%grada, k_a, k_b)
393 0 : delta(j) = eval_center(s%rmid, s%chiT, k_a, k_b)/eval_center(s%rmid, s%chiRho, k_a, k_b)
394 0 : nabla(j) = eval_center(s%r, s%gradT, k_a, k_b)
395 :
396 0 : kap(j) = eval_center(s%rmid, s%opacity, k_a, k_b)
397 0 : kap_kap_T(j) = eval_center(s%rmid, s%d_opacity_dlnT, k_a, k_b)
398 0 : kap_kap_rho(j) = eval_center(s%rmid, s%d_opacity_dlnd, k_a, k_b)
399 :
400 0 : eps_nuc(j) = eval_center(s%rmid, s%eps_nuc, k_a, k_b)
401 0 : eps_eps_T(j) = eval_center(s%rmid, s%d_epsnuc_dlnT, k_a, k_b)
402 0 : eps_eps_rho(j) = eval_center(s%rmid, s%d_epsnuc_dlnd, k_a, k_b)
403 0 : if (ASSOCIATED(eps_grav)) eps_grav(j) = eval_center(s%rmid, s%eps_grav_ad%val, k_a, k_b)
404 :
405 0 : if (s%rotation_flag) then
406 0 : Omega_rot(j) = eval_center(s%r, s%omega, k_a, k_b)
407 : else
408 0 : Omega_rot(j) = 0d0
409 : end if
410 :
411 :
412 0 : return
413 :
414 : end subroutine store_point_data_ctr
415 :
416 : end subroutine get_gyre_data
417 :
418 :
419 0 : subroutine write_gyre_data (id, filename, global_data, point_data, ierr)
420 :
421 : integer, intent(in) :: id
422 : character(*), intent(in) :: filename
423 : real(dp), intent(in) :: global_data(:)
424 : real(dp), intent(in) :: point_data(:,:)
425 : integer, intent(out) :: ierr
426 :
427 : type(star_info), pointer :: s
428 : integer :: iounit
429 : integer :: nn
430 : integer :: j
431 :
432 : ! Write GYRE data to file
433 :
434 0 : call get_star_ptr(id, s, ierr)
435 0 : if (ierr /= 0) then
436 0 : write(*,*) 'bad star id for write_gyre_data'
437 0 : return
438 : end if
439 :
440 0 : select case(s%gyre_data_schema)
441 : case(101,120)
442 : case default
443 0 : write(*,*) 'invalid gyre_data_schema'
444 0 : ierr = -1
445 0 : return
446 : end select
447 :
448 : ! Open the file
449 :
450 0 : open(newunit=iounit, file=TRIM(filename), status='REPLACE', iostat=ierr)
451 0 : if (ierr /= 0) then
452 0 : write(*,*) 'failed to open '//TRIM(filename)
453 0 : return
454 : end if
455 :
456 : ! Write the data
457 :
458 0 : nn = SIZE(point_data, 2)
459 :
460 0 : write(iounit, 100) nn, global_data, s%gyre_data_schema
461 : 100 format(I6, 3(1X,1PE26.16), 1X, I6)
462 :
463 0 : do j = 1, nn
464 0 : write(iounit, 110) j, point_data(:,j)
465 : 110 format(I6, 99(1X,1PE26.16))
466 : end do
467 :
468 : ! Close the file
469 :
470 0 : close(iounit)
471 :
472 0 : return
473 :
474 : end subroutine write_gyre_data
475 :
476 : end module pulse_gyre
|
