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_saio
21 :
22 : use star_private_def
23 : use const_def, only: dp, pi, lsun, rsun, clight, crad
24 : use utils_lib
25 : use atm_def
26 : use atm_support
27 : use pulse_utils
28 :
29 : implicit none
30 :
31 : private
32 : public :: get_saio_data
33 : public :: write_saio_data
34 :
35 : contains
36 :
37 0 : subroutine get_saio_data (id, &
38 : keep_surface_point, add_atmosphere, global_data, point_data, ierr)
39 :
40 : integer, intent(in) :: id
41 : logical, intent(in) :: keep_surface_point
42 : logical, intent(in) :: add_atmosphere
43 : real(dp), allocatable, intent(out) :: global_data(:)
44 : real(dp), allocatable, intent(out) :: point_data(:,:)
45 : integer, intent(out) :: ierr
46 :
47 : type(star_info), pointer :: s
48 0 : integer, allocatable :: k_a(:)
49 0 : integer, allocatable :: k_b(:)
50 : integer :: n_sg
51 : integer :: n_atm
52 : integer :: nn_atm
53 : integer :: n_env
54 : integer :: nn_env
55 : integer :: nn
56 0 : real(dp) :: r_outer
57 0 : real(dp) :: m_outer
58 : integer :: j
59 : integer :: k
60 : integer :: sg
61 :
62 : ! Get SAIO data
63 :
64 0 : call get_star_ptr(id, s, ierr)
65 0 : if (ierr /= 0) then
66 0 : write(*,*) 'bad star id for get_saio_data'
67 0 : return
68 : end if
69 :
70 : ! Set up segment indices
71 :
72 0 : call set_segment_indices(s, k_a, k_b, .FALSE.)
73 :
74 0 : n_sg = SIZE(k_a)
75 :
76 : ! Determine data dimensions
77 :
78 0 : if (add_atmosphere) then
79 0 : call build_atm(s, s%L(1), s%r(1), s%Teff, s%m_grav(1), s%cgrav(1), ierr)
80 0 : if (ierr /= 0) then
81 0 : write(*,*) 'failed in build_atm'
82 0 : return
83 : end if
84 0 : n_atm = s%atm_structure_num_pts
85 0 : nn_atm = n_atm - 1
86 : else
87 : n_atm = 0
88 : nn_atm = 0
89 : end if
90 :
91 0 : n_env = s%nz
92 :
93 0 : if (keep_surface_point) then
94 0 : nn_env = n_env + n_sg - 1
95 : else
96 0 : nn_env = n_env - 1 + n_sg - 1
97 : end if
98 :
99 0 : nn = nn_env + nn_atm
100 :
101 : ! Store global data
102 :
103 0 : allocate(global_data(4))
104 :
105 0 : r_outer = Rsun*s%photosphere_r
106 0 : m_outer = s%m_grav(1)
107 :
108 0 : global_data(1) = m_outer
109 0 : global_data(2) = log10(s%L(1)/Lsun)
110 0 : global_data(3) = log10(r_outer/Rsun)
111 0 : global_data(4) = s%star_age
112 :
113 : ! Store point data
114 :
115 0 : allocate(point_data(20,nn))
116 :
117 0 : j = 1
118 :
119 : ! Atmosphere (we skip the point at the base of the atm to
120 : ! smooth the transition)
121 :
122 0 : atm_loop : do k = 1, n_atm-1
123 0 : call store_saio_data_atm(j, k, k_a(1), k_b(1))
124 0 : j = j + 1
125 : end do atm_loop
126 :
127 : ! Envelope
128 :
129 0 : sg = 1
130 :
131 0 : env_loop : do k = 1, n_env
132 :
133 0 : if (k == 1 .AND. .NOT. keep_surface_point) cycle env_loop
134 :
135 0 : call store_saio_data_env(j, k, k_a(sg), k_b(sg))
136 0 : j = j + 1
137 :
138 0 : if (k == k_b(sg)+1) then
139 :
140 0 : sg = sg + 1
141 :
142 0 : call store_saio_data_env(j, k, k_a(sg), k_b(sg))
143 0 : j = j + 1
144 :
145 : end if
146 :
147 : end do env_loop
148 :
149 : ! Reverse the data ordering (SAIO format is center-to-surface)
150 :
151 0 : point_data = point_data(:,nn:1:-1)
152 :
153 : ! Tidy up
154 :
155 0 : if (ASSOCIATED(s%atm_structure)) then
156 0 : deallocate(s%atm_structure)
157 : end if
158 :
159 0 : return
160 :
161 : contains
162 :
163 0 : subroutine store_saio_data_atm (j, k, k_a, k_b)
164 :
165 : integer, intent(in) :: j
166 : integer, intent(in) :: k
167 : integer, intent(in) :: k_a
168 : integer, intent(in) :: k_b
169 :
170 : ! Store data associated with atmosphere point k into the
171 : ! point_data array at position j
172 :
173 : associate ( &
174 : r => point_data(1,j), &
175 : m => point_data(2,j), &
176 : Lrad => point_data(3,j), &
177 : T => point_data(4,j), &
178 : rho => point_data(5,j), &
179 : P => point_data(6,j), &
180 : eps => point_data(7,j), &
181 : kap => point_data(8,j), &
182 : c_V => point_data(9,j), &
183 : chi_rho => point_data(10,j), &
184 : chi_T => point_data(11,j), &
185 : eps_rho => point_data(12,j), &
186 : eps_T => point_data(13,j), &
187 : kap_rho => point_data(14,j), &
188 : kap_T => point_data(15,j), &
189 : nabla => point_data(16,j), &
190 : nabla_ad => point_data(17,j), &
191 : X => point_data(18,j), &
192 : Y => point_data(19,j), &
193 : L => point_data(20,j))
194 :
195 0 : r = s% r(1) + s% atm_structure(atm_delta_r,k)
196 0 : m = s% m_grav(1)
197 0 : T = exp(s% atm_structure(atm_lnT,k))
198 0 : rho = exp(s% atm_structure(atm_lnd,k))
199 0 : P = exp(s% atm_structure(atm_lnP,k))
200 0 : eps = 0d0
201 0 : kap = s% atm_structure(atm_kap,k)
202 0 : c_V = s% atm_structure(atm_cv,k)
203 0 : chi_rho = s% atm_structure(atm_chiRho,k)
204 0 : chi_T = s% atm_structure(atm_chiT,k)
205 0 : eps_rho = 0d0
206 0 : eps_T = 0d0
207 0 : kap_rho = s% atm_structure(atm_dlnkap_dlnd,k)
208 0 : kap_T = s% atm_structure(atm_dlnkap_dlnT,k)
209 0 : nabla = s% atm_structure(atm_gradT,k)
210 0 : nabla_ad = s% atm_structure(atm_grada,k)
211 0 : X = eval_face(s%dq, s%X, 1, k_a, k_b, v_lo=0d0, v_hi=1d0)
212 0 : Y = eval_face(s%dq, s%Y, 1, k_a, k_b, v_lo=0d0, v_hi=1d0)
213 0 : L = s% L(1)
214 :
215 0 : if (s%atm_structure(atm_tau,k) >= 2d0/3d0) then
216 0 : Lrad = (16*pi*clight*crad*s%cgrav(k)*m*T*T*T*T*nabla)/(3*P*kap)
217 : else
218 0 : Lrad = L
219 : end if
220 :
221 : end associate
222 :
223 0 : return
224 :
225 : end subroutine store_saio_data_atm
226 :
227 :
228 0 : subroutine store_saio_data_env (j, k, k_a, k_b)
229 :
230 : integer, intent(in) :: j
231 : integer, intent(in) :: k
232 : integer, intent(in) :: k_a
233 : integer, intent(in) :: k_b
234 :
235 : ! Store data for envelope face k into the point_data array at
236 : ! position j
237 :
238 : associate ( &
239 : r => point_data(1,j), &
240 : m => point_data(2,j), &
241 : Lrad => point_data(3,j), &
242 : T => point_data(4,j), &
243 : rho => point_data(5,j), &
244 : P => point_data(6,j), &
245 : eps => point_data(7,j), &
246 : kap => point_data(8,j), &
247 : c_V => point_data(9,j), &
248 : chi_rho => point_data(10,j), &
249 : chi_T => point_data(11,j), &
250 : eps_rho => point_data(12,j), &
251 : eps_T => point_data(13,j), &
252 : kap_rho => point_data(14,j), &
253 : kap_T => point_data(15,j), &
254 : nabla => point_data(16,j), &
255 : nabla_ad => point_data(17,j), &
256 : X => point_data(18,j), &
257 : Y => point_data(19,j), &
258 : L => point_data(20,j))
259 :
260 0 : r = s%r(k)
261 0 : m = s%m_grav(k)
262 0 : T = eval_face(s%dq, s%T, k, 1, s%nz)
263 0 : if (s%interpolate_rho_for_pulse_data) then
264 0 : rho = eval_face(s%dq, s%rho, k, k_a, k_b)
265 : else
266 0 : rho = eval_face_rho(s, k, k_a, k_b)
267 : end if
268 0 : P = eval_face(s%dq, s%Peos, k, 1, s%nz)
269 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)
270 0 : c_V = eval_face(s%dq, s%Cv, k, k_a, k_b)
271 0 : chi_rho = eval_face(s%dq, s%chiRho, k, k_a, k_b)
272 0 : chi_T = eval_face(s%dq, s%chiT, k, k_a, k_b)
273 0 : eps_rho = eval_face(s%dq, s%d_epsnuc_dlnd, k, k_a, k_b)
274 0 : eps_T = eval_face(s%dq, s%d_epsnuc_dlnT, k, k_a, k_b)
275 0 : kap_rho = eval_face(s%dq, s%d_opacity_dlnd, k, k_a, k_b)/kap
276 0 : kap_T = eval_face(s%dq, s%d_opacity_dlnT, k, k_a, k_b)/kap
277 0 : nabla = s%gradT(k) ! Not quite right; gradT can be discontinuous
278 0 : nabla_ad = eval_face(s%dq, s%grada, k, k_a, k_b)
279 0 : X = eval_face(s%dq, s%X, k, k_a, k_b, v_lo=0d0, v_hi=1d0)
280 0 : Y = eval_face(s%dq, s%Y, k, k_a, k_b, v_lo=0d0, v_hi=1d0)
281 0 : L = s%L(k)
282 :
283 0 : if (s%tau(k) >= 2d0/3d0) then
284 0 : Lrad = (16*pi*clight*crad*s%cgrav(k)*m*T*T*T*T*nabla)/(3*P*kap)
285 : else
286 0 : Lrad = L
287 : end if
288 :
289 : end associate
290 :
291 0 : return
292 :
293 : end subroutine store_saio_data_env
294 :
295 : end subroutine get_saio_data
296 :
297 :
298 0 : subroutine write_saio_data (id, filename, global_data, point_data, ierr)
299 :
300 : integer, intent(in) :: id
301 : character(*), intent(in) :: filename
302 : real(dp), intent(in) :: global_data(:)
303 : real(dp), intent(in) :: point_data(:,:)
304 : integer, intent(out) :: ierr
305 :
306 : type(star_info), pointer :: s
307 : integer :: iounit
308 : integer :: nn
309 : integer :: j
310 :
311 : ! Write SAIO data to file
312 :
313 0 : call get_star_ptr(id, s, ierr)
314 0 : if (ierr /= 0) then
315 0 : write(*,*) 'bad star id for get_fgong_info'
316 0 : return
317 : end if
318 :
319 : ! Open the file
320 :
321 0 : open(newunit=iounit, file=TRIM(filename), status='REPLACE', iostat=ierr)
322 0 : if (ierr /= 0) then
323 0 : write(*,*) 'failed to open '//TRIM(filename)
324 0 : return
325 : end if
326 :
327 : ! Write the data
328 :
329 0 : nn = SIZE(point_data, 2)
330 :
331 0 : write(iounit, 100) nn, global_data
332 : 100 format(I6,99(1PE16.9))
333 :
334 0 : do j = 1, nn
335 0 : write(iounit, 110) point_data(:,j)
336 : 110 format(5(1PE16.9))
337 : end do
338 :
339 : ! Close the file
340 :
341 0 : close(iounit)
342 :
343 0 : return
344 :
345 : end subroutine write_saio_data
346 :
347 : end module pulse_saio
|
