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