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