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_step
21 :
22 : use star_private_def
23 : use overshoot_utils
24 :
25 : implicit none
26 :
27 : private
28 : public :: eval_overshoot_step
29 :
30 : contains
31 :
32 0 : subroutine eval_overshoot_step (s, i, j, k_a, k_b, D, vc, ierr)
33 :
34 : type(star_info), pointer :: s
35 : integer, intent(in) :: i
36 : integer, intent(in) :: j
37 : integer, intent(out) :: k_a
38 : integer, intent(out) :: k_b
39 : real(dp), intent(out) :: D(:)
40 : real(dp), intent(out) :: vc(:)
41 : integer, intent(out) :: ierr
42 :
43 : real(dp) :: f, f2
44 : real(dp) :: f0
45 : real(dp) :: D0
46 : real(dp) :: Delta0
47 : real(dp) :: w
48 : real(dp) :: factor
49 : real(dp) :: Hp_cb
50 : integer :: k_ob
51 : real(dp) :: r_ob
52 : real(dp) :: D_ob
53 : real(dp) :: vc_ob
54 : logical :: outward
55 : integer :: dk
56 : integer :: k
57 : real(dp) :: r
58 : real(dp) :: dr
59 :
60 : ! Evaluate the overshoot diffusion coefficient D(k_a:k_b) and
61 : ! mixing velocity vc(k_a:k_b) at the i'th convective boundary,
62 : ! using the j'th set of overshoot parameters. The overshoot
63 : ! follows a simple step scheme, with an optional exponential tail.
64 :
65 0 : ierr = 0
66 :
67 : ! Extract parameters
68 :
69 0 : f = s%overshoot_f(j)
70 0 : f2 = 0._dp
71 0 : f0 = s%overshoot_f0(j)
72 :
73 0 : D0 = s%overshoot_D0(j)
74 0 : Delta0 = s%overshoot_Delta0(j)
75 :
76 0 : if (s%overshoot_scheme(j) == 'step+exponential') then
77 0 : f2 = s%overshoot_f2(j)
78 0 : if (f <= 0._dp .OR. f0 <= 0._dp .OR. f2 <= 0._dp) then
79 0 : write(*,*) 'ERROR: for step+exp overshooting, must set f, f2 and f0 > 0'
80 0 : write(*,*) 'see description of overshooting in star/defaults/control.defaults'
81 0 : ierr = -1
82 0 : return
83 : end if
84 0 : else if (f <= 0._dp .OR. f0 <= 0._dp) then
85 0 : write(*,*) 'ERROR: for step overshooting, must set f and f0 > 0'
86 0 : write(*,*) 'see description of overshooting in star/defaults/control.defaults'
87 0 : ierr = -1
88 0 : return
89 : end if
90 :
91 : ! Apply mass limits
92 :
93 0 : if (s%star_mass < s%overshoot_mass_full_on(j)) then
94 0 : if (s%star_mass > s%overshoot_mass_full_off(j)) then
95 : w = (s%star_mass - s%overshoot_mass_full_off(j)) / &
96 0 : (s%overshoot_mass_full_on(j) - s%overshoot_mass_full_off(j))
97 0 : factor = 0.5_dp*(1._dp - cospi(w))
98 0 : f = f*factor
99 0 : f0 = f0*factor
100 : else
101 0 : f = 0._dp
102 0 : f0 = 0._dp
103 : end if
104 : end if
105 :
106 : ! Evaluate convective boundary (_cb) parameters
107 :
108 0 : call eval_conv_bdy_Hp(s, i, Hp_cb, ierr)
109 0 : if (ierr /= 0) return
110 :
111 : ! Evaluate overshoot boundary (_ob) parameters
112 :
113 0 : call eval_over_bdy_params(s, i, f0, k_ob, r_ob, D_ob, vc_ob, ierr)
114 0 : if (ierr /= 0) return
115 :
116 : ! Loop over cell faces, adding overshoot until D <= overshoot_D_min
117 :
118 0 : outward = s%top_conv_bdy(i)
119 :
120 0 : if (outward) then
121 0 : k_a = k_ob
122 0 : k_b = 1
123 0 : dk = -1
124 : else
125 0 : k_a = k_ob+1
126 0 : k_b = s%nz
127 0 : dk = 1
128 : end if
129 :
130 0 : face_loop : do k = k_a, k_b, dk
131 :
132 0 : r = s%r(k)
133 :
134 0 : if (outward) then
135 0 : dr = r - r_ob
136 : else
137 0 : dr = r_ob - r
138 : end if
139 :
140 0 : if (dr < f*Hp_cb) then
141 : factor = 1._dp
142 0 : else if (f2 > 0._dp) then
143 0 : factor = exp(-2._dp*(dr - f*Hp_cb)/(f2*Hp_cb))
144 : else
145 : factor = 0._dp
146 : end if
147 :
148 : ! Store the diffusion coefficient and velocity
149 :
150 0 : D(k) = (D0 + Delta0*D_ob)*factor
151 0 : if(D_ob /= 0d0) then
152 0 : vc(k) = (D0/D_ob + Delta0)*vc_ob*factor
153 : else
154 0 : vc(k) = 0d0
155 : end if
156 :
157 : ! Check for early overshoot completion
158 :
159 0 : if (D(k) < s%overshoot_D_min) then
160 0 : k_b = k
161 0 : exit face_loop
162 : end if
163 :
164 : end do face_loop
165 :
166 0 : ierr = 0
167 :
168 0 : return
169 :
170 : end subroutine eval_overshoot_step
171 :
172 : end module overshoot_step
|
