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