Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2018-2019 Radek Smolec & 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 rsp_step
21 : use rsp_def
22 : use star_def, only: star_ptr, star_info
23 : use const_def, only: dp, qp, i8, crad
24 : use utils_lib, only: is_bad
25 :
26 : implicit none
27 :
28 : private
29 : public :: calculate_energies
30 : public :: init_HYD
31 : public :: HYD
32 : public :: turn_off_time_weighting
33 : public :: turn_on_time_weighting
34 : public :: eval_vars
35 : public :: eval_eqns
36 : public :: calc_equations
37 : public :: save_start_vars
38 : public :: do1_eos_and_kap
39 : public :: calc_Fr
40 : public :: do1_specific_volume
41 : public :: get_Psurf
42 : public :: set_f_Edd
43 : public :: calc_Prad
44 : public :: calc_Hp_face
45 : public :: calc_Y_face
46 : public :: calc_PII_face
47 : public :: calc_Pvsc
48 : public :: calc_Pturb
49 : public :: calc_Chi
50 : public :: calc_Eq
51 : public :: calc_source_sink
52 : public :: acceleration_eqn
53 : public :: calc_cell_equations
54 : public :: T_form_of_calc_Fr
55 : public :: calc_Lc, calc_Lt
56 : public :: rsp_set_Teff
57 :
58 : logical, parameter :: call_is_bad = .false.
59 : integer, parameter :: i_var_Vol = 99 ! for remeshing tests with dfridr
60 : integer, parameter :: &
61 : i_var_T = 2, i_var_w = 3, i_var_er = 4, i_var_Fr = 5, i_var_R = 6, & ! R must be last
62 :
63 : i_r_dr_in2 = 1, &
64 : i_r_dT_in = 2, &
65 : i_r_dw_in = 3, &
66 : i_r_der_in = 4, &
67 : i_r_dFr_in = 5, &
68 :
69 : i_r_dr_in = i_r_dr_in2 + NV, &
70 : i_r_dr_00 = i_r_dr_in + NV, &
71 : i_r_dr_out = i_r_dr_00 + NV, &
72 : i_r_dr_out2 = i_r_dr_out + NV, &
73 :
74 : i_r_dT_00 = i_r_dT_in + NV, &
75 : i_r_dT_out = i_r_dT_00 + NV, &
76 : i_r_dT_out2 = i_r_dT_out + NV, &
77 :
78 : i_r_dw_00 = i_r_dw_in + NV, &
79 : i_r_dw_out = i_r_dw_00 + NV, &
80 : i_r_dw_out2 = i_r_dw_out + NV, &
81 :
82 : i_r_der_00 = i_r_der_in + NV, &
83 : i_r_der_out = i_r_der_00 + NV, &
84 : i_r_der_out2 = i_r_der_out + NV, &
85 :
86 : i_r_dFr_00 = i_r_dFr_in + NV, &
87 : i_r_dFr_out = i_r_dFr_00 + NV, &
88 : i_r_dFr_out2 = i_r_dFr_out + NV, &
89 :
90 : i_T_dT_in2 = 1, &
91 : i_T_dw_in2 = 2, &
92 : i_T_der_in2 = 3, &
93 : i_T_dFr_in2 = 4, &
94 : i_T_dr_in2 = 5, &
95 :
96 : i_T_dT_in = i_T_dT_in2 + NV, &
97 : i_T_dT_00 = i_T_dT_in + NV, &
98 : i_T_dT_out = i_T_dT_00 + NV, &
99 : i_T_dT_out2 = i_T_dT_out + NV, &
100 :
101 : i_T_dw_in = i_T_dw_in2 + NV, &
102 : i_T_dw_00 = i_T_dw_in + NV, &
103 : i_T_dw_out = i_T_dw_00 + NV, &
104 :
105 : i_T_der_in = i_T_der_in2 + NV, &
106 : i_T_der_00 = i_T_der_in + NV, &
107 : i_T_der_out = i_T_der_00 + NV, &
108 :
109 : i_T_dFr_in = i_T_dFr_in2 + NV, &
110 : i_T_dFr_00 = i_T_dFr_in + NV, &
111 : i_T_dFr_out = i_T_dFr_00 + NV, &
112 :
113 : i_T_dr_in = i_T_dr_in2 + NV, &
114 : i_T_dr_00 = i_T_dr_in + NV, &
115 : i_T_dr_out = i_T_dr_00 + NV, &
116 :
117 : i_w_dw_in2 = 1, &
118 : i_w_der_in2 = 2, &
119 : i_w_dFr_in2 = 3, &
120 : i_w_dr_in2 = 4, &
121 : i_w_dT_in = 5, &
122 :
123 : i_w_dw_in = i_w_dw_in2 + NV, &
124 : i_w_dw_00 = i_w_dw_in + NV, &
125 : i_w_dw_out = i_w_dw_00 + NV, &
126 : i_w_dw_out2 = i_w_dw_out + NV, &
127 :
128 : i_w_der_in = i_w_der_in2 + NV, &
129 : i_w_der_00 = i_w_der_in + NV, &
130 : i_w_der_out = i_w_der_00 + NV, &
131 :
132 : i_w_dFr_in = i_w_dFr_in2 + NV, &
133 : i_w_dFr_00 = i_w_dFr_in + NV, &
134 : i_w_dFr_out = i_w_dFr_00 + NV, &
135 :
136 : i_w_dr_in = i_w_dr_in2 + NV, &
137 : i_w_dr_00 = i_w_dr_in + NV, &
138 : i_w_dr_out = i_w_dr_00 + NV, &
139 :
140 : i_w_dT_00 = i_w_dT_in + NV, &
141 : i_w_dT_out = i_w_dT_00 + NV, &
142 : i_w_dT_out2 = i_w_dT_out + NV, &
143 :
144 : i_er_der_in2 = 1, &
145 : i_er_dFr_in2 = 2, &
146 : i_er_dr_in2 = 3, &
147 : i_er_dT_in = 4, &
148 : i_er_dw_in = 5, &
149 :
150 : i_er_der_in = i_er_der_in2 + NV, &
151 : i_er_der_00 = i_er_der_in + NV, &
152 : i_er_der_out = i_er_der_00 + NV, &
153 : i_er_der_out2 = i_er_der_out + NV, &
154 :
155 : i_er_dFr_in = i_er_dFr_in2 + NV, &
156 : i_er_dFr_00 = i_er_dFr_in + NV, &
157 : i_er_dFr_out = i_er_dFr_00 + NV, &
158 :
159 : i_er_dr_in = i_er_dr_in2 + NV, &
160 : i_er_dr_00 = i_er_dr_in + NV, &
161 : i_er_dr_out = i_er_dr_00 + NV, &
162 :
163 : i_er_dT_00 = i_er_dT_in + NV, &
164 : i_er_dT_out = i_er_dT_00 + NV, &
165 : i_er_dT_out2 = i_er_dT_out + NV, &
166 :
167 : i_er_dw_00 = i_er_dw_in + NV, &
168 : i_er_dw_out = i_er_dw_00 + NV, &
169 : i_er_dw_out2 = i_er_dw_out + NV, &
170 :
171 : i_Fr_dFr_in2 = 1, &
172 : i_Fr_dr_in2 = 2, &
173 : i_Fr_dT_in = 3, &
174 : i_Fr_dw_in = 4, &
175 : i_Fr_der_in = 5, &
176 :
177 : i_Fr_dFr_in = i_Fr_dFr_in2 + NV, &
178 : i_Fr_dFr_00 = i_Fr_dFr_in + NV, &
179 : i_Fr_dFr_out = i_Fr_dFr_00 + NV, &
180 : i_Fr_dFr_out2 = i_Fr_dFr_out + NV, &
181 :
182 : i_Fr_dr_in = i_Fr_dr_in2 + NV, &
183 : i_Fr_dr_00 = i_Fr_dr_in + NV, &
184 : i_Fr_dr_out = i_Fr_dr_00 + NV, &
185 :
186 : i_Fr_dT_00 = i_Fr_dT_in + NV, &
187 : i_Fr_dT_out = i_Fr_dT_00 + NV, &
188 : i_Fr_dT_out2 = i_Fr_dT_out + NV, &
189 :
190 : i_Fr_dw_00 = i_Fr_dw_in + NV, &
191 : i_Fr_dw_out = i_Fr_dw_00 + NV, &
192 : i_Fr_dw_out2 = i_Fr_dw_out + NV, &
193 :
194 : i_Fr_der_00 = i_Fr_der_in + NV, &
195 : i_Fr_der_out = i_Fr_der_00 + NV, &
196 : i_Fr_der_out2 = i_Fr_der_out + NV
197 :
198 : integer :: iter, min_k_for_turbulent_flux
199 :
200 : contains
201 :
202 0 : subroutine HYD(s,ierr)
203 : use star_utils, only: start_time, update_time
204 : type (star_info), pointer :: s
205 : integer, intent(out) :: ierr
206 0 : real(dp) :: DXXT, DXXC, DXXE, DXXL, PREC2, DXH, EZH, &
207 0 : P_surf, R_center_start, dt_max, total
208 : integer :: &
209 : i_min, i_max, num_tries, max_retries, max_iters, k, nz, &
210 : kT_max, kW_max, kE_max, kL_max, iter_for_dfridr, test_partials_k
211 : integer(i8) :: time0
212 : logical :: converged, dbg_msg, trace
213 : include 'formats'
214 :
215 0 : ierr = 0
216 0 : dbg_msg = s% report_solver_progress
217 0 : trace = s% trace_evolve
218 0 : iter = 0
219 0 : iter_for_dfridr = - 1
220 0 : test_partials_k = s% solver_test_partials_k
221 0 : s% solver_test_partials_k = 0
222 : if (s% model_number == s% solver_test_partials_call_number &
223 : .and. s% solver_test_partials_dx_0 > 0 &
224 : .and. s% solver_test_partials_iter_number > 0 &
225 0 : .and. test_partials_k > 0) then
226 0 : iter_for_dfridr = s% solver_test_partials_iter_number
227 0 : s% solver_test_partials_var = 0
228 0 : s% solver_test_partials_val = 0
229 0 : s% solver_test_partials_dval_dx = 0
230 : end if
231 0 : nz = s% nz
232 0 : i_min = 1
233 0 : i_max = nz
234 0 : PREC2 = s% RSP_tol_max_corr
235 0 : DXH = 0.3d0
236 0 : max_retries = s% RSP_max_retries_per_step
237 0 : max_iters = s% RSP_max_iters_per_try
238 0 : R_center_start = s% R_center
239 0 : call save_start_vars(s)
240 0 : do k=1,s% nz
241 0 : s% L(k) = s% Fr(k)*4*pi*s% r(k)**2 + s% Lc(k) + s% Lt(k)
242 0 : s% L_start(k) = s% L(k)
243 : end do
244 0 : call set_f_Edd(s,ierr)
245 0 : if (ierr /= 0) return
246 0 : P_surf = get_Psurf(s,ierr)
247 0 : if (ierr /= 0) return
248 0 : if (ALFAT > 0d0) then
249 0 : call set_min_k_for_turbulent_flux
250 : else
251 0 : min_k_for_turbulent_flux = 0
252 : end if
253 :
254 0 : if (s% v_center > s% v(nz) .and. s% v_center > 0d0) then ! compressing innermost cell
255 0 : dt_max = 1d-2*(s% r(nz) - s% r_center)/(s% v_center - s% v(nz))
256 0 : if (s% dt > dt_max) then
257 : !write(*,2) 'reduce dt in HYD', s% model_number, s% dt, dt_max
258 : !write(*,2) 's% r(nz) - s% r_center', s% model_number, s% r(nz) - s% r_center
259 : !write(*,2) 's% v_center - s% v(nz)', s% model_number, s% v_center - s% v(nz)
260 0 : if (s% RSP_report_limit_dt) then
261 0 : write(*,4) 'limit dt to avoid over-compressing innermost cell', s% model_number
262 0 : write(*,2) 's% r(nz) - s% r_center', s% model_number, s% r(nz) - s% r_center
263 0 : write(*,2) 's% v_center - s% v(nz)', s% model_number, s% v_center - s% v(nz)
264 0 : write(*,2) 'old dt', s% model_number, s% dt
265 0 : write(*,2) 'reduced dt', s% model_number, dt_max
266 0 : write(*,'(A)')
267 : end if
268 0 : call mesa_error(__FILE__,__LINE__,'HYD compressing innermost cell')
269 0 : s% dt = dt_max
270 : if (call_is_bad) then
271 : if (is_bad(s% dt)) then
272 : write(*,1) 'dt', s% dt
273 : call mesa_error(__FILE__,__LINE__,'HYD compressing innermost cell')
274 : end if
275 : end if
276 : end if
277 : end if
278 0 : if (s% dt < s% force_timestep_min .and. s% force_timestep_min > 0) &
279 0 : s% dt = s% force_timestep_min
280 0 : if (s% force_timestep > 0) s% dt = s% force_timestep
281 :
282 0 : retry_loop: do num_tries = 1, max_retries+1
283 :
284 0 : converged = .false.
285 0 : call set_1st_iter_R_using_v_start(s)
286 0 : s% R_center = s% R_center + s% dt*s% v_center
287 : ! write(*,3) 'RSP HYD w', 22, 0, s% RSP_w(22)
288 :
289 0 : iter_loop: do iter = 1, max_iters
290 0 : s% solver_iter = iter
291 0 : if (iter == iter_for_dfridr) then
292 0 : s% solver_test_partials_k = test_partials_k
293 : end if
294 0 : call eval_vars(s,iter,i_min,i_max,ierr)
295 0 : if (ierr /= 0) return
296 0 : call eval_eqns(s,P_surf)
297 0 : if (iter == iter_for_dfridr) call check_partial()
298 0 : if (converged) then
299 0 : s% num_solver_iterations = iter - 1
300 0 : s% solver_test_partials_k = test_partials_k
301 0 : return
302 : end if
303 0 : if (s% doing_timing) call start_time(s, time0, total)
304 0 : call solve_for_corrections(s,iter)
305 : call apply_corrections(s, &
306 0 : DXH,DXXT,DXXC,DXXE,DXXL,EZH,kT_max,kW_max,kE_max,kL_max)
307 0 : if (s% doing_timing) call update_time(s, time0, total, s% time_solver_matrix)
308 0 : if (dbg_msg) call write_msg
309 : ! write(*,3) 'RSP HYD w', 22, iter, s% RSP_w(22)
310 0 : if (iter == 1) cycle iter_loop
311 0 : converged = (abs(DXXT) < PREC2 .and. abs(DXXC) < PREC2)
312 : end do iter_loop
313 0 : call doing_retry
314 : end do retry_loop
315 :
316 0 : write(*,*) ' NO CONVERGENCE IN HYD, TIME STEP: ',s% model_number
317 0 : stop
318 :
319 : contains
320 :
321 0 : subroutine set_min_k_for_turbulent_flux
322 0 : real(dp) :: tau, rmid
323 : integer :: k
324 0 : tau = 0
325 0 : min_k_for_turbulent_flux = nz
326 0 : do k = 1, nz-1
327 0 : rmid = 0.5d0*(s% r(k) + s% r(k+1))
328 0 : tau = tau + s% dm(k)*s% opacity(k)/(4*pi*rmid**2)
329 0 : if (tau >= s% RSP_min_tau_for_turbulent_flux) then
330 0 : min_k_for_turbulent_flux = k
331 0 : exit
332 : end if
333 : end do
334 0 : end subroutine set_min_k_for_turbulent_flux
335 :
336 0 : subroutine doing_retry
337 : include 'formats'
338 0 : if (num_tries == max_retries+1) then
339 0 : write(*,*) 'NO CONVERGENCE IN HYD, TIME STEP num tries, max allowed', &
340 0 : s% model_number, num_tries, max_retries+1
341 0 : call mesa_error(__FILE__,__LINE__,'RSP: step num_retries = RSP_max_retries_per_step')
342 : end if
343 0 : call restore_start_vars(s)
344 0 : s% R_center = R_center_start
345 0 : s% dt = s% dt/2.d0
346 0 : s% num_retries = s% num_retries + 1
347 0 : if (s% max_number_retries < 0) return
348 0 : if (s% num_retries > s% max_number_retries) then
349 0 : write(*,3) 'model max_number_retries', s% model_number, s% max_number_retries
350 0 : call mesa_error(__FILE__,__LINE__,'RSP: num_retries > max_number_retries')
351 : end if
352 : end subroutine doing_retry
353 :
354 0 : subroutine write_msg
355 : include 'formats'
356 : !if (EZH < 1d0) write(*,3) 'undercorrection factor', s% model_number, iter, EZH
357 : write(*,'(i6, 2x, i3, 4(4x, a, 1x, i4, 1x, 1pe11.4, 1x, 1pe11.4))') &
358 0 : s% model_number, iter, &
359 0 : 'T', kT_max, DXXT, s% T(max(1,kT_max)), &
360 0 : 'w', kW_max, DXXC, max(1d-99,s% RSP_w(max(1,kW_max))), &
361 0 : 'erad', kE_max, DXXE, s% erad(max(1,kE_max)), &
362 0 : 'Fr', kL_max, DXXL, s% Fr(max(1,kL_max))
363 0 : end subroutine write_msg
364 :
365 0 : subroutine check_partial
366 : integer :: i_var
367 0 : real(dp) :: dvardx_0, dx_0, dvardx, xdum, err
368 : include 'formats'
369 0 : i_var = s% solver_test_partials_var
370 0 : dvardx_0 = s% solver_test_partials_dval_dx ! analytic partial
371 0 : if (i_var <= 0) then
372 0 : write(*,2) 'need to set test_partials_var', i_var
373 0 : call mesa_error(__FILE__,__LINE__,'check_partial')
374 : end if
375 0 : dx_0 = get1_val(i_var, s% solver_test_partials_k)
376 0 : dx_0 = s% solver_test_partials_dx_0*max(1d-99, abs(dx_0))
377 0 : dvardx = dfridr(dx_0,err)
378 0 : xdum = (dvardx - dvardx_0)/max(abs(dvardx_0),1d-50)
379 0 : write(*,1) 'analytic numeric err rel_diff',dvardx_0,dvardx,err,xdum
380 : !write(*,*)
381 0 : call mesa_error(__FILE__,__LINE__,'check_partial')
382 0 : end subroutine check_partial
383 :
384 0 : real(dp) function get1_val(i_var,k) result(val)
385 : integer, intent(in) :: i_var, k
386 : include 'formats'
387 0 : if (i_var == i_var_w) then
388 0 : val = s% RSP_w(k)
389 : else if (i_var == i_var_R) then
390 0 : val = s% r(k)
391 : else if (i_var == i_var_T) then
392 0 : val = s% T(k)
393 : else if (i_var == i_var_er) then
394 0 : val = s% erad(k)
395 : else if (i_var == i_var_Fr) then
396 0 : val = s% Fr(k)
397 : else if (i_var == i_var_Vol) then
398 0 : val = s% Vol(k)
399 : else
400 0 : write(*,2) 'bad value for solver_test_partials_var', i_var
401 0 : call mesa_error(__FILE__,__LINE__,'solver_test_partials')
402 : end if
403 0 : end function get1_val
404 :
405 0 : subroutine store1_val(i_var, k, val)
406 : integer, intent(in) :: i_var, k
407 : real(dp), intent(in) :: val
408 : include 'formats'
409 0 : if (i_var == i_var_w) then
410 0 : s% RSP_w(k) = val
411 : else if (i_var == i_var_R) then
412 0 : s% r(k) = val
413 0 : s% v(k) = 2.d0*(s% r(k) - s% r_start(k))/s% dt - s% v_start(k)
414 : ! partials wrt R assume v changes along with R,
415 : ! so need to do it here.
416 : else if (i_var == i_var_T) then
417 0 : s% T(k) = val
418 : else if (i_var == i_var_er) then
419 0 : s% erad(k) = val
420 : else if (i_var == i_var_Fr) then
421 0 : s% Fr(k) = val
422 : else if (i_var == i_var_Vol) then
423 0 : s% Vol(k) = val
424 : else
425 0 : write(*,2) 'bad value for solver_test_partials_var', i_var
426 0 : call mesa_error(__FILE__,__LINE__,'solver_test_partials')
427 : end if
428 0 : end subroutine store1_val
429 :
430 0 : real(dp) function dfridr_func(delta_x) result(val)
431 : real(dp), intent(in) :: delta_x
432 : integer :: i_var, k, ierr
433 : real(dp) :: save1
434 : include 'formats'
435 0 : i_var = s% solver_test_partials_var
436 0 : k = s% solver_test_partials_k
437 0 : save1 = get1_val(i_var, k)
438 0 : call store1_val(i_var, k, save1 + delta_x)
439 0 : call eval_vars(s,0,i_min,i_max,ierr)
440 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'failed in eval_vars')
441 0 : call eval_eqns(s,P_surf)
442 0 : val = s% solver_test_partials_val
443 : !write(*,2) 'dfridr val', k, val
444 0 : call store1_val(i_var, k, save1)
445 0 : end function dfridr_func
446 :
447 0 : real(dp) function dfridr(hx,err)
448 : real(dp), intent(in) :: hx
449 : real(dp), intent(out) :: err
450 : ! this routine returns the first derivative of a function func(x)
451 : ! at the point x, by ridders method of polynomial extrapolation.
452 : ! value hx is the initial step size;
453 : ! it should be an increment for which func changes substantially.
454 : ! an estimate of the error in the first derivative is returned in err.
455 : integer, parameter :: ntab = 20
456 : integer :: i,j
457 0 : real(dp) :: errt,fac,hh,a(ntab,ntab),f1,f2
458 : real(dp), parameter :: con2 = 2d0, con = sqrt(con2), big = 1d199, safe = 2d0
459 : include 'formats'
460 0 : dfridr = 0d0
461 0 : hh = hx
462 : ! 2nd order central difference
463 0 : f1 = dfridr_func(hh)
464 : !write(*,2) 'f1', 1, f1, save_dx(i_var,k) + hh
465 0 : f2 = dfridr_func( - hh)
466 : !write(*,2) 'f2', 1, f2, save_dx(i_var,k) - hh
467 0 : a(1,1) = (f1 - f2)/(2d0*hh)
468 : !write(*,2) 'dfdx', 1, a(1,1), hh
469 0 : err = big
470 : ! successive columns in the neville tableu will go to smaller stepsizes
471 : ! and higher orders of extrapolation
472 0 : do i = 2,ntab
473 0 : hh = hh/con
474 0 : f1 = dfridr_func(hh)
475 : !write(*,2) 'f1', i, f1
476 0 : f2 = dfridr_func( - hh)
477 : !write(*,2) 'f2', i, f2
478 0 : a(1,i) = (f1 - f2)/(2d0*hh)
479 : !write(*,2) 'dfdx', i, a(1,i), hh
480 : ! compute extrapolations of various orders; the error strategy is to compare
481 : ! each new extrapolation to one order lower but both at the same stepsize
482 : ! and at the previous stepsize
483 0 : fac = con2
484 0 : do j = 2,i
485 0 : a(j,i) = (a(j - 1,i)*fac - a(j - 1,i - 1))/(fac - 1d0)
486 0 : fac = con2*fac
487 0 : errt = max(abs(a(j,i) - a(j - 1,i)),abs(a(j,i) - a(j - 1,i - 1)))
488 : !write(*,2) 'errt, err', j, errt, err
489 0 : if (errt <= err) then
490 0 : err = errt
491 0 : dfridr = a(j,i)
492 0 : write(*,3) 'dfridr err', i, j, dfridr, err
493 : end if
494 : end do
495 : ! if higher order is worse by a significant factor safe, then bail
496 0 : if (abs(a(i,i) - a(i - 1,i - 1)) >= safe*err) then
497 : !write(*,1) 'higher order is worse'
498 0 : write(*,'(A)')
499 0 : return
500 : end if
501 : end do
502 0 : end function dfridr
503 :
504 : end subroutine HYD
505 :
506 :
507 0 : real(dp) function get_Psurf(s,ierr) result(P_surf)
508 : use rsp_eval_eos_and_kap, only: get_surf_P_T_kap
509 : type (star_info), pointer :: s
510 : integer, intent(out) :: ierr
511 0 : real(dp) :: kap_guess, T_surf, kap_surf, Teff_atm
512 : include 'formats'
513 0 : ierr = 0
514 0 : if (s% RSP_fixed_Psurf) then
515 0 : P_surf = Psurf_from_atm
516 0 : else if (s% RSP_use_atm_grey_with_kap_for_Psurf) then
517 : ierr = 0
518 0 : kap_guess = 1d-2
519 : call get_surf_P_T_kap(s, &
520 : s% M(1), s% r(1), s% L(1), &
521 : (2d0/3d0)*s% tau_factor, kap_guess, &
522 0 : T_surf, P_surf, kap_surf, Teff_atm, ierr)
523 0 : if (ierr/= 0) then
524 0 : write(*,*) 'ierr from get_surf_P_T_kap'
525 0 : return
526 : end if
527 0 : else if (s% RSP_use_Prad_for_Psurf) then
528 0 : P_surf = crad*s% T(1)**4/3d0
529 0 : else if (s% RSP_Psurf >= 0d0) then
530 0 : P_surf = s% RSP_Psurf
531 : else
532 0 : P_surf = 0d0
533 : end if
534 0 : end function get_Psurf
535 :
536 :
537 0 : subroutine calculate_energies(s,total_radiation)
538 0 : use star_utils, only: cell_specific_KE, cell_specific_PE
539 : type (star_info), pointer :: s
540 : real(dp), intent(out) :: total_radiation
541 : integer :: i, k
542 0 : real(dp) :: d_dlnR00, d_dlnRp1, d_dv00, d_dvp1
543 : include 'formats'
544 0 : EGRV = 0.d0
545 0 : ETHE = 0.d0
546 0 : EKIN = 0.d0
547 0 : ECON = 0.d0
548 0 : do i=1,NZN
549 0 : k = NZN+1-i
550 0 : ETHE = ETHE + (s% egas(k)+s% erad(k))*s% dm(k)
551 0 : ECON = ECON + s% RSP_w(k)**2*s% dm(k)
552 0 : EKIN = EKIN + cell_specific_KE(s,k,d_dv00,d_dvp1)*s% dm(k)
553 0 : EGRV = EGRV + cell_specific_PE(s,k,d_dlnR00,d_dlnRp1)*s% dm(k)
554 :
555 : !EKIN = EKIN + 0.5d0*s% v(k)**2*s% dm_bar(k)
556 : !if (k < NZN) then
557 : ! EGRV = EGRV - s% cgrav(k) * (s%m(k)-0.5d0*s%dm(k))*s%dm(k)/(0.5d0*(s%r(k)+s%r(k+1)))
558 : !else
559 : ! EGRV = EGRV - s% cgrav(k) * (s%m(k)-0.5d0*s%dm(k))*s%dm(k)/(0.5d0*(s%r(k)+s%r_center))
560 : !end if
561 :
562 : end do
563 0 : if (s% RSP_hydro_only) then
564 0 : total_radiation = 0d0
565 : else
566 0 : total_radiation = s% dt*(s% L_center - (WTR*s% L(1) + WTR1*s% L_start(1)))
567 : end if
568 0 : ETOT = ETHE+EKIN+ECON+EGRV
569 0 : EDE_start = EDE_start-E0
570 0 : if (EKIN > EKMAX) EKMAX = EKIN
571 0 : if (EKIN < EKMIN) EKMIN = EKIN
572 : !write(*,1) 'ETHE', ETHE
573 : !write(*,1) 'ECON', ECON
574 : !write(*,1) 'EKIN', EKIN
575 : !write(*,1) 'EGRV', EGRV
576 : !write(*,1) 'ETOT', ETOT
577 0 : end subroutine calculate_energies
578 :
579 :
580 0 : subroutine init_HYD(s,ierr)
581 : type (star_info), pointer :: s
582 : integer, intent(out) :: ierr
583 : integer :: i_min, i_max
584 : include 'formats'
585 0 : i_min = 1
586 0 : i_max = s% nz
587 : ! setup so can save start vals in 1st call to HYD
588 0 : s% f_Edd(1:NZN) = f_Edd_isotropic ! fake it for 1st call on eval_vars
589 0 : call eval_vars(s,0,i_min,i_max,ierr)
590 0 : if (ierr /= 0) return
591 0 : call set_f_Edd(s,ierr) ! needs opacities
592 0 : if (ierr /= 0) return
593 0 : call save_start_vars(s) ! needed by eval_eqns
594 0 : call eval_eqns(s,0d0)
595 0 : end subroutine init_HYD
596 :
597 :
598 0 : subroutine set_1st_iter_R_using_v_start(s)
599 : type (star_info), pointer :: s
600 0 : real(dp) :: EH1,EHJT
601 : integer :: k
602 : include 'formats'
603 0 : EHJT = 1.d0
604 0 : do k = 1,NZN
605 0 : if (k /= NZN) then
606 : EH1 = - (s% v_start(k) - s% v_start(k+1))*s% dt/ &
607 0 : (s% r_start(k) - s% r_start(k+1))/0.8d0
608 : else
609 : EH1 = - (s% v_start(k) - s% v_center)*s% dt/ &
610 0 : (s% r_start(k) - s% R_center)/0.8d0
611 : end if
612 0 : EHJT = 1.d0/max(1.d0/EHJT,EH1,1d0)
613 : end do
614 0 : do k = 1,NZN
615 0 : s% r(k) = s% r_start(k) + EHJT*s% dt*s% v_start(k) ! initial guess for R
616 0 : s% v(k) = (2d0*EHJT - 1d0)*s% v_start(k) ! initial guess for v
617 : end do
618 0 : end subroutine set_1st_iter_R_using_v_start
619 :
620 :
621 0 : subroutine rsp_set_Teff(s)
622 : use star_utils, only: get_phot_info
623 : use atm_lib, only: atm_Teff
624 : type (star_info), pointer :: s
625 : real(dp) :: r, m, v, L, T_phot, cs, kap, logg, ysum
626 : integer :: k_phot
627 : include 'formats'
628 0 : call get_phot_info(s,r,m,v,L,T_phot,cs,kap,logg,ysum,k_phot)
629 0 : s% Teff = atm_Teff(L,r)
630 0 : end subroutine rsp_set_Teff
631 :
632 :
633 0 : subroutine set_f_Edd(s, ierr)
634 : type (star_info), pointer :: s
635 : integer, intent(out) :: ierr
636 : include 'formats'
637 0 : ierr = 0
638 0 : s% g_Edd = 0.5d0
639 0 : s% f_Edd(1:NZN) = f_Edd_isotropic
640 0 : end subroutine set_f_Edd
641 :
642 :
643 0 : subroutine save_start_vars(s)
644 : type (star_info), pointer :: s
645 : integer :: I, k
646 0 : do I = 1,NZN
647 0 : k = NZN+1-i
648 0 : s% T_start(k) = s% T(k)
649 0 : s% r_start(k) = s% r(k)
650 0 : s% Pgas_start(k) = s% Pgas(k)
651 0 : s% Prad_start(k) = s% Prad(k)
652 0 : s% Pvsc_start(k) = s% Pvsc(k)
653 0 : s% Vol_start(k) = s% Vol(k)
654 0 : s% csound_start(k) = s% csound(k)
655 0 : s% opacity_start(k) = s% opacity(k)
656 0 : s% egas_start(k) = s% egas(k)
657 0 : s% erad_start(k) = s% erad(k)
658 0 : s% RSP_w_start(k) = s% RSP_w(k)
659 0 : s% Ptrb_start(k) = s% Ptrb(k)
660 0 : s% Chi_start(k) = s% Chi(k)
661 0 : s% v_start(k) = s% v(k)
662 0 : s% Fr_start(k) = s% Fr(k)
663 0 : s% Lc_start(k) = s% Lc(k)
664 0 : s% Lt_start(k) = s% Lt(k)
665 0 : s% COUPL_start(k) = s% COUPL(k)
666 : end do
667 0 : end subroutine save_start_vars
668 :
669 :
670 0 : subroutine restore_start_vars(s)
671 : type (star_info), pointer :: s
672 : integer :: I, k
673 0 : do I = 1,NZN
674 0 : k = NZN+1-i
675 0 : s% T(k) = s% T_start(k)
676 0 : s% r(k) = s% r_start(k)
677 0 : s% Pgas(k) = s% Pgas_start(k)
678 0 : s% Prad(k) = s% Prad_start(k)
679 0 : s% Pvsc(k) = s% Pvsc_start(k)
680 0 : s% Vol(k) = s% Vol_start(k)
681 0 : s% csound(k) = s% csound_start(k)
682 0 : s% opacity(k) = s% opacity_start(k)
683 0 : s% egas(k) = s% egas_start(k)
684 0 : s% erad(k) = s% erad_start(k)
685 0 : s% RSP_w(k) = s% RSP_w_start(k)
686 0 : s% Ptrb(k) = s% Ptrb_start(k)
687 0 : s% Chi(k) = s% Chi_start(k)
688 0 : s% v(k) = s% v_start(k)
689 0 : s% L(k) = s% L_start(k)
690 0 : s% Fr(k) = s% Fr_start(k)
691 0 : s% Lc(k) = s% Lc_start(k)
692 0 : s% Lt(k) = s% Lt_start(k)
693 0 : s% COUPL(k) = s% COUPL_start(k)
694 : end do
695 0 : end subroutine restore_start_vars
696 :
697 :
698 0 : subroutine eval_vars(s,iter,i_min,i_max,ierr)
699 : use star_utils, only: start_time, update_time
700 : type (star_info), pointer :: s
701 : integer, intent(in) :: iter,i_min,i_max
702 : integer, intent(out) :: ierr
703 : integer :: i, op_err
704 : integer(i8) :: time0
705 0 : real(dp) :: total
706 : include 'formats'
707 0 : ierr = 0
708 0 : if (s% doing_timing) call start_time(s, time0, total)
709 0 : !$OMP PARALLEL DO PRIVATE(I,op_err) SCHEDULE(dynamic,2)
710 : do i = i_min,i_max
711 : call do1_specific_volume(s,i)
712 : call do1_eos_and_kap(s,i,op_err)
713 : if (op_err /= 0) ierr = op_err
714 : call calc_Prad(s,i)
715 : end do
716 : !$OMP END PARALLEL DO
717 0 : if (s% doing_timing) call update_time(s, time0, total, s% time_eos)
718 0 : if (ierr /= 0) return
719 0 : !$OMP PARALLEL DO PRIVATE(I) SCHEDULE(dynamic,2)
720 : do i = i_min,i_max
721 : call calc_Hp_face(s,i)
722 : call calc_Y_face(s,i)
723 : call calc_PII_face(s,i)
724 : call calc_Pvsc(s,i)
725 : end do
726 : !$OMP END PARALLEL DO
727 0 : if (iter == 1 .and. s% RSP_do_check_omega) then
728 0 : do i = i_min,i_max
729 0 : call check_omega(s,i)
730 : end do
731 : end if
732 0 : !$OMP PARALLEL DO PRIVATE(I) SCHEDULE(dynamic,2)
733 : do i = i_min,i_max
734 : call calc_Pturb(s,i)
735 : call calc_Chi(s,i)
736 : call calc_Eq(s,i)
737 : call calc_source_sink(s,i)
738 : end do
739 : !$OMP END PARALLEL DO
740 0 : call zero_boundaries(s)
741 0 : end subroutine eval_vars
742 :
743 :
744 0 : subroutine eval_eqns(s,P_surf)
745 : type (star_info), pointer :: s
746 : real(dp), intent(in) :: P_surf
747 : integer :: i
748 : include 'formats'
749 0 : !$OMP PARALLEL DO PRIVATE(I) SCHEDULE(dynamic,2)
750 : do i = 1,NZN
751 : call calc_equations(s, i, P_surf)
752 : end do
753 : !$OMP END PARALLEL DO
754 0 : end subroutine eval_eqns
755 :
756 :
757 0 : subroutine do1_specific_volume(s,i)
758 : type (star_info), pointer :: s
759 : integer, intent(in) :: i
760 : integer :: k
761 0 : real(dp) :: T1
762 0 : real(qp) :: q1, q2, q3, q4
763 : include 'formats'
764 0 : k = NZN+1-i
765 0 : T1 = P43/s% dm(k)
766 0 : if (i /= 1) then
767 0 : q1 = T1
768 0 : q2 = s% r(k)
769 0 : q3 = s% r(k+1)
770 0 : q4 = q1*(q2**3 - q3**3)
771 0 : s% Vol(k) = dble(q4)
772 0 : dVol_dr_in(i) = - 3.0d0*T1*s% r(k+1)**2
773 : else
774 0 : s% Vol(k) = T1*(s% r(k)**3 - s% R_center**3)
775 0 : dVol_dr_in(i) = 0d0
776 : end if
777 0 : if (s% Vol(k) <= 0d0) then
778 0 : write(*,2) 'bad Vol', k, s% Vol(k)
779 0 : call mesa_error(__FILE__,__LINE__,'do1_specific_volume')
780 : end if
781 0 : dVol_dr_00(i) = 3.d0*T1*s% r(k)**2
782 0 : end subroutine do1_specific_volume
783 :
784 :
785 0 : subroutine solve_for_corrections(s,iter)
786 : type (star_info), pointer :: s
787 : integer, intent(in) :: iter
788 : integer :: i, j, info, IR, IT, IW, IE, IL, N
789 : logical :: okay
790 : include 'formats'
791 :
792 : !if (s% model_number >= s% max_model_number-1 .and. iter == 1) then
793 : !if (.false. .and. s% model_number == s% max_model_number) then !.and. iter == 1) then
794 : if (.false.) then
795 : ! this can be useful for finding restart bugs.
796 : write(*,*) 'solve_for_corrections'
797 : okay = .true.
798 : do i = 1,NZN
799 : IR = i_var_R + NV*(i-1)
800 : IT = i_var_T + NV*(i-1)
801 : IW = i_var_w + NV*(i-1)
802 : IE = i_var_er + NV*(i-1)
803 : IL = i_var_Fr + NV*(i-1)
804 :
805 : if (is_bad(HR(IR))) then
806 : write(*,4) 'HR(IR)', iter, i, IR, HR(IR)
807 : okay = .false.
808 : end if
809 : if (is_bad(HR(IT))) then
810 : write(*,4) 'HR(IT)', iter, i, IT, HR(IT)
811 : okay = .false.
812 : end if
813 : if (is_bad(HR(IW))) then
814 : write(*,4) 'HR(IW)', iter, i, IW, HR(IW)
815 : okay = .false.
816 : end if
817 : if (is_bad(HR(IL))) then
818 : write(*,4) 'HR(IL)', iter, i, IL, HR(IL)
819 : okay = .false.
820 : end if
821 :
822 : if (.not. okay) call mesa_error(__FILE__,__LINE__,'solve_for_corrections')
823 :
824 : do j = 1,LD_HD
825 : if (is_bad(HD(j,IR))) then
826 : write(*,5) 'HD(j,IR)', iter, i, j, IR, HD(j,IR)
827 : okay = .false.
828 : end if
829 : if (is_bad(HD(j,IT))) then
830 : write(*,5) 'HD(j,IT)', iter, i, j, IT, HD(j,IT)
831 : okay = .false.
832 : end if
833 : if (is_bad(HD(j,IW))) then
834 : write(*,5) 'HD(j,IW)', iter, i, j, IW, HD(j,IW)
835 : okay = .false.
836 : end if
837 : if (is_bad(HD(j,IE))) then
838 : write(*,5) 'HD(j,IE)', iter, i, j, IE, HD(j,IE)
839 : okay = .false.
840 : end if
841 : if (is_bad(HD(j,IL))) then
842 : write(*,5) 'HD(j,IL)', iter, i, j, IL, HD(j,IL)
843 : okay = .false.
844 : end if
845 :
846 : if (.not. okay) call mesa_error(__FILE__,__LINE__,'solve_for_corrections')
847 :
848 : end do
849 :
850 : cycle
851 :
852 : write(*,3) 'HR(IR)', iter, i, HR(IR)
853 : write(*,3) 'HR(IT)', iter, i, HR(IT)
854 : write(*,3) 'HR(IW)', iter, i, HR(IW)
855 : write(*,3) 'HR(IL)', iter, i, HR(IL)
856 :
857 : do j = 1,13
858 : write(*,5) 'HD(j,IR)', iter, j, IR, i, HD(j,IR)
859 : end do
860 :
861 : do j = 1,13
862 : write(*,5) 'HD(j,IT)', iter, j, IT, i, HD(j,IT)
863 : end do
864 :
865 : do j = 1,13
866 : write(*,5) 'HD(j,IW)', iter, j, IW, i, HD(j,IW)
867 : end do
868 :
869 : write(*,3) 'HR(IE)', iter, i, HR(IE)
870 : do j = 1,13
871 : write(*,5) 'HD(j,IE)', iter, j, IE, i, HD(j,IE)
872 : end do
873 :
874 : do j = 1,13
875 : write(*,5) 'HD(j,IL)', iter, j, IL, i, HD(j,IL)
876 : end do
877 :
878 : end do
879 : !call mesa_error(__FILE__,__LINE__,'solve_for_corrections')
880 : end if
881 :
882 0 : N = NV*NZN+1
883 :
884 : if (.false.) then ! check HR and HD for NaN's
885 : do I = 1,N
886 : if (is_bad(HR(I))) then
887 : write(*,3) 'HR(I)', iter, I, HR(I)
888 : call mesa_error(__FILE__,__LINE__,'solve_for_corrections')
889 : end if
890 : end do
891 : do I = 1,N
892 : do j = 1,LD_HD
893 : if (is_bad(HD(j,i))) then
894 : write(*,4) 'HD(j,i)', iter, j, i, HD(j,i)
895 : call mesa_error(__FILE__,__LINE__,'solve_for_corrections')
896 : end if
897 : end do
898 : end do
899 : end if
900 :
901 0 : do J = 1,2*NV ! translate hd into band storage of LAPACK
902 0 : do I = 1,N
903 0 : ABB(J,I) = 0.0d0
904 : end do
905 : end do
906 0 : do J = 1,2*NV
907 0 : do I = 1,N - J
908 0 : ABB(LD_HD - J,I + J) = HD(HD_DIAG + J,I) ! upper diagonals
909 : end do
910 : end do
911 0 : do J = 1,2*NV
912 0 : do I = 1,N - J
913 0 : ABB(LD_HD + J,I) = HD(HD_DIAG - J,I + J) ! lower diagonals
914 : end do
915 : end do
916 0 : do I = 1,N
917 0 : ABB(LD_HD,I) = HD(HD_DIAG,I)
918 : end do
919 :
920 0 : call DGBTRF(N,N,2*NV,2*NV,ABB,LD_ABB,IPVT,INFO)
921 0 : if (INFO/= 0) then
922 0 : write(*,*) 'hyd: LAPACK/dgbtrf problem',INFO
923 0 : stop
924 : end if
925 :
926 0 : call DGBTRS('n',N,2*NV,2*NV,1,ABB,LD_ABB,IPVT,HR,N,INFO)
927 0 : if (INFO/= 0) then
928 0 : write(*,*) 'hyd: LAPACK/dgbtrs problem',INFO
929 0 : stop
930 : end if
931 :
932 0 : do I = 1,N
933 0 : DX(I) = HR(I)
934 0 : if (call_is_bad) then
935 : if (is_bad(DX(I))) then
936 : write(*,2) 'DX(I)', I, DX(I)
937 : call mesa_error(__FILE__,__LINE__,'solve_for_corrections')
938 : end if
939 : end if
940 : end do
941 :
942 0 : end subroutine solve_for_corrections
943 :
944 :
945 0 : subroutine apply_corrections(s, &
946 : DXH, XXT, XXC, XXE, XXL, EZH, &
947 : kT_max, kW_max, kE_max, kL_max)
948 : use star_utils, only: rand
949 : type (star_info), pointer :: s
950 : real(dp), intent(in) :: DXH
951 : real(dp), intent(out) :: XXT, XXC, XXE, XXL, EZH
952 : integer, intent(out) :: kT_max,kW_max,kE_max,kL_max
953 : integer :: i, k, IR, IT, IW, IE, IL, kEZH, &
954 : iTM, kTM, iRM, kRM, iEM, kEM, iCM, kCM, iLM, kLM
955 0 : real(dp) :: old_w, XXR, DXXT, DXXC, DXXE, DXXL, DXRM, DXR, &
956 0 : EZH1, XXTM, XXCM, XXRM, XXEM, XXLM, DXKT, DXKC, DXKE, DXKL
957 : include 'formats'
958 0 : EZH = 1.0d0; kEZH = 0
959 0 : XXTM = 0d0; kTM = 0; iTM = 0
960 0 : XXRM = 0d0; kRM = 0; iRM = 0
961 0 : XXEM = 0d0; kEM = 0; iEM = 0
962 0 : XXCM = 0d0; kCM = 0; iCM = 0
963 0 : XXLM = 0d0; kLM = 0; iLM = 0
964 0 : do I = 1,NZN
965 0 : k = NZN+1-i
966 0 : IR = i_var_R + NV*(i-1)
967 0 : IT = i_var_T + NV*(i-1)
968 0 : IW = i_var_w + NV*(i-1)
969 0 : IE = i_var_er + NV*(i-1)
970 0 : IL = i_var_Fr + NV*(i-1)
971 0 : if (s% RSP_w(k) > (1.d+2)*EFL0) then
972 0 : XXC = abs(DX(IW)/s% RSP_w(k))/DXH
973 0 : if (XXC > XXCM) then
974 0 : XXCM = XXC; kCM = k; iCM = IW
975 : end if
976 : else
977 0 : XXC = 0d0
978 : end if
979 0 : if (i == 1) then
980 0 : DXR = -DX(IR)
981 0 : XXR = (DXR/(s% r(k) - s% r_center))/DXH
982 0 : if (XXR > XXRM) then
983 0 : DXRM = DXR/(s% r(k) - s% r_center)
984 0 : XXRM = XXR; kRM = k; iRM = IR
985 : end if
986 : else
987 0 : DXR = DX(IR-NV) - DX(IR)
988 0 : XXR = (DXR/(s% r(k) - s% r(k+1)))/DXH
989 0 : if (XXR > XXRM) then
990 0 : DXRM = DXR/(s% r(k) - s% r(k+1))
991 0 : XXRM = XXR; kRM = k; iRM = IR
992 : end if
993 : end if
994 0 : XXE = abs(DX(IE)/s% erad(k))/DXH
995 0 : if (XXE > XXEM) then
996 0 : XXEM = XXE; kEM = k; iEM = IE
997 : end if
998 0 : XXL = abs(DX(IL)/s% Fr(k))/DXH
999 0 : if (XXL > XXLM) then
1000 0 : XXLM = XXL; kLM = k; iLM = IL
1001 : end if
1002 0 : XXT = abs(DX(IT)/s% T(k))/DXH
1003 0 : if (XXT > XXTM) then
1004 0 : XXTM = XXT; kTM = k; iTM = IT
1005 : end if
1006 0 : EZH1 = EZH
1007 0 : EZH = 1.d0/max(1.d0/EZH,XXR,XXT,XXC,XXE)
1008 0 : if (EZH1 /= EZH) kEZH = k
1009 : end do
1010 :
1011 0 : if (EZH < 1d0 .and. s% RSP_report_undercorrections) then
1012 : write(*,'(i6, 2x, i3, 6(4x, a, 1x, i4, 1x, 1pe10.3))') &
1013 0 : s% model_number, iter, &
1014 0 : 'EZH', kEZH, EZH, &
1015 0 : 'r', kRM, DXRM, &
1016 0 : 'T', kTM, DX(iTM)/s% T(kTM), &
1017 0 : 'w', kCM, DX(iCM)/max(1d-99,s% RSP_w(kCM)), &
1018 0 : 'erad', kEM, DX(iEM)/s% erad(kEM), &
1019 0 : 'Fr', kLM, DX(iLM)/s% Fr(kLM)
1020 : end if
1021 :
1022 0 : DXXT = 0.0d0
1023 0 : DXXC = 0.0d0
1024 0 : DXXE = 0.0d0
1025 0 : DXXL = 0.0d0
1026 0 : XXT = 0.0d0
1027 0 : XXC = 0.0d0
1028 0 : XXE = 0.0d0
1029 0 : XXL = 0.0d0
1030 0 : kT_max = 0
1031 0 : kW_max = 0
1032 0 : kE_max = 0
1033 0 : kL_max = 0
1034 :
1035 0 : do I = 1,NZN
1036 0 : k = NZN+1-i
1037 0 : IT = i_var_T + NV*(i-1)
1038 0 : IR = i_var_R + NV*(i-1)
1039 0 : IW = i_var_w + NV*(i-1)
1040 0 : IE = i_var_er + NV*(i-1)
1041 0 : IL = i_var_Fr + NV*(i-1)
1042 0 : s% T(k) = s% T(k) + EZH*DX(IT)
1043 0 : s% erad(k) = s% erad(k) + EZH*DX(IE)
1044 0 : if (I > IBOTOM .and. I < NZN)then
1045 0 : if ((s% RSP_w(k) + EZH*DX(IW)) <= 0d0)then
1046 0 : old_w = s% RSP_w(k)
1047 0 : s% RSP_w(k) = EFL0*rand(s)*1d-6 ! RSP NEEDS THIS to give seed for SOURCE
1048 : !write(*,3) 'rand RSP_w', k, iter, s% model_number, &
1049 : ! s% RSP_w(k), old_w, EZH*DX(IW), EZH
1050 : else
1051 0 : s% RSP_w(k) = s% RSP_w(k) + EZH*DX(IW)
1052 : end if
1053 : end if
1054 : s% v(k) = s% v(k) - &
1055 : (s% v(k) + s% v_start(k)) + &
1056 0 : 2.d0/s% dt*(EZH*DX(IR) + (s% r(k) - s% r_start(k)))
1057 0 : s% r(k) = s% r(k) + EZH*DX(IR)
1058 0 : s% Fr(k) = s% Fr(k) + EZH*DX(IL)
1059 0 : s% Lr(k) = s% Fr(k)*4d0*pi*s% r(k)**2
1060 0 : DXKT = DXXT
1061 0 : DXKC = DXXC
1062 0 : DXKE = DXXE
1063 0 : DXKL = DXXL
1064 0 : DXXT = max(DXXT,abs(DX(IT)/s% T(k)))
1065 0 : DXXE = max(DXXE,abs(DX(IE)/s% erad(k)))
1066 0 : DXXL = max(DXXL,abs(DX(IL)/s% Fr(k)))
1067 0 : if (s% RSP_w(k) > (1.d-2)*EFL0) &
1068 0 : DXXC = max(DXXC,abs(DX(IW)/s% RSP_w(k)))
1069 0 : if (DXXC > DXKC) then
1070 0 : kW_max = k; XXC = DX(IW)/max(1d-99,s% RSP_w(k))
1071 : end if
1072 0 : if (DXXT > DXKT) then
1073 0 : kT_max = k; XXT = DX(IT)/s% T(k)
1074 : end if
1075 0 : if (DXXE > DXKE) then
1076 0 : kE_max = k; XXE = DX(IE)/s% erad(k)
1077 : end if
1078 0 : if (DXXL > DXKL) then
1079 0 : kL_max = k; XXL = DX(IL)/s% Fr(k)
1080 : end if
1081 : end do
1082 :
1083 0 : end subroutine apply_corrections
1084 :
1085 :
1086 0 : subroutine do1_eos_and_kap(s,i,ierr)
1087 0 : use rsp_eval_eos_and_kap, only : eval_mesa_eos_and_kap
1088 : type (star_info), pointer :: s
1089 : integer, intent(in) :: i
1090 : integer, intent(out) :: ierr
1091 : real(dp) :: Prad, d_Pr_dT, erad, d_erad_dVol, d_erad_dT
1092 : integer :: k
1093 : include 'formats'
1094 0 : k = NZN + 1 - i
1095 : call eval_mesa_eos_and_kap(&
1096 : s, k, s% T(k), s% Vol(k), &
1097 : s% Pgas(k), d_Pg_dVol(I), d_Pg_dT(I), &
1098 : Prad, d_Pr_dT, &
1099 : s% egas(k), d_egas_dVol(I), d_egas_dT(I), &
1100 : erad, d_erad_dVol, d_erad_dT, &
1101 : s% csound(k), s% Cp(k), dCp_dVol(I), dCp_dT(I), &
1102 : s% QQ(k), dQQ_dVol(I), dQQ_dT(I), &
1103 0 : s% opacity(k), dK_dVol(I), dK_dT(I),ierr)
1104 0 : if (ierr /= 0) return
1105 0 : d_Pg_dr_00(I) = d_Pg_dVol(I)*dVol_dr_00(I)
1106 0 : d_Pg_dr_in(I) = d_Pg_dVol(I)*dVol_dr_in(I)
1107 0 : d_egas_dr_00(I) = d_egas_dVol(I)*dVol_dr_00(I)
1108 0 : d_egas_dr_in(I) = d_egas_dVol(I)*dVol_dr_in(I)
1109 0 : dCp_dr_00(I) = dCp_dVol(I)*dVol_dr_00(I)
1110 0 : dCp_dr_in(I) = dCp_dVol(I)*dVol_dr_in(I)
1111 0 : dQQ_dr_00(I) = dQQ_dVol(I)*dVol_dr_00(I)
1112 0 : dQQ_dr_in(I) = dQQ_dVol(I)*dVol_dr_in(I)
1113 0 : dK_dr_00(I) = dK_dVol(I)*dVol_dr_00(I)
1114 0 : dK_dr_in(I) = dK_dVol(I)*dVol_dr_in(I)
1115 : if (call_is_bad) then
1116 : if (is_bad(s% opacity(k))) then
1117 : !$omp critical (rsp_step_1)
1118 : write(*,2) 's% opacity(k)', k, s% opacity(k)
1119 : call mesa_error(__FILE__,__LINE__,'do1_eos_and_kap')
1120 : !$omp end critical (rsp_step_1)
1121 : end if
1122 : end if
1123 0 : end subroutine do1_eos_and_kap
1124 :
1125 :
1126 0 : subroutine calc_Prad(s,i)
1127 : type (star_info), pointer :: s
1128 : integer, intent(in) :: i
1129 : integer :: k
1130 0 : real(dp) :: V
1131 : logical :: test_partials
1132 : include 'formats'
1133 0 : k = NZN+1-i
1134 0 : V = s% Vol(k)
1135 0 : s% Prad(k) = s% f_Edd(k)*s% erad(k)/V
1136 0 : d_Pr_der(i) = s% f_Edd(k)/V
1137 0 : d_Pr_dVol(i) = -s% Prad(k)/V
1138 0 : d_Pr_dr_00(i) = d_Pr_dVol(i)*dVol_dr_00(i)
1139 0 : d_Pr_dr_in(i) = d_Pr_dVol(i)*dVol_dr_in(i)
1140 :
1141 : !test_partials = (k == s% solver_test_partials_k)
1142 0 : test_partials = .false.
1143 :
1144 : if (test_partials) then
1145 : s% solver_test_partials_val = s% Prad(k)
1146 : s% solver_test_partials_var = i_var_er
1147 : s% solver_test_partials_dval_dx = d_Pr_der(i)
1148 : write(*,*) 'calc_Prad', s% solver_test_partials_var
1149 : write(*,2) 'erad Prad f_Edd', k, s% erad(k), s% Prad(k), s% f_Edd(k)
1150 : end if
1151 :
1152 0 : end subroutine calc_Prad
1153 :
1154 :
1155 0 : subroutine calc_Hp_face(s,i)
1156 : type (star_info), pointer :: s
1157 : integer, intent(in) :: I
1158 0 : real(dp) :: POM
1159 : integer :: k
1160 : logical :: test_partials
1161 : include 'formats'
1162 0 : k = NZN+1-i
1163 0 : if (I < NZN) then
1164 :
1165 0 : POM = (s% r(k)**2)/(2.d0*s% cgrav(k)*s% m(k))
1166 : s% Hp_face(k) = POM*( &
1167 : (s% Pgas(k) + s% Prad(k))*s% Vol(k) &
1168 0 : + (s% Pgas(k-1) + s% Prad(k-1))*s% Vol(k-1))
1169 :
1170 : dHp_dVol_00(I) = POM*( &
1171 : + s% Vol(k)*(d_Pg_dVol(i) + d_Pr_dVol(i)) &
1172 0 : + s% Pgas(k) + s% Prad(k))
1173 : dHp_dVol_out(I) = POM*( &
1174 : + s% Vol(k-1)*(d_Pg_dVol(i+1) + d_Pr_dVol(i+1)) &
1175 0 : + s% Pgas(k-1) + s% Prad(k-1))
1176 :
1177 : dHp_dr_in(I) = POM*( &
1178 : (s% Pgas(k) + s% Prad(k))*dVol_dr_in(I) &
1179 0 : + (d_Pg_dr_in(I) + d_Pr_dr_in(I))*s% Vol(k)) !
1180 : dHp_dr_00(I) = 2.d0*s% Hp_face(k)/s% r(k) + POM*( &
1181 : (s% Pgas(k) + s% Prad(k))*dVol_dr_00(I) &
1182 : + (d_Pg_dr_00(I) + d_Pr_dr_00(I))*s% Vol(k) &
1183 : + (s% Pgas(k-1) + s% Prad(k-1))*dVol_dr_in(i+1) &
1184 0 : + (d_Pg_dr_in(i+1) + d_Pr_dr_in(i+1))*s% Vol(k-1)) !
1185 : dHp_dr_out(I) = POM*( &
1186 : (s% Pgas(k-1) + s% Prad(k-1))*dVol_dr_00(i+1) &
1187 0 : + (d_Pg_dr_00(i+1) + d_Pr_dr_00(i+1))*s% Vol(k-1)) !
1188 :
1189 0 : dHp_dT_00(I) = POM*s% Vol(k)*d_Pg_dT(I) !
1190 0 : dHp_dT_out(I) = POM*s% Vol(k-1)*d_Pg_dT(I+1) !
1191 :
1192 0 : dHp_der_00(I) = POM*s% Vol(k)*d_Pr_der(I) !
1193 0 : dHp_der_out(I) = POM*s% Vol(k-1)*d_Pr_der(I+1) !
1194 :
1195 : else ! surface
1196 :
1197 0 : POM = (s% r(k)**2)/(s% cgrav(k)*s% M(k))
1198 0 : s% Hp_face(k) = POM*(s% Pgas(k) + s% Prad(k))*s% Vol(k)
1199 :
1200 : dHp_dVol_00(I) = POM*( &
1201 : + s% Vol(k)*(d_Pg_dVol(i) + d_Pr_dVol(i)) &
1202 0 : + s% Pgas(k) + s% Prad(k))
1203 0 : dHp_dVol_out(I) = 0d0
1204 :
1205 : dHp_dr_in(i) = POM*( &
1206 : (s% Pgas(k) + s% Prad(k))*dVol_dr_in(i) &
1207 0 : + (d_Pg_dr_in(i) + d_Pr_dr_in(i))*s% Vol(k))
1208 : dHp_dr_00(i) = 2.d0*s% Hp_face(k)/s% r(k) + POM*( &
1209 : (s% Pgas(k) + s% Prad(k))*dVol_dr_00(i) &
1210 0 : + (d_Pg_dr_00(i) + d_Pr_dr_00(i))*s% Vol(k))
1211 0 : dHp_dr_out(i) = 0.d0
1212 0 : dHp_dT_00(i) = POM*s% Vol(k)*d_Pg_dT(i)
1213 0 : dHp_dT_out(i) = 0.d0
1214 0 : dHp_der_00(i) = POM*s% Vol(k)*d_Pr_der(i)
1215 0 : dHp_der_out(i) = 0.d0
1216 :
1217 : end if
1218 :
1219 : !test_partials = (k-1 == s% solver_test_partials_k)
1220 0 : test_partials = .false.
1221 : if (test_partials) then
1222 : s% solver_test_partials_val = s% Hp_face(k)
1223 : s% solver_test_partials_var = i_var_Vol
1224 : s% solver_test_partials_dval_dx = dHp_dVol_out(I)
1225 : write(*,*) 'calc_Hp_face', s% solver_test_partials_var
1226 : end if
1227 :
1228 0 : end subroutine calc_Hp_face
1229 :
1230 :
1231 0 : subroutine calc_Y_face(s,i)
1232 : type (star_info), pointer :: s
1233 : integer, intent(in) :: I
1234 0 : real(dp) :: POM, POM2, &
1235 0 : Y1, d_Y1_dr_00, d_Y1_dr_in, d_Y1_dr_out, &
1236 0 : d_Y1_dVol_00, d_Y1_dVol_out, &
1237 0 : d_Y1_dT_00, d_Y1_dT_out, &
1238 0 : d_Y1_der_00, d_Y1_der_out, &
1239 0 : Y2, d_Y2_dr_00, d_Y2_dr_in, d_Y2_dr_out, &
1240 0 : d_Y2_dVol_00, d_Y2_dVol_out, &
1241 0 : d_Y2_dT_00, d_Y2_dT_out, &
1242 0 : d_Y2_der_00, d_Y2_der_out
1243 : integer :: k
1244 : logical :: test_partials
1245 : include 'formats'
1246 0 : k = NZN+1-i
1247 0 : if (i < NZN .and. ALFA /= 0d0) then
1248 0 : POM = 0.5d0*(s% QQ(k)/s% Cp(k) + s% QQ(k-1)/s% Cp(k-1))
1249 0 : POM2 = 0.5d0*((s% Pgas(k-1) + s% Prad(k-1)) - (s% Pgas(k) + s% Prad(k)))
1250 : Y1 = &
1251 : 0.5d0*(s% QQ(k)/s% Cp(k)+ s% QQ(k-1)/s% Cp(k-1))* &
1252 : ((s% Pgas(k-1) + s% Prad(k-1)) - (s% Pgas(k) + s% Prad(k))) &
1253 0 : - (s% lnT(k-1) - s% lnT(k))
1254 :
1255 : d_Y1_dVol_00 = &
1256 : - POM*(d_Pg_dVol(i) + d_Pr_dVol(i)) &
1257 0 : + POM2*(dQQ_dVol(i) - s% QQ(k)*dCp_dVol(i)/s% Cp(k))/s% Cp(k)
1258 : d_Y1_dVol_out = &
1259 : + POM*(d_Pg_dVol(i+1) + d_Pr_dVol(i+1)) &
1260 0 : + POM2*(dQQ_dVol(i+1) - s% QQ(k-1)*dCp_dVol(i+1)/s% Cp(k-1))/s% Cp(k-1)
1261 :
1262 : d_Y1_dr_00 = &
1263 : POM2*( &
1264 : (dQQ_dr_00(I) - s% QQ(k)/s% Cp(k)*dCp_dr_00(I))/s% Cp(k)&
1265 : +(dQQ_dr_in(i+1) - s% QQ(k-1)/s% Cp(k-1)*dCp_dr_in(i+1))/s% Cp(k-1)) &
1266 0 : + POM*(d_Pg_dr_in(i+1) - d_Pg_dr_00(I) + d_Pr_dr_in(i+1) - d_Pr_dr_00(I))
1267 :
1268 : d_Y1_dr_in = &
1269 : POM2*(dQQ_dr_in(I) - s% QQ(k)/s% Cp(k)*dCp_dr_in(I))/s% Cp(k) &
1270 0 : + POM*(- d_Pg_dr_in(I) - d_Pr_dr_in(I))
1271 :
1272 : d_Y1_dr_out = &
1273 : POM2*(dQQ_dr_00(i+1) - s% QQ(k-1)/s% Cp(k-1)*dCp_dr_00(i+1))/s% Cp(k-1) &
1274 0 : + POM*(d_Pg_dr_00(i+1) + d_Pr_dr_00(i+1))
1275 :
1276 : d_Y1_dT_00 = &
1277 : POM2*(dQQ_dT(I) - s% QQ(k)/s% Cp(k)*dCp_dT(I))/s% Cp(k) &
1278 : - POM*d_Pg_dT(I) &
1279 0 : + 1.d0/s% T(k)
1280 :
1281 : d_Y1_dT_out = &
1282 : POM2*(dQQ_dT(i+1) - s% QQ(k-1)/s% Cp(k-1)*dCp_dT(i+1))/s% Cp(k-1) &
1283 : + POM*d_Pg_dT(I+1) &
1284 0 : - 1.d0/s% T(k-1)
1285 :
1286 0 : d_Y1_der_00 = -POM*d_Pr_der(I)
1287 :
1288 0 : d_Y1_der_out = POM*d_Pr_der(I+1)
1289 :
1290 0 : POM = 2.d0/(s% Vol(k) + s% Vol(k-1))
1291 0 : POM2 = 8.d0*PI*(s% r(k)**2)/s% dm_bar(k)*s% Hp_face(k)
1292 0 : Y2 = 4.d0*PI*(s% r(k)**2)*s% Hp_face(k)*POM/s% dm_bar(k)
1293 :
1294 : d_Y2_dVol_00 = &
1295 : + Y2/s% Hp_face(k)*dHp_dVol_00(i) &
1296 0 : - POM2/(s% Vol(k) + s% Vol(k-1))**2
1297 : d_Y2_dVol_out = &
1298 : + Y2/s% Hp_face(k)*dHp_dVol_out(i) &
1299 0 : - POM2/(s% Vol(k) + s% Vol(k-1))**2
1300 :
1301 : d_Y2_dr_00 = 2.d0*Y2/s% r(k) &
1302 : + Y2/s% Hp_face(k)*dHp_dr_00(I) &
1303 0 : - POM2/(s% Vol(k) + s% Vol(k-1))**2*(dVol_dr_00(I) + dVol_dr_in(i+1))
1304 :
1305 : d_Y2_dr_in = - POM2/(s% Vol(k) &
1306 : + s% Vol(k-1))**2*dVol_dr_in(I) &
1307 0 : + Y2/s% Hp_face(k)*dHp_dr_in(I)
1308 :
1309 : d_Y2_dr_out = - POM2/(s% Vol(k) &
1310 : + s% Vol(k-1))**2*dVol_dr_00(i+1) &
1311 0 : + Y2/s% Hp_face(k)*dHp_dr_out(I)
1312 :
1313 0 : d_Y2_dT_00 = Y2/s% Hp_face(k)*dHp_dT_00(I)
1314 :
1315 0 : d_Y2_dT_out = Y2/s% Hp_face(k)*dHp_dT_out(I)
1316 :
1317 0 : d_Y2_der_00 = Y2/s% Hp_face(k)*dHp_der_00(I)
1318 :
1319 0 : d_Y2_der_out = Y2/s% Hp_face(k)*dHp_der_out(I)
1320 :
1321 0 : s% Y_face(k) = Y1*Y2
1322 :
1323 0 : if (k==-35 .and. iter == 1) then
1324 0 : write(*,3) 'RSP Y_face Y1 Y2', k, s% solver_iter, s% Y_face(k), Y1, Y2
1325 0 : write(*,3) 'Peos', k, s% solver_iter, s% Pgas(k) + s% Prad(k)
1326 0 : write(*,3) 'Peos', k-1, s% solver_iter, s% Pgas(k-1) + s% Prad(k-1)
1327 0 : write(*,3) 'QQ', k, s% solver_iter, s% QQ(k)
1328 0 : write(*,3) 'QQ', k-1, s% solver_iter, s% QQ(k-1)
1329 0 : write(*,3) 'Cp', k, s% solver_iter, s% Cp(k)
1330 0 : write(*,3) 'Cp', k-1, s% solver_iter, s% Cp(k-1)
1331 0 : write(*,3) 'lgT', k, s% solver_iter, s% lnT(k)/ln10
1332 0 : write(*,3) 'lgT', k-1, s% solver_iter, s% lnT(k-1)/ln10
1333 0 : write(*,3) 'lgd', k, s% solver_iter, log(1d0/s% Vol(k))/ln10
1334 0 : write(*,3) 'lgd', k-1, s% solver_iter, log(1d0/s% Vol(k-1))/ln10
1335 : end if
1336 :
1337 0 : dY_dr_00(I) = Y1*d_Y2_dr_00 + Y2*d_Y1_dr_00 !
1338 0 : dY_dr_in(I) = Y1*d_Y2_dr_in + Y2*d_Y1_dr_in !
1339 0 : dY_dr_out(I) = Y1*d_Y2_dr_out + Y2*d_Y1_dr_out !
1340 0 : dY_dVol_00(I) = Y1*d_Y2_dVol_00 + Y2*d_Y1_dVol_00 !
1341 0 : dY_dVol_out(I) = Y1*d_Y2_dVol_out + Y2*d_Y1_dVol_out !
1342 0 : dY_dT_00(I) = Y1*d_Y2_dT_00 + Y2*d_Y1_dT_00 !
1343 0 : dY_dT_out(I) = Y1*d_Y2_dT_out + Y2*d_Y1_dT_out !
1344 0 : dY_der_00(I) = Y1*d_Y2_der_00 + Y2*d_Y1_der_00 !
1345 0 : dY_der_out(I) = Y1*d_Y2_der_out + Y2*d_Y1_der_out !
1346 :
1347 : if (call_is_bad) then
1348 0 : if (is_bad(s% Y_face(k))) then
1349 : !$omp critical (rsp_step_2)
1350 : write(*,2) 's% Y_face(k)', k, s% Y_face(k)
1351 : write(*,2) 'Y1', k, Y1
1352 : write(*,2) 'Y2', k, Y2
1353 : write(*,2) 's% QQ(k)', k, s% QQ(k)
1354 : write(*,2) 's% Cp(k)', k, s% Cp(k)
1355 : write(*,2) 's% Pgas(k)', k, s% Pgas(k)
1356 : write(*,2) 's% Prad(k)', k, s% Prad(k)
1357 : write(*,2) 's% T(k)', k, s% T(k)
1358 : call mesa_error(__FILE__,__LINE__,'calc_Y_face')
1359 : !$omp end critical (rsp_step_2)
1360 : end if
1361 : end if
1362 :
1363 : else
1364 0 : s% Y_face(k) = 0
1365 0 : dY_dr_00(I) = 0
1366 0 : dY_dr_in(I) = 0
1367 0 : dY_dr_out(I) = 0
1368 0 : dY_dVol_00(I) = 0
1369 0 : dY_dVol_out(I) = 0
1370 0 : dY_dT_00(I) = 0
1371 0 : dY_dT_out(I) = 0
1372 0 : dY_der_00(I) = 0
1373 0 : dY_der_out(I) = 0
1374 : end if
1375 :
1376 : !test_partials = (k == s% solver_test_partials_k)
1377 0 : test_partials = .false.
1378 : if (test_partials) then
1379 : s% solver_test_partials_val = s% Y_face(k)
1380 : s% solver_test_partials_var = i_var_Vol
1381 : s% solver_test_partials_dval_dx = dY_dVol_00(I)
1382 : write(*,*) 'calc_Y_face', s% solver_test_partials_var
1383 : end if
1384 :
1385 0 : end subroutine calc_Y_face
1386 :
1387 :
1388 0 : subroutine calc_PII_face(s,i)
1389 : type (star_info), pointer :: s
1390 : integer, intent(in) :: I
1391 0 : real(dp) :: POM, POM2
1392 : integer :: k
1393 : logical :: test_partials
1394 : include 'formats'
1395 0 : k = NZN+1-i
1396 0 : if (k == 1 .or. k == s% nz .or. ALFA == 0d0) then
1397 0 : s% PII(k) = 0
1398 0 : dPII_dr_00(I) = 0
1399 0 : dPII_dr_in(I) = 0
1400 0 : dPII_dr_out(I) = 0
1401 0 : dPII_dVol_00(I) = 0
1402 0 : dPII_dVol_out(I) = 0
1403 0 : dPII_dT_00(I) = 0
1404 0 : dPII_dT_out(I) = 0
1405 0 : dPII_der_00(I) = 0
1406 0 : dPII_der_out(I) = 0
1407 : else
1408 0 : POM = ALFAS*ALFA
1409 0 : POM2 = 0.5d0*(s% Cp(k) + s% Cp(k-1))
1410 0 : s% PII(k) = POM*POM2*s% Y_face(k)
1411 :
1412 : dPII_dVol_00(I) = &
1413 0 : POM*(POM2*dY_dVol_00(I) + s% Y_face(k)*0.5d0*dCp_dVol(I))
1414 : dPII_dVol_out(I) = &
1415 0 : POM*(POM2*dY_dVol_out(I) + s% Y_face(k)*0.5d0*dCp_dVol(I+1))
1416 : dPII_dr_in(I) = & !
1417 0 : POM*(POM2*dY_dr_in(I) + s% Y_face(k)*0.5d0*dCp_dr_in(I))
1418 : dPII_dr_00(I) = & !
1419 0 : POM*(POM2*dY_dr_00(I) + s% Y_face(k)*0.5d0*(dCp_dr_00(I) + dCp_dr_in(i+1)))
1420 : dPII_dr_out(I) = & !
1421 0 : POM*(POM2*dY_dr_out(I) + s% Y_face(k)*0.5d0*dCp_dr_00(i+1))
1422 : dPII_dT_00(I) = & !
1423 0 : POM*(POM2*dY_dT_00(I) + s% Y_face(k)*0.5d0*dCp_dT(I))
1424 : dPII_dT_out(I) = & !
1425 0 : POM*(POM2*dY_dT_out(I) + s% Y_face(k)*0.5d0*dCp_dT(i+1))
1426 0 : dPII_der_00(I) = POM*POM2*dY_der_00(I) !
1427 0 : dPII_der_out(I) = POM*POM2*dY_der_out(I) !
1428 : end if
1429 :
1430 : !test_partials = (k-1 == s% solver_test_partials_k)
1431 0 : test_partials = .false.
1432 : if (test_partials) then
1433 : s% solver_test_partials_val = s% PII(k)
1434 : s% solver_test_partials_var = i_var_Vol
1435 : s% solver_test_partials_dval_dx = dPII_dVol_out(I)
1436 : write(*,*) 'calc_PII_face', s% solver_test_partials_var
1437 : end if
1438 :
1439 0 : end subroutine calc_PII_face
1440 :
1441 :
1442 0 : subroutine calc_Pvsc(s,i)
1443 : type (star_info), pointer :: s
1444 : integer, intent(in) :: I
1445 : real(dp) :: &
1446 0 : dv, dv1, P, dP_dT, dP_der, dP_dr_in, dP_dr_00, V, sqrt_PV, &
1447 0 : d_PV_dT, d_PV_der, d_PV_dr_in, d_PV_dr_00, &
1448 0 : d_sqrt_PV_dT, d_sqrt_PV_der, d_sqrt_PV_dr_in, d_sqrt_PV_dr_00, &
1449 0 : d_dv_dT, d_dv_der, d_dv_dr_in, d_dv_dr_00, &
1450 0 : dP_dVol, d_PV_dVol, d_sqrt_PV_dVol, d_dv_dVol
1451 : integer :: k
1452 : logical :: test_partials
1453 : include 'formats'
1454 0 : k = NZN+1-i
1455 0 : if (CQ == 0d0) then
1456 0 : s% Pvsc(k) = 0d0
1457 0 : d_Pvsc_dVol(i) = 0d0
1458 0 : d_Pvsc_dT(i) = 0d0
1459 0 : d_Pvsc_der(i) = 0d0
1460 0 : d_Pvsc_dr_in(i) = 0d0
1461 0 : d_Pvsc_dr_00(i) = 0d0
1462 0 : return
1463 : end if
1464 0 : if (I > 1) then
1465 0 : dv = s% v(k+1) - s% v(k)
1466 : else
1467 0 : dv = s% v_center - s% v(NZN)
1468 : end if
1469 0 : P = s% Pgas(k) + s% Prad(k)
1470 0 : dP_dT = d_Pg_dT(I)
1471 0 : dP_dVol = d_Pg_dVol(I) + d_Pr_dVol(I)
1472 0 : dP_der = d_Pr_der(I)
1473 0 : dP_dr_in = d_Pg_dr_in(I) + d_Pr_dr_in(I)
1474 0 : dP_dr_00 = d_Pg_dr_00(I) + d_Pr_dr_00(I)
1475 0 : V = s% Vol(k)
1476 0 : sqrt_PV = sqrt(P*V)
1477 0 : if (dv <= ZSH*sqrt_PV) then
1478 0 : s% Pvsc(k) = 0d0
1479 0 : d_Pvsc_dVol(i) = 0d0
1480 0 : d_Pvsc_dT(i) = 0d0
1481 0 : d_Pvsc_der(i) = 0d0
1482 0 : d_Pvsc_dr_in(i) = 0d0
1483 0 : d_Pvsc_dr_00(i) = 0d0
1484 0 : return
1485 : end if
1486 :
1487 0 : d_PV_dT = dP_dT*V
1488 0 : d_PV_dVol = P + dP_dVol*V
1489 0 : d_PV_der = dP_der*V
1490 0 : d_PV_dr_in = dP_dr_in*V + P*dVol_dr_in(I)
1491 0 : d_PV_dr_00 = dP_dr_00*V + P*dVol_dr_00(I)
1492 :
1493 0 : d_sqrt_PV_dT = 0.5d0*d_PV_dT/sqrt_PV
1494 0 : d_sqrt_PV_dVol = 0.5d0*d_PV_dVol/sqrt_PV
1495 0 : d_sqrt_PV_der = 0.5d0*d_PV_der/sqrt_PV
1496 0 : d_sqrt_PV_dr_in = 0.5d0*d_PV_dr_in/sqrt_PV
1497 0 : d_sqrt_PV_dr_00 = 0.5d0*d_PV_dr_00/sqrt_PV
1498 :
1499 0 : dv1 = dv
1500 0 : dv = dv - ZSH*sqrt_PV
1501 :
1502 0 : d_dv_dT = - ZSH*d_sqrt_PV_dT
1503 0 : d_dv_dVol = - ZSH*d_sqrt_PV_dVol
1504 0 : d_dv_der = - ZSH*d_sqrt_PV_der
1505 :
1506 0 : d_dv_dr_in = 2d0/s% dt - ZSH*d_sqrt_PV_dr_in ! not used if I == 1
1507 0 : d_dv_dr_00 = -2d0/s% dt - ZSH*d_sqrt_PV_dr_00
1508 :
1509 : ! Pvsc = CQ*P*(dv1/sqrt_PV - cut)^2 eqn 3.60
1510 : ! = CQ*P*((dv1 - cut*sqrt_PV)/sqrt_PV)^2
1511 : ! = CQ*P/(P*V)*dv^2
1512 : ! = CQ/V*dv^2
1513 :
1514 0 : s% Pvsc(k) = CQ/V*dv**2
1515 0 : d_Pvsc_dVol(i) = -s% Pvsc(k)/V + 2d0*d_dv_dVol*CQ/V*dv
1516 0 : d_Pvsc_dT(i) = CQ/V*2d0*dv*d_dv_dT
1517 0 : d_Pvsc_der(i) = CQ/V*2d0*dv*d_dv_der
1518 0 : d_Pvsc_dr_in(i) = CQ/V*2d0*dv*d_dv_dr_in - CQ*dv**2*dVol_dr_in(I)/V**2
1519 0 : d_Pvsc_dr_00(i) = CQ/V*2d0*dv*d_dv_dr_00 - CQ*dv**2*dVol_dr_00(I)/V**2
1520 :
1521 : !test_partials = (k == s% solver_test_partials_k)
1522 0 : test_partials = .false.
1523 : if (test_partials) then
1524 : s% solver_test_partials_val = s% Pvsc(k)
1525 : s% solver_test_partials_var = i_var_T
1526 : s% solver_test_partials_dval_dx = d_Pvsc_dT(i)
1527 : write(*,*) 'calc_Pvsc', s% solver_test_partials_var
1528 : end if
1529 : end subroutine calc_Pvsc
1530 :
1531 :
1532 0 : subroutine check_omega(s,i)
1533 : type (star_info), pointer :: s
1534 : integer, intent(in) :: i
1535 0 : real(dp) :: SOURS, DAMPS, DAMPRS, DELTA, SOL, POM, POM2, POM3, POM4, w_start
1536 : integer :: k
1537 : include 'formats'
1538 0 : if (I > IBOTOM .and. I < NZN .and. ALFA /= 0d0) then
1539 : ! JAK OKRESLIC OMEGA DLA PIERWSZEJ ITERACJI
1540 0 : k = NZN+1-i
1541 0 : if (s% RSP_w(k) > EFL0) return
1542 0 : w_start = s% RSP_w(k)
1543 0 : POM = (s% PII(k)/s% Hp_face(k) + s% PII(k+1)/s% Hp_face(k+1))*0.5d0
1544 0 : POM2 = s% T(k)*(s% Pgas(k) + s% Prad(k))*s% QQ(k)/s% Cp(k)
1545 0 : SOURS = POM*POM2
1546 0 : DAMPS = (CEDE/ALFA)/((s% Hp_face(k) + s% Hp_face(k+1))*0.5d0)
1547 0 : POM3 = (GAMMAR**2)/(ALFA**2)*4.d0*SIG
1548 0 : POM4 = (s% T(k)**3)*(s% Vol(k)**2)/(s% Cp(k)*s% opacity(k))
1549 0 : DAMPRS = POM3*POM4/((s% Hp_face(k)**2 + s% Hp_face(k+1)**2)*0.5d0)
1550 0 : DELTA = DAMPRS**2 + 4.d0*DAMPS*SOURS
1551 : if (k==-35) then
1552 : write(*,2) 'del', k, DELTA
1553 : write(*,2) 'DAMPR', k, DAMPRS
1554 : write(*,2) 'DAMP', k, DAMPS
1555 : write(*,2) 'SOURCE', k, SOURS
1556 : write(*,2) 'POM', k, POM
1557 : write(*,2) 'POM2', k, POM2
1558 : write(*,2) 's% Hp_face(k)', k, s% Hp_face(k)
1559 : write(*,2) 's% Hp_face(k+1)', k+1, s% Hp_face(k+1)
1560 : write(*,2) 's% PII(k)', k, s% PII(k)
1561 : write(*,2) 's% PII(k+1)', k+1, s% PII(k+1)
1562 : write(*,2) 's% Y_face(k)', k, s% Y_face(k)
1563 : write(*,2) 's% Y_face(k+1)', k+1, s% Y_face(k+1)
1564 : end if
1565 0 : if (DELTA >= 0.d0) SOL = ( - DAMPRS + sqrt(DELTA))/(2.d0*DAMPS)
1566 0 : if (DELTA < 0.d0) SOL = - 99.99d0
1567 0 : if (SOL >= 0.d0) SOL = SOL**2
1568 0 : if (SOL > 0.d0) then
1569 0 : s% RSP_w(k) = sqrt(SOL)
1570 0 : if (s% RSP_report_check_omega_changes) &
1571 0 : write(*,3) 'RSP_w modified initial guess vs initial', k, s% model_number, &
1572 0 : s% RSP_w(k), w_start
1573 : end if
1574 : end if
1575 : end subroutine check_omega
1576 :
1577 :
1578 0 : subroutine calc_Pturb(s,i) ! TURBULENT PRESSURE (ZONE)
1579 : type (star_info), pointer :: s
1580 : integer, intent(in) :: I
1581 0 : real(dp) :: TEM1, Vol
1582 : integer :: k
1583 : logical :: test_partials
1584 : include 'formats'
1585 0 : k = NZN+1-i
1586 : if (ALFA == 0d0 .or. ALFAP == 0d0 .or. &
1587 0 : I <= IBOTOM .or. I >= NZN) then
1588 0 : s% Ptrb(k) = 0.d0
1589 0 : dPtrb_dVol_00(I) = 0.d0
1590 0 : dPtrb_dw_00(I) = 0.d0
1591 0 : dPtrb_dr_00(I) = 0.d0
1592 0 : dPtrb_dr_in(I) = 0.d0
1593 : else
1594 0 : Vol = s% Vol(k)
1595 0 : s% Ptrb(k) = ALFAP*s% RSP_w(k)**2/Vol
1596 0 : dPtrb_dVol_00(I) = -s% Ptrb(k)/Vol
1597 0 : dPtrb_dw_00(I) = 2.d0*ALFAP*s% RSP_w(k)/Vol
1598 0 : TEM1 = - ALFAP*s% RSP_w(k)**2/Vol**2
1599 0 : dPtrb_dr_00(I) = TEM1*dVol_dr_00(I)
1600 0 : dPtrb_dr_in(I) = TEM1*dVol_dr_in(I)
1601 : end if
1602 : !test_partials = (k == s% solver_test_partials_k)
1603 0 : test_partials = .false.
1604 : if (test_partials) then
1605 : s% solver_test_partials_val = 0 ! residual
1606 : s% solver_test_partials_var = i_var_R
1607 : s% solver_test_partials_dval_dx = 0 ! d_residual_dr_00
1608 : write(*,*) 'calc_Pturb', s% solver_test_partials_var
1609 : end if
1610 0 : end subroutine calc_Pturb
1611 :
1612 :
1613 0 : subroutine calc_Chi(s,i) ! eddy viscosity (Kuhfuss 1986)
1614 : type (star_info), pointer :: s
1615 : integer, intent(in) :: I
1616 0 : real(dp) :: POM, POM1, POM2, POM3, POM4, &
1617 0 : POMT1, POMT2, POMT3, POMT4, Vol
1618 : integer :: k
1619 : logical :: test_partials
1620 : include 'formats'
1621 0 : k = NZN+1-i
1622 0 : if (ALFA == 0d0 .or. I <= IBOTOM .or. I >= NZN) then
1623 0 : s% Chi(k) = 0
1624 0 : dChi_dT_out(I) = 0
1625 0 : dChi_dT_00(I) = 0
1626 0 : dChi_dT_in(I) = 0
1627 0 : dChi_der_out(I) = 0
1628 0 : dChi_der_00(I) = 0
1629 0 : dChi_der_in(I) = 0
1630 0 : dChi_dw_00(I) = 0
1631 0 : dChi_dr_out(I) = 0
1632 0 : dChi_dr_00(I) = 0
1633 0 : dChi_dr_in(I) = 0
1634 0 : dChi_dr_in2(I) = 0
1635 : !s% profile_extra(k,2) = 0
1636 : !s% profile_extra(k,3) = 0
1637 : !s% profile_extra(k,4) = 0
1638 : else
1639 0 : POM = (16.d0/3.d0)*PI*ALFA*ALFAM/s% dm(k)
1640 0 : Vol = s% Vol(k)
1641 0 : POM1 = s% RSP_w(k)/Vol**2
1642 0 : POM2 = 0.5d0*(s% r(k)**6 + s% r(k+1)**6)
1643 0 : POM4 = 0.5d0*(s% Hp_face(k) + s% Hp_face(k+1))
1644 0 : POM3 = s% v(k)/s% r(k) - s% v(k+1)/s% r(k+1)
1645 0 : POMT3 = POM*POM1*POM2*POM4
1646 0 : POMT1 = POM*POM2*POM3*POM4
1647 0 : POMT2 = POM*POM1*POM3*POM4
1648 0 : POMT4 = POM*POM1*POM2*POM3
1649 :
1650 0 : s% Chi(k) = POMT1*POM1
1651 :
1652 : if (call_is_bad) then
1653 : if (is_bad(s% Chi(k))) then
1654 : !$omp critical (rsp_step_3)
1655 : write(*,2) 'POM', k, POM
1656 : write(*,2) 'POM1', k, POM1
1657 : write(*,2) 'POM2', k, POM2
1658 : write(*,2) 'POM3', k, POM3
1659 : write(*,2) 'POM4', k, POM4
1660 : write(*,2) 's% RSP_w(k)', k, s% RSP_w(k)
1661 : write(*,2) 's% Volk)', k, s% Vol(k)
1662 : write(*,2) 's% r(k)', k, s% r(k)
1663 : write(*,2) 's% r(k+1)', k+1, s% r(k+1)
1664 : write(*,2) 's% v(k)', k, s% v(k)
1665 : write(*,2) 's% v(k+1)', k+1, s% v(k+1)
1666 : write(*,2) 's% Vol(k+1)', k+1, s% Vol(k+1)
1667 : write(*,2) 's% rho(k+1)', k+1, s% rho(k+1)
1668 : call mesa_error(__FILE__,__LINE__,'calc_Chi')
1669 : !$omp end critical (rsp_step_3)
1670 : end if
1671 : end if
1672 :
1673 : dChi_dVol_out(I) = &
1674 0 : POMT4*0.5d0*dHp_dVol_out(I)
1675 : dChi_dVol_00(I) = &
1676 : + POMT4*0.5d0*(dHp_dVol_00(I) + dHp_dVol_out(i-1)) &
1677 0 : - 2d0*s% Chi(k)/Vol
1678 : dChi_dVol_in(I) = &
1679 0 : POMT4*0.5d0*dHp_dVol_00(i-1)
1680 :
1681 0 : dChi_dT_out(I) = POMT4*0.5d0*dHp_dT_out(I) !
1682 0 : dChi_dT_00(I) = POMT4*0.5d0*(dHp_dT_00(I) + dHp_dT_out(i-1)) !
1683 0 : dChi_dT_in(I) = POMT4*0.5d0*dHp_dT_00(i-1) !
1684 :
1685 0 : dChi_der_out(I) = POMT4*0.5d0*(dHp_der_out(I)) !
1686 0 : dChi_der_00(I) = POMT4*0.5d0*(dHp_der_00(I) + dHp_der_out(i-1)) !
1687 0 : dChi_der_in(I) = POMT4*0.5d0*(dHp_der_00(i-1)) !
1688 :
1689 0 : dChi_dw_00(I) = POMT1/Vol**2 !
1690 :
1691 0 : dChi_dr_out(I) = POMT4*0.5d0*(dHp_dr_out(I)) !
1692 : dChi_dr_00(I) = &
1693 : - 2.d0*s% Chi(k)/Vol*dVol_dr_00(I) & !
1694 : + POMT2*3.d0*s% r(k)**5 &
1695 : + POMT3*(2.d0/s% dt/s% r(k) - s% v(k)/s% r(k)**2) &
1696 0 : + POMT4*0.5d0*(dHp_dr_00(I) + dHp_dr_out(i-1))
1697 : dChi_dr_in(I) = &
1698 : - 2.d0*s% Chi(k)/Vol*dVol_dr_in(I) & !
1699 : + POMT2*3.d0*s% r(k+1)**5 &
1700 : - POMT3*(2.d0/s% dt/s% r(k+1) - s% v(k+1)/s% r(k+1)**2) &
1701 0 : + POMT4*0.5d0*(dHp_dr_in(I) + dHp_dr_00(i-1))
1702 0 : dChi_dr_in2(I) = POMT4*0.5d0*dHp_dr_in(i-1) !
1703 :
1704 : if (call_is_bad) then
1705 : if (is_bad(dChi_dr_in(I))) then
1706 : !$omp critical (rsp_step_4)
1707 : write(*,2) 's% Chi(k)', k, s% Chi(k)
1708 : write(*,2) 'Vol', k, Vol
1709 : write(*,2) 'POMT2', k, POMT2
1710 : write(*,2) 'POMT3', k, POMT3
1711 : write(*,2) 'POMT4', k, POMT4
1712 : write(*,2) 'dVol_dr_in(I)', I, dVol_dr_in(I)
1713 : write(*,2) 'dChi_dr_in(I)', I, dChi_dr_in(I)
1714 : write(*,2) 'dHp_dr_00(i-1)', I-1, dHp_dr_00(i-1)
1715 : write(*,2) 'dHp_dr_in(I)', I, dHp_dr_in(I)
1716 : write(*,2) 's% r(k+1)', k+1, s% r(k+1)
1717 : write(*,2) 's% v(k+1)', k+1, s% v(k+1)
1718 : call mesa_error(__FILE__,__LINE__,'calc_Chi')
1719 : !$omp end critical (rsp_step_4)
1720 : end if
1721 : end if
1722 :
1723 : end if
1724 :
1725 : !if (k==194) then
1726 : ! write(*,2) 'RSP Chi', k, s% Chi(k)
1727 : !end if
1728 :
1729 : !test_partials = (k-1 == s% solver_test_partials_k)
1730 0 : test_partials = .false.
1731 : if (test_partials) then
1732 : s% solver_test_partials_val = s% Chi(k)
1733 : s% solver_test_partials_var = i_var_Vol
1734 : s% solver_test_partials_dval_dx = dChi_dVol_out(I)
1735 : write(*,*) 'calc_Chi', s% solver_test_partials_var
1736 : end if
1737 0 : end subroutine calc_Chi
1738 :
1739 :
1740 0 : subroutine calc_Eq(s,i)
1741 : type (star_info), pointer :: s
1742 : integer, intent(in) :: I
1743 0 : real(dp) :: RRI, RRM, UUI, UUM, POM, POM2
1744 : integer :: k
1745 : logical :: test_partials
1746 : include 'formats'
1747 0 : k = NZN+1-i
1748 0 : if (ALFA == 0d0 .or. I <= IBOTOM .or. I >= NZN) then
1749 0 : s% Eq(k) = 0
1750 0 : dEq_dr_out(I) = 0
1751 0 : dEq_dr_00(I) = 0
1752 0 : dEq_dr_in(I) = 0
1753 0 : dEq_dr_in2(I) = 0
1754 0 : dEq_dVol_out(I) = 0
1755 0 : dEq_dVol_00(I) = 0
1756 0 : dEq_dVol_in(I) = 0
1757 0 : dEq_dT_out(I) = 0
1758 0 : dEq_dT_00(I) = 0
1759 0 : dEq_dT_in(I) = 0
1760 0 : dEq_der_out(I) = 0
1761 0 : dEq_der_00(I) = 0
1762 0 : dEq_der_in(I) = 0
1763 0 : dEq_dw_00(I) = 0
1764 : else
1765 :
1766 0 : RRI = 0.5d0*(s% r(k) + s% r_start(k))
1767 0 : RRM = 0.5d0*(s% r(k+1) + s% r_start(k+1))
1768 0 : UUI = 0.5d0*(s% v(k) + s% v_start(k))
1769 0 : UUM = 0.5d0*(s% v(k+1) + s% v_start(k+1))
1770 :
1771 0 : POM = P4/s% dm(k)*(UUI/RRI - UUM/RRM)
1772 0 : POM2 = P4/s% dm(k)*(THETAU*s% Chi(k) + THETAU1*s% Chi_start(k))
1773 :
1774 0 : s% Eq(k) = POM*(THETAU*s% Chi(k) + THETAU1*s% Chi_start(k))
1775 :
1776 0 : POM = POM*THETAU
1777 :
1778 0 : dEq_dVol_out(I) = POM*THETAU*dChi_dVol_out(i)
1779 0 : dEq_dVol_00(I) = POM*THETAU*dChi_dVol_00(i)
1780 0 : dEq_dVol_in(I) = POM*THETAU*dChi_dVol_in(i)
1781 :
1782 0 : dEq_dr_out(I) = POM*dChi_dr_out(I) !
1783 : dEq_dr_00(I) = POM*dChi_dr_00(I) & !
1784 0 : + POM2*(1.d0/(RRI*s% dt) - 0.5d0*UUI/(RRI**2))
1785 : dEq_dr_in(I) = POM*dChi_dr_in(I) & !
1786 0 : - POM2*(1.d0/(RRM*s% dt) - 0.5d0*UUM/(RRM**2))
1787 0 : dEq_dr_in2(I) = POM*dChi_dr_in2(I) !
1788 0 : dEq_dT_out(I) = POM*dChi_dT_out(I) !
1789 0 : dEq_dT_00(I) = POM*dChi_dT_00(I) !
1790 0 : dEq_dT_in(I) = POM*dChi_dT_in(I) !
1791 0 : dEq_der_out(I) = POM*dChi_der_out(I) !
1792 0 : dEq_der_00(I) = POM*dChi_der_00(I) !
1793 0 : dEq_der_in(I) = POM*dChi_der_in(I) !
1794 0 : dEq_dw_00(I) = POM*dChi_dw_00(I) !
1795 : end if
1796 :
1797 : !test_partials = (k+1 == s% solver_test_partials_k)
1798 0 : test_partials = .false.
1799 : if (test_partials) then
1800 : s% solver_test_partials_val = s% Eq(k)
1801 : s% solver_test_partials_var = i_var_Vol
1802 : s% solver_test_partials_dval_dx = dEq_dVol_in(I)
1803 : write(*,*) 'calc_Eq', s% solver_test_partials_var
1804 : end if
1805 0 : end subroutine calc_Eq
1806 :
1807 :
1808 0 : subroutine calc_source_sink(s,i)
1809 : type (star_info), pointer :: s
1810 : integer, intent(in) :: I
1811 0 : real(dp) :: POM, POM2, POM3, POM4, TEM1, TEMI, TEMM, &
1812 0 : dsrc_dr_in2, dsrc_dr_in, dsrc_dr_00, dsrc_dr_out, &
1813 0 : dsrc_dVol_out, dsrc_dVol_00, dsrc_dVol_in, &
1814 0 : dsrc_dT_out, dsrc_dT_00, dsrc_dT_in, &
1815 0 : dsrc_der_out, dsrc_der_00, dsrc_der_in, &
1816 0 : dsrc_dw_00, &
1817 0 : d_damp_dr_in2, d_damp_dr_in, d_damp_dr_00, d_damp_dr_out, &
1818 0 : d_damp_dVol_out, d_damp_dVol_00, d_damp_dVol_in, &
1819 0 : d_damp_dT_out, d_damp_dT_00, d_damp_dT_in, &
1820 0 : d_damp_der_out, d_damp_der_00, d_damp_der_in, &
1821 0 : d_damp_dw_00, &
1822 0 : d_dampR_dr_in2, d_dampR_dr_in, d_dampR_dr_00, d_dampR_dr_out, &
1823 0 : d_dampR_dVol_out, d_dampR_dVol_00, d_dampR_dVol_in, &
1824 0 : d_dampR_dT_out, d_dampR_dT_00, d_dampR_dT_in, &
1825 0 : d_dampR_der_out, d_dampR_der_00, d_dampR_der_in, &
1826 0 : d_dampR_dw_00, d_POM_dVol_00, d_POM_dVol_in, &
1827 0 : d_QQ_div_Cp_d_Vol_00, d_POM2_dVol_00, QQ_div_Cp, &
1828 0 : POM2a, POM2b, d_POM2a_dVol_00, d_POM2b_dVol_00, &
1829 0 : d_POM4_dVol_00
1830 : integer :: k
1831 : logical :: test_partials
1832 : include 'formats'
1833 0 : k = NZN+1-i
1834 0 : if (ALFA == 0d0 .or. I <= IBOTOM .or. I >= NZN) then
1835 0 : s% SOURCE(k) = 0
1836 0 : s% DAMP(k) = 0
1837 0 : s% DAMPR(k) = 0
1838 0 : s% COUPL(k) = 0
1839 : else
1840 : ! SOURCE TERM
1841 0 : POM = 0.5d0*(s% PII(k)/s% Hp_face(k) + s% PII(k+1)/s% Hp_face(k+1))
1842 0 : QQ_div_Cp = s% QQ(k)/s% Cp(k)
1843 0 : POM2 = s% T(k)*(s% Pgas(k) + s% Prad(k))*QQ_div_Cp
1844 0 : POM3 = s% RSP_w(k)
1845 0 : s% SOURCE(k) = POM*POM2*POM3
1846 :
1847 : ! P*QQ/Cp = grad_ad
1848 0 : if (k==-109) write(*,3) 'w grada PII_00 PII_p1 SOURCE', k, s% solver_iter, &
1849 0 : s% RSP_w(k), (s% Pgas(k) + s% Prad(k))*QQ_div_Cp, s% PII(k), &
1850 0 : s% PII(k+1), s% SOURCE(k)
1851 :
1852 0 : TEM1 = POM2*POM3*0.5d0
1853 0 : TEMI = - s% PII(k)/s% Hp_face(k)**2
1854 0 : TEMM = - s% PII(k+1)/s% Hp_face(k+1)**2
1855 :
1856 : d_POM_dVol_00 = &
1857 : 0.5d0*(dPII_dVol_00(I) - s% PII(k)/s% Hp_face(k)*dHp_dVol_00(I))/s% Hp_face(k) + &
1858 0 : 0.5d0*(dPII_dVol_out(I-1) - s% PII(k+1)/s% Hp_face(k+1)*dHp_dVol_out(I-1))/s% Hp_face(k+1)
1859 : d_POM_dVol_in = 0.5d0*( &
1860 0 : dPII_dVol_00(I-1) - s% PII(k+1)/s% Hp_face(k+1)*dHp_dVol_00(I-1))/s% Hp_face(k+1)
1861 0 : d_QQ_div_Cp_d_Vol_00 = (dQQ_dVol(i) - dCp_dVol(i)*s% QQ(k)/s% Cp(k))/s% Cp(k)
1862 : d_POM2_dVol_00 = s% T(k)*( &
1863 : (d_Pg_dVol(i) + d_Pr_dVol(i))*QQ_div_Cp + &
1864 0 : (s% Pgas(k) + s% Prad(k))*d_QQ_div_Cp_d_Vol_00)
1865 : dsrc_dVol_out = TEM1*( & ! ok
1866 : dPII_dVol_out(I)/s% Hp_face(k) &
1867 0 : + TEMI*dHp_dVol_out(I))
1868 0 : dsrc_dVol_00 = POM3*(d_POM_dVol_00*POM2 + POM*d_POM2_dVol_00)
1869 : dsrc_dVol_in = TEM1*( & ! ok
1870 : dPII_dVol_00(i-1)/s% Hp_face(k+1) &
1871 0 : + TEMM*dHp_dVol_00(i-1))
1872 :
1873 : dsrc_dr_out = TEM1*(dPII_dr_out(I)/s% Hp_face(k) &
1874 0 : + TEMI*dHp_dr_out(I))
1875 : dsrc_dr_00 = TEM1*(dPII_dr_00(I)/s% Hp_face(k) + dPII_dr_out(i-1)/s% Hp_face(k+1) &
1876 0 : + TEMI*dHp_dr_00(I) + TEMM*dHp_dr_out(i-1))
1877 : dsrc_dr_in = TEM1*(dPII_dr_in(I)/s% Hp_face(k) + dPII_dr_00(i-1)/s% Hp_face(k+1) &
1878 0 : + TEMI*dHp_dr_in(I) + TEMM*dHp_dr_00(i-1))
1879 : dsrc_dr_in2 = TEM1*(dPII_dr_in(i-1)/s% Hp_face(k+1) &
1880 0 : + TEMM*dHp_dr_in(i-1))
1881 :
1882 : dsrc_dT_out = TEM1*( &
1883 : dPII_dT_out(I)/s% Hp_face(k) &
1884 0 : + TEMI*dHp_dT_out(I))
1885 : dsrc_dT_00 = TEM1*( &
1886 : dPII_dT_00(I)/s% Hp_face(k) &
1887 : + dPII_dT_out(i-1)/s% Hp_face(k+1) &
1888 : + TEMI*dHp_dT_00(I) &
1889 0 : + TEMM*dHp_dT_out(i-1))
1890 : dsrc_dT_in = TEM1*( &
1891 : dPII_dT_00(i-1)/s% Hp_face(k+1) &
1892 0 : + TEMM*dHp_dT_00(i-1))
1893 :
1894 : dsrc_der_out = TEM1*( &
1895 : dPII_der_out(I)/s% Hp_face(k) &
1896 0 : + TEMI*dHp_der_out(I))
1897 : dsrc_der_00 = TEM1*( &
1898 : dPII_der_00(I)/s% Hp_face(k) &
1899 : + dPII_der_out(i-1)/s% Hp_face(k+1) &
1900 : + TEMI*dHp_der_00(I) &
1901 0 : + TEMM*dHp_der_out(i-1))
1902 : dsrc_der_in = TEM1*( &
1903 : dPII_der_00(i-1)/s% Hp_face(k+1) &
1904 0 : + TEMM*dHp_der_00(i-1))
1905 :
1906 0 : dsrc_dw_00 = POM*POM2
1907 :
1908 0 : POM = POM*POM3
1909 :
1910 : dsrc_dT_00 = dsrc_dT_00 + POM/s% Cp(k)*( &
1911 : (s% Pgas(k) + s% Prad(k))*s% QQ(k) &
1912 : + s% T(k)*s% QQ(k)*d_Pg_dT(I) &
1913 : + s% T(k)*(s% Pgas(k) + s% Prad(k))*dQQ_dT(I) &
1914 0 : - s% T(k)*(s% Pgas(k) + s% Prad(k))*s% QQ(k)/s% Cp(k)*dCp_dT(I))
1915 :
1916 : dsrc_der_00 = dsrc_der_00 &
1917 0 : + POM/s% Cp(k)*s% T(k)*s% QQ(k)*d_Pr_der(I)
1918 :
1919 : dsrc_dr_00 = dsrc_dr_00 + POM*s% T(k)/s% Cp(k)*( &
1920 : s% QQ(k)*(d_Pg_dr_00(I) + d_Pr_dr_00(I)) &
1921 : + (s% Pgas(k) + s% Prad(k))*dQQ_dr_00(I) &
1922 0 : - (s% Pgas(k) + s% Prad(k))*s% QQ(k)/s% Cp(k)*dCp_dr_00(I))
1923 : dsrc_dr_in = dsrc_dr_in + POM*s% T(k)/s% Cp(k)*( &
1924 : s% QQ(k)*(d_Pg_dr_in(I) + d_Pr_dr_in(I)) &
1925 : + (s% Pgas(k) + s% Prad(k))*dQQ_dr_in(I) &
1926 0 : - (s% Pgas(k) + s% Prad(k))*s% QQ(k)/s% Cp(k)*dCp_dr_in(I))
1927 :
1928 : ! DAMP TERM
1929 0 : POM = (CEDE/ALFA)*(s% RSP_w(k)**3 - EFL0**3)
1930 0 : POM2 = 0.5d0*(s% Hp_face(k) + s% Hp_face(k+1))
1931 0 : s% DAMP(k) = POM/POM2
1932 :
1933 0 : TEM1 = - 0.5d0*POM/POM2**2
1934 :
1935 0 : d_damp_dVol_out = TEM1*dHp_dVol_out(I)
1936 0 : d_damp_dVol_00 = TEM1*(dHp_dVol_00(I) + dHp_dVol_out(i-1))
1937 0 : d_damp_dVol_in = TEM1*dHp_dVol_00(i-1)
1938 0 : d_damp_dr_out = TEM1*dHp_dr_out(I)
1939 0 : d_damp_dr_00 = TEM1*(dHp_dr_00(I) + dHp_dr_out(i-1))
1940 0 : d_damp_dr_in = TEM1*(dHp_dr_in(I) + dHp_dr_00(i-1))
1941 0 : d_damp_dr_in2 = TEM1*(dHp_dr_in(i-1))
1942 0 : d_damp_dT_00 = TEM1*(dHp_dT_00(I) + dHp_dT_out(i-1))
1943 0 : d_damp_dT_out = TEM1*dHp_dT_out(I)
1944 0 : d_damp_dT_in = TEM1*dHp_dT_00(i-1)
1945 0 : d_damp_der_00 = TEM1*(dHp_der_00(I) + dHp_der_out(i-1))
1946 0 : d_damp_der_out = TEM1*dHp_der_out(I)
1947 0 : d_damp_der_in = TEM1*dHp_der_00(i-1)
1948 0 : d_damp_dw_00 = 3.0d0*(CEDE/ALFA)/POM2*s% RSP_w(k)**2
1949 :
1950 : ! RADIATIVE DAMP TERM
1951 0 : if (GAMMAR == 0.d0)then
1952 0 : s% DAMPR(k) = 0.d0
1953 0 : d_dampR_dr_out = 0.d0
1954 0 : d_dampR_dr_00 = 0.d0
1955 0 : d_dampR_dr_in = 0.d0
1956 0 : d_dampR_dr_in2 = 0.d0
1957 0 : d_dampR_dVol_out = 0.d0
1958 0 : d_dampR_dVol_in = 0.d0
1959 0 : d_dampR_dVol_00 = 0.d0
1960 0 : d_dampR_dT_out = 0.d0
1961 0 : d_dampR_dT_in = 0.d0
1962 0 : d_dampR_dT_00 = 0.d0
1963 0 : d_dampR_der_out = 0.d0
1964 0 : d_dampR_der_in = 0.d0
1965 0 : d_dampR_der_00 = 0.d0
1966 0 : d_dampR_dw_00 = 0.d0
1967 : else
1968 0 : POM = (GAMMAR**2)/(ALFA**2)*4.d0*SIG
1969 0 : POM2a = s% T(k)**3*s% Vol(k)**2
1970 0 : POM2b = 1d0/(s% Cp(k)*s% opacity(k))
1971 0 : POM2 = POM2a*POM2b
1972 0 : POM3 = s% RSP_w(k)**2
1973 0 : POM4 = 0.5d0*(s% Hp_face(k)**2 + s% Hp_face(k+1)**2)
1974 0 : s% DAMPR(k) = POM*POM2*POM3/POM4
1975 :
1976 0 : TEM1 = - s% DAMPR(k)/POM4
1977 :
1978 0 : d_POM2a_dVol_00 = 2d0*s% T(k)**3*s% Vol(k)
1979 : d_POM2b_dVol_00 = &
1980 0 : -POM2b*(dCp_dVol(i)/s% Cp(k) + dK_dVol(i)/s% opacity(k))
1981 0 : d_POM2_dVol_00 = d_POM2a_dVol_00*POM2b + POM2a*d_POM2b_dVol_00
1982 : d_POM4_dVol_00 = &
1983 : s% Hp_face(k)*dHp_dVol_00(I) + &
1984 0 : s% Hp_face(k+1)*dHp_dVol_out(I-1)
1985 :
1986 0 : d_dampR_dVol_out = TEM1*s% Hp_face(k)*dHp_dVol_out(I)
1987 : d_dampR_dVol_00 = POM*POM3*( &
1988 0 : d_POM2_dVol_00 - POM2*d_POM4_dVol_00/POM4)/POM4
1989 :
1990 0 : d_dampR_dVol_in = TEM1*s% Hp_face(k+1)*dHp_dVol_00(i-1)
1991 :
1992 0 : d_dampR_dr_out = TEM1*s% Hp_face(k)*dHp_dr_out(I)
1993 : d_dampR_dr_00 = TEM1*(s% Hp_face(k)*dHp_dr_00(I) &
1994 0 : + s% Hp_face(k+1)*dHp_dr_out(i-1))
1995 : d_dampR_dr_in = TEM1*(s% Hp_face(k)*dHp_dr_in(I) &
1996 0 : + s% Hp_face(k+1)*dHp_dr_00(i-1))
1997 0 : d_dampR_dr_in2 = TEM1*s% Hp_face(k+1)*dHp_dr_in(i-1)
1998 :
1999 0 : d_dampR_dT_out = TEM1*s% Hp_face(k)*dHp_dT_out(I)
2000 : d_dampR_dT_00 = TEM1*(s% Hp_face(k)*dHp_dT_00(I) &
2001 0 : + s% Hp_face(k+1)*dHp_dT_out(i-1))
2002 0 : d_dampR_dT_in = TEM1*s% Hp_face(k+1)*dHp_dT_00(i-1)
2003 :
2004 0 : d_dampR_der_out = TEM1*s% Hp_face(k)*dHp_der_out(I)
2005 : d_dampR_der_00 = TEM1*(s% Hp_face(k)*dHp_der_00(I) &
2006 0 : + s% Hp_face(k+1)*dHp_der_out(i-1))
2007 0 : d_dampR_der_in = TEM1*s% Hp_face(k+1)*dHp_der_00(i-1)
2008 :
2009 0 : d_dampR_dw_00 = POM*POM2/POM4*2.d0*s% RSP_w(k)
2010 :
2011 0 : TEM1 = POM*POM3/POM4
2012 : d_dampR_dr_00 = d_dampR_dr_00 &
2013 : + TEM1*s% T(k)**3*(2.d0*s% Vol(k)*dVol_dr_00(I) &
2014 : - s% Vol(k)**2*(1.d0/s% Cp(k)*dCp_dr_00(I) &
2015 : + 1.d0/s% opacity(k)*dK_dr_00(I))) &
2016 0 : /(s% Cp(k)*s% opacity(k))
2017 :
2018 : d_dampR_dr_in = d_dampR_dr_in &
2019 : + TEM1*s% T(k)**3*(&
2020 : 2.d0*s% Vol(k)*dVol_dr_in(I) &
2021 : - s% Vol(k)**2* &
2022 : (dCp_dr_in(I)/s% Cp(k) + dK_dr_in(I)/s% opacity(k))) &
2023 0 : /(s% Cp(k)*s% opacity(k))
2024 :
2025 : d_dampR_dT_00 = d_dampR_dT_00 &
2026 : + TEM1*s% Vol(k)**2*(3.d0*s% T(k)**2 &
2027 : - s% T(k)**3*(1.d0/s% Cp(k)*dCp_dT(I) &
2028 : + 1.d0/s% opacity(k)*dK_dT(I))) &
2029 0 : /(s% Cp(k)*s% opacity(k))
2030 :
2031 : end if
2032 :
2033 0 : s% COUPL(k) = s% SOURCE(k) - s% DAMP(k) - s% DAMPR(k)
2034 0 : dC_dr_00(I) = dsrc_dr_00 - d_damp_dr_00 - d_dampR_dr_00 !
2035 0 : dC_dr_out(I) = dsrc_dr_out - d_damp_dr_out - d_dampR_dr_out !
2036 0 : dC_dr_in(I) = dsrc_dr_in - d_damp_dr_in - d_dampR_dr_in !
2037 0 : dC_dr_in2(I) = dsrc_dr_in2 - d_damp_dr_in2 - d_dampR_dr_in2 !
2038 0 : dC_dVol_00(I) = dsrc_dVol_00 - d_damp_dVol_00 - d_dampR_dVol_00 !
2039 0 : dC_dVol_out(I) = dsrc_dVol_out - d_damp_dVol_out - d_dampR_dVol_out !
2040 0 : dC_dVol_in(I) = dsrc_dVol_in - d_damp_dVol_in - d_dampR_dVol_in !
2041 0 : dC_dT_00(I) = dsrc_dT_00 - d_damp_dT_00 - d_dampR_dT_00 !
2042 0 : dC_dT_out(I) = dsrc_dT_out - d_damp_dT_out - d_dampR_dT_out !
2043 0 : dC_dT_in(I) = dsrc_dT_in - d_damp_dT_in - d_dampR_dT_in !
2044 0 : dC_der_00(I) = dsrc_der_00 - d_damp_der_00 - d_dampR_der_00 !
2045 0 : dC_der_out(I) = dsrc_der_out - d_damp_der_out - d_dampR_der_out !
2046 0 : dC_der_in(I) = dsrc_der_in - d_damp_der_in - d_dampR_der_in !
2047 0 : dC_dw_00(I) = dsrc_dw_00 - d_damp_dw_00 - d_dampR_dw_00 !
2048 :
2049 : !test_partials = (k == s% solver_test_partials_k)
2050 0 : test_partials = .false.
2051 : if (test_partials) then
2052 : s% solver_test_partials_val = s% COUPL(k)
2053 : s% solver_test_partials_var = i_var_Vol
2054 : s% solver_test_partials_dval_dx = dC_dVol_00(I)
2055 : write(*,*) 'calc_source_sink', s% solver_test_partials_var
2056 : end if
2057 :
2058 : end if
2059 :
2060 0 : end subroutine calc_source_sink
2061 :
2062 :
2063 0 : subroutine zero_boundaries(s)
2064 : type (star_info), pointer :: s
2065 : integer :: I, k
2066 0 : do I = 1,IBOTOM
2067 0 : k = NZN+1-i
2068 0 : s% Eq(k) = 0.d0
2069 0 : dEq_dr_in2(I) = 0.d0
2070 0 : dEq_dr_in(I) = 0.d0
2071 0 : dEq_dr_00(I) = 0.d0
2072 0 : dEq_dr_out(I) = 0.d0
2073 0 : dEq_dT_in(I) = 0.d0
2074 0 : dEq_dT_00(I) = 0.d0
2075 0 : dEq_dT_out(I) = 0.d0
2076 0 : dEq_dw_00(I) = 0.d0 ! -
2077 0 : s% Chi(k) = 0.d0
2078 0 : dChi_dr_in2(I) = 0.d0
2079 0 : dChi_dr_in(I) = 0.d0
2080 0 : dChi_dr_00(I) = 0.d0
2081 0 : dChi_dr_out(I) = 0.d0
2082 0 : dChi_dT_in(I) = 0.d0
2083 0 : dChi_dT_00(I) = 0.d0
2084 0 : dChi_dT_out(I) = 0.d0
2085 0 : dChi_dw_00(I) = 0.d0 ! -
2086 0 : s% COUPL(k) = 0.d0
2087 0 : dC_dr_00(I) = 0.d0
2088 0 : dC_dr_out(I) = 0.d0
2089 0 : dC_dr_in(I) = 0.d0
2090 0 : dC_dr_in2(I) = 0.d0
2091 0 : dC_dT_in(I) = 0.d0
2092 0 : dC_dT_00(I) = 0.d0
2093 0 : dC_dT_out(I) = 0.d0
2094 0 : dC_dw_00(I) = 0.d0 ! -
2095 0 : s% Ptrb(k) = 0.d0
2096 0 : dPtrb_dr_00(I) = 0.d0
2097 0 : dPtrb_dr_in(I) = 0.d0
2098 0 : dPtrb_dw_00(I) = 0.d0 ! -
2099 : end do
2100 0 : do I = NZN,NZN
2101 0 : k = NZN+1-i
2102 0 : s% Eq(k) = 0.d0
2103 0 : dEq_dr_in2(I) = 0.d0
2104 0 : dEq_dr_in(I) = 0.d0
2105 0 : dEq_dr_00(I) = 0.d0
2106 0 : dEq_dr_out(I) = 0.d0
2107 0 : dEq_dT_in(I) = 0.d0
2108 0 : dEq_dT_00(I) = 0.d0
2109 0 : dEq_dT_out(I) = 0.d0
2110 0 : dEq_dw_00(I) = 0.d0 ! -
2111 0 : s% Chi(k) = 0.d0
2112 0 : dChi_dr_in2(I) = 0.d0
2113 0 : dChi_dr_in(I) = 0.d0
2114 0 : dChi_dr_00(I) = 0.d0
2115 0 : dChi_dr_out(I) = 0.d0
2116 0 : dChi_dT_in(I) = 0.d0
2117 0 : dChi_dT_00(I) = 0.d0
2118 0 : dChi_dT_out(I) = 0.d0
2119 0 : dChi_dw_00(I) = 0.d0 ! -
2120 0 : s% COUPL(k) = 0.d0
2121 0 : dC_dr_00(I) = 0.d0
2122 0 : dC_dr_out(I) = 0.d0
2123 0 : dC_dr_in(I) = 0.d0
2124 0 : dC_dr_in2(I) = 0.d0
2125 0 : dC_dT_in(I) = 0.d0
2126 0 : dC_dT_00(I) = 0.d0
2127 0 : dC_dT_out(I) = 0.d0
2128 0 : dC_dw_00(I) = 0.d0 ! -
2129 0 : s% Ptrb(k) = 0.d0
2130 0 : dPtrb_dr_00(I) = 0.d0
2131 0 : dPtrb_dr_in(I) = 0.d0
2132 0 : dPtrb_dw_00(I) = 0.d0 ! -
2133 : end do
2134 0 : end subroutine zero_boundaries
2135 :
2136 :
2137 0 : subroutine calc_Lt(s,i,Lt_00, &
2138 : dLt_dr_00, dLt_dr_in, dLt_dr_out, &
2139 : dLt_dVol_00, dLt_dVol_out, &
2140 : dLt_dT_00, dLt_dT_out, &
2141 : dLt_der_00, dLt_der_out, &
2142 : dLt_dw_00, dLt_dw_out)
2143 : type (star_info), pointer :: s
2144 : integer, intent(in) :: I
2145 : real(dp), intent(out) :: &
2146 : Lt_00, &
2147 : dLt_dr_00, dLt_dr_in, dLt_dr_out, &
2148 : dLt_dVol_00, dLt_dVol_out, &
2149 : dLt_dT_00, dLt_dT_out, &
2150 : dLt_der_00, dLt_der_out, &
2151 : dLt_dw_00, dLt_dw_out
2152 0 : real(dp) :: POM, POM2, POM3, TEM1, TEM2, &
2153 0 : d_POM2_dVol_00, d_POM2_dVol_out, rho2_face
2154 : integer :: k
2155 : logical :: test_partials
2156 : include 'formats'
2157 0 : k = NZN+1-i
2158 : if (I <= IBOTOM .or. I == NZN .or. ALFA == 0d0 .or. &
2159 0 : ALFAT == 0.d0 .or. k < min_k_for_turbulent_flux) then
2160 0 : Lt_00 = 0.d0
2161 0 : dLt_dr_00 = 0.d0
2162 0 : dLt_dr_in = 0.d0
2163 0 : dLt_dr_out = 0.d0
2164 0 : dLt_dVol_00 = 0.d0
2165 0 : dLt_dVol_out = 0.d0
2166 0 : dLt_dT_00 = 0.d0
2167 0 : dLt_dT_out = 0.d0
2168 0 : dLt_der_00 = 0.d0
2169 0 : dLt_der_out = 0.d0
2170 0 : dLt_dw_00 = 0.d0
2171 0 : dLt_dw_out = 0.d0
2172 : else
2173 0 : POM3 = (s% RSP_w(k-1)**3 - s% RSP_w(k)**3)/s% dm_bar(k)
2174 0 : POM = - 2.d0/3.d0*ALFA*ALFAT*(P4*(s% r(k)**2))**2
2175 0 : rho2_face = 0.5d0*(1.d0/s% Vol(k)**2 + 1.d0/s% Vol(k-1)**2)
2176 0 : POM2 = s% Hp_face(k)*rho2_face
2177 0 : Lt_00 = POM*POM2*POM3
2178 :
2179 0 : TEM1 = Lt_00/s% Hp_face(k)
2180 0 : TEM2 = Lt_00/POM2*s% Hp_face(k)
2181 :
2182 0 : d_POM2_dVol_00 = dHp_dVol_00(i)*rho2_face - s% Hp_face(k)/s% Vol(k)**3
2183 0 : d_POM2_dVol_out = dHp_dVol_out(i)*rho2_face - s% Hp_face(k)/s% Vol(k-1)**3
2184 0 : dLt_dVol_00 = POM*POM3*d_POM2_dVol_00
2185 0 : dLt_dVol_out = POM*POM3*d_POM2_dVol_out
2186 :
2187 : dLt_dr_00 = 4.d0*Lt_00/s% r(k) & !
2188 : - TEM2/s% Vol(k)**3*dVol_dr_00(I) &
2189 : - TEM2/s% Vol(k-1)**3*dVol_dr_in(i+1) &
2190 0 : + TEM1*dHp_dr_00(I)
2191 : dLt_dr_in = &
2192 : - TEM2/s% Vol(k)**3*dVol_dr_in(I) & !
2193 0 : + TEM1*dHp_dr_in(I)
2194 : dLt_dr_out = &
2195 : - TEM2/s% Vol(k-1)**3*dVol_dr_00(i+1) & !
2196 0 : + TEM1*dHp_dr_out(I)
2197 0 : dLt_dT_00 = TEM1*dHp_dT_00(I) !
2198 0 : dLt_dT_out = TEM1*dHp_dT_out(I) !
2199 0 : dLt_der_00 = TEM1*dHp_der_00(I) !
2200 0 : dLt_der_out = TEM1*dHp_der_out(I) !
2201 0 : TEM1 = POM*POM2*3.0d0/s% dm_bar(k)
2202 0 : dLt_dw_00 = -TEM1*s% RSP_w(k)**2 !
2203 0 : dLt_dw_out = TEM1*s% RSP_w(k-1)**2 !
2204 : end if
2205 0 : if (i > 0) s% Lt(k) = Lt_00
2206 :
2207 : !test_partials = (k-1 == s% solver_test_partials_k)
2208 0 : test_partials = .false.
2209 : if (test_partials) then
2210 : s% solver_test_partials_val = Lt_00
2211 : s% solver_test_partials_var = i_var_Vol
2212 : s% solver_test_partials_dval_dx = dLt_dVol_out
2213 : write(*,*) 'calc_Lt', s% solver_test_partials_var
2214 : end if
2215 0 : end subroutine calc_Lt
2216 :
2217 :
2218 0 : subroutine calc_Lc(s,i,Lc_00, &
2219 : dLc_dr_in, dLc_dr_00, dLc_dr_out, &
2220 : dLc_dVol_00, dLc_dVol_out, &
2221 : dLc_dT_00, dLc_dT_out, &
2222 : dLc_der_00, dLc_der_out, &
2223 : dLc_dw_00, dLc_dw_out)
2224 : type (star_info), pointer :: s
2225 : integer, intent(in) :: I
2226 : real(dp), intent(out) :: &
2227 : Lc_00, dLc_dr_in, dLc_dr_00, dLc_dr_out, &
2228 : dLc_dVol_00, dLc_dVol_out, &
2229 : dLc_dT_00, dLc_dT_out, &
2230 : dLc_der_00, dLc_der_out, &
2231 : dLc_dw_00, dLc_dw_out
2232 0 : real(dp) :: POM,POM2,POM3
2233 : integer :: k
2234 : logical :: test_partials
2235 : include 'formats'
2236 0 : k = NZN+1-i
2237 0 : if (I <= IBOTOM .or. I == NZN .or. ALFA == 0d0)then
2238 0 : Lc_00 = 0.d0 !DODANE SMOLEC
2239 0 : dLc_dr_in = 0.d0
2240 0 : dLc_dr_00 = 0.d0
2241 0 : dLc_dr_out = 0.d0
2242 0 : dLc_dVol_00 = 0.d0
2243 0 : dLc_dVol_out = 0.d0
2244 0 : dLc_dT_00 = 0.d0
2245 0 : dLc_dT_out = 0.d0
2246 0 : dLc_der_00 = 0.d0
2247 0 : dLc_der_out = 0.d0
2248 0 : dLc_dw_00 = 0.d0
2249 0 : dLc_dw_out = 0.d0 ! -
2250 0 : else if (s% RSP_w(k) < EFL0*1d-8)then
2251 0 : Lc_00 = 0.d0
2252 0 : dLc_dr_00 = 0.d0
2253 0 : dLc_dr_in = 0.d0
2254 0 : dLc_dr_out = 0.d0
2255 0 : dLc_dVol_00 = 0.d0
2256 0 : dLc_dVol_out = 0.d0
2257 0 : dLc_dT_00 = 0.d0
2258 0 : dLc_dT_out = 0.d0
2259 0 : dLc_der_00 = 0.d0
2260 0 : dLc_der_out = 0.d0
2261 0 : dLc_dw_00 = 0.d0
2262 0 : dLc_dw_out = 0.d0
2263 : else
2264 :
2265 0 : POM3 = 0.5d0*(s% RSP_w(k) + s% RSP_w(k-1))
2266 :
2267 : POM = P4*(s% r(k)**2)*(ALFAC/ALFAS)* &
2268 0 : 0.5d0*(s% T(k)/s% Vol(k) + s% T(k-1)/s% Vol(k-1))
2269 0 : Lc_00 = POM*s% PII(k)*POM3
2270 :
2271 0 : dLc_dw_00 = POM*s% PII(k)*0.5d0 !
2272 0 : if (I >= NZN - 1) then
2273 0 : dLc_dw_out = 0.d0
2274 : else
2275 0 : dLc_dw_out = POM*s% PII(k)*0.5d0 !
2276 : end if
2277 :
2278 0 : POM2 = P4*(s% r(k)**2)*s% PII(k)*POM3*0.5d0*(ALFAC/ALFAS)
2279 0 : POM = POM*POM3
2280 :
2281 : dLc_dr_00 = & !
2282 : POM*dPII_dr_00(I) &
2283 : + 2.d0*Lc_00/s% r(k) &
2284 : - POM2*(s% T(k)/(s% Vol(k)**2)*dVol_dr_00(I) + &
2285 0 : s% T(k-1)/(s% Vol(k-1)**2)*dVol_dr_in(i+1))
2286 : dLc_dr_in = &
2287 : POM*dPII_dr_in(I) &
2288 0 : - POM2*s% T(k)/(s% Vol(k)**2)*dVol_dr_in(I) !
2289 : dLc_dr_out = &
2290 : POM*dPII_dr_out(I) &
2291 0 : - POM2*s% T(k-1)/(s% Vol(k-1)**2)*dVol_dr_00(i+1) !
2292 :
2293 : dLc_dVol_00 = &
2294 : POM*dPII_dVol_00(I) &
2295 0 : - POM2*s% T(k)/(s% Vol(k)**2) !
2296 : dLc_dVol_out = &
2297 : POM*dPII_dVol_out(I) &
2298 0 : - POM2*s% T(k-1)/(s% Vol(k-1)**2) !
2299 :
2300 0 : dLc_dT_00 = POM*dPII_dT_00(I) + POM2/s% Vol(k) !
2301 0 : dLc_dT_out = POM*dPII_dT_out(I) + POM2/s% Vol(k-1) !
2302 :
2303 0 : dLc_der_00 = POM*dPII_der_00(I) !
2304 0 : dLc_der_out = POM*dPII_der_out(I) !
2305 :
2306 : end if
2307 :
2308 0 : if (i > 0) s% Lc(k) = Lc_00
2309 :
2310 : !test_partials = (k-1 == s% solver_test_partials_k)
2311 0 : test_partials = .false.
2312 : if (test_partials) then
2313 : s% solver_test_partials_val = Lc_00
2314 : s% solver_test_partials_var = i_var_Vol
2315 : s% solver_test_partials_dval_dx = dLc_dVol_out
2316 : write(*,*) 'calc_Lc', s% solver_test_partials_var
2317 : end if
2318 0 : end subroutine calc_Lc
2319 :
2320 : ! in diffusion limit, radiative flux equation reduces to Fr calculated from d_erad_dm as below.
2321 : ! note that can have nonequilibrium diffusion regime with different T for gas and photons.
2322 : ! this happens when absorption mean opacity is different than Planck mean opacity.
2323 0 : subroutine calc_Fr(s, i, Fr_00, & !rs Stellingwerf 1975, Appendix A
2324 : dFr_dr_out, dFr_dr_00, dFr_dr_in, &
2325 : dFr_dVol_out, dFr_dVol_00, &
2326 : dFr_dT_out, dFr_dT_00, &
2327 : dFr_der_out, dFr_der_00)
2328 : type (star_info), pointer :: s
2329 : integer, intent(in) :: I
2330 : real(dp), intent(out) :: &
2331 : Fr_00, dFr_dr_out, dFr_dr_00, dFr_dr_in, &
2332 : dFr_dVol_out, dFr_dVol_00, &
2333 : dFr_dT_out, dFr_dT_00, &
2334 : dFr_der_out, dFr_der_00
2335 : real(dp) :: &
2336 0 : W_00, d_W_00_dr_in, d_W_00_dr_00, d_W_00_dVol_00, d_W_00_der_00, &
2337 0 : W_out, d_W_out_dr_00, d_W_out_dr_out, d_W_out_dVol_out, d_W_out_der_out, &
2338 0 : Prad_factor, Fr2a, d_Fr2a_dW_00, d_Fr2a_dW_out, &
2339 0 : Fr2b, d_Fr2b_dW_00, d_Fr2b_dW_out, &
2340 0 : BW, kap_00, kap_out, BK, Fr1, Fr2, Fr3, &
2341 0 : d_Fr_dK_00, d_Fr_dK_out, d_Fr_dW_00, d_Fr_dW_out
2342 : integer :: k
2343 : logical :: test_partials
2344 : include 'formats'
2345 :
2346 0 : k = NZN+1-i
2347 :
2348 0 : if (i < 1) then
2349 0 : if (s% RSP_hydro_only) then
2350 0 : Fr_00 = 0d0
2351 : else
2352 0 : Fr_00 = s% L_center
2353 : end if
2354 0 : Fr_00 = Fr_00/(4d0*pi*s% r_center**2)
2355 0 : dFr_dr_out = 0
2356 0 : dFr_dr_00 = 0
2357 0 : dFr_dr_in = 0
2358 0 : dFr_dVol_out = 0
2359 0 : dFr_dVol_00 = 0
2360 0 : dFr_dT_out = 0
2361 0 : dFr_dT_00 = 0
2362 0 : dFr_der_out = 0
2363 0 : dFr_der_00 = 0
2364 0 : return
2365 : end if
2366 :
2367 0 : Prad_factor = 3d0/crad ! 3d0 to cancel the 1/3d0 factor in CL below
2368 0 : W_00 = Prad_factor*s% Prad(k) ! replaces s% T(k)**4
2369 0 : d_W_00_dVol_00 = Prad_factor*d_Pr_dVol(i)
2370 0 : d_W_00_dr_in = Prad_factor*d_Pr_dr_in(i)
2371 0 : d_W_00_dr_00 = Prad_factor*d_Pr_dr_00(i)
2372 0 : d_W_00_der_00 = Prad_factor*d_Pr_der(i)
2373 :
2374 0 : if (k == 1) then ! surface
2375 0 : Fr1 = s% g_Edd*4d0*SIG
2376 0 : Fr_00 = Fr1*W_00 ! s% T(k)**4 => W_00
2377 0 : dFr_dr_out = 0
2378 0 : dFr_dr_in = Fr1*d_W_00_dr_in
2379 0 : dFr_dr_00 = Fr1*d_W_00_dr_00
2380 0 : dFr_dVol_out = 0
2381 0 : dFr_dVol_00 = Fr1*d_W_00_dVol_00
2382 0 : dFr_dT_out = 0
2383 0 : dFr_dT_00 = 0
2384 0 : dFr_der_out = 0
2385 0 : dFr_der_00 = Fr1*d_W_00_der_00
2386 0 : return
2387 : end if
2388 :
2389 0 : W_out = Prad_factor*s% Prad(k-1) ! replaces s% T(k-1)**4
2390 0 : d_W_out_dVol_out = Prad_factor*d_Pr_dVol(i+1)
2391 0 : d_W_out_dr_00 = Prad_factor*d_Pr_dr_in(i+1)
2392 0 : d_W_out_dr_out = Prad_factor*d_Pr_dr_00(i+1)
2393 0 : d_W_out_der_out = Prad_factor*d_Pr_der(i+1)
2394 :
2395 0 : BW = log(W_out/W_00)
2396 0 : if (abs(BW) < 1d-30) then
2397 0 : Fr_00 = 0
2398 0 : dFr_dr_out = 0
2399 0 : dFr_dr_in = 0
2400 0 : dFr_dr_00 = 0
2401 0 : dFr_dVol_out = 0
2402 0 : dFr_dVol_00 = 0
2403 0 : dFr_dT_out = 0
2404 0 : dFr_dT_00 = 0
2405 0 : dFr_der_out = 0
2406 0 : dFr_der_00 = 0
2407 0 : return
2408 : end if
2409 :
2410 0 : kap_00 = s% opacity(k)
2411 0 : kap_out = s% opacity(k-1)
2412 0 : BK = log(kap_out/kap_00)
2413 :
2414 0 : Fr1 = -CL*s% r(k)**2/(4d0*pi*s% dm_bar(k)) ! CL = 4d0*(4d0*PI)**2*SIG/3d0
2415 :
2416 0 : Fr2a = W_out/kap_out - W_00/kap_00
2417 0 : d_Fr2a_dW_00 = -1d0/kap_00
2418 0 : d_Fr2a_dW_out = 1d0/kap_out
2419 :
2420 0 : Fr2b = 1d0 - BK/BW
2421 0 : d_Fr2b_dW_00 = -BK/BW**2/W_00
2422 0 : d_Fr2b_dW_out = BK/BW**2/W_out
2423 :
2424 0 : Fr2 = Fr2a/Fr2b
2425 0 : Fr_00 = Fr1*Fr2
2426 0 : d_Fr_dW_00 = Fr1*(d_Fr2a_dW_00 - Fr2*d_Fr2b_dW_00)/Fr2b
2427 0 : d_Fr_dW_out = Fr1*(d_Fr2a_dW_out - Fr2*d_Fr2b_dW_out)/Fr2b
2428 :
2429 0 : Fr3 = Fr1/(BW - BK)
2430 0 : d_Fr_dK_00 = (Fr3/kap_00)*(W_00*BW/kap_00 - Fr2)
2431 0 : d_Fr_dK_out = -(Fr3/kap_out)*(W_out*BW/kap_out - Fr2)
2432 :
2433 : dFr_dr_in = & !
2434 : + d_Fr_dK_00*dK_dr_in(i) &
2435 0 : + d_Fr_dW_00*d_W_00_dr_in
2436 : dFr_dr_00 = 2d0*Fr_00/s% r(k) & !
2437 : + d_Fr_dK_00*dK_dr_00(i) &
2438 : + d_Fr_dK_out*dK_dr_in(i+1) &
2439 : + d_Fr_dW_00*d_W_00_dr_00 &
2440 0 : + d_Fr_dW_out*d_W_out_dr_00
2441 : dFr_dr_out = & !
2442 : + d_Fr_dK_out*dK_dr_00(i+1) &
2443 0 : + d_Fr_dW_out*d_W_out_dr_out
2444 :
2445 : dFr_dVol_00 = &
2446 : + d_Fr_dK_00*dK_dVol(i) &
2447 0 : + d_Fr_dW_00*d_W_00_dVol_00
2448 : dFr_dVol_out = &
2449 : + d_Fr_dK_out*dK_dVol(i+1) &
2450 0 : + d_Fr_dW_out*d_W_out_dVol_out
2451 :
2452 : dFr_dT_out = & !
2453 0 : + d_Fr_dK_out*dK_dT(i+1)
2454 : dFr_dT_00 = & !
2455 0 : + d_Fr_dK_00*dK_dT(i)
2456 :
2457 : dFr_der_out = & !
2458 0 : + d_Fr_dW_out*d_W_out_der_out
2459 : dFr_der_00 = & !
2460 0 : + d_Fr_dW_00*d_W_00_der_00
2461 :
2462 : !test_partials = (k-1 == s% solver_test_partials_k)
2463 0 : test_partials = .false.
2464 : if (test_partials) then
2465 : s% solver_test_partials_val = Fr_00
2466 : s% solver_test_partials_var = i_var_Vol
2467 : s% solver_test_partials_dval_dx = dFr_dVol_out
2468 : write(*,*) 'calc_Fr', s% solver_test_partials_var
2469 : end if
2470 :
2471 : if (call_is_bad) then
2472 : if (is_bad(Fr_00)) then
2473 : !$omp critical (rsp_step_5)
2474 : write(*,2) 'Fr_00', k, Fr_00
2475 : write(*,2) 'Fr1', k, Fr1
2476 : write(*,2) 'Fr2', k, Fr2
2477 : write(*,2) 'Fr2a', k, Fr2a
2478 : write(*,2) 'Fr2b', k, Fr2b
2479 : write(*,2) 'W_00', k, W_00
2480 : write(*,2) 'W_out', k, W_out
2481 : write(*,2) 'kap_00', k, kap_00
2482 : write(*,2) 'kap_out', k, kap_out
2483 : write(*,2) 'r(k)', k, s% r(k)
2484 : write(*,2) 'dm_bar(k)', k, s% dm_bar(k)
2485 : write(*,2) 'erad(k)', k, s% erad(k)
2486 : write(*,2) 'erad(k-1)', k-1, s% erad(k-1)
2487 : write(*,2) 'BK', k, BK
2488 : write(*,2) 'BW', k, BW
2489 : write(*,2) 'nz', s% nz
2490 : call mesa_error(__FILE__,__LINE__,'calc_Fr')
2491 : !$omp end critical (rsp_step_5)
2492 : end if
2493 : end if
2494 :
2495 : end subroutine calc_Fr
2496 :
2497 :
2498 0 : subroutine rsp_calc_XP(s, P_surf, i, with_Prad, & ! time weighted combined pressure
2499 : XP, d_XP_dVol_00, d_XP_dT_00, d_XP_der_00, &
2500 : d_XP_dw_00, d_XP_dr_in, d_XP_dr_00)
2501 : type (star_info), pointer :: s
2502 : real(dp), intent(in) :: P_surf
2503 : integer, intent(in) :: i
2504 : logical, intent(in) :: with_Prad
2505 : real(dp), intent(out) :: &
2506 : XP, d_XP_dVol_00, d_XP_dT_00, d_XP_der_00, &
2507 : d_XP_dw_00, d_XP_dr_in, d_XP_dr_00
2508 0 : real(dp) :: T_surf, Prad_factor
2509 : logical :: test_partials
2510 : integer :: k
2511 : include 'formats'
2512 :
2513 0 : k = NZN+1-i
2514 0 : if (k == 0) then ! pressure outside of surface
2515 0 : if (s% RSP_use_atm_grey_with_kap_for_Psurf) then
2516 0 : XP = P_surf
2517 0 : else if (s% RSP_use_Prad_for_Psurf) then
2518 0 : T_surf = s% T_start(1)
2519 0 : XP = crad*T_surf**4/3d0
2520 : else
2521 0 : XP = 0.0d0
2522 : end if
2523 0 : d_XP_dVol_00 = 0d0
2524 0 : d_XP_dT_00 = 0d0
2525 0 : d_XP_der_00 = 0d0
2526 0 : d_XP_dw_00 = 0d0
2527 0 : d_XP_dr_in = 0d0
2528 0 : d_XP_dr_00 = 0d0
2529 : return
2530 : end if
2531 0 : if (with_Prad) then
2532 : Prad_factor = 1d0
2533 : else
2534 0 : Prad_factor = 0d0
2535 : end if
2536 :
2537 : XP = THETA*(s% Pgas(k) + Prad_factor*s% Prad(k)) &
2538 : + THETA1*(s% Pgas_start(k) + Prad_factor*s% Prad_start(k)) &
2539 : + THETAQ*s% Pvsc(k) + THETAQ1*s% Pvsc_start(k) &
2540 0 : + THETAT*s% Ptrb(k) + THETAT1*s% Ptrb_start(k)
2541 : d_XP_dVol_00 = &
2542 : THETA*(d_Pg_dVol(i) + Prad_factor*d_Pr_dVol(i)) &
2543 : + THETAQ*d_Pvsc_dVol(i) &
2544 0 : + THETAT*dPtrb_dVol_00(i)
2545 0 : d_XP_dT_00 = THETA*d_Pg_dT(i) + THETAQ*d_Pvsc_dT(i) !
2546 0 : d_XP_der_00 = THETA*Prad_factor*d_Pr_der(i) !
2547 0 : d_XP_dw_00 = THETAT*dPtrb_dw_00(i) !
2548 : d_XP_dr_in = & !
2549 : THETA*(d_Pg_dr_in(i) + Prad_factor*d_Pr_dr_in(i)) &
2550 : + THETAQ*d_Pvsc_dr_in(i) &
2551 0 : + THETAT*dPtrb_dr_in(i)
2552 : d_XP_dr_00 = & !
2553 : THETA*(d_Pg_dr_00(i) + Prad_factor*d_Pr_dr_00(i)) &
2554 : + THETAQ*d_Pvsc_dr_00(i) &
2555 0 : + THETAT*dPtrb_dr_00(i)
2556 :
2557 : if (call_is_bad) then
2558 : if (is_bad(XP)) then
2559 : !$omp critical (rsp_step_6)
2560 : write(*,2) 'XP', k, XP
2561 : write(*,2) 's% Pgas(k)', k, s% Pgas(k)
2562 : write(*,2) 's% Prad(k)', k, s% Prad(k)
2563 : write(*,2) 's% Pvsc(k)', k, s% Pvsc(k)
2564 : write(*,2) 's% Ptrb(k)', k, s% Ptrb(k)
2565 : write(*,2) 'THETA', k, THETA
2566 : write(*,2) 'THETAQ', k, THETAQ
2567 : write(*,2) 'THETAT', k, THETAT
2568 : call mesa_error(__FILE__,__LINE__,'rsp_calc_XP')
2569 : !$omp end critical (rsp_step_6)
2570 : end if
2571 : end if
2572 :
2573 : !test_partials = (k == s% solver_test_partials_k)
2574 0 : test_partials = .false.
2575 :
2576 : if (test_partials) then
2577 : s% solver_test_partials_val = XP
2578 : s% solver_test_partials_var = i_var_Vol
2579 : s% solver_test_partials_dval_dx = d_XP_dVol_00
2580 : write(*,*) 'rsp_calc_XP', s% solver_test_partials_var
2581 : end if
2582 : end subroutine rsp_calc_XP
2583 :
2584 :
2585 0 : subroutine calc_equations(s,i,P_surf)
2586 : type (star_info), pointer :: s
2587 : integer, intent(in) :: i
2588 : real(dp), intent(in) :: P_surf
2589 0 : call calc_face_equations(s,i,P_surf)
2590 0 : call calc_cell_equations(s,i,P_surf,.true.,.true.,.true.)
2591 0 : end subroutine calc_equations
2592 :
2593 :
2594 0 : subroutine calc_face_equations(s,i,P_surf)
2595 : type (star_info), pointer :: s
2596 : integer, intent(in) :: i
2597 : real(dp), intent(in) :: P_surf
2598 0 : call acceleration_eqn(s,i,P_surf)
2599 0 : call Fr_eqn(s,i)
2600 0 : end subroutine calc_face_equations
2601 :
2602 :
2603 0 : subroutine calc_cell_equations( &
2604 : s, i, P_surf, do_etot, do_eturb, do_erad)
2605 : type (star_info), pointer :: s
2606 : integer, intent(in) :: i
2607 : real(dp), intent(in) :: P_surf
2608 : logical, intent(in) :: do_etot, do_eturb, do_erad
2609 : integer :: k
2610 : real(dp) :: &
2611 0 : Lt_00, Lt_00_start, Lt_in, Lt_in_start, &
2612 0 : dLt_00_dr_00, dLt_00_dr_in, dLt_00_dr_out, &
2613 0 : dLt_00_dVol_00, dLt_00_dVol_out, &
2614 0 : dLt_00_dT_00, dLt_00_dT_out, &
2615 0 : dLt_00_der_00, dLt_00_der_out, &
2616 0 : dLt_00_dw_00, dLt_00_dw_out, &
2617 0 : dLt_in_dr_in, dLt_in_dr_in2, dLt_in_dr_00, &
2618 0 : dLt_in_dVol_in, dLt_in_dVol_00, &
2619 0 : dLt_in_dT_in, dLt_in_dT_00, &
2620 0 : dLt_in_der_in, dLt_in_der_00, &
2621 0 : dLt_in_dw_in, dLt_in_dw_00
2622 :
2623 : include 'formats'
2624 :
2625 : ! HR = -residual
2626 : ! partials of residual go in HD
2627 :
2628 0 : k = NZN+1-i
2629 :
2630 0 : call get_Lt
2631 :
2632 0 : if (do_etot) call total_energy_eqn(s, i, P_surf, &
2633 : Lt_00, Lt_00_start, Lt_in, Lt_in_start, &
2634 : dLt_00_dr_00, dLt_00_dr_in, dLt_00_dr_out, &
2635 : dLt_00_dVol_00, dLt_00_dVol_out, &
2636 : dLt_00_dT_00, dLt_00_dT_out, &
2637 : dLt_00_der_00, dLt_00_der_out, &
2638 : dLt_00_dw_00, dLt_00_dw_out, &
2639 : dLt_in_dr_in, dLt_in_dr_in2, dLt_in_dr_00, &
2640 : dLt_in_dVol_in, dLt_in_dVol_00, &
2641 : dLt_in_dT_in, dLt_in_dT_00, &
2642 : dLt_in_der_in, dLt_in_der_00, &
2643 0 : dLt_in_dw_in, dLt_in_dw_00)
2644 :
2645 0 : if (do_eturb) call turbulent_energy_eqn(s, i, &
2646 : Lt_00, Lt_00_start, Lt_in, Lt_in_start, &
2647 : dLt_00_dr_00, dLt_00_dr_in, dLt_00_dr_out, &
2648 : dLt_00_dVol_00, dLt_00_dVol_out, &
2649 : dLt_00_dT_00, dLt_00_dT_out, &
2650 : dLt_00_der_00, dLt_00_der_out, &
2651 : dLt_00_dw_00, dLt_00_dw_out, &
2652 : dLt_in_dr_in, dLt_in_dr_in2, dLt_in_dr_00, &
2653 : dLt_in_dVol_in, dLt_in_dVol_00, &
2654 : dLt_in_dT_in, dLt_in_dT_00, &
2655 : dLt_in_der_in, dLt_in_der_00, &
2656 0 : dLt_in_dw_in, dLt_in_dw_00)
2657 :
2658 0 : if (do_erad) call erad_eqn(s,i)
2659 :
2660 0 : if (I == 1) call inner_boundary_eqn
2661 :
2662 : contains
2663 :
2664 0 : subroutine inner_boundary_eqn
2665 0 : HR(1) = 0.d0
2666 0 : HD(1:LD_HD,1) = 0.d0
2667 0 : HD(i_r_dr_00,1) = 1.d0
2668 0 : end subroutine inner_boundary_eqn
2669 :
2670 0 : subroutine get_Lt
2671 : integer :: k
2672 0 : k = NZN+1-i
2673 : call calc_Lt(s,i - 1,Lt_in, &
2674 : dLt_in_dr_in, dLt_in_dr_in2, dLt_in_dr_00, &
2675 : dLt_in_dVol_in, dLt_in_dVol_00, &
2676 : dLt_in_dT_in, dLt_in_dT_00, &
2677 : dLt_in_der_in, dLt_in_der_00, &
2678 0 : dLt_in_dw_in, dLt_in_dw_00)
2679 : call calc_Lt(s,i,Lt_00, &
2680 : dLt_00_dr_00, dLt_00_dr_in, dLt_00_dr_out, &
2681 : dLt_00_dVol_00, dLt_00_dVol_out, &
2682 : dLt_00_dT_00, dLt_00_dT_out, &
2683 : dLt_00_der_00, dLt_00_der_out, &
2684 : dLt_00_dw_00, dLt_00_dw_out)
2685 0 : if (I == NZN) then
2686 0 : Lt_00_start = 0.d0
2687 : else
2688 0 : Lt_00_start = s% Lt_start(k)
2689 : end if
2690 0 : if (i == 1) then
2691 0 : Lt_in_start = 0
2692 : else
2693 0 : Lt_in_start = s% Lt_start(k+1)
2694 : end if
2695 0 : end subroutine get_Lt
2696 :
2697 : end subroutine calc_cell_equations
2698 :
2699 :
2700 0 : subroutine acceleration_eqn(s, i, P_surf)
2701 : type (star_info), pointer :: s
2702 : integer, intent(in) :: i
2703 : real(dp), intent(in) :: P_surf
2704 : integer :: IR, k
2705 0 : real(dp) :: dt, dm_bar, residual, area, d_area_dr_00, &
2706 0 : grav, dXP_dm, Uq1, A_dm, R_00, dv_dr, dvdt_factor, &
2707 0 : XP_00, dXP_00_dVol_00, dXP_00_dT_00, dXP_00_der_00, &
2708 0 : dXP_00_dw_00, dXP_00_dr_in, dXP_00_dr_00, &
2709 0 : XP_out, dXP_out_dVol_out, dXP_out_dT_out, dXP_out_der_out, &
2710 : dXP_out_dw_out, dXP_out_dr_00, dXP_out_dr_out, &
2711 0 : Chi_00, Chi_out, d_Chi_out_dVol_00, &
2712 0 : d_Chi_out_dr_in, d_Chi_out_dr_00, d_Chi_out_dr_out, d_Chi_out_dr_out2, &
2713 0 : d_Chi_out_dT_00, d_Chi_out_dT_out, d_Chi_out_dT_out2, &
2714 0 : d_Chi_out_der_00, d_Chi_out_der_out, d_Chi_out_der_out2, &
2715 0 : d_Chi_out_dw_out, &
2716 0 : d_Chi_00_dr_in2, d_Chi_00_dr_in, d_Chi_00_dr_00, d_Chi_00_dr_out, &
2717 0 : d_Chi_00_dT_in, d_Chi_00_dT_00, d_Chi_00_dT_out, &
2718 0 : d_Chi_00_der_in, d_Chi_00_der_00, d_Chi_00_der_out, &
2719 0 : d_Chi_00_dw_00, d_Chi_00_dVol_00, &
2720 0 : d_Uq1_dr_00, d_Uq_dVol_00, d_Uq_dr_in2, d_Uq_dw_00, d_Uq_dw_out, &
2721 0 : d_Uq_dT_in, d_Uq_dT_00, d_Uq_dT_out, d_Uq_dT_out2, &
2722 0 : d_Uq_der_in, d_Uq_der_00, d_Uq_der_out, d_Uq_der_out2, &
2723 0 : d_Uq_dr_in, d_Uq_dr_00, d_Uq_dr_out, d_Uq_dr_out2, &
2724 0 : kap_face, d_kap_dVol_00, d_kap_dT_00, d_kap_dT_out, &
2725 0 : d_kap_dr_in, d_kap_dr_00, d_kap_dr_out, &
2726 0 : Fr_tw, Fr_term, d_Fr_term_dFr_00, d_Fr_term_dVol_00, d_Fr_term_dT_00, d_Fr_term_dT_out, &
2727 0 : d_Fr_term_dr_in, d_Fr_term_dr_00, d_Fr_term_dr_out
2728 : logical :: test_partials, use_Prad
2729 :
2730 : include 'formats'
2731 :
2732 0 : use_Prad = .true. ! s% RSP_accel_eqn_use_Prad_instead_of_Fr_term
2733 :
2734 0 : dvdt_factor = 1d0
2735 :
2736 0 : k = NZN+1-i
2737 0 : IR = i_var_R + NV*(i-1)
2738 0 : HD(1:LD_HD,IR) = 0.d0
2739 :
2740 0 : dt = s% dt
2741 : !if (s% use_compression_outer_BC .and. I == NZN) then
2742 : ! call mesa_error(__FILE__,__LINE__,'no rsp support for use_compression_outer_BC')
2743 : !end if
2744 :
2745 : ! XP doesn't include Prad for acceleration equation
2746 : ! instead introduce term using Fr
2747 :
2748 : call rsp_calc_XP(s, P_surf, i+1, use_Prad, &
2749 : XP_out, dXP_out_dVol_out, dXP_out_dT_out, dXP_out_der_out, &
2750 0 : dXP_out_dw_out, dXP_out_dr_00, dXP_out_dr_out)
2751 :
2752 : call rsp_calc_XP(s, P_surf, i, use_Prad, &
2753 : XP_00, dXP_00_dVol_00, dXP_00_dT_00, dXP_00_der_00, &
2754 0 : dXP_00_dw_00, dXP_00_dr_in, dXP_00_dr_00)
2755 :
2756 0 : area = P43*(s% r(k)**2 + s% r(k)*s% r_start(k) + s% r_start(k)**2)
2757 0 : d_area_dr_00 = P43*(2d0*s% r(k) + s% r_start(k))
2758 0 : dm_bar = s% dm_bar(k)
2759 0 : A_dm = area/dm_bar
2760 0 : dv_dr = 2d0/dt
2761 0 : dXP_dm = (XP_out - XP_00)/dm_bar
2762 0 : grav = - s% cgrav(k)*s% m(k)/(s% r(k)*s% r_start(k))
2763 :
2764 0 : R_00 = 0.5d0*(s% r(k) + s% r_start(k))
2765 0 : Uq1 = P4/(dm_bar*R_00)
2766 0 : d_Uq1_dr_00 = -Uq1*0.5d0/R_00
2767 :
2768 0 : Chi_00 = THETAU*s% Chi(k) + THETAU1*s% Chi_start(k)
2769 0 : d_Chi_00_dVol_00 = THETAU*dChi_dVol_00(i)
2770 0 : d_Chi_00_dT_in = THETAU*dChi_dT_in(i)
2771 0 : d_Chi_00_dT_00 = THETAU*dChi_dT_00(i)
2772 0 : d_Chi_00_dT_out = THETAU*dChi_dT_out(i)
2773 0 : d_Chi_00_der_in = THETAU*dChi_der_in(i)
2774 0 : d_Chi_00_der_00 = THETAU*dChi_der_00(i)
2775 0 : d_Chi_00_der_out = THETAU*dChi_der_out(i)
2776 0 : d_Chi_00_dw_00 = THETAU*dChi_dw_00(i)
2777 0 : d_Chi_00_dr_in2 = THETAU*dChi_dr_in2(i)
2778 0 : d_Chi_00_dr_in = THETAU*dChi_dr_in(i)
2779 0 : d_Chi_00_dr_00 = THETAU*dChi_dr_00(i)
2780 0 : d_Chi_00_dr_out = THETAU*dChi_dr_out(i)
2781 0 : if (I == NZN) then
2782 0 : Chi_out = 0d0
2783 0 : d_Chi_out_dVol_00 = 0d0
2784 0 : d_Chi_out_dT_00 = 0d0
2785 0 : d_Chi_out_dT_out = 0d0
2786 0 : d_Chi_out_dT_out2 = 0d0
2787 0 : d_Chi_out_der_00 = 0d0
2788 0 : d_Chi_out_der_out = 0d0
2789 0 : d_Chi_out_der_out2 = 0d0
2790 0 : d_Chi_out_dw_out = 0d0
2791 0 : d_Chi_out_dr_in = 0d0
2792 0 : d_Chi_out_dr_00 = 0d0
2793 0 : d_Chi_out_dr_out = 0d0
2794 0 : d_Chi_out_dr_out2 = 0d0
2795 : else
2796 0 : Chi_out = THETAU*s% Chi(k-1) + THETAU1*s% Chi_start(k-1)
2797 0 : d_Chi_out_dVol_00 = THETAU*dChi_dVol_in(i+1)
2798 0 : d_Chi_out_dT_00 = THETAU*dChi_dT_in(i+1)
2799 0 : d_Chi_out_dT_out = THETAU*dChi_dT_00(i+1)
2800 0 : d_Chi_out_dT_out2 = THETAU*dChi_dT_out(i+1)
2801 0 : d_Chi_out_der_00 = THETAU*dChi_der_in(i+1)
2802 0 : d_Chi_out_der_out = THETAU*dChi_der_00(i+1)
2803 0 : d_Chi_out_der_out2 = THETAU*dChi_der_out(i+1)
2804 0 : d_Chi_out_dw_out = THETAU*dChi_dw_00(i+1)
2805 0 : d_Chi_out_dr_in = THETAU*dChi_dr_in2(i+1)
2806 0 : d_Chi_out_dr_00 = THETAU*dChi_dr_in(i+1)
2807 0 : d_Chi_out_dr_out = THETAU*dChi_dr_00(i+1)
2808 0 : d_Chi_out_dr_out2 = THETAU*dChi_dr_out(i+1)
2809 : end if
2810 :
2811 0 : s% Uq(k) = Uq1*(Chi_out - Chi_00)
2812 0 : d_Uq_dVol_00 = Uq1*(d_Chi_out_dVol_00 - d_Chi_00_dVol_00)
2813 0 : d_Uq_dT_in = -Uq1*d_Chi_00_dT_in
2814 0 : d_Uq_dT_00 = Uq1*(d_Chi_out_dT_00 - d_Chi_00_dT_00)
2815 0 : d_Uq_dT_out = Uq1*(d_Chi_out_dT_out - d_Chi_00_dT_out)
2816 0 : d_Uq_dT_out2 = Uq1*d_Chi_out_dT_out2
2817 0 : d_Uq_der_in = -Uq1*d_Chi_00_der_in
2818 0 : d_Uq_der_00 = Uq1*(d_Chi_out_der_00 - d_Chi_00_der_00)
2819 0 : d_Uq_der_out = Uq1*(d_Chi_out_der_out - d_Chi_00_der_out)
2820 0 : d_Uq_der_out2 = Uq1*d_Chi_out_der_out2
2821 0 : d_Uq_dw_00 = -Uq1*d_Chi_00_dw_00
2822 0 : d_Uq_dw_out = Uq1*d_Chi_out_dw_out
2823 0 : d_Uq_dr_in2 = -Uq1*d_Chi_00_dr_in2
2824 0 : d_Uq_dr_in = Uq1*(d_Chi_out_dr_in - d_Chi_00_dr_in)
2825 : d_Uq_dr_00 = &
2826 : Uq1*(d_Chi_out_dr_00 - d_Chi_00_dr_00) &
2827 0 : + d_Uq1_dr_00*(Chi_out - Chi_00)
2828 0 : d_Uq_dr_out = Uq1*(d_Chi_out_dr_out - d_Chi_00_dr_out)
2829 0 : d_Uq_dr_out2 = Uq1*d_Chi_out_dr_out2
2830 :
2831 : if (use_Prad) then
2832 0 : Fr_term = 0
2833 0 : d_Fr_term_dFr_00 = 0
2834 0 : d_Fr_term_dVol_00 = 0
2835 0 : d_Fr_term_dT_00 = 0
2836 0 : d_Fr_term_dT_out = 0
2837 0 : d_Fr_term_dr_in = 0
2838 : d_Fr_term_dr_00 = 0
2839 : d_Fr_term_dr_out = 0
2840 : else ! include radiative force, Fr*kap_face/clight
2841 : if (k == 1) then
2842 : kap_face = s% opacity(k)
2843 : d_kap_dVol_00 = dK_dVol(i)
2844 : d_kap_dT_00 = dK_dT(i)
2845 : d_kap_dT_out = 0d0
2846 : d_kap_dr_in = dK_dr_in(i)
2847 : d_kap_dr_00 = dK_dr_00(i)
2848 : d_kap_dr_out = 0d0
2849 : else
2850 : kap_face = 0.5d0*(s% opacity(k) + s% opacity(k-1))
2851 : d_kap_dVol_00 = 0.5d0*dK_dVol(i)
2852 : d_kap_dT_00 = 0.5d0*dK_dT(i)
2853 : d_kap_dT_out = 0.5d0*dK_dT(i+1)
2854 : d_kap_dr_in = 0.5d0*dK_dr_in(i)
2855 : d_kap_dr_00 = 0.5d0*(dK_dr_00(i) + dK_dr_in(i+1))
2856 : d_kap_dr_out = 0.5d0*dK_dr_00(i+1)
2857 : end if
2858 : Fr_tw = WTR*s% Fr(k) + WTR1*s% Fr_start(k)
2859 : Fr_term = Fr_tw*kap_face/clight
2860 : d_Fr_term_dFr_00 = WTR*kap_face/clight
2861 : d_Fr_term_dVol_00 = Fr_term*d_kap_dVol_00/kap_face
2862 : d_Fr_term_dT_00 = Fr_term*d_kap_dT_00/kap_face
2863 : d_Fr_term_dT_out = Fr_term*d_kap_dT_out/kap_face
2864 : d_Fr_term_dr_in = Fr_term*d_kap_dr_in/kap_face
2865 : d_Fr_term_dr_00 = Fr_term*d_kap_dr_00/kap_face
2866 : d_Fr_term_dr_out = Fr_term*d_kap_dr_out/kap_face
2867 : end if
2868 :
2869 : residual = &
2870 : dvdt_factor*(s% v(k) - s% v_start(k))/dt &
2871 0 : + area*dXP_dm - grav - s% Uq(k) - Fr_term
2872 0 : HR(IR) = -residual
2873 :
2874 : !s% xtra1_array(k) = s% Pgas(k) + s% Prad(k)
2875 : !s% xtra2_array(k) = s% Vol(k)
2876 : !s% xtra3_array(k) = s% T(k)
2877 : !s% xtra4_array(k) = s% v(k)
2878 : !s% xtra5_array(k) = s% RSP_w(k)**2
2879 : !s% xtra6_array(k) = s% r(k)
2880 :
2881 0 : HD(i_r_dFr_00,IR) = - d_Fr_term_dFr_00
2882 :
2883 : HD(i_r_dT_in,IR) = &
2884 0 : - d_Uq_dT_in
2885 : HD(i_r_dT_00,IR) = &
2886 : + A_dm*(-dXP_00_dT_00) &
2887 : - d_Uq_dT_00 &
2888 0 : - d_Fr_term_dT_00
2889 : HD(i_r_dT_out,IR) = &
2890 : + A_dm*dXP_out_dT_out &
2891 : - d_Uq_dT_out &
2892 0 : - d_Fr_term_dT_out
2893 : HD(i_r_dT_out2,IR) = &
2894 0 : - d_Uq_dT_out2
2895 :
2896 0 : HD(i_r_dr_in2,IR) = - d_Uq_dr_in2
2897 : HD(i_r_dr_in,IR) = & !
2898 : + A_dm*(-dXP_00_dr_in) &
2899 : - d_Uq_dr_in &
2900 0 : - d_Fr_term_dr_in
2901 : HD(i_r_dr_00,IR) = &
2902 : dvdt_factor*dv_dr/dt &
2903 : + d_area_dr_00*dXP_dm &
2904 : + A_dm*(dXP_out_dr_00 - dXP_00_dr_00) &
2905 : + grav/s% r(k) &
2906 : - d_Uq_dr_00 &
2907 0 : - d_Fr_term_dr_00
2908 : HD(i_r_dr_00,IR) = &
2909 0 : HD(i_r_dr_00,IR) + 0.5d0*Uq1/R_00*(Chi_out - Chi_00)
2910 : HD(i_r_dr_out,IR) = &
2911 : + A_dm*dXP_out_dr_out &
2912 : - d_Uq_dr_out &
2913 0 : - d_Fr_term_dr_out
2914 : HD(i_r_dr_out2,IR) = & !
2915 0 : - d_Uq_dr_out2
2916 :
2917 : HD(i_r_der_in,IR) = &
2918 0 : - d_Uq_der_in
2919 : HD(i_r_der_00,IR) = &
2920 : + A_dm*(-dXP_00_der_00) &
2921 0 : - d_Uq_der_00
2922 : HD(i_r_der_out,IR) = &
2923 : + A_dm*dXP_out_der_out &
2924 0 : - d_Uq_der_out
2925 : HD(i_r_der_out2,IR) = &
2926 0 : - d_Uq_der_out2
2927 :
2928 0 : HD(i_r_dw_in,IR) = 0.d0
2929 0 : if (I <= IBOTOM .or. I == NZN) then
2930 0 : HD(i_r_dw_00,IR) = 0.d0
2931 : else
2932 : HD(i_r_dw_00,IR) = &
2933 : + A_dm*(-dXP_00_dw_00) &
2934 0 : - d_Uq_dw_00
2935 : end if
2936 0 : if (I <= IBOTOM - 1 .or. I >= NZN - 1) then
2937 0 : HD(i_r_dw_out,IR) = 0.d0
2938 : else
2939 : HD(i_r_dw_out,IR) = &
2940 : + A_dm*dXP_out_dw_out &
2941 0 : - d_Uq_dw_out
2942 : end if
2943 0 : HD(i_r_dw_out2,IR) = 0.0d0
2944 :
2945 : !call check_is_bad
2946 :
2947 :
2948 : !HD(i_r_dr_in2,IR) !
2949 : !HD(i_r_dr_in,IR) ! ok
2950 : !HD(i_r_dr_00,IR) ! ok
2951 : !HD(i_r_dr_out,IR) ! ok
2952 : !HD(i_r_dr_out2,IR) !
2953 :
2954 : !HD(i_r_dT_in,IR) ! 0
2955 : !HD(i_r_dT_00,IR) ! ok
2956 : !HD(i_r_dT_out,IR) ! need 1d-3 like for er
2957 : !HD(i_r_dT_out2,IR) ! 0
2958 :
2959 : ! NOTE: may need solver_test_partials_dx_0 = 1d-3 for er
2960 : !HD(i_r_der_in,IR) !
2961 : !HD(i_r_der_00,IR) ! 0
2962 : !HD(i_r_der_out,IR) ! 0
2963 : !HD(i_r_der_out2,IR) !
2964 :
2965 : !HD(i_r_dw_00,IR) !
2966 : !HD(i_r_dw_out,IR) !
2967 :
2968 : !test_partials = (k == s% solver_test_partials_k)
2969 0 : test_partials = .false.
2970 :
2971 0 : if (test_partials) then
2972 : s% solver_test_partials_val = residual
2973 : s% solver_test_partials_var = i_var_r
2974 : s% solver_test_partials_dval_dx = i_r_dr_00
2975 : write(*,*) 'acceleration_eqn', s% solver_test_partials_var
2976 : end if
2977 :
2978 : contains
2979 :
2980 : subroutine check_is_bad
2981 : include 'formats'
2982 : if (is_bad(residual)) then
2983 : !$omp critical (rsp_step_7)
2984 : write(*,2) 'residual', k, residual
2985 : write(*,2) 's% v(k)', k, s% v(k)
2986 : write(*,2) 's% v_start(k)', k, s% v_start(k)
2987 : write(*,2) 'area', k, area
2988 : write(*,2) 'dXP_dm', k, dXP_dm
2989 : write(*,2) 'XP_out', k, XP_out
2990 : write(*,2) 'XP_00', k, XP_00
2991 : write(*,2) 's% dm_bar(k)', k, s% dm_bar(k)
2992 : write(*,2) 'grav', k, grav
2993 : write(*,2) 's% Uq(k)', k, s% Uq(k)
2994 : write(*,2) 'Fr_term', k, Fr_term
2995 : write(*,2) 'dt', k, dt
2996 : call mesa_error(__FILE__,__LINE__,'acceleration_eqn')
2997 : !$omp end critical (rsp_step_7)
2998 : end if
2999 :
3000 : if (is_bad(HD(i_r_dr_in2,IR))) then
3001 : !$omp critical (rsp_step_8)
3002 : write(*,2) 'HD(i_r_dr_in2,IR)', k, HD(i_r_dr_in2,IR)
3003 : call mesa_error(__FILE__,__LINE__,'acceleration_eqn')
3004 : !$omp end critical (rsp_step_8)
3005 : end if
3006 :
3007 : if (is_bad(HD(i_r_dr_in,IR))) then
3008 : !$omp critical (rsp_step_9)
3009 : write(*,2) 'HD(i_r_dr_in,IR)', k, HD(i_r_dr_in,IR)
3010 : write(*,2) 'dXP_00_dr_in', k, dXP_00_dr_in
3011 : write(*,2) 'd_Uq_dr_in', k, d_Uq_dr_in
3012 : write(*,2) 'd_Fr_term_dr_in', k, d_Fr_term_dr_in
3013 : write(*,2) 'd_Chi_out_dr_in', k, d_Chi_out_dr_in
3014 : write(*,2) 'd_Chi_00_dr_in', k, d_Chi_00_dr_in
3015 : call mesa_error(__FILE__,__LINE__,'acceleration_eqn')
3016 : !$omp end critical (rsp_step_9)
3017 : end if
3018 :
3019 : if (is_bad(HD(i_r_dr_00,IR))) then
3020 : !$omp critical (rsp_step_10)
3021 : write(*,2) 'HD(i_r_dr_00,IR)', k, HD(i_r_dr_00,IR)
3022 : call mesa_error(__FILE__,__LINE__,'acceleration_eqn')
3023 : !$omp end critical (rsp_step_10)
3024 : end if
3025 :
3026 : if (is_bad(HD(i_r_dr_out,IR))) then
3027 : !$omp critical (rsp_step_11)
3028 : write(*,2) 'HD(i_r_dr_out,IR)', k, HD(i_r_dr_out,IR)
3029 : call mesa_error(__FILE__,__LINE__,'acceleration_eqn')
3030 : !$omp end critical (rsp_step_11)
3031 : end if
3032 :
3033 : if (is_bad(HD(i_r_dr_out2,IR))) then
3034 : !$omp critical (rsp_step_12)
3035 : write(*,2) 'HD(i_r_dr_out2,IR)', k, HD(i_r_dr_out2,IR)
3036 : call mesa_error(__FILE__,__LINE__,'acceleration_eqn')
3037 : !$omp end critical (rsp_step_12)
3038 : end if
3039 :
3040 : if (is_bad(HD(i_r_dT_in,IR))) then
3041 : !$omp critical (rsp_step_13)
3042 : write(*,2) 'HD(i_r_dT_in,IR)', k, HD(i_r_dT_in,IR)
3043 : call mesa_error(__FILE__,__LINE__,'acceleration_eqn')
3044 : !$omp end critical (rsp_step_13)
3045 : end if
3046 :
3047 : if (is_bad(HD(i_r_dT_00,IR))) then
3048 : !$omp critical (rsp_step_14)
3049 : write(*,2) 'HD(i_r_dT_00,IR)', k, HD(i_r_dT_00,IR)
3050 : call mesa_error(__FILE__,__LINE__,'acceleration_eqn')
3051 : !$omp end critical (rsp_step_14)
3052 : end if
3053 :
3054 : if (is_bad(HD(i_r_dT_out,IR))) then
3055 : !$omp critical (rsp_step_15)
3056 : write(*,2) 'HD(i_r_dT_out,IR)', k, HD(i_r_dT_out,IR)
3057 : call mesa_error(__FILE__,__LINE__,'acceleration_eqn')
3058 : !$omp end critical (rsp_step_15)
3059 : end if
3060 :
3061 : if (is_bad(HD(i_r_dT_out2,IR))) then
3062 : !$omp critical (rsp_step_16)
3063 : write(*,2) 'HD(i_r_dT_out2,IR)', k, HD(i_r_dT_out2,IR)
3064 : call mesa_error(__FILE__,__LINE__,'acceleration_eqn')
3065 : !$omp end critical (rsp_step_16)
3066 : end if
3067 :
3068 : if (is_bad(HD(i_r_der_in,IR))) then
3069 : !$omp critical (rsp_step_17)
3070 : write(*,2) 'HD(i_r_der_in,IR)', k, HD(i_r_der_in,IR)
3071 : call mesa_error(__FILE__,__LINE__,'acceleration_eqn')
3072 : !$omp end critical (rsp_step_17)
3073 : end if
3074 :
3075 : if (is_bad(HD(i_r_der_00,IR))) then
3076 : !$omp critical (rsp_step_18)
3077 : write(*,2) 'HD(i_r_der_00,IR)', k, HD(i_r_der_00,IR)
3078 : call mesa_error(__FILE__,__LINE__,'acceleration_eqn')
3079 : !$omp end critical (rsp_step_18)
3080 : end if
3081 :
3082 : if (is_bad(HD(i_r_der_out,IR))) then
3083 : !$omp critical (rsp_step_19)
3084 : write(*,2) 'HD(i_r_der_out,IR)', k, HD(i_r_der_out,IR)
3085 : call mesa_error(__FILE__,__LINE__,'acceleration_eqn')
3086 : !$omp end critical (rsp_step_19)
3087 : end if
3088 :
3089 : if (is_bad(HD(i_r_der_out2,IR))) then
3090 : !$omp critical (rsp_step_20)
3091 : write(*,2) 'HD(i_r_der_out2,IR)', k, HD(i_r_der_out2,IR)
3092 : call mesa_error(__FILE__,__LINE__,'acceleration_eqn')
3093 : !$omp end critical (rsp_step_20)
3094 : end if
3095 :
3096 : if (is_bad(HD(i_r_dw_00,IR))) then
3097 : !$omp critical (rsp_step_21)
3098 : write(*,2) 'HD(i_r_dw_00,IR)', k, HD(i_r_dw_00,IR)
3099 : call mesa_error(__FILE__,__LINE__,'acceleration_eqn')
3100 : !$omp end critical (rsp_step_21)
3101 : end if
3102 :
3103 : if (is_bad(HD(i_r_dw_out,IR))) then
3104 : !$omp critical (rsp_step_22)
3105 : write(*,2) 'HD(i_r_dw_out,IR)', k, HD(i_r_dw_out,IR)
3106 : call mesa_error(__FILE__,__LINE__,'acceleration_eqn')
3107 : !$omp end critical (rsp_step_22)
3108 : end if
3109 : end subroutine check_is_bad
3110 :
3111 : end subroutine acceleration_eqn
3112 :
3113 :
3114 0 : subroutine total_energy_eqn(s, i, P_surf, &
3115 : Lt_00, Lt_00_start, Lt_in, Lt_in_start, &
3116 : dLt_00_dr_00, dLt_00_dr_in, dLt_00_dr_out, &
3117 : dLt_00_dVol_00, dLt_00_dVol_out, &
3118 : dLt_00_dT_00, dLt_00_dT_out, &
3119 : dLt_00_der_00, dLt_00_der_out, &
3120 : dLt_00_dw_00, dLt_00_dw_out, &
3121 : dLt_in_dr_in, dLt_in_dr_in2, dLt_in_dr_00, &
3122 : dLt_in_dVol_in, dLt_in_dVol_00, &
3123 : dLt_in_dT_in, dLt_in_dT_00, &
3124 : dLt_in_der_in, dLt_in_der_00, &
3125 : dLt_in_dw_in, dLt_in_dw_00)
3126 : type (star_info), pointer :: s
3127 : integer, intent(in) :: i
3128 : real(dp), intent(in) :: P_surf, &
3129 : Lt_00, Lt_00_start, Lt_in, Lt_in_start, &
3130 : dLt_00_dr_00, dLt_00_dr_in, dLt_00_dr_out, &
3131 : dLt_00_dVol_00, dLt_00_dVol_out, &
3132 : dLt_00_dT_00, dLt_00_dT_out, &
3133 : dLt_00_der_00, dLt_00_der_out, &
3134 : dLt_00_dw_00, dLt_00_dw_out, &
3135 : dLt_in_dr_in, dLt_in_dr_in2, dLt_in_dr_00, &
3136 : dLt_in_dVol_in, dLt_in_dVol_00, &
3137 : dLt_in_dT_in, dLt_in_dT_00, &
3138 : dLt_in_der_in, dLt_in_der_00, &
3139 : dLt_in_dw_in, dLt_in_dw_00
3140 : integer :: IT, k
3141 0 : real(dp) :: dt, dm, residual, erad, erad_tw, DV, dt_div_dm, &
3142 0 : area_00, area_00_start, area_in, area_in_start, &
3143 0 : L_00, Lr_00, Lr_00_start, d_Lr_00_dFr_00, d_Lr_00_dr_00, &
3144 0 : L_in, Lr_in, Lr_in_start, d_Lr_in_dFr_in, d_Lr_in_dr_in, &
3145 0 : Lc_00, Lc_00_start, Lc_in, Lc_in_start, &
3146 : dLc_in_dr_in, dLc_in_dr_in2, dLc_in_dr_00, &
3147 : dLc_in_dVol_in, dLc_in_dVol_00, &
3148 : dLc_in_dT_in, dLc_in_dT_00, &
3149 : dLc_in_der_in, dLc_in_der_00, &
3150 : dLc_in_dw_in, dLc_in_dw_00, &
3151 : dLc_00_dr_00, dLc_00_dr_in, dLc_00_dr_out, &
3152 : dLc_00_dVol_00, dLc_00_dVol_out, &
3153 : dLc_00_dT_00, dLc_00_dT_out, &
3154 : dLc_00_der_00, dLc_00_der_out, &
3155 : dLc_00_dw_00, dLc_00_dw_out, &
3156 : XP_00, dXP_00_dr_00, dXP_00_dr_in, &
3157 : dXP_00_dVol_00, dXP_00_dT_00, dXP_00_der_00, dXP_00_dw_00, &
3158 0 : u_div_r, d_u_div_r_dr_00, d_u_div_r_dr_in, u_div_r_factor
3159 : logical :: test_partials
3160 :
3161 : include 'formats'
3162 :
3163 0 : k = NZN+1-i
3164 :
3165 0 : IT = i_var_T + NV*(i-1)
3166 0 : HD(1:LD_HD,IT) = 0.d0
3167 :
3168 : call calc_Lc(s, i-1, Lc_in, &
3169 : dLc_in_dr_in2, dLc_in_dr_in, dLc_in_dr_00, &
3170 : dLc_in_dVol_in, dLc_in_dVol_00, &
3171 : dLc_in_dT_in, dLc_in_dT_00, &
3172 : dLc_in_der_in, dLc_in_der_00, &
3173 0 : dLc_in_dw_in, dLc_in_dw_00)
3174 : call calc_Lc(s, i, Lc_00, &
3175 : dLc_00_dr_in, dLc_00_dr_00, dLc_00_dr_out, &
3176 : dLc_00_dVol_00, dLc_00_dVol_out, &
3177 : dLc_00_dT_00, dLc_00_dT_out, &
3178 : dLc_00_der_00, dLc_00_der_out, &
3179 0 : dLc_00_dw_00, dLc_00_dw_out)
3180 :
3181 0 : if (I == NZN) then
3182 : Lc_00_start = 0.d0
3183 : else
3184 0 : Lc_00_start = s% Lc_start(k)
3185 : end if
3186 0 : if (i == 1) then
3187 : Lc_in_start = 0
3188 : else
3189 0 : Lc_in_start = s% Lc_start(k+1)
3190 : end if
3191 :
3192 0 : area_00 = 4d0*pi*s% r(k)**2
3193 0 : Lr_00 = s% Fr(k)*area_00
3194 0 : d_Lr_00_dFr_00 = area_00
3195 0 : d_Lr_00_dr_00 = 2d0*Lr_00/s% r(k)
3196 0 : area_00_start = 4d0*pi*s% r_start(k)**2
3197 0 : Lr_00_start = s% Fr_start(k)*area_00_start
3198 0 : if (i == 1) then
3199 0 : if (s% RSP_hydro_only) then
3200 : Lr_in = 0d0
3201 : else
3202 0 : Lr_in = s% L_center
3203 : end if
3204 : d_Lr_in_dFr_in = 0d0
3205 : d_Lr_in_dr_in = 0d0
3206 : Lr_in_start = Lr_in
3207 : else
3208 0 : area_in = 4d0*pi*s% r(k+1)**2
3209 0 : Lr_in = s% Fr(k+1)*area_in
3210 0 : d_Lr_in_dFr_in = area_in
3211 0 : d_Lr_in_dr_in = 2d0*Lr_in/s% r(k+1)
3212 0 : area_in_start = 4d0*pi*s% r_start(k+1)**2
3213 0 : Lr_in_start = s% Fr_start(k+1)*area_in_start
3214 : end if
3215 :
3216 : L_00 = &
3217 : WTR*Lr_00 + WTC*Lc_00 + WTT*Lt_00 + &
3218 0 : WTR1*Lr_00_start + WTC1*Lc_00_start + WTT1*Lt_00_start
3219 : L_in = &
3220 : WTR*Lr_in + WTC*Lc_in + WTT*Lt_in + &
3221 0 : WTR1*Lr_in_start + WTC1*Lc_in_start + WTT1*Lt_in_start
3222 :
3223 0 : dt = s% dt
3224 0 : dm = s% dm(k)
3225 0 : dt_div_dm = dt/dm
3226 0 : DV = s% Vol(k) - s% Vol_start(k)
3227 :
3228 0 : if (s% f_Edd(k) == f_Edd_isotropic .or. k == NZN) then
3229 : u_div_r = 0d0
3230 : d_u_div_r_dr_00 = 0d0
3231 : d_u_div_r_dr_in = 0d0
3232 : u_div_r_factor = 0d0
3233 : else
3234 0 : u_div_r = 0.5d0*(s% v(k)/s% r(k) + s% v(k+1)/s% r(k+1))
3235 0 : d_u_div_r_dr_00 = 0.5d0*(2d0/dt - s% v(k)/s% r(k))/s% r(k)
3236 0 : d_u_div_r_dr_in = 0.5d0*(2d0/dt - s% v(k+1)/s% r(k+1))/s% r(k+1)
3237 0 : u_div_r_factor = dt*(1d0 - 3d0*s% f_Edd(k))
3238 : end if
3239 :
3240 0 : erad = s% erad(k)
3241 0 : erad_tw = THETAE*erad + THETAE1*s% erad_start(k)
3242 :
3243 : call rsp_calc_XP(s, P_surf, i, .true., &
3244 : XP_00, dXP_00_dVol_00, dXP_00_dT_00, dXP_00_der_00, &
3245 0 : dXP_00_dw_00, dXP_00_dr_in, dXP_00_dr_00)
3246 :
3247 : residual = &
3248 : s% egas(k) - s% egas_start(k) &
3249 : + erad - s% erad_start(k) &
3250 : + s% RSP_w(k)**2 - s% RSP_w_start(k)**2 &
3251 : + dt_div_dm*(L_00 - L_in) &
3252 : + XP_00*DV &
3253 : + erad_tw*u_div_r_factor*u_div_r &
3254 0 : - dt*s% Eq(k)
3255 :
3256 0 : s% ergs_error(k) = s% dm(k)*residual
3257 :
3258 0 : HR(IT) = -residual
3259 :
3260 0 : HD(i_T_dFr_in,IT) = -dt_div_dm*WTR*d_Lr_in_dFr_in
3261 0 : HD(i_T_dFr_00,IT) = dt_div_dm*WTR*d_Lr_00_dFr_00
3262 :
3263 : HD(i_T_dr_in2,IT) = &
3264 : - dt_div_dm*WTC*dLc_in_dr_in2 &
3265 : - dt_div_dm*WTT*dLt_in_dr_in2 &
3266 0 : - dt*dEq_dr_in2(I)
3267 : HD(i_T_dr_in,IT) = &
3268 : + d_egas_dr_in(i) &
3269 : - dt_div_dm*WTR*d_Lr_in_dr_in &
3270 : + dt_div_dm*WTC*(dLc_00_dr_in - dLc_in_dr_in) &
3271 : + dt_div_dm*WTT*(dLt_00_dr_in - dLt_in_dr_in) &
3272 : + dVol_dr_in(I)*XP_00 &
3273 : + DV*dXP_00_dr_in &
3274 : + erad_tw*u_div_r_factor*d_u_div_r_dr_in &
3275 0 : - dt*dEq_dr_in(I)
3276 : HD(i_T_dr_00,IT) = &
3277 : + d_egas_dr_00(i) &
3278 : + dt_div_dm*WTR*d_Lr_00_dr_00 &
3279 : + dt_div_dm*WTC*(dLc_00_dr_00 - dLc_in_dr_00) &
3280 : + dt_div_dm*WTT*(dLt_00_dr_00 - dLt_in_dr_00) &
3281 : + dVol_dr_00(I)*XP_00 &
3282 : + DV*dXP_00_dr_00 &
3283 : + erad_tw*u_div_r_factor*d_u_div_r_dr_00 &
3284 0 : - dt*dEq_dr_00(I)
3285 : HD(i_T_dr_out,IT) = &
3286 : + dt_div_dm*WTC*dLc_00_dr_out &
3287 : + dt_div_dm*WTT*dLt_00_dr_out &
3288 0 : - dt*dEq_dr_out(I)
3289 :
3290 : HD(i_T_dT_in,IT) = &
3291 : - dt_div_dm*WTC*dLc_in_dT_in &
3292 : - dt_div_dm*WTT*dLt_in_dT_in &
3293 0 : - dt*dEq_dT_in(I)
3294 : HD(i_T_dT_00,IT) = &
3295 : d_egas_dT(i) &
3296 : + DV*dXP_00_dT_00 &
3297 : + dt_div_dm*WTC*(dLc_00_dT_00 - dLc_in_dT_00) &
3298 : + dt_div_dm*WTT*(dLt_00_dT_00 - dLt_in_dT_00) &
3299 0 : - dt*dEq_dT_00(I)
3300 : HD(i_T_dT_out,IT) = & !
3301 : + dt_div_dm*WTC*dLc_00_dT_out &
3302 : + dt_div_dm*WTT*dLt_00_dT_out &
3303 0 : - dt*dEq_dT_out(I)
3304 :
3305 : HD(i_T_der_in,IT) = &
3306 : - dt_div_dm*WTC*dLc_in_der_in &
3307 : - dt_div_dm*WTT*dLt_in_der_in &
3308 0 : - dt*dEq_der_in(I)
3309 : HD(i_T_der_00,IT) = &
3310 : 1d0 &
3311 : + DV*dXP_00_der_00 &
3312 : + dt_div_dm*WTC*(dLc_00_der_00 - dLc_in_der_00) &
3313 : + dt_div_dm*WTT*(dLt_00_der_00 - dLt_in_der_00) &
3314 : + THETAE*u_div_r_factor*u_div_r &
3315 0 : - dt*dEq_der_00(I)
3316 : HD(i_T_der_out,IT) = & !
3317 : + dt_div_dm*WTC*dLc_00_der_out &
3318 : + dt_div_dm*WTT*dLt_00_der_out &
3319 0 : - dt*dEq_der_out(I)
3320 :
3321 0 : if (I <= IBOTOM + 1) then
3322 0 : HD(i_T_dw_in,IT) = 0.d0
3323 : else
3324 : HD(i_T_dw_in,IT) = &
3325 : - dt_div_dm*WTC*dLc_in_dw_in &
3326 0 : - dt_div_dm*WTT*dLt_in_dw_in
3327 : end if
3328 0 : if (I <= IBOTOM .or. I == NZN) then
3329 0 : HD(i_T_dw_00,IT) = 0.d0
3330 : else
3331 : HD(i_T_dw_00,IT) = &
3332 : 2.d0*s% RSP_w(k) &
3333 : + dt_div_dm*WTC*(dLc_00_dw_00 - dLc_in_dw_00) &
3334 : + dt_div_dm*WTT*(dLt_00_dw_00 - dLt_in_dw_00) &
3335 : + DV*dXP_00_dw_00 &
3336 0 : - dt*dEq_dw_00(I)
3337 : end if
3338 0 : if (I <= IBOTOM - 1 .or. I >= NZN - 1) then
3339 0 : HD(i_T_dw_out,IT) = 0.d0
3340 : else
3341 : HD(i_T_dw_out,IT) = &
3342 : dt_div_dm*WTT*dLt_00_dw_out &
3343 0 : + dt_div_dm*WTC*dLc_00_dw_out
3344 : end if
3345 :
3346 0 : if (i == -6 .and. s% model_number == s% max_model_number) then
3347 0 : write(*,5) 'HD(i_T_dw_00,IT)', k, i, iter, s% model_number, HD(i_T_dw_00,IT)
3348 0 : write(*,5) 's% RSP_w(k)', k, i, iter, s% model_number, s% RSP_w(k)
3349 0 : write(*,5) 'dt_div_dm', k, i, iter, s% model_number, dt_div_dm
3350 0 : write(*,5) 'WTC', k, i, iter, s% model_number, WTC
3351 0 : write(*,5) 'WTT', k, i, iter, s% model_number, WTT
3352 0 : write(*,5) 'DV', k, i, iter, s% model_number, DV
3353 0 : write(*,5) 'dt', k, i, iter, s% model_number, dt
3354 0 : write(*,5) 'dLc_00_dw_00', k, i, iter, s% model_number, dLc_00_dw_00
3355 0 : write(*,5) 'dLc_in_dw_00', k, i, iter, s% model_number, dLc_in_dw_00
3356 0 : write(*,5) 'dLt_00_dw_00', k, i, iter, s% model_number, dLt_00_dw_00
3357 0 : write(*,5) 'dLt_in_dw_00', k, i, iter, s% model_number, dLt_in_dw_00
3358 0 : write(*,5) 'dXP_00_dw_00', k, i, iter, s% model_number, dXP_00_dw_00
3359 0 : write(*,5) 'dEq_dw_00(I)', k, i, iter, s% model_number, dEq_dw_00(I)
3360 : end if
3361 :
3362 : !HD(i_T_dr_in2,IT) !
3363 : !HD(i_T_dr_in,IT) ! ok
3364 : !HD(i_T_dr_00,IT) ! ok
3365 : !HD(i_T_dr_out,IT) !
3366 :
3367 : !HD(i_T_dT_in,IT) ! 0
3368 : !HD(i_T_dT_00,IT) ! ok
3369 : !HD(i_T_dT_out,IT) ! 0
3370 :
3371 : !HD(i_T_der_in,IT) !
3372 : !HD(i_T_der_00,IT) ! ok
3373 : !HD(i_T_der_out,IT) !
3374 :
3375 : !HD(i_T_dw_in,IT) !
3376 : !HD(i_T_dw_00,IT) !
3377 : !HD(i_T_dw_out,IT) !
3378 :
3379 : !test_partials = (k+1 == s% solver_test_partials_k)
3380 0 : test_partials = .false.
3381 :
3382 : if (test_partials) then
3383 : s% solver_test_partials_val = residual
3384 : s% solver_test_partials_var = i_var_r
3385 : s% solver_test_partials_dval_dx = HD(i_T_dr_00,IT)
3386 : write(*,*) 'total_energy_eqn', s% solver_test_partials_var
3387 : end if
3388 :
3389 0 : end subroutine total_energy_eqn
3390 :
3391 :
3392 0 : subroutine turbulent_energy_eqn(s, i, &
3393 : Lt_00, Lt_00_start, Lt_in, Lt_in_start, &
3394 : dLt_00_dr_00, dLt_00_dr_in, dLt_00_dr_out, &
3395 : dLt_00_dVol_00, dLt_00_dVol_out, &
3396 : dLt_00_dT_00, dLt_00_dT_out, &
3397 : dLt_00_der_00, dLt_00_der_out, &
3398 : dLt_00_dw_00, dLt_00_dw_out, &
3399 : dLt_in_dr_in, dLt_in_dr_in2, dLt_in_dr_00, &
3400 : dLt_in_dVol_in, dLt_in_dVol_00, &
3401 : dLt_in_dT_in, dLt_in_dT_00, &
3402 : dLt_in_der_in, dLt_in_der_00, &
3403 : dLt_in_dw_in, dLt_in_dw_00)
3404 : type (star_info), pointer :: s
3405 : integer, intent(in) :: i
3406 : real(dp), intent(in) :: &
3407 : Lt_00, Lt_00_start, Lt_in, Lt_in_start, &
3408 : dLt_00_dr_00, dLt_00_dr_in, dLt_00_dr_out, &
3409 : dLt_00_dVol_00, dLt_00_dVol_out, &
3410 : dLt_00_dT_00, dLt_00_dT_out, &
3411 : dLt_00_der_00, dLt_00_der_out, &
3412 : dLt_00_dw_00, dLt_00_dw_out, &
3413 : dLt_in_dr_in, dLt_in_dr_in2, dLt_in_dr_00, &
3414 : dLt_in_dVol_in, dLt_in_dVol_00, &
3415 : dLt_in_dT_in, dLt_in_dT_00, &
3416 : dLt_in_der_in, dLt_in_der_00, &
3417 : dLt_in_dw_in, dLt_in_dw_00
3418 : integer :: IW, k
3419 0 : real(dp) :: dt, dm, residual, Ptrb_tw, DV, dt_div_dm, L_00, L_in
3420 : logical :: test_partials
3421 :
3422 : include 'formats'
3423 :
3424 0 : k = NZN+1-i
3425 :
3426 0 : IW = i_var_w + NV*(i-1)
3427 0 : HD(1:LD_HD,IW) = 0.d0
3428 :
3429 0 : if (ALFA == 0d0 .or. I <= IBOTOM .or. I == NZN) then
3430 0 : HD(i_w_dw_00,IW) = 1.d0
3431 0 : HR(IW) = 0.d0
3432 0 : return
3433 : end if
3434 :
3435 0 : L_00 = WTT*Lt_00 + WTT1*Lt_00_start
3436 0 : L_in = WTT*Lt_in + WTT1*Lt_in_start
3437 :
3438 0 : dt = s% dt
3439 0 : dm = s% dm(k)
3440 0 : dt_div_dm = dt/dm
3441 0 : DV = s% Vol(k) - s% Vol_start(k)
3442 0 : Ptrb_tw = THETAT*s% Ptrb(k) + THETAT1*s% Ptrb_start(k)
3443 :
3444 : residual = &
3445 : (s% RSP_w(k)**2 - s% RSP_w_start(k)**2) &
3446 : + dt_div_dm*(L_00 - L_in) &
3447 : + DV*Ptrb_tw &
3448 0 : - dt*(GAM*s% COUPL(k) + GAM1*s% COUPL_start(k) + s% Eq(k))
3449 0 : HR(IW) = -residual
3450 :
3451 0 : if (k==-35) then
3452 0 : write(*,3) 'RSP w dEt PdV dtC dtEq', k, iter, &
3453 0 : s% RSP_w(k), s% RSP_w(k)**2 - s% RSP_w_start(k)**2, DV*Ptrb_tw, &
3454 0 : dt*(GAM*s% COUPL(k) + GAM1*s% COUPL_start(k)), dt*s% Eq(k)
3455 : !write(*,2) 'RSP w COUPL SOURCE DAMP DAMPR', k, &
3456 : ! s% RSP_w(k), s% COUPL(k), s% SOURCE(k), s% DAMP(k), s% DAMPR(k)
3457 : !write(*,2) 'RSP w SOURCE PII/Hp P*QQ_div_Cp P T', k, &
3458 : ! s% RSP_w(k), s% SOURCE(k), &
3459 : ! 0.5d0*(s% PII(k)/s% Hp_face(k) + s% PII(k+1)/s% Hp_face(k+1)), &
3460 : ! (s% Pgas(k) + s% Prad(k))*s% QQ(k)/s% Cp(k), s% Pgas(k) + s% Prad(k), s% T(k)
3461 : !write(*,2) 'RSP PII_00 PII_p1 Hp_00 Hp_p1', k, &
3462 : ! s% PII(k), s% PII(k+1), s% Hp_face(k), s% Hp_face(k+1)
3463 : end if
3464 :
3465 0 : HD(i_w_dw_in2,IW) = 0.d0
3466 0 : HD(i_w_dw_in,IW) = - dt_div_dm*WTT*dLt_in_dw_in
3467 : HD(i_w_dw_00,IW) = &
3468 : 2.d0*s% RSP_w(k) &
3469 : - dt*GAM*dC_dw_00(I) &
3470 : - dt*dEq_dw_00(I) &
3471 : + dt_div_dm*WTT*(dLt_00_dw_00 - dLt_in_dw_00) &
3472 0 : + DV*THETAT*dPtrb_dw_00(I)
3473 0 : HD(i_w_dw_out,IW) = dt_div_dm*WTT*dLt_00_dw_out
3474 0 : HD(i_w_dw_out2,IW) = 0.d0
3475 :
3476 : HD(i_w_dr_in2,IW) = &
3477 : - dt*GAM*dC_dr_in2(I) &
3478 : - dt_div_dm*WTT*dLt_in_dr_in2 &
3479 0 : - dt*dEq_dr_in2(I)
3480 : HD(i_w_dr_in,IW) = &
3481 : - dt*GAM*dC_dr_in(I) &
3482 : - dt*dEq_dr_in(I) &
3483 : + dt_div_dm*WTT*(dLt_00_dr_in - dLt_in_dr_in) &
3484 : + DV*THETAT*dPtrb_dr_in(I) &
3485 0 : + (THETAT*s% Ptrb(k) + THETAT1*s% Ptrb_start(k))*dVol_dr_in(I)
3486 : HD(i_w_dr_00,IW) = &
3487 : - dt*GAM*dC_dr_00(I) &
3488 : - dt*dEq_dr_00(I) &
3489 : + dt_div_dm*WTT*(dLt_00_dr_00 - dLt_in_dr_00) &
3490 : + DV*THETAT*dPtrb_dr_00(I) &
3491 0 : + (THETAT*s% Ptrb(k) + THETAT1*s% Ptrb_start(k))*dVol_dr_00(I)
3492 : HD(i_w_dr_out,IW) = &
3493 : - dt*GAM*dC_dr_out(I) &
3494 : + dt_div_dm*WTT*dLt_00_dr_out &
3495 0 : - dt*dEq_dr_out(I)
3496 :
3497 0 : if (I <= IBOTOM + 1) then
3498 0 : HD(i_w_dT_in,IW) = 0.d0
3499 : else
3500 : HD(i_w_dT_in,IW) = &
3501 : - dt_div_dm*WTT*dLt_in_dT_in &
3502 : - dt*dEq_dT_in(I) &
3503 0 : - dt*GAM*dC_dT_in(I)
3504 : end if
3505 : HD(i_w_dT_00,IW) = &
3506 : - dt*GAM*dC_dT_00(I) &
3507 : - dt*dEq_dT_00(I) &
3508 0 : + dt_div_dm*WTT*(dLt_00_dT_00 - dLt_in_dT_00)
3509 0 : if (I <= IBOTOM - 1 .or. I >= NZN - 1) then
3510 0 : HD(i_w_dT_out,IW) = 0.d0
3511 : else
3512 : HD(i_w_dT_out,IW) = &
3513 : - dt*GAM*dC_dT_out(I) &
3514 : + dt_div_dm*WTT*dLt_00_dT_out &
3515 0 : - dt*dEq_dT_out(I)
3516 : end if
3517 :
3518 0 : if (I <= IBOTOM + 1) then
3519 0 : HD(i_w_der_in,IW) = 0.d0
3520 : else
3521 : HD(i_w_der_in,IW) = &
3522 : - dt_div_dm*WTT*dLt_in_der_in &
3523 : - dt*dEq_der_in(I) &
3524 0 : - dt*GAM*dC_der_in(I)
3525 : end if
3526 : HD(i_w_der_00,IW) = &
3527 : - dt*GAM*dC_der_00(I) &
3528 : - dt*dEq_der_00(I) &
3529 0 : + dt_div_dm*WTT*(dLt_00_der_00 - dLt_in_der_00)
3530 0 : if (I <= IBOTOM - 1 .or. I >= NZN - 1) then
3531 0 : HD(i_w_der_out,IW) = 0.d0
3532 : else
3533 : HD(i_w_der_out,IW) = &
3534 : - dt*GAM*dC_der_out(I) &
3535 : + dt_div_dm*WTT*dLt_00_der_out &
3536 0 : - dt*dEq_der_out(I)
3537 : end if
3538 :
3539 : !HD(i_w_dw_in,IW) !
3540 : !HD(i_w_dw_00,IW) !
3541 : !HD(i_w_dw_out,IW) !
3542 :
3543 : !HD(i_w_dr_in2,IW) !
3544 : !HD(i_w_dr_in,IW) !
3545 : !HD(i_w_dr_00,IW) !
3546 : !HD(i_w_dr_out,IW) !
3547 :
3548 : !HD(i_w_dT_in,IW) !
3549 : !HD(i_w_dT_00,IW) !
3550 : !HD(i_w_dT_out,IW) !
3551 :
3552 : !HD(i_w_der_in,IW) !
3553 : !HD(i_w_der_00,IW) !
3554 : !HD(i_w_der_out,IW) !
3555 :
3556 : !test_partials = (k+1 == s% solver_test_partials_k)
3557 0 : test_partials = .false.
3558 :
3559 : if (test_partials) then
3560 : s% solver_test_partials_val = residual
3561 : s% solver_test_partials_var = i_var_r
3562 : s% solver_test_partials_dval_dx = HD(i_w_dr_00,IW)
3563 : write(*,*) 'turbulent_energy_eqn', s% solver_test_partials_var
3564 : end if
3565 :
3566 : end subroutine turbulent_energy_eqn
3567 :
3568 :
3569 0 : subroutine erad_eqn(s, i)
3570 : use const_def, only: crad, clight
3571 : type (star_info), pointer :: s
3572 : integer, intent(in) :: i
3573 : include 'formats'
3574 :
3575 0 : call T_form_of_erad_eqn(s, i)
3576 :
3577 0 : end subroutine erad_eqn
3578 :
3579 :
3580 0 : subroutine Fr_eqn(s, i)
3581 : use const_def, only: clight
3582 : type (star_info), pointer :: s
3583 : integer, intent(in) :: i
3584 : integer :: IL, k
3585 0 : real(dp) :: residual
3586 :
3587 : include 'formats'
3588 :
3589 0 : if (s% RSP_hydro_only) then
3590 0 : k = NZN+1-i
3591 0 : IL = i_var_Fr + NV*(i-1)
3592 0 : HD(1:LD_HD,IL) = 0d0
3593 0 : residual = - s% Fr(k) ! want Fr = 0d0
3594 0 : HR(IL) = -residual
3595 0 : HD(i_Fr_dFr_00,IL) = -1d0
3596 0 : return
3597 : end if
3598 :
3599 0 : call T_form_of_Fr_eqn(s,i)
3600 :
3601 : end subroutine Fr_eqn
3602 :
3603 :
3604 : subroutine d_Prad_dm_Fr_eqn(s,i)
3605 : type (star_info), pointer :: s
3606 : integer, intent(in) :: i
3607 : integer :: IL, k
3608 : real(dp) :: residual, Fr_00, &
3609 : dFr_dr_out, dFr_dr_00, dFr_dr_in, &
3610 : dFr_dVol_out, dFr_dVol_00, &
3611 : dFr_dT_out, dFr_dT_00, &
3612 : dFr_der_out, dFr_der_00
3613 :
3614 : logical :: test_partials
3615 :
3616 : include 'formats'
3617 :
3618 : k = NZN+1-i
3619 : IL = i_var_Fr + NV*(i-1)
3620 : HD(1:LD_HD,IL) = 0d0
3621 :
3622 : call calc_Fr(s, i, Fr_00, &
3623 : dFr_dr_out, dFr_dr_00, dFr_dr_in, &
3624 : dFr_dVol_out, dFr_dVol_00, &
3625 : dFr_dT_out, dFr_dT_00, &
3626 : dFr_der_out, dFr_der_00)
3627 :
3628 : residual = Fr_00 - s% Fr(k)
3629 : HR(IL) = -residual
3630 :
3631 : HD(i_Fr_der_00,IL) = dFr_der_00
3632 : HD(i_Fr_der_out,IL) = dFr_der_out
3633 : HD(i_Fr_dr_in,IL) = dFr_dr_in !
3634 : HD(i_Fr_dr_00,IL) = dFr_dr_00 !
3635 : HD(i_Fr_dr_out,IL) = dFr_dr_out !
3636 : HD(i_Fr_dT_00,IL) = dFr_dT_00
3637 : HD(i_Fr_dT_out,IL) = dFr_dT_out
3638 : HD(i_Fr_dFr_00,IL) = -1d0
3639 :
3640 : !test_partials = (k == s% solver_test_partials_k)
3641 : test_partials = .false.
3642 :
3643 : if (test_partials) then
3644 : s% solver_test_partials_val = residual
3645 : s% solver_test_partials_var = i_var_er
3646 : s% solver_test_partials_dval_dx = dFr_der_00
3647 : write(*,*) 'd_Prad_dm_Fr_eqn', s% solver_test_partials_var
3648 : end if
3649 :
3650 : end subroutine d_Prad_dm_Fr_eqn
3651 :
3652 :
3653 0 : subroutine T_form_of_erad_eqn(s,i)
3654 : type (star_info), pointer :: s
3655 : integer, intent(in) :: i
3656 : integer :: IE, k
3657 0 : real(dp) :: T, V, erad_expected, residual
3658 : logical :: test_partials
3659 :
3660 : include 'formats'
3661 :
3662 0 : k = NZN+1-i
3663 0 : IE = i_var_er + NV*(i-1)
3664 0 : HD(1:LD_HD,IE) = 0.d0
3665 :
3666 0 : T = s% T(k)
3667 0 : V = s% Vol(k)
3668 0 : erad_expected = crad*T**4*V ! ergs/gm
3669 0 : residual = erad_expected - s% erad(k)
3670 :
3671 0 : HR(IE) = -residual
3672 :
3673 0 : HD(i_er_der_00,IE) = -1d0
3674 0 : HD(i_er_dT_00,IE) = 4d0*crad*T**3*V
3675 0 : HD(i_er_dr_00,IE) = crad*T**4*dVol_dr_00(i)
3676 0 : HD(i_er_dr_in,IE) = crad*T**4*dVol_dr_in(i)
3677 :
3678 : !test_partials = (k == s% solver_test_partials_k)
3679 0 : test_partials = .false.
3680 :
3681 : if (test_partials) then
3682 : s% solver_test_partials_val = residual
3683 : s% solver_test_partials_var = i_var_T
3684 : s% solver_test_partials_dval_dx = HD(i_er_dT_00,IE)
3685 : write(*,*) 'T_form_of_erad_eqn', s% solver_test_partials_var
3686 : end if
3687 0 : end subroutine T_form_of_erad_eqn
3688 :
3689 :
3690 0 : subroutine T_form_of_calc_Fr(s, i, Fr_00, & !rs Stellingwerf 1975, Appendix A
3691 : dFr_dr_out, dFr_dr_in, dFr_dr_00, &
3692 : dFr_dT_out, dFr_dT_00, dFr_dVol_00)
3693 : type (star_info), pointer :: s
3694 : integer, intent(in) :: I
3695 : real(dp), intent(out) :: &
3696 : Fr_00, dFr_dr_out, dFr_dr_in, &
3697 : dFr_dr_00, dFr_dT_out, dFr_dT_00, dFr_dVol_00
3698 0 : real(dp) :: dFr_dK_00, dFr_dK_out, W, WP, BW, BK, T1, T2, T3
3699 : integer :: k
3700 : logical :: test_partials
3701 : include 'formats'
3702 0 : k = NZN+1-i
3703 0 : if (i < 1) then
3704 0 : if (s% RSP_hydro_only) then
3705 0 : Fr_00 = 0d0
3706 : else
3707 0 : Fr_00 = s% L_center
3708 : end if
3709 0 : Fr_00 = Fr_00/(4d0*pi*s% r_center**2)
3710 0 : dFr_dr_00 = 0
3711 0 : dFr_dT_00 = 0
3712 0 : dFr_dK_00 = 0
3713 0 : dFr_dK_out = 0
3714 0 : dFr_dr_out = 0
3715 0 : dFr_dr_in = 0
3716 0 : dFr_dT_out = 0
3717 0 : else if (i == NZN) then
3718 0 : Fr_00 = 2d0*SIG*s% T(k)**4 !EDDI
3719 0 : dFr_dT_00 = 4.d0*Fr_00/s% T(k)
3720 0 : dFr_dK_00 = 0
3721 0 : dFr_dK_out = 0
3722 0 : dFr_dr_out = 0
3723 0 : dFr_dr_in = 0
3724 0 : dFr_dr_00 = 0
3725 0 : dFr_dT_out = 0
3726 : else
3727 0 : Fr_00 = 0d0
3728 0 : dFr_dr_00 = 0d0
3729 0 : dFr_dT_00 = 0d0
3730 0 : dFr_dK_00 = 0d0
3731 0 : dFr_dK_out = 0d0
3732 0 : dFr_dr_out = 0d0
3733 0 : dFr_dr_in = 0d0
3734 0 : dFr_dT_out = 0d0
3735 0 : W = s% T(k)**4
3736 0 : WP = s% T(k-1)**4
3737 0 : BW = 4d0*(s% lnT(k-1) - s% lnT(k)) ! log(WP/W)
3738 0 : if (abs(BW) < 1d-20) return
3739 0 : BK = log(s% opacity(k-1)/s% opacity(k))
3740 0 : if (abs(1.d0 - BK/BW) < 1d-15 .or. abs(BW - BK) < 1d-15) return
3741 0 : T1 = - CL*s% r(k)**2/(4d0*pi*s% dm_bar(k)) ! CL = 4d0*(4d0*PI)**2*SIG/3d0
3742 0 : T2 = (WP/s% opacity(k-1) - W/s% opacity(k))/(1.d0 - BK/BW)
3743 0 : Fr_00 = T1*T2
3744 0 : T3 = T1/(BW - BK)
3745 : !rs radiative luminosity derivatives
3746 0 : dFr_dK_00 = (T3/s% opacity(k))*(W*BW/s% opacity(k) - T2)
3747 0 : dFr_dK_out = -(T3/s% opacity(k-1))*(WP*BW/s% opacity(k-1) - T2)
3748 0 : dFr_dr_out = dFr_dK_out*dK_dr_00(i+1) !
3749 0 : dFr_dr_in = dFr_dK_00*dK_dr_in(I) !
3750 : dFr_dr_00 = 2d0*Fr_00/s% r(k) & !
3751 : + dFr_dK_00*dK_dr_00(I) &
3752 0 : + dFr_dK_out*dK_dr_in(i+1)
3753 : dFr_dT_out = & !
3754 : 4.d0*(T3/s% T(k-1))*(WP*BW/s% opacity(k-1) &
3755 0 : - T2*BK/BW) + dFr_dK_out*dK_dT(i+1)
3756 : dFr_dT_00 = & !
3757 : - 4.d0*(T3/s% T(k))*(W*BW/s% opacity(k) - T2*BK/BW) &
3758 0 : + dFr_dK_00*dK_dT(I)
3759 :
3760 : if (call_is_bad) then
3761 : if (is_bad(dFr_dT_out + dFr_dT_00)) then
3762 : write(*,3) 'dFr_dT_out', i, k, dFr_dT_out
3763 : write(*,3) 'dFr_dT_00', i, k, dFr_dT_00
3764 : write(*,3) 's% T(k-1)', i, k, s% T(k-1)
3765 : write(*,3) 's% T(k)', i, k, s% T(k)
3766 : write(*,3) 's% opacity(k-1)', i, k, s% opacity(k-1)
3767 : write(*,3) 's% opacity(k)', i, k, s% opacity(k)
3768 : write(*,3) 'dK_dT(I)', i, k, dK_dT(I)
3769 : write(*,3) 'dK_dT(i+1)', i, k, dK_dT(i+1)
3770 : write(*,3) 'T2', i, k, T2
3771 : write(*,3) 'T3', i, k, T3
3772 : write(*,3) 'BK', i, k, BK
3773 : write(*,3) 'BW', i, k, BW
3774 : write(*,3) '1.d0 - BK/BW', i, k, 1.d0 - BK/BW
3775 : write(*,3) 'W', i, k, W
3776 : write(*,3) 'WP', i, k, WP
3777 : call mesa_error(__FILE__,__LINE__,'T_form_of_calc_Fr')
3778 : end if
3779 : end if
3780 :
3781 0 : dFr_dVol_00 = dFr_dK_00*dK_dVol(I)
3782 : end if
3783 :
3784 : !test_partials = (k-1 == s% solver_test_partials_k)
3785 : test_partials = .false.
3786 : if (test_partials) then
3787 : s% solver_test_partials_val = Fr_00
3788 : s% solver_test_partials_var = i_var_T
3789 : s% solver_test_partials_dval_dx = dFr_dT_out
3790 : write(*,*) 'T_form_of_calc_Fr', s% solver_test_partials_var
3791 : end if
3792 : end subroutine T_form_of_calc_Fr
3793 :
3794 :
3795 0 : subroutine T_form_of_Fr_eqn(s,i)
3796 : type (star_info), pointer :: s
3797 : integer, intent(in) :: i
3798 : integer :: IL, k
3799 0 : real(dp) :: residual, Fr_00, &
3800 : dFr_00_dr_out, dFr_00_dr_in, dFr_00_dr_00, &
3801 : dFr_00_dT_out, dFr_00_dT_00, dFr_dVol_00
3802 : logical :: test_partials
3803 :
3804 : include 'formats'
3805 :
3806 0 : k = NZN+1-i
3807 0 : IL = i_var_Fr + NV*(i-1)
3808 0 : HD(1:LD_HD,IL) = 0d0
3809 :
3810 : call T_form_of_calc_Fr(s, i, Fr_00, &
3811 : dFr_00_dr_out, dFr_00_dr_in, dFr_00_dr_00, &
3812 0 : dFr_00_dT_out, dFr_00_dT_00, dFr_dVol_00)
3813 :
3814 0 : residual = Fr_00 - s% Fr(k)
3815 0 : HR(IL) = -residual
3816 :
3817 0 : HD(i_Fr_dFr_00,IL) = -1d0
3818 0 : HD(i_Fr_dr_in,IL) = dFr_00_dr_in
3819 0 : HD(i_Fr_dr_00,IL) = dFr_00_dr_00
3820 0 : HD(i_Fr_dr_out,IL) = dFr_00_dr_out
3821 0 : HD(i_Fr_dT_00,IL) = dFr_00_dT_00
3822 0 : HD(i_Fr_dT_out,IL) = dFr_00_dT_out
3823 : !test_partials = (k == s% solver_test_partials_k)
3824 0 : test_partials = .false.
3825 :
3826 : if (test_partials) then
3827 : s% solver_test_partials_val = residual
3828 : s% solver_test_partials_var = i_var_T
3829 : s% solver_test_partials_dval_dx = HD(i_Fr_dT_00,IL)
3830 : write(*,*) 'T_form_of_Fr_eqn', s% solver_test_partials_var
3831 : end if
3832 0 : end subroutine T_form_of_Fr_eqn
3833 :
3834 : end module rsp_step
|
