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