Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2011-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 mix_info
21 :
22 : use const_def, only: dp, i8, ln10, pi4, msun, no_mixing, &
23 : minimum_mixing, &
24 : rayleigh_taylor_mixing, &
25 : convective_mixing, &
26 : semiconvective_mixing, &
27 : overshoot_mixing, &
28 : thermohaline_mixing, &
29 : rotation_mixing
30 : use num_lib
31 : use utils_lib
32 : use star_private_def
33 :
34 : implicit none
35 :
36 : private
37 : public :: set_mixing_info
38 : public :: set_RTI_mixing_info
39 : public :: do_smoothing_by_mass
40 : public :: update_rotation_mixing_info
41 : public :: set_dPdr_dRhodr_info
42 : public :: get_convection_sigmas
43 : public :: set_dxdt_mix
44 : public :: set_cz_bdy_mass
45 :
46 : contains
47 :
48 33 : subroutine set_mixing_info(s, skip_set_cz_bdy_mass, ierr)
49 : ! set convection variables cdc and conv_vel starting from local MLT results.
50 : ! overshooting can also be added. and rotation mixing.
51 : use rates_def, only: i_rate
52 : use chem_def, only: ipp, icno, i3alf, ih1, ihe4, ic12
53 : use star_utils, only: start_time, update_time
54 : use overshoot, only: add_overshooting
55 : use predictive_mix, only: add_predictive_mixing
56 : use auto_diff_support, only: get_RSP2_conv_velocity
57 : type (star_info), pointer :: s
58 : logical, intent(in) :: skip_set_cz_bdy_mass
59 : integer, intent(out) :: ierr
60 :
61 : integer :: nz, k, max_conv_bdy, max_mix_bdy, k_Tmax, i_h1, i_he4, i_c12
62 : real(dp) :: rho_face, f, Tmax, min_conv_vel_for_convective_mixing_type, &
63 : region_bottom_q, region_top_q, L_val
64 33 : real(dp), allocatable, dimension(:) :: eps_h, eps_he, eps_z, cdc_factor
65 :
66 : logical :: RSP2_or_RSP
67 :
68 : integer(i8) :: time0
69 : real(dp) :: total
70 :
71 : include 'formats'
72 :
73 33 : ierr = 0
74 33 : nz = s% nz
75 :
76 33 : min_conv_vel_for_convective_mixing_type = 1d0 ! make this a control parameter
77 :
78 33 : RSP2_or_RSP = s% RSP_flag .or. s% RSP2_flag
79 :
80 33 : if (s% doing_timing) call start_time(s, time0, total)
81 :
82 33 : if (s% RTI_flag) then
83 0 : call set_RTI_mixing_info(s, ierr)
84 0 : if (failed('set_RTI_mixing_info')) return
85 0 : call set_dPdr_dRhodr_info(s, ierr)
86 0 : if (failed('set_dPdr_dRhodr_info')) return
87 : end if
88 :
89 33 : max_conv_bdy = 10 ! will automatically be increased if necessary
90 33 : max_mix_bdy = 10 ! will automatically be increased if necessary
91 :
92 33 : s% num_conv_boundaries = 0
93 33 : if (.not. associated(s% conv_bdy_loc)) allocate(s% conv_bdy_loc(max_conv_bdy))
94 33 : if (.not. associated(s% conv_bdy_q)) allocate(s% conv_bdy_q(max_conv_bdy))
95 33 : if (.not. associated(s% top_conv_bdy)) allocate(s% top_conv_bdy(max_conv_bdy))
96 33 : if (.not. associated(s% burn_h_conv_region)) allocate(s% burn_h_conv_region(max_conv_bdy))
97 33 : if (.not. associated(s% burn_he_conv_region)) allocate(s% burn_he_conv_region(max_conv_bdy))
98 33 : if (.not. associated(s% burn_z_conv_region)) allocate(s% burn_z_conv_region(max_conv_bdy))
99 :
100 33 : s% num_mix_boundaries = 0
101 33 : if (.not. associated(s% mix_bdy_loc)) allocate(s% mix_bdy_loc(max_mix_bdy))
102 33 : if (.not. associated(s% mix_bdy_q)) allocate(s% mix_bdy_q(max_mix_bdy))
103 33 : if (.not. associated(s% top_mix_bdy)) allocate(s% top_mix_bdy(max_mix_bdy))
104 33 : if (.not. associated(s% burn_h_mix_region)) allocate(s% burn_h_mix_region(max_mix_bdy))
105 33 : if (.not. associated(s% burn_he_mix_region)) allocate(s% burn_he_mix_region(max_mix_bdy))
106 33 : if (.not. associated(s% burn_z_mix_region)) allocate(s% burn_z_mix_region(max_mix_bdy))
107 :
108 198 : allocate(eps_h(nz), eps_he(nz), eps_z(nz), cdc_factor(nz))
109 :
110 33 : if (.not. RSP2_or_RSP) then
111 40766 : do k = 1, nz
112 40766 : s% mixing_type(k) = s% mlt_mixing_type(k)
113 : end do
114 : end if
115 :
116 33 : cdc_factor(1) = 1d0
117 40733 : do k = 2, nz
118 : rho_face = (s% dq(k-1)*s% rho(k) + s% dq(k)*s% rho(k-1))/&
119 40700 : (s% dq(k-1) + s% dq(k))
120 40700 : f = pi4*s% r(k)*s% r(k)*rho_face
121 40733 : cdc_factor(k) = f*f
122 : end do
123 :
124 33 : if (s% RSP_flag) then
125 0 : do k = 1, nz
126 0 : s% mixing_type(k) = no_mixing
127 0 : s% D_mix(k) = 0d0
128 0 : s% cdc(k) = 0d0
129 0 : s% conv_vel(k) = 0d0
130 : end do
131 33 : else if (s% RSP2_flag) then
132 0 : do k = 1, nz
133 0 : s% conv_vel(k) = get_RSP2_conv_velocity(s,k)
134 0 : s% D_mix(k) = s% conv_vel(k)*s% mixing_length_alpha*s% Hp_face(k)/3d0
135 0 : s% cdc(k) = cdc_factor(k)*s% D_mix(k)
136 0 : L_val = max(1d-99,abs(s% L(k)))
137 0 : if (abs(s% Lt(k)) > &
138 0 : L_val*s% RSP2_min_Lt_div_L_for_overshooting_mixing_type) then
139 0 : s% mixing_type(k) = overshoot_mixing
140 0 : else if (abs(s% Lc(k)) > &
141 : L_val*s% RSP2_min_Lc_div_L_for_convective_mixing_type) then
142 0 : s% mixing_type(k) = convective_mixing
143 : else
144 0 : s% mixing_type(k) = no_mixing
145 : end if
146 : end do
147 : else
148 40766 : do k = 1, nz
149 40733 : s% D_mix(k) = s% mlt_D(k)
150 40733 : s% cdc(k) = s% mlt_cdc(k)
151 40766 : s% conv_vel(k) = s% mlt_vc(k)
152 : end do
153 : end if
154 :
155 33 : call check('after get mlt_D')
156 :
157 33 : if (s% remove_mixing_glitches .and. .not. RSP2_or_RSP) then
158 :
159 33 : call remove_tiny_mixing(s, ierr)
160 33 : if (failed('remove_tiny_mixing')) return
161 :
162 33 : call remove_mixing_singletons(s, ierr)
163 33 : if (failed('remove_mixing_singletons')) return
164 :
165 33 : call close_convection_gaps(s, ierr)
166 33 : if (failed('close_convection_gaps')) return
167 :
168 33 : call close_thermohaline_gaps(s, ierr)
169 33 : if (failed('close_thermohaline_gaps')) return
170 :
171 33 : call remove_thermohaline_dropouts(s, ierr)
172 33 : if (failed('remove_thermohaline_dropouts')) return
173 :
174 33 : call close_semiconvection_gaps(s, ierr)
175 33 : if (failed('close_semiconvection_gaps')) return
176 :
177 33 : call remove_embedded_semiconvection(s, ierr)
178 33 : if (failed('remove_embedded_semiconvection')) return
179 :
180 : end if
181 :
182 33 : call check('after get remove_mixing_glitches')
183 :
184 33 : call do_mix_envelope(s)
185 :
186 : ! Keep convection velocities before adding other mixing sources.
187 33 : if (.not. s% RSP2_flag) then
188 40766 : do k=1,nz
189 40766 : if (s% mixing_type(k) /= convective_mixing) s% conv_vel(k) = 0d0
190 : end do
191 : end if
192 :
193 40766 : do k=1,s% nz
194 : eps_h(k) = s% eps_nuc_categories(ipp,k) + &
195 40733 : s% eps_nuc_categories(icno,k)
196 40733 : eps_he(k) = s% eps_nuc_categories(i3alf,k)
197 40766 : eps_z(k) = s% eps_nuc(k) - (eps_h(k) + eps_he(k))
198 : end do
199 :
200 33 : if (.not. s% RSP_flag) then
201 :
202 33 : call set_cz_boundary_info(s, ierr)
203 33 : if (failed('set_cz_boundary_info')) return
204 :
205 : call locate_convection_boundaries( &
206 : s, nz, eps_h, eps_he, eps_z, s% mstar, &
207 33 : s% q, s% cdc, ierr)
208 33 : if (failed('locate_convection_boundaries')) return
209 :
210 33 : call add_predictive_mixing(s, ierr)
211 33 : if (failed('add_predictive_mixing')) return
212 :
213 : end if
214 :
215 33 : call check('after add_predictive_mixing')
216 :
217 : ! NB: re-call locate_convection_boundaries to take into
218 : ! account changes from add_predictive_mixing
219 :
220 33 : if (.not. s% RSP_flag) then
221 :
222 : call locate_convection_boundaries( &
223 : s, nz, eps_h, eps_he, eps_z, s% mstar, &
224 33 : s% q, s% cdc, ierr)
225 33 : if (failed('locate_convection_boundaries')) return
226 :
227 33 : call locate_mixing_boundaries(s, eps_h, eps_he, eps_z, ierr)
228 33 : if (failed('locate_mixing_boundaries')) return
229 :
230 33 : call add_overshooting(s, ierr)
231 33 : if (failed('add_overshooting')) return
232 :
233 : end if
234 :
235 33 : call check('after add_overshooting')
236 :
237 33 : call add_RTI_turbulence(s, ierr)
238 33 : if (failed('add_RTI_turbulence')) return
239 :
240 33 : call s% other_after_set_mixing_info(s% id, ierr)
241 33 : if (failed('other_after_set_mixing_info')) return
242 :
243 33 : if (.not. (skip_set_cz_bdy_mass .or. RSP2_or_RSP)) then
244 11 : call set_cz_bdy_mass(s, ierr)
245 11 : if (failed('set_cz_bdy_mass')) return
246 : end if
247 :
248 33 : if (s% set_min_D_mix .and. s% ye(nz) >= s% min_center_Ye_for_min_D_mix) then
249 0 : do k=1,nz
250 0 : if (s% D_mix(k) >= s% min_D_mix) cycle
251 0 : if (s% m(k) > s% mass_upper_limit_for_min_D_mix*Msun) cycle
252 0 : if (s% m(k) < s% mass_lower_limit_for_min_D_mix*Msun) cycle
253 0 : s% D_mix(k) = s% min_D_mix
254 0 : s% mixing_type(k) = minimum_mixing
255 : end do
256 : end if
257 :
258 33 : if (s% set_min_D_mix_below_Tmax) then
259 0 : Tmax = -1
260 0 : k_Tmax = -1
261 0 : do k=1,nz
262 0 : if (s% T(k) > Tmax) then
263 0 : Tmax = s% T(k)
264 0 : k_Tmax = k
265 : end if
266 : end do
267 0 : do k=k_Tmax+1,nz
268 0 : if (s% D_mix(k) >= s% min_D_mix_below_Tmax) cycle
269 0 : s% D_mix(k) = s% min_D_mix_below_Tmax
270 0 : s% mixing_type(k) = minimum_mixing
271 : end do
272 : end if
273 :
274 33 : if (s% set_min_D_mix_in_H_He) then
275 0 : do k=1,nz
276 0 : if (s% X(k) + s% Y(k) < 0.5d0) exit
277 0 : if (s% D_mix(k) >= s% min_D_mix_in_H_He) cycle
278 0 : s% D_mix(k) = s% min_D_mix_in_H_He
279 0 : s% mixing_type(k) = minimum_mixing
280 : end do
281 : end if
282 :
283 33 : if (s% use_other_D_mix) then
284 0 : call s% other_D_mix(s% id, ierr)
285 0 : if (failed('other_D_mix')) return
286 : end if
287 :
288 40766 : do k=1,nz
289 40766 : s% D_mix_non_rotation(k) = s% D_mix(k)
290 : end do
291 :
292 33 : call check('before rotation_flag')
293 :
294 33 : if (s% rotation_flag) then
295 :
296 0 : call update_rotation_mixing_info(s,ierr)
297 0 : if (failed('update_rotation_mixing_info')) return
298 :
299 0 : do k = 2, nz
300 0 : if (s% D_mix(k) < 1d-10) s% D_mix(k) = 0d0
301 0 : s% cdc(k) = s% D_mix(k)*cdc_factor(k)
302 0 : if (s% D_mix(k) /= 0 .and. s% D_mix_non_rotation(k) < s% D_mix_rotation(k)) then
303 0 : s% mixing_type(k) = rotation_mixing
304 : end if
305 : end do
306 0 : s% cdc(1) = s% cdc(2)
307 :
308 : end if
309 :
310 33 : call check('after update_rotation_mixing_info')
311 :
312 33 : region_bottom_q = s% D_mix_zero_region_bottom_q
313 33 : region_top_q = s% D_mix_zero_region_top_q
314 :
315 33 : if (s% dq_D_mix_zero_at_H_He_crossover > 0d0) then
316 0 : i_h1 = s% net_iso(ih1)
317 0 : i_he4 = s% net_iso(ihe4)
318 0 : if (i_h1 > 0 .and. i_he4 > 0) then
319 0 : do k=2,nz
320 0 : if (s% xa(i_h1,k-1) > s% xa(i_he4,k-1) .and. &
321 0 : s% xa(i_h1,k) <= s% xa(i_he4,k)) then ! crossover
322 : region_bottom_q = &
323 0 : s% q(k) - 0.5d0*s% dq_D_mix_zero_at_H_He_crossover
324 : region_top_q = &
325 0 : s% q(k) + 0.5d0*s% dq_D_mix_zero_at_H_He_crossover
326 0 : exit
327 : end if
328 : end do
329 : end if
330 : end if
331 :
332 33 : if (region_bottom_q < region_top_q) then
333 0 : do k=2,nz
334 0 : if (s% q(k) >= region_bottom_q .and. s% q(k) <= region_top_q) then
335 0 : s% D_mix(k) = 0d0
336 0 : s% mixing_type(k) = no_mixing
337 : end if
338 : end do
339 : end if
340 :
341 33 : region_bottom_q = 1d99
342 33 : region_top_q = -1d99
343 33 : if (s% dq_D_mix_zero_at_H_C_crossover > 0d0) then
344 0 : i_h1 = s% net_iso(ih1)
345 0 : i_c12 = s% net_iso(ic12)
346 0 : if (i_h1 > 0 .and. i_c12 > 0) then
347 0 : do k=2,nz
348 0 : if (s% xa(i_h1,k-1) > s% xa(i_c12,k-1) .and. &
349 0 : s% xa(i_h1,k) <= s% xa(i_c12,k)) then ! crossover
350 : region_bottom_q = &
351 0 : s% q(k) - 0.5d0*s% dq_D_mix_zero_at_H_C_crossover
352 : region_top_q = &
353 0 : s% q(k) + 0.5d0*s% dq_D_mix_zero_at_H_C_crossover
354 0 : exit
355 : end if
356 : end do
357 : end if
358 : end if
359 :
360 0 : if (region_bottom_q < region_top_q) then
361 0 : do k=2,nz
362 0 : if (s% q(k) >= region_bottom_q .and. s% q(k) <= region_top_q) then
363 0 : s% D_mix(k) = 0d0
364 0 : s% mixing_type(k) = no_mixing
365 : end if
366 : end do
367 : end if
368 :
369 : ! set these just for plotting. not used.
370 33 : s% mixing_type(1) = s% mixing_type(2)
371 33 : s% D_mix(1) = s% D_mix(2)
372 33 : s% conv_vel(1) = 0d0
373 :
374 33 : call check('final')
375 33 : if (failed('set_mixing_info')) return
376 :
377 33 : if (s% doing_timing) &
378 0 : call update_time(s, time0, total, s% time_set_mixing_info)
379 :
380 99 : s% have_mixing_info = .true.
381 :
382 : contains
383 :
384 539 : logical function failed(str)
385 : character (len=*), intent(in) :: str
386 539 : if (ierr == 0) then
387 : failed = .false.
388 : return
389 : end if
390 0 : if (s% report_ierr) &
391 0 : write(*,*) 'set_mixing_info failed in call to ' // trim(str)
392 : failed = .true.
393 : end function failed
394 :
395 231 : subroutine check(str)
396 : character(len=*) :: str
397 : integer :: k
398 : include 'formats'
399 285362 : do k = 1, s% nz
400 285362 : if (is_bad_num(s% D_mix(k))) then
401 0 : ierr = -1
402 0 : if (s% report_ierr) then
403 0 : write(*,3) trim(str) // ' mixing_type, D_mix', k, s% mixing_type(k), s% D_mix(k)
404 0 : if (s% rotation_flag) then
405 0 : if (is_bad_num(s% D_mix_non_rotation(k))) &
406 0 : write(*,2) 's% D_mix_non_rotation(k)', k, s% D_mix_non_rotation(k)
407 0 : if (is_bad_num(s% D_visc(k))) write(*,2) 's% D_visc(k)', k, s% D_visc(k)
408 0 : if (is_bad_num(s% D_DSI(k))) write(*,2) 's% D_DSI(k)', k, s% D_DSI(k)
409 0 : if (is_bad_num(s% D_SH(k))) write(*,2) 's% D_SH(k)', k, s% D_SH(k)
410 0 : if (is_bad_num(s% D_SSI(k))) write(*,2) 's% D_SSI(k)', k, s% D_SSI(k)
411 0 : if (is_bad_num(s% D_ES(k))) write(*,2) 's% D_ES(k)', k, s% D_ES(k)
412 0 : if (is_bad_num(s% D_GSF(k))) write(*,2) 's% D_GSF(k)', k, s% D_GSF(k)
413 0 : if (is_bad_num(s% D_ST(k))) write(*,2) 's% D_ST(k)', k, s% D_ST(k)
414 : end if
415 : end if
416 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'set mixing info')
417 : end if
418 : end do
419 231 : end subroutine check
420 :
421 : end subroutine set_mixing_info
422 :
423 :
424 33 : subroutine set_cz_boundary_info(s, ierr)
425 : type (star_info), pointer :: s
426 : integer, intent(out) :: ierr
427 :
428 : integer :: k, mt1, mt2, nz
429 : real(dp) :: dg0, dg1
430 :
431 : include 'formats'
432 :
433 : ! NOTE: this routine must be called BEFORE overshooting is done.
434 :
435 : ! convective zone boundary is where gradL = gradr
436 : ! semi-convective zone boundary is where grada_face = gradr
437 :
438 33 : ierr = 0
439 33 : nz = s% nz
440 40766 : s% cz_bdy_dq(1:nz) = 0d0
441 :
442 40733 : do k = 2, nz
443 40700 : mt1 = s% mixing_type(k-1)
444 40700 : mt2 = s% mixing_type(k)
445 40700 : if (mt1 == mt2) cycle
446 105 : if (mt2 == convective_mixing .or. mt1 == convective_mixing) then
447 105 : dg0 = s% gradL(k-1) - s% gradr(k-1)
448 105 : dg1 = s% gradL(k) - s% gradr(k)
449 0 : else if (mt2 == semiconvective_mixing .or. mt1 == semiconvective_mixing) then
450 0 : dg0 = s% grada_face(k-1) - s% gradr(k-1)
451 0 : dg1 = s% grada_face(k) - s% gradr(k)
452 : else
453 : cycle
454 : end if
455 105 : if (dg0*dg1 >= 0) cycle
456 105 : s% cz_bdy_dq(k-1) = find0(0d0,dg0,s% dq(k-1),dg1)
457 138 : if (s% cz_bdy_dq(k-1) < 0d0 .or. s% cz_bdy_dq(k-1) > s% dq(k-1)) then
458 0 : write(*,2) 'bad cz_bdy_dq', k-1, s% cz_bdy_dq(k-1), s% dq(k-1)
459 0 : call mesa_error(__FILE__,__LINE__,'set_cz_boundary_info')
460 : ierr = -1
461 : return
462 : end if
463 : end do
464 :
465 : end subroutine set_cz_boundary_info
466 :
467 :
468 11 : subroutine set_cz_bdy_mass(s, ierr)
469 : type (star_info), pointer :: s
470 : integer, intent(out) :: ierr
471 :
472 : logical :: in_convective_region
473 : integer :: k, j, nz
474 : logical, parameter :: dbg = .false.
475 :
476 : include 'formats'
477 11 : ierr = 0
478 11 : nz = s% nz
479 :
480 1111 : s% cz_top_mass(:) = s% mstar
481 1111 : s% cz_bot_mass(:) = s% mstar
482 :
483 11 : s% n_conv_regions = 0
484 11 : in_convective_region = (s% mixing_type(nz) == convective_mixing)
485 11 : if (in_convective_region) then
486 11 : s% n_conv_regions = 1
487 11 : s% cz_bot_mass(1) = s% M_center
488 : end if
489 :
490 : if (dbg) write(*,*) 'initial in_convective_region', in_convective_region
491 :
492 13076 : do k=nz-1, 2, -1
493 13076 : if (in_convective_region) then
494 1894 : if (s% mixing_type(k) /= convective_mixing) then ! top of convective region
495 : s% cz_top_mass(s% n_conv_regions) = &
496 23 : s% M_center + (s% q(k) - s% cz_bdy_dq(k))*s% xmstar
497 23 : in_convective_region = .false.
498 : end if
499 : else
500 11171 : if (s% mixing_type(k) == convective_mixing) then ! bottom of convective region
501 12 : if (s% n_conv_regions < max_num_mixing_regions) then
502 12 : s% n_conv_regions = s% n_conv_regions + 1
503 : s% cz_bot_mass(s% n_conv_regions) = &
504 12 : s% M_center + (s% q(k) - s% cz_bdy_dq(k))*s% xmstar
505 : end if
506 : in_convective_region = .true.
507 : end if
508 : end if
509 : end do
510 11 : if (in_convective_region) then
511 0 : s% cz_top_mass(s% n_conv_regions) = s% mstar
512 : end if
513 :
514 11 : s% have_new_cz_bdy_info = .true.
515 :
516 : if (dbg) then
517 : write(*,'(A)')
518 : write(*,2) 'set_mixing_info s% n_conv_regions', s% n_conv_regions
519 : do j = 1, s% n_conv_regions
520 : write(*,2) 'conv region', j, s% cz_bot_mass(j)/Msun, s% cz_top_mass(j)/Msun
521 : end do
522 : write(*,'(A)')
523 : end if
524 :
525 11 : end subroutine set_cz_bdy_mass
526 :
527 :
528 66 : subroutine locate_convection_boundaries( &
529 66 : s, nz, eps_h, eps_he, eps_z, mstar, q, cdc, ierr)
530 : type (star_info), pointer :: s
531 : integer, intent(in) :: nz
532 : real(dp), dimension(:), intent(in) :: eps_h, eps_he, eps_z
533 : real(dp), intent(in) :: mstar
534 : real(dp), pointer, dimension(:) :: q, cdc
535 : integer, intent(out) :: ierr
536 :
537 : logical :: in_convective_region
538 : integer :: k, k_bot, i, j, iounit, max_conv_bdy
539 : real(dp) :: turnover_time, bot_Hp, bot_r, top_Hp, top_r, dr
540 :
541 : logical :: dbg
542 : logical, parameter :: write_debug = .false.
543 :
544 : include 'formats'
545 :
546 66 : ierr = 0
547 66 : dbg = .false.
548 :
549 : if (write_debug) then
550 : open(newunit=iounit, file=trim('debug.data'), action='write', iostat=ierr)
551 : if (ierr /= 0) then
552 : write(*, *) 'open debug.data failed'
553 : return
554 : end if
555 : write(*,*) 'write debug.data'
556 : write(iounit,*) 'nz', nz
557 : write(iounit,1) 'mstar', mstar
558 : do k=1,nz
559 : write(iounit,2) 'q', k, q(k)
560 : write(iounit,2) 'cdc', k, cdc(k)
561 : write(iounit,2) 'eps_h', k, eps_h(k)
562 : write(iounit,2) 'eps_he', k, eps_he(k)
563 : write(iounit,2) 'eps_z', k, eps_z(k)
564 : end do
565 : end if
566 :
567 66 : max_conv_bdy = size(s% conv_bdy_q, dim=1)
568 726 : s% conv_bdy_q(:) = 0
569 726 : s% conv_bdy_loc(:) = 0
570 726 : s% top_conv_bdy(:) = .false.
571 726 : s% burn_h_conv_region(:) = .false.
572 726 : s% burn_he_conv_region(:) = .false.
573 726 : s% burn_z_conv_region(:) = .false.
574 66 : bot_Hp = 0; bot_r = 0; top_Hp = 0; top_r = 0; dr = 0
575 :
576 66 : s% num_conv_boundaries = 0
577 66 : in_convective_region = (s% mixing_type(nz) == convective_mixing)
578 66 : k_bot = nz
579 66 : turnover_time = 0
580 :
581 81400 : do k=nz-1, 2, -1
582 81400 : if (in_convective_region) then
583 12002 : if (s% mixing_type(k) /= convective_mixing) then
584 138 : call end_of_convective_region
585 81334 : else if(s% conv_vel(k) /= 0d0) then
586 81334 : turnover_time = turnover_time + (s% rmid(k-1) - s% rmid(k))/s% conv_vel(k)
587 : end if
588 : else ! in non-convective region
589 69332 : if (s% mixing_type(k) == convective_mixing) then ! start of a convective region
590 72 : if (s% num_conv_boundaries == max_conv_bdy) then
591 0 : call realloc(ierr)
592 0 : if (ierr /= 0) then
593 : return
594 : end if
595 : end if
596 72 : s% num_conv_boundaries = s% num_conv_boundaries+1
597 72 : i = s% num_conv_boundaries
598 72 : k_bot = k+1
599 72 : if (k == 1) then
600 0 : s% conv_bdy_q(i) = 1
601 : else ! bottom of region is between k+1 and k
602 72 : s% conv_bdy_q(i) = s% q(k) - s% cz_bdy_dq(k)
603 : end if
604 72 : s% top_conv_bdy(i) = .false.
605 72 : s% conv_bdy_loc(i) = k_bot
606 : ! bottom of region is between k_bot and k_bot-1
607 72 : in_convective_region = .true.
608 72 : bot_r = s% r(k_bot)
609 : !bot_Hp = scale_height(s,k_bot,.true.)
610 72 : bot_Hp = s% scale_height(k_bot)
611 72 : turnover_time = 0
612 : end if
613 : end if
614 : end do
615 :
616 66 : if (in_convective_region) then
617 0 : k = 1 ! end at top
618 0 : call end_of_convective_region
619 : end if
620 :
621 : if (write_debug) then
622 : write(iounit,*) 's% num_conv_boundaries', s% num_conv_boundaries
623 : do j=1,s% num_conv_boundaries
624 : write(iounit,2) 's% conv_bdy_q', j, s% conv_bdy_q(j)
625 : write(iounit,*) 's% top_conv_bdy', j, s% top_conv_bdy(j)
626 : write(iounit,*) 's% burn_h_conv_region', j, s% burn_h_conv_region(j)
627 : write(iounit,*) 's% burn_he_conv_region', j, s% burn_he_conv_region(j)
628 : write(iounit,*) 's% burn_z_conv_region', j, s% burn_z_conv_region(j)
629 : write(iounit,*) 's% conv_bdy_loc', j, s% conv_bdy_loc(j)
630 : end do
631 : close(iounit)
632 : end if
633 :
634 : if (dbg) then
635 : write(*,*) 's% num_conv_boundaries', s% num_conv_boundaries
636 : do j=1,s% num_conv_boundaries
637 : write(*,*) 's% conv_bdy_q', j, s% conv_bdy_q(j)
638 : write(*,*) 's% top_conv_bdy', j, s% top_conv_bdy(j)
639 : write(*,*) 's% burn_h_conv_region', j, s% burn_h_conv_region(j)
640 : write(*,*) 's% burn_he_conv_region', j, s% burn_he_conv_region(j)
641 : write(*,*) 's% burn_z_conv_region', j, s% burn_z_conv_region(j)
642 : write(*,*) 's% conv_bdy_loc', j, s% conv_bdy_loc(j)
643 : write(*,*) 'mixing type', s% mixing_type(s% conv_bdy_loc(j)-3:s% conv_bdy_loc(j)+3)
644 : end do
645 : write(*,'(A)')
646 : write(*,3) 'mixing_type', 1152, s% mixing_type(1152)
647 : call mesa_error(__FILE__,__LINE__,'locate_convection_boundaries')
648 : end if
649 :
650 :
651 : contains
652 :
653 :
654 138 : subroutine end_of_convective_region()
655 : integer :: max_logT_loc, kk
656 : real(dp) :: max_logT, max_X, max_Y, Hp, max_eps
657 : logical :: end_dbg
658 :
659 : 9 format(a40, 3i7, 99(1pd26.16))
660 : include 'formats'
661 :
662 138 : in_convective_region = .false.
663 :
664 138 : end_dbg = .false.
665 :
666 138 : top_r = s% r(k)
667 138 : top_Hp = s% scale_height(k)
668 138 : dr = top_r - bot_r
669 138 : Hp = (bot_Hp + top_Hp)/2
670 :
671 138 : if (dr/Hp < s% prune_bad_cz_min_Hp_height .and. s% prune_bad_cz_min_Hp_height > 0) then
672 0 : max_eps = maxval(eps_h(k:k_bot) + eps_he(k:k_bot) + eps_z(k:k_bot))
673 : if (end_dbg) write(*,3) 'max_eps', k, k_bot, max_eps, &
674 : exp10(s% prune_bad_cz_min_log_eps_nuc)
675 : if (max_eps < exp10(s% prune_bad_cz_min_log_eps_nuc) &
676 0 : .and. all(s% mixing_type(k+1:k_bot-1) /= thermohaline_mixing)) then
677 0 : do kk = k, k_bot ! this includes the radiative points at boundaries
678 0 : call set_use_gradr(s,kk)
679 0 : s% cdc(kk) = 0
680 0 : s% D_mix(kk) = 0
681 0 : s% conv_vel(kk) = 0
682 0 : s% mixing_type(kk) = no_mixing
683 : end do
684 0 : if (s% num_conv_boundaries > 0) &
685 0 : s% num_conv_boundaries = s% num_conv_boundaries-1
686 0 : return
687 : end if
688 : end if
689 :
690 138 : if (s% num_conv_boundaries == max_conv_bdy) then
691 0 : call realloc(ierr)
692 0 : if (ierr /= 0) return
693 : end if
694 :
695 138 : s% num_conv_boundaries = s% num_conv_boundaries+1
696 138 : i = s% num_conv_boundaries
697 :
698 : ! check for burning in region
699 138 : max_logT = -1d99
700 138 : max_X = 0d0
701 138 : max_Y = 0d0
702 138 : max_logT_loc = 0
703 12422 : do kk=k,min(nz,k_bot+1)
704 12284 : if (s% lnT(kk)/ln10 > max_logT) then
705 : max_logT = s% lnT(kk)/ln10
706 : max_logT_loc = kk
707 : end if
708 12284 : if (s% X(kk) > max_X) then
709 : max_X = s% X(kk)
710 : end if
711 12422 : if (s% Y(kk) > max_Y) then
712 : max_Y = s% Y(kk)
713 : end if
714 : end do
715 : if (max_logT > s% burn_z_mix_region_logT &
716 138 : .and. max_Y < s% max_Y_for_burn_z_mix_region) then
717 0 : s% burn_z_conv_region(i) = .true.
718 0 : if (i > 1) s% burn_z_conv_region(i-1) = .true.
719 : !write(*,*) 'burn z mix region', i, &
720 : ! s% burn_z_mix_region_logT, max_logT, s% m(max_logT_loc)/Msun
721 : else if (max_logT > s% burn_he_mix_region_logT &
722 138 : .and. max_X < s% max_X_for_burn_he_mix_region) then
723 0 : s% burn_he_conv_region(i) = .true.
724 0 : if (i > 1) s% burn_he_conv_region(i-1) = .true.
725 : !write(*,*) 'burn he mix region', i, &
726 : ! s% burn_he_mix_region_logT, max_logT, s% m(max_logT_loc)/Msun
727 138 : else if (max_logT > s% burn_h_mix_region_logT) then
728 66 : s% burn_h_conv_region(i) = .true.
729 66 : if (i > 1) s% burn_h_conv_region(i-1) = .true.
730 : !write(*,*) 'burn h mix region', i, &
731 : ! s% burn_h_mix_region_logT, max_logT, s% m(max_logT_loc)/Msun
732 : end if
733 :
734 138 : if (k == 1) then
735 0 : s% conv_bdy_q(i) = 1
736 : else
737 : ! top of region is between k+1 and k
738 138 : s% conv_bdy_q(i) = s% q(k) - s% cz_bdy_dq(k)
739 : end if
740 138 : s% top_conv_bdy(i) = .true.
741 138 : s% conv_bdy_loc(i) = k
742 :
743 : end subroutine end_of_convective_region
744 :
745 :
746 0 : subroutine realloc(ierr)
747 : use utils_lib
748 : integer, intent(out) :: ierr
749 :
750 : integer :: sz
751 :
752 0 : sz = size(s% conv_bdy_q, dim=1)
753 :
754 : ierr = 0
755 0 : max_conv_bdy = 2*(10+max_conv_bdy)
756 :
757 0 : call realloc_double(s% conv_bdy_q,max_conv_bdy,ierr)
758 0 : if (ierr /= 0) return
759 :
760 0 : call realloc_integer(s% conv_bdy_loc,max_conv_bdy,ierr)
761 0 : if (ierr /= 0) return
762 :
763 0 : call realloc_logical(s% top_conv_bdy,max_conv_bdy,ierr)
764 0 : if (ierr /= 0) return
765 :
766 0 : call realloc_logical(s% burn_h_conv_region,max_conv_bdy,ierr)
767 0 : if (ierr /= 0) return
768 :
769 0 : call realloc_logical(s% burn_he_conv_region,max_conv_bdy,ierr)
770 0 : if (ierr /= 0) return
771 :
772 0 : call realloc_logical(s% burn_z_conv_region,max_conv_bdy,ierr)
773 0 : if (ierr /= 0) return
774 :
775 0 : s% conv_bdy_q(sz+1:max_conv_bdy) = 0
776 0 : s% conv_bdy_loc(sz+1:max_conv_bdy) = 0
777 0 : s% top_conv_bdy(sz+1:max_conv_bdy) = .false.
778 0 : s% burn_h_conv_region(sz+1:max_conv_bdy) = .false.
779 0 : s% burn_he_conv_region(sz+1:max_conv_bdy) = .false.
780 0 : s% burn_z_conv_region(sz+1:max_conv_bdy) = .false.
781 :
782 : end subroutine realloc
783 :
784 : end subroutine locate_convection_boundaries
785 :
786 :
787 30 : subroutine set_use_gradr(s,k)
788 : use turb_info, only: switch_to_radiative
789 : type (star_info), pointer :: s
790 : integer, intent(in) :: k
791 30 : call switch_to_radiative(s,k)
792 : end subroutine set_use_gradr
793 :
794 :
795 33 : subroutine locate_mixing_boundaries(s, eps_h, eps_he, eps_z, ierr)
796 : type (star_info), pointer :: s
797 : real(dp), dimension(:), intent(in) :: eps_h, eps_he, eps_z
798 : integer, intent(out) :: ierr
799 :
800 : logical :: in_mixing_region
801 : integer :: k, k_bot, i, max_mix_bdy, nz
802 :
803 : logical, parameter :: dbg = .false.
804 :
805 : include 'formats'
806 :
807 33 : ierr = 0
808 33 : nz = s% nz
809 :
810 33 : max_mix_bdy = size(s% mix_bdy_q, dim=1)
811 363 : s% mix_bdy_q(:) = 0
812 363 : s% mix_bdy_loc(:) = 0
813 363 : s% top_mix_bdy(:) = .false.
814 363 : s% burn_h_mix_region(:) = .false.
815 363 : s% burn_he_mix_region(:) = .false.
816 363 : s% burn_z_mix_region(:) = .false.
817 :
818 33 : s% num_mix_boundaries = 0
819 33 : s% num_mix_regions = 0
820 33 : in_mixing_region = (s% mixing_type(nz) /= no_mixing)
821 33 : k_bot = nz
822 40700 : do k=nz-1, 2, -1
823 40700 : if (in_mixing_region) then
824 6001 : if (s% mixing_type(k) == no_mixing) call end_of_mixing_region
825 : else ! in non-mixing region
826 34666 : if (s% mixing_type(k) /= no_mixing) then ! start of a mixing region
827 36 : if (s% num_mix_boundaries == max_mix_bdy) then
828 0 : call realloc(ierr)
829 0 : if (ierr /= 0) return
830 : end if
831 36 : s% num_mix_boundaries = s% num_mix_boundaries+1
832 36 : i = s% num_mix_boundaries
833 36 : k_bot = k+1
834 36 : if (k == 1) then
835 0 : s% mix_bdy_q(i) = 1
836 : else ! bottom of region is between k+1 and k
837 36 : s% mix_bdy_q(i) = s% q(k) - s% cz_bdy_dq(k)
838 : end if
839 36 : s% top_mix_bdy(i) = .false.
840 36 : s% mix_bdy_loc(i) = k_bot
841 36 : in_mixing_region = .true.
842 : end if
843 : end if
844 : end do
845 :
846 33 : if (in_mixing_region) then
847 0 : k = 1 ! end at top
848 0 : call end_of_mixing_region
849 : end if
850 :
851 :
852 : !do i=1,s% num_conv_boundaries
853 : ! write(*,*) 'locate_mixing_boundaries region burn_h he z', i, &
854 : ! s% burn_h_mix_region(i), s% burn_he_conv_region(i), s% burn_z_conv_region(i)
855 : !end do
856 :
857 : contains
858 :
859 :
860 69 : subroutine end_of_mixing_region()
861 : integer :: kk, max_logT_loc
862 : real(dp) :: max_logT, max_X, max_Y
863 :
864 : 9 format(a40, 3i7, 99(1pd26.16))
865 : include 'formats'
866 :
867 69 : in_mixing_region = .false.
868 :
869 69 : if (s% num_mix_boundaries == max_mix_bdy) then
870 0 : call realloc(ierr)
871 0 : if (ierr /= 0) return
872 : end if
873 :
874 69 : s% num_mix_regions = s% num_mix_regions+1
875 69 : s% num_mix_boundaries = s% num_mix_boundaries+1
876 69 : i = s% num_mix_boundaries
877 :
878 : ! check for burning in region
879 69 : max_logT = -1d99
880 69 : max_X = 0d0
881 69 : max_Y = 0d0
882 69 : max_logT_loc = 0
883 6211 : do kk=k,min(nz,k_bot+1)
884 6142 : if (s% lnT(kk)/ln10 > max_logT) then
885 : max_logT = s% lnT(kk)/ln10
886 : max_logT_loc = kk
887 : end if
888 6142 : if (s% X(kk) > max_X) then
889 : max_X = s% X(kk)
890 : end if
891 6211 : if (s% Y(kk) > max_Y) then
892 : max_Y = s% Y(kk)
893 : end if
894 : end do
895 : if (max_logT > s% burn_z_mix_region_logT &
896 69 : .and. max_Y < s% max_Y_for_burn_z_mix_region) then
897 0 : s% burn_z_mix_region(i) = .true.
898 0 : if (i > 1) s% burn_z_mix_region(i-1) = .true.
899 : !write(*,*) 'burn z mix region', i
900 : else if (max_logT > s% burn_he_mix_region_logT &
901 69 : .and. max_X < s% max_X_for_burn_he_mix_region) then
902 0 : s% burn_he_mix_region(i) = .true.
903 0 : if (i > 1) s% burn_he_mix_region(i-1) = .true.
904 : !write(*,*) 'burn he mix region', i
905 69 : else if (max_logT > s% burn_h_mix_region_logT) then
906 33 : s% burn_h_mix_region(i) = .true.
907 33 : if (i > 1) s% burn_h_mix_region(i-1) = .true.
908 : !write(*,*) 'burn h mix region', i
909 : end if
910 :
911 69 : if (k == 1) then
912 0 : s% mix_bdy_q(i) = 1
913 : else
914 : ! top of region is between k+1 and k
915 69 : s% mix_bdy_q(i) = s% q(k) - s% cz_bdy_dq(k)
916 : end if
917 69 : s% top_mix_bdy(i) = .true.
918 69 : s% mix_bdy_loc(i) = k
919 :
920 : end subroutine end_of_mixing_region
921 :
922 :
923 0 : subroutine realloc(ierr)
924 : use utils_lib
925 : integer, intent(out) :: ierr
926 :
927 : integer :: sz
928 :
929 0 : sz = size(s% mix_bdy_q, dim=1)
930 :
931 : ierr = 0
932 0 : max_mix_bdy = 2*(10+max_mix_bdy)
933 :
934 0 : call realloc_double(s% mix_bdy_q,max_mix_bdy,ierr)
935 0 : if (ierr /= 0) return
936 :
937 0 : call realloc_integer(s% mix_bdy_loc,max_mix_bdy,ierr)
938 0 : if (ierr /= 0) return
939 :
940 0 : call realloc_logical(s% top_mix_bdy,max_mix_bdy,ierr)
941 0 : if (ierr /= 0) return
942 :
943 0 : call realloc_logical(s% burn_h_mix_region,max_mix_bdy,ierr)
944 0 : if (ierr /= 0) return
945 :
946 0 : call realloc_logical(s% burn_he_mix_region,max_mix_bdy,ierr)
947 0 : if (ierr /= 0) return
948 :
949 0 : call realloc_logical(s% burn_z_mix_region,max_mix_bdy,ierr)
950 0 : if (ierr /= 0) return
951 :
952 0 : s% mix_bdy_q(sz+1:max_mix_bdy) = 0
953 0 : s% mix_bdy_loc(sz+1:max_mix_bdy) = 0
954 0 : s% top_mix_bdy(sz+1:max_mix_bdy) = .false.
955 0 : s% burn_h_mix_region(sz+1:max_mix_bdy) = .false.
956 0 : s% burn_he_mix_region(sz+1:max_mix_bdy) = .false.
957 0 : s% burn_z_mix_region(sz+1:max_mix_bdy) = .false.
958 :
959 : end subroutine realloc
960 :
961 :
962 : end subroutine locate_mixing_boundaries
963 :
964 :
965 33 : subroutine remove_tiny_mixing(s, ierr)
966 : type (star_info), pointer :: s
967 : integer, intent(out) :: ierr
968 :
969 : integer :: k, nz
970 : logical, parameter :: dbg = .false.
971 : real(dp) :: tiny
972 :
973 : include 'formats'
974 :
975 : if (dbg) write(*,*) 'remove_tiny_mixing'
976 :
977 33 : ierr = 0
978 33 : nz = s% nz
979 :
980 33 : tiny = s% clip_D_limit
981 40766 : do k=1,nz
982 40766 : if (s% D_mix(k) < tiny) then
983 0 : s% cdc(k) = 0
984 0 : s% D_mix(k) = 0
985 0 : s% conv_vel(k) = 0
986 0 : s% mixing_type(k) = no_mixing
987 : end if
988 : end do
989 :
990 33 : end subroutine remove_tiny_mixing
991 :
992 :
993 : ! remove single point mixing or non-mixing regions
994 : ! NOTE: does not remove thermohaline singletons
995 33 : subroutine remove_mixing_singletons(s, ierr)
996 : !use star_utils, only: scale_height
997 : type (star_info), pointer :: s
998 : integer, intent(out) :: ierr
999 :
1000 : integer :: k, nz
1001 : logical, parameter :: dbg = .false.
1002 : real(dp) :: lambda
1003 :
1004 : include 'formats'
1005 :
1006 : if (dbg) write(*,*) 'remove_mixing_singletons'
1007 :
1008 33 : ierr = 0
1009 33 : nz = s% nz
1010 :
1011 40700 : do k=2,nz-1
1012 40700 : if (s% cdc(k) == 0) then
1013 34669 : if (s% cdc(k-1) /= 0 .and. s% cdc(k+1) /= 0) then
1014 0 : s% cdc(k) = (s% cdc(k-1) + s% cdc(k+1))/2
1015 0 : s% D_mix(k) = s% cdc(k)/pow2(pi4*s% r(k)*s% r(k)*s% rho(k))
1016 0 : s% mixing_type(k) = max(s% mixing_type(k-1), s% mixing_type(k+1))
1017 0 : s% conv_vel(k) = 0d0
1018 0 : if (s% mixing_type(k) == convective_mixing) then
1019 : lambda = s% alpha_mlt(k)* &
1020 0 : (s% scale_height(k-1) + s% scale_height(k+1))/2
1021 0 : if (lambda > 0d0) s% conv_vel(k) = 3*s% D_mix(k)/lambda
1022 : end if
1023 : if (dbg) write(*,3) 'remove radiative singleton', k, nz
1024 : end if
1025 5998 : else if (s% okay_to_remove_mixing_singleton) then
1026 5998 : if (s% cdc(k-1) == 0 .and. s% cdc(k+1) == 0) then
1027 30 : call set_use_gradr(s,k)
1028 30 : s% cdc(k) = 0
1029 30 : s% D_mix(k) = 0
1030 30 : s% conv_vel(k) = 0
1031 30 : s% mixing_type(k) = no_mixing
1032 : if (dbg) write(*,3) 'remove mixing singleton', k, nz
1033 : end if
1034 : end if
1035 : end do
1036 :
1037 33 : if (s% cdc(1) == 0) then
1038 33 : if (s% cdc(2) /= 0) then
1039 0 : s% cdc(1) = s% cdc(2)
1040 0 : s% D_mix(1) = s% D_mix(2)
1041 0 : s% mixing_type(1) = s% mixing_type(2)
1042 0 : s% conv_vel(1) = 0d0
1043 0 : if (s% mixing_type(1) == convective_mixing) &
1044 0 : s% conv_vel(1) = s% conv_vel(2)
1045 : if (dbg) write(*,3) 'remove radiative singleton', 1, nz
1046 : end if
1047 : else
1048 0 : if (s% cdc(2) == 0) then
1049 0 : call set_use_gradr(s,1)
1050 0 : s% cdc(1) = 0
1051 0 : s% D_mix(1) = 0
1052 0 : s% conv_vel(1) = 0
1053 0 : s% mixing_type(1) = no_mixing
1054 : if (dbg) write(*,2) 'remove mixing singleton', 1
1055 : end if
1056 : end if
1057 :
1058 33 : if (s% cdc(nz) == 0) then
1059 0 : if (s% cdc(nz-1) /= 0) then
1060 0 : s% cdc(nz) = s% cdc(nz-1)
1061 0 : s% D_mix(nz) = s% D_mix(nz-1)
1062 0 : s% mixing_type(nz) = s% mixing_type(nz-1)
1063 0 : s% conv_vel(nz) = 0d0
1064 0 : if (s% mixing_type(nz) == convective_mixing) &
1065 0 : s% conv_vel(nz) = s% conv_vel(nz-1)
1066 : if (dbg) write(*,2) 'remove radiative singleton: s% cdc(nz-1)', nz, s% cdc(nz-1)
1067 : end if
1068 : else
1069 33 : if (s% cdc(nz-1) == 0) then
1070 0 : call set_use_gradr(s,nz)
1071 0 : s% cdc(nz) = 0
1072 0 : s% D_mix(nz) = 0
1073 0 : s% conv_vel(nz) = 0
1074 0 : s% mixing_type(nz) = no_mixing
1075 : if (dbg) write(*,2) 'remove mixing singleton: s% cdc(nz)', nz, s% cdc(nz)
1076 : end if
1077 : end if
1078 :
1079 33 : end subroutine remove_mixing_singletons
1080 :
1081 :
1082 33 : subroutine close_convection_gaps(s, ierr)
1083 : type (star_info), pointer :: s
1084 : integer, intent(out) :: ierr
1085 33 : if (s% use_other_close_gaps) then
1086 0 : call s% other_close_gaps(s% id, convective_mixing, s% min_convective_gap, ierr)
1087 : else
1088 33 : call close_gaps(s, convective_mixing, s% min_convective_gap, ierr)
1089 : end if
1090 33 : end subroutine close_convection_gaps
1091 :
1092 :
1093 33 : subroutine close_thermohaline_gaps(s, ierr)
1094 : type (star_info), pointer :: s
1095 : integer, intent(out) :: ierr
1096 33 : if (s% use_other_close_gaps) then
1097 0 : call s% other_close_gaps(s% id, thermohaline_mixing, s% min_thermohaline_gap, ierr)
1098 : else
1099 33 : call close_gaps(s, thermohaline_mixing, s% min_thermohaline_gap, ierr)
1100 : end if
1101 33 : end subroutine close_thermohaline_gaps
1102 :
1103 :
1104 33 : subroutine close_semiconvection_gaps(s, ierr)
1105 : type (star_info), pointer :: s
1106 : integer, intent(out) :: ierr
1107 33 : if (s% use_other_close_gaps) then
1108 0 : call s% other_close_gaps(s% id, semiconvective_mixing, s% min_semiconvection_gap, ierr)
1109 : else
1110 33 : call close_gaps(s, semiconvective_mixing, s% min_semiconvection_gap, ierr)
1111 : end if
1112 33 : end subroutine close_semiconvection_gaps
1113 :
1114 :
1115 99 : subroutine close_gaps(s, mix_type, min_gap, ierr)
1116 : type (star_info), pointer :: s
1117 : integer, intent(in) :: mix_type
1118 : real(dp), intent(in) :: min_gap
1119 : integer, intent(out) :: ierr
1120 :
1121 : integer :: k, nz
1122 : logical :: in_region, dbg
1123 : real(dp) :: rtop, rbot, Hp
1124 : integer :: ktop, kbot ! k's for gap
1125 : include 'formats'
1126 :
1127 99 : dbg = .false.
1128 : !dbg = (mix_type == convective_mixing)
1129 : if (dbg) write(*,*) 'close_gaps convective_mixing'
1130 : if (dbg) write(*,3) 'mixing_type', 1152, s% mixing_type(1152)
1131 99 : ierr = 0
1132 99 : if (min_gap < 0) return
1133 0 : nz = s% nz
1134 0 : in_region = (s% mixing_type(nz) == mix_type)
1135 0 : rbot = 0
1136 0 : kbot = nz
1137 0 : do k=nz-1, 2, -1
1138 0 : if (in_region) then
1139 0 : if (s% mixing_type(k) /= mix_type) then ! end of region
1140 0 : kbot = k+1
1141 0 : rbot = s% r(kbot)
1142 0 : in_region = .false.
1143 : if (dbg) write(*,2) 'end of region', kbot, rbot
1144 : end if
1145 : else
1146 0 : if (s% mixing_type(k) == mix_type) then ! start of region
1147 0 : ktop = k
1148 0 : rtop = s% r(ktop)
1149 0 : Hp = s% Peos(ktop)/(s% rho(ktop)*s% grav(ktop))
1150 : if (dbg) write(*,2) 'start of region', ktop, rtop
1151 : if (dbg) write(*,1) 'rtop - rbot < Hp*min_gap', (rtop - rbot) - Hp*min_gap, &
1152 : rtop - rbot, Hp*min_gap, Hp, min_gap, (rtop-rbot)/Hp
1153 0 : if (rtop - rbot < Hp*min_gap) then
1154 0 : if (kbot < nz) then
1155 0 : s% cdc(ktop+1:kbot-1) = (s% cdc(ktop) + s% cdc(kbot))/2
1156 : s% D_mix(ktop+1:kbot-1) = &
1157 0 : (s% D_mix(ktop) + s% D_mix(kbot))/2
1158 0 : s% conv_vel(ktop+1:kbot-1) = 0d0
1159 0 : if (mix_type == convective_mixing) &
1160 0 : s% conv_vel(ktop+1:kbot-1) = (s% conv_vel(ktop) + s% conv_vel(kbot))/2
1161 0 : s% mixing_type(ktop+1:kbot-1) = mix_type
1162 : if (dbg) write(*,3) 'close mixing gap', &
1163 : ktop+1, kbot-1, (rtop - rbot)/Hp, rtop - rbot, Hp
1164 : else
1165 0 : s% cdc(ktop+1:kbot) = s% cdc(ktop)
1166 0 : s% D_mix(ktop+1:kbot) = s% D_mix(ktop)
1167 0 : s% conv_vel(ktop+1:kbot) = 0d0
1168 0 : if (mix_type == convective_mixing) &
1169 0 : s% conv_vel(ktop+1:kbot) = s% conv_vel(ktop)
1170 0 : s% mixing_type(ktop+1:kbot) = mix_type
1171 : if (dbg) write(*,3) 'close mixing gap', &
1172 : ktop+1, kbot, (rtop - rbot)/Hp, rtop - rbot, Hp
1173 : end if
1174 : end if
1175 : in_region = .true.
1176 : end if
1177 : end if
1178 : end do
1179 : if (dbg) write(*,3) 'mixing_type', 1152, s% mixing_type(1152)
1180 : if (dbg) write(*,*) 'done close_gaps'
1181 : !if (dbg) call mesa_error(__FILE__,__LINE__,'done close_gaps')
1182 :
1183 : end subroutine close_gaps
1184 :
1185 :
1186 : ! if find radiative region embedded in thermohaline,
1187 : ! and max(gradL - grada) in region is < 1d-3
1188 : ! and region height is < min_thermohaline_dropout
1189 : ! then convert the region to thermohaline
1190 33 : subroutine remove_thermohaline_dropouts(s, ierr)
1191 : type (star_info), pointer :: s
1192 : integer, intent(out) :: ierr
1193 :
1194 : integer :: k, nz, j
1195 : logical :: in_region
1196 : real(dp) :: rtop, rbot, Hp, q_upper, q_lower, alfa, beta
1197 : integer :: ktop, kbot ! k's for gap
1198 : logical :: all_small
1199 : logical, parameter :: dbg = .false.
1200 : include 'formats'
1201 33 : ierr = 0
1202 33 : nz = s% nz
1203 33 : rbot = s% r(nz)
1204 33 : kbot = nz-1
1205 33 : in_region = (s% mixing_type(kbot) == no_mixing)
1206 33 : all_small = .false.
1207 40667 : do k=nz-2, 2, -1
1208 40667 : if (in_region) then
1209 34666 : if (s% mixing_type(k) == no_mixing) then ! check if okay
1210 34630 : if (s% gradL(k) - s% grada_face(k) > s% max_dropout_gradL_sub_grada) &
1211 0 : all_small = .false.
1212 : else ! end of radiative region
1213 36 : ktop = k+1
1214 36 : rtop = s% r(ktop)
1215 36 : Hp = s% Peos(ktop)/(s% rho(ktop)*s% grav(ktop))
1216 36 : q_upper = s% q(ktop-1)
1217 36 : q_lower = s% q(kbot+1)
1218 : if (rtop - rbot < Hp*s% min_thermohaline_dropout .and. &
1219 : s% mixing_type(ktop-1) == thermohaline_mixing .and. &
1220 : s% mixing_type(kbot+1) == thermohaline_mixing .and. &
1221 36 : q_upper - q_lower > 1d-20 .and. all_small) then
1222 0 : do j = ktop, kbot ! interpolate in q
1223 0 : alfa = (s% q(j) - q_lower)/(q_upper - q_lower)
1224 0 : beta = 1 - alfa
1225 0 : s% cdc(j) = alfa*s% cdc(ktop-1) + beta*s% cdc(kbot+1)
1226 0 : s% D_mix(j) = alfa*s% D_mix(ktop-1) + beta*s% D_mix(kbot+1)
1227 0 : s% conv_vel(j) = 0d0
1228 0 : s% mixing_type(j) = thermohaline_mixing
1229 : end do
1230 : end if
1231 : in_region = .false.
1232 : end if
1233 : else
1234 5968 : if (s% mixing_type(k) == no_mixing) then ! start of region
1235 69 : kbot = k
1236 69 : rbot = s% r(kbot)
1237 69 : in_region = .true.
1238 : all_small = &
1239 69 : (s% gradL(k) - s% grada_face(k) <= s% max_dropout_gradL_sub_grada)
1240 : end if
1241 : end if
1242 : end do
1243 :
1244 33 : end subroutine remove_thermohaline_dropouts
1245 :
1246 :
1247 33 : subroutine remove_embedded_semiconvection(s, ierr)
1248 : type (star_info), pointer :: s
1249 : integer, intent(out) :: ierr
1250 :
1251 : integer :: k, nz
1252 : logical :: in_region
1253 : integer :: kbot, ktop
1254 :
1255 : logical, parameter :: dbg = .false.
1256 :
1257 : include 'formats'
1258 :
1259 33 : ierr = 0
1260 33 : if (.not. s% remove_embedded_semiconvection) return
1261 :
1262 0 : nz = s% nz
1263 :
1264 0 : in_region = check(nz)
1265 0 : kbot = nz
1266 0 : do k=nz-1, 2, -1
1267 0 : if (in_region) then
1268 0 : if (.not. check(k)) then ! end of region
1269 0 : ktop = k+1
1270 0 : in_region = .false.
1271 0 : call clean_region
1272 : end if
1273 : else ! not in region
1274 0 : if (check(k)) then ! start of region
1275 0 : kbot = k
1276 0 : in_region = .true.
1277 : end if
1278 : end if
1279 : end do
1280 :
1281 0 : if (in_region) then
1282 0 : ktop = 1
1283 0 : call clean_region
1284 : end if
1285 :
1286 :
1287 : contains
1288 :
1289 :
1290 0 : subroutine clean_region
1291 : integer :: kbot1, ktop1
1292 : include 'formats'
1293 : if (dbg) write(*,3) 'clean_region semiconvective', kbot, ktop
1294 : ! move top to below top convective region
1295 0 : do while (s% mixing_type(ktop) == convective_mixing)
1296 0 : ktop = ktop + 1
1297 0 : if (ktop >= kbot) return
1298 : end do
1299 : if (dbg) write(*,2) 'new ktop 1', ktop
1300 : ! move top to below top semiconvective region
1301 0 : do while (s% mixing_type(ktop) == semiconvective_mixing)
1302 0 : ktop = ktop + 1
1303 0 : if (ktop >= kbot) return
1304 : end do
1305 : if (dbg) write(*,2) 'new ktop 2', ktop
1306 : ! move bot to above bottom convective region
1307 0 : do while (s% mixing_type(kbot) == convective_mixing)
1308 0 : kbot = kbot - 1
1309 0 : if (ktop >= kbot) return
1310 : end do
1311 : if (dbg) write(*,2) 'new kbot 1', kbot
1312 : ! move bot to above bottom semiconvective region
1313 0 : do while (s% mixing_type(kbot) == semiconvective_mixing)
1314 0 : kbot = kbot - 1
1315 0 : if (ktop >= kbot) return
1316 : end do
1317 : if (dbg) write(*,2) 'new kbot 2', kbot
1318 : ! convert any semiconvective region between kbot and ktop
1319 : kbot1 = kbot
1320 0 : do while (kbot1 > ktop)
1321 : ! move kbot1 to bottom of next semiconvective region
1322 0 : do while (s% mixing_type(kbot1) == convective_mixing)
1323 0 : kbot1 = kbot1 - 1
1324 0 : if (kbot1 <= ktop) return
1325 : end do
1326 : ktop1 = kbot1
1327 : ! move ktop1 to top of semiconvective region
1328 0 : do while (s% mixing_type(ktop1) == semiconvective_mixing)
1329 0 : ktop1 = ktop1 - 1
1330 0 : if (ktop1 <= ktop) return
1331 : end do
1332 0 : s% D_mix(ktop1+1:kbot1) = s% D_mix(ktop1)
1333 0 : s% mixing_type(ktop1+1:kbot1) = convective_mixing
1334 0 : s% conv_vel(ktop1+1:kbot1) = s% conv_vel(ktop1)
1335 : if (dbg) write(*,3) 'merge semiconvective island', kbot1, ktop1+1
1336 0 : kbot1 = ktop1
1337 : end do
1338 : end subroutine clean_region
1339 :
1340 :
1341 0 : logical function check(k)
1342 : integer, intent(in) :: k
1343 : check = &
1344 : (s% mixing_type(k) == semiconvective_mixing) .or. &
1345 0 : (s% mixing_type(k) == convective_mixing)
1346 0 : end function check
1347 :
1348 : end subroutine remove_embedded_semiconvection
1349 :
1350 :
1351 33 : subroutine do_mix_envelope(s)
1352 : type (star_info), pointer :: s
1353 : real(dp) :: T_mix_limit
1354 : integer :: j, k, i, nz
1355 :
1356 : include 'formats'
1357 :
1358 33 : T_mix_limit = s% T_mix_limit
1359 : !write(*,1) 'T_mix_limit', T_mix_limit
1360 33 : if (T_mix_limit <= 0) return
1361 0 : nz = s% nz
1362 0 : j = 0
1363 0 : do k = 1, nz ! search inward until find T >= T_mix_limit
1364 0 : if (s% T(k) >= T_mix_limit) then
1365 : j = k; exit
1366 : end if
1367 : end do
1368 0 : if (j==0) j=nz ! all T < T_mix_limit
1369 : ! find base of innermost convection that has T < T_mix_limit
1370 0 : i = 0
1371 0 : do k = j, 1, -1
1372 0 : if (s% mixing_type(k) == convective_mixing) then
1373 : i = k; exit
1374 : end if
1375 : end do
1376 0 : if (i == 0) then
1377 : return ! no convection in region with T < T_mix_limit
1378 : end if
1379 : ! extend convection region to surface
1380 0 : j = maxloc(s% cdc(1:i), dim=1, mask=s% mixing_type(1:i) == convective_mixing)
1381 0 : s% conv_vel(1:i) = s% conv_vel(j)
1382 0 : j = maxloc(s% cdc(1:i), dim=1)
1383 0 : s% cdc(1:i) = s% cdc(j)
1384 0 : do k = 1,i
1385 0 : s% D_mix(k) = s% D_mix(j)
1386 : end do
1387 0 : s% mixing_type(1:i) = convective_mixing
1388 :
1389 : end subroutine do_mix_envelope
1390 :
1391 :
1392 11 : subroutine get_convection_sigmas(s, dt, ierr)
1393 : type (star_info), pointer :: s
1394 : real(dp), intent(in) :: dt
1395 : integer, intent(out) :: ierr
1396 :
1397 : integer :: nz, k, j, species, ktop, kbot, bdy
1398 : real(dp) :: sig_term_limit ! sig_term_limit is used to help convergence
1399 :
1400 : real(dp) :: siglim, xmstar, dq00, dqm1, cdcterm, dmavg, rho_face, &
1401 : cdc, max_sig, D, xm1, x00, dm, dX, X, cushion, limit, &
1402 : Tc, full_off, full_on, qbot, qtop, f1, f, alfa
1403 11 : real(dp), dimension(:), pointer :: sig, D_mix
1404 :
1405 : include 'formats'
1406 :
1407 11 : sig_term_limit = s% sig_term_limit
1408 :
1409 11 : ierr = 0
1410 11 : nz = s% nz
1411 11 : xmstar = s% xmstar
1412 11 : sig => s% sig
1413 11 : D_mix => s% D_mix
1414 :
1415 11 : sig(1) = 0d0
1416 13087 : do k = 2, nz
1417 13076 : D = D_mix(k)
1418 13076 : if (D == 0) then
1419 11182 : sig(k) = 0d0
1420 11182 : cycle
1421 : end if
1422 : rho_face = (s% dq(k-1)*s% rho(k) + s% dq(k)*s% rho(k-1))/ &
1423 1894 : (s% dq(k-1) + s% dq(k))
1424 1894 : f1 = pi4*s% r(k)*s% r(k)*rho_face
1425 1894 : f = f1*f1
1426 1894 : cdc = D*f
1427 1894 : cdcterm = s% mix_factor*cdc
1428 1894 : dq00 = s% dq(k)
1429 1894 : dqm1 = s% dq(k-1)
1430 1894 : dmavg = xmstar*(dq00+dqm1)/2
1431 1894 : sig(k) = cdcterm/dmavg
1432 1905 : if (is_bad_num(sig(k))) then
1433 0 : if (s% stop_for_bad_nums) then
1434 0 : write(*,2) 'sig(k)', k, sig(k)
1435 0 : call mesa_error(__FILE__,__LINE__,'get_convection_sigmas')
1436 : end if
1437 0 : sig(k) = 1d99
1438 : end if
1439 : end do
1440 :
1441 : ! can get numerical problems unless limit sig
1442 13098 : max_sig = maxval(sig(1:nz))
1443 22 : if (max_sig < 1) return
1444 13098 : do k = 1, nz
1445 13087 : s% sig_raw(k) = sig(k)
1446 13087 : if (sig(k) == 0) cycle
1447 1894 : if (k > 1) then
1448 1894 : siglim = sig_term_limit*xmstar*min(s% dq(k),s% dq(k-1))/dt
1449 : else
1450 0 : siglim = sig_term_limit*xmstar*s% dq(k)/dt
1451 : end if
1452 1905 : if (sig(k) > siglim) sig(k) = siglim
1453 : end do
1454 :
1455 11 : species = s% species
1456 : ! limit sigma to avoid negative mass fractions
1457 11 : cushion = 10d0
1458 11 : Tc = s% T(s% nz)
1459 11 : full_off = s% Tcenter_max_for_sig_min_factor_full_off
1460 11 : if (Tc <= full_off) return
1461 0 : full_on = s% Tcenter_min_for_sig_min_factor_full_on
1462 0 : limit = s% sig_min_factor_for_high_Tcenter
1463 0 : if (Tc < full_on) then
1464 0 : alfa = (Tc - full_off)/(full_on - full_off)
1465 : ! limit is full on for alfa = 1 and limit is 1 for alfa = 0
1466 0 : limit = limit*alfa + (1d0 - alfa)
1467 : end if
1468 0 : if (limit >= 1d0 .or. s% num_mix_regions == 0) return
1469 : ! boundaries are in order from center to surface
1470 : ! no bottom boundary at loc=nz included if center is mixed
1471 : ! however, do include top boundary at loc=1 if surface is mixed
1472 0 : if (s% top_mix_bdy(1)) then
1473 0 : kbot = s% nz
1474 0 : qbot = 0d0
1475 0 : bdy = 1
1476 : else
1477 0 : kbot = s% mix_bdy_loc(1)
1478 0 : qbot = s% q(kbot)
1479 0 : bdy = 2
1480 : end if
1481 0 : if (.not. s% top_mix_bdy(bdy)) then
1482 0 : call mesa_error(__FILE__,__LINE__,'bad mix bdy info 1')
1483 : end if
1484 0 : ktop = s% mix_bdy_loc(bdy)
1485 0 : qtop = s% q(ktop)
1486 0 : call do1_region
1487 11 : do while (bdy < s% num_mix_boundaries)
1488 0 : bdy = bdy+1
1489 0 : if (s% top_mix_bdy(bdy)) then
1490 0 : call mesa_error(__FILE__,__LINE__,'bad mix bdy info 2')
1491 : end if
1492 0 : kbot = s% mix_bdy_loc(bdy)
1493 0 : qbot = s% q(kbot)
1494 0 : bdy = bdy+1
1495 0 : if (.not. s% top_mix_bdy(bdy)) then
1496 0 : call mesa_error(__FILE__,__LINE__,'bad mix bdy info 3')
1497 : end if
1498 0 : ktop = s% mix_bdy_loc(bdy)
1499 0 : qtop = s% q(ktop)
1500 0 : call do1_region
1501 : end do
1502 :
1503 :
1504 : contains
1505 :
1506 :
1507 0 : subroutine do1_region
1508 : real(dp) :: delta_m, max_lim, alfa, beta, delta_m_upper, delta_m_lower
1509 0 : delta_m = s% m(ktop) - s% m(kbot)
1510 0 : delta_m_upper = s% delta_m_upper_for_sig_min_factor
1511 0 : delta_m_lower = s% delta_m_lower_for_sig_min_factor
1512 0 : if (delta_m >= delta_m_upper) then
1513 0 : max_lim = 1d0
1514 0 : else if (delta_m <= delta_m_lower) then
1515 0 : max_lim = limit
1516 : else
1517 0 : alfa = (delta_m - delta_m_lower)/(delta_m_upper - delta_m_lower)
1518 0 : beta = 1d0 - alfa
1519 0 : max_lim = alfa + beta*limit
1520 : end if
1521 0 : do k=max(2,ktop),kbot
1522 : call do1(k, max_lim, &
1523 0 : min(s% q(k) - qbot, qtop - s% q(k))*s% xmstar/Msun)
1524 : end do
1525 0 : end subroutine do1_region
1526 :
1527 :
1528 0 : subroutine do1(k, max_lim, delta_m_to_bdy)
1529 : integer, intent(in) :: k
1530 : real(dp), intent(in) :: max_lim, delta_m_to_bdy
1531 : real(dp) :: lim, max_delta_m_to_bdy
1532 : include 'formats'
1533 0 : siglim = sig(k)
1534 0 : if (siglim == 0d0) return
1535 : ! okay to increase limit up to max_lim
1536 0 : max_delta_m_to_bdy = s% max_delta_m_to_bdy_for_sig_min_factor
1537 0 : if (delta_m_to_bdy >= max_delta_m_to_bdy) return ! no change in sig
1538 0 : lim = limit + (max_lim - limit)*delta_m_to_bdy/max_delta_m_to_bdy
1539 0 : if (lim >= 1d0) return
1540 0 : do j=1,species
1541 0 : xm1 = s% xa(j,k-1)
1542 0 : x00 = s% xa(j,k)
1543 0 : if (xm1 > x00) then
1544 0 : X = xm1
1545 0 : dX = xm1 - x00
1546 0 : dm = s% dm(k-1)
1547 : else
1548 0 : X = x00
1549 0 : dX = x00 - xm1
1550 0 : dm = s% dm(k)
1551 0 : dX = -dX
1552 : end if
1553 0 : if (X < 1d-5) cycle
1554 0 : if (cushion*dt*dX*siglim > dm*X) then
1555 0 : siglim = dm*X/(dt*dX*cushion)
1556 : end if
1557 : end do
1558 0 : if (siglim > lim*sig(k)) then
1559 0 : sig(k) = siglim
1560 : else
1561 0 : sig(k) = lim*sig(k)
1562 : end if
1563 : end subroutine do1
1564 :
1565 :
1566 : end subroutine get_convection_sigmas
1567 :
1568 :
1569 0 : subroutine update_rotation_mixing_info(s, ierr)
1570 : use rotation_mix_info, only: set_rotation_mixing_info
1571 : type (star_info), pointer :: s
1572 : integer, intent(out) :: ierr
1573 :
1574 : integer :: k, nz
1575 : logical :: set_min_am_nu_non_rot, okay
1576 : real(dp) :: &
1577 : am_nu_visc_factor, &
1578 : am_nu_DSI_factor, &
1579 : am_nu_SH_factor, &
1580 : am_nu_SSI_factor, &
1581 : am_nu_ES_factor, &
1582 : am_nu_GSF_factor, &
1583 : am_nu_ST_factor, &
1584 : f, lgT, full_off, full_on
1585 : real(dp), dimension(:), allocatable :: & ! work vectors for tridiagonal solve
1586 0 : sig, rhs, d, du, dl, bp, vp, xp, x
1587 :
1588 : include 'formats'
1589 :
1590 0 : ierr = 0
1591 0 : nz = s% nz
1592 :
1593 0 : call set_rotation_mixing_info(s, ierr)
1594 0 : if (ierr /= 0) then
1595 0 : if (s% report_ierr) &
1596 0 : write(*,*) 'update_rotation_mixing_info failed in call to set_rotation_mixing_info'
1597 0 : return
1598 : end if
1599 :
1600 0 : call check('after set_rotation_mixing_info')
1601 0 : if (s% D_omega_flag) call check_D_omega('check_D_omega after set_rotation_mixing_info')
1602 :
1603 : ! include rotation part for mixing abundances
1604 0 : full_on = s% D_mix_rotation_max_logT_full_on
1605 0 : full_off = s% D_mix_rotation_min_logT_full_off
1606 :
1607 0 : do k = 2, nz
1608 0 : lgT = s% lnT(k)/ln10
1609 0 : if (lgT <= full_on) then
1610 : f = 1d0
1611 0 : else if (lgT >= full_off) then
1612 : f = 0d0
1613 : else ! lgT > full_on and < full_off
1614 0 : f = (full_off - lgT) / (full_off - full_on)
1615 : end if
1616 :
1617 0 : if (s% D_omega_flag) then
1618 0 : s% D_mix_rotation(k) = f*s% am_D_mix_factor*s% D_omega(k)
1619 : else
1620 : s% D_mix_rotation(k) = &
1621 : f*s% am_D_mix_factor *( &
1622 : ! note: have dropped viscosity from mixing
1623 : s% D_DSI_factor * s% D_DSI(k) + &
1624 : s% D_SH_factor * s% D_SH(k) + &
1625 : s% D_SSI_factor * s% D_SSI(k) + &
1626 : s% D_ES_factor * s% D_ES(k) + &
1627 : s% D_GSF_factor * s% D_GSF(k) + &
1628 0 : s% D_ST_factor * s% D_ST(k))
1629 : end if
1630 0 : s% D_mix(k) = s% D_mix_non_rotation(k) + s% D_mix_rotation(k)
1631 : end do
1632 :
1633 0 : call check('after include rotation part for mixing abundances')
1634 :
1635 0 : am_nu_DSI_factor = s% am_nu_DSI_factor
1636 0 : am_nu_SH_factor = s% am_nu_SH_factor
1637 0 : am_nu_SSI_factor = s% am_nu_SSI_factor
1638 0 : am_nu_ES_factor = s% am_nu_ES_factor
1639 0 : am_nu_GSF_factor = s% am_nu_GSF_factor
1640 0 : am_nu_ST_factor = s% am_nu_ST_factor
1641 0 : am_nu_visc_factor = s% am_nu_visc_factor
1642 :
1643 0 : if ((.not. s% am_nu_rot_flag) .and. &
1644 : (s% D_omega_flag .and. .not. s% job% use_D_omega_for_am_nu_rot)) then
1645 : ! check for any am_nu factors > 0 and not same as for D_omega
1646 0 : okay = .true.
1647 0 : if (am_nu_DSI_factor >= 0 .and. am_nu_DSI_factor /= s% D_DSI_factor) okay = .false.
1648 0 : if (am_nu_SH_factor >= 0 .and. am_nu_SH_factor /= s% D_SH_factor) okay = .false.
1649 0 : if (am_nu_SSI_factor >= 0 .and. am_nu_SSI_factor /= s% D_SSI_factor) okay = .false.
1650 0 : if (am_nu_DSI_factor >= 0 .and. am_nu_DSI_factor /= s% D_DSI_factor) okay = .false.
1651 0 : if (am_nu_ES_factor >= 0 .and. am_nu_ES_factor /= s% D_ES_factor) okay = .false.
1652 0 : if (am_nu_GSF_factor >= 0 .and. am_nu_GSF_factor /= s% D_GSF_factor) okay = .false.
1653 0 : if (am_nu_ST_factor >= 0 .and. am_nu_ST_factor /= s% D_ST_factor) okay = .false.
1654 0 : if (.not. okay) then
1655 0 : write(*,*) 'have an am_nu factor >= 0 and not same as corresponding D_omega factor'
1656 0 : write(*,*) 'so if want smoothing like D_omega, must set am_nu_rot_flag true'
1657 0 : write(*,*) 'please fix this'
1658 0 : ierr = -1
1659 0 : return
1660 : end if
1661 : end if
1662 :
1663 : ! If am_nu_..._factor < -1, use the D_..._factor
1664 0 : if (am_nu_DSI_factor < 0) am_nu_DSI_factor = s% D_DSI_factor
1665 0 : if (am_nu_SH_factor < 0) am_nu_SH_factor = s% D_SH_factor
1666 0 : if (am_nu_SSI_factor < 0) am_nu_SSI_factor = s% D_SSI_factor
1667 0 : if (am_nu_ES_factor < 0) am_nu_ES_factor = s% D_ES_factor
1668 0 : if (am_nu_GSF_factor < 0) am_nu_GSF_factor = s% D_GSF_factor
1669 0 : if (am_nu_ST_factor < 0) am_nu_ST_factor = s% D_ST_factor
1670 0 : if (am_nu_visc_factor < 0) am_nu_visc_factor = s% D_visc_factor
1671 :
1672 : ! set set_min_am_nu_non_rot to s% set_min_am_nu_non_rot if
1673 : ! ye > min_center_ye_for_min_am_nu_non_rot
1674 : ! and s% set_min_am_nu_non_rot
1675 : set_min_am_nu_non_rot = &
1676 : s% set_min_am_nu_non_rot .and. &
1677 : s% ye(nz) >= s% min_center_Ye_for_min_am_nu_non_rot .and. &
1678 0 : (.not. s% set_uniform_am_nu_non_rot)
1679 :
1680 0 : if (s% am_nu_rot_flag) then
1681 0 : call set_am_nu_rot(ierr)
1682 0 : if (ierr /= 0) return
1683 : end if
1684 :
1685 0 : do k=1,nz
1686 0 : if (s% set_uniform_am_nu_non_rot) then
1687 0 : s% am_nu_non_rot(k) = s% uniform_am_nu_non_rot
1688 : else
1689 : s% am_nu_non_rot(k) = s% am_nu_factor* &
1690 0 : s% am_nu_non_rotation_factor*s% D_mix_non_rotation(k)
1691 : end if
1692 0 : if (set_min_am_nu_non_rot) &
1693 0 : s% am_nu_non_rot(k) = max(s% min_am_nu_non_rot, s% am_nu_non_rot(k))
1694 0 : if (s% am_nu_rot_flag) then
1695 : ! already have am_nu_rot(k) from calling set_am_nu_rot above
1696 0 : else if (s% D_omega_flag .and. s% job% use_D_omega_for_am_nu_rot) then
1697 0 : s% am_nu_rot(k) = s% am_nu_factor*s% D_omega(k)
1698 : else
1699 : s% am_nu_rot(k) = s% am_nu_factor * ( &
1700 : am_nu_visc_factor*s% D_visc(k) + &
1701 : am_nu_DSI_factor*s% D_DSI(k) + &
1702 : am_nu_SH_factor*s% D_SH(k) + &
1703 : am_nu_SSI_factor*s% D_SSI(k) + &
1704 : am_nu_ES_factor*s% D_ES(k) + &
1705 : am_nu_GSF_factor*s% D_GSF(k) + &
1706 0 : am_nu_ST_factor*s% nu_ST(k))
1707 : end if
1708 0 : if (s% am_nu_rot(k) < 0d0) s% am_nu_rot(k) = 0d0
1709 : s% am_nu_omega(k) = &
1710 : s% am_nu_omega_non_rot_factor*s% am_nu_non_rot(k) + &
1711 0 : s% am_nu_omega_rot_factor*s% am_nu_rot(k)
1712 0 : if (s% am_nu_omega(k) < 0) then
1713 0 : ierr = -1
1714 0 : if (s% report_ierr) write(*,3) &
1715 0 : 'update_rotation_mixing_info: am_nu_omega(k) < 0', k, nz, s% am_nu_omega(k), &
1716 0 : s% am_nu_non_rot(k), s% am_nu_rot(k), s% D_omega(k)
1717 0 : return
1718 : end if
1719 : s% am_nu_j(k) = &
1720 : s% am_nu_j_non_rot_factor*s% am_nu_non_rot(k) + &
1721 0 : s% am_nu_j_rot_factor*s% am_nu_rot(k)
1722 0 : if (s% am_nu_j(k) < 0) then
1723 0 : ierr = -1
1724 0 : if (s% report_ierr) write(*,3) &
1725 0 : 'update_rotation_mixing_info: am_nu_j(k) < 0', k, nz, s% am_nu_j(k), &
1726 0 : s% am_nu_non_rot(k), s% am_nu_rot(k), s% D_omega(k)
1727 0 : return
1728 : end if
1729 : end do
1730 :
1731 0 : if (s% use_other_am_mixing) then
1732 0 : call s% other_am_mixing(s% id, ierr)
1733 0 : if (ierr /= 0) return
1734 : end if
1735 :
1736 : contains
1737 :
1738 0 : subroutine check(str)
1739 : character(len=*) :: str
1740 : integer :: k
1741 : include 'formats'
1742 0 : do k = 2, s% nz
1743 0 : if (is_bad_num(s% D_mix(k))) then
1744 0 : ierr = -1
1745 0 : if (s% report_ierr) then
1746 0 : write(*,3) trim(str) // ' mixing_type, D_mix', k, s% mixing_type(k), s% D_mix(k)
1747 0 : if (s% rotation_flag) then
1748 0 : if (is_bad_num(s% D_mix_non_rotation(k))) &
1749 0 : write(*,2) 's% D_mix_non_rotation(k)', k, s% D_mix_non_rotation(k)
1750 0 : if (is_bad_num(s% D_visc(k))) write(*,2) 's% D_visc(k)', k, s% D_visc(k)
1751 0 : if (is_bad_num(s% D_DSI(k))) write(*,2) 's% D_DSI(k)', k, s% D_DSI(k)
1752 0 : if (is_bad_num(s% D_SH(k))) write(*,2) 's% D_SH(k)', k, s% D_SH(k)
1753 0 : if (is_bad_num(s% D_SSI(k))) write(*,2) 's% D_SSI(k)', k, s% D_SSI(k)
1754 0 : if (is_bad_num(s% D_ES(k))) write(*,2) 's% D_ES(k)', k, s% D_ES(k)
1755 0 : if (is_bad_num(s% D_GSF(k))) write(*,2) 's% D_GSF(k)', k, s% D_GSF(k)
1756 0 : if (is_bad_num(s% D_ST(k))) write(*,2) 's% D_ST(k)', k, s% D_ST(k)
1757 : end if
1758 : end if
1759 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'set mixing info')
1760 : end if
1761 : end do
1762 0 : end subroutine check
1763 :
1764 0 : subroutine check_D_omega(str)
1765 : character(len=*) :: str
1766 : integer :: k
1767 : include 'formats'
1768 0 : do k = 2, s% nz
1769 0 : if (is_bad_num(s% D_omega(k))) then
1770 0 : ierr = -1
1771 0 : if (s% report_ierr) then
1772 0 : write(*,3) trim(str) // ' mixing_type, D_omega', k, s% mixing_type(k), s% D_omega(k)
1773 0 : write(*,*) 's% doing_finish_load_model', s% doing_finish_load_model
1774 0 : if (s% rotation_flag) then
1775 0 : if (is_bad_num(s% D_mix_non_rotation(k))) &
1776 0 : write(*,2) 's% D_mix_non_rotation(k)', k, s% D_mix_non_rotation(k)
1777 0 : if (is_bad_num(s% D_visc(k))) write(*,2) 's% D_visc(k)', k, s% D_visc(k)
1778 0 : if (is_bad_num(s% D_DSI(k))) write(*,2) 's% D_DSI(k)', k, s% D_DSI(k)
1779 0 : if (is_bad_num(s% D_SH(k))) write(*,2) 's% D_SH(k)', k, s% D_SH(k)
1780 0 : if (is_bad_num(s% D_SSI(k))) write(*,2) 's% D_SSI(k)', k, s% D_SSI(k)
1781 0 : if (is_bad_num(s% D_ES(k))) write(*,2) 's% D_ES(k)', k, s% D_ES(k)
1782 0 : if (is_bad_num(s% D_GSF(k))) write(*,2) 's% D_GSF(k)', k, s% D_GSF(k)
1783 0 : if (is_bad_num(s% D_ST(k))) write(*,2) 's% D_ST(k)', k, s% D_ST(k)
1784 : end if
1785 : end if
1786 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'set mixing info')
1787 : end if
1788 : end do
1789 0 : end subroutine check_D_omega
1790 :
1791 0 : subroutine set_am_nu_rot(ierr)
1792 : use alloc
1793 : use rotation_mix_info, only: smooth_for_rotation
1794 : integer, intent(out) :: ierr
1795 : integer :: i, k, nz
1796 : real(dp) :: &
1797 : dt, rate, d_ddt_dm1, d_ddt_d00, d_ddt_dp1, m, &
1798 : d_dt, d_dt_in, d_dt_out, am_nu_rot_source
1799 : include 'formats'
1800 :
1801 0 : ierr = 0
1802 0 : nz = s% nz
1803 0 : dt = s% dt
1804 :
1805 0 : if (s% am_nu_rot_flag .and. s% doing_finish_load_model) then
1806 0 : do k=1,nz
1807 0 : s% am_nu_rot(k) = 0d0
1808 : end do
1809 0 : else if (s% am_nu_rot_flag) then
1810 :
1811 0 : do k=1,nz
1812 0 : if (s% q(k) <= s% max_q_for_nu_omega_zero_in_convection_region .and. &
1813 : s% mixing_type(k) == convective_mixing) then
1814 0 : s% am_nu_rot(k) = 0d0
1815 0 : cycle
1816 : end if
1817 : am_nu_rot_source = s% am_nu_factor * ( &
1818 : am_nu_visc_factor*s% D_visc(k) + &
1819 : am_nu_DSI_factor*s% D_DSI(k) + &
1820 : am_nu_SH_factor*s% D_SH(k) + &
1821 : am_nu_SSI_factor*s% D_SSI(k) + &
1822 : am_nu_ES_factor*s% D_ES(k) + &
1823 : am_nu_GSF_factor*s% D_GSF(k) + &
1824 0 : am_nu_ST_factor*s% nu_ST(k))
1825 0 : if (is_bad(am_nu_rot_source)) then
1826 0 : write(*,2) 'am_nu_rot_source', k, am_nu_rot_source
1827 0 : call mesa_error(__FILE__,__LINE__,'set am_nu_rot')
1828 : end if
1829 0 : s% am_nu_rot(k) = am_nu_rot_source
1830 0 : if (is_bad(s% am_nu_rot(k))) then
1831 0 : write(*,2) 's% am_nu_rot(k)', k, s% am_nu_rot(k)
1832 0 : write(*,2) 'am_nu_rot_source', k, am_nu_rot_source
1833 0 : call mesa_error(__FILE__,__LINE__,'set am_nu_rot')
1834 : end if
1835 : end do
1836 :
1837 0 : if (s% smooth_am_nu_rot > 0 .or. &
1838 : (s% nu_omega_mixing_rate > 0d0 .and. s% dt > 0)) then
1839 :
1840 0 : allocate(sig(nz), rhs(nz), d(nz), du(nz), dl(nz), bp(nz), vp(nz), xp(nz), x(nz))
1841 :
1842 0 : if (s% smooth_am_nu_rot > 0) then
1843 0 : call smooth_for_rotation(s, s% am_nu_rot, s% smooth_am_nu_rot, sig)
1844 : end if
1845 :
1846 0 : if (s% nu_omega_mixing_rate > 0d0 .and. s% dt > 0) then ! mix am_nu_rot
1847 :
1848 0 : rate = min(s% nu_omega_mixing_rate, 1d0/dt)
1849 0 : do k=2,nz-1
1850 0 : if (s% am_nu_rot(k) == 0 .or. s% am_nu_rot(k+1) == 0) then
1851 0 : sig(k) = 0
1852 : else if ((.not. s% nu_omega_mixing_across_convection_boundary) .and. &
1853 0 : s% mixing_type(k) /= convective_mixing .and. &
1854 : (s% mixing_type(k-1) == convective_mixing .or. &
1855 : s% mixing_type(k+1) == convective_mixing)) then
1856 0 : sig(k) = 0
1857 : else
1858 0 : sig(k) = rate*dt
1859 : end if
1860 : end do
1861 0 : sig(1) = 0
1862 0 : sig(nz) = 0
1863 :
1864 0 : do k=1,nz
1865 0 : if (k < nz) then
1866 0 : d_dt_in = sig(k)*(s% am_nu_rot(k+1) - s% am_nu_rot(k))
1867 : else
1868 0 : d_dt_in = -sig(k)*s% am_nu_rot(k)
1869 : end if
1870 0 : if (k > 1) then
1871 0 : d_dt_out = sig(k-1)*(s% am_nu_rot(k) - s% am_nu_rot(k-1))
1872 0 : d_ddt_dm1 = sig(k-1)
1873 0 : d_ddt_d00 = -(sig(k-1) + sig(k))
1874 : else
1875 0 : d_dt_out = 0
1876 0 : d_ddt_dm1 = 0
1877 0 : d_ddt_d00 = -sig(k)
1878 : end if
1879 0 : d_dt = d_dt_in - d_dt_out
1880 0 : d_ddt_dp1 = sig(k)
1881 0 : rhs(k) = d_dt
1882 0 : d(k) = 1d0 - d_ddt_d00
1883 0 : if (k < nz) then
1884 0 : du(k) = -d_ddt_dp1
1885 : else
1886 0 : du(k) = 0
1887 : end if
1888 0 : if (k > 1) dl(k-1) = -d_ddt_dm1
1889 : end do
1890 0 : dl(nz) = 0
1891 :
1892 : ! solve tridiagonal
1893 0 : bp(1) = d(1)
1894 0 : vp(1) = rhs(1)
1895 0 : do i = 2,nz
1896 0 : if (bp(i-1) == 0) then
1897 0 : write(*,*) 'failed in set_am_nu_rot', s% model_number
1898 0 : call mesa_error(__FILE__,__LINE__,'mix_am_nu_rot')
1899 : ierr = -1
1900 0 : return
1901 : end if
1902 0 : m = dl(i-1)/bp(i-1)
1903 0 : bp(i) = d(i) - m*du(i-1)
1904 0 : vp(i) = rhs(i) - m*vp(i-1)
1905 : end do
1906 0 : xp(nz) = vp(nz)/bp(nz)
1907 0 : x(nz) = xp(nz)
1908 0 : do i = nz-1, 1, -1
1909 0 : xp(i) = (vp(i) - du(i)*xp(i+1))/bp(i)
1910 0 : x(i) = xp(i)
1911 : end do
1912 :
1913 0 : do k=2,nz
1914 0 : if (is_bad(x(k))) then
1915 : return
1916 : write(*,3) 's% am_nu_rot(k) prev, x', k, &
1917 : s% model_number, s% am_nu_rot(k), x(k), bp(i)
1918 : call mesa_error(__FILE__,__LINE__,'mix_am_nu_rot')
1919 : end if
1920 : end do
1921 :
1922 : ! update am_nu_rot
1923 0 : do k=2,nz
1924 0 : s% am_nu_rot(k) = s% am_nu_rot(k) + x(k)
1925 0 : if (is_bad(s% am_nu_rot(k))) then
1926 0 : write(*,3) 's% am_nu_rot(k)', k, s% model_number, s% am_nu_rot(k)
1927 0 : call mesa_error(__FILE__,__LINE__,'mix_am_nu_rot')
1928 : end if
1929 0 : if (s% am_nu_rot(k) < 0d0) s% am_nu_rot(k) = 0d0
1930 : end do
1931 0 : s% am_nu_rot(1) = 0d0
1932 :
1933 : end if
1934 :
1935 : end if
1936 :
1937 : end if
1938 :
1939 0 : if (s% am_nu_rot_flag) then ! check
1940 0 : do k=1,nz
1941 0 : if (is_bad(s% am_nu_rot(k))) then
1942 0 : write(*,2) 'before return s% am_nu_rot(k)', k, s% am_nu_rot(k)
1943 0 : call mesa_error(__FILE__,__LINE__,'set_am_nu_rot')
1944 : end if
1945 0 : if (s% am_nu_rot(k) < 0d0) s% am_nu_rot(k) = 0d0
1946 : end do
1947 : end if
1948 :
1949 : end subroutine set_am_nu_rot
1950 :
1951 : end subroutine update_rotation_mixing_info
1952 :
1953 :
1954 0 : subroutine set_RTI_mixing_info(s, ierr)
1955 : use chem_def, only: ih1
1956 : use star_utils, only: get_shock_info
1957 : type (star_info), pointer :: s
1958 : integer, intent(out) :: ierr
1959 : real(dp) :: &
1960 : C, alpha_face, f, v, &
1961 : min_dm, alfa, cs, r, shock_mass_start, &
1962 : log_max_boost, m_full_boost, m_no_boost, max_boost, &
1963 : dm_for_center_eta_nondecreasing, min_eta
1964 : integer :: k, nz, i_h1
1965 : include 'formats'
1966 0 : ierr = 0
1967 0 : if (.not. s% RTI_flag) return
1968 :
1969 0 : nz = s% nz
1970 :
1971 0 : s% eta_RTI(1:nz) = 0d0
1972 0 : s% etamid_RTI(1:nz) = 0d0
1973 0 : s% boost_for_eta_RTI(1:nz) = 0d0
1974 0 : s% sig_RTI(1:nz) = 0d0
1975 0 : s% sigmid_RTI(1:nz) = 0d0
1976 :
1977 0 : if (s% RTI_C <= 0d0) return
1978 :
1979 0 : i_h1 = s% net_iso(ih1)
1980 :
1981 0 : shock_mass_start = 1d99
1982 0 : do k = 1, nz
1983 0 : if (s% u_flag) then
1984 0 : v = s% u(k)
1985 : else
1986 0 : v = s% v(k)
1987 : end if
1988 0 : if (v > s% csound(k)) then
1989 0 : if (k > 1) shock_mass_start = s% m(k) ! skip this after breakout
1990 : exit
1991 : end if
1992 : end do
1993 :
1994 0 : min_dm = s% RTI_min_dm_behind_shock_for_full_on*Msun
1995 0 : log_max_boost = s% RTI_log_max_boost
1996 0 : max_boost = exp10(log_max_boost)
1997 0 : m_full_boost = s% RTI_m_full_boost*Msun
1998 0 : m_no_boost = s% RTI_m_no_boost*Msun
1999 :
2000 0 : min_eta = -1d0
2001 0 : dm_for_center_eta_nondecreasing = Msun*s% RTI_dm_for_center_eta_nondecreasing
2002 :
2003 0 : do k=1,nz
2004 :
2005 0 : f = max(0d0, s% X(k) - s% RTI_C_X0_frac*s% surface_h1)
2006 0 : C = s% RTI_C*(1d0 + f*f*s% RTI_C_X_factor)
2007 :
2008 0 : if (s% m(k) < m_no_boost) then
2009 0 : if (s% m(k) <= m_full_boost) then
2010 0 : C = C*max_boost
2011 : else
2012 0 : alfa = (m_no_boost - s% m(k))/(m_no_boost - m_full_boost)
2013 0 : C = C*exp10(alfa*log_max_boost)
2014 : end if
2015 : end if
2016 :
2017 0 : if (shock_mass_start - s% m(k) < min_dm) &
2018 0 : C = C*(shock_mass_start - s% m(k))/min_dm
2019 :
2020 0 : if (k > 1) then
2021 : alpha_face = &
2022 : (s% dq(k-1)*s% alpha_RTI(k) + s% dq(k)*s% alpha_RTI(k-1))/ &
2023 0 : (s% dq(k-1) + s% dq(k))
2024 0 : cs = s% csound_face(k)
2025 0 : r = s% r(k)
2026 0 : s% eta_RTI(k) = C*alpha_face*cs*r
2027 :
2028 0 : if (is_bad(s% eta_RTI(k)) .and. s% q(k) <= s% alpha_RTI_src_max_q) then
2029 0 : ierr = -1
2030 0 : return
2031 : if (s% stop_for_bad_nums) then
2032 : write(*,2) 's% eta_RTI(k)', k, s% eta_RTI(k)
2033 : call mesa_error(__FILE__,__LINE__,'set_RTI_mixing_info')
2034 : end if
2035 : end if
2036 :
2037 0 : if (s% m(k) - s% M_center <= dm_for_center_eta_nondecreasing) then
2038 0 : if (min_eta < 0d0) then
2039 : min_eta = s% eta_RTI(k)
2040 0 : else if (min_eta > s% eta_RTI(k)) then
2041 0 : s% eta_RTI(k) = min_eta
2042 : end if
2043 : end if
2044 :
2045 : end if
2046 :
2047 0 : s% etamid_RTI(k) = max(min_eta, C*s% alpha_RTI(k)*s% csound(k)*s% rmid(k))
2048 0 : s% boost_for_eta_RTI(k) = C/s% RTI_C
2049 :
2050 0 : if (is_bad(s% etamid_RTI(k))) then
2051 0 : ierr = -1
2052 0 : return
2053 : end if
2054 :
2055 : end do
2056 :
2057 0 : call get_shock_info(s)
2058 :
2059 : ! sig_RTI(k) is mixing flow across face k in (gm sec^1)
2060 : call get_RTI_sigmas(s, s% sig_RTI, s% eta_RTI, &
2061 0 : s% rho_face, s% r, s% dm_bar, s% dt, ierr)
2062 0 : if (ierr /= 0) return
2063 :
2064 0 : if (s% v_flag) then
2065 : ! sigmid_RTI(k) is mixing flow at center k in (gm sec^1)
2066 : call get_RTI_sigmas(s, s% sigmid_RTI, s% etamid_RTI, &
2067 0 : s% rho, s% rmid, s% dm, s% dt, ierr)
2068 0 : if (ierr /= 0) return
2069 : end if
2070 :
2071 : end subroutine set_RTI_mixing_info
2072 :
2073 :
2074 0 : subroutine get_RTI_sigmas(s, sig, eta, rho, r, dm_bar, dt, ierr)
2075 : type (star_info), pointer :: s
2076 : real(dp), dimension(:) :: sig, eta, rho, r, dm_bar
2077 : real(dp), intent(in) :: dt
2078 : integer, intent(out) :: ierr
2079 :
2080 : integer :: nz, k
2081 : real(dp) :: D, f1, f, cdcterm, max_sig, &
2082 : siglim, sig_term_limit, xmstar
2083 :
2084 : include 'formats'
2085 :
2086 0 : ierr = 0
2087 0 : nz = s% nz
2088 0 : xmstar = s% xmstar
2089 0 : sig_term_limit = s% sig_term_limit
2090 :
2091 0 : sig(1) = 0d0
2092 0 : max_sig = 0d0
2093 0 : do k = 2, nz
2094 0 : D = eta(k)
2095 0 : if (D <= 0 .or. &
2096 : (s% q(k) > s% shock_q .and. s% shock_q > 0d0)) then
2097 0 : sig(k) = 0d0
2098 0 : cycle
2099 : end if
2100 0 : f1 = pi4*r(k)*r(k)*rho(k)
2101 0 : f = f1*f1
2102 0 : cdcterm = s% mix_factor*D*f
2103 0 : sig(k) = cdcterm/dm_bar(k)
2104 0 : if (is_bad(sig(k))) then
2105 0 : if (s% stop_for_bad_nums) then
2106 0 : write(*,2) 'sig(k)', k, sig(k)
2107 0 : call mesa_error(__FILE__,__LINE__,'get_RTI_sigmas')
2108 : end if
2109 0 : sig(k) = 1d99
2110 : end if
2111 0 : if (sig(k) > max_sig) max_sig = sig(k)
2112 : end do
2113 :
2114 0 : if (max_sig < 1) return
2115 0 : do k = 1, nz
2116 0 : if (sig(k) == 0) cycle
2117 0 : siglim = sig_term_limit*dm_bar(k)/dt
2118 0 : if (sig(k) > siglim) sig(k) = siglim
2119 : end do
2120 :
2121 : end subroutine get_RTI_sigmas
2122 :
2123 :
2124 0 : subroutine set_dPdr_dRhodr_info(s, ierr)
2125 : type (star_info), pointer :: s
2126 : integer, intent(out) :: ierr
2127 : real(dp) :: rho, r00, alfa00, beta00, &
2128 : dr_m1, dr_00, c, d, am1, a00, ap1, v, rmid
2129 0 : real(dp), allocatable, dimension(:) :: dPdr, drhodr, P_face, rho_face
2130 : integer :: k, nz
2131 : logical, parameter :: do_slope_limiting = .false.
2132 : include 'formats'
2133 0 : ierr = 0
2134 0 : if (.not. s% RTI_flag) return
2135 0 : if (s% dPdr_dRhodr_info(1) >= 0d0) then
2136 0 : if (is_bad(s% dPdr_dRhodr_info(1))) call mesa_error(__FILE__,__LINE__,'set_dPdr_dRhodr_info')
2137 : return ! already set this step
2138 : end if
2139 :
2140 0 : nz = s% nz
2141 :
2142 0 : allocate(dPdr(nz), drhodr(nz), P_face(nz), rho_face(nz))
2143 :
2144 0 : do k=2,nz
2145 0 : rho = s% rho(k)
2146 0 : r00 = s% r(k)
2147 0 : alfa00 = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
2148 0 : beta00 = 1d0 - alfa00
2149 0 : P_face(k) = alfa00*s% Peos(k) + beta00*s% Peos(k-1)
2150 0 : rho_face(k) = alfa00*s% rho(k) + beta00*s% rho(k-1)
2151 : end do
2152 : P_face(1) = s% Peos(1)
2153 0 : rho_face(1) = s% rho(1)
2154 :
2155 0 : do k=1,nz
2156 0 : if (s% u_flag) then
2157 0 : v = s% u(k)
2158 : else
2159 0 : v = s% v(k)
2160 : end if
2161 0 : if (k == nz .or. k == 1 .or. &
2162 : (s% alpha_RTI_start(k) == 0d0 .and. &
2163 0 : v/s% csound(k) < s% alpha_RTI_src_min_v_div_cs)) then
2164 0 : drhodr(k) = 0d0
2165 0 : dPdr(k) = 0d0
2166 : else if (do_slope_limiting) then
2167 : dr_m1 = s% r(k-1) - s% r(k)
2168 : dr_00 = s% r(k) - s% r(k+1)
2169 : dPdr(k) = slope_limit(P_face, k, dr_m1, dr_00)
2170 : drhodr(k) = slope_limit(rho_face, k, dr_m1, dr_00)
2171 : else
2172 : !dr_00 = s% r(k) - s% r(k+1)
2173 0 : rmid = 0.5d0*(s% r(k) + s% r(k+1))
2174 0 : dr_00 = s% dm(k)/(pi4*rmid*rmid*s% rho(k)) ! don't subtract r's to get dr
2175 0 : dPdr(k) = (P_face(k) - P_face(k+1))/dr_00
2176 0 : drhodr(k) = (rho_face(k) - rho_face(k+1))/dr_00
2177 : end if
2178 : end do
2179 :
2180 : ! smooth dPdr and drhodr
2181 0 : if (s% RTI_smooth_mass > 0d0 .and. s% RTI_smooth_iterations > 0) then
2182 : call do_smoothing_by_mass( &
2183 0 : s, s% RTI_smooth_mass, s% RTI_smooth_iterations, drhodr, ierr)
2184 0 : if (ierr /= 0) return
2185 : call do_smoothing_by_mass( &
2186 0 : s, s% RTI_smooth_mass, s% RTI_smooth_iterations, dPdr, ierr)
2187 0 : if (ierr /= 0) return
2188 0 : else if (s% RTI_smooth_fraction < 1d0) then
2189 0 : c = s% RTI_smooth_fraction
2190 0 : d = 0.5d0*(1d0 - c)
2191 0 : am1 = 0
2192 0 : a00 = dPdr(1)
2193 0 : ap1 = dPdr(2)
2194 0 : do k=2,nz-1
2195 0 : am1 = a00
2196 0 : a00 = ap1
2197 0 : ap1 = dPdr(k+1)
2198 0 : dPdr(k) = c*a00 + d*(am1 + ap1)
2199 : end do
2200 0 : dPdr(nz) = c*a00 + 2*d*am1
2201 0 : a00 = drhodr(1)
2202 0 : ap1 = drhodr(2)
2203 0 : do k=2,nz-1
2204 0 : am1 = a00
2205 0 : a00 = ap1
2206 0 : ap1 = drhodr(k+1)
2207 0 : if (am1 == 0d0 .and. a00 == 0d0) cycle
2208 0 : drhodr(k) = c*a00 + d*(am1 + ap1)
2209 : end do
2210 0 : drhodr(nz) = c*a00 + 2*d*am1
2211 : end if
2212 :
2213 : ! set
2214 0 : do k=1,nz
2215 0 : s% dPdr_info(k) = dPdr(k)
2216 0 : s% dRhodr_info(k) = drhodr(k)
2217 0 : s% dPdr_dRhodr_info(k) = min(0d0,dPdr(k)*drhodr(k))
2218 : end do
2219 :
2220 : contains
2221 :
2222 : real(dp) function slope_limit(v, k, dr_m1, dr_00)
2223 : real(dp), intent(in) :: v(:), dr_m1, dr_00
2224 : integer, intent(in) :: k
2225 : real(dp) :: s_00, s_p1
2226 : s_00 = (v(k-1) - v(k))/dr_m1
2227 : s_p1 = (v(k) - v(k+1))/dr_00
2228 : if (s_00*s_p1 < 0d0) then
2229 : slope_limit = 0d0
2230 : else if(abs(s_00) > abs(s_p1)) then
2231 : slope_limit = s_p1
2232 : else
2233 : slope_limit = s_00
2234 : end if
2235 : end function slope_limit
2236 :
2237 : end subroutine set_dPdr_dRhodr_info
2238 :
2239 :
2240 0 : subroutine do_smoothing_by_mass( &
2241 0 : s, smooth_mass, number_iterations, val, ierr)
2242 : type (star_info), pointer :: s
2243 : real(dp), intent(in) :: smooth_mass
2244 : integer, intent(in) :: number_iterations
2245 : real(dp) :: val(:)
2246 : integer, intent(out) :: ierr
2247 : integer :: nz, iter, k_center, k_inner, k_outer, k
2248 : real(dp) :: mlo, mhi, mmid, smooth_m, v, dm_half, mtotal, mass_from_cell
2249 0 : real(dp), allocatable :: work(:)
2250 : include 'formats'
2251 0 : ierr = 0
2252 0 : nz = s% nz
2253 0 : allocate(work(nz))
2254 0 : do k = 1, nz
2255 0 : work(k) = val(k)
2256 : end do
2257 0 : smooth_m = smooth_mass*Msun
2258 0 : dm_half = 0.5d0*smooth_m
2259 0 : do iter = 1, number_iterations
2260 0 : do k_center = nz-1, 2, -1
2261 0 : mmid = s% m(k_center) - 0.5d0*s% dm(k_center)
2262 0 : mlo = mmid - dm_half
2263 0 : mhi = mmid + dm_half
2264 0 : k_inner = k_center
2265 0 : v = 0d0
2266 0 : mtotal = 0d0
2267 0 : do k=k_center, nz, -1
2268 0 : if (s% m(k) - 0.5d0*s% dm(k) <= mlo) then
2269 0 : k_inner = k
2270 0 : mass_from_cell = s% m(k) - mlo
2271 0 : v = v + val(k)*mass_from_cell
2272 0 : mtotal = mtotal + mass_from_cell
2273 0 : exit
2274 : end if
2275 : end do
2276 0 : k_outer = k_center
2277 0 : do k=k_center,1,-1
2278 0 : if (s% m(k) >= mhi) then
2279 0 : k_outer = k
2280 0 : mass_from_cell = mhi - (s% m(k) - s% dm(k))
2281 0 : v = v + val(k)*mass_from_cell
2282 0 : mtotal = mtotal + mass_from_cell
2283 0 : exit
2284 : end if
2285 : end do
2286 0 : do k=k_outer+1,k_inner-1
2287 0 : v = v + val(k)*s% dm(k)
2288 0 : mtotal = mtotal + s% dm(k)
2289 : end do
2290 0 : if(mtotal>0d0) then
2291 0 : work(k_center) = v/mtotal
2292 : else
2293 0 : work(k_center) = 0d0
2294 : end if
2295 : end do
2296 0 : do k = 2, nz-1
2297 0 : val(k) = work(k)
2298 : end do
2299 : end do
2300 0 : end subroutine do_smoothing_by_mass
2301 :
2302 :
2303 44 : subroutine set_dxdt_mix(s)
2304 : type (star_info), pointer :: s
2305 : real(dp) :: x00, xp1, xm1, dx00, dxp1, dm, sig00, sigp1, &
2306 : flux00, dflux00_dxm1, dflux00_dx00, &
2307 : fluxp1, dfluxp1_dx00, dfluxp1_dxp1
2308 : integer :: j, k
2309 : include 'formats'
2310 :
2311 52392 : do k = 1, s% nz
2312 :
2313 52348 : dm = s% dm(k)
2314 52348 : sig00 = s% sig(k)
2315 :
2316 52348 : if (k < s% nz) then
2317 52304 : sigp1 = s% sig(k+1)
2318 : else
2319 : sigp1 = 0
2320 : end if
2321 :
2322 52348 : if (k > 1) then
2323 52304 : dflux00_dxm1 = -sig00
2324 52304 : dflux00_dx00 = sig00
2325 : else
2326 : dflux00_dxm1 = 0
2327 : dflux00_dx00 = 0
2328 : end if
2329 :
2330 52348 : if (k < s% nz) then
2331 52304 : dfluxp1_dx00 = -sigp1
2332 52304 : dfluxp1_dxp1 = sigp1
2333 : else
2334 : dfluxp1_dx00 = 0
2335 : dfluxp1_dxp1 = 0
2336 : end if
2337 :
2338 52348 : s% d_dxdt_mix_dx00(k) = (dfluxp1_dx00 - dflux00_dx00)/dm
2339 52348 : s% d_dxdt_mix_dxm1(k) = -dflux00_dxm1/dm
2340 52348 : s% d_dxdt_mix_dxp1(k) = dfluxp1_dxp1/dm
2341 :
2342 :
2343 471176 : do j=1, s% species
2344 418784 : x00 = s% xa(j,k)
2345 418784 : if (k > 1) then
2346 418432 : xm1 = s% xa(j,k-1)
2347 418432 : dx00 = xm1 - x00
2348 418432 : flux00 = -sig00*dx00
2349 : else
2350 : flux00 = 0
2351 : end if
2352 418784 : if (k < s% nz) then
2353 418432 : xp1 = s% xa(j,k+1)
2354 418432 : dxp1 = x00 - xp1
2355 418432 : fluxp1 = -sigp1*dxp1
2356 : else
2357 : fluxp1 = 0
2358 : end if
2359 471132 : s% dxdt_mix(j,k) = (fluxp1 - flux00)/dm
2360 : end do
2361 :
2362 : end do
2363 :
2364 44 : end subroutine set_dxdt_mix
2365 :
2366 :
2367 :
2368 33 : subroutine add_RTI_turbulence(s, ierr)
2369 : type (star_info), pointer :: s
2370 : integer, intent(out) :: ierr
2371 :
2372 : integer :: nz, k
2373 : real(dp) :: coeff, D, cdc, alfa, beta, rho_face, lg_coeff, full_off, full_on, min_m
2374 : logical :: dbg
2375 :
2376 : include 'formats'
2377 :
2378 33 : dbg = .false.
2379 :
2380 33 : ierr = 0
2381 33 : nz = s% nz
2382 :
2383 33 : if (.not. s% RTI_flag) return
2384 :
2385 0 : coeff = s% composition_RTI_diffusion_factor
2386 0 : if (coeff <= 0) return
2387 0 : lg_coeff = log10(coeff)
2388 0 : full_off = s% min_M_RTI_factors_full_off*Msun
2389 0 : full_on = s% max_M_RTI_factors_full_on*Msun
2390 0 : min_m = s% RTI_min_m_for_D_mix_floor*Msun
2391 0 : do k = 2, nz
2392 0 : D = s% eta_RTI(k)
2393 0 : if (s% m(k) <= full_on) then
2394 0 : D = D*coeff
2395 0 : else if (s% m(k) < full_off) then
2396 : D = D*exp10( &
2397 0 : lg_coeff*(full_off - s% m(k))/(full_off - full_on))
2398 : ! else full off, i.e., coeff = 1
2399 : end if
2400 0 : if (is_bad(D)) then
2401 0 : write(*,2) 'eta_RTI', k, D
2402 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'add_RTI_turbulence')
2403 0 : ierr = -1
2404 0 : cycle
2405 : end if
2406 0 : if (D < s% RTI_D_mix_floor .and. s% m(k) >= min_m) &
2407 0 : D = s% RTI_D_mix_floor
2408 0 : if (D > s% D_mix(k)) &
2409 0 : s% mixing_type(k) = rayleigh_taylor_mixing
2410 0 : D = D + s% D_mix(k)
2411 0 : alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
2412 0 : beta = 1 - alfa
2413 0 : rho_face = alfa*s% rho(k) + beta*s% rho(k-1)
2414 0 : cdc = pow2(pi4*s% r(k)*s% r(k)*rho_face)*D ! gm^2/sec
2415 0 : s% cdc(k) = cdc
2416 0 : s% D_mix(k) = D
2417 : end do
2418 :
2419 : end subroutine add_RTI_turbulence
2420 :
2421 : end module mix_info
|
