Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2021 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 :
21 : module turb_info
22 :
23 : use star_private_def
24 : use const_def, only: dp, i8, ln10, pi4, no_mixing, convective_mixing, crystallized, phase_separation_mixing
25 : use num_lib
26 : use utils_lib
27 : use auto_diff_support
28 :
29 : implicit none
30 :
31 : private
32 : public :: set_mlt_vars ! for hydro_vars and conv_premix
33 : public :: do1_mlt_2 ! for predictive_mix
34 : public :: switch_to_radiative ! mix_info
35 : public :: check_for_redo_MLT ! for hydro_vars
36 : public :: set_gradT_excess_alpha ! for evolve
37 :
38 : contains
39 :
40 66 : subroutine set_mlt_vars(s, nzlo, nzhi, ierr)
41 : use star_utils, only: start_time, update_time
42 : type (star_info), pointer :: s
43 : integer, intent(in) :: nzlo, nzhi
44 : integer, intent(out) :: ierr
45 : integer :: k, op_err
46 : integer(i8) :: time0
47 66 : real(dp) :: total
48 : logical :: make_gradr_sticky_in_solver_iters
49 : include 'formats'
50 66 : ierr = 0
51 66 : if (s% doing_timing) call start_time(s, time0, total)
52 66 : !$OMP PARALLEL DO PRIVATE(k,op_err,make_gradr_sticky_in_solver_iters) SCHEDULE(dynamic,2)
53 : do k = nzlo, nzhi
54 : op_err = 0
55 : call do1_mlt_2(s, k, make_gradr_sticky_in_solver_iters, op_err)
56 : if (make_gradr_sticky_in_solver_iters .and. s% solver_iter > 3) then
57 : if (.not. s% fixed_gradr_for_rest_of_solver_iters(k)) then
58 : s% fixed_gradr_for_rest_of_solver_iters(k) = &
59 : (s% mlt_mixing_type(k) == no_mixing)
60 : end if
61 : end if
62 : if (op_err /= 0) ierr = op_err
63 : end do
64 : !$OMP END PARALLEL DO
65 66 : if (s% doing_timing) call update_time(s, time0, total, s% time_mlt)
66 :
67 66 : end subroutine set_mlt_vars
68 :
69 :
70 79994 : subroutine do1_mlt_2(s, k, &
71 : make_gradr_sticky_in_solver_iters, ierr, &
72 : mixing_length_alpha_in, gradL_composition_term_in)
73 : ! get convection info for point k
74 66 : use star_utils
75 : use turb_support, only: do1_mlt_eval
76 : use eos_def
77 : use auto_diff_support
78 : type (star_info), pointer :: s
79 : integer, intent(in) :: k
80 : logical, intent(out) :: make_gradr_sticky_in_solver_iters
81 : integer, intent(out) :: ierr
82 : real(dp), intent(in), optional :: &
83 : mixing_length_alpha_in, gradL_composition_term_in
84 :
85 : type(auto_diff_real_star_order1) :: gradr_factor
86 79994 : real(dp) :: f, gradL_composition_term, abs_du_div_cs, cs, mixing_length_alpha
87 79994 : real(dp), pointer :: vel(:)
88 : integer :: i, mixing_type, nz, k_T_max
89 : real(dp), parameter :: conv_vel_mach_limit = 0.9d0
90 79994 : real(dp) :: crystal_pad
91 : logical :: no_mix
92 : type(auto_diff_real_star_order1) :: &
93 : grada_face_ad, scale_height_ad, gradr_ad, rho_face_ad, &
94 : gradT_ad, Y_face_ad, mlt_vc_ad, D_ad, Gamma_ad
95 : include 'formats'
96 :
97 79994 : ierr = 0
98 79994 : nz = s% nz
99 :
100 79994 : if (k < 1 .or. k > nz) then
101 0 : write(*,3) 'bad k for do1_mlt', k, nz
102 0 : ierr = -1
103 0 : return
104 : call mesa_error(__FILE__,__LINE__)
105 : end if
106 :
107 79994 : if (present(mixing_length_alpha_in)) then
108 0 : mixing_length_alpha = mixing_length_alpha_in
109 : else
110 79994 : mixing_length_alpha = s% mixing_length_alpha
111 : end if
112 :
113 79994 : if (present(gradL_composition_term_in)) then
114 0 : gradL_composition_term = gradL_composition_term_in
115 79994 : else if (s% use_Ledoux_criterion) then
116 0 : gradL_composition_term = s% gradL_composition_term(k)
117 : else
118 79994 : gradL_composition_term = 0d0
119 : end if
120 :
121 79994 : grada_face_ad = get_grada_face(s,k)
122 79994 : scale_height_ad = get_scale_height_face(s,k)
123 79994 : gradr_ad = get_gradr_face(s,k)
124 :
125 79994 : if (s% rotation_flag .and. s% mlt_use_rotation_correction) then
126 0 : gradr_factor = s% ft_rot(k)/s% fp_rot(k)*s% gradr_factor(k)
127 : else
128 79994 : gradr_factor = s% gradr_factor(k)
129 : end if
130 79994 : if (is_bad_num(gradr_factor% val)) then
131 0 : ierr = -1
132 0 : return
133 : end if
134 79994 : gradr_ad = gradr_ad*gradr_factor
135 :
136 : ! now can call set_no_mixing if necessary
137 :
138 79994 : if (k == 1 .and. s% mlt_make_surface_no_mixing) then
139 0 : call set_no_mixing('surface_no_mixing')
140 0 : return
141 : end if
142 :
143 79994 : crystal_pad = s% min_dq * s% m(1) * 0.5d0
144 : if ((s% phase(k) > 0.5d0 .and. s% mu(k) > 1.7d0) &
145 79994 : .or. s% crystal_core_boundary_mass + crystal_pad > s% m(k)) then
146 : ! mu(k) check is so that we only evaluate this in C/O dominated material or heavier.
147 : ! Helium can return bad phase info on Skye, so we don't want it to shut off
148 : ! convection because of wrong phase information.
149 0 : call set_no_mixing('solid_no_mixing')
150 0 : s% mlt_mixing_type(k) = crystallized
151 0 : return
152 : end if
153 :
154 79994 : if (s% m(k) <= s% phase_sep_mixing_mass) then
155 : ! Treat as radiative for MLT purposes, and label as already mixed by phase separation
156 0 : call set_no_mixing('phase_separation_mixing')
157 0 : s% mlt_mixing_type(k) = phase_separation_mixing
158 0 : return
159 : end if
160 :
161 79994 : if (s% lnT_start(k)/ln10 > s% max_logT_for_mlt) then
162 0 : call set_no_mixing('max_logT')
163 0 : return
164 : end if
165 :
166 79994 : if (s% no_MLT_below_shock .and. (s%u_flag .or. s%v_flag)) then ! check for outward shock above k
167 0 : if (s% u_flag) then
168 0 : vel => s% u
169 : else
170 0 : vel => s% v
171 : end if
172 0 : do i=k-1,1,-1
173 0 : cs = s% csound(i)
174 79994 : if (vel(i+1) >= cs .and. vel(i) < cs) then
175 0 : call set_no_mixing('below_shock')
176 0 : return
177 : end if
178 : end do
179 : end if
180 :
181 79994 : if (s% csound_start(k) > 0d0 .and. (s% u_flag .or. s% v_flag)) then
182 0 : no_mix = .false.
183 0 : if (s% u_flag) then
184 0 : vel => s% u_start
185 : else
186 0 : vel => s% v_start
187 : end if
188 0 : abs_du_div_cs = 0d0
189 0 : if (vel(k)/1d5 > s% max_v_for_convection) then
190 : no_mix = .true.
191 0 : else if (s% q(k) > s% max_q_for_convection_with_hydro_on) then
192 : no_mix = .true.
193 0 : else if ((abs(vel(k))) >= &
194 : s% csound_start(k)*s% max_v_div_cs_for_convection) then
195 : no_mix = .true.
196 0 : else if (s% u_flag) then
197 0 : if (k == 1) then
198 : abs_du_div_cs = 1d99
199 0 : else if (k < nz) then
200 : abs_du_div_cs = max(abs(vel(k) - vel(k+1)), &
201 0 : abs(vel(k) - vel(k-1))) / s% csound_start(k)
202 : end if
203 0 : if (abs_du_div_cs > s% max_abs_du_div_cs_for_convection) then
204 : no_mix = .true.
205 : end if
206 : end if
207 : if (no_mix) then
208 0 : call set_no_mixing('no_mix')
209 0 : return
210 : end if
211 : end if
212 :
213 79994 : make_gradr_sticky_in_solver_iters = s% make_gradr_sticky_in_solver_iters
214 79994 : if (.not. make_gradr_sticky_in_solver_iters .and. &
215 : s% min_logT_for_make_gradr_sticky_in_solver_iters < 1d20) then
216 0 : k_T_max = maxloc(s% lnT_start(1:nz),dim=1)
217 : make_gradr_sticky_in_solver_iters = &
218 0 : (s% lnT_start(k_T_max)/ln10 >= s% min_logT_for_make_gradr_sticky_in_solver_iters)
219 : end if
220 79994 : if (make_gradr_sticky_in_solver_iters .and. s% fixed_gradr_for_rest_of_solver_iters(k)) then
221 0 : call set_no_mixing('gradr_sticky')
222 0 : return
223 : end if
224 :
225 : call do1_mlt_eval(s, k, s% MLT_option, gradL_composition_term, &
226 : gradr_ad, grada_face_ad, scale_height_ad, mixing_length_alpha, &
227 79994 : mixing_type, gradT_ad, Y_face_ad, mlt_vc_ad, D_ad, Gamma_ad, ierr)
228 79994 : if (ierr /= 0) then
229 0 : if (s% report_ierr) then
230 0 : write(*,*) 'ierr in do1_mlt_eval for k', k
231 : end if
232 0 : return
233 : end if
234 :
235 79994 : call store_results
236 :
237 79994 : if (s% mlt_gradT_fraction >= 0d0 .and. s% mlt_gradT_fraction <= 1d0) then
238 0 : f = s% mlt_gradT_fraction
239 : else
240 79994 : f = s% adjust_mlt_gradT_fraction(k)
241 : end if
242 79994 : call adjust_gradT_fraction(s, k, f)
243 :
244 239982 : if (s% mlt_mixing_type(k) == no_mixing .or. abs(s% gradr(k)) < 1d-20) then
245 68285 : s% L_conv(k) = 0d0
246 : else
247 11709 : s% L_conv(k) = s% L(k) * (1d0 - s% gradT(k)/s% gradr(k)) ! C&G 14.109
248 : end if
249 :
250 : contains
251 :
252 79994 : subroutine store_results
253 79994 : s% mlt_mixing_type(k) = mixing_type
254 :
255 79994 : s% grada_face_ad(k) = grada_face_ad
256 79994 : s% grada_face(k) = grada_face_ad%val
257 :
258 79994 : s% gradT_ad(k) = gradT_ad
259 79994 : s% gradT(k) = s% gradT_ad(k)%val
260 79994 : s% mlt_gradT(k) = s% gradT(k) ! prior to adjustments
261 :
262 79994 : s% Y_face_ad(k) = Y_face_ad
263 79994 : s% Y_face(k) = s% Y_face_ad(k)%val
264 :
265 79994 : s% mlt_vc_ad(k) = mlt_vc_ad
266 79994 : if (s% okay_to_set_mlt_vc) s% mlt_vc(k) = s% mlt_vc_ad(k)%val
267 :
268 79994 : s% mlt_D_ad(k) = D_ad
269 79994 : s% mlt_D(k) = D_ad%val
270 :
271 79994 : rho_face_ad = get_rho_face(s,k)
272 79994 : s% mlt_cdc(k) = s% mlt_D(k)*pow2(pi4*pow2(s%r(k))*rho_face_ad%val)
273 :
274 79994 : s% mlt_Gamma_ad(k) = Gamma_ad
275 79994 : s% mlt_Gamma(k) = Gamma_ad%val
276 :
277 79994 : s% gradr_ad(k) = gradr_ad
278 79994 : s% gradr(k) = s% gradr_ad(k)%val
279 :
280 79994 : s% gradL_ad(k) = s% grada_face_ad(k) + gradL_composition_term
281 79994 : s% gradL(k) = s% gradL_ad(k)%val
282 :
283 79994 : s% scale_height_ad(k) = scale_height_ad
284 79994 : s% scale_height(k) = scale_height_ad%val
285 :
286 79994 : s% Lambda_ad(k) = mixing_length_alpha*scale_height_ad
287 79994 : s% mlt_mixing_length(k) = s% Lambda_ad(k)%val
288 :
289 79994 : end subroutine store_results
290 :
291 0 : subroutine set_no_mixing(str)
292 : character (len=*) :: str
293 : include 'formats'
294 :
295 0 : s% mlt_mixing_type(k) = no_mixing
296 :
297 0 : s% grada_face_ad(k) = grada_face_ad
298 0 : s% grada_face(k) = grada_face_ad%val
299 :
300 : gradT_ad = gradr_ad
301 0 : s% gradT_ad(k) = gradT_ad
302 0 : s% gradT(k) = s% gradT_ad(k)%val
303 :
304 0 : Y_face_ad = gradT_ad - grada_face_ad
305 0 : s% Y_face_ad(k) = Y_face_ad
306 0 : s% Y_face(k) = s% Y_face_ad(k)%val
307 :
308 0 : s% mlt_vc_ad(k) = 0d0
309 0 : if (s% okay_to_set_mlt_vc) s% mlt_vc(k) = 0d0
310 :
311 0 : s% mlt_D_ad(k) = 0d0
312 0 : s% mlt_D(k) = 0d0
313 0 : s% mlt_cdc(k) = 0d0
314 :
315 0 : s% mlt_Gamma_ad(k) = 0d0
316 0 : s% mlt_Gamma(k) = 0d0
317 :
318 0 : s% gradr_ad(k) = gradr_ad
319 0 : s% gradr(k) = s% gradr_ad(k)%val
320 :
321 0 : s% gradL_ad(k) = 0d0
322 0 : s% gradL(k) = 0d0
323 :
324 0 : s% scale_height_ad(k) = scale_height_ad
325 0 : s% scale_height(k) = scale_height_ad%val
326 :
327 0 : s% Lambda_ad(k) = mixing_length_alpha*scale_height_ad
328 0 : s% mlt_mixing_length(k) = s% Lambda_ad(k)%val
329 :
330 0 : s% L_conv(k) = 0d0
331 :
332 0 : end subroutine set_no_mixing
333 :
334 : end subroutine do1_mlt_2
335 :
336 :
337 79994 : subroutine adjust_gradT_fraction(s,k,f)
338 : ! replace gradT by combo of grada_face and gradr
339 : ! then check excess
340 : use eos_def
341 : type (star_info), pointer :: s
342 : real(dp), intent(in) :: f
343 : integer, intent(in) :: k
344 : include 'formats'
345 79994 : if (f >= 0.0d0 .and. f <= 1.0d0) then
346 0 : if (f == 0d0) then
347 0 : s% gradT_ad(k) = s% gradr_ad(k)
348 : else ! mix
349 0 : s% gradT_ad(k) = f*s% grada_face_ad(k) + (1.0d0 - f)*s% gradr_ad(k)
350 : end if
351 0 : s% gradT(k) = s% gradT_ad(k)%val
352 : end if
353 79994 : call adjust_gradT_excess(s, k)
354 79994 : s% gradT_sub_grada(k) = s% gradT(k) - s% grada_face(k)
355 79994 : end subroutine adjust_gradT_fraction
356 :
357 :
358 79994 : subroutine adjust_gradT_excess(s, k)
359 79994 : use eos_def
360 : type (star_info), pointer :: s
361 : integer, intent(in) :: k
362 79994 : real(dp) :: alfa, log_tau, gradT_excess_alpha, gradT_sub_grada
363 : include 'formats'
364 : !s% gradT_excess_alpha is calculated at start of step and held constant during iterations
365 : ! gradT_excess_alpha = 0 means no efficiency boost; = 1 means full efficiency boost
366 79994 : gradT_excess_alpha = s% gradT_excess_alpha
367 79994 : s% gradT_excess_effect(k) = 0.0d0
368 79994 : gradT_sub_grada = s% gradT(k) - s% grada_face(k)
369 79994 : if (gradT_excess_alpha <= 0.0d0 .or. &
370 : gradT_sub_grada <= s% gradT_excess_f1) return
371 0 : if (s% lnT(k)/ln10 > s% gradT_excess_max_logT) return
372 0 : log_tau = log10(s% tau(k))
373 0 : if (log_tau < s% gradT_excess_max_log_tau_full_off) return
374 0 : if (log_tau < s% gradT_excess_min_log_tau_full_on) &
375 : gradT_excess_alpha = gradT_excess_alpha* &
376 : (log_tau - s% gradT_excess_max_log_tau_full_off)/ &
377 0 : (s% gradT_excess_min_log_tau_full_on - s% gradT_excess_max_log_tau_full_off)
378 0 : alfa = s% gradT_excess_f2 ! for full boost, use this fraction of gradT
379 0 : if (gradT_excess_alpha < 1) & ! only partial boost, so increase alfa
380 : ! alfa goes to 1 as gradT_excess_alpha goes to 0
381 : ! alfa unchanged as gradT_excess_alpha goes to 1
382 0 : alfa = alfa + (1d0 - alfa)*(1d0 - gradT_excess_alpha)
383 0 : s% gradT_ad(k) = alfa*s% gradT_ad(k) + (1d0 - alfa)*s% grada_face_ad(k)
384 0 : s% gradT(k) = s% gradT_ad(k)%val
385 0 : s% gradT_excess_effect(k) = 1d0 - alfa
386 79994 : end subroutine adjust_gradT_excess
387 :
388 :
389 30 : subroutine switch_to_radiative(s,k)
390 : type (star_info), pointer :: s
391 : integer, intent(in) :: k
392 30 : s% mlt_mixing_type(k) = no_mixing
393 30 : s% mlt_mixing_length(k) = 0
394 30 : s% mlt_D(k) = 0
395 30 : s% mlt_cdc(k) = 0d0
396 30 : s% mlt_vc(k) = 0
397 30 : s% gradT_ad(k) = s% gradr_ad(k)
398 30 : s% gradT(k) = s% gradT_ad(k)%val
399 79994 : end subroutine switch_to_radiative
400 :
401 :
402 : subroutine switch_to_adiabatic(s,k)
403 : use eos_def, only: i_grad_ad
404 : type (star_info), pointer :: s
405 : integer, intent(in) :: k
406 : s% gradT_ad(k) = s% grada_face_ad(k)
407 : s% gradT(k) = s% gradT_ad(k)%val
408 : end subroutine switch_to_adiabatic
409 :
410 :
411 21 : subroutine set_gradT_excess_alpha(s, ierr)
412 : use alloc
413 : use star_utils, only: get_Lrad_div_Ledd, after_C_burn
414 : use chem_def, only: ih1, ihe4
415 : type (star_info), pointer :: s
416 : integer, intent(out) :: ierr
417 21 : real(dp) :: beta, lambda, tmp, alpha, &
418 21 : beta_limit, lambda1, beta1, lambda2, beta2, dlambda, dbeta
419 : integer :: k, k_beta, k_lambda, nz, h1, he4
420 : include 'formats'
421 21 : ierr = 0
422 21 : if (.not. s% okay_to_reduce_gradT_excess) then
423 21 : s% gradT_excess_alpha = 0
424 21 : return
425 : end if
426 0 : nz = s% nz
427 0 : h1 = s% net_iso(ih1)
428 0 : if (h1 /= 0) then
429 0 : if (s% xa(h1,nz) > s% gradT_excess_max_center_h1) then
430 0 : s% gradT_excess_alpha = 0
431 0 : return
432 : end if
433 : end if
434 0 : he4 = s% net_iso(ihe4)
435 0 : if (he4 /= 0) then
436 0 : if (s% xa(he4,nz) < s% gradT_excess_min_center_he4) then
437 0 : s% gradT_excess_alpha = 0
438 0 : return
439 : end if
440 : end if
441 0 : beta = 1d0 ! beta = min over k of Pgas(k)/Peos(k)
442 0 : k_beta = 0
443 0 : do k=1,nz
444 0 : tmp = s% Pgas(k)/s% Peos(k)
445 0 : if (tmp < beta) then
446 0 : k_beta = k
447 0 : beta = tmp
448 : end if
449 : end do
450 0 : beta = beta*(1d0 + s% xa(1,nz))
451 0 : s% gradT_excess_min_beta = beta
452 0 : lambda = 0d0 ! lambda = max over k of Lrad(k)/Ledd(k)
453 0 : do k=2,k_beta
454 0 : tmp = get_Lrad_div_Ledd(s,k)
455 0 : if (tmp > lambda) then
456 0 : k_lambda = k
457 0 : lambda = tmp
458 : end if
459 : end do
460 0 : lambda = min(1d0,lambda)
461 0 : s% gradT_excess_max_lambda = lambda
462 0 : lambda1 = s% gradT_excess_lambda1
463 0 : beta1 = s% gradT_excess_beta1
464 0 : lambda2 = s% gradT_excess_lambda2
465 0 : beta2 = s% gradT_excess_beta2
466 0 : dlambda = s% gradT_excess_dlambda
467 0 : dbeta = s% gradT_excess_dbeta
468 : ! alpha is fraction of full boost to apply
469 : ! depends on location in (beta,lambda) plane
470 0 : if (lambda1 < 0) then
471 : alpha = 1
472 0 : else if (lambda >= lambda1) then
473 0 : if (beta <= beta1) then
474 : alpha = 1
475 0 : else if (beta < beta1 + dbeta) then
476 0 : alpha = (beta1 + dbeta - beta)/dbeta
477 : else ! beta >= beta1 + dbeta
478 : alpha = 0
479 : end if
480 0 : else if (lambda >= lambda2) then
481 : beta_limit = beta2 + &
482 0 : (lambda - lambda2)*(beta1 - beta2)/(lambda1 - lambda2)
483 0 : if (beta <= beta_limit) then
484 : alpha = 1
485 0 : else if (beta < beta_limit + dbeta) then
486 0 : alpha = (beta_limit + dbeta - beta)/dbeta
487 : else
488 : alpha = 0
489 : end if
490 0 : else if (lambda > lambda2 - dlambda) then
491 0 : if (beta <= beta2) then
492 : alpha = 1
493 0 : else if (beta < beta2 + dbeta) then
494 0 : alpha = (lambda - (lambda2 - dlambda))/dlambda
495 : else ! beta >= beta2 + dbeta
496 : alpha = 0
497 : end if
498 : else ! lambda <= lambda2 - dlambda
499 : alpha = 0
500 : end if
501 0 : if (s% generations > 1 .and. lambda1 >= 0) then ! time smoothing
502 : s% gradT_excess_alpha = &
503 : (1d0 - s% gradT_excess_age_fraction)*alpha + &
504 0 : s% gradT_excess_age_fraction*s% gradT_excess_alpha_old
505 0 : if (s% gradT_excess_max_change > 0d0) then
506 0 : if (s% gradT_excess_alpha > s% gradT_excess_alpha_old) then
507 : s% gradT_excess_alpha = min(s% gradT_excess_alpha, s% gradT_excess_alpha_old + &
508 0 : s% gradT_excess_max_change)
509 : else
510 : s% gradT_excess_alpha = max(s% gradT_excess_alpha, s% gradT_excess_alpha_old - &
511 0 : s% gradT_excess_max_change)
512 : end if
513 : end if
514 : else
515 0 : s% gradT_excess_alpha = alpha
516 : end if
517 0 : if (s% gradT_excess_alpha < 1d-4) s% gradT_excess_alpha = 0d0
518 0 : if (s% gradT_excess_alpha > 0.9999d0) s% gradT_excess_alpha = 1d0
519 21 : end subroutine set_gradT_excess_alpha
520 :
521 :
522 66 : subroutine check_for_redo_MLT(s, nzlo, nzhi, ierr)
523 : type (star_info), pointer :: s
524 : integer, intent(in) :: nzlo, nzhi
525 : integer, intent(out) :: ierr
526 : logical :: in_convective_region
527 : integer :: k, k_bot
528 66 : real(dp) :: bot_Hp, bot_r, top_Hp, top_r, dr
529 : logical :: dbg
530 : include 'formats'
531 : ! check_for_redo_MLT assumes that nzlo = 1, nzhi = nz
532 : ! that is presently true; make sure that assumption doesn't change
533 66 : if (.not. ((nzlo==1).and.(nzhi==s%nz))) then
534 0 : write(*,*) 'nzlo != 1 or nzhi != nz'
535 0 : call mesa_error(__FILE__,__LINE__)
536 : end if
537 66 : ierr = 0
538 66 : dbg = .false.
539 66 : bot_Hp = 0; bot_r = 0; top_Hp = 0; top_r = 0; dr = 0
540 66 : in_convective_region = (s% mlt_mixing_type(nzhi) == convective_mixing)
541 66 : k_bot = nzhi
542 66 : bot_r = s% r(k_bot)
543 66 : bot_Hp = s% scale_height(k_bot)
544 79928 : do k=nzhi-1, nzlo+1, -1
545 79928 : if (in_convective_region) then
546 11709 : if (s% mlt_mixing_type(k) /= convective_mixing) then
547 198 : call end_of_convective_region
548 : end if
549 : else ! in non-convective region
550 68153 : if (s% mlt_mixing_type(k) == convective_mixing) then
551 : ! start of a convective region
552 132 : k_bot = k+1
553 132 : in_convective_region = .true.
554 132 : bot_r = s% r(k_bot)
555 132 : bot_Hp = s% scale_height(k_bot)
556 : end if
557 : end if
558 : end do
559 87 : if (in_convective_region) then
560 0 : k = 1 ! end at top
561 0 : call end_of_convective_region
562 : end if
563 :
564 : contains
565 :
566 198 : subroutine end_of_convective_region()
567 : integer :: kk, op_err
568 198 : real(dp) :: Hp
569 : logical :: end_dbg
570 : 9 format(a40, 3i7, 99(1pd26.16))
571 : include 'formats'
572 198 : in_convective_region = .false.
573 198 : end_dbg = .false.
574 198 : top_r = s% r(k)
575 198 : top_Hp = s% scale_height(k)
576 198 : dr = top_r - bot_r
577 198 : Hp = (bot_Hp + top_Hp)/2
578 198 : if (dr < s% alpha_mlt(k)*min(top_Hp, bot_Hp) .and. &
579 : s% redo_conv_for_dr_lt_mixing_length) then
580 0 : !$OMP PARALLEL DO PRIVATE(kk,op_err) SCHEDULE(dynamic,2)
581 : do kk = k, k_bot
582 : op_err = 0
583 : call redo1_mlt(s,kk,dr,op_err)
584 : if (op_err /= 0) ierr = op_err
585 : end do
586 : !$OMP END PARALLEL DO
587 : end if
588 198 : end subroutine end_of_convective_region
589 :
590 0 : subroutine redo1_mlt(s, k, dr, ierr)
591 : type (star_info), pointer :: s
592 : integer, intent(in) :: k
593 : real(dp), intent(in) :: dr
594 : integer, intent(out) :: ierr
595 : logical :: make_gradr_sticky_in_solver_iters
596 : include 'formats'
597 0 : ierr = 0
598 0 : if (dr >= s% mlt_mixing_length(k)) return
599 : ! if convection zone is smaller than mixing length
600 : ! redo MLT with reduced alpha so mixing_length = dr
601 : call do1_mlt_2(s, k, make_gradr_sticky_in_solver_iters, ierr, &
602 0 : mixing_length_alpha_in = dr/s% scale_height(k))
603 : end subroutine redo1_mlt
604 :
605 : end subroutine check_for_redo_MLT
606 :
607 :
608 : end module turb_info
|
