MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_ionization_degree.t
Go to the documentation of this file.
1 !> module ionization degree - get ionization degree for given temperature
3  implicit none
4  double precision, private :: Te_H_CL(1:100)
5  double precision, private :: iz_H_CL(1:100)
6  double precision, dimension(:), allocatable :: te_h_table
7  double precision, dimension(:), allocatable :: iz_h_table
9  ! Carlsson, M., & Leenaarts, J. 2012, A&A, 539, A39
10  data te_h_cl /2.1000005d+03, 2.2349495d+03, 2.3785708d+03, 2.5314199d+03, 2.6940928d+03, &
11  2.8672190d+03, 3.0514692d+03, 3.2475613d+03, 3.4562544d+03, 3.6783584d+03, &
12  3.9147351d+03, 4.1662998d+03, 4.4340322d+03, 4.7189697d+03, 5.0222153d+03, &
13  5.3449502d+03, 5.6884248d+03, 6.0539712d+03, 6.4430083d+03, 6.8570420d+03, &
14  7.2976860d+03, 7.7666460d+03, 8.2657383d+03, 8.7969062d+03, 9.3622090d+03, &
15  9.9638330d+03, 1.0604130d+04, 1.1285561d+04, 1.2010781d+04, 1.2782619d+04, &
16  1.3604041d+04, 1.4478250d+04, 1.5408652d+04, 1.6398826d+04, 1.7452648d+04, &
17  1.8574172d+04, 1.9767766d+04, 2.1038082d+04, 2.2390010d+04, 2.3828838d+04, &
18  2.5360100d+04, 2.6989764d+04, 2.8724182d+04, 3.0570023d+04, 3.2534518d+04, &
19  3.4625215d+04, 3.6850262d+04, 3.9218293d+04, 4.1738543d+04, 4.4420699d+04, &
20  4.7275266d+04, 5.0313219d+04, 5.3546391d+04, 5.6987395d+04, 6.0649453d+04, &
21  6.4546914d+04, 6.8694758d+04, 7.3109148d+04, 7.7807289d+04, 8.2807258d+04, &
22  8.8128625d+04, 9.3791852d+04, 9.9819000d+04, 1.0623346d+05, 1.1306024d+05, &
23  1.2032560d+05, 1.2805797d+05, 1.3628709d+05, 1.4504503d+05, 1.5436592d+05, &
24  1.6428561d+05, 1.7484295d+05, 1.8607852d+05, 1.9803609d+05, 2.1076208d+05, &
25  2.2430608d+05, 2.3872045d+05, 2.5406084d+05, 2.7038703d+05, 2.8776234d+05, &
26  3.0625456d+05, 3.2593475d+05, 3.4688000d+05, 3.6917081d+05, 3.9289406d+05, &
27  4.1814181d+05, 4.4501247d+05, 4.7360988d+05, 5.0404450d+05, 5.3643488d+05, &
28  5.7090725d+05, 6.0759431d+05, 6.4663956d+05, 6.8819319d+05, 7.3241712d+05, &
29  7.7948294d+05, 8.2957412d+05, 8.8288331d+05, 9.3961919d+05, 1.0000000d+06 /
30 
31  ! Carlsson, M., & Leenaarts, J. 2012, A&A, 539, A39
32  data iz_h_cl /1.0000000d+00, 1.0000000d+00, 1.0000000d+00, 1.0000000d+00, 1.0000000d+00, &
33  1.0000000d+00, 1.0000000d+00, 9.9999928d-01, 9.9999774d-01, 9.9999624d-01, &
34  9.9999428d-01, 9.9999219d-01, 9.9998862d-01, 9.9998313d-01, 9.9997944d-01, &
35  9.9995774d-01, 1.0000324d+00, 9.9822283d-01, 9.6391600d-01, 8.8995951d-01, &
36  8.0025935d-01, 7.4219799d-01, 7.1378660d-01, 6.9565988d-01, 6.7326778d-01, &
37  6.5032905d-01, 6.2487972d-01, 5.9925246d-01, 5.7028764d-01, 5.3673893d-01, &
38  4.9144468d-01, 4.4332290d-01, 3.9158964d-01, 3.4621361d-01, 2.7826238d-01, &
39  1.9633104d-01, 1.2593637d-01, 8.1383914d-02, 5.0944358d-02, 2.8770180d-02, &
40  1.4449068d-02, 7.6071853d-03, 4.3617180d-03, 2.7694483d-03, 1.8341533d-03, &
41  1.2467797d-03, 8.6159195d-04, 6.0232455d-04, 4.2862241d-04, 3.0666622d-04, &
42  2.2256430d-04, 1.6328457d-04, 1.2179965d-04, 9.2065704d-05, 7.0357783d-05, &
43  5.4168402d-05, 4.2026968d-05, 3.2807184d-05, 2.5801779d-05, 2.0414085d-05, &
44  1.6271379d-05, 1.3041737d-05, 1.0538992d-05, 8.5329248d-06, 6.9645471d-06, &
45  5.7069005d-06, 4.7155136d-06, 3.9203474d-06, 3.2734958d-06, 2.7415674d-06, &
46  2.2911979d-06, 1.9045801d-06, 1.5692771d-06, 1.2782705d-06, 1.0281080d-06, &
47  8.1567373d-07, 6.4200310d-07, 5.0668103d-07, 4.0577518d-07, 3.3245340d-07, &
48  2.8006056d-07, 2.4119515d-07, 2.1033340d-07, 1.8429903d-07, 1.6189058d-07, &
49  1.4264889d-07, 1.2616209d-07, 1.1193421d-07, 9.9408616d-08, 8.8335305d-08, &
50  7.8593857d-08, 7.0042269d-08, 6.2532941d-08, 5.5898361d-08, 5.0049657d-08, &
51  4.4888207d-08, 4.0317172d-08, 3.6192560d-08, 3.2330732d-08, 2.9387850d-08 /
52  double precision :: te_low_iz_he=5413.d0
53  contains
54 
57 
58  double precision :: fact1, fact2, fact3, dL1, dL2
59  integer :: i, j, ntable, n_interpolated_table
60  logical :: jump
61 
62  ntable=100
63  n_interpolated_table=4000
64 
65  ! change neutral fraction to ionization degree
66  iz_h_cl=1.d0-iz_h_cl
67 
68  allocate(te_h_table(1:n_interpolated_table))
69  allocate(iz_h_table(1:n_interpolated_table))
70  te_h_table=zero
71  iz_h_table=zero
72 
73  ! create cooling table(s) for use in amrvac
74  te_table_step = (te_h_cl(ntable)-te_h_cl(1))/dble(n_interpolated_table-1)
75 
76  te_h_table(1) = te_h_cl(1)
77  iz_h_table(1) = iz_h_cl(1)
78 
79  te_h_table(n_interpolated_table) = te_h_cl(ntable)
80  iz_h_table(n_interpolated_table) = iz_h_cl(ntable)
81 
82  do i=2,n_interpolated_table ! loop to create one table
84  do j=1,ntable-1 ! loop to create one spot on a table
85  ! Second order polynomial interpolation, except at the outer edge,
86  ! or in case of a large jump.
87  if(te_h_table(i) < te_h_cl(j+1)) then
88  if(j.eq. ntable-1 )then
89  fact1 = (te_h_table(i)-te_h_cl(j+1)) &
90  /(te_h_cl(j)-te_h_cl(j+1))
91 
92  fact2 = (te_h_table(i)-te_h_cl(j)) &
93  /(te_h_cl(j+1)-te_h_cl(j))
94 
95  iz_h_table(i) = iz_h_cl(j)*fact1 + iz_h_cl(j+1)*fact2
96  exit
97  else
98  dl1 = iz_h_cl(j+1)-iz_h_cl(j)
99  dl2 = iz_h_cl(j+2)-iz_h_cl(j+1)
100  jump =(max(dabs(dl1),dabs(dl2)) > 2.d0*min(dabs(dl1),dabs(dl2)))
101  end if
102 
103  if( jump ) then
104  fact1 = (te_h_table(i)-te_h_cl(j+1)) &
105  /(te_h_cl(j)-te_h_cl(j+1))
106 
107  fact2 = (te_h_table(i)-te_h_cl(j)) &
108  /(te_h_cl(j+1)-te_h_cl(j))
109 
110  iz_h_table(i) = iz_h_cl(j)*fact1 + iz_h_cl(j+1)*fact2
111  exit
112  else
113  fact1 = ((te_h_table(i)-te_h_cl(j+1)) &
114  * (te_h_table(i)-te_h_cl(j+2))) &
115  / ((te_h_cl(j)-te_h_cl(j+1)) &
116  * (te_h_cl(j)-te_h_cl(j+2)))
117 
118  fact2 = ((te_h_table(i)-te_h_cl(j)) &
119  * (te_h_table(i)-te_h_cl(j+2))) &
120  / ((te_h_cl(j+1)-te_h_cl(j)) &
121  * (te_h_cl(j+1)-te_h_cl(j+2)))
122 
123  fact3 = ((te_h_table(i)-te_h_cl(j)) &
124  * (te_h_table(i)-te_h_cl(j+1))) &
125  / ((te_h_cl(j+2)-te_h_cl(j)) &
126  * (te_h_cl(j+2)-te_h_cl(j+1)))
127 
128  iz_h_table(i) = iz_h_cl(j)*fact1 + iz_h_cl(j+1)*fact2 &
129  + iz_h_cl(j+2)*fact3
130  exit
131  endif
132  endif
133  enddo ! end loop to find create one spot on a table
134  enddo ! end loop to create one table
135 
136  ! nondimensionalize temperature
139  te_table_max=te_h_table(n_interpolated_table)
141 
142  ! transition temperature
144 
145  end subroutine ionization_degree_init
146 
147  subroutine ionization_degree_from_temperature(ixI^L,ixO^L,Te,iz_H,iz_He)
149  integer, intent(in) :: ixI^L, ixO^L
150  double precision, intent(in) :: Te(ixI^S)
151  double precision, intent(out) :: iz_H(ixO^S),iz_He(ixO^S)
152 
153  integer :: ix^D, i
154 
155  {do ix^db=ixomin^db,ixomax^db\}
156  if(te(ix^d)<te_table_min) then
157  iz_h(ix^d)=0.d0
158  else if (te(ix^d)>te_table_max) then
159  iz_h(ix^d)=1.d0
160  else
161  i=int((te(ix^d)-te_table_min)*inv_te_table_step)+1
162  iz_h(ix^d)=iz_h_table(i)+(te(ix^d)-te_h_table(i))&
164  end if
165  ! Ni, L. et al. A&A, 665, A116
166  if(te(ix^d)<te_low_iz_he) then
167  iz_he(ix^d)=1.0084814d-4
168  else
169  iz_he(ix^d)=1.d0-10**(0.322571d0-5.96d-5*te(ix^d)*unit_temperature)
170  end if
171  {end do\}
172 
174 
175  ! gas constant R in ideal gas law for solar plasma
178  double precision, intent(in) :: te
179  double precision :: r_ideal_gas_law_partial_ionization
180 
181  double precision :: iz_h, iz_he
182  integer :: i
183 
184  if(te<te_table_min) then
185  iz_h=0.d0
186  else if (te>te_table_max) then
187  iz_h=1.d0
188  else
189  i=int((te-te_table_min)*inv_te_table_step)+1
190  iz_h=iz_h_table(i)+(te-te_h_table(i))&
192  end if
193  if(te<te_low_iz_he) then
194  iz_he=1.0084814d-4
195  else
196  iz_he=1.d0-10**(0.322571d0-5.96d-5*te*unit_temperature)
197  end if
198  ! dimensionless: kB and mp are included in units
199  ! assume the first and second ionization of Helium have the same degree
200  r_ideal_gas_law_partial_ionization=(1.d0+iz_h+0.1d0*(1.d0+iz_he*(1.d0+iz_he)))/2.3d0
201  return
202 
204 
205 end module mod_ionization_degree
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer, dimension(:), allocatable, parameter d
double precision unit_temperature
Physical scaling factor for temperature.
module ionization degree - get ionization degree for given temperature
double precision, dimension(:), allocatable te_h_table
double precision te_table_max
double precision, dimension(:), allocatable iz_h_table
double precision function r_ideal_gas_law_partial_ionization(Te)
subroutine ionization_degree_from_temperature(ixIL, ixOL, Te, iz_H, iz_He)
double precision inv_te_table_step
double precision te_table_min
double precision te_low_iz_he
double precision te_table_step