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