Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2018-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_temperature
21 :
22 : use star_private_def
23 : use const_def, only: dp, ln10, pi4, crad, clight, convective_mixing
24 : use reconstructed_face_support, only: get_reconstructed_face_eos_kap_ad
25 : use utils_lib, only: mesa_error, is_bad
26 : use auto_diff
27 : use auto_diff_support
28 :
29 : implicit none
30 :
31 : private
32 : public :: do1_alt_dlnT_dm_eqn
33 : public :: do1_gradT_eqn
34 : public :: do1_dlnT_dm_eqn
35 :
36 : contains
37 :
38 : ! just relate L_rad to T gradient.
39 : ! d_P_rad/dm = -<opacity_face>*L_rad/(clight*area^2) -- see, e.g., K&W (5.12)
40 : ! P_rad = (1/3)*crad*T^4
41 : ! d_P_rad/dm = (crad/3)*(T(k-1)^4 - T(k)^4)/dm_bar
42 : ! L_rad = L - L_non_rad, L_non_rad = L_start - L_rad_start
43 : ! L_rad_start = (-d_P_rad/dm_bar*clight*area^2/<opacity_face>)_start
44 0 : subroutine do1_alt_dlnT_dm_eqn(s, k, nvar, ierr)
45 : use eos_def
46 : use star_utils, only: save_eqn_residual_info, get_face_weights
47 : type (star_info), pointer :: s
48 : integer, intent(in) :: k, nvar
49 : integer, intent(out) :: ierr
50 :
51 : real(dp) :: alfa, beta, scale, dm_bar
52 : type(auto_diff_real_star_order1) :: L_ad, r_00, area, area2, Lrad_ad, &
53 : kap_00, kap_m1, kap_face, d_P_rad_expected_ad, T_m1, T4_m1, T_00, T4_00, &
54 : P_rad_m1, P_rad_00, d_P_rad_actual_ad, resid
55 : type(auto_diff_real_star_order1) :: T_face, rho_face, P_face, Cp_face, ChiRho_face, ChiT_face, grada_face
56 : type(auto_diff_real_star_order1) :: flxR, flxLambda
57 :
58 : integer :: i_equL
59 : logical :: dbg
60 : logical :: test_partials
61 :
62 : include 'formats'
63 0 : ierr = 0
64 0 : i_equL = s% i_equL
65 0 : if (i_equL == 0) return
66 :
67 0 : if (.not. s% use_dPrad_dm_form_of_T_gradient_eqn) then
68 0 : ierr = -1
69 0 : return
70 : end if
71 :
72 : !test_partials = (k == s% solver_test_partials_k)
73 0 : test_partials = .false.
74 :
75 0 : dbg = .false.
76 :
77 0 : call get_face_weights(s, k, alfa, beta)
78 :
79 0 : scale = s% energy_start(k)*s% rho_start(k)
80 0 : dm_bar = s% dm_bar(k)
81 0 : L_ad = wrap_L_00(s,k)
82 0 : r_00 = wrap_r_00(s,k)
83 0 : area = pi4*pow2(r_00); area2 = pow2(area)
84 :
85 : if (s% lnT(k)/ln10 <= s% max_logT_for_mlt &
86 : .and. s% mixing_type(k) == convective_mixing .and. s% gradr(k) > 0d0 &
87 0 : .and. abs(s% gradr(k) - s% gradT(k)) > abs(s% gradr(k))*1d-5) then
88 0 : Lrad_ad = L_ad*s% gradT_ad(k)/s% gradr_ad(k) ! C&G 14.109
89 : else
90 0 : Lrad_ad = L_ad
91 : end if
92 :
93 0 : if (s% use_face_reconstruction) then
94 0 : if (s% reconstructed_face_state_valid(k)) then
95 0 : kap_face = s% reconstructed_opacity_face_ad(k)
96 : else
97 : call get_reconstructed_face_eos_kap_ad( &
98 0 : s, k, T_face, rho_face, P_face, Cp_face, ChiRho_face, ChiT_face, grada_face, kap_face, ierr)
99 0 : if (ierr /= 0) return
100 : end if
101 : else
102 0 : kap_00 = wrap_kap_00(s,k)
103 0 : kap_m1 = wrap_kap_m1(s,k)
104 0 : kap_face = alfa*kap_00 + beta*kap_m1
105 : end if
106 0 : if (kap_face%val < s% min_kap_for_dPrad_dm_eqn) &
107 0 : kap_face = s% min_kap_for_dPrad_dm_eqn
108 :
109 : ! calculate expected d_P_rad from current L_rad
110 0 : d_P_rad_expected_ad = -dm_bar*kap_face*Lrad_ad/(clight*area2)
111 :
112 : ! calculate actual d_P_rad in current model
113 0 : T_m1 = wrap_T_m1(s,k); T4_m1 = pow4(T_m1)
114 0 : T_00 = wrap_T_00(s,k); T4_00 = pow4(T_00)
115 :
116 : !d_P_rad_expected = d_P_rad_expected*s% gradr_factor(k) !TODO(Pablo): check this
117 :
118 0 : P_rad_m1 = (crad/3._dp)*T4_m1
119 0 : P_rad_00 = (crad/3._dp)*T4_00
120 0 : d_P_rad_actual_ad = P_rad_m1 - P_rad_00
121 :
122 : ! enable flux-limited radiation transport derived by Levermore & Pomraning 1981
123 0 : s% flux_limit_R(k) = 0._dp
124 0 : s% flux_limit_lambda(k) =0._dp
125 0 : if (s% use_flux_limiting_with_dPrad_dm_form) then
126 : ! calculate the flux ratio R
127 : flxR = area * abs(T4_m1 - T4_00) / dm_bar / &
128 0 : (kap_face * 0.5_dp * (T4_m1 + T4_00))
129 :
130 0 : s% flux_limit_R(k) = flxR%val
131 :
132 : ! calculate the flux limiter lambda
133 0 : flxLambda = (6._dp + 3._dp*flxR) / (6._dp + (3._dp + flxR)*flxR)
134 :
135 0 : s% flux_limit_lambda(k) = flxLambda%val
136 :
137 : ! calculate d_P_rad given the flux limiter
138 0 : d_P_rad_expected_ad = d_P_rad_expected_ad / flxLambda
139 : end if
140 :
141 : ! residual
142 0 : resid = (d_P_rad_expected_ad - d_P_rad_actual_ad)/scale
143 0 : s% equ(i_equL, k) = resid%val
144 :
145 0 : if (is_bad(resid%val)) then
146 0 : !$OMP critical (star_alt_dlntdm_bad_num)
147 0 : write(*,2) 'resid%val', k, resid%val
148 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'do1_alt_dlnT_dm_eqn')
149 : !$OMP end critical (star_alt_dlntdm_bad_num)
150 : end if
151 :
152 : if (test_partials) then
153 : s% solver_test_partials_val = s% gradT(k)
154 : end if
155 :
156 : call save_eqn_residual_info( &
157 0 : s, k, nvar, i_equL, resid, 'do1_alt_dlnT_dm_eqn', ierr)
158 :
159 : if (test_partials) then
160 : s% solver_test_partials_var = 0
161 : s% solver_test_partials_dval_dx = 0
162 : write(*,*) 'do1_alt_dlnT_dm_eqn', s% solver_test_partials_var
163 : end if
164 :
165 : contains
166 :
167 : end subroutine do1_alt_dlnT_dm_eqn
168 :
169 :
170 0 : subroutine do1_gradT_eqn(s, k, nvar, ierr)
171 : use eos_def
172 : use star_utils, only: save_eqn_residual_info
173 : type (star_info), pointer :: s
174 : integer, intent(in) :: k, nvar
175 : integer, intent(out) :: ierr
176 :
177 : type(auto_diff_real_star_order1) :: &
178 : resid, gradT, dlnT, dlnP
179 : integer :: i_equL
180 : logical :: test_partials
181 :
182 : include 'formats'
183 0 : ierr = 0
184 :
185 : !test_partials = (k == s% solver_test_partials_k)
186 0 : test_partials = .false.
187 :
188 0 : i_equL = s% i_equL
189 0 : if (i_equL == 0) return
190 :
191 0 : gradT = s% gradT_ad(k)
192 0 : dlnT = wrap_lnT_m1(s,k) - wrap_lnT_00(s,k)
193 0 : dlnP = wrap_lnPeos_m1(s,k) - wrap_lnPeos_00(s,k)
194 :
195 0 : resid = gradT*dlnP - dlnT
196 0 : s% equ(i_equL, k) = resid%val
197 :
198 0 : if (is_bad(s% equ(i_equL, k))) then
199 0 : ierr = -1
200 0 : if (s% report_ierr) write(*,2) 'equ(i_equL, k)', k, s% equ(i_equL, k)
201 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'do1_gradT_eqn')
202 : return
203 : write(*,2) 'equ(i_equL, k)', k, s% equ(i_equL, k)
204 : write(*,2) 'gradT', k, gradT
205 : write(*,2) 'dlnT', k, dlnT
206 : write(*,2) 'dlnP', k, dlnP
207 : call mesa_error(__FILE__,__LINE__,'do1_gradT_eqn')
208 : end if
209 :
210 : if (test_partials) then
211 : s% solver_test_partials_val = s% equ(i_equL,k)
212 : end if
213 :
214 : call save_eqn_residual_info( &
215 0 : s, k, nvar, i_equL, resid, 'do1_gradT_eqn', ierr)
216 :
217 : !call set_xtras
218 :
219 : contains
220 :
221 : subroutine set_xtras
222 : use auto_diff_support
223 : use star_utils, only: get_Lrad
224 : type(auto_diff_real_star_order1) :: &
225 : T4m1, T400, kap_m1, kap_00, alfa, beta, kap_face, &
226 : diff_T4_div_kap
227 : T4m1 = pow4(wrap_T_m1(s,k))
228 : T400 = pow4(wrap_T_00(s,k))
229 : kap_m1 = wrap_kap_m1(s,k)
230 : kap_00 = wrap_kap_00(s,k)
231 : alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
232 : beta = 1d0 - alfa
233 : kap_face = alfa*kap_00 + beta*kap_m1
234 : diff_T4_div_kap = (T4m1 - T400)/kap_face
235 : s% xtra1_array(k) = s% T_start(k)
236 : s% xtra2_array(k) = T4m1%val - T400%val
237 : s% xtra3_array(k) = kap_face%val
238 : s% xtra4_array(k) = diff_T4_div_kap%val
239 : s% xtra5_array(k) = get_Lrad(s,k)
240 : s% xtra6_array(k) = 1
241 : end subroutine set_xtras
242 :
243 : end subroutine do1_gradT_eqn
244 :
245 :
246 209216 : subroutine do1_dlnT_dm_eqn(s, k, nvar, ierr)
247 : use eos_def
248 : use star_utils, only: save_eqn_residual_info
249 : type (star_info), pointer :: s
250 : integer, intent(in) :: k, nvar
251 : integer, intent(out) :: ierr
252 :
253 : type(auto_diff_real_star_order1) :: resid, &
254 : dlnPdm, Ppoint, gradT, dlnTdm, T00, Tm1, dT, Tpoint, lnTdiff
255 : real(dp) :: delm, alfa
256 : integer :: i_equL
257 : logical :: test_partials
258 :
259 : include 'formats'
260 52304 : ierr = 0
261 :
262 : !test_partials = (k == s% solver_test_partials_k)
263 52304 : test_partials = .false.
264 :
265 52304 : i_equL = s% i_equL
266 52304 : if (i_equL == 0) return
267 :
268 52304 : if (k ==1 .and. s% use_RSP_L_eqn_outer_BC) then
269 0 : call set_RSP_Lsurf_BC(s, nvar, ierr)
270 0 : return
271 : end if
272 :
273 52304 : if (s% use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn) then
274 0 : call do1_gradT_eqn(s, k, nvar, ierr)
275 0 : return
276 : end if
277 :
278 52304 : if (s% use_dPrad_dm_form_of_T_gradient_eqn) then
279 0 : call do1_alt_dlnT_dm_eqn(s, k, nvar, ierr)
280 0 : return
281 : end if
282 :
283 : ! dT/dm = dP/dm * T/P * grad_T, grad_T = dlnT/dlnP from MLT.
284 : ! but use hydrostatic value for dP/dm in this.
285 : ! this is because of limitations of MLT for calculating grad_T.
286 : ! (MLT assumes hydrostatic equilibrium)
287 : ! see comment in K&W chpt 9.1.
288 :
289 52304 : call eval_dlnPdm_qhse(s, k, dlnPdm, Ppoint, ierr)
290 52304 : if (ierr /= 0) return
291 :
292 52304 : gradT = s% gradT_ad(k)
293 52304 : dlnTdm = dlnPdm*gradT
294 :
295 52304 : Tm1 = wrap_T_m1(s,k)
296 52304 : T00 = wrap_T_00(s,k)
297 52304 : dT = Tm1 - T00
298 52304 : alfa = s% dm(k-1)/(s% dm(k-1) + s% dm(k))
299 52304 : Tpoint = alfa*T00 + (1d0 - alfa)*Tm1
300 52304 : lnTdiff = dT/Tpoint ! use this in place of lnT(k-1)-lnT(k)
301 52304 : delm = (s% dm(k) + s% dm(k-1))/2
302 :
303 52304 : resid = delm*dlnTdm - lnTdiff
304 52304 : s% equ(i_equL, k) = resid%val
305 :
306 52304 : if (is_bad(s% equ(i_equL, k))) then
307 0 : ierr = -1
308 0 : if (s% report_ierr) write(*,2) 'equ(i_equL, k)', k, s% equ(i_equL, k)
309 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'hydro eqns')
310 : return
311 : write(*,2) 'equ(i_equL, k)', k, s% equ(i_equL, k)
312 : write(*,2) 'lnTdiff', k, lnTdiff
313 : write(*,2) 'delm', k, delm
314 : write(*,2) 'dlnPdm', k, dlnPdm
315 : write(*,2) 'gradT', k, gradT
316 : call mesa_error(__FILE__,__LINE__,'i_equL')
317 : end if
318 :
319 : if (test_partials) then
320 : s% solver_test_partials_val = s% equ(i_equL,k)
321 : end if
322 :
323 : call save_eqn_residual_info( &
324 52304 : s, k, nvar, i_equL, resid, 'do1_dlnT_dm_eqn', ierr)
325 :
326 : end subroutine do1_dlnT_dm_eqn
327 :
328 :
329 :
330 0 : subroutine set_RSP_Lsurf_BC(s, nvar, ierr)
331 : use const_def, only: crad, clight, pi4
332 : use eos_def
333 : use star_utils, only: save_eqn_residual_info, get_area_info_opt_time_center
334 : use auto_diff_support
335 : implicit none
336 :
337 : type(star_info), pointer :: s
338 : integer, intent(out) :: ierr
339 : integer, intent(in) :: nvar
340 :
341 : type(auto_diff_real_star_order1) :: L1_ad, r1_ad, area_ad, rhs_ad, lhs_ad, resid_ad, inv_R2
342 : type(auto_diff_real_star_order1) :: T_surf, Erad_ad
343 : integer :: i_equL
344 : real(dp) :: factor, scale, L_theta
345 : logical :: debug
346 :
347 0 : ierr = 0
348 0 : debug = .false.
349 :
350 0 : i_equL = s% i_equL
351 :
352 0 : if (s%nz < 1) then
353 0 : write(*,*) 'ERROR: Insufficient zones (nz < 1)'
354 0 : ierr = -1
355 0 : return
356 : end if
357 :
358 : if (debug) write(*,*) 'RSP zone 1 surface BC being set'
359 :
360 0 : call get_area_info_opt_time_center(s, 1, area_ad, inv_R2, ierr)
361 : ! no time centering the surface equations.
362 0 : L1_ad = wrap_L_00(s, 1)
363 0 : T_surf = wrap_T_00(s,1)
364 :
365 : if (debug) then
366 : write(*,*) 'T_surf =', T_surf%val, ' r_surf =', r1_ad%val, ' area =', area_ad%val
367 : end if
368 :
369 : ! rsp equation, zone 1
370 0 : rhs_ad = s%RSP2_Lsurf_factor * area_ad * clight * (crad * pow4(T_surf)) ! missing Lc at the moment, so only radiative surface
371 :
372 : if (debug) then
373 : write(*,*) 'RSP_Lsurf_factor =', s%RSP2_Lsurf_factor
374 : write(*,*) 'rhs_ad (RSP BC) =', rhs_ad%val
375 : end if
376 :
377 : ! residual
378 0 : lhs_ad = L1_ad
379 0 : resid_ad = lhs_ad - rhs_ad
380 :
381 0 : scale =maxval(s% L_start(1:s% nz))
382 0 : resid_ad = resid_ad / scale
383 :
384 : if (debug) then
385 : write(*,*) 'lhs (L1) =', lhs_ad%val
386 : write(*,*) 'scaled residual =', resid_ad%val
387 : end if
388 :
389 0 : s%equ(i_equL,1) = resid_ad%val
390 :
391 0 : if (is_bad(resid_ad%val)) then
392 0 : write(*,*) 'ERROR: NaN or Inf residual:', resid_ad%val
393 : ierr = -1
394 : end if
395 :
396 : call save_eqn_residual_info( &
397 0 : s, 1, nvar, i_equL, resid_ad, 'do1_dlnT_dm_eqn', ierr)
398 :
399 :
400 : end subroutine set_RSP_Lsurf_BC
401 :
402 : ! only used for dlnT_dm equation
403 52304 : subroutine eval_dlnPdm_qhse(s, k, & ! calculate the expected dlnPdm for HSE
404 : dlnPdm_qhse, Ppoint, ierr)
405 : use hydro_momentum, only: expected_HSE_grav_term
406 : type (star_info), pointer :: s
407 : integer, intent(in) :: k
408 : type(auto_diff_real_star_order1), intent(out) :: dlnPdm_qhse, Ppoint
409 : integer, intent(out) :: ierr
410 :
411 : real(dp) :: alfa, P_theta
412 : type(auto_diff_real_star_order1) :: grav, area, P00, Pm1, inv_R2, mlt_Ptrb00, mlt_Ptrbm1, mlt_Ptrb_face
413 : type(auto_diff_real_star_order1) :: T_face, rho_face, P_face, Cp_face, ChiRho_face, ChiT_face, grada_face, opacity_face
414 : include 'formats'
415 :
416 : ierr = 0
417 :
418 : ! basic eqn is dP/dm = -G m / (4 pi r^4)
419 : ! divide by Ppoint to make it unitless
420 :
421 : ! for rotation, multiply gravity by factor fp. MESA 2, eqn 22.
422 52304 : call expected_HSE_grav_term(s, k, grav, area, ierr) ! note that expected_HSE_grav_term is negative
423 :
424 : if (s% using_velocity_time_centering .and. &
425 52304 : s% include_P_in_velocity_time_centering .and. &
426 : s% lnT(k)/ln10 <= s% max_logT_for_include_P_and_L_in_velocity_time_centering) then
427 0 : P_theta = s% P_theta_for_velocity_time_centering
428 : else
429 52304 : P_theta = 1d0
430 : end if
431 :
432 52304 : if (s% use_face_reconstruction) then
433 0 : if (s% reconstructed_face_state_valid(k)) then
434 0 : rho_face = s% reconstructed_rho_face_ad(k)
435 0 : Ppoint = s% reconstructed_P_face_ad(k)
436 : else
437 : call get_reconstructed_face_eos_kap_ad( &
438 0 : s, k, T_face, rho_face, P_face, Cp_face, ChiRho_face, ChiT_face, grada_face, opacity_face, ierr)
439 0 : if (ierr /= 0) return
440 0 : Ppoint = P_face
441 : end if
442 0 : if (P_theta /= 1d0) then
443 0 : Ppoint = P_theta*Ppoint + (1d0 - P_theta)*s% reconstructed_P_face_start(k)
444 : end if
445 : if (s% have_mlt_vc .and. s% okay_to_set_mlt_vc .and. s% include_mlt_Pturb_in_thermodynamic_gradients &
446 0 : .and. s% mlt_Pturb_factor > 0d0) then
447 : ! Keep the lagged convective velocity, but form the pressure term from the same
448 : ! face density used by the reconstructed face thermodynamic quantities.
449 0 : mlt_Ptrb_face = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*rho_face/3d0
450 0 : Ppoint = Ppoint + mlt_Ptrb_face
451 : end if
452 : else
453 : ! mlt_pturb in thermodynamic gradients does not currently support time centering because it is timelagged.
454 : ! replace mlt_vc check with s% mlt_vc_old(k) >0 check.
455 : if ((s% have_mlt_vc .and. s% okay_to_set_mlt_vc) .and. s% include_mlt_Pturb_in_thermodynamic_gradients &
456 52304 : .and. s% mlt_Pturb_factor > 0d0) then
457 0 : if (k ==1) then
458 0 : mlt_Ptrb00 = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*wrap_d_00(s,k)/3d0
459 0 : mlt_Ptrbm1 = 0d0
460 : else
461 0 : mlt_Ptrb00 = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*wrap_d_00(s,k)/3d0
462 0 : mlt_Ptrbm1 = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*wrap_d_m1(s,k)/3d0
463 : end if
464 : else ! no mlt_pturb
465 52304 : mlt_Ptrb00 = 0d0
466 52304 : mlt_Ptrbm1 = 0d0
467 : end if
468 :
469 52304 : P00 = wrap_Peos_00(s,k)
470 :
471 : ! mlt Pturb doesn't support time centering yet.
472 52304 : if (P_theta /= 1d0) P00 = P_theta*P00 + (1d0 - P_theta)*s% Peos_start(k)
473 :
474 52304 : if (k == 1) then
475 0 : Pm1 = 0d0
476 0 : Ppoint = P00 + mlt_Ptrb00
477 : else
478 52304 : Pm1 = wrap_Peos_m1(s,k)
479 52304 : if (P_theta /= 1d0) Pm1 = P_theta*Pm1 + (1d0 - P_theta)*s% Peos_start(k-1)
480 52304 : Pm1 = Pm1 + mlt_Ptrbm1 ! include mlt Ptrb in k-1
481 52304 : P00 = P00 + mlt_Ptrb00 ! include mlt Ptrb in k
482 52304 : alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
483 52304 : Ppoint = alfa*P00 + (1d0-alfa)*Pm1
484 : end if
485 : end if
486 :
487 52304 : dlnPdm_qhse = grav/(area*Ppoint) ! note that expected_HSE_grav_term is negative
488 :
489 52304 : if (is_bad(dlnPdm_qhse%val)) then
490 0 : ierr = -1
491 0 : s% retry_message = 'eval_dlnPdm_qhse: is_bad(dlnPdm_qhse)'
492 0 : if (s% report_ierr) then
493 0 : !$OMP critical (hydro_vars_crit1)
494 0 : write(*,*) 'eval_dlnPdm_qhse: is_bad(dlnPdm_qhse)'
495 0 : stop
496 : !$OMP end critical (hydro_vars_crit1)
497 : end if
498 0 : if (s% stop_for_bad_nums) then
499 0 : write(*,2) 'dlnPdm_qhse', k, dlnPdm_qhse
500 0 : call mesa_error(__FILE__,__LINE__,'eval_dlnPdm_qhse')
501 : end if
502 : return
503 : end if
504 :
505 : end subroutine eval_dlnPdm_qhse
506 :
507 : end module hydro_temperature
|
