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