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