MPI-AMRVAC 3.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_ionization_degree.t
Go to the documentation of this file.
1!> module ionization degree - get ionization degree for given temperature
2module mod_ionization_degree
3 implicit none
4 private
5 integer, parameter :: n_Tonly = 100
6 double precision :: HI_fraction_chromosphere(n_Tonly)
7 double precision :: T_H_fraction_K(n_Tonly)
8 double precision :: ionization_He_abundance = 0.1d0
9 double precision :: ionization_Rfactor_norm = 2.3d0
10 double precision, allocatable :: prom_T_table(:)
11 double precision, allocatable :: prom_p_table(:)
12 double precision, allocatable :: prom_T_eos_table(:)
13 double precision, allocatable :: prom_Y_eos_table(:,:)
14 double precision :: Te_table_min, Te_table_max, inv_Te_table_step
15 double precision :: Te_low_iz_He=5413.d0
16 double precision :: Y_table_min, Y_table_max
17 double precision :: ionization_H_energy_unit = 0.d0
18 double precision, dimension(:), allocatable :: Rfactor_table
19 double precision, dimension(:), allocatable :: Y_table
20 double precision, dimension(:), allocatable :: eps_ion_H_table
21 double precision, dimension(:), allocatable :: Te_H_table
22 double precision, dimension(:), allocatable :: iz_H_table
23 ! chromosphere: T-only ionization table based on
24 ! Carlsson & Leenaarts (2012, A&A, 539, A39).
25 integer, parameter :: ionization_mode_chromosphere = 1
26 ! prominence: pressure-dependent prominence ionization table based on
27 ! Heinzel et al. (2015, A&A, 579, A16) and
28 ! Gunár et al. (2025, A&A, 699, A89). The bundled table uses model A,
29 ! side illumination, and the L = 500 km voxel from the latter work.
30 integer, parameter :: ionization_mode_prominence = 2
31 ! flare: T-only hydrogen H I/H fraction table for flaring chromosphere
32 ! conditions from Hong, Carlsson & Ding (2022, A&A, 661, A77). It is
33 ! handled by the same T-only H-fraction machinery as the chromosphere
34 ! table.
35 integer, parameter :: ionization_mode_flare = 3
36 integer, parameter :: n_prom_T = 10
37 integer, parameter :: n_prom_p = 7
38 integer, parameter :: n_prom_eos = 4000
39 integer :: ionization_table_mode = ionization_mode_chromosphere
40 ! Initialize
41 public :: ionization_degree_init
42 ! Current rho,p -> complete EOS state
43 public :: ionization_get_state
44 public :: ionization_get_state_scalar
45 ! Prescribed T,p -> Rfactor and optional ionization degrees
46 public :: ionization_state_tp
47 ! Prescribed rho,T -> consistent pthermal and Rfactor
48 public :: ionization_solve_p_rfactor
49 ! T-only cooling callback and mode query
50 public :: ionization_get_rfactor_from_temperature
51 public :: ionization_is_temperature_only
52 public :: ionization_includes_energy
53 public :: ionization_get_state_from_eint
54 public :: ionization_get_p_eint_from_rho_t
55 public :: ionization_get_eps_derivative_t
56 public :: ionization_get_csound2_t
57 public :: ionization_check_eint_table
58 ! Compatibility interfaces
59 public :: ionization_state_from_temperature_scalar
60 ! Legacy compatibility API
61 public :: ionization_degree_from_temperature
62 logical :: ionization_initialized = .false.
63 logical :: ionization_include_energy = .false.
64
65 double precision, parameter :: prom_T_K(n_prom_T) = (/ &
66 6000.d0, 7000.d0, 8000.d0, 9000.d0, 10000.d0, &
67 12000.d0, 15000.d0, 20000.d0, 25000.d0, 30000.d0 /)
68
69 double precision, parameter :: prom_p_cgs(n_prom_p) = (/ &
70 0.03d0, 0.05d0, 0.10d0, 0.30d0, 0.50d0, 0.70d0, 1.00d0 /)
71
72 ! First index: temperature. Second index: pressure.
73 ! Data are ordered one complete temperature column per pressure.
74 double precision, parameter :: &
75 prom_iz_H_data(n_prom_T,n_prom_p) = reshape((/ &
76 ! p = 0.03 dyn cm^-2
77 0.53d0, 0.59d0, 0.65d0, 0.73d0, 0.80d0, &
78 0.90d0, 0.95d0, 0.98d0, 0.99d0, 1.00d0, &
79 ! p = 0.05 dyn cm^-2
80 0.40d0, 0.47d0, 0.56d0, 0.69d0, 0.79d0, &
81 0.91d0, 0.96d0, 0.98d0, 0.99d0, 1.00d0, &
82 ! p = 0.10 dyn cm^-2
83 0.25d0, 0.32d0, 0.47d0, 0.67d0, 0.81d0, &
84 0.93d0, 0.97d0, 0.99d0, 0.99d0, 1.00d0, &
85 ! p = 0.30 dyn cm^-2
86 0.10d0, 0.17d0, 0.41d0, 0.69d0, 0.85d0, &
87 0.96d0, 0.98d0, 0.99d0, 1.00d0, 1.00d0, &
88 ! p = 0.50 dyn cm^-2
89 0.06d0, 0.14d0, 0.38d0, 0.70d0, 0.87d0, &
90 0.97d0, 0.99d0, 1.00d0, 1.00d0, 1.00d0, &
91 ! p = 0.70 dyn cm^-2
92 0.04d0, 0.12d0, 0.36d0, 0.70d0, 0.88d0, &
93 0.97d0, 0.99d0, 1.00d0, 1.00d0, 1.00d0, &
94 ! p = 1.00 dyn cm^-2
95 0.03d0, 0.10d0, 0.34d0, 0.71d0, 0.89d0, &
96 0.98d0, 0.99d0, 1.00d0, 1.00d0, 1.00d0 /), &
97 (/ n_prom_t, n_prom_p /))
98
99 ! H I/H population fraction from Hong, Carlsson & Ding
100 ! (2022, A&A, 661, A77). Its temperature grid is identical to
101 ! T_H_fraction_K.
102 double precision, parameter :: HI_fraction_flare(n_Tonly) = (/ &
103 9.9999877d-1, 9.9999877d-1, 9.9999877d-1, 9.9999877d-1, 9.9999877d-1, &
104 9.9999877d-1, 9.9999877d-1, 9.9999877d-1, 9.9999877d-1, 9.9999877d-1, &
105 9.9999877d-1, 9.9999890d-1, 9.9998024d-1, 9.9992560d-1, 9.9968205d-1, &
106 9.9899492d-1, 9.9682125d-1, 9.9069732d-1, 9.8273722d-1, 9.5165553d-1, &
107 8.8401513d-1, 7.5808722d-1, 6.3272201d-1, 4.7164997d-1, 1.6497967d-1, &
108 8.6151552d-2, 4.9787498d-2, 2.3629016d-2, 1.1010071d-2, 7.7663095d-3, &
109 5.7739222d-3, 4.1811404d-3, 3.4427999d-3, 2.9108183d-3, 2.2295092d-3, &
110 1.9451997d-3, 1.4754925d-3, 1.1849689d-3, 9.1381928d-4, 7.1369446d-4, &
111 6.3154203d-4, 5.4204957d-4, 4.6800799d-4, 4.1090740d-4, 3.4051796d-4, &
112 2.9397681d-4, 2.5497294d-4, 2.0900410d-4, 1.7267453d-4, 1.3474650d-4, &
113 1.0894726d-4, 8.6230407d-5, 7.0545733d-5, 5.8882317d-5, 4.9311732d-5, &
114 3.9687148d-5, 3.1954017d-5, 2.6460585d-5, 2.1367807d-5, 1.7836169d-5, &
115 1.4207015d-5, 1.2022741d-5, 9.9158901d-6, 8.5561393d-6, 7.1994135d-6, &
116 6.1326905d-6, 5.0944881d-6, 4.3131760d-6, 3.6080551d-6, 3.0571391d-6, &
117 2.5778758d-6, 2.1649433d-6, 1.7960119d-6, 1.4711565d-6, 1.1832833d-6, &
118 9.3858978d-7, 7.4227325d-7, 5.8536296d-7, 4.6590828d-7, 3.8339449d-7, &
119 3.2851748d-7, 2.8778528d-7, 2.5213506d-7, 2.2042234d-7, 1.9440791d-7, &
120 1.7116454d-7, 1.5143655d-7, 1.3390303d-7, 1.1879927d-7, 1.0569364d-7, &
121 9.4140342d-8, 8.3767251d-8, 7.4822295d-8, 6.6863923d-8, 5.9836421d-8, &
122 5.3685738d-8, 4.8177839d-8, 4.3423310d-8, 3.9109303d-8, 1.0188512d-8 /)
123
124 ! Common temperature grid for the T-only H-fraction tables. The
125 ! Carlsson & Leenaarts (2012) and Hong et al. (2022) grids coincide.
126 data t_h_fraction_k /2.1000005d+03, 2.2349495d+03, 2.3785708d+03, 2.5314199d+03, 2.6940928d+03, &
127 2.8672190d+03, 3.0514692d+03, 3.2475613d+03, 3.4562544d+03, 3.6783584d+03, &
128 3.9147351d+03, 4.1662998d+03, 4.4340322d+03, 4.7189697d+03, 5.0222153d+03, &
129 5.3449502d+03, 5.6884248d+03, 6.0539712d+03, 6.4430083d+03, 6.8570420d+03, &
130 7.2976860d+03, 7.7666460d+03, 8.2657383d+03, 8.7969062d+03, 9.3622090d+03, &
131 9.9638330d+03, 1.0604130d+04, 1.1285561d+04, 1.2010781d+04, 1.2782619d+04, &
132 1.3604041d+04, 1.4478250d+04, 1.5408652d+04, 1.6398826d+04, 1.7452648d+04, &
133 1.8574172d+04, 1.9767766d+04, 2.1038082d+04, 2.2390010d+04, 2.3828838d+04, &
134 2.5360100d+04, 2.6989764d+04, 2.8724182d+04, 3.0570023d+04, 3.2534518d+04, &
135 3.4625215d+04, 3.6850262d+04, 3.9218293d+04, 4.1738543d+04, 4.4420699d+04, &
136 4.7275266d+04, 5.0313219d+04, 5.3546391d+04, 5.6987395d+04, 6.0649453d+04, &
137 6.4546914d+04, 6.8694758d+04, 7.3109148d+04, 7.7807289d+04, 8.2807258d+04, &
138 8.8128625d+04, 9.3791852d+04, 9.9819000d+04, 1.0623346d+05, 1.1306024d+05, &
139 1.2032560d+05, 1.2805797d+05, 1.3628709d+05, 1.4504503d+05, 1.5436592d+05, &
140 1.6428561d+05, 1.7484295d+05, 1.8607852d+05, 1.9803609d+05, 2.1076208d+05, &
141 2.2430608d+05, 2.3872045d+05, 2.5406084d+05, 2.7038703d+05, 2.8776234d+05, &
142 3.0625456d+05, 3.2593475d+05, 3.4688000d+05, 3.6917081d+05, 3.9289406d+05, &
143 4.1814181d+05, 4.4501247d+05, 4.7360988d+05, 5.0404450d+05, 5.3643488d+05, &
144 5.7090725d+05, 6.0759431d+05, 6.4663956d+05, 6.8819319d+05, 7.3241712d+05, &
145 7.7948294d+05, 8.2957412d+05, 8.8288331d+05, 9.3961919d+05, 1.0000000d+06 /
146
147 ! Chromospheric H I/H fraction from
148 ! Carlsson, M., & Leenaarts, J. 2012, A&A, 539, A39.
149 data hi_fraction_chromosphere /1.0000000d+00, 1.0000000d+00, 1.0000000d+00, 1.0000000d+00, 1.0000000d+00, &
150 1.0000000d+00, 1.0000000d+00, 9.9999928d-01, 9.9999774d-01, 9.9999624d-01, &
151 9.9999428d-01, 9.9999219d-01, 9.9998862d-01, 9.9998313d-01, 9.9997944d-01, &
152 9.9995774d-01, 1.0000324d+00, 9.9822283d-01, 9.6391600d-01, 8.8995951d-01, &
153 8.0025935d-01, 7.4219799d-01, 7.1378660d-01, 6.9565988d-01, 6.7326778d-01, &
154 6.5032905d-01, 6.2487972d-01, 5.9925246d-01, 5.7028764d-01, 5.3673893d-01, &
155 4.9144468d-01, 4.4332290d-01, 3.9158964d-01, 3.4621361d-01, 2.7826238d-01, &
156 1.9633104d-01, 1.2593637d-01, 8.1383914d-02, 5.0944358d-02, 2.8770180d-02, &
157 1.4449068d-02, 7.6071853d-03, 4.3617180d-03, 2.7694483d-03, 1.8341533d-03, &
158 1.2467797d-03, 8.6159195d-04, 6.0232455d-04, 4.2862241d-04, 3.0666622d-04, &
159 2.2256430d-04, 1.6328457d-04, 1.2179965d-04, 9.2065704d-05, 7.0357783d-05, &
160 5.4168402d-05, 4.2026968d-05, 3.2807184d-05, 2.5801779d-05, 2.0414085d-05, &
161 1.6271379d-05, 1.3041737d-05, 1.0538992d-05, 8.5329248d-06, 6.9645471d-06, &
162 5.7069005d-06, 4.7155136d-06, 3.9203474d-06, 3.2734958d-06, 2.7415674d-06, &
163 2.2911979d-06, 1.9045801d-06, 1.5692771d-06, 1.2782705d-06, 1.0281080d-06, &
164 8.1567373d-07, 6.4200310d-07, 5.0668103d-07, 4.0577518d-07, 3.3245340d-07, &
165 2.8006056d-07, 2.4119515d-07, 2.1033340d-07, 1.8429903d-07, 1.6189058d-07, &
166 1.4264889d-07, 1.2616209d-07, 1.1193421d-07, 9.9408616d-08, 8.8335305d-08, &
167 7.8593857d-08, 7.0042269d-08, 6.2532941d-08, 5.5898361d-08, 5.0049657d-08, &
168 4.4888207d-08, 4.0317172d-08, 3.6192560d-08, 3.2330732d-08, 2.9387850d-08 /
169
170 contains
171
172 subroutine ionization_degree_init(He_abundance, Rfactor_norm, table_name, include_energy)
173 use mod_global_parameters, only: unit_temperature, unit_numberdensity, &
174 unit_pressure, si_unit
175 use mod_constants, only: const_ev
176 use mod_comm_lib, only: mpistop
177 double precision, intent(in) :: he_abundance, rfactor_norm
178 character(len=*), intent(in), optional :: table_name
179 character(len=32) :: selected_table
180 logical, intent(in), optional :: include_energy
181 double precision, parameter :: erg_to_joule = 1.0d-7
182
183 if (ionization_initialized) then
184 call mpistop("ionization_degree_init called more than once")
185 end if
186
187 if (he_abundance < 0.d0) call mpistop("negative He abundance")
188 if (rfactor_norm <= 0.d0) call mpistop("invalid Rfactor normalization")
189
190 ionization_he_abundance = he_abundance
191 ionization_rfactor_norm = rfactor_norm
192 te_low_iz_he = 5413.d0/unit_temperature
193 ionization_include_energy = .false.
194 if (present(include_energy)) ionization_include_energy = include_energy
195
196 selected_table = "chromosphere"
197 if (present(table_name)) selected_table = trim(adjustl(table_name))
198
199 select case (trim(selected_table))
200 case ("chromosphere")
201 ionization_table_mode = ionization_mode_chromosphere
202 case ("flare")
203 ionization_table_mode = ionization_mode_flare
204 case ("prominence")
205 ionization_table_mode = ionization_mode_prominence
206 case default
207 call mpistop( &
208 "Unknown ionization table. Use 'chromosphere', 'flare', "// &
209 "or 'prominence'.")
210 end select
211
212 if (ionization_include_energy .and. &
213 .not. ionization_is_temperature_only()) then
214 call mpistop("Ionization energy is supported only for T-only tables "// &
215 "('chromosphere' and 'flare'), not for the pressure-dependent "// &
216 "'prominence' table.")
217 end if
218
219 if (ionization_include_energy) then
220 if (unit_numberdensity <= 0.d0 .or. unit_pressure <= 0.d0) then
221 call mpistop("invalid units for hydrogen ionization energy")
222 end if
223 if (si_unit) then
224 ! const_eV is in erg/eV; convert to joule for SI pressure units.
225 ionization_h_energy_unit = &
226 13.6d0*const_ev*erg_to_joule * &
227 unit_numberdensity/unit_pressure
228 else
229 ionization_h_energy_unit = &
230 13.6d0*const_ev*unit_numberdensity/unit_pressure
231 end if
232 else
233 ionization_h_energy_unit = 0.d0
234 end if
235
236 select case (ionization_table_mode)
237 case (ionization_mode_chromosphere)
238 call ionization_build_temperature_table(hi_fraction_chromosphere)
239 case (ionization_mode_flare)
240 call ionization_build_temperature_table(hi_fraction_flare)
241 case (ionization_mode_prominence)
242 call ionization_build_prominence_table()
243 end select
244 ionization_initialized = .true.
245 end subroutine ionization_degree_init
246
247 subroutine ionization_build_temperature_table(HI_fraction)
248 use mod_global_parameters, only: unit_temperature, zero
249
250 double precision, intent(in) :: hi_fraction(n_tonly)
251 double precision :: fact1, fact2, fact3
252 double precision :: dl1, dl2, te_table_step
253 double precision :: iz_h_source(n_tonly)
254 integer :: i, j, n_interpolated_table
255 logical :: jump
256
257 n_interpolated_table=4000
258
259 ! The source tables store H I/H; the EOS uses H II/H.
260 iz_h_source=1.d0-hi_fraction
261
262 allocate(te_h_table(1:n_interpolated_table))
263 allocate(iz_h_table(1:n_interpolated_table))
264 te_h_table=zero
265 iz_h_table=zero
266
267 te_table_step = (t_h_fraction_k(n_tonly)-t_h_fraction_k(1)) &
268 /dble(n_interpolated_table-1)
269
270 te_h_table(1) = t_h_fraction_k(1)
271 iz_h_table(1) = iz_h_source(1)
272
273 te_h_table(n_interpolated_table) = t_h_fraction_k(n_tonly)
274 iz_h_table(n_interpolated_table) = iz_h_source(n_tonly)
275
276 do i=2,n_interpolated_table
277 te_h_table(i) = te_h_table(i-1)+te_table_step
278 do j=1,n_tonly-1
279 ! Second order interpolation except at the outer edge or a jump.
280 if(te_h_table(i) < t_h_fraction_k(j+1)) then
281 if(j.eq. n_tonly-1) then
282 fact1 = (te_h_table(i)-t_h_fraction_k(j+1)) &
283 /(t_h_fraction_k(j)-t_h_fraction_k(j+1))
284
285 fact2 = (te_h_table(i)-t_h_fraction_k(j)) &
286 /(t_h_fraction_k(j+1)-t_h_fraction_k(j))
287
288 iz_h_table(i) = iz_h_source(j)*fact1 + &
289 iz_h_source(j+1)*fact2
290 exit
291 else
292 dl1 = iz_h_source(j+1)-iz_h_source(j)
293 dl2 = iz_h_source(j+2)-iz_h_source(j+1)
294 jump = (max(dabs(dl1),dabs(dl2)) > &
295 2.d0*min(dabs(dl1),dabs(dl2)))
296 end if
297
298 if(jump) then
299 fact1 = (te_h_table(i)-t_h_fraction_k(j+1)) &
300 /(t_h_fraction_k(j)-t_h_fraction_k(j+1))
301
302 fact2 = (te_h_table(i)-t_h_fraction_k(j)) &
303 /(t_h_fraction_k(j+1)-t_h_fraction_k(j))
304
305 iz_h_table(i) = iz_h_source(j)*fact1 + &
306 iz_h_source(j+1)*fact2
307 exit
308 else
309 fact1 = ((te_h_table(i)-t_h_fraction_k(j+1)) &
310 * (te_h_table(i)-t_h_fraction_k(j+2))) &
311 / ((t_h_fraction_k(j)-t_h_fraction_k(j+1)) &
312 * (t_h_fraction_k(j)-t_h_fraction_k(j+2)))
313
314 fact2 = ((te_h_table(i)-t_h_fraction_k(j)) &
315 * (te_h_table(i)-t_h_fraction_k(j+2))) &
316 / ((t_h_fraction_k(j+1)-t_h_fraction_k(j)) &
317 * (t_h_fraction_k(j+1)-t_h_fraction_k(j+2)))
318
319 fact3 = ((te_h_table(i)-t_h_fraction_k(j)) &
320 * (te_h_table(i)-t_h_fraction_k(j+1))) &
321 / ((t_h_fraction_k(j+2)-t_h_fraction_k(j)) &
322 * (t_h_fraction_k(j+2)-t_h_fraction_k(j+1)))
323
324 iz_h_table(i) = iz_h_source(j)*fact1 + &
325 iz_h_source(j+1)*fact2 + iz_h_source(j+2)*fact3
326 exit
327 end if
328 end if
329 end do
330 end do
331
332 te_h_table=te_h_table/unit_temperature
333 te_table_min=te_h_table(1)
334 te_table_max=te_h_table(n_interpolated_table)
335 inv_te_table_step = &
336 dble(n_interpolated_table-1)/(te_table_max-te_table_min)
337
338 call ionization_build_eos_table()
339 end subroutine ionization_build_temperature_table
340
341 subroutine ionization_build_prominence_table()
342 use mod_global_parameters, only: unit_temperature, unit_pressure, si_unit
343 use mod_comm_lib, only: mpistop
344
345 integer :: it, ip
346 double precision :: pressure_unit_cgs
347 double precision :: dt, rfactor
348
349 if (unit_temperature <= 0.d0 .or. unit_pressure <= 0.d0) then
350 call mpistop("invalid units for prominence ionization table")
351 end if
352
353 if (si_unit) then
354 ! unit_pressure is Pa; 1 Pa = 10 dyn cm^-2.
355 pressure_unit_cgs = 10.d0*unit_pressure
356 else
357 ! In cgs mode unit_pressure is already dyn cm^-2.
358 pressure_unit_cgs = unit_pressure
359 end if
360
361 allocate(prom_t_table(n_prom_t))
362 allocate(prom_p_table(n_prom_p))
363 allocate(prom_t_eos_table(n_prom_eos))
364 allocate(prom_y_eos_table(n_prom_eos,n_prom_p))
365
366 prom_t_table = prom_t_k/unit_temperature
367 prom_p_table = prom_p_cgs/pressure_unit_cgs
368
369 dt = (1.d6 - 2.1000005d3) / dble(n_prom_eos-1)
370 do it = 1, n_prom_eos
371 prom_t_eos_table(it) = &
372 (2.1000005d3 + dble(it-1)*dt)/unit_temperature
373 end do
374
375 do ip = 1, n_prom_p
376 do it = 1, n_prom_eos
377 call ionization_state_tp(prom_t_eos_table(it), prom_p_table(ip), rfactor)
378 prom_y_eos_table(it,ip) = &
379 rfactor*prom_t_eos_table(it)
380 end do
381 end do
382
383 call ionization_check_prominence_y_table()
384 end subroutine ionization_build_prominence_table
385
386 subroutine ionization_prominence_pressure_bracket(p, ip, weight)
387 double precision, intent(in) :: p
388 integer, intent(out) :: ip
389 double precision, intent(out) :: weight
390
391 integer :: j
392 double precision :: p_use
393
394 p_use = min(prom_p_table(n_prom_p), &
395 max(prom_p_table(1), p))
396
397 if (p_use <= prom_p_table(1)) then
398 ip = 1
399 weight = 0.d0
400 return
401 else if (p_use >= prom_p_table(n_prom_p)) then
402 ip = n_prom_p-1
403 weight = 1.d0
404 return
405 end if
406
407 do j = 1, n_prom_p-1
408 if (p_use <= prom_p_table(j+1)) then
409 ip = j
410 weight = (p_use-prom_p_table(j)) / &
411 (prom_p_table(j+1)-prom_p_table(j))
412 return
413 end if
414 end do
415 end subroutine ionization_prominence_pressure_bracket
416
417 !> Legacy array interface for T -> ionization degrees.
418 !> Not used anymore; kept for compatibility and may be removed later.
419 subroutine ionization_degree_from_temperature(ixI^L,ixO^L,Te,iz_H,iz_He)
420 use mod_global_parameters
421 integer, intent(in) :: ixi^l, ixo^l
422 double precision, intent(in) :: te(ixi^s)
423 double precision, intent(out) :: iz_h(ixo^s), iz_he(ixo^s)
424 integer :: ix^d
425 double precision :: rtmp
426
427 {do ix^db=ixomin^db,ixomax^db\}
428 call ionization_state_from_temperature_scalar( &
429 te(ix^d), rtmp, iz_h(ix^d), iz_he(ix^d))
430 {end do\}
431 end subroutine ionization_degree_from_temperature
432
433 subroutine ionization_h_degree_temperature_scalar(T, iz_H)
434 use mod_comm_lib, only: mpistop
435
436 double precision, intent(in) :: t
437 double precision, intent(out) :: iz_h
438
439 integer :: i
440
441 if (.not. allocated(te_h_table) .or. &
442 .not. allocated(iz_h_table)) then
443 call mpistop("temperature-only ionization table is not initialized")
444 end if
445
446 if (t < te_table_min) then
447 iz_h = 0.d0
448 else if (t >= te_table_max) then
449 iz_h = 1.d0
450 else
451 i = int((t-te_table_min)*inv_te_table_step) + 1
452 i = max(1, min(size(te_h_table)-1, i))
453
454 iz_h = iz_h_table(i) + &
455 (t-te_h_table(i)) * &
456 (iz_h_table(i+1)-iz_h_table(i)) * &
457 inv_te_table_step
458 end if
459
460 iz_h = min(1.d0, max(0.d0, iz_h))
461 end subroutine ionization_h_degree_temperature_scalar
462
463 subroutine ionization_h_degree_prominence_scalar(T, p, iz_H)
464 use mod_comm_lib, only: mpistop
465
466 double precision, intent(in) :: t, p
467 double precision, intent(out) :: iz_h
468
469 integer :: it, ip, j
470 double precision :: weight_t, weight_p
471 double precision :: iz_low, iz_high
472
473 if (.not. allocated(prom_t_table) .or. &
474 .not. allocated(prom_p_table)) then
475 call mpistop("prominence ionization table is not initialized")
476 end if
477
478 call ionization_prominence_pressure_bracket(p, ip, weight_p)
479
480 if (t <= prom_t_table(1)) then
481 it = 1
482 weight_t = 0.d0
483 else if (t >= prom_t_table(n_prom_t)) then
484 iz_h = 1.d0
485 return
486 else
487 it = 1
488 do j = 1, n_prom_t-1
489 if (t <= prom_t_table(j+1)) then
490 it = j
491 exit
492 end if
493 end do
494
495 weight_t = (t-prom_t_table(it)) / &
496 (prom_t_table(it+1)-prom_t_table(it))
497 end if
498
499 if (it == 1 .and. weight_t == 0.d0) then
500 iz_h = (1.d0-weight_p)*prom_iz_h_data(1,ip) + &
501 weight_p*prom_iz_h_data(1,ip+1)
502 else
503 iz_low = (1.d0-weight_p)*prom_iz_h_data(it,ip) + &
504 weight_p*prom_iz_h_data(it,ip+1)
505
506 iz_high = (1.d0-weight_p)*prom_iz_h_data(it+1,ip) + &
507 weight_p*prom_iz_h_data(it+1,ip+1)
508
509 iz_h = (1.d0-weight_t)*iz_low + weight_t*iz_high
510 end if
511
512 iz_h = min(1.d0, max(0.d0, iz_h))
513 end subroutine ionization_h_degree_prominence_scalar
514
515 subroutine ionization_get_state_scalar(rho, p, T, Rfactor, iz_H, iz_He)
516 use mod_comm_lib, only: mpistop
517 double precision, intent(in) :: rho, p
518 double precision, intent(out) :: t, rfactor
519 double precision, intent(out), optional :: iz_h, iz_he
520
521 double precision :: y
522 double precision :: iz_h_loc, iz_he_loc
523
524 if (.not. ionization_initialized) then
525 call mpistop("ionization tables are not initialized")
526 end if
527
528 if (rho <= 0.d0 .or. p <= 0.d0) then
529 select case (ionization_table_mode)
530 case (ionization_mode_chromosphere, ionization_mode_flare)
531 t = te_h_table(1)
532 case (ionization_mode_prominence)
533 t = prom_t_eos_table(1)
534 end select
535 else
536 y = p/rho
537
538 select case (ionization_table_mode)
539 case (ionization_mode_chromosphere, ionization_mode_flare)
540 call ionization_invert_y_scalar(y, t)
541
542 case (ionization_mode_prominence)
543 call ionization_invert_y_prominence_scalar(y, p, t)
544
545 case default
546 call mpistop("invalid ionization table mode")
547 end select
548 end if
549
550 call ionization_state_tp( t, p, rfactor, iz_h_loc, iz_he_loc)
551
552 if (present(iz_h)) iz_h = iz_h_loc
553 if (present(iz_he)) iz_he = iz_he_loc
554 end subroutine ionization_get_state_scalar
555
556 subroutine ionization_he_degree_scalar(T, iz_He)
557 use mod_global_parameters, only: unit_temperature
558
559 double precision, intent(in) :: t
560 double precision, intent(out) :: iz_he
561
562 if (t < te_low_iz_he) then
563 iz_he = 1.0084814d-4
564 else
565 iz_he = 1.d0 - &
566 10.d0**(0.322571d0-5.96d-5*t*unit_temperature)
567 end if
568
569 iz_he = min(1.d0, max(0.d0, iz_he))
570 end subroutine ionization_he_degree_scalar
571
572 subroutine ionization_state_tp(T, p, Rfactor, iz_H, iz_He)
573 use mod_comm_lib, only: mpistop
574
575 double precision, intent(in) :: t, p
576 double precision, intent(out) :: rfactor
577 double precision, intent(out), optional :: iz_h, iz_he
578
579 double precision :: iz_h_loc, iz_he_loc
580
581 select case (ionization_table_mode)
582 case (ionization_mode_chromosphere, ionization_mode_flare)
583 call ionization_h_degree_temperature_scalar(t, iz_h_loc)
584
585 case (ionization_mode_prominence)
586 call ionization_h_degree_prominence_scalar( &
587 t, p, iz_h_loc)
588
589 case default
590 call mpistop("invalid ionization table mode")
591 end select
592
593 call ionization_he_degree_scalar(t, iz_he_loc)
594
595 iz_h_loc = min(1.d0, max(0.d0, iz_h_loc))
596 iz_he_loc = min(1.d0, max(0.d0, iz_he_loc))
597
598 rfactor = ionization_rfactor_from_degrees( &
599 iz_h_loc, iz_he_loc)
600
601 if (present(iz_h)) iz_h = iz_h_loc
602 if (present(iz_he)) iz_he = iz_he_loc
603 end subroutine ionization_state_tp
604
605 subroutine ionization_solve_p_rfactor(rho, T, p, Rfactor)
606 use mod_comm_lib, only: mpistop
607
608 integer, parameter :: max_expand = 32, max_iter = 80
609 double precision, intent(in) :: rho, t
610 double precision, intent(out) :: p, rfactor
611
612 integer :: iter
613 double precision :: plo, phi, pmid
614 double precision :: rlo, rhi, rmid
615 double precision :: fhi, fmid, scale
616
617 if (.not. ionization_initialized) then
618 call mpistop("ionization pressure solver called before initialization")
619 end if
620 if (rho <= 0.d0 .or. t <= 0.d0) then
621 call mpistop("invalid rho or T in ionization pressure solver")
622 end if
623
624 if (ionization_is_temperature_only()) then
625 call ionization_state_from_temperature_scalar(t, rfactor)
626 p = rho*rfactor*t
627 return
628 end if
629
630 ! Low-pressure clamp branch
631 plo = prom_p_table(1)
632 call ionization_state_tp(t, plo, rlo)
633 p = rho*rlo*t
634 if (p <= plo) then
635 rfactor = rlo
636 return
637 end if
638
639 ! High-pressure clamp branch
640 phi = prom_p_table(n_prom_p)
641 call ionization_state_tp(t, phi, rhi)
642 p = rho*rhi*t
643 if (p >= phi) then
644 rfactor = rhi
645 return
646 end if
647
648 ! Now the root must lie inside [plo, phi]
649 do iter = 1, max_iter
650 pmid = 0.5d0*(plo+phi)
651 call ionization_state_tp(t, pmid, rmid)
652 fmid = pmid-rho*rmid*t
653 scale = max(abs(pmid), abs(rho*rmid*t), tiny(1.d0))
654
655 if (abs(fmid) <= 1.d-12*scale) then
656 p = pmid
657 rfactor = rmid
658 return
659 end if
660
661 if (fmid > 0.d0) then
662 phi = pmid
663 else
664 plo = pmid
665 end if
666 end do
667
668 call mpistop("pressure-dependent EOS solve did not converge")
669 end subroutine ionization_solve_p_rfactor
670
671 subroutine ionization_get_state(ixI^L, ixO^L, rho, p, T, Rfactor, iz_H, iz_He)
672 integer, intent(in) :: ixi^l, ixo^l
673 double precision, intent(in) :: rho(ixi^s), p(ixi^s)
674 double precision, intent(out) :: t(ixi^s), rfactor(ixi^s)
675 double precision, intent(out), optional :: iz_h(ixi^s), iz_he(ixi^s)
676
677 double precision :: iz_h_loc, iz_he_loc
678 integer :: ix^d
679
680 {do ix^db=ixomin^db,ixomax^db\}
681 call ionization_get_state_scalar(rho(ix^d), p(ix^d), t(ix^d), &
682 rfactor(ix^d), iz_h_loc, iz_he_loc)
683
684 if (present(iz_h)) iz_h(ix^d) = iz_h_loc
685 if (present(iz_he)) iz_he(ix^d) = iz_he_loc
686 {end do\}
687 end subroutine ionization_get_state
688
689 subroutine ionization_invert_y_scalar(y, T)
690 use mod_comm_lib, only: mpistop
691 double precision, intent(in) :: y
692 double precision, intent(out) :: t
693
694 integer :: ilo, ihi, imid, n
695 double precision :: denom
696
697 if (.not. allocated(y_table) .or. .not. allocated(rfactor_table)) then
698 call mpistop("ionization_invert_y_scalar: EOS tables are not initialized")
699 end if
700
701 n = size(y_table)
702
703 if (n < 2) then
704 call mpistop("ionization_invert_y_scalar: EOS table has fewer than 2 points")
705 end if
706
707 if (y <= 0.d0) then
708 t = te_h_table(1)
709 return
710 end if
711
712 if (y < y_table_min) then
713 if (rfactor_table(1) > 0.d0) then
714 t = max(tiny(1.d0), y/rfactor_table(1))
715 else
716 t = te_h_table(1)
717 end if
718 return
719 end if
720
721 if (y >= y_table_max) then
722 if (rfactor_table(n) > 0.d0) then
723 t = y/rfactor_table(n)
724 else
725 t = te_h_table(n)
726 end if
727 return
728 end if
729
730 ilo = 1
731 ihi = n
732
733 do while (ihi - ilo > 1)
734 imid = (ilo + ihi)/2
735 if (y >= y_table(imid)) then
736 ilo = imid
737 else
738 ihi = imid
739 end if
740 end do
741
742 denom = y_table(ilo+1) - y_table(ilo)
743
744 if (denom <= max(tiny(denom), 1.d-12*dabs(y_table(ilo)))) then
745 call mpistop("ionization_invert_y_scalar: non-increasing Y_table interval")
746 end if
747
748 t = te_h_table(ilo) + (y - y_table(ilo)) * &
749 (te_h_table(ilo+1) - te_h_table(ilo))/denom
750 end subroutine ionization_invert_y_scalar
751
752 subroutine ionization_state_from_temperature_scalar(T, Rfactor, iz_H, iz_He)
753 use mod_comm_lib, only: mpistop
754
755 double precision, intent(in) :: t
756 double precision, intent(out) :: rfactor
757 double precision, intent(out), optional :: iz_h, iz_he
758 double precision :: iz_h_loc, iz_he_loc
759
760 if (.not. ionization_is_temperature_only()) then
761 call mpistop("temperature-only ionization state requested "// &
762 "for a pressure-dependent table")
763 end if
764
765 select case (ionization_table_mode)
766 case (ionization_mode_chromosphere, ionization_mode_flare)
767 call ionization_h_degree_temperature_scalar(t, iz_h_loc)
768 end select
769 call ionization_he_degree_scalar(t, iz_he_loc)
770
771 rfactor = ionization_rfactor_from_degrees( &
772 iz_h_loc, iz_he_loc)
773
774 if (present(iz_h)) iz_h = iz_h_loc
775 if (present(iz_he)) iz_he = iz_he_loc
776 end subroutine ionization_state_from_temperature_scalar
777
778 subroutine ionization_check_prominence_y_table()
779 use mod_comm_lib, only: mpistop
780 integer :: it, ip
781
782 do ip = 1, n_prom_p
783 do it = 1, n_prom_eos-1
784 if (prom_y_eos_table(it+1,ip) <= &
785 prom_y_eos_table(it,ip)) then
786 call mpistop("prominence EOS Y(T,p) is not monotonic")
787 end if
788 end do
789 end do
790 end subroutine ionization_check_prominence_y_table
791
792 subroutine ionization_invert_y_prominence_scalar(y, p, T)
793 use mod_comm_lib, only: mpistop
794
795 double precision, intent(in) :: y, p
796 double precision, intent(out) :: t
797
798 integer :: ip, ilo, ihi, imid
799 double precision :: wp, ylo, yhi, ymid
800 double precision :: ymin, ymax, rend, denom
801
802 if (.not. allocated(prom_y_eos_table)) then
803 call mpistop("prominence EOS inversion table is missing")
804 end if
805
806 call ionization_prominence_pressure_bracket(p, ip, wp)
807
808 ymin = (1.d0-wp)*prom_y_eos_table(1,ip) + &
809 wp*prom_y_eos_table(1,ip+1)
810
811 ymax = (1.d0-wp)*prom_y_eos_table(n_prom_eos,ip) + &
812 wp*prom_y_eos_table(n_prom_eos,ip+1)
813
814 if (y <= 0.d0) then
815 t = prom_t_eos_table(1)
816 return
817 end if
818
819 if (y < ymin) then
820 rend = ymin/prom_t_eos_table(1)
821 t = max(tiny(1.d0), y/rend)
822 return
823 end if
824
825 if (y >= ymax) then
826 rend = ymax/prom_t_eos_table(n_prom_eos)
827 t = y/rend
828 return
829 end if
830
831 ilo = 1
832 ihi = n_prom_eos
833
834 do while (ihi-ilo > 1)
835 imid = (ilo+ihi)/2
836 ymid = (1.d0-wp)*prom_y_eos_table(imid,ip) + &
837 wp*prom_y_eos_table(imid,ip+1)
838
839 if (y >= ymid) then
840 ilo = imid
841 else
842 ihi = imid
843 end if
844 end do
845
846 ylo = (1.d0-wp)*prom_y_eos_table(ilo,ip) + &
847 wp*prom_y_eos_table(ilo,ip+1)
848
849 yhi = (1.d0-wp)*prom_y_eos_table(ilo+1,ip) + &
850 wp*prom_y_eos_table(ilo+1,ip+1)
851
852 denom = yhi-ylo
853 if (denom <= max(tiny(denom),1.d-12*abs(ylo))) then
854 call mpistop("invalid prominence EOS inversion interval")
855 end if
856
857 t = prom_t_eos_table(ilo) + (y-ylo) * &
858 (prom_t_eos_table(ilo+1)-prom_t_eos_table(ilo)) / &
859 denom
860 end subroutine ionization_invert_y_prominence_scalar
861
862 subroutine ionization_build_eos_table()
863 use mod_comm_lib, only: mpistop
864 integer :: i, n
865 double precision :: iz_h
866
867 if (.not. allocated(te_h_table) .or. .not. allocated(iz_h_table)) then
868 call mpistop("ionization_build_eos_table: ionization tables are not initialized")
869 end if
870
871 if (allocated(rfactor_table)) deallocate(rfactor_table)
872 if (allocated(y_table)) deallocate(y_table)
873 if (allocated(eps_ion_h_table)) deallocate(eps_ion_h_table)
874
875 n = size(te_h_table)
876
877 if (n < 2) then
878 call mpistop("ionization_build_eos_table: Te_H_table has fewer than 2 points")
879 end if
880
881 allocate(rfactor_table(1:n))
882 allocate(y_table(1:n))
883 if (ionization_include_energy) allocate(eps_ion_h_table(1:n))
884
885 do i = 1, n
886 call ionization_state_from_temperature_scalar( &
887 te_h_table(i), rfactor_table(i), iz_h)
888 y_table(i) = rfactor_table(i)*te_h_table(i)
889 if (ionization_include_energy) then
890 eps_ion_h_table(i) = ionization_h_energy_unit*iz_h
891 end if
892 end do
893
894 y_table_min = y_table(1)
895 y_table_max = y_table(n)
896
897 call ionization_check_y_table()
898 end subroutine ionization_build_eos_table
899
900 subroutine ionization_check_y_table()
901 use mod_comm_lib, only: mpistop
902 integer :: i
903
904 if (.not. allocated(y_table)) then
905 call mpistop("ionization_check_y_table: Y_table is not allocated")
906 end if
907
908 if (size(y_table) < 2) then
909 call mpistop("ionization_check_y_table: Y_table has fewer than 2 points")
910 end if
911
912 do i = 1, size(y_table)-1
913 if (y_table(i+1) <= y_table(i)) then
914 call mpistop("ionization_check_y_table: Y(T)=Rfactor(T)*T is not strictly increasing")
915 end if
916 end do
917 end subroutine ionization_check_y_table
918
919 double precision function ionization_rfactor_from_degrees(iz_H, iz_He)
920 double precision, intent(in) :: iz_h, iz_he
921
922 ionization_rfactor_from_degrees = &
923 (1.d0 + iz_h + ionization_he_abundance * &
924 (1.d0 + iz_he*(1.d0 + iz_he))) / ionization_rfactor_norm
925 end function ionization_rfactor_from_degrees
926
927 subroutine ionization_get_rfactor_from_temperature(T, Rfactor)
928 double precision, intent(in) :: t
929 double precision, intent(out) :: rfactor
930
931 call ionization_state_from_temperature_scalar(t, rfactor)
932 end subroutine ionization_get_rfactor_from_temperature
933
934 logical function ionization_is_temperature_only()
935 ionization_is_temperature_only = &
936 ionization_table_mode == ionization_mode_chromosphere .or. &
937 ionization_table_mode == ionization_mode_flare
938 end function ionization_is_temperature_only
939
940 logical function ionization_includes_energy()
941 ionization_includes_energy = ionization_include_energy
942 end function ionization_includes_energy
943
944 double precision function ionization_eint_table_value(i, invgam)
945 integer, intent(in) :: i
946 double precision, intent(in) :: invgam
947
948 ionization_eint_table_value = &
949 y_table(i)*invgam + eps_ion_h_table(i)
950 end function ionization_eint_table_value
951
952 subroutine ionization_check_eint_table(invgam)
953 use mod_comm_lib, only: mpistop
954 double precision, intent(in) :: invgam
955 integer :: i
956 double precision :: eps1, eps2
957
958 if (.not. ionization_include_energy) return
959 if (.not. allocated(eps_ion_h_table)) then
960 call mpistop("ionization-energy table is not initialized")
961 end if
962 if (invgam <= 0.d0) then
963 call mpistop("invalid inverse gamma minus one for ionization EOS")
964 end if
965
966 do i = 1, size(y_table)-1
967 eps1 = ionization_eint_table_value(i, invgam)
968 eps2 = ionization_eint_table_value(i+1, invgam)
969 if (eps2 <= eps1) then
970 call mpistop("ionization specific-energy table is not monotonic")
971 end if
972 end do
973 end subroutine ionization_check_eint_table
974
975 subroutine ionization_invert_eint_scalar(eps, invgam, T)
976 use mod_comm_lib, only: mpistop
977 double precision, intent(in) :: eps, invgam
978 double precision, intent(out) :: t
979
980 integer :: ilo, ihi, imid, n, iter
981 double precision :: eps_min, eps_max, eps_lo, eps_hi, denom
982 double precision :: eps_eval, slope, rfactor, iz_h, tnew
983
984 if (.not. ionization_include_energy .or. &
985 .not. allocated(eps_ion_h_table)) then
986 call mpistop("ionization-energy EOS is not initialized")
987 end if
988 if (invgam <= 0.d0) then
989 call mpistop("invalid inverse gamma minus one in energy inversion")
990 end if
991
992 n = size(y_table)
993 eps_min = ionization_eint_table_value(1, invgam)
994 eps_max = ionization_eint_table_value(n, invgam)
995
996 if (eps < eps_min) then
997 ! Below the tabulated eion energy range this is floor/illegal-state
998 ! recovery, not a strictly conservative EOS inversion.
999 t = max(tiny(1.d0), &
1000 (eps-eps_ion_h_table(1))/(rfactor_table(1)*invgam))
1001 return
1002 end if
1003
1004 if (eps >= eps_max) then
1005 t = max(tiny(1.d0), &
1006 (eps-eps_ion_h_table(n))/(rfactor_table(n)*invgam))
1007 return
1008 end if
1009
1010 ilo = 1
1011 ihi = n
1012 do while (ihi-ilo > 1)
1013 imid = (ilo+ihi)/2
1014 if (eps >= ionization_eint_table_value(imid, invgam)) then
1015 ilo = imid
1016 else
1017 ihi = imid
1018 end if
1019 end do
1020
1021 eps_lo = ionization_eint_table_value(ilo, invgam)
1022 eps_hi = ionization_eint_table_value(ilo+1, invgam)
1023 denom = eps_hi-eps_lo
1024 if (denom <= max(tiny(denom), 1.d-12*abs(eps_lo))) then
1025 call mpistop("invalid ionization-energy inversion interval")
1026 end if
1027
1028 t = te_h_table(ilo) + (eps-eps_lo) * &
1029 (te_h_table(ilo+1)-te_h_table(ilo))/denom
1030
1031 ! The tabulated value supplies a close initial guess. Correct it against
1032 ! the actual interpolated ionization state so rho,eint -> T and
1033 ! rho,T -> eint remain mutually consistent inside each table interval.
1034 slope = denom/(te_h_table(ilo+1)-te_h_table(ilo))
1035 do iter = 1, 3
1036 call ionization_state_from_temperature_scalar( &
1037 t, rfactor, iz_h)
1038 eps_eval = rfactor*t*invgam + ionization_h_energy_unit*iz_h
1039 if (abs(eps_eval-eps) <= &
1040 1.d-12*max(abs(eps), tiny(1.d0))) exit
1041 tnew = t-(eps_eval-eps)/slope
1042 t = max(te_h_table(ilo), min(te_h_table(ilo+1), tnew))
1043 end do
1044 end subroutine ionization_invert_eint_scalar
1045
1046 subroutine ionization_get_state_from_eint( &
1047 rho, eint, invgam, T, p, Rfactor, iz_H, iz_He)
1048 use mod_comm_lib, only: mpistop
1049 double precision, intent(in) :: rho, eint, invgam
1050 double precision, intent(out) :: t, p, rfactor
1051 double precision, intent(out), optional :: iz_h, iz_he
1052 double precision :: iz_h_loc, iz_he_loc
1053
1054 if (.not. ionization_include_energy) &
1055 call mpistop("ionization-energy EOS is not enabled")
1056 if (.not. ionization_is_temperature_only()) &
1057 call mpistop("ionization-energy EOS requires a T-only table")
1058 if (rho <= 0.d0 .or. eint <= 0.d0) &
1059 call mpistop("invalid rho or eint in ionization EOS")
1060 call ionization_invert_eint_scalar(eint/rho, invgam, t)
1061 call ionization_state_from_temperature_scalar( &
1062 t, rfactor, iz_h_loc, iz_he_loc)
1063 p = rho*rfactor*t
1064 if (present(iz_h)) iz_h = iz_h_loc
1065 if (present(iz_he)) iz_he = iz_he_loc
1066 end subroutine ionization_get_state_from_eint
1067
1068 subroutine ionization_get_p_eint_from_rho_t( &
1069 rho, T, invgam, p, eint, Rfactor, iz_H, iz_He)
1070 use mod_comm_lib, only: mpistop
1071 double precision, intent(in) :: rho, t, invgam
1072 double precision, intent(out) :: p, eint, rfactor
1073 double precision, intent(out), optional :: iz_h, iz_he
1074 double precision :: iz_h_loc, iz_he_loc
1075
1076 if (rho <= 0.d0 .or. t <= 0.d0 .or. invgam <= 0.d0) then
1077 call mpistop("invalid state in rho,T ionization EOS mapping")
1078 end if
1079
1080 if (ionization_include_energy) then
1081 if (.not. ionization_is_temperature_only()) then
1082 call mpistop("ionization energy requires a T-only ionization table")
1083 end if
1084 call ionization_state_from_temperature_scalar( &
1085 t, rfactor, iz_h_loc, iz_he_loc)
1086 p = rho*rfactor*t
1087 ! Only hydrogen ionization energy is included in this first version;
1088 ! helium ionization energy is intentionally deferred.
1089 eint = rho*(rfactor*t*invgam + &
1090 ionization_h_energy_unit*iz_h_loc)
1091 else
1092 call ionization_solve_p_rfactor(rho, t, p, rfactor)
1093 call ionization_state_tp(t, p, rfactor, iz_h_loc, iz_he_loc)
1094 eint = p*invgam
1095 end if
1096
1097 if (present(iz_h)) iz_h = iz_h_loc
1098 if (present(iz_he)) iz_he = iz_he_loc
1099 end subroutine ionization_get_p_eint_from_rho_t
1100
1101 subroutine ionization_get_eps_derivative_t( &
1102 T, invgam, eps, deps_dT, dq_dT_out)
1103 use mod_comm_lib, only: mpistop
1104 double precision, intent(in) :: t, invgam
1105 double precision, intent(out) :: eps, deps_dt
1106 double precision, intent(out), optional :: dq_dt_out
1107
1108 integer :: i, n
1109 double precision :: rfactor, iz_h, q, dq_dt, depsion_dt
1110
1111 if (.not. ionization_include_energy) then
1112 call mpistop("ionization-energy derivative requested while disabled")
1113 end if
1114 if (.not. ionization_is_temperature_only()) then
1115 call mpistop("ionization-energy derivative requires a T-only table")
1116 end if
1117 if (t <= 0.d0 .or. invgam <= 0.d0) then
1118 call mpistop("invalid state in ionization-energy derivative")
1119 end if
1120
1121 n = size(te_h_table)
1122 call ionization_state_from_temperature_scalar(t, rfactor, iz_h)
1123 q = rfactor*t
1124 eps = q*invgam + ionization_h_energy_unit*iz_h
1125
1126 if (t < te_table_min) then
1127 dq_dt = rfactor_table(1)
1128 depsion_dt = 0.d0
1129 else if (t >= te_table_max) then
1130 dq_dt = rfactor_table(n)
1131 depsion_dt = 0.d0
1132 else
1133 i = int((t-te_table_min)*inv_te_table_step) + 1
1134 i = max(1, min(n-1, i))
1135 dq_dt = (y_table(i+1)-y_table(i)) / &
1136 (te_h_table(i+1)-te_h_table(i))
1137 depsion_dt = (eps_ion_h_table(i+1)-eps_ion_h_table(i)) / &
1138 (te_h_table(i+1)-te_h_table(i))
1139 end if
1140
1141 deps_dt = dq_dt*invgam + depsion_dt
1142 if (deps_dt <= 0.d0) then
1143 call mpistop("non-positive heat capacity in ionization EOS")
1144 end if
1145 if (present(dq_dt_out)) dq_dt_out = dq_dt
1146 end subroutine ionization_get_eps_derivative_t
1147
1148 subroutine ionization_get_csound2_t(T, invgam, csound2)
1149 double precision, intent(in) :: t, invgam
1150 double precision, intent(out) :: csound2
1151
1152 double precision :: rfactor, q, eps, deps_dt, dq_dt
1153
1154 call ionization_get_eps_derivative_t( &
1155 t, invgam, eps, deps_dt, dq_dt)
1156 call ionization_state_from_temperature_scalar(t, rfactor)
1157 q = rfactor*t
1158 csound2 = q*(1.d0+dq_dt/deps_dt)
1159 end subroutine ionization_get_csound2_t
1160
1161end module mod_ionization_degree
mod_comm_lib
Definition mod_comm_lib.t:1
mod_comm_lib::mpistop
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition mod_comm_lib.t:208
mod_constants
Module for physical and numeric constants.
Definition mod_constants.t:2
mod_constants::const_ev
double precision, parameter const_ev
Definition mod_constants.t:72
mod_global_parameters
This module contains definitions of global parameters and variables and some generic functions/subrou...
Definition mod_global_parameters.t:4
mod_global_parameters::it
integer it
Number of time steps taken.
Definition mod_global_parameters.t:444
mod_global_parameters::unit_numberdensity
double precision unit_numberdensity
Physical scaling factor for number density.
Definition mod_global_parameters.t:93
mod_global_parameters::unit_pressure
double precision unit_pressure
Physical scaling factor for pressure.
Definition mod_global_parameters.t:87
mod_global_parameters::d
double precision, dimension(:), allocatable, parameter d
Definition mod_global_parameters.t:23
mod_global_parameters::dt
double precision dt
global time step
Definition mod_global_parameters.t:32
mod_global_parameters::l
double precision l
Definition mod_global_parameters.t:16
mod_global_parameters::unit_temperature
double precision unit_temperature
Physical scaling factor for temperature.
Definition mod_global_parameters.t:84
mod_global_parameters::si_unit
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
Definition mod_global_parameters.t:666
mod_ionization_degree
module ionization degree - get ionization degree for given temperature
Definition mod_ionization_degree.t:2
mod_ionization_degree::ionization_get_state
subroutine, public ionization_get_state(ixil, ixol, rho, p, t, rfactor, iz_h, iz_he)
Definition mod_ionization_degree.t:672
mod_ionization_degree::ionization_get_csound2_t
subroutine, public ionization_get_csound2_t(t, invgam, csound2)
Definition mod_ionization_degree.t:1149
mod_ionization_degree::ionization_get_state_from_eint
subroutine, public ionization_get_state_from_eint(rho, eint, invgam, t, p, rfactor, iz_h, iz_he)
Definition mod_ionization_degree.t:1048
mod_ionization_degree::ionization_solve_p_rfactor
subroutine, public ionization_solve_p_rfactor(rho, t, p, rfactor)
Definition mod_ionization_degree.t:606
mod_ionization_degree::ionization_get_eps_derivative_t
subroutine, public ionization_get_eps_derivative_t(t, invgam, eps, deps_dt, dq_dt_out)
Definition mod_ionization_degree.t:1103
mod_ionization_degree::ionization_is_temperature_only
logical function, public ionization_is_temperature_only()
Definition mod_ionization_degree.t:935
mod_ionization_degree::ionization_check_eint_table
subroutine, public ionization_check_eint_table(invgam)
Definition mod_ionization_degree.t:953
mod_ionization_degree::ionization_get_p_eint_from_rho_t
subroutine, public ionization_get_p_eint_from_rho_t(rho, t, invgam, p, eint, rfactor, iz_h, iz_he)
Definition mod_ionization_degree.t:1070
mod_ionization_degree::ionization_get_state_scalar
subroutine, public ionization_get_state_scalar(rho, p, t, rfactor, iz_h, iz_he)
Definition mod_ionization_degree.t:516
mod_ionization_degree::ionization_state_tp
subroutine, public ionization_state_tp(t, p, rfactor, iz_h, iz_he)
Definition mod_ionization_degree.t:573
mod_ionization_degree::ionization_degree_from_temperature
subroutine, public ionization_degree_from_temperature(ixil, ixol, te, iz_h, iz_he)
Legacy array interface for T -> ionization degrees. Not used anymore; kept for compatibility and may ...
Definition mod_ionization_degree.t:420
mod_ionization_degree::ionization_state_from_temperature_scalar
subroutine, public ionization_state_from_temperature_scalar(t, rfactor, iz_h, iz_he)
Definition mod_ionization_degree.t:753
mod_ionization_degree::ionization_includes_energy
logical function, public ionization_includes_energy()
Definition mod_ionization_degree.t:941
mod_ionization_degree::ionization_degree_init
subroutine, public ionization_degree_init(he_abundance, rfactor_norm, table_name, include_energy)
Definition mod_ionization_degree.t:173
mod_ionization_degree::ionization_get_rfactor_from_temperature
subroutine, public ionization_get_rfactor_from_temperature(t, rfactor)
Definition mod_ionization_degree.t:928
Generated by doxygen 1.9.8