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 hydro_vars
21 :
22 : use star_private_def
23 : use const_def, only: dp, i8, pi, ln10, boltz_sigma, clight, standard_cgrav, two_thirds, four_thirds, lsun, rsun, no_mixing
24 : use chem_def, only: chem_isos
25 : use utils_lib, only: mesa_error, is_bad
26 :
27 : implicit none
28 :
29 : private
30 : public :: set_vars_if_needed
31 : public :: set_vars
32 : public :: set_final_vars
33 : public :: update_vars
34 : public :: set_cgrav
35 : public :: set_hydro_vars
36 : public :: unpack_xh
37 : public :: set_Teff_info_for_eqns
38 : public :: set_Teff
39 : public :: get_surf_PT
40 : public :: set_grads
41 :
42 : logical, parameter :: dbg = .false.
43 : logical, parameter :: trace_setvars = .false.
44 :
45 : contains
46 :
47 119 : subroutine set_vars_if_needed(s, dt, str, ierr)
48 : type (star_info), pointer :: s
49 : real(dp), intent(in) :: dt
50 : character (len=*), intent(in) :: str
51 : integer, intent(out) :: ierr
52 87 : if (.not. s% need_to_setvars) then
53 : if (trace_setvars) write(*,*) '** skip set_vars ' // trim(str)
54 55 : s% num_skipped_setvars = s% num_skipped_setvars + 1
55 55 : return
56 : end if
57 : if (trace_setvars) write(*,*) 'set_vars ' // trim(str)
58 32 : call set_vars(s, dt, ierr)
59 : end subroutine set_vars_if_needed
60 :
61 :
62 33 : subroutine set_vars(s, dt, ierr)
63 : use star_utils, only: reset_starting_vectors
64 : type (star_info), pointer :: s
65 : real(dp), intent(in) :: dt
66 : integer, intent(out) :: ierr
67 : logical, parameter :: &
68 : skip_m_grav_and_grav = .false., &
69 : skip_net = .false., &
70 : skip_neu = .false., &
71 : skip_kap = .false., &
72 : skip_grads = .false., &
73 : skip_rotation = .false., &
74 : skip_brunt = .false., &
75 : skip_irradiation_heat = .false.
76 : logical :: skip_set_cz_bdy_mass, skip_mixing_info
77 33 : ierr = 0
78 33 : s% num_setvars = s% num_setvars + 1
79 : if (dbg) write(*,*) 'set_vars'
80 33 : call reset_starting_vectors(s)
81 33 : skip_set_cz_bdy_mass = .not. s% have_new_generation
82 33 : skip_mixing_info = .not. s% okay_to_set_mixing_info
83 : call set_some_vars( &
84 : s, skip_m_grav_and_grav, &
85 : skip_net, skip_neu, skip_kap, skip_grads, skip_rotation, &
86 : skip_brunt, skip_mixing_info, skip_set_cz_bdy_mass, skip_irradiation_heat, &
87 33 : dt, ierr)
88 33 : end subroutine set_vars
89 :
90 :
91 11 : subroutine set_final_vars(s, dt, ierr)
92 : use rates_def, only: num_rvs
93 : type (star_info), pointer :: s
94 : real(dp), intent(in) :: dt
95 : integer, intent(out) :: ierr
96 :
97 : logical :: &
98 : skip_basic_vars, &
99 : skip_micro_vars, &
100 : skip_m_grav_and_grav, &
101 : skip_eos, &
102 : skip_net, &
103 : skip_neu, &
104 : skip_kap, &
105 : skip_grads, &
106 : skip_rotation, &
107 : skip_brunt, &
108 : skip_other_cgrav, &
109 : skip_mixing_info, &
110 : skip_set_cz_bdy_mass, &
111 : skip_mlt
112 : integer :: nz, k
113 :
114 : include 'formats'
115 :
116 11 : ierr = 0
117 11 : nz = s% nz
118 :
119 11 : skip_grads = .false.
120 11 : skip_rotation = .false.
121 11 : skip_brunt = .false.
122 11 : skip_other_cgrav = .false.
123 11 : skip_set_cz_bdy_mass = .false.
124 11 : skip_m_grav_and_grav = .false.
125 11 : skip_mixing_info = .not. s% recalc_mix_info_after_evolve
126 :
127 : ! only need to do things that were skipped in set_vars_for_solver
128 : ! i.e., skip what it already did for the last solver iteration
129 11 : skip_basic_vars = .not. s% need_to_setvars
130 11 : skip_micro_vars = .not. s% need_to_setvars
131 11 : skip_kap = .not. s% need_to_setvars
132 11 : skip_neu = .not. s% need_to_setvars
133 11 : skip_net = .not. s% need_to_setvars
134 11 : skip_eos = .not. s% need_to_setvars
135 11 : skip_mlt = .not. s% need_to_setvars
136 :
137 11 : if (s% need_to_setvars) then
138 0 : s% num_setvars = s% num_setvars + 1
139 : if (trace_setvars) write(*,*) 'set_vars in set_final_vars'
140 : else
141 11 : s% num_skipped_setvars = s% num_skipped_setvars + 1
142 : if (trace_setvars) write(*,*) '** skip set_vars in set_final_vars'
143 : end if
144 : if (trace_setvars) write(*,*)
145 :
146 : call set_hydro_vars( &
147 : s, 1, nz, skip_basic_vars, &
148 : skip_micro_vars, skip_m_grav_and_grav, skip_eos, skip_net, skip_neu, &
149 : skip_kap, skip_grads, skip_rotation, skip_brunt, skip_other_cgrav, &
150 11 : skip_mixing_info, skip_set_cz_bdy_mass, skip_mlt, ierr)
151 11 : if (ierr /= 0) then
152 0 : if (s% report_ierr .or. dbg) &
153 0 : write(*,*) 'update_vars: set_hydro_vars returned ierr', ierr
154 0 : return
155 : end if
156 :
157 11 : if (s% op_split_burn) then
158 0 : do k = 1, nz
159 0 : if (s% T_start(k) >= s% op_split_burn_min_T) &
160 0 : s% eps_nuc(k) = s% burn_avg_epsnuc(k)
161 : end do
162 : end if
163 :
164 : end subroutine set_final_vars
165 :
166 :
167 33 : subroutine set_some_vars( &
168 : s, skip_m_grav_and_grav, &
169 : skip_net, skip_neu, skip_kap, skip_grads, skip_rotation, &
170 : skip_brunt, skip_mixing_info, skip_set_cz_bdy_mass, skip_irradiation_heat, &
171 : dt, ierr)
172 : use star_utils, only: start_time, update_time
173 : type (star_info), pointer :: s
174 : logical, intent(in) :: &
175 : skip_m_grav_and_grav, &
176 : skip_net, skip_neu, skip_kap, skip_grads, &
177 : skip_rotation, skip_brunt, &
178 : skip_mixing_info, skip_set_cz_bdy_mass, skip_irradiation_heat
179 : real(dp), intent(in) :: dt
180 : integer, intent(out) :: ierr
181 :
182 : logical, parameter :: skip_other_cgrav = .false.
183 : logical, parameter :: skip_basic_vars = .false.
184 : logical, parameter :: skip_micro_vars = .false.
185 : logical, parameter :: skip_mlt = .false.
186 : logical, parameter :: skip_eos = .false.
187 :
188 : include 'formats'
189 :
190 : call update_vars(s, &
191 : skip_basic_vars, skip_micro_vars, &
192 : skip_m_grav_and_grav, skip_net, skip_neu, skip_kap, &
193 : skip_grads, skip_rotation, skip_brunt, skip_other_cgrav, &
194 : skip_mixing_info, skip_set_cz_bdy_mass, skip_irradiation_heat, &
195 33 : skip_mlt, skip_eos, dt, ierr)
196 33 : if (ierr /= 0) then
197 0 : if (s% report_ierr .or. dbg) &
198 0 : write(*,*) 'set_some_vars: update_vars returned ierr', ierr
199 0 : return
200 : end if
201 :
202 : end subroutine set_some_vars
203 :
204 :
205 33 : subroutine update_vars(s, &
206 : skip_basic_vars, skip_micro_vars, &
207 : skip_m_grav_and_grav, skip_net, skip_neu, skip_kap, &
208 : skip_grads, skip_rotation, skip_brunt, skip_other_cgrav, &
209 : skip_mixing_info, skip_set_cz_bdy_mass, skip_irradiation_heat, &
210 : skip_mlt, skip_eos, dt, ierr)
211 : use star_utils, only: eval_irradiation_heat
212 : type (star_info), pointer :: s
213 : logical, intent(in) :: &
214 : skip_basic_vars, skip_micro_vars, &
215 : skip_m_grav_and_grav, skip_net, skip_neu, skip_kap, skip_grads, &
216 : skip_rotation, skip_brunt, skip_other_cgrav, &
217 : skip_mixing_info, skip_set_cz_bdy_mass, skip_irradiation_heat, &
218 : skip_mlt, skip_eos
219 : real(dp), intent(in) :: dt
220 : integer, intent(out) :: ierr
221 : integer :: k, nz
222 : include 'formats'
223 : ierr = 0
224 33 : nz = s% nz
225 33 : if (s% doing_finish_load_model .or. .not. s% RSP_flag) then
226 33 : call unpack_xh(s,ierr)
227 33 : if (ierr /= 0) then
228 0 : if (s% report_ierr .or. dbg) &
229 0 : write(*,*) 'after unpack ierr', ierr
230 0 : return
231 : end if
232 33 : if (.not. skip_mixing_info) then
233 40766 : s% mixing_type(1:nz) = no_mixing
234 40766 : s% adjust_mlt_gradT_fraction(1:nz) = -1
235 : end if
236 : end if
237 :
238 : call set_hydro_vars( &
239 : s, 1, nz, skip_basic_vars, &
240 : skip_micro_vars, skip_m_grav_and_grav, skip_eos, skip_net, skip_neu, &
241 : skip_kap, skip_grads, skip_rotation, skip_brunt, skip_other_cgrav, &
242 33 : skip_mixing_info, skip_set_cz_bdy_mass, skip_mlt, ierr)
243 33 : if (ierr /= 0) then
244 0 : if (s% report_ierr .or. dbg) &
245 0 : write(*,*) 'update_vars: set_hydro_vars returned ierr', ierr
246 0 : return
247 : end if
248 :
249 33 : if (s% use_other_momentum) then
250 0 : call s% other_momentum(s% id, ierr)
251 0 : if (ierr /= 0) then
252 0 : if (s% report_ierr .or. dbg) &
253 0 : write(*,*) 'update_vars: other_momentum returned ierr', ierr
254 0 : return
255 : end if
256 : end if
257 :
258 33 : if (.not. skip_irradiation_heat) then
259 33 : if (s% irradiation_flux /= 0) then
260 0 : do k=1,nz
261 0 : s% irradiation_heat(k) = eval_irradiation_heat(s,k)
262 : end do
263 : else
264 40766 : s% irradiation_heat(1:nz) = 0
265 : end if
266 : end if
267 :
268 : end subroutine update_vars
269 :
270 :
271 66 : subroutine unpack_xh(s,ierr)
272 : use star_utils, only: set_qs, set_dm_bar, set_m_and_dm
273 : type (star_info), pointer :: s
274 : integer, intent(out) :: ierr
275 : integer :: i_lnd, i_lnT, i_lnR, i_w, i_Hp, &
276 : i_lum, i_v, i_u, i_alpha_RTI, i_Et_RSP, &
277 : j, k, species, nvar_chem, nz, k_below_just_added
278 : include 'formats'
279 33 : ierr = 0
280 33 : nz = s% nz
281 33 : k_below_just_added = 1
282 33 : species = s% species
283 33 : nvar_chem = s% nvar_chem
284 33 : i_lnd = s% i_lnd
285 33 : i_lnT = s% i_lnT
286 33 : i_lnR = s% i_lnR
287 33 : i_lum = s% i_lum
288 33 : i_w = s% i_w
289 33 : i_Hp = s% i_Hp
290 33 : i_v = s% i_v
291 33 : i_u = s% i_u
292 33 : i_alpha_RTI = s% i_alpha_RTI
293 33 : i_Et_RSP = s% i_Et_RSP
294 :
295 165 : do j=1,s% nvar_hydro
296 165 : if (j == i_lnd) then
297 40766 : do k=1,nz
298 40733 : s% lnd(k) = s% xh(i_lnd,k)
299 40766 : s% rho(k) = exp(s% lnd(k))
300 : end do
301 40766 : s% dxh_lnd(1:nz) = 0d0
302 99 : else if (j == i_lnT) then
303 40766 : do k=1,nz
304 40733 : s% lnT(k) = s% xh(i_lnT,k)
305 40766 : s% T(k) = exp(s% lnT(k))
306 : end do
307 40766 : s% dxh_lnT(1:nz) = 0d0
308 66 : else if (j == i_lnR) then
309 40766 : do k=1,nz
310 40733 : s% lnR(k) = s% xh(i_lnR,k)
311 40766 : s% r(k) = exp(s% lnR(k))
312 : end do
313 40766 : s% dxh_lnR(1:nz) = 0d0
314 33 : else if (j == i_w) then
315 0 : do k=1,nz
316 0 : s% w(k) = s% xh(i_w, k)
317 0 : if (s% w(k) < 0d0) then
318 : !write(*,4) 'unpack: fix w < 0', k, &
319 : ! s% solver_iter, s% model_number, s% w(k)
320 0 : s% w(k) = s% RSP2_w_fix_if_neg
321 : end if
322 : end do
323 33 : else if (j == i_Hp) then
324 0 : do k=1,nz
325 0 : s% Hp_face(k) = s% xh(i_Hp, k)
326 : end do
327 33 : else if (j == i_lum) then
328 40766 : do k=1,nz
329 40766 : s% L(k) = s% xh(i_lum, k)
330 : end do
331 0 : else if (j == i_v) then
332 0 : do k=1,nz
333 0 : s% v(k) = s% xh(i_v,k)
334 0 : s% dxh_v(k) = 0d0
335 : end do
336 0 : else if (j == i_u) then
337 0 : do k=1,nz
338 0 : s% u(k) = s% xh(i_u,k)
339 0 : s% dxh_u(k) = 0d0
340 : end do
341 0 : else if (j == i_alpha_RTI) then
342 0 : do k=1,nz
343 0 : s% alpha_RTI(k) = max(0d0, s% xh(i_alpha_RTI,k))
344 0 : s% dxh_alpha_RTI(k) = 0d0
345 : end do
346 0 : else if (j == i_Et_RSP) then
347 0 : do k=1,nz
348 0 : s% RSP_Et(k) = max(0d0,s% xh(i_Et_RSP,k))
349 : end do
350 : end if
351 : end do
352 :
353 33 : if (i_lum == 0 .and. .not. s% RSP_flag) s% L(1:nz) = 0d0
354 40766 : if (i_v == 0) s% v(1:nz) = 0d0
355 40766 : if (i_u == 0) s% u(1:nz) = 0d0
356 40766 : if (i_w == 0) s% w(1:nz) = 0d0
357 40766 : if (i_Hp == 0) s% Hp_face(1:nz) = 0d0
358 :
359 33 : call set_qs(s, nz, s% q, s% dq, ierr)
360 33 : if (ierr /= 0) then
361 0 : write(*,*) 'update_vars failed in set_qs'
362 0 : return
363 : end if
364 33 : call set_m_and_dm(s)
365 33 : call set_dm_bar(s, s% nz, s% dm, s% dm_bar)
366 :
367 : end subroutine unpack_xh
368 :
369 :
370 0 : subroutine set_Teff(s, ierr)
371 : type (star_info), pointer :: s
372 : integer, intent(out) :: ierr
373 : real(dp) :: r_phot, L_surf
374 : logical, parameter :: skip_partials = .false., &
375 : need_atm_Psurf = .false., need_atm_Tsurf = .false.
376 : real(dp) :: Teff, &
377 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
378 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap
379 : ierr = 0
380 : call set_Teff_info_for_eqns(s, skip_partials, &
381 : need_atm_Psurf, need_atm_Tsurf, r_phot, L_surf, Teff, &
382 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
383 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
384 0 : ierr)
385 0 : end subroutine set_Teff
386 :
387 :
388 44 : subroutine set_Teff_info_for_eqns(s, skip_partials, &
389 : need_atm_Psurf_in, need_atm_Tsurf_in, r_surf, L_surf, Teff, &
390 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
391 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
392 : ierr)
393 : use star_utils, only: set_phot_info
394 : use atm_lib, only: atm_Teff
395 : use starspots, only: starspot_tweak_PT, starspot_restore_PT
396 : type (star_info), pointer :: s
397 : logical, intent(in) :: skip_partials, &
398 : need_atm_Psurf_in, need_atm_Tsurf_in
399 : real(dp), intent(out) :: r_surf, L_surf, Teff, &
400 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
401 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap
402 : integer, intent(out) :: ierr
403 : logical :: need_atm_Psurf, need_atm_Tsurf
404 :
405 : include 'formats'
406 :
407 44 : need_atm_Psurf = need_atm_Psurf_in
408 44 : need_atm_Tsurf = need_atm_Tsurf_in
409 :
410 44 : ierr = 0
411 :
412 44 : r_surf = s% r(1)
413 44 : L_surf = s% L(1)
414 :
415 44 : s% P_surf = s% Peos(1)
416 44 : s% T_surf = s% T(1)
417 :
418 44 : call set_phot_info(s) ! sets Teff using L_phot and R_phot
419 44 : Teff = s% Teff
420 :
421 44 : if (s% RSP_flag) then
422 0 : lnT_surf = s% lnT(1)
423 0 : dlnT_dL = 0d0
424 0 : dlnT_dlnR = 0d0
425 0 : dlnT_dlnM = 0d0
426 0 : dlnT_dlnkap = 0d0
427 0 : lnP_surf = s% lnPeos(1)
428 0 : dlnP_dL = 0d0
429 0 : dlnP_dlnR = 0d0
430 0 : dlnP_dlnM = 0d0
431 0 : dlnP_dlnkap = 0d0
432 0 : return
433 : end if
434 :
435 44 : if (s% use_other_surface_PT) then
436 : call s% other_surface_PT( &
437 : s% id, skip_partials, &
438 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
439 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
440 0 : ierr)
441 : else
442 : ! starspot YREC routine
443 44 : if (s% do_starspots) then
444 0 : need_atm_Psurf = .true.
445 0 : need_atm_Tsurf = .true.
446 0 : call starspot_tweak_PT(s)
447 : end if
448 : call get_surf_PT( &
449 : s, skip_partials, need_atm_Psurf, need_atm_Tsurf, &
450 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
451 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
452 44 : ierr)
453 : ! starspot YREC routine
454 44 : if (s% do_starspots) then
455 0 : call starspot_restore_PT(s)
456 : end if
457 : end if
458 44 : if (ierr /= 0) then
459 0 : if (s% report_ierr) then
460 0 : write(*,*) 'error in get_surf_PT'
461 : end if
462 0 : return
463 : end if
464 44 : s% T_surf = exp(lnT_surf)
465 44 : s% P_surf = exp(lnP_surf)
466 :
467 44 : call set_phot_info(s) ! s% T_surf might have changed so call again
468 :
469 : end subroutine set_Teff_info_for_eqns
470 :
471 :
472 154 : subroutine set_hydro_vars( &
473 : s, nzlo, nzhi, skip_basic_vars, skip_micro_vars, skip_m_grav_and_grav, &
474 : skip_eos, skip_net, skip_neu, skip_kap, skip_grads, skip_rotation, &
475 : skip_brunt, skip_other_cgrav, &
476 : skip_mixing_info, skip_set_cz_bdy_mass, skip_mlt, ierr)
477 : use atm_lib, only: atm_eval_T_tau_dq_dtau
478 : use atm_support, only: get_T_tau_id
479 : use micro, only: set_micro_vars
480 : use turb_info, only: set_mlt_vars, check_for_redo_MLT
481 : use star_utils, only: start_time, update_time, &
482 : set_m_grav_and_grav, set_scale_height, get_tau, &
483 : set_abs_du_div_cs, set_conv_time_scales
484 : use hydro_rotation, only: set_rotation_info, compute_j_fluxes_and_extra_jdot
485 : use brunt, only: do_brunt_B, do_brunt_N2
486 : use mix_info, only: set_mixing_info
487 : use hydro_rsp2, only: set_RSP2_vars
488 : use tdc_hydro, only: set_viscosity_vars_TDC
489 :
490 : type (star_info), pointer :: s
491 : integer, intent(in) :: nzlo, nzhi
492 : logical, intent(in) :: &
493 : skip_basic_vars, skip_micro_vars, skip_m_grav_and_grav, &
494 : skip_eos, skip_net, skip_neu, skip_kap, skip_brunt, &
495 : skip_grads, skip_rotation, skip_other_cgrav, &
496 : skip_mixing_info, skip_set_cz_bdy_mass, skip_mlt
497 : integer, intent(out) :: ierr
498 :
499 : integer :: nz, k, T_tau_id
500 : integer(i8) :: time0
501 : logical, parameter :: dbg = .false.
502 : real(dp) :: total
503 :
504 : include 'formats'
505 :
506 : if (dbg) write(*,*) 'set_hydro_vars', nzlo, nzhi
507 77 : if (s% doing_timing) call start_time(s, time0, total)
508 :
509 77 : ierr = 0
510 77 : nz = s% nz
511 :
512 77 : if (.not. skip_basic_vars) then
513 : if (dbg) write(*,*) 'call set_basic_vars'
514 33 : call set_basic_vars(s, nzlo, nzhi, ierr)
515 33 : if (failed('set_basic_vars')) return
516 : end if
517 :
518 77 : if (.not. skip_micro_vars) then
519 : if (dbg) write(*,*) 'call set_micro_vars'
520 : call set_micro_vars(s, nzlo, nzhi, &
521 66 : skip_eos, skip_net, skip_neu, skip_kap, ierr)
522 66 : if (failed('set_micro_vars')) return
523 : end if
524 :
525 77 : if (.not. skip_other_cgrav) call set_cgrav(s, ierr)
526 :
527 77 : call get_tau(s, ierr)
528 77 : if (failed('get_tau')) return
529 :
530 77 : if (.not. skip_m_grav_and_grav) then
531 : ! don't change m_grav or grav during solver iterations
532 : if (dbg) write(*,*) 'call set_m_grav_and_grav'
533 44 : call set_m_grav_and_grav(s)
534 : end if
535 :
536 : if (dbg) write(*,*) 'call set_scale_height'
537 77 : call set_scale_height(s)
538 :
539 77 : if (s% rotation_flag .and. (.not. skip_rotation .or. s% w_div_wc_flag)) then
540 : if (dbg) write(*,*) 'call set_rotation_info'
541 0 : call set_rotation_info(s, .true., ierr)
542 0 : if (failed('set_rotation_info')) return
543 : end if
544 :
545 77 : if (.not. skip_grads) then
546 : if (dbg) write(*,*) 'call do_brunt_B'
547 44 : call do_brunt_B(s, nzlo, nzhi, ierr) ! for unsmoothed_brunt_B
548 44 : if (failed('do_brunt_B')) return
549 : if (dbg) write(*,*) 'call set_grads'
550 44 : call set_grads(s, ierr)
551 44 : if (failed('set_grads')) return
552 44 : call set_conv_time_scales(s) ! uses brunt_B
553 : end if
554 :
555 77 : if (.not. skip_mixing_info) then
556 33 : if (.not. s% RSP2_flag) then
557 : if (dbg) write(*,*) 'call other_adjust_mlt_gradT_fraction'
558 33 : call s% other_adjust_mlt_gradT_fraction(s% id,ierr)
559 33 : if (failed('other_adjust_mlt_gradT_fraction')) return
560 : end if
561 : if (dbg) write(*,*) 'call set_abs_du_div_cs'
562 33 : call set_abs_du_div_cs(s)
563 : end if
564 :
565 77 : if (.not. skip_mlt .and. .not. s% RSP_flag) then
566 :
567 80060 : s% reconstructed_face_state_valid(1:s%nz) = .false.
568 :
569 66 : if (.not. skip_mixing_info) then
570 33 : if (s% make_gradr_sticky_in_solver_iters) then
571 0 : s% fixed_gradr_for_rest_of_solver_iters(nzlo:nzhi) = .false.
572 : end if
573 40766 : s% alpha_mlt(nzlo:nzhi) = s% mixing_length_alpha
574 33 : if (s% use_other_alpha_mlt) then
575 0 : call s% other_alpha_mlt(s% id, ierr)
576 0 : if (ierr /= 0) then
577 0 : if (s% report_ierr .or. dbg) &
578 0 : write(*,*) 'other_alpha_mlt returned ierr', ierr
579 0 : return
580 : end if
581 : end if
582 : end if
583 :
584 66 : if (s% use_other_gradr_factor) then
585 : if (dbg) write(*,*) 'call other_gradr_factor'
586 0 : call s% other_gradr_factor(s% id, ierr)
587 0 : if (failed('other_gradr_factor')) return
588 66 : else if (s% use_T_tau_gradr_factor) then
589 : if (dbg) write(*,*) 'call use_T_tau_gradr_factor'
590 0 : call get_T_tau_id(s% atm_T_tau_relation, T_tau_id, ierr)
591 0 : do k = nzlo, nzhi
592 : ! slightly inconsistent to use s% tau, defined at
593 : ! faces, with s% gradr, which is a cell average
594 : ! but this performs better than using
595 : ! taumid = (tau(k)+tau(k+1))/2
596 0 : call atm_eval_T_tau_dq_dtau(T_tau_id, s% tau(k), s% gradr_factor(k), ierr)
597 0 : s% gradr_factor(k) = 1.0_dp + s% gradr_factor(k)
598 : end do
599 0 : if (failed('use_T_tau_gradr_factor')) return
600 : else
601 80060 : s% gradr_factor(nzlo:nzhi) = 1d0
602 : end if
603 :
604 : if (s% TDC_alpha_M > 0 .and. s% MLT_option == 'TDC' &
605 66 : .and. .not. (s% RSP2_flag .or. s% RSP_flag)) then
606 0 : call set_viscosity_vars_TDC(s,ierr)
607 0 : if (ierr /= 0) then
608 0 : if (len_trim(s% retry_message) == 0) s% retry_message = 'set_viscosity_vars_TDC failed'
609 0 : if (s% report_ierr) write(*,*) 'ierr from set_viscosity_vars_TDC'
610 0 : return
611 : end if
612 : end if
613 :
614 66 : call set_mlt_vars(s, nzlo, nzhi, ierr)
615 66 : if (failed('set_mlt_vars')) return
616 : if (dbg) write(*,*) 'call check_for_redo_MLT'
617 :
618 66 : call check_for_redo_MLT(s, nzlo, nzhi, ierr)
619 66 : if (failed('check_for_redo_MLT')) return
620 :
621 : end if
622 :
623 77 : if (.not. skip_brunt) then ! skip_brunt during solver iterations
624 : if (dbg) write(*,*) 'call do_brunt_N2'
625 44 : call do_brunt_N2(s, nzlo, nzhi, ierr)
626 44 : if (failed('do_brunt_N2')) return
627 : end if
628 :
629 77 : if (.not. skip_mixing_info) then
630 : if (dbg) write(*,*) 'call set_mixing_info'
631 33 : call set_mixing_info(s, skip_set_cz_bdy_mass, ierr)
632 33 : if (ierr /= 0) return
633 : end if
634 :
635 77 : if (s% j_rot_flag) then
636 0 : call compute_j_fluxes_and_extra_jdot(s% id, ierr)
637 0 : if (ierr /= 0) then
638 0 : write(*,*) 'failed in compute_j_fluxes'
639 : end if
640 : end if
641 :
642 77 : if (s% RSP2_flag) then
643 0 : call set_RSP2_vars(s,ierr)
644 0 : if (ierr /= 0) then
645 0 : if (len_trim(s% retry_message) == 0) s% retry_message = 'set_RSP2_vars failed'
646 0 : if (s% report_ierr) write(*,*) 'ierr from set_RSP2_vars'
647 0 : return
648 : end if
649 : end if
650 :
651 :
652 77 : if (s% doing_timing) &
653 0 : call update_time(s, time0, total, s% time_set_hydro_vars)
654 :
655 77 : s% need_to_setvars = .false.
656 :
657 :
658 : contains
659 :
660 473 : logical function failed(str)
661 : character (len=*), intent(in) :: str
662 473 : if (ierr == 0) then
663 : failed = .false.
664 : return
665 : end if
666 0 : if (s% report_ierr .or. dbg) &
667 0 : write(*,*) 'set_hydro_vars failed in call to ' // trim(str)
668 : failed = .true.
669 : end function failed
670 :
671 : end subroutine set_hydro_vars
672 :
673 :
674 : subroutine check_rs(s, ierr)
675 : type (star_info), pointer :: s
676 : integer, intent(out) :: ierr
677 : integer :: k
678 : logical :: okay
679 : include 'formats'
680 : ierr = 0
681 : okay = .true.
682 : do k=2, s% nz
683 : if (s% r(k) > s% r(k-1)) then
684 : if (s% report_ierr) then
685 : write(*,2) 's% r(k) > s% r(k-1)', k, &
686 : s% r(k)/Rsun, s% r(k-1)/Rsun, s% r(k)/Rsun-s% r(k-1)/Rsun
687 : end if
688 : okay = .false.
689 : end if
690 : end do
691 : if (okay) return
692 : s% retry_message = 'invalid values for r'
693 : ierr = -1
694 : if (s% report_ierr) write(*,*)
695 : end subroutine check_rs
696 :
697 :
698 66 : subroutine set_basic_vars(s, nzlo, nzhi, ierr)
699 : use star_utils, only: set_rv_info, set_rmid
700 : type (star_info), pointer :: s
701 : integer, intent(in) :: nzlo, nzhi
702 : integer, intent(out) :: ierr
703 : integer :: j, k, species, nz
704 : real(dp) :: twoGmrc2, sum_xa
705 :
706 : include 'formats'
707 :
708 : if (dbg) write(*,4) 'enter set_basic_vars: nzlo, nzhi, nz', &
709 : nzlo, nzhi, s% nz
710 33 : ierr = 0
711 33 : nz = s% nz
712 33 : species = s% species
713 33 : s% L_phot = s% L(1)/Lsun
714 33 : if (.not. s% using_velocity_time_centering) then
715 33 : s% d_vc_dv = 1d0
716 : else
717 0 : s% d_vc_dv = 0.5d0
718 : end if
719 :
720 33 : !$OMP PARALLEL DO PRIVATE(j,k,twoGmrc2,sum_xa) SCHEDULE(dynamic,2)
721 : do k=nzlo, nzhi
722 : if (s% lnd_start(k) < -1d90) then
723 : s% lnd_start(k) = s% lnd(k)
724 : s% rho_start(k) = s% rho(k)
725 : end if
726 : if (s% T_start(k) < 0d0) then
727 : s% T_start(k) = s% T(k)
728 : s% lnT_start(k) = s% lnT(k)
729 : end if
730 : if (s% v_flag) then
731 : if (s% v_start(k) < -1d90) s% v_start(k) = s% v(k)
732 : end if
733 : if (s% u_flag) then
734 : if (s% u_start(k) < -1d90) s% u_start(k) = s% u(k)
735 : end if
736 : if (s% RTI_flag) then
737 : if (s% alpha_RTI_start(k) < -1d90) &
738 : s% alpha_RTI_start(k) = s% alpha_RTI(k)
739 : end if
740 : if (s% RSP_flag) then
741 : s% RSP_w(k) = sqrt(s% RSP_Et(k))
742 : if (s% RSP_w_start(k) < -1d90) then
743 : s% RSP_w_start(k) = s% RSP_w(k)
744 : end if
745 : end if
746 : if (s% r_start(k) < 0) s% r_start(k) = s% r(k)
747 : call set_rv_info(s,k)
748 : do j=1,species
749 : s% xa(j,k) = max(0d0, min(1d0, s% xa(j,k)))
750 : end do
751 : sum_xa = sum(s% xa(1:species,k))
752 : if (abs(sum_xa - 1d0) > 1d-12) then
753 : do j=1,species
754 : s% xa(j,k) = s% xa(j,k)/sum_xa
755 : end do
756 : end if
757 : end do
758 : !$OMP END PARALLEL DO
759 :
760 33 : call set_rmid(s, nzlo, nzhi, ierr)
761 :
762 33 : end subroutine set_basic_vars
763 :
764 :
765 44 : subroutine set_cgrav(s, ierr)
766 : type (star_info), pointer :: s
767 : integer, intent(out) :: ierr
768 44 : if (s% use_other_cgrav) then
769 0 : call s% other_cgrav(s% id, ierr)
770 0 : if (ierr /= 0) then
771 0 : if (s% report_ierr .or. dbg) &
772 0 : write(*,*) 'other_cgrav returned ierr', ierr
773 0 : return
774 : end if
775 : else
776 53864 : s% cgrav(1:s% nz) = standard_cgrav
777 : end if
778 : end subroutine set_cgrav
779 :
780 :
781 44 : subroutine get_surf_PT( &
782 : s, skip_partials, &
783 : need_atm_Psurf, need_atm_Tsurf, &
784 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
785 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
786 : ierr)
787 :
788 : use atm_support, only: get_atm_PT
789 : use atm_def, only: star_debugging_atm_flag, &
790 : atm_test_partials_val, atm_test_partials_dval_dx
791 : use chem_def
792 : use eos_lib, only: Radiation_Pressure
793 :
794 : type (star_info), pointer :: s
795 : logical, intent(in) :: skip_partials, &
796 : need_atm_Psurf, need_atm_Tsurf
797 : real(dp), intent(out) :: &
798 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
799 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap
800 : integer, intent(out) :: ierr
801 :
802 : real(dp) :: L_surf
803 : real(dp) :: R_surf
804 : real(dp) :: Teff
805 : real(dp) :: tau_surf
806 : real(dp) :: Teff4
807 : real(dp) :: T_surf4
808 : real(dp) :: T_surf
809 : real(dp) :: P_surf_atm
810 : real(dp) :: P_surf
811 : real(dp) :: Pextra
812 : real(dp) :: kap_surf
813 : real(dp) :: M_surf
814 :
815 : include 'formats'
816 :
817 : ! Set up stellar surface parameters
818 :
819 44 : L_surf = s% L(1)
820 44 : R_surf = s% r(1)
821 44 : kap_surf = s% opacity(1)
822 44 : M_surf = s% m(1)
823 44 : Teff = s% Teff
824 :
825 : ! Initialize partials
826 44 : dlnT_dL = 0._dp; dlnT_dlnR = 0._dp; dlnT_dlnM = 0._dp; dlnT_dlnkap = 0._dp
827 44 : dlnP_dL = 0._dp; dlnP_dlnR = 0._dp; dlnP_dlnM = 0._dp; dlnP_dlnkap = 0._dp
828 :
829 : ! Evaluate the surface optical depth
830 :
831 44 : tau_surf = s% tau_factor*s% tau_base ! tau at outer edge of cell 1
832 44 : if (is_bad(tau_surf)) then
833 0 : write(*,1) 's% tau_factor', s% tau_factor
834 0 : write(*,1) 's% tau_base', s% tau_base
835 0 : write(*,1) 'tau_surf', tau_surf
836 0 : call mesa_error(__FILE__,__LINE__,'bad tau_surf in get_surf_PT')
837 : end if
838 :
839 : ! Evaluate surface temperature and pressure
840 :
841 44 : if (.not. (need_atm_Psurf .or. need_atm_Tsurf)) then
842 :
843 : ! Special-case boundary condition
844 :
845 0 : lnP_surf = s% lnPeos_start(1)
846 0 : if (is_bad(lnP_surf)) lnP_surf = 0._dp
847 0 : T_surf = s% T_start(1)
848 0 : lnT_surf = log(T_surf)
849 0 : T_surf4 = T_surf*T_surf*T_surf*T_surf
850 0 : Teff = pow(four_thirds*T_surf4/(tau_surf + two_thirds), 0.25_dp)
851 :
852 0 : if (.not. skip_partials) then
853 0 : dlnT_dL = 0._dp; dlnT_dlnR = 0._dp; dlnT_dlnM = 0._dp; dlnT_dlnkap = 0._dp
854 0 : dlnP_dL = 0._dp; dlnP_dlnR = 0._dp; dlnP_dlnM = 0._dp; dlnP_dlnkap = 0._dp
855 : end if
856 :
857 : else
858 : ! Evaluate temperature and pressure based on atm_option
859 : ! (yes, we might want to translate atm_option into an integer,
860 : ! but these string comparisons are cheap)
861 :
862 : ! The first few are special, 'trivial-atmosphere' options
863 :
864 0 : select case (s% atm_option)
865 :
866 : case ('fixed_Teff')
867 :
868 : ! set Tsurf from Eddington T-tau relation
869 : ! for current surface tau and Teff = `atm_fixed_Teff`.
870 : ! set Psurf = Radiation_Pressure(Tsurf)
871 :
872 0 : Teff = s% atm_fixed_Teff
873 0 : Teff4 = Teff*Teff*Teff*Teff
874 0 : T_surf = pow(3._dp/4._dp*Teff4*(tau_surf + 2._dp/3._dp), 0.25_dp)
875 0 : lnT_surf = log(T_surf)
876 0 : lnP_surf = Radiation_Pressure(T_surf)
877 :
878 0 : if (.not. skip_partials) then
879 0 : dlnT_dL = 0._dp; dlnT_dlnR = 0._dp; dlnT_dlnM = 0._dp; dlnT_dlnkap = 0._dp
880 0 : dlnP_dL = 0._dp; dlnP_dlnR = 0._dp; dlnP_dlnM = 0._dp; dlnP_dlnkap = 0._dp
881 : end if
882 :
883 : case ('fixed_Tsurf')
884 :
885 : ! set Teff from Eddington T-tau relation for given
886 : ! Tsurf and tau=2/3 set Psurf =
887 : ! Radiation_Pressure(Tsurf)
888 :
889 0 : T_surf = s% atm_fixed_Tsurf
890 0 : lnT_surf = log(T_surf)
891 0 : T_surf4 = T_surf*T_surf*T_surf*T_surf
892 0 : Teff = pow(4._dp/3._dp*T_surf4/(tau_surf + 2._dp/3._dp), 0.25_dp)
893 0 : lnP_surf = Radiation_Pressure(T_surf)
894 :
895 0 : if (.not. skip_partials) then
896 0 : dlnT_dL = 0._dp; dlnT_dlnR = 0._dp; dlnT_dlnM = 0._dp; dlnT_dlnkap = 0._dp
897 0 : dlnP_dL = 0._dp; dlnP_dlnR = 0._dp; dlnP_dlnM = 0._dp; dlnP_dlnkap = 0._dp
898 : end if
899 :
900 : case ('fixed_Psurf')
901 :
902 : ! set Teff from L and R using L = 4*pi*R^2*boltz_sigma*T^4.
903 : ! set Tsurf using Eddington T-tau relation
904 :
905 0 : if (L_surf < 0._dp) then
906 0 : if (s% report_ierr) then
907 0 : write(*,2) 'get_surf_PT: L_surf <= 0', s% model_number, L_surf
908 : end if
909 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__)
910 0 : s% retry_message = 'L_surf < 0'
911 0 : ierr = -1
912 0 : return
913 : end if
914 :
915 0 : lnP_surf = safe_log(s% atm_fixed_Psurf)
916 0 : if (R_surf <= 0._dp) then
917 0 : T_surf4 = 1._dp
918 0 : T_surf = 1._dp
919 0 : lnT_surf = 0._dp
920 0 : if (.not. skip_partials) then
921 0 : dlnT_dlnR = 0._dp
922 0 : dlnT_dL = 0._dp
923 : end if
924 0 : Teff = s% T(1)
925 : else
926 0 : Teff = pow(L_surf/(4._dp*pi*R_surf*R_surf*boltz_sigma), 0.25_dp)
927 0 : T_surf4 = 3d0/4d0*pow(Teff, 4.d0)*(tau_surf + 2._dp/3._dp)
928 0 : T_surf = pow(T_surf4, 0.25_dp)
929 0 : lnT_surf = log(T_surf)
930 0 : if (.not. skip_partials) then
931 0 : dlnT_dlnR = -0.5_dp
932 0 : dlnT_dL = 1._dp/(4._dp*L_surf)
933 : end if
934 : end if
935 :
936 0 : if (.not. skip_partials) then
937 0 : dlnT_dlnM = 0._dp; dlnT_dlnkap = 0._dp
938 0 : dlnP_dL = 0._dp; dlnP_dlnR = 0._dp; dlnP_dlnM = 0._dp; dlnP_dlnkap = 0._dp
939 : end if
940 :
941 : case ('fixed_Psurf_and_Tsurf')
942 :
943 0 : lnP_surf = safe_log(s% atm_fixed_Psurf)
944 0 : T_surf = s% atm_fixed_Tsurf
945 0 : lnT_surf = log(T_surf)
946 0 : T_surf4 = T_surf*T_surf*T_surf*T_surf
947 0 : Teff = pow(4d0/3d0*T_surf4/(tau_surf + 2d0/3d0), 0.25d0)
948 :
949 0 : if (.not. skip_partials) then
950 0 : dlnT_dL = 0; dlnT_dlnR = 0; dlnT_dlnM = 0; dlnT_dlnkap = 0
951 0 : dlnP_dL = 0; dlnP_dlnR = 0; dlnP_dlnM = 0; dlnP_dlnkap = 0
952 : end if
953 :
954 : case default
955 :
956 : ! Everything else -- the 'non-trivial atmospheres' ---
957 : ! gets passed to atm_support
958 :
959 : if (.false. .and. 1 == s% solver_test_partials_k .and. &
960 : s% solver_iter == s% solver_test_partials_iter_number) then
961 : star_debugging_atm_flag = .true.
962 : end if
963 :
964 : call get_atm_PT( &
965 : s, tau_surf, L_surf, R_surf, s% m(1), s% cgrav(1), skip_partials, &
966 : Teff, lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
967 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
968 44 : ierr)
969 44 : if (ierr /= 0) then
970 0 : if (s% report_ierr) then
971 0 : write(*,1) 'tau_surf', tau_surf
972 0 : write(*,1) 'L_surf', L_surf
973 0 : write(*,1) 'R_surf', R_surf
974 0 : write(*,1) 'Teff', Teff
975 0 : write(*,1) 's% m(1)', s% m(1)
976 0 : write(*,1) 's% cgrav(1)', s% cgrav(1)
977 0 : write(*,*) 'failed in get_atm_PT'
978 : end if
979 0 : return
980 : end if
981 :
982 : if (.false. .and. 1 == s% solver_test_partials_k .and. &
983 44 : s% solver_iter == s% solver_test_partials_iter_number) then
984 : s% solver_test_partials_val = atm_test_partials_val
985 : s% solver_test_partials_dval_dx = atm_test_partials_dval_dx
986 : end if
987 : end select
988 : end if
989 :
990 : ! if using fixed surface, calculate Pextra.
991 : if (s% atm_option == 'fixed_Tsurf' .or. s% atm_option == 'fixed_Psurf_and_Tsurf' &
992 44 : .or. s% atm_option == 'fixed_Psurf' .or. s% atm_option == 'fixed_Teff') then
993 : ! add extra pressure for fixed atmosphere options
994 0 : if (s% Pextra_factor /= 0._dp) then
995 0 : P_surf_atm = exp(lnP_surf)
996 0 : Pextra = s% Pextra_factor*(kap_surf/tau_surf)*(L_surf/M_surf)/(6._dp*pi*clight*s% cgrav(1))
997 0 : P_surf = P_surf_atm + Pextra
998 0 : if (P_surf < 1E-50_dp) then
999 0 : lnP_surf = -50*ln10
1000 0 : if (.not. skip_partials) then
1001 0 : dlnP_dL = 0._dp
1002 0 : dlnP_dlnR = 0._dp
1003 0 : dlnP_dlnM = 0._dp
1004 0 : dlnP_dlnkap = 0._dp
1005 : end if
1006 : else
1007 0 : lnP_surf = log(P_surf)
1008 0 : if (.not. skip_partials) then
1009 0 : dlnP_dL = dlnP_dL*P_surf_atm/P_surf
1010 0 : dlnP_dlnR = dlnP_dlnR*P_surf_atm/P_surf
1011 0 : dlnP_dlnM = dlnP_dlnM*P_surf_atm/P_surf
1012 0 : dlnP_dlnkap = dlnP_dlnkap*P_surf_atm/P_surf
1013 : end if
1014 : end if
1015 : end if
1016 : end if
1017 :
1018 : ! Check outputs
1019 :
1020 44 : if (is_bad(lnT_surf) .or. is_bad(lnP_surf)) then
1021 0 : if (len_trim(s% retry_message) == 0) s% retry_message = 'bad logT surf or logP surf'
1022 0 : ierr = -1
1023 0 : write(*,1) 'bad outputs in get_surf_PT'
1024 0 : write(*,1) 'lnT_surf', lnT_surf
1025 0 : write(*,1) 'lnP_surf', lnP_surf
1026 0 : write(*,*) 'atm_option = ', trim(s% atm_option)
1027 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'get PT surf')
1028 : end if
1029 :
1030 : return
1031 :
1032 : end subroutine get_surf_PT
1033 :
1034 :
1035 88 : subroutine set_grads(s, ierr)
1036 : use chem_def, only: chem_isos
1037 : use star_utils, only: smooth, safe_div_val
1038 : type (star_info), pointer :: s
1039 : integer, intent(out) :: ierr
1040 :
1041 : integer :: k, nz, j, cid, max_cid
1042 : real(dp) :: val, max_val, A, Z
1043 44 : real(dp), pointer, dimension(:) :: dlnP, dlnd, dlnT
1044 :
1045 : include 'formats'
1046 :
1047 44 : ierr = 0
1048 44 : nz = s% nz
1049 : call do_alloc(ierr)
1050 44 : if (ierr /= 0) return
1051 :
1052 53820 : do k = 2, nz
1053 53776 : dlnP(k) = s% lnPeos(k-1) - s% lnPeos(k)
1054 53776 : dlnd(k) = s% lnd(k-1) - s% lnd(k)
1055 53820 : dlnT(k) = s% lnT(k-1) - s% lnT(k)
1056 : end do
1057 44 : dlnP(1) = dlnP(2)
1058 44 : dlnd(1) = dlnd(2)
1059 44 : dlnT(1) = dlnT(2)
1060 :
1061 44 : call smooth(dlnP,nz)
1062 44 : call smooth(dlnd,nz)
1063 44 : call smooth(dlnT,nz)
1064 :
1065 44 : s% grad_density(1) = 0
1066 44 : s% grad_temperature(1) = 0
1067 53820 : do k = 2, nz
1068 53820 : if (dlnP(k) >= 0) then
1069 2 : s% grad_density(k) = 0
1070 2 : s% grad_temperature(k) = 0
1071 : else
1072 53774 : s% grad_density(k) = safe_div_val(s, dlnd(k), dlnP(k))
1073 53774 : s% grad_temperature(k) = safe_div_val(s, dlnT(k), dlnP(k))
1074 : end if
1075 : end do
1076 :
1077 44 : call smooth(s% grad_density,nz)
1078 44 : call smooth(s% grad_temperature,nz)
1079 :
1080 : ! this will compute the values of s% smoothed_brunt_B
1081 44 : if (s% calculate_Brunt_B) then
1082 53864 : do k=1,nz
1083 53864 : s% smoothed_brunt_B(k) = s% unsmoothed_brunt_B(k)
1084 : end do
1085 44 : call compute_smoothed_brunt_B
1086 : else
1087 0 : s% smoothed_brunt_B(:) = 0d0
1088 : end if
1089 44 : if (s% use_Ledoux_criterion .and. s% calculate_Brunt_B) then
1090 0 : do k=1,nz
1091 0 : s% gradL_composition_term(k) = s% smoothed_brunt_B(k)
1092 : end do
1093 : else
1094 53864 : do k=1,nz
1095 53864 : s% gradL_composition_term(k) = 0d0
1096 : end do
1097 : end if
1098 :
1099 : call dealloc
1100 :
1101 53688 : do k=3,nz-2
1102 53644 : max_cid = 0
1103 53644 : max_val = -1d99
1104 482796 : do j=1,s% species
1105 429152 : cid = s% chem_id(j)
1106 429152 : A = dble(chem_isos% Z_plus_N(cid))
1107 429152 : Z = dble(chem_isos% Z(cid))
1108 429152 : val = (s% xa(j,k-2) + s% xa(j,k-1) - s% xa(j,k) - s% xa(j,k+1))*(1d0 + Z)/A
1109 482796 : if (val > max_val) then
1110 66336 : max_val = val
1111 66336 : max_cid = cid
1112 : end if
1113 : end do
1114 53688 : s% dominant_iso_for_thermohaline(k) = max_cid
1115 : end do
1116 : s% dominant_iso_for_thermohaline(1:2) = &
1117 132 : s% dominant_iso_for_thermohaline(3)
1118 : s% dominant_iso_for_thermohaline(nz-1:nz) = &
1119 132 : s% dominant_iso_for_thermohaline(nz-2)
1120 :
1121 :
1122 : contains
1123 :
1124 44 : subroutine compute_smoothed_brunt_B
1125 : use star_utils, only: weighted_smoothing, threshold_smoothing
1126 : logical, parameter :: preserve_sign = .false.
1127 : real(dp), pointer, dimension(:) :: work
1128 : include 'formats'
1129 44 : ierr = 0
1130 44 : work => dlnd
1131 44 : if (s% num_cells_for_smooth_gradL_composition_term <= 0) return
1132 : call threshold_smoothing( &
1133 : s% smoothed_brunt_B, s% threshold_for_smooth_gradL_composition_term, s% nz, &
1134 44 : s% num_cells_for_smooth_gradL_composition_term, preserve_sign, work)
1135 : end subroutine compute_smoothed_brunt_B
1136 :
1137 44 : subroutine do_alloc(ierr)
1138 : integer, intent(out) :: ierr
1139 44 : call do_work_arrays(.true.,ierr)
1140 : end subroutine do_alloc
1141 :
1142 44 : subroutine dealloc
1143 44 : call do_work_arrays(.false.,ierr)
1144 : end subroutine dealloc
1145 :
1146 88 : subroutine do_work_arrays(alloc_flag, ierr)
1147 : use alloc, only: work_array
1148 : logical, intent(in) :: alloc_flag
1149 : integer, intent(out) :: ierr
1150 : logical, parameter :: crit = .false.
1151 : ierr = 0
1152 : call work_array(s, alloc_flag, crit, &
1153 88 : dlnP, nz, nz_alloc_extra, 'mlt', ierr)
1154 88 : if (ierr /= 0) return
1155 : call work_array(s, alloc_flag, crit, &
1156 88 : dlnd, nz, nz_alloc_extra, 'mlt', ierr)
1157 88 : if (ierr /= 0) return
1158 : call work_array(s, alloc_flag, crit, &
1159 88 : dlnT, nz, nz_alloc_extra, 'mlt', ierr)
1160 : if (ierr /= 0) return
1161 : end subroutine do_work_arrays
1162 :
1163 : end subroutine set_grads
1164 :
1165 : end module hydro_vars
|
