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