Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-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 rotation_mix_info
21 :
22 : use const_def, only: dp, pi, pi4, qe, clight, crad, boltz_sigma, lsun, msun, rsun, one_third, &
23 : convective_mixing, overshoot_mixing
24 : use num_lib
25 : use utils_lib
26 : use star_private_def
27 : use magnetic_diffusion
28 :
29 : implicit none
30 :
31 : private
32 : public :: set_rotation_mixing_info, smooth_for_rotation
33 :
34 : real(dp), parameter :: Ri_crit = 0.25d0 ! critical Richardson number
35 : real(dp), parameter :: R_crit = 2500d0 ! critical Reynolds number
36 :
37 : integer, parameter :: i_DSI = 1
38 : integer, parameter :: i_SH = i_DSI + 1
39 : integer, parameter :: i_SSI = i_SH + 1
40 : integer, parameter :: i_ES = i_SSI + 1
41 : integer, parameter :: i_GSF = i_ES + 1
42 : integer, parameter :: i_ST = i_GSF + 1
43 : integer, parameter :: num_instabilities = i_ST
44 :
45 : contains
46 :
47 0 : subroutine set_rotation_mixing_info(s, ierr)
48 : use star_utils, only: weighed_smoothing
49 :
50 : type (star_info), pointer :: s
51 : integer, intent(out) :: ierr
52 :
53 0 : real(dp) :: f_mu, q, D_ST_old, nu_ST_old
54 : ! the following are all defined at cell boundaries
55 : real(dp), dimension(:), pointer :: & ! just copies of pointers
56 0 : r, m, L, j_rot, gradT, grada, grav, visc, Ri
57 : real(dp), dimension(:), allocatable :: & ! allocated temporary storage
58 0 : csound, rho, T, P, cp, cv, chiRho, abar, zbar, gradT_sub_grada, &
59 0 : opacity, kap_cond, gamma1, mu_alt, eps_nuc, enu, L_neu, delta, &
60 0 : scale_height, omega, cell_dr, &
61 0 : dRho, dr, dPressure, domega, d_mu, d_j_rot, &
62 0 : dRho_dr, dRho_dr_ad, dr2omega, H_T, &
63 0 : domega_dlnR, Hj, dlnR_domega, &
64 0 : t_dyn, t_kh, &
65 0 : Ri_mu, Ri_T, &
66 0 : ve0, ve_mu, &
67 0 : v_ssi, h_ssi, Ris_1, Ris_2, &
68 0 : v_es, H_es, &
69 0 : v_gsf, H_gsf, &
70 0 : N2, N2_mu
71 0 : real(dp), allocatable :: smooth_work(:,:), saved(:,:)
72 0 : logical, allocatable :: unstable(:,:) ! (num_instabilities, nz)
73 :
74 : integer :: nz, k, k0, which, op_err
75 0 : real(dp) :: alfa, growth_limit, age_fraction, &
76 0 : D_omega_source, dgtau, angsml, angsmt, prev_out_q, prev_out_q_m1, prev_q, prev_q_m1
77 :
78 : include 'formats'
79 :
80 0 : ierr = 0
81 0 : nz = s% nz
82 :
83 0 : s% D_visc(1:nz) = 0
84 0 : s% D_DSI(1:nz) = 0
85 0 : s% D_SH(1:nz) = 0
86 0 : s% D_SSI(1:nz) = 0
87 0 : s% D_ES(1:nz) = 0
88 0 : s% D_GSF(1:nz) = 0
89 0 : s% D_ST(1:nz) = 0
90 0 : s% nu_ST(1:nz) = 0
91 0 : s% dynamo_B_r(1:nz) = 0
92 0 : s% dynamo_B_phi(1:nz) = 0
93 0 : s% omega_shear(1:nz) = 0
94 0 : s% D_omega(1:nz) = 0
95 :
96 0 : if (all(s% omega(1:nz) == 0d0)) return
97 :
98 0 : call setup(ierr)
99 0 : if (failed('setup for set_rotation_mixing_info', ierr)) return
100 :
101 0 : unstable(:,1:nz) = .false.
102 0 : growth_limit = 1d-10
103 :
104 0 : !$OMP PARALLEL DO PRIVATE(which, k, q, age_fraction, op_err) SCHEDULE(dynamic,2)
105 : do which = 1, num_instabilities
106 :
107 : if (ierr /= 0) cycle
108 : op_err = 0
109 :
110 : select case (which)
111 :
112 : case (i_DSI)
113 :
114 : if (s% D_DSI_factor > 0 .or. s% am_nu_DSI_factor > 0) then
115 : call set_D_DSI(op_err)
116 : if (failed('set_D_DSI', op_err)) then
117 : ierr = -1; cycle
118 : end if
119 : call smooth_for_rotation(s, s% D_DSI, s% smooth_D_DSI, smooth_work(1:nz,which))
120 : if (s% skip_rotation_in_convection_zones) &
121 : call zero_if_convective(nz, s% mixing_type, s% D_mix, s% D_DSI)
122 : call zero_if_tiny(s,s% D_DSI)
123 : end if
124 :
125 : case (i_SH)
126 :
127 : if (s% D_SH_factor > 0 .or. s% am_nu_SH_factor > 0) then
128 : call set_D_SH(op_err)
129 : if (failed('set_D_SH', op_err)) then
130 : ierr = -1; cycle
131 : end if
132 : call smooth_for_rotation(s, s% D_SH, s% smooth_D_SH, smooth_work(1:nz,which))
133 : if (s% skip_rotation_in_convection_zones) &
134 : call zero_if_convective(nz, s% mixing_type, s% D_mix, s% D_SH)
135 : call zero_if_tiny(s,s% D_SH)
136 : end if
137 :
138 : case (i_SSI)
139 :
140 : if (s% D_SSI_factor > 0 .or. s% am_nu_SSI_factor > 0) then
141 : call set_D_SSI(op_err)
142 : if (failed('set_D_SSI', op_err)) then
143 : ierr = -1; cycle
144 : end if
145 : call smooth_for_rotation(s, s% D_SSI, s% smooth_D_SSI, smooth_work(1:nz,which))
146 : if (s% skip_rotation_in_convection_zones) &
147 : call zero_if_convective(nz, s% mixing_type, s% D_mix, s% D_SSI)
148 : call zero_if_tiny(s,s% D_SSI)
149 : end if
150 :
151 : case (i_ES)
152 :
153 : if (s% D_ES_factor > 0 .or. s% am_nu_ES_factor > 0) then
154 : call set_D_ES(op_err)
155 : if (failed('set_D_ES', op_err)) then
156 : ierr = -1; cycle
157 : end if
158 : call smooth_for_rotation(s, s% D_ES, s% smooth_D_ES, smooth_work(1:nz,which))
159 : if (s% skip_rotation_in_convection_zones) &
160 : call zero_if_convective(nz, s% mixing_type, s% D_mix, s% D_ES)
161 : call zero_if_tiny(s,s% D_ES)
162 : end if
163 :
164 : case (i_GSF)
165 :
166 : if (s% D_GSF_factor > 0 .or. s% am_nu_GSF_factor > 0) then
167 : call set_D_GSF(op_err)
168 : if (failed('set_D_GSF', op_err)) then
169 : ierr = -1; cycle
170 : end if
171 : call smooth_for_rotation(s, s% D_GSF, s% smooth_D_GSF, smooth_work(1:nz,which))
172 : if (s% skip_rotation_in_convection_zones) &
173 : call zero_if_convective(nz, s% mixing_type, s% D_mix, s% D_GSF)
174 : call zero_if_tiny(s,s% D_GSF)
175 : end if
176 :
177 : case (i_ST)
178 :
179 : if (s% D_ST_factor > 0 .or. s% am_nu_ST_factor > 0) then
180 :
181 : call set_ST(s, &
182 : rho, T, r, L, omega, Cp, abar, zbar, delta, grav, &
183 : N2, N2_mu, opacity, kap_cond, scale_height, op_err)
184 : if (failed('set_ST', op_err)) then
185 : ierr = -1; cycle
186 : end if
187 :
188 : call smooth_for_rotation(s, s% D_ST, s% smooth_D_ST, smooth_work(1:nz,which))
189 : call smooth_for_rotation(s, s% nu_ST, s% smooth_nu_ST, smooth_work(1:nz,which))
190 :
191 : ! time smooth for ST only
192 : if (s% prev_mesh_have_ST_start_info .and. .not. s% doing_relax &
193 : .and. .not. s% doing_first_model_of_run .and. s% ST_angsmt>0d0) then
194 :
195 : k0 = 1
196 : angsml = s% ST_angsml
197 : angsmt = s% ST_angsmt
198 : prev_out_q = 0d0
199 : prev_q = 1d0
200 : do k=2, nz-1 ! ignore k=1 or nz edge case
201 : if (s% m(k) > s% mstar_old) cycle
202 : do while (k0 < s% prev_mesh_nz)
203 : if (s% m(k)/s% mstar_old > prev_q) exit
204 : prev_out_q = prev_out_q + s% prev_mesh_dq(k0)
205 : prev_q = 1d0 - prev_out_q
206 : k0 = k0+1
207 : end do
208 : if (s% m(k)/s% mstar_old < prev_q) exit
209 : prev_out_q_m1 = prev_out_q - s% prev_mesh_dq(k0-1)
210 : prev_q_m1 = 1 - prev_out_q_m1
211 :
212 : ! linear interpolation
213 : alfa = (s% m(k)/s% mstar_old - prev_q)/(prev_q_m1-prev_q)
214 : D_ST_old = (1d0-alfa)*s% prev_mesh_D_ST_start(k0) + alfa*s% prev_mesh_D_ST_start(k0-1)
215 : nu_ST_old = (1d0-alfa)*s% prev_mesh_nu_ST_start(k0) + alfa*s% prev_mesh_nu_ST_start(k0-1)
216 : dgtau = angsml*(s% r(k)-s% r(k+1))*(s% r(k-1)*s% r(k))/s% dt
217 : s% D_ST(k) = max(D_ST_old/(1d0 + angsmt), &
218 : min(s% D_ST(k), max(D_ST_old*(1d0 + angsmt), D_ST_old + dgtau)))
219 : s% nu_ST(k) = max(nu_ST_old/(1d0 + angsmt), &
220 : min(s% nu_ST(k), max(nu_ST_old*(1d0 + angsmt), nu_ST_old + dgtau)))
221 : end do
222 : end if
223 :
224 : ! calculate B_r and B_phi
225 : do k = 1, nz
226 : q = s% omega_shear(k)
227 : s% dynamo_B_r(k) = & ! eqn 11, Heger et al. 2005
228 : pow(pow2(pi4*rho(k)*s% nu_ST(k)*q/r(k))*abs(omega(k))*s% nu_ST(k),0.25D0)
229 : s% dynamo_B_phi(k) = & ! eqn 12, Heger et al. 2005
230 : pow(pow2(pi4*rho(k)*omega(k)*q*r(k))*abs(omega(k))*s% nu_ST(k),0.25d0)
231 : end do
232 :
233 : if (s% skip_rotation_in_convection_zones) &
234 : call zero_if_convective(nz, s% mixing_type, s% D_mix, s% D_ST)
235 : if (s% skip_rotation_in_convection_zones) &
236 : call zero_if_convective(nz, s% mixing_type, s% D_mix, s% nu_ST)
237 : if (s% skip_rotation_in_convection_zones) &
238 : call zero_if_convective(nz, s% mixing_type, s% D_mix, s% dynamo_B_r)
239 : if (s% skip_rotation_in_convection_zones) &
240 : call zero_if_convective(nz, s% mixing_type, s% D_mix, s% dynamo_B_phi)
241 : call zero_if_tiny(s,s% D_ST)
242 : call zero_if_tiny(s,s% nu_ST)
243 :
244 : end if
245 :
246 : case default
247 : call mesa_error(__FILE__,__LINE__,'bad case for rotation_mix_info')
248 :
249 : end select
250 :
251 : end do
252 : !$OMP END PARALLEL DO
253 0 : if (failed('set_rotation_mixing_info instabilities', ierr)) return
254 :
255 0 : if (s% D_omega_flag .and. s% doing_finish_load_model) then
256 0 : do k=1,nz
257 0 : s% D_omega(k) = 0d0
258 : end do
259 0 : else if (s% D_omega_flag) then
260 :
261 0 : do k=1,nz
262 0 : if (s% q(k) <= s% max_q_for_D_omega_zero_in_convection_region .and. &
263 : s% mixing_type(k) == convective_mixing) then
264 0 : s% D_omega(k) = 0d0
265 0 : cycle
266 : end if
267 : D_omega_source = &
268 : s% D_DSI_factor * s% D_DSI(k) + &
269 : s% D_SH_factor * s% D_SH(k) + &
270 : s% D_SSI_factor * s% D_SSI(k) + &
271 : s% D_ES_factor * s% D_ES(k) + &
272 : s% D_GSF_factor * s% D_GSF(k) + &
273 0 : s% D_ST_factor * s% D_ST(k)
274 0 : if (is_bad(D_omega_source)) then
275 0 : write(*,2) 'D_omega_source', k, D_omega_source
276 0 : write(*,2) 's% D_DSI_factor * s% D_DSI(k)', k, s% D_DSI_factor * s% D_DSI(k), &
277 0 : s% D_DSI_factor, s% D_DSI(k)
278 0 : write(*,2) 's% D_SH_factor * s% D_SH(k)', k, s% D_SH_factor * s% D_SH(k), &
279 0 : s% D_SH_factor, s% D_SH(k)
280 0 : write(*,2) 's% D_SSI_factor * s% D_SSI(k)', k, s% D_SSI_factor * s% D_SSI(k), &
281 0 : s% D_SSI_factor, s% D_SSI(k)
282 0 : write(*,2) 's% D_ES_factor * s% D_ES(k)', k, s% D_ES_factor * s% D_ES(k), &
283 0 : s% D_ES_factor, s% D_ES(k)
284 0 : write(*,2) 's% D_GSF_factor * s% D_GSF(k)', k, s% D_GSF_factor * s% D_GSF(k), &
285 0 : s% D_GSF_factor, s% D_GSF(k)
286 0 : write(*,2) 's% D_ST_factor * s% D_ST(k)', k, s% D_ST_factor * s% D_ST(k), &
287 0 : s% D_ST_factor, s% D_ST(k)
288 0 : call mesa_error(__FILE__,__LINE__,'rotation mix')
289 : end if
290 0 : s% D_omega(k) = D_omega_source
291 0 : if (is_bad(s% D_omega(k))) then
292 0 : write(*,2) 's% D_omega(k)', k, s% D_omega(k)
293 0 : write(*,2) 'D_omega_source', k, D_omega_source
294 0 : call mesa_error(__FILE__,__LINE__,'rotation mix')
295 : end if
296 : end do
297 :
298 0 : if (s% smooth_D_omega > 0) then
299 0 : call smooth_for_rotation(s, s% D_omega, s% smooth_D_omega, smooth_work(1:nz,1))
300 0 : do k=1,nz
301 0 : if (is_bad(s% D_omega(k))) then
302 0 : write(*,2) 'after smooth_for_rotation s% D_omega(k)', k, s% D_omega(k)
303 0 : call mesa_error(__FILE__,__LINE__,'rotation mix')
304 : end if
305 : end do
306 : end if
307 :
308 0 : if (s% D_omega_mixing_rate > 0d0 .and. s% dt > 0) then
309 0 : call mix_D_omega
310 0 : do k=1,nz
311 0 : if (is_bad(s% D_omega(k))) then
312 0 : write(*,2) 'after mix_D_omega s% D_omega(k)', k, s% D_omega(k)
313 0 : call mesa_error(__FILE__,__LINE__,'rotation mix')
314 : end if
315 : end do
316 : end if
317 :
318 : end if
319 :
320 0 : if (s% D_omega_flag) then
321 0 : do k=1,nz
322 0 : if (is_bad(s% D_omega(k))) then
323 0 : write(*,2) 'before return s% D_omega(k)', k, s% D_omega(k)
324 0 : call mesa_error(__FILE__,__LINE__,'rotation mix')
325 : end if
326 0 : if (s% D_omega(k) < 0d0) s% D_omega(k) = 0d0
327 : end do
328 : end if
329 :
330 :
331 : contains
332 :
333 :
334 0 : subroutine mix_D_omega
335 : integer :: i, k, nz
336 : real(dp), dimension(:), allocatable :: & ! work vectors
337 0 : sig, rhs, d, du, dl, bp, vp, xp, x
338 : real(dp) :: &
339 0 : dt, rate, d_ddt_dm1, d_ddt_d00, d_ddt_dp1, m, &
340 0 : d_dt, d_dt_in, d_dt_out
341 : include 'formats'
342 :
343 0 : nz = s% nz
344 0 : dt = s% dt
345 0 : if (dt == 0) return
346 :
347 0 : allocate(sig(nz), rhs(nz), d(nz), du(nz), dl(nz), bp(nz), vp(nz), xp(nz), x(nz))
348 :
349 0 : rate = min(s% D_omega_mixing_rate, 1d0/dt)
350 0 : do k=2,nz-1
351 0 : if (s% D_omega(k) == 0 .or. s% D_omega(k+1) == 0) then
352 0 : sig(k) = 0
353 : else if ((.not. s% D_omega_mixing_across_convection_boundary) .and. &
354 0 : s% mixing_type(k) /= convective_mixing .and. &
355 : (s% mixing_type(k-1) == convective_mixing .or. &
356 : s% mixing_type(k+1) == convective_mixing)) then
357 0 : sig(k) = 0
358 : else
359 0 : sig(k) = rate*dt
360 : end if
361 : end do
362 0 : sig(1) = 0
363 0 : sig(nz) = 0
364 :
365 0 : do k=1,nz
366 0 : if (k < nz) then
367 0 : d_dt_in = sig(k)*(s% D_omega(k+1) - s% D_omega(k))
368 : else
369 0 : d_dt_in = -sig(k)*s% D_omega(k)
370 : end if
371 0 : if (k > 1) then
372 0 : d_dt_out = sig(k-1)*(s% D_omega(k) - s% D_omega(k-1))
373 0 : d_ddt_dm1 = sig(k-1)
374 0 : d_ddt_d00 = -(sig(k-1) + sig(k))
375 : else
376 0 : d_dt_out = 0
377 0 : d_ddt_dm1 = 0
378 0 : d_ddt_d00 = -sig(k)
379 : end if
380 0 : d_dt = d_dt_in - d_dt_out
381 0 : d_ddt_dp1 = sig(k)
382 0 : rhs(k) = d_dt
383 0 : d(k) = 1d0 - d_ddt_d00
384 0 : if (k < nz) then
385 0 : du(k) = -d_ddt_dp1
386 : else
387 0 : du(k) = 0
388 : end if
389 0 : if (k > 1) dl(k-1) = -d_ddt_dm1
390 : end do
391 0 : dl(nz) = 0
392 :
393 : ! solve tridiagonal
394 0 : bp(1) = d(1)
395 0 : vp(1) = rhs(1)
396 0 : do i = 2,nz
397 0 : m = dl(i-1)/bp(i-1)
398 0 : bp(i) = d(i) - m*du(i-1)
399 0 : vp(i) = rhs(i) - m*vp(i-1)
400 : end do
401 0 : xp(nz) = vp(nz)/bp(nz)
402 0 : x(nz) = xp(nz)
403 0 : do i = nz-1, 1, -1
404 0 : xp(i) = (vp(i) - du(i)*xp(i+1))/bp(i)
405 0 : x(i) = xp(i)
406 : end do
407 :
408 0 : do k=2,nz
409 0 : if (is_bad(x(k))) then
410 : return
411 : write(*,3) 's% D_omega(k) prev, x', k, s% model_number, s% D_omega(k), x(k), bp(i)
412 : call mesa_error(__FILE__,__LINE__,'mix_D_omega')
413 : end if
414 : end do
415 :
416 : ! update D_omega
417 :
418 0 : do k=2,nz
419 0 : s% D_omega(k) = s% D_omega(k) + x(k)
420 0 : if (is_bad(s% D_omega(k))) then
421 0 : write(*,3) 's% D_omega(k)', k, s% model_number, s% D_omega(k)
422 0 : call mesa_error(__FILE__,__LINE__,'mix_D_omega')
423 : end if
424 0 : if (s% D_omega(k) < 0d0) s% D_omega(k) = 0d0
425 : end do
426 0 : s% D_omega(1) = 0d0
427 :
428 0 : end subroutine mix_D_omega
429 :
430 :
431 0 : subroutine setup(ierr)
432 : use kap_lib, only: kap_get_elect_cond_opacity
433 : integer, intent(out) :: ierr
434 :
435 : real(dp) :: &
436 0 : bracket_term, ri0, alfa, beta, enum1, enu00, &
437 0 : rho6, gamma, mu_e, rm23, ctmp, xi2, dynvisc, denom, &
438 0 : eps_nucm1, eps_nuc00, scale_height2, dlnRho_dlnP, dlnT_dlnP, &
439 0 : kap, dlnkap_dlnRho, dlnkap_dlnT
440 :
441 : integer :: i, k
442 :
443 : include 'formats'
444 :
445 0 : ierr = 0
446 0 : nz = s% nz
447 :
448 0 : f_mu = s% am_gradmu_factor
449 :
450 : ! copy some pointers
451 0 : r => s% r
452 0 : m => s% m
453 0 : L => s% L
454 0 : j_rot => s% j_rot
455 0 : gradT => s% gradT
456 0 : grada => s% grada_face
457 0 : grav => s% grav
458 0 : visc => s% D_visc
459 0 : Ri => s% richardson_number
460 :
461 : allocate( &
462 0 : csound(nz), rho(nz), T(nz), P(nz), cp(nz), cv(nz), chiRho(nz), abar(nz), zbar(nz), &
463 0 : opacity(nz), kap_cond(nz), gamma1(nz), mu_alt(nz), omega(nz), cell_dr(nz), eps_nuc(nz), enu(nz), L_neu(nz), &
464 0 : gradT_sub_grada(nz), delta(nz), scale_height(nz), &
465 0 : dRho(nz), dr(nz), dPressure(nz), domega(nz), d_mu(nz), d_j_rot(nz), &
466 0 : dRho_dr(nz), dRho_dr_ad(nz), dr2omega(nz), H_T(nz), &
467 0 : domega_dlnR(nz), Hj(nz), dlnR_domega(nz), t_dyn(nz), t_kh(nz), &
468 0 : Ri_mu(nz), Ri_T(nz), ve0(nz), ve_mu(nz), &
469 0 : v_ssi(nz), h_ssi(nz), Ris_1(nz), Ris_2(nz), &
470 0 : v_es(nz), H_es(nz), v_gsf(nz), H_gsf(nz), N2(nz), N2_mu(nz), &
471 0 : smooth_work(nz,num_instabilities), saved(nz,num_instabilities), &
472 0 : unstable(num_instabilities,nz))
473 :
474 : ! interpolate by mass to get values at cell boundaries
475 0 : enu00 = s% eps_nuc_neu_total(1) + s% non_nuc_neu(1)
476 0 : enu(1) = enu00
477 0 : eps_nuc00 = s% eps_nuc(1)
478 0 : eps_nuc(1) = eps_nuc00
479 0 : csound(1) = s% csound(1)
480 0 : rho(1) = s% rho(1)
481 0 : T(1) = s% T(1)
482 0 : P(1) = s% Peos(1)
483 0 : cp(1) = s% cp(1)
484 0 : cv(1) = s% Cv(1)
485 0 : chiRho(1) = s% chiRho(1)
486 0 : abar(1) = s% abar(1)
487 0 : zbar(1) = s% zbar(1)
488 0 : opacity(1) = s% opacity(1)
489 0 : gamma1(1) = s% gamma1(1)
490 0 : mu_alt(1) = s% abar(1)/(1 + s% zbar(1))
491 0 : delta(1) = s% chiT(1)/s% chiRho(1)
492 0 : L_neu(1) = enu00*s% dm(1)
493 0 : do k = 2, nz
494 0 : alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
495 0 : beta = 1 - alfa
496 0 : enum1 = enu00
497 0 : enu00 = s% eps_nuc_neu_total(k) + s% non_nuc_neu(k)
498 0 : enu(k) = alfa*enu00 + beta*enum1
499 0 : eps_nucm1 = eps_nuc00
500 0 : eps_nuc00 = s% eps_nuc(k)
501 0 : eps_nuc(k) = alfa*eps_nuc00 + beta*eps_nucm1
502 0 : csound(k) = alfa*s% csound(k) + beta*s% csound(k-1)
503 0 : rho(k) = alfa*s% rho(k) + beta*s% rho(k-1)
504 0 : T(k) = alfa*s% T(k) + beta*s% T(k-1)
505 0 : P(k) = alfa*s% Peos(k) + beta*s% Peos(k-1)
506 0 : cp(k) = alfa*s% cp(k) + beta*s% cp(k-1)
507 0 : cv(k) = alfa*s% Cv(k) + beta*s% Cv(k-1)
508 0 : chiRho(k) = alfa*s% chiRho(k) + beta*s% chiRho(k-1)
509 0 : abar(k) = alfa*s% abar(k) + beta*s% abar(k-1)
510 0 : zbar(k) = alfa*s% zbar(k) + beta*s% zbar(k-1)
511 0 : opacity(k) = alfa*s% opacity(k) + beta*s% opacity(k-1)
512 0 : gamma1(k) = alfa*s% gamma1(k) + beta*s% gamma1(k-1)
513 0 : mu_alt(k) = alfa*s% abar(k)/(1 + s% zbar(k)) + beta*s% abar(k-1)/(1 + s% zbar(k-1))
514 0 : delta(k) = alfa*s% chiT(k)/s% chiRho(k) + beta*s% chiT(k-1)/s% chiRho(k-1)
515 0 : L_neu(k) = enu00*s% dm(k) + L_neu(k-1)
516 0 : cell_dr(k-1) = s% rmid(k-1) - s% rmid(k)
517 : end do
518 0 : cell_dr(nz) = s% rmid(nz) - s% R_center
519 :
520 0 : do i = 1, nz
521 0 : gradT_sub_grada(i) = s% gradT(i) - s% grada_face(i)
522 : gradT_sub_grada(i) = & ! make sure it isn't too close to 0
523 0 : sign(max(abs(gradT_sub_grada(i)),1d-99),gradT_sub_grada(i))
524 0 : scale_height(i) = P(i)*r(i)*r(i)/(s% cgrav(i)*m(i)*rho(i))
525 0 : scale_height2 = sqrt(P(i)/s% cgrav(i))/rho(i)
526 0 : if (scale_height2 < scale_height(i)) scale_height(i) = scale_height2
527 0 : omega(i) = s% omega(i)
528 : end do
529 :
530 : ! differences (at cell boundaries)
531 0 : do i = 2, nz-1
532 0 : dRho(i) = rho(i-1) - rho(i)
533 0 : dr(i) = max(1d0, 0.5D0*(r(i-1) - r(i+1)))
534 0 : dPressure(i) = min(-1d-10, P(i-1) - P(i))
535 0 : d_mu(i) = mu_alt(i-1) - mu_alt(i)
536 0 : d_j_rot(i) = j_rot(i-1) - j_rot(i)
537 0 : domega(i) = 0.5D0*(omega(i-1) - omega(i+1))
538 : end do
539 :
540 0 : dRho(1) = 0
541 0 : dRho(nz) = 0
542 :
543 0 : dr(1) = max(1d0,r(1)-r(2))
544 0 : dr(nz) = r(nz) - s% R_center
545 :
546 0 : dPressure(1) = 0
547 0 : dPressure(nz) = 0
548 :
549 0 : d_mu(1) = 0
550 0 : d_mu(nz) = 0
551 :
552 0 : d_j_rot(1) = 0
553 0 : d_j_rot(nz) = 0
554 :
555 0 : domega(1) = 0
556 0 : domega(nz) = 0
557 :
558 0 : do i = 2, nz-1
559 0 : dRho_dr(i) = dRho(i)/dr(i)
560 0 : dRho_dr_ad(i) = rho(i)*dPressure(i)/(P(i)*gamma1(i)*dr(i))
561 0 : dr2omega(i) = 4.5d0*j_rot(i)*d_j_rot(i)/dr(i) ! d(r^2 omega)^2/dr using i = (2/3)*r^2
562 0 : domega_dlnR(i) = domega(i)*r(i)/dr(i)
563 0 : if (gradT(i) > 1d-20) then
564 0 : H_T(i) = scale_height(i)/gradT(i) ! -dr/dlnT, scale height of temperature
565 : else
566 0 : H_T(i) = scale_height(i)
567 : end if
568 0 : if(abs(d_j_rot(i))>0.d0)then
569 0 : Hj(i) = j_rot(i)*dr(i)/d_j_rot(i)
570 : else
571 0 : Hj(i) =0d0
572 : end if
573 : ! dr/dlnj, scale height of angular momentum
574 : end do
575 :
576 0 : dRho_dr(1) = 0; dRho_dr(nz) = 0
577 0 : dRho_dr_ad(1) = 0; dRho_dr_ad(nz) = 0
578 0 : dr2omega(1) = 0; dr2omega(nz) = 0
579 0 : domega_dlnR(1) = 0; domega_dlnR(nz) = 0
580 0 : H_T(1) = H_T(2); H_T(nz) = H_T(nz-1)
581 0 : Hj(1) = Hj(2); Hj(nz) = Hj(nz-1)
582 :
583 0 : do i = 1, nz
584 0 : dlnRho_dlnP = s% grad_density(i)
585 0 : dlnT_dlnP = s% grad_temperature(i)
586 0 : N2(i) = -grav(i)*(1/gamma1(i) - dlnRho_dlnP)/scale_height(i)
587 0 : N2_mu(i) = -grav(i)/scale_height(i)*(1/chiRho(i) - delta(i)*dlnT_dlnP - dlnRho_dlnP)
588 : end do
589 :
590 0 : do k=1,nz
591 0 : s% domega_dlnR(k) = domega_dlnR(k)
592 : end do
593 :
594 : ! safe inverse of domega/dlnR
595 0 : do i = 2, nz-1
596 0 : dlnR_domega(i) = sign(1d0/max(abs(domega_dlnR(i)),1d-30),domega_dlnR(i))
597 : end do
598 0 : dlnR_domega(1) = 0; dlnR_domega(nz) = 0
599 :
600 0 : do k = 2, nz - 1 ! shear
601 : s% omega_shear(k) = &
602 0 : (omega(k-1) - omega(k+1))/(r(k-1) - r(k+1))*(r(k)/omega(k))
603 0 : s% omega_shear(k) = max(1d-30,min(1d30,abs(s% omega_shear(k))))
604 : end do
605 0 : s% omega_shear(1) = 0
606 0 : s% omega_shear(nz) = 0
607 :
608 : ! timescales
609 0 : do i = 1, nz
610 0 : t_dyn(i) = sqrt(r(i)*r(i)*r(i)/(s% cgrav(i)*m(i)))
611 0 : t_kh(i) = s% cgrav(i)*m(i)*m(i)/(r(i)*max(1d0,L(i)+L_neu(i)))
612 : end do
613 :
614 : ! Richardson numbers (Heger 2000, eqn 20)
615 0 : do i = 2, nz-1
616 0 : ri0 = (rho(i)*delta(i)/P(i))*pow2(dlnR_domega(i)*grav(i))
617 0 : Ri_T(i) = ri0*max(0d0,-gradT_sub_grada(i)) ! turn off Ri_T in convection zones
618 0 : Ri_mu(i) = ri0*f_mu*s% smoothed_brunt_B(i)
619 : end do
620 0 : Ri_T(1) = 0; Ri_T(nz) = 0
621 0 : Ri_mu(1) = 0; Ri_mu(nz) = 0
622 0 : do i=1,nz
623 0 : if (N2(i) < 0d0) then
624 0 : Ri(i) = 1d0 ! disable in convective region
625 : else
626 0 : Ri(i) = Ri_T(i) + Ri_mu(i)
627 : end if
628 : end do
629 :
630 : ! dynamic viscosity
631 0 : do i=1,nz
632 0 : rho6 = rho(i)*1d-6
633 0 : gamma = 0.2275d0*zbar(i)*zbar(i)*pow(rho6/abar(i),one_third)*1.d8/T(i)
634 : ! gamma => eq (5) of Itoh et al 1987 ApJ 317,733
635 : ! electron viscosity according to Nandkumar & Pethick 1984 MNRAS
636 0 : mu_e = abar(i)/zbar(i)
637 0 : rm23 = pow(rho6/mu_e,2d0/3d0)
638 0 : ctmp = 1d0 + 1.018d0*rm23
639 : xi2 = sqrt(pi/3d0)*log(zbar(i))/3d0 + 2d0*log(1.32d0+2.33d0/sqrt(gamma)) - &
640 0 : 0.475d0*(1d0+2.036d0*rm23)/ctmp + 0.276d0*rm23/ctmp
641 : ! dynamic shear viscosity according to Wallenborn and Bauss (1978)
642 : ! and also Itoh et al 1987 ApJ 317,733
643 : ! fitting formula for eta* in Eq (12) of Itoh et al. 1987
644 : ctmp = -0.016321227d0+1.0198850d0*pow(gamma,-1.9217970d0) + &
645 0 : 0.024113535d0*pow(gamma,0.49999098d0)
646 : ! dynamic shear viscosity
647 0 : dynvisc = 5.53d3*zbar(i)*pow(rho6,5d0/6d0)*ctmp/pow(abar(i),one_third)
648 : ! add contribution of radiation
649 0 : dynvisc = dynvisc + 4.D0*crad*pow4(T(i))/(15.D0*clight*opacity(i)*rho(i))
650 : ! add contribution of electrons
651 0 : dynvisc = dynvisc + 1.893d6*pow(rm23,2.5d0)/(zbar(i)*ctmp*xi2)
652 : ! kinematic shear viscosity
653 0 : visc(i) = dynvisc/rho(i)
654 : end do
655 :
656 : ! velocities for ES and GSF
657 0 : if (s% D_ES_factor > 0 .or. s% D_GSF_factor > 0) then
658 0 : do i = 1, nz ! Heger 2000, eqns 35 and 36
659 : ! the bracket_term blows up at center since r^2/L and r^2/m can both -> Inf
660 : ! so bullet proof by including lower bounds
661 : bracket_term = &
662 : 2*r(i)*r(i)*(eps_nuc(i)/max(1d-3*Lsun,abs(L(i))) - 1/max(1d-3*Msun,m(i))) - &
663 0 : 3/(pi4*rho(i)*max(1d-3*Rsun,r(i)))
664 0 : if (abs(gradT_sub_grada(i)) < 1d-50) then
665 0 : ve0(i) = 1d99
666 0 : ve_mu(i) = 1d99
667 : else
668 0 : denom = (-gradT_sub_grada(i)*delta(i)*pow2(s% cgrav(i)*m(i)))
669 0 : if (abs(denom) < 1d-50 .or. is_bad(denom)) then
670 0 : if (s% stop_for_bad_nums) then
671 0 : write(*,2) 'denom', i, denom
672 0 : call mesa_error(__FILE__,__LINE__,'rotation mix info: velocities for ES and GSF')
673 : end if
674 0 : ve0(i) = 1d99
675 : else
676 0 : ve0(i) = grada(i)*omega(i)*omega(i)*r(i)*r(i)*r(i)*L(i)*bracket_term/denom
677 : end if
678 : ve_mu(i) = (scale_height(i)/t_kh(i))* &
679 0 : (f_mu*s% smoothed_brunt_B(i))/(gradT_sub_grada(i))
680 : end if
681 0 : if (is_bad(ve0(i))) then
682 0 : if (s% stop_for_bad_nums) then
683 0 : write(*,2) 've0(i)', i, ve0(i)
684 0 : call mesa_error(__FILE__,__LINE__,'rotation mix info')
685 : end if
686 0 : ve0(i) = 1d99
687 : end if
688 0 : if (is_bad(ve_mu(i))) then
689 0 : if (s% stop_for_bad_nums) then
690 0 : write(*,2) 've_mu(i)', i, ve_mu(i)
691 0 : call mesa_error(__FILE__,__LINE__,'rotation mix info')
692 : end if
693 0 : ve_mu(i) = 1d99
694 : end if
695 :
696 0 : if (s% model_number == -1 .and. i == 316) then
697 0 : write(*,2) 'grada(i)', i, grada(i)
698 0 : write(*,2) 'gradT(i)', i, gradT(i)
699 0 : write(*,2) 'gradT_sub_grada(i)', i, gradT_sub_grada(i)
700 0 : write(*,2) 's% smoothed_brunt_B(i)', i, s% smoothed_brunt_B(i)
701 0 : write(*,2) 'omega(i)', i, omega(i)
702 0 : write(*,2) 's% omega(i)', i, s% omega(i)
703 0 : write(*,2) '2*r**2*eps_nuc/L', i, 2*r(i)*r(i)*eps_nuc(i)/max(1d0,L(i))
704 0 : write(*,2) '2*r**2/m', i, 2*r(i)*r(i)/m(i)
705 0 : write(*,2) '3/(4*pi*rho*r)', i, 3/(4*pi*rho(i)*r(i))
706 0 : write(*,2) 've0(i)', i, ve0(i)
707 : end if
708 :
709 : end do
710 :
711 : ! conductive opacities for ST
712 0 : if (s% D_ST_factor > 0d0 .or. s% am_nu_factor > 0d0) then
713 0 : do i = 1, nz
714 : call kap_get_elect_cond_opacity( &
715 : s% kap_handle, zbar(i), log10(rho(i)), log10(T(i)), &
716 0 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
717 0 : if (ierr /= 0) return
718 0 : kap_cond(i) = kap
719 : end do
720 : end if
721 :
722 : end if
723 :
724 0 : end subroutine setup
725 :
726 :
727 0 : subroutine set_D_DSI(ierr)
728 : integer, intent(out) :: ierr
729 : integer :: i, k, kbot, ktop
730 0 : real(dp) :: instability_height, height, D
731 : logical, parameter :: dbg = .false.
732 : include 'formats'
733 :
734 0 : ierr = 0
735 0 : do i = 1, nz
736 0 : unstable(i_DSI,i) = (Ri(i) < Ri_crit) .and. (gradT_sub_grada(i) < 0)
737 : ! stable in convective regions where gradT >= grada
738 : end do
739 :
740 0 : kbot = nz
741 0 : do i = nz-1, 1, -1
742 :
743 0 : if (unstable(i_DSI,i) .and. .not. unstable(i_DSI,i+1)) kbot = i
744 :
745 : if (unstable(i_DSI,i+1) .and. &
746 0 : (i == 1 .or. .not. unstable(i_DSI,i)) .and. kbot > 1) then
747 :
748 0 : if (unstable(i_DSI,i)) then
749 : ktop = i
750 : else
751 0 : ktop = i+1
752 : end if
753 :
754 0 : if (ktop >= kbot) cycle
755 :
756 0 : instability_height = r(ktop) - r(kbot)
757 : if (dbg) write(*,3) 'DSI: ktop, kbot', ktop, kbot
758 0 : do k = ktop, kbot
759 0 : if (.not. unstable(i_DSI,k)) then
760 0 : write(*,2) 'D_DSI where stable?', k, s% q(k), &
761 0 : s% D_DSI(k), D, scale_height(k)*csound(k), &
762 0 : Ri(k), Ri_crit, height, t_dyn(k), &
763 0 : instability_height, scale_height(k), csound(k)
764 0 : call mesa_error(__FILE__,__LINE__,'set_D_DSI')
765 : end if
766 0 : height = min(instability_height, scale_height(k))
767 0 : D = height*height/t_dyn(k)
768 0 : s% D_DSI(k) = min(D, scale_height(k)*csound(k))
769 : if (dbg) write(*,2) 'D_DSI', k, s% q(k), &
770 : s% D_DSI(k), D, scale_height(k)*csound(k), &
771 : Ri(k), Ri_crit, height, t_dyn(k), &
772 0 : instability_height, scale_height(k), csound(k)
773 : end do
774 :
775 : if (dbg) write(*,*)
776 :
777 : end if
778 :
779 : end do
780 : if (dbg) call mesa_error(__FILE__,__LINE__,'set_D_DSI')
781 0 : end subroutine set_D_DSI
782 :
783 :
784 0 : subroutine set_D_SH(ierr) ! comment in Langer code says "DO NOT USE"
785 : integer, intent(out) :: ierr
786 : integer :: i, k, kbot, ktop
787 0 : real(dp) :: instability_height, height, D
788 0 : ierr = 0
789 :
790 0 : do i = 1, nz
791 0 : D = grav(i)/rho(i)*(dRho_dr_ad(i)-dRho_dr(i))+dr2omega(i)/(r(i)*r(i)*r(i))
792 0 : if (D < 0) then
793 0 : unstable(i_SH,i) = .true.
794 0 : s% D_SH(i) = D ! save for later
795 : else
796 0 : s% D_SH(i) = 0
797 : end if
798 : end do
799 :
800 0 : kbot = nz
801 0 : do i = nz-1, 1, -1
802 :
803 0 : if (unstable(i_SH,i) .and. .not. unstable(i_SH,i+1)) kbot = i
804 :
805 : if (unstable(i_SH,i+1) .and. &
806 0 : (i == 1 .or. .not. unstable(i_SH,i)) .and. kbot > 1) then
807 0 : if (unstable(i_SH,i)) then
808 : ktop = i
809 : else
810 0 : ktop = i+1
811 : end if
812 0 : if (ktop >= kbot) cycle
813 0 : instability_height = r(ktop) - r(kbot)
814 0 : do k = ktop, kbot
815 0 : height = min(instability_height, scale_height(k))
816 : ! use the previously calculated value saved in D_SH
817 0 : D = pow2(height*s% D_SH(k)*r(k)/grav(k))/t_dyn(k)
818 0 : s% D_SH(k) = min(D, scale_height(k)*csound(k))
819 : end do
820 : end if
821 :
822 : end do
823 0 : end subroutine set_D_SH
824 :
825 :
826 0 : subroutine set_D_SSI(ierr)
827 : use chem_def
828 : integer, intent(out) :: ierr
829 : integer :: i, k, kbot, ktop
830 0 : real(dp) :: qe3, qe4, dynvisc, Prandtl, radcon, D
831 : include 'formats'
832 :
833 0 : ierr = 0
834 0 : qe3 = qe*qe*qe
835 0 : qe4 = qe3*qe
836 :
837 0 : do i=1,nz
838 0 : unstable(i_SSI,i) = .false.
839 : ! thermal conductivity
840 0 : radcon = 4.D0*crad*clight*T(i)*T(i)*T(i)/(3.D0*opacity(i)*rho(i)) ! erg / (K cm sec)
841 : ! Prandtl-number according to Tassoul
842 0 : dynvisc = visc(i)*rho(i)
843 0 : Prandtl = dynvisc*Cv(i)/radcon
844 0 : Ris_1(i) = 0.125D0*R_crit*Prandtl*Ri_T(i)
845 0 : if (Ris_1(i) <= Ri_crit) then
846 0 : Ris_2(i) = Ri_mu(i)
847 0 : if (Ris_2(i) <= Ri_crit) unstable(i_SSI,i) = .true.
848 : else
849 0 : Ris_2(i) = 0
850 : end if
851 : end do
852 :
853 0 : kbot = nz
854 0 : do i = nz-1, 1, -1
855 :
856 0 : if (unstable(i_SSI,i) .and. .not. unstable(i_SSI,i+1)) kbot = i
857 :
858 : if (unstable(i_SSI,i+1) .and. &
859 0 : (i == 1 .or. .not. unstable(i_SSI,i)) .and. kbot > 1) then
860 0 : if (unstable(i_SSI,i)) then
861 : ktop = i
862 : else
863 0 : ktop = i+1
864 : end if
865 :
866 0 : if (ktop >= kbot) cycle
867 :
868 0 : do k = ktop, kbot ! Heger 2000, eqn 31
869 0 : v_ssi(k)=sqrt(visc(k)/R_crit*abs(domega_dlnR(k)))
870 : end do
871 :
872 : H_ssi(kbot) = v_ssi(kbot)*(r(kbot-1) - r(kbot))/ &
873 0 : max(1d-99,abs(v_ssi(kbot-1) - v_ssi(kbot)))
874 0 : do k = kbot-1, ktop+1, -1
875 : H_ssi(k) = v_ssi(k)*dr(k)/ &
876 0 : max(1d-99,0.5d0*abs(v_ssi(k+1) - v_ssi(k-1)))
877 : end do
878 : H_ssi(ktop) = v_ssi(ktop)*(r(ktop) - r(ktop+1))/ &
879 0 : max(1d-99,abs(v_ssi(ktop) - v_ssi(ktop+1)))
880 :
881 0 : do k = ktop, kbot ! Heger 2000, eqn 34
882 0 : H_ssi(k) = min(H_ssi(k),scale_height(k))
883 0 : v_ssi(k) = min(v_ssi(k),csound(k))
884 : D = H_ssi(k)*v_ssi(k)* &
885 0 : pow2(1d0-max(0d0,max(Ris_1(k),Ris_2(k))/Ri_crit))
886 0 : s% D_SSI(k) = min(D, scale_height(k)*csound(k))
887 :
888 0 : if (s% D_SSI(k) > 1d100) then ! bug
889 0 : write(*,2) 's% D_SSI(k)', k, s% D_SSI(k)
890 0 : write(*,2) 'D', k, D
891 0 : write(*,2) 'scale_height(k)', k, scale_height(k)
892 0 : write(*,2) 'csound(k)', k, csound(k)
893 0 : write(*,2) 'H_ssi(k)', k, H_ssi(k)
894 0 : write(*,2) 'v_ssi(k)', k, v_ssi(k)
895 0 : write(*,2) 'Ris_1(k)', k, Ris_1(k)
896 0 : write(*,2) 'Ris_2(k)', k, Ris_2(k)
897 0 : write(*,2) 'dr(k)', k, dr(k)
898 0 : write(*,2) 'visc(k)', k, visc(k)
899 0 : write(*,2) 'domega_dlnR(k)', k, domega_dlnR(k)
900 0 : call mesa_error(__FILE__,__LINE__,'set_D_SSI')
901 : end if
902 :
903 : end do
904 :
905 : end if
906 :
907 : end do
908 :
909 0 : end subroutine set_D_SSI
910 :
911 :
912 0 : subroutine set_D_ES(ierr)
913 : integer, intent(out) :: ierr
914 : integer :: i, k, kbot, ktop
915 0 : real(dp) :: instability_height, D, v
916 : include 'formats'
917 0 : ierr = 0
918 :
919 0 : do i = 1, nz
920 0 : v = abs(ve0(i)) - abs(ve_mu(i)) ! heger 2000, eqn 38
921 0 : if (v > 0) then
922 0 : unstable(i_ES,i) = .true.
923 0 : v_es(i) = v
924 : else
925 0 : v_es(i) = 0
926 : end if
927 : end do
928 :
929 0 : kbot = nz
930 0 : do i = nz-1, 1, -1
931 :
932 0 : if (unstable(i_ES,i) .and. .not. unstable(i_ES,i+1)) kbot = i
933 :
934 : if (unstable(i_ES,i+1) .and. &
935 0 : (i == 1 .or. .not. unstable(i_ES,i)).and. kbot > 1) then
936 0 : if (unstable(i_ES,i)) then
937 : ktop = i
938 : else
939 0 : ktop = i+1
940 : end if
941 :
942 0 : if (ktop >= kbot) cycle
943 :
944 0 : instability_height = r(ktop) - r(kbot)
945 :
946 : ! heger 2000, eqn 39
947 : H_es(kbot) = v_es(kbot)*(r(kbot-1) - r(kbot))/ &
948 0 : max(1d-99,abs(v_es(kbot-1) - v_es(kbot)))
949 0 : do k = kbot-1, ktop+1, -1
950 : H_es(k) = v_es(k)*dr(k)/ &
951 0 : max(1d-99,0.5d0*abs(v_es(k+1) - v_es(k-1)))
952 : end do
953 : H_es(ktop) = v_es(ktop)*(r(ktop) - r(ktop+1))/ &
954 0 : max(1d-99,abs(v_es(ktop) - v_es(ktop+1)))
955 :
956 0 : do k = ktop, kbot
957 0 : H_es(k) = min(instability_height, H_es(k), scale_height(k))
958 0 : v_es(k) = min(v_es(k), csound(k))
959 0 : D = H_es(k)*v_es(k)
960 0 : s% D_ES(k) = min(D, scale_height(k)*csound(k))
961 : end do
962 :
963 : end if
964 :
965 : end do
966 :
967 0 : end subroutine set_D_ES
968 :
969 :
970 0 : subroutine set_D_GSF(ierr)
971 : integer, intent(out) :: ierr
972 : integer :: i, k, kbot, ktop
973 0 : real(dp) :: instability_height, D, v, v_diff
974 : include 'formats'
975 0 : ierr = 0
976 :
977 0 : do i = 1, nz
978 :
979 : ! heger 2000, eqn 42
980 0 : if(abs(Hj(i))>0d0) then
981 0 : v = ve0(i)*2*H_T(i)*r(i)/(Hj(i)*Hj(i))/(1 + 2*omega(i)*dlnR_domega(i))
982 : else
983 0 : v = 0d0
984 : end if
985 0 : if (is_bad(v)) then
986 0 : write(*,2) 'bad v for GSF', i, v
987 0 : write(*,2) 've0(i)', i, ve0(i)
988 0 : write(*,2) 'H_T(i)', i, H_T(i)
989 0 : write(*,2) 'r(i)', i, r(i)
990 0 : write(*,2) 'Hj(i)', i, Hj(i)
991 0 : write(*,2) 'omega(i)', i, omega(i)
992 0 : write(*,2) 'dlnR_domega(i)', i, dlnR_domega(i)
993 0 : call mesa_error(__FILE__,__LINE__,'set_D_GSF')
994 0 : v = 0
995 : end if
996 0 : v_diff = abs(v) - abs(ve_mu(i)) ! heger 2000, eqn 43
997 0 : if (v_diff > 0) then
998 0 : unstable(i_GSF,i) = .true.
999 0 : v_gsf(i) = v_diff
1000 : else
1001 0 : v_gsf(i) = 0
1002 : end if
1003 :
1004 0 : if (s% model_number == -3 .and. i == -1) then
1005 0 : write(*,2) 've0(i)', i, ve0(i)
1006 0 : write(*,2) 'H_T(i)', i, H_T(i)
1007 0 : write(*,2) 'r(i)', i, r(i)
1008 0 : write(*,2) 'Hj(i)', i, Hj(i)
1009 0 : write(*,2) 'omega(i)', i, omega(i)
1010 0 : write(*,2) 'dlnR_domega(i)', i, dlnR_domega(i)
1011 0 : write(*,2) 'v', i, v
1012 0 : write(*,2) 've_mu(i)', i, ve_mu(i)
1013 0 : write(*,2) 'v_gsf(i)', i, v_gsf(i)
1014 : end if
1015 :
1016 : end do
1017 :
1018 0 : kbot = nz
1019 0 : do i = nz-1, 1, -1
1020 :
1021 0 : if (unstable(i_GSF,i) .and. .not. unstable(i_GSF,i+1)) kbot = i
1022 :
1023 : if (unstable(i_GSF,i+1) .and. &
1024 0 : (i == 1 .or. .not. unstable(i_GSF,i)) .and. kbot > 1) then
1025 0 : if (unstable(i_GSF,i)) then
1026 0 : ktop = i
1027 : else
1028 0 : ktop = i+1
1029 : end if
1030 :
1031 0 : if (ktop >= kbot) cycle
1032 :
1033 0 : instability_height = r(ktop) - r(kbot)
1034 :
1035 : ! heger 2000, eqn 45
1036 : if (kbot == 1) then
1037 : H_gsf(kbot) = v_gsf(kbot)*(s% R_center - r(kbot))/ &
1038 : max(1d-99,abs(v_gsf(kbot)))
1039 : else
1040 : H_gsf(kbot) = v_gsf(kbot)*(r(kbot-1) - r(kbot))/ &
1041 0 : max(1d-99,abs(v_gsf(kbot-1) - v_gsf(kbot)))
1042 : end if
1043 :
1044 : if (kbot == -1) then
1045 : write(*,2) 'r(kbot)', kbot, r(kbot)
1046 : write(*,2) 'r(kbot-1)', kbot-1, r(kbot-1)
1047 : write(*,2) 'v_gsf(kbot)', kbot, v_gsf(kbot)
1048 : write(*,2) 'v_gsf(kbot-1)', kbot-1, v_gsf(kbot-1)
1049 : write(*,2) 'H_gsf(kbot)', kbot, H_gsf(kbot)
1050 : end if
1051 :
1052 0 : do k = max(2,kbot-1), min(nz-1,ktop+1), -1
1053 : H_gsf(k) = v_gsf(k)*dr(k)/ &
1054 0 : max(1d-99,0.5d0*abs(v_gsf(k+1) - v_gsf(k-1)))
1055 0 : if (k == -1) then
1056 0 : write(*,2) 'dr(k)', k, dr(k)
1057 0 : write(*,2) 'v_gsf(k-1)', k-1, v_gsf(k-1)
1058 0 : write(*,2) 'v_gsf(k)', k, v_gsf(k)
1059 0 : write(*,2) 'v_gsf(k+1)', k+1, v_gsf(k+1)
1060 0 : write(*,2) 'H_gsf(k)', k, H_gsf(k)
1061 : end if
1062 : end do
1063 :
1064 0 : if (ktop == nz) then
1065 : H_gsf(ktop) = v_gsf(ktop)*(r(ktop) - s% R_center)/ &
1066 0 : max(1d-99,abs(v_gsf(ktop)))
1067 : else
1068 : H_gsf(ktop) = v_gsf(ktop)*(r(ktop) - r(ktop+1))/ &
1069 0 : max(1d-99,abs(v_gsf(ktop) - v_gsf(ktop+1)))
1070 : end if
1071 :
1072 0 : if (ktop == -1) then
1073 0 : write(*,2) 'r(ktop)', ktop, r(ktop)
1074 0 : write(*,2) 'r(ktop+1)', ktop+1, r(ktop+1)
1075 0 : write(*,2) 'v_gsf(ktop)', ktop, v_gsf(ktop)
1076 0 : write(*,2) 'v_gsf(ktop+1)', ktop+1, v_gsf(ktop+1)
1077 0 : write(*,2) 'H_gsf(ktop)', ktop, H_gsf(ktop)
1078 : end if
1079 :
1080 0 : do k = ktop, kbot
1081 0 : H_gsf(k) = min(instability_height, H_gsf(k), scale_height(k))
1082 0 : v_gsf(k) = min(v_gsf(k), csound(k))
1083 0 : D = H_gsf(k)*v_gsf(k)
1084 0 : s% D_GSF(k) = min(D, scale_height(k)*csound(k))
1085 : end do
1086 :
1087 : end if
1088 :
1089 : end do
1090 :
1091 0 : if (s% model_number == -1) then
1092 0 : k = -1
1093 0 : write(*,2) 've0(k)', k, ve0(k)
1094 0 : write(*,2) 'H_T(k)', k, H_T(k)
1095 0 : write(*,2) 'r(k)', k, r(k)
1096 0 : write(*,2) 'Hj(k)', k, Hj(k)
1097 0 : write(*,2) 'omega(k)', k, omega(k)
1098 0 : write(*,2) 'dlnR_domega(k)', k, dlnR_domega(k)
1099 0 : write(*,2) 'dr(k)', k, dr(k)
1100 0 : write(*,2) 'H_gsf(k)', k, H_gsf(k)
1101 0 : write(*,2) 'v_gsf(k)', k, v_gsf(k)
1102 0 : write(*,2) 'csound(k)', k, csound(k)
1103 0 : write(*,2) 'scale_height(k)', k, scale_height(k)
1104 0 : write(*,2) 's% D_GSF(k)', k, s% D_GSF(k)
1105 : end if
1106 :
1107 0 : end subroutine set_D_GSF
1108 :
1109 :
1110 0 : logical function failed(str, ierr)
1111 : character (len=*), intent(in) :: str
1112 : integer, intent(in) :: ierr
1113 0 : if (ierr == 0) then
1114 : failed = .false.
1115 : return
1116 : end if
1117 0 : write(s% retry_message,*) 'set_rotation_mixing_info failed in call to ' // trim(str)
1118 0 : if (s% report_ierr) write(*, *) s% retry_message
1119 : failed = .true.
1120 : end function failed
1121 :
1122 :
1123 : end subroutine set_rotation_mixing_info
1124 :
1125 :
1126 0 : subroutine smooth_for_rotation(s, v, width, work)
1127 : use star_utils, only: weighed_smoothing
1128 : type (star_info), pointer :: s
1129 : real(dp), dimension(:) :: v, work
1130 : integer :: width
1131 : logical, parameter :: preserve_sign = .false.
1132 0 : if (width <= 0) return
1133 0 : call weighed_smoothing(v, s% nz, width, preserve_sign, work)
1134 0 : end subroutine smooth_for_rotation
1135 :
1136 :
1137 0 : subroutine set_ST(s, &
1138 0 : rho, T, r, L, omega, Cp, abar, zbar, delta, grav, &
1139 0 : N2, N2_mu, opacity, kap_cond, scale_height, &
1140 : ierr)
1141 : ! with modifications by S.-C. Yoon, July 2003
1142 : type (star_info), pointer :: s
1143 : real(dp), dimension(:) :: & ! allocated temporary storage
1144 : rho, T, r, L, omega, Cp, abar, zbar, delta, grav, &
1145 : N2, N2_mu, opacity, kap_cond, scale_height
1146 : integer, intent(out) :: ierr
1147 :
1148 : integer :: nz, k, kk
1149 0 : real(dp) :: xmagfmu, xmagft, xmagfdif, xmagfnu, &
1150 0 : xkap, xgamma, xsig, &
1151 0 : xeta, xmagn, xmagnmu, xmagnt, xmagw, xmagdn, xmagtn, xmagrn, xmag4pd, &
1152 0 : dlnomega_dlnr, xmagq, xmager2w, &
1153 0 : xmagnn, xmagwn, xmagq0, xmagwa0, xmags0, xmagbphi0, xmagbr0, xmageta0, &
1154 0 : xmagkr2n, xmagq1, xmagwa1, xmags1a, xmags1b, xmags1, &
1155 0 : xmagbphi1, xmagbr1, xmageta1a, xmageta1b, xmageta1, &
1156 0 : xmagsm, xmagqm, xmagetam, xmagsf, xmags, xmagnu, xmagdif
1157 :
1158 : include 'formats'
1159 :
1160 0 : ierr = 0
1161 0 : nz = s% nz
1162 :
1163 0 : s% D_ST(1:nz) = 0
1164 0 : s% nu_ST(1:nz) = 0
1165 0 : s% dynamo_B_r(1:nz) = 0
1166 0 : s% dynamo_B_phi(1:nz) = 0
1167 :
1168 0 : xmagfmu = 1
1169 0 : xmagft = 1
1170 0 : xmagfdif = 1
1171 0 : xmagfnu = 1
1172 :
1173 0 : xmageta0 = 0; xmageta1 = 0; xmagetam = 0
1174 0 : xmagq0 = 0; xmagq1 = 0; xmagqm = 0
1175 0 : xmags0 = 0; xmags1 = 0; xmagsm = 0;
1176 :
1177 0 : do k = 2, nz-1
1178 :
1179 : xkap = 16d0*boltz_sigma*T(k)*T(k)*T(k)/ &
1180 0 : (3d0*opacity(k)*rho(k)*rho(k)*Cp(k)) ! thermal diffusivity
1181 :
1182 0 : xsig = calc_sige(abar(k), zbar(k), rho(k), T(k), Cp(k), kap_cond(k), opacity(k))
1183 0 : xeta = calc_eta(xsig) ! magnetic diffusivity
1184 :
1185 0 : xmagn = N2(k)
1186 0 : xmagnmu = N2_mu(k)
1187 0 : xmagnt = xmagn - xmagnmu ! N2_T
1188 0 : xmagw = abs(omega(k))
1189 0 : xmagdn = rho(k)
1190 0 : xmagtn = T(k)
1191 0 : xmagrn = r(k)
1192 0 : xmag4pd = sqrt(pi4*xmagdn)
1193 :
1194 : dlnomega_dlnr = &
1195 0 : (omega(k-1) - omega(k+1))/(r(k-1) - r(k+1))*(r(k)/omega(k))
1196 0 : xmagq = max(1d-30,min(1d30,abs(dlnomega_dlnr))) ! shear
1197 0 : xmager2w = xeta/(r(k)*r(k)*abs(omega(k)))
1198 :
1199 : ! magnetic quantities
1200 0 : if (xmagnmu > 0.0D0) then
1201 0 : xmagnn = xmagnmu*xmagfmu ! N^2
1202 0 : xmagwn = xmagw/sqrt(xmagnn) ! omega/N
1203 0 : xmagq0 = pow(xmagwn,-1.5D0)*pow(xmager2w,0.25D0) ! q0
1204 0 : xmagwa0 = xmagq*xmagwn*xmagw ! omega_A
1205 0 : xmags0 = xmagdn*xmagw*xmagw*xmagrn*xmagrn*pow3(xmagq)*pow4(xmagwn) ! S_0
1206 0 : if (is_bad(xmags0)) then
1207 0 : do kk=1,s% nz
1208 0 : write(*,2) 'omega(kk)', kk, omega(kk)
1209 : end do
1210 0 : write(*,2) 'dlnomega_dlnr', k, dlnomega_dlnr
1211 0 : write(*,2) 'xmagfmu', k, xmagfmu
1212 0 : write(*,2) 'xmagnmu', k, xmagnmu
1213 0 : write(*,2) 'xmagnn', k, xmagnn
1214 0 : write(*,2) 'xmagwn', k, xmagwn
1215 0 : write(*,2) 'xmagq', k, xmagq
1216 0 : write(*,2) 'xmagrn', k, xmagrn
1217 0 : write(*,2) 'xmagw', k, xmagw
1218 0 : write(*,2) 'xmagdn', k, xmagdn
1219 0 : write(*,2) 'xmags0', k, xmags0
1220 0 : call mesa_error(__FILE__,__LINE__,'set_ST')
1221 : end if
1222 0 : xmagbphi0 = xmagwa0*xmag4pd*xmagrn ! B_\phi
1223 0 : xmagbr0 = xmagbphi0*xmagq*xmagwn*xmagwn ! B_r
1224 0 : xmageta0 = pow4(xmagq)*pow6(xmagwn)*xmagrn*xmagrn*xmagw ! eta_e
1225 : end if
1226 :
1227 0 : if (xmagnt > 0.0D0) then
1228 0 : xmagnn = xmagnt*xmagft ! N^2
1229 0 : xmagwn = xmagw/sqrt(xmagnn) ! omega/N
1230 0 : xmagkr2n = xkap/(xmagrn*xmagrn*sqrt(xmagnn)) ! kappa/(r^2 N)
1231 :
1232 : !xmagq1 = pow(xmager2w*xmagwn*pow3(xeta/xkap)/pow7(xmagwn),0.25D0) ! q_1
1233 0 : if(xmagwn > 1d-42) then ! fix from rob
1234 0 : xmagq1 = pow(xmager2w*xmagwn*pow3(xeta/xkap)/pow7(xmagwn),0.25D0) ! q_1
1235 : else
1236 : xmagq1 = 0d0
1237 : end if
1238 :
1239 0 : xmagwa1 = sqrt(xmagq)*xmagw*pow(xmagwn*xmagkr2n,0.125D0) ! \omega_A
1240 0 : xmags1a = xmagdn*pow2(xmagw*xmagrn)*xmagq*sqrt(xmagwn*xmagkr2n) ! S_1a
1241 0 : xmags1b = xmagdn*pow2(xmagw)*pow2(xmagrn)*pow3(xmagq)*pow4(xmagwn) ! S_1b
1242 0 : xmags1 = max(xmags1a,xmags1b) ! S_1
1243 0 : xmagbphi1 = xmagwa1*xmag4pd*xmagrn ! B_\phi
1244 0 : xmagbr1 = xmagbphi1*pow(xmagwn*xmagkr2n,0.25D0) ! B_r
1245 0 : xmageta1a = xmagrn*xmagrn*xmagw*xmagq*pow(xmagwn*xmagkr2n,0.75D0) ! eta_e1a
1246 0 : xmageta1b = pow4(xmagq)*pow6(xmagwn)*pow2(xmagrn)*xmagw ! eta_e1b
1247 0 : xmageta1 = max(xmageta1a,xmageta1b) ! eta_e
1248 : end if
1249 :
1250 0 : if ((xmagnt > 0.D0) .and. (xmagnmu > 0.D0) .and. xmags0+xmags1 /= 0d0) then
1251 0 : xmagsm=xmags0*xmags1/(xmags0+xmags1) ! S_m
1252 0 : if (is_bad(xmagsm)) then
1253 0 : write(*,2) 'xmags0', k, xmags0
1254 0 : write(*,2) 'xmags1', k, xmags1
1255 0 : write(*,2) 'xmags0+xmags1', k, xmags0+xmags1
1256 0 : write(*,2) 'xmags1', k, xmags1
1257 0 : write(*,2) 'xmags0', k, xmags0
1258 0 : write(*,2) 'xmagsm', k, xmagsm
1259 0 : call mesa_error(__FILE__,__LINE__,'set_ST')
1260 : end if
1261 0 : xmagqm=xmagq0+xmagq1 ! q_m
1262 0 : xmagetam=xmageta0*xmageta1/(xmageta0+xmageta1) ! eta_m
1263 0 : else if ((xmagnt <= 0.D0) .and. (xmagnmu > 0.D0)) then
1264 0 : xmagsm=xmags0
1265 0 : xmagqm=xmagq0
1266 0 : xmagetam=xmageta0
1267 0 : else if ((xmagnt > 0.D0) .and. (xmagnmu <= 0.D0)) then
1268 0 : xmagsm=xmags1
1269 0 : xmagqm=xmagq1
1270 0 : xmagetam=xmageta1
1271 0 : else if ((xmagnt <= 0.D0) .and. (xmagnmu <= 0.D0)) then
1272 0 : xmagsm=0.0D0
1273 0 : xmagqm=0.0D0
1274 0 : xmagetam=0.0D0
1275 : end if
1276 :
1277 0 : if (xmagn < 0.0D0) then
1278 0 : xmagsm=0.0D0
1279 0 : xmagqm=0.0D0
1280 0 : xmagetam=0.0D0
1281 : end if
1282 :
1283 0 : if (s% mixing_type(k) == convective_mixing .and. (xmagn > 0.D0)) then
1284 0 : xmagsm=0.0D0
1285 0 : xmagqm=0.0D0
1286 0 : xmagetam=0.0D0
1287 : end if
1288 :
1289 0 : if (s% mixing_type(k) == overshoot_mixing) then
1290 0 : xmagsm=0.0D0
1291 0 : xmagqm=0.0D0
1292 0 : xmagetam=0.0D0
1293 : end if
1294 :
1295 0 : xmagsf=1.D0-MIN(1.D0,xmagqm/xmagq)
1296 0 : xmags=xmagsf*xmagsm ! S
1297 0 : xmagnu=xmags/(xmagw*xmagq*xmagdn) ! \nu_e
1298 0 : xmagdif=xmagsf*xmagetam ! \kappa_c (molecular diffusion)
1299 : ! (set equal to eta_e and compute "mean")
1300 :
1301 : !..... Special treatment for semiconvective regions:
1302 : !..... Geometric mean between semiconvective "effective viscosity"
1303 : !..... and \nu_e as resulting from \mu-dominated case
1304 : !..... (the T-dominated case is indefined in these regions).
1305 : !..... Assume flux is dominated by convective flux and given by
1306 : !..... the luminosity of the star at this point.
1307 : !..... From Kippenhahn&Weigert, Eqs. (7.6) and (7.7), the convective
1308 : !..... velocity can be computed. Now assume \nu=1/3 l_mix v and assume
1309 : !..... that l_mix = H_p.
1310 0 : if ((xmagnt <= 0.D0) .AND. (xmagnmu > 0.D0) .AND. (xmagn > 0.D0)) &
1311 : xmagnu = &
1312 : sqrt(xmagnu*scale_height(k)*(one_third)* &
1313 : pow(grav(k)*delta(k)*scale_height(k)*MAX(0.D0,L(k))/ &
1314 0 : (64.D0*pi*xmagdn*Cp(k)*xmagtn*xmagrn*xmagrn),one_third))
1315 :
1316 0 : s% D_ST(k) = min(xmagdif,xmagnu)*xmagfdif
1317 0 : s% nu_ST(k) = xmagnu*xmagfnu
1318 :
1319 0 : if (is_bad(s% D_ST(k))) then
1320 0 : write(*,4) 'set_ST mixing_type', k, s% model_number, s% mixing_type(k)
1321 0 : write(*,3) 'set_ST xeta', k, s% model_number, xeta
1322 0 : write(*,3) 'set_ST xsig', k, s% model_number, xsig
1323 0 : write(*,3) 'set_ST dlnomega_dlnr', k, s% model_number, dlnomega_dlnr
1324 0 : write(*,3) 'set_ST xmager2w', k, s% model_number, xmager2w
1325 0 : write(*,3) 'set_ST xkap', k, s% model_number, xkap
1326 0 : write(*,3) 'set_ST xgamma', k, s% model_number, xgamma
1327 0 : write(*,3) 'set_ST xmagq', k, s% model_number, xmagq
1328 0 : write(*,3) 'set_ST xmagqm', k, s% model_number, xmagqm
1329 0 : write(*,3) 'set_ST xmagsf', k, s% model_number, xmagsf
1330 0 : write(*,3) 'set_ST xmagsm', k, s% model_number, xmagsm
1331 0 : write(*,3) 'set_ST xmags', k, s% model_number, xmags
1332 0 : write(*,3) 'set_ST xmagnu', k, s% model_number, xmagnu
1333 0 : write(*,3) 'set_ST xmagfdif', k, s% model_number, xmagfdif
1334 0 : write(*,3) 'set_ST D_ST(k)', k, s% model_number, s% D_ST(k)
1335 0 : call mesa_error(__FILE__,__LINE__,'set_ST')
1336 : end if
1337 :
1338 : end do
1339 :
1340 0 : s% D_ST(1) = s% D_ST(2)
1341 0 : s% D_ST(nz) = s% D_ST(nz-1)
1342 :
1343 0 : s% nu_ST(1) = s% nu_ST(2)
1344 0 : s% nu_ST(nz) = s% nu_ST(nz-1)
1345 :
1346 0 : end subroutine set_ST
1347 :
1348 :
1349 0 : subroutine zero_if_convective(nz, mixing_type, D_mix, dc)
1350 : integer, intent(in) :: nz
1351 : integer, dimension(:) :: mixing_type
1352 : real(dp), dimension(:) :: D_mix, dc
1353 : integer :: k
1354 0 : do k=1,nz
1355 0 : if (mixing_type(k) == convective_mixing) dc(k) = 0
1356 : end do
1357 0 : end subroutine zero_if_convective
1358 :
1359 :
1360 0 : subroutine zero_if_tiny(s, dc)
1361 : type (star_info), pointer :: s
1362 : real(dp), dimension(:) :: dc
1363 : integer :: k
1364 0 : real(dp) :: tiny
1365 0 : tiny = s% clip_D_limit
1366 0 : do k=1,s% nz
1367 0 : if (dc(k) < tiny) dc(k) = 0
1368 : end do
1369 0 : end subroutine zero_if_tiny
1370 :
1371 : end module rotation_mix_info
|
