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 33 : 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 : 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 11 : 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 11 : 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 : return
200 : end if
201 :
202 33 : 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 33 : 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 33 : end subroutine update_vars
269 :
270 :
271 33 : subroutine unpack_xh(s,ierr)
272 33 : 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 : 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 33 : 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 33 : 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 44 : end subroutine set_Teff_info_for_eqns
470 :
471 :
472 77 : 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 44 : 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 :
489 : type (star_info), pointer :: s
490 : integer, intent(in) :: nzlo, nzhi
491 : logical, intent(in) :: &
492 : skip_basic_vars, skip_micro_vars, skip_m_grav_and_grav, &
493 : skip_eos, skip_net, skip_neu, skip_kap, skip_brunt, &
494 : skip_grads, skip_rotation, skip_other_cgrav, &
495 : skip_mixing_info, skip_set_cz_bdy_mass, skip_mlt
496 : integer, intent(out) :: ierr
497 :
498 : integer :: nz, k, T_tau_id
499 : integer(i8) :: time0
500 : logical, parameter :: dbg = .false.
501 77 : real(dp) :: total
502 :
503 : include 'formats'
504 :
505 : if (dbg) write(*,*) 'set_hydro_vars', nzlo, nzhi
506 77 : if (s% doing_timing) call start_time(s, time0, total)
507 :
508 77 : ierr = 0
509 77 : nz = s% nz
510 :
511 77 : if (.not. skip_basic_vars) then
512 : if (dbg) write(*,*) 'call set_basic_vars'
513 33 : call set_basic_vars(s, nzlo, nzhi, ierr)
514 33 : if (failed('set_basic_vars')) return
515 : end if
516 :
517 77 : if (.not. skip_micro_vars) then
518 : if (dbg) write(*,*) 'call set_micro_vars'
519 : call set_micro_vars(s, nzlo, nzhi, &
520 66 : skip_eos, skip_net, skip_neu, skip_kap, ierr)
521 66 : if (failed('set_micro_vars')) return
522 : end if
523 :
524 77 : if (.not. skip_other_cgrav) call set_cgrav(s, ierr)
525 :
526 77 : call get_tau(s, ierr)
527 77 : if (failed('get_tau')) return
528 :
529 77 : if (.not. skip_m_grav_and_grav) then
530 : ! don't change m_grav or grav during solver iterations
531 : if (dbg) write(*,*) 'call set_m_grav_and_grav'
532 44 : call set_m_grav_and_grav(s)
533 : end if
534 :
535 : if (dbg) write(*,*) 'call set_scale_height'
536 77 : call set_scale_height(s)
537 :
538 77 : if (s% rotation_flag .and. (.not. skip_rotation .or. s% w_div_wc_flag)) then
539 : if (dbg) write(*,*) 'call set_rotation_info'
540 0 : call set_rotation_info(s, .true., ierr)
541 0 : if (failed('set_rotation_info')) return
542 : end if
543 :
544 77 : if (.not. skip_grads) then
545 : if (dbg) write(*,*) 'call do_brunt_B'
546 44 : call do_brunt_B(s, nzlo, nzhi, ierr) ! for unsmoothed_brunt_B
547 44 : if (failed('do_brunt_B')) return
548 : if (dbg) write(*,*) 'call set_grads'
549 44 : call set_grads(s, ierr)
550 44 : if (failed('set_grads')) return
551 44 : call set_conv_time_scales(s) ! uses brunt_B
552 : end if
553 :
554 77 : if (.not. skip_mixing_info) then
555 33 : if (.not. s% RSP2_flag) then
556 : if (dbg) write(*,*) 'call other_adjust_mlt_gradT_fraction'
557 33 : call s% other_adjust_mlt_gradT_fraction(s% id,ierr)
558 33 : if (failed('other_adjust_mlt_gradT_fraction')) return
559 : end if
560 : if (dbg) write(*,*) 'call set_abs_du_div_cs'
561 33 : call set_abs_du_div_cs(s)
562 : end if
563 :
564 77 : if (.not. skip_mlt .and. .not. s% RSP_flag) then
565 :
566 66 : if (.not. skip_mixing_info) then
567 33 : if (s% make_gradr_sticky_in_solver_iters) then
568 0 : s% fixed_gradr_for_rest_of_solver_iters(nzlo:nzhi) = .false.
569 : end if
570 40766 : s% alpha_mlt(nzlo:nzhi) = s% mixing_length_alpha
571 33 : if (s% use_other_alpha_mlt) then
572 0 : call s% other_alpha_mlt(s% id, ierr)
573 0 : if (ierr /= 0) then
574 0 : if (s% report_ierr .or. dbg) &
575 0 : write(*,*) 'other_alpha_mlt returned ierr', ierr
576 0 : return
577 : end if
578 : end if
579 : end if
580 :
581 66 : if (s% use_other_gradr_factor) then
582 : if (dbg) write(*,*) 'call other_gradr_factor'
583 0 : call s% other_gradr_factor(s% id, ierr)
584 0 : if (failed('other_gradr_factor')) return
585 66 : else if (s% use_T_tau_gradr_factor) then
586 : if (dbg) write(*,*) 'call use_T_tau_gradr_factor'
587 0 : call get_T_tau_id(s% atm_T_tau_relation, T_tau_id, ierr)
588 0 : do k = nzlo, nzhi
589 : ! slightly inconsistent to use s% tau, defined at
590 : ! faces, with s% gradr, which is a cell average
591 : ! but this performs better than using
592 : ! taumid = (tau(k)+tau(k+1))/2
593 0 : call atm_eval_T_tau_dq_dtau(T_tau_id, s% tau(k), s% gradr_factor(k), ierr)
594 0 : s% gradr_factor(k) = 1.0_dp + s% gradr_factor(k)
595 : end do
596 0 : if (failed('use_T_tau_gradr_factor')) return
597 : else
598 80060 : s% gradr_factor(nzlo:nzhi) = 1d0
599 : end if
600 :
601 66 : call set_mlt_vars(s, nzlo, nzhi, ierr)
602 66 : if (failed('set_mlt_vars')) return
603 : if (dbg) write(*,*) 'call check_for_redo_MLT'
604 :
605 66 : call check_for_redo_MLT(s, nzlo, nzhi, ierr)
606 66 : if (failed('check_for_redo_MLT')) return
607 :
608 : end if
609 :
610 77 : if (.not. skip_brunt) then ! skip_brunt during solver iterations
611 : if (dbg) write(*,*) 'call do_brunt_N2'
612 44 : call do_brunt_N2(s, nzlo, nzhi, ierr)
613 44 : if (failed('do_brunt_N2')) return
614 : end if
615 :
616 77 : if (.not. skip_mixing_info) then
617 : if (dbg) write(*,*) 'call set_mixing_info'
618 33 : call set_mixing_info(s, skip_set_cz_bdy_mass, ierr)
619 33 : if (ierr /= 0) return
620 : end if
621 :
622 77 : if (s% j_rot_flag) then
623 0 : call compute_j_fluxes_and_extra_jdot(s% id, ierr)
624 0 : if (ierr /= 0) then
625 0 : write(*,*) 'failed in compute_j_fluxes'
626 : end if
627 : end if
628 :
629 77 : if (s% RSP2_flag) then
630 0 : call set_RSP2_vars(s,ierr)
631 0 : if (ierr /= 0) then
632 0 : if (len_trim(s% retry_message) == 0) s% retry_message = 'set_RSP2_vars failed'
633 0 : if (s% report_ierr) write(*,*) 'ierr from set_RSP2_vars'
634 0 : return
635 : end if
636 : end if
637 :
638 77 : if (s% doing_timing) &
639 0 : call update_time(s, time0, total, s% time_set_hydro_vars)
640 :
641 77 : s% need_to_setvars = .false.
642 :
643 :
644 : contains
645 :
646 473 : logical function failed(str)
647 : character (len=*), intent(in) :: str
648 473 : if (ierr == 0) then
649 : failed = .false.
650 : return
651 : end if
652 0 : if (s% report_ierr .or. dbg) &
653 0 : write(*,*) 'set_hydro_vars failed in call to ' // trim(str)
654 : failed = .true.
655 77 : end function failed
656 :
657 : end subroutine set_hydro_vars
658 :
659 :
660 : subroutine check_rs(s, ierr)
661 : type (star_info), pointer :: s
662 : integer, intent(out) :: ierr
663 : integer :: k
664 : logical :: okay
665 : include 'formats'
666 : ierr = 0
667 : okay = .true.
668 : do k=2, s% nz
669 : if (s% r(k) > s% r(k-1)) then
670 : if (s% report_ierr) then
671 : write(*,2) 's% r(k) > s% r(k-1)', k, &
672 : s% r(k)/Rsun, s% r(k-1)/Rsun, s% r(k)/Rsun-s% r(k-1)/Rsun
673 : end if
674 : okay = .false.
675 : end if
676 : end do
677 : if (okay) return
678 : s% retry_message = 'invalid values for r'
679 : ierr = -1
680 : if (s% report_ierr) write(*,*)
681 : end subroutine check_rs
682 :
683 :
684 33 : subroutine set_basic_vars(s, nzlo, nzhi, ierr)
685 : use star_utils, only: set_rv_info, set_rmid
686 : type (star_info), pointer :: s
687 : integer, intent(in) :: nzlo, nzhi
688 : integer, intent(out) :: ierr
689 : integer :: j, k, species, nz
690 33 : real(dp) :: twoGmrc2, sum_xa
691 :
692 : include 'formats'
693 :
694 : if (dbg) write(*,4) 'enter set_basic_vars: nzlo, nzhi, nz', &
695 : nzlo, nzhi, s% nz
696 33 : ierr = 0
697 33 : nz = s% nz
698 33 : species = s% species
699 33 : s% L_phot = s% L(1)/Lsun
700 33 : if (.not. s% using_velocity_time_centering) then
701 33 : s% d_vc_dv = 1d0
702 : else
703 0 : s% d_vc_dv = 0.5d0
704 : end if
705 :
706 33 : !$OMP PARALLEL DO PRIVATE(j,k,twoGmrc2,sum_xa) SCHEDULE(dynamic,2)
707 : do k=nzlo, nzhi
708 : if (s% lnd_start(k) < -1d90) then
709 : s% lnd_start(k) = s% lnd(k)
710 : s% rho_start(k) = s% rho(k)
711 : end if
712 : if (s% T_start(k) < 0d0) then
713 : s% T_start(k) = s% T(k)
714 : s% lnT_start(k) = s% lnT(k)
715 : end if
716 : if (s% v_flag) then
717 : if (s% v_start(k) < -1d90) s% v_start(k) = s% v(k)
718 : end if
719 : if (s% u_flag) then
720 : if (s% u_start(k) < -1d90) s% u_start(k) = s% u(k)
721 : end if
722 : if (s% RTI_flag) then
723 : if (s% alpha_RTI_start(k) < -1d90) &
724 : s% alpha_RTI_start(k) = s% alpha_RTI(k)
725 : end if
726 : if (s% RSP_flag) then
727 : s% RSP_w(k) = sqrt(s% RSP_Et(k))
728 : if (s% RSP_w_start(k) < -1d90) then
729 : s% RSP_w_start(k) = s% RSP_w(k)
730 : end if
731 : end if
732 : if (s% r_start(k) < 0) s% r_start(k) = s% r(k)
733 : call set_rv_info(s,k)
734 : do j=1,species
735 : s% xa(j,k) = max(0d0, min(1d0, s% xa(j,k)))
736 : end do
737 : sum_xa = sum(s% xa(1:species,k))
738 : if (abs(sum_xa - 1d0) > 1d-12) then
739 : do j=1,species
740 : s% xa(j,k) = s% xa(j,k)/sum_xa
741 : end do
742 : end if
743 : end do
744 : !$OMP END PARALLEL DO
745 :
746 33 : call set_rmid(s, nzlo, nzhi, ierr)
747 :
748 33 : end subroutine set_basic_vars
749 :
750 :
751 44 : subroutine set_cgrav(s, ierr)
752 : type (star_info), pointer :: s
753 : integer, intent(out) :: ierr
754 44 : if (s% use_other_cgrav) then
755 0 : call s% other_cgrav(s% id, ierr)
756 0 : if (ierr /= 0) then
757 0 : if (s% report_ierr .or. dbg) &
758 0 : write(*,*) 'other_cgrav returned ierr', ierr
759 0 : return
760 : end if
761 : else
762 53864 : s% cgrav(1:s% nz) = standard_cgrav
763 : end if
764 33 : end subroutine set_cgrav
765 :
766 :
767 44 : subroutine get_surf_PT( &
768 : s, skip_partials, &
769 : need_atm_Psurf, need_atm_Tsurf, &
770 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
771 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
772 : ierr)
773 :
774 : use atm_support, only: get_atm_PT
775 : use atm_def, only: star_debugging_atm_flag, &
776 : atm_test_partials_val, atm_test_partials_dval_dx
777 : use chem_def
778 : use eos_lib, only: Radiation_Pressure
779 :
780 : type (star_info), pointer :: s
781 : logical, intent(in) :: skip_partials, &
782 : need_atm_Psurf, need_atm_Tsurf
783 : real(dp), intent(out) :: &
784 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
785 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap
786 : integer, intent(out) :: ierr
787 :
788 44 : real(dp) :: L_surf
789 44 : real(dp) :: R_surf
790 44 : real(dp) :: Teff
791 44 : real(dp) :: tau_surf
792 44 : real(dp) :: Teff4
793 44 : real(dp) :: T_surf4
794 44 : real(dp) :: T_surf
795 44 : real(dp) :: P_surf_atm
796 44 : real(dp) :: P_surf
797 44 : real(dp) :: Pextra
798 44 : real(dp) :: kap_surf
799 44 : real(dp) :: M_surf
800 :
801 : include 'formats'
802 :
803 : ! Set up stellar surface parameters
804 :
805 44 : L_surf = s% L(1)
806 44 : R_surf = s% r(1)
807 44 : kap_surf = s% opacity(1)
808 44 : M_surf = s% m(1)
809 44 : Teff = s% Teff
810 :
811 : ! Initialize partials
812 44 : dlnT_dL = 0._dp; dlnT_dlnR = 0._dp; dlnT_dlnM = 0._dp; dlnT_dlnkap = 0._dp
813 44 : dlnP_dL = 0._dp; dlnP_dlnR = 0._dp; dlnP_dlnM = 0._dp; dlnP_dlnkap = 0._dp
814 :
815 : ! Evaluate the surface optical depth
816 :
817 44 : tau_surf = s% tau_factor*s% tau_base ! tau at outer edge of cell 1
818 44 : if (is_bad(tau_surf)) then
819 0 : write(*,1) 's% tau_factor', s% tau_factor
820 0 : write(*,1) 's% tau_base', s% tau_base
821 0 : write(*,1) 'tau_surf', tau_surf
822 0 : call mesa_error(__FILE__,__LINE__,'bad tau_surf in get_surf_PT')
823 : end if
824 :
825 : ! Evaluate surface temperature and pressure
826 :
827 44 : if (.not. (need_atm_Psurf .or. need_atm_Tsurf)) then
828 :
829 : ! Special-case boundary condition
830 :
831 0 : lnP_surf = s% lnPeos_start(1)
832 0 : if (is_bad(lnP_surf)) lnP_surf = 0._dp
833 0 : T_surf = s% T_start(1)
834 0 : lnT_surf = log(T_surf)
835 0 : T_surf4 = T_surf*T_surf*T_surf*T_surf
836 0 : Teff = pow(four_thirds*T_surf4/(tau_surf + two_thirds), 0.25_dp)
837 :
838 0 : if (.not. skip_partials) then
839 0 : dlnT_dL = 0._dp; dlnT_dlnR = 0._dp; dlnT_dlnM = 0._dp; dlnT_dlnkap = 0._dp
840 0 : dlnP_dL = 0._dp; dlnP_dlnR = 0._dp; dlnP_dlnM = 0._dp; dlnP_dlnkap = 0._dp
841 : end if
842 :
843 : else
844 : ! Evaluate temperature and pressure based on atm_option
845 : ! (yes, we might want to translate atm_option into an integer,
846 : ! but these string comparisons are cheap)
847 :
848 : ! The first few are special, 'trivial-atmosphere' options
849 :
850 44 : select case (s% atm_option)
851 :
852 : case ('fixed_Teff')
853 :
854 : ! set Tsurf from Eddington T-tau relation
855 : ! for current surface tau and Teff = `atm_fixed_Teff`.
856 : ! set Psurf = Radiation_Pressure(Tsurf)
857 :
858 0 : Teff = s% atm_fixed_Teff
859 0 : Teff4 = Teff*Teff*Teff*Teff
860 0 : T_surf = pow(3._dp/4._dp*Teff4*(tau_surf + 2._dp/3._dp), 0.25_dp)
861 0 : lnT_surf = log(T_surf)
862 0 : lnP_surf = Radiation_Pressure(T_surf)
863 :
864 0 : if (.not. skip_partials) then
865 0 : dlnT_dL = 0._dp; dlnT_dlnR = 0._dp; dlnT_dlnM = 0._dp; dlnT_dlnkap = 0._dp
866 0 : dlnP_dL = 0._dp; dlnP_dlnR = 0._dp; dlnP_dlnM = 0._dp; dlnP_dlnkap = 0._dp
867 : end if
868 :
869 : case ('fixed_Tsurf')
870 :
871 : ! set Teff from Eddington T-tau relation for given
872 : ! Tsurf and tau=2/3 set Psurf =
873 : ! Radiation_Pressure(Tsurf)
874 :
875 0 : T_surf = s% atm_fixed_Tsurf
876 0 : lnT_surf = log(T_surf)
877 0 : T_surf4 = T_surf*T_surf*T_surf*T_surf
878 0 : Teff = pow(4._dp/3._dp*T_surf4/(tau_surf + 2._dp/3._dp), 0.25_dp)
879 0 : lnP_surf = Radiation_Pressure(T_surf)
880 :
881 0 : if (.not. skip_partials) then
882 0 : dlnT_dL = 0._dp; dlnT_dlnR = 0._dp; dlnT_dlnM = 0._dp; dlnT_dlnkap = 0._dp
883 0 : dlnP_dL = 0._dp; dlnP_dlnR = 0._dp; dlnP_dlnM = 0._dp; dlnP_dlnkap = 0._dp
884 : end if
885 :
886 : case ('fixed_Psurf')
887 :
888 : ! set Teff from L and R using L = 4*pi*R^2*boltz_sigma*T^4.
889 : ! set Tsurf using Eddington T-tau relation
890 :
891 0 : if (L_surf < 0._dp) then
892 0 : if (s% report_ierr) then
893 0 : write(*,2) 'get_surf_PT: L_surf <= 0', s% model_number, L_surf
894 : end if
895 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__)
896 0 : s% retry_message = 'L_surf < 0'
897 0 : ierr = -1
898 0 : return
899 : end if
900 :
901 0 : lnP_surf = safe_log(s% atm_fixed_Psurf)
902 0 : if (R_surf <= 0._dp) then
903 0 : T_surf4 = 1._dp
904 0 : T_surf = 1._dp
905 0 : lnT_surf = 0._dp
906 0 : if (.not. skip_partials) then
907 0 : dlnT_dlnR = 0._dp
908 0 : dlnT_dL = 0._dp
909 : end if
910 0 : Teff = s% T(1)
911 : else
912 0 : Teff = pow(L_surf/(4._dp*pi*R_surf*R_surf*boltz_sigma), 0.25_dp)
913 0 : T_surf4 = 3d0/4d0*pow(Teff, 4.d0)*(tau_surf + 2._dp/3._dp)
914 0 : T_surf = pow(T_surf4, 0.25_dp)
915 0 : lnT_surf = log(T_surf)
916 0 : if (.not. skip_partials) then
917 0 : dlnT_dlnR = -0.5_dp
918 0 : dlnT_dL = 1._dp/(4._dp*L_surf)
919 : end if
920 : end if
921 :
922 0 : if (.not. skip_partials) then
923 0 : dlnT_dlnM = 0._dp; dlnT_dlnkap = 0._dp
924 0 : dlnP_dL = 0._dp; dlnP_dlnR = 0._dp; dlnP_dlnM = 0._dp; dlnP_dlnkap = 0._dp
925 : end if
926 :
927 : case ('fixed_Psurf_and_Tsurf')
928 :
929 0 : lnP_surf = safe_log(s% atm_fixed_Psurf)
930 0 : T_surf = s% atm_fixed_Tsurf
931 0 : lnT_surf = log(T_surf)
932 0 : T_surf4 = T_surf*T_surf*T_surf*T_surf
933 0 : Teff = pow(4d0/3d0*T_surf4/(tau_surf + 2d0/3d0), 0.25d0)
934 :
935 0 : if (.not. skip_partials) then
936 0 : dlnT_dL = 0; dlnT_dlnR = 0; dlnT_dlnM = 0; dlnT_dlnkap = 0
937 0 : dlnP_dL = 0; dlnP_dlnR = 0; dlnP_dlnM = 0; dlnP_dlnkap = 0
938 : end if
939 :
940 : case default
941 :
942 : ! Everything else -- the 'non-trivial atmospheres' ---
943 : ! gets passed to atm_support
944 :
945 : if (.false. .and. 1 == s% solver_test_partials_k .and. &
946 : s% solver_iter == s% solver_test_partials_iter_number) then
947 : star_debugging_atm_flag = .true.
948 : end if
949 :
950 : call get_atm_PT( &
951 : s, tau_surf, L_surf, R_surf, s% m(1), s% cgrav(1), skip_partials, &
952 : Teff, lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
953 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
954 44 : ierr)
955 44 : if (ierr /= 0) then
956 0 : if (s% report_ierr) then
957 0 : write(*,1) 'tau_surf', tau_surf
958 0 : write(*,1) 'L_surf', L_surf
959 0 : write(*,1) 'R_surf', R_surf
960 0 : write(*,1) 'Teff', Teff
961 0 : write(*,1) 's% m(1)', s% m(1)
962 0 : write(*,1) 's% cgrav(1)', s% cgrav(1)
963 0 : write(*,*) 'failed in get_atm_PT'
964 : end if
965 0 : return
966 : end if
967 :
968 : if (.false. .and. 1 == s% solver_test_partials_k .and. &
969 44 : s% solver_iter == s% solver_test_partials_iter_number) then
970 : s% solver_test_partials_val = atm_test_partials_val
971 : s% solver_test_partials_dval_dx = atm_test_partials_dval_dx
972 : end if
973 : end select
974 : end if
975 :
976 : ! if using fixed surface, calculate Pextra.
977 : if (s% atm_option == 'fixed_Tsurf' .or. s% atm_option == 'fixed_Psurf_and_Tsurf' &
978 44 : .or. s% atm_option == 'fixed_Psurf' .or. s% atm_option == 'fixed_Teff') then
979 : ! add extra pressure for fixed atmosphere options
980 0 : if (s% Pextra_factor /= 0._dp) then
981 0 : P_surf_atm = exp(lnP_surf)
982 0 : Pextra = s% Pextra_factor*(kap_surf/tau_surf)*(L_surf/M_surf)/(6._dp*pi*clight*s% cgrav(1))
983 0 : P_surf = P_surf_atm + Pextra
984 0 : if (P_surf < 1E-50_dp) then
985 0 : lnP_surf = -50*ln10
986 0 : if (.not. skip_partials) then
987 0 : dlnP_dL = 0._dp
988 0 : dlnP_dlnR = 0._dp
989 0 : dlnP_dlnM = 0._dp
990 0 : dlnP_dlnkap = 0._dp
991 : end if
992 : else
993 0 : lnP_surf = log(P_surf)
994 0 : if (.not. skip_partials) then
995 0 : dlnP_dL = dlnP_dL*P_surf_atm/P_surf
996 0 : dlnP_dlnR = dlnP_dlnR*P_surf_atm/P_surf
997 0 : dlnP_dlnM = dlnP_dlnM*P_surf_atm/P_surf
998 0 : dlnP_dlnkap = dlnP_dlnkap*P_surf_atm/P_surf
999 : end if
1000 : end if
1001 : end if
1002 : end if
1003 :
1004 : ! Check outputs
1005 :
1006 44 : if (is_bad(lnT_surf) .or. is_bad(lnP_surf)) then
1007 0 : if (len_trim(s% retry_message) == 0) s% retry_message = 'bad logT surf or logP surf'
1008 0 : ierr = -1
1009 0 : write(*,1) 'bad outputs in get_surf_PT'
1010 0 : write(*,1) 'lnT_surf', lnT_surf
1011 0 : write(*,1) 'lnP_surf', lnP_surf
1012 0 : write(*,*) 'atm_option = ', trim(s% atm_option)
1013 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'get PT surf')
1014 : end if
1015 :
1016 : return
1017 :
1018 44 : end subroutine get_surf_PT
1019 :
1020 :
1021 44 : subroutine set_grads(s, ierr)
1022 44 : use chem_def, only: chem_isos
1023 : use star_utils, only: smooth, safe_div_val
1024 : type (star_info), pointer :: s
1025 : integer, intent(out) :: ierr
1026 :
1027 : integer :: k, nz, j, cid, max_cid
1028 44 : real(dp) :: val, max_val, A, Z
1029 44 : real(dp), pointer, dimension(:) :: dlnP, dlnd, dlnT
1030 :
1031 : include 'formats'
1032 :
1033 44 : ierr = 0
1034 44 : nz = s% nz
1035 44 : call do_alloc(ierr)
1036 44 : if (ierr /= 0) return
1037 :
1038 53820 : do k = 2, nz
1039 53776 : dlnP(k) = s% lnPeos(k-1) - s% lnPeos(k)
1040 53776 : dlnd(k) = s% lnd(k-1) - s% lnd(k)
1041 53820 : dlnT(k) = s% lnT(k-1) - s% lnT(k)
1042 : end do
1043 44 : dlnP(1) = dlnP(2)
1044 44 : dlnd(1) = dlnd(2)
1045 44 : dlnT(1) = dlnT(2)
1046 :
1047 44 : call smooth(dlnP,nz)
1048 44 : call smooth(dlnd,nz)
1049 44 : call smooth(dlnT,nz)
1050 :
1051 44 : s% grad_density(1) = 0
1052 44 : s% grad_temperature(1) = 0
1053 53820 : do k = 2, nz
1054 53820 : if (dlnP(k) >= 0) then
1055 2 : s% grad_density(k) = 0
1056 2 : s% grad_temperature(k) = 0
1057 : else
1058 53774 : s% grad_density(k) = safe_div_val(s, dlnd(k), dlnP(k))
1059 53774 : s% grad_temperature(k) = safe_div_val(s, dlnT(k), dlnP(k))
1060 : end if
1061 : end do
1062 :
1063 44 : call smooth(s% grad_density,nz)
1064 44 : call smooth(s% grad_temperature,nz)
1065 :
1066 : ! this will compute the values of s% smoothed_brunt_B
1067 44 : if (s% calculate_Brunt_B) then
1068 53864 : do k=1,nz
1069 53864 : s% smoothed_brunt_B(k) = s% unsmoothed_brunt_B(k)
1070 : end do
1071 44 : call compute_smoothed_brunt_B
1072 : else
1073 0 : s% smoothed_brunt_B(:) = 0d0
1074 : end if
1075 44 : if (s% use_Ledoux_criterion .and. s% calculate_Brunt_B) then
1076 0 : do k=1,nz
1077 0 : s% gradL_composition_term(k) = s% smoothed_brunt_B(k)
1078 : end do
1079 : else
1080 53864 : do k=1,nz
1081 53864 : s% gradL_composition_term(k) = 0d0
1082 : end do
1083 : end if
1084 :
1085 : call dealloc
1086 :
1087 53688 : do k=3,nz-2
1088 53644 : max_cid = 0
1089 53644 : max_val = -1d99
1090 482796 : do j=1,s% species
1091 429152 : cid = s% chem_id(j)
1092 429152 : A = dble(chem_isos% Z_plus_N(cid))
1093 429152 : Z = dble(chem_isos% Z(cid))
1094 429152 : val = (s% xa(j,k-2) + s% xa(j,k-1) - s% xa(j,k) - s% xa(j,k+1))*(1d0 + Z)/A
1095 482796 : if (val > max_val) then
1096 66337 : max_val = val
1097 66337 : max_cid = cid
1098 : end if
1099 : end do
1100 53688 : s% dominant_iso_for_thermohaline(k) = max_cid
1101 : end do
1102 : s% dominant_iso_for_thermohaline(1:2) = &
1103 132 : s% dominant_iso_for_thermohaline(3)
1104 : s% dominant_iso_for_thermohaline(nz-1:nz) = &
1105 132 : s% dominant_iso_for_thermohaline(nz-2)
1106 :
1107 :
1108 : contains
1109 :
1110 132 : subroutine compute_smoothed_brunt_B
1111 44 : use star_utils, only: weighed_smoothing, threshold_smoothing
1112 : logical, parameter :: preserve_sign = .false.
1113 : real(dp), pointer, dimension(:) :: work
1114 : include 'formats'
1115 44 : ierr = 0
1116 44 : work => dlnd
1117 44 : if (s% num_cells_for_smooth_gradL_composition_term <= 0) return
1118 : call threshold_smoothing( &
1119 : s% smoothed_brunt_B, s% threshold_for_smooth_gradL_composition_term, s% nz, &
1120 44 : s% num_cells_for_smooth_gradL_composition_term, preserve_sign, work)
1121 44 : end subroutine compute_smoothed_brunt_B
1122 :
1123 44 : subroutine do_alloc(ierr)
1124 : integer, intent(out) :: ierr
1125 44 : call do_work_arrays(.true.,ierr)
1126 44 : end subroutine do_alloc
1127 :
1128 44 : subroutine dealloc
1129 44 : call do_work_arrays(.false.,ierr)
1130 : end subroutine dealloc
1131 :
1132 88 : subroutine do_work_arrays(alloc_flag, ierr)
1133 : use alloc, only: work_array
1134 : logical, intent(in) :: alloc_flag
1135 : integer, intent(out) :: ierr
1136 : logical, parameter :: crit = .false.
1137 : ierr = 0
1138 : call work_array(s, alloc_flag, crit, &
1139 88 : dlnP, nz, nz_alloc_extra, 'mlt', ierr)
1140 88 : if (ierr /= 0) return
1141 : call work_array(s, alloc_flag, crit, &
1142 88 : dlnd, nz, nz_alloc_extra, 'mlt', ierr)
1143 88 : if (ierr /= 0) return
1144 : call work_array(s, alloc_flag, crit, &
1145 88 : dlnT, nz, nz_alloc_extra, 'mlt', ierr)
1146 88 : if (ierr /= 0) return
1147 88 : end subroutine do_work_arrays
1148 :
1149 : end subroutine set_grads
1150 :
1151 : end module hydro_vars
|
