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