Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2021 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 eos_def
21 :
22 : use const_def, only: dp, use_mesa_temp_cache, strlen
23 : use chem_def, only: max_el_z
24 :
25 : implicit none
26 :
27 : ! interfaces for procedure pointers
28 : abstract interface
29 :
30 : subroutine other_eos_frac_interface( &
31 : handle, &
32 : species, chem_id, net_iso, xa, &
33 : Rho, logRho, T, logT, &
34 : frac, dfrac_dlogRho, dfrac_dlogT, &
35 : ierr)
36 : use const_def, only: dp
37 : implicit none
38 : integer, intent(in) :: handle ! eos handle; from star, pass s% eos_handle
39 : integer, intent(in) :: species
40 : integer, pointer :: chem_id(:) ! maps species to chem id
41 : ! index from 1 to species
42 : ! value is between 1 and num_chem_isos
43 : integer, pointer :: net_iso(:) ! maps chem id to species number
44 : ! index from 1 to num_chem_isos (defined in chem_def)
45 : ! value is 0 if the iso is not in the current net
46 : ! else is value between 1 and number of species in current net
47 : real(dp), intent(in) :: xa(:) ! mass fractions
48 :
49 : real(dp), intent(in) :: Rho, logRho ! the density
50 : real(dp), intent(in) :: T, logT ! the temperature
51 :
52 : real(dp), intent(out) :: frac ! fraction of other_eos to use
53 : real(dp), intent(out) :: dfrac_dlogRho ! its partial derivative at constant T
54 : real(dp), intent(out) :: dfrac_dlogT ! its partial derivative at constant Rho
55 :
56 : integer, intent(out) :: ierr ! 0 means AOK.
57 : end subroutine other_eos_frac_interface
58 :
59 : subroutine other_eos_interface( &
60 : handle, &
61 : species, chem_id, net_iso, xa, &
62 : Rho, logRho, T, logT, &
63 : res, d_dlnd, d_dlnT, d_dxa, ierr)
64 : use const_def, only: dp
65 : implicit none
66 : integer, intent(in) :: handle ! eos handle; from star, pass s% eos_handle
67 :
68 : integer, intent(in) :: species
69 : integer, pointer :: chem_id(:) ! maps species to chem id
70 : integer, pointer :: net_iso(:) ! maps chem id to species number
71 : real(dp), intent(in) :: xa(:) ! mass fractions
72 :
73 : real(dp), intent(in) :: Rho, logRho ! the density
74 : real(dp), intent(in) :: T, logT ! the temperature
75 :
76 : real(dp), intent(inout) :: res(:) ! (num_eos_basic_results)
77 : real(dp), intent(inout) :: d_dlnd(:) ! (num_eos_basic_results)
78 : real(dp), intent(inout) :: d_dlnT(:) ! (num_eos_basic_results)
79 : real(dp), intent(inout) :: d_dxa(:,:) ! (num_eos_basic_results,species)
80 :
81 : integer, intent(out) :: ierr ! 0 means AOK.
82 : end subroutine other_eos_interface
83 :
84 : end interface
85 :
86 :
87 : logical, parameter :: show_allocations = .false. ! for debugging memory usage
88 : integer, parameter :: eos_name_length = 20 ! String length for storing EOS variable names
89 :
90 :
91 : ! cgs units
92 :
93 : ! the basic eos results
94 :
95 : integer, parameter :: i_lnPgas = 1
96 : ! gas pressure (total pressure minus radiation pressure)
97 : integer, parameter :: i_lnE = i_lnPgas+1
98 : ! internal energy per gram
99 : integer, parameter :: i_lnS = i_lnE+1
100 : ! entropy per gram
101 : integer, parameter :: i_mu = i_lnS+1
102 : ! mean molecular weight per gas particle (ions + free electrons)
103 : integer, parameter :: i_lnfree_e = i_mu+1
104 : ! free_e := total combined number per nucleon of free electrons
105 : integer, parameter :: i_eta = i_lnfree_e+1
106 : ! electron degeneracy parameter (eta > 1 for significant degeneracy)
107 : ! eta = ratio of electron chemical potential to kT
108 : integer, parameter :: i_grad_ad = i_eta+1
109 : ! dlnT_dlnP at constant S
110 : integer, parameter :: i_chiRho = i_grad_ad+1
111 : ! dlnP_dlnRho at constant T
112 : integer, parameter :: i_chiT = i_chiRho+1
113 : ! dlnP_dlnT at constant Rho
114 : integer, parameter :: i_Cp = i_chiT+1
115 : ! dh_dT at constant P, specific heat at constant total pressure
116 : ! where h is enthalpy, h = E + P/Rho
117 : integer, parameter :: i_Cv = i_Cp+1
118 : ! dE_dT at constant Rho, specific heat at constant volume
119 : integer, parameter :: i_dE_dRho = i_Cv+1
120 : ! at constant T
121 : integer, parameter :: i_dS_dT = i_dE_dRho+1
122 : ! at constant Rho
123 : integer, parameter :: i_dS_dRho = i_dS_dT+1
124 : ! at constant T
125 : integer, parameter :: i_gamma1 = i_dS_dRho+1
126 : ! dlnP_dlnRho at constant S
127 : integer, parameter :: i_gamma3 = i_gamma1+1
128 : ! gamma3 - 1 = dlnT_dlnRho at constant S
129 : integer, parameter :: i_phase = i_gamma3+1
130 : ! phase: 1 for solid, 0 for liquid, in-between for blend
131 : integer, parameter :: i_latent_ddlnT = i_phase+1
132 : ! T*dS/dlnT from phase transition
133 : integer, parameter :: i_latent_ddlnRho = i_latent_ddlnT+1
134 : ! T*dS/dlnRho from phase transition
135 :
136 : ! blend information
137 : integer, parameter :: i_eos_HELM = 1
138 : integer, parameter :: i_frac_HELM = i_latent_ddlnRho+1
139 : ! fraction HELM
140 : integer, parameter :: i_eos_OPAL_SCVH = 2
141 : integer, parameter :: i_frac_OPAL_SCVH = i_frac_HELM+1
142 : ! fraction OPAL/SCVH
143 : integer, parameter :: i_eos_FreeEOS = 3
144 : integer, parameter :: i_frac_FreeEOS = i_frac_OPAL_SCVH+1
145 : ! fraction FreeEOS
146 : integer, parameter :: i_eos_PC = 4
147 : integer, parameter :: i_frac_PC = i_frac_FreeEOS+1
148 : ! fraction PC
149 : integer, parameter :: i_eos_Skye = 5
150 : integer, parameter :: i_frac_Skye = i_frac_PC+1
151 : ! fraction Skye
152 : integer, parameter :: i_eos_CMS = 6
153 : integer, parameter :: i_frac_CMS = i_frac_Skye+1
154 : ! fraction CMS
155 : integer, parameter :: i_eos_ideal = 7
156 : integer, parameter :: i_frac_ideal = i_frac_CMS+1
157 : ! fraction ideal ion+radiation
158 :
159 : ! eos frac entries correspond to the slice
160 : ! i_frac:i_frac+num_eos_frac_results-1
161 : integer, parameter :: i_frac = i_frac_HELM ! first frac entry
162 : integer, parameter :: num_eos_frac_results = 7
163 :
164 : integer, parameter :: num_eos_basic_results = i_frac_ideal
165 : integer, parameter :: nv = num_eos_basic_results
166 :
167 : ! only return d_dxa of lnE and lnPgas to star
168 : integer, parameter :: num_eos_d_dxa_results = 2
169 :
170 : character (len=eos_name_length) :: eosDT_result_names(nv)
171 :
172 : ! non-positive indices for special EOS quantities
173 : ! these are supported in the root-finds (eosDT_get_{Rho,T})
174 : ! even though their values are not basic EOS results
175 :
176 : integer, parameter :: i_logPtot = 0 ! log10 total pressure (gas + radiation)
177 : integer, parameter :: i_egas = -1 ! gas specific energy density (no radiation)
178 :
179 :
180 : ! NOTE: the calculation of eta is based on the following equation for ne,
181 : ! the mean number of free electrons per cm^3,
182 : ! assuming non-relativistic electrons (okay for T < 10^9 at least)
183 : !
184 : ! ne = 4 Pi / h^3 (2 me k T)^1.5 F[1/2,eta] -- see, for example, Clayton, eqn 2-57
185 : ! where F is the fermi-dirac integral:
186 : ! F[beta,eta] := Integrate[(u^beta)/(1+E^(u-eta)),{u,0,Infinity}]
187 : !
188 : ! CAVEAT: when free_e, the mean number of free electrons per nucleon gets really small,
189 : ! eta isn't very interesting because there aren't a lot of free electrons to be degenerate!
190 : ! our calculation of eta gets flaky at this point as well.
191 : ! we sweep this problem under the rug by making eta tend to a fairly large negative value
192 : ! when free_e < 0.01 or so. this roughly corresponds to T < 10^4 or less.
193 :
194 : integer, parameter :: eosdt_OPAL_SCVH = 1
195 : integer, parameter :: eosdt_max_FreeEOS = 2
196 :
197 : integer, parameter :: max_num_DT_Zs = 15
198 : integer, parameter :: max_num_DT_Xs = 12
199 :
200 : integer, parameter :: num_eosDT_Zs = 3
201 : integer, parameter :: num_eosDT_Xs = 6
202 :
203 : integer, parameter :: num_FreeEOS_Zs = 15
204 : integer, parameter :: num_FreeEOS_Xs = 12
205 :
206 : type DT_XZ_Info
207 : integer :: nZs
208 : real(dp) :: Zs(max_num_DT_Zs)
209 : integer :: nXs_for_Z(max_num_DT_Zs)
210 : real(dp) :: Xs_for_Z(max_num_DT_Xs, max_num_DT_Zs)
211 : end type DT_XZ_Info
212 :
213 : type (DT_XZ_Info), target :: eosDT_XZ_struct, FreeEOS_XZ_struct
214 :
215 : integer, parameter :: sz_per_eos_point = 4 ! for bicubic spline interpolation
216 :
217 : type EoS_General_Info
218 :
219 : ! limits for HELM
220 : real(dp) :: Z_all_HELM ! all HELM for Z >= this unless eos_use_FreeEOS
221 : real(dp) :: logT_all_HELM ! all HELM for lgT >= this
222 : real(dp) :: logT_low_all_HELM ! all HELM for lgT <= this
223 : real(dp) :: coulomb_temp_cut_HELM, coulomb_den_cut_HELM
224 :
225 : ! limits for OPAL_SCVH
226 : logical :: use_OPAL_SCVH
227 : real(dp) :: logT_low_all_SCVH ! SCVH for lgT >= this
228 : real(dp) :: logT_all_OPAL ! OPAL for lgT <= this
229 : real(dp) :: logRho1_OPAL_SCVH_limit ! don't use OPAL_SCVH for logRho > this
230 : real(dp) :: logRho2_OPAL_SCVH_limit ! full OPAL_SCVH okay for logRho < this
231 : real(dp) :: logRho_min_OPAL_SCVH_limit ! no OPAL/SCVH for logRho < this
232 : real(dp) :: logQ_max_OPAL_SCVH ! no OPAL/SCVH for logQ > this
233 : real(dp) :: logQ_min_OPAL_SCVH ! no OPAL/SCVH for logQ <= this.
234 : real(dp) :: Z_all_OPAL ! all OPAL for Z <= this
235 :
236 : ! limits for FreeEOS
237 : logical :: use_FreeEOS
238 : real(dp) :: logQ_max_FreeEOS_hi
239 : real(dp) :: logQ_max_FreeEOS_lo
240 : real(dp) :: logQ_min_FreeEOS_hi
241 : real(dp) :: logQ_min_FreeEOS_lo
242 : real(dp) :: logRho_min_FreeEOS_hi
243 : real(dp) :: logRho_min_FreeEOS_lo
244 : real(dp) :: logRho_max_FreeEOS_hi
245 : real(dp) :: logRho_max_FreeEOS_lo
246 : real(dp) :: logT_min_FreeEOS_hi
247 : real(dp) :: logT_min_FreeEOS_lo
248 : real(dp) :: logT_max_FreeEOS_hi
249 : real(dp) :: logT_max_FreeEOS_lo
250 : real(dp) :: logQ_cut_FreeEOS_lo_Z_max
251 : real(dp) :: logQ_cut_lo_Z_FreeEOS_hi
252 : real(dp) :: logQ_cut_lo_Z_FreeEOS_lo
253 : real(dp) :: logQ_cut_hi_Z_FreeEOS_hi
254 : real(dp) :: logQ_cut_hi_Z_FreeEOS_lo
255 : real(dp) :: logRho_cut_FreeEOS_hi
256 : real(dp) :: logRho_cut_FreeEOS_lo
257 : real(dp) :: logT_cut_FreeEOS_hi
258 : real(dp) :: logT_cut_FreeEOS_lo
259 : character (len=30) :: suffix_for_FreeEOS_Z(num_FreeEOS_Zs)
260 :
261 : ! limits for CMS
262 : logical :: use_CMS, CMS_use_fixed_composition
263 : integer :: CMS_fixed_composition_index ! in [0,10]
264 : real(dp) :: max_Z_for_any_CMS, max_Z_for_all_CMS ! set to -1 to disable CMS
265 : real(dp) :: logQ_max_for_any_CMS, logQ_max_for_all_CMS ! for upper blend zone in logQ = logRho - 2*logT + 12
266 : real(dp) :: logQ_min_for_any_CMS, logQ_min_for_all_CMS ! for lower blend zone in logQ
267 : real(dp) :: logRho_max_for_all_CMS, logRho_max_for_any_CMS ! for upper blend zone in logRho
268 : real(dp) :: logRho_min_for_all_CMS, logRho_min_for_any_CMS ! for lower blend zone in logRho
269 : real(dp) :: logT_max_for_all_CMS, logT_max_for_any_CMS ! for upper blend zone in logT
270 : real(dp) :: logT_min_for_all_CMS, logT_min_for_any_CMS ! for lower blend zone in logT
271 : real(dp) :: logT_max_for_all_CMS_pure_He, logT_max_for_any_CMS_pure_He ! upper logT blend zone is different for pure He
272 :
273 : ! limits for PC
274 : logical :: use_PC
275 : real(dp) :: mass_fraction_limit_for_PC ! skip any species with abundance < this
276 : real(dp) :: logRho1_PC_limit ! okay for pure PC for logRho > this
277 : real(dp) :: logRho2_PC_limit ! don't use PC for logRho < this (>= 2.8)
278 : logical :: PC_use_Gamma_limit_instead_of_T
279 : real(dp) :: logT1_PC_limit ! okay for pure PC for logT < this (like logT_all_OPAL)
280 : real(dp) :: logT2_PC_limit ! don't use PC for logT > this (like logT_all_HELM)
281 : real(dp) :: log_Gamma_e_all_HELM ! HELM for log_Gamma_e <= this
282 : real(dp) :: Gamma_e_all_HELM ! 10**log_Gamma_e_all_HELM
283 : real(dp) :: log_Gamma_e_all_PC ! PC for log_Gamma_e >= this
284 : real(dp) :: tiny_fuzz
285 : ! crystallization boundaries
286 : real(dp) :: PC_Gamma_start_crystal ! Begin releasing latent heat of crystallization
287 : real(dp) :: PC_Gamma_full_crystal ! Fully into the solid phase
288 :
289 : ! limits for Skye
290 : logical :: use_Skye
291 : logical :: Skye_use_ion_offsets
292 : real(dp) :: mass_fraction_limit_for_Skye
293 : real(dp) :: Skye_min_gamma_for_solid ! The minimum Gamma_i at which to use the solid free energy fit (below this, extrapolate).
294 : real(dp) :: Skye_max_gamma_for_liquid ! The maximum Gamma_i at which to use the liquid free energy fit (above this, extrapolate).
295 : character(len=128) :: Skye_solid_mixing_rule ! Currently support 'Ogata' or 'PC'
296 :
297 : logical :: use_simple_Skye_blends
298 : real(dp) :: logRho_min_for_any_Skye, logRho_min_for_all_Skye
299 : real(dp) :: logT_min_for_any_Skye, logT_min_for_all_Skye
300 :
301 : ! misc
302 : logical :: include_radiation, include_elec_pos
303 : logical :: eosDT_use_linear_interp_for_X
304 : logical :: eosDT_use_linear_interp_to_HELM
305 :
306 : character(len=128) :: eosDT_file_prefix
307 :
308 : logical :: okay_to_convert_ierr_to_skip
309 :
310 : ! other eos
311 : logical :: use_other_eos_component
312 : procedure (other_eos_frac_interface), pointer, nopass :: &
313 : other_eos_frac => null()
314 : procedure (other_eos_interface), pointer, nopass :: &
315 : other_eos_component => null()
316 : logical :: use_other_eos_results
317 : procedure (other_eos_interface), pointer, nopass :: &
318 : other_eos_results => null()
319 :
320 : ! debugging
321 : logical :: dbg
322 : real(dp) :: logT_lo, logT_hi
323 : real(dp) :: logRho_lo, logRho_hi
324 : real(dp) :: X_lo, X_hi
325 : real(dp) :: Z_lo, Z_hi
326 :
327 : ! bookkeeping
328 : integer :: handle
329 : logical :: in_use
330 :
331 :
332 : ! User supplied inputs
333 : real(dp) :: eos_ctrl(10)
334 : integer :: eos_integer_ctrl(10)
335 : logical :: eos_logical_ctrl(10)
336 : character(len=strlen) :: eos_character_ctrl(10)
337 :
338 : end type EoS_General_Info
339 :
340 :
341 : include 'helm_def.dek'
342 :
343 :
344 : ! THE FOLLOWING ARE PRIVATE DEFS -- NOT FOR USE BY CLIENTS
345 :
346 :
347 : ! data table types
348 :
349 : type (HELM_Table), pointer :: eos_ht
350 :
351 : ! for mesa (logQ,logT) tables
352 : type EosDT_XZ_Info
353 : real(dp) :: logQ_min ! logQ = logRho - 2*logT + 12
354 : real(dp) :: logQ_max
355 : real(dp) :: del_logQ ! spacing for the logQs
356 : integer :: num_logQs
357 : real(dp) :: logT_min
358 : real(dp) :: logT_max
359 : real(dp) :: del_logT ! spacing for the logTs
360 : integer :: num_logTs
361 : real(dp), pointer :: logQs(:), logTs(:)
362 : real(dp), pointer :: tbl1(:)
363 : integer :: version
364 : end type EosDT_XZ_Info
365 :
366 :
367 : ! data table variables
368 : type (EosDT_XZ_Info), dimension(num_eosDT_Xs, num_eosDT_Zs), target :: &
369 : eosDT_XZ_data, eosSCVH_XZ_data, eosCMS_XZ_data
370 : type (EosDT_XZ_Info), dimension(num_FreeEOS_Xs, num_FreeEOS_Zs), target :: &
371 : FreeEOS_XZ_data
372 :
373 : logical, dimension(num_eosDT_Xs, num_eosDT_Zs) :: &
374 : eosDT_XZ_loaded, eosSCVH_XZ_loaded, eosCMS_XZ_loaded
375 : logical, dimension(num_FreeEOS_Xs, num_FreeEOS_Zs) :: FreeEOS_XZ_loaded
376 :
377 :
378 : ! interpolation info for eosPC support tables FITION9
379 : type FITION_Info
380 : real(dp), pointer :: f1(:), f(:,:)
381 : end type FITION_Info
382 : integer, parameter :: FITION_vals = 6
383 : real(dp), pointer :: FITION9_lnGAMIs(:)
384 : type (FITION_Info), target :: FITION_data(FITION_vals)
385 : logical :: FITION9_loaded ! all get created at once
386 :
387 : ! interpolation info for eosPC support tables, FSCRliq8 and EXCOR7
388 : type eosPC_Support_Info
389 : integer :: Zion, nlnRS, nlnGAME, nvals
390 : real(dp) :: lnRS_min, lnRS_max, dlnRS
391 : real(dp) :: lnGAME_min, lnGAME_max, dlnGAME
392 : real(dp), pointer :: lnRSs(:) ! (nlnRS)
393 : real(dp), pointer :: lnGAMEs(:) ! (nlnGAME)
394 : real(dp), pointer :: tbl1(:) ! (4,nvals,nlnRS,nlnGAME)
395 : real(dp), pointer :: tbl(:,:,:,:) ! => tbl1(:)
396 : end type eosPC_Support_Info
397 :
398 : integer, parameter :: max_FSCRliq8_Zion = max_el_z
399 : type (eosPC_Support_Info), target :: FSCRliq8_data(max_FSCRliq8_Zion)
400 : logical, dimension(max_FSCRliq8_Zion) :: FSCRliq8_Zion_loaded
401 :
402 : type (eosPC_Support_Info), target :: EXCOR7_data
403 : logical :: EXCOR7_table_loaded
404 :
405 : integer, parameter :: max_eos_handles = 10
406 : type (EoS_General_Info), target :: eos_handles(max_eos_handles)
407 :
408 : logical :: use_cache_for_eos = .false.
409 : logical :: eos_root_is_initialized = .false.
410 :
411 : character(len=1000) :: eosDT_cache_dir
412 : character(len=1000) :: eosDT_temp_cache_dir
413 :
414 : logical :: eos_test_partials
415 : real(dp) :: eos_test_partials_val, eos_test_partials_dval_dx ! for dfridr from star
416 :
417 :
418 : contains
419 :
420 :
421 13 : subroutine eos_def_init
422 : integer :: i
423 : type (DT_XZ_Info), pointer :: eosDT_XZ_ptr, FreeEOS_XZ_ptr
424 : include 'formats'
425 :
426 13 : use_cache_for_eos = .false.
427 13 : eos_root_is_initialized = .false.
428 13 : eos_test_partials = .false.
429 :
430 13 : EXCOR7_table_loaded = .false.
431 13 : FSCRliq8_Zion_loaded(:) = .false.
432 13 : FITION9_loaded = .false.
433 :
434 143 : do i=1,max_eos_handles
435 130 : eos_handles(i)% handle = i
436 143 : eos_handles(i)% in_use = .false.
437 : end do
438 :
439 13 : eosDT_XZ_ptr => eosDT_XZ_struct
440 13 : eosDT_XZ_ptr% nZs = num_eosDT_Zs
441 52 : eosDT_XZ_ptr% Zs(1:num_eosDT_Zs) = [ 0.00d0, 0.02d0, 0.04d0 ]
442 52 : eosDT_XZ_ptr% nXs_for_Z(1:num_eosDT_Zs) = [ 6, 5, 5 ]
443 52 : do i=1,num_eosDT_Zs
444 : eosDT_XZ_ptr% Xs_for_Z(1:num_eosDT_Xs,i) = &
445 286 : [ 0.0d0, 0.2d0, 0.4d0, 0.6d0, 0.8d0, 1.0d0 ]
446 : end do
447 :
448 13 : FreeEOS_XZ_ptr => FreeEOS_XZ_struct
449 13 : FreeEOS_XZ_ptr% nZs = num_FreeEOS_Zs
450 : FreeEOS_XZ_ptr% Zs(1:num_FreeEOS_Zs) = &
451 : [ 0.00d0, 0.02d0, 0.04d0, 0.06d0, 0.08d0, 0.1d0, 0.2d0, &
452 208 : 0.3d0, 0.4d0, 0.5d0, 0.6d0, 0.7d0, 0.8d0, 0.9d0, 1.0d0 ]
453 : FreeEOS_XZ_ptr% nXs_for_Z(1:num_FreeEOS_Zs) = &
454 208 : [ 11, 11, 11, 11, 11, 4, 3, 3, 3, 3, 3, 3, 3, 2, 1 ]
455 78 : do i=1,5 ! 0.0 to 0.08
456 : FreeEOS_XZ_ptr% Xs_for_Z(1:11,i) = &
457 : [ 0.0d0, 0.1d0, 0.2d0, 0.3d0, 0.4d0, 0.5d0, &
458 793 : 0.6d0, 0.7d0, 0.8d0, 0.9d0, 1.0d0-FreeEOS_XZ_ptr% Zs(i) ]
459 : end do
460 65 : FreeEOS_XZ_ptr% Xs_for_Z(1:4,6) = [ 0.0d0, 0.1d0, 0.4d0, 0.9d0 ] ! 0.1
461 104 : do i=7,13 ! 0.2, 0.8
462 : FreeEOS_XZ_ptr% Xs_for_Z(1:3,i) = &
463 377 : [ 0.0d0, 0.1d0, 1.0d0-FreeEOS_XZ_ptr% Zs(i) ]
464 : end do
465 39 : FreeEOS_XZ_ptr% Xs_for_Z(1:2,14) = [ 0.0d0, 0.1d0 ] ! 0.9
466 13 : FreeEOS_XZ_ptr% Xs_for_Z(1,15) = 0.0d0 ! 1.0
467 :
468 13 : eosDT_XZ_loaded(:,:) = .false.
469 13 : eosSCVH_XZ_loaded(:,:)=.false.
470 13 : FreeEOS_XZ_loaded(:,:)=.false.
471 13 : eosCMS_XZ_loaded(:,:)=.false.
472 :
473 13 : eosDT_result_names(i_lnPgas) = 'lnPgas'
474 13 : eosDT_result_names(i_lnE) = 'lnE'
475 13 : eosDT_result_names(i_lnS) = 'lnS'
476 13 : eosDT_result_names(i_grad_ad) = 'grad_ad'
477 13 : eosDT_result_names(i_chiRho) = 'chiRho'
478 13 : eosDT_result_names(i_chiT) = 'chiT'
479 13 : eosDT_result_names(i_Cp) = 'Cp'
480 13 : eosDT_result_names(i_Cv) = 'Cv'
481 13 : eosDT_result_names(i_dE_dRho) = 'dE_dRho'
482 13 : eosDT_result_names(i_dS_dT) = 'dS_dT'
483 13 : eosDT_result_names(i_dS_dRho) = 'dS_dRho'
484 13 : eosDT_result_names(i_mu) = 'mu'
485 13 : eosDT_result_names(i_lnfree_e) = 'lnfree_e'
486 13 : eosDT_result_names(i_gamma1) = 'gamma1'
487 13 : eosDT_result_names(i_gamma3) = 'gamma3'
488 13 : eosDT_result_names(i_eta) = 'eta'
489 13 : eosDT_result_names(i_phase) = 'phase'
490 13 : eosDT_result_names(i_latent_ddlnT) = 'latent_ddlnT'
491 13 : eosDT_result_names(i_latent_ddlnRho) = 'latent_ddlnRho'
492 13 : eosDT_result_names(i_frac_OPAL_SCVH) = 'OPAL/SCVH'
493 13 : eosDT_result_names(i_frac_HELM) = 'HELM'
494 13 : eosDT_result_names(i_frac_Skye) = 'Skye'
495 13 : eosDT_result_names(i_frac_PC) = 'PC'
496 13 : eosDT_result_names(i_frac_FreeEOS) = 'FreeEOS'
497 13 : eosDT_result_names(i_frac_CMS) = 'CMS'
498 13 : eosDT_result_names(i_frac_ideal) = 'ideal'
499 :
500 13 : end subroutine eos_def_init
501 :
502 :
503 19 : integer function do_alloc_eos(ierr) result(handle)
504 : integer, intent(out) :: ierr
505 : integer :: i
506 19 : ierr = 0
507 19 : handle = -1
508 38 : !$omp critical (eos_handle)
509 27 : do i = 1, max_eos_handles
510 27 : if (.not. eos_handles(i)% in_use) then
511 19 : eos_handles(i)% in_use = .true.
512 19 : handle = i
513 19 : exit
514 : end if
515 : end do
516 : !$omp end critical (eos_handle)
517 19 : if (handle == -1) then
518 0 : ierr = -1
519 0 : return
520 : end if
521 19 : if (eos_handles(handle)% handle /= handle) then
522 0 : ierr = -1
523 0 : return
524 : end if
525 19 : call init_eos_handle_data(handle)
526 19 : end function do_alloc_eos
527 :
528 :
529 19 : subroutine init_eos_handle_data(handle)
530 : use other_eos
531 : use math_lib
532 : integer, intent(in) :: handle
533 : type (EoS_General_Info), pointer :: rq
534 19 : rq => eos_handles(handle)
535 19 : rq% in_use = .true.
536 19 : rq% handle = handle
537 :
538 19 : rq% other_eos_frac => null_other_eos_frac
539 19 : rq% other_eos_component => null_other_eos_component
540 19 : rq% other_eos_results => null_other_eos_results
541 :
542 19 : end subroutine init_eos_handle_data
543 :
544 :
545 5 : subroutine do_free_eos_handle(handle)
546 : integer, intent(in) :: handle
547 : type (EoS_General_Info), pointer :: rq
548 5 : if (handle >= 1 .and. handle <= max_eos_handles) then
549 5 : rq => eos_handles(handle)
550 5 : eos_handles(handle)% in_use = .false.
551 : end if
552 19 : end subroutine do_free_eos_handle
553 :
554 :
555 221488 : subroutine get_eos_ptr(handle,rq,ierr)
556 : integer, intent(in) :: handle
557 : type (EoS_General_Info), pointer :: rq
558 : integer, intent(out):: ierr
559 221488 : if (handle < 1 .or. handle > max_eos_handles) then
560 0 : ierr = -1
561 0 : return
562 : end if
563 221488 : rq => eos_handles(handle)
564 221488 : ierr = 0
565 : end subroutine get_eos_ptr
566 :
567 :
568 5 : subroutine eos_def_shutdown
569 : type (eosPC_Support_Info), pointer :: fq
570 : type (FITION_Info), pointer :: fi
571 : integer :: iz
572 :
573 : ! DT tables
574 : call free_EosDT_XZ_Info( &
575 5 : eosDT_XZ_data, eosDT_XZ_loaded, num_eosDT_Xs, num_eosDT_Zs)
576 : call free_EosDT_XZ_Info( &
577 5 : eosSCVH_XZ_data, eosSCVH_XZ_loaded, num_eosDT_Xs, num_eosDT_Zs)
578 : call free_EosDT_XZ_Info( &
579 5 : eosCMS_XZ_data, eosCMS_XZ_loaded, num_eosDT_Xs, num_eosDT_Zs)
580 : call free_EosDT_XZ_Info( &
581 5 : FreeEOS_XZ_data, FreeEOS_XZ_loaded, num_FreeEOS_Xs, num_FreeEOS_Zs)
582 :
583 : ! Misc. stuff
584 :
585 5 : if (FITION9_loaded) then
586 0 : do iz = 1, FITION_vals
587 0 : fi => FITION_data(iz)
588 0 : if (ASSOCIATED(fi% f1)) deallocate(fi% f1)
589 0 : nullify(fi% f)
590 : end do
591 0 : if (ASSOCIATED(FITION9_lnGAMIs)) deallocate(FITION9_lnGAMIs)
592 0 : FITION9_loaded = .false.
593 : end if
594 :
595 5 : fq => EXCOR7_data
596 5 : call free_eosPC_support_Info(fq)
597 5 : EXCOR7_table_loaded = .false.
598 :
599 565 : do iz = 1, max_FSCRliq8_Zion
600 560 : if (.not. FSCRliq8_Zion_loaded(iz)) cycle
601 0 : fq => FSCRliq8_data(iz)
602 565 : call free_eosPC_support_Info(fq)
603 : end do
604 5 : FSCRliq8_Zion_loaded(:) = .false.
605 :
606 5 : eos_root_is_initialized = .false.
607 :
608 :
609 : contains
610 :
611 5 : subroutine free_eosPC_support_Info(fq)
612 : type (eosPC_Support_Info), pointer :: fq
613 5 : if (ASSOCIATED(fq% lnRSs)) deallocate(fq% lnRSs)
614 5 : if (ASSOCIATED(fq% lnGAMEs)) deallocate(fq% lnGAMEs)
615 5 : if (ASSOCIATED(fq% tbl1)) deallocate(fq% tbl1)
616 5 : nullify(fq% tbl)
617 5 : end subroutine free_eosPC_support_Info
618 :
619 20 : subroutine free_EosDT_XZ_Info(d, flgs, numXs, numZs)
620 : integer, intent(in) :: numXs, numZs
621 : type (EosDT_XZ_Info), dimension(numXs, numZs) :: d
622 : logical, dimension(numXs, numZs) :: flgs
623 : integer :: ix, iz
624 170 : do ix = 1, numXs
625 1340 : do iz = 1, numZs
626 1170 : if (ASSOCIATED(d(ix,iz)% tbl1)) then
627 12 : deallocate(d(ix,iz)% tbl1)
628 : end if
629 1170 : if (ASSOCIATED(d(ix,iz)% logQs)) then
630 12 : deallocate(d(ix,iz)% logQs)
631 : end if
632 1320 : if (ASSOCIATED(d(ix,iz)% logTs)) then
633 12 : deallocate(d(ix,iz)% logTs)
634 : end if
635 : end do
636 : end do
637 1310 : flgs(:,:) = .false.
638 20 : end subroutine free_EosDT_XZ_Info
639 :
640 : end subroutine eos_def_shutdown
641 :
642 0 : end module eos_def
|
