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_cafein
21 :
22 : use star_private_def
23 : use const_def, only: dp, pi, pi4, crad, clight, msun, rsun, cgas
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_cafein_data
33 : public :: write_cafein_data
34 :
35 : integer, parameter :: NCOL = 35
36 :
37 : contains
38 :
39 0 : subroutine get_cafein_data (id, &
40 : add_center_point, keep_surface_point, add_atmosphere, global_data, point_data, ierr)
41 :
42 : integer, intent(in) :: id
43 : logical, intent(in) :: add_center_point
44 : logical, intent(in) :: keep_surface_point
45 : logical, intent(in) :: add_atmosphere
46 : real(dp), allocatable, intent(out) :: global_data(:)
47 : real(dp), allocatable, intent(out) :: point_data(:,:)
48 : integer, intent(out) :: ierr
49 :
50 : type(star_info), pointer :: s
51 : integer :: n_atm
52 : integer :: nn_atm
53 : integer :: n_env
54 : integer :: nn_env
55 : integer :: nn
56 0 : real(dp) :: R_star
57 0 : real(dp) :: M_star
58 0 : real(dp) :: P_c
59 0 : real(dp) :: T_c
60 0 : real(dp) :: rho_c
61 0 : real(dp) :: s_units
62 0 : real(dp) :: g_units
63 0 : real(dp) :: cm_units
64 0 : real(dp) :: erg_units
65 0 : real(dp) :: K_units
66 0 : real(dp), allocatable :: l_rad(:)
67 : integer :: j
68 : integer :: k
69 :
70 : ! Get model data for CAFein output
71 :
72 0 : call get_star_ptr(id, s, ierr)
73 0 : if (ierr /= 0) then
74 0 : write(*,*) 'bad star id for get_cafein_data'
75 0 : return
76 : end if
77 :
78 : ! Determine data dimensions
79 :
80 0 : if (add_atmosphere) then
81 0 : call build_atm(s, s%L(1), s%r(1), s%Teff, s%m_grav(1), s%cgrav(1), ierr)
82 0 : if (ierr /= 0) then
83 0 : write(*,*) 'failed in build_atm'
84 0 : return
85 : end if
86 0 : n_atm = s%atm_structure_num_pts
87 0 : nn_atm = n_atm - 1
88 : else
89 : n_atm = 0
90 : nn_atm = 0
91 : end if
92 :
93 0 : n_env = s%nz
94 :
95 0 : if (keep_surface_point) then
96 : nn_env = n_env
97 : else
98 0 : nn_env = n_env - 1
99 : end if
100 :
101 0 : if (add_center_point) then
102 0 : nn = nn_env + nn_atm + 1
103 : else
104 0 : nn = nn_env + nn_atm
105 : end if
106 :
107 : ! Store global data
108 :
109 0 : allocate(global_data(13))
110 :
111 0 : M_star = s%m_grav(1)
112 0 : R_star = Rsun*s%photosphere_r
113 :
114 0 : P_c = eval_center(s%rmid, s%Peos, 1, s%nz)
115 0 : T_c = eval_center(s%rmid, s%T, 1, s%nz)
116 0 : rho_c = eval_center(s%rmid, s%rho, 1, s%nz)
117 :
118 0 : s_units = SQRT(pow3(R_star)/(s%cgrav(1)*M_star))
119 0 : g_units = M_star
120 0 : cm_units = R_star
121 0 : erg_units = s%cgrav(1)*M_star**2/R_star
122 0 : K_units = s%cgrav(1)*M_star**2/(R_star*cgas)
123 :
124 0 : global_data(1) = M_star/Msun
125 0 : global_data(2) = R_star/Rsun
126 0 : global_data(3) = s%L(1)
127 0 : global_data(4) = P_c
128 0 : global_data(5) = T_c
129 0 : global_data(6) = rho_c
130 0 : global_data(7) = s%Teff
131 0 : global_data(8) = s_units
132 0 : global_data(9) = g_units
133 0 : global_data(10) = cm_units
134 0 : global_data(11) = erg_units
135 0 : global_data(12) = K_units
136 0 : global_data(13) = 0.d0
137 :
138 : ! Store point data. IMPORTANT NOTE: CAFein is ambiguous about
139 : ! which luminosity (total l or radiative l_rad) should be stored
140 : ! in output files (both as the Lr field, and for calculation of
141 : ! the dlnLr/dlnr, c_3 and c_4 fields). From context (i.e., reading
142 : ! the CAFein code), both dlnL/dlnr and the c_3/c_4 coefficients
143 : ! should be calculated using l_rad; however, it seems more useful
144 : ! for the Lr field to contain l. This is what the following code
145 : ! does.
146 :
147 0 : allocate(point_data(NCOL,nn))
148 :
149 0 : allocate(l_rad(nn))
150 :
151 0 : j = 1
152 :
153 : ! Atmosphere (we skip the point at the base of the atm to smooth
154 : ! the transition)
155 :
156 0 : atm_loop : do k = 1, n_atm-1
157 0 : call store_point_data_atm(j, k)
158 0 : j = j + 1
159 : end do atm_loop
160 :
161 : ! Envelope
162 :
163 0 : env_loop : do k = 1, n_env
164 :
165 0 : if (k == 1 .AND. .NOT. keep_surface_point) cycle env_loop
166 :
167 0 : call store_point_data_env(j, k)
168 0 : j = j + 1
169 :
170 : end do env_loop
171 :
172 : ! Center
173 :
174 0 : if (add_center_point) then
175 0 : call store_point_data_ctr(j)
176 0 : j = j + 1
177 : end if
178 :
179 : ! Check that all point data has correctly been calculated
180 :
181 0 : if (j /= nn+1) call mesa_error(__FILE__,__LINE__,'Invalid cell index in get_cafein_data')
182 :
183 : ! Reverse the data ordering (CAFein format is center-to-surface)
184 :
185 0 : point_data = point_data(:,nn:1:-1)
186 :
187 0 : l_rad = l_rad(nn:1:-1)
188 :
189 : ! Set remaining point data
190 :
191 : associate ( &
192 : r => point_data(5,:), &
193 : U => point_data(8,:), &
194 : N2 => point_data(10,:), &
195 : V => point_data(14,:), &
196 : nabla_ad => point_data(15,:), &
197 : c_2_fit => point_data(22,:), &
198 : dlnLr_dlnr => point_data(25,:), &
199 : surf_r_rad_beg => point_data(27,:), &
200 : surf_r_rad_end => point_data(28,:), &
201 : U_U_surf => point_data(34,:), &
202 : c_2 => point_data(35,:))
203 :
204 0 : c_2 = c_2 + nabla_ad*(log_deriv(r, nabla_ad, dy_a=0d0) + V)
205 0 : c_2_fit = c_2
206 :
207 0 : dlnLr_dlnr = log_deriv(r, l_rad, dy_a=3d0)
208 :
209 0 : U_U_surf = U/U(nn)
210 :
211 0 : surf_r_rad_beg = R_star
212 :
213 0 : do j = 1, nn
214 0 : surf_r_rad_end = r(j)
215 0 : if (N2(j) < 0d0) exit
216 : end do
217 :
218 : end associate
219 :
220 : ! Convert dimensioned data to dimensionless
221 :
222 : associate ( &
223 : m => point_data(1,:), &
224 : logrho => point_data(2,:) , &
225 : logP => point_data(3,:), &
226 : logT => point_data(4,:), &
227 : r => point_data(5,:), &
228 : N2 => point_data(10,:), &
229 : L2_ll1 => point_data(11,:), &
230 : g => point_data(12,:), &
231 : l => point_data(13,:), &
232 : C_P => point_data(17,:), &
233 : P_scale => point_data(26,:), &
234 : surf_r_rad_beg => point_data(27,:), &
235 : surf_r_rad_end => point_data(28,:), &
236 : entropy => point_data(30,:), &
237 : kap => point_data(31,:))
238 :
239 0 : m = m/g_units
240 :
241 0 : r = r/cm_units
242 0 : l = l/(erg_units/s_units)
243 :
244 0 : logrho = logrho - log10(g_units/cm_units**3)
245 0 : logP = logP - log10(erg_units/cm_units**3)
246 0 : logT = logT - log10(K_units)
247 :
248 0 : N2 = N2*s_units**2
249 0 : L2_ll1 = L2_ll1*s_units**2
250 :
251 0 : g = g/(cm_units/s_units**2)
252 :
253 0 : C_p = C_p/(erg_units/(K_units*g_units))
254 :
255 0 : P_scale = P_scale/cm_units
256 :
257 0 : surf_r_rad_beg = surf_r_rad_beg/cm_units
258 0 : surf_r_rad_end = surf_r_rad_end/cm_units
259 :
260 0 : entropy = entropy/(erg_units/(K_units*g_units))
261 :
262 0 : kap = kap/(cm_units**2/g_units)
263 :
264 : end associate
265 :
266 : ! Tidy up
267 :
268 0 : if (ASSOCIATED(s%atm_structure)) then
269 0 : deallocate(s%atm_structure)
270 : end if
271 :
272 0 : return
273 :
274 : contains
275 :
276 0 : subroutine store_point_data_atm (j, k)
277 :
278 : integer, intent(in) :: j
279 : integer, intent(in) :: k
280 :
281 0 : real(dp) :: rho
282 0 : real(dp) :: P
283 0 : real(dp) :: T
284 0 : real(dp) :: kap_rho
285 0 : real(dp) :: kap_T
286 0 : real(dp) :: kap_ad
287 0 : real(dp) :: eps
288 :
289 : ! Store data associated with atmosphere point k into the
290 : ! point_data array at position j
291 :
292 : associate ( &
293 : m => point_data(1,j), &
294 : logrho => point_data(2,j) , &
295 : logP => point_data(3,j), &
296 : logT => point_data(4,j), &
297 : r => point_data(5,j), &
298 : V_g => point_data(6,j), &
299 : As => point_data(7,j), &
300 : U => point_data(8,j), &
301 : c_1 => point_data(9,j), &
302 : N2 => point_data(10,j), &
303 : L2_ll1 => point_data(11,j), &
304 : g => point_data(12,j), &
305 : l => point_data(13,j), &
306 : V => point_data(14,j), &
307 : nabla_ad => point_data(15,j), &
308 : nabla => point_data(16,j), &
309 : C_P => point_data(17,j), &
310 : delta => point_data(18,j), &
311 : kap_S => point_data(19,j), &
312 : eps_ad => point_data(20,j), &
313 : eps_S => point_data(21,j), &
314 : c_3 => point_data(23,j), &
315 : c_4 => point_data(24,j), &
316 : P_scale => point_data(26,j), &
317 : Gamma_1 => point_data(29,j), &
318 : entropy => point_data(30,j), &
319 : kap => point_data(31,j), &
320 : chi_rho => point_data(32,j), &
321 : chi_T => point_data(33,j), &
322 : c_2 => point_data(35,j))
323 :
324 0 : m = s%m_grav(1) !+ s%atm_structure(atm_delta_m,k)
325 0 : r = s%r(1) + s%atm_structure(atm_delta_r,k)
326 0 : l = s%L(1)
327 :
328 0 : logrho = s%atm_structure(atm_lnd,k)/log(10d0)
329 0 : logP = s%atm_structure(atm_lnP,k)/log(10d0)
330 0 : logT = s%atm_structure(atm_lnT,k)/log(10d0)
331 :
332 0 : rho = pow(10d0, logrho)
333 0 : P = pow(10d0, logP)
334 0 : T = pow(10d0, logT)
335 :
336 0 : nabla_ad = s%atm_structure(atm_grada,k)
337 0 : nabla = s%atm_structure(atm_gradT,k)
338 0 : C_P = s%atm_structure(atm_cp,k)
339 0 : delta = s%atm_structure(atm_chiT,k)/s%atm_structure(atm_chiRho,k)
340 0 : Gamma_1 = s%atm_structure(atm_gamma1,k)
341 0 : entropy = 0d0 ! No entropy data in surface layers
342 0 : chi_rho = s%atm_structure(atm_chiRho,k)
343 0 : chi_T = s%atm_structure(atm_chiT,k)
344 :
345 0 : kap = s%atm_structure(atm_kap,k)
346 0 : kap_rho = s%atm_structure(atm_dlnkap_dlnd,k)
347 0 : kap_T = s%atm_structure(atm_dlnkap_dlnT,k)
348 0 : kap_ad = nabla_ad*kap_T + kap_rho/Gamma_1
349 0 : kap_S = kap_T - delta*kap_rho
350 :
351 0 : eps = 0d0
352 0 : eps_ad = 0d0
353 0 : eps_S = 0d0
354 :
355 0 : g = s%cgrav(1)*m/(r*r)
356 0 : N2 = g*g*(rho/P)*delta*(nabla_ad - nabla)
357 0 : L2_ll1 = Gamma_1*P/(rho*r**2)
358 0 : P_scale = P/(rho*g)
359 :
360 0 : V = rho*g*r/P
361 0 : V_g = V/Gamma_1
362 0 : As = N2*r/g
363 0 : U = pi4*rho*r**3/m
364 :
365 0 : l_rad(j) = l
366 :
367 0 : c_1 = (r/R_star)**3*(M_star/m)
368 0 : c_2 = (kap_ad - 4d0*nabla_ad)*V*nabla ! Note -- we omit the nabla_ad*(dnabla_ad + V) term for now
369 0 : c_3 = 0d0
370 0 : c_4 = pi4*r**3*rho*T*c_P/l_rad(j)*SQRT(s%cgrav(1)*M_star/R_star**3)
371 :
372 : end associate
373 :
374 :
375 0 : return
376 :
377 : end subroutine store_point_data_atm
378 :
379 :
380 0 : subroutine store_point_data_env (j, k)
381 :
382 : integer, intent(in) :: j
383 : integer, intent(in) :: k
384 :
385 0 : real(dp) :: rho
386 0 : real(dp) :: P
387 0 : real(dp) :: T
388 0 : real(dp) :: kap_rho
389 0 : real(dp) :: kap_T
390 0 : real(dp) :: kap_ad
391 0 : real(dp) :: eps
392 0 : real(dp) :: eps_rho
393 0 : real(dp) :: eps_T
394 :
395 : ! Store data associated with envelope face k into the point_data
396 : ! array at position j
397 :
398 : ! Store data associated with atmosphere point k into the
399 : ! point_data array at position j
400 :
401 : associate ( &
402 : m => point_data(1,j), &
403 : logrho => point_data(2,j) , &
404 : logP => point_data(3,j), &
405 : logT => point_data(4,j), &
406 : r => point_data(5,j), &
407 : V_g => point_data(6,j), &
408 : As => point_data(7,j), &
409 : U => point_data(8,j), &
410 : c_1 => point_data(9,j), &
411 : N2 => point_data(10,j), &
412 : L2_ll1 => point_data(11,j), &
413 : g => point_data(12,j), &
414 : l => point_data(13,j), &
415 : V => point_data(14,j), &
416 : nabla_ad => point_data(15,j), &
417 : nabla => point_data(16,j), &
418 : C_P => point_data(17,j), &
419 : delta => point_data(18,j), &
420 : kap_S => point_data(19,j), &
421 : eps_ad => point_data(20,j), &
422 : eps_S => point_data(21,j), &
423 : c_3 => point_data(23,j), &
424 : c_4 => point_data(24,j), &
425 : P_scale => point_data(26,j), &
426 : Gamma_1 => point_data(29,j), &
427 : entropy => point_data(30,j), &
428 : kap => point_data(31,j), &
429 : chi_rho => point_data(32,j), &
430 : chi_T => point_data(33,j), &
431 : c_2 => point_data(35,j))
432 :
433 0 : m = s%m_grav(k)
434 0 : r = s%r(k)
435 0 : l = s%L(k)
436 :
437 0 : if (s%interpolate_rho_for_pulse_data) then
438 0 : rho = eval_face(s%dq, s%rho, k, 1, s%nz)
439 : else
440 0 : rho = eval_face_rho(s, k, 1, s%nz)
441 : end if
442 0 : P = eval_face(s%dq, s%Peos, k, 1, s%nz)
443 0 : T = eval_face(s%dq, s%T, k, 1, s%nz)
444 :
445 0 : logrho = log10(rho)
446 0 : logP = log10(P)
447 0 : logT = log10(T)
448 :
449 0 : nabla_ad = eval_face(s%dq, s%grada, k, 1, s%nz)
450 0 : nabla = s%gradT(k)
451 0 : C_P = eval_face(s%dq, s%cp, k, 1, s%nz)
452 0 : delta = eval_face(s%dq, s%chiT, k, 1, s%nz)/eval_face(s%dq, s%chiRho, k, 1, s%nz)
453 0 : Gamma_1 = eval_face(s%dq, s%gamma1, k, 1, s%nz)
454 0 : entropy = exp(eval_face(s%dq, s%lnS, k, 1, s%nz))
455 0 : chi_rho = eval_face(s%dq, s%chiRho, k, 1, s%nz)
456 0 : chi_T = eval_face(s%dq, s%chiT, k, 1, s%nz)
457 :
458 0 : kap = eval_face(s%dq, s%opacity, k, 1, s%nz)
459 0 : kap_rho = eval_face(s%dq, s%d_opacity_dlnd, k, 1, s%nz)/kap
460 0 : kap_T = eval_face(s%dq, s%d_opacity_dlnT, k, 1, s%nz)/kap
461 0 : kap_ad = nabla_ad*kap_T + kap_rho/Gamma_1
462 0 : kap_S = kap_T - delta*kap_rho
463 :
464 0 : eps = eval_face(s%dq, s%eps_nuc, k, 1, s%nz)
465 0 : if (ABS(eps) > 1D-99) then
466 0 : eps_rho = eval_face(s%dq, s%d_epsnuc_dlnd, k, 1, s%nz)/eps
467 0 : eps_T = eval_face(s%dq, s%d_epsnuc_dlnT, k, 1, s%nz)/eps
468 : else
469 : eps_rho = 0d0
470 : eps_T = 0d0
471 : end if
472 0 : eps_ad = nabla_ad*eps_T + eps_rho/Gamma_1
473 0 : eps_S = eps_T - delta*eps_rho
474 :
475 0 : g = s%grav(k)
476 0 : N2 = eval_face_A_ast(s, k, 1, s%nz)*g/r
477 0 : L2_ll1 = Gamma_1*P/(rho*r**2)
478 0 : P_scale = s%scale_height(k)
479 :
480 0 : V = rho*g*r/P
481 0 : V_g = V/Gamma_1
482 0 : As = N2*r/g
483 0 : U = pi4*rho*r**3/m
484 :
485 0 : l_rad(j) = 16d0*pi*r*crad*clight*T**4*nabla*V/(3d0*kap*rho)
486 :
487 0 : c_1 = (r/R_star)**3*(M_star/m)
488 0 : c_2 = (kap_ad - 4d0*nabla_ad)*V*nabla ! Note -- we omit the nabla_ad*(dnabla_ad + V) term for now
489 0 : c_3 = pi4*r**3*rho*eps/l_rad(j)
490 0 : c_4 = pi4*r**3*rho*T*c_P/l_rad(j)*SQRT(s%cgrav(1)*M_star/R_star**3)
491 :
492 : end associate
493 :
494 :
495 0 : return
496 :
497 : end subroutine store_point_data_env
498 :
499 :
500 0 : subroutine store_point_data_ctr (j)
501 :
502 : integer, intent(in) :: j
503 :
504 0 : real(dp) :: rho
505 0 : real(dp) :: P
506 0 : real(dp) :: T
507 0 : real(dp) :: kap_rho
508 0 : real(dp) :: kap_T
509 0 : real(dp) :: kap_ad
510 0 : real(dp) :: eps
511 0 : real(dp) :: eps_rho
512 0 : real(dp) :: eps_T
513 0 : real(dp) :: cgrav
514 :
515 : ! Store data for the center into the point_data array at position j
516 :
517 : associate ( &
518 : m => point_data(1,j), &
519 : logrho => point_data(2,j) , &
520 : logP => point_data(3,j), &
521 : logT => point_data(4,j), &
522 : r => point_data(5,j), &
523 : V_g => point_data(6,j), &
524 : As => point_data(7,j), &
525 : U => point_data(8,j), &
526 : c_1 => point_data(9,j), &
527 : N2 => point_data(10,j), &
528 : L2_ll1 => point_data(11,j), &
529 : g => point_data(12,j), &
530 : l => point_data(13,j), &
531 : V => point_data(14,j), &
532 : nabla_ad => point_data(15,j), &
533 : nabla => point_data(16,j), &
534 : C_P => point_data(17,j), &
535 : delta => point_data(18,j), &
536 : kap_S => point_data(19,j), &
537 : eps_ad => point_data(20,j), &
538 : eps_S => point_data(21,j), &
539 : c_3 => point_data(23,j), &
540 : c_4 => point_data(24,j), &
541 : P_scale => point_data(26,j), &
542 : Gamma_1 => point_data(29,j), &
543 : entropy => point_data(30,j), &
544 : kap => point_data(31,j), &
545 : chi_rho => point_data(32,j), &
546 : chi_T => point_data(33,j), &
547 : c_2 => point_data(35,j))
548 :
549 0 : m = 0d0
550 0 : r = 0d0
551 0 : l = 0d0
552 :
553 0 : if (s%interpolate_rho_for_pulse_data) then
554 0 : rho = eval_center(s%rmid, s%rho, 1, s%nz)
555 : else
556 0 : rho = eval_center_rho(s, s%nz)
557 : end if
558 0 : P = eval_center(s%rmid, s%Peos, 1, s%nz)
559 0 : T = eval_center(s%rmid, s%T, 1, s%nz)
560 :
561 0 : logrho = log10(rho)
562 0 : logP = log10(P)
563 0 : logT = log10(T)
564 :
565 0 : nabla_ad = eval_center(s%rmid, s%grada, 1, s%nz)
566 0 : nabla = eval_center(s%r, s%gradT, 1, s%nz)
567 0 : C_P = eval_center(s%rmid, s%cp, 1, s%nz)
568 0 : delta = eval_center(s%rmid, s%chiT, 1, s%nz)/eval_center(s%rmid, s%chiRho, 1, s%nz)
569 0 : Gamma_1 = eval_center(s%rmid, s%gamma1, 1, s%nz)
570 0 : entropy = exp(eval_center(s%rmid, s%lnS, 1, s%nz))
571 0 : chi_rho = eval_center(s%rmid, s%chiRho, 1, s%nz)
572 0 : chi_T = eval_center(s%rmid, s%chiT, 1, s%nz)
573 :
574 0 : kap = eval_center(s%rmid, s%opacity, 1, s%nz)
575 0 : kap_rho = eval_center(s%rmid, s%d_opacity_dlnd, 1, s%nz)/kap
576 0 : kap_T = eval_center(s%rmid, s%d_opacity_dlnT, 1, s%nz)/kap
577 0 : kap_ad = nabla_ad*kap_T + kap_rho/Gamma_1
578 0 : kap_S = kap_T - delta*kap_rho
579 :
580 0 : eps = eval_center(s%rmid, s%eps_nuc, 1, s%nz)
581 0 : if (ABS(eps) > 1D-99) then
582 0 : eps_rho = eval_center(s%rmid, s%d_epsnuc_dlnd, 1, s%nz)/eps
583 0 : eps_T = eval_center(s%rmid, s%d_epsnuc_dlnT, 1, s%nz)/eps
584 : else
585 : eps_rho = 0d0
586 : eps_T = 0d0
587 : end if
588 0 : eps_ad = nabla_ad*eps_T + eps_rho/Gamma_1
589 0 : eps_S = eps_T - delta*eps_rho
590 :
591 0 : g = 0d0
592 0 : N2 = 0d0
593 0 : L2_ll1 = HUGE(0d0)
594 0 : P_scale = eval_center(s%r, s%scale_height, 1, s%nz)
595 :
596 0 : V = 0d0
597 0 : V_g = 0d0
598 0 : As = 0d0
599 0 : U = 3d0
600 :
601 0 : l_rad(j) = 0d0
602 :
603 0 : cgrav = eval_center(s%r, s%cgrav, 1, s%nz)
604 :
605 0 : c_1 = 3d0*(M_star/R_star**3)/(pi4*rho)
606 0 : c_2 = 0d0
607 0 : c_3 = 9d0*eps*kap*P/(16*pi*crad*clight*T**4*nabla*cgrav)
608 0 : c_4 = 9d0*T*c_P*kap*P/(16*pi*crad*clight*T**4*nabla*cgrav)*SQRT(s%cgrav(s%nz)*M_star/R_star**3)
609 :
610 : end associate
611 :
612 :
613 0 : return
614 :
615 : end subroutine store_point_data_ctr
616 :
617 :
618 0 : function log_deriv (x, y, dy_a, dy_b) result (dy)
619 :
620 : real(dp), intent(in) :: x(:)
621 : real(dp), intent(in) :: y(:)
622 : real(dp), intent(in), optional :: dy_a
623 : real(dp), intent(in), optional :: dy_b
624 : real(dp) :: dy(SIZE(x))
625 :
626 : integer :: n
627 : integer :: j
628 :
629 : ! Evaluate the logarithmic derivative of y wrt x, using simple
630 : ! finite differences
631 :
632 0 : n = SIZE(x)
633 :
634 0 : if (PRESENT(dy_a)) then
635 0 : dy(1) = dy_a
636 : else
637 0 : dy(1) = x(1)/y(1) * (y(2) - y(1))/(x(2) - x(1))
638 : end if
639 :
640 0 : do j = 2, n-1
641 0 : dy(j) = x(j)/y(j) * (y(j+1) - y(j-1))/(x(j+1) - x(j-1))
642 : end do
643 :
644 0 : if (PRESENT(dy_b)) then
645 0 : dy(n) = dy_b
646 : else
647 0 : dy(n) = x(n)/y(n) * (y(n) - y(n-1))/(x(n) - x(n-1))
648 : end if
649 :
650 0 : return
651 :
652 : end function log_deriv
653 :
654 : end subroutine get_cafein_data
655 :
656 :
657 0 : subroutine write_cafein_data (id, filename, global_data, point_data, ierr)
658 :
659 : integer, intent(in) :: id
660 : character(*), intent(in) :: filename
661 : real(dp), intent(in) :: global_data(:)
662 : real(dp), intent(in) :: point_data(:,:)
663 : integer, intent(out) :: ierr
664 :
665 : integer :: iounit
666 : integer :: nn
667 : integer :: j
668 :
669 : ! Write CAFein data to file
670 :
671 : ! Open the file
672 :
673 0 : open(newunit=iounit, file=TRIM(filename), status='REPLACE', iostat=ierr)
674 0 : if (ierr /= 0) then
675 0 : write(*,*) 'failed to open '//TRIM(filename)
676 0 : return
677 : end if
678 :
679 : ! Write the data
680 :
681 0 : nn = SIZE(point_data, 2)
682 :
683 : write(iounit, 100) &
684 : '..................M(Msun)' // '..................R(Rsun)' // '..................L(Lsun)' // &
685 : '..............Pc(erg/cm3)' // '....................Tc(K)' // '..............rhoC(g/cm3)' // &
686 : '..................Teff(K)' // '................sec_units' // '..............grams_units' // &
687 : '.................cm_units' // '................erg_units' // '..................K_units' // &
688 0 : '...................unused'
689 : 100 format(A)
690 :
691 0 : write(iounit, 110) global_data
692 : 110 format(13(1X,E24.14E3))
693 :
694 0 : write(iounit, *)
695 : write(iounit, 100) &
696 : '........................m' // '...................logrho' // '.....................logP' // &
697 : '.....................logT' // '........................r' // '.......................Vg' // &
698 : '....................Astar' // '........................U' // '.......................c1' // &
699 : '.......................N2' // '..............L2/(l(l+1))' // '........................g' // &
700 : '.......................Lr' // '........................V' // '............del_ad(noDim)' // &
701 : '...............del(noDim)' // '.......................Cp' // '......................v_t' // &
702 : '...........kappa_s(noDim)' // '............eps_ad(noDim)' // '.............eps_s(noDim)' // &
703 : '................c2_fitted' // '.......................c3' // '.......................c4' // &
704 : '...............dlnLr/dlnr' // '..................P_scale' // '........StartRadSurfLayer' // &
705 : '..........EndRadSurfLayer' // '...................Gamma1' // '..................entropy' // &
706 : '....................kappa' // '...................chiRho' // '.....................chiT' // &
707 0 : '..................U/Usurf' // '.......................c2'
708 :
709 0 : do j = 1, nn
710 0 : write(iounit, 120) point_data(:,j)
711 : 120 format(35(1X,E24.14E3))
712 : end do
713 :
714 : ! Close the file
715 :
716 0 : close(iounit)
717 :
718 0 : return
719 :
720 : end subroutine write_cafein_data
721 :
722 : end module pulse_cafein
|
