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 :
21 : module winds
22 :
23 : use star_private_def
24 : use const_def, only: dp, ln10, rsun, msun, lsun, clight, crad, pi, pi4, secyer
25 : use chem_def, only: ih1, ihe4, ic12, ic13, in14, io16
26 : use utils_lib, only: is_bad
27 :
28 : implicit none
29 :
30 : private
31 : public :: set_mdot
32 :
33 : contains
34 :
35 11 : subroutine set_mdot(s, L_phot, M_phot, T_phot, ierr)
36 : use chem_def
37 : type (star_info), pointer :: s
38 : real(dp), intent(in) :: L_phot, M_phot, T_phot ! photosphere values (cgs)
39 : integer, intent(out) :: ierr
40 : include 'formats'
41 : ierr = 0
42 11 : call do_set_mdot(s, L_phot, M_phot, T_phot, ierr)
43 11 : if (ierr /= 0) return
44 11 : if (s% use_other_adjust_mdot) then
45 0 : call s% other_adjust_mdot(s% id, ierr)
46 0 : if (ierr /= 0) then
47 0 : if (s% report_ierr) write(*, *) 'other_adjust_mdot'
48 0 : return
49 : end if
50 : end if
51 11 : end subroutine set_mdot
52 :
53 :
54 11 : subroutine do_set_mdot(s, L_phot, M_phot, T_phot, ierr)
55 11 : use chem_def
56 : type (star_info), pointer :: s
57 : real(dp), intent(in) :: L_phot, M_phot, T_phot ! photosphere values (cgs)
58 : integer, intent(out) :: ierr
59 :
60 : integer :: h1, he4, nz
61 11 : real(dp) :: wind_mdot, wind, alfa, &
62 11 : X, Y, Z, w1, w2, T_high, T_low, L1, M1, R1, T1, &
63 11 : beta, divisor, &
64 11 : center_h1, center_he4, surface_h1, surface_he4, mdot, xfer_ratio, &
65 11 : L_div_Ledd, full_off, full_on, max_boost, super_eddington_boost, &
66 11 : hot_wind, cool_wind, H_env_mass, H_He_env_mass, He_layer_mass
67 : character (len=strlen) :: scheme
68 : logical :: using_wind_scheme_mdot, use_other
69 : real(dp), parameter :: Zsolar = 0.019d0 ! for Vink et al formula
70 :
71 : logical, parameter :: dbg = .false.
72 :
73 : include 'formats'
74 :
75 : if (dbg) write(*,1) 'enter set_mdot mass_change', s% mass_change
76 :
77 11 : ierr = 0
78 :
79 11 : xfer_ratio = 1d0
80 :
81 11 : L1 = L_phot
82 11 : M1 = M_phot
83 11 : T1 = T_phot
84 11 : R1 = sqrt(L1/(pi*crad*clight*T1*T1*T1*T1)) ! assume L1 and T1 for photosphere
85 :
86 11 : h1 = s% net_iso(ih1)
87 11 : he4 = s% net_iso(ihe4)
88 11 : nz = s% nz
89 11 : wind_mdot = 0
90 11 : using_wind_scheme_mdot = .false.
91 :
92 11 : call eval_super_eddington_wind(s, L1, M1, R1, ierr)
93 11 : if (ierr /= 0) then
94 0 : if (dbg .or. s% report_ierr) write(*, *) 'set_mdot: eval_super_eddington_wind ierr'
95 0 : return
96 : end if
97 :
98 11 : if (s% super_eddington_wind_mdot > wind_mdot) then
99 0 : wind_mdot = s% super_eddington_wind_mdot
100 : if (dbg) write(*,1) 'super eddington wind lg_Mdot', &
101 : log10(s% super_eddington_wind_mdot/(Msun/secyer))
102 : end if
103 :
104 11 : mdot = eval_rlo_wind(s, L1/Lsun, R1/Rsun, T1, xfer_ratio, ierr) ! Msun/year
105 11 : mdot = mdot*Msun/secyer
106 11 : if (ierr /= 0) then
107 0 : if (dbg .or. s% report_ierr) write(*, *) 'set_mdot: eval_rlo_wind ierr'
108 0 : return
109 : end if
110 11 : s% doing_rlo_wind = (mdot /= 0)
111 : if (dbg) write(*,*) 's% doing_rlo_wind', s% doing_rlo_wind, mdot, wind_mdot
112 :
113 11 : if (s% doing_rlo_wind .and. mdot > wind_mdot) then
114 11 : wind_mdot = mdot
115 : if (dbg) write(*,1) 's% doing_rlo_wind mdot', wind_mdot
116 : end if
117 :
118 11 : if (h1 > 0) then
119 11 : center_h1 = s% xa(h1,nz)
120 11 : surface_h1 = s% xa(h1,1)
121 : else
122 0 : center_h1 = 0
123 0 : surface_h1 = 0
124 : end if
125 11 : if (he4 > 0) then
126 11 : center_he4 = s% xa(he4,nz)
127 11 : surface_he4 = s% xa(he4,1)
128 : else
129 0 : center_he4 = 0
130 0 : surface_he4 = 0
131 : end if
132 :
133 11 : if ((s% mass_change > 0 .and. wind_mdot == 0) .or. &
134 : (s% mass_change < 0 .and. -s% mass_change*Msun/secyer > wind_mdot)) then
135 : if (dbg) write(*,*) 'mass_change mdot', s% mass_change
136 0 : wind_mdot = -s% mass_change*Msun/secyer
137 0 : if (s% mass_change > 0 .and. xfer_ratio < 1d0) then
138 0 : write(*,1) 'almost full roche lobe: reduce mdot by xfer_ratio', xfer_ratio
139 0 : wind_mdot = wind_mdot*xfer_ratio
140 : end if
141 : end if
142 :
143 : !separate calls for hot and cool winds
144 11 : if(s% hot_wind_full_on_T < s% cool_wind_full_on_T)then
145 0 : ierr = -1
146 0 : write(*,*) ' *** set_mdot error: hot_wind_full_on_T < cool_wind_full_on_T '
147 0 : return
148 : end if
149 :
150 11 : if(T1 >= s% cool_wind_full_on_T)then !do hot_wind calculation
151 11 : call eval_wind_for_scheme(s% hot_wind_scheme,hot_wind)
152 : else
153 0 : hot_wind = 0d0
154 : end if
155 :
156 11 : if(T1 <= s% hot_wind_full_on_T)then
157 0 : if (center_h1 < 0.01d0 .and. center_he4 < s% RGB_to_AGB_wind_switch) then
158 0 : scheme = s% cool_wind_AGB_scheme
159 : if (dbg) &
160 : write(*,1) 'using cool_wind_AGB_scheme: "' // trim(scheme) // '"', &
161 0 : center_h1, center_he4, s% RGB_to_AGB_wind_switch
162 : else
163 0 : scheme = s% cool_wind_RGB_scheme
164 : if (dbg) &
165 : write(*,*) 'using cool_wind_RGB_scheme: "' // trim(scheme) // '"'
166 : end if
167 0 : call eval_wind_for_scheme(scheme, cool_wind)
168 : else
169 11 : cool_wind = 0d0
170 : end if
171 :
172 : !now combine the contributions of hot and cool winds
173 11 : if(T1 >= s% hot_wind_full_on_T)then
174 11 : wind = hot_wind
175 0 : else if(T1 <= s% cool_wind_full_on_T)then
176 0 : wind = cool_wind
177 0 : else if(s% hot_wind_full_on_T == s% cool_wind_full_on_T)then
178 0 : wind = 0.5d0*(hot_wind + cool_wind)
179 : else ! blend
180 0 : divisor = s% hot_wind_full_on_T - s% cool_wind_full_on_T
181 0 : beta = min( (s% hot_wind_full_on_T - T1) / divisor, 1d0)
182 0 : alfa = 1d0 - beta
183 0 : wind = alfa*hot_wind + beta*cool_wind
184 : end if
185 :
186 11 : if (wind*Msun/secyer > abs(wind_mdot)) then
187 : using_wind_scheme_mdot = .true.
188 : if (dbg) write(*,1) 'use wind scheme mdot', wind*Msun/secyer, wind_mdot
189 : wind_mdot = wind*Msun/secyer
190 : end if
191 :
192 : if (dbg) write(*,1) 'wind_mdot 1', wind_mdot
193 :
194 : if (using_wind_scheme_mdot) then
195 0 : if (s% no_wind_if_no_rotation .and. .not. s% rotation_flag) then
196 0 : s% mstar_dot = 0
197 0 : if (s% trace_dt_control_mass_change) &
198 0 : write(*,1) 'no_wind_if_no_rotation'
199 0 : return
200 : end if
201 0 : if (s% dt > 0 .and. s% dt < s% mass_change_full_on_dt) then
202 0 : if (s% dt <= s% mass_change_full_off_dt) then
203 0 : s% mstar_dot = 0
204 0 : if (s% trace_dt_control_mass_change .or. dbg) &
205 0 : write(*,1) 'no wind: dt <= mass_change_full_off_dt'
206 0 : return
207 : end if
208 : alfa = (s% dt - s% mass_change_full_off_dt)/ &
209 0 : (s% mass_change_full_on_dt - s% mass_change_full_off_dt)
210 0 : if (s% trace_dt_control_mass_change .or. dbg) &
211 0 : write(*,1) 'reduce wind: dt <= mass_change_full_on_dt', alfa
212 0 : wind_mdot = wind_mdot*alfa
213 : end if
214 : end if
215 :
216 11 : if (wind_mdot >= 0 .and. s% super_eddington_scaling_factor <= 0) then
217 : ! check for super eddington boost to wind
218 11 : L_div_Ledd = L1 / s% prev_Ledd
219 11 : full_off = s% wind_boost_full_off_L_div_Ledd
220 11 : if (L_div_Ledd > full_off) then
221 0 : full_on = s% wind_boost_full_on_L_div_Ledd
222 0 : max_boost = s% super_eddington_wind_max_boost
223 0 : if (L_div_Ledd >= full_on) then
224 0 : super_eddington_boost = max_boost
225 : else
226 : super_eddington_boost = &
227 0 : 1 + (max_boost-1)*(L_div_Ledd - full_off)/(full_on - full_off)
228 : end if
229 0 : wind_mdot = wind_mdot*super_eddington_boost
230 0 : if (s% trace_super_eddington_wind_boost .or. dbg) then
231 0 : write(*,1) 'super eddington wind boost factor, L_div_Ledd', &
232 0 : super_eddington_boost, L_div_Ledd
233 0 : write(*,'(A)')
234 : end if
235 : end if
236 : end if
237 :
238 : if (dbg) write(*,1) 'wind_mdot 2', wind_mdot
239 :
240 11 : if (wind_mdot >= 0 .and. s% min_wind > 0 .and. &
241 : wind_mdot < s% min_wind*Msun/secyer) then
242 : if (dbg) write(*,1) 'use s% min_wind', s% min_wind
243 11 : wind_mdot = s% min_wind*Msun/secyer
244 : end if
245 :
246 : if (dbg) write(*,1) 'wind_mdot 3', wind_mdot
247 :
248 11 : if (wind_mdot >= 0 .and. s% max_wind > 0 .and. &
249 : wind_mdot > s% max_wind*Msun/secyer) then
250 : if (dbg) write(*,1) 'use s% max_wind', s% max_wind
251 11 : wind_mdot = s% max_wind*Msun/secyer
252 : end if
253 : if (dbg) write(*,1) 'wind_mdot 4', wind_mdot
254 :
255 : if (wind_mdot >= 0) then
256 11 : if (s% starting_T_center > s% max_T_center_for_any_mass_loss) then
257 : if (dbg) write(*,1) 'starting_T_center > max_T_center_for_any_mass_loss', &
258 : s% starting_T_center, s% max_T_center_for_any_mass_loss
259 : wind_mdot = 0
260 11 : else if (s% starting_T_center > s% max_T_center_for_full_mass_loss) then
261 : if (dbg) write(*,1) 'starting_T_center > max_T_center_for_full_mass_loss', &
262 : s% starting_T_center, s% max_T_center_for_full_mass_loss
263 : wind_mdot = wind_mdot* &
264 : (s% max_T_center_for_any_mass_loss - s% starting_T_center)/ &
265 : (s% max_T_center_for_any_mass_loss - &
266 0 : s% max_T_center_for_full_mass_loss)
267 : end if
268 : end if
269 : if (dbg) write(*,1) 'wind_mdot 5', wind_mdot
270 :
271 11 : if (wind_mdot >= 0) then
272 11 : H_env_mass = s% star_mass - s% he_core_mass
273 11 : H_He_env_mass = s% star_mass - s% co_core_mass
274 11 : He_layer_mass = s% he_core_mass - s% co_core_mass
275 11 : if (s% wind_H_envelope_limit > 0 .and. &
276 : H_env_mass < s% wind_H_envelope_limit) then
277 : wind_mdot = 0
278 11 : else if (s% wind_H_He_envelope_limit > 0 .and. &
279 : H_He_env_mass < s% wind_H_He_envelope_limit) then
280 : wind_mdot = 0
281 11 : else if (s% wind_He_layer_limit > 0 .and. &
282 : He_layer_mass < s% wind_He_layer_limit) then
283 0 : wind_mdot = 0
284 : end if
285 : end if
286 :
287 11 : s% mstar_dot = -wind_mdot
288 : if (dbg) write(*,1) 'mstar_dot', s% mstar_dot
289 11 : if (s% mstar_dot < 0 .and. &
290 : (s% min_wind > 0 .or. &
291 : using_wind_scheme_mdot .or. &
292 : s% v_div_v_crit_avg_surf > 0.8d0)) then
293 0 : call rotation_enhancement(ierr)
294 0 : if (is_bad(s% rotational_mdot_boost)) then
295 0 : write(*,2) 'is_bad(s% rotational_mdot_boost)', s% model_number
296 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'winds: rotation_enhancement')
297 : end if
298 0 : if (ierr /= 0) then
299 0 : if (dbg .or. s% report_ierr) write(*, *) 'set_mdot: rotation_enhancement ierr'
300 0 : return
301 : end if
302 : end if
303 :
304 11 : s% explicit_mstar_dot = s% mstar_dot
305 :
306 11 : if (dbg) then
307 : write(*,1) 'final star_mdot', s% mstar_dot/(Msun/secyer)
308 : write(*,1) 'final lg abs s% mstar_dot/(Msun/secyer)', safe_log10(abs(s% mstar_dot/(Msun/secyer)))
309 : write(*,'(A)')
310 : end if
311 :
312 : contains
313 :
314 11 : subroutine eval_wind_for_scheme(scheme,wind)
315 : character(len=strlen) :: scheme
316 : real(dp), intent(out) :: wind
317 : include 'formats'
318 :
319 11 : use_other = (s% use_other_wind .or. scheme == 'other')
320 11 : if ((.not. use_other) .and. len_trim(scheme) == 0) then
321 11 : wind = 0
322 : else
323 0 : wind = 4d-13*(L1*R1/M1)/(Lsun*Rsun/Msun) ! in Msun/year
324 : if (dbg) write(*,1) 'wind', wind
325 0 : if (wind <= 0 .or. is_bad_num(wind)) then
326 0 : ierr = -1
327 0 : write(*,*) 'bad value for wind :', wind,L1,R1,M1
328 : if (dbg) call mesa_error(__FILE__,__LINE__,'debug: bad value for wind')
329 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'winds')
330 0 : return
331 : end if
332 0 : X = surface_h1
333 0 : Y = surface_he4
334 : ! use the Zbase as a proxy for Fe-group abundances
335 : ! unlikely to change during the stellar evolution
336 0 : if (s% kap_rq% Zbase >= 0) then
337 0 : Z = s% kap_rq% Zbase
338 : else
339 : if (dbg) write(*,*) 'Zbase not set, using present-day Z for wind'
340 0 : Z = 1 - (X + Y)
341 : end if
342 :
343 0 : if (use_other) then
344 : if (dbg) write(*,*) 'call other_wind'
345 0 : call s% other_wind(s% id, L1, M1, R1, T1, X, Y, Z, wind, ierr)
346 : if (ierr /= 0) return
347 0 : else if (scheme == 'Dutch') then
348 0 : T_high = 11000
349 0 : T_low = 10000
350 0 : if (s% Dutch_scaling_factor == 0) then
351 0 : wind = 0
352 0 : else if (T1 <= T_low) then
353 0 : call eval_lowT_Dutch(wind)
354 0 : else if (T1 >= T_high) then
355 0 : call eval_highT_Dutch(wind)
356 : else ! transition
357 : call eval_lowT_Dutch(w1)
358 : call eval_highT_Dutch(w2)
359 0 : alfa = (T1 - T_low)/(T_high - T_low)
360 0 : wind = (1-alfa)*w1 + alfa*w2
361 : end if
362 0 : wind = s% Dutch_scaling_factor * wind
363 0 : else if (scheme == 'Reimers') then
364 0 : wind = wind * s% Reimers_scaling_factor
365 : if (dbg) then
366 : write(*,1) 's% Reimers_scaling_factor', s% Reimers_scaling_factor
367 : write(*,1) 'Reimers_wind', wind
368 : write(*,1) 'L1/Lsun', L1/Lsun
369 : write(*,1) 'R1/Rsun', R1/Rsun
370 : write(*,1) 'M1/Msun', M1/Msun
371 : write(*,1) 'Reimers_scaling_factorReimers_scaling_factor', s% Reimers_scaling_factor
372 : write(*,1) 'wind', wind
373 : write(*,1) 'log10 wind', log10(wind)
374 : write(*,'(A)')
375 : call mesa_error(__FILE__,__LINE__,'debug: winds')
376 : end if
377 0 : else if (scheme == 'Vink') then
378 0 : call eval_Vink_wind(wind)
379 0 : wind = wind * s% Vink_scaling_factor
380 : if (dbg) write(*,1) 'Vink_wind', wind
381 : if (dbg) write(*,1) 'Grafener_wind', wind
382 0 : else if (scheme == 'Blocker') then
383 0 : call eval_blocker_wind(wind)
384 : if (dbg) write(*,1) 'Blocker_wind', wind
385 0 : else if (scheme == 'de Jager') then
386 0 : call eval_de_Jager_wind(wind)
387 0 : wind = s% de_Jager_scaling_factor * wind
388 : if (dbg) write(*,1) 'de_Jager_wind', wind
389 0 : else if (scheme == 'van Loon') then
390 0 : call eval_van_Loon_wind(wind)
391 0 : wind = s% van_Loon_scaling_factor * wind
392 : if (dbg) write(*,1) 'van_Loon_wind', wind
393 0 : else if (scheme == 'Nieuwenhuijzen') then
394 0 : call eval_Nieuwenhuijzen_wind(wind)
395 0 : wind = s% Nieuwenhuijzen_scaling_factor * wind
396 : if (dbg) write(*,1) 'Nieuwenhuijzen_wind', wind
397 0 : else if (scheme == 'Bjorklund') then
398 0 : call eval_Bjorklund_wind(wind)
399 0 : wind = s% Bjorklund_scaling_factor * wind
400 : if (dbg) write(*,1) 'Bjorklund_wind', wind
401 : else
402 0 : ierr = -1
403 0 : write(*,*) 'unknown name for wind scheme : ' // trim(scheme)
404 : if (dbg) call mesa_error(__FILE__,__LINE__,'debug: bad value for wind scheme')
405 0 : return
406 : end if
407 : end if
408 :
409 11 : end subroutine eval_wind_for_scheme
410 :
411 :
412 0 : subroutine rotation_enhancement(ierr)
413 : use star_utils, only: eval_kh_timescale
414 : integer, intent(out) :: ierr
415 : ! as in Heger, Langer, and Woosley, 2000, ApJ, 528:368-396. section 2.6
416 : ! Mdot = Mdot_no_rotation/(1 - Osurf/Osurf_crit)^mdot_omega_power
417 : ! where Osurf = angular velocity at surface
418 : ! Osurf_crit^2 = (1 - Gamma_edd)*G*M/R_equatorial^3
419 : ! Gamma_edd = kappa*L/(4 pi c G M), Eddington factor
420 0 : real(dp) :: enhancement, wind_mdot, wind_mdot_lim, &
421 0 : kh_timescale, wind_mdot_prev, dmsfac, dmskhf
422 :
423 : include 'formats'
424 :
425 0 : ierr = 0
426 :
427 0 : if (.not. s% rotation_flag) return
428 0 : if (s% mdot_omega_power <= 0) return
429 0 : if (s% mstar_dot >= 0) return
430 :
431 0 : wind_mdot = -s% mstar_dot
432 :
433 0 : kh_timescale = eval_kh_timescale(s% cgrav(1), M1, R1, L1)
434 0 : dmskhf = s% rotational_mdot_kh_fac
435 0 : dmsfac = s% rotational_mdot_boost_fac
436 0 : wind_mdot_lim = min(dmskhf*M1/kh_timescale, wind_mdot*dmsfac)
437 :
438 0 : enhancement = pow(max(1d-22, 1d0 - s% v_div_v_crit_avg_surf), -s% mdot_omega_power)
439 0 : if (s% max_rotational_mdot_boost > 0 .and. &
440 : enhancement > s% max_rotational_mdot_boost) then
441 0 : enhancement = s% max_rotational_mdot_boost
442 : end if
443 :
444 : if (enhancement > s% lim_trace_rotational_mdot_boost) then
445 : if (dbg) write(*,1) &
446 : 'mdot rotation enhancement factor for mdot, v_div_v_crit_avg_surf', &
447 : enhancement, s% v_div_v_crit_avg_surf
448 : end if
449 :
450 0 : if (wind_mdot*enhancement < wind_mdot_lim) then
451 : wind_mdot = wind_mdot*enhancement
452 : if (dbg) write(*,2) 'wind_mdot = wind_mdot*enhancement', &
453 : s% model_number, enhancement, &
454 : log10(wind_mdot/(Msun/secyer)), log10(wind_mdot_lim/(Msun/secyer))
455 : else
456 0 : enhancement = wind_mdot_lim/wind_mdot
457 0 : wind_mdot = wind_mdot_lim
458 : if (dbg) write(*,2) 'wind_mdot = wind_mdot_lim', &
459 : s% model_number, log10(wind_mdot/(Msun/secyer))
460 : end if
461 :
462 0 : wind_mdot_prev = -s% mstar_dot_old
463 0 : if (wind_mdot > 0 .and. wind_mdot_prev > 0) &
464 0 : wind_mdot = min(wind_mdot, s% max_mdot_jump_for_rotation*wind_mdot_prev)
465 :
466 : ! recheck max_wind
467 0 : if (wind_mdot >= 0 .and. s% max_wind > 0 .and. &
468 : wind_mdot > s% max_wind*Msun/secyer) then
469 : if (dbg) write(*,1) 'use s% max_wind', s% max_wind
470 0 : wind_mdot = s% max_wind*Msun/secyer
471 : end if
472 :
473 0 : s% mstar_dot = -wind_mdot
474 0 : s% rotational_mdot_boost = enhancement
475 :
476 0 : end subroutine rotation_enhancement
477 :
478 :
479 0 : subroutine eval_Bjorklund_wind(w)
480 : real(dp), intent(inout) :: w
481 : real(dp), parameter :: Zbjork = 0.014d0
482 0 : real(dp) :: logw, Meff
483 :
484 : ! electron opacity is constant 0.34 in their models (Eq. 6)
485 0 : Meff = M1*(1d0 - 0.34d0*L1/(pi4*clight*s% cgrav(1)*M1)) ! effective mass
486 :
487 : ! eq 7 from Björklund et al, 2022, arXiv:2203.08218, accepted to A&A
488 : logw = - 5.52d0 &
489 : + 2.39d0 * log10(L1/(1d6*Lsun)) &
490 : - 1.48d0 * log10(Meff/(4.5d1*Msun)) &
491 : + 2.12d0 * log10(T1/4.5d4) &
492 0 : + (0.75d0 - 1.87d0 * log10(T1/4.5d4)) * log10(Z/Zbjork)
493 0 : w = exp10(logw)
494 :
495 0 : end subroutine eval_Bjorklund_wind
496 :
497 0 : subroutine eval_Vink_wind(w)
498 : real(dp), intent(inout) :: w
499 0 : real(dp) :: alfa, w1, w2, Teff_jump, logMdot, dT, vinf_div_vesc
500 :
501 : ! alfa = 1 for hot side, = 0 for cool side
502 0 : if (T1 > 27500d0) then
503 : alfa = 1
504 0 : else if (T1 < 22500d0) then
505 : alfa = 0
506 : else ! use Vink et al 2001, eqns 14 and 15 to set "jump" temperature
507 0 : Teff_jump = 1d3*(61.2d0 + 2.59d0*(-13.636d0 + 0.889d0*log10(Z/Zsolar)))
508 0 : dT = 100d0
509 0 : if (T1 > Teff_jump + dT) then
510 : alfa = 1
511 0 : else if (T1 < Teff_jump - dT) then
512 : alfa = 0
513 : else
514 0 : alfa = 0.5d0*(T1 - (Teff_jump - dT)) / dT
515 : end if
516 : end if
517 :
518 0 : if (alfa > 0) then ! eval hot side wind (eqn 24)
519 0 : vinf_div_vesc = 2.6d0 ! this is the hot side galactic value
520 0 : vinf_div_vesc = vinf_div_vesc*pow(Z/Zsolar,0.13d0) ! corrected for Z
521 : logMdot = &
522 : - 6.697d0 &
523 : + 2.194d0*log10(L1/Lsun/1d5) &
524 : - 1.313d0*log10(M1/Msun/30d0) &
525 : - 1.226d0*log10(vinf_div_vesc/2d0) &
526 : + 0.933d0*log10(T1/4d4) &
527 : - 10.92d0*pow2(log10(T1/4d4)) &
528 0 : + 0.85d0*log10(Z/Zsolar)
529 0 : w1 = exp10(logMdot)
530 : else
531 : w1 = 0
532 : end if
533 :
534 0 : if (alfa < 1) then ! eval cool side wind (eqn 25)
535 0 : vinf_div_vesc = 1.3d0 ! this is the cool side galactic value
536 0 : vinf_div_vesc = vinf_div_vesc*pow(Z/Zsolar,0.13d0) ! corrected for Z
537 : logMdot = &
538 : - 6.688d0 &
539 : + 2.210d0*log10(L1/Lsun/1d5) &
540 : - 1.339d0*log10(M1/Msun/30d0) &
541 : - 1.601d0*log10(vinf_div_vesc/2d0) &
542 : + 1.07d0*log10(T1/2d4) &
543 0 : + 0.85d0*log10(Z/Zsolar)
544 0 : w2 = exp10(logMdot)
545 : else
546 : w2 = 0
547 : end if
548 :
549 0 : w = alfa*w1 + (1 - alfa)*w2
550 :
551 : if (dbg) write(*,*) 'vink wind', w
552 :
553 0 : end subroutine eval_Vink_wind
554 :
555 :
556 0 : subroutine eval_blocker_wind(w)
557 : real(dp), intent(inout) :: w
558 : w = w * s% Blocker_scaling_factor * &
559 0 : 4.83d-9 * pow(M1/Msun,-2.1d0) * pow(L1/Lsun,2.7d0)
560 : if (dbg) write(*,*) 'blocker wind', w
561 0 : end subroutine eval_blocker_wind
562 :
563 :
564 0 : subroutine eval_highT_Dutch(w)
565 : real(dp), intent(out) :: w
566 : include 'formats'
567 0 : if (surface_h1 < 0.4d0) then ! helium rich Wolf-Rayet star: Nugis & Lamers
568 0 : w = 1d-11 * pow(L1/Lsun,1.29d0) * pow(Y,1.7d0) * sqrt(Z)
569 : if (dbg) write(*,1) 'Dutch_wind = Nugis & Lamers', log10(wind)
570 : else
571 0 : call eval_Vink_wind(w)
572 : end if
573 0 : end subroutine eval_highT_Dutch
574 :
575 :
576 0 : subroutine eval_lowT_Dutch(w)
577 : real(dp), intent(out) :: w
578 : include 'formats'
579 0 : if (s% Dutch_wind_lowT_scheme == 'de Jager') then
580 0 : call eval_de_Jager_wind(w)
581 : if (dbg) write(*,1) 'Dutch_wind = de Jager', safe_log10(wind), T1, T_low, T_high
582 0 : else if (s% Dutch_wind_lowT_scheme == 'van Loon') then
583 0 : call eval_van_Loon_wind(w)
584 : if (dbg) write(*,1) 'Dutch_wind = van Loon', safe_log10(wind), T1, T_low, T_high
585 0 : else if (s% Dutch_wind_lowT_scheme == 'Nieuwenhuijzen') then
586 0 : call eval_Nieuwenhuijzen_wind(w)
587 : if (dbg) write(*,1) 'Dutch_wind = Nieuwenhuijzen', safe_log10(wind), T1, T_low, T_high
588 : else
589 : write(*,*) 'unknown value for Dutch_wind_lowT_scheme ' // &
590 0 : trim(s% Dutch_wind_lowT_scheme)
591 0 : w = 0
592 : end if
593 0 : end subroutine eval_lowT_Dutch
594 :
595 :
596 0 : subroutine eval_de_Jager_wind(w)
597 : ! de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259.
598 : real(dp), intent(out) :: w
599 0 : real(dp) :: log10w
600 : include 'formats'
601 0 : log10w = 1.769d0*log10(L1/Lsun) - 1.676d0*log10(T1) - 8.158d0
602 0 : w = exp10(log10w)
603 : if (dbg) then
604 : write(*,1) 'de_Jager log10 wind', log10w
605 : end if
606 0 : end subroutine eval_de_Jager_wind
607 :
608 :
609 0 : subroutine eval_van_Loon_wind(w)
610 : ! van Loon et al. 2005, A&A, 438, 273
611 : real(dp), intent(out) :: w
612 0 : real(dp) :: log10w
613 : include 'formats'
614 0 : log10w = -5.65d0 + 1.05d0*log10(L1/(1d4*Lsun)) - 6.3d0*log10(T1/35d2)
615 0 : w = exp10(log10w)
616 0 : end subroutine eval_van_Loon_wind
617 :
618 :
619 0 : subroutine eval_Nieuwenhuijzen_wind(w)
620 : ! Nieuwenhuijzen, H.; de Jager, C. 1990, A&A, 231, 134 (eqn 2)
621 : real(dp), intent(out) :: w
622 0 : real(dp) :: log10w
623 : include 'formats'
624 : log10w = -14.02d0 + &
625 : 1.24d0*log10(L1/Lsun) + &
626 : 0.16d0*log10(M1/Msun) + &
627 0 : 0.81d0*log10(R1/Rsun)
628 0 : w = exp10(log10w)
629 : if (dbg) then
630 : write(*,1) 'Nieuwenhuijzen log10 wind', log10w
631 : end if
632 0 : end subroutine eval_Nieuwenhuijzen_wind
633 :
634 : end subroutine do_set_mdot
635 :
636 :
637 11 : subroutine eval_super_eddington_wind(s, L, M, R, ierr)
638 : type (star_info), pointer :: s
639 : real(dp), intent(in) :: L, M, R
640 : integer, intent(out) :: ierr
641 :
642 11 : real(dp) :: Ledd, Leff, vesc2
643 : include 'formats'
644 :
645 11 : ierr = 0
646 11 : s% super_eddington_wind_mdot = 0
647 11 : if (s% super_eddington_scaling_factor <= 0) return
648 :
649 0 : Ledd = s% prev_Ledd
650 0 : Leff = L/s% super_eddington_wind_Ledd_factor
651 0 : if (Leff <= Ledd) return
652 0 : vesc2 = 2d0 * s% cgrav(1)*M/R
653 0 : s% super_eddington_wind_mdot = s% super_eddington_scaling_factor*(Leff - Ledd)/(0.5d0 * vesc2)
654 0 : if (mod(s% model_number, s% terminal_interval) == 0) &
655 0 : write(*,'(a60,i12,1p2e12.4)') 'super eddington wind: lg_Mdot, L/Ledd', &
656 0 : s% model_number, log10(s% super_eddington_wind_mdot/(Msun/secyer)), L/Ledd
657 :
658 : end subroutine eval_super_eddington_wind
659 :
660 :
661 11 : real(dp) function eval_rlo_wind(s, L_surf, R, Teff, xfer_ratio, ierr) ! value in Msun/year
662 : type (star_info), pointer :: s
663 : real(dp), intent(in) :: L_surf, R, Teff ! Lsun, Rsun, K
664 : real(dp), intent(inout) :: xfer_ratio
665 : integer, intent(out) :: ierr
666 11 : real(dp) :: roche_lobe_radius ! Rsun
667 11 : real(dp) :: ratio, scale_height, mdot
668 : include 'formats'
669 11 : ierr = 0
670 11 : eval_rlo_wind = 0
671 11 : if (s% rlo_scaling_factor <= 0) return
672 0 : if (L_surf < s% rlo_wind_min_L) return
673 0 : if (Teff > s% rlo_wind_max_Teff) return
674 0 : scale_height = s% rlo_wind_scale_height
675 0 : if (scale_height <= 0) then
676 0 : scale_height = s% Peos(1) / (s% cgrav(1)*s% m(1)*s% rho(1) / (s% r(1)**2)) / Rsun
677 : end if
678 0 : roche_lobe_radius = s% rlo_wind_roche_lobe_radius
679 0 : ratio = R/roche_lobe_radius
680 : !write(*,2) 'R/roche_lobe_radius', s% model_number, ratio
681 0 : if (ratio < 1) then
682 : ! check for reduction in transfer ratio for almost full Roche lobe
683 0 : if (ratio < s% roche_lobe_xfer_full_on) return
684 0 : if (ratio > s% roche_lobe_xfer_full_off) then
685 0 : xfer_ratio = 0
686 0 : return
687 : end if
688 : xfer_ratio = (s% roche_lobe_xfer_full_off - ratio) / &
689 0 : (s% roche_lobe_xfer_full_off - s% roche_lobe_xfer_full_on)
690 0 : xfer_ratio = 0.5d0*(1 - cospi(xfer_ratio))
691 0 : return
692 : end if
693 : mdot = s% rlo_wind_base_mdot* &
694 0 : exp(min(6*ln10,(R - roche_lobe_radius)/scale_height))
695 0 : eval_rlo_wind = s% rlo_scaling_factor*mdot ! Msun/year
696 :
697 : !write(*,1) 's% rlo_wind_base_mdot', s% rlo_wind_base_mdot
698 : !write(*,1) 'R - roche_lobe_radius', R - roche_lobe_radius, R, roche_lobe_radius
699 : !write(*,1) 'scale_height', scale_height
700 : !write(*,1) 's% rlo_scaling_factor', s% rlo_scaling_factor
701 : !write(*,1) 'eval_rlo_wind, log eval_rlo_wind, R/R_L', log10(eval_rlo_wind), R/roche_lobe_radius
702 :
703 0 : end function eval_rlo_wind
704 :
705 :
706 : end module winds
|
