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 overshoot_utils
21 :
22 : use num_lib
23 : use star_private_def
24 :
25 : implicit none
26 :
27 : private
28 : public :: eval_conv_bdy_k
29 : public :: eval_conv_bdy_r
30 : public :: eval_conv_bdy_Hp
31 : public :: eval_over_bdy_params
32 :
33 : contains
34 :
35 0 : subroutine eval_conv_bdy_k (s, i, k, ierr)
36 :
37 : type(star_info), pointer :: s
38 : integer, intent(in) :: i
39 : integer, intent(out) :: k
40 : integer, intent(out) :: ierr
41 :
42 : ! Evaluate the index k of the cell containing the i'th convective
43 : ! boundary
44 :
45 0 : ierr = 0
46 :
47 0 : if (s%top_conv_bdy(i)) then
48 0 : k = s%conv_bdy_loc(i)
49 : else
50 0 : k = s%conv_bdy_loc(i) - 1
51 : end if
52 :
53 0 : if (k >= s%nz .OR. k < 1) then
54 0 : write(*,*) 'Invalid cell for convective boundary: i, k, nz=', i, k, s%nz
55 0 : ierr = -1
56 0 : return
57 : end if
58 :
59 : return
60 :
61 : end subroutine eval_conv_bdy_k
62 :
63 :
64 0 : subroutine eval_conv_bdy_r (s, i, r, ierr)
65 :
66 : type(star_info), pointer :: s
67 : integer, intent(in) :: i
68 : real(dp), intent(out) :: r
69 : integer, intent(out) :: ierr
70 :
71 : integer :: k
72 : real(dp) :: w
73 :
74 : ! Evaluate the radius r at the i'th convective boundary
75 :
76 : ! Find the convective boundary cell
77 :
78 : ierr = 0
79 :
80 0 : call eval_conv_bdy_k(s, i, k, ierr)
81 0 : if (ierr /= 0) return
82 :
83 : ! Interpolate r based on the fact that r^3 varies linearly with q
84 : ! across the (constant-density) cell
85 :
86 0 : w = s%cz_bdy_dq(k)/s%dq(k)
87 :
88 0 : if (w < 0._dp .OR. w > 1._dp) then
89 0 : write(*,*) 'Invalid weight for convective boundary: i, w=', i, w
90 0 : ierr = -1
91 0 : return
92 : end if
93 :
94 : associate (k_o => k, &
95 : k_i => k+1)
96 :
97 : r = pow(( w)*s%r(k_i)*s%r(k_i)*s%r(k_i) + &
98 0 : (1._dp-w)*s%r(k_o)*s%r(k_o)*s%r(k_o), 1._dp/3._dp)
99 :
100 : end associate
101 :
102 0 : return
103 :
104 0 : end subroutine eval_conv_bdy_r
105 :
106 :
107 0 : subroutine eval_conv_bdy_Hp (s, i, Hp, ierr)
108 :
109 : type(star_info), pointer :: s
110 : integer, intent(in) :: i
111 : real(dp), intent(out) :: Hp
112 : integer, intent(out) :: ierr
113 :
114 : integer :: k
115 : real(dp) :: r
116 0 : real(dp) :: x0
117 0 : real(dp) :: x1
118 0 : real(dp) :: x2
119 0 : real(dp) :: x
120 0 : real(dp) :: a0
121 0 : real(dp) :: a1
122 0 : real(dp) :: a2
123 0 : real(dp) :: P
124 0 : real(dp) :: rho
125 0 : real(dp) :: r_top
126 0 : real(dp) :: r_bot
127 0 : real(dp) :: dr
128 :
129 : ! Evaluate the pressure scale height Hp at the i'th convective boundary
130 :
131 : ! Find the convective boundary cell
132 :
133 : ierr = 0
134 :
135 0 : call eval_conv_bdy_k(s, i, k, ierr)
136 0 : if (ierr /= 0) return
137 :
138 : ! Evaluate the radius at the convective boundary
139 :
140 0 : call eval_conv_bdy_r(s, i, r, ierr)
141 0 : if (ierr /= 0) return
142 :
143 : ! Interpolate the pressure and density at the boundary, using a
144 : ! quadratic fit across the boundary cell and its neighbors (the
145 : ! x's are fractional mass distances from the outer edge of cell
146 : ! k-1); then, evaluate the pressure scale height
147 :
148 : associate (k_o => k-1, &
149 : k_m => k, &
150 : k_i => k+1)
151 :
152 0 : x0 = s%dq(k_o)/2._dp
153 0 : x1 = s%dq(k_o) + s%dq(k_m)/2._dp
154 0 : x2 = s%dq(k_o) + s%dq(k_m) + s%dq(k_i)/2._dp
155 :
156 0 : x = s%dq(k_o) + s%cz_bdy_dq(k)
157 :
158 0 : call two_piece_linear_coeffs(x, x0, x1, x2, a0, a1, a2, ierr)
159 0 : if (ierr /= 0) return
160 :
161 0 : P = exp(a0*s%lnPeos(k_o) + a1*s%lnPeos(k_m) + a2*s%lnPeos(k_i))
162 0 : rho = exp(a0*s%lnd(k_o) + a1*s%lnd(k_m) + a2*s%lnd(k_i))
163 :
164 : ! Evaluate the pressure scale height
165 :
166 : Hp = P/(rho*s%cgrav(k_m)* &
167 0 : (s%M_center + s%xmstar*s%conv_bdy_q(i))/(r*r))
168 :
169 : end associate
170 :
171 : ! (Possibly) limit the scale height using the size of the
172 : ! convection zone
173 :
174 0 : if (s%limit_overshoot_Hp_using_size_of_convection_zone) then
175 :
176 : ! Determine the radial extent of the convection zone (note that
177 : ! r_top/r_bot don't coincide exactly with the r calculated
178 : ! above)
179 :
180 0 : if (s%top_conv_bdy(i)) then
181 :
182 0 : if (i == 1) then
183 0 : r_bot = s%R_center
184 : else
185 0 : if (s%top_conv_bdy(i-1)) then
186 0 : write(*,*) 'Double top boundary in overshoot; i=', i
187 0 : ierr = -1
188 0 : return
189 : end if
190 0 : r_bot = s%r(s%conv_bdy_loc(i-1))
191 : end if
192 :
193 0 : r_top = s%r(k)
194 :
195 : else
196 :
197 0 : r_bot = s%r(k+1)
198 :
199 0 : if (i == s%num_conv_boundaries) then
200 0 : r_top = s%r(1)
201 : else
202 0 : if (.NOT. s%top_conv_bdy(i+1)) then
203 0 : write(*,*) 'Double bottom boundary in overshoot; i=', i
204 0 : ierr = -1
205 0 : return
206 : end if
207 0 : r_top = s%r(s%conv_bdy_loc(i+1))
208 : end if
209 :
210 : end if
211 :
212 0 : dr = r_top - r_bot
213 :
214 : ! Apply the limit
215 :
216 0 : if (s%overshoot_alpha > 0d0) then
217 0 : if (s%overshoot_alpha*Hp > dr) Hp = dr/s%overshoot_alpha
218 : else
219 0 : if (s%alpha_mlt(k)*Hp > dr) Hp = dr/s%mixing_length_alpha
220 : end if
221 :
222 : end if
223 :
224 : return
225 :
226 0 : end subroutine eval_conv_bdy_Hp
227 :
228 :
229 0 : subroutine eval_over_bdy_params (s, i, f0, k, r, D, vc, ierr)
230 :
231 : type(star_info), pointer :: s
232 : integer, intent(in) :: i
233 : real(dp), intent(in) :: f0
234 : integer, intent(out) :: k
235 : real(dp), intent(out) :: r
236 : real(dp), intent(out) :: D
237 : real(dp), intent(out) :: vc
238 : integer, intent(out) :: ierr
239 :
240 : integer :: k_cb
241 : real(dp) :: r_cb
242 0 : real(dp) :: Hp_cb
243 0 : real(dp) :: w
244 0 : real(dp) :: lambda
245 :
246 : ! Evaluate parameters (cell index k, radius r, diffusion
247 : ! coefficients D and cdc) for the overshoot boundary associated
248 : ! with the i'th convective boundary
249 :
250 : ! Find the convective boundary cell
251 :
252 : ierr = 0
253 :
254 0 : call eval_conv_bdy_k(s, i, k_cb, ierr)
255 0 : if (ierr /= 0) return
256 :
257 : ! Evaluate the radius at the convective boundary
258 :
259 0 : call eval_conv_bdy_r(s, i, r_cb, ierr)
260 0 : if (ierr /= 0) return
261 :
262 : ! Evaluate the pressure scale height at the convective boundary
263 :
264 0 : call eval_conv_bdy_Hp(s, i, Hp_cb, ierr)
265 0 : if (ierr /= 0) return
266 :
267 : ! Search for the overshoot boundary cell
268 :
269 0 : ierr = 0
270 :
271 0 : if (s%top_conv_bdy(i)) then
272 :
273 : ! Overshooting outward -- search inward
274 :
275 0 : r = r_cb - f0*Hp_cb
276 :
277 0 : if (r <= s%r(s%nz)) then
278 :
279 0 : r = s%r(s%nz)
280 0 : k = s%nz - 1
281 :
282 : else
283 :
284 0 : search_in_loop: do k = k_cb, s%nz-1
285 0 : if (s%r(k+1) <= r) exit search_in_loop
286 : end do search_in_loop
287 :
288 : end if
289 :
290 : else
291 :
292 : ! Overshooting inward -- search outward
293 :
294 0 : r = r_cb + f0*Hp_cb
295 :
296 0 : if (r >= s%r(1)) then
297 :
298 0 : r = s%r(1)
299 0 : k = 1
300 :
301 : else
302 :
303 0 : search_out_loop : do k = k_cb, 1, -1
304 0 : if (s%r(k) > r) exit search_out_loop
305 : end do search_out_loop
306 :
307 : end if
308 :
309 : end if
310 :
311 0 : if (.NOT. (s%r(k+1) <= r .AND. s%r(k) >= r)) then
312 0 : write(*,*) 'r_ob not correctly bracketed: r(k+1), r, r(k)=', s%r(k+1), r, s%r(k)
313 0 : ierr = -1
314 0 : return
315 : end if
316 :
317 : ! Interpolate mixing parameters
318 :
319 : w = (s%r(k)*s%r(k)*s%r(k) - r*r*r)/ &
320 0 : (s%r(k)*s%r(k)*s%r(k) - s%r(k+1)*s%r(k+1)*s%r(k+1))
321 :
322 0 : lambda = (1._dp-w)*s%mlt_mixing_length(k) + w*s%mlt_mixing_length(k+1)
323 :
324 0 : if (s%conv_vel(k) /= 0._dp .AND. s%conv_vel(k+1) /= 0._dp) then
325 :
326 : ! Both faces of cell have non-zero mixing; interpolate vc between faces
327 :
328 0 : vc = (1._dp-w)*s%conv_vel(k) + w*s%conv_vel(k+1)
329 :
330 0 : elseif (s%conv_vel(k) /= 0._dp .AND. s%conv_vel(k+1) == 0._dp) then
331 :
332 : ! Outer face of cell has non-zero mixing; interpolate vc
333 : ! between this face and r_cb, assuming vc = 0 at the latter
334 :
335 0 : if(s%r(k) /= r_cb) then
336 : w = (s%r(k)*s%r(k)*s%r(k) - r*r*r)/ &
337 0 : (s%r(k)*s%r(k)*s%r(k) - r_cb*r_cb*r_cb)
338 : else
339 : w = 0d0
340 : end if
341 :
342 0 : vc = (1._dp-w)*s%conv_vel(k)
343 :
344 0 : elseif (s%conv_vel(k) == 0._dp .AND. s%conv_vel(k+1) /= 0._dp) then
345 :
346 : ! Inner face of cell has non-zero mixing; interpolate vc
347 : ! between this face and r_cb, assuming vc = 0 at the latter
348 :
349 0 : if(s%r(k+1) /= r_cb) then
350 : w = (r_cb*r_cb*r_cb - r*r*r)/ &
351 0 : (r_cb*r_cb*r_cb - s%r(k+1)*s%r(k+1)*s%r(k+1))
352 : else
353 : w = 0d0
354 : end if
355 :
356 0 : vc = w*s%conv_vel(k+1)
357 :
358 : else
359 :
360 : ! Neither face of cell has non-zero mixing; return
361 :
362 0 : vc = 0._dp
363 :
364 : end if
365 :
366 : ! Evaluate the diffusion coefficient
367 :
368 0 : D = vc*lambda/3._dp
369 :
370 : ierr = 0
371 :
372 0 : return
373 :
374 : end subroutine eval_over_bdy_params
375 :
376 : end module overshoot_utils
|
