Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2022 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 pc_support
21 : use utils_lib, only: is_bad, mesa_error, mv, switch_str
22 : use eos_def
23 : use const_def, only: ln10, mesa_temp_caches_dir, one_sixth
24 : use math_lib
25 :
26 : implicit none
27 :
28 : public :: get_FITION9, get_EXCOR7, get_FSCRliq8
29 : private
30 :
31 : ! the file data for EXCOR7
32 : integer, parameter :: jFXC = 1
33 : integer, parameter :: jUXC = 2
34 : integer, parameter :: jPXC = 3
35 : integer, parameter :: jCVXC = 4
36 : integer, parameter :: jSXC = 5
37 : integer, parameter :: jPDTXC = 6
38 : integer, parameter :: jPDRXC = 7
39 : integer, parameter :: nvals_EXCOR7 = 7
40 :
41 : ! the file data for FSCRliq8 FSCR, USCR, PSCR, CVSCR, PDTSCR, PDRSCR
42 : integer, parameter :: iFSCR = 1
43 : integer, parameter :: iUSCR = 2
44 : integer, parameter :: iPSCR = 3
45 : integer, parameter :: iCVSCR = 4
46 : integer, parameter :: iPDTSCR = 5
47 : integer, parameter :: iPDRSCR = 6
48 : integer, parameter :: nvals_FSCRliq8 = 6
49 :
50 :
51 : contains
52 :
53 :
54 0 : subroutine get_FITION9(GAMI, &
55 : FION, dFION_dlnGAMI, &
56 : UION, dUION_dlnGAMI, &
57 : PION, dPION_dlnGAMI, &
58 : CVii, dCVii_dlnGAMI, &
59 : PDTii, dPDTii_dlnGAMI, &
60 : PDRii, dPDRii_dlnGAMI, &
61 : skip, ierr)
62 : real(dp), intent(in) :: GAMI
63 : real(dp), intent(out) :: &
64 : FION, dFION_dlnGAMI, &
65 : UION, dUION_dlnGAMI, &
66 : PION, dPION_dlnGAMI, &
67 : CVii, dCVii_dlnGAMI, &
68 : PDTii, dPDTii_dlnGAMI, &
69 : PDRii, dPDRii_dlnGAMI
70 : logical, intent(out) :: skip
71 : integer, intent(out) :: ierr
72 :
73 : integer, parameter :: iFION = 1
74 : integer, parameter :: iUION = 2
75 : integer, parameter :: iPION = 3
76 : integer, parameter :: iCVii = 4
77 : integer, parameter :: iPDTii = 5
78 : integer, parameter :: iPDRii = 6
79 0 : real(dp) :: lnGAMI_min, lnGAMI_max, dlnGAMI, lnGAMI
80 : integer :: nlnGAMI, k_old
81 0 : real(dp) :: xk_old, xkp1_old, xk_new, delta
82 0 : real(dp), pointer, dimension(:,:) :: f_1, f_2, f_3, f_4, f_5, f_6
83 :
84 : include 'formats'
85 :
86 0 : ierr = 0
87 0 : skip = .false.
88 :
89 0 : lnGAMI_min = -2.0d0*ln10
90 0 : lnGAMI_max = 6.5d0*ln10
91 0 : dlnGAMI = 1d-2*ln10
92 0 : nlnGAMI = (lnGAMI_max - lnGAMI_min)/dlnGAMI + 1
93 :
94 0 : if (.not. FITION9_loaded) then
95 0 : !$OMP CRITICAL(eosPC_support_FITION9_load)
96 0 : if (.not. FITION9_loaded) then
97 0 : call build_FITION9_data(ierr)
98 0 : if (ierr == 0) FITION9_loaded = .true.
99 : end if
100 : !$OMP END CRITICAL(eosPC_support_FITION9_load)
101 0 : if (ierr /= 0) return
102 : end if
103 :
104 0 : lnGAMI = log(GAMI)
105 0 : if (lnGAMI < lnGAMI_min .or. lnGAMI > lnGAMI_max) then
106 0 : skip = .true.
107 0 : return
108 : end if
109 :
110 0 : f_1(1:4,1:nlnGAMI) => FITION_data(iFION)% f1
111 0 : f_2(1:4,1:nlnGAMI) => FITION_data(iUION)% f1
112 0 : f_3(1:4,1:nlnGAMI) => FITION_data(iPION)% f1
113 0 : f_4(1:4,1:nlnGAMI) => FITION_data(iCVii)% f1
114 0 : f_5(1:4,1:nlnGAMI) => FITION_data(iPDTii)% f1
115 0 : f_6(1:4,1:nlnGAMI) => FITION_data(iPDRii)% f1
116 :
117 0 : xk_new = lnGAMI
118 0 : k_old = int((lnGAMI - lnGAMI_min)/dlnGAMI + 1d-4) + 1
119 0 : if (k_old < 1 .or. k_old >= nlnGAMI) then
120 : if (k_old < 1) then
121 0 : k_old = 1
122 0 : xk_old = lnGAMI_min
123 0 : xkp1_old = xk_old + dlnGAMI
124 : else
125 0 : k_old = nlnGAMI-1
126 0 : xk_old = lnGAMI_min + (k_old-1)*dlnGAMI
127 0 : xkp1_old = xk_old + dlnGAMI
128 : end if
129 : else
130 0 : xk_old = lnGAMI_min + (k_old-1)*dlnGAMI
131 0 : xkp1_old = xk_old + dlnGAMI
132 : end if
133 :
134 0 : delta = xk_new - xk_old
135 :
136 : FION = &
137 : f_1(1, k_old) + delta*(f_1(2, k_old) &
138 0 : + delta*(f_1(3, k_old) + delta*f_1(4, k_old)))
139 : dFION_dlnGAMI = &
140 0 : f_1(2, k_old) + 2*delta*(f_1(3, k_old) + 1.5d0*delta*f_1(4, k_old))
141 :
142 : UION = &
143 : f_2(1, k_old) + delta*(f_2(2, k_old) &
144 0 : + delta*(f_2(3, k_old) + delta*f_2(4, k_old)))
145 : dUION_dlnGAMI = &
146 0 : f_2(2, k_old) + 2*delta*(f_2(3, k_old) + 1.5d0*delta*f_2(4, k_old))
147 :
148 : PION = &
149 : f_3(1, k_old) + delta*(f_3(2, k_old) &
150 0 : + delta*(f_3(3, k_old) + delta*f_3(4, k_old)))
151 : dPION_dlnGAMI = &
152 0 : f_3(2, k_old) + 2*delta*(f_3(3, k_old) + 1.5d0*delta*f_3(4, k_old))
153 :
154 : CVii = &
155 : f_4(1, k_old) + delta*(f_4(2, k_old) &
156 0 : + delta*(f_4(3, k_old) + delta*f_4(4, k_old)))
157 : dCVii_dlnGAMI = &
158 0 : f_4(2, k_old) + 2*delta*(f_4(3, k_old) + 1.5d0*delta*f_4(4, k_old))
159 :
160 : PDTii = &
161 : f_5(1, k_old) + delta*(f_5(2, k_old) &
162 0 : + delta*(f_5(3, k_old) + delta*f_5(4, k_old)))
163 : dPDTii_dlnGAMI = &
164 0 : f_5(2, k_old) + 2*delta*(f_5(3, k_old) + 1.5d0*delta*f_5(4, k_old))
165 :
166 : PDRii = &
167 : f_6(1, k_old) + delta*(f_6(2, k_old) &
168 0 : + delta*(f_6(3, k_old) + delta*f_6(4, k_old)))
169 : dPDRii_dlnGAMI = &
170 0 : f_6(2, k_old) + 2*delta*(f_6(3, k_old) + 1.5d0*delta*f_6(4, k_old))
171 :
172 : contains
173 :
174 :
175 0 : subroutine build_FITION9_data(ierr)
176 : use interp_1d_def, only : pm_work_size
177 : use interp_1d_lib, only : interp_pm
178 : integer, intent(out) :: ierr
179 : type (FITION_Info), pointer :: fi
180 0 : real(dp), pointer :: work(:)
181 0 : real(dp) :: lnGAMI, GAMI, FION, UION, PION, CVii, PDTii, PDRii
182 : integer :: j
183 : include 'formats'
184 :
185 0 : ierr = 0
186 : !write(*,'(a)') 'eosPC build_FITION9_data'
187 :
188 0 : allocate(FITION9_lnGAMIs(nlnGAMI), work(pm_work_size*nlnGAMI))
189 0 : do j=1,FITION_vals
190 0 : fi => FITION_data(j)
191 0 : allocate(fi% f1(4*nlnGAMI))
192 0 : fi% f(1:4, 1:nlnGAMI) => fi% f1(1:4*nlnGAMI)
193 : end do
194 :
195 0 : do j=1,nlnGAMI-1
196 0 : FITION9_lnGAMIs(j) = lnGAMI_min + (j-1)*dlnGAMI
197 : end do
198 0 : FITION9_lnGAMIs(nlnGAMI) = lnGAMI_max
199 :
200 0 : do j=1,nlnGAMI
201 0 : lnGAMI = FITION9_lnGAMIs(j)
202 0 : GAMI = exp(lnGAMI)
203 0 : call FITION9(GAMI,FION,UION,PION,CVii,PDTii,PDRii)
204 0 : fi => FITION_data(iFION); fi% f(1,j) = FION
205 0 : fi => FITION_data(iUION); fi% f(1,j) = UION
206 0 : fi => FITION_data(iPION); fi% f(1,j) = PION
207 0 : fi => FITION_data(iCVii); fi% f(1,j) = CVii
208 0 : fi => FITION_data(iPDTii); fi% f(1,j) = PDTii
209 0 : fi => FITION_data(iPDRii); fi% f(1,j) = PDRii
210 : end do
211 :
212 0 : do j=1,FITION_vals
213 0 : fi => FITION_data(j)
214 : call interp_pm(FITION9_lnGAMIs, nlnGAMI, &
215 0 : fi% f1, pm_work_size, work, 'build_FITION9_data', ierr)
216 0 : if (ierr /= 0) return
217 : end do
218 :
219 0 : deallocate(work)
220 :
221 : !write(*,'(a)') 'done eosPC build_FITION9_data'
222 :
223 0 : end subroutine build_FITION9_data
224 :
225 : end subroutine get_FITION9
226 :
227 :
228 0 : subroutine get_FSCRliq8(iZion, RS, GAME, &
229 : FSCR, dFSCR_dlnRS, dFSCR_dlnGAME, &
230 : USCR, dUSCR_dlnRS, dUSCR_dlnGAME, &
231 : PSCR, dPSCR_dlnRS, dPSCR_dlnGAME, &
232 : CVSCR, dCVSCR_dlnRS, dCVSCR_dlnGAME, &
233 : PDTSCR, dPDTSCR_dlnRS, dPDTSCR_dlnGAME, &
234 : PDRSCR, dPDRSCR_dlnRS, dPDRSCR_dlnGAME, &
235 : skip, ierr)
236 : integer, intent(in) :: iZion
237 : real(dp), intent(in) :: RS, GAME
238 : real(dp), intent(out) :: &
239 : FSCR, dFSCR_dlnRS, dFSCR_dlnGAME, &
240 : USCR, dUSCR_dlnRS, dUSCR_dlnGAME, &
241 : PSCR, dPSCR_dlnRS, dPSCR_dlnGAME, &
242 : CVSCR, dCVSCR_dlnRS, dCVSCR_dlnGAME, &
243 : PDTSCR, dPDTSCR_dlnRS, dPDTSCR_dlnGAME, &
244 : PDRSCR, dPDRSCR_dlnRS, dPDRSCR_dlnGAME
245 : logical, intent(out) :: skip
246 : integer, intent(out) :: ierr
247 :
248 : integer :: iRS, jGAME
249 0 : real(dp) :: lnRS, lnRS0, lnRS1, lnGAME, lnGAME0, lnGAME1
250 0 : real(dp), dimension(nvals_FSCRliq8) :: fval, df_dlnRS, df_dlnGAME
251 : type (eosPC_Support_Info), pointer :: fq
252 : include 'formats'
253 0 : ierr = 0
254 0 : skip = .false.
255 :
256 0 : if (iZion < 1 .or. iZion > max_FSCRliq8_Zion) then
257 0 : write(*,2) 'invalid value for Z ion in get_FSCRliq8', iZion
258 0 : ierr = -1
259 0 : return
260 : end if
261 :
262 0 : fq => FSCRliq8_data(iZion)
263 :
264 0 : if (.not. FSCRliq8_Zion_loaded(iZion)) then
265 0 : !$OMP CRITICAL(eosPC_support_FSCRliq8_load)
266 0 : if (.not. FSCRliq8_Zion_loaded(iZion)) then
267 0 : call load_eosPC_support_Info(fq,iZion,ierr)
268 0 : if (ierr == 0) FSCRliq8_Zion_loaded(iZion) = .true.
269 : end if
270 : !$OMP END CRITICAL(eosPC_support_FSCRliq8_load)
271 0 : if (ierr /= 0) return
272 : end if
273 :
274 : ! get interpolation results
275 0 : lnRS = log(RS)
276 0 : lnGAME = log(GAME)
277 : if (lnRS < fq% lnRS_min .or. lnRS > fq% lnRS_max .or. &
278 0 : lnGAME < fq% lnGAME_min .or. lnGAME > fq% lnGAME_max) then
279 0 : skip = .true.
280 0 : return
281 : end if
282 :
283 : call Locate_lnRS(lnRS, &
284 : fq% nlnRS, fq% lnRS_min, fq% dlnRS, &
285 0 : iRS, lnRS0, lnRS1)
286 : call Locate_lnGAME(lnGAME, &
287 : fq% nlnGAME, fq% lnGAME_min, fq% dlnGAME, &
288 0 : jGAME, lnGAME0, lnGAME1)
289 :
290 : call Do_Interpolations( &
291 : 1, nvals_FSCRliq8, nvals_FSCRliq8, nvals_FSCRliq8, &
292 : fq% nlnRS, fq% lnRSs, fq% nlnGAME, fq% lnGAMEs, fq% tbl1, &
293 : iRS, jGAME, lnRS0, lnRS, lnRS1, lnGAME0, lnGAME, lnGAME1, &
294 0 : fval, df_dlnRS, df_dlnGAME, ierr)
295 0 : if (ierr /= 0) then
296 0 : write(*,1) 'Do_Interpolations failed in get_FSCRliq8'
297 0 : return
298 : end if
299 :
300 0 : FSCR = fval(iFSCR); dFSCR_dlnRS = df_dlnRS(iFSCR); dFSCR_dlnGAME = df_dlnGAME(iFSCR)
301 0 : USCR = fval(iUSCR); dUSCR_dlnRS = df_dlnRS(iUSCR); dUSCR_dlnGAME = df_dlnGAME(iUSCR)
302 0 : PSCR = fval(iPSCR); dPSCR_dlnRS = df_dlnRS(iPSCR); dPSCR_dlnGAME = df_dlnGAME(iPSCR)
303 0 : CVSCR = fval(iCVSCR); dCVSCR_dlnRS = df_dlnRS(iCVSCR); dCVSCR_dlnGAME = df_dlnGAME(iCVSCR)
304 0 : PDTSCR = fval(iPDTSCR); dPDTSCR_dlnRS = df_dlnRS(iPDTSCR); dPDTSCR_dlnGAME = df_dlnGAME(iPDTSCR)
305 0 : PDRSCR = fval(iPDRSCR); dPDRSCR_dlnRS = df_dlnRS(iPDRSCR); dPDRSCR_dlnGAME = df_dlnGAME(iPDRSCR)
306 :
307 : end subroutine get_FSCRliq8
308 :
309 :
310 0 : subroutine get_EXCOR7(RS, GAME, &
311 : FXC, dFXC_dlnRS, dFXC_dlnGAME, &
312 : UXC, dUXC_dlnRS, dUXC_dlnGAME, &
313 : PXC, dPXC_dlnRS, dPXC_dlnGAME, &
314 : CVXC, dCVXC_dlnRS, dCVXC_dlnGAME, &
315 : SXC, dSXC_dlnRS, dSXC_dlnGAME, &
316 : PDTXC, dPDTXC_dlnRS, dPDTXC_dlnGAME, &
317 : PDRXC, dPDRXC_dlnRS, dPDRXC_dlnGAME, &
318 : skip, ierr)
319 : real(dp), intent(in) :: RS, GAME
320 : real(dp), intent(out) :: &
321 : FXC, dFXC_dlnRS, dFXC_dlnGAME, &
322 : UXC, dUXC_dlnRS, dUXC_dlnGAME, &
323 : PXC, dPXC_dlnRS, dPXC_dlnGAME, &
324 : CVXC, dCVXC_dlnRS, dCVXC_dlnGAME, &
325 : SXC, dSXC_dlnRS, dSXC_dlnGAME, &
326 : PDTXC, dPDTXC_dlnRS, dPDTXC_dlnGAME, &
327 : PDRXC, dPDRXC_dlnRS, dPDRXC_dlnGAME
328 : logical, intent(out) :: skip
329 : integer, intent(out) :: ierr
330 :
331 : integer :: iRS, jGAME
332 0 : real(dp) :: lnRS, lnRS0, lnRS1, lnGAME, lnGAME0, lnGAME1
333 0 : real(dp), dimension(nvals_EXCOR7) :: fval, df_dlnRS, df_dlnGAME
334 : type (eosPC_Support_Info), pointer :: fq
335 : include 'formats'
336 0 : ierr = 0
337 0 : skip = .false.
338 :
339 0 : fq => EXCOR7_data
340 :
341 0 : if (.not. EXCOR7_table_loaded) then
342 0 : !$OMP CRITICAL(eosPC_EXCOR7_support_load)
343 0 : if (.not. EXCOR7_table_loaded) then
344 0 : call load_eosPC_support_Info(fq,0,ierr)
345 0 : if (ierr == 0) EXCOR7_table_loaded = .true.
346 : end if
347 : !$OMP END CRITICAL(eosPC_EXCOR7_support_load)
348 0 : if (ierr /= 0) return
349 : end if
350 :
351 : ! get interpolation results
352 0 : lnRS = log(RS)
353 0 : lnGAME = log(GAME)
354 : if (lnRS < fq% lnRS_min .or. lnRS > fq% lnRS_max .or. &
355 0 : lnGAME < fq% lnGAME_min .or. lnGAME > fq% lnGAME_max) then
356 0 : skip = .true.
357 0 : return
358 : end if
359 :
360 : call Locate_lnRS(lnRS, &
361 0 : fq% nlnRS, fq% lnRS_min, fq% dlnRS, iRS, lnRS0, lnRS1)
362 : call Locate_lnGAME(lnGAME, &
363 0 : fq% nlnGAME, fq% lnGAME_min, fq% dlnGAME, jGAME, lnGAME0, lnGAME1)
364 :
365 : call Do_Interpolations( &
366 : 1, nvals_EXCOR7, nvals_EXCOR7, nvals_EXCOR7, &
367 : fq% nlnRS, fq% lnRSs, fq% nlnGAME, fq% lnGAMEs, &
368 : fq% tbl1, iRS, jGAME, lnRS0, lnRS, lnRS1, lnGAME0, lnGAME, lnGAME1, &
369 0 : fval, df_dlnRS, df_dlnGAME, ierr)
370 0 : if (ierr /= 0) then
371 0 : write(*,1) 'Do_Interpolations failed in get_EXCOR7'
372 0 : return
373 : end if
374 :
375 0 : FXC = fval(jFXC); dFXC_dlnRS = df_dlnRS(jFXC); dFXC_dlnGAME = df_dlnGAME(jFXC)
376 0 : UXC = fval(jUXC); dUXC_dlnRS = df_dlnRS(jUXC); dUXC_dlnGAME = df_dlnGAME(jUXC)
377 0 : PXC = fval(jPXC); dPXC_dlnRS = df_dlnRS(jPXC); dPXC_dlnGAME = df_dlnGAME(jPXC)
378 0 : CVXC = fval(jCVXC); dCVXC_dlnRS = df_dlnRS(jCVXC); dCVXC_dlnGAME = df_dlnGAME(jCVXC)
379 0 : SXC = fval(jSXC); dSXC_dlnRS = df_dlnRS(jSXC); dSXC_dlnGAME = df_dlnGAME(jSXC)
380 0 : PDTXC = fval(jPDTXC); dPDTXC_dlnRS = df_dlnRS(jPDTXC); dPDTXC_dlnGAME = df_dlnGAME(jPDTXC)
381 0 : PDRXC = fval(jPDRXC); dPDRXC_dlnRS = df_dlnRS(jPDRXC); dPDRXC_dlnGAME = df_dlnGAME(jPDRXC)
382 :
383 : end subroutine get_EXCOR7
384 :
385 :
386 0 : subroutine Locate_lnRS( &
387 : lnRS, nlnRS, lnRS_min, dlnRS, iQ, lnRS0, lnRS1)
388 : real(dp), intent(inout) :: lnRS
389 : integer, intent(in) :: nlnRS
390 : real(dp), intent(in) :: lnRS_min, dlnRS
391 : integer, intent(out) :: iQ
392 : real(dp), intent(out) :: lnRS0, lnRS1
393 0 : iQ = int((lnRS - lnRS_min)/dlnRS + 1d-4) + 1
394 0 : if (iQ < 1 .or. iQ >= nlnRS) then
395 0 : if (iQ < 1) then
396 0 : iQ = 1
397 0 : lnRS0 = lnRS_min
398 0 : lnRS1 = lnRS0 + dlnRS
399 0 : lnRS = lnRS0
400 : else
401 0 : iQ = nlnRS-1
402 0 : lnRS0 = lnRS_min + (iQ-1)*dlnRS
403 0 : lnRS1 = lnRS0 + dlnRS
404 0 : lnRS = lnRS1
405 : end if
406 : else
407 0 : lnRS0 = lnRS_min + (iQ-1)*dlnRS
408 0 : lnRS1 = lnRS0 + dlnRS
409 : end if
410 0 : end subroutine Locate_lnRS
411 :
412 :
413 0 : subroutine Locate_lnGAME( &
414 : lnGAME, nlnGAME, lnGAME_min, dlnGAME, iQ, lnGAME0, lnGAME1)
415 : real(dp), intent(inout) :: lnGAME
416 : integer, intent(in) :: nlnGAME
417 : real(dp), intent(in) :: lnGAME_min, dlnGAME
418 : integer, intent(out) :: iQ
419 : real(dp), intent(out) :: lnGAME0, lnGAME1
420 0 : iQ = int((lnGAME - lnGAME_min)/dlnGAME + 1d-4) + 1
421 0 : if (iQ < 1 .or. iQ >= nlnGAME) then
422 0 : if (iQ < 1) then
423 0 : iQ = 1
424 0 : lnGAME0 = lnGAME_min
425 0 : lnGAME1 = lnGAME0 + dlnGAME
426 0 : lnGAME = lnGAME0
427 : else
428 0 : iQ = nlnGAME-1
429 0 : lnGAME0 = lnGAME_min + (iQ-1)*dlnGAME
430 0 : lnGAME1 = lnGAME0 + dlnGAME
431 0 : lnGAME = lnGAME1
432 : end if
433 : else
434 0 : lnGAME0 = lnGAME_min + (iQ-1)*dlnGAME
435 0 : lnGAME1 = lnGAME0 + dlnGAME
436 : end if
437 0 : end subroutine Locate_lnGAME
438 :
439 :
440 0 : subroutine Do_Interpolations( &
441 : nvlo, nvhi, nvals, n, nlnGAMI, x, ny, y, fin1, i, j, &
442 : x0, xget, x1, y0, yget, y1, &
443 0 : fval, df_dx, df_dy, ierr)
444 : integer, intent(in) :: nvlo, nvhi, nvals, n, nlnGAMI, ny
445 : real(dp), intent(in) :: x(:) ! (nlnGAMI)
446 : real(dp), intent(in) :: y(:) ! (ny)
447 : real(dp), intent(in), pointer :: fin1(:) ! = (4,n,nlnGAMI,ny)
448 : integer, intent(in) :: i, j ! target cell in f
449 : real(dp), intent(in) :: x0, xget, x1 ! x0 <= xget <= x1; x0 = xs(i), x1 = xs(i+1)
450 : real(dp), intent(in) :: y0, yget, y1 ! y0 <= yget <= y1; y0 = ys(j), y1 = ys(j+1)
451 : real(dp), intent(inout), dimension(nvals) :: fval, df_dx, df_dy
452 : integer, intent(out) :: ierr
453 :
454 0 : real(dp) :: xp, xpi, xp2, xpi2, ax, axbar, bx, bxbar, cx, cxi, hx2, cxd, cxdi, hx, hxi
455 0 : real(dp) :: yp, ypi, yp2, ypi2, ay, aybar, by, bybar, cy, cyi, hy2, cyd, cydi, hy, hyi
456 0 : real(dp) :: sixth_hx2, sixth_hy2, z36th_hx2_hy2
457 0 : real(dp) :: sixth_hx, sixth_hxi_hy2, z36th_hx_hy2
458 0 : real(dp) :: sixth_hx2_hyi, sixth_hy, z36th_hx2_hy
459 : integer :: k, ip1, jp1
460 0 : real(dp), pointer :: fin(:,:,:,:)
461 :
462 : include 'formats'
463 :
464 0 : ierr = 0
465 :
466 : fin(1:4,1:n,1:nlnGAMI,1:ny) => &
467 0 : fin1(1:4*n*nlnGAMI*ny)
468 :
469 0 : hx=x1-x0
470 0 : hxi=1d0/hx
471 0 : hx2=hx*hx
472 :
473 0 : xp=(xget-x0)*hxi
474 :
475 0 : xpi=1d0-xp
476 0 : xp2=xp*xp
477 0 : xpi2=xpi*xpi
478 :
479 0 : ax=xp2*(3d0-2d0*xp)
480 0 : axbar=1d0-ax
481 :
482 0 : bx=-xp2*xpi
483 0 : bxbar=xpi2*xp
484 :
485 0 : cx=xp*(xp2-1d0)
486 0 : cxi=xpi*(xpi2-1d0)
487 0 : cxd=3d0*xp2-1d0
488 0 : cxdi=-3d0*xpi2+1d0
489 :
490 0 : hy=y1-y0
491 0 : hyi=1d0/hy
492 0 : hy2=hy*hy
493 :
494 0 : yp=(yget-y0)*hyi
495 :
496 0 : ypi=1d0-yp
497 0 : yp2=yp*yp
498 0 : ypi2=ypi*ypi
499 :
500 0 : ay=yp2*(3d0-2d0*yp)
501 0 : aybar=1d0-ay
502 :
503 0 : by=-yp2*ypi
504 0 : bybar=ypi2*yp
505 :
506 0 : cy=yp*(yp2-1d0)
507 0 : cyi=ypi*(ypi2-1d0)
508 0 : cyd=3d0*yp2-1d0
509 0 : cydi=-3d0*ypi2+1d0
510 :
511 0 : sixth_hx2 = one_sixth*hx2
512 0 : sixth_hy2 = one_sixth*hy2
513 0 : z36th_hx2_hy2 = sixth_hx2*sixth_hy2
514 :
515 0 : sixth_hx = one_sixth*hx
516 0 : sixth_hxi_hy2 = one_sixth*hxi*hy2
517 0 : z36th_hx_hy2 = sixth_hx*sixth_hy2
518 :
519 0 : sixth_hx2_hyi = sixth_hx2*hyi
520 0 : sixth_hy = one_sixth*hy
521 0 : z36th_hx2_hy = sixth_hx2*sixth_hy
522 :
523 0 : ip1 = i+1
524 0 : jp1 = j+1
525 :
526 0 : !$omp simd
527 : do k = nvlo, nvhi
528 : ! bicubic spline interpolation
529 :
530 : ! f(1,i,j) = f(x(i),y(j))
531 : ! f(2,i,j) = d2f/dx2(x(i),y(j))
532 : ! f(3,i,j) = d2f/dy2(x(i),y(j))
533 : ! f(4,i,j) = d4f/dx2dy2(x(i),y(j))
534 :
535 : fval(k) = &
536 : xpi*( &
537 : ypi*fin(1,k,i,j) +yp*fin(1,k,i,jp1)) &
538 : +xp*(ypi*fin(1,k,ip1,j)+yp*fin(1,k,ip1,jp1)) &
539 : +sixth_hx2*( &
540 : cxi*(ypi*fin(2,k,i,j) +yp*fin(2,k,i,jp1))+ &
541 : cx*(ypi*fin(2,k,ip1,j)+yp*fin(2,k,ip1,jp1))) &
542 : +sixth_hy2*( &
543 : xpi*(cyi*fin(3,k,i,j) +cy*fin(3,k,i,jp1))+ &
544 : xp*(cyi*fin(3,k,ip1,j)+cy*fin(3,k,ip1,jp1))) &
545 : +z36th_hx2_hy2*( &
546 : cxi*(cyi*fin(4,k,i,j) +cy*fin(4,k,i,jp1))+ &
547 0 : cx*(cyi*fin(4,k,ip1,j)+cy*fin(4,k,ip1,jp1)))
548 :
549 : ! derivatives of bicubic splines
550 : df_dx(k) = &
551 : hxi*( &
552 : -(ypi*fin(1,k,i,j) +yp*fin(1,k,i,jp1)) &
553 : +(ypi*fin(1,k,ip1,j)+yp*fin(1,k,ip1,jp1))) &
554 : +sixth_hx*( &
555 : cxdi*(ypi*fin(2,k,i,j) +yp*fin(2,k,i,jp1))+ &
556 : cxd*(ypi*fin(2,k,ip1,j)+yp*fin(2,k,ip1,jp1))) &
557 : +sixth_hxi_hy2*( &
558 : -(cyi*fin(3,k,i,j) +cy*fin(3,k,i,jp1)) &
559 : +(cyi*fin(3,k,ip1,j)+cy*fin(3,k,ip1,jp1))) &
560 : +z36th_hx_hy2*( &
561 : cxdi*(cyi*fin(4,k,i,j) +cy*fin(4,k,i,jp1))+ &
562 0 : cxd*(cyi*fin(4,k,ip1,j)+cy*fin(4,k,ip1,jp1)))
563 :
564 : df_dy(k) = &
565 : hyi*( &
566 : xpi*(-fin(1,k,i,j) +fin(1,k,i,jp1))+ &
567 : xp*(-fin(1,k,ip1,j)+fin(1,k,ip1,jp1))) &
568 : +sixth_hx2_hyi*( &
569 : cxi*(-fin(2,k,i,j) +fin(2,k,i,jp1))+ &
570 : cx*(-fin(2,k,ip1,j)+fin(2,k,ip1,jp1))) &
571 : +sixth_hy*( &
572 : xpi*(cydi*fin(3,k,i,j) +cyd*fin(3,k,i,jp1))+ &
573 : xp*(cydi*fin(3,k,ip1,j)+cyd*fin(3,k,ip1,jp1))) &
574 : +z36th_hx2_hy*( &
575 : cxi*(cydi*fin(4,k,i,j) +cyd*fin(4,k,i,jp1))+ &
576 0 : cx*(cydi*fin(4,k,ip1,j)+cyd*fin(4,k,ip1,jp1)))
577 :
578 : end do
579 :
580 0 : end subroutine Do_Interpolations
581 :
582 :
583 0 : subroutine load_eosPC_support_Info(fq, iZion, ierr)
584 : use const_def, only: mesa_data_dir
585 : use utils_lib, only: alloc_iounit, free_iounit
586 : use create_EXCOR7_table, only: do_create_EXCOR7_table
587 : use create_FSCRliq8_table, only: do_create_FSCRliq8_table
588 : type (eosPC_Support_Info), pointer :: fq
589 : integer, intent(in) :: iZion ! 0 means EXCOR7
590 : integer, intent(out) :: ierr
591 :
592 : integer :: io_unit, n, j, i, iQ, nparams, nvals
593 : character (len=256) :: filename, fname, cache_fname, temp_cache_fname
594 : character (len=1000) :: message
595 0 : real(dp), pointer :: tbl2_1(:), tbl2(:,:,:)
596 0 : real(dp), target :: vec_ary(50)
597 : real(dp), pointer :: vec(:)
598 0 : real(dp) :: lnGAME, lnRS
599 :
600 : include 'formats'
601 0 : ierr = 0
602 0 : vec => vec_ary
603 :
604 0 : if (iZion == 0) then
605 0 : filename = 'EXCOR7'
606 0 : else if (iZion < 10) then
607 0 : write(filename,'(a,i1)') 'FSCRliq8_Zion_0', iZion
608 : else
609 0 : write(filename,'(a,i2)') 'FSCRliq8_Zion_', iZion
610 : end if
611 0 : fname = trim(mesa_data_dir) // '/eosPC_support_data/' // trim(filename) // '.data'
612 0 : cache_fname = trim(mesa_data_dir) // '/eosPC_support_data/cache/' // trim(filename) // '.bin'
613 0 : temp_cache_fname = trim(mesa_temp_caches_dir) // '/' // trim(filename) // '.bin'
614 :
615 0 : io_unit = alloc_iounit(ierr)
616 0 : if (ierr /= 0) then
617 0 : write(*,'(a)') 'failed to alloc_iounit'
618 0 : call mesa_error(__FILE__,__LINE__)
619 : end if
620 0 : open(unit=io_unit, FILE=trim(fname), ACTION='READ', STATUS='OLD', IOSTAT=ierr)
621 0 : if (ierr /= 0) then ! need to open the table file
622 0 : ierr = 0
623 0 : if (iZion == 0) then
624 0 : call do_create_EXCOR7_table(fname)
625 : else
626 0 : call do_create_FSCRliq8_table(fname,iZion)
627 : end if
628 0 : open(unit=io_unit, FILE=trim(fname), ACTION='READ', STATUS='OLD', IOSTAT=ierr)
629 0 : if (ierr /= 0) then
630 0 : write(*,'(a)') 'failed to create ' // trim(fname)
631 0 : write(*,'(A)')
632 0 : call mesa_error(__FILE__,__LINE__)
633 : end if
634 : end if
635 :
636 0 : read(io_unit,*,iostat=ierr) ! skip line 1
637 0 : if (ierr /= 0) return
638 :
639 0 : if (iZion == 0) then
640 0 : nparams = 8
641 0 : nvals = nvals_EXCOR7
642 : else
643 0 : nparams = 9
644 0 : nvals = nvals_FSCRliq8
645 : end if
646 0 : read(io_unit,'(a)',iostat=ierr) message ! get parameters line
647 0 : if (ierr == 0) call str_to_vector(message, vec, n, ierr)
648 0 : if (ierr /= 0 .or. n < nparams) then
649 0 : write(*,'(a)') 'failed while reading ' // trim(fname)
650 0 : close(io_unit)
651 0 : ierr = -1
652 0 : return
653 : end if
654 :
655 0 : read(io_unit,*,iostat=ierr) ! line 3
656 0 : if (ierr /= 0) return
657 :
658 0 : fq% Zion = iZion
659 0 : fq% nvals = nvals
660 0 : fq% nlnRS = int(vec(1))
661 0 : fq% lnRS_min = vec(2)*ln10
662 0 : fq% lnRS_max = vec(3)*ln10
663 0 : fq% dlnRS = vec(4)*ln10
664 :
665 0 : fq% nlnGAME = int(vec(5))
666 0 : fq% lnGAME_min = vec(6)*ln10
667 0 : fq% lnGAME_max = vec(7)*ln10
668 0 : fq% dlnGAME = vec(8)*ln10
669 :
670 0 : if (iZion > 0 .and. int(vec(9)) /= iZion) then
671 0 : write(*,*) 'bad value for Zion in file', int(vec(9)), iZion
672 0 : close(io_unit)
673 0 : ierr = -1
674 0 : return
675 : end if
676 :
677 : if (show_allocations) write(*,2) 'allocate pc_support', &
678 : 4*nvals*fq% nlnRS*fq% nlnGAME + fq% nlnRS + fq% nlnGAME
679 : allocate(fq% tbl1(4*nvals*fq% nlnRS*fq% nlnGAME), &
680 0 : fq% lnRSs(fq% nlnRS), fq% lnGAMEs(fq% nlnGAME), STAT=ierr)
681 0 : if (ierr /= 0) return
682 : fq% tbl(1:4, 1:nvals, 1:fq% nlnRS, 1:fq% nlnGAME) => &
683 0 : fq% tbl1(1:4*nvals*fq% nlnRS*fq% nlnGAME)
684 :
685 0 : fq% lnRSs(1) = fq% lnRS_min
686 0 : do i = 2, fq% nlnRS-1
687 0 : fq% lnRSs(i) = fq% lnRSs(i-1) + fq% dlnRS
688 : end do
689 0 : fq% lnRSs(fq% nlnRS) = fq% lnRS_max
690 :
691 0 : fq% lnGAMEs(1) = fq% lnGAME_min
692 0 : do i = 2, fq% nlnGAME-1
693 0 : fq% lnGAMEs(i) = fq% lnGAMEs(i-1) + fq% dlnGAME
694 : end do
695 0 : fq% lnGAMEs(fq% nlnGAME) = fq% lnGAME_max
696 :
697 0 : if (use_cache_for_eos) then
698 0 : call Read_Cache(fq, cache_fname, ierr)
699 0 : if (ierr == 0) then ! got it from the cache
700 0 : close(io_unit)
701 0 : call free_iounit(io_unit)
702 0 : return
703 : end if
704 0 : ierr = 0
705 : end if
706 :
707 0 : allocate(tbl2_1(nvals*fq% nlnRS*fq% nlnGAME), STAT=ierr)
708 0 : if (ierr /= 0) return
709 :
710 : tbl2(1:nvals, 1:fq% nlnRS, 1:fq% nlnGAME) => &
711 0 : tbl2_1(1:nvals*fq% nlnRS*fq% nlnGAME)
712 :
713 0 : do iQ=1,fq% nlnRS
714 :
715 0 : read(io_unit,*,iostat=ierr)
716 0 : if (failed('skip line1')) return
717 :
718 0 : read(io_unit,'(a)',iostat=ierr) message
719 0 : if (ierr == 0) call str_to_double(message, vec(1), ierr)
720 0 : if (failed('read lnRS')) return
721 0 : lnRS = vec(1)
722 :
723 0 : read(io_unit,*,iostat=ierr)
724 0 : if (failed('skip line2')) return
725 :
726 0 : read(io_unit,*,iostat=ierr)
727 0 : if (failed('skip line3')) return
728 :
729 0 : do i=1,fq% nlnGAME
730 :
731 0 : read(io_unit,'(a)',iostat=ierr) message
732 0 : if (failed('read line')) then
733 0 : write(*,'(a)') trim(message)
734 0 : write(*,*) trim(fname)
735 0 : write(*,*) 'iQ, i', iQ, i
736 0 : write(*,*) 'lnRS', lnRS
737 0 : write(*,*) 'bad input line?'
738 0 : call mesa_error(__FILE__,__LINE__)
739 : end if
740 :
741 0 : call str_to_vector(message, vec, n, ierr)
742 0 : if (ierr /= 0 .or. n < 1+nvals) then
743 0 : write(*,'(a)') trim(message)
744 0 : write(*,*) trim(fname)
745 0 : write(*,*) 'iQ, i', iQ, i
746 0 : write(*,*) 'lnRS', lnRS
747 0 : write(*,*) 'bad input line?'
748 0 : call mesa_error(__FILE__,__LINE__)
749 : end if
750 0 : lnGAME = vec(1)
751 0 : do j=1,nvals
752 0 : tbl2(j,iQ,i) = vec(1+j)
753 : end do
754 :
755 : end do
756 :
757 0 : if(iQ == fq% nlnRS) exit
758 0 : read(io_unit,*,iostat=ierr)
759 0 : if (failed('skip line4')) return
760 0 : read(io_unit,*,iostat=ierr)
761 0 : if (failed('skip line5')) return
762 :
763 : end do
764 :
765 0 : close(io_unit)
766 :
767 0 : call Make_Interpolation_Data(fq, nvals, tbl2_1, ierr)
768 0 : deallocate(tbl2_1)
769 0 : if (failed('Make_Interpolation_Data')) return
770 :
771 0 : call Check_Interpolation_Data
772 :
773 0 : if (.not. use_cache_for_eos) then
774 0 : call free_iounit(io_unit)
775 0 : return
776 : end if
777 :
778 : open(unit=io_unit, &
779 : file=trim(switch_str(temp_cache_fname, cache_fname, use_mesa_temp_cache)), &
780 0 : iostat=ierr, action='write', form='unformatted')
781 0 : if (ierr == 0) then
782 0 : write(*,'(a)') ' write ' // trim(cache_fname)
783 : write(io_unit) &
784 0 : fq% Zion, fq% nlnGAME, fq% lnGAME_min, fq% lnGAME_max, fq% dlnGAME, &
785 0 : fq% nlnRS, fq% lnRS_min, fq% lnRS_max, fq% dlnRS
786 0 : write(io_unit) fq% tbl1(1:4*nvals*fq% nlnRS*fq% nlnGAME)
787 0 : close(io_unit)
788 0 : if (use_mesa_temp_cache) call mv(temp_cache_fname, cache_fname, .true.)
789 : end if
790 0 : call free_iounit(io_unit)
791 :
792 0 : if (iZion == 0) then
793 0 : EXCOR7_table_loaded = .true.
794 : else
795 0 : FSCRliq8_Zion_loaded(iZion) = .true.
796 : end if
797 :
798 : contains
799 :
800 0 : subroutine Check_Interpolation_Data
801 : integer :: i, j, iQ, jtemp
802 0 : do i = 1, 4
803 0 : do j = 1, nvals
804 0 : do iQ = 1, fq% nlnRS
805 0 : do jtemp = 1, fq% nlnGAME
806 0 : if (is_bad(fq% tbl(i,j,iQ,jtemp))) then
807 0 : fq% tbl(i,j,iQ,jtemp) = 0
808 : end if
809 : end do
810 : end do
811 : end do
812 : end do
813 0 : end subroutine Check_Interpolation_Data
814 :
815 0 : logical function failed(str)
816 : character (len=*), intent(in) :: str
817 0 : failed = (ierr /= 0)
818 0 : if (failed) write(*,*) 'load_eosPC_support_Info failed: ' // trim(str)
819 0 : if (failed) call mesa_error(__FILE__,__LINE__,'load_eosPC_support_Info')
820 0 : end function failed
821 :
822 : end subroutine load_eosPC_support_Info
823 :
824 :
825 0 : subroutine Make_Interpolation_Data(fq, nvals, tbl2_1, ierr)
826 : use interp_2d_lib_db
827 : type (eosPC_Support_Info), pointer :: fq
828 : integer, intent(in) :: nvals
829 : real(dp), pointer :: tbl2_1(:) ! =(nvals, nlnRS, nlnGAME)
830 : integer, intent(out) :: ierr
831 :
832 0 : real(dp), allocatable, target :: f1_ary(:) ! data & spline coefficients
833 0 : real(dp), pointer :: f1(:), f(:,:,:), tbl2(:,:,:)
834 : integer :: ibcxmin ! bc flag for x=xmin
835 0 : real(dp) :: bcxmin(fq% nlnGAME) ! bc data vs. y at x=xmin
836 : integer :: ibcxmax ! bc flag for x=xmax
837 0 : real(dp) :: bcxmax(fq% nlnGAME) ! bc data vs. y at x=xmax
838 : integer :: ibcymin ! bc flag for y=ymin
839 0 : real(dp) :: bcymin(fq% nlnRS) ! bc data vs. x at y=ymin
840 : integer :: ibcymax ! bc flag for y=ymax
841 0 : real(dp) :: bcymax(fq% nlnRS) ! bc data vs. x at y=ymax
842 : integer :: ili_lnRSs ! =1: logRho grid is "nearly" equally spaced
843 : integer :: ili_lnGAMEs ! =1: lnGAME grid is "nearly" equally spaced
844 : integer :: ier ! =0 on exit if there is no error.
845 :
846 : integer :: v, i, j
847 :
848 : include 'formats'
849 :
850 0 : ierr = 0
851 :
852 : ! just use "not a knot" bc's at edges of tables
853 0 : ibcxmin = 0; bcxmin(:) = 0
854 0 : ibcxmax = 0; bcxmax(:) = 0
855 0 : ibcymin = 0; bcymin(:) = 0
856 0 : ibcymax = 0; bcymax(:) = 0
857 :
858 : tbl2(1:nvals, 1:fq% nlnRS, 1:fq% nlnGAME) => &
859 0 : tbl2_1(1:nvals*fq% nlnRS*fq% nlnGAME)
860 : ! copy file variables to internal interpolation tables
861 0 : do j=1,fq% nlnGAME
862 0 : do i=1,fq% nlnRS
863 0 : do v = 1, nvals
864 0 : fq% tbl(1,v,i,j) = tbl2(v,i,j)
865 : end do
866 : end do
867 : end do
868 :
869 0 : allocate(f1_ary(4 * fq% nlnRS * fq% nlnGAME))
870 :
871 0 : f1 => f1_ary
872 : f(1:4, 1:fq% nlnRS, 1:fq% nlnGAME) => &
873 0 : f1_ary(1:4*fq% nlnRS*fq% nlnGAME)
874 :
875 : ! create tables for bicubic spline interpolation
876 0 : do v = 1, nvals
877 0 : do i=1,fq% nlnRS
878 0 : do j=1,fq% nlnGAME
879 0 : f(1,i,j) = fq% tbl(1,v,i,j)
880 : end do
881 : end do
882 : call interp_mkbicub_db( &
883 : fq% lnRSs,fq% nlnRS,fq% lnGAMEs,fq% nlnGAME,f1,fq% nlnRS, &
884 : ibcxmin,bcxmin,ibcxmax,bcxmax, &
885 : ibcymin,bcymin,ibcymax,bcymax, &
886 0 : ili_lnRSs,ili_lnGAMEs,ier)
887 0 : if (ier /= 0) then
888 0 : write(*,*) 'interp_mkbicub_db error happened for EXCOR7 value', v
889 0 : ierr = 3
890 0 : return
891 : end if
892 0 : do i=1,fq% nlnRS
893 0 : do j=1,fq% nlnGAME
894 0 : fq% tbl(2,v,i,j) = f(2,i,j)
895 0 : fq% tbl(3,v,i,j) = f(3,i,j)
896 0 : fq% tbl(4,v,i,j) = f(4,i,j)
897 : end do
898 : end do
899 : end do
900 :
901 0 : end subroutine Make_Interpolation_Data
902 :
903 :
904 0 : subroutine Read_Cache(fq, cache_fname, ierr)
905 : use utils_lib, only: alloc_iounit, free_iounit
906 : type (eosPC_Support_Info), pointer :: fq
907 : character (*), intent(in) :: cache_fname
908 : integer, intent(out) :: ierr
909 :
910 0 : real(dp) :: lnGAME_min_in, lnGAME_max_in, dlnGAME_in, &
911 0 : lnRS_min_in, lnRS_max_in, dlnRS_in
912 : integer :: Zion_in, nlnRS_in, nlnGAME_in, io_unit
913 : real(dp), parameter :: tiny = 1d-10
914 :
915 : include 'formats'
916 :
917 : ierr = 0
918 :
919 0 : io_unit = alloc_iounit(ierr)
920 0 : if (ierr /= 0) then
921 0 : write(*,'(a)') 'failed to alloc_iounit'
922 0 : call mesa_error(__FILE__,__LINE__)
923 : end if
924 : open(unit=io_unit, file=trim(cache_fname), action='read', &
925 0 : status='old', iostat=ierr, form='unformatted')
926 0 : if (ierr /= 0) then
927 0 : call free_iounit(io_unit)
928 0 : return
929 : end if
930 :
931 : !write(*,'(a)') 'read ' // trim(cache_fname)
932 :
933 : read(io_unit, iostat=ierr) &
934 0 : Zion_in, nlnGAME_in, lnGAME_min_in, lnGAME_max_in, dlnGAME_in, &
935 0 : nlnRS_in, lnRS_min_in, lnRS_max_in, dlnRS_in
936 0 : if (ierr /= 0) then
937 0 : write(*,*) 'read cache failed'
938 : end if
939 :
940 0 : if (fq% Zion /= Zion_in) then
941 0 : ierr = 1
942 0 : write(*,*) 'read cache failed for Zion'
943 : end if
944 0 : if (fq% nlnRS /= nlnRS_in) then
945 0 : ierr = 1
946 0 : write(*,*) 'read cache failed for nlnRS'
947 : end if
948 0 : if (fq% nlnGAME /= nlnGAME_in) then
949 0 : ierr = 1
950 0 : write(*,*) 'read cache failed for nlnGAME'
951 : end if
952 0 : if (abs(fq% lnGAME_min-lnGAME_min_in) > tiny) then
953 0 : ierr = 1
954 0 : write(*,*) 'read cache failed for eos_lnGAME_min'
955 : end if
956 0 : if (abs(fq% lnGAME_max-lnGAME_max_in) > tiny) then
957 0 : ierr = 1
958 0 : write(*,*) 'read cache failed for eos_lnGAME_max'
959 : end if
960 0 : if (abs(fq% dlnGAME-dlnGAME_in) > tiny) then
961 0 : ierr = 1
962 0 : write(*,*) 'read cache failed for eos_dlnGAME'
963 : end if
964 0 : if (abs(fq% lnRS_min-lnRS_min_in) > tiny) then
965 0 : ierr = 1
966 0 : write(*,*) 'read cache failed for eos_lnRS_min'
967 : end if
968 0 : if (abs(fq% lnRS_max-lnRS_max_in) > tiny) then
969 0 : ierr = 1
970 0 : write(*,*) 'read cache failed for eos_lnRS_max'
971 : end if
972 0 : if (abs(fq% dlnRS-dlnRS_in) > tiny) then
973 0 : ierr = 1
974 0 : write(*,*) 'read cache failed for eos_dlnRS'
975 : end if
976 :
977 0 : if (ierr /= 0) then
978 0 : write(*,*) 'read cache file failed 1'
979 0 : call mesa_error(__FILE__,__LINE__,'Read_Cache')
980 0 : close(io_unit); return
981 : end if
982 :
983 : read(io_unit, iostat=ierr) &
984 0 : fq% tbl1(1:4*fq% nvals*fq% nlnRS*fq% nlnGAME)
985 0 : if (ierr /= 0) then
986 0 : write(*,*) 'read cache file failed 2'
987 0 : call mesa_error(__FILE__,__LINE__,'Read_Cache')
988 : end if
989 :
990 0 : close(io_unit)
991 0 : call free_iounit(io_unit)
992 :
993 : end subroutine Read_Cache
994 :
995 :
996 : ! ================== ELECTRON-ION COULOMB LIQUID =================== *
997 0 : subroutine FITION9(GAMI, &
998 : FION,UION,PION,CVii,PDTii,PDRii)
999 : ! Version 11.09.08
1000 : ! Non-ideal contributions to thermodynamic functions of classical OCP,
1001 : ! corrected at small density for a mixture.
1002 : ! Stems from FITION00 v.24.05.00.
1003 : ! Input: GAMI - ion coupling parameter
1004 : ! Output: FION - ii free energy / N_i kT
1005 : ! UION - ii internal energy / N_i kT
1006 : ! PION - ii pressure / n_i kT
1007 : ! CVii - ii heat capacity / N_i k
1008 : ! PDTii = PION + d(PION)/d ln T = (1/N_i kT)*(d P_{ii}/d ln T)
1009 : ! PDRii = PION + d(PION)/d ln\rho
1010 : ! Parameters adjusted to Caillol (1999).
1011 : real(dp), intent(in) :: GAMI
1012 : real(dp), intent(out) :: FION, UION, PION, CVii, PDTii, PDRii
1013 0 : real(dp) :: F0, U0
1014 :
1015 : real(dp), parameter :: A1=-.907347d0
1016 : real(dp), parameter :: A2=.62849d0
1017 : real(dp) :: A3
1018 : real(dp), parameter :: C1=.004500d0
1019 : real(dp), parameter :: G1=170.0d0
1020 : real(dp), parameter :: C2=-8.4d-5
1021 : real(dp), parameter :: G2=.0037d0
1022 : real(dp), parameter :: SQ32=.8660254038d0 ! SQ32=sqrt(3)/2
1023 :
1024 0 : A3=-SQ32-A1/sqrt(A2)
1025 : F0=A1*(sqrt(GAMI*(A2+GAMI)) &
1026 : - A2*log(sqrt(GAMI/A2)+sqrt(1d0+GAMI/A2))) &
1027 0 : + 2d0*A3*(sqrt(GAMI)-atan(sqrt(GAMI)))
1028 0 : U0=pow3(sqrt(GAMI))*(A1/sqrt(A2+GAMI)+A3/(1.d0+GAMI))
1029 : ! This is the zeroth approximation. Correction:
1030 0 : UION=U0+C1*GAMI*GAMI/(G1+GAMI)+C2*GAMI*GAMI/(G2+GAMI*GAMI)
1031 : FION=F0+C1*(GAMI-G1*log(1.d0+GAMI/G1)) &
1032 0 : + C2/2d0*log(1.d0+GAMI*GAMI/G2)
1033 : CVii=-0.5d0*pow3(sqrt(GAMI))*(A1*A2/pow3(sqrt(A2+GAMI)) &
1034 : + A3*(1.d0-GAMI)/pow2(1.d0+GAMI)) &
1035 0 : - GAMI*GAMI*(C1*G1/pow2(G1+GAMI)+C2*(G2-GAMI*GAMI)/pow2(G2+GAMI*GAMI))
1036 0 : PION=UION/3.0d0
1037 0 : PDRii=(4.0d0*UION-CVii)/9.0d0 ! p_{ii} + d p_{ii} / d ln\rho
1038 0 : PDTii=CVii/3.0d0 ! p_{ii} + d p_{ii} / d ln T
1039 0 : return
1040 : end subroutine FITION9
1041 :
1042 :
1043 : end module pc_support
|
