-
Notifications
You must be signed in to change notification settings - Fork 37
Expand file tree
/
Copy pathfang.f90
More file actions
132 lines (96 loc) · 3.81 KB
/
Copy pathfang.f90
File metadata and controls
132 lines (96 loc) · 3.81 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
module ionize_fang
use, intrinsic :: iso_fortran_env, only: sp=>real32, dp=>real64
use phys_consts, only: wp, kb
implicit none
private
real(wp), parameter :: deps = 0.035_wp
real(wp), parameter :: erg2kev = 624150648._wp
!! keV, kinetic energy lost per ion-electron pair produced
public :: fang2008, fang2010, gravity_accel, erg2kev
contains
elemental real(wp) function fang2010(Q0_keV, Emono_keV, Tn, massden_gcm3, meanmass_g, g_ms2) result(Qtot)
real(wp), intent(in) :: Q0_keV, Emono_keV, Tn, massden_gcm3, meanmass_g, g_ms2
real(wp) :: y, H_cm, f
integer :: i, j
real(wp), dimension(8) :: C
real(wp), parameter :: P(8,4) = reshape( &
[1.24616_wp, 1.45903_wp, -2.42269e-1_wp, 5.95459e-2_wp, &
2.23976_wp, -4.22918e-7_wp, 1.36458e-2_wp, 2.53332e-3_wp, &
1.41754_wp, 1.44597e-1_wp, 1.70433e-2_wp, 6.39717e-4_wp, &
2.48775e-1_wp, -1.50890e-1_wp, 6.30894e-9_wp, 1.23707e-3_wp, &
-4.65119e-1_wp, -1.05081e-1_wp, -8.95701e-2_wp, 1.22450e-2_wp, &
3.86019e-1_wp, 1.75430e-3_wp, -7.42960e-4_wp, 4.60881e-4_wp, &
-6.45454e-1_wp, 8.49555e-4_wp, -4.28581e-2_wp, -2.99302e-3_wp, &
9.48930e-1_wp, 1.97385e-1_wp, -2.50660e-3_wp, -2.06938e-3_wp], shape(P), order=[2,1])
!! scale height
!! Equation (2)
H_cm = 100 * kb * Tn / (meanmass_g/1000) / abs(g_ms2)
!! normalized atmospheric column mass
!! Equation (1)
y = 2/Emono_keV * (massden_gcm3 * H_cm / 6e-6_wp)**0.7_wp
!! Equation (5)
C = 0
do i=1,size(P,1)
do j=1,size(P,2)
C(i) = C(i) + P(i,j) * log(Emono_keV)**(j-1)
end do
end do
C = exp(C)
!! Equation (4)
!! Energy deposition "f"
f = C(1)*y**C(2)*exp(-1*C(3)*y**C(4)) + C(5)*y**C(6) * exp(-1*C(7)*y**C(8))
!! Equation (3)
!! total ionization rate "qtot" [cm^-3 s^-1]
Qtot = f * Q0_keV / deps / H_cm
end function fang2010
elemental real(wp) function fang2008(Q0_keV, E0_keV, Tn, massden_gcm3, meanmass_g, g_ms2) result(Qtot)
!! COMPUTE IONIZATION RATES PER THE FANG 2008 SEMI-EMPIRICAL METHOD.
!! https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2008JA013384
!! Total Ionization Rate by Precipitating Electrons With a
!! Maxwellian Energy and Isotropic Pitch Angle Distribution
real(wp), intent(in) :: Q0_keV, E0_keV, Tn, massden_gcm3, meanmass_g, g_ms2
!! Q0: [keV cm^-2 s^-1]
!! massden_gcm3: [g cm^-3]
!! meanmass_g: [g]
!! g_ms2: [m s^-2]
real(wp) :: y, H_cm, f
integer :: i, j
real(wp), dimension(8) :: C
real(wp), parameter :: P(8,4) = reshape( &
[3.49979e-1_wp, -6.18200e-2_wp, -4.08124e-2_wp, 1.65414e-2_wp, &
5.85425e-1_wp, -5.00793e-2_wp, 5.69309e-2_wp, -4.02491e-3_wp, &
1.69692e-1_wp, -2.58981e-2_wp, 1.96822e-2_wp, 1.20505e-3_wp, &
-1.22271e-1_wp, -1.15532e-2_wp, 5.37951e-6_wp, 1.20189e-3_wp, &
1.57018_wp, 2.87896e-1_wp, -4.14857e-1_wp, 5.18158e-2_wp, &
8.83195e-1_wp, 4.31402e-2_wp, -8.33599e-2_wp, 1.02515e-2_wp, &
1.90953_wp, -4.74704e-2_wp, -1.80200e-1_wp, 2.46652e-2_wp, &
-1.29566_wp, -2.10952e-1_wp, 2.73106e-1_wp, -2.92752e-2_wp], shape(P), order=[2,1])
!! Equation (3)
!! scale height
!! kb [m^2 kg s^-2 K^-1]
H_cm = 100 * kb * Tn / (meanmass_g/1000) / abs(g_ms2)
!! H_cm [centimeters]
!! Equation (4)
y = 1 / E0_keV*(massden_gcm3 * H_cm / 4e-6_wp)**(0.606_wp)
!! Ci COEFFS and SHAPE FUNCTION
!! Equation (7)
C = 0
do i=1,size(P,1)
do j=1,size(P,2)
C(i) = C(i) + P(i,j) * log(E0_keV)**(j-1)
end do
end do
C = exp(C)
!! Equation (6) energy deposition "f"
f = C(1)*y**C(2)*exp(-1*C(3)*y**C(4))+C(5)*y**C(6)*exp(-1*C(7)*y**C(8))
!! Equation (2) total electron impact ionization rate
Qtot = Q0_keV / 2._wp / deps / H_cm * f
!! [cm^-3 s^-1]
end function fang2008
elemental real(wp) function gravity_accel(alt_km)
!! computes gravitational acceleration normal to Earth, where up is positive acceleration
real(wp), intent(in) :: alt_km
gravity_accel = -1 * 6.67408e-11_wp * 5.9722e24_wp / (6371.e3_wp + alt_km*1000)**2
!! [m s^-2]
end function gravity_accel
end module ionize_fang