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 pulse_osc
21 :
22 : use star_private_def
23 : use const_def, only: dp, pi, four_thirds, rsun
24 : use utils_lib
25 : use chem_def
26 : use atm_def
27 : use atm_support
28 :
29 : use pulse_utils
30 :
31 : implicit none
32 :
33 : ! Parameter definitions (these values are for OSC format version 2K,
34 : ! with 14 abundances)
35 :
36 : integer, parameter :: ICONST = 15
37 : integer, parameter :: IVAR = 22
38 : integer, parameter :: IABUND = 14
39 :
40 : private
41 : public :: get_osc_data
42 : public :: write_osc_data
43 :
44 : contains
45 :
46 0 : subroutine get_osc_data (id, &
47 : add_center_point, keep_surface_point, add_atmosphere, global_data, point_data, ierr)
48 :
49 : integer, intent(in) :: id
50 : logical, intent(in) :: add_center_point
51 : logical, intent(in) :: keep_surface_point
52 : logical, intent(in) :: add_atmosphere
53 : real(dp), allocatable, intent(out) :: global_data(:)
54 : real(dp), allocatable, intent(out) :: point_data(:,:)
55 : integer, intent(out) :: ierr
56 :
57 : type(star_info), pointer :: s
58 0 : integer, allocatable :: k_a(:)
59 0 : integer, allocatable :: k_b(:)
60 : integer :: n_sg
61 : integer :: h1
62 : integer :: h2
63 : integer :: he3
64 : integer :: he4
65 : integer :: li7
66 : integer :: be7
67 : integer :: be9
68 : integer :: c12
69 : integer :: c13
70 : integer :: n14
71 : integer :: n15
72 : integer :: o16
73 : integer :: o17
74 : integer :: o18
75 : integer :: si28
76 : integer :: n_atm
77 : integer :: nn_atm
78 : integer :: n_env
79 : integer :: nn_env
80 : integer :: nn
81 0 : real(dp) :: r_outer
82 0 : real(dp) :: m_outer
83 0 : real(dp) :: P_c
84 0 : real(dp) :: rho_c
85 0 : real(dp) :: d2P_dr2_c
86 : integer :: j
87 : integer :: k
88 : integer :: sg
89 :
90 : ! Get model data for OSC output
91 :
92 0 : call get_star_ptr(id, s, ierr)
93 0 : if (ierr /= 0) then
94 0 : write(*,*) 'bad star id for get_osc_data'
95 0 : return
96 : end if
97 :
98 : ! Set up segment indices
99 :
100 : call set_segment_indices(s, k_a, k_b, add_center_point)
101 :
102 0 : n_sg = SIZE(k_a)
103 :
104 : ! Set up specicies indices
105 :
106 0 : h1 = s%net_iso(ih1)
107 0 : h2 = s%net_iso(ih2)
108 0 : he3 = s%net_iso(ihe3)
109 0 : he4 = s%net_iso(ihe4)
110 0 : li7 = s%net_iso(ili7)
111 0 : be7 = s%net_iso(ibe7)
112 0 : be9 = s%net_iso(ibe9)
113 0 : c12 = s%net_iso(ic12)
114 0 : c13 = s%net_iso(ic13)
115 0 : n14 = s%net_iso(in14)
116 0 : n15 = s%net_iso(in15)
117 0 : o16 = s%net_iso(io16)
118 0 : o17 = s%net_iso(io17)
119 0 : o18 = s%net_iso(io18)
120 0 : si28 = s%net_iso(isi28)
121 :
122 : ! Determine data dimensions
123 :
124 0 : if (add_atmosphere) then
125 0 : call build_atm(s, s%L(1), s%r(1), s%Teff, s%m_grav(1), s%cgrav(1), ierr)
126 0 : if (ierr /= 0) then
127 0 : write(*,*) 'failed in build_atm'
128 0 : return
129 : end if
130 0 : n_atm = s%atm_structure_num_pts
131 0 : nn_atm = n_atm - 1
132 : else
133 : n_atm = 0
134 : nn_atm = 0
135 : end if
136 :
137 0 : n_env = s%nz
138 :
139 0 : if (keep_surface_point) then
140 0 : nn_env = n_env + n_sg - 1
141 : else
142 0 : nn_env = n_env - 1 + n_sg - 1
143 : end if
144 :
145 0 : if (add_center_point) then
146 0 : nn = nn_env + nn_atm + 1
147 : else
148 0 : nn = nn_env + nn_atm
149 : end if
150 :
151 : ! Store global data
152 :
153 0 : allocate(global_data(ICONST))
154 0 : global_data = 0d0
155 :
156 0 : r_outer = Rsun*s%photosphere_r
157 0 : m_outer = s%m_grav(1)
158 :
159 0 : global_data(1) = m_outer
160 0 : global_data(2) = r_outer
161 0 : global_data(3) = s%L(1)
162 0 : global_data(4) = s%initial_z
163 0 : global_data(5) = 1d0 - (s%initial_y + s%initial_z)
164 0 : global_data(6) = s%mixing_length_alpha
165 0 : if (s%largest_conv_mixing_region /= 0) then
166 0 : k = s%mixing_region_bottom(s%largest_conv_mixing_region)
167 0 : global_data(7) = s%X(k)
168 0 : global_data(8) = s%Y(k)
169 : end if
170 :
171 0 : if (s%interpolate_rho_for_pulse_data) then
172 0 : rho_c = eval_center_rho(s, k_b(n_sg))
173 : else
174 0 : rho_c = eval_center(s%rmid, s%rho, k_a(n_sg), k_b(n_sg))
175 : end if
176 :
177 : ! at the centre d²P/dr² = -4πGρ²/3
178 0 : d2P_dr2_c = -four_thirds*pi*s% cgrav(s% nz)*rho_c**2
179 0 : P_c = s%Peos(s% nz) - 0.5d0*d2P_dr2_c*s% rmid(s% nz)**2
180 0 : global_data(9) = r_outer**2*d2P_dr2_c/P_c
181 0 : global_data(10) = r_outer**2*eval_center_d2(s%rmid, s%rho, k_a(n_sg), k_b(n_sg)) / rho_c
182 :
183 0 : global_data(11) = s%star_age
184 0 : if (s%rotation_flag) then
185 0 : global_data(12) = s%omega(1)
186 : else
187 0 : global_data(12) = 0d0
188 : end if
189 :
190 : ! global_data(13) should be the initial rotation rate, but we lack
191 : ! that datum
192 :
193 : ! Store point data
194 :
195 0 : allocate(point_data(IVAR+IABUND,nn))
196 0 : point_data = 0d0
197 :
198 0 : j = 1
199 :
200 : ! Atmosphere (we skip the point at the base of the atm to
201 : ! smooth the transition)
202 :
203 0 : atm_loop : do k = 1, n_atm-1
204 0 : call store_point_data_atm(j, k, k_a(1), k_b(1))
205 0 : j = j + 1
206 : end do atm_loop
207 :
208 : ! Envelope
209 :
210 0 : sg = 1
211 :
212 0 : env_loop : do k = 1, n_env
213 :
214 0 : if (k == 1 .AND. .NOT. keep_surface_point) cycle env_loop
215 :
216 0 : call store_point_data_env(j, k, k_a(sg), k_b(sg))
217 0 : j = j + 1
218 :
219 0 : if (k == k_b(sg)+1) then
220 :
221 0 : sg = sg + 1
222 :
223 0 : call store_point_data_env(j, k, k_a(sg), k_b(sg))
224 0 : j = j + 1
225 :
226 : end if
227 :
228 : end do env_loop
229 :
230 : ! Center
231 :
232 0 : if (add_center_point) then
233 0 : call store_point_data_ctr(j, k_a(n_sg), k_b(n_sg))
234 0 : j = j + 1
235 : end if
236 :
237 : ! Tidy up
238 :
239 0 : if (ASSOCIATED(s%atm_structure)) then
240 0 : deallocate(s%atm_structure)
241 : end if
242 :
243 0 : return
244 :
245 : contains
246 :
247 0 : subroutine store_point_data_atm (j, k, k_a, k_b)
248 :
249 : integer, intent(in) :: j
250 : integer, intent(in) :: k
251 : integer, intent(in) :: k_a
252 : integer, intent(in) :: k_b
253 :
254 0 : real(dp) :: grav
255 0 : real(dp) :: N2
256 :
257 : ! Store data associated with atmosphere point k into the point_data array at
258 : ! position j
259 :
260 : associate ( &
261 : r => point_data(1,j), &
262 : lnq => point_data(2,j), &
263 : T => point_data(3,j), &
264 : P => point_data(4,j), &
265 : rho => point_data(5,j), &
266 : nabla => point_data(6,j), &
267 : L => point_data(7,j), &
268 : kap => point_data(8,j), &
269 : eps => point_data(9,j), &
270 : Gamma_1 => point_data(10,j), &
271 : nabla_ad => point_data(11,j), &
272 : delta => point_data(12,j), &
273 : c_P => point_data(13, j), &
274 : rec_mu_e => point_data(14,j), &
275 : A_ast => point_data(15,j), &
276 : omega => point_data(16,j), &
277 : kap_T => point_data(17,j), &
278 : kap_rho => point_data(18,j), &
279 : eps_T => point_data(19,j), &
280 : eps_rho => point_data(20,j), &
281 : P_tot_gas => point_data(21,j), &
282 : nabla_rad => point_data(22,j), &
283 : X_H1 => point_data(23,j), &
284 : X_H2 => point_data(24,j), &
285 : X_He3 => point_data(25,j), &
286 : X_He4 => point_data(26,j), &
287 : X_Li7 => point_data(27,j), &
288 : X_Be7 => point_data(28,j), &
289 : X_C12 => point_data(29,j), &
290 : X_C13 => point_data(30,j), &
291 : X_N14=> point_data(31,j), &
292 : X_N15 => point_data(32,j), &
293 : X_O16 => point_data(33,j), &
294 : X_O17 => point_data(34,j), &
295 : X_Be9 => point_data(35,j), &
296 : X_Si28 => point_data(36,j))
297 :
298 0 : r = s%r(1) + s%atm_structure(atm_delta_r,k)
299 0 : lnq = log(s%m_grav(1)/m_outer)
300 0 : T = exp(s%atm_structure(atm_lnT,k))
301 0 : P = exp(s%atm_structure(atm_lnP,k))
302 0 : rho = exp(s%atm_structure(atm_lnd,k))
303 0 : nabla = s%atm_structure(atm_gradT,k)
304 0 : L = s%L(1)
305 0 : kap = s%atm_structure(atm_kap,k)
306 0 : eps = 0d0
307 0 : Gamma_1 = s%atm_structure(atm_gamma1,k)
308 0 : nabla_ad = s%atm_structure(atm_grada,k)
309 0 : delta = s%atm_structure(atm_chiT,k)/s%atm_structure(atm_chiRho,k)
310 0 : c_P = s%atm_structure(atm_cp,k)
311 0 : rec_mu_e = exp(s%atm_structure(atm_lnfree_e,k)) ! check
312 :
313 0 : grav = s%cgrav(1)*s%m_grav(1)/(r*r)
314 0 : N2 = grav*grav*(rho/P)*delta*(nabla_ad - nabla)
315 0 : A_ast = N2*r/grav
316 :
317 0 : if (s%rotation_flag) then
318 0 : omega = s%omega(1)
319 : else
320 0 : omega = 0d0
321 : end if
322 0 : kap_T = s%atm_structure(atm_dlnkap_dlnT,k)
323 0 : kap_rho = s%atm_structure(atm_dlnkap_dlnd,k)
324 0 : eps_T = 0d0
325 0 : eps_rho = 0d0
326 0 : P_tot_gas = P/exp(s%atm_structure(atm_lnPgas,k))
327 0 : nabla_rad = s%atm_structure(atm_gradr,k)
328 0 : X_H1 = eval_face_X(s, h1, 1, k_a, k_b)
329 0 : X_H2 = eval_face_X(s, h2, 1, k_a, k_b)
330 0 : X_He3 = eval_face_X(s, he3, 1, k_a, k_b)
331 0 : X_He4 = eval_face_X(s, he4, 1, k_a, k_b)
332 0 : X_Li7 = eval_face_X(s, li7, 1, k_a, k_b)
333 0 : X_Be7 = eval_face_X(s, be7, 1, k_a, k_b)
334 0 : X_C12 = eval_face_X(s, c12, 1, k_a, k_b)
335 0 : X_C13 = eval_face_X(s, c13, 1, k_a, k_b)
336 0 : X_N14 = eval_face_X(s, n14, 1, k_a, k_b)
337 0 : X_N15 = eval_face_X(s, n15, 1, k_a, k_b)
338 0 : X_O16 = eval_face_X(s, o16, 1, k_a, k_b)
339 0 : X_O17 = eval_face_X(s, o17, 1, k_a, k_b)
340 0 : X_Be9 = eval_face_X(s, be9, 1, k_a, k_b)
341 0 : X_Si28 = eval_face_X(s, si28, 1, k_a, k_b)
342 :
343 : end associate
344 :
345 0 : return
346 :
347 : end subroutine store_point_data_atm
348 :
349 :
350 0 : subroutine store_point_data_env (j, k, k_a, k_b)
351 :
352 : integer, intent(in) :: j
353 : integer, intent(in) :: k
354 : integer, intent(in) :: k_a
355 : integer, intent(in) :: k_b
356 :
357 : ! Store data associated with envelope face k into the point_data array at
358 : ! position j
359 :
360 : associate ( &
361 : r => point_data(1,j), &
362 : lnq => point_data(2,j), &
363 : T => point_data(3,j), &
364 : P => point_data(4,j), &
365 : rho => point_data(5,j), &
366 : nabla => point_data(6,j), &
367 : L => point_data(7,j), &
368 : kap => point_data(8,j), &
369 : eps => point_data(9,j), &
370 : Gamma_1 => point_data(10,j), &
371 : nabla_ad => point_data(11,j), &
372 : delta => point_data(12,j), &
373 : c_P => point_data(13, j), &
374 : rec_mu_e => point_data(14,j), &
375 : A_ast => point_data(15,j), &
376 : omega => point_data(16,j), &
377 : kap_T => point_data(17,j), &
378 : kap_rho => point_data(18,j), &
379 : eps_T => point_data(19,j), &
380 : eps_rho => point_data(20,j), &
381 : P_tot_gas => point_data(21,j), &
382 : nabla_rad => point_data(22,j), &
383 : X_H1 => point_data(23,j), &
384 : X_H2 => point_data(24,j), &
385 : X_He3 => point_data(25,j), &
386 : X_He4 => point_data(26,j), &
387 : X_Li7 => point_data(27,j), &
388 : X_Be7 => point_data(28,j), &
389 : X_C12 => point_data(29,j), &
390 : X_C13 => point_data(30,j), &
391 : X_N14=> point_data(31,j), &
392 : X_N15 => point_data(32,j), &
393 : X_O16 => point_data(33,j), &
394 : X_O17 => point_data(34,j), &
395 : X_Be9 => point_data(35,j), &
396 : X_Si28 => point_data(36,j))
397 :
398 0 : r = s%r(k)
399 0 : lnq = log(s%m_grav(k)/m_outer)
400 0 : T = eval_face(s%dq, s%T, k, 1, s%nz)
401 0 : P = eval_face(s%dq, s%Peos, k, 1, s%nz)
402 0 : if (s%interpolate_rho_for_pulse_data) then
403 0 : rho = eval_face(s%dq, s%rho, k, k_a, k_b)
404 : else
405 0 : rho = eval_face_rho(s, k, k_a, k_b)
406 : end if
407 0 : nabla = s%gradT(k) ! Not quite right; gradT can be discontinuous
408 0 : L = s%L(k)
409 0 : kap = eval_face(s%dq, s%opacity, k, k_a, k_b)
410 0 : eps = eval_face(s%dq, s%eps_nuc, k, k_a, k_b) + eval_face(s%dq, s%eps_grav_ad%val, k, k_a, k_b)
411 0 : Gamma_1 = eval_face(s%dq, s%gamma1, k, k_a, k_b)
412 0 : nabla_ad = eval_face(s%dq, s%grada, k, k_a, k_b)
413 0 : delta = eval_face(s%dq, s%chiT, k, k_a, k_b)/eval_face(s%dq, s%chiRho, k, k_a, k_b)
414 0 : c_P = eval_face(s%dq, s%cp, k, k_a, k_b)
415 0 : rec_mu_e = exp(eval_face(s%dq, s%lnfree_e, k, k_a, k_b)) ! check
416 0 : A_ast = eval_face_A_ast(s, k, k_a, k_b)
417 0 : if (s%rotation_flag) then
418 0 : omega = s%omega(k) ! Not quite right; omega can be discontinuous
419 : else
420 0 : omega = 0d0
421 : end if
422 0 : kap_T = eval_face(s%dq, s%d_opacity_dlnT, k, k_a, k_b)/kap
423 0 : kap_rho = eval_face(s%dq, s%d_opacity_dlnd, k, k_a, k_b)/kap
424 0 : eps_T = eval_face(s%dq, s%d_epsnuc_dlnT, k, k_a, k_b)
425 0 : eps_rho = eval_face(s%dq, s%d_epsnuc_dlnd, k, k_a, k_b)
426 0 : P_tot_gas = eval_face(s%dq, s%Peos, k, 1, s%nz)/eval_face(s%dq, s%Pgas, k, k_a, k_b)
427 0 : nabla_rad = s%gradr(k) ! Not quite right; gradr can be discontinuous
428 0 : X_H1 = eval_face_X(s, h1, k, k_a, k_b)
429 0 : X_H2 = eval_face_X(s, h2, k, k_a, k_b)
430 0 : X_He3 = eval_face_X(s, he3, k, k_a, k_b)
431 0 : X_He4 = eval_face_X(s, he4, k, k_a, k_b)
432 0 : X_Li7 = eval_face_X(s, li7, k, k_a, k_b)
433 0 : X_Be7 = eval_face_X(s, be7, k, k_a, k_b)
434 0 : X_C12 = eval_face_X(s, c12, k, k_a, k_b)
435 0 : X_C13 = eval_face_X(s, c13, k, k_a, k_b)
436 0 : X_N14 = eval_face_X(s, n14, k, k_a, k_b)
437 0 : X_N15 = eval_face_X(s, n15, k, k_a, k_b)
438 0 : X_O16 = eval_face_X(s, o16, k, k_a, k_b)
439 0 : X_O17 = eval_face_X(s, o17, k, k_a, k_b)
440 0 : X_Be9 = eval_face_X(s, be9, k, k_a, k_b)
441 0 : X_Si28 = eval_face_X(s, si28, k, k_a, k_b)
442 :
443 : end associate
444 :
445 0 : return
446 :
447 : end subroutine store_point_data_env
448 :
449 :
450 0 : subroutine store_point_data_ctr (j, k_a, k_b)
451 :
452 : integer, intent(in) :: j
453 : integer, intent(in) :: k_a
454 : integer, intent(in) :: k_b
455 :
456 : ! Store data for the center into the point_data array at position j
457 :
458 : associate ( &
459 : r => point_data(1,j), &
460 : lnq => point_data(2,j), &
461 : T => point_data(3,j), &
462 : P => point_data(4,j), &
463 : rho => point_data(5,j), &
464 : nabla => point_data(6,j), &
465 : L => point_data(7,j), &
466 : kap => point_data(8,j), &
467 : eps => point_data(9,j), &
468 : Gamma_1 => point_data(10,j), &
469 : nabla_ad => point_data(11,j), &
470 : delta => point_data(12,j), &
471 : c_P => point_data(13, j), &
472 : rec_mu_e => point_data(14,j), &
473 : A_ast => point_data(15,j), &
474 : omega => point_data(16,j), &
475 : kap_T => point_data(17,j), &
476 : kap_rho => point_data(18,j), &
477 : eps_T => point_data(19,j), &
478 : eps_rho => point_data(20,j), &
479 : P_tot_gas => point_data(21,j), &
480 : nabla_rad => point_data(22,j), &
481 : X_H1 => point_data(23,j), &
482 : X_H2 => point_data(24,j), &
483 : X_He3 => point_data(25,j), &
484 : X_He4 => point_data(26,j), &
485 : X_Li7 => point_data(27,j), &
486 : X_Be7 => point_data(28,j), &
487 : X_C12 => point_data(29,j), &
488 : X_C13 => point_data(30,j), &
489 : X_N14=> point_data(31,j), &
490 : X_N15 => point_data(32,j), &
491 : X_O16 => point_data(33,j), &
492 : X_O17 => point_data(34,j), &
493 : X_Be9 => point_data(35,j), &
494 : X_Si28 => point_data(36,j))
495 :
496 0 : r = 0d0
497 0 : lnq = log(TINY(0d0))
498 0 : T = eval_center(s%rmid, s%T, 1, s%nz)
499 0 : P = P_c
500 0 : rho = rho_c
501 0 : nabla = eval_center(s%r, s%gradT, k_a, k_b)
502 0 : L = 0d0
503 0 : kap = eval_center(s%rmid, s%opacity, k_a, k_b)
504 0 : eps = eval_center(s%rmid, s%eps_nuc, k_a, k_b) + eval_center(s%rmid, s%eps_grav_ad%val, k_a, k_b)
505 0 : Gamma_1 = eval_center(s%rmid, s%gamma1, k_a, k_b)
506 0 : nabla_ad = eval_center(s%rmid, s%grada, k_a, k_b)
507 0 : delta = eval_center(s%rmid, s%chiT, k_a, k_b)/eval_center(s%rmid, s%chiRho, k_a, k_b)
508 0 : c_P = eval_center(s%rmid, s%cp, k_a, k_b)
509 0 : rec_mu_e = exp(eval_center(s%rmid, s%lnfree_e, k_a, k_b)) ! check
510 0 : A_ast = point_data(15,j)
511 0 : if (s%rotation_flag) then
512 0 : omega = eval_center(s%r, s%omega, k_a, k_b)
513 : else
514 0 : omega = 0d0
515 : end if
516 0 : kap_T = eval_center(s%rmid, s%d_opacity_dlnT, k_a, k_b)/kap
517 0 : kap_rho = eval_center(s%rmid, s%d_opacity_dlnd, k_a, k_b)/kap
518 0 : eps_T = eval_center(s%rmid, s%d_epsnuc_dlnT, k_a, k_b)
519 0 : eps_rho = eval_center(s%rmid, s%d_epsnuc_dlnd, k_a, k_b)
520 0 : P_tot_gas = eval_center(s%rmid, s%Peos, 1, s%nz)/eval_center(s%rmid, s%Pgas, k_a, k_b)
521 0 : nabla_rad = eval_center(s%r, s%gradr, k_a, k_b)
522 0 : X_H1 = eval_center_X(s, h1, k_a, k_b)
523 0 : X_H2 = eval_center_X(s, h2, k_a, k_b)
524 0 : X_He3 = eval_center_X(s, he3, k_a, k_b)
525 0 : X_He4 = eval_center_X(s, he4, k_a, k_b)
526 0 : X_Li7 = eval_center_X(s, li7, k_a, k_b)
527 0 : X_Be7 = eval_center_X(s, be7, k_a, k_b)
528 0 : X_C12 = eval_center_X(s, c12, k_a, k_b)
529 0 : X_C13 = eval_center_X(s, c13, k_a, k_b)
530 0 : X_N14 = eval_center_X(s, n14, k_a, k_b)
531 0 : X_N15 = eval_center_X(s, n15, k_a, k_b)
532 0 : X_O16 = eval_center_X(s, o16, k_a, k_b)
533 0 : X_O17 = eval_center_X(s, o17, k_a, k_b)
534 0 : X_Be9 = eval_center_X(s, be9, k_a, k_b)
535 0 : X_Si28 = eval_center_X(s, si28, k_a, k_b)
536 :
537 : end associate
538 :
539 0 : return
540 :
541 : end subroutine store_point_data_ctr
542 :
543 : end subroutine get_osc_data
544 :
545 :
546 0 : subroutine write_osc_data (id, filename, global_data, point_data, ierr)
547 :
548 : integer, intent(in) :: id
549 : character(*), intent(in) :: filename
550 : real(dp), intent(in) :: global_data(:)
551 : real(dp), intent(in) :: point_data(:,:)
552 : integer, intent(out) :: ierr
553 :
554 : type(star_info), pointer :: s
555 : integer :: iounit
556 : integer :: n_global
557 : integer :: n_point
558 : integer :: nn
559 : integer :: i
560 : integer :: j
561 : integer :: k
562 :
563 : ! Write OSC data to file
564 :
565 0 : call get_star_ptr(id, s, ierr)
566 0 : if (ierr /= 0) then
567 0 : write(*,*) 'bad star id for write_osc_data'
568 0 : return
569 : end if
570 :
571 : ! Open the file
572 :
573 0 : open(newunit=iounit, file=TRIM(filename), status='REPLACE', iostat=ierr)
574 0 : if (ierr /= 0) then
575 0 : write(*,*) 'failed to open '//TRIM(filename)
576 0 : return
577 : end if
578 :
579 : ! Write the data
580 :
581 0 : n_global = SIZE(global_data)
582 0 : n_point = SIZE(point_data, 1)
583 :
584 0 : nn = SIZE(point_data, 2)
585 :
586 0 : write(iounit, *) 'OSC file'
587 0 : write(iounit, *) 'Created by MESAstar version', version_number
588 0 : write(iounit, *)
589 0 : write(iounit, *)
590 :
591 0 : write(iounit, '(a)') '14 H1 H2 He3 He4 Li7 Be7 C12 C13 N14 N15 O16 O17 Be9 Si28'
592 0 : write(iounit, '(5I10)') nn, ICONST, IVAR, IABUND, 2000
593 :
594 0 : write(iounit, s%format_for_osc_data) (global_data(i), i=1,n_global)
595 :
596 0 : do k = 1, nn
597 0 : write(iounit, s%format_for_osc_data) (point_data(j,k), j=1,n_point)
598 : end do
599 :
600 : ! Close the file
601 :
602 0 : close(iounit)
603 :
604 0 : return
605 :
606 : end subroutine write_osc_data
607 :
608 : end module pulse_osc
|
