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_fgong
21 :
22 : use star_private_def
23 : use const_def, only: dp, pi, rsun, four_thirds, no_mixing, standard_cgrav
24 : use utils_lib
25 : use chem_def
26 : use atm_def
27 : use eos_def, only: i_Gamma1, i_lnPgas, i_chiRho, i_chiT
28 : use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results
29 : use atm_support
30 : use eos_support
31 : use pulse_utils
32 :
33 : implicit none
34 :
35 : ! Parameter definitions (these values are for FGONG format versions
36 : ! 300 & 1300)
37 :
38 : integer, parameter :: ICONST = 15
39 : integer, parameter :: IVAR = 40
40 :
41 : private
42 : public :: get_fgong_data
43 : public :: write_fgong_data
44 :
45 : contains
46 :
47 0 : subroutine get_fgong_data (id, &
48 : add_center_point, keep_surface_point, add_atmosphere, global_data, point_data, ierr)
49 :
50 : integer, intent(in) :: id
51 : logical, intent(in) :: add_center_point
52 : logical, intent(in) :: keep_surface_point
53 : logical, intent(in) :: add_atmosphere
54 : real(dp), allocatable, intent(out) :: global_data(:)
55 : real(dp), allocatable, intent(out) :: point_data(:,:)
56 : integer, intent(out) :: ierr
57 :
58 : type(star_info), pointer :: s
59 0 : integer, allocatable :: k_a(:)
60 0 : integer, allocatable :: k_b(:)
61 : integer :: n_sg
62 : integer :: h1
63 : integer :: h2
64 : integer :: he3
65 : integer :: he4
66 : integer :: li7
67 : integer :: be7
68 : integer :: c12
69 : integer :: c13
70 : integer :: n14
71 : integer :: n15
72 : integer :: o16
73 : integer :: o17
74 : integer :: o18
75 : integer :: ne20
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 : real(dp), dimension(num_eos_basic_results) :: &
91 0 : res, res_loX, res_hiX, dres_dlnRho, dres_dlnT
92 :
93 : ! Rather than allocating these arrays before use in each region,
94 : ! allocate them in outer calling subroutine.
95 : ! ∂lnΓ₁/∂… blocks in each region can probably be refactored.
96 0 : real(dp), allocatable :: dres_dxa(:,:)
97 0 : real(dp), allocatable :: xa(:)
98 : real(dp), parameter :: dxa = 1d-3 ! perturbation of abundances for ∂lnΓ₁/∂Y
99 :
100 : ! Get FGONG data
101 :
102 0 : call get_star_ptr(id, s, ierr)
103 0 : if (ierr /= 0) then
104 0 : write(*,*) 'bad star id for get_fgong_data'
105 0 : return
106 : end if
107 :
108 : ! Set up segment indices
109 :
110 : call set_segment_indices(s, k_a, k_b, add_center_point)
111 :
112 0 : n_sg = SIZE(k_a)
113 :
114 : ! Set up specicies indices
115 :
116 0 : h1 = s%net_iso(ih1)
117 0 : h2 = s%net_iso(ih2)
118 0 : he3 = s%net_iso(ihe3)
119 0 : he4 = s%net_iso(ihe4)
120 0 : li7 = s%net_iso(ili7)
121 0 : be7 = s%net_iso(ibe7)
122 0 : c12 = s%net_iso(ic12)
123 0 : c13 = s%net_iso(ic13)
124 0 : n14 = s%net_iso(in14)
125 0 : n15 = s%net_iso(in15)
126 0 : o16 = s%net_iso(io16)
127 0 : o17 = s%net_iso(io17)
128 0 : o18 = s%net_iso(io18)
129 0 : ne20 = s%net_iso(ine20)
130 :
131 : ! Determine data dimensions
132 :
133 0 : allocate(dres_dxa(num_eos_basic_results, s% species))
134 0 : allocate(xa(s% species))
135 :
136 0 : if (add_atmosphere) then
137 0 : call build_atm(s, s%L(1), s%r(1), s%Teff, s%m_grav(1), s%cgrav(1), ierr)
138 0 : if (ierr /= 0) then
139 0 : write(*,*) 'failed in build_atm'
140 0 : return
141 : end if
142 0 : n_atm = s%atm_structure_num_pts
143 0 : nn_atm = n_atm - 1
144 : else
145 : n_atm = 0
146 : nn_atm = 0
147 : end if
148 :
149 0 : n_env = s%nz
150 :
151 0 : if (keep_surface_point) then
152 0 : nn_env = n_env + n_sg - 1
153 : else
154 0 : nn_env = n_env - 1 + n_sg - 1
155 : end if
156 :
157 0 : if (add_center_point) then
158 0 : nn = nn_env + nn_atm + 1
159 : else
160 0 : nn = nn_env + nn_atm
161 : end if
162 :
163 : ! Store global data
164 :
165 0 : allocate(global_data(ICONST))
166 0 : global_data = 0d0
167 :
168 0 : r_outer = Rsun*s%photosphere_r
169 0 : m_outer = s%m_grav(1)
170 :
171 0 : global_data(1) = m_outer
172 0 : global_data(2) = r_outer
173 0 : global_data(3) = s%L(1)
174 0 : global_data(4) = s%initial_z
175 0 : global_data(5) = 1d0 - (s%initial_y + s%initial_z)
176 0 : global_data(6) = s%mixing_length_alpha
177 :
178 : ! global_data(7)-global_data(10) are not used
179 :
180 0 : if (s%interpolate_rho_for_pulse_data) then
181 0 : rho_c = eval_center_rho(s, k_b(n_sg))
182 : else
183 0 : rho_c = eval_center(s%rmid, s%rho, k_a(n_sg), k_b(n_sg))
184 : end if
185 :
186 : ! at the centre d²P/dr² = -4πGρ²/3
187 0 : d2P_dr2_c = -four_thirds*pi*s% cgrav(s% nz)*rho_c**2
188 0 : P_c = s%Peos(s% nz) - 0.5d0*d2P_dr2_c*s% rmid(s% nz)**2
189 0 : global_data(11) = r_outer**2*d2P_dr2_c/P_c
190 0 : global_data(12) = r_outer**2*eval_center_d2(s%rmid, s%rho, k_a(n_sg), k_b(n_sg)) / rho_c
191 0 : global_data(13) = s%star_age
192 0 : global_data(14) = s%Teff
193 0 : global_data(15) = standard_cgrav
194 :
195 : ! Store point data
196 :
197 0 : allocate(point_data(IVAR,nn))
198 0 : point_data = 0d0
199 :
200 0 : j = 1
201 :
202 : ! Atmosphere (we skip the point at the base of the atm to smooth
203 : ! the transition)
204 :
205 0 : atm_loop : do k = 1, n_atm-1
206 0 : call store_point_data_atm(j, k, k_a(1), k_b(1))
207 0 : j = j + 1
208 : end do atm_loop
209 :
210 : ! Envelope
211 :
212 0 : sg = 1
213 :
214 0 : env_loop : do k = 1, n_env
215 :
216 0 : if (k == 1 .AND. .NOT. keep_surface_point) cycle env_loop
217 :
218 0 : call store_point_data_env(j, k, k_a(sg), k_b(sg))
219 0 : j = j + 1
220 :
221 0 : if (k == k_b(sg)+1) then
222 :
223 0 : sg = sg + 1
224 :
225 0 : call store_point_data_env(j, k, k_a(sg), k_b(sg))
226 0 : j = j + 1
227 :
228 : end if
229 :
230 : end do env_loop
231 :
232 : ! Center
233 :
234 0 : if (add_center_point) then
235 0 : call store_point_data_ctr(j, k_a(n_sg), k_b(n_sg))
236 0 : j = j + 1
237 : end if
238 :
239 : ! Tidy up
240 :
241 0 : if (ASSOCIATED(s%atm_structure)) then
242 0 : deallocate(s%atm_structure)
243 : end if
244 :
245 0 : deallocate(dres_dxa)
246 0 : deallocate(xa)
247 :
248 0 : return
249 :
250 : contains
251 :
252 0 : subroutine store_point_data_atm (j, k, k_a, k_b)
253 :
254 : integer, intent(in) :: j
255 : integer, intent(in) :: k
256 : integer, intent(in) :: k_a
257 : integer, intent(in) :: k_b
258 :
259 0 : real(dp) :: grav
260 0 : real(dp) :: N2
261 :
262 : ! Store data associated with atmosphere point k into the
263 : ! point_data array at position j
264 :
265 : associate ( &
266 : r => point_data(1,j), &
267 : lnq => point_data(2,j), &
268 : T => point_data(3,j), &
269 : P => point_data(4,j), &
270 : rho => point_data(5,j), &
271 : X => point_data(6,j), &
272 : L => point_data(7,j), &
273 : kap => point_data(8,j), &
274 : eps => point_data(9,j), &
275 : Gamma_1 => point_data(10,j), &
276 : nabla_ad => point_data(11,j), &
277 : delta => point_data(12,j), &
278 : c_P => point_data(13,j), &
279 : rec_mu_e => point_data(14,j), &
280 : A_ast => point_data(15,j), &
281 : r_X => point_data(16,j), &
282 : Z => point_data(17,j), &
283 : R_r => point_data(18,j), &
284 : eps_grav => point_data(19,j), &
285 : X_He3 => point_data(21,j), &
286 : X_C12 => point_data(22,j), &
287 : X_C13 => point_data(23,j), &
288 : X_N14 => point_data(24,j), &
289 : X_O16 => point_data(25,j), &
290 : dlnGamma1_dlnrho => point_data(26,j), &
291 : dlnGamma1_dlnP => point_data(27,j), &
292 : dlnGamma1_dY => point_data(28,j), &
293 : X_H2 => point_data(29,j), &
294 : X_He4 => point_data(30,j), &
295 : X_Li7 => point_data(31,j), &
296 : X_Be7 => point_data(32,j), &
297 : X_N15 => point_data(33,j), &
298 : X_O17 => point_data(34,j), &
299 : X_O18 => point_data(35,j), &
300 : X_Ne20 => point_data(36,j))
301 :
302 0 : r = s%r(1) + s%atm_structure(atm_delta_r,k)
303 0 : lnq = log(s%m_grav(1)/m_outer)
304 0 : T = exp(s%atm_structure(atm_lnT,k))
305 0 : P = exp(s%atm_structure(atm_lnP,k))
306 0 : rho = exp(s%atm_structure(atm_lnd,k))
307 0 : X = eval_face(s%dq, s%X, 1, k_a, k_b, v_lo=0d0, v_hi=1d0)
308 0 : L = s%L(1)
309 0 : kap = s%atm_structure(atm_kap,k)
310 0 : eps = 0d0
311 0 : Gamma_1 = s%atm_structure(atm_gamma1,k)
312 0 : nabla_ad = s%atm_structure(atm_grada,k)
313 0 : delta = s%atm_structure(atm_chiT,k)/s%atm_structure(atm_chiRho,k)
314 0 : c_P = s%atm_structure(atm_cp,k)
315 0 : rec_mu_e = exp(s%atm_structure(atm_lnfree_e,k)) ! check
316 :
317 0 : grav = s%cgrav(1)*s%m_grav(1)/(r*r)
318 0 : N2 = grav*grav*(rho/P)*delta*(nabla_ad - s%atm_structure(atm_gradT,k))
319 0 : A_ast = N2*r/grav
320 :
321 0 : r_X = 0d0 ! dxdt_nuc_h1
322 : Z = MIN(1d0, MAX(0d0, 1d0 - &
323 : eval_face(s%dq, s%X, 1, k_a, k_b) - &
324 0 : eval_face(s%dq, s%Y, 1, k_a, k_b)))
325 0 : R_r = r_outer - r
326 0 : eps_grav = 0d0
327 :
328 0 : X_He3 = eval_face_X(s, he3, 1, k_a, k_b)
329 0 : X_C12 = eval_face_X(s, c12, 1, k_a, k_b)
330 0 : X_C13 = eval_face_X(s, c13, 1, k_a, k_b)
331 0 : X_N14 = eval_face_X(s, n14, 1, k_a, k_b)
332 0 : X_O16 = eval_face_X(s, o16, 1, k_a, k_b)
333 0 : X_H2 = eval_face_X(s, h2, 1, k_a, k_b)
334 0 : X_He4 = eval_face_X(s, he4, 1, k_a, k_b)
335 0 : X_Li7 = eval_face_X(s, li7, 1, k_a, k_b)
336 0 : X_Be7 = eval_face_X(s, be7, 1, k_a, k_b)
337 0 : X_N15 = eval_face_X(s, n15, 1, k_a, k_b)
338 0 : X_O17 = eval_face_X(s, o17, 1, k_a, k_b)
339 0 : X_O18 = eval_face_X(s, o18, 1, k_a, k_b)
340 0 : X_Ne20 = eval_face_X(s, ne20, 1, k_a, k_b)
341 :
342 : ! Γ₁ derivatives
343 0 : xa(:) = s% xa(:,1)
344 :
345 : ! evaluate EoS at low X
346 0 : xa(h1) = s% xa(h1,1) - dxa
347 0 : xa(he4) = s% xa(he4,1) + dxa
348 : call get_eos(s, k, xa(:), rho, log10(rho), T, log10(T), &
349 0 : res_loX, dres_dlnRho, dres_dlnT, dres_dxa, ierr)
350 :
351 : ! evaluate EoS at high X
352 0 : xa(h1) = s% xa(h1,1) + dxa
353 0 : xa(he4) = s% xa(he4,1) - dxa
354 : call get_eos(s, k, xa(:), rho, log10(rho), T, log10(T), &
355 0 : res_hiX, dres_dlnRho, dres_dlnT, dres_dxa, ierr)
356 :
357 : ! evaluate EoS at actual X
358 : call get_eos(s, k, s% xa(:,1), rho, log10(rho), T, log10(T), &
359 0 : res, dres_dlnRho, dres_dlnT, dres_dxa, ierr)
360 :
361 : ! use dres_dxa(:,1) to store finite difference,
362 : ! which is evaluated at constant (ρ,T)
363 0 : dres_dxa(:,1) = (res_hiX-res_loX)/2/dxa
364 :
365 : dlnGamma1_dlnRho = dres_dlnRho(i_Gamma1) &
366 0 : - dres_dlnT(i_Gamma1)*res(i_chiRho)/res(i_chiT)
367 0 : dlnGamma1_dlnRho = dlnGamma1_dlnRho/res(i_Gamma1)
368 :
369 0 : dlnGamma1_dlnP = dres_dlnT(i_Gamma1)/res(i_chiT)/res(i_Gamma1)
370 :
371 : dlnGamma1_dY = dres_dxa(i_Gamma1,1) &
372 0 : - dres_dlnT(i_Gamma1)*dres_dxa(i_lnPgas,1)/res(i_chiT)
373 0 : dlnGamma1_dY = -dlnGamma1_dY/res(i_Gamma1) ! d/dY ~ -d/dX
374 :
375 : end associate
376 :
377 0 : return
378 :
379 : end subroutine store_point_data_atm
380 :
381 :
382 0 : subroutine store_point_data_env (j, k, k_a, k_b)
383 :
384 : integer, intent(in) :: j
385 : integer, intent(in) :: k
386 : integer, intent(in) :: k_a
387 : integer, intent(in) :: k_b
388 :
389 : ! Store data for envelope face k into the point_data array at
390 : ! position j
391 :
392 : associate ( &
393 : r => point_data(1,j), &
394 : lnq => point_data(2,j), &
395 : T => point_data(3,j), &
396 : P => point_data(4,j), &
397 : rho => point_data(5,j), &
398 : X => point_data(6,j), &
399 : L => point_data(7,j), &
400 : kap => point_data(8,j), &
401 : eps => point_data(9,j), &
402 : Gamma_1 => point_data(10,j), &
403 : nabla_ad => point_data(11,j), &
404 : delta => point_data(12,j), &
405 : c_P => point_data(13,j), &
406 : rec_mu_e => point_data(14,j), &
407 : A_ast => point_data(15,j), &
408 : r_X => point_data(16,j), &
409 : Z => point_data(17,j), &
410 : R_r => point_data(18,j), &
411 : eps_grav => point_data(19,j), &
412 : X_He3 => point_data(21,j), &
413 : X_C12 => point_data(22,j), &
414 : X_C13 => point_data(23,j), &
415 : X_N14 => point_data(24,j), &
416 : X_O16 => point_data(25,j), &
417 : dlnGamma1_dlnrho => point_data(26,j), &
418 : dlnGamma1_dlnP => point_data(27,j), &
419 : dlnGamma1_dY => point_data(28,j), &
420 : X_H2 => point_data(29,j), &
421 : X_He4 => point_data(30,j), &
422 : X_Li7 => point_data(31,j), &
423 : X_Be7 => point_data(32,j), &
424 : X_N15 => point_data(33,j), &
425 : X_O17 => point_data(34,j), &
426 : X_O18 => point_data(35,j), &
427 : X_Ne20 => point_data(36,j))
428 :
429 0 : r = s%r(k)
430 0 : lnq = log(s%m_grav(k)/m_outer)
431 0 : T = eval_face(s%dq, s%T, k, 1, s%nz)
432 0 : P = eval_face(s%dq, s%Peos, k, 1, s%nz)
433 0 : if (s%interpolate_rho_for_pulse_data) then
434 0 : rho = eval_face(s%dq, s%rho, k, k_a, k_b)
435 : else
436 0 : rho = eval_face_rho(s, k, k_a, k_b)
437 : end if
438 0 : X = eval_face(s%dq, s%X, k, k_a, k_b, v_lo=0d0, v_hi=1d0)
439 0 : L = s%L(k)
440 0 : kap = eval_face(s%dq, s%opacity, k, k_a, k_b)
441 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)
442 0 : Gamma_1 = eval_face(s%dq, s%gamma1, k, k_a, k_b)
443 0 : nabla_ad = eval_face(s%dq, s%grada, k, k_a, k_b)
444 0 : delta = eval_face(s%dq, s%chiT, k, k_a, k_b)/eval_face(s%dq, s%chiRho, k, k_a, k_b)
445 0 : c_P = eval_face(s%dq, s%cp, k, k_a, k_b)
446 0 : rec_mu_e = exp(eval_face(s%dq, s%lnfree_e, k, k_a, k_b))
447 0 : if (r <= s%fgong_zero_A_inside_r*Rsun .AND. s%mixing_type(k) /= no_mixing) then
448 0 : A_ast = 0d0
449 : else
450 0 : A_ast = eval_face_A_ast(s, k, k_a, k_b)
451 : end if
452 0 : r_X = eval_face(s%dq, s%dxdt_nuc(h1,:), k, k_a, k_b)
453 : Z = MIN(1d0, MAX(0d0, 1d0 - &
454 : eval_face(s%dq, s%X, k, k_a, k_b) - &
455 0 : eval_face(s%dq, s%Y, k, k_a, k_b)))
456 0 : R_r = r_outer - s%r(k)
457 0 : eps_grav = eval_face(s%dq, s%eps_grav_ad%val, k, k_a, k_b)
458 0 : X_He3 = eval_face_X(s, he3, k, k_a, k_b)
459 0 : X_C12 = eval_face_X(s, c12, k, k_a, k_b)
460 0 : X_C13 = eval_face_X(s, c13, k, k_a, k_b)
461 0 : X_N14 = eval_face_X(s, n14, k, k_a, k_b)
462 0 : X_O16 = eval_face_X(s, o16, k, k_a, k_b)
463 0 : X_H2 = eval_face_X(s, h2, k, k_a, k_b)
464 0 : X_He4 = eval_face_X(s, he4, k, k_a, k_b)
465 0 : X_Li7 = eval_face_X(s, li7, k, k_a, k_b)
466 0 : X_Be7 = eval_face_X(s, be7, k, k_a, k_b)
467 0 : X_N15 = eval_face_X(s, n15, k, k_a, k_b)
468 0 : X_O17 = eval_face_X(s, o17, k, k_a, k_b)
469 0 : X_O18 = eval_face_X(s, o18, k, k_a, k_b)
470 0 : X_Ne20 = eval_face_X(s, ne20, k, k_a, k_b)
471 :
472 : ! Γ₁ derivatives
473 0 : xa(:) = s% xa(:,k)
474 :
475 : ! evaluate EoS at low X
476 0 : xa(h1) = s% xa(h1,k) - dxa
477 0 : xa(he4) = s% xa(he4,k) + dxa
478 : call get_eos(s, k, xa(:), rho, log10(rho), T, log10(T), &
479 0 : res_loX, dres_dlnRho, dres_dlnT, dres_dxa, ierr)
480 :
481 : ! evaluate EoS at high X
482 0 : xa(h1) = s% xa(h1,k) + dxa
483 0 : xa(he4) = s% xa(he4,k) - dxa
484 : call get_eos(s, k, xa(:), rho, log10(rho), T, log10(T), &
485 0 : res_hiX, dres_dlnRho, dres_dlnT, dres_dxa, ierr)
486 :
487 : ! evaluate EoS at actual X
488 : call get_eos(s, k, s% xa(:,k), rho, log10(rho), T, log10(T), &
489 0 : res, dres_dlnRho, dres_dlnT, dres_dxa, ierr)
490 :
491 : ! use dres_dxa(:,1) to store finite difference,
492 : ! which is evaluated at constant (ρ,T)
493 0 : dres_dxa(:,1) = (res_hiX-res_loX)/2/dxa
494 :
495 : dlnGamma1_dlnRho = dres_dlnRho(i_Gamma1) &
496 0 : - dres_dlnT(i_Gamma1)*res(i_chiRho)/res(i_chiT)
497 0 : dlnGamma1_dlnRho = dlnGamma1_dlnRho/res(i_Gamma1)
498 :
499 0 : dlnGamma1_dlnP = dres_dlnT(i_Gamma1)/res(i_chiT)/res(i_Gamma1)
500 :
501 : dlnGamma1_dY = dres_dxa(i_Gamma1,1) &
502 0 : - dres_dlnT(i_Gamma1)*dres_dxa(i_lnPgas,1)/res(i_chiT)
503 0 : dlnGamma1_dY = -dlnGamma1_dY/res(i_Gamma1) ! d/dY ~ -d/dX
504 :
505 : end associate
506 :
507 0 : return
508 :
509 : end subroutine store_point_data_env
510 :
511 :
512 0 : subroutine store_point_data_ctr (j, k_a, k_b)
513 :
514 : integer, intent(in) :: j
515 : integer, intent(in) :: k_a
516 : integer, intent(in) :: k_b
517 :
518 : ! Store data for the center into the point_data array at position j
519 :
520 : associate ( &
521 : r => point_data(1,j), &
522 : lnq => point_data(2,j), &
523 : T => point_data(3,j), &
524 : P => point_data(4,j), &
525 : rho => point_data(5,j), &
526 : X => point_data(6,j), &
527 : L => point_data(7,j), &
528 : kap => point_data(8,j), &
529 : eps => point_data(9,j), &
530 : Gamma_1 => point_data(10,j), &
531 : nabla_ad => point_data(11,j), &
532 : delta => point_data(12,j), &
533 : c_P => point_data(13,j), &
534 : rec_mu_e => point_data(14,j), &
535 : A_ast => point_data(15,j), &
536 : r_X => point_data(16,j), &
537 : Z => point_data(17,j), &
538 : R_r => point_data(18,j), &
539 : eps_grav => point_data(19,j), &
540 : X_He3 => point_data(21,j), &
541 : X_C12 => point_data(22,j), &
542 : X_C13 => point_data(23,j), &
543 : X_N14 => point_data(24,j), &
544 : X_O16 => point_data(25,j), &
545 : dlnGamma1_dlnrho => point_data(26,j), &
546 : dlnGamma1_dlnP => point_data(27,j), &
547 : dlnGamma1_dY => point_data(28,j), &
548 : X_H2 => point_data(29,j), &
549 : X_He4 => point_data(30,j), &
550 : X_Li7 => point_data(31,j), &
551 : X_Be7 => point_data(32,j), &
552 : X_N15 => point_data(33,j), &
553 : X_O17 => point_data(34,j), &
554 : X_O18 => point_data(35,j), &
555 : X_Ne20 => point_data(36,j))
556 :
557 0 : r = 0d0
558 0 : lnq = log(TINY(0d0))
559 0 : T = eval_center(s%rmid, s%T, 1, s%nz)
560 0 : P = P_c
561 0 : rho = rho_c
562 0 : X = eval_center(s%rmid, s%X, k_a, k_b, v_lo=0d0, v_hi=1d0)
563 0 : L = 0d0
564 0 : kap = eval_center(s%rmid, s%opacity, k_a, k_b)
565 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)
566 0 : Gamma_1 = eval_center(s%rmid, s%gamma1, k_a, k_b)
567 0 : nabla_ad = eval_center(s%rmid, s%grada, k_a, k_b)
568 0 : delta = eval_center(s%rmid, s%chiT, k_a, k_b)/eval_center(s%rmid, s%chiRho, k_a, k_b)
569 0 : c_P = eval_center(s%rmid, s%cp, k_a, k_b)
570 0 : rec_mu_e = exp(eval_center(s%rmid, s%lnfree_e, k_a, k_b))
571 0 : A_ast = 0d0
572 0 : r_X = eval_center(s%rmid, s%dxdt_nuc(h1,:), k_a, k_b)
573 : Z = MIN(1d0, MAX(0d0, 1d0 - &
574 : eval_center(s%rmid, s%X, k_a, k_b) - &
575 0 : eval_center(s%rmid, s%Y, k_a, k_b)))
576 0 : R_r = r_outer
577 0 : eps_grav = eval_center(s%rmid, s%eps_grav_ad%val, k_a, k_b)
578 0 : X_He3 = eval_center_X(s, he3, k_a, k_b)
579 0 : X_C12 = eval_center_X(s, c12, k_a, k_b)
580 0 : X_C13 = eval_center_X(s, c13, k_a, k_b)
581 0 : X_N14 = eval_center_X(s, n14, k_a, k_b)
582 0 : X_O16 = eval_center_X(s, o16, k_a, k_b)
583 0 : X_H2 = eval_center_X(s, h2, k_a, k_b)
584 0 : X_He4 = eval_center_X(s, he4, k_a, k_b)
585 0 : X_Li7 = eval_center_X(s, li7, k_a, k_b)
586 0 : X_Be7 = eval_center_X(s, be7, k_a, k_b)
587 0 : X_N15 = eval_center_X(s, n15, k_a, k_b)
588 0 : X_O17 = eval_center_X(s, o17, k_a, k_b)
589 0 : X_O18 = eval_center_X(s, o18, k_a, k_b)
590 0 : X_Ne20 = eval_center_X(s, ne20, k_a, k_b)
591 :
592 : ! Γ₁ derivatives
593 : ! I assume the abundances at r=0 are the same as in the innermost cell
594 : ! This could be improved using eval_center_X
595 0 : xa(:) = s% xa(:, s% nz)
596 :
597 : ! evaluate EoS at low X
598 0 : xa(h1) = s% xa(h1, s% nz) - dxa
599 0 : xa(he4) = s% xa(he4, s% nz) + dxa
600 : call get_eos(s, k, xa(:), rho, log10(rho), T, log10(T), &
601 0 : res_loX, dres_dlnRho, dres_dlnT, dres_dxa, ierr)
602 :
603 : ! evaluate EoS at high X
604 0 : xa(h1) = s% xa(h1, s% nz) + dxa
605 0 : xa(he4) = s% xa(he4, s% nz) - dxa
606 : call get_eos(s, k, xa(:), rho, log10(rho), T, log10(T), &
607 0 : res_hiX, dres_dlnRho, dres_dlnT, dres_dxa, ierr)
608 :
609 : ! evaluate EoS at actual X
610 : call get_eos(s, k, s% xa(:, s% nz), rho, log10(rho), T, log10(T), &
611 0 : res, dres_dlnRho, dres_dlnT, dres_dxa, ierr)
612 :
613 : ! use dres_dxa(:,1) to store finite difference,
614 : ! which is evaluated at constant (ρ,T)
615 0 : dres_dxa(:,1) = (res_hiX-res_loX)/2/dxa
616 :
617 : dlnGamma1_dlnRho = dres_dlnRho(i_Gamma1) &
618 0 : - dres_dlnT(i_Gamma1)*res(i_chiRho)/res(i_chiT)
619 0 : dlnGamma1_dlnRho = dlnGamma1_dlnRho/res(i_Gamma1)
620 :
621 0 : dlnGamma1_dlnP = dres_dlnT(i_Gamma1)/res(i_chiT)/res(i_Gamma1)
622 :
623 : dlnGamma1_dY = dres_dxa(i_Gamma1,1) &
624 0 : - dres_dlnT(i_Gamma1)*dres_dxa(i_lnPgas,1)/res(i_chiT)
625 0 : dlnGamma1_dY = -dlnGamma1_dY/res(i_Gamma1) ! d/dY ~ -d/dX
626 :
627 : end associate
628 :
629 0 : return
630 :
631 : end subroutine store_point_data_ctr
632 :
633 : end subroutine get_fgong_data
634 :
635 :
636 0 : subroutine write_fgong_data (id, filename, global_data, point_data, ierr)
637 :
638 : integer, intent(in) :: id
639 : character(*), intent(in) :: filename
640 : real(dp), intent(in) :: global_data(:)
641 : real(dp), intent(in) :: point_data(:,:)
642 : integer, intent(out) :: ierr
643 :
644 : type(star_info), pointer :: s
645 : integer :: iounit
646 : integer :: n_global
647 : integer :: n_point
648 : integer :: nn
649 : integer :: i
650 : integer :: j
651 : integer :: k
652 : character(len=strlen) :: format_for_fgong_data
653 :
654 : ! Write FGONG data to file
655 :
656 0 : call get_star_ptr(id, s, ierr)
657 0 : if (ierr /= 0) then
658 0 : write(*,*) 'bad star id for write_fgong_info'
659 0 : return
660 : end if
661 :
662 : ! Set float format from version number (ivers)
663 :
664 0 : if (s% fgong_ivers == 300) then
665 0 : format_for_fgong_data = '(1P5E16.9,x)'
666 0 : else if (s% fgong_ivers == 1300) then
667 0 : format_for_fgong_data = '(1P,5(X,E26.18E3))'
668 : else
669 0 : write(*,*) ''
670 0 : write(*,'(a,i4)') 'bad fgong_ivers: must be 300 or 1300, not ', s% fgong_ivers
671 0 : ierr = 1
672 0 : return
673 : end if
674 :
675 : ! Open the file
676 :
677 0 : open(newunit=iounit, file=TRIM(filename), status='REPLACE', iostat=ierr)
678 0 : if (ierr /= 0) then
679 0 : write(*,*) 'failed to open '//TRIM(filename)
680 0 : return
681 : end if
682 :
683 : ! Write the data
684 :
685 0 : n_global = SIZE(global_data)
686 0 : n_point = SIZE(point_data, 1)
687 :
688 0 : nn = SIZE(point_data, 2)
689 :
690 0 : do k = 1, 4
691 0 : write(iounit, *) trim(s% fgong_header(k))
692 : end do
693 :
694 0 : write(iounit,'(4I10)') nn, ICONST, IVAR, s% fgong_ivers
695 :
696 0 : write(iounit, format_for_fgong_data) (global_data(i), i=1,n_global)
697 :
698 0 : do k = 1, nn
699 0 : write(iounit, format_for_fgong_data) (point_data(j,k), j=1,n_point)
700 : end do
701 :
702 : ! Close the file
703 :
704 0 : close(iounit)
705 :
706 0 : return
707 :
708 : end subroutine write_fgong_data
709 :
710 : end module pulse_fgong
|
