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