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 pre_ms_model
21 :
22 : use star_private_def
23 : use const_def, only: dp, pi, pi4, ln10, clight, standard_cgrav, crad, lsun, msun, rsun, one_third, two_thirds, four_thirds_pi
24 :
25 : implicit none
26 :
27 : private
28 : public :: build_pre_ms_model
29 :
30 : logical, parameter :: dbg = .false.
31 :
32 : contains
33 :
34 0 : subroutine build_pre_ms_model(id, s, nvar_hydro, species, ierr)
35 : use chem_def
36 : use chem_lib, only: basic_composition_info, chem_Xsol
37 : use adjust_xyz, only: get_xa_for_standard_metals
38 : use alloc, only: allocate_star_info_arrays
39 : use num_lib, only: look_for_brackets, safe_root
40 :
41 : type (star_info), pointer :: s
42 : integer, intent(in) :: id, nvar_hydro, species
43 : integer, intent(out) :: ierr
44 :
45 : real(dp) :: &
46 0 : initial_z, x, y, z, xa(species), mstar, mstar1, lgM, rstar, rho_c, &
47 0 : abar, zbar, z53bar, mass_correction, z2bar, ye, sumx, &
48 0 : lgL, eps_grav, lnd, dlnd, lnd1, lnd3, y1, y3, epsx, epsy, &
49 0 : T_c, guess_rho_c, d_log10_P
50 : integer :: i, j, k, nz, pre_ms_lrpar
51 0 : real(dp), pointer :: xh(:,:), q(:), dq(:)
52 : integer, parameter :: pre_ms_lipar = 1, imax = 100, rpar_init = 8
53 : integer, target :: ipar_ary(pre_ms_lipar)
54 : integer, pointer :: ipar(:)
55 0 : real(dp), target :: rpar_ary(rpar_init+species) ! (pre_ms_lrpar)
56 : real(dp), pointer :: rpar(:)
57 0 : real(dp) :: mu_eff ! effective mean molecular weight for gas particles
58 0 : real(dp) :: initial_y, initial_h1, initial_h2, initial_he3, initial_he4, &
59 0 : xsol_he3, xsol_he4
60 : integer :: initial_zfracs
61 : real(dp), parameter :: max_mass_to_create = 90, min_mass_to_create = 0.03d0
62 :
63 : include 'formats'
64 :
65 0 : ipar => ipar_ary
66 0 : rpar => rpar_ary
67 :
68 0 : pre_ms_lrpar = rpar_init+species
69 :
70 0 : if (nvar_hydro > 4) then
71 0 : write(*,*) 'sorry, build_pre_ms_model only supports the basic 4 vars.'
72 0 : ierr = -1
73 0 : return
74 : end if
75 0 : mstar = min(max_mass_to_create*Msun, max(min_mass_to_create*Msun, s% mstar))
76 0 : s% mstar = mstar
77 0 : s% star_mass = mstar/Msun
78 0 : s% xmstar = mstar
79 :
80 0 : T_c = s% pre_ms_T_c
81 0 : if (T_c <= 0) T_c = 9d5
82 0 : if (T_c >= 1d6) then
83 0 : write(*,1) 'log T_c', log10(T_c)
84 0 : write(*,*) 'center temperature is too high for a pre main sequence model.'
85 0 : write(*,*) 'please pick T_c below 10^6 so code can ignore nuclear reactions.'
86 0 : ierr = -1
87 0 : return
88 : end if
89 :
90 0 : ierr = 0
91 0 : initial_z = s% initial_z
92 0 : initial_y = s% initial_y
93 0 : s% M_center = 0
94 0 : s% L_center = 0
95 0 : s% R_center = 0
96 0 : s% v_center = 0
97 :
98 0 : initial_h1 = max(0d0, min(1d0, 1d0 - (initial_z + initial_y)))
99 0 : initial_h2 = 0
100 0 : if (s% initial_he3 < 0d0) then
101 0 : xsol_he3 = chem_Xsol('he3')
102 0 : xsol_he4 = chem_Xsol('he4')
103 0 : initial_he3 = initial_y*xsol_he3/(xsol_he3 + xsol_he4)
104 0 : initial_he4 = initial_y*xsol_he4/(xsol_he3 + xsol_he4)
105 0 : else if (s% initial_he3 < initial_y) then
106 0 : initial_he3 = s% initial_he3
107 0 : initial_he4 = initial_y - s% initial_he3
108 : else
109 0 : write(*,*) "ERROR: initial_he3 is larger than initial_y"
110 : ierr = -1
111 : end if
112 0 : initial_zfracs = s% pre_ms_initial_zfracs
113 :
114 : call get_xa_for_standard_metals(s, &
115 : species, s% chem_id, s% net_iso, &
116 : initial_h1, initial_h2, initial_he3, initial_he4, &
117 : initial_zfracs, s% pre_ms_dump_missing_heaviest, &
118 0 : xa, ierr)
119 0 : if (ierr /= 0) return
120 :
121 : if (dbg) then
122 : write(*,*) 'abundances'
123 : do i=1,species
124 : write(*,1) trim(chem_isos% name(s% chem_id(i))), xa(i)
125 : end do
126 : write(*,'(A)')
127 : end if
128 :
129 0 : if (abs(1-sum(xa(:))) > 1d-8) then
130 0 : write(*,1) 'initial_h1', initial_h1
131 0 : write(*,1) 'initial_h2', initial_h2
132 0 : write(*,1) 'initial_he3', initial_he3
133 0 : write(*,1) 'initial_he4', initial_he4
134 0 : write(*,1) 'initial_y', initial_y
135 0 : write(*,1) 'initial_z', initial_z
136 0 : write(*,1) 'initial_h1+h2+he3+he4+z', &
137 0 : initial_h1 + initial_h2 + initial_he3 + initial_he4 + initial_z
138 0 : write(*,1) 'sum(xa(:))', sum(xa(:))
139 0 : write(*,*) 'build_pre_ms_model'
140 0 : ierr = -1
141 0 : return
142 : end if
143 :
144 : call basic_composition_info( &
145 : species, s% chem_id, xa(:), x, y, z, abar, zbar, z2bar, z53bar, ye, &
146 0 : mass_correction, sumx)
147 :
148 0 : mu_eff = 4 / (3 + 5*x)
149 : ! estimate mu_eff assuming complete ionization and Z << 1
150 :
151 0 : guess_rho_c = s% pre_ms_guess_rho_c
152 0 : if (guess_rho_c <= 0) then ! use n=3/2 polytrope
153 0 : rstar = Rsun*7.41d6*(mu_eff/0.6d0)*(mstar/Msun)/T_c ! Ushomirsky et al, 6
154 0 : rho_c = 8.44d0*(mstar/Msun)*pow3(Rsun/rstar) ! Ushomirsky et al, 5
155 0 : rho_c = 1.2d0*rho_c ! polytrope value is usually too small
156 : else
157 0 : rho_c = guess_rho_c
158 : end if
159 :
160 : ! pick a luminosity that is above the zams level
161 0 : lgM = log10(mstar/Msun)
162 0 : if (lgM > 1) then
163 0 : lgL = 4.5d0 + 2.5d0*(lgM - 1)
164 0 : else if (lgM > 0.5d0) then
165 0 : lgL = 3d0 + 3*(lgM - 0.5d0)
166 0 : else if (lgM > 0) then
167 0 : lgL = 4d0*lgM + 0.5d0
168 : else
169 0 : lgL = 0.5d0
170 : end if
171 :
172 : ! use uniform eps_grav to give that luminosity
173 0 : eps_grav = exp10(lgL)*Lsun/mstar
174 :
175 : if (dbg) then
176 : write(*,1) 'initial_z', initial_z
177 : write(*,1) 'T_c', T_c
178 : write(*,1) 'rho_c', rho_c
179 : write(*,1) 'lgL', lgL
180 : write(*,1) 'eps_grav', eps_grav
181 : write(*,1) 'mstar/Msun', mstar/Msun
182 : end if
183 :
184 0 : if (ASSOCIATED(s% xh)) deallocate(s% xh)
185 0 : nullify(s% q)
186 0 : nullify(s% dq)
187 :
188 0 : d_log10_P = s% pre_ms_d_log10_P
189 :
190 : i = 1 ! rpar(1) for mstar result
191 0 : rpar(i+1) = T_c; i = i+1
192 0 : rpar(i+1) = eps_grav; i = i+1
193 0 : rpar(i+1) = x; i = i+1
194 0 : rpar(i+1) = initial_z; i = i+1
195 0 : rpar(i+1) = abar; i = i+1
196 0 : rpar(i+1) = zbar; i = i+1
197 0 : rpar(i+1) = d_log10_P; i = i+1
198 :
199 0 : rpar(i+1:i+species) = xa(1:species); i = i+species
200 :
201 : if (i /= pre_ms_lrpar) then
202 : write(*,*) 'i /= pre_ms_lrpar', i, pre_ms_lrpar
203 : write(*,*) 'pre ms'
204 : ierr = -1
205 : return
206 : end if
207 :
208 0 : ipar(1) = id
209 :
210 0 : lnd = log(rho_c)
211 0 : dlnd = 0.01d0
212 :
213 : call look_for_brackets(lnd, dlnd, lnd1, lnd3, pre_ms_f, y1, y3, &
214 0 : imax, pre_ms_lrpar, rpar, pre_ms_lipar, ipar, ierr)
215 0 : if (ierr /= 0) then
216 : if (dbg) then
217 : if (dbg) write(*,*) 'look_for_brackets ierr', ierr
218 : write(*,1) 'lnd1', lnd1
219 : write(*,1) 'lnd3', lnd3
220 : write(*,1) 'y1', y1
221 : write(*,1) 'y3', y3
222 : end if
223 : return
224 : end if
225 :
226 0 : epsx = 1d-3 ! limit for variation in lnd
227 0 : epsy = 1d-3 ! limit for matching desired mass as fraction of total mass
228 :
229 : lnd = safe_root(pre_ms_f, lnd1, lnd3, y1, y3, imax, epsx, epsy, &
230 0 : pre_ms_lrpar, rpar, pre_ms_lipar, ipar, ierr)
231 0 : if (ierr /= 0) then
232 : if (dbg) write(*,*) 'safe_root ierr', ierr
233 : return
234 : end if
235 :
236 0 : mstar1 = rpar(1)
237 :
238 : ! these pointers are set, but none of these vars are set.
239 0 : xh => s% xh
240 0 : q => s% q
241 0 : dq => s% dq
242 0 : nz = s% nz
243 :
244 : ! this debug statement will cause a backtrace because the above vars are not set.
245 : if (dbg) then
246 : write(*,'(A)')
247 : write(*,*) 'finished pre-MS model'
248 : write(*,1) 'mstar1/Msun', mstar1/Msun
249 : write(*,1) '(mstar-mstar1)/mstar', (mstar-mstar1)/mstar
250 : write(*,*) 'nz', nz
251 : write(*,'(A)')
252 : call mesa_error(__FILE__,__LINE__,'debug: pre ms')
253 : end if
254 :
255 : ! The following deallocations deal with arrays which were
256 : ! needed in the root bracket/solve above, but will get
257 : ! overwritten during the call to allocate_star_info_arrays
258 :
259 0 : if (ASSOCIATED(s% xh_old)) deallocate(s% xh_old)
260 0 : if (ASSOCIATED(s% xh_start)) deallocate(s% xh_start)
261 :
262 0 : call allocate_star_info_arrays(s, ierr)
263 0 : if (ierr /= 0) then
264 0 : call dealloc
265 0 : return
266 : end if
267 :
268 0 : do k=1,nz
269 0 : do j=1,nvar_hydro
270 0 : s% xh(j,k) = xh(j,k)
271 : !write(*,3) trim(s% nameofvar(j)), j, k, xh(j,k)
272 : end do
273 0 : do j=1,species
274 0 : s% xa(j,k) = xa(j)
275 : end do
276 0 : s% q(k) = q(k)
277 0 : s% dq(k) = dq(k)
278 : end do
279 :
280 0 : call dealloc
281 :
282 : contains
283 :
284 0 : subroutine dealloc
285 0 : deallocate(xh, q, dq)
286 0 : end subroutine dealloc
287 :
288 : end subroutine build_pre_ms_model
289 :
290 :
291 0 : real(dp) function pre_ms_f(lnd, dfdx, lrpar, rpar, lipar, ipar, ierr)
292 : integer, intent(in) :: lrpar, lipar
293 : real(dp), intent(in) :: lnd
294 : real(dp), intent(out) :: dfdx
295 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
296 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
297 : integer, intent(out) :: ierr
298 :
299 : type (star_info), pointer :: s
300 0 : real(dp) :: rho_c, T_c, eps_grav, x, z, abar, zbar, d_log10_P
301 0 : real(dp), pointer :: xa(:)
302 : integer :: i, nz, species
303 0 : real(dp) :: mstar, mstar1
304 :
305 : logical, parameter :: dbg = .false.
306 :
307 : include 'formats'
308 :
309 0 : ierr = 0
310 0 : pre_ms_f = 0
311 0 : if (lipar <= 0) then
312 0 : write(*,*) 'lipar', lipar
313 0 : write(*,*) 'pre_ms f'
314 0 : ierr = -1
315 0 : return
316 : end if
317 :
318 0 : call get_star_ptr(ipar(1), s, ierr)
319 0 : if (ierr /= 0) return
320 :
321 0 : species = s% species
322 :
323 0 : if (associated(s% xh)) deallocate(s% xh)
324 0 : if (associated(s% q)) deallocate(s% q)
325 0 : if (associated(s% dq)) deallocate(s% dq)
326 :
327 0 : rho_c = exp(lnd)
328 :
329 : i = 1 ! rpar(1) for mstar result
330 0 : T_c = rpar(i+1); i = i+1
331 0 : eps_grav = rpar(i+1); i = i+1
332 0 : x = rpar(i+1); i = i+1
333 0 : z = rpar(i+1); i = i+1
334 0 : abar = rpar(i+1); i = i+1
335 0 : zbar = rpar(i+1); i = i+1
336 0 : d_log10_P = rpar(i+1); i = i+1
337 0 : xa => rpar(i+1:i+species); i = i+species
338 0 : if (i > lrpar) then
339 0 : write(*,*) 'i > lrpar', i, lrpar
340 0 : write(*,*) 'pre_ms f'
341 0 : ierr = -1
342 0 : return
343 : end if
344 :
345 0 : mstar = s% mstar ! desired value
346 : mstar1 = mstar ! to keep gfortran quiet
347 :
348 : call build1_pre_ms_model( &
349 : s, T_c, rho_c, d_log10_P, eps_grav, &
350 : x, s% initial_z, abar, zbar, &
351 0 : xa, nz, mstar1, ierr)
352 0 : if (ierr /= 0) then
353 0 : write(*,*) 'failed in build1_pre_ms_model'
354 0 : return
355 : end if
356 :
357 0 : s% nz = nz
358 :
359 0 : rpar(1) = mstar1 ! return the actual mass
360 :
361 0 : pre_ms_f = (mstar - mstar1) / mstar
362 0 : dfdx = 0
363 :
364 : if (dbg) then
365 : write(*,1) 'rho_c', rho_c
366 : write(*,1) 'pre_ms_f', pre_ms_f
367 : write(*,'(A)')
368 : end if
369 :
370 0 : end function pre_ms_f
371 :
372 :
373 0 : subroutine build1_pre_ms_model( &
374 : s, T_c, rho_c, d_log10_P_in, eps_grav_in, &
375 0 : x, z, abar, zbar, xa, nz, mstar, ierr)
376 : use chem_def
377 : use eos_def
378 : use kap_lib
379 : use chem_lib
380 : use eos_lib, only: Radiation_Pressure
381 : use eos_support, only: get_eos, solve_eos_given_PgasT_auto
382 : use star_utils, only: normalize_dqs, set_qs, &
383 : store_r_in_xh, store_lnT_in_xh, store_lnd_in_xh
384 : type (star_info), pointer :: s
385 : real(dp), intent(in) :: &
386 : T_c, rho_c, d_log10_P_in, eps_grav_in, &
387 : x, z, abar, zbar, xa(:)
388 : integer, intent(out) :: nz
389 : real(dp), intent(out) :: mstar ! the mass of the constructed model
390 : integer, intent(out) :: ierr
391 :
392 : real(dp), parameter :: LOGRHO_TOL = 1E-6_dp
393 : real(dp), parameter :: LOGPGAS_TOL = 1E-6_dp
394 :
395 : integer :: i, ii, k, j, prune, max_retries
396 : real(dp), parameter :: &
397 : delta_logPgas = 0.004d0, q_at_nz = 1d-5
398 : real(dp) :: &
399 0 : P_surf_limit, y, dlogPgas, logPgas, Prad, Pgas, try_dlogPgas, logPgas0, &
400 0 : res(num_eos_basic_results), eps_grav, P_c, logP, m, &
401 0 : d_eos_dlnd(num_eos_basic_results), d_eos_dlnT(num_eos_basic_results), &
402 0 : d_eos_dxa(num_eos_basic_results,s% species), &
403 0 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
404 0 : eta, d_eta_dlnRho, d_eta_dlnT, &
405 0 : cgrav, r, rmid, rho, logRho, T, lnT, L, P, P0, dm, m0, L0, r0, lnT0, T0, &
406 0 : rho0, rho_mid, Pmid, chiRho0, chiRho_mid, chiT0, chiT_mid, Cp0, Cp_mid, &
407 0 : grada0, grada_mid, mmid, Tmid, Lmid, &
408 0 : chiRho, chiT, Cp, grada, gradT, logT_surf_limit, logP_surf_limit
409 0 : real(dp), pointer :: xh(:,:), q(:), dq(:) ! model structure info
410 :
411 : logical, parameter :: dbg = .false.
412 :
413 : include 'formats'
414 :
415 0 : ierr = 0
416 :
417 0 : logP_surf_limit = s% pre_ms_logP_surf_limit
418 0 : if (logP_surf_limit <= 0) logP_surf_limit = 3.5d0
419 0 : P_surf_limit = exp10(logP_surf_limit)
420 :
421 0 : logT_surf_limit = s% pre_ms_logT_surf_limit
422 0 : if (logT_surf_limit <= 0) logT_surf_limit = 3.7d0
423 :
424 : if (dbg) write(*,1) 'logT_surf_limit', logT_surf_limit
425 :
426 0 : cgrav = standard_cgrav
427 :
428 0 : eps_grav = eps_grav_in
429 : if (dbg) write(*,1) 'eps_grav', eps_grav
430 :
431 0 : if (d_log10_P_in == 0) then
432 : dlogPgas = delta_logPgas
433 : else
434 0 : dlogPgas = abs(d_log10_P_in)
435 : end if
436 : if (dbg) write(*,1) 'dlogPgas', dlogPgas
437 :
438 : call get_eos( &
439 : s, 0, xa, &
440 : rho_c, log10(rho_c), T_c, log10(T_c), &
441 : res, d_eos_dlnd, d_eos_dlnT, &
442 0 : d_eos_dxa, ierr)
443 0 : if (ierr /= 0) return
444 0 : call unpack_eos_results
445 :
446 0 : logPgas = res(i_lnPgas)/ln10
447 0 : Pgas = exp10(logPgas)
448 0 : P_c = Pgas + Radiation_Pressure(T_c) ! center pressure
449 :
450 0 : mstar = s% mstar ! desired total mass
451 0 : m = q_at_nz*mstar ! mass at nz
452 : ! pressure at innermost point using K&W 10.6
453 0 : P = P_c - 3*cgrav/(8*pi)*pow(pi4*rho_c/3,4d0/3d0)*pow(m,two_thirds)
454 0 : logP = log10(P)
455 :
456 : ! estimate nz from lgP
457 0 : nz = 1 + (logP - s% pre_ms_logP_surf_limit)/dlogPgas
458 :
459 : ! temperature at nz using K&W 10.9 assuming convective core
460 : lnT = log(T_c) - &
461 0 : pow(pi/6,one_third)*cgrav*grada*pow(rho_c*rho_c*m,two_thirds)/P_c
462 0 : T = exp(lnT)
463 :
464 : ! density at nz
465 : call solve_eos_given_PgasT_auto( &
466 : s, 0, xa, &
467 : lnT/ln10, log10(Pgas), LOGRHO_TOL, LOGPGAS_TOL, &
468 : logRho, res, d_eos_dlnd, d_eos_dlnT, d_eos_dxa, &
469 0 : ierr)
470 0 : if (ierr /= 0) return
471 0 : rho = exp10(logRho)
472 0 : call unpack_eos_results
473 :
474 0 : r = pow(m/(pi4*rho/3),one_third) ! radius at nz
475 :
476 0 : y = 1 - (x+z)
477 :
478 0 : do
479 :
480 0 : L = eps_grav*m ! L at nz
481 :
482 : ! check for convective core
483 : call eval_gradT( &
484 : s, zbar, x, y, xa, rho, m, mstar, r, T, lnT, L, P, &
485 : chiRho, chiT, Cp, grada, &
486 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
487 : eta, d_eta_dlnRho, d_eta_dlnT, &
488 0 : gradT, ierr )
489 0 : if (ierr /= 0) return
490 :
491 0 : if (gradT >= grada) exit
492 :
493 0 : eps_grav = 1.1d0*eps_grav
494 :
495 : end do
496 :
497 0 : allocate(xh(s% nvar_hydro,nz), q(nz), dq(nz), stat=ierr)
498 0 : if (ierr /= 0) return
499 0 : s% xh => xh
500 0 : s% dq => dq
501 0 : s% q => q
502 :
503 0 : call store_lnd_in_xh(s, nz, logRho*ln10, xh)
504 0 : call store_lnT_in_xh(s, nz, lnT, xh)
505 0 : call store_r_in_xh(s, nz, r, xh)
506 0 : if (s% i_lum /= 0) xh(s% i_lum,nz) = L
507 :
508 0 : q(nz) = q_at_nz
509 0 : dq(nz) = q_at_nz
510 :
511 : if (dbg) write(*,*) 'nz', nz
512 :
513 0 : max_retries = 10
514 0 : prune = 0
515 0 : step_loop: do k = nz-1, 1, -1
516 :
517 0 : try_dlogPgas = dlogPgas
518 0 : logPgas0 = logPgas
519 0 : P0 = P
520 0 : m0 = m
521 0 : L0 = L
522 0 : r0 = r
523 0 : lnT0 = lnT
524 0 : T0 = T
525 0 : rho0 = rho
526 0 : chiRho0 = chiRho
527 0 : chiT0 = chiT
528 0 : Cp0 = Cp
529 0 : grada0 = grada
530 0 : dm = 0 ! for gfortran
531 :
532 : if (dbg) write(*,3) 'step', k, nz, logPgas0
533 :
534 0 : retry_loop: do j = 1, max_retries
535 :
536 0 : logPgas = logPgas0 - try_dlogPgas
537 0 : Pgas = exp10(logPgas)
538 :
539 : if (j > 1) write(*,2) 'retry', j, logPgas
540 :
541 0 : do i = 1, 2
542 :
543 0 : Prad = Radiation_Pressure(T)
544 0 : P = Pgas + Prad
545 :
546 0 : rho_mid = (rho+rho0)/2
547 :
548 0 : do ii = 1, 10 ! repeat to get hydrostatic balance
549 0 : rmid = pow((r*r*r + r0*r0*r0)/2,one_third)
550 0 : mmid = (m + m0)/2
551 0 : if (ii == 10) exit
552 0 : dm = -pi4*pow4(rmid)*(P-P0)/(cgrav*mmid)
553 0 : m = m0 + dm ! mass at point k
554 0 : r = pow(r0*r0*r0 + dm/(four_thirds_pi*rho_mid),one_third)
555 0 : if (dbg) write(*,2) 'r', ii, r, m, dm
556 : end do
557 :
558 0 : L = L0 + dm*eps_grav ! luminosity at point k
559 0 : Lmid = (L0+L)/2
560 :
561 0 : Pmid = (P+P0)/2
562 :
563 0 : chiRho_mid = (chiRho0 + chiRho)/2
564 0 : chiT_mid = (chiT0 + chiT)/2
565 0 : Cp_mid = (Cp0 + Cp)/2
566 0 : grada_mid = (grada0 + grada)/2
567 :
568 0 : do ii = 1, 2
569 0 : Tmid = (T+T0)/2
570 : call eval_gradT( &
571 : s, zbar, x, y, xa, rho_mid, mmid, mstar, rmid, Tmid, log(Tmid), Lmid, Pmid, &
572 : chiRho_mid, chiT_mid, Cp_mid, grada_mid, &
573 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
574 : eta, d_eta_dlnRho, d_eta_dlnT, &
575 0 : gradT, ierr )
576 0 : if (ierr /= 0) return
577 0 : T = T0 + Tmid*gradT*(P-P0)/Pmid
578 0 : lnT = log(T)
579 0 : if (dbg) write(*,2) 'T', ii, T
580 : end do
581 :
582 0 : if (i == 2) exit
583 :
584 : call solve_eos_given_PgasT_auto( &
585 : s, 0, xa, &
586 : lnT/ln10, logPgas, LOGRHO_TOL, LOGPGAS_TOL, &
587 : logRho, res, d_eos_dlnd, d_eos_dlnT, d_eos_dxa, &
588 0 : ierr)
589 0 : rho = exp10(logRho)
590 0 : if (ierr /= 0) return
591 0 : call unpack_eos_results
592 :
593 : end do
594 :
595 0 : if (lnT <= logT_surf_limit*ln10) then
596 : if (dbg) write(*,*) 'have reached lgT_surf_limit', lnT/ln10, logT_surf_limit
597 : prune = k
598 : exit step_loop
599 : end if
600 :
601 0 : if (P <= P_surf_limit) then
602 : if (dbg) write(*,1) 'have reached P_surf limit', P, P_surf_limit
603 : prune = k
604 : exit step_loop
605 : end if
606 :
607 0 : call store_lnd_in_xh(s, k, logRho*ln10, xh)
608 0 : call store_lnT_in_xh(s, k, lnT, xh)
609 0 : call store_r_in_xh(s, k, r, xh)
610 0 : if (s% i_lum /= 0) xh(s% i_lum,k) = L
611 0 : q(k) = m/mstar
612 0 : dq(k) = dm/mstar
613 :
614 : if (dbg) then
615 : write(*,2) 'L', k, L
616 : write(*,2) 'q(k)', k, q(k)
617 : write(*,2) 'dq(k)', k, dq(k)
618 : end if
619 :
620 : exit retry_loop
621 :
622 : end do retry_loop
623 :
624 : end do step_loop
625 :
626 0 : if (prune > 0) then ! move stuff and reduce nz
627 : if (dbg) write(*,*) 'prune', prune
628 0 : do k=1,nz-prune
629 0 : xh(:,k) = xh(:,k+prune)
630 0 : q(k) = q(k+prune)
631 0 : dq(k) = dq(k+prune)
632 : end do
633 0 : m = mstar*q(1)
634 0 : nz = nz-prune
635 : if (dbg) write(*,*) 'final nz', nz
636 : end if
637 :
638 0 : mstar = m ! actual total mass
639 :
640 0 : if (.not. s% do_normalize_dqs_as_part_of_set_qs) then
641 0 : call normalize_dqs(s, nz, dq, ierr)
642 0 : if (ierr /= 0) then
643 0 : if (s% report_ierr) write(*,*) 'normalize_dqs failed in pre ms model'
644 0 : return
645 : end if
646 : end if
647 0 : call set_qs(s, nz, q, dq, ierr)
648 0 : if (ierr /= 0) then
649 0 : if (s% report_ierr) write(*,*) 'set_qs failed in pre ms model'
650 0 : return
651 : end if
652 :
653 : contains
654 :
655 :
656 0 : subroutine unpack_eos_results
657 0 : chiRho = res(i_chiRho)
658 0 : chiT = res(i_chiT)
659 0 : Cp = res(i_cp)
660 0 : grada = res(i_grad_ad)
661 0 : lnfree_e = res(i_lnfree_e)
662 0 : d_lnfree_e_dlnRho = d_eos_dlnd(i_lnfree_e)
663 0 : d_lnfree_e_dlnT = d_eos_dlnT(i_lnfree_e)
664 0 : eta = res(i_eta)
665 0 : d_eta_dlnRho = d_eos_dlnd(i_eta)
666 0 : d_eta_dlnT = d_eos_dlnT(i_eta)
667 0 : end subroutine unpack_eos_results
668 :
669 :
670 : end subroutine build1_pre_ms_model
671 :
672 :
673 0 : subroutine eval_gradT( &
674 0 : s, zbar, x, y, xa, rho, m, mstar, r, T, lnT, L, P, &
675 : chiRho, chiT, Cp, grada, &
676 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
677 : eta, d_eta_dlnRho, d_eta_dlnT, &
678 : gradT, ierr )
679 : use chem_def, only: ih1
680 : use turb_support, only: get_gradT
681 : use kap_def, only : num_kap_fracs
682 : use kap_lib, only : kap_get
683 :
684 : type (star_info), pointer :: s
685 : real(dp), intent(in) :: &
686 : zbar, x, y, xa(:), rho, m, mstar, r, T, lnT, L, P, &
687 : chiRho, chiT, Cp, grada, &
688 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
689 : eta, d_eta_dlnRho, d_eta_dlnT
690 : real(dp), intent(out) :: gradT
691 : integer, intent(out) :: ierr
692 :
693 :
694 0 : real(dp) :: dlnkap_dlnd, dlnkap_dlnT, gradL_composition_term, &
695 0 : opacity, grav, scale_height, scale_height2, gradr, cgrav
696 0 : real(dp) :: kap_fracs(num_kap_fracs), dlnkap_dxa(s% species)
697 0 : real(dp) :: Y_face, conv_vel, D, Gamma ! Not used
698 : integer :: mixing_type
699 :
700 0 : ierr = 0
701 :
702 0 : if (s% use_simple_es_for_kap) then
703 0 : opacity = 0.2d0*(1 + x)
704 : dlnkap_dlnd = 0
705 : dlnkap_dlnT = 0
706 : else
707 : call kap_get( &
708 : s% kap_handle, s% species, s% chem_id, s% net_iso, xa, &
709 : log10(Rho), lnT/ln10, &
710 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
711 : eta, d_eta_dlnRho, d_eta_dlnT, &
712 0 : kap_fracs, opacity, dlnkap_dlnd, dlnkap_dlnT, dlnkap_dxa, ierr)
713 0 : if (ierr /= 0) then
714 0 : write(*,*) 'failed in kap_get in eval_gradT'
715 : return
716 : end if
717 : end if
718 :
719 0 : gradL_composition_term = 0d0
720 0 : cgrav = standard_cgrav
721 0 : grav = cgrav*m/pow2(r)
722 0 : scale_height = P/(grav*rho) ! this assumes HSE
723 0 : if (s% alt_scale_height_flag) then
724 0 : scale_height2 = sqrt(P/cgrav)/rho
725 0 : if (scale_height2 < scale_height) then
726 0 : scale_height = scale_height2
727 : end if
728 : end if
729 0 : gradr = P*opacity*L/(16d0*pi*clight*m*cgrav*crad*pow4(T)/3d0)
730 :
731 : call get_gradT(s, s% MLT_option, & ! used to create models
732 : r, L, T, P, opacity, rho, chiRho, chiT, Cp, gradr, grada, scale_height, &
733 : s% net_iso(ih1), x, standard_cgrav, m, gradL_composition_term, s% mixing_length_alpha, &
734 0 : mixing_type, gradT, Y_face, conv_vel, D, Gamma, ierr)
735 :
736 0 : end subroutine eval_gradT
737 :
738 : end module pre_ms_model
|
