Line data Source code
1 : !------------------------------------------------------------------------------
2 : ! Utility routines for medium energy electron ionization base on Ap
3 : ! geomagnetic activity index or observed electron fluxes
4 : !------------------------------------------------------------------------------
5 : module mee_ap_util_mod
6 : use shr_kind_mod, only: r8 => shr_kind_r8
7 : use shr_const_mod, only: pi => shr_const_pi
8 : use mee_fluxes, only: mee_fluxes_nenergy, mee_fluxes_energy, mee_fluxes_denergy
9 : use mee_fluxes, only: mee_fluxes_active, mee_fluxes_extract
10 : use cam_abortutils, only : endrun
11 :
12 : implicit none
13 :
14 : private
15 : public :: mee_ap_iprs
16 : public :: mee_ap_init
17 : public :: mee_ap_final
18 : public :: mee_ap_error
19 : public :: mee_ap_noerror
20 :
21 : integer, parameter :: mee_ap_error=-1
22 : integer, parameter :: mee_ap_noerror=0
23 :
24 : integer :: nbins=0
25 :
26 : real(r8),pointer :: energies(:) => null() ! energy bins
27 : real(r8),pointer :: denergies(:) => null() ! width of each energy bin
28 : real(r8) :: solid_angle_factor = -huge(1._r8)
29 : real(r8),allocatable :: fang_coefs(:,:)
30 :
31 : contains
32 :
33 : !-----------------------------------------------------------------------------
34 : ! Sets up energy spectrum grid for medium energy electrons in earth's
35 : ! radiation belt
36 : !-----------------------------------------------------------------------------
37 0 : subroutine mee_ap_init(loss_cone_angle, status)
38 :
39 : real(r8), intent(in) :: loss_cone_angle ! Bounce Loss Cone angle (degrees)
40 : integer, intent(out) :: status ! error status
41 :
42 : integer :: ierr
43 : character(len=*), parameter :: subname = 'mee_ap_init: '
44 :
45 0 : status = mee_ap_noerror
46 0 : if ( loss_cone_angle<0._r8 .or. loss_cone_angle>90._r8 ) then
47 0 : status = mee_ap_error
48 0 : return
49 : endif
50 :
51 : ! The area of the BLC in sr is 2pi(1-cos(BLC))
52 0 : solid_angle_factor = 2._r8*pi*(1._r8-cos(loss_cone_angle*pi/180._r8))
53 :
54 0 : if (mee_fluxes_active) then
55 0 : nbins=mee_fluxes_nenergy
56 0 : energies=>mee_fluxes_energy
57 0 : denergies=>mee_fluxes_denergy
58 : else
59 0 : nbins=100
60 0 : allocate(energies(nbins), stat=ierr)
61 0 : if (ierr/=0) call endrun(subname//'not able to allocate energies')
62 0 : allocate(denergies(nbins), stat=ierr)
63 0 : if (ierr/=0) call endrun(subname//'not able to allocate denergies')
64 0 : call gen_energy_grid(nbins, energies, denergies)
65 : endif
66 :
67 0 : call init_fang_coefs()
68 :
69 : end subroutine mee_ap_init
70 :
71 : !-----------------------------------------------------------------------------
72 : !-----------------------------------------------------------------------------
73 0 : subroutine mee_ap_final()
74 :
75 0 : deallocate(fang_coefs)
76 :
77 0 : if (.not.mee_fluxes_active) then
78 0 : deallocate(energies)
79 0 : deallocate(denergies)
80 : endif
81 0 : nullify(energies)
82 0 : nullify(denergies)
83 :
84 0 : end subroutine mee_ap_final
85 :
86 : !-----------------------------------------------------------------------------
87 : ! Computes ion pair production rates base on Ap geomagnetic activity index
88 : ! or observed electron fluxes
89 : !-----------------------------------------------------------------------------
90 0 : subroutine mee_ap_iprs( ncols, nlyrs, airden, scaleh, Ap, ionpairs, status, maglat, lshell)
91 :
92 : integer ,intent(in) :: ncols, nlyrs
93 : real(r8),intent(in) :: airden(ncols,nlyrs) ! g/cm3
94 : real(r8),intent(in) :: scaleh(ncols,nlyrs) ! cm
95 : real(r8),intent(in) :: Ap ! geomagnetic activity index
96 : real(r8),intent(out) :: ionpairs(ncols,nlyrs) ! pairs/cm3/sec
97 : integer, intent(out) :: status ! error status
98 : real(r8),intent(in), optional :: maglat(ncols) ! magnetic latitude (radians)
99 : real(r8),intent(in), optional :: lshell(ncols) ! magnetosphere L-Shells
100 :
101 : integer :: i,k
102 0 : real(r8) :: flux_sd(nbins)
103 0 : logical :: valid(nbins)
104 0 : real(r8) :: flux(nbins)
105 0 : real(r8) :: ipr(nbins,nlyrs)
106 0 : real(r8) :: l_shells(ncols)
107 :
108 0 : status = 0
109 :
110 0 : if (present(lshell)) then
111 0 : l_shells(:) = lshell(:)
112 0 : elseif (present(maglat)) then
113 : ! get L-shell values corresponeding to the column geo-mag latitudes
114 0 : l_shells = maglat2lshell( maglat )
115 : else
116 0 : ionpairs(:,:) = -huge(1._r8)
117 0 : status = mee_ap_error
118 : return
119 : endif
120 :
121 0 : ionpairs(:,:) = 0._r8
122 :
123 0 : do i = 1,ncols
124 :
125 0 : if ( l_shells(i) >= 2._r8 .and. l_shells(i) <= 10._r8 ) then
126 :
127 0 : valid(:) = .false.
128 0 : flux_sd(:) = 0._r8
129 :
130 0 : if (mee_fluxes_active) then
131 : ! use prescribed top-of-atmosphere fluxes
132 : call mee_fluxes_extract( l_shells(i), flux_sd, valid )
133 : end if
134 :
135 0 : where ( (.not.valid(:)) .and. (energies(:)>=30._r8) .and. (energies(:)<=1000._r8) )
136 : ! calculate the top of the atmosphere energetic electron energy spectrum
137 : ! van de Kamp is per steradian (electrons / (cm2 sr s keV))
138 : flux_sd(:) = FluxSpectrum(energies(:), l_shells(i), Ap)
139 : end where
140 :
141 : ! assume flux is isotropic inside a nominal bounce loss cone (BLC) angle.
142 : ! The area of the BLC in sr is 2pi(1-cosd(BLC))
143 0 : flux(:) = solid_angle_factor*flux_sd(:)
144 :
145 : ! calculate the IPR as a function of height and energy
146 0 : ipr(:,:) = iprmono(energies, flux, airden(i,:), scaleh(i,:))
147 :
148 : ! integrate across the energy range to get total IPR
149 0 : do k=1,nlyrs
150 0 : ionpairs(i,k) = sum(ipr(:,k)*denergies(:))
151 : end do
152 :
153 : end if
154 :
155 : end do
156 :
157 0 : end subroutine mee_ap_iprs
158 :
159 : !------------------------------------------------------------------------------
160 : ! Electron fluxes for a specific L-shell and Ap
161 : !
162 : ! Based on:
163 : !
164 : ! van de Kamp, M., Seppala, A., Clilverd, M. A., Rodger, C. J., Verronen, P. T., and Whittaker, I. C. (2016),
165 : ! A model providing long‐term datasets of energetic electron precipitation during geomagnetic storms,
166 : ! J. Geophys. Res. Atmos., 121, 12,520– 12,540,
167 : ! [doi:10.1002/2015JD024212](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2015JD024212)
168 : !------------------------------------------------------------------------------
169 0 : pure function FluxSpectrum( en, lshell, Ap ) result(flux)
170 :
171 : real(r8), intent(in) :: en(:) ! electron energy bins (keV)
172 : real(r8), intent(in) :: lshell ! magnetosphere L-Shell number
173 : real(r8), intent(in) :: Ap ! geomagnetic activity index
174 :
175 : real(r8) :: flux(size(en))
176 :
177 : real(r8) :: lpp
178 : real(r8) :: Spp
179 : real(r8) :: A
180 : real(r8) :: b
181 : real(r8) :: c
182 : real(r8) :: s
183 : real(r8) :: d
184 : real(r8) :: F30
185 : real(r8) :: E
186 : real(r8) :: bk
187 : real(r8) :: sk
188 : real(r8) :: k
189 : real(r8) :: x
190 :
191 0 : lpp = -0.7430_r8*log(Ap) + 6.5257_r8
192 0 : Spp = lshell - lpp
193 :
194 : ! vdK2016 eqn.(8)
195 :
196 0 : A = 8.2091_r8*Ap**(0.16255_r8)
197 0 : b = 1.3754_r8*Ap**(0.33042_r8)
198 0 : c = 0.13334_r8*Ap**(0.42616_r8)
199 0 : s = 2.2833_r8*Ap**(-0.22990_r8)
200 0 : d = 2.7563e-4_r8*Ap**(2.6116_r8)
201 :
202 0 : F30 = exp(A) / (exp(-b*(Spp-s)) + exp(c*(Spp-s)) + d)
203 :
204 : ! vdK2016 eqn.(9)
205 :
206 0 : E = 3.3777_r8*Ap**(-1.7038_r8) + 0.15_r8
207 0 : bk = 3.7632_r8*Ap**(-0.16034_r8)
208 0 : sk = 12.184_r8*Ap**(-0.30111_r8)
209 :
210 0 : k = -1.0_r8 / (E*exp(-bk*Spp) + 0.30450_r8*cosh(0.20098_r8*(Spp-sk))) - 1._r8
211 :
212 0 : x=k+1
213 0 : c = F30*x / ((1e3_r8**x) - (30._r8**x))
214 0 : flux(:) = c * (en(:)**k)
215 :
216 0 : end function FluxSpectrum
217 :
218 : !------------------------------------------------------------------------------
219 : ! Calculate how energy from top of atmosphere is deposited in rest of atmosphere
220 : !
221 : ! The function takes the energy spectrum at the top of the atmosphere and
222 : ! calculates how that energy is deposited in the atmosphere using the parameterization
223 : ! described in [Fang et al., (2010)](https://opensky.ucar.edu/islandora/object/articles:10653)
224 : !
225 : ! Fang, X., C. E. Randall, D. Lummerzheim, W. Wang, G. Lu, S. C. Solomon,
226 : ! and R. A. Frahm (2010), Parameterization of monoenergetic electron impact
227 : ! ionization, Geophys. Res. Lett., 37, L22106, [doi:10.1029/2010GL045406.]
228 : ! (https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2010GL045406)
229 : !
230 : ! Application of the new parameterization requires the following steps:
231 : !
232 : ! 1. Fang coefficients using equation (5) and Table 1. are calculated at initialization time
233 : ! 2. Calculate the y values throughout the atmosphere using equation (1).
234 : ! 3. Calculate the normalized energy dissipation f values using equation (4).
235 : ! 4. Obtain the altitude profile of qtot by substituting the f values into equation (3).
236 : !
237 : !------------------------------------------------------------------------------
238 0 : function iprmono(e, flux, rho, scaleh) result(ipr)
239 : real(r8), intent(in) :: e(:) ! electron energy bins (keV)
240 : real(r8), intent(in) :: flux(:) ! top of atmos electron fluxes (electrons / (cm2 sr s keV))
241 : real(r8), intent(in) :: rho(:) ! density of atmos (g/cm3)
242 : real(r8), intent(in) :: scaleh(:) ! scale height (cm)
243 :
244 : real(r8) :: ipr(size(e),size(rho))
245 : integer :: nenergies, n
246 : integer :: nlayers, k
247 :
248 0 : real(r8) :: y(size(rho))
249 0 : real(r8) :: f(size(rho))
250 : real(r8) :: Qmono
251 :
252 : ! assign constants
253 : real(r8), parameter :: epsilon = 0.035_r8 ! keV energy loss per ion pair produced
254 :
255 0 : ipr = 0._r8
256 0 : nenergies = size(e)
257 0 : nlayers = size(rho)
258 :
259 0 : do n = 1,nenergies
260 :
261 : ! step 1. (eq. 1)
262 0 : y(:) = (2/e(n))*(rho(:)*scaleh(:)/6.0e-6_r8)**0.7_r8
263 :
264 0 : do k = 1,nlayers
265 0 : f(k) = fang(y(k), n)
266 : end do
267 :
268 : ! calculate ipr (qtot) using eq. (3) for a specified flux at ea. energy
269 0 : Qmono = flux(n)*e(n) ! (keV cm−2 s−1)
270 0 : ipr(n,:) = f(:)*Qmono/(epsilon*scaleh(:))
271 : end do
272 :
273 : contains
274 :
275 0 : pure function fang(y,n) result(f)
276 : real(r8), intent(in) :: y
277 : integer, intent(in) :: n ! energy ndx
278 :
279 : real(r8) :: f
280 : ! Input:
281 : ! y - normalized atmospheric column mass as a function of vertical location (z)
282 : ! Emono - is incident electron energy (keV)
283 : ! Output:
284 : ! f - quanity calculated by eqn. (4)
285 :
286 :
287 : ! eq. (4) - Normalized energy deposition
288 0 : f = fang_coefs(1,n) * (y**fang_coefs(2,n)) * exp(-fang_coefs(3,n)*y**fang_coefs(4,n)) &
289 0 : + fang_coefs(5,n) * (y**fang_coefs(6,n)) * exp(-fang_coefs(7,n)*y**fang_coefs(8,n))
290 :
291 0 : end function fang
292 :
293 : end function iprmono
294 :
295 : !------------------------------------------------------------------------------
296 : ! initializes the coeffs used in the fang function in iprmono
297 : !------------------------------------------------------------------------------
298 0 : subroutine init_fang_coefs
299 :
300 : integer :: n,i, ierr
301 :
302 : real(r8) :: lne, lne2, lne3
303 : ! Table 1. of Fang et al., 2010
304 : ! (https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2010GL045406)
305 : real(r8), parameter :: p1(8,4) = reshape( &
306 : (/ 1.24616E+0_r8, 1.45903E+0_r8, -2.42269E-1_r8, 5.95459E-2_r8, &
307 : 2.23976E+0_r8, -4.22918E-7_r8, 1.36458E-2_r8, 2.53332E-3_r8, &
308 : 1.41754E+0_r8, 1.44597E-1_r8, 1.70433E-2_r8, 6.39717E-4_r8, &
309 : 2.48775E-1_r8, -1.50890E-1_r8, 6.30894E-9_r8, 1.23707E-3_r8, &
310 : -4.65119E-1_r8, -1.05081E-1_r8, -8.95701E-2_r8, 1.22450E-2_r8, &
311 : 3.86019E-1_r8, 1.75430E-3_r8, -7.42960E-4_r8, 4.60881E-4_r8, &
312 : -6.45454E-1_r8, 8.49555E-4_r8, -4.28581E-2_r8, -2.99302E-3_r8, &
313 : 9.48930E-1_r8, 1.97385E-1_r8, -2.50660E-3_r8, -2.06938E-3_r8 /), &
314 : shape=(/8,4/),order=(/2,1/))
315 :
316 0 : allocate(fang_coefs(8,nbins), stat=ierr)
317 0 : if (ierr/=0) call endrun('init_fang_coefs: not able to allocate fang_coefs')
318 :
319 0 : do n = 1,nbins
320 : ! terms in eq. (5)
321 0 : lne = log(energies(n))
322 0 : lne2 = lne*lne
323 0 : lne3 = lne*lne2
324 :
325 : ! step 2. calculate the C array in (5)
326 0 : do i = 1,8
327 0 : fang_coefs(i,n) = exp(p1(i,1) + p1(i,2)*lne + p1(i,3)*lne2 + p1(i,4)*lne3)
328 : end do
329 :
330 : end do
331 0 : end subroutine init_fang_coefs
332 :
333 : !------------------------------------------------------------------------------
334 : ! Generate a grid of energy bins for the flux spectrum.
335 : ! The energy range of the spectrum is 30-1000 keV,
336 : ! with nbins of logarithmically spaced grid points.
337 : !------------------------------------------------------------------------------
338 0 : subroutine gen_energy_grid(nbins, energies, deltas)
339 : integer, intent(in) :: nbins
340 : real(r8),intent(out) :: energies(nbins)
341 : real(r8),intent(out) :: deltas(nbins)
342 :
343 : integer :: i
344 : real(r8) :: low,med,hig
345 :
346 : ! for energy bins ranging from 30 keV to 1000 keV
347 : real(r8), parameter :: e1 = log10(30._r8)
348 : real(r8), parameter :: e2 = log10(1000._r8)
349 :
350 0 : do i = 1,nbins
351 0 : low = e1 + (e2-e1)*(i-1.0_r8)/nbins
352 0 : med = e1 + (e2-e1)*(i-0.5_r8)/nbins
353 0 : hig = e1 + (e2-e1)*(i)/nbins
354 :
355 0 : energies(i) = 10**med
356 0 : deltas(i) = (10**hig)-(10**low)
357 : end do
358 :
359 0 : end subroutine gen_energy_grid
360 :
361 : !------------------------------------------------------------------------------
362 : ! returns L-Shell number for a given magnetic latitude (radians)
363 : !------------------------------------------------------------------------------
364 0 : pure elemental function maglat2lshell( rmaglat ) result(lshell)
365 : real(r8), intent(in) :: rmaglat ! mag latitude in radians
366 : real(r8) :: lshell ! magnetosphere L-Shell number
367 :
368 : ! lshell = r/(cos(rmaglat)**2) (https://en.wikipedia.org/wiki/L-shell)
369 : ! where r is the radial distance (in planetary radii) to a point on the line.
370 : ! r = 1.01 corresponds to an altitude in the lower mesosphere (~64 km)
371 : ! where medium-energy electrons typically deposit their energy.
372 0 : lshell = 1.01_r8/(cos(rmaglat)**2)
373 :
374 0 : end function maglat2lshell
375 :
376 : end module mee_ap_util_mod
|
