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 report
21 :
22 : use star_private_def
23 : use chem_def
24 : use utils_lib
25 : use star_utils
26 : use num_lib, only: find0
27 : use const_def, only: avo, kerg, pi, amu, clight, crad, Rsun, Lsun, Msun, &
28 : secday, secyer, ln10, mev_amu, ev2erg, one_third, two_thirds, four_thirds_pi, &
29 : no_mixing, convective_mixing, semiconvective_mixing
30 :
31 : implicit none
32 :
33 : contains
34 :
35 33 : subroutine set_min_gamma1(s)
36 : type (star_info), pointer :: s
37 : integer :: k
38 33 : real(dp) :: vesc
39 : include 'formats'
40 33 : s% min_gamma1 = 1d99
41 18256 : do k = s% nz, 1, -1
42 18256 : if (s% q(k) > s% gamma1_limit_max_q) exit
43 18223 : vesc = sqrt(2*s% cgrav(k)*s% m(k)/(s% r(k)))
44 18223 : if (s% u_flag) then
45 0 : if (s% u(k) > vesc*s% gamma1_limit_max_v_div_vesc) exit
46 18223 : else if (s% v_flag) then
47 0 : if (s% v(k) > vesc*s% gamma1_limit_max_v_div_vesc) exit
48 : end if
49 18256 : if (s% gamma1(k) < s% min_gamma1) s% min_gamma1 = s% gamma1(k)
50 : end do
51 33 : end subroutine set_min_gamma1
52 :
53 :
54 11 : subroutine set_power_info(s)
55 : use chem_def, only: category_name
56 : type (star_info), pointer :: s
57 : integer :: j, k, nz
58 11 : real(dp) :: eps_nuc
59 : include 'formats'
60 11 : nz = s% nz
61 :
62 275 : do j=1,num_categories
63 : s% L_by_category(j) = &
64 314352 : dot_product(s% dm(1:nz), s% eps_nuc_categories(j,1:nz))/Lsun
65 264 : s% center_eps_burn(j) = center_value_eps_burn(j)
66 275 : if (is_bad(s% L_by_category(j))) then
67 0 : do k=1,nz
68 0 : if (is_bad(s% eps_nuc_categories(j,k))) then
69 0 : write(*,2) trim(category_name(j)) // ' eps_nuc logT', k, s% eps_nuc_categories(j,k), s% lnT(k)/ln10
70 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'set_power_info')
71 : end if
72 : end do
73 : end if
74 : end do
75 :
76 11 : if (s% eps_nuc_factor == 0d0) then
77 0 : s% power_nuc_burn = 0d0
78 0 : s% power_nuc_neutrinos = 0d0
79 0 : s% power_nonnuc_neutrinos = 0d0
80 0 : s% power_neutrinos = 0d0
81 0 : s% power_h_burn = 0d0
82 0 : s% power_he_burn = 0d0
83 0 : s% power_z_burn = 0d0
84 0 : s% power_photo = 0d0
85 : else
86 : ! better if set power_nuc_burn using eps_nuc instead of categories
87 : ! categories can be subject to numerical jitters at very high temperatures
88 11 : s% power_nuc_burn = 0d0
89 13098 : do k=1,nz
90 13087 : if (s% op_split_burn .and. s% T_start(k) >= s% op_split_burn_min_T) then
91 0 : eps_nuc = s% burn_avg_epsnuc(k)
92 : else
93 13087 : eps_nuc = s% eps_nuc(k)
94 : end if
95 13098 : s% power_nuc_burn = s% power_nuc_burn + eps_nuc*s% dm(k)
96 : end do
97 11 : s% power_nuc_burn = s% power_nuc_burn/Lsun
98 13098 : s% power_nuc_neutrinos = dot_product(s% dm(1:nz),s% eps_nuc_neu_total(1:nz))/Lsun
99 11 : s% power_h_burn = s% L_by_category(ipp) + s% L_by_category(icno)
100 11 : s% power_he_burn = s% L_by_category(i3alf)
101 11 : s% power_z_burn = s% power_nuc_burn - (s% power_h_burn + s% power_he_burn)
102 11 : s% power_photo = s% L_by_category(iphoto)
103 : end if
104 :
105 11 : if (s% non_nuc_neu_factor == 0d0) then
106 0 : s% power_nonnuc_neutrinos = 0d0
107 : else
108 : s% power_nonnuc_neutrinos = &
109 13098 : dot_product(s% dm(1:nz),s% non_nuc_neu(1:nz))/Lsun
110 : end if
111 11 : s% power_neutrinos = s% power_nuc_neutrinos + s% power_nonnuc_neutrinos
112 11 : s% L_nuc_burn_total = s% power_nuc_burn
113 :
114 : contains
115 :
116 264 : real(dp) function center_value_eps_burn(j)
117 : integer, intent(in) :: j
118 264 : real(dp) :: sum_x, sum_dq, dx, dq
119 : integer :: k
120 264 : sum_x = 0
121 264 : sum_dq = 0
122 264 : do k = s% nz, 1, -1
123 264 : dq = s% dq(k)
124 264 : dx = s% eps_nuc_categories(j,k)*dq
125 264 : if (sum_dq+dq >= s% center_avg_value_dq) then
126 264 : sum_x = sum_x + dx*(s% center_avg_value_dq - sum_dq)/dq
127 264 : sum_dq = s% center_avg_value_dq
128 264 : exit
129 : end if
130 0 : sum_x = sum_x + dx
131 0 : sum_dq = sum_dq + dq
132 : end do
133 264 : center_value_eps_burn = sum_x/sum_dq
134 11 : end function center_value_eps_burn
135 :
136 : end subroutine set_power_info
137 :
138 :
139 33 : subroutine do_report(s, ierr)
140 : use rates_def, only: &
141 : i_rate, i_rate_dRho, i_rate_dT
142 : use star_utils, only: get_phot_info
143 : use hydro_rotation, only: set_surf_avg_rotation_info
144 : type (star_info), pointer :: s
145 : integer, intent(out) :: ierr
146 :
147 : integer :: k, nz, h1, h2, he3, he4, c12, n14, o16, ne20, si28, co56, ni56, k_min, k_fe
148 33 : real(dp) :: radius, dr, non_fe_core_mass, nu_for_delta_Pg, v, mstar, luminosity, mass_sum
149 33 : integer, pointer :: net_iso(:)
150 : real(dp), pointer :: velocity(:) => null()
151 :
152 : include 'formats'
153 :
154 33 : ierr = 0
155 33 : nz = s% nz
156 33 : net_iso => s% net_iso
157 :
158 33 : h1 = net_iso(ih1)
159 33 : h2 = net_iso(ih2)
160 33 : he3 = net_iso(ihe3)
161 33 : he4 = net_iso(ihe4)
162 33 : c12 = net_iso(ic12)
163 33 : n14 = net_iso(in14)
164 33 : o16 = net_iso(io16)
165 33 : ne20 = net_iso(ine20)
166 33 : si28 = net_iso(isi28)
167 33 : co56 = net_iso(ico56)
168 33 : ni56 = net_iso(ini56)
169 :
170 33 : radius = s% r(1) ! radius in cm
171 33 : s% log_surface_radius = log10(radius/Rsun)
172 : ! log10(stellar radius in solar units)
173 33 : s% log_center_density = center_value(s, s% lnd)/ln10
174 40799 : s% log_max_temperature = maxval(s% lnT(1:nz))/ln10
175 33 : s% log_center_temperature = center_value(s, s% lnT)/ln10
176 33 : s% log_center_pressure = center_value(s, s% lnPeos)/ln10
177 33 : s% center_degeneracy = center_value(s, s% eta)
178 :
179 33 : s% center_eps_nuc = center_value(s, s% eps_nuc)
180 :
181 33 : s% d_center_eps_nuc_dlnT = center_value(s, s% d_epsnuc_dlnT)
182 33 : s% d_center_eps_nuc_dlnd = center_value(s, s% d_epsnuc_dlnd)
183 33 : s% center_non_nuc_neu = center_value(s, s% non_nuc_neu)
184 :
185 33 : s% center_abar = center_value(s, s% abar)
186 33 : s% center_zbar = center_value(s, s% zbar)
187 33 : s% center_mu = center_value(s, s% mu)
188 33 : s% center_ye = center_value(s, s% ye)
189 33 : s% center_entropy = exp(center_value(s, s% lnS))*amu/kerg
190 40799 : s% max_entropy = exp(maxval(s% lnS(1:nz)))*amu/kerg
191 :
192 33 : if (.not. s% rotation_flag) then
193 40766 : s% omega(1:nz) = 0
194 33 : s% center_omega = 0
195 33 : s% center_omega_div_omega_crit = 0
196 : else
197 0 : s% center_omega = center_value(s, s% omega)
198 0 : s% center_omega_div_omega_crit = center_omega_div_omega_crit()
199 : end if
200 :
201 33 : luminosity = s% L(1)
202 :
203 33 : if (s% u_flag) then
204 0 : s% v_surf = s% u(1)
205 33 : else if (s% v_flag) then
206 0 : s% v_surf = s% v(1)
207 : else
208 33 : s% v_surf = 0d0
209 : end if
210 :
211 33 : call set_surf_avg_rotation_info(s)
212 :
213 33 : call set_min_gamma1(s)
214 :
215 : ! s% time is in seconds
216 33 : s% star_age = s% time/secyer
217 33 : s% time_years = s% time/secyer
218 33 : s% time_days = s% time/secday
219 33 : if ( s% model_number <= 0 ) then
220 1 : s% star_age = 0d0
221 1 : s% time_days = 0d0
222 1 : s% time_years = 0d0
223 : end if
224 :
225 : ! s% dt is in seconds
226 33 : s% time_step = s% dt/secyer ! timestep in years
227 33 : s% dt_years = s% dt/secyer
228 33 : s% dt_days = s% dt/secday
229 :
230 33 : mstar = s% mstar
231 33 : s% star_mass = mstar/Msun ! stellar mass in solar units
232 :
233 33 : s% kh_timescale = eval_kh_timescale(s% cgrav(1), mstar, radius, luminosity)/secyer
234 : ! kelvin-helmholtz timescale in years (about 1.6x10^7 for the sun)
235 :
236 33 : if (is_bad(s% kh_timescale)) then
237 0 : write(*,1) 's% kh_timescale', s% kh_timescale
238 0 : write(*,1) 's% cgrav(1)', s% cgrav(1)
239 0 : write(*,1) 'mstar', mstar
240 0 : write(*,1) 'radius', radius
241 0 : write(*,1) 'luminosity', luminosity
242 0 : stop
243 : end if
244 :
245 33 : if (luminosity > 0d0) then
246 33 : s% nuc_timescale = 1d10*s% star_mass/(luminosity/Lsun)
247 : else
248 0 : s% nuc_timescale = 1d99
249 : end if
250 : ! nuclear timescale in years (e.g., about 10^10 years for sun)
251 33 : if (s% cgrav(1) <= 0d0) then
252 0 : s% dynamic_timescale = 1d99
253 : else
254 33 : s% dynamic_timescale = 2*pi*sqrt(radius*radius*radius/(s% cgrav(1)*mstar))
255 : end if
256 :
257 : ! center_avg_x and surface_avg_x check if the species is not in the net
258 : ! and set's the values to 0 if so. So dont check species here.
259 33 : s% center_h1 = center_avg_x(s,h1)
260 33 : s% surface_h1 = surface_avg_x(s,h1)
261 33 : s% center_he3 = center_avg_x(s,he3)
262 33 : s% surface_he3 = surface_avg_x(s,he3)
263 33 : s% center_he4 = center_avg_x(s,he4)
264 33 : s% surface_he4 = surface_avg_x(s,he4)
265 33 : s% center_c12 = center_avg_x(s,c12)
266 33 : s% surface_c12 = surface_avg_x(s,c12)
267 33 : s% center_n14 = center_avg_x(s,n14)
268 33 : s% surface_n14 = surface_avg_x(s,n14)
269 33 : s% center_o16 = center_avg_x(s,o16)
270 33 : s% surface_o16 = surface_avg_x(s,o16)
271 33 : s% center_ne20 = center_avg_x(s,ne20)
272 33 : s% surface_ne20 = surface_avg_x(s,ne20)
273 33 : s% center_si28 = center_avg_x(s,si28)
274 :
275 : ! FYI profile stuff
276 40766 : do k=1,nz
277 :
278 40733 : s% entropy(k) = exp(s% lnS(k))/(avo*kerg)
279 40733 : if (is_bad(s% entropy(k))) then
280 0 : ierr = -1
281 0 : write(*,2) 'report: s% entropy(k)', k, s% entropy(k)
282 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'report')
283 0 : return
284 : end if
285 :
286 40733 : if (k == nz) then
287 33 : dr = s% r(k) - s% R_center
288 : else
289 40700 : dr = s% r(k) - s% r(k+1)
290 : end if
291 :
292 40733 : s% dr_div_csound(k) = dr/s% csound(k)
293 40733 : if (is_bad(s% dr_div_csound(k))) then
294 0 : ierr = -1
295 0 : write(*,2) 'report: s% dr_div_csound(k)', &
296 0 : k, s% dr_div_csound(k), dr, s% csound(k)
297 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'report')
298 0 : return
299 : end if
300 :
301 40733 : if (s% u_flag) then
302 0 : v = s% u(k)
303 40733 : else if (s% v_flag) then
304 0 : v = s% v(k)
305 : else
306 40733 : v = 0d0
307 : end if
308 :
309 40733 : if (is_bad(v)) then
310 0 : ierr = -1
311 0 : write(*,2) 'report: v', k, v
312 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'report')
313 0 : return
314 : end if
315 :
316 40733 : s% v_div_csound(k) = v/s% csound_face(k)
317 40766 : if (is_bad(s% v_div_csound(k))) then
318 0 : ierr = -1
319 0 : write(*,2) 'report: s% v_div_csound(k)', k, s% v_div_csound(k), &
320 0 : v, s% csound_face(k)
321 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'report')
322 0 : return
323 : end if
324 :
325 : end do
326 :
327 33 : call set_phot_info(s)
328 :
329 33 : if (s% photosphere_r*Rsun >= s% r(1)) then
330 : s% photosphere_acoustic_r = sum(s% dr_div_csound(1:nz)) + &
331 40766 : (s% photosphere_r*Rsun - s% r(1))/s% csound(1)
332 : else
333 0 : do k=2,nz
334 0 : if (s% photosphere_r*Rsun > s% r(k)) then
335 : s% photosphere_acoustic_r = sum(s% dr_div_csound(k:nz)) + &
336 0 : (s% photosphere_r*Rsun - s% r(k))/s% csound(k-1)
337 0 : exit
338 : end if
339 : end do
340 : end if
341 :
342 33 : if (.not. s% get_delta_nu_from_scaled_solar) then
343 33 : s% delta_nu = 1d6/(2*s% photosphere_acoustic_r) ! microHz
344 : else
345 : s% delta_nu = &
346 : s% delta_nu_sun*sqrt(s% star_mass)*pow3(s% Teff/s% astero_Teff_sun) / &
347 0 : pow(s% L_phot,0.75d0)
348 : end if
349 :
350 33 : call get_mass_info(s, s% dm, ierr)
351 33 : if (failed('get_mass_info')) return
352 :
353 : s% nu_max = s% nu_max_sun*s% star_mass/ &
354 33 : (pow2(s% photosphere_r)*sqrt(max(0d0,s% Teff)/s% astero_Teff_sun))
355 : s% acoustic_cutoff = &
356 33 : 0.25d6/pi*s% grav(1)*sqrt(s% gamma1(1)*s% rho(1)/s% Peos(1))
357 33 : nu_for_delta_Pg = s% nu_max
358 33 : if (s% delta_Pg_mode_freq > 0) nu_for_delta_Pg = s% delta_Pg_mode_freq
359 33 : if ( .not. s% delta_Pg_traditional) then
360 0 : call get_delta_Pg_bildsten2012(s, nu_for_delta_Pg, s% delta_Pg)
361 : else
362 33 : call get_delta_Pg_traditional(s, s% delta_Pg)
363 : end if
364 :
365 33 : if (s% rsp_flag) return
366 :
367 33 : call get_mixing_regions(s, ierr)
368 33 : if (failed('get_mixing_regions')) return
369 :
370 33 : call set_mass_conv_core
371 33 : call set_mass_semiconv_core
372 33 : call find_conv_mx_regions
373 33 : call find_mx_regions
374 33 : call get_burn_zone_info(s, ierr)
375 33 : if (failed('get_burn_zone_info')) return
376 :
377 :
378 33 : s% fe_core_infall = 0d0
379 33 : s% non_fe_core_infall = 0d0
380 33 : s% non_fe_core_rebound = 0d0
381 33 : s% max_infall_speed_mass = 0d0
382 33 : k_fe = -1 ! initialize, in case u_flag & v_flag = .false.
383 33 : k_min = -1 ! initialize, in case u_flag & v_flag = .false.
384 :
385 66 : if(s% u_flag .or. s% v_flag) then
386 :
387 0 : if(s% u_flag) then
388 0 : velocity => s% u
389 : else
390 0 : velocity => s% v
391 : end if
392 0 : k_min = minloc(velocity(1:nz), dim=1)
393 :
394 0 : s% max_infall_speed_mass = s% m(k_min)/Msun
395 :
396 0 : mass_sum = 0d0
397 0 : if (s% fe_core_mass > 0) then
398 : ! check if [> fe_core_infall_mass] of core is infalling
399 : ! instead of a for loop, we move inside out and only check regions inside the fe_core
400 0 : k = nz
401 0 : do while (k > 1)
402 0 : if (s% m(k) > Msun * s% fe_core_mass) exit ! exit when outside Fe core
403 0 : if (-velocity(k) > s% fe_core_infall) then
404 0 : mass_sum = mass_sum + s% dm(k)
405 : end if
406 0 : k = k-1 ! loop outwards
407 : end do
408 0 : k_fe = k ! mark fe_core_mass boundary cell location, k_fe = nz if no fe_core
409 :
410 0 : if (s% report_max_infall_inside_fe_core) then ! report peak infall velocity inside fe_core (not necessarily the maximum, since infall tends to start outside in)
411 0 : if (mass_sum > s% fe_core_infall_mass*msun) then ! prevents toggling when k == nz
412 0 : s% fe_core_infall = - minval(s%v(k_fe:nz))
413 : end if
414 : else !(default in r24.03.1 prior) report max infall velocity anywhere
415 0 : if(mass_sum > s% fe_core_infall_mass*msun) then
416 0 : s% fe_core_infall = -velocity(k_min)
417 : end if
418 : end if
419 : end if
420 :
421 0 : non_fe_core_mass = s% he_core_mass
422 0 : mass_sum = 0d0
423 :
424 0 : if (non_fe_core_mass > 0) then
425 0 : do k=1, nz
426 0 : if (s% m(k) > Msun * non_fe_core_mass) cycle
427 0 : if (s% m(k) < Msun * s% fe_core_mass) exit
428 0 : if(-velocity(k) > s% non_fe_core_infall) mass_sum = mass_sum + s% dm(k)
429 : end do
430 :
431 0 : if ((mass_sum > s% non_fe_core_infall_mass*msun) .and. &
432 : (s%m(k_min) <= s% he_core_mass * msun)) then
433 0 : s% non_fe_core_infall = -velocity(k_min)
434 : end if
435 :
436 0 : s% non_fe_core_rebound = maxval(velocity(s%he_core_k:s%fe_core_k),dim=1)
437 :
438 : end if
439 :
440 0 : nullify(velocity)
441 : end if
442 :
443 :
444 : contains
445 :
446 :
447 99 : logical function failed(str)
448 : character (len=*), intent(in) :: str
449 99 : failed = (ierr /= 0)
450 99 : if (failed) then
451 0 : write(*, *) trim(str) // ' ierr', ierr
452 : end if
453 33 : end function failed
454 :
455 :
456 33 : subroutine set_mass_conv_core
457 : integer :: j, nz, k
458 33 : real(dp) :: dm_limit
459 : include 'formats'
460 33 : s% mass_conv_core = 0
461 33 : dm_limit = s% conv_core_gap_dq_limit*s% xmstar
462 33 : nz = s% nz
463 33 : do j = 1, s% n_conv_regions
464 : ! ignore possible small gap at center
465 33 : if (s% cz_bot_mass(j) <= s% m(nz) + dm_limit) then
466 32 : s% mass_conv_core = s% cz_top_mass(j)/Msun
467 : ! jump over small gaps
468 32 : do k = j+1, s% n_conv_regions
469 32 : if (s% cz_bot_mass(k) - s% cz_top_mass(k-1) >= dm_limit) exit
470 32 : s% mass_conv_core = s% cz_top_mass(k)/Msun
471 : end do
472 : exit
473 : end if
474 : end do
475 33 : end subroutine set_mass_conv_core
476 :
477 :
478 33 : subroutine set_mass_semiconv_core
479 : integer :: k, ktop, nz
480 33 : real(dp) :: qb
481 : include 'formats'
482 33 : s% mass_semiconv_core = s% mass_conv_core
483 33 : nz = s% nz
484 33 : if (s% mixing_type(nz) /= convective_mixing) return
485 33 : ktop = 1
486 33 : do k=nz-1,1,-1
487 33 : if (s% mixing_type(k) == convective_mixing) then
488 : ktop = k
489 : exit
490 : end if
491 : end do
492 33 : if (ktop == 1 .or. s% mixing_type(ktop) /= semiconvective_mixing) return
493 0 : do k=ktop-1,1,-1
494 0 : if (s% mixing_type(k) /= semiconvective_mixing) then
495 0 : qb = s% q(k) - s% cz_bdy_dq(k)
496 0 : s% mass_semiconv_core = (qb*s% xmstar + s% M_center)/Msun
497 0 : return
498 : end if
499 : end do
500 0 : s% mass_semiconv_core = s% star_mass
501 : end subroutine set_mass_semiconv_core
502 :
503 :
504 : real(dp) function volume_at_q(q)
505 : use interp_1d_def
506 : use interp_1d_lib
507 : real(dp), intent(in) :: q
508 : integer, parameter :: n_old = 4, n_new = 1, nwork = pm_work_size
509 : real(dp) :: qlo, x_old(n_old), v_old(n_old), x_new(n_new), v_new(n_new)
510 : integer :: k, nz, k00, ierr
511 : real(dp), target :: work_ary(n_old*nwork)
512 : real(dp), pointer :: work(:)
513 : work => work_ary
514 :
515 : include 'formats'
516 : nz = s% nz
517 :
518 : if (q == 0d0) then
519 : volume_at_q = 0
520 : return
521 : end if
522 :
523 : if (q <= s% q(nz)) then
524 : volume_at_q = (q/s% q(nz))*four_thirds_pi*s% r(nz)*s% r(nz)*s% r(nz)
525 : return
526 : end if
527 : k00 = 1
528 : do k=nz-1,2,-1
529 : if (s% q(k) >= q) then
530 : k00 = k; exit
531 : end if
532 : end do
533 : if (k00 == 1) then
534 : volume_at_q = four_thirds_pi*q*s% r(1)*s% r(1)*s% r(1)
535 : write(*,1) 'volume_at_q', volume_at_q
536 : return
537 : end if
538 :
539 : x_old(1) = 0
540 : if (k00+1 == nz) then
541 : v_old(1) = s% R_center*s% R_center*s% R_center
542 : qlo = 0
543 : else
544 : v_old(1) = s% r(k00+2)*s% r(k00+2)*s% r(k00+2)
545 : qlo = s% q(k00+2)
546 : end if
547 :
548 : x_old(2) = s% dq(k00+1)
549 : v_old(2) = s% r(k00+1)*s% r(k00+1)*s% r(k00+1)
550 :
551 : x_old(3) = x_old(2) + s% dq(k00)
552 : v_old(3) = s% r(k00)*s% r(k00)*s% r(k00)
553 :
554 : x_old(4) = x_old(3) + s% dq(k00-1)
555 : v_old(4) = s% r(k00-1)*s% r(k00-1)*s% r(k00-1)
556 :
557 : ierr = 0
558 : x_new(1) = q-qlo
559 : call interpolate_vector( &
560 : n_old, x_old, n_new, x_new, v_old, v_new, interp_pm, nwork, work, &
561 : 'report volume_at_q', ierr)
562 : if (ierr /= 0) then
563 : write(*,*) 'volume_at_q: failed in interpolate_vector'
564 : volume_at_q = (v_old(2) + v_old(3))/2
565 : return
566 : end if
567 :
568 : volume_at_q = four_thirds_pi*v_new(1)
569 :
570 : end function volume_at_q
571 :
572 :
573 0 : real(dp) function center_omega_div_omega_crit()
574 0 : real(dp) :: sum_x, sum_dq, dx, dq
575 : integer :: k
576 0 : center_omega_div_omega_crit = 0
577 0 : if (.not. s% rotation_flag) return
578 0 : sum_x = 0
579 0 : sum_dq = 0
580 0 : do k = s% nz, 1, -1
581 0 : dq = s% dq(k)
582 0 : dx = dq*s% omega(k)/omega_crit(s,k)
583 0 : if (sum_dq+dq >= s% center_avg_value_dq) then
584 0 : sum_x = sum_x+ dx*(s% center_avg_value_dq - sum_dq)/dq
585 0 : sum_dq = s% center_avg_value_dq
586 0 : exit
587 : end if
588 0 : sum_x = sum_x + dx
589 0 : sum_dq = sum_dq + dq
590 : end do
591 0 : center_omega_div_omega_crit = min(1d0,sum_x/sum_dq)
592 0 : end function center_omega_div_omega_crit
593 :
594 :
595 33 : subroutine find_conv_mx_regions
596 33 : real(dp) :: conv_mx1_dq, conv_mx2_dq, mx_dq
597 : integer :: i, ktop, kbot, conv_mx1, conv_mx2
598 :
599 : include 'formats'
600 :
601 33 : s% largest_conv_mixing_region = 0
602 :
603 33 : s% conv_mx1_top = 0
604 33 : s% conv_mx1_bot = 0
605 33 : s% conv_mx2_top = 0
606 33 : s% conv_mx2_bot = 0
607 :
608 33 : s% conv_mx1_top_r = 0
609 33 : s% conv_mx1_bot_r = 0
610 33 : s% conv_mx2_top_r = 0
611 33 : s% conv_mx2_bot_r = 0
612 :
613 33 : if (s% num_mixing_regions == 0) return
614 :
615 : conv_mx1 = 0
616 : conv_mx2 = 0
617 : conv_mx1_dq = 0
618 : conv_mx2_dq = 0
619 :
620 102 : do i = 1, s% num_mixing_regions
621 69 : if (s% mixing_region_type(i) /= convective_mixing) cycle
622 : ! mx_dq = s% q(s% mixing_region_top(i)) - s% q(s% mixing_region_bottom(i))
623 6070 : mx_dq = sum(s% dq(s% mixing_region_top(i):s% mixing_region_bottom(i)))
624 102 : if (mx_dq > conv_mx1_dq) then
625 : conv_mx2_dq = conv_mx1_dq; conv_mx1_dq = mx_dq
626 : conv_mx2 = conv_mx1; conv_mx1 = i
627 0 : else if (mx_dq > conv_mx2_dq) then
628 69 : conv_mx2_dq = mx_dq
629 69 : conv_mx2 = i
630 : end if
631 : end do
632 :
633 33 : if (conv_mx1 > 0) then
634 33 : s% largest_conv_mixing_region = conv_mx1
635 33 : ktop = s% mixing_region_top(conv_mx1)
636 33 : kbot = s% mixing_region_bottom(conv_mx1)
637 33 : s% conv_mx1_top = s% q(ktop) - s% cz_bdy_dq(ktop)
638 33 : s% conv_mx1_bot = s% q(kbot) - s% cz_bdy_dq(kbot)
639 33 : s% conv_mx1_top_r = s% r(ktop)/Rsun
640 33 : s% conv_mx1_bot_r = s% r(kbot)/Rsun
641 : end if
642 :
643 33 : if (conv_mx2 > 0) then
644 33 : ktop = s% mixing_region_top(conv_mx2)
645 33 : kbot = s% mixing_region_bottom(conv_mx2)
646 33 : s% conv_mx2_top = s% q(ktop) - s% cz_bdy_dq(ktop)
647 33 : s% conv_mx2_bot = s% q(kbot) - s% cz_bdy_dq(kbot)
648 33 : s% conv_mx2_top_r = s% r(ktop)/Rsun
649 33 : s% conv_mx2_bot_r = s% r(kbot)/Rsun
650 : end if
651 :
652 : end subroutine find_conv_mx_regions
653 :
654 :
655 33 : subroutine find_mx_regions
656 33 : real(dp) :: mx1_dq, mx2_dq, mx_dq
657 : integer :: i, ktop, kbot, mx1_top_region, mx1_bottom_region, &
658 : mx2_top_region, mx2_bottom_region, &
659 : current_top_region, current_bottom_region, &
660 : current_top_point, current_bottom_point
661 :
662 : include 'formats'
663 :
664 33 : s% mx1_top = 0
665 33 : s% mx1_bot = 0
666 33 : s% mx2_top = 0
667 33 : s% mx2_bot = 0
668 :
669 33 : s% mx1_top_r = 0
670 33 : s% mx1_bot_r = 0
671 33 : s% mx2_top_r = 0
672 33 : s% mx2_bot_r = 0
673 :
674 33 : if (s% num_mixing_regions == 0) return
675 :
676 : mx1_top_region = 0
677 : mx1_bottom_region = 0
678 : mx1_dq = 0
679 :
680 : mx2_top_region = 0
681 : mx2_bottom_region = 0
682 : mx2_dq = 0
683 :
684 : i = 1
685 : do
686 102 : if (i > s% num_mixing_regions) exit
687 69 : current_top_region = i
688 69 : current_top_point = s% mixing_region_top(current_top_region)
689 : do
690 69 : current_bottom_region = i
691 69 : current_bottom_point = s% mixing_region_bottom(current_bottom_region)
692 69 : i = i+1
693 69 : if (i > s% num_mixing_regions) exit
694 69 : if (s% mixing_region_top(i) /= current_bottom_point+1) exit
695 : end do
696 : ! mx_dq = s% q(current_top_point) - s% q(current_bottom_point)
697 6070 : mx_dq = sum(s% dq(current_top_point:current_bottom_point))
698 102 : if (mx_dq > mx1_dq) then
699 : mx2_dq = mx1_dq; mx1_dq = mx_dq
700 : mx2_top_region = mx1_top_region
701 : mx1_top_region = current_top_region
702 : mx2_bottom_region = mx1_bottom_region
703 : mx1_bottom_region = current_bottom_region
704 0 : else if (mx_dq > mx2_dq) then
705 0 : mx2_dq = mx_dq
706 0 : mx2_top_region = current_top_region
707 0 : mx2_bottom_region = current_bottom_region
708 : end if
709 : end do
710 :
711 33 : if (mx1_top_region > 0) then
712 33 : ktop = s% mixing_region_top(mx1_top_region)
713 33 : kbot = s% mixing_region_bottom(mx1_bottom_region)
714 33 : s% mx1_top = s% q(ktop) - s% cz_bdy_dq(ktop)
715 33 : s% mx1_bot = s% q(kbot) - s% cz_bdy_dq(kbot)
716 33 : s% mx1_top_r = s% r(ktop)/Rsun
717 33 : s% mx1_bot_r = s% r(kbot)/Rsun
718 : end if
719 :
720 33 : if (mx2_top_region > 0) then
721 33 : ktop = s% mixing_region_top(mx2_top_region)
722 33 : kbot = s% mixing_region_bottom(mx2_bottom_region)
723 33 : s% mx2_top = s% q(ktop) - s% cz_bdy_dq(ktop)
724 33 : s% mx2_bot = s% q(kbot) - s% cz_bdy_dq(kbot)
725 33 : s% mx2_top_r = s% r(ktop)/Rsun
726 33 : s% mx2_bot_r = s% r(kbot)/Rsun
727 : end if
728 :
729 : end subroutine find_mx_regions
730 :
731 :
732 : end subroutine do_report
733 :
734 :
735 33 : subroutine get_burn_zone_info(s, ierr)
736 : type (star_info), pointer :: s
737 : integer, intent(out) :: ierr
738 33 : real(dp) :: burn_min1, burn_min2
739 : integer :: i, i_start
740 : include 'formats'
741 33 : burn_min1 = s% burn_min1; burn_min2 = s% burn_min2
742 : ! up to 3 zones where eps_nuc > burn_min1 erg/g/s
743 : ! for each zone have 4 numbers: start1, start2, end2, end1
744 : ! start1 is mass of inner edge where first goes > burn_min1 (or -20 if none such)
745 : ! start2 is mass of inner edge where first zone reaches burn_min2 erg/g/sec (or -20 if none such)
746 : ! end2 is mass of outer edge where first zone drops back below burn_min2 erg/g/s
747 : ! end1 is mass of outer edge where first zone ends (i.e. eps_nuc < burn_min1)
748 : ! similar for second and third zones
749 33 : i_start = s% nz
750 132 : do i=1,3
751 : call find_epsnuc_zone(s, i_start, &
752 : s% burn_zone_mass(1,i), s% burn_zone_mass(2,i), &
753 : s% burn_zone_mass(3,i), s% burn_zone_mass(4,i), &
754 132 : s% burn_min1, s% burn_min2, ierr)
755 : end do
756 33 : end subroutine get_burn_zone_info
757 :
758 :
759 99 : subroutine find_epsnuc_zone( &
760 : s, i_start, bzm_1, bzm_2, bzm_3, bzm_4, burn_min1, burn_min2, ierr)
761 : type (star_info), pointer :: s
762 : integer, intent(inout) :: i_start
763 : real(dp), intent(out) :: bzm_1, bzm_2, bzm_3, bzm_4
764 : real(dp), intent(in) :: burn_min1, burn_min2
765 : integer, intent(out) :: ierr
766 :
767 : real(dp), parameter :: null_zone = -20
768 : integer :: i, burn_zone
769 99 : real(dp) :: prev_m, prev_x, cur_m, cur_x
770 99 : ierr = 0
771 99 : bzm_1 = null_zone; bzm_2 = null_zone; bzm_3 = null_zone; bzm_4 = null_zone
772 99 : burn_zone = 0 ! haven't entered the zone yet
773 99 : if (i_start /= s% nz) then
774 66 : i = i_start+1
775 66 : prev_m = s% m(i)
776 66 : prev_x = s% eps_nuc(i)
777 : else ! keep the compiler happy
778 33 : prev_m = 0
779 33 : prev_x = 0
780 : end if
781 40832 : do i = i_start, 1, -1
782 40766 : cur_m = s% m(i)
783 40766 : cur_x = s% eps_nuc(i)
784 36981 : select case (burn_zone)
785 : case (0)
786 36981 : if ( cur_x > burn_min2 ) then
787 33 : if ( i == s% nz ) then ! use star center as start of zone
788 33 : bzm_2 = 0d0
789 : else ! interpolate to estimate where rate reached burn_min1
790 0 : bzm_2 = find0(prev_m, prev_x-burn_min2, cur_m, cur_x-burn_min2)
791 : end if
792 33 : bzm_1 = bzm_2
793 33 : burn_zone = 2
794 36948 : elseif ( cur_x > burn_min1 ) then
795 0 : if ( i == s% nz ) then ! use star center as start of zone
796 0 : bzm_1 = 0d0
797 : else ! interpolate to estimate where rate reached burn_min1
798 0 : bzm_1 = find0(prev_m, prev_x-burn_min1, cur_m, cur_x-burn_min1)
799 : end if
800 : burn_zone = 1
801 : end if
802 : case (1) ! in the initial eps > burn_min1 region
803 0 : if ( cur_x > burn_min2 ) then
804 0 : bzm_2 = find0(prev_m, prev_x-burn_min2, cur_m, cur_x-burn_min2)
805 0 : burn_zone = 2
806 0 : else if ( cur_x < burn_min1 ) then
807 0 : bzm_4 = find0(prev_m, prev_x-burn_min1, cur_m, cur_x-burn_min1)
808 0 : i_start = i
809 99 : return
810 : end if
811 : case (2) ! in the initial eps > burn_min2 region
812 1429 : if ( cur_x < burn_min1 ) then
813 0 : bzm_4 = find0(prev_m, prev_x-burn_min1, cur_m, cur_x-burn_min1)
814 0 : bzm_3 = bzm_4
815 0 : i_start = i
816 0 : return
817 : end if
818 1429 : if ( cur_x < burn_min2 ) then
819 33 : bzm_3 = find0(prev_m, prev_x-burn_min2, cur_m, cur_x-burn_min2)
820 33 : burn_zone = 3
821 : end if
822 : case (3) ! in the final eps > burn_min1 region
823 2356 : if ( cur_x < burn_min1 ) then
824 33 : bzm_4 = find0(prev_m, prev_x-burn_min1, cur_m, cur_x-burn_min1)
825 33 : i_start = i
826 33 : return
827 : end if
828 : case default
829 : ierr = -1
830 : write(*,*) 'error in find_eps_nuc_zone'
831 40766 : return
832 : end select
833 40799 : prev_m = cur_m; prev_x = cur_x
834 : end do
835 66 : i_start = 0
836 : select case (burn_zone)
837 : case (0)
838 0 : return
839 : case (1)
840 0 : bzm_4 = cur_m
841 : case (2)
842 0 : bzm_3 = cur_m
843 0 : bzm_4 = cur_m
844 : case (3)
845 0 : bzm_4 = cur_m
846 : case default
847 : ierr = -1
848 : write(*,*) 'error in find_eps_nuc_zone'
849 66 : return
850 : end select
851 : end subroutine find_epsnuc_zone
852 :
853 :
854 33 : subroutine get_mass_info(s, cell_masses, ierr)
855 : type (star_info), pointer :: s
856 : real(dp), pointer, intent(in) :: cell_masses(:)
857 : integer, intent(out) :: ierr
858 :
859 : integer :: k, nz, j
860 33 : real(dp) :: cell_mass
861 33 : integer, pointer :: net_iso(:)
862 :
863 : include 'formats'
864 :
865 33 : ierr = 0
866 33 : nz = s% nz
867 33 : net_iso => s% net_iso
868 :
869 33 : s% star_mass_h1=0d0
870 33 : s% star_mass_he3=0d0
871 33 : s% star_mass_he4=0d0
872 33 : s% star_mass_c12 = 0d0
873 33 : s% star_mass_n14 = 0d0
874 33 : s% star_mass_o16 = 0d0
875 33 : s% star_mass_ne20 = 0d0
876 :
877 40766 : do k = 1, nz
878 40733 : cell_mass = cell_masses(k)
879 366630 : do j=1, s% species
880 366597 : if (s% chem_id(j) == ih1) then
881 40733 : s% star_mass_h1 = s% star_mass_h1 + cell_mass*s% xa(j, k)
882 285131 : else if (s% chem_id(j) == ihe3) then
883 40733 : s% star_mass_he3 = s% star_mass_he3 + cell_mass*s% xa(j, k)
884 244398 : else if (s% chem_id(j) == ihe4) then
885 40733 : s% star_mass_he4 = s% star_mass_he4 + cell_mass*s% xa(j, k)
886 203665 : else if (s% chem_id(j) == ic12) then
887 40733 : s% star_mass_c12 = s% star_mass_c12 + cell_mass*s% xa(j, k)
888 162932 : else if (s% chem_id(j) == in14) then
889 40733 : s% star_mass_n14 = s% star_mass_n14 + cell_mass*s% xa(j, k)
890 122199 : else if (s% chem_id(j) == io16) then
891 40733 : s% star_mass_o16 = s% star_mass_o16 + cell_mass*s% xa(j, k)
892 81466 : else if (s% chem_id(j) == ine20) then
893 40733 : s% star_mass_ne20 = s% star_mass_ne20 + cell_mass*s% xa(j, k)
894 : end if
895 : end do
896 : end do
897 :
898 33 : s% star_mass_h1 = s% star_mass_h1 / Msun
899 33 : s% star_mass_he3 = s% star_mass_he3 / Msun
900 33 : s% star_mass_he4 = s% star_mass_he4 / Msun
901 33 : s% star_mass_c12 = s% star_mass_c12 / Msun
902 33 : s% star_mass_n14 = s% star_mass_n14 / Msun
903 33 : s% star_mass_o16 = s% star_mass_o16 / Msun
904 33 : s% star_mass_ne20 = s% star_mass_ne20 / Msun
905 :
906 33 : call get_core_info(s)
907 33 : call get_shock_info(s)
908 :
909 33 : end subroutine get_mass_info
910 :
911 :
912 0 : subroutine get_info_at_surface( &
913 : s, bdy_m, bdy_r, bdy_lgT, bdy_lgRho, bdy_L, bdy_v, bdy_time, &
914 : bdy_omega, bdy_omega_div_omega_crit)
915 : type (star_info), pointer :: s
916 : real(dp), intent(out) :: &
917 : bdy_m, bdy_r, bdy_lgT, bdy_lgRho, bdy_L, bdy_v, &
918 : bdy_omega, bdy_omega_div_omega_crit, bdy_time
919 :
920 0 : real(dp) :: bdy_omega_crit
921 :
922 0 : bdy_time = 0
923 0 : bdy_m = s% star_mass
924 0 : bdy_r = s% r(1)/Rsun
925 0 : bdy_lgT = s% lnT(1)/ln10
926 0 : bdy_lgRho = s% lnd(1)/ln10
927 0 : bdy_L = s% L(1)/Lsun
928 0 : if (s% v_flag) then
929 0 : bdy_v = s% v(1)
930 0 : else if (s% u_flag) then
931 0 : bdy_v = s% u_face_ad(1)%val
932 : else
933 0 : bdy_v = 0d0
934 : end if
935 0 : if (s% rotation_flag) then
936 0 : bdy_omega = s% omega_avg_surf
937 0 : bdy_omega_crit = s% omega_crit_avg_surf
938 0 : bdy_omega_div_omega_crit = s% w_div_w_crit_avg_surf
939 : else
940 0 : bdy_omega = 0
941 0 : bdy_omega_div_omega_crit = 0
942 : end if
943 :
944 0 : end subroutine get_info_at_surface
945 :
946 :
947 0 : subroutine get_mach1_location_info( &
948 : s, dbg, k_shock, v, r, &
949 : mach1_mass, &
950 : mach1_q, &
951 : mach1_radius, &
952 : mach1_velocity, &
953 : mach1_csound, &
954 : mach1_lgT, &
955 : mach1_lgRho, &
956 : mach1_lgP, &
957 : mach1_gamma1, &
958 : mach1_entropy, &
959 : mach1_tau, &
960 : mach1_k)
961 : type (star_info), pointer :: s
962 : integer, intent(in) :: k_shock
963 : logical, intent(in) :: dbg
964 : real(dp), intent(in), pointer :: v(:)
965 : real(dp), intent(in) :: r
966 : real(dp), intent(out) :: &
967 : mach1_mass, &
968 : mach1_q, &
969 : mach1_radius, &
970 : mach1_velocity, &
971 : mach1_csound, &
972 : mach1_lgT, &
973 : mach1_lgRho, &
974 : mach1_lgP, &
975 : mach1_gamma1, &
976 : mach1_entropy, &
977 : mach1_tau
978 : integer, intent(out) :: mach1_k
979 :
980 : integer :: k
981 0 : real(dp) :: alfa, beta
982 :
983 : include 'formats'
984 :
985 0 : k = k_shock
986 : if (r < s% R_center .or. r > s% r(1) .or. &
987 : k < 1 .or. k > s% nz .or. &
988 0 : .not. associated(v) .or. &
989 : .not. (s% v_flag .or. s% u_flag)) then
990 0 : mach1_mass = 0
991 0 : mach1_q = 0
992 0 : mach1_radius = 0
993 0 : mach1_velocity = 0
994 0 : mach1_csound = 0
995 0 : mach1_lgT = 0
996 0 : mach1_lgRho = 0
997 0 : mach1_lgP = 0
998 0 : mach1_gamma1 = 0
999 0 : mach1_entropy = 0
1000 0 : mach1_tau = 0
1001 0 : mach1_k = 0
1002 0 : return
1003 : end if
1004 :
1005 0 : mach1_radius = r/Rsun
1006 0 : mach1_k = k
1007 0 : if (k < s% nz) then
1008 0 : alfa = (r - s% r(k))/(s% r(k+1) - s% r(k))
1009 0 : beta = 1d0 - alfa
1010 0 : mach1_mass = (alfa*s% m(k+1) + beta*s% m(k))/Msun
1011 0 : mach1_q = alfa*s% q(k+1) + beta*s% q(k)
1012 0 : mach1_velocity = alfa*v(k+1) + beta*v(k)
1013 : else
1014 0 : mach1_mass = s% m(k)/Msun
1015 0 : mach1_q = s% q(k)
1016 0 : mach1_velocity = v(k)
1017 : end if
1018 0 : mach1_csound = s% csound(k)
1019 0 : mach1_lgT = s% lnT(k)/ln10
1020 0 : mach1_lgRho = s% lnd(k)/ln10
1021 0 : mach1_lgP = s% lnPeos(k)/ln10
1022 0 : mach1_gamma1 = s% gamma1(k)
1023 0 : mach1_entropy = s% entropy(k)
1024 0 : mach1_tau = s% tau(k)
1025 :
1026 : end subroutine get_mach1_location_info
1027 :
1028 :
1029 33 : subroutine get_core_info(s)
1030 : type (star_info), pointer :: s
1031 :
1032 : integer :: k, j, A_max, h1, he4, c12, o16, ne20, si28, species, nz
1033 33 : integer, pointer :: net_iso(:)
1034 : logical :: have_he, have_co, have_one, have_fe
1035 33 : real(dp) :: sumx, min_x
1036 : integer, parameter :: &
1037 : A_max_fe = 47, &
1038 : A_max_si = 28, &
1039 : A_max_one = 20, &
1040 : A_max_co = 16, &
1041 : A_max_he = 4
1042 :
1043 : include 'formats'
1044 :
1045 33 : net_iso => s% net_iso
1046 33 : species = s% species
1047 33 : nz = s% nz
1048 33 : h1 = net_iso(ih1)
1049 33 : he4 = net_iso(ihe4)
1050 33 : c12 = net_iso(ic12)
1051 33 : o16 = net_iso(io16)
1052 33 : ne20 = net_iso(ine20)
1053 33 : si28 = net_iso(isi28)
1054 33 : min_x = s% min_boundary_fraction
1055 :
1056 : call clear_core_info(s% neutron_rich_core_k, &
1057 : s% neutron_rich_core_mass, s% neutron_rich_core_radius, s% neutron_rich_core_lgT, &
1058 : s% neutron_rich_core_lgRho, s% neutron_rich_core_L, s% neutron_rich_core_v, &
1059 33 : s% neutron_rich_core_omega, s% neutron_rich_core_omega_div_omega_crit)
1060 :
1061 40766 : do k=1,nz
1062 40766 : if (s% Ye(k) <= s% neutron_rich_core_boundary_Ye_max) then
1063 : call set_core_info(s, k, s% neutron_rich_core_k, &
1064 : s% neutron_rich_core_mass, s% neutron_rich_core_radius, s% neutron_rich_core_lgT, &
1065 : s% neutron_rich_core_lgRho, s% neutron_rich_core_L, s% neutron_rich_core_v, &
1066 0 : s% neutron_rich_core_omega, s% neutron_rich_core_omega_div_omega_crit)
1067 0 : exit
1068 : end if
1069 : end do
1070 :
1071 : call clear_core_info(s% fe_core_k, &
1072 : s% fe_core_mass, s% fe_core_radius, s% fe_core_lgT, &
1073 : s% fe_core_lgRho, s% fe_core_L, s% fe_core_v, &
1074 33 : s% fe_core_omega, s% fe_core_omega_div_omega_crit)
1075 : call clear_core_info(s% co_core_k, &
1076 : s% co_core_mass, s% co_core_radius, s% co_core_lgT, &
1077 : s% co_core_lgRho, s% co_core_L, s% co_core_v, &
1078 33 : s% co_core_omega, s% co_core_omega_div_omega_crit)
1079 : call clear_core_info(s% one_core_k, &
1080 : s% one_core_mass, s% one_core_radius, s% one_core_lgT, &
1081 : s% one_core_lgRho, s% one_core_L, s% one_core_v, &
1082 33 : s% one_core_omega, s% one_core_omega_div_omega_crit)
1083 : call clear_core_info(s% he_core_k, &
1084 : s% he_core_mass, s% he_core_radius, s% he_core_lgT, &
1085 : s% he_core_lgRho, s% he_core_L, s% he_core_v, &
1086 33 : s% he_core_omega, s% he_core_omega_div_omega_crit)
1087 :
1088 33 : have_he = .false.
1089 33 : have_co = .false.
1090 33 : have_one = .false.
1091 33 : have_fe = .false.
1092 :
1093 40766 : do k=1,nz
1094 :
1095 407330 : j = maxloc(s% xa(1:species,k),dim=1)
1096 40733 : A_max = chem_isos% Z_plus_N(s% chem_id(j))
1097 :
1098 : if (.not. have_fe) then
1099 40733 : if (s% fe_core_boundary_si28_fraction < 0) then
1100 0 : if (A_max >= A_max_fe) then
1101 : call set_core_info(s, k, s% fe_core_k, &
1102 : s% fe_core_mass, s% fe_core_radius, s% fe_core_lgT, &
1103 : s% fe_core_lgRho, s% fe_core_L, s% fe_core_v, &
1104 0 : s% fe_core_omega, s% fe_core_omega_div_omega_crit)
1105 0 : exit
1106 : end if
1107 40733 : else if (si28 /= 0) then
1108 0 : if (s% xa(si28,k) <= s% fe_core_boundary_si28_fraction) then
1109 : sumx = 0
1110 0 : do j=1,species
1111 0 : if (chem_isos% Z_plus_N(s% chem_id(j)) >= A_max_fe) &
1112 0 : sumx = sumx + s% xa(j,k)
1113 : end do
1114 0 : if (sumx >= min_x) then
1115 : call set_core_info(s, k, s% fe_core_k, &
1116 : s% fe_core_mass, s% fe_core_radius, s% fe_core_lgT, &
1117 : s% fe_core_lgRho, s% fe_core_L, s% fe_core_v, &
1118 0 : s% fe_core_omega, s% fe_core_omega_div_omega_crit)
1119 0 : exit
1120 : end if
1121 : end if
1122 : end if
1123 : end if
1124 :
1125 40733 : if (.not. have_one) then
1126 40733 : if (s% one_core_boundary_he4_c12_fraction < 0) then
1127 0 : if (A_max >= A_max_one) then
1128 : call set_core_info(s, k, s% one_core_k, &
1129 : s% one_core_mass, s% one_core_radius, s% one_core_lgT, &
1130 : s% one_core_lgRho, s% one_core_L, s% one_core_v, &
1131 0 : s% one_core_omega, s% one_core_omega_div_omega_crit)
1132 0 : have_one = .true.
1133 : end if
1134 40733 : else if (he4 /= 0 .and. c12 /= 0 .and. o16 /= 0 .and. ne20 /=0) then
1135 40733 : if (s% xa(he4,k)+s% xa(c12,k) <= s% one_core_boundary_he4_c12_fraction .and. &
1136 : s% xa(o16,k)+s% xa(ne20,k) >= min_x) then
1137 : call set_core_info(s, k, s% one_core_k, &
1138 : s% one_core_mass, s% one_core_radius, s% one_core_lgT, &
1139 : s% one_core_lgRho, s% one_core_L, s% one_core_v, &
1140 0 : s% one_core_omega, s% one_core_omega_div_omega_crit)
1141 0 : have_one = .true.
1142 : end if
1143 : end if
1144 : end if
1145 :
1146 40733 : if (.not. have_co) then
1147 40733 : if (s% co_core_boundary_he4_fraction < 0) then
1148 0 : if (A_max >= A_max_co) then
1149 : call set_core_info(s, k, s% co_core_k, &
1150 : s% co_core_mass, s% co_core_radius, s% co_core_lgT, &
1151 : s% co_core_lgRho, s% co_core_L, s% co_core_v, &
1152 0 : s% co_core_omega, s% co_core_omega_div_omega_crit)
1153 0 : have_co = .true.
1154 : end if
1155 40733 : else if (he4 /= 0 .and. c12 /= 0 .and. o16 /= 0) then
1156 40733 : if (s% xa(he4,k) <= s% co_core_boundary_he4_fraction .and. &
1157 : s% xa(c12,k)+s% xa(o16,k) >= min_x) then
1158 : call set_core_info(s, k, s% co_core_k, &
1159 : s% co_core_mass, s% co_core_radius, s% co_core_lgT, &
1160 : s% co_core_lgRho, s% co_core_L, s% co_core_v, &
1161 0 : s% co_core_omega, s% co_core_omega_div_omega_crit)
1162 0 : have_co = .true.
1163 : end if
1164 : end if
1165 : end if
1166 :
1167 40766 : if (.not. have_he) then
1168 40733 : if (s% he_core_boundary_h1_fraction < 0) then
1169 0 : if (A_max >= A_max_he) then
1170 : call set_core_info(s, k, s% he_core_k, &
1171 : s% he_core_mass, s% he_core_radius, s% he_core_lgT, &
1172 : s% he_core_lgRho, s% he_core_L, s% he_core_v, &
1173 0 : s% he_core_omega, s% he_core_omega_div_omega_crit)
1174 0 : have_he = .true.
1175 : end if
1176 40733 : else if (h1 /= 0 .and. he4 /= 0) then
1177 40733 : if (s% xa(h1,k) <= s% he_core_boundary_h1_fraction .and. &
1178 : s% xa(he4,k) >= min_x) then
1179 : call set_core_info(s, k, s% he_core_k, &
1180 : s% he_core_mass, s% he_core_radius, s% he_core_lgT, &
1181 : s% he_core_lgRho, s% he_core_L, s% he_core_v, &
1182 0 : s% he_core_omega, s% he_core_omega_div_omega_crit)
1183 0 : have_he = .true.
1184 : end if
1185 : end if
1186 : end if
1187 :
1188 : end do
1189 :
1190 66 : end subroutine get_core_info
1191 :
1192 :
1193 66 : subroutine clear_core_info( &
1194 : core_k, core_m, core_r, core_lgT, core_lgRho, core_L, core_v, &
1195 : core_omega, core_omega_div_omega_crit)
1196 : integer, intent(out) :: core_k
1197 : real(dp), intent(out) :: &
1198 : core_m, core_r, core_lgT, core_lgRho, core_L, core_v, &
1199 : core_omega, core_omega_div_omega_crit
1200 66 : core_k = 0
1201 66 : core_m = 0
1202 66 : core_r = 0
1203 66 : core_lgT = 0
1204 66 : core_lgRho = 0
1205 66 : core_L = 0
1206 66 : core_v = 0
1207 66 : core_omega = 0
1208 66 : core_omega_div_omega_crit = 0
1209 0 : end subroutine clear_core_info
1210 :
1211 :
1212 0 : subroutine set_core_info(s, k, &
1213 : core_k, core_m, core_r, core_lgT, core_lgRho, core_L, core_v, &
1214 : core_omega, core_omega_div_omega_crit)
1215 : type (star_info), pointer :: s
1216 : integer, intent(in) :: k
1217 : integer, intent(out) :: core_k
1218 : real(dp), intent(out) :: &
1219 : core_m, core_r, core_lgT, core_lgRho, core_L, core_v, &
1220 : core_omega, core_omega_div_omega_crit
1221 :
1222 : integer :: jm1, j00
1223 0 : real(dp) :: dm1, d00, qm1, q00, core_q, &
1224 0 : core_lgP, core_g, core_X, core_Y, &
1225 0 : core_scale_height, core_dlnX_dr, core_dlnY_dr, core_dlnRho_dr
1226 :
1227 : include 'formats'
1228 :
1229 0 : if (k == 1) then
1230 0 : core_q = 1d0
1231 : else
1232 0 : jm1 = maxloc(s% xa(:,k-1), dim=1)
1233 0 : j00 = maxloc(s% xa(:,k), dim=1)
1234 0 : qm1 = s% q(k-1) - 0.5d0*s% dq(k-1) ! center of k-1
1235 0 : q00 = s% q(k) - 0.5d0*s% dq(k) ! center of k
1236 0 : dm1 = s% xa(j00,k-1) - s% xa(jm1,k-1)
1237 0 : d00 = s% xa(j00,k) - s% xa(jm1,k)
1238 0 : if (dm1*d00 > 0d0) then
1239 0 : write(*,2) 'bad args for set_core_info', k, dm1, d00
1240 0 : call mesa_error(__FILE__,__LINE__)
1241 0 : core_q = 0.5d0*(qm1 + q00)
1242 0 : else if (dm1 == 0d0 .and. d00 == 0d0) then
1243 0 : core_q = 0.5d0*(qm1 + q00)
1244 0 : else if (dm1 == 0d0) then
1245 0 : core_q = qm1
1246 0 : else if (d00 == 0d0) then
1247 0 : core_q = q00
1248 : else
1249 0 : core_q = find0(qm1, dm1, q00, d00)
1250 : end if
1251 : end if
1252 :
1253 : call get_info_at_q(s, core_q, &
1254 : core_k, core_m, core_r, core_lgT, core_lgRho, core_L, core_v, &
1255 : core_lgP, core_g, core_X, core_Y, &
1256 : core_scale_height, core_dlnX_dr, core_dlnY_dr, core_dlnRho_dr, &
1257 0 : core_omega, core_omega_div_omega_crit)
1258 :
1259 0 : end subroutine set_core_info
1260 :
1261 :
1262 0 : subroutine get_info_at_q(s, bdy_q, &
1263 : kbdy, bdy_m, bdy_r, bdy_lgT, bdy_lgRho, bdy_L, bdy_v, &
1264 : bdy_lgP, bdy_g, bdy_X, bdy_Y, &
1265 : bdy_scale_height, bdy_dlnX_dr, bdy_dlnY_dr, bdy_dlnRho_dr, &
1266 : bdy_omega, bdy_omega_div_omega_crit)
1267 :
1268 : type (star_info), pointer :: s
1269 : real(dp), intent(in) :: bdy_q
1270 : integer, intent(out) :: kbdy
1271 : real(dp), intent(out) :: &
1272 : bdy_m, bdy_r, bdy_lgT, bdy_lgRho, bdy_L, bdy_v, &
1273 : bdy_lgP, bdy_g, bdy_X, bdy_Y, &
1274 : bdy_scale_height, bdy_dlnX_dr, bdy_dlnY_dr, bdy_dlnRho_dr, &
1275 : bdy_omega, bdy_omega_div_omega_crit
1276 :
1277 0 : real(dp) :: x, x0, x1, x2, alfa, bdy_omega_crit
1278 : integer :: k, klo, khi
1279 :
1280 : include 'formats'
1281 :
1282 0 : bdy_m=0; bdy_r=0; bdy_lgT=0; bdy_lgRho=0; bdy_L=0; bdy_v=0
1283 0 : bdy_lgP=0; bdy_g=0; bdy_X=0; bdy_Y=0
1284 0 : bdy_scale_height=0; bdy_dlnX_dr=0; bdy_dlnY_dr=0; bdy_dlnRho_dr=0
1285 0 : bdy_omega=0; bdy_omega_div_omega_crit=0
1286 0 : kbdy = 0
1287 :
1288 0 : if (bdy_q <= 0) return
1289 0 : k = k_for_q(s,bdy_q)
1290 0 : if (k >= s% nz) then
1291 0 : kbdy = s% nz
1292 0 : return
1293 : end if
1294 0 : if (k <= 1) then
1295 0 : bdy_m = s% star_mass
1296 0 : bdy_r = s% r(1)/Rsun
1297 0 : bdy_lgT = s% lnT(1)/ln10
1298 0 : bdy_lgRho = s% lnd(1)/ln10
1299 0 : bdy_L = s% L(1)/Lsun
1300 0 : if (s% u_flag) then
1301 0 : bdy_v = s% u(1)
1302 0 : else if (s% v_flag) then
1303 0 : bdy_v = s% v(1)
1304 : end if
1305 0 : bdy_lgP = s% lnPeos(1)/ln10
1306 0 : bdy_g = s% grav(1)
1307 0 : bdy_X = s% X(1)
1308 0 : bdy_Y = s% Y(1)
1309 0 : bdy_scale_height = s% scale_height(1)
1310 0 : bdy_omega = s% omega(k)
1311 0 : if (s% rotation_flag) then
1312 0 : bdy_omega_crit = omega_crit(s,1)
1313 0 : bdy_omega_div_omega_crit = min(1d0,bdy_omega/bdy_omega_crit)
1314 : else
1315 0 : bdy_omega_crit = 0
1316 : bdy_omega_div_omega_crit = 0
1317 : end if
1318 0 : kbdy = 1
1319 0 : return
1320 : end if
1321 :
1322 0 : kbdy = k+1
1323 :
1324 0 : bdy_m = (s% M_center + s% xmstar*bdy_q)/Msun
1325 :
1326 0 : x = s% q(k-1) - bdy_q
1327 0 : x0 = s% dq(k-1)/2
1328 0 : x1 = s% dq(k)/2 + s% dq(k-1)
1329 0 : x2 = s% dq(k+1)/2 + s% dq(k) + s% dq(k-1)
1330 :
1331 0 : alfa = max(0d0, min(1d0, (bdy_q - s% q(k+1))/s% dq(k)))
1332 :
1333 0 : bdy_lgT = interp3(s% lnT(k-1), s% lnT(k), s% lnT(k+1))/ln10
1334 0 : bdy_lgRho = interp3(s% lnd(k-1), s% lnd(k), s% lnd(k+1))/ln10
1335 0 : bdy_lgP = interp3(s% lnPeos(k-1), s% lnPeos(k), s% lnPeos(k+1))/ln10
1336 0 : bdy_X = interp3(s% X(k-1), s% X(k), s% X(k+1))
1337 0 : bdy_Y = interp3(s% Y(k-1), s% Y(k), s% Y(k+1))
1338 :
1339 : bdy_r = pow( &
1340 0 : interp2(s% r(k)*s% r(k)*s% r(k), s% r(k+1)*s% r(k+1)*s% r(k+1)),one_third)/Rsun
1341 0 : bdy_L = interp2(s% L(k), s% L(k+1))/Lsun
1342 0 : bdy_g = interp2(s% grav(k), s% grav(k+1))
1343 0 : bdy_scale_height = interp2(s% scale_height(k), s% scale_height(k+1))
1344 :
1345 0 : klo = k-1
1346 0 : khi = k+1
1347 : bdy_dlnX_dr = log(max(1d-99,max(1d-99,s% X(klo))/max(1d-99,s% X(khi)))) &
1348 0 : / (s% rmid(klo) - s% rmid(khi))
1349 : bdy_dlnY_dr = log(max(1d-99,max(1d-99,s% Y(klo))/max(1d-99,s% Y(khi)))) &
1350 0 : / (s% rmid(klo) - s% rmid(khi))
1351 0 : bdy_dlnRho_dr = (s% lnd(klo) - s% lnd(khi))/(s% rmid(klo) - s% rmid(khi))
1352 :
1353 0 : if (s% u_flag) then
1354 0 : bdy_v = interp2(s% u(k), s% u(k+1))
1355 0 : else if (s% v_flag) then
1356 0 : bdy_v = interp2(s% v(k), s% v(k+1))
1357 : end if
1358 0 : if (s% rotation_flag) then
1359 0 : bdy_omega = interp2(s% omega(k), s% omega(k+1))
1360 0 : bdy_omega_crit = interp2(omega_crit(s,k), omega_crit(s,k+1))
1361 0 : bdy_omega_div_omega_crit = min(1d0,bdy_omega/bdy_omega_crit)
1362 : else
1363 : bdy_omega = 0
1364 0 : bdy_omega_crit = 0
1365 : bdy_omega_div_omega_crit = 0
1366 : end if
1367 :
1368 : contains
1369 :
1370 0 : real(dp) function interp3(f0, f1, f2)
1371 : real(dp), intent(in) :: f0, f1, f2
1372 0 : real(dp) :: fmin, fmax
1373 0 : fmin = min(f0,f1,f2)
1374 0 : fmax = max(f0,f1,f2)
1375 : interp3 = (f1*(x-x0)*(x-x2)*(x0-x2)-&
1376 : (x-x1)*(f0*(x-x2)*(x1-x2) + (x-x0)*(x0-x1)*x2))/ &
1377 0 : ((x0-x1)*(x0-x2)*(-x1+x2))
1378 0 : interp3 = min(fmax, max(fmin, interp3))
1379 0 : end function interp3
1380 :
1381 0 : real(dp) function interp2(f0, f1)
1382 : real(dp), intent(in) :: f0, f1
1383 0 : interp2 = alfa*f0 + (1-alfa)*f1
1384 : end function interp2
1385 :
1386 : end subroutine get_info_at_q
1387 :
1388 :
1389 33 : subroutine get_mixing_regions(s, ierr)
1390 : type (star_info), pointer :: s
1391 : integer, intent(out) :: ierr
1392 : integer :: prev_type, cur_type, cur_top, n, k
1393 : include 'formats'
1394 33 : ierr = 0
1395 33 : cur_type = s% mixing_type(1)
1396 33 : cur_top = 1
1397 33 : n = 0
1398 40733 : do k = 2, s% nz
1399 40700 : prev_type = cur_type
1400 40700 : cur_type = s% mixing_type(k)
1401 40700 : if (cur_type == prev_type .and. k < s% nz) cycle
1402 : ! change of type from k-1 to k
1403 138 : if (prev_type /= no_mixing) then
1404 69 : n = n + 1
1405 69 : s% mixing_region_type(n) = prev_type
1406 69 : s% mixing_region_top(n) = cur_top
1407 69 : if (k == s% nz) then
1408 33 : s% mixing_region_bottom(n) = k
1409 : else
1410 36 : s% mixing_region_bottom(n) = k-1
1411 : end if
1412 69 : if (n == max_num_mixing_regions) exit
1413 : end if
1414 40595 : cur_top = k
1415 : end do
1416 33 : s% num_mixing_regions = n
1417 33 : end subroutine get_mixing_regions
1418 :
1419 : end module report
|
