Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2018-2019 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module rsp_def
21 : use const_def, only: dp, qp, pi, crad, clight, ln10
22 : use math_lib
23 : use utils_lib, only: is_bad, mesa_error
24 : use star_def, only: star_info
25 : use star_utils, only: normalize_dqs, set_qs, set_m_and_dm, set_dm_bar, &
26 : store_T_in_xh, get_T_and_lnT_from_xh, store_r_in_xh, get_r_and_lnR_from_xh, &
27 : store_rho_in_xh, get_rho_and_lnd_from_xh
28 :
29 : implicit none
30 :
31 : integer, parameter :: MAX_NZN = 1501
32 :
33 : real(dp), parameter :: f_Edd_isotropic = 1d0/3d0, f_Edd_free_stream = 1d0
34 :
35 : real(dp) :: rsp_tau_factor, rsp_min_dr_div_cs, rsp_min_rad_diff_time, Psurf_from_atm
36 : integer :: i_min_dr_div_cs, i_min_rad_diff_time
37 :
38 : real(dp), pointer, dimension(:) :: &
39 : dVol_dr_00, dVol_dr_in, &
40 : d_egas_dVol, d_egas_dT, d_egas_dr_00, d_egas_dr_in, &
41 : d_Pg_dVol, d_Pg_dT, d_Pg_dr_00, d_Pg_dr_in, &
42 : dK_dVol, dK_dT, dK_dr_00, dK_dr_in, &
43 : dQQ_dVol, dQQ_dT, dQQ_dr_00, dQQ_dr_in, &
44 : dCp_dVol, dCp_dT, dCp_dr_00, dCp_dr_in, &
45 : d_Pr_dVol, d_Pr_der, d_Pr_dr_00, d_Pr_dr_in, &
46 :
47 : dHp_dr_out, dHp_dr_00, dHp_dr_in, &
48 : dHp_dVol_00, dHp_dVol_out, &
49 : dHp_dT_00, dHp_dT_out, &
50 : dHp_der_00, dHp_der_out, &
51 :
52 : dY_dr_in, dY_dr_00, dY_dr_out, &
53 : dY_dVol_00, dY_dVol_out, &
54 : dY_dT_00, dY_dT_out, &
55 : dY_der_00, dY_der_out, &
56 :
57 : dPII_dr_in, dPII_dr_00, dPII_dr_out, &
58 : dPII_dVol_00, dPII_dVol_out, &
59 : dPII_dT_00, dPII_dT_out, &
60 : dPII_der_00, dPII_der_out, &
61 :
62 : d_Pvsc_dr_00, d_Pvsc_dr_in, &
63 : d_Pvsc_dVol, d_Pvsc_dT, d_Pvsc_der, &
64 :
65 : dPtrb_dr_00, dPtrb_dr_in, dPtrb_dVol_00, dPtrb_dw_00, &
66 :
67 : dChi_dr_in2, dChi_dr_in, dChi_dr_00, dChi_dr_out, &
68 : dChi_dVol_in, dChi_dVol_00, dChi_dVol_out, &
69 : dChi_dT_in, dChi_dT_00, dChi_dT_out, &
70 : dChi_der_in, dChi_der_00, dChi_der_out, &
71 : dChi_dw_00, &
72 :
73 : dEq_dr_out, dEq_dr_00, dEq_dr_in, dEq_dr_in2, &
74 : dEq_dVol_out, dEq_dVol_00, dEq_dVol_in, &
75 : dEq_dT_out, dEq_dT_00, dEq_dT_in, &
76 : dEq_der_out, dEq_der_00, dEq_der_in, &
77 : dEq_dw_00, &
78 :
79 : dC_dr_in2, dC_dr_in, dC_dr_00, dC_dr_out, &
80 : dC_dVol_in, dC_dVol_00, dC_dVol_out, &
81 : dC_dT_in, dC_dT_00, dC_dT_out, &
82 : dC_der_in, dC_der_00, dC_der_out, &
83 : dC_dw_00, &
84 :
85 : photo_T, photo_r, photo_Vol, photo_w, photo_opacity, photo_QQ, &
86 : photo_Pgas, photo_Prad, photo_egas, photo_erad, photo_Cp, &
87 : photo_v, photo_M, photo_dm, photo_dm_bar, photo_csound
88 :
89 : ! arrays for LAPACK must be declared at compile time. legacy fortran issue.
90 : integer, parameter :: NV=5
91 : integer, parameter :: HD_DIAG=2*NV+1, LD_HD=4*NV+1, LD_ABB=6*NV+1, LPSZ=NV*MAX_NZN+1
92 : real(dp) :: DX(LPSZ), HR(LPSZ), HD(LD_HD, LPSZ), ABB(LD_ABB, LPSZ)
93 : integer :: IPVT(LPSZ)
94 :
95 : integer, parameter :: LD_LLL = 4*MAX_NZN
96 : real(dp), dimension(LD_LLL, LD_LLL) :: LLL, VLx, VRx
97 : integer :: ISORTx(LD_LLL)
98 : real(dp) :: WORKx(4*LD_LLL), WRx(LD_LLL), WIx(LD_LLL)
99 :
100 : ! for rsp_eval_eos_and_kap
101 : real(dp), pointer :: xa(:)
102 : real(dp) :: X, Z, Y, abar, zbar, z53bar, XC, XN, XO, Xne
103 :
104 : ! these for for rsp.f90 period and work calculations
105 : real(dp) :: ETOT, EGRV, ETHE, EKIN, EDE_start, ECON, &
106 : TE, ELSTA, TEFF, E0, TT1, TE_start, T0, UN, ULL, &
107 : VMAX, RMAX, LMAX, LMIN, EKMAX, EKMIN, EKMAXL, EKDEL, &
108 : RSTA, RMIN, PERIODL, PERIODLIN, &
109 : PDVWORK, FASE0
110 : real(dp), pointer, dimension(:) :: &
111 : PPP0, PPQ0, PPT0, PPC0, VV0, &
112 : WORK, WORKQ, WORKT, WORKC
113 : integer :: INSIDE, IWORK, ID, NSTART, FIRST, &
114 : run_num_retries_prev_period, prev_cycle_run_num_steps, &
115 : run_num_iters_prev_period
116 :
117 : ! for maps
118 : logical :: writing_map, done_writing_map
119 : integer, parameter :: max_map_cols = 200
120 : integer :: map_ids(max_map_cols), num_map_cols
121 : character(256) :: map_col_names(max_map_cols)
122 :
123 : ! marsaglia and zaman random number generator. period is 2**43 with
124 : ! 900 million different sequences. the state of the generator (for restarts)
125 : integer, parameter :: rn_u_len=97
126 : integer :: rn_i97, rn_j97
127 : real(dp) :: rn_u(rn_u_len), rn_c, rn_cd, rn_cm
128 :
129 : ! from const
130 : real(dp) :: G, SIG, SUNL, SUNM, SUNR, CL, P43, P4
131 :
132 : ! these are set from inlist
133 : real(dp) :: ALFA, ALFAP, ALFAM, ALFAT, ALFAS, ALFAC, CEDE, GAMMAR
134 : real(dp) :: THETA, THETAT, THETAQ, THETAU, THETAE, WTR, WTC, WTT, GAM
135 : real(dp) :: THETA1, THETAT1, THETAQ1, THETAU1, THETAE1, WTR1, WTC1, WTT1, GAM1
136 : real(dp) :: EFL0, CQ, ZSH, kapE_factor, kapP_factor
137 : integer :: NZN, IBOTOM
138 :
139 : contains
140 :
141 0 : subroutine init_def(s)
142 : use const_def, only: standard_cgrav, boltz_sigma, &
143 : Lsun, Msun, Rsun
144 : use utils_lib, only: mkdir, folder_exists
145 : type (star_info), pointer :: s
146 :
147 0 : P4=4.d0*PI
148 0 : P43=P4/3.d0
149 :
150 0 : G=standard_cgrav
151 0 : SIG=boltz_sigma
152 0 : SUNL=Lsun
153 0 : SUNM=Msun
154 0 : SUNR=Rsun
155 0 : CL=4d0*(4d0*PI)**2*SIG/3d0
156 :
157 0 : ALFA = s% RSP_alfa
158 0 : ALFAP = s% RSP_alfap
159 0 : ALFAM = s% RSP_alfam
160 0 : ALFAT = s% RSP_alfat
161 0 : ALFAS = s% RSP_alfas
162 0 : ALFAC = s% RSP_alfac
163 0 : CEDE = s% RSP_alfad
164 0 : GAMMAR = s% RSP_gammar
165 :
166 0 : ALFAP = ALFAP*2.d0/3.d0
167 0 : ALFAS = ALFAS*(1.d0/2.d0)*sqrt(2.d0/3.d0)
168 0 : ALFAC = ALFAC*(1.d0/2.d0)*sqrt(2.d0/3.d0)
169 0 : CEDE = CEDE*(8.d0/3.d0)*sqrt(2.d0/3.d0)
170 0 : GAMMAR = GAMMAR*2.d0*sqrt(3.d0)
171 :
172 0 : call turn_on_time_weighting(s)
173 :
174 0 : CQ = s% RSP_cq
175 0 : ZSH = s% RSP_zsh
176 0 : EFL0 = s% RSP_efl0
177 0 : kapE_factor = 1d0 ! s% RSP_kapE_factor
178 0 : kapP_factor = 1d0 ! s% RSP_kapP_factor
179 :
180 0 : if (ALFA == 0.d0) EFL0=0.d0
181 :
182 0 : writing_map = .false.
183 :
184 0 : if(.not. folder_exists(trim(s% log_directory))) call mkdir(trim(s% log_directory))
185 :
186 0 : end subroutine init_def
187 :
188 :
189 0 : subroutine turn_off_time_weighting(s)
190 : type (star_info), pointer :: s
191 0 : THETA = 1d0
192 0 : THETAT = 1d0
193 0 : THETAQ = 1d0
194 0 : THETAE = 1d0
195 0 : THETAU = 1d0
196 0 : WTR = 1d0
197 0 : WTC = 1d0
198 0 : WTT = 1d0
199 0 : GAM = 1d0
200 0 : GAM1=0d0
201 0 : WTR1=0d0
202 0 : WTC1=0d0
203 0 : WTT1=0d0
204 0 : THETA1 =0d0
205 0 : THETAT1=0d0
206 0 : THETAQ1=0d0
207 0 : THETAE1=0d0
208 0 : THETAU1=0d0
209 0 : end subroutine turn_off_time_weighting
210 :
211 :
212 0 : subroutine turn_on_time_weighting(s)
213 : type (star_info), pointer :: s
214 0 : THETA = s% RSP_theta
215 0 : THETAT = s% RSP_thetat
216 0 : THETAQ = s% RSP_thetaq
217 0 : THETAE = s% RSP_thetae
218 0 : THETAU = s% RSP_thetau
219 0 : WTR = s% RSP_wtr
220 0 : WTC = s% RSP_wtc
221 0 : WTT = s% RSP_wtt
222 0 : GAM = s% RSP_gam
223 0 : GAM1=1.d0-GAM
224 0 : WTR1=1.d0-WTR
225 0 : WTC1=1.d0-WTC
226 0 : WTT1=1.d0-WTT
227 0 : THETA1 =1.d0-THETA
228 0 : THETAT1=1.d0-THETAT
229 0 : THETAQ1=1.d0-THETAQ
230 0 : THETAE1=1.d0-THETAE
231 0 : THETAU1=1.d0-THETAU
232 0 : end subroutine turn_on_time_weighting
233 :
234 :
235 0 : subroutine init_allocate(s,nz)
236 : !use rsp_eddfac, only: eddfac_allocate
237 : type (star_info), pointer :: s
238 : integer, intent(in) :: nz
239 : integer :: n
240 0 : NZN = nz
241 0 : if (NZN > MAX_NZN) then
242 0 : write(*,*) 'NZN > MAX_NZN', NZN, MAX_NZN
243 0 : call mesa_error(__FILE__,__LINE__,'rsp init_allocate')
244 : end if
245 0 : IBOTOM = NZN/s% RSP_nz_div_IBOTOM
246 0 : n = NZN + 1 ! room for ghost cell
247 : allocate(xa(s% species), &
248 : dVol_dr_00(n), dVol_dr_in(n), &
249 : d_egas_dVol(n), d_egas_dT(n), d_egas_dr_00(n), d_egas_dr_in(n), &
250 : d_Pg_dVol(n), d_Pg_dT(n), d_Pg_dr_00(n), d_Pg_dr_in(n), &
251 : d_Pr_dVol(n), d_Pr_der(n), d_Pr_dr_00(n), d_Pr_dr_in(n), &
252 : dK_dVol(n), dK_dT(n), dK_dr_00(n), dK_dr_in(n), &
253 : dQQ_dVol(n), dQQ_dT(n), dQQ_dr_00(n), dQQ_dr_in(n), &
254 : dCp_dVol(n), dCp_dT(n), dCp_dr_00(n), dCp_dr_in(n), &
255 : dHp_dr_out(n), dHp_dr_00(n), dHp_dr_in(n), &
256 : dHp_dVol_00(n), dHp_dVol_out(n), &
257 : dHp_dT_00(n), dHp_dT_out(n), dHp_der_00(n), dHp_der_out(n), &
258 : dY_dr_in(n), dY_dr_00(n), dY_dr_out(n), &
259 : dY_dVol_00(n), dY_dVol_out(n), &
260 : dY_dT_00(n), dY_dT_out(n), dY_der_00(n), dY_der_out(n), &
261 : dPII_dr_in(n), dPII_dr_00(n), dPII_dr_out(n), &
262 : dPII_dVol_00(n), dPII_dVol_out(n), &
263 : dPII_dT_00(n), dPII_dT_out(n), dPII_der_00(n), dPII_der_out(n), &
264 : d_Pvsc_dr_00(n), d_Pvsc_dr_in(n), d_Pvsc_dVol(n), d_Pvsc_dT(n), d_Pvsc_der(n), &
265 : dPtrb_dr_00(n), dPtrb_dr_in(n), dPtrb_dVol_00(n), dPtrb_dw_00(n), &
266 : dChi_dr_in2(n), dChi_dr_in(n), dChi_dr_00(n), dChi_dr_out(n), &
267 : dChi_dVol_in(n), dChi_dVol_00(n), dChi_dVol_out(n), &
268 : dChi_dT_in(n), dChi_dT_00(n), dChi_dT_out(n), &
269 : dChi_der_in(n), dChi_der_00(n), dChi_der_out(n), &
270 : dChi_dw_00(n), &
271 : dEq_dr_out(n), dEq_dr_00(n), dEq_dr_in(n), dEq_dr_in2(n), &
272 : dEq_dVol_out(n), dEq_dVol_00(n), dEq_dVol_in(n), &
273 : dEq_dT_out(n), dEq_dT_00(n), dEq_dT_in(n), &
274 : dEq_der_out(n), dEq_der_00(n), dEq_der_in(n), &
275 : dEq_dw_00(n), &
276 : dC_dr_in2(n), dC_dr_in(n), dC_dr_00(n), dC_dr_out(n), &
277 : dC_dVol_in(n), dC_dVol_00(n), dC_dVol_out(n), &
278 : dC_dT_in(n), dC_dT_00(n), dC_dT_out(n), &
279 : dC_der_in(n), dC_der_00(n), dC_der_out(n), dC_dw_00(n), &
280 : photo_T(n), photo_r(n), photo_Vol(n), photo_w(n), photo_opacity(n), photo_QQ(n), &
281 : photo_Pgas(n), photo_Prad(n), photo_egas(n), photo_erad(n), photo_Cp(n), &
282 : photo_v(n), photo_M(n), photo_dm(n), photo_dm_bar(n), photo_csound(n), &
283 : PPP0(n), PPQ0(n), PPT0(n), PPC0(n), VV0(n), &
284 0 : WORK(n), WORKQ(n), WORKT(n), WORKC(n))
285 0 : end subroutine init_allocate
286 :
287 :
288 0 : subroutine init_free(s)
289 : type (star_info), pointer :: s
290 0 : deallocate(xa, &
291 0 : dVol_dr_00, dVol_dr_in, &
292 0 : d_egas_dVol, d_egas_dT, d_egas_dr_00, d_egas_dr_in, &
293 0 : d_Pg_dVol, d_Pg_dT, d_Pg_dr_00, d_Pg_dr_in, &
294 0 : d_Pr_dVol, d_Pr_der, d_Pr_dr_00, d_Pr_dr_in, &
295 0 : dK_dVol, dK_dT, dK_dr_00, dK_dr_in, &
296 0 : dQQ_dVol, dQQ_dT, dQQ_dr_00, dQQ_dr_in, &
297 0 : dCp_dVol, dCp_dT, dCp_dr_00, dCp_dr_in, &
298 0 : dHp_dr_out, dHp_dr_00, dHp_dr_in, &
299 0 : dHp_dVol_00, dHp_dVol_out, &
300 0 : dHp_dT_00, dHp_dT_out, &
301 0 : dHp_der_00, dHp_der_out, &
302 0 : dY_dr_in, dY_dr_00, dY_dr_out, &
303 0 : dY_dVol_00, dY_dVol_out, &
304 0 : dY_dT_00, dY_dT_out, dY_der_00, dY_der_out, &
305 0 : dPII_dr_in, dPII_dr_00, dPII_dr_out, &
306 0 : dPII_dVol_00, dPII_dVol_out, &
307 0 : dPII_dT_00, dPII_dT_out, dPII_der_00, dPII_der_out, &
308 0 : d_Pvsc_dr_00, d_Pvsc_dr_in, d_Pvsc_dVol, d_Pvsc_dT, d_Pvsc_der, &
309 0 : dPtrb_dr_00, dPtrb_dr_in, dPtrb_dVol_00, dPtrb_dw_00, &
310 0 : dChi_dr_in2, dChi_dr_in, dChi_dr_00, dChi_dr_out, &
311 0 : dChi_dVol_in, dChi_dVol_00, dChi_dVol_out, &
312 0 : dChi_dT_in, dChi_dT_00, dChi_dT_out, &
313 0 : dChi_der_in, dChi_der_00, dChi_der_out, &
314 0 : dChi_dw_00, &
315 0 : dEq_dr_out, dEq_dr_00, dEq_dr_in, dEq_dr_in2, &
316 0 : dEq_dVol_out, dEq_dVol_00, dEq_dVol_in, &
317 0 : dEq_dT_out, dEq_dT_00, dEq_dT_in, &
318 0 : dEq_der_out, dEq_der_00, dEq_der_in, &
319 0 : dEq_dw_00, &
320 0 : dC_dr_in2, dC_dr_in, dC_dr_00, dC_dr_out, &
321 0 : dC_dVol_in, dC_dVol_00, dC_dVol_out, &
322 0 : dC_dT_in, dC_dT_00, dC_dT_out, &
323 0 : dC_der_in, dC_der_00, dC_der_out, dC_dw_00, &
324 0 : photo_T, photo_r, photo_Vol, photo_w, photo_csound, &
325 0 : photo_Pgas, photo_Prad, photo_egas, photo_erad, photo_Cp, photo_QQ, &
326 0 : photo_v, photo_M, photo_dm, photo_dm_bar, photo_opacity, &
327 0 : PPP0, PPQ0, PPT0, PPC0, VV0, &
328 0 : WORK, WORKQ, WORKT, WORKC)
329 0 : end subroutine init_free
330 :
331 :
332 0 : real(dp) function rsp_phase_time0()
333 0 : rsp_phase_time0 = TT1
334 0 : end function rsp_phase_time0
335 :
336 :
337 0 : real(dp) function rsp_WORK(s, k)
338 : type (star_info), pointer :: s
339 : integer, intent(in) :: k
340 0 : rsp_WORK = WORK(k)
341 0 : end function rsp_WORK
342 :
343 :
344 0 : real(dp) function rsp_WORKQ(s, k)
345 : type (star_info), pointer :: s
346 : integer, intent(in) :: k
347 0 : rsp_WORKQ = WORKQ(k)
348 0 : end function rsp_WORKQ
349 :
350 :
351 0 : real(dp) function rsp_WORKT(s, k)
352 : type (star_info), pointer :: s
353 : integer, intent(in) :: k
354 0 : rsp_WORKT = WORKT(k)
355 0 : end function rsp_WORKT
356 :
357 :
358 0 : real(dp) function rsp_WORKC(s, k)
359 : type (star_info), pointer :: s
360 : integer, intent(in) :: k
361 0 : rsp_WORKC = WORKC(k)
362 0 : end function rsp_WORKC
363 :
364 :
365 0 : subroutine rsp_photo_out(s, iounit)
366 : type (star_info), pointer :: s
367 : integer, intent(in) :: iounit
368 : integer :: n
369 : include 'formats'
370 0 : n = NZN + 1
371 0 : write(iounit) NZN
372 0 : write(iounit) xa(1:s% species), &
373 0 : X, Z, Y, abar, zbar, z53bar, XC, XN, XO, Xne, IBOTOM, &
374 0 : s% RSP_num_periods, s% RSP_dt, s% RSP_period, s% rsp_DeltaR, &
375 0 : s% rsp_DeltaMag, s% rsp_GRPDV, s% rsp_GREKM, s% rsp_GREKM_avg_abs, &
376 0 : rsp_tau_factor, rsp_min_dr_div_cs, rsp_min_rad_diff_time, &
377 0 : i_min_dr_div_cs, i_min_rad_diff_time, Psurf_from_atm, &
378 0 : s% Fr(1:n), s% Lc(1:n), s% Lt(1:n), s% Y_face(1:n), &
379 0 : s% Ptrb(1:n), s% Chi(1:n), s% COUPL(1:n), s% Pvsc(1:n), &
380 0 : s% T(1:n), s% r(1:n), s% Vol(1:n), s% RSP_w(1:n), &
381 0 : s% Pgas(1:n), s% Prad(1:n), s% csound(1:n), s% Cp(1:n), &
382 0 : s% egas(1:n), s% erad(1:n), s% opacity(1:n), s% QQ(1:n), &
383 0 : s% v(1:n), s% M(1:n), s% dm(1:n), s% dm_bar(1:n), &
384 0 : ETOT, EGRV, ETHE, EKIN, EDE_start, ECON, &
385 0 : TE, ELSTA, TEFF, E0, TT1, TE_start, T0, UN, ULL, &
386 0 : VMAX, RMAX, LMAX, LMIN, EKMAX, EKMIN, EKMAXL, EKDEL, &
387 0 : RSTA, RMIN, PERIODL, PERIODLIN, &
388 0 : PDVWORK, FASE0, INSIDE, IWORK, ID, NSTART, FIRST, &
389 0 : s% rsp_LINA_periods(1:3), s% rsp_LINA_growth_rates(1:3), &
390 0 : run_num_retries_prev_period, prev_cycle_run_num_steps, &
391 0 : run_num_iters_prev_period, writing_map, &
392 0 : PPP0(1:n), PPQ0(1:n), PPT0(1:n), PPC0(1:n), VV0(1:n), &
393 0 : WORK(1:n), WORKQ(1:n), WORKT(1:n), WORKC(1:n)
394 0 : end subroutine rsp_photo_out
395 :
396 :
397 0 : subroutine rsp_photo_in(s, iounit, ierr)
398 : type (star_info), pointer :: s
399 : integer, intent(in) :: iounit
400 : integer, intent(out) :: ierr
401 : integer :: n
402 : include 'formats'
403 0 : call init_def(s)
404 0 : ierr = 0
405 0 : read(iounit, iostat=ierr) NZN
406 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'read failed in rsp_photo_in')
407 0 : s% nz = NZN
408 0 : call init_allocate(s,NZN)
409 0 : n = NZN + 1
410 0 : read(iounit, iostat=ierr) xa(1:s% species), &
411 0 : X, Z, Y, abar, zbar, z53bar, XC, XN, XO, Xne, IBOTOM, &
412 0 : s% RSP_num_periods, s% RSP_dt, s% RSP_period, s% rsp_DeltaR, &
413 0 : s% rsp_DeltaMag, s% rsp_GRPDV, s% rsp_GREKM, s% rsp_GREKM_avg_abs, &
414 0 : rsp_tau_factor, rsp_min_dr_div_cs, rsp_min_rad_diff_time, &
415 0 : i_min_dr_div_cs, i_min_rad_diff_time, Psurf_from_atm, &
416 0 : s% Fr(1:n), s% Lc(1:n), s% Lt(1:n), s% Y_face(1:n), &
417 0 : s% Ptrb(1:n), s% Chi(1:n), s% COUPL(1:n), s% Pvsc(1:n), &
418 0 : photo_T(1:n), photo_r(1:n), photo_Vol(1:n), photo_w(1:n), &
419 0 : photo_Pgas(1:n), photo_Prad(1:n), photo_csound(1:n), photo_Cp(1:n), &
420 0 : photo_egas(1:n), photo_erad(1:n), photo_opacity(1:n), photo_QQ(1:n), &
421 0 : photo_v(1:n), photo_M(1:n), photo_dm(1:n), photo_dm_bar(1:n), &
422 0 : ETOT, EGRV, ETHE, EKIN, EDE_start, ECON, &
423 0 : TE, ELSTA, TEFF, E0, TT1, TE_start, T0, UN, ULL, &
424 0 : VMAX, RMAX, LMAX, LMIN, EKMAX, EKMIN, EKMAXL, EKDEL, &
425 0 : RSTA, RMIN, PERIODL, PERIODLIN, &
426 0 : PDVWORK, FASE0, INSIDE, IWORK, ID, NSTART, FIRST, &
427 0 : s% rsp_LINA_periods(1:3), s% rsp_LINA_growth_rates(1:3), &
428 0 : run_num_retries_prev_period, prev_cycle_run_num_steps, &
429 0 : run_num_iters_prev_period, writing_map, &
430 0 : PPP0(1:n), PPQ0(1:n), PPT0(1:n), PPC0(1:n), VV0(1:n), &
431 0 : WORK(1:n), WORKQ(1:n), WORKT(1:n), WORKC(1:n)
432 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'read failed in rsp_photo_in')
433 0 : if (writing_map) then
434 0 : write(*,*) 'sorry, cannot use photo to resume writing map data file'
435 0 : writing_map = .false.
436 : end if
437 0 : end subroutine rsp_photo_in
438 :
439 :
440 0 : subroutine finish_after_build_model(s)
441 : type (star_info), pointer :: s
442 : integer :: k
443 : ! restore bit-for-bit erad = crad*T**4*Vol
444 : include 'formats'
445 0 : do k = 1, NZN
446 0 : s% erad(k) = crad*s% T(k)**4*s% Vol(k)
447 0 : s% Prad(k) = s% f_Edd(k)*s% erad(k)/s% Vol(k)
448 : end do
449 0 : end subroutine finish_after_build_model
450 :
451 :
452 0 : subroutine finish_rsp_photo_in(s)
453 : use star_utils, only: set_rmid
454 : type (star_info), pointer :: s
455 : integer :: k, ierr
456 : ! restore bit-for-bit same as before photo
457 : include 'formats'
458 0 : do k = 1, NZN
459 0 : s% T(k) = photo_T(k)
460 0 : s% Pgas(k) = photo_Pgas(k)
461 0 : s% Prad(k) = photo_Prad(k)
462 0 : s% Peos(k) = s% Pgas(k) + s% Prad(k)
463 0 : s% egas(k) = photo_egas(k)
464 0 : s% erad(k) = photo_erad(k)
465 0 : s% opacity(k) = photo_opacity(k)
466 0 : s% csound(k) = photo_csound(k)
467 0 : s% Cp(k) = photo_Cp(k)
468 0 : s% QQ(k) = photo_QQ(k)
469 0 : s% r(k) = photo_r(k)
470 0 : s% Vol(k) = photo_Vol(k)
471 0 : s% RSP_w(k) = photo_w(k)
472 0 : s% v(k) = photo_v(k)
473 0 : s% M(k) = photo_M(k)
474 0 : s% dm(k) = photo_dm(k)
475 0 : s% dm_bar(k) = photo_dm_bar(k)
476 : end do
477 0 : call set_rmid(s, 1, NZN, ierr)
478 0 : if (ierr /= 0) then
479 0 : write(*,*) 'finish_rsp_photo_in failed in set_rmid'
480 0 : call mesa_error(__FILE__,__LINE__,'finish_rsp_photo_in')
481 : end if
482 0 : end subroutine finish_rsp_photo_in
483 :
484 :
485 0 : subroutine set_build_vars(s, &
486 0 : m, dm, dm_bar, r, Vol, T, RSP_Et, Lr, Lc)
487 : type (star_info), pointer :: s
488 : real(dp), dimension(:), intent(in) :: &
489 : m, dm, dm_bar, r, Vol, T, RSP_Et, Lr, Lc
490 : integer :: k, i
491 : include 'formats'
492 0 : do i=1, NZN
493 0 : k = NZN+1 - i
494 0 : s% m(k) = m(i)
495 0 : s% dm(k) = dm(i)
496 0 : s% dm_bar(k) = dm_bar(i)
497 0 : s% r(k) = r(i)
498 0 : s% Vol(k) = Vol(i)
499 0 : s% T(k) = T(i)
500 0 : s% RSP_Et(k) = RSP_Et(i)
501 0 : s% RSP_w(k) = sqrt(s% RSP_Et(k))
502 0 : s% Fr(k) = Lr(i)/(4d0*pi*s% r(k)**2)
503 0 : s% erad(k) = crad*s% T(k)**4*s% Vol(k)
504 0 : s% L(k) = Lr(i) + Lc(i)
505 0 : if (is_bad(s% L(k))) then
506 0 : write(*,2) 'L Lr Lc', k, s% L(k), Lr(i), Lc(i)
507 0 : call mesa_error(__FILE__,__LINE__,'set_build_vars')
508 : end if
509 : end do
510 0 : end subroutine set_build_vars
511 :
512 :
513 0 : subroutine set_star_vars(s, ierr)
514 : type (star_info), pointer :: s
515 : integer, intent(out) :: ierr
516 0 : real(dp) :: sum_dm
517 : integer :: k
518 : include 'formats'
519 0 : ierr = 0
520 0 : sum_dm = 0d0
521 0 : do k=1, NZN
522 0 : sum_dm = sum_dm + s% dm(k)
523 : end do
524 0 : do k=1, NZN
525 0 : s% dq(k) = s% dm(k)/sum_dm
526 : end do
527 0 : s% xmstar = sum_dm
528 0 : if (s% nz /= NZN) then
529 0 : write(*,2) 'NZN', NZN
530 0 : write(*,2) 's% nz', s% nz
531 0 : call mesa_error(__FILE__,__LINE__,'bad nz')
532 : end if
533 0 : if (.not. s% do_normalize_dqs_as_part_of_set_qs) then
534 0 : call normalize_dqs(s, s% nz, s% dq, ierr)
535 0 : if (ierr /= 0) then
536 0 : if (s% report_ierr) write(*,*) 'normalize_dqs failed in rsp_def set_star_vars'
537 0 : return
538 : end if
539 : end if
540 0 : call set_qs(s, s% nz, s% q, s% dq, ierr)
541 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'failed in set_qs')
542 0 : s% m(1) = s% mstar
543 0 : do k=2,s% nz
544 0 : s% dm(k-1) = s% dq(k-1)*s% xmstar
545 0 : s% m(k) = s% m(k-1) - s% dm(k-1)
546 : end do
547 : k = s% nz
548 0 : s% dm(k) = s% m(k) - s% m_center
549 0 : call set_dm_bar(s, s% nz, s% dm, s% dm_bar)
550 0 : do k=1, NZN
551 :
552 0 : if (k==NZN) then
553 0 : s% Vol(k)=P43/s% dm(k)*(s% r(k)**3 - s% R_center**3)
554 0 : s% rmid(k) = 0.5d0*(s% r(k) + s% R_center)
555 : else
556 0 : s% Vol(k)=P43/s% dm(k)*(s% r(k)**3 - s% r(k+1)**3)
557 0 : s% rmid(k) = 0.5d0*(s% r(k) + s% r(k+1))
558 : end if
559 0 : if (is_bad(s% Vol(k)))then
560 0 : write(*, 2) 's% Vol(k)', k, s% Vol(k)
561 0 : call mesa_error(__FILE__,__LINE__,'set_star_vars')
562 : end if
563 :
564 0 : s% rho(k) = 1d0/s% Vol(k)
565 0 : call store_rho_in_xh(s, k, s% rho(k))
566 0 : call get_rho_and_lnd_from_xh(s, k, s% rho(k), s% lnd(k))
567 0 : s% Vol(k) = 1d0/s% rho(k)
568 :
569 0 : call store_T_in_xh(s, k, s% T(k))
570 0 : call get_T_and_lnT_from_xh(s, k, s% T(k), s% lnT(k))
571 :
572 0 : call store_r_in_xh(s, k, s% r(k))
573 0 : call get_r_and_lnR_from_xh(s, k, s% r(k), s% lnR(k))
574 :
575 0 : s% RSP_Et(k) = s% RSP_w(k)*s% RSP_w(k)
576 0 : s% xh(s% i_Et_RSP, k) = s% RSP_Et(k)
577 0 : s% xh(s% i_v, k) = s% v(k)
578 : end do
579 : end subroutine set_star_vars
580 :
581 :
582 0 : subroutine copy_from_xh_to_rsp(s, nz_new) ! do this when load a file
583 : use star_utils, only: get_T_and_lnT_from_xh, get_r_and_lnR_from_xh
584 : type (star_info), pointer :: s
585 : integer, intent(in) :: nz_new
586 : integer :: k
587 0 : real(qp) :: q1, q2, q3, q4
588 : include 'formats'
589 0 : if (nz_new > 0) NZN = nz_new
590 0 : do k=NZN,1,-1
591 0 : call get_rho_and_lnd_from_xh(s, k, s% rho(k), s% lnd(k))
592 0 : call get_T_and_lnT_from_xh(s, k, s% T(k), s% lnT(k))
593 0 : call get_r_and_lnR_from_xh(s, k, s% r(k), s% lnR(k))
594 0 : s% RSP_Et(k) = s% xh(s% i_Et_RSP,k)
595 0 : s% RSP_w(k) = sqrt(s% RSP_Et(k))
596 0 : s% Fr(k) = s% xh(s% i_Fr_RSP,k)
597 0 : s% v(k) = s% xh(s% i_v,k)
598 0 : if (k == NZN) then ! center
599 0 : s% Vol(k)=P43/s% dm(k)*(s% r(k)**3 - s% R_center**3)
600 0 : s% rmid(k) = 0.5d0*(s% r(k) + s% R_center)
601 0 : if (s% Vol(k) <= 0d0 .or. is_bad(s% Vol(k))) then
602 0 : write(*,2) 'copy from xh to rsp s% Vol(k) r00 r_center dm', k, &
603 0 : s% Vol(k), s% r(k), s% R_center, s% dm(k)
604 0 : call mesa_error(__FILE__,__LINE__,'copy_from_xh_to_rsp')
605 : end if
606 : else
607 0 : q1 = P43/s% dm(k)
608 0 : q2 = s% r(k)
609 0 : q3 = s% r(k+1)
610 0 : q4 = q1*(q2**3 - q3**3)
611 0 : s% Vol(k) = dble(q4)
612 0 : s% rmid(k) = 0.5d0*(s% r(k) + s% r(k+1))
613 0 : if (s% Vol(k) <= 0d0 .or. is_bad(s% Vol(k))) then
614 0 : write(*,2) 'copy from xh to rsp s% Vol(k) r00 rp1 dm', k, &
615 0 : s% Vol(k), s% r(k), s% r(k+1), s% dm(k)
616 0 : call mesa_error(__FILE__,__LINE__,'copy_from_xh_to_rsp')
617 : end if
618 : end if
619 0 : s% erad(k) = s% xh(s% i_erad_RSP,k)
620 0 : s% Prad(k) = f_Edd_isotropic*s% erad(k)/s% Vol(k)
621 : end do
622 0 : end subroutine copy_from_xh_to_rsp
623 :
624 :
625 0 : subroutine check_for_T_or_P_inversions(s,str)
626 : type (star_info), pointer :: s
627 : character (len=*), intent(in) :: str
628 : integer :: k
629 : logical :: okay
630 : include 'formats'
631 0 : okay = .true.
632 0 : do k=2, s% nz
633 0 : if (s% T(k) <= s% T(k-1)) then
634 0 : write(*,3) trim(str) // ' T inversion', k, s% model_number, s% T(k), s% T(k-1)
635 0 : okay = .false.
636 : end if
637 : end do
638 0 : do k=2, s% nz
639 0 : if (s% Peos(k) <= s% Peos(k-1)) then
640 0 : write(*,3) trim(str) // ' Peos inversion', k, s% model_number, s% Peos(k), s% Peos(k-1), &
641 0 : s% Ptrb(k), s% Ptrb(k-1), s% Pvsc(k), s% Pvsc(k-1), s% v(k+1), s% v(k), s% v(k-1)
642 0 : okay = .false.
643 : end if
644 : end do
645 0 : if (.not. okay) call mesa_error(__FILE__,__LINE__,'rsp_one_step check_for_T_or_P_inversions')
646 0 : end subroutine check_for_T_or_P_inversions
647 :
648 :
649 0 : subroutine check_R(s,str)
650 : type (star_info), pointer :: s
651 : character (len=*), intent(in) :: str
652 : integer :: k
653 : include 'formats'
654 0 : do k=s% nz,1,-1
655 0 : if (k == s% nz) then
656 0 : if (s% r(k) <= s% r_center) then
657 0 : write(*,3) trim(str) // ' bad r', k, s% model_number, s% r(k), s% r_center
658 0 : call mesa_error(__FILE__,__LINE__,'rsp_one_step check_R')
659 : end if
660 : else
661 0 : if (s% r(k) <= s% r(k+1)) then
662 0 : write(*,3) trim(str) // ' bad r', k, s% model_number, s% r(k), s% r(k+1)
663 0 : call mesa_error(__FILE__,__LINE__,'rsp_one_step check_R')
664 : end if
665 : end if
666 : end do
667 0 : end subroutine check_R
668 :
669 :
670 0 : subroutine rsp_dump_for_debug(s)
671 : type (star_info), pointer :: s
672 : integer :: k
673 : include 'formats'
674 0 : write(*,*) 'rsp_dump_for_debug'
675 0 : write(*,2) 'R_center', s% model_number, s% R_center
676 0 : write(*,2) 'xmstar', s% model_number, s% xmstar
677 0 : write(*,2) 'M/Msun', s% model_number, s% star_mass
678 0 : write(*,2) 's% R_center', s% model_number, s% R_center
679 0 : write(*,4) 'species', s% model_number, s% species
680 0 : write(*,2) 'm_center', s% model_number, s% M_center
681 0 : write(*,2) 'mstar', s% model_number, s% mstar
682 0 : write(*,2) 'L_center', s% model_number, s% L_center
683 0 : write(*,2) 'X', s% model_number, X
684 0 : write(*,2) 'Z', s% model_number, Z
685 0 : write(*,2) 'Y', s% model_number, Y
686 0 : write(*,2) 'abar', s% model_number, abar
687 0 : write(*,2) 'zbar', s% model_number, zbar
688 0 : write(*,2) 'z53bar', s% model_number, z53bar
689 0 : write(*,2) 'XC', s% model_number, XC
690 0 : write(*,2) 'XN', s% model_number, XN
691 0 : write(*,2) 'XO', s% model_number, XO
692 0 : write(*,2) 'Xne', s% model_number, Xne
693 0 : write(*,2) 'dt dt_next', s% model_number, s% dt
694 0 : write(*,2) 's% rsp_period', s% model_number, s% rsp_period
695 0 : write(*,2) 'dt', s% model_number, s% dt
696 0 : do k=1,s% nz
697 0 : write(*,2) '% Vol(k)', k, s% Vol(k)
698 0 : write(*,2) 's% v(k)', k, s% v(k)
699 0 : write(*,2) 's% r(k)', k, s% r(k)
700 0 : write(*,2) 's% dm(k)', k, s% dm(k)
701 0 : write(*,2) 's% RSP_w(k)', k, s% RSP_w(k)
702 0 : write(*,2) 's% T(k)', k, s% T(k)
703 0 : write(*,2) 's% erad(k)', k, s% erad(k)
704 0 : write(*,2) 's% Prad(k)', k, s% Prad(k)
705 0 : write(*,2) 's% Fr(k)', k, s% Fr(k)
706 : !write(*,2) '', k,
707 : end do
708 : !call mesa_error(__FILE__,__LINE__,'rsp_dump_for_debug')
709 0 : end subroutine rsp_dump_for_debug
710 :
711 :
712 0 : subroutine cleanup_for_LINA( &
713 0 : s, M, DM, DM_BAR, R, Vol, T, RSP_Et, Peos, ierr)
714 : use star_utils, only: normalize_dqs, set_qs, set_m_and_dm, set_dm_bar
715 : type (star_info), pointer :: s
716 : real(dp), intent(inout), dimension(:) :: &
717 : M, DM, DM_BAR, R, Vol, T, RSP_Et, Peos
718 : integer, intent(out) :: ierr
719 :
720 : integer :: I, k
721 :
722 : include 'formats'
723 :
724 : ! get
725 0 : do i=1,NZN
726 0 : k = NZN+1-i
727 0 : s% m(k) = M(i)
728 0 : s% dq(k) = DM(i)/s% xmstar
729 0 : s% r(k) = R(i)
730 0 : s% Vol(k) = Vol(i)
731 0 : s% T(k) = T(i)
732 0 : s% RSP_w(k) = sqrt(RSP_Et(i))
733 0 : s% Peos(k) = Peos(i)
734 0 : s% Prad(k) = crad*s% T(k)**4/3d0
735 0 : s% Pgas(k) = s% Peos(k) - s% Prad(k)
736 : end do
737 0 : s% dq(s% nz) = (s% m(NZN) - s% M_center)/s% xmstar
738 :
739 : ! fix
740 0 : if (.not. s% do_normalize_dqs_as_part_of_set_qs) then
741 0 : call normalize_dqs(s, NZN, s% dq, ierr)
742 0 : if (ierr /= 0) then
743 0 : if (s% report_ierr) write(*,*) 'normalize_dqs failed in rsp_def cleanup_for_LINA'
744 : return
745 : end if
746 : end if
747 0 : call set_qs(s, NZN, s% q, s% dq, ierr)
748 0 : if (ierr /= 0) then
749 0 : write(*,*) 'failed in set_qs'
750 0 : call mesa_error(__FILE__,__LINE__,'build_rsp_model')
751 : end if
752 0 : call set_m_and_dm(s)
753 0 : call set_dm_bar(s, s% nz, s% dm, s% dm_bar)
754 0 : s% Y_face(1:s% nz) = 0d0
755 0 : s% grada(1:s% nz) = 0d0
756 0 : s% rmid_start(1:s% nz) = -1d0
757 0 : s% Fr(1:s% nz) = 0d0
758 0 : s% Lc(1:s% nz) = 0d0
759 0 : s% Lt(1:s% nz) = 0d0
760 0 : s% csound(1:s% nz) = 0d0
761 0 : call copy_results(s)
762 :
763 : ! put back
764 0 : do i=1,NZN
765 0 : k = NZN+1-i
766 0 : M(i) = s% m(k)
767 0 : DM(i) = s% dm(k)
768 0 : DM_BAR(i) = s% dm_bar(k)
769 0 : R(i) = s% r(k)
770 0 : Vol(i) = s% Vol(k)
771 0 : T(i) = s% T(k)
772 0 : RSP_Et(i) = s% RSP_w(k)**2
773 : end do
774 :
775 0 : end subroutine cleanup_for_LINA
776 :
777 :
778 0 : subroutine copy_results(s)
779 : use star_utils, only: set_rmid, store_r_in_xh, &
780 : get_r_and_lnR_from_xh, store_r_in_xh, &
781 0 : get_T_and_lnT_from_xh, store_T_in_xh
782 : use const_def, only: convective_mixing, no_mixing, qp
783 : type (star_info), pointer :: s
784 : integer :: i, k, ierr
785 : real(dp) :: RSP_efl0_2
786 0 : real(qp) :: q1, q2, q3, q4
787 :
788 0 : RSP_efl0_2 = EFL0**2
789 0 : do i=1, NZN
790 :
791 0 : k = NZN+1 - i
792 0 : s% xh(s% i_v,k) = s% v(k)
793 0 : s% xh(s% i_erad_RSP,k) = s% erad(k)
794 0 : s% xh(s% i_Fr_RSP,k) = s% Fr(k)
795 :
796 : ! sqrt(w**2) /= original w, so need to redo
797 0 : s% RSP_Et(k) = s% RSP_w(k)**2
798 0 : s% xh(s% i_Et_RSP,k) = s% RSP_Et(k)
799 0 : s% RSP_w(k) = sqrt(s% xh(s% i_Et_RSP,k))
800 :
801 : ! exp(log(r)) /= original r, so need to redo
802 0 : call store_r_in_xh(s, k, s% r(k))
803 0 : call get_r_and_lnR_from_xh(s, k, s% r(k), s% lnR(k))
804 :
805 : ! exp(log(T)) /= original T, so need to redo
806 0 : call store_T_in_xh(s, k, s% T(k))
807 0 : call get_T_and_lnT_from_xh(s, k, s% T(k), s% lnT(k))
808 :
809 0 : s% Peos(k) = s% Pgas(k) + s% Prad(k)
810 0 : if (k > 1) s% gradT(k) = &
811 0 : s% Y_face(k) + 0.5d0*(s% grada(k-1) + s% grada(k))
812 :
813 : end do
814 0 : s% gradT(1) = s% gradT(2)
815 :
816 0 : call set_rmid(s, 1, NZN, ierr)
817 0 : if (ierr /= 0) then
818 0 : write(*,*) 'copy_results failed in set_rmid'
819 0 : call mesa_error(__FILE__,__LINE__,'copy_results')
820 : end if
821 :
822 0 : do i=1, NZN
823 0 : k = NZN+1 - i
824 : ! revise Vol and rho using revised r
825 0 : if (i==1) then
826 0 : s% Vol(k)=P43/s% dm(k)*(s% r(k)**3 - s% R_center**3)
827 : else
828 0 : q1 = P43/s% dm(k)
829 0 : q2 = s% r(k)
830 0 : q3 = s% r(k+1)
831 0 : q4 = q1*(q2**3 - q3**3)
832 0 : s% Vol(k) = dble(q4)
833 : end if
834 0 : s% rho(k) = 1d0/s% Vol(k)
835 0 : call store_rho_in_xh(s, k, s% rho(k))
836 0 : call get_rho_and_lnd_from_xh(s, k, s% rho(k), s% lnd(k))
837 0 : s% Vol(k) = 1d0/s% rho(k)
838 0 : s% L(k) = 4d0*pi*s% r(k)**2*s% Fr(k) + s% Lc(k) + s% Lt(k)
839 0 : if (s% RSP_w(k) > 1d4) then ! arbitrary cut
840 0 : s% mixing_type(k) = convective_mixing
841 : else
842 0 : s% mixing_type(k) = no_mixing
843 : end if
844 : end do
845 :
846 : ! set some things for mesa output reporting
847 0 : i = 1
848 0 : s% rho_face(i) = s% rho(i)
849 0 : s% P_face_ad(i)%val = s% Peos(i)
850 0 : s% csound_face(i) = s% csound(i)
851 0 : do i = 2,NZN
852 0 : s% rho_face(i) = 0.5d0*(s% rho(i) + s% rho(i-1))
853 0 : s% P_face_ad(i)%val = 0.5d0*(s% Peos(i) + s% Peos(i-1))
854 0 : s% csound_face(i) = 0.5d0*(s% csound(i) + s% csound(i-1))
855 : end do
856 :
857 : ! these are necessary to make files consistent with photos.
858 0 : s% R_center = pow(s% r(NZN)**3 - s% Vol(NZN)*s% dm(NZN)/P43, 1d0/3d0)
859 0 : s% M_center = s% mstar - s% xmstar
860 :
861 0 : end subroutine copy_results
862 :
863 : end module rsp_def
|
