Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2021 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 phase_separation
21 :
22 : use star_private_def
23 : use const_def
24 : use forum_m, only: hdf5io_t, OPEN_FILE_RO
25 :
26 : implicit none
27 :
28 : private
29 : public :: do_phase_separation
30 :
31 : logical, parameter :: dbg = .false.
32 :
33 : ! offset to higher phase than 0.5 to avoid interference
34 : ! between phase separation mixing and latent heat for Skye.
35 : real(dp), parameter :: eos_phase_boundary = 0.9d0
36 :
37 : contains
38 :
39 0 : subroutine do_phase_separation(s, dt, ierr)
40 : type (star_info), pointer :: s
41 : real(dp), intent(in) :: dt
42 : integer, intent(out) :: ierr
43 :
44 : ! 'CO' or 'ONe' will implement 2-species phase separation, for 'ONe' 22Ne is included
45 0 : if(s% phase_separation_option == 'CO') then
46 0 : call separate_mix_and_heat(s, dt, 'CO', ierr)
47 0 : else if(s% phase_separation_option == 'ONe') then
48 0 : call separate_mix_and_heat(s, dt, 'ONe', ierr)
49 0 : else if(s% phase_separation_option == '3c') then
50 0 : call separate_mix_and_heat(s, dt, '3c', ierr)
51 : else
52 0 : write(*,*) 'invalid phase_separation_option'
53 0 : stop
54 : end if
55 :
56 0 : if(ierr/=0) return
57 :
58 0 : if(s% smooth_phase_separation_heating) then
59 : ! Redistribute energy associated with phase separation evenly through inner half of star.
60 : ! This can help with small timesteps where there may be too much localized heating.
61 0 : call smooth_eps_phase_sep(s, dt, ierr)
62 : end if
63 :
64 : end subroutine do_phase_separation
65 :
66 0 : subroutine separate_mix_and_heat(s, dt, components, ierr)
67 : use chem_def, only: ic12, io16, ine20, ine22
68 : use chem_lib, only: chem_get_iso_id
69 : type (star_info), pointer :: s
70 : real(dp), intent(in) :: dt
71 : character (len=*), intent(in) :: components
72 : integer, intent(out) :: ierr
73 : real(dp) :: XNe20, XNe22, XO, XC, pad
74 : integer :: k, k_bound, kstart, net_ic12, net_io16, net_ine20, net_ine22
75 : logical :: save_Skye_use_ion_offsets
76 :
77 : ! Set phase separation mixing mass negative at beginning of phase separation
78 0 : s% phase_sep_mixing_mass = -1d0
79 0 : s% eps_phase_separation(1:s%nz) = 0d0
80 :
81 0 : if(s% phase(s% nz) < eos_phase_boundary) then
82 0 : s% crystal_core_boundary_mass = 0d0
83 0 : return
84 : end if
85 :
86 0 : net_ic12 = s% net_iso(ic12)
87 0 : net_io16 = s% net_iso(io16)
88 0 : net_ine20 = s% net_iso(ine20)
89 0 : net_ine22 = s% net_iso(ine22)
90 :
91 : ! Find zone of phase transition from liquid to solid
92 0 : k_bound = -1
93 0 : do k = s%nz,1,-1
94 0 : if(s% phase(k-1) <= eos_phase_boundary .and. s% phase(k) > eos_phase_boundary) then
95 : k_bound = k
96 : exit
97 : end if
98 : end do
99 :
100 0 : XC = s% xa(net_ic12,k_bound)
101 0 : XO = s% xa(net_io16,k_bound)
102 0 : XNe20 = s% xa(net_ine20,k_bound)
103 0 : XNe22 = s% xa(net_ine22,k_bound)
104 :
105 : ! Check that we're still in C/O or O/Ne dominated material as appropriate,
106 : ! otherwise skip phase separation
107 0 : if(components == 'CO' .and. XO + XC < 0.9d0) return
108 0 : if(components == 'ONe' .and. XNe20 + XNe22 + XO < 0.8d0) return ! O/Ne mixtures tend to have more byproducts of burning mixed in
109 :
110 : ! If there is a phase transition, reset the composition at the boundary
111 0 : if(k_bound > 0) then
112 :
113 : ! core boundary needs to be padded by a minimal amount (less than a zone worth of mass)
114 : ! to account for loss of precision during remeshing.
115 0 : pad = s% min_dq * s% m(1) * 0.5d0
116 0 : do k = s%nz,1,-1
117 0 : if(s% m(k) > s% crystal_core_boundary_mass + pad) then
118 : kstart = k
119 : exit
120 : end if
121 : end do
122 :
123 : ! calculate energy associated with phase separation, ignoring the ionization
124 : ! energy term that Skye sometimes calculates
125 0 : save_Skye_use_ion_offsets = s% eos_rq% Skye_use_ion_offsets
126 0 : s% eos_rq% Skye_use_ion_offsets = .false.
127 0 : call update_model_(s,1,s%nz,.false.)
128 0 : do k=1,s% nz
129 0 : s% eps_phase_separation(k) = s% energy(k)
130 : end do
131 :
132 : ! loop runs outward starting at previous crystallization boundary
133 0 : do k = kstart,1,-1
134 : ! Start by checking if this material should be crystallizing
135 0 : if(s% phase(k) <= eos_phase_boundary) then
136 0 : s% crystal_core_boundary_mass = s% m(k+1)
137 0 : exit
138 : end if
139 :
140 0 : if(components == '3c') then
141 0 : call move_one_zone_MCT_3comp(s,k,components)
142 : else
143 0 : call move_one_zone_for_2comp(s,k,components)
144 : end if
145 :
146 : ! crystallized out to k now, liquid starts at k-1.
147 : ! now mix the liquid material outward until stably stratified
148 0 : call mix_outward(s, k-1, 0)
149 :
150 : end do
151 :
152 0 : call update_model_(s,1,s%nz,.false.)
153 :
154 : ! phase separation heating term for use by energy equation
155 0 : do k=1,s% nz
156 0 : s% eps_phase_separation(k) = (s% eps_phase_separation(k) - s% energy(k)) / dt
157 : end do
158 0 : s% eos_rq% Skye_use_ion_offsets = save_Skye_use_ion_offsets
159 0 : s% need_to_setvars = .true.
160 : end if
161 :
162 0 : ierr = 0
163 : end subroutine separate_mix_and_heat
164 :
165 0 : subroutine move_one_zone_for_2comp(s,k,components)
166 : use chem_def, only: ic12, io16, ine20, ine22
167 : use chem_lib, only: chem_get_iso_id
168 : type(star_info), pointer :: s
169 : integer, intent(in) :: k
170 : character (len=*), intent(in) :: components
171 : real(dp) :: XC, XO, XNe20, XNe22, XC1, XO1, XNe120, XNe122, dXO, Xfac
172 : real(dp), dimension(2) :: dXNe
173 : integer :: net_ic12, net_io16, net_ine20, net_ine22
174 :
175 0 : net_ic12 = s% net_iso(ic12)
176 0 : net_io16 = s% net_iso(io16)
177 0 : net_ine20 = s% net_iso(ine20)
178 0 : net_ine22 = s% net_iso(ine22)
179 :
180 0 : if(components == 'CO') then
181 0 : XO = s% xa(net_io16,k)
182 0 : XC = s% xa(net_ic12,k)
183 :
184 : ! Call Blouin phase diagram.
185 : ! Need to rescale temporarily because phase diagram assumes XO + XC = 1
186 0 : Xfac = XO + XC
187 0 : XO = XO/Xfac
188 0 : XC = XC/Xfac
189 :
190 0 : dXO = blouin_delta_xo(XO)
191 :
192 0 : s% xa(net_io16,k) = Xfac*(XO + dXO)
193 0 : s% xa(net_ic12,k) = Xfac*(XC - dXO)
194 :
195 : ! Redistribute change in C,O into zone k-1,
196 : ! conserving total mass of C,O
197 0 : XC1 = s% xa(net_ic12,k-1)
198 0 : XO1 = s% xa(net_io16,k-1)
199 0 : s% xa(net_ic12,k-1) = XC1 + Xfac*dXO * s% dq(k) / s% dq(k-1)
200 0 : s% xa(net_io16,k-1) = XO1 - Xfac*dXO * s% dq(k) / s% dq(k-1)
201 0 : else if(components == 'ONe') then
202 0 : XNe20 = s% xa(net_ine20,k)
203 0 : XNe22 = s% xa(net_ine22,k)
204 0 : XO = s% xa(net_io16,k)
205 :
206 : ! Call Blouin phase diagram.
207 : ! Need to rescale temporarily because phase diagram assumes XO + XNe = 1
208 0 : Xfac = XO + XNe20 + XNe22
209 0 : XO = XO/Xfac
210 0 : XNe20 = XNe20/Xfac
211 0 : XNe22 = XNe22/Xfac
212 :
213 0 : dXNe = blouin_delta_xne(XNe20,XNe22)
214 :
215 0 : s% xa(net_ine20,k) = Xfac*(XNe20 + dXNe(1))
216 0 : s% xa(net_ine22,k) = Xfac*(XNe22 + dXNe(2))
217 0 : s% xa(net_io16,k) = Xfac*(XO - dXNe(1) - dXNe(2))
218 :
219 : ! Redistribute change in Ne,O into zone k-1,
220 : ! conserving total mass of Ne,O
221 0 : XO1 = s% xa(net_io16,k-1)
222 0 : XNe120 = s% xa(net_ine20,k-1)
223 0 : XNe122 = s% xa(net_ine22,k-1)
224 0 : s% xa(net_io16,k-1) = XO1 + Xfac*(dXNe(1) + dXNe(2)) * s% dq(k) / s% dq(k-1)
225 0 : s% xa(net_ine20,k-1) = XNe120 - Xfac*dXNe(1) * s% dq(k) / s% dq(k-1)
226 0 : s% xa(net_ine22,k-1) = XNe122 - Xfac*dXNe(2) * s% dq(k) / s% dq(k-1)
227 : else
228 0 : write(*,*) 'invalid components option in 2-component phase separation'
229 0 : stop
230 : end if
231 :
232 0 : call update_model_(s,k-1,s%nz,.true.)
233 :
234 0 : end subroutine move_one_zone_for_2comp
235 :
236 0 : subroutine move_one_zone_MCT_3comp(s,k,components)
237 : use chem_def, only: ic12, io16, ine20, ine22, ina23, img24
238 : use chem_lib, only: chem_get_iso_id
239 : type(star_info), pointer :: s
240 : integer, intent(in) :: k
241 : real(dp), dimension(4) :: Dd
242 : character (len=*), intent(in) :: components
243 : real(dp) :: XC, XO, XNe20, XNe22, XNa, XMg, XC1, XO1, XNe120, XNe122, XNa1, XMg1, Xfac
244 : integer :: net_ic12, net_io16, net_ine20, net_ine22, net_ina23, net_img24
245 :
246 0 : net_ic12 = s% net_iso(ic12)
247 0 : net_io16 = s% net_iso(io16)
248 0 : net_ine20 = s% net_iso(ine20)
249 0 : net_ine22 = s% net_iso(ine22)
250 0 : net_ina23 = s% net_iso(ina23)
251 0 : net_img24 = s% net_iso(img24)
252 0 : XO = s% xa(net_io16,k)
253 0 : XC = s% xa(net_ic12,k)
254 0 : XNe20 = s% xa(net_ine20,k)
255 0 : XNe22 = s% xa(net_ine22,k)
256 0 : XNa = s% xa(net_ina23,k)
257 0 : XMg = s% xa(net_img24,k)
258 :
259 : ! check the abundances to decide which table use for interpolation
260 0 : if (XO + XC + XNe20 + XNe22 > 0.7d0 .and. XC > XMg .and. XC > XNa) then
261 0 : Xfac = XO + XC + XNe20 + XNe22
262 0 : XO = XO/Xfac
263 0 : XC = XC/Xfac
264 0 : XNe20 = XNe20/Xfac
265 0 : XNe22 = XNe22/Xfac
266 : ! call the deltas resulting from interpolation (in mass fraction)
267 0 : call medin_cumming_3p_d_cone(XC,XO,XNe20,XNe22,Dd)
268 : ! apply fractionation as given by the deltas from interpolation
269 0 : s% xa(net_ic12,k) = Xfac*(XC + Dd(1))
270 0 : s% xa(net_io16,k) = Xfac*(XO + Dd(2))
271 0 : s% xa(net_ine20,k) = Xfac*(XNe20 + Dd(3))
272 0 : s% xa(net_ine22,k) = Xfac*(XNe22 + Dd(4))
273 0 : XC1 = s% xa(net_ic12,k-1)
274 0 : XO1 = s% xa(net_io16,k-1)
275 0 : XNe120 = s% xa(net_ine20,k-1)
276 0 : XNe122 = s% xa(net_ine22,k-1)
277 0 : s% xa(net_ic12,k-1) = XC1 - Xfac*Dd(1) * s% dq(k) / s% dq(k-1)
278 0 : s% xa(net_io16,k-1) = XO1 - Xfac*Dd(2) * s% dq(k) / s% dq(k-1)
279 0 : s% xa(net_ine20,k-1) = XNe120 - Xfac*(Dd(3)) * s% dq(k) / s% dq(k-1)
280 0 : s% xa(net_ine22,k-1) = XNe122 - Xfac*(Dd(4)) * s% dq(k) / s% dq(k-1)
281 : ! write(*,*) 'phase 3 CONe abundances',XC,XO,XNe20+XNe22
282 0 : else if (XO + XNe20 + XNe22 + XMg > 0.7d0 .and. XMg > XC .and. XMg > XNa) then
283 0 : Xfac = XO + XNe20 + XNe22 + XMg
284 0 : XMg = XMg/Xfac
285 0 : XO = XO/Xfac
286 0 : XNe20 = XNe20/Xfac
287 0 : XNe22 = XNe22/Xfac
288 : ! call the deltas resulting from interpolation (in mass fraction)
289 0 : call medin_cumming_3p_d_neomg(XMg,XO,XNe20,XNe22,Dd)
290 : ! apply fractionation as given by the deltas from interpolation
291 0 : s% xa(net_img24,k) = Xfac*(XMg + Dd(1))
292 0 : s% xa(net_io16,k) = Xfac*(XO + Dd(2))
293 0 : s% xa(net_ine20,k) = Xfac*(XNe20 + Dd(3))
294 0 : s% xa(net_ine22,k) = Xfac*(XNe22 + Dd(4))
295 0 : XMg1 = s% xa(net_img24,k-1)
296 0 : XO1 = s% xa(net_io16,k-1)
297 0 : XNe120 = s% xa(net_ine20,k-1)
298 0 : XNe122 = s% xa(net_ine22,k-1)
299 0 : s% xa(net_img24,k-1) = XMg1 - Xfac*Dd(1) * s% dq(k) / s% dq(k-1)
300 0 : s% xa(net_io16,k-1) = XO1 - Xfac*Dd(2) * s% dq(k) / s% dq(k-1)
301 0 : s% xa(net_ine20,k-1) = XNe120 - Xfac*Dd(3) * s% dq(k) / s% dq(k-1)
302 0 : s% xa(net_ine22,k-1) = XNe122 - Xfac*Dd(4) * s% dq(k) / s% dq(k-1)
303 : ! write(*,*) 'phase 3 ONeMg abundances',XO,XNe20+XNe22,XMg
304 0 : else if (XO + XNe20 + XNe22 + XNa > 0.7d0 .and. XNa > XC .and. XNa > XMg) then
305 0 : Xfac = XO + XNe20 + XNe22 + XNa
306 0 : XNa = XNa/Xfac
307 0 : XO = XO/Xfac
308 0 : XNe20 = XNe20/Xfac
309 0 : XNe22 = XNe22/Xfac
310 : ! call the deltas resulting from interpolation (in mass fraction)
311 0 : call medin_cumming_3p_d_onena(XNa,XO,XNe20,XNe22,Dd)
312 : ! apply fractionation as given by the deltas from interpolation
313 0 : s% xa(net_ina23,k) = Xfac*(XNa + Dd(1))
314 0 : s% xa(net_io16,k) = Xfac*(XO + Dd(2))
315 0 : s% xa(net_ine20,k) = Xfac*(XNe20 + Dd(3))
316 0 : s% xa(net_ine22,k) = Xfac*(XNe22 + Dd(4))
317 0 : XNa1 = s% xa(net_ina23,k-1)
318 0 : XO1 = s% xa(net_io16,k-1)
319 0 : XNe120 = s% xa(net_ine20,k-1)
320 0 : XNe122 = s% xa(net_ine22,k-1)
321 0 : s% xa(net_ina23,k-1) = XNa1 - Xfac*Dd(1) * s% dq(k) / s% dq(k-1)
322 0 : s% xa(net_io16,k-1) = XO1 - Xfac*Dd(2) * s% dq(k) / s% dq(k-1)
323 0 : s% xa(net_ine20,k-1) = XNe120 - Xfac*Dd(3) * s% dq(k) / s% dq(k-1)
324 0 : s% xa(net_ine22,k-1) = XNe122 - Xfac*Dd(4) * s% dq(k) / s% dq(k-1)
325 : ! write(*,*) 'phase 3 ONeNa abundances',XO,XNe20+XNe22,XNa
326 0 : else if (XC + XO + XMg > 0.7d0 .and. XMg > XNa .and. XMg > XNe20+XNe22) then
327 0 : Xfac = XC + XO + XMg
328 0 : XC = XC/Xfac
329 0 : XO = XO/Xfac
330 0 : XMg = XMg/Xfac
331 : ! call the deltas resulting from interpolation (in mass fraction)
332 0 : call medin_cumming_3p_d_comg(XC,XMg,XO,Dd)
333 : ! apply fractionation as given by the deltas from interpolation
334 0 : s% xa(net_ic12,k) = Xfac*(XC + Dd(1))
335 0 : s% xa(net_img24,k) = Xfac*(XMg + Dd(2))
336 0 : s% xa(net_io16,k) = Xfac*(XO - (Dd(1) + Dd(2)))
337 0 : XC1 = s% xa(net_ic12,k-1)
338 0 : XO1 = s% xa(net_io16,k-1)
339 0 : XMg1 = s% xa(net_img24,k-1)
340 0 : s% xa(net_ic12,k-1) = XC1 - Xfac*Dd(1) * s% dq(k) / s% dq(k-1)
341 0 : s% xa(net_img24,k-1) = XMg1 - Xfac*Dd(2) * s% dq(k) / s% dq(k-1)
342 0 : s% xa(net_io16,k-1) = XO1 + Xfac*(Dd(1)+Dd(2)) * s% dq(k) / s% dq(k-1)
343 : ! write(*,*) 'phase 3 COMg abundances',XC,XO,XMg
344 : end if
345 :
346 0 : call update_model_(s,k-1,s%nz,.true.)
347 :
348 0 : end subroutine move_one_zone_MCT_3comp
349 :
350 : ! mix composition outward until reaching stable composition profile
351 0 : subroutine mix_outward(s,kbot,min_mix_zones)
352 : type(star_info), pointer :: s
353 : integer, intent(in) :: kbot, min_mix_zones
354 :
355 0 : real(dp) :: avg_xa(s%species)
356 : real(dp) :: mass, B_term, grada, gradr
357 : integer :: k, l, ktop
358 : logical :: use_brunt
359 :
360 0 : use_brunt = s% phase_separation_mixing_use_brunt
361 :
362 0 : do k=kbot-min_mix_zones,1,-1
363 0 : ktop = k
364 :
365 0 : if (s% m(ktop) > s% phase_sep_mixing_mass) then
366 0 : s% phase_sep_mixing_mass = s% m(ktop)
367 : end if
368 :
369 0 : mass = SUM(s%dm(ktop:kbot))
370 0 : do l = 1, s%species
371 0 : avg_xa(l) = SUM(s%dm(ktop:kbot)*s%xa(l,ktop:kbot))/mass
372 : end do
373 :
374 : ! some potential safeguards from conv_premix
375 : ! avg_xa = MAX(MIN(avg_xa, 1._dp), 0._dp)
376 : ! avg_xa = avg_xa/SUM(avg_xa)
377 :
378 0 : do l = 1, s%species
379 0 : s%xa(l,ktop:kbot) = avg_xa(l)
380 : end do
381 :
382 : ! updates, eos, opacities, mu, etc now that abundances have changed,
383 : ! but only in the cells near the boundary where we need to check here.
384 : ! Will call full update over mixed region after exiting loop.
385 0 : call update_model_(s, ktop-1, ktop+1, use_brunt)
386 :
387 0 : if(use_brunt) then
388 0 : B_term = s% unsmoothed_brunt_B(ktop)
389 0 : grada = s% grada_face(ktop)
390 0 : gradr = s% gradr(ktop)
391 0 : if(B_term + grada - gradr > 0d0) then
392 : ! stable against further mixing, so exit loop
393 : exit
394 : end if
395 : else ! simpler calculation based on mu gradient
396 0 : if(s% mu(ktop) >= s% mu(ktop-1)) then
397 : ! stable against further mixing, so exit loop
398 : exit
399 : end if
400 : end if
401 :
402 : end do
403 :
404 : ! Call a final update over all mixed cells now.
405 0 : call update_model_(s, ktop, kbot+1, .true.)
406 :
407 0 : end subroutine mix_outward
408 :
409 0 : real(dp) function blouin_delta_xo(Xin)
410 : real(dp), intent(in) :: Xin ! mass fraction
411 : real(dp) :: Xnew ! mass fraction
412 : real(dp) :: xo, dxo ! number fractions
413 : real(dp) :: a0, a1, a2, a3, a4, a5
414 :
415 : ! Convert input mass fraction to number fraction, assuming C/O mixture
416 0 : xo = (Xin/16d0)/(Xin/16d0 + (1d0 - Xin)/12d0)
417 :
418 0 : a0 = 0d0
419 0 : a1 = -0.311540d0
420 0 : a2 = 2.114743d0
421 0 : a3 = -1.661095d0
422 0 : a4 = -1.406005d0
423 0 : a5 = 1.263897d0
424 :
425 : dxo = &
426 : a0 + &
427 : a1*xo + &
428 : a2*xo*xo + &
429 : a3*xo*xo*xo + &
430 : a4*xo*xo*xo*xo + &
431 0 : a5*xo*xo*xo*xo*xo
432 :
433 0 : xo = xo + dxo
434 :
435 : ! Convert back to mass fraction
436 0 : Xnew = 16d0*xo/(16d0*xo + 12d0*(1d0-xo))
437 :
438 0 : blouin_delta_xo = Xnew - Xin
439 0 : end function blouin_delta_xo
440 :
441 0 : function blouin_delta_xne(Xin20,Xin22)
442 : real(dp), intent(in) :: Xin20, Xin22 ! mass fraction
443 : real(dp) :: Xnew20, Xnew22 ! mass fraction
444 : real(dp) :: xne, dxne, xne1, xne2 ! number fractions
445 : real(dp) :: a0, a1, a2, a3, a4, a5
446 : real(dp), dimension(2) :: blouin_delta_xne
447 :
448 : ! Convert input mass fraction to number fraction, assuming O/Ne mixture
449 0 : xne1 =(Xin20/20d0)/(Xin20/20d0 + Xin22/22d0 + (1d0 - Xin20 - Xin22)/16d0)
450 0 : xne2 =(Xin22/22d0)/(Xin20/20d0 + Xin22/22d0 + (1d0 - Xin20 - Xin22)/16d0)
451 : ! isotope 22Ne is added to the Ne separation along with 20Ne
452 0 : xne = xne1 + xne2
453 :
454 0 : a0 = 0d0
455 0 : a1 = -0.120299d0
456 0 : a2 = 1.304399d0
457 0 : a3 = -1.722625d0
458 0 : a4 = 0.393996d0
459 0 : a5 = 0.144529d0
460 :
461 : dxne = &
462 : a0 + &
463 : a1*xne + &
464 : a2*xne*xne + &
465 : a3*xne*xne*xne + &
466 : a4*xne*xne*xne*xne + &
467 0 : a5*xne*xne*xne*xne*xne
468 :
469 0 : xne1 = xne1 + dxne*xne1/xne
470 0 : xne2 = xne2 + dxne*xne2/xne
471 0 : xne = xne1 + xne2
472 :
473 : ! Convert back to mass fraction
474 0 : Xnew20 = (20d0*xne1)/(20d0*xne1 + 22d0*xne2 + 16d0*(1d0-xne))
475 0 : Xnew22 = (22d0*xne2)/(20d0*xne1 + 22d0*xne2 + 16d0*(1d0-xne))
476 :
477 0 : blouin_delta_xne(1) = Xnew20 - Xin20
478 0 : blouin_delta_xne(2) = Xnew22 - Xin22
479 : end function blouin_delta_xne
480 :
481 0 : subroutine tab_interp_medin_cumming_dx1(x1_,x2_,components,dx1_)
482 : use interp_2D_lib_db, only: interp_mkbicub_db, interp_evbicub_db
483 : use utils_lib, only: mesa_error, mkdir, is_bad
484 : implicit none
485 : integer :: ilinx,iliny,ibcxmin,ibcxmax,ibcymin,ibcymax,ict(6),ierr,i,j,k
486 : integer :: num_x1, num_x2
487 0 : real(dp), allocatable :: bcxmin(:), bcxmax(:)
488 0 : real(dp), allocatable :: bcymin(:), bcymax(:)
489 0 : real(dp), pointer, dimension(:) :: x1_l, x2_l, deltax1_sob_f1
490 0 : real(dp), pointer :: deltax1_sob_f(:,:,:)
491 0 : real(dp), allocatable :: delta_grid(:,:)
492 : real(dp), intent(in) :: x1_,x2_ ! target of this interpolation
493 : character (len=*), intent(in) :: components
494 : character(len=256) :: filename, phase_sep_data_dir
495 : real(dp) :: fval(6) ! output data
496 : real(dp), intent(out) :: dx1_
497 : integer :: ier
498 0 : type(hdf5io_t) :: hi
499 :
500 0 : ict = 0
501 0 : ict(1) = 1
502 :
503 0 : phase_sep_data_dir = trim(mesa_data_dir) // '/star_data/phase_separation'
504 :
505 : ! setup interpolation table for x1 x2 dx1
506 0 : if (components=='CONe') then
507 0 : filename = trim(phase_sep_data_dir) // '/CONe_deltaC.h5'
508 0 : else if (components=='NeOMg') then
509 0 : filename = trim(phase_sep_data_dir) // '/NeOMg_deltaMg.h5'
510 0 : else if (components=='ONeNa') then
511 0 : filename = trim(phase_sep_data_dir) // '/ONeNa_deltaNa.h5'
512 0 : else if (components=='COMg') then
513 0 : filename = trim(phase_sep_data_dir) // '/COMg_deltaC.h5'
514 : end if
515 :
516 : ! Open HDF5 file
517 0 : hi = hdf5io_t(filename, OPEN_FILE_RO)
518 :
519 : ! Read grid dimensions
520 0 : call hi%read_attr('num_x1', num_x1)
521 0 : call hi%read_attr('num_x2', num_x2)
522 :
523 : ! Allocate arrays
524 : allocate(x1_l(num_x1), x2_l(num_x2), &
525 0 : deltax1_sob_f1(4*num_x1*num_x2))
526 0 : allocate(delta_grid(num_x1, num_x2))
527 0 : allocate(bcxmin(num_x1), bcxmax(num_x1))
528 0 : allocate(bcymin(num_x2), bcymax(num_x2))
529 :
530 : deltax1_sob_f(1:4,1:num_x1,1:num_x2) => &
531 0 : deltax1_sob_f1(1:4*num_x1*num_x2)
532 :
533 : ! Read data from HDF5
534 0 : call hi%read_dset('x1', x1_l)
535 0 : call hi%read_dset('x2', x2_l)
536 0 : call hi%read_dset('delta', delta_grid)
537 :
538 : ! Copy delta data to interpolation array
539 0 : do j=1,num_x1
540 0 : do i=1,num_x2
541 0 : deltax1_sob_f(1,j,i) = delta_grid(j,i)
542 : end do
543 : end do
544 :
545 : ! Close HDF5 file
546 0 : call hi%final()
547 :
548 0 : deallocate(delta_grid)
549 :
550 : ! just use "not a knot" bc's at edges of tables
551 0 : ibcxmin = 0; bcxmin(1:num_x1) = 0
552 0 : ibcxmax = 0; bcxmax(1:num_x1) = 0
553 0 : ibcymin = 0; bcymin(1:num_x2) = 0
554 0 : ibcymax = 0; bcymax(1:num_x2) = 0
555 : call interp_mkbicub_db( &
556 : x1_l, num_x1, x2_l, num_x2, deltax1_sob_f1, num_x1, &
557 : ibcxmin,bcxmin,ibcxmax,bcxmax, &
558 : ibcymin,bcymin,ibcymax,bcymax, &
559 0 : ilinx,iliny,ierr)
560 0 : if (ierr /= 0) then
561 0 : write(*,*) 'interp_mkbicub_db error'
562 0 : ierr = -1
563 0 : call mesa_error(__FILE__,__LINE__)
564 : end if
565 0 : do j=1,num_x1
566 0 : do i=1,num_x2
567 0 : do k=1,4
568 0 : if (is_bad(deltax1_sob_f(k,j,i))) then
569 0 : write(*,*) 'deltax1_sob_f', i, j, k, deltax1_sob_f(k,j,i)
570 : end if
571 : end do
572 : end do
573 : end do
574 : call interp_evbicub_db( &
575 : x1_, x2_, x1_l, num_x1, x2_l, num_x2, &
576 0 : ilinx, iliny, deltax1_sob_f1, num_x1, ict, fval, ier)
577 0 : dx1_=fval(1) ! delta_x1 from 2d interpolation
578 :
579 0 : deallocate(x1_l, x2_l, deltax1_sob_f1)
580 0 : deallocate(bcxmin, bcxmax, bcymin, bcymax)
581 0 : end subroutine tab_interp_medin_cumming_dx1
582 :
583 :
584 0 : subroutine tab_interp_medin_cumming_dx2(x1_,x2_,components,dx2_)
585 : use interp_2D_lib_db, only: interp_mkbicub_db, interp_evbicub_db
586 : use utils_lib, only: mesa_error, mkdir, is_bad
587 : implicit none
588 : integer :: ilinx,iliny,ibcxmin,ibcxmax,ibcymin,ibcymax,ict(6),ierr,i,j,k
589 : integer :: num_x1, num_x2
590 0 : real(dp), allocatable :: bcxmin(:), bcxmax(:)
591 0 : real(dp), allocatable :: bcymin(:), bcymax(:)
592 0 : real(dp), pointer, dimension(:) :: x1_l, x2_l, deltax1_sob_f1
593 0 : real(dp), pointer :: deltax1_sob_f(:,:,:)
594 0 : real(dp), allocatable :: delta_grid(:,:)
595 : real(dp), intent(in) :: x1_,x2_ ! target of this interpolation
596 : character (len=*), intent(in) :: components
597 : character(len=256) :: filename, phase_sep_data_dir
598 : real(dp) :: fval(6) ! output data
599 : real(dp), intent(out) :: dx2_
600 : integer :: ier
601 0 : type(hdf5io_t) :: hi
602 :
603 0 : ict = 0
604 0 : ict(1) = 1
605 :
606 0 : phase_sep_data_dir = trim(mesa_data_dir) // '/star_data/phase_separation'
607 :
608 : ! setup interpolation table for tau sob eta
609 0 : if (components=='CONe') then
610 0 : filename = trim(phase_sep_data_dir) // '/CONe_deltaO.h5'
611 0 : else if (components=='NeOMg') then
612 0 : filename = trim(phase_sep_data_dir) // '/NeOMg_deltaO.h5'
613 0 : else if (components=='ONeNa') then
614 0 : filename = trim(phase_sep_data_dir) // '/ONeNa_deltaO.h5'
615 0 : else if (components=='COMg') then
616 0 : filename = trim(phase_sep_data_dir) // '/COMg_deltaMg.h5'
617 : end if
618 :
619 : ! Open HDF5 file
620 0 : hi = hdf5io_t(filename, OPEN_FILE_RO)
621 :
622 : ! Read grid dimensions
623 0 : call hi%read_attr('num_x1', num_x1)
624 0 : call hi%read_attr('num_x2', num_x2)
625 :
626 : ! Allocate arrays
627 : allocate(x1_l(num_x1), x2_l(num_x2), &
628 0 : deltax1_sob_f1(4*num_x1*num_x2))
629 0 : allocate(delta_grid(num_x1, num_x2))
630 0 : allocate(bcxmin(num_x1), bcxmax(num_x1))
631 0 : allocate(bcymin(num_x2), bcymax(num_x2))
632 :
633 : deltax1_sob_f(1:4,1:num_x1,1:num_x2) => &
634 0 : deltax1_sob_f1(1:4*num_x1*num_x2)
635 :
636 : ! Read data from HDF5
637 0 : call hi%read_dset('x1', x1_l)
638 0 : call hi%read_dset('x2', x2_l)
639 0 : call hi%read_dset('delta', delta_grid)
640 :
641 : ! Copy delta data to interpolation array
642 0 : do j=1,num_x1
643 0 : do i=1,num_x2
644 0 : deltax1_sob_f(1,j,i) = delta_grid(j,i)
645 : end do
646 : end do
647 :
648 : ! Close HDF5 file
649 0 : call hi%final()
650 :
651 0 : deallocate(delta_grid)
652 :
653 : ! just use "not a knot" bc's at edges of tables
654 0 : ibcxmin = 0; bcxmin(1:num_x1) = 0
655 0 : ibcxmax = 0; bcxmax(1:num_x1) = 0
656 0 : ibcymin = 0; bcymin(1:num_x2) = 0
657 0 : ibcymax = 0; bcymax(1:num_x2) = 0
658 : call interp_mkbicub_db( &
659 : x1_l, num_x1, x2_l, num_x2, deltax1_sob_f1, num_x1, &
660 : ibcxmin,bcxmin,ibcxmax,bcxmax, &
661 : ibcymin,bcymin,ibcymax,bcymax, &
662 0 : ilinx,iliny,ierr)
663 0 : if (ierr /= 0) then
664 0 : write(*,*) 'interp_mkbicub_db error'
665 0 : ierr = -1
666 0 : call mesa_error(__FILE__,__LINE__)
667 : end if
668 0 : do j=1,num_x1
669 0 : do i=1,num_x2
670 0 : do k=1,4
671 0 : if (is_bad(deltax1_sob_f(k,j,i))) then
672 0 : write(*,*) 'deltax1_sob_f', i, j, k, deltax1_sob_f(k,j,i)
673 : end if
674 : end do
675 : end do
676 : end do
677 : call interp_evbicub_db( &
678 : x1_, x2_, x1_l, num_x1, x2_l, num_x2, &
679 0 : ilinx, iliny, deltax1_sob_f1, num_x1, ict, fval, ier)
680 0 : dx2_=fval(1) ! delta_x2 from 2d interpolation
681 :
682 0 : deallocate(x1_l, x2_l, deltax1_sob_f1)
683 0 : deallocate(bcxmin, bcxmax, bcymin, bcymax)
684 0 : end subroutine tab_interp_medin_cumming_dx2
685 :
686 :
687 0 : subroutine medin_cumming_3p_d_cone(X1,X2,X3_1,X3_2,Dd)
688 : real(dp), intent(in) :: X1, X2, X3_1, X3_2 ! mass fraction
689 : real(dp), dimension(4),intent(out) :: Dd
690 : real(dp) :: Xnew1, Xnew2, Xnew3_1, Xnew3_2, Xfac ! mass fraction
691 : real(dp) :: xc, dxc, xo, dxo, xne1, xne2 ! number fractions
692 : real(dp) :: dx1_,dx2_
693 :
694 0 : Xfac = X1 + X2 + X3_1 + X3_2
695 0 : xc = (X1/12)/(X1/12 + X2/16 + X3_1/20 + X3_2/22)
696 0 : xo = (X2/16)/(X1/12 + X2/16 + X3_1/20 + X3_2/22)
697 0 : xne1 = (X3_1/20)/(X1/12 + X2/16 + X3_1/20 + X3_2/22)
698 0 : xne2 = (X3_2/22)/(X1/12 + X2/16 + X3_1/20 + X3_2/22)
699 0 : call tab_interp_medin_cumming_dx1(xc,xo,'CONe',dx1_)
700 0 : call tab_interp_medin_cumming_dx2(xc,xo,'CONe',dx2_)
701 0 : dxc=dx1_
702 0 : dxo=dx2_
703 : !write(*,*) 'delta_xc: ',dxc,' delta_xo: ', dxo
704 0 : xc = xc + dxc
705 0 : xo = xo + dxo
706 : ! convert deltas in number fraction to mass fraction
707 0 : Xnew1 = 12*xc/(12*xc + 16*xo + 20*(1-xc-xo)*(xne1)/(xne1+xne2)+22*(1-xc-xo)*(xne2)/(xne1+xne2))
708 0 : Xnew2 = 16*xo/(12*xc + 16*xo + 20*(1-xc-xo)*(xne1)/(xne1+xne2)+22*(1-xc-xo)*(xne2)/(xne1+xne2))
709 : Xnew3_1 = (20*(1-xc-xo)*(xne1)/(xne1+xne2))/(12*xc + 16*xo + 20*(1-xc-xo)*(xne1) &
710 0 : /(xne1+xne2)+22*(1-xc-xo)*(xne2)/(xne1+xne2))
711 : Xnew3_2 = (22*(1-xc-xo)*(xne2)/(xne1+xne2))/(12*xc + 16*xo + 20*(1-xc-xo)*(xne1) &
712 0 : /(xne1+xne2)+22*(1-xc-xo)*(xne2)/(xne1+xne2))
713 : Dd=[0,0,0,0]
714 0 : Dd(1)= Xnew1 - X1
715 0 : Dd(2)= Xnew2 - X2
716 0 : Dd(3)= Xnew3_1 - X3_1
717 0 : Dd(4)= Xnew3_2 - X3_2
718 : !write(*,*) 'delta_XC: ',Dd(1),' delta_XO: ', Dd(2), 'delta_XNe:', Dd(3)+Dd(4)
719 0 : end subroutine medin_cumming_3p_d_cone
720 :
721 0 : subroutine medin_cumming_3p_d_neomg(X1,X2,X3_1,X3_2,Dd)
722 : real(dp), intent(in) :: X1, X2, X3_1, X3_2 ! mass fraction
723 : real(dp), dimension(4),intent(out) :: Dd
724 : real(dp) :: Xnew1, Xnew2, Xnew3_1, Xnew3_2, Xfac ! mass fraction
725 : real(dp) :: xmg, dxmg, xo, dxo, xne1, xne2 ! number fractions
726 : real(dp) :: dx1_,dx2_
727 :
728 0 : Xfac = X1 + X2 + X3_1 + X3_2
729 0 : xmg = (X1/24)/(X1/24 + X2/16 + X3_1/20 + X3_2/22)
730 0 : xo = (X2/16)/(X1/24 + X2/16 + X3_1/20 + X3_2/22)
731 0 : xne1 = (X3_1/20)/(X1/24 + X2/16 + X3_1/20 + X3_2/22)
732 0 : xne2 = (X3_2/22)/(X1/24 + X2/16 + X3_1/20 + X3_2/22)
733 0 : call tab_interp_medin_cumming_dx1(xmg,xo,'NeOMg',dx1_)
734 0 : call tab_interp_medin_cumming_dx2(xmg,xo,'NeOMg',dx2_)
735 0 : dxmg=dx1_
736 0 : dxo=dx2_
737 0 : xmg = xmg + dxmg
738 0 : xo = xo + dxo
739 : ! convert deltas in number fraction to mass fraction
740 0 : Xnew1 = 24*xmg/(24*xmg + 16*xo + 20*(1-xmg-xo)*(xne1)/(xne1+xne2)+22*(1-xmg-xo)*(xne2)/(xne1+xne2))
741 0 : Xnew2 = 16*xo/(24*xmg + 16*xo + 20*(1-xmg-xo)*(xne1)/(xne1+xne2)+22*(1-xmg-xo)*(xne2)/(xne1+xne2))
742 : Xnew3_1 = (20*(1-xmg-xo)*(xne1)/(xne1+xne2))/(24*xmg + 16*xo + 20*(1-xmg-xo)*(xne1) &
743 0 : /(xne1+xne2)+22*(1-xmg-xo)*(xne2)/(xne1+xne2))
744 : Xnew3_2 = (22*(1-xmg-xo)*(xne2)/(xne1+xne2))/(24*xmg + 16*xo + 20*(1-xmg-xo)*(xne1) &
745 0 : /(xne1+xne2)+22*(1-xmg-xo)*(xne2)/(xne1+xne2))
746 : Dd=[0,0,0,0]
747 0 : Dd(1)= Xnew1 - X1
748 0 : Dd(2)= Xnew2 - X2
749 0 : Dd(3)= Xnew3_1 - X3_1
750 0 : Dd(4)= Xnew3_2 - X3_2
751 : !write(*,*) 'delta_XMg: ',Dd(1),' delta_XO: ', Dd(2), 'delta_XNe:', Dd(3)+Dd(4)
752 0 : end subroutine medin_cumming_3p_d_neomg
753 :
754 0 : subroutine medin_cumming_3p_d_onena(X1,X2,X3_1,X3_2,Dd)
755 : real(dp), intent(in) :: X1, X2, X3_1, X3_2 ! mass fraction
756 : real(dp), dimension(4),intent(out) :: Dd
757 : real(dp) :: Xnew1, Xnew2, Xnew3_1, Xnew3_2, Xfac ! mass fraction
758 : real(dp) :: xna, dxna, xo, dxo, xne1, xne2 ! number fractions
759 : real(dp) :: dx1_,dx2_
760 :
761 0 : Xfac = X1 + X2 + X3_1 + X3_2
762 0 : xna = (X1/23)/(X1/23 + X2/16 + X3_1/20 + X3_2/22)
763 0 : xo = (X2/16)/(X1/23 + X2/16 + X3_1/20 + X3_2/22)
764 0 : xne1 = (X3_1/20)/(X1/23 + X2/16 + X3_1/20 + X3_2/22)
765 0 : xne2 = (X3_2/22)/(X1/23 + X2/16 + X3_1/20 + X3_2/22)
766 0 : call tab_interp_medin_cumming_dx1(xna,xo,'ONeNa',dx1_)
767 0 : call tab_interp_medin_cumming_dx2(xna,xo,'ONeNa',dx2_)
768 0 : dxna=dx1_
769 0 : dxo=dx2_
770 : !write(*,*) xna,xo
771 : !write(*,*) 'delta_xna: ',dxna,' delta_xo: ', dxo
772 0 : xna = xna + dxna
773 0 : xo = xo + dxo
774 : ! convert deltas in number fraction to mass fraction
775 0 : Xnew1 = 23*xna/(23*xna + 16*xo + 20*(1-xna-xo)*(xne1)/(xne1+xne2)+22*(1-xna-xo)*(xne2)/(xne1+xne2))
776 0 : Xnew2 = 16*xo/(23*xna + 16*xo + 20*(1-xna-xo)*(xne1)/(xne1+xne2)+22*(1-xna-xo)*(xne2)/(xne1+xne2))
777 : Xnew3_1 = (20*(1-xna-xo)*(xne1)/(xne1+xne2))/(23*xna + 16*xo + 20*(1-xna-xo)*(xne1) &
778 0 : /(xne1+xne2)+22*(1-xna-xo)*(xne2)/(xne1+xne2))
779 : Xnew3_2 = (22*(1-xna-xo)*(xne2)/(xne1+xne2))/(23*xna + 16*xo + 20*(1-xna-xo)*(xne1) &
780 0 : /(xne1+xne2)+22*(1-xna-xo)*(xne2)/(xne1+xne2))
781 : Dd=[0,0,0,0]
782 0 : Dd(1)= Xnew1 - X1
783 0 : Dd(2)= Xnew2 - X2
784 0 : Dd(3)= Xnew3_1 - X3_1
785 0 : Dd(4)= Xnew3_2 - X3_2
786 : !write(*,*) 'delta_XNa: ',Dd(1),' delta_XO: ', Dd(2), 'delta_XNe:', Dd(3)+Dd(4)
787 0 : end subroutine medin_cumming_3p_d_onena
788 :
789 0 : subroutine medin_cumming_3p_d_comg(X1,X2,X3,Dd)
790 : real(dp), intent(in) :: X1, X2, X3 ! mass fraction
791 : real(dp), dimension(4),intent(out) :: Dd
792 : real(dp) :: Xnew1, Xnew2, Xfac ! mass fraction
793 : real(dp) :: xc, dxc, xmg, dxmg, xo ! number fractions
794 : real(dp) :: dx1_,dx2_
795 :
796 0 : Xfac = X1 + X2 + X3
797 0 : xc = (X1/12)/(X1/12 + X2/24 + X3/16)
798 0 : xmg = (X2/24)/(X1/12 + X2/24 + X3/16)
799 0 : xo = (X3/16)/(X1/12 + X2/24 + X3/16)
800 0 : call tab_interp_medin_cumming_dx1(xc,xmg,'COMg',dx1_)
801 0 : call tab_interp_medin_cumming_dx2(xc,xmg,'COMg',dx2_)
802 0 : dxc=dx1_
803 0 : dxmg=dx2_
804 0 : xc = xc + dxc
805 0 : xmg = xmg + dxmg
806 : ! convert deltas in number fraction to mass fraction
807 0 : Xnew1 = 12*xc/(12*xc + 24*xmg + 16*(1-xc-xmg))
808 0 : Xnew2 = 24*xmg/(12*xc + 24*xmg + 16*(1-xc-xmg))
809 0 : Dd=[0,0,0,0]
810 0 : Dd(1)= Xnew1 - X1
811 0 : Dd(2)= Xnew2 - X2
812 0 : end subroutine medin_cumming_3p_d_comg
813 :
814 0 : subroutine update_model_ (s, kc_t, kc_b, do_brunt)
815 :
816 : use turb_info, only: set_mlt_vars
817 : use brunt, only: do_brunt_B
818 : use micro
819 :
820 : type(star_info), pointer :: s
821 : integer, intent(in) :: kc_t
822 : integer, intent(in) :: kc_b
823 : logical, intent(in) :: do_brunt
824 :
825 : integer :: ierr
826 : integer :: kf_t
827 : integer :: kf_b
828 :
829 0 : logical :: mask(s%nz)
830 :
831 0 : mask(:) = .true.
832 :
833 : ! Update the model to reflect changes in the abundances across
834 : ! cells kc_t:kc_b (the mask part of this call is unused, mask=true for all zones).
835 : ! Do updates at constant (P,T) rather than constant (rho,T).
836 0 : s%fix_Pgas = .true.
837 0 : call set_eos_with_mask(s, kc_t, kc_b, mask, ierr)
838 0 : if (ierr /= 0) then
839 0 : write(*,*) 'phase_separation: error from call to set_eos_with_mask'
840 0 : stop
841 : end if
842 0 : s%fix_Pgas = .false.
843 :
844 : ! Update opacities across cells kc_t:kc_b (this also sets rho_face
845 : ! and related quantities on faces kc_t:kc_b)
846 : call set_micro_vars(s, kc_t, kc_b, &
847 0 : skip_eos=.TRUE., skip_net=.TRUE., skip_neu=.TRUE., skip_kap=.FALSE., ierr=ierr)
848 0 : if (ierr /= 0) then
849 0 : write(*,*) 'phase_separation: error from call to set_micro_vars'
850 0 : stop
851 : end if
852 :
853 : ! This is expensive, so only do it if we really need to.
854 0 : if(do_brunt) then
855 : ! Need to make sure we can set brunt for mix_outward calculation.
856 0 : if(.not. s% calculate_Brunt_B) then
857 0 : stop "phase separation requires s% calculate_Brunt_B = .true."
858 : end if
859 0 : call do_brunt_B(s, kc_t, kc_b, ierr) ! for unsmoothed_brunt_B
860 0 : if (ierr /= 0) then
861 0 : write(*,*) 'phase_separation: error from call to do_brunt_B'
862 0 : stop
863 : end if
864 : end if
865 :
866 : ! Finally update MLT for interior faces
867 :
868 0 : kf_t = kc_t
869 0 : kf_b = kc_b + 1
870 :
871 0 : if (s% use_face_reconstruction) then
872 : ! update_model_ changed the local composition and refreshed EOS and
873 : ! kap on kc_t:kc_b, so the cached face thermo bundle must be rebuilt
874 : ! for the interior faces before set_mlt_vars uses it.
875 0 : s% reconstructed_face_state_valid(kf_t+1:kf_b-1) = .false.
876 : end if
877 :
878 0 : call set_mlt_vars(s, kf_t+1, kf_b-1, ierr)
879 0 : if (ierr /= 0) then
880 0 : write(*,*) 'phase_separation: failed in call to set_mlt_vars during update_model_'
881 0 : stop
882 : end if
883 :
884 0 : return
885 :
886 : end subroutine update_model_
887 :
888 0 : subroutine smooth_eps_phase_sep(s,dt,ierr)
889 : type (star_info), pointer :: s
890 : real(dp), intent(in) :: dt
891 : integer, intent(out) :: ierr
892 :
893 : real(dp) :: integrated_luminosity
894 : integer :: k, kmid
895 :
896 0 : integrated_luminosity = dot_product(s% dm(1:s%nz), s% eps_phase_separation(1:s%nz))
897 :
898 : ! redistribute evenly through the inner half of the star
899 0 : kmid = s%nz / 2
900 0 : do k = 1,s%nz
901 0 : if(s% q(k) < 0.5d0) then
902 : kmid = k
903 : exit
904 : end if
905 : end do
906 :
907 0 : s% eps_phase_separation(:) = 0d0
908 0 : s% eps_phase_separation(kmid:s%nz) = integrated_luminosity/s% m(kmid)
909 :
910 0 : end subroutine smooth_eps_phase_sep
911 :
912 : end module phase_separation
|
