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