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 : module turb_support
21 :
22 : use star_private_def
23 : use const_def, only: dp, crad, no_mixing
24 : use num_lib
25 : use utils_lib
26 : use auto_diff_support
27 : use star_utils
28 : use turb
29 :
30 : implicit none
31 :
32 : private
33 : public :: get_gradT
34 : public :: do1_mlt_eval
35 : public :: Get_results
36 :
37 : contains
38 :
39 : !> Determines if it is safe (physically) to use TDC instead of MLT.
40 : !!
41 : !! Currently we only know we have to fall back to MLT in cells that get touched
42 : !! by adjust_mass, because there the convection speeds at the start of the
43 : !! step can be badly out of whack. This can be disabled with TDC_adjust_mass_fallback_to_mlt
44 : !! to let those cells use TDC.
45 : !!
46 : !! @param s star pointer
47 : !! @param k face index
48 : !! @param fallback False if we can use TDC, True if we can fall back to MLT.
49 65210 : logical function check_if_must_fall_back_to_MLT(s, k) result(fallback)
50 : type (star_info), pointer :: s
51 : integer, intent(in) :: k
52 :
53 65210 : fallback = .false.
54 65210 : if (s% TDC_adjust_mass_fallback_to_mlt .and. abs(s%mstar_dot) > 1d-99 .and. k < s% k_const_mass) then
55 65210 : fallback = .true.
56 : end if
57 65210 : end function check_if_must_fall_back_to_MLT
58 :
59 0 : subroutine get_gradT(s, MLT_option, & ! used to create models
60 : r, L, T, P, opacity, rho, chiRho, chiT, Cp, gradr, grada, scale_height, &
61 : iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, &
62 : mixing_type, gradT, Y_face, conv_vel, D, Gamma, ierr)
63 : type (star_info), pointer :: s
64 : character (len=*), intent(in) :: MLT_option
65 : real(dp), intent(in) :: &
66 : r, L, T, P, opacity, rho, chiRho, chiT, Cp, gradr, grada, scale_height, &
67 : XH1, cgrav, m, gradL_composition_term, mixing_length_alpha
68 : integer, intent(in) :: iso
69 : real(dp), intent(out) :: gradT, Y_face, conv_vel, D, Gamma
70 : integer, intent(out) :: mixing_type, ierr
71 : type(auto_diff_real_star_order1) :: &
72 : gradr_ad, grada_ad, scale_height_ad, gradT_ad, Y_face_ad, mlt_vc_ad, D_ad, &
73 : Gamma_ad, r_ad, L_ad, T_ad, P_ad, opacity_ad, rho_ad, dV_ad, chiRho_ad, chiT_ad, Cp_ad, energy_ad
74 0 : ierr = 0
75 0 : r_ad = r
76 0 : L_ad = L
77 0 : T_ad = T
78 0 : P_ad = P
79 0 : opacity_ad = opacity
80 0 : rho_ad = rho
81 0 : dV_ad = 0d0
82 0 : chiRho_ad = chiRho
83 0 : chiT_ad = chiT
84 0 : Cp_ad = Cp
85 0 : gradr_ad = gradr
86 0 : grada_ad = grada
87 0 : energy_ad = 0d0 ! correct to a value
88 0 : scale_height_ad = scale_height
89 0 : if (s% use_other_mlt_results) then
90 : call s% other_mlt_results(s% id, 0, MLT_option, &
91 : r_ad, L_ad, T_ad, P_ad, opacity_ad, rho_ad, dV_ad, chiRho_ad, &
92 : chiT_ad, Cp_ad, gradr_ad, grada_ad, scale_height_ad, &
93 : iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, &
94 : s% alpha_semiconvection, s% thermohaline_coeff, &
95 0 : mixing_type, gradT_ad, Y_face_ad, mlt_vc_ad, D_ad, Gamma_ad, energy_ad, ierr)
96 : else
97 : call Get_results(s, 0, MLT_option, &
98 : r_ad, L_ad, T_ad, P_ad, opacity_ad, rho_ad, dV_ad, chiRho_ad, &
99 : chiT_ad, Cp_ad, gradr_ad, grada_ad, scale_height_ad, &
100 : iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, &
101 : s% alpha_semiconvection, s% thermohaline_coeff, &
102 0 : mixing_type, gradT_ad, Y_face_ad, mlt_vc_ad, D_ad, Gamma_ad, energy_ad, ierr)
103 : end if
104 0 : gradT = gradT_ad%val
105 0 : Y_face = Y_face_ad%val
106 0 : conv_vel = mlt_vc_ad%val
107 0 : D = D_ad%val
108 0 : Gamma = Gamma_ad%val
109 0 : end subroutine get_gradT
110 :
111 :
112 719946 : subroutine do1_mlt_eval( &
113 : s, k, MLT_option, gradL_composition_term, &
114 : T_in, P_in, energy_in, opacity_in, rho_in, chiRho_in, chiT_in, Cp_in, &
115 : gradr_in, grada, scale_height, mixing_length_alpha, &
116 : mixing_type, gradT, Y_face, mlt_vc, D, Gamma, ierr)
117 : use chem_def, only: ih1
118 : use const_def, only: ln10, convective_mixing
119 : use starspots, only: starspot_tweak_gradr
120 : type (star_info), pointer :: s
121 : integer, intent(in) :: k
122 : character (len=*), intent(in) :: MLT_option
123 : type(auto_diff_real_star_order1), intent(in) :: &
124 : T_in, P_in, energy_in, opacity_in, rho_in, chiRho_in, chiT_in, Cp_in, &
125 : gradr_in, grada, scale_height
126 : real(dp), intent(in) :: gradL_composition_term, mixing_length_alpha
127 : integer, intent(out) :: mixing_type
128 : type(auto_diff_real_star_order1), intent(out) :: &
129 : gradT, Y_face, mlt_vc, D, Gamma
130 : integer, intent(out) :: ierr
131 :
132 : real(dp) :: cgrav, m, XH1, P_theta, L_theta
133 : integer :: iso
134 : type(auto_diff_real_star_order1) :: gradr, r, L, T, P, opacity, rho, dV, &
135 : chiRho, chiT, Cp, rho_start, energy
136 : include 'formats'
137 79994 : ierr = 0
138 :
139 :
140 79994 : P = P_in ! if u_flag, should this be P_face_ad? (time centered in riemann)
141 79994 : if (s% include_mlt_in_velocity_time_centering) then
142 : ! could be cleaner with a wrapper for time_centered P and L
143 : if (s% using_velocity_time_centering .and. &
144 0 : s% include_P_in_velocity_time_centering .and. &
145 : s% lnT(k)/ln10 <= s% max_logT_for_include_P_and_L_in_velocity_time_centering) then
146 0 : P_theta = s% P_theta_for_velocity_time_centering
147 : else
148 0 : P_theta = 1d0
149 : end if
150 : ! consder building a wrapper : wrap_opt_time_center_L_00(s,k)
151 : if (s% using_velocity_time_centering .and. &
152 0 : s% include_L_in_velocity_time_centering .and. &
153 : s% lnT(k)/ln10 <= s% max_logT_for_include_P_and_L_in_velocity_time_centering) then
154 0 : L_theta = s% L_theta_for_velocity_time_centering
155 : else
156 0 : L_theta = 1d0
157 : end if
158 0 : L = L_theta*wrap_L_00(s, k) + (1d0 - L_theta)*s% L_start(k)
159 0 : if (s% use_face_reconstruction) then
160 0 : P = P_theta*P + (1d0-P_theta)*s% reconstructed_P_face_start(k)
161 : else
162 0 : P = P_theta*P + (1d0-P_theta)*s% Peos_face_start(k)
163 : end if
164 0 : r = wrap_opt_time_center_r_00(s,k)
165 : else
166 79994 : L = wrap_L_00(s,k)
167 79994 : r = wrap_r_00(s,k)
168 : end if
169 79994 : gradr = gradr_in
170 79994 : cgrav = s% cgrav(k)
171 79994 : m = s% m_grav(k)
172 79994 : T = T_in
173 79994 : opacity = opacity_in
174 79994 : rho = rho_in
175 79994 : rho_start = get_rho_start_face(s,k)
176 79994 : dV = 1d0/rho - 1d0/rho_start ! both variables are face wrapped.
177 79994 : chiRho = chiRho_in
178 79994 : chiT = chiT_in
179 79994 : Cp = Cp_in
180 79994 : energy = energy_in
181 79994 : iso = s% dominant_iso_for_thermohaline(k)
182 79994 : XH1 = s% xa(s% net_iso(ih1),k)
183 :
184 79994 : if (s% use_other_mlt_results) then
185 : call s% other_mlt_results(s% id, k, MLT_option, &
186 : r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, &
187 : iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, &
188 : s% alpha_semiconvection, s% thermohaline_coeff, &
189 0 : mixing_type, gradT, Y_face, mlt_vc, D, Gamma, energy, ierr)
190 : else
191 : ! starspot YREC routine
192 79994 : if (s% do_starspots) then
193 : !dV = 0d0 ! dV = 1/rho - 1/rho_start and we assume rho = rho_start.
194 0 : call starspot_tweak_gradr(s, P, gradr_in, gradr)
195 : end if
196 : call Get_results(s, k, MLT_option, &
197 : r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, &
198 : iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, &
199 : s% alpha_semiconvection, s% thermohaline_coeff, &
200 79994 : mixing_type, gradT, Y_face, mlt_vc, D, Gamma, energy, ierr)
201 : end if
202 79994 : if (ierr /= 0) return
203 :
204 : ! Do not use diffusion velocities as MLT/TDC turbulent velocity.
205 79994 : if (mixing_type /= convective_mixing) mlt_vc = 0d0
206 :
207 : end subroutine do1_mlt_eval
208 :
209 :
210 159988 : subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg
211 : r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, &
212 : iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, &
213 : alpha_semiconvection, thermohaline_coeff, &
214 : mixing_type, gradT, Y_face, conv_vel, D, Gamma, energy, ierr)
215 : use star_utils
216 : use tdc_hydro, only: compute_tdc_Eq_div_w_face
217 : type (star_info), pointer :: s
218 : integer, intent(in) :: k
219 : character (len=*), intent(in) :: MLT_option
220 : type(auto_diff_real_star_order1), intent(in) :: &
221 : r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, energy
222 : integer, intent(in) :: iso
223 : real(dp), intent(in) :: &
224 : XH1, cgrav, m, gradL_composition_term, &
225 : mixing_length_alpha, alpha_semiconvection, thermohaline_coeff
226 : integer, intent(out) :: mixing_type
227 : type(auto_diff_real_star_order1), intent(out) :: gradT, Y_face, conv_vel, D, Gamma
228 : integer, intent(out) :: ierr
229 :
230 : type(auto_diff_real_star_order1) :: Pr, Pg, grav, Lambda, gradL, beta
231 : real(dp) :: conv_vel_start, scale, max_conv_vel, Y_face_guess
232 :
233 : ! these are used by use_superad_reduction
234 : real(dp) :: Gamma_limit, scale_value1, scale_value2, diff_grads_limit, reduction_limit, lambda_limit
235 : type(auto_diff_real_star_order1) :: Lrad_div_Ledd, Gamma_inv_threshold, Gamma_factor, alfa0, &
236 : diff_grads_factor, Gamma_term, exp_limit, grad_scale, gradr_scaled, Eq_div_w, check_Eq, mlt_Pturb, Ptot
237 : logical :: test_partials, using_TDC, have_Y_face_guess
238 : logical, parameter :: report = .false.
239 : include 'formats'
240 :
241 : ! check if this particular k can be done with TDC
242 79994 : using_TDC = .false.
243 79994 : if (s% MLT_option == 'TDC') using_TDC = .true.
244 79994 : if (.not. s% have_mlt_vc) using_TDC = .false.
245 79994 : if (k <= 0 .or. s%dt <= 0d0) using_TDC = .false.
246 75066 : if (using_TDC) using_TDC = .not. check_if_must_fall_back_to_MLT(s, k)
247 :
248 : ! Pre-calculate some things.
249 79994 : Eq_div_w = 0d0
250 79994 : if ((s% v_flag .or. s% u_flag) .and. k > 0 ) then ! only include Eq_div_w if v_flag or u_flag is true.
251 0 : if (using_TDC .and. s% TDC_alpha_M > 0) then
252 0 : check_Eq = compute_tdc_Eq_div_w_face(s, k, ierr)
253 0 : Eq_div_w = check_Eq
254 : end if
255 : end if
256 :
257 : ! Wrap Pturb into P
258 79994 : if (s% okay_to_set_mlt_vc .and. s% include_mlt_Pturb_in_thermodynamic_gradients .and. k > 0) then
259 0 : mlt_Pturb = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*rho/3d0
260 0 : Ptot = P + mlt_Pturb
261 : else
262 79994 : Ptot = P
263 : end if
264 :
265 79994 : Pr = crad*pow4(T)/3d0
266 79994 : Pg = Ptot - Pr
267 79994 : beta = Pg / Ptot
268 79994 : Lambda = mixing_length_alpha*scale_height
269 :
270 79994 : if (k == 0) then
271 0 : grav = cgrav*m/pow2(r)
272 : else
273 79994 : grav = cgrav*m/pow2(r) !try replacing with wrap_geff_face(s,k)
274 : end if
275 :
276 79994 : if (s% use_Ledoux_criterion) then
277 0 : gradL = grada + gradL_composition_term ! Ledoux temperature gradient
278 : else
279 : gradL = grada
280 : end if
281 :
282 : ! maximum convection velocity.
283 79994 : if (k > 0) then
284 79994 : if (s% q(k) <= s% max_conv_vel_div_csound_maxq) then
285 0 : if (s% use_face_reconstruction) then
286 0 : max_conv_vel = s% reconstructed_csound_face(k)*s% max_conv_vel_div_csound
287 : else
288 0 : max_conv_vel = s% csound_face(k)*s% max_conv_vel_div_csound
289 : end if
290 : else
291 79994 : max_conv_vel = 1d99
292 : end if
293 : else ! if k == 0
294 0 : max_conv_vel = 1d99
295 : end if
296 :
297 :
298 : ! Initialize with no mixing
299 79994 : mixing_type = no_mixing
300 79994 : gradT = gradr
301 79994 : Y_face = gradT - gradL
302 79994 : conv_vel = 0d0
303 79994 : D = 0d0
304 79994 : Gamma = 0d0
305 79994 : if (k /= 0) s% superad_reduction_factor(k) = 1d0
306 :
307 : ! Bail if we asked for no mixing, or if parameters are bad.
308 : if (MLT_option == 'none' .or. beta < 1d-10 .or. mixing_length_alpha <= 0d0 .or. &
309 : opacity%val < 1d-10 .or. P%val < 1d-20 .or. T%val < 1d-10 .or. Rho%val < 1d-20 &
310 79994 : .or. m < 1d-10 .or. r%val < 1d-10 .or. cgrav < 1d-10) return
311 :
312 : !test_partials = (k == s% solver_test_partials_k)
313 79994 : test_partials = .false.
314 79994 : ierr = 0
315 79994 : if (k > 0) then
316 79994 : s% tdc_num_iters(k) = 0
317 : end if
318 :
319 : if (report) then
320 : write(*,'(A)')
321 : write(*,4) 'enter Get_results k slvr_itr model gradr grada scale_height ' // trim(MLT_option), &
322 : k, s% solver_iter, s% model_number, gradr%val, grada%val, scale_height%val
323 : end if
324 :
325 79994 : if (k >= 1) then
326 79994 : s% dvc_dt_TDC(k) = 0d0
327 : end if
328 79994 : if (using_TDC) then
329 : if (report) write(*,3) 'call set_TDC', k, s% solver_iter
330 65210 : if (s% okay_to_set_mlt_vc) then
331 42492 : conv_vel_start = s% mlt_vc_old(k)
332 : else
333 22718 : conv_vel_start = s% mlt_vc(k)
334 : end if
335 :
336 : ! Set scale for judging the TDC luminosity equation Q(Y)=0.
337 : ! Q has units of a luminosity, so the scale should be a luminosity.
338 65210 : if (s% solver_iter == 0) then
339 39833456 : scale = max(abs(s% L(k)), 1d-3*maxval(s% L(1:s% nz)))
340 : else
341 34744752 : scale = max(abs(s% L_start(k)), 1d-3*maxval(s% L_start(1:s% nz)))
342 : end if
343 :
344 65210 : have_Y_face_guess = s% use_TDC_Y_face_seeded_newton .and. s% doing_solver_iterations
345 : if (have_Y_face_guess) then
346 0 : Y_face_guess = s% Y_face(k)
347 : else
348 : ! Non-positive Y_face_guess values are ignored by the TDC seeded bracket.
349 65210 : Y_face_guess = 0d0
350 : end if
351 :
352 : call set_TDC(&
353 : conv_vel_start, mixing_length_alpha, s%TDC_alpha_D, s%TDC_alpha_R, s%TDC_alpha_Pt, &
354 : s%dt, cgrav, m, report, &
355 : mixing_type, scale, chiT, chiRho, gradr, r, Ptot, T, rho, dV, Cp, opacity, &
356 : scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), max_conv_vel, &
357 : Eq_div_w, grav, &
358 : s% include_mlt_corr_to_TDC, s% TDC_alpha_C, s% TDC_alpha_S, s% use_TDC_enthalpy_flux_limiter, energy, &
359 65210 : Y_face_guess, ierr)
360 65210 : s% dvc_dt_TDC(k) = (conv_vel%val - conv_vel_start) / s%dt
361 :
362 65210 : if (ierr /= 0) then
363 0 : if (s% report_ierr) write(*,*) 'ierr from set_TDC'
364 0 : return
365 : end if
366 :
367 : ! Experimental method to lower superadiabaticity. Call TDC again with an artificially reduced
368 : ! gradr if the resulting gradT would lead to the radiative luminosity approaching the Eddington
369 : ! limit, or when a density inversion is expected to happen.
370 : ! This is meant as an implicit alternative to okay_to_reduce_gradT_excess
371 65210 : if (s% use_superad_reduction) then
372 0 : call set_superad_reduction
373 0 : if (Gamma_factor > 1d0) then
374 : call set_TDC(&
375 : conv_vel_start, mixing_length_alpha, s%TDC_alpha_D, s%TDC_alpha_R, s%TDC_alpha_Pt, &
376 : s%dt, cgrav, m, report, &
377 : mixing_type, scale, chiT, chiRho, gradr_scaled, r, Ptot, T, rho, dV, Cp, opacity, &
378 : scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), max_conv_vel, &
379 : Eq_div_w, grav, &
380 : s% include_mlt_corr_to_TDC, s% TDC_alpha_C, s% TDC_alpha_S, s% use_TDC_enthalpy_flux_limiter, energy, &
381 0 : Y_face_guess, ierr)
382 0 : s% dvc_dt_TDC(k) = (conv_vel%val - conv_vel_start) / s%dt
383 0 : if (ierr /= 0) then
384 0 : if (s% report_ierr) write(*,*) 'ierr from set_TDC when using superad_reduction'
385 0 : return
386 : end if
387 : end if
388 : end if
389 :
390 14784 : else if (gradr > gradL) then
391 : if (report) write(*,3) 'call set_MLT', k, s% solver_iter
392 : call set_MLT(MLT_option, mixing_length_alpha, s% Henyey_MLT_nu_param, s% Henyey_MLT_y_param, &
393 : chiT, chiRho, Cp, grav, Lambda, rho, Ptot, T, opacity, &
394 : gradr, grada, gradL, &
395 2729 : Gamma, gradT, Y_face, conv_vel, D, mixing_type, max_conv_vel, ierr)
396 :
397 :
398 2729 : if (ierr /= 0) then
399 0 : if (s% report_ierr) write(*,*) 'ierr from set_MLT'
400 0 : return
401 : end if
402 :
403 : ! Experimental method to lower superadiabaticity. Call MLT again with an artificially reduced
404 : ! gradr if the resulting gradT would lead to the radiative luminosity approaching the Eddington
405 : ! limit, or when a density inversion is expected to happen.
406 : ! This is meant as an implicit alternative to okay_to_reduce_gradT_excess
407 2729 : if (s% use_superad_reduction) then
408 0 : call set_superad_reduction
409 0 : if (Gamma_factor > 1d0) then
410 : call set_MLT(MLT_option, mixing_length_alpha, s% Henyey_MLT_nu_param, s% Henyey_MLT_y_param, &
411 : chiT, chiRho, Cp, grav, Lambda, rho, Ptot, T, opacity, &
412 : gradr_scaled, grada, gradL, &
413 0 : Gamma, gradT, Y_face, conv_vel, D, mixing_type, max_conv_vel, ierr)
414 :
415 0 : if (ierr /= 0) then
416 0 : if (s% report_ierr) write(*,*) 'ierr from set_MLT when using superad_reduction'
417 0 : return
418 : end if
419 : end if
420 : end if
421 : end if
422 :
423 : ! If we're not convecting, try thermohaline and semiconvection.
424 79994 : if (mixing_type == no_mixing) then
425 68285 : if (gradL_composition_term < 0) then
426 : if (report) write(*,3) 'call set_thermohaline', k, s% solver_iter
427 : call set_thermohaline(s%thermohaline_option, Lambda, grada, gradr, T, opacity, rho, Cp, gradL_composition_term, &
428 : iso, XH1, thermohaline_coeff, &
429 0 : D, gradT, Y_face, conv_vel, mixing_type, ierr)
430 0 : if (ierr /= 0) then
431 0 : if (s% report_ierr) write(*,*) 'ierr from set_thermohaline'
432 0 : return
433 : end if
434 68285 : else if (gradr > grada) then
435 : if (report) write(*,3) 'call set_semiconvection', k, s% solver_iter
436 : call set_semiconvection(L, Lambda, m, T, Ptot, Pr, beta, opacity, rho, alpha_semiconvection, &
437 : s% semiconvection_option, cgrav, Cp, gradr, grada, gradL, &
438 : gradL_composition_term, &
439 0 : gradT, Y_face, conv_vel, D, mixing_type, ierr)
440 0 : if (ierr /= 0) then
441 0 : if (s% report_ierr) write(*,*) 'ierr from set_semiconvection'
442 0 : return
443 : end if
444 : end if
445 : end if
446 :
447 : ! If there's too-little mixing to bother, or we hit a bad value, fall back on no mixing.
448 79994 : if (D%val < s% remove_small_D_limit .or. is_bad(D%val)) then
449 : if (report) write(*,2) 'D < s% remove_small_D_limit', k, D%val, s% remove_small_D_limit
450 68285 : mixing_type = no_mixing
451 68285 : gradT = gradr
452 68285 : Y_face = gradT - gradL
453 68285 : conv_vel = 0d0
454 68285 : D = 0d0
455 68285 : Gamma = 0d0
456 : end if
457 :
458 : ! Prevent convection near center of model for MLT or TDC pulsations
459 : ! We don't check for the using_TDC flag, because mlt is sometimes called when using TDC
460 79994 : if (k > s% nz - s% TDC_num_innermost_cells_forced_nonturbulent .or. &
461 : k < s% TDC_num_outermost_cells_forced_nonturbulent) then
462 : if (report) write(*,2) 'make TDC center cells non-turbulent', k
463 0 : mixing_type = no_mixing
464 0 : gradT = gradr
465 0 : Y_face = gradT - gradL
466 0 : conv_vel = 0d0
467 0 : D = 0d0
468 0 : Gamma = 0d0
469 : end if
470 :
471 :
472 : contains
473 :
474 0 : subroutine set_superad_reduction()
475 0 : Gamma_limit = s% superad_reduction_Gamma_limit
476 0 : scale_value1 = s% superad_reduction_Gamma_limit_scale
477 0 : scale_value2 = s% superad_reduction_Gamma_inv_scale
478 0 : diff_grads_limit = s% superad_reduction_diff_grads_limit
479 0 : reduction_limit = s% superad_reduction_limit
480 0 : Lrad_div_Ledd = 4d0*crad/3d0*pow4(T)/P*gradT
481 0 : Gamma_inv_threshold = 4d0*(1d0-beta)/(4d0-3*beta)
482 :
483 0 : Gamma_factor = 1d0
484 0 : if (gradT > gradL) then
485 0 : if (Lrad_div_Ledd > Gamma_limit .or. Lrad_div_Ledd > Gamma_inv_threshold) then
486 0 : alfa0 = (gradT-gradL)/diff_grads_limit
487 0 : if (alfa0 < 1d0) then
488 0 : diff_grads_factor = -alfa0*alfa0*alfa0*(-10d0 + alfa0*(15d0 - 6d0*alfa0))
489 : else
490 0 : diff_grads_factor = 1d0
491 : end if
492 :
493 0 : Gamma_term = 0d0
494 : !if (Lrad_div_Ledd > Gamma_limit) then
495 : ! Gamma_term = Gamma_term + scale_value1*pow2(Lrad_div_Ledd/Gamma_limit-1d0)
496 : !end if
497 : !if (Lrad_div_Ledd% val > Gamma_inv_threshold) then
498 : ! Gamma_term = Gamma_term + scale_value2*pow2(Lrad_div_Ledd/Gamma_inv_threshold-1d0)
499 : !end if
500 0 : if (Lrad_div_Ledd > Gamma_limit) then
501 0 : alfa0 = Lrad_div_Ledd/Gamma_limit-1d0
502 0 : if (alfa0 < 1d0) then
503 0 : Gamma_term = Gamma_term + scale_value1*(0.5d0*alfa0*alfa0)
504 : else
505 0 : Gamma_term = Gamma_term + scale_value1*(alfa0-0.5d0)
506 : end if
507 : !Gamma_term = Gamma_term + scale_value1*pow2(Lrad_div_Ledd/Gamma_limit-1d0)
508 : end if
509 0 : if (Lrad_div_Ledd% val > Gamma_inv_threshold) then
510 0 : alfa0 = Lrad_div_Ledd/Gamma_inv_threshold-1d0
511 0 : if (alfa0 < 1d0) then
512 0 : Gamma_term = Gamma_term + scale_value1*(0.5d0*alfa0*alfa0)
513 : else
514 0 : Gamma_term = Gamma_term + scale_value1*(alfa0-0.5d0)
515 : end if
516 : !Gamma_term = Gamma_term + scale_value2*pow2(Lrad_div_Ledd/Gamma_inv_threshold-1d0)
517 : end if
518 :
519 0 : if (Gamma_term > 0d0) then
520 0 : Gamma_factor = Gamma_term/pow(beta,0.5d0)*diff_grads_factor
521 0 : Gamma_factor = Gamma_factor + 1d0
522 0 : if (reduction_limit > 1d0) then
523 0 : lambda_limit = 2d0/(reduction_limit-1d0)
524 0 : exp_limit = exp(-lambda_limit*(Gamma_factor-1d0))
525 0 : Gamma_factor = 2d0*(reduction_limit-1d0)*(1d0/(1d0+exp_limit)-0.5d0)+1d0
526 : end if
527 : end if
528 : end if
529 : end if
530 0 : if (k /= 0) s% superad_reduction_factor(k) = Gamma_factor% val
531 0 : if (Gamma_factor > 1d0) then
532 0 : grad_scale = (gradr-gradL)/(Gamma_factor*gradr) + gradL/gradr
533 0 : gradr_scaled = grad_scale*gradr
534 : end if
535 0 : end subroutine set_superad_reduction
536 : end subroutine Get_results
537 :
538 : end module turb_support
|
