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 star_utils
21 :
22 : use star_private_def
23 : use const_def, only: dp, i8, pi, pi4, ln10, clight, crad, msun, rsun, lsun, one_third, four_thirds_pi
24 : use num_lib
25 : use utils_lib
26 : use auto_diff_support
27 :
28 : implicit none
29 :
30 : private
31 : public :: eval_min_cell_collapse_time
32 : public :: min_dr_div_cs
33 : public :: get_XYZ
34 : public :: lookup_nameofvar
35 : public :: start_time
36 : public :: update_time
37 : public :: foreach_cell
38 : public :: write_eos_call_info
39 : public :: eval_csound
40 : public :: eval_current_abundance
41 : public :: eval_current_z
42 : public :: eval_ledd
43 : public :: normalize_dqs
44 : public :: set_qs
45 : public :: set_m_and_dm
46 : public :: set_dm_bar
47 : public :: set_rmid
48 : public :: get_r_from_xh
49 : public :: get_r_and_lnr_from_xh
50 : public :: get_t_and_lnt_from_xh
51 : public :: get_lnt_from_xh
52 : public :: get_lnr_from_xh
53 : public :: get_rho_and_lnd_from_xh
54 : public :: store_t_in_xh
55 : public :: store_r_in_xh
56 : public :: store_rho_in_xh
57 : public :: store_lnd_in_xh
58 : public :: store_lnt_in_xh
59 : public :: cell_specific_ke
60 : public :: cell_specific_pe
61 : public :: get_phot_info
62 : public :: init_random
63 : public :: rand
64 : public :: interp_val_to_pt
65 : public :: get_log_concentration
66 : public :: get_lconv
67 : public :: get_effective_mixing_speed
68 : public :: get_ladv
69 : public :: get_lrad_div_ledd
70 : public :: get_lrad
71 : public :: get_ledd
72 : public :: tau_eff
73 : public :: get_rho_face_val
74 : public :: omega_crit
75 : public :: eval_cell_section_total_energy
76 : public :: conv_time_scale
77 : public :: qhse_time_scale
78 : public :: eps_nuc_time_scale
79 : public :: cooling_time_scale
80 : public :: get_face_values
81 : public :: threshold_smoothing
82 : public :: save_eqn_residual_info
83 : public :: calc_ptrb_ad_tw
84 : public :: set_energy_eqn_scal
85 : public :: get_scale_height_face_val
86 : public :: weighted_smoothing
87 : public :: get_kap_face
88 : public :: get_rho_face
89 : public :: get_rho_start_face
90 : public :: get_e_face
91 : public :: get_ChiRho_face
92 : public :: get_ChiT_face
93 : public :: get_T_face
94 : public :: get_Peos_face
95 : public :: get_Peos_face_val
96 : public :: get_cp_face
97 : public :: get_grada_face
98 : public :: get_gradr_face
99 : public :: get_scale_height_face
100 : public :: get_tau
101 : public :: after_c_burn
102 : public :: get_shock_info
103 : public :: reset_starting_vectors
104 : public :: eval_irradiation_heat
105 : public :: set_phot_info
106 : public :: set_m_grav_and_grav
107 : public :: set_scale_height
108 : public :: set_abs_du_div_cs
109 : public :: set_conv_time_scales
110 : public :: set_rv_info
111 : public :: smooth
112 : public :: safe_div_val
113 : public :: center_avg_x
114 : public :: surface_avg_x
115 : public :: get_remnant_mass
116 : public :: get_ejecta_mass
117 : public :: get_phi_joss
118 : public :: center_value
119 : public :: eval_kh_timescale
120 : public :: k_for_q
121 : public :: total_angular_momentum
122 : public :: reset_epsnuc_vectors
123 : public :: yrs_for_init_timestep
124 : public :: set_phase_of_evolution
125 : public :: get_string_for_model_number
126 : public :: e00
127 : public :: em1
128 : public :: ep1
129 : public :: store_partials
130 : public :: save_eqn_dxa_partials
131 : public :: unpack_residual_partials
132 : public :: get_area_info_opt_time_center
133 : public :: calc_ptot_ad_tw
134 : public :: get_face_weights
135 : public :: get_dke_dt_dpe_dt
136 : public :: get_pvsc_ad
137 : public :: show_matrix
138 : public :: no_extra_profile_columns
139 : public :: no_data_for_extra_profile_columns
140 : public :: total_times
141 : public :: current_min_xa_hard_limit
142 : public :: current_sum_xa_hard_limit
143 : public :: lookup_nameofequ
144 : public :: eval_total_energy_integrals
145 : public :: set_luminosity_by_category
146 : public :: eval_deltam_total_energy_integrals
147 : public :: report_xa_bad_nums
148 : public :: eval_current_y
149 : public :: get_name_for_restart_file
150 : public :: eval_integrated_total_energy_profile
151 : public :: use_xh_to_set_rho_to_dm_div_dv
152 : public :: cell_specific_total_energy
153 : public :: store_lnr_in_xh
154 : public :: interp_q
155 : public :: set_zero_alpha_rti
156 : public :: after_he_burn
157 : public :: interp_xa_to_pt
158 : public :: set_xqs
159 : public :: get_tau_at_r
160 : public :: smooth_abundances
161 : public :: do_boxcar_mixing
162 : public :: eval_total_energy_profile
163 : public :: get_delta_pg_traditional
164 : public :: get_delta_pg_bildsten2012
165 : public :: write_to_extra_terminal_output_file
166 :
167 : logical, parameter :: mdb = .false.
168 :
169 : contains
170 :
171 66 : subroutine foreach_cell(s,nzlo,nzhi,use_omp,do1,ierr)
172 : type (star_info), pointer :: s
173 : integer, intent(in) :: nzlo, nzhi
174 : logical, intent(in) :: use_omp
175 : interface
176 : subroutine do1(s,k,ierr)
177 : use star_private_def
178 : implicit none
179 : type (star_info), pointer :: s
180 : integer, intent(in) :: k
181 : integer, intent(out) :: ierr
182 : end subroutine do1
183 : end interface
184 : integer, intent(out) :: ierr
185 :
186 : integer :: k, op_err
187 : logical :: okay
188 66 : ierr = 0
189 :
190 66 : if (nzlo == nzhi) then
191 0 : call do1(s,nzlo,ierr)
192 0 : return
193 : end if
194 :
195 66 : if (use_omp) then
196 66 : okay = .true.
197 66 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
198 : do k = nzlo, nzhi
199 : if (.not. okay) cycle
200 : op_err = 0
201 : call do1(s,k,op_err)
202 : if (op_err /= 0) okay = .false. ! cannot just exit from a parallel loop
203 : end do
204 : !$OMP END PARALLEL DO
205 66 : if (.not. okay) ierr = -1
206 : else
207 0 : do k = nzlo, nzhi
208 0 : call do1(s,k,ierr)
209 0 : if (ierr /= 0) exit
210 : end do
211 : end if
212 :
213 : end subroutine foreach_cell
214 :
215 :
216 0 : subroutine get_average_Y_and_Z(s, nzlo, nzhi, y_avg, z_avg, ierr)
217 : use chem_def
218 : type (star_info), pointer :: s
219 : integer, intent(in) :: nzlo, nzhi
220 : real(dp), intent(out) :: y_avg, z_avg
221 : integer, intent(out) :: ierr
222 : integer :: k, nz, h1, h2, he3, he4
223 : real(dp) :: total_mass_h, total_mass_he, total_mass_z, cell_mass, total_mass
224 0 : ierr = 0
225 0 : nz = s% nz
226 0 : h1 = s% net_iso(ih1)
227 0 : h2 = s% net_iso(ih2)
228 0 : he3 = s% net_iso(ihe3)
229 0 : he4 = s% net_iso(ihe4)
230 0 : total_mass=0; total_mass_h=0; total_mass_he=0; total_mass_z=0
231 0 : do k=nzlo, nzhi
232 0 : cell_mass = s% dm(k)
233 0 : total_mass = total_mass + cell_mass
234 0 : total_mass_h = total_mass_h + cell_mass*s% xa(h1, k)
235 0 : if (h2 /= 0) total_mass_h = total_mass_h + cell_mass*s% xa(h2, k)
236 0 : total_mass_he = total_mass_he + cell_mass*s% xa(he4, k)
237 0 : if (he3 /= 0) total_mass_he = total_mass_he + cell_mass*s% xa(he3, k)
238 : end do
239 0 : total_mass_z = total_mass - (total_mass_h + total_mass_he)
240 0 : z_avg = total_mass_z / total_mass
241 0 : y_avg = total_mass_he / total_mass
242 0 : end subroutine get_average_Y_and_Z
243 :
244 :
245 0 : real(dp) function eval_current_y(s, nzlo, nzhi, ierr)
246 : type (star_info), pointer :: s
247 : integer, intent(in) :: nzlo, nzhi
248 : integer, intent(out) :: ierr
249 : real(dp) :: y_avg, z_avg
250 0 : call get_average_Y_and_Z(s, nzlo, nzhi, y_avg, z_avg, ierr)
251 0 : eval_current_y = y_avg
252 0 : end function eval_current_y
253 :
254 :
255 0 : real(dp) function eval_current_z(s, nzlo, nzhi, ierr)
256 : type (star_info), pointer :: s
257 : integer, intent(in) :: nzlo, nzhi
258 : integer, intent(out) :: ierr
259 : real(dp) :: y_avg, z_avg
260 0 : call get_average_Y_and_Z(s, nzlo, nzhi, y_avg, z_avg, ierr)
261 0 : eval_current_z = z_avg
262 0 : end function eval_current_z
263 :
264 :
265 11 : real(dp) function eval_current_abundance(s, j, nzlo, nzhi, ierr)
266 : type (star_info), pointer :: s
267 : integer, intent(in) :: j, nzlo, nzhi
268 : integer, intent(out) :: ierr
269 : integer :: k, nz
270 : real(dp) :: cell_mass, jmass, total_mass
271 :
272 11 : ierr = 0
273 :
274 11 : if (j == 0) then
275 11 : eval_current_abundance = 0
276 : return
277 : end if
278 :
279 11 : nz = s% nz
280 11 : total_mass=0; jmass=0
281 22 : do k=nzlo, nzhi
282 11 : cell_mass = s% dm(k)
283 11 : total_mass = total_mass + cell_mass
284 22 : jmass = jmass + cell_mass*s% xa(j, k)
285 : end do
286 11 : eval_current_abundance = jmass / total_mass
287 :
288 11 : end function eval_current_abundance
289 :
290 :
291 0 : subroutine smooth_abundances(s, cnt, nzlo, nzhi, ierr)
292 : type (star_info), pointer :: s
293 : integer, intent(in) :: cnt ! make this many passes
294 : integer, intent(in) :: nzlo, nzhi ! only smooth zones nzlo to nzhi inclusive
295 : integer, intent(out) :: ierr
296 : integer :: k, j, nz
297 0 : ierr = 0
298 0 : nz = s% nz
299 0 : do j = 1, cnt
300 0 : do k = max(nzlo,2), min(nzhi, nz)
301 0 : s% xa(:,k) = (s% xa(:,k-1) + s% xa(:,k) + s% xa(:,k+1))/3
302 : end do
303 0 : if (nzhi == nz) s% xa(:,nz) = (s% xa(:,nz-1) + s% xa(:,nz) + s% xa(:,nz))/3
304 0 : if (nzlo == 1) s% xa(:,1) = (s% xa(:,2) + s% xa(:,1) + s% xa(:,1))/3
305 : end do
306 0 : s% need_to_setvars = .true.
307 0 : end subroutine smooth_abundances
308 :
309 :
310 0 : integer function k_for_q(s, q)
311 : ! return k s.t. q(k) >= q > q(k)-dq(k)
312 : type (star_info), pointer :: s
313 : real(dp), intent(in) :: q
314 : integer :: k, nz
315 0 : nz = s% nz
316 0 : if (q >= 1) then
317 : k_for_q = 1; return
318 0 : else if (q <= s% q(nz)) then
319 : k_for_q = nz; return
320 : end if
321 0 : do k = 1, nz-1
322 0 : if (q > s% q(k+1)) then
323 : k_for_q = k; return
324 : end if
325 : end do
326 : k_for_q = nz
327 : end function k_for_q
328 :
329 :
330 1 : subroutine get_name_for_restart_file(n, num_digits, num)
331 : integer, intent(in) :: n, num_digits
332 : character (len=*), intent(out) :: num
333 1 : call get_string_for_model_number('x', n, num_digits, num)
334 1 : end subroutine get_name_for_restart_file
335 :
336 :
337 1 : subroutine get_string_for_model_number(prefix, n, num_digits, num)
338 : character (len=*), intent(in) :: prefix
339 : integer, intent(in) :: n, num_digits
340 : character (len=*), intent(out) :: num
341 : integer :: val
342 : character (len=32) :: fstring
343 : include 'formats'
344 1 : val = mod(n, 10**num_digits) ! wrap around
345 1 : if (val == 0) then
346 0 : write(num,*) n
347 0 : num = adjustl(num)
348 0 : return
349 : end if
350 1 : write(fstring,'( "(a,i",i2.2,".",i2.2,")" )') num_digits, num_digits
351 1 : write(num,fstring) trim(prefix), val
352 1 : end subroutine get_string_for_model_number
353 :
354 :
355 0 : subroutine report_xa_bad_nums(s,ierr)
356 : type (star_info), pointer :: s
357 : integer, intent(out) :: ierr
358 : integer :: k, j
359 0 : ierr = 0
360 0 : do k=1,s% nz
361 0 : do j=1,s% species
362 0 : if (is_bad(s% xa(j,k))) then
363 0 : ierr = -1
364 0 : write(*,*) j, k, s% xa(j,k)
365 : end if
366 : end do
367 : end do
368 0 : end subroutine report_xa_bad_nums
369 :
370 :
371 79998 : real(dp) function eval_csound(s,k,ierr) result(cs)
372 : type (star_info), pointer :: s
373 : integer, intent(in) :: k
374 : integer, intent(out) :: ierr
375 : real(dp) :: cs2
376 : include 'formats'
377 79998 : ierr = 0
378 79998 : cs2 = s% gamma1(k)*s% Peos(k)/s% rho(k)
379 79998 : if (cs2 < 0d0) then
380 0 : cs = 0d0
381 0 : ierr = -1
382 0 : return
383 : end if
384 79998 : cs = sqrt(cs2)
385 79998 : end function eval_csound
386 :
387 :
388 55 : subroutine set_m_grav_and_grav(s) ! using mass_corrections
389 : type (star_info), pointer :: s
390 : real(dp) :: r, lnR
391 : integer :: k, nz
392 : include 'formats'
393 55 : nz = s% nz
394 55 : if (.not. s% use_mass_corrections) then
395 66962 : do k=1,nz
396 66962 : s% m_grav(k) = s% m(k)
397 : end do
398 : else
399 : s% m_grav(nz) = &
400 0 : s% M_center + s% dm(nz)*s% mass_correction(nz)
401 0 : do k=nz-1,1,-1
402 : s% m_grav(k) = &
403 0 : s% m_grav(k+1) + s% dm(k)*s% mass_correction(k)
404 : end do
405 : end if
406 66962 : do k=1,nz
407 : ! We need to call set_m_grav_and_grav during model loading before we have set all vars
408 66907 : call get_r_and_lnR_from_xh(s, k, r, lnR)
409 66962 : s% grav(k) = s% cgrav(k)*s% m_grav(k)/(r*r)
410 : end do
411 55 : end subroutine set_m_grav_and_grav
412 :
413 :
414 66907 : subroutine get_r_and_lnR_from_xh(s, k, r, lnR, xh_in)
415 : type (star_info), pointer :: s
416 : integer, intent(in) :: k
417 : real(dp), intent(out) :: r, lnR
418 : real(dp), intent(in), pointer, optional :: xh_in(:,:)
419 66907 : real(dp), pointer :: xh(:,:)
420 66907 : if (present(xh_in)) then
421 0 : xh => xh_in
422 : else
423 66907 : xh => s% xh
424 : end if
425 66907 : lnR = xh(s% i_lnR,k)
426 66907 : r = exp(lnR)
427 66907 : end subroutine get_r_and_lnR_from_xh
428 :
429 :
430 14559 : real(dp) function get_r_from_xh(s, k, xh_in) result(r)
431 : type (star_info), pointer :: s
432 : integer, intent(in) :: k
433 : real(dp), intent(in), pointer, optional :: xh_in(:,:)
434 14559 : real(dp), pointer :: xh(:,:)
435 14559 : if (present(xh_in)) then
436 0 : xh => xh_in
437 : else
438 14559 : xh => s% xh
439 : end if
440 14559 : r = exp(xh(s% i_lnR,k))
441 14559 : end function get_r_from_xh
442 :
443 :
444 0 : real(dp) function get_lnR_from_xh(s, k, xh_in) result(lnR)
445 : type (star_info), pointer :: s
446 : integer, intent(in) :: k
447 : real(dp), intent(in), pointer, optional :: xh_in(:,:)
448 0 : real(dp), pointer :: xh(:,:)
449 0 : if (present(xh_in)) then
450 0 : xh => xh_in
451 : else
452 0 : xh => s% xh
453 : end if
454 0 : lnR = xh(s% i_lnR,k)
455 0 : end function get_lnR_from_xh
456 :
457 1434 : subroutine store_r_in_xh(s, k, r, xh_in)
458 : type (star_info), pointer :: s
459 : integer, intent(in) :: k
460 : real(dp), intent(in) :: r
461 : real(dp), intent(in), pointer, optional :: xh_in(:,:)
462 1434 : real(dp), pointer :: xh(:,:)
463 1434 : if (present(xh_in)) then
464 1434 : xh => xh_in
465 : else
466 0 : xh => s% xh
467 : end if
468 1434 : xh(s% i_lnR,k) = log(r)
469 1434 : end subroutine store_r_in_xh
470 :
471 :
472 0 : subroutine store_lnR_in_xh(s, k, lnR, xh_in)
473 : type (star_info), pointer :: s
474 : integer, intent(in) :: k
475 : real(dp), intent(in) :: lnR
476 : real(dp), intent(in), pointer, optional :: xh_in(:,:)
477 0 : real(dp), pointer :: xh(:,:)
478 0 : if (present(xh_in)) then
479 0 : xh => xh_in
480 : else
481 0 : xh => s% xh
482 : end if
483 0 : xh(s% i_lnR,k) = lnR
484 0 : end subroutine store_lnR_in_xh
485 :
486 :
487 0 : subroutine get_T_and_lnT_from_xh(s, k, T, lnT, xh_in)
488 : type (star_info), pointer :: s
489 : integer, intent(in) :: k
490 : real(dp), intent(out) :: T, lnT
491 : real(dp), intent(in), pointer, optional :: xh_in(:,:)
492 0 : real(dp), pointer :: xh(:,:)
493 0 : if (present(xh_in)) then
494 0 : xh => xh_in
495 : else
496 0 : xh => s% xh
497 : end if
498 0 : lnT = xh(s% i_lnT,k)
499 0 : T = exp(lnT)
500 0 : end subroutine get_T_and_lnT_from_xh
501 :
502 :
503 : real(dp) function get_T_from_xh(s, k, xh_in) result(T)
504 : type (star_info), pointer :: s
505 : integer, intent(in) :: k
506 : real(dp), intent(in), pointer, optional :: xh_in(:,:)
507 : real(dp), pointer :: xh(:,:)
508 : if (present(xh_in)) then
509 : xh => xh_in
510 : else
511 : xh => s% xh
512 : end if
513 : T = exp(xh(s% i_lnT,k))
514 : end function get_T_from_xh
515 :
516 :
517 0 : real(dp) function get_lnT_from_xh(s, k, xh_in) result(lnT)
518 : type (star_info), pointer :: s
519 : integer, intent(in) :: k
520 : real(dp), intent(in), pointer, optional :: xh_in(:,:)
521 0 : real(dp), pointer :: xh(:,:)
522 0 : if (present(xh_in)) then
523 0 : xh => xh_in
524 : else
525 0 : xh => s% xh
526 : end if
527 0 : lnT = xh(s% i_lnT,k)
528 0 : end function get_lnT_from_xh
529 :
530 :
531 0 : subroutine store_T_in_xh(s, k, T, xh_in)
532 : type (star_info), pointer :: s
533 : integer, intent(in) :: k
534 : real(dp), intent(in) :: T
535 : real(dp), intent(in), pointer, optional :: xh_in(:,:)
536 0 : real(dp), pointer :: xh(:,:)
537 0 : if (present(xh_in)) then
538 0 : xh => xh_in
539 : else
540 0 : xh => s% xh
541 : end if
542 0 : xh(s% i_lnT,k) = log(T)
543 0 : end subroutine store_T_in_xh
544 :
545 :
546 1472 : subroutine store_lnT_in_xh(s, k, lnT, xh_in)
547 : type (star_info), pointer :: s
548 : integer, intent(in) :: k
549 : real(dp), intent(in) :: lnT
550 : real(dp), intent(in), pointer, optional :: xh_in(:,:)
551 1472 : real(dp), pointer :: xh(:,:)
552 1472 : if (present(xh_in)) then
553 1472 : xh => xh_in
554 : else
555 0 : xh => s% xh
556 : end if
557 1472 : xh(s% i_lnT,k) = lnT
558 1472 : end subroutine store_lnT_in_xh
559 :
560 :
561 0 : subroutine get_rho_and_lnd_from_xh(s, k, rho, lnd, xh_in)
562 : type (star_info), pointer :: s
563 : integer, intent(in) :: k
564 : real(dp), intent(out) :: rho, lnd
565 : real(dp), intent(in), pointer, optional :: xh_in(:,:)
566 0 : real(dp), pointer :: xh(:,:)
567 0 : if (present(xh_in)) then
568 0 : xh => xh_in
569 : else
570 0 : xh => s% xh
571 : end if
572 0 : lnd = xh(s% i_lnd,k)
573 0 : rho = exp(lnd)
574 0 : end subroutine get_rho_and_lnd_from_xh
575 :
576 :
577 17509 : subroutine store_rho_in_xh(s, k, rho, xh_in)
578 : type (star_info), pointer :: s
579 : integer, intent(in) :: k
580 : real(dp), intent(in) :: rho
581 : real(dp), intent(in), pointer, optional :: xh_in(:,:)
582 17509 : real(dp), pointer :: xh(:,:)
583 17509 : if (present(xh_in)) then
584 2950 : xh => xh_in
585 : else
586 14559 : xh => s% xh
587 : end if
588 17509 : xh(s% i_lnd,k) = log(rho)
589 17509 : end subroutine store_rho_in_xh
590 :
591 :
592 493 : subroutine store_lnd_in_xh(s, k, lnd, xh_in)
593 : type (star_info), pointer :: s
594 : integer, intent(in) :: k
595 : real(dp), intent(in) :: lnd
596 : real(dp), intent(in), pointer, optional :: xh_in(:,:)
597 493 : real(dp), pointer :: xh(:,:)
598 493 : if (present(xh_in)) then
599 493 : xh => xh_in
600 : else
601 0 : xh => s% xh
602 : end if
603 493 : xh(s% i_lnd,k) = lnd
604 493 : end subroutine store_lnd_in_xh
605 :
606 :
607 11 : subroutine use_xh_to_set_rho_to_dm_div_dV(s, ierr)
608 : type (star_info), pointer :: s
609 : integer, intent(out) :: ierr
610 : integer :: k, nz, i_lnR, i_lnd
611 : real(dp) :: rL, rR, dm, dV, rho
612 : include 'formats'
613 11 : ierr = 0
614 11 : i_lnR = s% i_lnR
615 11 : i_lnd = s% i_lnd
616 11 : if (i_lnd == 0 .or. i_lnR == 0) return
617 11 : nz = s% nz
618 11 : rR = s% R_center
619 14570 : do k = nz, 1, -1
620 14559 : rL = rR
621 14559 : rR = get_r_from_xh(s,k)
622 14559 : dm = s% dm(k)
623 14559 : dV = four_thirds_pi*(rR*rR*rR - rL*rL*rL)
624 14559 : rho = dm/dV
625 14559 : if (rho <= 0d0) then
626 0 : write(*,3) 'set_rho_to_dm_div_dV: rho <= 0', &
627 0 : k, nz, rho, dm, dV, rR, rL
628 0 : ierr = -1
629 0 : return
630 : end if
631 14570 : call store_rho_in_xh(s, k, rho)
632 : end do
633 : end subroutine use_xh_to_set_rho_to_dm_div_dV
634 :
635 :
636 66 : subroutine set_m_and_dm(s)
637 : type (star_info), pointer :: s
638 : integer :: k
639 : include 'formats'
640 81532 : do k = 1, s% nz
641 81466 : s% m(k) = s% M_center + s% q(k)*s% xmstar
642 81466 : s% dm(k) = s% dq(k)*s% xmstar
643 81532 : if (s% dm(k) <= 0d0 .or. is_bad(s% m(k) + s% dm(k))) then
644 0 : write(*,2) 'dm m dq q M_center', k, &
645 0 : s% dm(k), s% m(k), s% dq(k), s% q(k), s% M_center
646 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'set_m_and_dm')
647 : end if
648 : end do
649 66 : end subroutine set_m_and_dm
650 :
651 :
652 66 : subroutine set_dm_bar(s, nz, dm, dm_bar)
653 : type (star_info), pointer :: s
654 : integer, intent(in) :: nz
655 : real(dp), intent(in) :: dm(:) ! (nz)
656 : real(dp), intent(inout) :: dm_bar(:) ! (nz)
657 : integer :: k
658 81400 : do k=2,nz-1
659 81400 : dm_bar(k) = 0.5d0*(dm(k-1) + dm(k))
660 : end do
661 66 : dm_bar(1) = 0.5d0*dm(1)
662 66 : if (s% rsp_flag .or. s% RSP2_flag) then ! rsp and RSP2 use this definition
663 0 : dm_bar(nz) = 0.5d0*(dm(nz-1) + dm(nz))
664 : else
665 66 : dm_bar(nz) = 0.5d0*dm(nz-1) + dm(nz)
666 : end if
667 66 : end subroutine set_dm_bar
668 :
669 :
670 0 : subroutine normalize_dqs(s, nz, dq, ierr)
671 : ! rescale dq's so that add to 1.000
672 : ! work in from boundaries to meet at largest dq
673 : type (star_info), pointer :: s
674 : integer, intent(in) :: nz
675 : real(dp), intent(inout) :: dq(:) ! (nz)
676 : integer, intent(out) :: ierr
677 : integer :: k, midq
678 : real(dp) :: dqsum1, dqsum2, dq_min
679 : include 'formats'
680 0 : k = minloc(dq(1:nz),dim=1)
681 0 : dq_min = dq(k)
682 0 : if (dq_min <= 0d0) then
683 0 : write(*,2) 'bad dq', k, dq(k)
684 0 : ierr = -1
685 0 : stop
686 0 : return
687 : end if
688 0 : midq = maxloc(dq(1:nz),dim=1)
689 : ! surface inward
690 0 : dqsum1 = 0
691 0 : do k=1, midq
692 0 : dqsum1 = dqsum1 + dq(k)
693 0 : if (dq(k) <= 0) then
694 0 : ierr = -1
695 0 : return
696 : end if
697 : end do
698 : ! center outward
699 0 : dqsum2 = 0
700 0 : do k=nz, midq+1, -1
701 0 : dqsum2 = dqsum2 + dq(k)
702 0 : if (dq(k) <= 0) then
703 0 : ierr = -1
704 0 : return
705 : end if
706 : end do
707 0 : do k=1,nz
708 0 : dq(k) = dq(k)/(dqsum1 + dqsum2)
709 : end do
710 : end subroutine normalize_dqs
711 :
712 :
713 68 : subroutine set_qs(s, nz, q, dq, ierr) ! set q's using normalized dq's
714 : type (star_info), pointer :: s
715 : integer, intent(in) :: nz
716 : real(dp), intent(inout) :: dq(:) ! (nz)
717 : real(dp), intent(inout) :: q(:) ! (nz)
718 : integer, intent(out) :: ierr
719 : integer :: k, midq
720 : real(dp) :: dqsum1, dqsum2
721 : logical :: okay
722 : include 'formats'
723 68 : ierr = 0
724 : ! normalize_dqs destroys bit-for-bit read as inverse of write for models.
725 : ! ok for create pre ms etc., but not for read model
726 68 : if (s% do_normalize_dqs_as_part_of_set_qs) then
727 0 : call normalize_dqs(s, nz, dq, ierr)
728 0 : if (ierr /= 0) return
729 : end if
730 68 : q(1) = 1d0
731 68 : okay = .true.
732 86396 : do k=2,nz
733 86328 : q(k) = q(k-1) - dq(k-1)
734 86396 : if (q(k) < 0d0 .or. q(k) > 1d0) then
735 : okay = .false.
736 : exit
737 : end if
738 : end do
739 68 : if (okay) return
740 0 : midq = maxloc(dq(1:nz),dim=1)
741 : ! surface inward
742 0 : dqsum1 = 0
743 0 : do k=1, midq
744 0 : q(k) = 1d0 - dqsum1
745 0 : dqsum1 = dqsum1 + dq(k)
746 : end do
747 : ! center outward
748 : dqsum2 = 0
749 0 : do k=nz, midq+1, -1
750 0 : dqsum2 = dqsum2 + dq(k)
751 0 : q(k) = dqsum2
752 : end do
753 : end subroutine set_qs
754 :
755 :
756 20 : subroutine set_xqs(nz, xq, dq, ierr) ! set xq's using dq's
757 : integer, intent(in) :: nz
758 : real(dp), intent(inout) :: dq(:) ! (nz)
759 : real(dp), intent(inout) :: xq(:) ! (nz)
760 : integer, intent(out) :: ierr
761 : integer :: k
762 : include 'formats'
763 20 : ierr = 0
764 20 : xq(1) = 0
765 22698 : do k=2,nz-1
766 22698 : xq(k) = xq(k-1) + dq(k-1)
767 : end do
768 20 : xq(nz) = 1 - dq(nz)
769 20 : if (xq(nz) < xq(nz-1)) then
770 0 : xq(nz) = xq(nz-1) + dq(nz-1)
771 0 : dq(nz) = 1 - xq(nz)
772 0 : if (dq(nz) <= 0) then
773 0 : ierr = -1
774 0 : return
775 : end if
776 : end if
777 : end subroutine set_xqs
778 :
779 :
780 0 : real(dp) function interp_val_to_pt(v,k,sz,dq,str)
781 : use interp_1d_lib, only: interp_4_to_1
782 : integer, intent(in) :: k, sz
783 : real(dp), intent(in) :: v(:), dq(:)
784 : character (len=*), intent(in) :: str
785 : integer :: ierr
786 : include 'formats'
787 0 : if (k == 1) then
788 0 : interp_val_to_pt = v(k)
789 0 : return
790 : end if
791 0 : if (k > 2 .and. k < sz) then
792 : ierr = 0
793 : call interp_4_to_1( &
794 : 0.5d0*(dq(k-2)+dq(k-1)), &
795 : 0.5d0*(dq(k-1)+dq(k)), &
796 : 0.5d0*(dq(k)+dq(k+1)), &
797 : 0.5d0*dq(k-2)+dq(k-1), &
798 : v(k-2), v(k-1), v(k), v(k+1), &
799 0 : interp_val_to_pt, str, ierr)
800 0 : if (ierr == 0) return
801 0 : write(*,1) '0.5d0*(dq(k-2)+dq(k-1))', 0.5d0*(dq(k-2)+dq(k-1))
802 0 : write(*,1) '0.5d0*(dq(k-1)+dq(k))', 0.5d0*(dq(k-1)+dq(k))
803 0 : write(*,1) '0.5d0*(dq(k)+dq(k+1))', 0.5d0*(dq(k)+dq(k+1))
804 0 : write(*,2) 'dq(k-2)', k-2, dq(k-2)
805 0 : write(*,2) 'dq(k-1)', k-1, dq(k-1)
806 0 : write(*,2) 'dq(k)', k, dq(k)
807 0 : write(*,2) 'dq(k+1)', k+1, dq(k+1)
808 :
809 0 : call mesa_error(__FILE__,__LINE__,'interp_val_to_pt')
810 : end if
811 0 : interp_val_to_pt = (v(k)*dq(k-1) + v(k-1)*dq(k))/(dq(k-1) + dq(k))
812 0 : end function interp_val_to_pt
813 :
814 :
815 0 : real(dp) function interp_xa_to_pt(xa,j,k,sz,dq,str)
816 : use interp_1d_lib, only: interp_4_to_1
817 : real(dp), intent(in) :: xa(:,:), dq(:)
818 : character (len=*), intent(in) :: str
819 : integer, intent(in) :: j, k, sz
820 : integer :: ierr
821 : include 'formats'
822 0 : if (j == 0) then
823 : interp_xa_to_pt = 0
824 : return
825 : end if
826 0 : if (k == 1) then
827 0 : interp_xa_to_pt = xa(j,k)
828 0 : return
829 : end if
830 0 : if (k > 2 .and. k < sz) then
831 : ierr = 0
832 : call interp_4_to_1( &
833 : 0.5d0*(dq(k-2)+dq(k-1)), &
834 : 0.5d0*(dq(k-1)+dq(k)), &
835 : 0.5d0*(dq(k)+dq(k+1)), &
836 : 0.5d0*dq(k-2)+dq(k-1), &
837 : xa(j,k-2), xa(j,k-1), xa(j,k), xa(j,k+1), &
838 0 : interp_xa_to_pt, str, ierr)
839 0 : interp_xa_to_pt = min(1d0,max(0d0,interp_xa_to_pt))
840 0 : if (ierr == 0) return
841 : end if
842 0 : interp_xa_to_pt = (xa(j,k)*dq(k-1) + xa(j,k-1)*dq(k))/(dq(k-1) + dq(k))
843 0 : interp_xa_to_pt = min(1d0,max(0d0,interp_xa_to_pt))
844 0 : end function interp_xa_to_pt
845 :
846 :
847 98 : real(dp) function get_dtau1(s, ierr)
848 : type (star_info), pointer :: s
849 : integer, intent(out) :: ierr
850 : integer :: k
851 : real(dp) :: kap
852 : include 'formats'
853 98 : ierr = 0
854 98 : k = 1
855 98 : kap = s% opacity(1)
856 98 : get_dtau1 = s% dm(1)*kap/(pi4*s% rmid(1)*s% rmid(1))
857 98 : if (is_bad(get_dtau1)) then
858 0 : ierr = -1
859 0 : if (.not. s% report_ierr) return
860 : k = 1
861 0 : write(*,2) 'get_dtau1', k, get_dtau1
862 0 : write(*,2) 's% dm(1)', k, s% dm(k)
863 0 : write(*,2) 's% opacity(1)', k, s% opacity(k)
864 0 : write(*,2) 's% rmid(1)', k, s% rmid(k)
865 0 : write(*,2) 's% r(1)', k, s% r(k)
866 0 : write(*,2) 's% r(2)', 2, s% r(2)
867 0 : call mesa_error(__FILE__,__LINE__,'get_dtau1')
868 : end if
869 : end function get_dtau1
870 :
871 :
872 77 : subroutine get_tau(s, ierr)
873 : type (star_info), pointer :: s
874 : integer, intent(out) :: ierr
875 : ! tau(k) is optical depth at outer boundary of cell k
876 : real(dp) :: dtau, kap, dm_sum, L_sum
877 : integer :: k
878 : logical, parameter :: dbg = .true.
879 : include 'formats'
880 : ierr = 0
881 77 : dtau = get_dtau1(s, ierr)
882 77 : if (ierr /= 0) return
883 77 : s% tau(1) = s% tau_factor*s% tau_base
884 77 : dm_sum = 0
885 77 : L_sum = 0
886 93081 : do k = 2, s% nz
887 93004 : s% tau(k) = s% tau(k-1) + dtau
888 93004 : kap = s% opacity(k)
889 93004 : dtau = s% dm(k)*kap/(pi4*s% rmid(k)*s% rmid(k))
890 93004 : if (is_bad(dtau)) then
891 0 : ierr = -1
892 0 : if (.not. s% report_ierr) return
893 0 : write(*,2) 'dtau', k, dtau
894 0 : write(*,2) 's% dm(k)', k, s% dm(k)
895 0 : write(*,2) 's% opacity(k)', k, s% opacity(k)
896 0 : write(*,2) 's% rmid(k)', k, s% rmid(k)
897 0 : call mesa_error(__FILE__,__LINE__,'get_tau')
898 : end if
899 93081 : if (k == s% nz) s% tau_center = s% tau(k) + dtau
900 : !write(*,*) 'dtau, dlogtau', k, tau(k) - tau(k-1), &
901 : ! log10(tau(k)/tau(k-1))
902 : end do
903 : end subroutine get_tau
904 :
905 :
906 : integer function find_cell_for_mass(s, m)
907 : type (star_info), pointer :: s
908 : real(dp), intent(in) :: m
909 : integer :: k
910 : find_cell_for_mass = s% nz
911 : do k = 1, s% nz-1
912 : if (s% m(k) >= m .and. m > s% m(k+1)) then
913 : find_cell_for_mass = k
914 : return
915 : end if
916 : end do
917 : end function find_cell_for_mass
918 :
919 :
920 33 : subroutine get_delta_Pg_traditional(s, delta_Pg)
921 : type (star_info), pointer :: s
922 : real(dp), intent(out) :: delta_Pg ! seconds
923 : ! g-mode period spacing for l=1
924 : real(dp) :: dr, N2, integral, r
925 : integer :: k
926 : logical, parameter :: dbg = .false.
927 : include 'formats'
928 : if (dbg) then
929 : write(*,2) 's% star_mass', s% model_number, s% star_mass
930 : write(*,2) 's% photosphere_r', s% model_number, s% photosphere_r
931 : write(*,2) 's% Teff', s% model_number, s% Teff
932 : end if
933 33 : delta_Pg = 0._dp
934 33 : if (.not. s% calculate_Brunt_N2) return
935 33 : k = 0
936 33 : r = 0._dp
937 33 : N2 = 0._dp
938 33 : dr = 0._dp
939 33 : integral = 0._dp
940 :
941 : ! we integrate at cell edges.
942 40733 : do k = 2, s% nz
943 40700 : N2 = s% brunt_N2(k) ! brunt_N2 at cell_face
944 40700 : r = s% r(k) ! r evaluated at cell_face
945 40700 : dr = s% rmid(k-1) - s% rmid(k) ! dr evaluated at cell face.
946 40733 : if (N2 > 0d0) integral = integral + sqrt(N2)*dr/r
947 : end do
948 :
949 : if (dbg) write(*,2) ' integral ', &
950 : s% model_number,integral
951 :
952 33 : if (integral == 0) return
953 33 : delta_Pg = sqrt(2._dp)*pi*pi/integral
954 33 : if (is_bad(delta_Pg)) delta_Pg = 0._dp
955 :
956 : if (dbg) write(*,2) 'delta_Pg', s% model_number, delta_Pg
957 :
958 : end subroutine get_delta_Pg_traditional
959 :
960 0 : subroutine get_delta_Pg_bildsten2012(s, nu_max, delta_Pg)
961 : type (star_info), pointer :: s
962 : real(dp), intent(in) :: nu_max ! microHz
963 : real(dp), intent(out) :: delta_Pg ! seconds
964 : ! g-mode period spacing for l=1
965 : real(dp) :: integral, N2, omega2, kr2, L2, el, &
966 : dr, r, r2, cs2, sl2, I_integral, I_integral_limit
967 : integer :: k, k_sl2
968 : logical, parameter :: dbg = .false.
969 : include 'formats'
970 : if (dbg) then
971 : write(*,2) 'nu_max', s% model_number, nu_max
972 : write(*,2) 's% star_mass', s% model_number, s% star_mass
973 : write(*,2) 's% photosphere_r', s% model_number, s% photosphere_r
974 : write(*,2) 's% Teff', s% model_number, s% Teff
975 : end if
976 0 : delta_Pg = 0
977 0 : if (.not. s% calculate_Brunt_N2) return
978 0 : integral = 0
979 0 : I_integral = 0
980 0 : I_integral_limit = 0.5d0
981 0 : omega2 = pow2(2*pi*nu_max/1d6)
982 : if (dbg) write(*,1) 'log omega2', log10(omega2)
983 0 : el = 1
984 0 : L2 = el*(el+1)
985 0 : k_sl2 = 0
986 0 : do k = 2, s% nz
987 0 : N2 = s% brunt_N2(k)
988 0 : r = s% r(k)
989 0 : r2 = r*r
990 0 : cs2 = s% csound_face(k)*s% csound_face(k)
991 0 : sl2 = L2*cs2/r2
992 0 : dr = s% rmid(k-1) - s% rmid(k)
993 0 : if (omega2 >= sl2) then
994 : cycle
995 : end if
996 : if (k_sl2 == 0) then
997 : k_sl2 = k
998 : if (dbg) write(*,2) 'k_sl2', k
999 : end if
1000 0 : if (N2 > omega2) then ! in g-cavity
1001 : if (dbg .and. integral == 0) write(*,2) 'enter g-cavity', k
1002 0 : integral = integral + sqrt(N2)*dr/r
1003 : else ! in decay region
1004 0 : if (integral == 0) cycle ! ! haven't been in g-cavity yet
1005 : if (dbg .and. I_integral == 0) write(*,2) 'enter decay', k
1006 : ! in decay region below g-cavity; I_integral estimates decay
1007 0 : kr2 = (1 - n2/omega2)*(1 - Sl2/omega2)*omega2/cs2
1008 0 : I_integral = I_integral + sqrt(-kr2)*dr
1009 0 : if (I_integral > I_integral_limit) exit
1010 : end if
1011 : end do
1012 :
1013 : if (dbg) write(*,2) 'omega2 nu_max integral I_integral', &
1014 : s% model_number, omega2, nu_max, integral, I_integral
1015 :
1016 0 : if (integral == 0) return
1017 0 : delta_Pg = sqrt(2d0)*pi*pi/integral
1018 0 : if (is_bad(delta_Pg)) delta_Pg = 0
1019 :
1020 : if (dbg) write(*,2) 'delta_Pg', s% model_number, delta_Pg
1021 :
1022 : end subroutine get_delta_Pg_bildsten2012
1023 :
1024 :
1025 66 : subroutine set_rmid(s, nzlo, nzhi, ierr)
1026 : type (star_info), pointer :: s
1027 : integer, intent(in) :: nzlo, nzhi
1028 : integer, intent(out) :: ierr
1029 : integer :: k, nz
1030 : real(dp) :: r003, rp13, rmid, rmid2
1031 : logical :: dbg
1032 : include 'formats'
1033 66 : ierr = 0
1034 66 : nz = s% nz
1035 66 : dbg = .false.
1036 66 : if (s% RSP_flag) then ! mid = avg r
1037 : !dbg = s% model_number >= s% max_model_number - 1
1038 0 : do k=nzlo, nzhi
1039 0 : if (k < nz) then
1040 0 : rmid = 0.5d0*(s% r(k) + s% r(k+1))
1041 : else
1042 0 : rmid = 0.5d0*(s% r(k) + s% R_center)
1043 : end if
1044 0 : s% rmid(k) = rmid
1045 : if (dbg) write(*,3) 'set_rmid s% r(k)', k, s% model_number, s% r(k)
1046 : if (dbg) write(*,3) 'set_rmid s% rmid(k)', k, s% model_number, s% rmid(k)
1047 0 : if (s% rmid_start(k) < 0) s% rmid_start(k) = s% rmid(k)
1048 0 : rmid2 = rmid*rmid
1049 : end do
1050 66 : return
1051 : end if
1052 : ! mid = middle by volume
1053 80060 : do k=nzlo, nzhi
1054 79994 : r003 = s% r(k)*s% r(k)*s% r(k)
1055 79994 : if (k < nz) then
1056 79928 : rp13 = s% r(k+1)*s% r(k+1)*s% r(k+1)
1057 : else
1058 66 : rp13 = s% R_center*s% R_center*s% R_center
1059 : end if
1060 79994 : rmid = pow(0.5d0*(r003 + rp13),one_third)
1061 79994 : s% rmid(k) = rmid
1062 79994 : if (s% rmid_start(k) < 0) s% rmid_start(k) = s% rmid(k)
1063 80060 : rmid2 = rmid*rmid
1064 : end do
1065 : end subroutine set_rmid
1066 :
1067 :
1068 0 : real(dp) function get_tau_at_r(s, r, ierr)
1069 : type (star_info), pointer :: s
1070 : real(dp), intent(in) :: r
1071 : integer, intent(out) :: ierr
1072 : real(dp) :: dtau, tau_m1, tau_00
1073 : integer :: k
1074 : logical, parameter :: dbg = .false.
1075 : include 'formats'
1076 : ierr = 0
1077 0 : dtau = get_dtau1(s, ierr)
1078 0 : if (ierr /= 0) return
1079 0 : tau_00 = s% tau_factor*s% tau_base
1080 0 : get_tau_at_r = tau_00
1081 0 : if (r >= s% r(1)) return
1082 0 : do k = 2, s% nz
1083 0 : tau_m1 = tau_00
1084 0 : tau_00 = tau_m1 + dtau
1085 0 : if (r < s% r(k-1) .and. r >= s% r(k)) then
1086 : get_tau_at_r = &
1087 0 : (tau_00*(s% r(k-1)-r) + tau_m1*(r-s% r(k)))/(s% r(k-1)-s% r(k))
1088 0 : return
1089 : end if
1090 0 : dtau = s% dm(k)*s% opacity(k)/(pi4*s% rmid(k)*s% rmid(k))
1091 : end do
1092 : end function get_tau_at_r
1093 :
1094 :
1095 132 : subroutine set_phot_info(s)
1096 : use atm_lib, only: atm_black_body_T
1097 : type (star_info), pointer :: s
1098 : real(dp) :: luminosity
1099 : include 'formats'
1100 : call get_phot_info(s, &
1101 : s% photosphere_r, s% photosphere_m, s% photosphere_v, &
1102 : s% photosphere_L, s% photosphere_T, s% photosphere_csound, &
1103 : s% photosphere_opacity, s% photosphere_logg, &
1104 132 : s% photosphere_column_density, s% photosphere_cell_k)
1105 : s% photosphere_black_body_T = &
1106 132 : atm_black_body_T(s% photosphere_L, s% photosphere_r)
1107 132 : s% Teff = s% photosphere_black_body_T
1108 132 : s% photosphere_r = s% photosphere_r/Rsun
1109 132 : s% photosphere_m = s% photosphere_m/Msun
1110 132 : s% photosphere_L = s% photosphere_L/Lsun
1111 132 : s% L_phot = s% photosphere_L
1112 132 : luminosity = s% L(1)
1113 132 : if (is_bad(luminosity)) then
1114 0 : write(*,2) 's% L(1)', s% model_number, s% L(1)
1115 0 : write(*,2) 's% xh(s% i_lum,1)', s% model_number, s% xh(s% i_lum,1)
1116 0 : call mesa_error(__FILE__,__LINE__,'set_phot_info')
1117 : luminosity = 0d0
1118 : end if
1119 132 : s% L_surf = luminosity/Lsun
1120 132 : s% log_surface_luminosity = log10(max(1d-99,luminosity/Lsun))
1121 : ! log10(stellar luminosity in solar units)
1122 132 : if (is_bad(s% L_surf)) then
1123 0 : write(*,2) 's% L_surf', s% model_number, s% L_surf
1124 0 : call mesa_error(__FILE__,__LINE__,'set_phot_info')
1125 : end if
1126 132 : end subroutine set_phot_info
1127 :
1128 :
1129 132 : subroutine get_phot_info(s,r,m,v,L,T_phot,cs,kap,logg,ysum,k_phot)
1130 : type (star_info), pointer :: s
1131 : real(dp), intent(out) :: r, m, v, L, T_phot, cs, kap, logg, ysum
1132 : integer, intent(out) :: k_phot
1133 :
1134 : integer :: k
1135 : real(dp) :: tau00, taup1, dtau, r003, rp13, r3, tau_phot, &
1136 : Tface_0, Tface_1
1137 :
1138 : include 'formats'
1139 :
1140 : ! set values for surface as defaults in case phot not in model
1141 132 : r = s% r(1)
1142 132 : m = s% m(1)
1143 132 : if (s% u_flag) then
1144 0 : v = s% u(1)
1145 132 : else if (s% v_flag) then
1146 0 : v = s% v(1)
1147 : else
1148 132 : v = 0d0
1149 : end if
1150 132 : L = max(1d0, s% L(1)) ! don't use negative L(1)
1151 132 : T_phot = s% T(1)
1152 132 : cs = s% csound(1)
1153 132 : kap = s% opacity(1)
1154 132 : logg = safe_log10(s% cgrav(1)*m/(r*r))
1155 132 : k_phot = 1
1156 132 : if (s% tau_factor >= 1) then
1157 : ! This is always true for tables, regardless of tau_base.
1158 132 : return ! just use surface values
1159 : end if
1160 0 : tau_phot = s% tau_base ! this holds for case of tau_factor < 1
1161 0 : tau00 = s% tau_factor*s% tau_base ! start at tau_surf < tau_phot and go inward
1162 0 : taup1 = 0
1163 0 : ysum = 0
1164 0 : do k = 1, s% nz-1
1165 0 : dtau = s% dm(k)*s% opacity(k)/(pi4*s% rmid(k)*s% rmid(k))
1166 0 : taup1 = tau00 + dtau
1167 0 : ysum = ysum + s% rho(k)*(s% r(k) - s% r(k+1))
1168 0 : if (taup1 >= tau_phot .and. dtau > 0d0) then
1169 0 : if (k == 1) then
1170 0 : Tface_0 = s% T_surf
1171 : else
1172 0 : Tface_0 = 0.5d0*(s% T(k) + s% T(k-1))
1173 : end if
1174 0 : Tface_1 = 0.5d0*(s% T(k) + s% T(k+1))
1175 0 : T_phot = Tface_0 + (Tface_1 - Tface_0)*(tau_phot - tau00)/dtau
1176 0 : r003 = s% r(k)*s% r(k)*s% r(k)
1177 0 : rp13 = s% r(k+1)*s% r(k+1)*s% r(k+1)
1178 0 : r3 = r003 + (rp13 - r003)*(tau_phot - tau00)/dtau
1179 0 : r = pow(r3,one_third)
1180 0 : m = s% m(k) - s% dm(k)*(tau_phot - tau00)/dtau
1181 0 : if (s% u_flag) then
1182 0 : v = s% v_center
1183 : ! skip it since get_phot_info can be called before u_face has been set
1184 0 : else if (s% v_flag) then
1185 0 : v = s% v(k) + (s% v(k+1) - s% v(k))*(tau_phot - tau00)/dtau
1186 : end if
1187 : L = s% L(k) + (s% L(k+1) - s% L(k))*(tau_phot - tau00)/dtau
1188 0 : L = max(1d0, s% L(1)) ! don't use negative L(1)
1189 0 : logg = safe_log10(s% cgrav(k_phot)*m/(r*r))
1190 0 : k_phot = k
1191 : ! don't bother interpolating these.
1192 0 : cs = s% csound(k_phot)
1193 0 : kap = s% opacity(k_phot)
1194 0 : return
1195 : end if
1196 0 : tau00 = taup1
1197 : end do
1198 : !write(*,*) 'get_phot_info failed to find photosphere'
1199 0 : k_phot = s% nz
1200 0 : r = s% R_center
1201 0 : m = s% m_center
1202 0 : v = s% v_center
1203 0 : T_phot = s% T(k_phot)
1204 0 : L = max(1d0, s% L_center)
1205 0 : cs = s% csound(k_phot)
1206 0 : kap = s% opacity(k_phot)
1207 0 : logg = safe_log10(s% cgrav(k_phot)*m/(r*r))
1208 : end subroutine get_phot_info
1209 :
1210 :
1211 440 : real(dp) function center_value(s, p)
1212 : type (star_info), pointer :: s
1213 : real(dp), intent(in) :: p(:)
1214 : real(dp) :: sum_x, sum_dq, dx, dq
1215 : integer :: k
1216 440 : sum_x = 0
1217 440 : sum_dq = 0
1218 440 : do k = s% nz, 1, -1
1219 440 : dq = s% dq(k)
1220 440 : dx = p(k)*dq
1221 440 : if (sum_dq+dq >= s% center_avg_value_dq) then
1222 440 : sum_x = sum_x + dx*(s% center_avg_value_dq - sum_dq)/dq
1223 440 : sum_dq = s% center_avg_value_dq
1224 440 : exit
1225 : end if
1226 0 : sum_x = sum_x + dx
1227 0 : sum_dq = sum_dq + dq
1228 : end do
1229 440 : center_value = sum_x/sum_dq
1230 440 : end function center_value
1231 :
1232 :
1233 2464 : subroutine interp_q( &
1234 2464 : nz2, nvar_hydro, species, qval, xh, xa, q, dq, struct, comp, ierr)
1235 : use num_lib, only: binary_search
1236 : integer, intent(in) :: nz2, nvar_hydro, species
1237 : real(dp), intent(in) :: qval
1238 : real(dp), intent(in) :: xh(:,:), xa(:,:), q(:), dq(:)
1239 : real(dp), intent(inout) :: struct(:), comp(:)
1240 : integer, intent(out) :: ierr
1241 : integer :: k
1242 : real(dp) :: alfa
1243 2464 : ierr = 0
1244 2464 : if (qval <= q(nz2)) then
1245 1 : if (nvar_hydro > 0) &
1246 5 : struct(1:nvar_hydro) = xh(1:nvar_hydro,nz2)
1247 1 : if (species > 0) &
1248 9 : comp(1:species) = xa(1:species,nz2)
1249 1 : return
1250 : end if
1251 2463 : k = binary_search(nz2, q, 0, qval)
1252 2463 : if (k < 1 .or. k >= nz2) then
1253 0 : ierr = -1
1254 0 : return
1255 : end if
1256 2463 : if (qval <= q(k) .and. qval > q(k+1)) then
1257 2463 : alfa = (qval - q(k+1)) / dq(k)
1258 2463 : if (nvar_hydro > 0) &
1259 : struct(1:nvar_hydro) = &
1260 12315 : alfa*xh(1:nvar_hydro,k) + (1-alfa)*xh(1:nvar_hydro,k+1)
1261 2463 : if (species > 0) &
1262 22167 : comp(1:species) = alfa*xa(1:species,k) + (1-alfa)*xa(1:species,k+1)
1263 2463 : return
1264 : end if
1265 0 : ierr = -1
1266 : end subroutine interp_q
1267 :
1268 :
1269 33 : subroutine set_abs_du_div_cs(s)
1270 : type (star_info), pointer :: s
1271 :
1272 : integer :: k, nz, j
1273 : real(dp) :: abs_du, cs
1274 : include 'formats'
1275 33 : nz = s% nz
1276 :
1277 33 : if (s% v_flag) then
1278 0 : do k=2,nz
1279 0 : abs_du = abs(s% v_start(k) - s% v_start(k-1))
1280 0 : cs = maxval(s% csound(max(1,k-5):min(nz,k+5)))
1281 0 : s% abs_du_plus_cs(k) = abs_du + cs
1282 0 : s% abs_du_div_cs(k) = abs_du/cs
1283 : end do
1284 0 : k = 1
1285 0 : s% abs_du_plus_cs(k) = s% abs_du_plus_cs(k+1)
1286 0 : s% abs_du_div_cs(k) = s% abs_du_div_cs(k+1)
1287 0 : do j = 1,3
1288 0 : do k=2,nz-1
1289 0 : s% abs_du_div_cs(k) = sum(s% abs_du_div_cs(k-1:k+1))/3d0
1290 : end do
1291 : end do
1292 33 : else if (s% u_flag) then
1293 0 : do k=2,nz-1
1294 : abs_du = &
1295 : max(abs(s% u_start(k) - s% u_start(k+1)), &
1296 0 : abs(s% u_start(k) - s% u_start(k-1)))
1297 0 : cs = maxval(s% csound(max(1,k-5):min(nz,k+5)))
1298 0 : s% abs_du_plus_cs(k) = abs_du + cs
1299 0 : s% abs_du_div_cs(k) = abs_du/cs
1300 : end do
1301 0 : k = nz
1302 : s% abs_du_plus_cs(k) = &
1303 0 : abs(s% u_start(k) - s% u_start(k-1)) + s% csound_start(k)
1304 : s% abs_du_div_cs(k) = &
1305 0 : abs(s% u_start(k) - s% u_start(k-1))/s% csound_start(k)
1306 0 : k = 2
1307 : s% abs_du_plus_cs(k) = &
1308 0 : abs(s% u_start(k) - s% u_start(k+1)) + s% csound_start(k)
1309 : s% abs_du_div_cs(k) = &
1310 0 : abs(s% u_start(k) - s% u_start(k+1))/s% csound_start(k)
1311 0 : k = 1
1312 0 : s% abs_du_plus_cs(k) = s% abs_du_plus_cs(k+1)
1313 0 : s% abs_du_div_cs(k) = s% abs_du_div_cs(k+1)
1314 0 : do j = 1,3
1315 0 : do k=2,nz-1
1316 0 : s% abs_du_div_cs(k) = sum(s% abs_du_div_cs(k-1:k+1))/3d0
1317 : end do
1318 : end do
1319 : else
1320 40766 : do k=1,nz
1321 40733 : s% abs_du_plus_cs(k) = 1d99
1322 40766 : s% abs_du_div_cs(k) = 1d99
1323 : end do
1324 : end if
1325 :
1326 33 : end subroutine set_abs_du_div_cs
1327 :
1328 :
1329 33 : subroutine get_shock_info(s)
1330 : type (star_info), pointer :: s
1331 : integer :: k, nz
1332 : real(dp) :: v_div_cs_00, v_div_cs_m1, v_div_cs_min, v_div_cs_max, shock_radius
1333 : real(dp), pointer :: v(:)
1334 :
1335 : include 'formats'
1336 :
1337 33 : s% shock_mass = 0d0
1338 33 : s% shock_q = 0d0
1339 33 : s% shock_radius = 0d0
1340 33 : s% shock_velocity = 0d0
1341 33 : s% shock_csound = 0d0
1342 33 : s% shock_lgT = 0d0
1343 33 : s% shock_lgRho = 0d0
1344 33 : s% shock_lgP = 0d0
1345 33 : s% shock_gamma1 = 0d0
1346 33 : s% shock_entropy = 0d0
1347 33 : s% shock_tau = 0d0
1348 33 : s% shock_k = 0
1349 :
1350 33 : if (s% u_flag) then
1351 0 : v => s% u
1352 33 : else if (s% v_flag) then
1353 0 : v => s% v
1354 : else
1355 33 : return
1356 : end if
1357 :
1358 0 : nz = s% nz
1359 0 : shock_radius = -1
1360 0 : v_div_cs_00 = v(1)/s% csound(1)
1361 0 : do k = 2,nz-1
1362 0 : v_div_cs_m1 = v_div_cs_00
1363 0 : v_div_cs_00 = v(k)/s% csound(k)
1364 0 : v_div_cs_max = max(v_div_cs_00, v_div_cs_m1)
1365 0 : v_div_cs_min = min(v_div_cs_00, v_div_cs_m1)
1366 0 : if (v_div_cs_max >= 1d0 .and. v_div_cs_min < 1d0) then
1367 0 : if (v(k+1) > s% csound(k+1)) then ! skip single point glitches
1368 : shock_radius = &
1369 0 : find0(s% r(k), v_div_cs_00-1d0, s% r(k-1), v_div_cs_m1-1d0)
1370 0 : if (shock_radius <= 0d0) then
1371 0 : call mesa_error(__FILE__,__LINE__,'get_shock_info 1')
1372 : end if
1373 : exit
1374 : end if
1375 : end if
1376 0 : if (v_div_cs_min <= -1d0 .and. v_div_cs_max > -1d0) then
1377 0 : if (v(k+1) < -s% csound(k+1)) then ! skip single point glitches
1378 : shock_radius = &
1379 0 : find0(s% r(k), v_div_cs_00+1d0, s% r(k-1), v_div_cs_m1+1d0)
1380 0 : if (shock_radius <= 0d0) then
1381 0 : call mesa_error(__FILE__,__LINE__,'get_shock_info 2')
1382 : end if
1383 : exit
1384 : end if
1385 : end if
1386 : end do
1387 0 : if (shock_radius < 0d0) return
1388 :
1389 : call get_shock_location_info( &
1390 : s, .false., k-1, v, shock_radius, &
1391 : s% shock_mass, &
1392 : s% shock_q, &
1393 : s% shock_radius, &
1394 : s% shock_velocity, &
1395 : s% shock_csound, &
1396 : s% shock_lgT, &
1397 : s% shock_lgRho, &
1398 : s% shock_lgP, &
1399 : s% shock_gamma1, &
1400 : s% shock_entropy, &
1401 : s% shock_tau, &
1402 0 : s% shock_k)
1403 :
1404 : end subroutine get_shock_info
1405 :
1406 :
1407 0 : subroutine get_shock_location_info( &
1408 : s, dbg, k_shock, v, r, &
1409 : shock_mass, &
1410 : shock_q, &
1411 : shock_radius, &
1412 : shock_velocity, &
1413 : shock_csound, &
1414 : shock_lgT, &
1415 : shock_lgRho, &
1416 : shock_lgP, &
1417 : shock_gamma1, &
1418 : shock_entropy, &
1419 : shock_tau, &
1420 : shock_k)
1421 : type (star_info), pointer :: s
1422 : integer, intent(in) :: k_shock
1423 : logical, intent(in) :: dbg
1424 : real(dp), intent(in), pointer :: v(:)
1425 : real(dp), intent(in) :: r
1426 : real(dp), intent(out) :: &
1427 : shock_mass, &
1428 : shock_q, &
1429 : shock_radius, &
1430 : shock_velocity, &
1431 : shock_csound, &
1432 : shock_lgT, &
1433 : shock_lgRho, &
1434 : shock_lgP, &
1435 : shock_gamma1, &
1436 : shock_entropy, &
1437 : shock_tau
1438 : integer, intent(out) :: shock_k
1439 :
1440 : integer :: k
1441 : real(dp) :: alfa, beta
1442 :
1443 : include 'formats'
1444 :
1445 0 : k = k_shock
1446 : if (r < s% R_center .or. r > s% r(1) .or. &
1447 : k < 1 .or. k > s% nz .or. &
1448 0 : .not. associated(v) .or. &
1449 : .not. (s% v_flag .or. s% u_flag)) then
1450 0 : shock_mass = 0
1451 0 : shock_q = 0
1452 0 : shock_radius = 0
1453 0 : shock_velocity = 0
1454 0 : shock_csound = 0
1455 0 : shock_lgT = 0
1456 0 : shock_lgRho = 0
1457 0 : shock_lgP = 0
1458 0 : shock_gamma1 = 0
1459 0 : shock_entropy = 0
1460 0 : shock_tau = 0
1461 0 : shock_k = 0
1462 0 : return
1463 : end if
1464 :
1465 0 : shock_radius = r/Rsun
1466 0 : shock_k = k
1467 0 : if (k < s% nz) then
1468 0 : alfa = (r - s% r(k))/(s% r(k+1) - s% r(k))
1469 0 : beta = 1d0 - alfa
1470 0 : shock_mass = (alfa*s% m(k+1) + beta*s% m(k))/Msun
1471 0 : shock_q = alfa*s% q(k+1) + beta*s% q(k)
1472 0 : shock_velocity = alfa*v(k+1) + beta*v(k)
1473 : else
1474 0 : shock_mass = s% m(k)/Msun
1475 0 : shock_q = s% q(k)
1476 0 : shock_velocity = v(k)
1477 : end if
1478 0 : shock_csound = s% csound(k)
1479 0 : shock_lgT = s% lnT(k)/ln10
1480 0 : shock_lgRho = s% lnd(k)/ln10
1481 0 : shock_lgP = s% lnPeos(k)/ln10
1482 0 : shock_gamma1 = s% gamma1(k)
1483 0 : shock_entropy = s% entropy(k)
1484 0 : shock_tau = s% tau(k)
1485 :
1486 : end subroutine get_shock_location_info
1487 :
1488 :
1489 0 : real(dp) function min_dr_div_cs(s,min_k) ! seconds
1490 : type (star_info), pointer :: s
1491 : integer, intent(out) :: min_k
1492 : integer :: k, nz, k_min
1493 : real(dp) :: dr, dt, min_q, max_q, min_abs_u_div_cs, min_abs_du_div_cs, r00, rp1, dr_div_cs, remnant_mass
1494 : include 'formats'
1495 0 : nz = s% nz
1496 0 : min_k = nz
1497 0 : min_dr_div_cs = 1d99
1498 0 : min_q = s% min_q_for_dt_div_min_dr_div_cs_limit
1499 0 : max_q = s% max_q_for_dt_div_min_dr_div_cs_limit
1500 0 : k_min = max(1, s% min_k_for_dt_div_min_dr_div_cs_limit)
1501 0 : if (s% check_remnant_only_for_dt_div_min_dr_div_cs_limit) then
1502 0 : remnant_mass = get_remnant_mass(s)
1503 : else
1504 0 : remnant_mass = s% m(1)
1505 : end if
1506 : min_abs_u_div_cs = &
1507 0 : s% min_abs_u_div_cs_for_dt_div_min_dr_div_cs_limit
1508 : min_abs_du_div_cs = &
1509 0 : s% min_abs_du_div_cs_for_dt_div_min_dr_div_cs_limit
1510 0 : if (s% v_flag) then
1511 0 : do k = k_min, nz-1
1512 0 : if (s% m(k) > remnant_mass) cycle
1513 0 : if (s% q(k) > max_q) cycle
1514 0 : if (s% q(k) < min_q) exit
1515 0 : if (abs(s% v_start(k))/s% csound(k) < min_abs_u_div_cs) cycle
1516 0 : if (s% abs_du_div_cs(k) < min_abs_du_div_cs) cycle
1517 0 : r00 = s% r(k)
1518 0 : rp1 = s% r(k+1)
1519 0 : dr_div_cs = (r00 - rp1)/s% csound(k)
1520 0 : if (dr_div_cs < min_dr_div_cs) then
1521 0 : min_dr_div_cs = dr_div_cs
1522 0 : min_k = k
1523 : end if
1524 : end do
1525 : !write(*,3) 'min_dr_div_cs', min_k, s% model_number, min_dr_div_cs
1526 : return
1527 : end if
1528 0 : if (.not. s% u_flag) return
1529 0 : do k = k_min, nz-1
1530 0 : if (s% m(k) > remnant_mass) cycle
1531 0 : if (s% q(k) > max_q) cycle
1532 0 : if (s% q(k) < min_q) exit
1533 0 : if (abs(s% u_start(k))/s% csound(k) < min_abs_u_div_cs) cycle
1534 0 : if (s% abs_du_div_cs(k) < min_abs_du_div_cs) cycle
1535 0 : dr = s% r(k) - s% r(k+1)
1536 0 : dt = dr/s% abs_du_plus_cs(k)
1537 0 : if (dt < min_dr_div_cs) then
1538 0 : min_dr_div_cs = dt
1539 0 : min_k = k
1540 : end if
1541 : end do
1542 : end function min_dr_div_cs
1543 :
1544 :
1545 12 : subroutine reset_epsnuc_vectors(s)
1546 : type (star_info), pointer :: s
1547 : integer :: k, nz
1548 12 : nz = s% nz
1549 15563 : do k=1,s% nz
1550 15551 : s% eps_nuc(k) = 0d0
1551 15551 : s% d_epsnuc_dlnd(k) = 0d0
1552 15551 : s% d_epsnuc_dlnT(k) = 0d0
1553 139959 : s% d_epsnuc_dx(:,k) = 0d0
1554 388775 : s% eps_nuc_categories(:,k) = 0d0
1555 139959 : s% dxdt_nuc(:,k) = 0d0
1556 139959 : s% d_dxdt_nuc_dRho(:,k) = 0d0
1557 139959 : s% d_dxdt_nuc_dT(:,k) = 0d0
1558 1135223 : s% d_dxdt_nuc_dx(:,:,k) = 0d0
1559 15563 : s% eps_nuc_neu_total(k) = 0d0
1560 : end do
1561 12 : end subroutine reset_epsnuc_vectors
1562 :
1563 :
1564 44 : subroutine reset_starting_vectors(s)
1565 : type (star_info), pointer :: s
1566 : integer :: k, nz
1567 44 : nz = s% nz
1568 55336 : do k=1,s% nz
1569 55292 : s% T_start(k) = -1d99
1570 55292 : s% r_start(k) = -1d99
1571 55292 : s% rmid_start(k) = -1d99
1572 55292 : s% v_start(k) = -1d99
1573 55292 : s% u_start(k) = -1d99
1574 55292 : s% lnd_start(k) = -1d99
1575 55292 : s% lnT_start(k) = -1d99
1576 55292 : s% csound_start(k) = -1d99
1577 55292 : s% rho_start(k) = -1d99
1578 55292 : s% erad_start(k) = -1d99
1579 55292 : s% alpha_RTI_start(k) = -1d99
1580 55292 : s% opacity_start(k) = -1d99
1581 55292 : s% w_start(k) = -1d99
1582 55336 : s% dPdr_dRhodr_info(k) = -1d99
1583 : end do
1584 44 : end subroutine reset_starting_vectors
1585 :
1586 :
1587 104652 : subroutine save_eqn_dxa_partials(&
1588 104652 : s, k, nvar, i_eqn, species, dxam1, dxa00, dxap1, str, ierr)
1589 : type (star_info), pointer :: s
1590 : integer, intent(in) :: k, nvar, i_eqn, species
1591 : real(dp), intent(in), dimension(species) :: dxam1, dxa00, dxap1
1592 : character (len=*), intent(in) :: str
1593 : integer, intent(out) :: ierr
1594 : integer :: j
1595 104652 : ierr = 0
1596 941868 : do j=1,species
1597 837216 : call em1(s, i_eqn, j+s% nvar_hydro, k, nvar, dxam1(j))
1598 837216 : call e00(s, i_eqn, j+s% nvar_hydro, k, nvar, dxa00(j))
1599 941868 : call ep1(s, i_eqn, j+s% nvar_hydro, k, nvar, dxap1(j))
1600 : end do
1601 104652 : end subroutine save_eqn_dxa_partials
1602 :
1603 :
1604 104740 : subroutine save_eqn_residual_info(s, k, nvar, i_eqn, resid, str, ierr)
1605 : type (star_info), pointer :: s
1606 : integer, intent(in) :: k, nvar, i_eqn
1607 : type(auto_diff_real_star_order1), intent(in) :: resid
1608 : character (len=*), intent(in) :: str
1609 : integer, intent(out) :: ierr
1610 209480 : real(dp) :: d_dm1(nvar), d_d00(nvar), d_dp1(nvar)
1611 : call unpack_residual_partials(s, k, nvar, i_eqn, &
1612 104740 : resid, d_dm1, d_d00, d_dp1)
1613 : call store_partials( &
1614 104740 : s, k, i_eqn, nvar, d_dm1, d_d00, d_dp1, str, ierr)
1615 104740 : end subroutine save_eqn_residual_info
1616 :
1617 :
1618 209392 : subroutine unpack_residual_partials(s, k, nvar, i_eqn, &
1619 209392 : residual, d_dm1, d_d00, d_dp1)
1620 : use auto_diff
1621 : use auto_diff_support
1622 : type (star_info), pointer :: s
1623 : integer, intent(in) :: k, nvar, i_eqn
1624 : type(auto_diff_real_star_order1) :: residual
1625 : real(dp) :: d_dm1(nvar), d_d00(nvar), d_dp1(nvar)
1626 :
1627 : real(dp) :: val, dlnd_m1, dlnd_00, dlnd_p1, dlnT_m1, dlnT_00, dlnT_p1, &
1628 : dw_m1, dw_00, dw_p1, &
1629 : dlnR_m1, dlnR_00, dlnR_p1, &
1630 : dv_m1, dv_00, dv_p1, dL_m1, dL_00, dL_p1, &
1631 : dHp_m1, dHp_00, dHp_p1, &
1632 : dw_div_wc_m1, dw_div_wc_00, dw_div_wc_p1, &
1633 : djrot_m1, djrot_00, djrot_p1, &
1634 : dxtra1_m1, dxtra1_00, dxtra1_p1, &
1635 : dxtra2_m1, dxtra2_00, dxtra2_p1
1636 :
1637 : include 'formats'
1638 :
1639 : call unwrap(residual, val, &
1640 : dlnd_m1, dlnd_00, dlnd_p1, dlnT_m1, dlnT_00, dlnT_p1, &
1641 : dw_m1, dw_00, dw_p1, dlnR_m1, dlnR_00, dlnR_p1, &
1642 : dv_m1, dv_00, dv_p1, dL_m1, dL_00, dL_p1, &
1643 : dHp_m1, dHp_00, dHp_p1, &
1644 : dw_div_wc_m1, dw_div_wc_00, dw_div_wc_p1, &
1645 : djrot_m1, djrot_00, djrot_p1, &
1646 : dxtra1_m1, dxtra1_00, dxtra1_p1, &
1647 209392 : dxtra2_m1, dxtra2_00, dxtra2_p1)
1648 :
1649 8166288 : d_dm1 = 0; d_d00 = 0; d_dp1 = 0
1650 209392 : call unpack1(s% i_lnd, dlnd_m1, dlnd_00, dlnd_p1)
1651 209392 : call unpack1(s% i_lnT, dlnT_m1, dlnT_00, dlnT_p1)
1652 209392 : call unpack1(s% i_lnR, dlnR_m1, dlnR_00, dlnR_p1)
1653 209392 : if (s% i_v /= 0) call unpack1(s% i_v, dv_m1, dv_00, dv_p1)
1654 209392 : if (s% i_u /= 0) call unpack1(s% i_u, dv_m1, dv_00, dv_p1)
1655 209392 : if (s% i_lum /= 0) call unpack1(s% i_lum, dL_m1, dL_00, dL_p1)
1656 209392 : if (s% i_w /= 0) call unpack1(s% i_w, dw_m1, dw_00, dw_p1)
1657 209392 : if (s% i_Hp /= 0) call unpack1(s% i_Hp, dHp_m1, dHp_00, dHp_p1)
1658 209392 : if (s% i_w_div_wc /= 0) call unpack1(s% i_w_div_wc, dw_div_wc_m1, dw_div_wc_00, dw_div_wc_p1)
1659 209392 : if (s% i_j_rot /= 0) call unpack1(s% i_j_rot, djrot_m1, djrot_00, djrot_p1)
1660 :
1661 : contains
1662 :
1663 418784 : subroutine unpack1(j, dvar_m1, dvar_00, dvar_p1)
1664 : integer, intent(in) :: j
1665 : real(dp), intent(in) :: dvar_m1, dvar_00, dvar_p1
1666 418784 : d_dm1(j) = dvar_m1
1667 418784 : d_d00(j) = dvar_00
1668 209392 : d_dp1(j) = dvar_p1
1669 209392 : end subroutine unpack1
1670 :
1671 : end subroutine unpack_residual_partials
1672 :
1673 209392 : subroutine store_partials(s, k, i_eqn, nvar, d_dm1, d_d00, d_dp1, str, ierr)
1674 : type (star_info), pointer :: s
1675 : integer, intent(in) :: k, i_eqn, nvar
1676 : real(dp), intent(in) :: d_dm1(nvar), d_d00(nvar), d_dp1(nvar)
1677 : character (len=*), intent(in) :: str
1678 : integer, intent(out) :: ierr
1679 : integer :: nz, j
1680 : logical, parameter :: checking = .true.
1681 209392 : ierr = 0
1682 209392 : nz = s% nz
1683 2722096 : do j=1,nvar
1684 2512704 : if (k > 1) then
1685 2510592 : if (checking) call check_dequ(d_dm1(j),trim(str) // ' d_dm1')
1686 2510592 : call em1(s, i_eqn, j, k, nvar, d_dm1(j))
1687 : end if
1688 2512704 : if (checking) call check_dequ(d_d00(j),trim(str) // ' d_d00')
1689 2512704 : call e00(s, i_eqn, j, k, nvar, d_d00(j))
1690 2722096 : if (k < nz) then
1691 2510592 : if (checking) call check_dequ(d_dp1(j),trim(str) // ' d_dp1')
1692 2510592 : call ep1(s, i_eqn, j, k, nvar, d_dp1(j))
1693 : end if
1694 : end do
1695 :
1696 : contains
1697 :
1698 7533888 : subroutine check_dequ(dequ, str)
1699 : real(dp), intent(in) :: dequ
1700 : character (len=*), intent(in) :: str
1701 : include 'formats'
1702 7533888 : if (is_bad(dequ)) then
1703 0 : !$omp critical (store_partials_crit)
1704 0 : ierr = -1
1705 0 : if (s% report_ierr) then
1706 0 : write(*,2) 'store_partials: bad ' // trim(str), k, dequ
1707 : end if
1708 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'store_partials')
1709 : !$omp end critical (store_partials_crit)
1710 0 : return
1711 : end if
1712 : end subroutine check_dequ
1713 :
1714 : end subroutine store_partials
1715 :
1716 :
1717 77 : subroutine set_scale_height(s)
1718 : type (star_info), pointer :: s
1719 : real(dp) :: Hp, alt_Hp, alfa, beta, rho_face, Peos_face
1720 : integer :: k
1721 : include 'formats'
1722 93158 : do k=1,s% nz
1723 93081 : if (s% cgrav(k) == 0) then
1724 0 : s% scale_height(k) = 0
1725 0 : cycle
1726 : end if
1727 93081 : if (k == 1) then
1728 : alfa = 1
1729 : else
1730 93004 : alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
1731 : end if
1732 93081 : beta = 1 - alfa
1733 93004 : if (alfa == 1) then
1734 77 : rho_face = s% rho(k)
1735 77 : Peos_face = s% Peos(k)
1736 : else
1737 93004 : rho_face = alfa*s% rho(k) + beta*s% rho(k-1)
1738 93004 : Peos_face = alfa*s% Peos(k) + beta*s% Peos(k-1)
1739 : end if
1740 93081 : Hp = Peos_face/(rho_face*s% grav(k))
1741 93081 : if (s% cgrav(k) <= 0) then
1742 0 : alt_Hp = s% r(k)
1743 : else
1744 93081 : alt_Hp = sqrt(Peos_face / s% cgrav(k)) / rho_face
1745 : end if
1746 93158 : s% scale_height(k) = min(Hp, alt_Hp)
1747 : end do
1748 77 : end subroutine set_scale_height
1749 :
1750 :
1751 0 : real(dp) function tau_eff(s,k)
1752 : ! tau_eff = tau that gives the local P == P_atm if this location at surface
1753 : type (star_info), pointer :: s
1754 : integer, intent(in) :: k
1755 : real(dp) :: Peos, g, Pextra_factor
1756 0 : if (k == 1) then
1757 0 : tau_eff = s% tau(1)
1758 0 : return
1759 : end if
1760 0 : if (s% cgrav(k) <= 0d0) then
1761 0 : tau_eff = 0d0
1762 : return
1763 : end if
1764 0 : Peos = (s% dq(k-1)*s% Peos(k) + s% dq(k)*s% Peos(k-1))/(s% dq(k-1) + s% dq(k))
1765 0 : g = s% cgrav(k)*s% m_grav(k)/(s% r(k)*s% r(k))
1766 0 : Pextra_factor = s% Pextra_factor
1767 : tau_eff = s% opacity(k)*(Peos/g - &
1768 0 : Pextra_factor*(s% L(k)/s% m_grav(k))/(6d0*pi*clight*s% cgrav(k)))
1769 0 : end function tau_eff
1770 :
1771 :
1772 21 : real(dp) function eval_Ledd(s, ierr)
1773 : type (star_info), pointer :: s
1774 : integer, intent(out) :: ierr
1775 : real(dp) :: dtau1, dtau, tau, dqsum, Ledd_sum
1776 : integer :: k
1777 : logical, parameter :: dbg = .false.
1778 : include 'formats'
1779 21 : ierr = 0
1780 21 : eval_Ledd = 0d0
1781 21 : if (s% cgrav(1) <= 0d0) return
1782 21 : dtau1 = get_dtau1(s, ierr)
1783 21 : if (ierr /= 0) return
1784 21 : dtau = dtau1
1785 21 : tau = s% tau_factor*s% tau_base
1786 21 : dqsum = s% dq(1)
1787 21 : Ledd_sum = s% dq(1)*pi4*clight*s% cgrav(1)*s% m_grav(1)/s% opacity(1)
1788 2415 : do k = 2, s% nz
1789 2415 : tau = tau + dtau
1790 2415 : if (tau > s% surf_avg_tau) exit
1791 2394 : dtau = s% dm(k)*s% opacity(k)/(pi4*s% rmid(k)*s% rmid(k))
1792 2394 : dqsum = dqsum + s% dq(k)
1793 : Ledd_sum = Ledd_sum + &
1794 2415 : s% dq(k)*pi4*clight*s% cgrav(1)*s% m_grav(1)/s% opacity(k)
1795 : end do
1796 21 : eval_Ledd = Ledd_sum/dqsum
1797 21 : end function eval_Ledd
1798 :
1799 :
1800 10 : real(dp) function eval_min_cell_collapse_time(s,k_lo,k_hi,min_collapse_k,ierr) &
1801 : result(min_collapse_time)
1802 : type (star_info), pointer :: s
1803 : integer, intent(in) :: k_lo, k_hi
1804 : integer, intent(out) :: min_collapse_k, ierr
1805 : real(dp) :: rp1, vp1, r00, v00, time
1806 : integer :: k
1807 : logical, parameter :: dbg = .false.
1808 : include 'formats'
1809 10 : ierr = 0
1810 10 : rp1 = s% R_center
1811 10 : vp1 = s% v_center
1812 10 : min_collapse_time = 1d99
1813 10 : min_collapse_k = -1
1814 10 : if (.not. s% v_flag) return
1815 0 : do k = k_hi, k_lo, -1
1816 0 : v00 = s% v(k)
1817 0 : r00 = s% r(k)
1818 0 : if (r00 <= rp1) then
1819 : !ierr = -1
1820 0 : min_collapse_time = -1
1821 0 : min_collapse_k = -1
1822 0 : return
1823 : write(*,2) 'bad radii', k, r00, rp1
1824 : call mesa_error(__FILE__,__LINE__,'eval_min_cell_collapse_time')
1825 : end if
1826 0 : if (vp1 > v00) then
1827 0 : time = (r00 - rp1)/(vp1 - v00)
1828 0 : if (time < min_collapse_time) then
1829 0 : min_collapse_time = time
1830 0 : min_collapse_k = k
1831 : end if
1832 : end if
1833 0 : rp1 = r00
1834 0 : vp1 = v00
1835 : end do
1836 : end function eval_min_cell_collapse_time
1837 :
1838 :
1839 42 : real(dp) function total_angular_momentum(s) result(J)
1840 : type (star_info), pointer :: s
1841 : include 'formats'
1842 42 : if (.not. s% rotation_flag) then
1843 : J = 0
1844 : else
1845 0 : J = dot_product(s% dm_bar(1:s% nz), s% j_rot(1:s% nz))
1846 : end if
1847 42 : end function total_angular_momentum
1848 :
1849 :
1850 0 : real(dp) function eval_irradiation_heat(s,k)
1851 : type (star_info), pointer :: s
1852 : integer, intent(in) :: k
1853 : real(dp) :: irradiation_dq, xq, eps
1854 0 : eval_irradiation_heat = 0
1855 0 : if (s% irradiation_flux /= 0) then
1856 0 : irradiation_dq = pi4*s% r(1)*s% r(1)*s% column_depth_for_irradiation/s% xmstar
1857 0 : xq = 1 - s% q(k)
1858 0 : if (irradiation_dq > xq) then ! add irradiation heat for cell k
1859 0 : eps = 0.25d0 * s% irradiation_flux / s% column_depth_for_irradiation
1860 0 : if (irradiation_dq < xq + s% dq(k)) then ! only part of cell gets heated
1861 0 : eval_irradiation_heat = eps*(irradiation_dq - xq)/s% dq(k)
1862 : else ! all of cell gets heated
1863 : eval_irradiation_heat = eps
1864 : end if
1865 : end if
1866 : end if
1867 0 : end function eval_irradiation_heat
1868 :
1869 :
1870 0 : subroutine start_time(s, time0, total_all_before)
1871 : type (star_info), pointer :: s
1872 : integer(i8), intent(out) :: time0
1873 : real(dp), intent(out) :: total_all_before
1874 : integer(i8) :: clock_rate
1875 0 : if (.not. s% doing_timing) return
1876 0 : total_all_before = total_times(s)
1877 0 : call system_clock(time0,clock_rate)
1878 : end subroutine start_time
1879 :
1880 :
1881 0 : subroutine update_time(s, time0, total_all_before, total)
1882 : type (star_info), pointer :: s
1883 : integer(i8), intent(in) :: time0
1884 : real(dp), intent(in) :: total_all_before
1885 : real(dp), intent(inout) :: total
1886 : real(dp) :: total_all_after, other_stuff
1887 : integer(i8) :: time1, clock_rate
1888 0 : if (.not. s% doing_timing) return
1889 0 : call system_clock(time1,clock_rate)
1890 0 : total_all_after = total_times(s)
1891 0 : other_stuff = total_all_after - total_all_before
1892 : ! don't double count any other stuff
1893 0 : total = total + (dble(time1-time0)/clock_rate - other_stuff)
1894 : end subroutine update_time
1895 :
1896 :
1897 0 : real(dp) function total_times(s)
1898 : type (star_info), pointer :: s
1899 : total_times = &
1900 : s% time_evolve_step + &
1901 : s% time_remesh + &
1902 : s% time_adjust_mass + &
1903 : s% time_conv_premix + &
1904 : s% time_element_diffusion + &
1905 : s% time_struct_burn_mix + &
1906 : s% time_solver_matrix + &
1907 : s% time_solve_mix + &
1908 : s% time_solve_burn + &
1909 : s% time_solve_omega_mix + &
1910 : s% time_eos + &
1911 : s% time_neu_kap + &
1912 : s% time_nonburn_net + &
1913 : s% time_mlt + &
1914 : s% time_set_hydro_vars + &
1915 0 : s% time_set_mixing_info
1916 0 : end function total_times
1917 :
1918 :
1919 11 : real(dp) function get_remnant_mass(s)
1920 : type (star_info), pointer :: s
1921 11 : get_remnant_mass = s% m(1) - get_ejecta_mass(s)
1922 11 : end function get_remnant_mass
1923 :
1924 :
1925 22 : real(dp) function get_ejecta_mass(s)
1926 : use num_lib, only: find0
1927 : type (star_info), pointer :: s
1928 : integer :: k
1929 : real(dp) :: v, vesc, v_div_vesc, v_div_vesc_prev, dm
1930 : include 'formats'
1931 22 : get_ejecta_mass = 0d0
1932 22 : if (.not. (s% u_flag .or. s% v_flag)) return
1933 0 : v_div_vesc_prev = 0d0
1934 0 : do k=1,s% nz
1935 0 : if (s% u_flag) then
1936 : !v = s% u_face_ad(k)%val ! CANNOT USE u_face for this
1937 : ! approximate value is good enough for this estimate
1938 0 : if (k == 1) then
1939 0 : v = s% u(k)
1940 : else
1941 0 : v = 0.5d0*(s% u(k-1) + s% u(k))
1942 : end if
1943 : else
1944 0 : v = s% v(k)
1945 : end if
1946 0 : vesc = sqrt(2d0*s% cgrav(k)*s% m(k)/s% r(k))
1947 0 : v_div_vesc = v/vesc
1948 0 : if (v_div_vesc < 1d0) then
1949 0 : if (k == 1) return
1950 0 : dm = find0(0d0, v_div_vesc_prev-1d0, s% dm(k-1), v_div_vesc-1d0)
1951 0 : if (dm < 0d0) then
1952 0 : write(*,2) 'v_div_vesc_prev-1d0', k, v_div_vesc_prev-1d0
1953 0 : write(*,2) 'v_div_vesc-1d0', k, v_div_vesc-1d0
1954 0 : write(*,2) 's% dm(k-1)', k, s% dm(k-1)
1955 0 : write(*,2) 'dm', k, dm
1956 0 : call mesa_error(__FILE__,__LINE__,'get_ejecta_mass')
1957 : end if
1958 0 : if (k == 2) then
1959 : get_ejecta_mass = dm
1960 : else
1961 0 : get_ejecta_mass = sum(s% dm(1:k-2)) + dm
1962 : end if
1963 0 : return
1964 : end if
1965 0 : v_div_vesc_prev = v_div_vesc
1966 : end do
1967 : end function get_ejecta_mass
1968 :
1969 :
1970 220 : subroutine smooth(dc, sz)
1971 : real(dp), intent(inout) :: dc(:)
1972 : integer, intent(in) :: sz
1973 : integer :: k
1974 220 : k = 1
1975 220 : dc(k) = (3*dc(k) + dc(k+1))/4
1976 268880 : do k=2,sz-1
1977 268880 : dc(k) = (dc(k-1) + 2*dc(k) + dc(k+1))/4
1978 : end do
1979 220 : k = sz
1980 220 : dc(k) = (dc(k-1) + 3*dc(k))/4
1981 220 : end subroutine smooth
1982 :
1983 :
1984 0 : subroutine do_boxcar_mixing( &
1985 : s, min_mass, max_mass, boxcar_nominal_mass, number_iterations, ierr)
1986 : ! code based on SNEC
1987 : type (star_info), pointer :: s
1988 : real(dp), intent(in) :: min_mass, max_mass, boxcar_nominal_mass
1989 : integer, intent(in) :: number_iterations
1990 : integer, intent(out) :: ierr
1991 : integer :: nz, iter, k_inner, k_outer, j, k
1992 : real(dp) :: dm_sum, target_mass, actual_mass, xa, min_m, max_m
1993 : include 'formats'
1994 0 : ierr = 0
1995 0 : nz = s% nz
1996 0 : target_mass = boxcar_nominal_mass*Msun
1997 0 : min_m = min_mass*Msun
1998 0 : max_m = max_mass*Msun
1999 0 : write(*,2) 'iters min max box', number_iterations, min_mass, max_mass, boxcar_nominal_mass
2000 0 : do iter = 1, number_iterations
2001 0 : do k_inner = nz, 1, -1
2002 0 : if (s% m(k_inner) < min_m) cycle
2003 : dm_sum = 0
2004 0 : k_outer = 1
2005 0 : do k = k_inner, 1, -1
2006 0 : dm_sum = dm_sum + s% dm(k)
2007 0 : if (dm_sum > target_mass) then
2008 : k_outer = k
2009 : exit
2010 : end if
2011 : end do
2012 0 : if (s% m(k_outer) > max_m) cycle
2013 0 : actual_mass = dm_sum
2014 0 : do j=1,s% species
2015 : xa = 0d0
2016 0 : do k=k_outer,k_inner
2017 0 : xa = xa + s% xa(j,k)*s% dm(k)
2018 : end do
2019 0 : do k=k_outer,k_inner
2020 0 : s% xa(j,k) = xa/dm_sum
2021 : end do
2022 : end do
2023 0 : if (actual_mass < target_mass) exit
2024 : end do
2025 : end do
2026 0 : s% need_to_setvars = .true.
2027 0 : end subroutine do_boxcar_mixing
2028 :
2029 :
2030 79994 : subroutine get_XYZ(s, xa, X, Y, Z)
2031 : use chem_def, only: ih1, ih2, ihe3, ihe4
2032 : type (star_info), pointer :: s
2033 : real(dp), intent(in) :: xa(:)
2034 : real(dp), intent(out) :: X, Y, Z
2035 79994 : X = 0d0
2036 79994 : if (s% net_iso(ih1) /= 0) X = X + xa(s% net_iso(ih1))
2037 79994 : if (s% net_iso(ih2) /= 0) X = X + xa(s% net_iso(ih2))
2038 79994 : X = min(1d0, max(0d0, X))
2039 79994 : Y = 0d0
2040 79994 : if (s% net_iso(ihe3) /= 0) Y = Y + xa(s% net_iso(ihe3))
2041 79994 : if (s% net_iso(ihe4) /= 0) Y = Y + xa(s% net_iso(ihe4))
2042 79994 : Y = min(1d0, max(0d0, Y))
2043 79994 : Z = min(1d0, max(0d0, 1d0 - (X + Y)))
2044 79994 : end subroutine get_XYZ
2045 :
2046 :
2047 225 : subroutine get_face_values(s, v_mid, v_face, ierr)
2048 : ! simple interpolation by mass
2049 : type (star_info), pointer :: s
2050 : real(dp), intent(in) :: v_mid(:)
2051 : real(dp), intent(inout) :: v_face(:)
2052 : integer, intent(out) :: ierr
2053 : integer :: k
2054 : real(dp) :: dq_sum
2055 225 : ierr = 0
2056 225 : v_face(1) = v_mid(1)
2057 281420 : do k=2, s% nz
2058 281195 : dq_sum = s% dq(k-1) + s% dq(k)
2059 281420 : v_face(k) = (v_mid(k)*s% dq(k-1) + v_mid(k-1)*s% dq(k))/dq_sum
2060 : end do
2061 225 : end subroutine get_face_values
2062 :
2063 :
2064 0 : real(dp) function get_Ledd(s,k) result(Ledd)
2065 : type (star_info), pointer :: s
2066 : integer, intent(in) :: k
2067 : real(dp) :: kap_face
2068 : integer :: j
2069 0 : if (k == 1) then
2070 0 : j = 2
2071 : else
2072 0 : j = k
2073 : end if
2074 0 : kap_face = interp_val_to_pt(s% opacity,j,s% nz,s% dq,'get_Ledd')
2075 0 : Ledd = pi4*clight*s% cgrav(j)*s% m_grav(j)/kap_face
2076 0 : end function get_Ledd
2077 :
2078 :
2079 0 : real(dp) function get_Lrad(s,k)
2080 : type (star_info), pointer :: s
2081 : integer, intent(in) :: k
2082 0 : get_Lrad = s% L(k) - get_Lconv(s,k)
2083 0 : end function get_Lrad
2084 :
2085 :
2086 0 : real(dp) function get_effective_mixing_speed(s,k) result(v_mix)
2087 : type (star_info), pointer :: s
2088 : integer, intent(in) :: k
2089 : real(dp) :: lambda
2090 :
2091 0 : v_mix = 0d0
2092 0 : if (k == 1) return ! retain the zero surface value used for plotting
2093 0 : lambda = s% alpha_mlt(k)*s% scale_height(k)
2094 0 : if (lambda > 0d0) v_mix = 3d0*s% D_mix(k)/lambda
2095 : end function get_effective_mixing_speed
2096 :
2097 :
2098 0 : real(dp) function get_Lconv(s,k)
2099 : type (star_info), pointer :: s
2100 : integer, intent(in) :: k
2101 0 : if (k == 1) then
2102 0 : get_Lconv = 0d0
2103 : return
2104 : end if
2105 0 : if (s% RSP2_flag .or. s% RSP_flag) then
2106 0 : get_Lconv = s% Lc(k)
2107 : else
2108 0 : get_Lconv = s% L_conv(k) ! L_conv set by last call on mlt
2109 : end if
2110 : end function get_Lconv
2111 :
2112 :
2113 0 : real(dp) function get_Ladv(s,k)
2114 : type (star_info), pointer :: s
2115 : integer, intent(in) :: k
2116 : real(dp) :: T, Erad, v, r
2117 0 : T = s% T(k)
2118 0 : Erad = crad*T*T*T*T
2119 0 : if (s% u_flag) then
2120 0 : v = s% u(k)
2121 0 : else if (s% v_flag) then
2122 0 : v = s% v(k)
2123 : else
2124 : v = 0d0
2125 : end if
2126 0 : r = s% rmid(k)
2127 0 : get_Ladv = pi4*r*r*v*Erad
2128 0 : end function get_Ladv
2129 :
2130 :
2131 0 : real(dp) function get_Lrad_div_Ledd(s,k) result(L_rad_div_Ledd)
2132 : type (star_info), pointer :: s
2133 : integer, intent(in) :: k
2134 : integer :: j
2135 : real(dp) :: del_m, del_T4, area
2136 0 : if (s% cgrav(k) <= 0) then
2137 0 : L_rad_div_Ledd = 0d0
2138 : return
2139 : end if
2140 0 : if (k == 1) then
2141 : j = 2
2142 : else
2143 0 : j = k
2144 : end if
2145 0 : del_m = s% dm_bar(j)
2146 0 : del_T4 = pow4(s% T(j-1)) - pow4(s% T(j))
2147 0 : area = pi4*s% r(j)*s% r(j)
2148 : L_rad_div_Ledd = &
2149 0 : -(area*area*crad*(del_T4/del_m)/3)/(pi4*s% cgrav(j)*s% m_grav(j))
2150 0 : end function get_Lrad_div_Ledd
2151 :
2152 :
2153 : real(dp) function cell_start_specific_KE(s,k)
2154 : type (star_info), pointer :: s
2155 : integer, intent(in) :: k
2156 : cell_start_specific_KE = cell_start_specific_KE_qp(s,k)
2157 : end function cell_start_specific_KE
2158 :
2159 :
2160 52348 : real(qp) function cell_start_specific_KE_qp(s,k)
2161 : ! for consistency with dual cells at faces, use <v**2> instead of <v>**2
2162 : type (star_info), pointer :: s
2163 : integer, intent(in) :: k
2164 : real(qp) :: qhalf, v0, v1, v2, Mbar
2165 52348 : qhalf = 0.5d0
2166 52348 : if (s% use_mass_corrections) then
2167 0 : Mbar = s% mass_correction(k)
2168 : else
2169 : Mbar = 1.0_qp
2170 : end if
2171 52348 : if (s% u_flag) then
2172 0 : v0 = s% u_start(k)
2173 0 : cell_start_specific_KE_qp = qhalf*Mbar*v0**2
2174 52348 : else if (s% v_flag) then
2175 0 : v0 = s% v_start(k)
2176 0 : if (k < s% nz) then
2177 0 : v1 = s% v_start(k+1)
2178 : else
2179 0 : v1 = s% v_center
2180 : end if
2181 0 : v2 = qhalf*(v0**2 + v1**2)
2182 0 : cell_start_specific_KE_qp = qhalf*Mbar*v2
2183 : else ! ignore kinetic energy if no velocity variables
2184 : cell_start_specific_KE_qp = 0d0
2185 : end if
2186 52348 : end function cell_start_specific_KE_qp
2187 :
2188 :
2189 13567 : real(dp) function cell_specific_KE(s,k,d_dv00,d_dvp1)
2190 : type (star_info), pointer :: s
2191 : integer, intent(in) :: k
2192 : real(dp), intent(out) :: d_dv00, d_dvp1
2193 13567 : cell_specific_KE = cell_specific_KE_qp(s,k,d_dv00,d_dvp1)
2194 13567 : end function cell_specific_KE
2195 :
2196 :
2197 65915 : real(qp) function cell_specific_KE_qp(s,k,d_dv00,d_dvp1)
2198 : ! for consistency with dual cells at faces, use <v**2> instead of <v>**2
2199 : type (star_info), pointer :: s
2200 : integer, intent(in) :: k
2201 : real(dp), intent(out) :: d_dv00, d_dvp1
2202 : real(dp) :: dv2_dv00, dv2_dvp1
2203 : real(qp) :: qhalf, v0, v1, v2, Mbar
2204 65915 : qhalf = 0.5d0
2205 65915 : if (s% use_mass_corrections) then
2206 0 : Mbar = s% mass_correction(k)
2207 : else
2208 : Mbar = 1.0_qp
2209 : end if
2210 65915 : if (s% u_flag) then
2211 0 : v0 = s% u(k)
2212 0 : cell_specific_KE_qp = qhalf*Mbar*v0**2
2213 0 : d_dv00 = Mbar*s% u(k)
2214 0 : d_dvp1 = 0d0
2215 65915 : else if (s% v_flag) then
2216 0 : v0 = s% v(k)
2217 0 : if (k < s% nz) then
2218 0 : v1 = s% v(k+1)
2219 0 : dv2_dvp1 = s% v(k+1)
2220 : else
2221 0 : v1 = s% v_center
2222 0 : dv2_dvp1 = 0d0
2223 : end if
2224 0 : v2 = qhalf*(v0**2 + v1**2)
2225 0 : dv2_dv00 = s% v(k)
2226 0 : cell_specific_KE_qp = qhalf*Mbar*v2
2227 0 : d_dv00 = qhalf*Mbar*dv2_dv00
2228 0 : d_dvp1 = qhalf*Mbar*dv2_dvp1
2229 : else ! ignore kinetic energy if no velocity variables
2230 65915 : cell_specific_KE_qp = 0d0
2231 65915 : d_dv00 = 0d0
2232 65915 : d_dvp1 = 0d0
2233 : end if
2234 65915 : end function cell_specific_KE_qp
2235 :
2236 :
2237 68390 : real(dp) function cell_specific_PE(s,k,d_dlnR00,d_dlnRp1)
2238 : type (star_info), pointer :: s
2239 : integer, intent(in) :: k
2240 : real(dp), intent(out) :: d_dlnR00,d_dlnRp1
2241 68390 : cell_specific_PE = cell_specific_PE_qp(s,k,d_dlnR00,d_dlnRp1)
2242 68390 : end function cell_specific_PE
2243 :
2244 :
2245 120738 : real(qp) function cell_specific_PE_qp(s,k,d_dlnR00,d_dlnRp1)
2246 : ! for consistency with dual cells at faces, <m/r>_cntr => (m(k)/r(k) + m(k+1)/r(k+1))/2 /= m_cntr/r_cntr
2247 : ! i.e., use avg of m/r at faces of cell rather than ratio of cell center mass over cell center r.
2248 : type (star_info), pointer :: s
2249 : integer, intent(in) :: k
2250 : real(dp), intent(out) :: d_dlnR00,d_dlnRp1
2251 : real(qp) :: qhalf, rp1, r00, mp1, m00, Gp1, G00, gravp1, grav00, Mbar
2252 : real(dp) :: d_grav00_dlnR00, d_gravp1_dlnRp1
2253 : include 'formats'
2254 120738 : qhalf = 0.5d0
2255 120738 : if (s% use_mass_corrections) then
2256 0 : Mbar = s% mass_correction(k)
2257 : else
2258 : Mbar = 1.0_qp
2259 : end if
2260 120738 : if (k == s% nz) then
2261 99 : rp1 = s% R_center
2262 99 : mp1 = s% m_center
2263 99 : Gp1 = s% cgrav(s% nz)
2264 : else
2265 120639 : rp1 = s% r(k+1)
2266 120639 : mp1 = s% m_grav(k+1)
2267 120639 : Gp1 = s% cgrav(k+1)
2268 : end if
2269 120738 : if (rp1 <= 0d0) then
2270 99 : gravp1 = 0d0
2271 99 : d_gravp1_dlnRp1 = 0d0
2272 : else
2273 120639 : gravp1 = -Gp1*mp1/rp1
2274 120639 : d_gravp1_dlnRp1 = -gravp1
2275 : end if
2276 120738 : r00 = s% r(k)
2277 120738 : m00 = s% m_grav(k)
2278 120738 : G00 = s% cgrav(k)
2279 120738 : grav00 = -G00*m00/r00
2280 120738 : d_grav00_dlnR00 = -grav00
2281 120738 : cell_specific_PE_qp = qhalf*Mbar*(gravp1 + grav00)
2282 120738 : d_dlnR00 = qhalf*Mbar*d_grav00_dlnR00
2283 120738 : d_dlnRp1 = qhalf*Mbar*d_gravp1_dlnRp1
2284 120738 : if (is_bad(cell_specific_PE_qp)) then
2285 0 : write(*,2) 'cell_specific_PE_qp', k, cell_specific_PE_qp
2286 0 : write(*,2) 'gravp1', k, gravp1
2287 0 : write(*,2) 'grav00', k, grav00
2288 0 : call mesa_error(__FILE__,__LINE__,'cell_specific_PE')
2289 : end if
2290 120738 : end function cell_specific_PE_qp
2291 :
2292 :
2293 : real(dp) function cell_start_specific_PE(s,k)
2294 : type (star_info), pointer :: s
2295 : integer, intent(in) :: k
2296 : cell_start_specific_PE = cell_start_specific_PE_qp(s,k)
2297 : end function cell_start_specific_PE
2298 :
2299 :
2300 52348 : real(dp) function cell_start_specific_PE_qp(s,k)
2301 : ! for consistency with dual cells at faces, <m/r>_cntr => (m(k)/r(k) + m(k+1)/r(k+1))/2 /= m_cntr/r_cntr
2302 : ! i.e., use avg of m/r at faces of cell rather than ratio of cell center mass over cell center r.
2303 : type (star_info), pointer :: s
2304 : integer, intent(in) :: k
2305 : real(qp) :: qhalf, rp1, r00, mp1, m00, Gp1, G00, gravp1, grav00, Mbar
2306 : include 'formats'
2307 52348 : qhalf = 0.5d0
2308 52348 : if (s% use_mass_corrections) then
2309 0 : Mbar = s% mass_correction_start(k)
2310 : else
2311 : Mbar = 1.0_qp
2312 : end if
2313 52348 : if (k == s% nz) then
2314 44 : rp1 = s% R_center
2315 44 : mp1 = s% m_center
2316 44 : Gp1 = s% cgrav(s% nz)
2317 : else
2318 52304 : rp1 = s% r_start(k+1)
2319 52304 : mp1 = s% m_grav_start(k+1)
2320 52304 : Gp1 = s% cgrav(k+1)
2321 : end if
2322 52348 : if (rp1 <= 0d0) then
2323 44 : gravp1 = 0d0
2324 : else
2325 52304 : gravp1 = -Gp1*mp1/rp1
2326 : end if
2327 52348 : r00 = s% r_start(k)
2328 52348 : m00 = s% m_grav_start(k)
2329 52348 : G00 = s% cgrav(k)
2330 52348 : grav00 = -G00*m00/r00
2331 52348 : cell_start_specific_PE_qp = qhalf*Mbar*(gravp1 + grav00)
2332 52348 : if (is_bad(cell_start_specific_PE_qp)) then
2333 0 : write(*,2) 'cell_start_specific_PE_qp', k, cell_start_specific_PE_qp
2334 0 : write(*,2) 'gravp1', k, gravp1
2335 0 : write(*,2) 'grav00', k, grav00
2336 0 : call mesa_error(__FILE__,__LINE__,'cell_start_specific_PE_qp')
2337 : end if
2338 52348 : end function cell_start_specific_PE_qp
2339 :
2340 :
2341 0 : real(dp) function cell_specific_rotational_energy(s,k)
2342 : type (star_info), pointer :: s
2343 : integer, intent(in) :: k
2344 : real(dp) :: e_00, e_p1
2345 0 : e_00 = s% i_rot(k)% val*s% omega(k)*s% omega(k)
2346 0 : if (k < s% nz) then
2347 0 : e_p1 = s% i_rot(k+1)% val*s% omega(k+1)*s% omega(k+1)
2348 : else
2349 : e_p1 = 0
2350 : end if
2351 0 : cell_specific_rotational_energy = 0.5d0*(e_p1 + e_00)
2352 0 : end function cell_specific_rotational_energy
2353 :
2354 :
2355 52348 : subroutine get_dke_dt_dpe_dt(s, k, dt, &
2356 : dke_dt, d_dkedt_dv00, d_dkedt_dvp1, &
2357 : dpe_dt, d_dpedt_dlnR00, d_dpedt_dlnRp1, ierr)
2358 : type (star_info), pointer :: s
2359 : integer, intent(in) :: k
2360 : real(dp), intent(in) :: dt
2361 : real(dp), intent(out) :: &
2362 : dke_dt, d_dkedt_dv00, d_dkedt_dvp1, &
2363 : dpe_dt, d_dpedt_dlnR00, d_dpedt_dlnRp1
2364 : integer, intent(out) :: ierr
2365 : real(dp) :: PE_start, PE_new, KE_start, KE_new, q1
2366 : real(dp) :: dpe_dlnR00, dpe_dlnRp1, dke_dv00, dke_dvp1
2367 : include 'formats'
2368 52348 : ierr = 0
2369 : ! rate of change in specific PE (erg/g/s)
2370 52348 : PE_start = cell_start_specific_PE_qp(s,k)
2371 52348 : PE_new = cell_specific_PE_qp(s,k,dpe_dlnR00,dpe_dlnRp1)
2372 52348 : q1 = PE_new - PE_start
2373 52348 : dpe_dt = q1/dt ! erg/g/s
2374 52348 : d_dpedt_dlnR00 = dpe_dlnR00/dt
2375 52348 : d_dpedt_dlnRp1 = dpe_dlnRp1/dt
2376 : ! rate of change in specific KE (erg/g/s)
2377 52348 : KE_start = cell_start_specific_KE_qp(s,k)
2378 52348 : KE_new = cell_specific_KE_qp(s,k,dke_dv00,dke_dvp1)
2379 52348 : q1 = KE_new - KE_start
2380 52348 : dke_dt = q1/dt ! erg/g/s
2381 52348 : d_dkedt_dv00 = dke_dv00/dt
2382 52348 : d_dkedt_dvp1 = dke_dvp1/dt
2383 52348 : end subroutine get_dke_dt_dpe_dt
2384 :
2385 :
2386 : real(dp) function eval_deltaM_total_from_profile( &
2387 : deltaM, premesh_dm, profile)
2388 : real(dp), intent(in) :: deltaM
2389 : real(dp), intent(in) :: premesh_dm(:), profile(:)
2390 : real(dp) :: sum_dm, dm, total, f
2391 : integer :: k
2392 : include 'formats'
2393 : total = 0
2394 : sum_dm = 0
2395 : do k=1,size(premesh_dm,dim=1)
2396 : if (sum_dm >= deltaM) exit
2397 : dm = premesh_dm(k)
2398 : if (sum_dm + dm > deltaM) then
2399 : f = (deltaM - sum_dm)/dm
2400 : total = total + f*profile(k)
2401 : exit
2402 : end if
2403 : total = total + profile(k)
2404 : sum_dm = sum_dm + dm
2405 : end do
2406 : eval_deltaM_total_from_profile = total
2407 : end function eval_deltaM_total_from_profile
2408 :
2409 :
2410 11 : real(dp) function cell_specific_total_energy(s, k) result(cell_total)
2411 : type (star_info), pointer :: s
2412 : integer, intent(in) :: k
2413 : real(dp) :: d_dv00,d_dvp1,d_dlnR00,d_dlnRp1
2414 : include 'formats'
2415 11 : cell_total = s% energy(k)
2416 11 : if (s% v_flag .or. s% u_flag) &
2417 0 : cell_total = cell_total + cell_specific_KE(s,k,d_dv00,d_dvp1)
2418 11 : cell_total = cell_total + cell_specific_PE(s,k,d_dlnR00,d_dlnRp1)
2419 11 : if (s% rotation_flag .and. s% include_rotation_in_total_energy) &
2420 0 : cell_total = cell_total + cell_specific_rotational_energy(s,k)
2421 11 : if (s% RSP2_flag) cell_total = cell_total + pow2(s% w(k))
2422 11 : if (s% rsp_flag) cell_total = cell_total + s% RSP_Et(k)
2423 11 : end function cell_specific_total_energy
2424 :
2425 :
2426 0 : subroutine eval_integrated_total_energy_profile(s, arr, direction, ierr)
2427 : type (star_info), pointer :: s
2428 : integer, intent(in) :: direction
2429 : integer, intent(out) :: ierr
2430 : real(dp), intent(out) :: arr(:)
2431 : integer :: k, start, finish
2432 :
2433 0 : ierr = 0
2434 0 : if (direction == 1) then
2435 0 : start = 1
2436 0 : finish = s%nz
2437 0 : else if (direction == -1) then
2438 0 : start = s%nz
2439 0 : finish = 1
2440 : else
2441 0 : ierr = 1
2442 0 : call mesa_error(__FILE__,__LINE__)
2443 : end if
2444 :
2445 0 : arr(start) = cell_specific_total_energy(s,start) * s%dm(start)
2446 0 : do k=start+direction, finish, direction
2447 0 : arr(k) = arr(k-direction) + cell_specific_total_energy(s,k) * s%dm(k)
2448 : end do
2449 :
2450 0 : end subroutine eval_integrated_total_energy_profile
2451 :
2452 45 : subroutine eval_deltaM_total_energy_integrals( &
2453 : s, klo, khi, deltaM, save_profiles, &
2454 45 : total_energy_profile, &
2455 : total_internal_energy, total_gravitational_energy, &
2456 : total_radial_kinetic_energy, total_rotational_kinetic_energy, &
2457 : total_turbulent_energy, sum_total)
2458 : type (star_info), pointer :: s
2459 : integer, intent(in) :: klo, khi ! sum from klo to khi
2460 : real(dp), intent(in) :: deltaM
2461 : logical, intent(in) :: save_profiles
2462 : real(dp), intent(out), dimension(:) :: total_energy_profile
2463 : real(dp), intent(out) :: &
2464 : total_internal_energy, total_gravitational_energy, &
2465 : total_radial_kinetic_energy, total_rotational_kinetic_energy, &
2466 : total_turbulent_energy, sum_total
2467 : integer :: k
2468 : real(dp) :: dm, sum_dm, cell_total, cell1, d_dv00, d_dvp1, d_dlnR00, d_dlnRp1, alfa, beta,TDC_eturb_cell
2469 : include 'formats'
2470 :
2471 45 : total_internal_energy = 0d0
2472 45 : total_gravitational_energy = 0d0
2473 45 : total_radial_kinetic_energy = 0d0
2474 45 : total_rotational_kinetic_energy = 0d0
2475 45 : total_turbulent_energy = 0d0
2476 45 : sum_total = 0d0
2477 :
2478 45 : if (klo < 1 .or. khi > s% nz .or. klo > khi) return
2479 :
2480 45 : sum_dm = 0
2481 54857 : do k=klo,khi
2482 54812 : if (sum_dm >= deltaM) exit
2483 54812 : cell_total = 0
2484 54812 : dm = s% dm(k)
2485 54812 : if (sum_dm + dm > deltaM) dm = deltaM - sum_dm
2486 54812 : cell1 = dm*s% energy(k)
2487 54812 : cell_total = cell_total + cell1
2488 54812 : total_internal_energy = total_internal_energy + cell1
2489 54812 : if (s% v_flag .or. s% u_flag) then
2490 0 : cell1 = dm*cell_specific_KE(s,k,d_dv00,d_dvp1)
2491 0 : cell_total = cell_total + cell1
2492 0 : total_radial_kinetic_energy = total_radial_kinetic_energy + cell1
2493 : end if
2494 54812 : cell1 = dm*cell_specific_PE(s,k,d_dlnR00,d_dlnRp1)
2495 54812 : cell_total = cell_total + cell1
2496 54812 : total_gravitational_energy = total_gravitational_energy + cell1
2497 54812 : if (s% rotation_flag) then
2498 0 : cell1 = dm*cell_specific_rotational_energy(s,k)
2499 0 : total_rotational_kinetic_energy = total_rotational_kinetic_energy + cell1
2500 0 : if (s% include_rotation_in_total_energy) &
2501 0 : cell_total = cell_total + cell1
2502 : end if
2503 54812 : if (s% RSP2_flag) then
2504 0 : cell1 = dm*pow2(s% w(k))
2505 0 : cell_total = cell_total + cell1
2506 0 : total_turbulent_energy = total_turbulent_energy + cell1
2507 54812 : else if ( s% MLT_option == 'TDC' .and. &
2508 : s% TDC_include_eturb_in_energy_equation) then ! needs corrected s% mlt_vc(k) >= 0 causes failures.
2509 0 : if (k < s% nz) then
2510 : TDC_eturb_cell = 0.75d0*(pow2(s% mlt_vc(k)) + &
2511 0 : pow2(s% mlt_vc(k+1)))
2512 : else ! k == s% nz
2513 0 : TDC_eturb_cell = 0.75d0*pow2(s% mlt_vc(k))
2514 : end if
2515 0 : cell1 = dm*TDC_eturb_cell
2516 0 : cell_total = cell_total + cell1
2517 0 : total_turbulent_energy = total_turbulent_energy + cell1
2518 : end if
2519 54812 : if (s% rsp_flag) then
2520 0 : cell1 = dm*s% RSP_Et(k)
2521 0 : cell_total = cell_total + cell1
2522 0 : total_turbulent_energy = total_turbulent_energy + cell1
2523 : end if
2524 54857 : if (save_profiles) then
2525 13087 : total_energy_profile(k) = cell_total
2526 : end if
2527 : end do
2528 :
2529 : sum_total = total_internal_energy + total_gravitational_energy + &
2530 45 : total_radial_kinetic_energy + total_turbulent_energy
2531 :
2532 45 : if (s% include_rotation_in_total_energy) &
2533 0 : sum_total = sum_total + total_rotational_kinetic_energy
2534 :
2535 : end subroutine eval_deltaM_total_energy_integrals
2536 :
2537 :
2538 0 : subroutine eval_total_energy_profile(s, total_energy_profile)
2539 : type (star_info), pointer :: s
2540 : real(dp), intent(out), dimension(:) :: total_energy_profile
2541 :
2542 : integer :: k
2543 : real(dp) :: dm, cell_total, cell1, d_dv00, d_dvp1, d_dlnR00, d_dlnRp1
2544 : include 'formats'
2545 :
2546 0 : do k=1, s%nz
2547 0 : cell_total = 0
2548 0 : dm = s% dm(k)
2549 0 : cell1 = dm*s% energy(k)
2550 0 : cell_total = cell_total + cell1
2551 0 : if (s% v_flag .or. s% u_flag) then
2552 0 : cell1 = dm*cell_specific_KE(s,k,d_dv00,d_dvp1)
2553 0 : cell_total = cell_total + cell1
2554 : end if
2555 0 : cell1 = dm*cell_specific_PE(s,k,d_dlnR00,d_dlnRp1)
2556 0 : cell_total = cell_total + cell1
2557 0 : if (s% rotation_flag) then
2558 0 : cell1 = dm*cell_specific_rotational_energy(s,k)
2559 0 : if (s% include_rotation_in_total_energy) &
2560 0 : cell_total = cell_total + cell1
2561 : end if
2562 0 : if (s% RSP2_flag) then
2563 0 : cell1 = dm*pow2(s% w(k))
2564 0 : cell_total = cell_total + cell1
2565 : end if
2566 0 : if (s% rsp_flag) then
2567 0 : cell1 = dm*s% RSP_Et(k)
2568 0 : cell_total = cell_total + cell1
2569 : end if
2570 0 : total_energy_profile(k) = cell_total
2571 : end do
2572 :
2573 0 : end subroutine eval_total_energy_profile
2574 :
2575 :
2576 0 : real(dp) function eval_cell_section_total_energy( &
2577 : s, klo, khi) result(sum_total)
2578 : type (star_info), pointer :: s
2579 : integer, intent(in) :: klo, khi ! sum from klo to khi
2580 : real(dp) :: &
2581 : total_internal_energy, total_gravitational_energy, &
2582 : total_radial_kinetic_energy, total_rotational_kinetic_energy, &
2583 : total_turbulent_energy
2584 0 : real(dp), allocatable, dimension(:) :: total_energy_profile
2585 0 : allocate(total_energy_profile(1:s% nz))
2586 : call eval_deltaM_total_energy_integrals( &
2587 : s, klo, khi, s% mstar, .false., &
2588 : total_energy_profile, &
2589 : total_internal_energy, total_gravitational_energy, &
2590 : total_radial_kinetic_energy, total_rotational_kinetic_energy, &
2591 0 : total_turbulent_energy, sum_total)
2592 0 : end function eval_cell_section_total_energy
2593 :
2594 :
2595 33 : subroutine eval_total_energy_integrals(s, &
2596 : total_internal_energy, total_gravitational_energy, &
2597 : total_radial_kinetic_energy, total_rotational_kinetic_energy, &
2598 : total_turbulent_energy, sum_total)
2599 : type (star_info), pointer :: s
2600 : real(dp), intent(out) :: &
2601 : total_internal_energy, total_gravitational_energy, &
2602 : total_radial_kinetic_energy, total_rotational_kinetic_energy, &
2603 : total_turbulent_energy, sum_total
2604 33 : real(dp), allocatable, dimension(:) :: total_energy_profile
2605 99 : allocate(total_energy_profile(1:s% nz))
2606 : call eval_deltaM_total_energy_integrals( &
2607 : s, 1, s% nz, s% mstar, .false., &
2608 : total_energy_profile, &
2609 : total_internal_energy, total_gravitational_energy, &
2610 : total_radial_kinetic_energy, total_rotational_kinetic_energy, &
2611 33 : total_turbulent_energy, sum_total)
2612 33 : end subroutine eval_total_energy_integrals
2613 :
2614 :
2615 : subroutine eval_total_energy_integrals_and_profile(s, &
2616 : total_energy_profile, &
2617 : total_internal_energy, total_gravitational_energy, &
2618 : total_radial_kinetic_energy, total_rotational_kinetic_energy, &
2619 : total_turbulent_energy, sum_total)
2620 : type (star_info), pointer :: s
2621 : real(dp), intent(out) :: total_energy_profile(:), &
2622 : total_internal_energy, total_gravitational_energy, &
2623 : total_radial_kinetic_energy, total_rotational_kinetic_energy, &
2624 : total_turbulent_energy, sum_total
2625 : call eval_deltaM_total_energy_integrals( &
2626 : s, 1, s% nz, s% mstar, .true., &
2627 : total_energy_profile, &
2628 : total_internal_energy, total_gravitational_energy, &
2629 : total_radial_kinetic_energy, total_rotational_kinetic_energy, &
2630 : total_turbulent_energy, sum_total)
2631 : end subroutine eval_total_energy_integrals_and_profile
2632 :
2633 :
2634 : real(dp) function get_total_energy_integral(s,k) result(sum_total)
2635 : ! from surface down to k
2636 : type (star_info), pointer :: s
2637 : integer, intent(in) :: k
2638 : real(dp) :: &
2639 : total_internal_energy, total_gravitational_energy, &
2640 : total_radial_kinetic_energy, total_rotational_kinetic_energy, &
2641 : total_turbulent_energy
2642 : real(dp), allocatable, dimension(:) :: total_energy_profile
2643 : allocate(total_energy_profile(1:s% nz))
2644 : call eval_deltaM_total_energy_integrals( &
2645 : s, 1, k, s% mstar, .false., &
2646 : total_energy_profile, &
2647 : total_internal_energy, total_gravitational_energy, &
2648 : total_radial_kinetic_energy, total_rotational_kinetic_energy, &
2649 : total_turbulent_energy, sum_total)
2650 : end function get_total_energy_integral
2651 :
2652 :
2653 : real(dp) function get_total_energy_integral_outward(s,k) result(sum_total)
2654 : ! from surface down to k
2655 : type (star_info), pointer :: s
2656 : integer, intent(in) :: k
2657 : real(dp) :: &
2658 : total_internal_energy, total_gravitational_energy, &
2659 : total_radial_kinetic_energy, total_rotational_kinetic_energy, &
2660 : total_turbulent_energy
2661 : real(dp), allocatable, dimension(:) :: total_energy_profile
2662 : allocate(total_energy_profile(1:s% nz))
2663 : call eval_deltaM_total_energy_integrals( &
2664 : s, k, s% nz, s% mstar, .false., &
2665 : total_energy_profile, &
2666 : total_internal_energy, total_gravitational_energy, &
2667 : total_radial_kinetic_energy, total_rotational_kinetic_energy, &
2668 : total_turbulent_energy, sum_total)
2669 : end function get_total_energy_integral_outward
2670 :
2671 :
2672 0 : real(dp) function get_log_concentration(s,j,k) result(log_c)
2673 : use chem_def, only: chem_isos
2674 : type (star_info), pointer :: s
2675 : integer, intent(in) :: j, k
2676 : ! concentration = number density / number density of electrons
2677 : ! Ci = (Xi/Ai) / sum(Zi*Xi/Ai) [see Thoul et al, ApJ 421:828-842, 1994]
2678 : integer :: i, cid, species
2679 : real(dp) :: tmp, c
2680 0 : log_c = -1d99
2681 0 : if (s% chem_id(j) == 0) return
2682 0 : species = s% species
2683 0 : tmp = 0d0
2684 0 : do i=1,species
2685 0 : cid = s% chem_id(i)
2686 0 : tmp = tmp + chem_isos% Z(cid)*s% xa(i,k)/chem_isos% Z_plus_N(cid)
2687 : end do
2688 0 : cid = s% chem_id(j)
2689 0 : c = (s% xa(j,k)/chem_isos% Z_plus_N(cid))/tmp
2690 0 : log_c = safe_log10(c)
2691 0 : end function get_log_concentration
2692 :
2693 :
2694 0 : real(dp) function get_phi_Joss(s,k) result(phi)
2695 : use eos_def, only: i_lnPgas
2696 : ! Joss, Salpeter, Ostriker, 1973. density inversion when Lrad/Ledd > phi.
2697 : type (star_info), pointer :: s
2698 : integer, intent(in) :: k
2699 0 : phi = 1d0/(1d0 + (s% Pgas(k)/(4* s% Prad(k)))*s% d_eos_dlnT(i_lnPgas,k))
2700 0 : end function get_phi_Joss
2701 :
2702 :
2703 0 : logical function after_He_burn(s, he4_limit)
2704 : use chem_def
2705 : type (star_info), pointer :: s
2706 : real(dp), intent(in) :: he4_limit
2707 : integer :: nz, h1, he4
2708 : real(dp), parameter :: small = 1d-4
2709 0 : after_He_burn = .false.
2710 0 : nz = s% nz
2711 0 : h1 = s% net_iso(ih1)
2712 0 : he4 = s% net_iso(ihe4)
2713 0 : if (h1 == 0 .or. he4 == 0) return
2714 0 : if (s% xa(h1,nz) > small .or. s% xa(he4,nz) > he4_limit) return
2715 0 : after_He_burn = .true.
2716 : end function after_He_burn
2717 :
2718 :
2719 0 : logical function after_C_burn(s, c12_limit)
2720 : use chem_def
2721 : type (star_info), pointer :: s
2722 : real(dp), intent(in) :: c12_limit
2723 : integer :: nz, h1, he4, c12
2724 : real(dp), parameter :: small = 1d-4
2725 0 : after_C_burn = .false.
2726 0 : nz = s% nz
2727 0 : h1 = s% net_iso(ih1)
2728 0 : he4 = s% net_iso(ihe4)
2729 0 : c12 = s% net_iso(ic12)
2730 0 : if (h1 == 0 .or. he4 == 0 .or. c12 == 0) return
2731 0 : if (s% xa(h1,nz) > small .or. s% xa(he4,nz) > small .or. &
2732 0 : s% xa(c12,nz) > c12_limit) return
2733 0 : after_C_burn = .true.
2734 : end function after_C_burn
2735 :
2736 :
2737 0 : integer function lookup_nameofvar(s, namestr)
2738 : type (star_info), pointer :: s
2739 : character (len=*), intent(in) :: namestr
2740 : integer :: i
2741 0 : lookup_nameofvar = 0
2742 0 : do i=1,s% nvar_total
2743 0 : if (namestr == s% nameofvar(i)) then
2744 0 : lookup_nameofvar = i
2745 : return
2746 : end if
2747 : end do
2748 : end function lookup_nameofvar
2749 :
2750 :
2751 0 : integer function lookup_nameofequ(s, namestr)
2752 : type (star_info), pointer :: s
2753 : character (len=*), intent(in) :: namestr
2754 : integer :: i
2755 0 : lookup_nameofequ = 0
2756 0 : do i=1,s% nvar_total
2757 0 : if (namestr == s% nameofequ(i)) then
2758 0 : lookup_nameofequ = i
2759 : return
2760 : end if
2761 : end do
2762 : end function lookup_nameofequ
2763 :
2764 :
2765 0 : real(dp) function omega_crit(s, k) ! result always is > 0
2766 : type (star_info), pointer :: s
2767 : integer, intent(in) :: k
2768 : real(dp) :: Ledd, gamma_factor, Lrad_div_Ledd, rmid, r003, rp13
2769 : include 'formats'
2770 0 : if (s% rotation_flag) then
2771 : ! Use equatorial radius (at center of cell by volume)
2772 0 : r003 = s% r_equatorial(k)*s% r_equatorial(k)*s% r_equatorial(k)
2773 0 : if (k < s% nz) then
2774 0 : rp13 = s% r_equatorial(k+1)*s% r_equatorial(k+1)*s% r_equatorial(k+1)
2775 : else
2776 0 : rp13 = s% R_center*s% R_center*s% R_center
2777 : end if
2778 0 : rmid = pow(0.5d0*(r003 + rp13),one_third)
2779 : else
2780 0 : rmid = s% rmid(k)
2781 : end if
2782 :
2783 0 : Ledd = pi4*clight*s% cgrav(k)*s% m_grav(k)/s% opacity(k)
2784 0 : Lrad_div_Ledd = get_Lrad_div_Ledd(s,k)
2785 0 : gamma_factor = 1d0 - min(Lrad_div_Ledd, 0.9999d0)
2786 0 : omega_crit = sqrt(gamma_factor*s% cgrav(k)*s% m_grav(k)/pow3(rmid))
2787 0 : end function omega_crit
2788 :
2789 :
2790 89 : subroutine weighted_smoothing(dd, n, ns, preserve_sign, ddold)
2791 : ! based on routine written by S.-C. Yoon, 18 Sept. 2002
2792 : ! for smoothing any variable (dd) with size n over 2*ns+1 cells.
2793 : real(dp), intent(inout) :: dd(:) ! (n)
2794 : integer, intent(in) :: n, ns
2795 : logical, intent(in) :: preserve_sign
2796 : real(dp), intent(inout) :: ddold(:) ! (n) work array
2797 :
2798 : integer :: nweight, mweight, i, j, k
2799 178 : real(dp) :: weight(2*ns+1), sweight, v0
2800 :
2801 : include 'formats'
2802 :
2803 110193 : do i = 1,n
2804 110193 : ddold(i) = dd(i)
2805 : end do
2806 :
2807 : !--preparation for smoothing --------
2808 : nweight = ns
2809 : mweight = 2*nweight+1
2810 622 : do i = 1,mweight
2811 622 : weight(i) = 0d0
2812 : end do
2813 89 : weight(1) = 1d0
2814 533 : do i = 1,mweight-1
2815 1907 : do j = i+1,2,-1
2816 1818 : weight(j) = weight(j) + weight(j-1)
2817 : end do
2818 : end do
2819 :
2820 : !--smoothing ------------------------
2821 110015 : do i=2,n-1
2822 109926 : sweight=0d0
2823 109926 : dd(i)=0d0
2824 109926 : v0 = ddold(i)
2825 493259 : do j = i, max(1,i-nweight), -1
2826 383333 : k=j-i+nweight+1
2827 383333 : if (preserve_sign .and. v0*ddold(j) <= 0) exit
2828 383333 : sweight = sweight+weight(k)
2829 493259 : dd(i) = dd(i)+ddold(j)*weight(k)
2830 : end do
2831 383333 : do j = i+1, min(n,i+nweight)
2832 273407 : k=j-i+nweight+1
2833 273407 : if (preserve_sign .and. v0*ddold(j) <= 0) exit
2834 273407 : sweight = sweight+weight(k)
2835 383333 : dd(i) = dd(i)+ddold(j)*weight(k)
2836 : end do
2837 110015 : if (sweight > 0) then
2838 109926 : sweight = 1d0/sweight
2839 109926 : dd(i) = dd(i)*sweight
2840 : end if
2841 : end do
2842 :
2843 89 : end subroutine weighted_smoothing
2844 :
2845 :
2846 89 : subroutine threshold_smoothing(dd, dd_thresh, n, ns, preserve_sign, ddold)
2847 :
2848 : ! Same as weighted_smoothing, but only smooth contiguous regions where |dd| >= dd_thresh
2849 :
2850 : real(dp), intent(inout) :: dd(:) ! (n)
2851 : real(dp), intent(in) :: dd_thresh
2852 : integer, intent(in) :: n
2853 : integer, intent(in) :: ns
2854 : logical, intent(in) :: preserve_sign
2855 : real(dp), intent(inout) :: ddold(:) ! (n) work array
2856 :
2857 : logical :: in_region
2858 : integer :: i
2859 : integer :: i_a
2860 : integer :: i_b
2861 :
2862 : include 'formats'
2863 :
2864 : ! Process regions
2865 :
2866 89 : in_region = .FALSE.
2867 :
2868 89 : i_a = 1
2869 110193 : do i = 1, n
2870 :
2871 110193 : if (in_region) then
2872 :
2873 110015 : if (ABS(dd(i)) < dd_thresh) then
2874 0 : i_b = i-1
2875 0 : if (i_b > i_a) call weighted_smoothing(dd(i_a:i_b), i_b-i_a+1, ns, preserve_sign, ddold(i_a:i_b))
2876 : in_region = .FALSE.
2877 : end if
2878 :
2879 : else
2880 89 : if (ABS(dd(i)) >= dd_thresh) then
2881 110104 : i_a = i
2882 110104 : in_region = .TRUE.
2883 : end if
2884 :
2885 : end if
2886 :
2887 : end do
2888 :
2889 : ! Handle the final region
2890 :
2891 89 : if (in_region) then
2892 :
2893 89 : i_b = n
2894 89 : if (i_b > i_a) call weighted_smoothing(dd(i_a:i_b), i_b-i_a+1, ns, preserve_sign, ddold(i_a:i_b))
2895 :
2896 : end if
2897 :
2898 89 : return
2899 :
2900 : end subroutine threshold_smoothing
2901 :
2902 :
2903 33 : real(dp) function eval_kh_timescale(G,M,R,L) result(kh)
2904 : real(dp), intent(in) :: G,M,R,L
2905 33 : if (L <= 0) then
2906 : kh = 0d0
2907 : else
2908 33 : kh = 0.75d0*G*M*M/(R*L) ! 0.75 is based on sun. Hansen & Kawaler eqn 1.30
2909 : end if
2910 33 : end function eval_kh_timescale
2911 :
2912 :
2913 1 : real(dp) function yrs_for_init_timestep(s)
2914 : type (star_info), pointer :: s
2915 1 : if (s% initial_mass <= 1) then
2916 : yrs_for_init_timestep = 1d5
2917 : else
2918 1 : yrs_for_init_timestep = 1d5 / pow(s% initial_mass,2.5d0)
2919 : end if
2920 1 : end function yrs_for_init_timestep
2921 :
2922 :
2923 12 : subroutine set_phase_of_evolution(s) ! from evolve after call do_report
2924 : use rates_def, only: i_rate
2925 : use chem_def
2926 : type (star_info), pointer :: s
2927 : real(dp) :: center_h1, center_he4
2928 : integer :: nz, j
2929 : include 'formats'
2930 12 : nz = s% nz
2931 :
2932 12 : s% phase_of_evolution = phase_starting
2933 12 : if (s%doing_first_model_of_run) return
2934 :
2935 10 : j = s% net_iso(ih1)
2936 10 : if (j > 0) then
2937 10 : center_h1 = center_avg_x(s,j)
2938 : else
2939 : center_h1 = 1d99
2940 : end if
2941 10 : j = s% net_iso(ihe4)
2942 10 : if (j > 0) then
2943 10 : center_he4 = center_avg_x(s,j)
2944 : else
2945 : center_he4 = 1d99
2946 : end if
2947 :
2948 10 : if (s% doing_relax) then
2949 0 : s% phase_of_evolution = phase_relax
2950 10 : else if (s% photosphere_logg > 6d0) then
2951 0 : s% phase_of_evolution = phase_WDCS
2952 10 : else if (s% L_by_category(i_burn_si) > 1d2) then
2953 0 : s% phase_of_evolution = phase_Si_Burn
2954 10 : else if ( (s% L_by_category(ioo) + s% L_by_category(ico)) > 1d2 .and. center_he4 < 1d-4) then
2955 0 : s% phase_of_evolution = phase_O_Burn
2956 10 : else if (s% L_by_category(i_burn_ne) > 1d2 .and. center_he4 < 1d-4) then
2957 0 : s% phase_of_evolution = phase_Ne_Burn
2958 10 : else if (s% L_by_category(icc) > 1d2 .and. center_he4 < 1d-4) then
2959 0 : s% phase_of_evolution = phase_C_Burn
2960 : else if (center_he4 < 1d-4 .and. &
2961 40 : s% he_core_mass - s% co_core_mass <= 0.1d0 .and. &
2962 : any(s% burn_he_conv_region(1:s% num_conv_boundaries))) then
2963 0 : s% phase_of_evolution = phase_TP_AGB
2964 10 : else if (center_he4 <= 1d-4) then
2965 0 : s% phase_of_evolution = phase_TACHeB
2966 10 : else if (s% center_eps_burn(i3alf) > 1d2) then
2967 0 : s% phase_of_evolution = phase_ZACHeB
2968 10 : else if (s% L_by_category(i3alf) > 1d2) then
2969 0 : s% phase_of_evolution = phase_He_Burn
2970 10 : else if (center_h1 <= 1d-6) then
2971 0 : s% phase_of_evolution = phase_TAMS
2972 10 : else if (center_h1 <= 0.3d0) then
2973 0 : s% phase_of_evolution = phase_IAMS
2974 10 : else if (s% L_nuc_burn_total >= s% L_phot*s% Lnuc_div_L_zams_limit) then
2975 10 : s% phase_of_evolution = phase_ZAMS
2976 0 : else if (s% log_center_temperature > 5d0) then
2977 0 : s% phase_of_evolution = phase_PreMS
2978 : else
2979 : s% phase_of_evolution = phase_starting
2980 : end if
2981 :
2982 : end subroutine set_phase_of_evolution
2983 :
2984 :
2985 79994 : subroutine set_rv_info(s,k)
2986 : type (star_info), pointer :: s
2987 : integer, intent(in) :: k
2988 : include 'formats'
2989 79994 : if (s% v_flag) then
2990 0 : if (s% using_velocity_time_centering) then
2991 0 : s% vc(k) = 0.5d0*(s% v_start(k) + s% v(k))
2992 : else
2993 0 : s% vc(k) = s% v(k)
2994 : end if
2995 : else
2996 79994 : s% vc(k) = 0d0
2997 : end if
2998 79994 : end subroutine set_rv_info
2999 :
3000 :
3001 0 : subroutine show_matrix(s, dmat, nvar)
3002 : type (star_info), pointer :: s
3003 : integer, intent(in) :: nvar
3004 : real(dp) :: dmat(nvar,nvar)
3005 : integer :: i, j
3006 0 : write(*,'(A)')
3007 0 : write(*,'(18x)', advance = 'no')
3008 0 : do j = 1, nvar
3009 0 : write(*,'(a15)', advance = 'no') s% nameofvar(j)
3010 : end do
3011 0 : write(*,'(A)')
3012 0 : do i = 1, nvar
3013 0 : write(*,'(a15)', advance = 'no') s% nameofequ(i)
3014 0 : do j = 1, nvar
3015 0 : write(*,'(e13.4,2x)', advance = 'no') dmat(i,j)
3016 : end do
3017 0 : write(*,'(A)')
3018 : end do
3019 0 : write(*,'(A)')
3020 0 : end subroutine show_matrix
3021 :
3022 :
3023 : ! e00(i,j,k) is partial of equ(i,k) wrt var(j,k)
3024 8375328 : subroutine e00(s,i,j,k,nvar,v)
3025 : type (star_info), pointer :: s
3026 : integer, intent(in) :: i, j, k, nvar
3027 : real(dp), intent(in) :: v
3028 : logical, parameter :: dbg = .false.
3029 : include 'formats'
3030 :
3031 : if (mdb .and. k==397 .and. v /= 0d0) &
3032 : write(*,4) 'e00(i,j,k) ' // &
3033 : trim(s% nameofequ(i)) // ' ' // trim(s% nameofvar(j)), i, j, k, v
3034 :
3035 : !if (j == s% i_lnd .and. k /= s% nz) return ! this variable is being held constant
3036 :
3037 8375328 : if (v == 0d0) return
3038 :
3039 : if (.false. .and. j == s% i_lnT .and. k == 30) then
3040 : write(*,4) 'e00(i,j,k) ' // &
3041 : trim(s% nameofequ(i)) // ' ' // trim(s% nameofvar(j)), i, j, k, v, s% x_scale(j,k)
3042 : end if
3043 :
3044 3641728 : if (is_bad(v)) then
3045 0 : !$omp critical (star_utils_e00_crit1)
3046 : write(*,4) 'e00(i,j,k) ' // &
3047 0 : trim(s% nameofequ(i)) // ' ' // trim(s% nameofvar(j)), i, j, k, v
3048 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'1 e00')
3049 : !$omp end critical (star_utils_e00_crit1)
3050 : end if
3051 :
3052 3641728 : if (i <= 0 .or. j <= 0 .or. k <= 0 .or. k > s% nz) then
3053 : write(*,4) 'bad i,j,k e00(i,j,k) ' // &
3054 0 : trim(s% nameofequ(i)) // ' ' // trim(s% nameofvar(j)), i, j, k, v
3055 0 : call mesa_error(__FILE__,__LINE__,'2 e00')
3056 : end if
3057 :
3058 3641728 : if (j > nvar) return ! hybrid
3059 :
3060 3641728 : if (i > nvar) then
3061 0 : !$omp critical (star_utils_e00_crit2)
3062 : write(*,5) 'bad i e00(i,j,k) ' // &
3063 0 : trim(s% nameofequ(i)) // ' ' // trim(s% nameofvar(j)), &
3064 0 : s% solver_iter, i, j, k, v
3065 0 : call mesa_error(__FILE__,__LINE__,'3 e00')
3066 : !$omp end critical (star_utils_e00_crit2)
3067 : end if
3068 :
3069 3641728 : if (abs(v) < 1d-250) return
3070 :
3071 3641728 : s% dblk(i,j,k) = s% dblk(i,j,k) + v*s% x_scale(j,k)
3072 :
3073 : end subroutine e00
3074 :
3075 :
3076 : ! em1(i,j,k) is partial of equ(i,k) wrt var(j,k-1)
3077 3766240 : subroutine em1(s,i,j,k,nvar,v)
3078 : type (star_info), pointer :: s
3079 : integer, intent(in) :: i, j, k, nvar
3080 : real(dp), intent(in) :: v
3081 : logical, parameter :: dbg = .false.
3082 3766240 : if (k == 1) return
3083 : include 'formats'
3084 :
3085 : if (mdb .and. k==397 .and. v /= 0d0) &
3086 : write(*,4) 'em1(i,j,k) ' // &
3087 : trim(s% nameofequ(i)) // ' ' // trim(s% nameofvar(j)), i, j, k, v
3088 :
3089 3765888 : if (v == 0d0) return
3090 :
3091 : if (.false. .and. j == s% i_lnT .and. k == 31) then
3092 : write(*,4) 'em1(i,j,k) ' // &
3093 : trim(s% nameofequ(i)) // ' ' // trim(s% nameofvar(j)), i, j, k, v, s% x_scale(j,k-1)
3094 : end if
3095 :
3096 887140 : if (is_bad(v)) then
3097 0 : !$omp critical (star_utils_em1_crit1)
3098 : write(*,4) 'em1(i,j,k) ' // &
3099 0 : trim(s% nameofequ(i)) // ' ' // trim(s% nameofvar(j)), i, j, k, v
3100 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'em1')
3101 : !$omp end critical (star_utils_em1_crit1)
3102 : end if
3103 :
3104 887140 : if (i <= 0 .or. j <= 0 .or. k <= 0 .or. k > s% nz) then
3105 : write(*,4) 'bad i,j,k em1(i,j,k) ' // &
3106 0 : trim(s% nameofequ(i)) // ' ' // trim(s% nameofvar(j)), i, j, k, v
3107 0 : call mesa_error(__FILE__,__LINE__,'em1')
3108 : end if
3109 :
3110 887140 : if (j > nvar) return ! hybrid
3111 :
3112 887140 : if (i > nvar) then
3113 : write(*,5) 'bad i em1(i,j,k) ' // &
3114 0 : trim(s% nameofequ(i)) // ' ' // trim(s% nameofvar(j)), &
3115 0 : s% solver_iter, i, j, k, v
3116 0 : call mesa_error(__FILE__,__LINE__,'em1')
3117 : end if
3118 :
3119 887140 : if (abs(v) < 1d-250) return
3120 :
3121 887140 : s% lblk(i,j,k) = s% lblk(i,j,k) + v*s% x_scale(j,k-1)
3122 :
3123 : end subroutine em1
3124 :
3125 :
3126 : ! ep1(i,j,k) is partial of equ(i,k) wrt var(j,k+1)
3127 3766240 : subroutine ep1(s,i,j,k,nvar,v)
3128 : type (star_info), pointer :: s
3129 : integer, intent(in) :: i, j, k, nvar
3130 : real(dp), intent(in) :: v
3131 : logical, parameter :: dbg = .false.
3132 : include 'formats'
3133 :
3134 : if (mdb .and. k==397 .and. v /= 0d0) &
3135 : write(*,4) 'ep1(i,j,k) ' // &
3136 : trim(s% nameofequ(i)) // ' ' // trim(s% nameofvar(j)), i, j, k, v
3137 :
3138 3766240 : if (v == 0d0) return
3139 :
3140 : if (.false. .and. j == s% i_lnT .and. k == 29) then
3141 : write(*,4) 'ep1(i,j,k) ' // &
3142 : trim(s% nameofequ(i)) // ' ' // trim(s% nameofvar(j)), i, j, k, v, s% x_scale(j,k+1)
3143 : end if
3144 :
3145 526718 : if (is_bad(v)) then
3146 0 : !$omp critical (star_utils_ep1_crit1)
3147 : write(*,4) 'ep1(i,j,k) ' // &
3148 0 : trim(s% nameofequ(i)) // ' ' // trim(s% nameofvar(j)), i, j, k, v
3149 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'ep1')
3150 : !$omp end critical (star_utils_ep1_crit1)
3151 : end if
3152 :
3153 526718 : if (i <= 0 .or. j <= 0 .or. k <= 0 .or. k > s% nz) then
3154 : write(*,4) 'bad i,j,k ep1(i,j,k) ' // &
3155 0 : trim(s% nameofequ(i)) // ' ' // trim(s% nameofvar(j)), i, j, k, v
3156 0 : call mesa_error(__FILE__,__LINE__,'ep1')
3157 : end if
3158 :
3159 526718 : if (j > nvar) return
3160 :
3161 526718 : if (i > nvar) then
3162 : write(*,5) 'bad i ep1(i,j,k) ' // &
3163 0 : trim(s% nameofequ(i)) // ' ' // trim(s% nameofvar(j)), &
3164 0 : s% solver_iter, i, j, k, v
3165 0 : call mesa_error(__FILE__,__LINE__,'ep1')
3166 : end if
3167 :
3168 526718 : if (abs(v) < 1d-250) return
3169 :
3170 526718 : s% ublk(i,j,k) = s% ublk(i,j,k) + v*s% x_scale(j,k+1)
3171 :
3172 : end subroutine ep1
3173 :
3174 :
3175 66 : real(dp) function current_min_xa_hard_limit(s) result(min_xa_hard_limit)
3176 : type (star_info), pointer :: s
3177 : real(dp) :: logTc, alfa
3178 66 : logTc = s% lnT(s% nz)/ln10
3179 66 : if (logTc <= s% logT_max_for_min_xa_hard_limit) then
3180 66 : min_xa_hard_limit = s% min_xa_hard_limit
3181 0 : else if (logTc >= s% logT_min_for_min_xa_hard_limit_for_highT) then
3182 0 : min_xa_hard_limit = s% min_xa_hard_limit_for_highT
3183 : else
3184 : alfa = (logTc - s% logT_max_for_min_xa_hard_limit) / &
3185 0 : (s% logT_min_for_min_xa_hard_limit_for_highT - s% logT_max_for_min_xa_hard_limit)
3186 : min_xa_hard_limit = &
3187 0 : alfa*s% min_xa_hard_limit_for_highT + (1d0 - alfa)*s% min_xa_hard_limit
3188 : end if
3189 66 : end function current_min_xa_hard_limit
3190 :
3191 :
3192 33 : real(dp) function current_sum_xa_hard_limit(s) result(sum_xa_hard_limit)
3193 : type (star_info), pointer :: s
3194 : real(dp) :: logTc, alfa
3195 : include 'formats'
3196 33 : logTc = s% lnT(s% nz)/ln10
3197 33 : if (logTc <= s% logT_max_for_sum_xa_hard_limit) then
3198 33 : sum_xa_hard_limit = s% sum_xa_hard_limit
3199 0 : else if (logTc >= s% logT_min_for_sum_xa_hard_limit_for_highT) then
3200 0 : sum_xa_hard_limit = s% sum_xa_hard_limit_for_highT
3201 : else
3202 : alfa = (logTc - s% logT_max_for_sum_xa_hard_limit) / &
3203 0 : (s% logT_min_for_sum_xa_hard_limit_for_highT - s% logT_max_for_sum_xa_hard_limit)
3204 : sum_xa_hard_limit = &
3205 0 : alfa*s% sum_xa_hard_limit_for_highT + (1d0 - alfa)*s% sum_xa_hard_limit
3206 : end if
3207 : !write(*,1) 'logTc hard_limit', logTc, sum_xa_hard_limit
3208 33 : end function current_sum_xa_hard_limit
3209 :
3210 :
3211 0 : integer function no_extra_profile_columns(id)
3212 : use star_def, only: star_info
3213 : integer, intent(in) :: id
3214 0 : no_extra_profile_columns = 0
3215 0 : end function no_extra_profile_columns
3216 :
3217 :
3218 0 : subroutine no_data_for_extra_profile_columns(id, n, nz, names, vals, ierr)
3219 : use star_def, only: maxlen_profile_column_name, star_info
3220 : integer, intent(in) :: id, n, nz
3221 : character (len=maxlen_profile_column_name) :: names(n)
3222 : real(dp) :: vals(nz,n)
3223 : integer, intent(out) :: ierr
3224 0 : ierr = 0
3225 0 : end subroutine no_data_for_extra_profile_columns
3226 :
3227 :
3228 : ! Stiriba, Youssef, Appl, Numer. Math. 45, 499-511. 2003.
3229 :
3230 : ! LPP-HARMOD -- local piecewise parabolic reconstruction
3231 :
3232 : ! interpolant is derived to conserve integral of v in cell k
3233 : ! interpolant slope at cell midpoint is harmonic mean of slopes between adjacent cells
3234 : ! where these slopes between cells are defined as the difference between cell averages
3235 : ! divided by the distance between cell midpoints.
3236 : ! interpolant curvature based on difference between the midpoint slope
3237 : ! and the smaller in magnitude of the slopes between adjacent cells.
3238 :
3239 : ! interpolant f(dq) = a + b*dq + (c/2)*dq^2, with dq = q - q_midpoint
3240 : ! c0 holds a's, c1 holds b's, and c2 holds c's.
3241 :
3242 :
3243 : subroutine get1_lpp(k, nz, &
3244 : dm1, d00, dp1, vm1, v00, vp1, quadratic, monotonic, dbg, c0, c1, c2)
3245 : integer, intent(in) :: k, nz
3246 : real(dp), intent(in) :: dm1, d00, dp1, vm1, v00, vp1
3247 : logical, intent(in) :: quadratic, monotonic
3248 : logical, intent(in) :: dbg
3249 : real(dp), intent(out) :: c0, c1, c2
3250 : real(dp) :: vbdy1, vbdy2, dhalf, sm1, s00, sp1, sprod, dv1, dv2
3251 : logical :: okay
3252 : include 'formats'
3253 :
3254 : c0 = v00
3255 : c2 = 0d0
3256 : if (k == 1) then
3257 : sm1 = 0d0
3258 : sp1 = (v00 - vp1) / (0.5d0*(d00 + dp1))
3259 : c1 = sp1
3260 : else if (k == nz) then
3261 : sp1 = 0d0
3262 : sm1 = (vm1 - v00) / (0.5d0*(dm1 + d00))
3263 : c1 = sm1
3264 : if (dbg) then
3265 : write(*,2) 'get1_lpp c1', k, c1
3266 : end if
3267 : else
3268 : sm1 = (vm1 - v00) / (0.5d0*(dm1 + d00))
3269 : sp1 = (v00 - vp1) / (0.5d0*(d00 + dp1))
3270 : sprod = sm1*sp1
3271 : if (sprod <= 0d0) then
3272 : ! at local min or max, so set slope and curvature to 0.
3273 : c1 = 0d0
3274 : return
3275 : end if
3276 : if (.not. quadratic) then ! minmod
3277 : s00 = (vm1 - vp1)/(d00 + (dm1 + dp1)/2)
3278 : c1 = sm1
3279 : if (sm1*s00 < 0d0) c1 = 0d0
3280 : if (abs(s00) < abs(c1)) c1 = s00
3281 : if (s00*sp1 < 0d0) c1 = 0d0
3282 : if (abs(sp1) < abs(c1)) c1 = sp1
3283 : else
3284 : c1 = sprod*2/(sp1 + sm1) ! harmonic mean slope
3285 : if (abs(sm1) <= abs(sp1)) then
3286 : c2 = (sm1 - c1)/(2d0*d00)
3287 : else
3288 : c2 = (c1 - sp1)/(2d0*d00)
3289 : end if
3290 : c0 = v00 - c2*d00*d00/24d0
3291 : end if
3292 : end if
3293 :
3294 : if (.not. monotonic) return
3295 :
3296 : dhalf = 0.5d0*d00
3297 : dv1 = c1*dhalf
3298 : dv2 = 0.5d0*c2*dhalf*dhalf
3299 : vbdy1 = c0 + dv1 + dv2 ! value at face(k)
3300 : vbdy2 = c0 - dv1 + dv2 ! value at face(k+1)
3301 :
3302 : okay = .true.
3303 :
3304 : if (k > 1) then ! check v00 <= vbdy1 <= vm1 or v00 >= vbdy1 >= vm1
3305 : if ((vm1 - vbdy1)*(vbdy1 - v00) < 0d0) okay = .false.
3306 : end if
3307 :
3308 : if (k < nz) then ! check vp1 <= vdby2 <= v00 or vp1 >= vdby2 >= v00
3309 : if ((v00 - vbdy2)*(vbdy2 - vp1) < 0d0) okay = .false.
3310 : end if
3311 :
3312 : if (okay) return
3313 :
3314 : c0 = v00
3315 : c2 = 0d0
3316 :
3317 : if (k == 1) then
3318 : c1 = (v00 - vp1) / (0.5d0*(d00 + dp1))
3319 : else if (k == nz) then
3320 : c1 = (vm1 - v00) / (0.5d0*(dm1 + d00))
3321 : else if (abs(sm1) <= abs(sp1)) then
3322 : c1 = sm1
3323 : else
3324 : c1 = sp1
3325 : end if
3326 :
3327 : end subroutine get1_lpp
3328 :
3329 :
3330 0 : subroutine calc_Ptrb_ad_tw(s, k, Ptrb, Ptrb_div_etrb, ierr)
3331 : ! note: Ptrb_div_etrb is not time weighted
3332 : ! erg cm^-3 = g cm^2 s^-2 cm^-3 = g cm^-1 s^-2
3333 : use auto_diff
3334 : use auto_diff_support
3335 : type (star_info), pointer :: s
3336 : integer, intent(in) :: k
3337 : type(auto_diff_real_star_order1), intent(out) :: Ptrb, Ptrb_div_etrb
3338 : integer, intent(out) :: ierr
3339 : type(auto_diff_real_star_order1) :: etrb, rho
3340 : real(dp) :: Ptrb_start
3341 : real(dp), parameter :: x_ALFAP = 2.d0/3.d0
3342 : logical :: time_center, test_partials
3343 : include 'formats'
3344 0 : ierr = 0
3345 : if (s% RSP2_alfap == 0 .or. s% mixing_length_alpha == 0 .or. &
3346 0 : k <= s% RSP2_num_outermost_cells_forced_nonturbulent .or. &
3347 : k > s% nz - int(s% nz/s% RSP_nz_div_IBOTOM)) then
3348 0 : Ptrb_div_etrb = 0d0
3349 0 : Ptrb = 0d0
3350 0 : return
3351 : end if
3352 0 : rho = wrap_d_00(s,k)
3353 0 : etrb = wrap_etrb_00(s,k)
3354 0 : Ptrb_div_etrb = s% RSP2_alfap*x_ALFAP*etrb*rho
3355 0 : Ptrb = Ptrb_div_etrb*etrb ! cm^2 s^-2 g cm^-3 = erg cm^-3
3356 : time_center = (s% using_velocity_time_centering .and. &
3357 0 : s% include_P_in_velocity_time_centering)
3358 : if (time_center) then
3359 0 : Ptrb_start = s% RSP2_alfap*get_etrb_start(s,k)*s% rho_start(k)
3360 : Ptrb = s% P_theta_for_velocity_time_centering*Ptrb + &
3361 0 : (1d0 - s% P_theta_for_velocity_time_centering)*Ptrb_start
3362 : end if
3363 :
3364 0 : if (is_bad(Ptrb%val)) then
3365 0 : !$omp critical (calc_Ptrb_ad_tw_crit)
3366 0 : write(*,2) 'Ptrb', k, Ptrb%val
3367 0 : call mesa_error(__FILE__,__LINE__,'calc_Ptrb_tw')
3368 : !$omp end critical (calc_Ptrb_ad_tw_crit)
3369 : end if
3370 :
3371 : !test_partials = (k == s% solver_test_partials_k)
3372 0 : test_partials = .false.
3373 : if (test_partials) then
3374 : s% solver_test_partials_val = Ptrb%val
3375 : !s% solver_test_partials_var = i_var_R
3376 : !s% solver_test_partials_dval_dx = 0 ! d_residual_dr_00
3377 : write(*,*) 'calc_Ptrb_ad_tw', s% solver_test_partials_var
3378 : end if
3379 :
3380 : end subroutine calc_Ptrb_ad_tw
3381 :
3382 :
3383 : ! Ptot_ad = Peos_ad + Pvsc_ad + Ptrb_ad + mlt_Pturb_ad with time weighting
3384 104608 : subroutine calc_Ptot_ad_tw( &
3385 104608 : s, k, skip_Peos, skip_mlt_Pturb, Ptot_ad, d_Ptot_dxa, ierr)
3386 : use auto_diff_support
3387 : type (star_info), pointer :: s
3388 : integer, intent(in) :: k
3389 : logical, intent(in) :: skip_Peos, skip_mlt_Pturb
3390 : type(auto_diff_real_star_order1), intent(out) :: Ptot_ad
3391 : real(dp), dimension(s% species), intent(out) :: d_Ptot_dxa
3392 : integer, intent(out) :: ierr
3393 : integer :: j
3394 : real(dp) :: mlt_Pturb_start, alfa, beta
3395 : type(auto_diff_real_star_order1) :: &
3396 : Peos_ad, Pvsc_ad, Ptrb_ad, mlt_Pturb_ad, Ptrb_ad_div_etrb
3397 : logical :: time_center
3398 : include 'formats'
3399 :
3400 104608 : ierr = 0
3401 941472 : d_Ptot_dxa = 0d0
3402 :
3403 : time_center = (s% using_velocity_time_centering .and. &
3404 : s% include_P_in_velocity_time_centering .and. &
3405 104608 : s% lnT(k)/ln10 <= s% max_logT_for_include_P_and_L_in_velocity_time_centering)
3406 : if (time_center) then
3407 0 : alfa = s% P_theta_for_velocity_time_centering
3408 : else
3409 104608 : alfa = 1d0
3410 : end if
3411 104608 : beta = 1d0 - alfa
3412 :
3413 104608 : Peos_ad = 0d0
3414 104608 : if (.not. skip_Peos) then
3415 104608 : Peos_ad = wrap_peos_00(s, k)
3416 104608 : Peos_ad = alfa*Peos_ad + beta*s% Peos_start(k)
3417 941472 : do j=1,s% species
3418 836864 : d_Ptot_dxa(j) = s% Peos(k)*s% dlnPeos_dxa_for_partials(j,k)
3419 941472 : d_Ptot_dxa(j) = alfa*d_Ptot_dxa(j)
3420 : end do
3421 : end if
3422 :
3423 104608 : Pvsc_ad = 0d0
3424 104608 : if (s% use_Pvsc_art_visc) then
3425 0 : call get_Pvsc_ad(s, k, Pvsc_ad, ierr) ! no time centering for Pvsc
3426 0 : if (ierr /= 0) return
3427 : ! NO TIME CENTERING FOR Pvsc: Pvsc_ad = alfa*Pvsc_ad + beta*s% Pvsc_start(k)
3428 : end if
3429 :
3430 104608 : Ptrb_ad = 0d0
3431 104608 : if (s% RSP2_flag) then
3432 0 : call calc_Ptrb_ad_tw(s, k, Ptrb_ad, Ptrb_ad_div_etrb, ierr)
3433 0 : if (ierr /= 0) return
3434 : ! note that Ptrb_ad is already time weighted
3435 : end if
3436 :
3437 104608 : mlt_Pturb_ad = 0d0
3438 104608 : if ((.not. skip_mlt_Pturb) .and. s% mlt_Pturb_factor > 0d0 .and. k > 1) then
3439 0 : mlt_Pturb_ad = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*get_rho_face(s,k)/3d0
3440 0 : if (time_center) then
3441 : mlt_Pturb_start = &
3442 0 : s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*(s% rho_start(k-1) + s% rho_start(k))/6d0
3443 0 : mlt_Pturb_ad = alfa*mlt_Pturb_ad + beta*mlt_Pturb_start
3444 : end if
3445 : end if
3446 :
3447 104608 : Ptot_ad = Peos_ad + Pvsc_ad + Ptrb_ad + mlt_Pturb_ad
3448 :
3449 104608 : if (s% use_other_pressure) Ptot_ad = Ptot_ad + s% extra_pressure(k)
3450 :
3451 : end subroutine calc_Ptot_ad_tw
3452 :
3453 :
3454 0 : subroutine get_Pvsc_ad(s, k, Pvsc, ierr)
3455 : use auto_diff
3456 : use auto_diff_support
3457 : type (star_info), pointer :: s
3458 : integer, intent(in) :: k
3459 : type(auto_diff_real_star_order1), intent(out) :: Pvsc
3460 : integer, intent(out) :: ierr
3461 : type(auto_diff_real_star_order1) :: v00, vp1, Peos, rho, &
3462 : Peos_div_rho, dv
3463 : real(dp) :: Pvsc_start, cq, zsh
3464 0 : Pvsc = 0
3465 0 : s% Pvsc(k) = 0d0
3466 0 : Pvsc_start = s% Pvsc_start(k)
3467 0 : if (Pvsc_start < 0d0) s% Pvsc_start(k) = 0d0
3468 0 : if (.not. (s% v_flag .and. s% use_Pvsc_art_visc)) return
3469 0 : cq = s% Pvsc_cq
3470 0 : if (cq == 0d0) return
3471 0 : zsh = s% Pvsc_zsh
3472 0 : v00 = wrap_v_00(s,k)
3473 0 : vp1 = wrap_v_p1(s,k)
3474 0 : Peos = wrap_Peos_00(s,k)
3475 0 : rho = wrap_d_00(s,k)
3476 0 : Peos_div_rho = Peos/rho
3477 0 : dv = (vp1 - v00) - zsh*sqrt(Peos_div_rho)
3478 0 : if (dv%val <= 0d0) return
3479 0 : Pvsc = cq*rho*pow2(dv)
3480 0 : s% Pvsc(k) = Pvsc%val
3481 0 : if (Pvsc_start < 0d0) s% Pvsc_start(k) = s% Pvsc(k)
3482 : end subroutine get_Pvsc_ad
3483 :
3484 :
3485 : ! marsaglia and zaman random number generator. period is 2**43 with
3486 : ! 900 million different sequences. the state of the generator (for restarts)
3487 2 : subroutine init_random(s)
3488 : type (star_info), pointer :: s
3489 : integer :: ijkl,ij,kl,i,j,k,l,ii,jj,m
3490 : real(dp) :: x,t
3491 2 : ijkl = 54217137
3492 2 : ij = ijkl/30082
3493 2 : kl = ijkl - 30082 * ij
3494 2 : i = mod(ij/177,177) + 2
3495 2 : j = mod(ij,177) + 2
3496 2 : k = mod(kl/169,178) + 1
3497 2 : l = mod(kl,169)
3498 196 : do ii=1,97
3499 : x = 0.0d0
3500 : t = 0.5d0
3501 4850 : do jj =1,24
3502 4656 : m = mod(mod(i*j,179)*k,179)
3503 4656 : i = j
3504 4656 : j = k
3505 4656 : k = m
3506 4656 : l = mod(53*l + 1, 169)
3507 4656 : if (mod(l*m,64) >= 32) x = x + t
3508 4850 : t = 0.5d0 * t
3509 : end do
3510 196 : s% rand_u(ii) = x
3511 : end do
3512 2 : s% rand_c = 362436.0d0/16777216.0d0
3513 2 : s% rand_cd = 7654321.0d0/16777216.0d0
3514 2 : s% rand_cm = 16777213.0d0/16777216.0d0
3515 2 : s% rand_i97 = 97
3516 2 : s% rand_j97 = 33
3517 2 : end subroutine init_random
3518 :
3519 :
3520 0 : real(dp) function rand(s)
3521 : type (star_info), pointer :: s
3522 : real(dp) :: uni
3523 0 : uni = s% rand_u(s% rand_i97) - s% rand_u(s% rand_j97)
3524 0 : if (uni < 0.0d0) uni = uni + 1.0d0
3525 0 : s% rand_u(s% rand_i97) = uni
3526 0 : s% rand_i97 = s% rand_i97 - 1
3527 0 : if (s% rand_i97 == 0) s% rand_i97 = 97
3528 0 : s% rand_j97 = s% rand_j97 - 1
3529 0 : if (s% rand_j97 == 0) s% rand_j97 = 97
3530 0 : s% rand_c = s% rand_c - s% rand_cd
3531 0 : if (s% rand_c < 0.0d0) s% rand_c = s% rand_c + s% rand_cm
3532 0 : uni = uni - s% rand_c
3533 0 : if (uni < 0.0d0) uni = uni + 1.0d0
3534 0 : rand = uni
3535 0 : end function rand
3536 :
3537 :
3538 4 : subroutine write_to_extra_terminal_output_file(s, str, advance)
3539 : type (star_info), pointer :: s
3540 : character(len=*), intent(in) :: str
3541 : logical, intent(in) :: advance
3542 : integer :: ierr, io
3543 4 : if (len_trim(s% extra_terminal_output_file) > 0) then
3544 0 : ierr = 0
3545 : open(newunit=io, file=trim(s% extra_terminal_output_file), &
3546 0 : action='write', position='append', iostat=ierr)
3547 0 : if (ierr == 0) then
3548 0 : if (advance) then
3549 0 : write(io,'(a)') trim(str)
3550 : else
3551 0 : write(io,'(a)',advance='no') trim(str)
3552 : end if
3553 0 : close(io)
3554 : else
3555 : write(*,*) 'failed to open extra_terminal_output_file ' // &
3556 0 : trim(s% extra_terminal_output_file)
3557 : end if
3558 : end if
3559 4 : end subroutine write_to_extra_terminal_output_file
3560 :
3561 :
3562 0 : subroutine write_eos_call_info(s,k)
3563 : use chem_def
3564 : type (star_info), pointer :: s
3565 : integer, intent(in) :: k ! 0 means not being called for a particular cell
3566 : integer :: j
3567 : include 'formats'
3568 0 : !$OMP critical (omp_write_eos_call_info)
3569 0 : write(*,'(A)')
3570 0 : do j=1,s% species
3571 0 : write(*,4) 'xa(j,k) ' // trim(chem_isos% name(s% chem_id(j))), j, j+s% nvar_hydro, k, s% xa(j,k)
3572 : end do
3573 0 : write(*,1) 'sum(xa) =', sum(s% xa(:,k))
3574 0 : write(*,1) 'Pgas = ', s% Pgas(k)
3575 0 : write(*,1) 'logPgas = ', s% lnPgas(k)/ln10
3576 0 : write(*,1) 'rho = ', s% rho(k)
3577 0 : write(*,1) 'T = ', s% T(k)
3578 0 : write(*,1) 'logQ = ', s% lnd(k)/ln10 - 2*s% lnT(k)/ln10 + 12
3579 0 : write(*,'(A)')
3580 0 : write(*,1) 'eos_frac_OPAL_SCVH', s% eos_frac_OPAL_SCVH(k)
3581 0 : write(*,1) 'eos_frac_HELM', s% eos_frac_HELM(k)
3582 0 : write(*,1) 'eos_frac_Skye', s% eos_frac_Skye(k)
3583 0 : write(*,1) 'eos_frac_PC', s% eos_frac_PC(k)
3584 0 : write(*,1) 'eos_frac_FreeEOS', s% eos_frac_FreeEOS(k)
3585 0 : write(*,1) 'eos_frac_CMS', s% eos_frac_CMS(k)
3586 0 : write(*,1) 'eos_frac_ideal', s% eos_frac_ideal(k)
3587 0 : write(*,'(A)')
3588 0 : write(*,1) 'Peos = ', s% Peos(k)
3589 0 : write(*,1) 'Prad = ', s% Prad(k)
3590 0 : write(*,1) 'logPeos = ', s% lnPeos(k)/ln10
3591 0 : write(*,1) 'logS = ', s% lnS(k)/ln10
3592 0 : write(*,1) 'logE = ', s% lnE(k)/ln10
3593 0 : write(*,1) 'energy = ', s% energy(k)
3594 0 : write(*,1) 'grada = ', s% grada(k)
3595 0 : write(*,1) 'Cv = ', s% Cv(k)
3596 0 : write(*,1) 'dE_dRho = ', s% dE_dRho(k)
3597 0 : write(*,1) 'cp = ', s% cp(k)
3598 0 : write(*,1) 'gamma1 = ', s% gamma1(k)
3599 0 : write(*,1) 'gamma3 = ', s% gamma3(k)
3600 0 : write(*,1) 'eta = ', s% eta(k)
3601 0 : write(*,1) 'gam = ', s% gam(k)
3602 0 : write(*,1) 'mu = ', s% mu(k)
3603 0 : write(*,1) 'log_free_e = ', s% lnfree_e(k)/ln10
3604 0 : write(*,1) 'chiRho = ', s% chiRho(k)
3605 0 : write(*,1) 'chiT = ', s% chiT(k)
3606 0 : write(*,'(A)')
3607 0 : write(*,*) 'do_eos_for_cell k, nz', k, s% nz
3608 0 : write(*,1) 'logRho = ', s% lnd(k)/ln10
3609 0 : write(*,1) 'logT = ', s% lnT(k)/ln10
3610 0 : write(*,1) 'z = ', s% Z(k)
3611 0 : write(*,1) 'x = ', s% X(k)
3612 0 : write(*,1) 'abar = ', s% abar(k)
3613 0 : write(*,1) 'zbar = ', s% zbar(k)
3614 0 : write(*,'(A)')
3615 0 : write(*,1) 'tau = ', s% tau(k)
3616 0 : write(*,'(A)')
3617 0 : write(*,*) 's% eos_rq% tiny_fuzz', s% eos_rq% tiny_fuzz
3618 0 : write(*,'(A)')
3619 : !call mesa_error(__FILE__,__LINE__,'write_eos_call_info')
3620 : !$OMP end critical (omp_write_eos_call_info)
3621 0 : end subroutine write_eos_call_info
3622 :
3623 :
3624 239 : real(dp) function surface_avg_x(s,j)
3625 : type (star_info), pointer :: s
3626 : integer, intent(in) :: j
3627 : real(dp) :: sum_x, sum_dq
3628 : integer :: k
3629 : include 'formats'
3630 239 : if (j == 0) then
3631 239 : surface_avg_x = 0d0
3632 : return
3633 : end if
3634 239 : sum_x = 0
3635 239 : sum_dq = 0
3636 50668 : do k = 1, s% nz
3637 50668 : sum_x = sum_x + s% xa(j,k)*s% dq(k)
3638 50668 : sum_dq = sum_dq + s% dq(k)
3639 50668 : if (sum_dq >= s% surface_avg_abundance_dq) exit
3640 : end do
3641 239 : surface_avg_x = sum_x/sum_dq
3642 239 : end function surface_avg_x
3643 :
3644 :
3645 300 : real(dp) function center_avg_x(s,j)
3646 : type (star_info), pointer :: s
3647 : integer, intent(in) :: j
3648 : real(dp) :: sum_x, sum_dq, dx, dq
3649 : integer :: k
3650 300 : if (j == 0) then
3651 300 : center_avg_x = 0d0
3652 : return
3653 : end if
3654 267 : sum_x = 0
3655 267 : sum_dq = 0
3656 267 : do k = s% nz, 1, -1
3657 267 : dq = s% dq(k)
3658 267 : dx = s% xa(j,k)*dq
3659 267 : if (sum_dq+dq >= s% center_avg_value_dq) then
3660 267 : sum_x = sum_x+ dx*(s% center_avg_value_dq - sum_dq)/dq
3661 267 : sum_dq = s% center_avg_value_dq
3662 267 : exit
3663 : end if
3664 0 : sum_x = sum_x + dx
3665 0 : sum_dq = sum_dq + dq
3666 : end do
3667 267 : center_avg_x = sum_x/sum_dq
3668 267 : end function center_avg_x
3669 :
3670 :
3671 209260 : subroutine get_area_info_opt_time_center(s, k, area, inv_R2, ierr)
3672 : use auto_diff_support
3673 : type (star_info), pointer :: s
3674 : integer, intent(in) :: k
3675 : type(auto_diff_real_star_order1), intent(out) :: area, inv_R2
3676 : integer, intent(out) :: ierr
3677 : type(auto_diff_real_star_order1) :: r_00, r2_00
3678 209260 : ierr = 0
3679 209260 : r_00 = wrap_r_00(s,k)
3680 209260 : r2_00 = pow2(r_00)
3681 209260 : if (s% using_velocity_time_centering) then
3682 0 : area = 4d0*pi*(r2_00 + r_00*s% r_start(k) + s% r_start(k)**2)/3d0
3683 0 : inv_R2 = 1d0/(r_00*s% r_start(k))
3684 : else
3685 209260 : area = 4d0*pi*r2_00
3686 209260 : inv_R2 = 1d0/r2_00
3687 : end if
3688 209260 : end subroutine get_area_info_opt_time_center
3689 :
3690 :
3691 52348 : subroutine set_energy_eqn_scal(s, k, scal, ierr) ! 1/(erg g^-1 s^-1)
3692 : type (star_info), pointer :: s
3693 : integer, intent(in) :: k
3694 : real(dp), intent(out) :: scal
3695 : integer, intent(out) :: ierr
3696 : real(dp) :: cell_energy_fraction_start
3697 : include 'formats'
3698 52348 : ierr = 0
3699 52348 : if (k > 1) then
3700 52304 : scal = 1d0
3701 : else
3702 44 : scal = 1d-6
3703 : end if
3704 52348 : if (s% dedt_eqn_r_scale > 0d0) then
3705 : cell_energy_fraction_start = &
3706 52348 : s% energy_start(k)*s% dm(k)/s% total_internal_energy_old
3707 52348 : scal = min(scal, cell_energy_fraction_start*s% dedt_eqn_r_scale)
3708 : end if
3709 52348 : scal = scal*s% dt/s% energy_start(k)
3710 52348 : end subroutine set_energy_eqn_scal
3711 :
3712 :
3713 107640 : real(dp) function conv_time_scale(s,k_in) result(tau_conv)
3714 : type (star_info), pointer :: s
3715 : integer, intent(in) :: k_in
3716 : integer :: k
3717 : real(dp) :: brunt_B, alfa, beta, rho_face, Peos_face, chiT_face, chiRho_face, &
3718 : f, dlnP, dlnT, grada_face, gradT_actual, brunt_N2
3719 53820 : if (.not. s% calculate_Brunt_B) then
3720 19764 : tau_conv = 0d0
3721 : return
3722 : end if
3723 53820 : k = max(2,k_in)
3724 53820 : brunt_B = s% brunt_B(k)
3725 53820 : call get_face_weights(s, k, alfa, beta)
3726 53820 : rho_face = alfa*s% rho(k) + beta*s% rho(k-1)
3727 53820 : Peos_face = alfa*s% Peos(k) + beta*s% Peos(k-1)
3728 53820 : chiT_face = alfa*s% chiT(k) + beta*s% chiT(k-1)
3729 53820 : chiRho_face = alfa*s% chiRho(k) + beta*s% chiRho(k-1)
3730 53820 : f = pow2(s% grav(k))*rho_face/Peos_face*chiT_face/chiRho_face
3731 53820 : dlnP = s% lnPeos(k-1) - s% lnPeos(k)
3732 53820 : dlnT = s% lnT(k-1) - s% lnT(k)
3733 53820 : grada_face = alfa*s% grada(k) + beta*s% grada(k-1)
3734 53820 : gradT_actual = safe_div_val(s, dlnT, dlnP) ! mlt has not been called yet when doing this
3735 53820 : brunt_N2 = f*(brunt_B - (gradT_actual - grada_face))
3736 53820 : if(abs(brunt_B) > 0d0) then
3737 34056 : tau_conv = 1d0/sqrt(abs(brunt_N2))
3738 : else
3739 : tau_conv = 0d0
3740 : end if
3741 : end function conv_time_scale
3742 :
3743 :
3744 44 : subroutine set_conv_time_scales(s)
3745 : type (star_info), pointer :: s
3746 : integer :: k
3747 : real(dp) :: tau_conv
3748 : include 'formats'
3749 44 : s% min_conv_time_scale = 1d99
3750 44 : s% max_conv_time_scale = 0d0
3751 53864 : do k=1,s%nz
3752 53820 : if (s% X(k) > s% max_X_for_conv_timescale) cycle
3753 53820 : if (s% X(k) < s% min_X_for_conv_timescale) cycle
3754 53820 : if (s% q(k) > s% max_q_for_conv_timescale) cycle
3755 53820 : if (s% q(k) < s% min_q_for_conv_timescale) exit
3756 53820 : tau_conv = conv_time_scale(s,k)
3757 53820 : if (tau_conv < s% min_conv_time_scale) &
3758 44 : s% min_conv_time_scale = tau_conv
3759 53820 : if (tau_conv > s% max_conv_time_scale) &
3760 2783 : s% max_conv_time_scale = tau_conv
3761 : end do
3762 44 : if (s% max_conv_time_scale == 0d0) s% max_conv_time_scale = 1d99
3763 44 : if (s% min_conv_time_scale == 1d99) s% min_conv_time_scale = 0d0
3764 44 : end subroutine set_conv_time_scales
3765 :
3766 :
3767 0 : real(dp) function QHSE_time_scale(s,k) result(tau_qhse)
3768 : type (star_info), pointer :: s
3769 : integer, intent(in) :: k
3770 : real(dp) :: abs_dv
3771 0 : if (s% v_flag) then
3772 0 : abs_dv = abs(s% v(k) - s% v_start(k))
3773 0 : else if (s% u_flag) then
3774 0 : abs_dv = abs(s% u_face_ad(k)%val - s% u_face_start(k))
3775 : else
3776 : abs_dv = 0d0
3777 : end if
3778 0 : tau_qhse = abs_dv/(s% cgrav(k)*s% m_grav(k)/pow2(s% r(k)))
3779 0 : end function QHSE_time_scale
3780 :
3781 :
3782 0 : real(dp) function eps_nuc_time_scale(s,k) result(tau_epsnuc)
3783 : type (star_info), pointer :: s
3784 : integer, intent(in) :: k
3785 0 : tau_epsnuc = s% Cp(k)*s% T(k)/max(1d-10,abs(s% eps_nuc(k)))
3786 0 : end function eps_nuc_time_scale
3787 :
3788 :
3789 0 : real(dp) function cooling_time_scale(s,k) result(tau_cool)
3790 : type (star_info), pointer :: s
3791 : integer, intent(in) :: k
3792 : real(dp) :: thermal_conductivity
3793 0 : thermal_conductivity = (4d0*crad*clight*pow3(s% T(k)))/(3d0*s% opacity(k)*s% rho(k)*s% Cp(k))
3794 0 : tau_cool = pow2(s% scale_height(k)) / thermal_conductivity
3795 0 : end function cooling_time_scale
3796 :
3797 :
3798 319844 : function get_rho_face(s,k) result(rho_face)
3799 : type (star_info), pointer :: s
3800 : integer, intent(in) :: k
3801 : type(auto_diff_real_star_order1) :: rho_face
3802 : real(dp) :: alfa, beta
3803 159988 : if (k == 1) then
3804 132 : rho_face = wrap_d_00(s,k)
3805 132 : return
3806 : end if
3807 159856 : call get_face_weights(s, k, alfa, beta)
3808 159856 : rho_face = alfa*wrap_d_00(s,k) + beta*wrap_d_m1(s,k)
3809 159856 : end function get_rho_face
3810 :
3811 :
3812 0 : real(dp) function get_rho_face_val(s,k) result(rho_face)
3813 : type (star_info), pointer :: s
3814 : integer, intent(in) :: k
3815 : real(dp) :: alfa, beta
3816 0 : if (k == 1) then
3817 0 : rho_face = s% rho(1)
3818 0 : return
3819 : end if
3820 0 : call get_face_weights(s, k, alfa, beta)
3821 0 : rho_face = alfa*s% rho(k) + beta*s% rho(k-1)
3822 0 : end function get_rho_face_val
3823 :
3824 159922 : function get_rho_start_face(s,k) result(rho_face)
3825 : type (star_info), pointer :: s
3826 : integer, intent(in) :: k
3827 : type(auto_diff_real_star_order1) :: rho_face
3828 : real(dp) :: alfa, beta
3829 79994 : if (k == 1) then
3830 66 : rho_face = wrap_d_00_start(s,k)
3831 66 : return
3832 : end if
3833 79928 : call get_face_weights(s, k, alfa, beta)
3834 79928 : rho_face = alfa*wrap_d_00_start(s,k) + beta*wrap_d_m1_start(s,k)
3835 79928 : end function get_rho_start_face
3836 :
3837 319844 : function get_T_face(s,k) result(T_face)
3838 : type (star_info), pointer :: s
3839 : integer, intent(in) :: k
3840 : type(auto_diff_real_star_order1) :: T_face
3841 : real(dp) :: alfa, beta
3842 159988 : if (k == 1) then
3843 132 : T_face = wrap_T_00(s,k)
3844 132 : return
3845 : end if
3846 159856 : call get_face_weights(s, k, alfa, beta)
3847 159856 : T_face = alfa*wrap_T_00(s,k) + beta*wrap_T_m1(s,k)
3848 159856 : end function get_T_face
3849 :
3850 :
3851 79994 : function get_Prad_face(s,k) result(Prad_face)
3852 : type (star_info), pointer :: s
3853 : integer, intent(in) :: k
3854 : type(auto_diff_real_star_order1) :: Prad_face
3855 79994 : Prad_face = crad*pow4(get_T_face(s,k))/3d0
3856 79994 : end function get_Prad_face
3857 :
3858 :
3859 639688 : function get_Peos_face(s,k) result(Peos_face)
3860 : type (star_info), pointer :: s
3861 : integer, intent(in) :: k
3862 : type(auto_diff_real_star_order1) :: Peos_face
3863 : real(dp) :: alfa, beta
3864 319976 : if (k == 1) then
3865 264 : Peos_face = wrap_Peos_00(s,k)
3866 264 : return
3867 : end if
3868 319712 : call get_face_weights(s, k, alfa, beta)
3869 319712 : Peos_face = alfa*wrap_Peos_00(s,k) + beta*wrap_Peos_m1(s,k)
3870 319712 : end function get_Peos_face
3871 :
3872 26163 : function get_Peos_face_val(s,k) result(Peos_face)
3873 : type (star_info), pointer :: s
3874 : integer, intent(in) :: k
3875 : real(dp) :: Peos_face
3876 : real(dp) :: alfa, beta
3877 13087 : if (k == 1) then
3878 11 : Peos_face = s% Peos(k)
3879 11 : return
3880 : end if
3881 13076 : call get_face_weights(s, k, alfa, beta)
3882 13076 : Peos_face = alfa*s% Peos(k) + beta*s% Peos(k-1)
3883 13076 : end function get_Peos_face_val
3884 :
3885 159922 : function get_e_face(s,k) result(e_face) ! specific energy on face
3886 : type (star_info), pointer :: s
3887 : integer, intent(in) :: k
3888 : type(auto_diff_real_star_order1) :: e_face
3889 : real(dp) :: alfa, beta
3890 79994 : if (k == 1) then
3891 66 : e_face = wrap_e_00(s,k)
3892 66 : return
3893 : end if
3894 79928 : call get_face_weights(s, k, alfa, beta)
3895 79928 : e_face = alfa*wrap_e_00(s,k) + beta*wrap_e_m1(s,k)
3896 79928 : end function get_e_face
3897 :
3898 159922 : function get_Cp_face(s,k) result(Cp_face)
3899 : type (star_info), pointer :: s
3900 : integer, intent(in) :: k
3901 : type(auto_diff_real_star_order1) :: Cp_face
3902 : real(dp) :: alfa, beta
3903 79994 : if (k == 1) then
3904 66 : Cp_face = wrap_Cp_00(s,k)
3905 66 : return
3906 : end if
3907 79928 : call get_face_weights(s, k, alfa, beta)
3908 79928 : Cp_face = alfa*wrap_Cp_00(s,k) + beta*wrap_Cp_m1(s,k)
3909 79928 : end function get_Cp_face
3910 :
3911 :
3912 159922 : function get_ChiRho_face(s,k) result(ChiRho_face)
3913 : type (star_info), pointer :: s
3914 : integer, intent(in) :: k
3915 : type(auto_diff_real_star_order1) :: ChiRho_face
3916 : real(dp) :: alfa, beta
3917 79994 : if (k == 1) then
3918 66 : ChiRho_face = wrap_ChiRho_00(s,k)
3919 66 : return
3920 : end if
3921 79928 : call get_face_weights(s, k, alfa, beta)
3922 79928 : ChiRho_face = alfa*wrap_ChiRho_00(s,k) + beta*wrap_ChiRho_m1(s,k)
3923 79928 : end function get_ChiRho_face
3924 :
3925 :
3926 159922 : function get_ChiT_face(s,k) result(ChiT_face)
3927 : type (star_info), pointer :: s
3928 : integer, intent(in) :: k
3929 : type(auto_diff_real_star_order1) :: ChiT_face
3930 : real(dp) :: alfa, beta
3931 79994 : if (k == 1) then
3932 66 : ChiT_face = wrap_ChiT_00(s,k)
3933 66 : return
3934 : end if
3935 79928 : call get_face_weights(s, k, alfa, beta)
3936 79928 : ChiT_face = alfa*wrap_ChiT_00(s,k) + beta*wrap_ChiT_m1(s,k)
3937 79928 : end function get_ChiT_face
3938 :
3939 :
3940 319844 : function get_kap_face(s,k) result(kap_face)
3941 : type (star_info), pointer :: s
3942 : integer, intent(in) :: k
3943 : type(auto_diff_real_star_order1) :: kap_face
3944 : real(dp) :: alfa, beta
3945 159988 : if (k == 1) then
3946 132 : kap_face = wrap_kap_00(s,k)
3947 132 : return
3948 : end if
3949 159856 : call get_face_weights(s, k, alfa, beta)
3950 159856 : kap_face = alfa*wrap_kap_00(s,k) + beta*wrap_kap_m1(s,k)
3951 159856 : end function get_kap_face
3952 :
3953 :
3954 159922 : function get_grada_face(s,k) result(grada_face)
3955 : type (star_info), pointer :: s
3956 : integer, intent(in) :: k
3957 : type(auto_diff_real_star_order1) :: grada_face, mlt_Pturb_ad, P, gamma1_face, alpha
3958 : real(dp) :: alfa, beta
3959 79994 : P = get_Peos_face(s,k)
3960 79994 : if (k == 1) then
3961 66 : grada_face = wrap_grad_ad_00(s,k)
3962 66 : return
3963 : end if
3964 79928 : call get_face_weights(s, k, alfa, beta)
3965 79928 : grada_face = alfa*wrap_grad_ad_00(s,k) + beta*wrap_grad_ad_m1(s,k)
3966 : if (s% have_mlt_vc .and. s% okay_to_set_mlt_vc .and. s% include_mlt_Pturb_in_thermodynamic_gradients &
3967 79928 : .and. s% mlt_Pturb_factor > 0d0 .and. k > 1) then
3968 0 : gamma1_face = alfa*wrap_gamma1_00(s,k) + beta*wrap_gamma1_m1(s,k)
3969 0 : mlt_Pturb_ad = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*get_rho_face(s,k)/3d0
3970 0 : alpha = mlt_Pturb_ad/ (P*gamma1_face)
3971 0 : grada_face = grada_face*(P+mlt_Pturb_ad) / (P * (1d0 + alpha))
3972 : end if
3973 79928 : end function get_grada_face
3974 :
3975 :
3976 79994 : function get_gradr_face(s,k) result(gradr)
3977 : type (star_info), pointer :: s
3978 : integer, intent(in) :: k
3979 : type(auto_diff_real_star_order1) :: gradr
3980 : type(auto_diff_real_star_order1) :: P, opacity, L, Pr, r
3981 : real(dp) :: L_theta
3982 : !include 'formats'
3983 79994 : P = get_Peos_face(s,k)
3984 79994 : opacity = get_kap_face(s,k)
3985 :
3986 79994 : if (s% include_mlt_in_velocity_time_centering) then
3987 : ! consider building a wrapper : wrap_opt_time_center_L_00(s,k)
3988 : if (s% using_velocity_time_centering .and. &
3989 0 : s% include_L_in_velocity_time_centering .and. &
3990 : s% lnT(k)/ln10 <= s% max_logT_for_include_P_and_L_in_velocity_time_centering) then
3991 0 : L_theta = s% L_theta_for_velocity_time_centering
3992 : else
3993 0 : L_theta = 1d0
3994 : end if
3995 0 : L = L_theta*wrap_L_00(s, k) + (1d0 - L_theta)*s% L_start(k)
3996 0 : r = wrap_opt_time_center_r_00(s,k)
3997 : else
3998 79994 : L = wrap_L_00(s,k)
3999 79994 : r = wrap_r_00(s,k)
4000 : end if
4001 79994 : Pr = get_Prad_face(s,k)
4002 79994 : gradr = P*opacity*L/(16d0*pi*clight*s% m_grav(k)*s% cgrav(k)*Pr)
4003 79994 : end function get_gradr_face
4004 :
4005 :
4006 79994 : function get_scale_height_face(s,k) result(scale_height)
4007 : type (star_info), pointer :: s
4008 : integer, intent(in) :: k
4009 : type(auto_diff_real_star_order1) :: scale_height
4010 : type(auto_diff_real_star_order1) :: grav, scale_height2, P, rho
4011 : real(dp) :: G, alfa, beta
4012 : include 'formats'
4013 79994 : G = s% cgrav(k)
4014 79994 : grav = G*s% m_grav(k)/pow2(wrap_r_00(s,k)) ! try geff later.
4015 79994 : if (s% use_rsp_form_of_scale_height .and. k >1) then ! use rsp form of Hp, assumes HSE, wraps P/rho together.
4016 0 : call get_face_weights(s, k, alfa, beta)
4017 0 : scale_height = (alfa*(wrap_Peos_00(s,k))/wrap_d_00(s,k) + beta*(wrap_Peos_m1(s,k))/wrap_d_m1(s,k))/grav
4018 : else
4019 79994 : P = get_Peos_face(s,k)
4020 79994 : rho = get_rho_face(s,k)
4021 79994 : scale_height = P/(grav*rho) ! this assumes HSE
4022 : end if
4023 79994 : if (s% alt_scale_height_flag) then
4024 : ! consider sound speed*hydro time scale as an alternative scale height
4025 : ! (this comes from Eggleton's code.)
4026 79994 : scale_height2 = sqrt(P/G)/rho
4027 79994 : if (scale_height2 < scale_height) then
4028 1905 : scale_height = scale_height2
4029 : end if
4030 : end if
4031 79994 : end function get_scale_height_face
4032 :
4033 :
4034 0 : real(dp) function get_scale_height_face_val(s,k) result(scale_height)
4035 : type (star_info), pointer :: s
4036 : integer, intent(in) :: k
4037 : real(dp) :: G, scale_height2, P, rho, alfa, beta
4038 : type(auto_diff_real_star_order1) :: P_face, rho_face, grav
4039 0 : G = s% cgrav(k)
4040 0 : grav = G*s% m_grav(k)/pow2(wrap_r_00(s,k))! try geff later
4041 0 : if (s% use_rsp_form_of_scale_height .and. k >1) then ! use rsp form of Hp, assumes HSE, wraps P/rho together.
4042 0 : call get_face_weights(s, k, alfa, beta)
4043 0 : scale_height = (alfa*(s% Peos(k)/s% rho(k)) + beta*(s% Peos(k-1)/s% rho(k-1)))/grav%val
4044 : else
4045 0 : P_face = get_Peos_face(s,k)
4046 0 : P = P_face%val
4047 0 : rho_face = get_rho_face(s,k)
4048 0 : rho = rho_face%val
4049 0 : scale_height = P/(grav%val*rho) ! this assumes HSE
4050 : end if
4051 0 : if (s% alt_scale_height_flag) then
4052 : ! consider sound speed*hydro time scale as an alternative scale height
4053 : ! (this comes from Eggleton's code.)
4054 0 : scale_height2 = sqrt(P/G)/rho
4055 : if (scale_height2 < scale_height) then
4056 0 : scale_height = scale_height2
4057 : end if
4058 : end if
4059 0 : end function get_scale_height_face_val
4060 :
4061 :
4062 : function get_QQ_cell(s,k) result(QQ_cell)
4063 : type (star_info), pointer :: s
4064 : integer, intent(in) :: k
4065 : type(auto_diff_real_star_order1) :: QQ_cell
4066 : type(auto_diff_real_star_order1) :: &
4067 : T_00, d_00, chiT_00, chiRho_00
4068 : T_00 = wrap_T_00(s,k)
4069 : d_00 = wrap_d_00(s,k)
4070 : chiT_00 = wrap_chiT_00(s,k)
4071 : chiRho_00 = wrap_chiRho_00(s,k)
4072 : QQ_cell = chiT_00/(d_00*T_00*chiRho_00)
4073 : end function get_QQ_cell
4074 :
4075 :
4076 1345744 : subroutine get_face_weights(s, k, alfa, beta)
4077 : type (star_info), pointer :: s
4078 : integer, intent(in) :: k
4079 : real(dp), intent(out) :: alfa, beta
4080 : ! face_value(k) = alfa*cell_value(k) + beta*cell_value(k-1)
4081 1345744 : if (k == 1) call mesa_error(__FILE__,__LINE__,'bad k==1 for get_face_weights')
4082 1345744 : alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
4083 1345744 : beta = 1d0 - alfa
4084 1345744 : end subroutine get_face_weights
4085 :
4086 :
4087 161368 : real(dp) function safe_div_val(s, x, y, lim) result(x_div_y)
4088 : type (star_info), pointer :: s
4089 : real(dp), intent(in) :: x, y, lim
4090 : optional :: lim
4091 : real(dp) :: limit
4092 161368 : if (present(lim)) then
4093 0 : limit = lim
4094 : else
4095 : limit = 1d-20
4096 : end if
4097 161368 : if (abs(y) < limit) then
4098 : x_div_y = 0d0
4099 : else
4100 161368 : x_div_y = x/y
4101 : end if
4102 161368 : end function safe_div_val
4103 :
4104 :
4105 : function safe_div(s, x, y, lim) result(x_div_y)
4106 : type (star_info), pointer :: s
4107 : type(auto_diff_real_star_order1), intent(in) :: x, y
4108 : type(auto_diff_real_star_order1) :: x_div_y
4109 : real(dp), intent(in) :: lim
4110 : optional :: lim
4111 : real(dp) :: limit
4112 : if (present(lim)) then
4113 : limit = lim
4114 : else
4115 : limit = 1d-20
4116 : end if
4117 : if (abs(y) < limit) then
4118 : x_div_y = 0d0
4119 : else
4120 : x_div_y = x/y
4121 : end if
4122 : end function safe_div
4123 :
4124 :
4125 43 : subroutine set_luminosity_by_category(s) ! integral by mass from center out
4126 : use chem_def, only: category_name
4127 : use rates_def, only: i_rate
4128 : use utils_lib, only: is_bad
4129 : type (star_info), pointer :: s
4130 : integer :: k, j
4131 : real(dp) :: L_burn_by_category(num_categories)
4132 : include 'formats'
4133 43 : L_burn_by_category(:) = 0
4134 51399 : do k = s% nz, 1, -1
4135 1283943 : do j = 1, num_categories
4136 : L_burn_by_category(j) = &
4137 1232544 : L_burn_by_category(j) + s% dm(k)*s% eps_nuc_categories(j, k)
4138 1232544 : if (is_bad(L_burn_by_category(j))) then
4139 0 : write(*,2) trim(category_name(j)) // ' eps_nuc logT', k, s% eps_nuc_categories(j,k), s% lnT(k)/ln10
4140 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'set_luminosity_by_category')
4141 : end if
4142 1283900 : s% luminosity_by_category(j,k) = L_burn_by_category(j)
4143 : end do
4144 : end do
4145 43 : end subroutine set_luminosity_by_category
4146 :
4147 :
4148 0 : subroutine set_zero_alpha_RTI(id, ierr)
4149 : integer, intent(in) :: id
4150 : integer, intent(out) :: ierr
4151 : type (star_info), pointer :: s
4152 : include 'formats'
4153 : ierr = 0
4154 0 : call get_star_ptr(id, s, ierr)
4155 0 : if (ierr /= 0) return
4156 0 : if (.not. s% u_flag) return
4157 0 : s% xh(s% i_alpha_RTI,1:s% nz) = 0d0
4158 0 : s% alpha_RTI(1:s% nz) = 0d0
4159 0 : s% need_to_setvars = .true.
4160 : end subroutine set_zero_alpha_RTI
4161 :
4162 : end module star_utils
|
