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