Line data Source code
1 :
2 : module mo_aurora
3 : !-----------------------------------------------------------------------
4 : !
5 : ! Auroral oval parameterization. See reference:
6 : ! R.G. Roble, E.C. Ridley
7 : ! An auroral model for the NCAR thermospheric general circulation model (TGCM)
8 : ! Annales Geophysicae,5A, (6), 369-382, 1987.
9 : !
10 : ! The aurora oval is a circle in auroral circle coordinates. Auroral circle
11 : ! coordinates are offset from magnetic coordinates by offa degrees (radians)
12 : ! towards 0 MLT and by dskofa degrees (radians) towards dusk (18 MLT).
13 : ! The aurora assumes a Maxwellian in energy, so that the characteristic
14 : ! energy is half of the mean energy (or mean energy = 2*alfa, where alfa
15 : ! is the characteristic energy). The Maxwellian is approximated in the
16 : ! aion subroutine.
17 : ! The aurora oval is assumed to be a Gaussian in auroral latitude, with
18 : ! peak values on the day (=1) and night (=2) sides that change from one to
19 : ! the other using cosines of the auroral longitude coordinate.
20 : ! There is provision for a low energy (~75 eV) aurora at the location of the
21 : ! regular (~1-6 keV) aurora in order to simulate the energy flux found
22 : ! at higher altitudes that is non-Maxwellian, but the flux is usually
23 : ! set to zero (1.e-80).
24 : ! There is provision for a proton (MeV) aurora, but the flux is usually
25 : ! set to zero (1.e-20).
26 : ! The drizzle is a constant low energy electron flux over the polar cap,
27 : ! which goes to 1/e over twice the half-width of the aurora at the
28 : ! radius of the aurora.
29 : ! The cusp is a low energy electron flux centered over the dayside convection
30 : ! entrance at phid at the convection reversal boundary theta0. The cusp
31 : ! falls off over 5 degrees in latitude and over 20 degrees in longitude
32 : ! to 1/e values of the peak at the center.
33 : ! 1.e-20 and 1.e-80 are used to give a near zero answer.
34 : !
35 : ! The polar drizzle and cusp electron energies are low, and soft particles
36 : ! have great influence on the over-all thermospheric and ionospheric
37 : ! structure, especially on the electron density profiles at mid-latitudes
38 : ! and in winter since low energy electrons produce ionization at high
39 : ! altitudes where loss rates are very low. (Comment by Wenbin Wang.)
40 : ! The original energies for drizzle and cusp were alfad=0.75, alfac=0.5 keV.
41 : ! The original guess at energy fluxes were: ed=0.1+2.0*power/100.,ec=0.1+0.9*power/100.
42 : ! The next guess at energy fluxes were: ed=0.01+0.2*power/100., ec=0.01+0.09*power/100.
43 : ! The values below reflect higher estimates for the electron energy (lower alt)
44 : !
45 : ! Calling sequence (all subs in mo_aurora, mo_aurora.F):
46 : ! 1) sub aurora_cons called once per time step from advance.
47 : ! 2) sub aurora called from dynamics, inside parallel latitude scan.
48 : ! 3) subs aurora_cusp and aurora_heat called from sub aurora.
49 : ! 4) sub aurora_ions called from sub aurora.
50 : !
51 : !-----------------------------------------------------------------------
52 :
53 : use shr_kind_mod, only: r8 => shr_kind_r8
54 : use mo_constants, only: pi, gask => rgas_cgs
55 : use cam_logfile, only: iulog
56 : use spmd_utils, only: masterproc
57 : use aurora_params, only: power=>hpower, plevel, aurora_params_set
58 : use aurora_params, only: ctpoten, theta0, dskofa, offa, phid, rrad
59 : use aurora_params, only: prescribed_period
60 :
61 : implicit none
62 :
63 : interface aurora
64 : module procedure aurora_prod
65 : module procedure aurora_hrate
66 : end interface
67 :
68 : save
69 :
70 :
71 : private
72 : public :: aurora_inti, aurora_timestep_init, aurora
73 : public :: aurora_register
74 :
75 : integer, parameter :: isouth = 1
76 : integer, parameter :: inorth = 2
77 :
78 : ! g = 8.7 m/s^2? Because this is 400 km up?
79 : real(r8), parameter :: grav = 870._r8 ! (cm/s^2)
80 :
81 : integer :: lev1 = 1
82 : real(r8) :: rmass_o1
83 : real(r8) :: rmass_o2
84 : real(r8) :: rmass_n2
85 : real(r8) :: rmassinv_o1
86 : real(r8) :: rmassinv_o2
87 : real(r8) :: rmassinv_n2
88 : real(r8), parameter :: twopi = 2._r8*pi
89 : real(r8), parameter :: d2r = pi/180._r8
90 : real(r8), parameter :: r2d = 180._r8/pi
91 :
92 : !-----------------------------------------------------------------------
93 : ! ... polar drizzle parameters:
94 : ! alfad: Characteristic Maxwellian energy of drizzle electrons (keV)
95 : ! ed : Column energy input of drizzle electrons (ergs/cm**2/s)
96 : ! fd : Electron particle flux of drizzle electrons (particles/cm**2/s)
97 : !-----------------------------------------------------------------------
98 : real(r8), parameter :: alfad = 0.5_r8
99 : real(r8) :: ed
100 : real(r8) :: fd ! set in sub aurora_ions
101 :
102 : !-----------------------------------------------------------------------
103 : ! ... polar cusp parameters:
104 : ! alfac: Characteristic Maxwellian energy of polar cusp electons (keV)
105 : ! ec : Column energy input of polar cusp electrons (ergs/cm**2/s)
106 : ! fc : Electron particle flux of polar cusp electrons (particles/cm**2/s)
107 : !-----------------------------------------------------------------------
108 : real(r8), parameter :: alfac = 0.1_r8
109 : real(r8) :: ec
110 : real(r8) :: fc ! set in sub aurora_ions
111 :
112 : !-----------------------------------------------------------------------
113 : ! e1: Peak energy flux in noon sector of the aurora (ergs/cm**2/s)
114 : ! e2: Peak energy flux in midnight sector of the aurora (ergs/cm**2/s)
115 : ! h1: Gaussian half-width of the noon auroral oval in degrees
116 : ! h2: Gaussian half-width of the midnight auroral oval in degrees
117 : !-----------------------------------------------------------------------
118 : real(r8) :: &
119 : e1, e2, & ! set in sub aurora_cons (function of hem power)
120 : h1, h2 ! set in sub aurora_cons (function of hem power)
121 :
122 : !-----------------------------------------------------------------------
123 : ! ... additional auroral parameters
124 : !-----------------------------------------------------------------------
125 : real(r8) :: &
126 : alfa0, & ! average of noon and midnight characteristic Maxw energies
127 : ralfa,ralfa2, & ! difference ratios of characteristic energies
128 : rrote, & ! clockwise rotation from noon of peak dayside energy flux (e1)
129 : rroth, & ! clockwise rotation from noon of dayside h1 Gaussian half-width
130 : h0, & ! average of noon and midnight Gaussian half-widths
131 : rh, & ! difference ratio of half-widths (rh=(h2-h1)/(h2+h1))
132 : e0,e20, & ! e0 = average of noon and midnight electrons
133 : ree,re2, & ! difference ratios of peak energy fluxes (ree=(e2-e1)/(e2+e1))
134 : alfa20 ! average of noon and midnight char energies for high alt aurora
135 :
136 : logical :: aurora_active = .false.
137 : integer :: indxAIPRS = -1
138 : integer :: indxQTe = -1
139 : integer :: indxEfx = -1
140 : integer :: indxKev = -1
141 :
142 : real(r8), parameter :: h2deg = 15._r8 ! hour to degree
143 :
144 : contains
145 :
146 :
147 : !----------------------------------------------------------------------
148 : !----------------------------------------------------------------------
149 0 : subroutine aurora_register
150 : use ppgrid, only : pver,pcols
151 : use physics_buffer, only : pbuf_add_field, dtype_r8
152 :
153 : ! add ionization rates to phys buffer for waccmx ionosphere module
154 :
155 : ! Sum of ion auroral production rates for O2
156 0 : call pbuf_add_field('AurIPRateSum', 'physpkg', dtype_r8, (/pcols,pver/), indxAIPRS)
157 0 : call pbuf_add_field('QTeAur', 'physpkg', dtype_r8, (/pcols/), indxQTe) ! for electron temperature
158 :
159 0 : endsubroutine aurora_register
160 :
161 1536 : subroutine aurora_inti(pbuf2d)
162 : !-----------------------------------------------------------------------
163 : ! ... initialize aurora module
164 : !-----------------------------------------------------------------------
165 :
166 0 : use ppgrid, only : pver
167 : use constituents, only : cnst_get_ind, cnst_mw
168 : use ref_pres, only : pref_mid
169 : use mo_chem_utls, only : get_spc_ndx
170 : use cam_history, only : addfld, horiz_only
171 : use physics_buffer,only: pbuf_get_index
172 : use infnan, only : nan, assignment(=)
173 : use physics_buffer, only: physics_buffer_desc, pbuf_set_field
174 :
175 : implicit none
176 :
177 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
178 : !-----------------------------------------------------------------------
179 : ! ... local variables
180 : !-----------------------------------------------------------------------
181 : integer :: k, m
182 : real(r8), parameter :: e = 1.e-10_r8
183 :
184 : real(r8) :: plb
185 : real(r8) :: alfa_1, alfa_2, alfa21, alfa22
186 : real(r8) :: e21, e22
187 : integer :: op_ndx,o2p_ndx,np_ndx,n2p_ndx,e_ndx
188 : integer :: ierr
189 : real(r8) :: x_nan
190 :
191 1536 : indxEfx = pbuf_get_index('AUREFX', errcode=ierr)
192 1536 : indxKev = pbuf_get_index('AURKEV', errcode=ierr)
193 :
194 1536 : if (indxEfx>0 .and. indxKev>0) then
195 0 : x_nan = nan
196 0 : call pbuf_set_field(pbuf2d, indxEfx, x_nan)
197 0 : call pbuf_set_field(pbuf2d, indxKev, x_nan)
198 : endif
199 :
200 1536 : if (indxAIPRS>0) then
201 0 : call pbuf_set_field(pbuf2d, indxAIPRS, 0._r8)
202 : endif
203 :
204 1536 : theta0(:) = nan
205 1536 : offa(:) = nan
206 1536 : dskofa(:) = nan
207 1536 : phid(:) = nan
208 1536 : rrad(:) = nan
209 1536 : ctpoten = nan
210 1536 : power = nan
211 :
212 1536 : op_ndx = get_spc_ndx( 'Op' )
213 1536 : o2p_ndx = get_spc_ndx( 'O2p' )
214 1536 : np_ndx = get_spc_ndx( 'Np' )
215 1536 : n2p_ndx = get_spc_ndx( 'N2p' )
216 1536 : e_ndx = get_spc_ndx( 'e' )
217 :
218 : aurora_active = op_ndx > 0 .and. o2p_ndx > 0 .and. np_ndx > 0 .and. n2p_ndx > 0 .and. e_ndx > 0 &
219 1536 : .and. pref_mid(1) < 0.1_r8 ! need high-top
220 :
221 1536 : if (.not. aurora_active) return
222 :
223 : !-----------------------------------------------------------------------
224 : ! ... initialize module variables
225 : !-----------------------------------------------------------------------
226 :
227 : !-----------------------------------------------------------------------
228 : ! ... set molecular weights
229 : !-----------------------------------------------------------------------
230 1536 : call cnst_get_ind( 'O2', m )
231 1536 : rmass_o2 = cnst_mw(m)
232 1536 : rmassinv_o2 = 1._r8/rmass_o2
233 1536 : call cnst_get_ind( 'O', m )
234 1536 : rmass_o1 = cnst_mw(m)
235 1536 : rmassinv_o1 = 1._r8/rmass_o1
236 1536 : call cnst_get_ind( 'N', m )
237 1536 : rmass_n2 = 2._r8*cnst_mw(m)
238 1536 : rmassinv_n2 = 1._r8/rmass_n2
239 :
240 1536 : offa(isouth) = 1.0_r8*d2r
241 1536 : offa(inorth) = 1.0_r8*d2r
242 :
243 1536 : alfa_1 = 1.5_r8
244 1536 : alfa_2 = 2._r8
245 :
246 : !-----------------------------------------------------------------------
247 : ! Values from 10/05/94 HPI estimates (50% or more higher than old estimates):
248 : ! alfa_1 = amin1(1.5,1.25+0.05*plevel)
249 : ! alfa_2 = 1.2 + 0.095*plevel
250 : !-----------------------------------------------------------------------
251 1536 : alfa0 = 0.5_r8*(alfa_1 + alfa_2)
252 1536 : ralfa = (alfa_2 - alfa_1) / (alfa_1 + alfa_2 + e)
253 1536 : alfa21 = 0.075_r8
254 1536 : alfa22 = 0.075_r8
255 1536 : alfa20 = 0.5_r8 * (alfa21 + alfa22)
256 1536 : ralfa2 = (alfa22 - alfa21) / (alfa21 + alfa22 + e)
257 1536 : e21 = 1.e-80_r8
258 1536 : e22 = 1.e-80_r8
259 1536 : e20 = 0.5_r8 * (e21 + e22)
260 1536 : re2 = (e22 - e21) / (e21 + e22)
261 :
262 : !-----------------------------------------------------------------------
263 : ! ... set auroral lower bndy index
264 : !-----------------------------------------------------------------------
265 1536 : plb = 5.e-4_r8*exp( 7._r8 ) * .1_r8 ! Pa
266 16896 : do k = 1,pver
267 16896 : if( pref_mid(k) >= plb ) then
268 1536 : lev1 = k-1
269 1536 : exit
270 : end if
271 : end do
272 :
273 1536 : if (masterproc) write(iulog,*) ' '
274 1536 : if (masterproc) write(iulog,*) 'aurora_inti: aurora will go down to lev,p = ',lev1,pref_mid(lev1)
275 1536 : if (masterproc) write(iulog,*) ' '
276 :
277 : !-----------------------------------------------------------------------
278 : ! Report to stdout:
279 : !-----------------------------------------------------------------------
280 : #ifdef AURORA_DIAGS
281 : write(iulog,"(/,'aurora_cons:')")
282 : ! write(iulog,"(' cusp: alfac=',f8.3,' ec=',f8.3,' fc=',e10.4)") &
283 : ! alfac,ec,fc
284 : ! write(iulog,"(' drizzle: alfad=',f8.3,' ed=',f8.3,' fd=',e10.4)") &
285 : ! alfad,ed,fd
286 : write(iulog,"(' half-widths = h1,h2=',2f10.3)") h1,h2
287 : write(iulog,"(' energy flux = e1,e2=',2f10.3)") e1,e2
288 : write(iulog,"(' add_sproton = ',l1)") add_sproton
289 : write(iulog,"(' ')")
290 : #endif
291 : call addfld('ALATM', horiz_only, 'I','degrees', &
292 1536 : 'Magnetic latitude at each geographic coordinate')
293 : call addfld('ALONM', horiz_only, 'I','degrees', &
294 1536 : 'Magnetic longitude at each geographic coordinate')
295 : call addfld( 'QSUM', (/ 'lev' /), 'I','/s', &
296 3072 : 'total ion production' )
297 :
298 1536 : end subroutine aurora_inti
299 :
300 16128 : subroutine aurora_timestep_init( )
301 : !-----------------------------------------------------------------------
302 : ! ... per timestep initialization
303 : !-----------------------------------------------------------------------
304 :
305 1536 : use heelis_mod, only : heelis_update
306 :
307 : !-----------------------------------------------------------------------
308 : ! ... local variables
309 : !-----------------------------------------------------------------------
310 : real(r8) :: roth, rote
311 : real(r8), parameter :: convert = 3.1211e8_r8
312 :
313 16128 : if (.not. aurora_active) return
314 :
315 16128 : if (.not.aurora_params_set) then
316 16128 : call heelis_update() ! use heelis if aurora parameters are not already set
317 : endif
318 :
319 16128 : if( power >= 1.0_r8 ) then
320 16128 : plevel = 2.09_r8*log( power )
321 : else
322 0 : plevel = 0._r8
323 : end if
324 :
325 : #ifdef AURORA_DIAGS
326 : if( masterproc ) then
327 : write(iulog,*) '----------------------------------------'
328 : write(iulog,*) 'aurora_timestep_init: power,ctpoten = ',power,ctpoten
329 : write(iulog,*) '----------------------------------------'
330 : end if
331 : #endif
332 :
333 : !-----------------------------------------------------------------------
334 : ! h1 = Gaussian half-width of the noon auroral oval in degrees
335 : ! h2 = Gaussian half-width of the midnight auroral oval in degrees
336 : !-----------------------------------------------------------------------
337 : ! produce realistic oval compared to NOAA empirical auroral oval and TIMED/GUVI
338 : ! h1 formula given by Wenbin base on POLARVIS image;
339 : ! h2 formula based on Emery et al original auroral parameterization report
340 16128 : h1 = min(2.35_r8, 0.83_r8 + 0.33_r8*plevel)
341 16128 : h2 = 2.5_r8+0.025_r8*max(power,55._r8)+0.01_r8*min(0._r8,power-55._r8)
342 :
343 : !-----------------------------------------------------------------------
344 : ! Values from corrections to Emery et al Parameterization report:
345 : ! h1 = amin1(2.35, 0.83 + 0.33*plevel)
346 : ! h2 = 2.87 + 0.15*plevel
347 : !-----------------------------------------------------------------------
348 :
349 16128 : rh = (h2 - h1) / (h1 + h2)
350 16128 : h0 = 0.5_r8 * (h1 + h2) * d2r
351 :
352 :
353 : ! roth = MLT of max width of aurora in hours
354 : ! rote = MLT of max energy flux of aurora in hours
355 :
356 16128 : roth = 0.81_r8 - 0.06_r8 * plevel
357 16128 : rote = 0.17_r8 - 0.04_r8 * plevel
358 :
359 : ! Convert MLT from hours to degrees to radians
360 :
361 16128 : rroth = roth * h2deg * d2r
362 16128 : rrote = rote * h2deg * d2r
363 :
364 : !-----------------------------------------------------------------------
365 : ! e1 = energy flux in the noon sector of the aurora (ergs/cm**2/s)
366 : ! e2 = energy flux in the midnight sector of the aurora (ergs/cm**2/s)
367 : !-----------------------------------------------------------------------
368 : ! produce realistic oval compared to NOAA empirical auroral oval and TIMED/GUVI
369 : ! e1 formula given by Wenbin base on POLARVIS image;
370 : ! e2 formula based on Emery et al original auroral parameterization report
371 : !-----------------------------------------------------------------------
372 16128 : e1 = max(0.50_r8, -2.15_r8 + 0.62_r8*plevel)
373 16128 : e2=1._r8+0.11_r8*power
374 :
375 : !-----------------------------------------------------------------------
376 : ! ed : Column energy input of drizzle electrons (ergs/cm**2/s)
377 : ! ec : Column energy input of polar cusp electrons (ergs/cm**2/s)
378 : !-----------------------------------------------------------------------
379 16128 : ed = .0012_r8+.0006_r8*power
380 16128 : ec = (0.24_r8+0.0067_r8*power)/5._r8
381 :
382 : !-----------------------------------------------------------------------
383 : ! Set cusp and drizzle parameters:
384 : ! (conversion between particle number density and characteristic
385 : ! energy and column energy input)
386 : !-----------------------------------------------------------------------
387 16128 : fc = convert * ec / alfac
388 16128 : fd = convert * ed / alfad
389 :
390 : !-----------------------------------------------------------------------
391 : ! Values from corrections to Emery et al Parameterization report:
392 : !-----------------------------------------------------------------------
393 16128 : e0 = 0.5_r8 * (e1 + e2)
394 16128 : ree = (e2 - e1) / (e1 + e2)
395 :
396 : end subroutine aurora_timestep_init
397 :
398 0 : subroutine aurora_prod( tn, o2, o1, mbar, rlats, &
399 0 : qo2p, qop, qn2p, qnp, pmid, &
400 : lchnk, calday, ncol, rlons, pbuf )
401 : !-----------------------------------------------------------------------
402 : ! ... auroral parameterization driver
403 : !-----------------------------------------------------------------------
404 :
405 : use mo_apex, only : alatm, alonm ! magnetic latitude,longitude grid (radians)
406 : use mo_apex, only : maglon0
407 : use ppgrid, only : pcols, pver
408 : use cam_history, only : outfld
409 : use physics_buffer,only: physics_buffer_desc,pbuf_get_field
410 :
411 : implicit none
412 :
413 : !-----------------------------------------------------------------------
414 : ! ... dummy arguments
415 : !-----------------------------------------------------------------------
416 : integer, intent(in) :: &
417 : ncol, & ! column count
418 : lchnk ! chunk index
419 : real(r8), intent(in) :: &
420 : calday ! calendar day of year
421 : real(r8), intent(in) :: &
422 : tn(pcols,pver), & ! neutral gas temperature (K)
423 : o2(ncol,pver), & ! O2 concentration (kg/kg)
424 : o1(ncol,pver), & ! O concentration (kg/kg)
425 : mbar(ncol,pver) ! mean molecular weight (g/mole)
426 : real(r8), intent(in) :: &
427 : pmid(pcols,pver) ! midpoint pressure (Pa)
428 : real(r8), intent(in) :: &
429 : rlats(ncol), & ! column latitudes (radians)
430 : rlons(ncol)
431 : real(r8), intent(out) :: &
432 : qo2p(ncol,pver), & ! o2+ production
433 : qop(ncol,pver), & ! o+ production
434 : qn2p(ncol,pver), & ! n2+ production
435 : qnp(ncol,pver) ! n+ production
436 :
437 : type(physics_buffer_desc),pointer :: pbuf(:)
438 :
439 : !-----------------------------------------------------------------------
440 : ! ... local variables
441 : !-----------------------------------------------------------------------
442 : integer :: i, k
443 0 : integer :: hemis(ncol)
444 : real(r8) :: ofda, cosofa, sinofa, aslona
445 0 : real(r8) :: dlat_aur(ncol)
446 0 : real(r8) :: dlon_aur(ncol)
447 0 : real(r8) :: colat(ncol)
448 0 : real(r8) :: sinlat(ncol)
449 0 : real(r8) :: coslat(ncol)
450 0 : real(r8) :: coslon(ncol)
451 0 : real(r8) :: sinlon(ncol)
452 0 : real(r8) :: alon(ncol)
453 0 : real(r8) :: cusp(ncol)
454 0 : real(r8) :: alfa(ncol)
455 0 : real(r8) :: alfa2(ncol)
456 0 : real(r8) :: flux(ncol)
457 0 : real(r8) :: flux2(ncol)
458 0 : real(r8) :: drizl(ncol)
459 0 : logical :: do_aurora(ncol)
460 :
461 : real(r8) :: dayfrac, rotation
462 0 : real(r8), pointer :: qteaur(:) ! for electron temperature
463 :
464 0 : if (.not. aurora_active) return
465 :
466 : !-----------------------------------------------------------------------
467 : ! ... initialize ion production
468 : !-----------------------------------------------------------------------
469 0 : do k = 1,pver
470 0 : qo2p(:,k) = 0._r8
471 0 : qop(:,k) = 0._r8
472 0 : qn2p(:,k) = 0._r8
473 0 : qnp(:,k) = 0._r8
474 : end do
475 :
476 : !-----------------------------------------------------------------------
477 : ! ... output mag lons, lats
478 : !-----------------------------------------------------------------------
479 0 : call outfld( 'ALONM', r2d*alonm(:ncol,lchnk), ncol, lchnk )
480 0 : call outfld( 'ALATM', r2d*alatm(:ncol,lchnk), ncol, lchnk )
481 :
482 0 : if (indxQTe>0) then
483 0 : call pbuf_get_field(pbuf, indxQTe, qteaur)
484 0 : qteaur(:) = 0._r8
485 : endif
486 :
487 : !-----------------------------------------------------------------------
488 : ! aurora is active for columns poleward of 30 deg
489 : !-----------------------------------------------------------------------
490 0 : do_aurora(:) = abs( rlats(:) ) > pi/6._r8
491 0 : if( all( .not. do_aurora(:) ) ) then
492 : return
493 : end if
494 :
495 : !-----------------------------------------------------------------------
496 : ! ... set rotation based on sun location
497 : !-----------------------------------------------------------------------
498 0 : dayfrac = (calday - int(calday))
499 0 : rotation = maglon0 + dayfrac*twopi-pi
500 :
501 0 : do i = 1,ncol
502 0 : if( do_aurora(i) ) then
503 0 : dlat_aur(i) = alatm(i,lchnk)
504 0 : dlon_aur(i) = alonm(i,lchnk) + rotation ! rotate it
505 0 : if( dlon_aur(i) > pi ) then
506 0 : dlon_aur(i) = dlon_aur(i) - twopi
507 0 : else if( dlon_aur(i) < -pi ) then
508 0 : dlon_aur(i) = dlon_aur(i) + twopi
509 : end if
510 0 : if( dlat_aur(i) > 0._r8 ) then
511 0 : hemis(i) = 2
512 : else
513 0 : hemis(i) = 1
514 : end if
515 : !-----------------------------------------------------------------------
516 : ! ... find auroral circle coordinates
517 : !-----------------------------------------------------------------------
518 0 : ofda = sqrt( offa(hemis(i))**2 + dskofa(hemis(i))**2)
519 0 : cosofa = cos( ofda )
520 0 : sinofa = sin( ofda )
521 0 : aslona = asin( dskofa(hemis(i))/ofda )
522 0 : sinlat(i) = sin( abs( dlat_aur(i) ) )
523 0 : coslat(i) = cos( dlat_aur(i) )
524 0 : sinlon(i) = sin( dlon_aur(i) + aslona )
525 0 : coslon(i) = cos( dlon_aur(i) + aslona )
526 0 : colat(i) = acos( cosofa*sinlat(i) - sinofa*coslat(i)*coslon(i))
527 : alon(i) = mod( atan2( sinlon(i)*coslat(i),sinlat(i)*sinofa &
528 0 : + cosofa*coslat(i)*coslon(i) ) - aslona + 3._r8*pi,twopi) - pi
529 : end if
530 : end do
531 : #ifdef AURORA_DIAGS
532 : write(iulog,*) '-----------------------------------------------------'
533 : write(iulog,*) 'aurora: diagnostics for lchnk = ',lchnk
534 : write(iulog,*) ' geo lats'
535 : write(iulog,'(1p,5g15.7)') r2d*rlats(:ncol)
536 : write(iulog,*) ' geo lons'
537 : write(iulog,'(1p,5g15.7)') r2d*rlons(:ncol)
538 : write(iulog,*) ' mag lats'
539 : write(iulog,'(1p,5g15.7)') r2d*dlat_aur(:ncol)
540 : write(iulog,*) ' mag lons'
541 : write(iulog,'(1p,5g15.7)') r2d*alonm(:ncol,lchnk)
542 : write(iulog,*) ' mag table lons'
543 : write(iulog,'(1p,5g15.7)') r2d*dlon_aur(:ncol)
544 : write(iulog,*) ' min,max mag lons = ',r2d*minval(alonm(:ncol,lchnk)),r2d*maxval(alonm(:ncol,lchnk))
545 : write(iulog,*) '-----------------------------------------------------'
546 : #endif
547 :
548 : !-----------------------------------------------------------------------
549 : ! ... make cusp
550 : !-----------------------------------------------------------------------
551 0 : call aurora_cusp( cusp, do_aurora, hemis, colat, alon, ncol )
552 :
553 : !-----------------------------------------------------------------------
554 : ! ... make alfa, flux, and drizzle
555 : !-----------------------------------------------------------------------
556 : call aurora_heat( flux, flux2, alfa, alfa2, &
557 : drizl, do_aurora, hemis, &
558 0 : alon, colat, ncol, pbuf )
559 :
560 : !-----------------------------------------------------------------------
561 : ! ... auroral additions to ionization rates
562 : !-----------------------------------------------------------------------
563 : call aurora_ions( drizl, cusp, alfa, alfa2, &
564 : flux, flux2, tn, o2, &
565 : o1, mbar, qo2p, qop, qn2p, &
566 0 : qnp, pmid, do_aurora, ncol, lchnk, pbuf )
567 :
568 0 : end subroutine aurora_prod
569 :
570 161280 : subroutine aurora_hrate( tn, mbar, rlats, &
571 80640 : aur_hrate, cpair, pmid, lchnk, calday, &
572 : ncol, rlons, pbuf )
573 : !-----------------------------------------------------------------------
574 : ! ... auroral parameterization driver
575 : !-----------------------------------------------------------------------
576 :
577 0 : use mo_apex, only : alatm, alonm ! magnetic latitude,longitude grid (radians)
578 : use mo_apex, only : maglon0
579 : use ppgrid, only : pcols, pver
580 : use physics_buffer,only: physics_buffer_desc
581 :
582 : implicit none
583 :
584 : !-----------------------------------------------------------------------
585 : ! ... dummy arguments
586 : !-----------------------------------------------------------------------
587 : integer, intent(in) :: &
588 : ncol, & ! column count
589 : lchnk ! chunk index
590 : real(r8), intent(in) :: &
591 : calday ! calendar day of year
592 : real(r8), intent(in) :: &
593 : tn(pcols,pver), & ! neutral gas temperature (K)
594 : mbar(ncol,pver) ! mean molecular weight (g/mole)
595 : real(r8), intent(in) :: &
596 : cpair(ncol,pver) ! specific heat capacity (J/K/kg)
597 : real(r8), intent(in) :: &
598 : pmid(pcols,pver) ! midpoint pressure (Pa)
599 : real(r8), intent(in) :: &
600 : rlats(ncol), & ! column latitudes (radians)
601 : rlons(ncol)
602 : real(r8), intent(out) :: &
603 : aur_hrate(ncol,pver) ! auroral heating rate
604 : type(physics_buffer_desc),pointer :: pbuf(:)
605 :
606 : !-----------------------------------------------------------------------
607 : ! ... local variables
608 : !-----------------------------------------------------------------------
609 : real(r8), parameter :: aur_therm = 807._r8
610 : real(r8), parameter :: jkcal = 4184._r8
611 : real(r8), parameter :: aur_heat_eff = .05_r8
612 : real(r8), parameter :: aur_hconst = 1.e3_r8*jkcal*aur_therm*aur_heat_eff
613 :
614 : integer :: i, k
615 161280 : integer :: hemis(ncol)
616 : real(r8) :: ofda, cosofa, sinofa, aslona
617 161280 : real(r8) :: dlat_aur(ncol)
618 161280 : real(r8) :: dlon_aur(ncol)
619 161280 : real(r8) :: colat(ncol)
620 161280 : real(r8) :: sinlat(ncol)
621 161280 : real(r8) :: coslat(ncol)
622 161280 : real(r8) :: coslon(ncol)
623 161280 : real(r8) :: sinlon(ncol)
624 161280 : real(r8) :: alon(ncol)
625 161280 : real(r8) :: cusp(ncol)
626 161280 : real(r8) :: alfa(ncol)
627 161280 : real(r8) :: alfa2(ncol)
628 161280 : real(r8) :: flux(ncol)
629 161280 : real(r8) :: flux2(ncol)
630 161280 : real(r8) :: drizl(ncol)
631 161280 : real(r8) :: qsum(ncol,pver) ! total ion production (1/s)
632 161280 : logical :: do_aurora(ncol)
633 :
634 : real(r8) :: dayfrac, rotation
635 :
636 : !-----------------------------------------------------------------------
637 : ! ... initialize ion production
638 : !-----------------------------------------------------------------------
639 5725440 : do k = 1,pver
640 87010560 : aur_hrate(:,k) = 0._r8
641 : end do
642 :
643 80640 : if (.not. aurora_active) return
644 :
645 : !-----------------------------------------------------------------------
646 : ! ... check latitudes, and return if all below 32.5 deg
647 : !-----------------------------------------------------------------------
648 1241856 : do_aurora(:) = abs( rlats(:) ) > pi/6._r8
649 80640 : if( all( .not. do_aurora(:) ) ) then
650 : return
651 : end if
652 :
653 : !-----------------------------------------------------------------------
654 : ! ... set rotation based on sun location
655 : !-----------------------------------------------------------------------
656 80640 : dayfrac = (calday - int(calday))
657 80640 : rotation = maglon0 + dayfrac*twopi-pi
658 :
659 1241856 : do i = 1,ncol
660 1241856 : if( do_aurora(i) ) then
661 774144 : dlat_aur(i) = alatm(i,lchnk)
662 774144 : dlon_aur(i) = alonm(i,lchnk) + rotation ! rotate it
663 774144 : if( dlon_aur(i) > pi ) then
664 0 : dlon_aur(i) = dlon_aur(i) - twopi
665 774144 : else if( dlon_aur(i) < -pi ) then
666 353910 : dlon_aur(i) = dlon_aur(i) + twopi
667 : end if
668 774144 : if( dlat_aur(i) > 0._r8 ) then
669 387072 : hemis(i) = 2
670 : else
671 387072 : hemis(i) = 1
672 : end if
673 : !-----------------------------------------------------------------------
674 : ! ... find auroral circle coordinates
675 : !-----------------------------------------------------------------------
676 774144 : ofda = sqrt( offa(hemis(i))**2 + dskofa(hemis(i))**2)
677 774144 : cosofa = cos( ofda )
678 774144 : sinofa = sin( ofda )
679 774144 : aslona = asin( dskofa(hemis(i))/ofda )
680 774144 : sinlat(i) = sin( abs( dlat_aur(i) ) )
681 774144 : coslat(i) = cos( dlat_aur(i) )
682 774144 : sinlon(i) = sin( dlon_aur(i) + aslona )
683 774144 : coslon(i) = cos( dlon_aur(i) + aslona )
684 774144 : colat(i) = acos( cosofa*sinlat(i) - sinofa*coslat(i)*coslon(i))
685 : alon(i) = mod( atan2( sinlon(i)*coslat(i),sinlat(i)*sinofa &
686 774144 : + cosofa*coslat(i)*coslon(i) ) - aslona + 3._r8*pi,twopi) - pi
687 : end if
688 : end do
689 : #ifdef AURORA_DIAGS
690 : write(iulog,*) '-----------------------------------------------------'
691 : write(iulog,*) 'aurora: diagnostics for lchnk = ',lchnk
692 : write(iulog,*) ' geo lats'
693 : write(iulog,'(1p,5g15.7)') r2d*rlats(:ncol)
694 : write(iulog,*) ' geo lons'
695 : write(iulog,'(1p,5g15.7)') r2d*rlons(:ncol)
696 : write(iulog,*) ' mag lats'
697 : write(iulog,'(1p,5g15.7)') r2d*dlat_aur(:ncol)
698 : write(iulog,*) ' mag lons'
699 : write(iulog,'(1p,5g15.7)') r2d*alonm(:ncol,lchnk)
700 : write(iulog,*) ' mag table lons'
701 : write(iulog,'(1p,5g15.7)') r2d*dlon_aur(:ncol)
702 : write(iulog,*) ' min,max mag lons = ',r2d*minval(alonm(:ncol,lchnk)),r2d*maxval(alonm(:ncol,lchnk))
703 : write(iulog,*) '-----------------------------------------------------'
704 : #endif
705 :
706 : !-----------------------------------------------------------------------
707 : ! ... make cusp
708 : !-----------------------------------------------------------------------
709 80640 : call aurora_cusp( cusp, do_aurora, hemis, colat, alon, ncol )
710 :
711 : !-----------------------------------------------------------------------
712 : ! ... make alfa, flux, and drizzle
713 : !-----------------------------------------------------------------------
714 : call aurora_heat( flux, flux2, alfa, alfa2, &
715 : drizl, do_aurora, hemis, &
716 80640 : alon, colat, ncol, pbuf )
717 :
718 : !-----------------------------------------------------------------------
719 : ! ... auroral additions to ionization rates
720 : !-----------------------------------------------------------------------
721 : call total_ion_prod( drizl, cusp, alfa, alfa2, &
722 : flux, flux2, tn, &
723 : mbar, qsum, pmid, do_aurora, &
724 80640 : ncol )
725 :
726 : !-----------------------------------------------------------------------
727 : ! ... form auroral heating rate
728 : !-----------------------------------------------------------------------
729 5725440 : do k = 1,pver
730 87010560 : aur_hrate(:,k) = aur_hconst * qsum(:,k) / (cpair(:,k) * mbar(:,k))
731 : end do
732 :
733 80640 : end subroutine aurora_hrate
734 :
735 80640 : subroutine aurora_cusp( cusp, do_aurora, hemis, colat, alon, ncol )
736 : !-----------------------------------------------------------------------
737 : ! ... calculate horizontal variation of polar cusp heating
738 : !-----------------------------------------------------------------------
739 :
740 : implicit none
741 :
742 : !-----------------------------------------------------------------------
743 : ! ... dummy arguments
744 : !-----------------------------------------------------------------------
745 : integer, intent(in) :: ncol
746 : integer, intent(in) :: hemis(ncol)
747 : real(r8), intent(in) :: colat(ncol)
748 : real(r8), intent(in) :: alon(ncol)
749 : real(r8), intent(out) :: cusp(ncol)
750 : logical, intent(in) :: do_aurora(ncol)
751 :
752 : !-----------------------------------------------------------------------
753 : ! ... local variables
754 : !-----------------------------------------------------------------------
755 : real(r8), parameter :: s5 =.08726646_r8, &
756 : s20 =.34906585_r8
757 :
758 1241856 : where( do_aurora(:) )
759 774144 : cusp(:) = (exp( -((theta0(hemis(:)) - colat(:))/s5)**2 ) &
760 : + exp( -((pi - theta0(hemis(:)) - colat(:))/s5)**2) ) &
761 : *exp( -(atan2( sin(alon(:) - phid(hemis(:))), cos(alon(:) - phid(hemis(:))) )/s20)**2 )
762 : elsewhere
763 : cusp(:) = 0._r8
764 : endwhere
765 :
766 80640 : end subroutine aurora_cusp
767 :
768 161280 : subroutine aurora_heat( flux, flux2, alfa, alfa2, &
769 80640 : drizl, do_aurora, hemis, &
770 80640 : alon, colat, ncol, pbuf )
771 : !-----------------------------------------------------------------------
772 : ! ... calculate alfa, flux, and drizzle
773 : !-----------------------------------------------------------------------
774 : use physics_buffer,only: physics_buffer_desc,pbuf_get_field
775 :
776 : implicit none
777 :
778 : !-----------------------------------------------------------------------
779 : ! ... dummy arguments
780 : !-----------------------------------------------------------------------
781 : integer, intent(in) :: ncol
782 : integer, intent(in) :: hemis(ncol)
783 : real(r8), intent(in) :: colat(ncol)
784 : real(r8), intent(in) :: alon(ncol)
785 : real(r8), intent(inout) :: flux(ncol)
786 : real(r8), intent(inout) :: flux2(ncol)
787 : real(r8), intent(inout) :: drizl(ncol)
788 : real(r8), intent(inout) :: alfa(ncol)
789 : real(r8), intent(inout) :: alfa2(ncol)
790 : logical, intent(in) :: do_aurora(ncol)
791 : type(physics_buffer_desc),pointer :: pbuf(:)
792 :
793 : !-----------------------------------------------------------------------
794 : ! ... local variables
795 : !-----------------------------------------------------------------------
796 : real(r8), dimension(ncol) :: &
797 161280 : coslamda, & ! cos(angle from throat)
798 161280 : halfwidth, & ! oval half-width
799 161280 : wrk, & ! temp wrk array
800 161280 : dtheta ! latitudinal variation (Gaussian)
801 : real(r8) :: ekev
802 80640 : real(r8), pointer :: pr_efx(:) ! Pointer to pbuf prescribed energy flux (mW m-2)
803 80640 : real(r8), pointer :: pr_kev(:) ! Pointer to pbuf prescribed mean energy (keV)
804 80640 : real(r8), pointer :: qteaur(:) ! for electron temperature
805 : integer :: n
806 :
807 : !-----------------------------------------------------------------------
808 : ! Low-energy protons:
809 : !
810 : ! alfap0 = 0.5*(alfap1+alfap2)
811 : ! e0p = 0.5*(pe1+pe2)
812 : !
813 : ! coslamda = cos(lamda)
814 : ! halfwidth = auroral half width
815 : ! dtheta = colat-theta0(ihem)
816 : ! alfa = electron energy
817 : !-----------------------------------------------------------------------
818 4725504 : where( do_aurora(:) )
819 : coslamda(:) = cos( atan2( sin( alon(:) - rrote ),cos( alon(:) - rrote ) ) )
820 : !-----------------------------------------------------------------------
821 : ! ... auroral oval half-width (equation (1) in Roble,1987):
822 : !-----------------------------------------------------------------------
823 : halfwidth(:) = h0*(1._r8 - rh*cos( atan2( sin(alon(:) - rroth),cos( alon(:) - rroth ) ) ) )
824 774144 : dtheta(:) = colat(:) - rrad(hemis(:))
825 : endwhere
826 : !-----------------------------------------------------------------------
827 : ! ... characteristic energy (equation (2) in Roble,1987):
828 : !-----------------------------------------------------------------------
829 80640 : if( alfa0 > .01_r8 ) then
830 1241856 : where( do_aurora(:) )
831 : alfa(:) = alfa0*(1._r8 - ralfa*coslamda(:))
832 : endwhere
833 : else
834 0 : alfa(:) = 0._r8
835 : end if
836 :
837 1241856 : wrk = 0._r8
838 7047936 : where( do_aurora(:) )
839 : !-----------------------------------------------------------------------
840 : ! ... flux, drizzle, alfa2, flux2
841 : !-----------------------------------------------------------------------
842 : !-----------------------------------------------------------------------
843 : ! ... energy flux (equation (3) in Roble,1987):
844 : !-----------------------------------------------------------------------
845 : wrk(:) = exp( -(dtheta(:)/halfwidth(:))**2 )
846 : flux(:) = e0*(1._r8 - ree*coslamda(:))*wrk(:) / (2._r8*alfa(:)*1.602e-9_r8)
847 : drizl(:) = exp( -((dtheta(:) + abs(dtheta(:)))/(2._r8*h0))**2 )
848 : alfa2(:) = alfa20*(1._r8 - ralfa2*coslamda(:))
849 : flux2(:) = e20*(1._r8 - re2*coslamda(:))*wrk(:) / (2._r8*alfa2(:)*1.602e-9_r8)
850 : endwhere
851 :
852 : !-----------------------------------------------------------------------
853 : ! ... for electron temperature (used in settei):
854 : !-----------------------------------------------------------------------
855 80640 : if (indxQTe>0) then
856 0 : call pbuf_get_field(pbuf, indxQTe, qteaur)
857 : ! The factor of -7e8 is based on the energy per ion pair and heating efficiency from Roble (1987)
858 0 : qteaur(:ncol) = -7.e8_r8*wrk(:ncol)
859 : end if
860 :
861 : !----------------------------------------------------------------------------------------------
862 : ! ... If turned on, use amie energy flux and mean energy to replace flux(:) and alfa(:)
863 : !----------------------------------------------------------------------------------------------
864 80640 : if (prescribed_period .and. indxEfx>0 .and. indxKev>0) then
865 : !---------------------------------------------------------------------------
866 : ! Overwrite with prescribed mean energy and energy flux in physics buffer
867 : !---------------------------------------------------------------------------
868 0 : call pbuf_get_field(pbuf, indxEfx, pr_efx)
869 0 : call pbuf_get_field(pbuf, indxKev, pr_kev)
870 0 : do n=1,ncol
871 0 : ekev = max(pr_kev(n),1._r8)
872 0 : alfa(n) = ekev/2._r8
873 0 : flux(n) = max(pr_efx(n)/(ekev*1.602e-9_r8),1.e-20_r8)
874 : enddo
875 : endif
876 :
877 161280 : end subroutine aurora_heat
878 :
879 0 : subroutine aurora_ions( drizl, cusp, alfa1, alfa2, &
880 0 : flux1, flux2, tn, o2, &
881 0 : o1, mbar, qo2p, qop, qn2p, &
882 0 : qnp, pmid, do_aurora, ncol, lchnk, pbuf )
883 : !-----------------------------------------------------------------------
884 : ! ... calculate auroral additions to ionization rates
885 : !-----------------------------------------------------------------------
886 :
887 80640 : use ppgrid, only : pcols, pver
888 : use cam_history, only : outfld
889 :
890 : use physics_buffer,only: physics_buffer_desc, pbuf_get_field
891 :
892 : implicit none
893 :
894 : !-----------------------------------------------------------------------
895 : ! ... dummy arguments
896 : !-----------------------------------------------------------------------
897 : integer, intent(in) :: ncol
898 : integer, intent(in) :: lchnk
899 : real(r8), intent(in), dimension(ncol) :: &
900 : drizl, &
901 : cusp, &
902 : alfa1, &
903 : alfa2, &
904 : flux1, &
905 : flux2
906 : real(r8), dimension(pcols,pver), intent(in) :: &
907 : tn, & ! midpoint neutral temperature (K)
908 : pmid ! midpoint pressure (Pa)
909 : real(r8), dimension(ncol,pver), intent(in) :: &
910 : o2, & ! midpoint o2 concentration (kg/kg)
911 : o1, & ! midpoint o concentration (kg/kg)
912 : mbar ! mean molecular mass (g/mole)
913 : real(r8), dimension(ncol,pver), intent(inout) :: &
914 : qo2p, & ! o2p prod from aurora (molecules/cm^3/s)
915 : qop, & ! op prod from aurora (molecules/cm^3/s)
916 : qn2p, & ! n2p prod from aurora (molecules/cm^3/s)
917 : qnp ! np prod from aurora (molecules/cm^3/s)
918 : logical, intent(in) :: do_aurora(ncol)
919 :
920 : type(physics_buffer_desc),pointer :: pbuf(:)
921 :
922 : !-----------------------------------------------------------------------
923 : ! ... local variables
924 : !-----------------------------------------------------------------------
925 : real(r8), parameter :: const0 = 1.e-20_r8
926 :
927 : integer :: k
928 : real(r8), dimension(ncol) :: &
929 0 : p0ez, &
930 0 : press, & ! pressure at interface levels (dyne/cm^2)
931 0 : tempi, & ! temperature at interface levels (K)
932 0 : xalfa1, &
933 0 : xalfa2, &
934 0 : xcusp, &
935 0 : xdrizl, & ! input to sub aion
936 0 : cusp_ion, &
937 0 : drizl_ion, & ! output from sub aion
938 0 : alfa1_ion, &
939 0 : alfa2_ion, &
940 0 : barm_t, &
941 0 : qsum, &
942 0 : denom, &
943 0 : barm, &
944 0 : falfa1, &
945 0 : falfa2, &
946 0 : fcusp, &
947 0 : fdrizl, &
948 0 : xn2
949 : real(r8), dimension(ncol) :: &
950 0 : qo2p_aur, &
951 0 : qop_aur, &
952 0 : qn2p_aur ! auroral ionization for O2+, O+, N2+
953 : real(r8) :: qia(5) ! low energy proton source (not in use, 1/02)
954 0 : real(r8) :: wrk(ncol,pver)
955 :
956 0 : real(r8), pointer :: aurIPRateSum(:,:) ! Pointer to pbuf auroral ion production sum for O2+,O+,N2+ (s-1 cm-3)
957 :
958 0 : qia(:) = 0._r8
959 0 : wrk(:,:) = 0._r8
960 :
961 : !-----------------------------------------------------------
962 : ! Point to production rates array in physics buffer where
963 : ! rates will be stored for ionosphere module access. Also,
964 : ! initialize rates to zero before column loop since only
965 : ! daylight values are filled
966 : !-----------------------------------------------------------
967 0 : if (indxAIPRS>0) then
968 0 : call pbuf_get_field(pbuf, indxAIPRS, aurIPRateSum)
969 0 : aurIPRateSum(:,:) = 0._r8
970 : endif
971 :
972 : level_loop : &
973 0 : do k = 1,lev1
974 0 : where( do_aurora(:) )
975 0 : press(:ncol) = 10._r8*pmid(:ncol,k) ! from Pa to dyne/cm^2
976 : tempi(:ncol) = tn(:ncol,k)
977 : barm(:) = mbar(:,k)
978 : p0ez(:) = (press(:)/(grav*4.e-6_r8))**.606_r8
979 : xalfa1(:) = p0ez(:)/alfa1(:)
980 : xalfa2(:) = p0ez(:)/alfa2(:)
981 : xcusp (:) = p0ez(:)/alfac
982 : xdrizl(:) = p0ez(:)/alfad
983 :
984 : !-----------------------------------------------------------------------
985 : ! ... initialize (whole array operations):
986 : !-----------------------------------------------------------------------
987 : alfa1_ion(:) = const0
988 : alfa2_ion(:) = const0
989 : cusp_ion(:) = const0
990 : drizl_ion(:) = const0
991 : endwhere
992 : !-----------------------------------------------------------------------
993 : ! ... auroral electrons
994 : !-----------------------------------------------------------------------
995 0 : call aion( xalfa1, alfa1_ion, do_aurora, ncol )
996 0 : call aion( xalfa2, alfa2_ion, do_aurora, ncol )
997 0 : call aion( xcusp , cusp_ion, do_aurora, ncol )
998 0 : call aion( xdrizl, drizl_ion, do_aurora, ncol )
999 0 : where( do_aurora(:) )
1000 : falfa1(:) = alfa1(:)*flux1(:) ! s7
1001 : falfa2(:) = alfa2(:)*flux2(:) ! s8
1002 : fcusp (:) = cusp(:)*alfac*fc ! s9
1003 : fdrizl(:) = drizl(:)*alfad*fd ! s10
1004 : qsum(:) = falfa1(:)*alfa1_ion(:) & ! s7*s3
1005 : + falfa2(:)*alfa2_ion(:) & ! s8*s4
1006 : + fcusp(:)*cusp_ion (:) & ! s9*s5
1007 : + fdrizl(:)*drizl_ion(:) ! s10*s6
1008 : endwhere
1009 :
1010 : !-----------------------------------------------------------------------
1011 : ! ... form production
1012 : !-----------------------------------------------------------------------
1013 0 : where( do_aurora(:) )
1014 : barm_t(:) = grav*barm(:)/(35.e-3_r8*gask*tempi(:))
1015 : qsum(:) = qsum(:)*barm_t(:) ! s1 = s1*s11
1016 : wrk(:,k) = qsum(:)
1017 : !-----------------------------------------------------------------------
1018 : ! ... denominator of equations (13-16) in Roble,1987.
1019 : !-----------------------------------------------------------------------
1020 : xn2(:) = max( (1._r8 - o2(:,k) - o1(:,k)),1.e-8_r8 )
1021 : denom(:) = 0.92_r8*xn2(:)*rmassinv_n2 &
1022 : + 1.5_r8*o2(:,k) *rmassinv_o2 + 0.56_r8*o1(:,k) *rmassinv_o1
1023 : !-----------------------------------------------------------------------
1024 : ! ... production of O2+ (equation (15) in Roble,1987):
1025 : !-----------------------------------------------------------------------
1026 : qo2p_aur(:) = qsum(:)*o2(:,k)/(rmass_o2*denom(:)) + qia(2)
1027 : !-----------------------------------------------------------------------
1028 : ! ... production of O+ (equation (16) in Roble,1987):
1029 : !-----------------------------------------------------------------------
1030 : qop_aur(:) = qsum(:)*(.5_r8 *o2(:,k)*rmassinv_o2 &
1031 : + .56_r8*o1(:,k)*rmassinv_o1)/denom(:) + qia(3)
1032 : !-----------------------------------------------------------------------
1033 : ! ... production of N2+ (equation (13) in Roble,1987)
1034 : !-----------------------------------------------------------------------
1035 : qn2p_aur(:) = qsum(:)*.7_r8*xn2(:)/(rmass_n2*denom(:)) + qia(1)
1036 : qo2p(:,k) = qo2p(:,k) + qo2p_aur(:)
1037 : qop(:,k) = qop(:,k) + qop_aur(:)
1038 : qn2p(:,k) = qn2p(:,k) + qn2p_aur(:)
1039 : qnp(:,k) = qnp (:,k) + .22_r8/.7_r8 * qn2p_aur(:)
1040 : endwhere
1041 : end do level_loop
1042 :
1043 : !----------------------------------------------------------------
1044 : ! Store the sum of the ion production rates in pbuf to be used
1045 : ! in the ionosx module
1046 : !----------------------------------------------------------------
1047 0 : if (indxAIPRS>0) then
1048 :
1049 0 : aurIPRateSum(1:ncol,1:pver) = wrk(1:ncol,1:pver)
1050 :
1051 : endif
1052 :
1053 0 : call outfld( 'QSUM', wrk, ncol, lchnk )
1054 :
1055 0 : end subroutine aurora_ions
1056 :
1057 80640 : subroutine total_ion_prod( drizl, cusp, alfa1, alfa2, &
1058 80640 : flux1, flux2, tn, &
1059 80640 : mbar, tpions, pmid, do_aurora, &
1060 : ncol )
1061 : !-----------------------------------------------------------------------
1062 : ! ... calculate auroral additions to ionization rates
1063 : !-----------------------------------------------------------------------
1064 :
1065 0 : use ppgrid, only : pcols, pver
1066 :
1067 : implicit none
1068 :
1069 : !-----------------------------------------------------------------------
1070 : ! ... dummy arguments
1071 : !-----------------------------------------------------------------------
1072 : integer, intent(in) :: ncol
1073 : real(r8), intent(in), dimension(ncol) :: &
1074 : drizl, &
1075 : cusp, &
1076 : alfa1, &
1077 : alfa2, &
1078 : flux1, &
1079 : flux2
1080 : real(r8), dimension(pcols,pver), intent(in) :: &
1081 : tn, & ! midpoint neutral temperature (K)
1082 : pmid ! midpoint pressure (Pa)
1083 : real(r8), dimension(ncol,pver), intent(in) :: &
1084 : mbar ! mean molecular mass (g/mole)
1085 : real(r8), dimension(ncol,pver), intent(inout) :: &
1086 : tpions ! total ion production (1/s)
1087 : logical, intent(in) :: do_aurora(ncol)
1088 :
1089 : !-----------------------------------------------------------------------
1090 : ! ... local variables
1091 : !-----------------------------------------------------------------------
1092 : real(r8), parameter :: const0 = 1.e-20_r8
1093 :
1094 : integer :: k
1095 : real(r8), dimension(ncol) :: &
1096 161280 : p0ez, &
1097 161280 : press, & ! pressure at interface levels (dyne/cm^2)
1098 161280 : tempi, & ! temperature at interface levels (K)
1099 161280 : xalfa1, &
1100 161280 : xalfa2, &
1101 161280 : xcusp, &
1102 161280 : xdrizl, & ! input to sub aion
1103 161280 : cusp_ion, &
1104 161280 : drizl_ion, & ! output from sub aion
1105 161280 : alfa1_ion, &
1106 161280 : alfa2_ion, &
1107 161280 : barm_t, &
1108 161280 : qsum, &
1109 161280 : barm, &
1110 161280 : falfa1, &
1111 161280 : falfa2, &
1112 80640 : fcusp
1113 :
1114 87010560 : tpions(:,:) = 0._r8
1115 :
1116 : level_loop : &
1117 887040 : do k = 1,lev1
1118 152570880 : where( do_aurora(:) )
1119 806400 : press(:ncol) = 10._r8*pmid(:ncol,k) ! from Pa to dyne/cm^2
1120 : tempi(:ncol) = tn(:ncol,k)
1121 : barm(:) = mbar(:,k)
1122 : p0ez(:) = (press(:)/(grav*4.e-6_r8))**.606_r8
1123 : xalfa1(:) = p0ez(:)/alfa1(:)
1124 : xalfa2(:) = p0ez(:)/alfa2(:)
1125 : xcusp (:) = p0ez(:)/alfac
1126 : xdrizl(:) = p0ez(:)/alfad
1127 :
1128 : !-----------------------------------------------------------------------
1129 : ! ... initiliaze (whole array operations):
1130 : !-----------------------------------------------------------------------
1131 : alfa1_ion(:) = const0
1132 : alfa2_ion(:) = const0
1133 : cusp_ion(:) = const0
1134 : drizl_ion(:) = const0
1135 : endwhere
1136 : !-----------------------------------------------------------------------
1137 : ! ... auroral electrons
1138 : !-----------------------------------------------------------------------
1139 806400 : call aion( xalfa1, alfa1_ion, do_aurora, ncol )
1140 806400 : call aion( xalfa2, alfa2_ion, do_aurora, ncol )
1141 806400 : call aion( xcusp , cusp_ion, do_aurora, ncol )
1142 806400 : call aion( xdrizl, drizl_ion, do_aurora, ncol )
1143 58867200 : where( do_aurora(:) )
1144 : falfa1(:) = alfa1(:)*flux1(:) ! s7
1145 : falfa2(:) = alfa2(:)*flux2(:) ! s8
1146 : fcusp (:) = cusp(:)*alfac*fc ! s9
1147 : qsum(:) = falfa1(:)*alfa1_ion(:) & ! s7*s3
1148 : + falfa2(:)*alfa2_ion(:) & ! s8*s4
1149 : + fcusp(:)*cusp_ion (:) & ! s9*s5
1150 : + drizl(:)*drizl_ion(:) ! s10*s6
1151 : endwhere
1152 :
1153 : !-----------------------------------------------------------------------
1154 : ! ... form production
1155 : !-----------------------------------------------------------------------
1156 35723520 : where( do_aurora(:) )
1157 : barm_t(:) = grav*barm(:)/(35.e-3_r8*gask*tempi(:))
1158 : tpions(:,k) = qsum(:)*barm_t(:) ! s1 = s1*s11
1159 : endwhere
1160 : end do level_loop
1161 :
1162 80640 : end subroutine total_ion_prod
1163 :
1164 3225600 : subroutine aion( si, so, do_aurora, ncol )
1165 : !-----------------------------------------------------------------------
1166 : ! Calculates integrated f(x) needed for total auroral ionization.
1167 : ! See equations (10-12) in Roble,1987.
1168 : ! Coefficients for equation (12) of Roble,1987 are in variable cc
1169 : ! (revised since 1987):
1170 : ! Uses the identity x**y = exp(y*ln(x)) for performance
1171 : ! (fewer (1/2) trancendental functions are required).
1172 : !------------------------------------------------------------------------
1173 :
1174 : implicit none
1175 :
1176 : !------------------------------------------------------------------------
1177 : ! ... dummy arguments
1178 : !------------------------------------------------------------------------
1179 : integer, intent(in) :: ncol
1180 : real(r8), intent(in) :: si(ncol)
1181 : real(r8), intent(out) :: so(ncol)
1182 : logical, intent(in) :: do_aurora(ncol)
1183 :
1184 : !------------------------------------------------------------------------
1185 : ! ... local variables
1186 : !------------------------------------------------------------------------
1187 : real(r8), parameter :: cc(8) = &
1188 : (/ 3.2333134511131_r8 , 2.5658873458085_r8 , 2.2540957232641_r8 , &
1189 : 0.72971983372673_r8, 1.1069072431948_r8 , 1.7134937681128_r8 , &
1190 : 1.8835442312993_r8 , 0.86472135072090_r8 /)
1191 :
1192 6451200 : real(r8) :: xlog(ncol)
1193 :
1194 189020160 : where( do_aurora(:) )
1195 : xlog(:) = log( si(:) )
1196 : so(:) = cc(1)*exp( cc(2)*xlog(:) - cc(3)*exp( cc(4)*xlog(:) ) ) &
1197 : + cc(5)*exp( cc(6)*xlog(:) - cc(7)*exp( cc(8)*xlog(:) ) )
1198 : elsewhere
1199 : so(:) = 0._r8
1200 : endwhere
1201 :
1202 3225600 : end subroutine aion
1203 :
1204 : end module mo_aurora
|
