Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2025 Ebraheem Farag & 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 tdc_hydro_support
21 :
22 : use star_private_def
23 : use const_def, only: dp, ln10, pi, lsun, rsun, msun
24 : use utils_lib, only: is_bad
25 : use auto_diff
26 : use auto_diff_support
27 : use accurate_sum_auto_diff_star_order1
28 : use star_utils
29 :
30 : implicit none
31 :
32 : private
33 : public :: remesh_for_TDC_pulsations
34 :
35 : contains
36 :
37 0 : subroutine remesh_for_TDC_pulsations(s, ierr)
38 : ! uses these controls
39 : ! TDC_hydro_nz = 150
40 : ! TDC_hydro_nz_outer = 40
41 : ! TDC_hydro_T_anchor = 11d3
42 : ! TDC_hydro_dq_1_factor = 2d0
43 : use interp_1d_def, only: pm_work_size
44 : use interp_1d_lib, only: interpolate_vector_pm
45 : type(star_info), pointer :: s
46 : integer, intent(out) :: ierr
47 : integer :: k, j, nz_old, nz
48 0 : real(dp) :: xm_anchor, P_surf, T_surf, old_L1, old_r1
49 : real(dp), allocatable, dimension(:) :: &
50 0 : xm_old, xm, xm_mid_old, xm_mid, v_old, v_new
51 0 : real(dp), pointer :: work1(:) ! =(nz_old+1, pm_work_size)
52 : include 'formats'
53 0 : ierr = 0
54 0 : nz_old = s%nz
55 0 : nz = s%TDC_hydro_nz
56 0 : if (nz == nz_old) return ! assume have already done remesh for RSP2
57 0 : if (nz > nz_old) call mesa_error(__FILE__, __LINE__, 'remesh_for_RSP2 cannot increase nz')
58 0 : call setvars2(ierr)
59 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'remesh_for_RSP2 failed in setvars')
60 0 : old_L1 = s%L(1)
61 0 : old_r1 = s%r(1)
62 0 : call set_phot_info(s) ! sets Teff
63 0 : call get_PT_surf2(P_surf, T_surf, ierr)
64 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'remesh_for_RSP2 failed in get_PT_surf')
65 : allocate ( &
66 0 : xm_old(nz_old + 1), xm_mid_old(nz_old), v_old(nz_old + 1), &
67 0 : xm(nz + 1), xm_mid(nz), v_new(nz + 1), work1((nz_old + 1)*pm_work_size))
68 0 : call set_xm_old2
69 0 : call find_xm_anchor2
70 0 : call set_xm_new2
71 0 : call interpolate1_face_val2(s%i_lnR, log(max(1d0, s%r_center)))
72 0 : call check_new_lnR2
73 0 : call interpolate1_face_val2(s%i_lum, s%L_center)
74 0 : if (s%i_v /= 0) call interpolate1_face_val2(s%i_v, s%v_center)
75 0 : call set_new_lnd2
76 0 : call interpolate1_cell_val2(s%i_lnT)
77 0 : do j = 1, s%species
78 0 : call interpolate1_xa2(j)
79 : end do
80 0 : call rescale_xa2
81 0 : call revise_lnT_for_QHSE2(P_surf, ierr)
82 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'remesh_for_RSP2 failed in revise_lnT_for_QHSE')
83 0 : do k = 1, nz
84 0 : call set_Hp_face2(k)
85 : end do
86 0 : deallocate (work1)
87 0 : s%nz = nz
88 0 : write (*, 1) 'new old L_surf/Lsun', s%xh(s%i_lum, 1)/Lsun, old_L1/Lsun
89 0 : write (*, 1) 'new old R_surf/Rsun', exp(s%xh(s%i_lnR, 1))/Rsun, old_r1/Rsun
90 0 : write (*, '(A)')
91 : !call mesa_error(__FILE__,__LINE__,'remesh_for_RSP2')
92 :
93 : contains
94 :
95 0 : subroutine setvars2(ierr)
96 0 : use hydro_vars, only: unpack_xh, set_hydro_vars
97 : integer, intent(out) :: ierr
98 : logical, parameter :: &
99 : skip_basic_vars = .false., &
100 : skip_micro_vars = .false., &
101 : skip_m_grav_and_grav = .false., &
102 : skip_net = .true., &
103 : skip_neu = .true., &
104 : skip_kap = .false., &
105 : skip_grads = .true., &
106 : skip_rotation = .true., &
107 : skip_brunt = .true., &
108 : skip_other_cgrav = .true., &
109 : skip_mixing_info = .true., &
110 : skip_set_cz_bdy_mass = .true., &
111 : skip_mlt = .true., &
112 : skip_eos = .false.
113 : ierr = 0
114 0 : call unpack_xh(s, ierr)
115 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'remesh_for_RSP2 failed in unpack_xh')
116 : call set_hydro_vars( &
117 : s, 1, nz_old, skip_basic_vars, &
118 : skip_micro_vars, skip_m_grav_and_grav, skip_eos, skip_net, skip_neu, &
119 : skip_kap, skip_grads, skip_rotation, skip_brunt, skip_other_cgrav, &
120 0 : skip_mixing_info, skip_set_cz_bdy_mass, skip_mlt, ierr)
121 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'remesh_for_RSP2 failed in set_hydro_vars')
122 0 : end subroutine setvars2
123 :
124 0 : subroutine get_PT_surf2(P_surf, T_surf, ierr)
125 0 : use atm_support, only: get_atm_PT
126 : real(dp), intent(out) :: P_surf, T_surf
127 : integer, intent(out) :: ierr
128 : real(dp) :: &
129 0 : Teff, lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
130 0 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap
131 : logical, parameter :: skip_partials = .true.
132 : include 'formats'
133 0 : ierr = 0
134 0 : call set_phot_info(s) ! sets s% Teff
135 0 : Teff = s%Teff
136 : call get_atm_PT( & ! this uses s% opacity(1)
137 : s, s%tau_factor*s%tau_base, s%L(1), s%r(1), s%m(1), s%cgrav(1), skip_partials, &
138 : Teff, lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
139 0 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, ierr)
140 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'get_P_surf failed in get_atm_PT')
141 0 : P_surf = exp(lnP_surf)
142 0 : T_surf = exp(lnT_surf)
143 : return
144 :
145 : write (*, 1) 'get_PT_surf P_surf', P_surf
146 : write (*, 1) 'get_PT_surf T_surf', T_surf
147 : write (*, 1) 'get_PT_surf Teff', Teff
148 : write (*, 1) 'get_PT_surf opacity(1)', s%opacity(1)
149 : write (*, 1)
150 : !call mesa_error(__FILE__,__LINE__,'get_PT_surf')
151 0 : end subroutine get_PT_surf2
152 :
153 0 : subroutine set_xm_old2
154 0 : xm_old(1) = 0d0
155 0 : do k = 2, nz_old
156 0 : xm_old(k) = xm_old(k - 1) + s%dm(k - 1)
157 : end do
158 0 : xm_old(nz_old + 1) = s%xmstar
159 0 : do k = 1, nz_old
160 0 : xm_mid_old(k) = xm_old(k) + 0.5d0*s%dm(k)
161 : end do
162 0 : end subroutine set_xm_old2
163 :
164 0 : subroutine find_xm_anchor2
165 0 : real(dp) :: lnT_anchor, xmm1, xm00, lnTm1, lnT00
166 : include 'formats'
167 0 : lnT_anchor = log(s%TDC_hydro_T_anchor)
168 0 : if (lnT_anchor <= s%xh(s%i_lnT, 1)) then
169 0 : write (*, 1) 'T_anchor < T_surf', s%TDC_hydro_T_anchor, exp(s%xh(s%i_lnT, 1))
170 0 : call mesa_error(__FILE__, __LINE__, 'find_xm_anchor')
171 : end if
172 0 : xm_anchor = xm_old(nz_old)
173 0 : do k = 2, nz_old
174 0 : if (s%xh(s%i_lnT, k) >= lnT_anchor) then
175 0 : xmm1 = xm_old(k - 1)
176 0 : xm00 = xm_old(k)
177 0 : lnTm1 = s%xh(s%i_lnT, k - 1)
178 0 : lnT00 = s%xh(s%i_lnT, k)
179 : xm_anchor = xmm1 + &
180 0 : (xm00 - xmm1)*(lnT_anchor - lnTm1)/(lnT00 - lnTm1)
181 0 : if (is_bad(xm_anchor) .or. xm_anchor <= 0d0) then
182 0 : write (*, 2) 'bad xm_anchor', k, xm_anchor, xmm1, xm00, lnTm1, lnT00, lnT_anchor, s%lnT(1)
183 0 : call mesa_error(__FILE__, __LINE__, 'find_xm_anchor')
184 : end if
185 0 : return
186 : end if
187 : end do
188 : end subroutine find_xm_anchor2
189 :
190 0 : subroutine set_xm_new2 ! sets xm, dm, m, dq, q
191 : integer :: nz_outer, k, n_inner
192 0 : real(dp) :: dq_1_factor, dxm_outer, lnx, dlnx, base_dm, rem_mass, H
193 0 : real(dp) :: H_low, H_high, H_mid, f_low, f_high, f_mid
194 : integer :: iter
195 : include 'formats'
196 0 : nz_outer = s%TDC_hydro_nz_outer
197 0 : dq_1_factor = s%TDC_hydro_dq_1_factor
198 0 : dxm_outer = xm_anchor/(nz_outer - 1d0 + dq_1_factor)
199 0 : xm(1) = 0d0
200 0 : xm(2) = dxm_outer*dq_1_factor
201 0 : s%dm(1) = xm(2)
202 0 : do k = 3, nz_outer + 1
203 0 : xm(k) = xm(k - 1) + dxm_outer
204 0 : s%dm(k - 1) = dxm_outer
205 : end do
206 :
207 0 : if (.not. s%remesh_for_TDC_pulsations_log_core_zoning) then
208 : ! do rsp style core zoning with a power law on dq
209 :
210 : ! solve for a smooth ramp factor H via bisection
211 0 : n_inner = nz - nz_outer
212 0 : rem_mass = s%xmstar - xm(nz_outer + 1)
213 0 : base_dm = dxm_outer !first dm equals outer spacing
214 :
215 : ! define function f(H) = base_dm*(sum_{j=1..n_inner-1}H^j) - rem_mass
216 :
217 0 : H_low = 1.001! Heuristics
218 0 : H_high = 1.40! Heuristics
219 : ! compute f at bounds
220 0 : f_low = base_dm*((H_low*(1d0 - H_low**(n_inner - 1))/(1d0 - H_low))) - rem_mass
221 0 : f_high = base_dm*((H_high*(1d0 - H_high**(n_inner - 1))/(1d0 - H_high))) - rem_mass
222 0 : do iter = 1, 1000
223 0 : H_mid = 0.5d0*(H_low + H_high)
224 0 : f_mid = base_dm*((H_mid*(1d0 - H_mid**(n_inner - 1))/(1d0 - H_mid))) - rem_mass
225 0 : if (abs(f_mid) < 1d-12*rem_mass) exit
226 0 : if (f_low*f_mid <= 0d0) then
227 : H_high = H_mid
228 0 : f_high = f_mid
229 : else
230 0 : H_low = H_mid
231 0 : f_low = f_mid
232 : end if
233 : end do
234 0 : H = H_mid
235 :
236 : ! first interior cell:
237 0 : s%dm(nz_outer + 1) = base_dm
238 0 : xm(nz_outer + 2) = xm(nz_outer + 1) + s%dm(nz_outer + 1)
239 :
240 : ! subsequent interior cells: ramp by H per zone (except final)
241 0 : do k = nz_outer + 2, nz - 1
242 0 : s%dm(k) = H**(k - nz_outer - 1)*base_dm
243 0 : xm(k + 1) = xm(k) + s%dm(k)
244 : end do
245 :
246 : ! final interior cell absorbs any remaining mass
247 0 : s%dm(nz) = s%xmstar - xm(nz)
248 0 : xm(nz + 1) = s%xmstar
249 :
250 : else ! use log zoning inward from anchor to core.
251 0 : lnx = log(xm(nz_outer + 1))
252 0 : if (is_bad(lnx)) then
253 0 : write (*, 2) 'bad lnx', nz_outer + 1, lnx, xm(nz_outer + 1)
254 0 : call mesa_error(__FILE__, __LINE__, 'set_xm_new')
255 : end if
256 0 : dlnx = (log(s%xmstar) - lnx)/(nz - nz_outer)
257 0 : do k = nz_outer + 2, nz
258 0 : lnx = lnx + dlnx
259 0 : xm(k) = exp(lnx)
260 0 : s%dm(k - 1) = xm(k) - xm(k - 1)
261 : end do
262 0 : s%dm(nz) = s%xmstar - xm(nz)
263 :
264 : ! — enforce the last boundary at total mass
265 0 : xm(nz + 1) = s%xmstar
266 :
267 : ! — recompute cell masses
268 0 : do k = nz_outer + 1, nz
269 0 : s%dm(k) = xm(k + 1) - xm(k)
270 : end do
271 :
272 : end if
273 :
274 0 : do k = 1, nz - 1
275 0 : xm_mid(k) = 0.5d0*(xm(k) + xm(k + 1))
276 : end do
277 0 : xm_mid(nz) = 0.5d0*(xm(nz) + s%xmstar)
278 0 : s%m(1) = s%mstar
279 0 : s%q(1) = 1d0
280 0 : s%dq(1) = s%dm(1)/s%xmstar
281 0 : do k = 2, nz
282 0 : s%m(k) = s%m(k - 1) - s%dm(k - 1)
283 0 : s%dq(k) = s%dm(k)/s%xmstar
284 0 : s%q(k) = s%q(k - 1) - s%dq(k - 1)
285 : end do
286 0 : call set_dm_bar(s, s%nz, s%dm, s%dm_bar)
287 0 : return
288 :
289 : do k = 2, nz
290 : write (*, 2) 'dm(k)/dm(k-1) m(k)', k, s%dm(k)/s%dm(k - 1), s%m(k)/Msun
291 : end do
292 : write (*, 1) 'm_center', s%m_center/msun
293 : call mesa_error(__FILE__, __LINE__, 'set_xm_new')
294 : end subroutine set_xm_new2
295 :
296 0 : subroutine interpolate1_face_val2(i, cntr_val)
297 : integer, intent(in) :: i
298 : real(dp), intent(in) :: cntr_val
299 0 : do k = 1, nz_old
300 0 : v_old(k) = s%xh(i, k)
301 : end do
302 0 : v_old(nz_old + 1) = cntr_val
303 : call interpolate_vector_pm( &
304 0 : nz_old + 1, xm_old, nz + 1, xm, v_old, v_new, work1, 'remesh_for_RSP2', ierr)
305 0 : do k = 1, nz
306 0 : s%xh(i, k) = v_new(k)
307 : end do
308 0 : end subroutine interpolate1_face_val2
309 :
310 0 : subroutine check_new_lnR2
311 : include 'formats'
312 0 : do k = 1, nz
313 0 : s%lnR(k) = s%xh(s%i_lnR, k)
314 0 : s%r(k) = exp(s%lnR(k))
315 : end do
316 0 : do k = 1, nz - 1
317 0 : if (s%r(k) <= s%r(k + 1)) then
318 0 : write (*, 2) 'bad r', k, s%r(k), s%r(k + 1)
319 0 : call mesa_error(__FILE__, __LINE__, 'check_new_lnR remesh rsp2')
320 : end if
321 : end do
322 0 : if (s%r(nz) <= s%r_center) then
323 0 : write (*, 2) 'bad r center', nz, s%r(nz), s%r_center
324 0 : call mesa_error(__FILE__, __LINE__, 'check_new_lnR remesh rsp2')
325 : end if
326 0 : end subroutine check_new_lnR2
327 :
328 0 : subroutine set_new_lnd2
329 0 : real(dp) :: vol, r300, r3p1
330 : include 'formats'
331 0 : do k = 1, nz
332 0 : r300 = pow3(s%r(k))
333 0 : if (k < nz) then
334 0 : r3p1 = pow3(s%r(k + 1))
335 : else
336 0 : r3p1 = pow3(s%r_center)
337 : end if
338 0 : vol = (4d0*pi/3d0)*(r300 - r3p1)
339 0 : s%rho(k) = s%dm(k)/vol
340 0 : s%lnd(k) = log(s%rho(k))
341 0 : s%xh(s%i_lnd, k) = s%lnd(k)
342 0 : if (is_bad(s%lnd(k))) then
343 0 : write (*, 2) 'bad lnd vol dm r300 r3p1', k, s%lnd(k), vol, s%dm(k), r300, r3p1
344 0 : call mesa_error(__FILE__, __LINE__, 'remesh for rsp2')
345 : end if
346 : end do
347 0 : end subroutine set_new_lnd2
348 :
349 0 : subroutine interpolate1_cell_val2(i)
350 : integer, intent(in) :: i
351 0 : do k = 1, nz_old
352 0 : v_old(k) = s%xh(i, k)
353 : end do
354 : call interpolate_vector_pm( &
355 0 : nz_old, xm_mid_old, nz, xm_mid, v_old, v_new, work1, 'remesh_for_RSP2', ierr)
356 0 : do k = 1, nz
357 0 : s%xh(i, k) = v_new(k)
358 : end do
359 0 : end subroutine interpolate1_cell_val2
360 :
361 0 : subroutine interpolate1_xa2(j)
362 : integer, intent(in) :: j
363 0 : do k = 1, nz_old
364 0 : v_old(k) = s%xa(j, k)
365 : end do
366 : call interpolate_vector_pm( &
367 0 : nz_old, xm_mid_old, nz, xm_mid, v_old, v_new, work1, 'remesh_for_RSP2', ierr)
368 0 : do k = 1, nz
369 0 : s%xa(j, k) = v_new(k)
370 : end do
371 0 : end subroutine interpolate1_xa2
372 :
373 0 : subroutine rescale_xa2
374 : integer :: k, j
375 0 : real(dp) :: sum_xa
376 0 : do k = 1, nz
377 0 : sum_xa = sum(s%xa(1:s%species, k))
378 0 : do j = 1, s%species
379 0 : s%xa(j, k) = s%xa(j, k)/sum_xa
380 : end do
381 : end do
382 0 : end subroutine rescale_xa2
383 :
384 0 : subroutine revise_lnT_for_QHSE2(P_surf, ierr)
385 : use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results
386 : use chem_def, only: chem_isos
387 : use eos_support, only: solve_eos_given_DP
388 : use eos_def, only: i_eta, i_lnfree_e
389 : use kap_def, only: num_kap_fracs
390 : use kap_support, only: get_kap
391 : real(dp), intent(in) :: P_surf
392 : integer, intent(out) :: ierr
393 0 : real(dp) :: logRho, logP, logT_guess, &
394 0 : logT_tol, logP_tol, logT, P_m1, P_00, dm_face, &
395 0 : kap_fracs(num_kap_fracs), kap, dlnkap_dlnRho, dlnkap_dlnT, &
396 0 : old_kap, new_P_surf, new_T_surf
397 : real(dp), dimension(num_eos_basic_results) :: &
398 0 : res, d_dlnd, d_dlnT
399 0 : real(dp) :: dres_dxa(num_eos_d_dxa_results, s%species)
400 : include 'formats'
401 0 : ierr = 0
402 0 : P_m1 = P_surf
403 0 : do k = 1, nz
404 0 : s%lnT(k) = s%xh(s%i_lnT, k)
405 0 : s%lnR(k) = s%xh(s%i_lnR, k)
406 0 : s%r(k) = exp(s%lnR(k))
407 : end do
408 : !write(*,1) 'before revise_lnT_for_QHSE: logT cntr', s% lnT(nz)/ln10
409 0 : do k = 1, nz
410 0 : if (k < nz) then
411 0 : dm_face = s%dm_bar(k)
412 : else
413 0 : dm_face = 0.5d0*(s%dm(k - 1) + s%dm(k))
414 : end if
415 0 : P_00 = P_m1 + s%cgrav(k)*s%m(k)*dm_face/(4d0*pi*pow4(s%r(k)))
416 0 : logP = log10(P_00) ! value for QHSE
417 0 : s%lnPeos(k) = logP/ln10
418 0 : s%Peos(k) = P_00
419 0 : logRho = s%lnd(k)/ln10
420 0 : logT_guess = s%lnT(k)/ln10
421 0 : logT_tol = 1d-11
422 0 : logP_tol = 1d-11
423 : call solve_eos_given_DP( &
424 : s, k, s%xa(:, k), &
425 : logRho, logP, logT_guess, logT_tol, logP_tol, &
426 0 : logT, res, d_dlnd, d_dlnT, dres_dxa, ierr)
427 0 : if (ierr /= 0) then
428 0 : write (*, 2) 'solve_eos_given_DP failed', k
429 0 : write (*, '(A)')
430 0 : write (*, 1) 'sum(xa)', sum(s%xa(:, k))
431 0 : do j = 1, s%species
432 0 : write (*, 4) 'xa(j,k) '//trim(chem_isos%name(s%chem_id(j))), j, j + s%nvar_hydro, k, s%xa(j, k)
433 : end do
434 0 : write (*, 1) 'logRho', logRho
435 0 : write (*, 1) 'logP', logP
436 0 : write (*, 1) 'logT_guess', logT_guess
437 0 : write (*, 1) 'logT_tol', logT_tol
438 0 : write (*, 1) 'logP_tol', logP_tol
439 0 : write (*, '(A)')
440 0 : call mesa_error(__FILE__, __LINE__, 'revise_lnT_for_QHSE')
441 : end if
442 0 : s%lnT(k) = logT*ln10
443 0 : s%xh(s%i_lnT, k) = s%lnT(k)
444 : !write(*,2) 'logP dlogT logT logT_guess logRho', k, &
445 : ! logP, logT - logT_guess, logT, logT_guess, logRho
446 0 : P_m1 = P_00
447 :
448 0 : if (k == 1) then ! get opacity and recheck surf BCs
449 : call get_kap( & ! assume zbar is set
450 : s, k, s%zbar(k), s%xa(:, k), logRho, logT, &
451 : res(i_lnfree_e), d_dlnd(i_lnfree_e), d_dlnT(i_lnfree_e), &
452 : res(i_eta), d_dlnd(i_eta), d_dlnT(i_eta), &
453 : kap_fracs, kap, dlnkap_dlnRho, dlnkap_dlnT, &
454 0 : ierr)
455 0 : if (ierr /= 0) then
456 0 : write (*, 2) 'get_kap failed', k
457 0 : call mesa_error(__FILE__, __LINE__, 'revise_lnT_for_QHSE')
458 : end if
459 0 : old_kap = s%opacity(1)
460 0 : s%opacity(1) = kap ! for use by atm surf PT
461 0 : call get_PT_surf2(new_P_surf, new_T_surf, ierr)
462 0 : if (ierr /= 0) then
463 0 : write (*, 2) 'get_PT_surf failed', k
464 0 : call mesa_error(__FILE__, __LINE__, 'revise_lnT_for_QHSE')
465 : end if
466 0 : write (*, 1) 'new old T_surf', new_T_surf, T_surf
467 0 : write (*, 1) 'new old P_surf', new_P_surf, P_surf
468 0 : write (*, 1) 'new old kap(1)', kap, old_kap
469 : !call mesa_error(__FILE__,__LINE__,'revise_lnT_for_QHSE')
470 : end if
471 :
472 : end do
473 : !write(*,1) 'after revise_lnT_for_QHSE: logT cntr', s% lnT(nz)/ln10
474 : !stop
475 0 : end subroutine revise_lnT_for_QHSE2
476 :
477 0 : subroutine set_Hp_face2(k)
478 0 : use tdc_hydro, only: get_RSP2_alfa_beta_face_weights
479 : integer, intent(in) :: k
480 0 : real(dp) :: r_00, d_00, Peos_00, Peos_div_rho, Hp_face, &
481 0 : d_m1, Peos_m1, alfa, beta
482 0 : r_00 = s%r(k)
483 0 : d_00 = s%rho(k)
484 0 : Peos_00 = s%Peos(k)
485 0 : if (k == 1) then
486 0 : Peos_div_rho = Peos_00/d_00
487 0 : Hp_face = pow2(r_00)*Peos_div_rho/(s%cgrav(k)*s%m(k))
488 : else
489 0 : d_m1 = s%rho(k - 1)
490 0 : Peos_m1 = s%Peos(k - 1)
491 0 : call get_RSP2_alfa_beta_face_weights(s, k, alfa, beta)
492 0 : Peos_div_rho = alfa*Peos_00/d_00 + beta*Peos_m1/d_m1
493 0 : Hp_face = pow2(r_00)*Peos_div_rho/(s%cgrav(k)*s%m(k))
494 : end if
495 0 : s%Hp_face(k) = get_scale_height_face_val(s, k)!Hp_face
496 : !s% xh(s% i_Hp, k) = Hp_face
497 0 : end subroutine set_Hp_face2
498 :
499 : end subroutine remesh_for_TDC_pulsations
500 :
501 : end module tdc_hydro_support
|
