Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2017 Rich Townsend & 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
21 :
22 : use const_def, only: dp, pi4, no_mixing, overshoot_mixing
23 : use num_lib
24 : use star_private_def
25 : use overshoot_utils
26 : use overshoot_exp
27 : use overshoot_step
28 :
29 : implicit none
30 :
31 : private
32 : public :: add_overshooting
33 :
34 : contains
35 :
36 33 : subroutine add_overshooting (s, ierr)
37 :
38 : type(star_info), pointer :: s
39 : integer, intent(out) :: ierr
40 :
41 : logical, parameter :: dbg = .false.
42 :
43 : integer :: i
44 : integer :: j
45 : logical :: is_core
46 : logical :: match_zone_type
47 : logical :: match_zone_loc
48 : logical :: match_bdy_loc
49 : integer :: k_cb
50 : integer :: k_a
51 : integer :: k_b
52 66 : real(dp) :: D(s%nz)
53 66 : real(dp) :: vc(s%nz)
54 : integer :: k
55 : integer :: dk
56 : real(dp) :: rho
57 : real(dp) :: cdc
58 : real(dp) :: r_cb
59 :
60 : include 'formats'
61 :
62 : ! Initialize
63 :
64 33 : ierr = 0
65 :
66 : if (dbg) then
67 : write(*, 3) 'add_overshooting; model, n_conv_bdy=', s%model_number, s%num_conv_boundaries
68 : end if
69 :
70 : ! Keep the convective reference velocities fixed during this loop.
71 : ! Loop over convective boundaries, from center to surface
72 :
73 138 : conv_bdy_loop : do i = 1, s%num_conv_boundaries
74 :
75 : ! Skip this boundary if it's too close to the center
76 :
77 105 : if (s%conv_bdy_q(i) < s%min_overshoot_q) then
78 : if (dbg) then
79 : write(*,*) 'skip since s%conv_bdy_q(i) < min_overshoot_q', i
80 : end if
81 : cycle conv_bdy_loop
82 : end if
83 :
84 : ! Skip this boundary if it's at the surface, since we don't
85 : ! overshoot there
86 :
87 105 : if (s%conv_bdy_loc(i) == 1) then
88 : if (dbg) then
89 : write(*,*) 'skip since s%conv_bdy_loc(i) == 1', i
90 : end if
91 : cycle conv_bdy_loop
92 : end if
93 :
94 : ! Loop over overshoot criteria
95 :
96 1818 : criteria_loop : do j = 1, NUM_OVERSHOOT_PARAM_SETS
97 :
98 1680 : if (s%overshoot_scheme(j) == '') cycle criteria_loop
99 :
100 : ! Check if the criteria match the current boundary
101 :
102 0 : select case (s%overshoot_zone_type(j))
103 : case ('burn_H')
104 0 : match_zone_type = s%burn_h_conv_region(i)
105 : case ('burn_He')
106 0 : match_zone_type = s%burn_he_conv_region(i)
107 : case ('burn_Z')
108 0 : match_zone_type = s%burn_z_conv_region(i)
109 : case ('nonburn')
110 : match_zone_type = .NOT. ( &
111 : s%burn_h_conv_region(i) .OR. &
112 : s%burn_he_conv_region(i) .OR. &
113 0 : s%burn_z_conv_region(i) )
114 : case ('any')
115 0 : match_zone_type = .true.
116 : case default
117 0 : write(*,*) 'Invalid overshoot_zone_type: j, s%overshoot_zone_type(j)=', j, s%overshoot_zone_type(j)
118 0 : ierr = -1
119 0 : return
120 : end select
121 :
122 0 : is_core = (i == 1 .AND. s%R_center == 0d0 .AND. s%top_conv_bdy(i))
123 :
124 : select case (s%overshoot_zone_loc(j))
125 : case ('core')
126 0 : match_zone_loc = is_core
127 : case ('shell')
128 0 : match_zone_loc = .NOT. is_core
129 : case ('any')
130 0 : match_zone_loc = .true.
131 : case default
132 0 : write(*,*) 'Invalid overshoot_zone_loc: j, s%overshoot_zone_loc(j)=', j, s%overshoot_zone_loc(j)
133 0 : ierr = -1
134 0 : return
135 : end select
136 :
137 0 : select case (s%overshoot_bdy_loc(j))
138 : case ('bottom')
139 0 : match_bdy_loc = .NOT. s%top_conv_bdy(i)
140 : case ('top')
141 0 : match_bdy_loc = s%top_conv_bdy(i)
142 : case ('any')
143 0 : match_bdy_loc = .true.
144 : case default
145 0 : write(*,*) 'Invalid overshoot_bdy_loc: j, s%overshoot_bdy_loc(j)=', j, s%overshoot_bdy_loc(j)
146 0 : ierr = -1
147 0 : return
148 : end select
149 :
150 0 : if (.NOT. (match_zone_type .AND. match_zone_loc .AND. match_bdy_loc)) cycle criteria_loop
151 :
152 : if (dbg) then
153 : write(*,*) 'Overshooting at convective boundary: i, j=', i, j
154 : write(*,*) ' s%overshoot_scheme=', TRIM(s%overshoot_scheme(j))
155 : write(*,*) ' s%overshoot_zone_type=', TRIM(s%overshoot_zone_type(j))
156 : write(*,*) ' s%overshoot_zone_loc=', TRIM(s%overshoot_zone_loc(j))
157 : write(*,*) ' s%overshoot_bdy_loc=', TRIM(s%overshoot_bdy_loc(j))
158 : end if
159 :
160 : ! Special-case check for an overshoot scheme of 'none' (this can be used
161 : ! to turn *off* overshoot for specific boundary configurations)
162 :
163 0 : if (s%overshoot_scheme(j) == 'none') exit criteria_loop
164 :
165 : ! Evaluate convective boundary (_cb) parameters
166 :
167 0 : call eval_conv_bdy_k(s, i, k_cb, ierr)
168 0 : if (ierr /= 0) return
169 :
170 : ! Evaluate the overshoot diffusion coefficient and velocity
171 : ! using the appropriate scheme-dependent routine
172 :
173 0 : select case (s%overshoot_scheme(j))
174 : case ('exponential')
175 0 : call eval_overshoot_exp(s, i, j, k_a, k_b, D, vc, ierr)
176 : case ('step', 'step+exponential')
177 0 : call eval_overshoot_step(s, i, j, k_a, k_b, D, vc, ierr)
178 : case ('other')
179 0 : call s% other_overshooting_scheme(s% id, i, j, k_a, k_b, D, vc, ierr)
180 : case default
181 0 : write(*,*) 'Invalid overshoot_scheme:', s%overshoot_scheme(j)
182 0 : ierr = -1
183 : end select
184 :
185 0 : if (ierr /= 0) return
186 :
187 : ! Update the model
188 :
189 0 : if (s%top_conv_bdy(i)) then
190 : dk = -1
191 : else
192 0 : dk = 1
193 : end if
194 :
195 0 : face_loop : do k = k_a, k_b, dk
196 :
197 : ! Check if the overshoot will be stabilized by the stratification
198 :
199 0 : if (s%overshoot_brunt_B_max > 0._dp .and. s% calculate_Brunt_B) then
200 :
201 0 : if (.not. s% calculate_Brunt_N2) &
202 : call mesa_error(__FILE__,__LINE__, &
203 0 : 'add_overshooting: when overshoot_brunt_B_max > 0, must have calculate_Brunt_N2 = .true.')
204 :
205 : ! (NB: we examine B(k+dk) rather than B(k), as the latter
206 : ! would allow the overshoot region to eat into a composition
207 : ! gradient). Special case for k == 1 or k == nz
208 :
209 0 : if (k > 1 .AND. k < s%nz) then
210 0 : if (s%unsmoothed_brunt_B(k+dk) > s%overshoot_brunt_B_max) exit face_loop
211 : else
212 0 : if (s%unsmoothed_brunt_B(k) > s%overshoot_brunt_B_max) exit face_loop
213 : end if
214 :
215 : end if
216 :
217 : ! Check whether D has dropped below the minimum
218 :
219 0 : if (D(k) < s%overshoot_D_min) then
220 :
221 : ! Update conv_bdy_dq to reflect where D drops below the minimum
222 : ! Convective regions can happen to be entirely below s%overshoot_D_min,
223 : ! in which case we ignore this correction.
224 0 : if (s%top_conv_bdy(i)) then
225 0 : if (s%D_mix(k+1) > s%overshoot_D_min) then
226 0 : s%cz_bdy_dq(k) = find0(0._dp, D(k)-s%overshoot_D_min, s%dq(k), s%D_mix(k+1)-s%overshoot_D_min)
227 0 : if (s%cz_bdy_dq(k) < 0._dp .OR. s%cz_bdy_dq(k) > s%dq(k)) then
228 0 : write(*,*) 'k, k_a, k_b', k, k_a, k_b
229 0 : write(*,*) 's%top_conv_bdy(i)=', s%top_conv_bdy(i)
230 0 : write(*,*) 'D(k)', D(k)
231 0 : write(*,*) 's%D_mix(k+1)', s%D_mix(k+1)
232 0 : write(*,*) 's%overshoot_D_min', s%overshoot_D_min
233 0 : write(*,*) 'Invalid location for overshoot boundary: cz_bdy_dq, dq=', s%cz_bdy_dq(k), s%dq(k)
234 0 : ierr = -1
235 0 : return
236 : end if
237 : end if
238 : else
239 0 : if (s%D_mix(k-1) > s%overshoot_D_min) then
240 0 : s%cz_bdy_dq(k-1) = find0(0._dp, s%D_mix(k-1)-s%overshoot_D_min, s%dq(k-1), D(k)-s%overshoot_D_min)
241 0 : if (s%cz_bdy_dq(k-1) < 0._dp .OR. s%cz_bdy_dq(k-1) > s%dq(k-1)) then
242 0 : write(*,*) 'k, k_a, k_b', k, k_a, k_b
243 0 : write(*,*) 's%top_conv_bdy(i)=', s%top_conv_bdy(i)
244 0 : write(*,*) 'D(k)', D(k)
245 0 : write(*,*) 's%D_mix(k-1)', s%D_mix(k-1)
246 0 : write(*,*) 's%overshoot_D_min', s%overshoot_D_min
247 0 : write(*,*) 'Invalid location for overshoot boundary: cz_bdy_dq, dq=', s%cz_bdy_dq(k-1), s%dq(k)
248 0 : ierr = -1
249 0 : return
250 : end if
251 : end if
252 : end if
253 :
254 : exit face_loop
255 :
256 : end if
257 :
258 : ! Revise mixing coefficients
259 :
260 0 : if (k > 1) then
261 : rho = (s%dq(k-1)*s%rho(k) + s%dq(k)*s%rho(k-1))/ &
262 0 : (s%dq(k-1) + s%dq(k))
263 : else
264 0 : rho = s%rho(k)
265 : end if
266 :
267 0 : cdc = (pi4*s%r(k)*s%r(k)*rho)*(pi4*s%r(k)*s%r(k)*rho)*D(k) ! gm^2/sec
268 :
269 0 : call eval_conv_bdy_r(s, i, r_cb, ierr)
270 0 : if (ierr /= 0) then
271 0 : write(*,*) 'error calling eval_conv_bdy_r in overshoot:add_overshooting'
272 0 : return
273 : end if
274 :
275 0 : if (s% r(k)*dk >= r_cb*dk) then
276 :
277 0 : s%cdc(k) = cdc
278 0 : s%D_mix(k) = D(k)
279 :
280 0 : elseif (D(k) > s%D_mix(k)) then
281 :
282 0 : s%cdc(k) = cdc
283 0 : s%D_mix(k) = D(k)
284 0 : s%mixing_type(k) = overshoot_mixing
285 :
286 : end if
287 :
288 : end do face_loop
289 :
290 : ! If we're still in the convection zone, set the rest of the
291 : ! zone to be non-convective
292 :
293 0 : s%cdc(k:k_cb:dk) = 0._dp
294 0 : s%D_mix(k:k_cb:dk) = 0._dp
295 0 : s%mixing_type(k:k_cb:dk) = no_mixing
296 :
297 : ! Finish (we apply at most a single overshoot scheme to each boundary)
298 :
299 105 : exit criteria_loop
300 :
301 : end do criteria_loop
302 :
303 : end do conv_bdy_loop
304 :
305 : ! Perform a sanity check on D_mix
306 :
307 40766 : check_loop : do k = 1, s%nz
308 :
309 40766 : if (is_bad_num(s%D_mix(k))) then
310 :
311 0 : ierr = -1
312 :
313 0 : if (s%report_ierr) then
314 0 : write(*,3) 'mixing_type, D_mix', k, s%mixing_type(k), s%D_mix(k)
315 : end if
316 :
317 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'add_overshooting')
318 :
319 : end if
320 :
321 : end do check_loop
322 :
323 : return
324 :
325 : end subroutine add_overshooting
326 :
327 : end module overshoot
|
