!
      module aurora_module
!
! Auroral oval parameterization. See reference:
! R.G. Roble, E.C. Ridley
! An auroral model for the NCAR thermospheric general circulation model (TGCM)
! Annales Geophysicae,5A, (6), 369-382, 1987. 
!
! The aurora oval is a circle in auroral circle coordinates.  Auroral circle
!  coordinates are offset from magnetic coordinates by offa degrees (radians)
!  towards 0 MLT and by dskofa degrees (radians) towards dusk (18 MLT).
! The aurora assumes a Maxwellian in energy, so that the characteristic
!  energy is half of the mean energy (or mean energy = 2*alfa, where alfa
!  is the characteristic energy).  The Maxwellian is approximated in the
!  aion and bion subroutines.
! The aurora oval is assumed to be a Gaussian in auroral latitude, with
!  peak values on the day (=1) and night (=2) sides that change from one to
!  the other using cosines of the auroral longitude coordinate.
! There is provision for a low energy (~75 eV) aurora at the location of the
!  regular (~1-6 keV) aurora in order to simulate the energy flux found
!  at higher altitudes that is non-Maxwellian, but the flux is usually
!  set to zero (1.e-80).
! There is provision for a proton (MeV) aurora, but the flux is usually
!  set to zero (1.e-20).
! The drizzle is a constant low energy electron flux over the polar cap,
!  which goes to 1/e over twice the half-width of the aurora at the
!  radius of the aurora.
! The cusp is a low energy electron flux centered over the dayside convection
!  entrance at phid at the convection reversal boundary theta0.  The cusp
!  falls off over 5 degrees in latitude and over 20 degrees in longitude
!  to 1/e values of the peak at the center.
! 1.e-20 and 1.e-80 are used to give a near zero answer.
!
! The polar drizzle and cusp electron energies are low, and soft particles
!  have great influence on the over-all thermospheric and ionospheric
!  structure, especially on the electron density profiles at mid-latitudes
!  and in winter since low energy electrons produce ionization at high
!  altitudes where loss rates are very low.  (Comment by Wenbin Wang.)
! The original energies for drizzle and cusp were alfad=0.75, alfac=0.5 keV.
! The original guess at energy fluxes were: ed=0.1+2.0*power/100.,ec=0.1+0.9*power/100.
! The next guess at energy fluxes were: ed=0.01+0.2*power/100., ec=0.01+0.09*power/100.
! The values below reflect higher estimates for the electron energy (lower alt)
!
! Calling sequence (all subs in aurora_module, aurora.F):
!   1) sub aurora_cons called once per time step from advance.
!   2) sub aurora called from dynamics, inside parallel latitude scan.
!   3) subs aurora_cusp and aurora_heat called from sub aurora.
!   4) sub aurora_ions called from sub aurora. 
!
      use cons_module,only: ! see cons.F
     |  dphi, ! delta latitude geographic (radians)
     |  dtr,  ! degrees to radians (pi/180)
     |  pi    ! 4.*atan(1.)
      use params_module,only: nlon,nlonp4,nlevp1,nlat
      use protons_module,only: get_protons
      use meped,only: get_meped
      use hist_module,only: nstep
      use init_module,only: istep
      implicit none
!
! VT vampir tracing:
!
#ifdef VT
#include <VT.inc>
#endif
!
! Grid dimensions:
      logical :: do_aurora ! flag to calculate aurora at current lat
      integer :: ihem      ! hemisphere at current latitude (1=SH,2=NH)
      integer,parameter :: isouth=1, inorth=2
!
! Polar drizzle parameters:
!   alfad: Characteristic Maxwellian energy of drizzle electrons (keV)
!   ed   : Column energy input of drizzle electrons (ergs/cm**2/s)
!   fd   : Electron particle flux of drizzle electrons (particles/cm**2/s)
! 4/25/08 btf: reduce alfad to 0.5, as per Wenbin:
!
      real,parameter :: alfad = 0.5
      real :: ed ! set in sub aurora_cons as function of hem power
      real :: fd ! set in sub aurora_ions
!
! Polar cusp parameters:
!   alfac: Characteristic Maxwellian energy of polar cusp electons (keV)
!   ec   : Column energy input of polar cusp electrons (ergs/cm**2/s)
!   fc   : Electron particle flux of polar cusp electrons (particles/cm**2/s)
!
! 1/24/08 btf: As per Wenbin: alfac = 0.1 (old tiegcm1.8 value was 1.0),
!              and ec is set as function of hemispheric power in aurora_cons
!              (old tiegcm1.8 value for ec was 0.5)
!
      real,parameter :: alfac = 0.1
      real :: ec ! set in sub aurora_cons (function of hem power)
      real :: fc ! set in sub aurora_ions
!
! e1: Peak energy flux in noon sector of the aurora (ergs/cm**2/s)
! e2: Peak energy flux in midnight sector of the aurora (ergs/cm**2/s)
! h1: Gaussian half-width of the noon auroral oval in degrees
! h2: Gaussian half-width of the midnight auroral oval in degrees
!
      real ::
     |  e1,e2, ! set in sub aurora_cons (function of hem power)
     |  h1,h2  ! set in sub aurora_cons (function of hem power)
!
! Solar proton parameters for the polar cap (time-gcm only)
!   alfa_sp: Characteristic Maxwellian energy of solar protons (MeV) (was alfad2)
!   e_sp   : Column energy input of solar protons (ergs/cm**2/s) (was ed2)
!   flx_sp : e_sp/1.602e-6, for input to sub bion                (was fd2)
!
      real,parameter ::
     |  alfa_sp_default = 10., ! in tgcm24 this is user namelist ED2
     |  e_sp_default = 1.e-20
      real,save :: alfa_sp, e_sp
      real :: flx_sp
!
! High energy electron parameters in the auroral oval (time-gcm only):
      logical :: add_helectron = .true. ! true in tgcm24 (see s13*s15 in qsum)
      real ::
     |  alfa30 = 40., ! Characteristic energy of auroral electrons (Mev)
     |  e30 = .05     ! Column energy of auroral electrons (ergs/cm**2/s)
!    |  e30 = 1.e-20  ! Column energy of auroral electrons (ergs/cm**2/s)
!
! 1/02 btf: low energy protons not included at this time.
!     real,parameter :: ! low energy protons:
!    |  alfap1 = 5.,
!    |  alfap2 = 15.,
!    |  pe1 = 1.e-20,
!    |  pe2 = 1.e-20
!     real :: alfap0,e0p
!
! Additional auroral parameters (see sub aurora_cons):
! (dimension 2 is for south, north hemispheres)
      real ::
     |  theta0(2),    ! convection reversal boundary in radians
     |  offa(2),      ! offset of oval towards 0 MLT relative to magnetic pole (rad)
     |  dskofa(2),    ! offset of oval in radians towards 18 MLT (f(By))
     |  phid(2),      ! dayside convection entrance in MLT converted to radians (f(By))
     |  phin(2),      ! night convection entrance in MLT converted to radians (f(By))
     |  rrad(2),      ! radius of auroral circle in radians
     |  alfa0,        ! average of noon and midnight characteristic Maxw energies
     |  ralfa,ralfa2, ! difference ratios of characteristic energies
     |  rrote,        ! clockwise rotation from noon of peak dayside energy flux (e1)
     |  rroth,        ! clockwise rotation from noon of dayside h1 Gaussian half-width
     |  h0,           ! average of noon and midnight Gaussian half-widths
     |  rh,           ! difference ratio of half-widths (rh=(h2-h1)/(h2+h1))
     |  e0,e20,       ! e0 = average of noon and midnight electrons
     |  ree,re2,      ! difference ratios of peak energy fluxes (ree=(e2-e1)/(e2+e1))
     |  alfa20        ! average of noon and midnight char energies for high alt aurora
!
! The following parameters are used only by heelis module for dynamo:
      real ::
     |  offc(2),      !
     |  dskofc(2),    !
     |  psim(2),      ! 
     |  psie(2),      !
     |  pcen(2),      !
     |  phidp0(2),    !
     |  phidm0(2),    !
     |  phinp0(2),    !
     |  phinm0(2),    !
     |  rr1(2)
!
! Auroral circle coordinates (see sub aurora):
      real ::
     |  dlat_aur(nlon), dlon_aur(nlon),
     |  colat(nlon),  sinlat(nlon), coslat(nlon),
     |  coslon(nlon), sinlon(nlon), alon(nlon)
!
! Local-time dependent energy and particle fluxes:
! (see subs aurora_heat and aurora_cusp)
      real ::
     |  alfa(nlonp4), alfa2(nlonp4), alfa3(nlonp4), ! characteristic energy
     |  flux(nlonp4), flux2(nlonp4), flux3(nlonp4), ! particle flux
     |  drizl(nlonp4),cusp(nlonp4),                 ! polar drizzle and cusp
     |  qteaur(nlonp4,nlat)                         ! for electron temperature
! TEMP
      real :: cusp2d(nlonp4,nlat),drzl2d(nlonp4,nlat),
     |        alfa2d(nlonp4,nlat),nflx2d(nlonp4,nlat)
!
! aureff = vertical profile for auroral heating efficiency (used by qjion)
! 
      real :: aureff(nlevp1)
      real :: bdriz(nlevp1)   ! vertical profile for background drizzle
!
! MEPED electron characteristic energy and hemespheric power:
      real :: echar_meped, epower_meped
!
! Subroutines internal to aurora module.
      contains
!-----------------------------------------------------------------------
#include <defs.h>
      subroutine aurora_cons(iprint,iyear,iday,secs)
!
! Set auroral constants (called from sub advance once per time-step).
!
      use params_module,only: dz,nlev
      use input_module,only:  ! from user input
     |  byimf,   ! BY component of IMF (nT)      (e.g., 0.)
     |  ctpoten, ! cross-cap potential (kV)      (e.g., 45.)
     |  power,   ! hemispheric power   (GW)      (e.g., 16.)
     |  solar_protons, ! 0/1 flag to get solar protons
     |  meped_file     ! if non-blank, use MEPED data (see meped.F)
!
! Args:
      integer,intent(in) :: iprint,iday,iyear
      real,intent(in) :: secs
!
! Local:
      real,parameter :: convert=3.1211e+8
      real :: alfa_1,alfa_2,plevel,roth,rote,rcp,rhp,
     |  alfa21,alfa22,e21,e22
      real :: arad(2),aurpower,hpgwsh
      real,parameter :: e=1.e-10
      logical :: found_protons
!
! Auroral heating efficiency (nlevp1) according to vertical resolution 
!
! 7/17/03: aureff and bdriz apparently not used (aureff in qjion is
!          local to qjion), so for now am not defining them for time-gcm.
!
! tiegcm: NLEV==28 for dz==0.5
!
#if (NLEV==28)
      aureff =
     |(/0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,
     |  0.55,0.50,0.45,0.37,0.30,0.20,0.16,0.13,0.10,0.08,0.07,0.06,
     |  0.05,0.05,0.05,0.05,0.05/)
      bdriz =
     |(/2.17E-4,4.48E-3,3.78E-2,1.55E-1,3.26E-1,3.85E-1,3.50E-1,2.81E-1,
     |  2.04E-1,1.38E-1,8.74E-2,5.23E-2,3.01E-2,1.68E-2,9.27E-3,5.08E-3,
     |  2.79E-3,1.55E-3,8.62E-4,4.85E-4,2.76E-4,1.58E-4,9.11E-5,5.28E-5,
     |  3.08E-5,1.79E-5,1.04E-5,6.01E-6,3.49E-6/)
!
! tiegcm: NLEV==56 for dz==0.25
!
#elif (NLEV==56)
        aureff =
     |  (/0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00,
     |    0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00,
     |    0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00,
     |    0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00,
     |    0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00,
     |    0.5250E+00, 0.5000E+00, 0.4750E+00, 0.4500E+00, 0.4100E+00,
     |    0.3700E+00, 0.3350E+00, 0.3000E+00, 0.2500E+00, 0.2000E+00,
     |    0.1800E+00, 0.1600E+00, 0.1450E+00, 0.1300E+00, 0.1150E+00,
     |    0.1000E+00, 0.9000E-01, 0.8000E-01, 0.7500E-01, 0.7000E-01,
     |    0.6500E-01, 0.6000E-01, 0.5500E-01, 0.5000E-01, 0.5000E-01,
     |    0.5000E-01, 0.5000E-01, 0.5000E-01, 0.5000E-01, 0.5000E-01,
     |    0.5000E-01, 0.5000E-01/)
        bdriz =
     |  (/0.2170E-03, 0.2349E-02, 0.4480E-02, 0.2114E-01, 0.3780E-01,
     |    0.9640E-01, 0.1550E+00, 0.2405E+00, 0.3260E+00, 0.3555E+00,
     |    0.3850E+00, 0.3675E+00, 0.3500E+00, 0.3155E+00, 0.2810E+00,
     |    0.2425E+00, 0.2040E+00, 0.1710E+00, 0.1380E+00, 0.1127E+00,
     |    0.8740E-01, 0.6985E-01, 0.5230E-01, 0.4120E-01, 0.3010E-01,
     |    0.2345E-01, 0.1680E-01, 0.1303E-01, 0.9270E-02, 0.7175E-02,
     |    0.5080E-02, 0.3935E-02, 0.2790E-02, 0.2170E-02, 0.1550E-02,
     |    0.1206E-02, 0.8620E-03, 0.6735E-03, 0.4850E-03, 0.3805E-03,
     |    0.2760E-03, 0.2170E-03, 0.1580E-03, 0.1246E-03, 0.9110E-04,
     |    0.7195E-04, 0.5280E-04, 0.4180E-04, 0.3080E-04, 0.2435E-04,
     |    0.1790E-04, 0.1415E-04, 0.1040E-04, 0.8205E-05, 0.6010E-05,
     |    0.4750E-05, 0.3490E-05/)
!
! time-gcm: NLEV==44 for dz==0.5
!
#elif (NLEV==44)

!
! time-gcm: NLEV==48 for dz==0.5 with ubc extension to +7
!
#elif (NLEV==48)

!
! time-gcm: NLEV==88 for dz==0.25
!

#elif (NLEV==96)
!
! time-gcm: NLEV==96 for dz==0.25 with ubc extension to +7
!

#else
        write(6,"('>>> aurora: unsupported NLEV=',i4,' dz=',f10.4)")
     |    nlev,dz
        call shutdown('aurora dz')
#endif
!
      plevel = 0.
! 8/1/06 btf: change conditional as per Richmond suggestion:
!     if (power >= 0.01) plevel = 2.09*alog(power)
      if (power >= 1.00) plevel = 2.09*alog(power)
!
! e1 = energy flux in the noon sector of the aurora (ergs/cm**2/s)
! e2 = energy flux in the midnight sector of the aurora (ergs/cm**2/s)
!
! e1 = energy flux in the noon sector of the aurora (ergs/cm**2/s)
! e2 = energy flux in the midnight sector of the aurora (ergs/cm**2/s)
! modified by LQIAN, 2007
! produce realistic oval compared to NOAA empirical auroral oval and TIMED/GUVI
! e1 formula given by Wenbin base on POLARVIS image;
! e2 formula based on Emery et al original auroral parameterization report
      e1 = max(0.50, -2.15 + 0.62*plevel)
      e2=1.+0.11*power
      e0 = 0.5 * (e1 + e2)
      ree = (e2 - e1) / (e1 + e2)
!
! 4/25/08 btf: Set ec and ed as function of power, as per Wenbin:
      ec=(0.24+0.0067*power)/5.
      ed=0.0012+0.0006*power
!
! h1 = Gaussian half-width of the noon auroral oval in degrees
! h2 = Gaussian half-width of the midnight auroral oval in degrees
! modified by LQIAN, 2007
! produce realistic oval compared to NOAA empirical auroral oval and TIMED/GUVI
! h1 formula given by Wenbin base on POLARVIS image;
! h2 formula based on Emery et al original auroral parameterization report
      h1 = min(2.35, 0.83 + 0.33*plevel)
      h2=2.5+0.025*max(power,55.)+0.01*min(0.,power-55.)
      rh = (h2 - h1) / (h1 + h2)
!
!     if (iamie==0) then
        theta0(isouth) = (-3.80+8.48*(ctpoten**0.1875))*dtr
        theta0(inorth) = (-3.80+8.48*(ctpoten**0.1875))*dtr
!     else
!       theta0(isouth) = (-1.92+8.10*(ctpoten**0.1875))*dtr
!       theta0(inorth) = (-1.92+8.10*(ctpoten**0.1875))*dtr
!     endif
      offa(isouth) = 1.0*dtr
      offa(inorth) = 1.0*dtr
      dskofa(isouth) = 0.
      dskofa(inorth) = 0.
!
      rhp = 14.20 + 0.96*plevel
      rcp = -0.43 + 9.69 * (ctpoten**0.1875)
      arad(isouth) = max(rcp,rhp)
      arad(inorth) = max(rcp,rhp)
      rrad(isouth) = arad(isouth)*dtr 
      rrad(inorth) = arad(inorth)*dtr
!
! 4/25/08 btf: New values for alfa_1 and alfa_2, as per Wenbin
      alfa_1 = 1.5
      alfa_2 = 2.
      alfa0 = 0.5*(alfa_1+alfa_2)
      ralfa = (alfa_2-alfa_1) / (alfa_1+alfa_2+e)
!     roth = (12.18 - 0.89 * plevel) * 15. ! tiegcm
!     rote = ( 2.62 - 0.55 * plevel) * 15. ! tiegcm
      roth = (12.18 - 0.89 * plevel)       ! tgcm24
      rote = ( 2.62 - 0.55 * plevel)       ! tgcm24
      rroth = roth * dtr
      rrote = rote * dtr
      h0 = 0.5 * (h1 + h2) * dtr
      alfa21 = 0.075
      alfa22 = 0.075
      alfa20 = 0.5 * (alfa21 + alfa22)
      ralfa2 = (alfa22 - alfa21) / (alfa21 + alfa22 + e)
 !    e21 = 1.e-80 ! tiegcm
 !    e22 = 1.e-80 ! tiegcm
      e21 = 1.e-20 ! timegcm
      e22 = 1.e-20 ! timegcm
      e20 = 0.5 * (e21 + e22)
      re2 = (e22 - e21) / (e21 + e22)
!
! Set cusp and drizzle parameters:
! (conversion between particle number density and characteristic
!  energy and column energy input)
!
      fc = convert * ec / alfac
      fd = convert * ed / alfad
!
! Solar proton flux.
! If user namelist solar_protons is set, get Jackman solar protons.
!   (see module protons.F)
! If model dates correspond to dates in which proton data is available,
!   characteristic energy and energy flux are returned by get_protons:
!   alfa_sp (alfad2) ! characteristic energy (MeV), interpolated from Eo
!   e_sp    (ed2)    ! energy flux (ergs/cm2/s), interpolated from Fe
! Note that alfa_sp and e_sp were set to constants in module data
!   above. If current model date is not within available proton data
!   periods, then get_protons returns silently, and the constant
!   values are used.
!
      alfa_sp = alfa_sp_default
      e_sp = e_sp_default
      if (solar_protons > 0) then
        call get_protons(iyear,iday,secs,istep,nstep,alfa_sp,e_sp,
     |    found_protons)
      endif
      flx_sp = e_sp/1.602e-6 ! tgcm24: fd = ed2/1.602E-6
!
! Get MEPED data:
      if (len_trim(meped_file) > 0) then
        call get_meped(iyear,iday,secs,
     |    epower_meped,echar_meped)
!       write(6,"('aurora: epower_meped=',e12.4,' echar_meped=',e12.4)")
!    |    epower_meped,echar_meped
!
! Use MEPED epower(GW) data in E30 calculation:
        hpgwsh = epower_meped/(1+arad(inorth)/arad(isouth))
        aurpower = hpgwsh/(1.+0.5*ree*rh*cos(rroth-rrote))/1.e+16
        e30 = aurpower/(2.*pi**1.5)*((6.37e8+1.e7)**2)*
     |        ((arad(isouth)*pi)/180.)*h0
        write(6,"(/,'aurora_cons: e1,e2=',2e12.4,' arad(south,north)=',
     |    2e12.4)") e1,e2,arad(isouth),arad(inorth)
        write(6,"('  h0=',e12.4,' ree=',e12.4,' rh=',
     |    e12.4,' rroth=',e12.4,' rrote=',e12.4)") h0,ree,rh,rroth,rrote
        write(6,"('  epower_meped=',e12.4,' hpgwsh=',e12.4,' aurpower=',
     |    e12.4,' e30=',e12.4)") epower_meped,hpgwsh,aurpower,e30
!
! Use MEPED characteristic energy (KeV) in alfa30:
        alfa30 = echar_meped
        write(6,"('aurora_cons: alfa30 = echar_meped = ',e12.4)")
     |    alfa30 
      endif
!
! The following parameters (offc through rr1) are used only in heelis 
!   potential calculation for the dynamo (see heelis.F)
!
      offc(isouth) = 1.*dtr
      offc(inorth) = 1.*dtr
      dskofc(isouth) = 0.
      dskofc(inorth) = 0.
      phid(isouth) = 0.
      phid(inorth) = 0.
      phin(isouth) = 180.*dtr
      phin(inorth) = 180.*dtr
      psim(:) =  0.50 * ctpoten * 1000.
      psie(:) = -0.50 * ctpoten * 1000.
!     pcen(isouth) = 0.
!     pcen(inorth) = 0.
      pcen(:) = 0. * 1000. ! tgcm24
      phidp0(:) = 90.*dtr
      phidm0(:) = 90.*dtr
      phinp0(:) = 90.*dtr
      phinm0(:) = 90.*dtr
      rr1(:) = -2.6
!
! Report to stdout:
      if (iprint > 0) then
        write(6,"(/,'aurora_cons:')")
        write(6,"('  cusp:    alfac=',f8.3,' ec=',f8.3,' fc=',e10.4)")
     |    alfac,ec,fc
        write(6,"('  drizzle: alfad=',f8.3,' ed=',f8.3,' fd=',e10.4)")
     |    alfad,ed,fd
        write(6,"('  half-widths = h1,h2=',2f10.3)") h1,h2
        write(6,"('  energy flux = e1,e2=',2f10.3)") e1,e2
        if (solar_protons > 0.and.found_protons) then
          write(6,"('  Time-dependent solar proton flux: alfa_sp=',
     |    f8.3,' e_sp = ',e10.4,' flx_sp = ',e10.4)") 
     |    alfa_sp,e_sp,flx_sp
        else
          write(6,"('  Constant solar proton flux: alfa_sp=',
     |    f8.3,' e_sp = ',e10.4,' flx_sp = ',e10.4)") 
     |    alfa_sp,e_sp,flx_sp
        endif
!       write(6,"('Aurora init: solar protons:')")
!       write(6,"('  alfa_sp=',e12.4)") alfa_sp
!       write(6,"('  e_sp   =',e12.4)") e_sp   
!       write(6,"('  flx_sp =',e12.4)") flx_sp   
        if (add_helectron)
     |    write(6,"('  high-energy electrons: alfa30 = ',f8.3,
     |    ' e30 = ',f8.3)") alfa30,e30
        if (len_trim(meped_file) > 0)
     |    write(6,"('  MEPED epower=',f10.2,' echar=',f10.2)")
     |      epower_meped,echar_meped
        write(6,"(' ')")
      endif
      end subroutine aurora_cons
!-----------------------------------------------------------------------
      subroutine aurora(tn,o2,o1,barm,lev0,lev1,lon0,lon1,lat)
!
! Auroral parameterization driver.
! (called by sub dynamics, from inside parallel latitude scan)
!
      use magfield_module,only: 
     |  rlatm,rlonm, ! magnetic latitude,longitude grid (radians)
     |  sunlons      ! sun's mag longitudes (see sub sunloc, magfield.F) 
      use params_module,only: nlat,zibot,dz
!
      use qrj_module,only:      ! for addfld debug only
     |  qo2p,  ! o2+ ionization
     |  qop,   ! o+  ionization
     |  qn2p,  ! n2+ ionization
     |  qnp    ! n+  ionization
!
! Args:
      integer,intent(in) :: 
     |  lev0,lev1,  ! bottom,top of pressure column (1,nlevp1) 
     |  lon0,lon1,  ! first,last longitude indices of current domain
     |  lat         ! current latitude index 
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  tn,    ! neutral gas temperature (deg K)
     |  o2,    ! O2 mass mixing ratio
     |  o1,    ! O mass mixing ratio
     |  barm   ! mean molecular weight (see sub addiag)
!
! Local:
      integer :: i,ier
      integer,parameter :: nlond2 = max(1,nlon/2)
      real :: aurlat  ! boundary lat below which auroral is not calculated
                      ! (32.5 degrees at 5 degree resolution)
      real :: ofda,cosofa,sinofa,aslona
!
#ifdef VT
!     code = 130 ; state = 'aurora' ; activity='ModelCode'
      call vtbegin(130,ier)
#endif
!
! Save ionization rates on entry to secondary history:
!     call addfld('QO2P_PRE',' ',' ',qo2p(lev0:lev1,lon0:lon1,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('QOP_PRE' ,' ',' ',qop(lev0:lev1,lon0:lon1,lat) ,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('QN2P_PRE',' ',' ',qn2p(lev0:lev1,lon0:lon1,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Check latitude, and return if below aurlat (32.5 deg at 5 deg res):
!
      aurlat = (float(lat-nlat)-0.5)*dphi+pi/2.
      do_aurora = .true.
      if (abs(aurlat) < pi/6.) do_aurora = .false.
      if (.not.do_aurora) then
#ifdef VT
        call vtend(130,ier)
#endif
        return
      endif
!
      do i=1,nlon
        dlat_aur(i) = rlatm(i+2,lat)
        dlon_aur(i) = rlonm(i+2,lat)-sunlons(lat)
      enddo
!
! Establish which hemisphere current latitude is in:
! (ihem will be either isouth (1) or inorth (2))
      ihem = int(dlat_aur(nlond2)*2./pi+2.)
!
! Find auroral circle coordinates (formerly in flowv3):
!
      ofda = sqrt(offa(ihem)**2+dskofa(ihem)**2)
      cosofa = cos(ofda)
      sinofa = sin(ofda)
      aslona = asin(dskofa(ihem)/ofda)
      do i=1,nlon
        sinlat(i) = sin(abs(dlat_aur(i)))
        coslat(i) = cos(dlat_aur(i))
        sinlon(i) = sin(dlon_aur(i)+aslona)
        coslon(i) = cos(dlon_aur(i)+aslona)
        colat(i) = acos(cosofa*sinlat(i)-sinofa*coslat(i)*coslon(i))
        alon(i) = amod(atan2(sinlon(i)*coslat(i),sinlat(i)*sinofa+
     |    cosofa*coslat(i)*coslon(i))-aslona+3.*pi,(pi*2.))-pi
      enddo
!
!     write(6,"('aurora: lat=',i2,' sinlat=',/,(6e12.4))") lat,sinlat
!     write(6,"('aurora: lat=',i2,' coslat=',/,(6e12.4))") lat,coslat
!     write(6,"('aurora: lat=',i2,' coslon=',/,(6e12.4))") lat,coslon
!     write(6,"('aurora: lat=',i2,' colat=',/,(6e12.4))") lat,colat
!     write(6,"('aurora: lat=',i2,' alon =',/,(6e12.4))") lat,alon
!
! Make cusp:
      call aurora_cusp(lat)
!
! Make alfa, flux, and drizzle:
      call aurora_heat(lat)

! TEMP
! Save aurora oval on exit to full arrays so can save later to secondary history
      do i=lon0,lon1
        cusp2d(i,lat)=cusp(i)
        drzl2d(i,lat)=drizl(i)
        alfa2d(i,lat)=alfa(i)
        nflx2d(i,lat)=flux(i)
      enddo
!
! Auroral additions to ionization rates:
! Sub aurora_ions is a standalone subroutine, i.e., not in aurora_module.
! Pass module data to aurora_ions through argument list.
!
!     subroutine aurora_ions(drizl,cusp,alfa1,alfa2,alfa3,flux1,flux2,
!    |  flux3,tn,o2,o1,barm,zpbot,dzp,lev0,lev1,lon0,lon1,lat)
!
      call aurora_ions(
     |  drizl(lon0:lon1), cusp (lon0:lon1),
     |  alfa (lon0:lon1), alfa2(lon0:lon1), alfa3(lon0:lon1),
     |  flux (lon0:lon1), flux2(lon0:lon1), flux3(lon0:lon1),
     |  tn,o2,o1,barm,zibot,dz,lev0,lev1,lon0,lon1,lat)
!
#ifdef VT
!     code = 130 ; state = 'aurora' ; activity='ModelCode'
      call vtend(130,ier)
#endif
!
      end subroutine aurora
!-----------------------------------------------------------------------
      subroutine aurora_cusp(lat)
      use params_module,only: nlonp2
      use cons_module,only: check_exp
!
! Calculate horizontal variation of polar cusp heating (was polht.F)
!
! Args:
      integer,intent(in) :: lat
!
! Local:
      integer :: i
      real,parameter :: 
     |  s5=.08726646, 
     |  s20=.34906585,
     |  pi_cusp = 3.1415927 ! hard-wired by tgcm15
!
! expo() (util.F) is used only if check_exp is true. This will avoid
! NaNS fpe, but will degrade performance. Check_exp is in cons.F.
!
      real,external :: expo
!
      do i=1,nlon
        if (.not.check_exp) then
          cusp(i+2) = (exp(-((theta0(ihem)-colat(i))/s5)**2)+
     |                 exp(-((pi_cusp-theta0(ihem)-colat(i))/s5)**2))*
     |                 exp(-(atan2(sin(alon(i)-phid(ihem)),
     |                             cos(alon(i)-phid(ihem)))/s20)**2)
        else
          cusp(i+2) = (expo(-((theta0(ihem)-colat(i))/s5)**2)+
     |                expo(-((pi_cusp-theta0(ihem)-colat(i))/s5)**2))*
     |                 expo(-(atan2(sin(alon(i)-phid(ihem)),
     |                             cos(alon(i)-phid(ihem)))/s20)**2)
        endif
      enddo ! i=3,nlonp2
!
! Periodic points:
      do i=1,2
        cusp(i) = cusp(nlon+i)     ! 1:2 <- 73:74 at 5 deg res
        cusp(nlonp2+i) = cusp(i+2) ! 75:76 <- 3:4
      enddo
!     write(6,"('aurora_cusp: lat=',i2,' cusp=',/,(6e12.4))") 
!    |  lat,cusp(3:nlon+2)
      end subroutine aurora_cusp 
!-----------------------------------------------------------------------
      subroutine aurora_heat(lat)
!
! Calculate alfa, flux, and drizzle.
!
      use params_module,only: nlon,nlonp2
!
! Args:
      integer,intent(in) :: lat
!
! Local:
      integer :: i
      real,dimension(nlon) :: 
     |  coslamda,  ! cos(angle from throat)
     |  halfwidth, ! oval half-width
     |  dtheta     ! latitudinal variation (Gaussian)
!
! Low-energy protons:
!     alfap0 = 0.5*(alfap1+alfap2)
!     e0p = 0.5*(pe1+pe2)
!
! coslamda  = cos(lamda)
! halfwidth = auroral half width
! dtheta    = colat-theta0(ihem)
! alfa      = electron energy
!
      if (alfa0 > 0.01) then
        do i=1,nlon
          coslamda(i) = cos(atan2(sin(alon(i)-rrote),
     |                            cos(alon(i)-rrote)))
!
! Auroral oval half-width (equation (1) in Roble,1987):
          halfwidth(i)=h0*(1.-rh*cos(atan2(sin(alon(i)-rroth),
     |      cos(alon(i)-rroth))))
          dtheta(i) = colat(i)-rrad(ihem)
!
! Characteristic energy (equation (2) in Roble,1987):
          alfa(i+2) =  alfa0*(1.-ralfa*coslamda(i))
        enddo ! i=1,nlon
      else
        do i=1,nlon
          coslamda(i)=cos(atan2(sin(alon(i)-rrote),cos(alon(i)-rrote))) ! wk1
          halfwidth(i)=h0*(1.-rh*cos(atan2(sin(alon(i)-rroth),          ! wk2
     |      cos(alon(i)-rroth))))
          dtheta(i) = colat(i)-rrad(ihem)                               ! wk3
!
! 1/02 btf: this alfa not currently in use:
!         alfa(i+2) = alfk+
!    |      alf6*exp(-((dtheta(i)-rd6-rd6v*cos(alon(i)) )/rh6)**2)*
!    |      exp(-(atan2(sin(alon(i)-rrot6),cos(alon(i)-rrot6))/rt6)**2)+
!    |      alf21*exp(-((dtheta(i)-rd21)/rh21)**2)*
!    |      exp(-(atan2(sin(alon(i)-rrot21),cos(alon(i)-rrot21))/
!    |      rt21)**2)
        enddo ! i=1,nlon
      endif
!     write(6,"('aurora_heat: lat=',i3,' alfa=',/,(6e12.4))")
!    |  lat,alfa(3:nlon+2)
!
! Flux, drizzle, alfa2, flux2:
      do i=1,nlon
!
! Energy flux (equation (3) in Roble,1987):
        flux(i+2) = e0*(1.-ree*coslamda(i))*exp(-(dtheta(i)/
     |    halfwidth(i))**2) / (2.*alfa(i+2)*1.602e-9)
        drizl(i+2) = exp(-((dtheta(i)+abs(dtheta(i)))/(2.*h0))**2)
        alfa2(i+2) = alfa20*(1.-ralfa2*coslamda(i))
        flux2(i+2) = e20*(1.-re2*coslamda(i))*exp(-(dtheta(i)/
     |    halfwidth(i))**2) / (2.*alfa2(i+2)*1.602E-9)
!
! Alfa3, flux3 for high energy electrons:
        alfa3(i+2) = alfa30
	flux3(i+2) = e30*exp(-(dtheta(i)/halfwidth(i))**2)
     |                      / (2.*alfa3(i+2)*1.602E-9)
!
! For electron temperature (used in settei):  
        qteaur(i+2,lat) = -7.e+8*exp(-(dtheta(i)/halfwidth(i))**2)
      enddo ! i=1,nlon
!
! Periodic points:
      do i=1,2
        alfa  (i) = alfa (nlon+i)
        flux  (i) = flux (nlon+i)
        drizl (i) = drizl(nlon+i)
        alfa2 (i) = alfa2(nlon+i)
        alfa3 (i) = alfa3(nlon+i)
        flux2 (i) = flux2(nlon+i)
        flux3 (i) = flux3(nlon+i)
        qteaur(i,lat)= qteaur(nlon+i,lat)
        alfa  (nlonp2+i) = alfa (i+2)
        flux  (nlonp2+i) = flux (i+2)
        drizl (nlonp2+i) = drizl(i+2)
        alfa2 (nlonp2+i) = alfa2(i+2)
        alfa3 (nlonp2+i) = alfa3(i+2)
        flux2 (nlonp2+i) = flux2(i+2)
        flux3 (nlonp2+i) = flux3(i+2)
        qteaur(nlonp2+i,lat)= qteaur(i+2,lat)
      enddo
      end subroutine aurora_heat
!-----------------------------------------------------------------------
      end module aurora_module
!-----------------------------------------------------------------------
      subroutine aurora_ions(drizl,cusp,alfa1,alfa2,alfa3,flux1,flux2,
     |  flux3,tn,o2,o1,barm,zpbot,dzp,lev0,lev1,lon0,lon1,lat)
!
! Calculate auroral additions to ionization rates.
! Note this routine is standalone, not in aurora_module.
!
      use cons_module,only: ! see cons.F
     |  p0,    ! standard pressure (5.0e-4)
     |  expz,  ! exp(-zp) where zp is log pressure (ln(p0/p))
     |  grav,  ! gravity
     |  gask,  ! gas constant (8.314e7)
     |  boltz, ! boltzman's constant (1.38e-16)
     |  rmass_o2, rmass_o1, rmass_n2,        ! molelcular weights of o2,o,n2
     |  rmassinv_o2,rmassinv_o1,rmassinv_n2, ! 1./rmass
     |  avo    ! avogadro number (6.023e23)
!
! Total ionizations. Auroral ionization is added to these:
      use qrj_module,only:  ! (nlevp1,nlonp4) see qrj.F
     |  qo2p,  ! o2+ ionization
     |  qop,   ! o+  ionization
     |  qn2p,  ! n2+ ionization
     |  qnp    ! n+  ionization
      use chemrates_module,only: disn2p
!
! Auroral parameters:
      use aurora_module,only: alfac, alfad, alfa_sp, ec,ed, fc,fd,
     |  flx_sp, e30, add_helectron
!    |  flx_sp, e30, add_sproton, add_helectron
      implicit none
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,intent(in),dimension(lon0:lon1) :: 
     |  drizl,cusp,alfa1,alfa2,alfa3,flux1,flux2,flux3
      real,intent(in) :: zpbot,dzp
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  tn,o2,o1,barm
!
! Local:
      integer :: i,k,klev
      real :: zlev,xn2
      real,dimension(lev0:lev1,lon0:lon1) :: 
     |  p0ez, xalfa1, xalfa2, xcusp, xdrizl,       ! input to sub aion
     |  xalfa_sp, xalfa3,
     |  flux1_ion, flux2_ion, cusp_ion, drizl_ion, ! output from sub aion
     |  alfa1_ion, alfa2_ion, alfa3_ion,           ! output from sub aion
     |  alfasp_bion,                               ! output from sub bion
     |  barm_t, qsum, denom, p0ez_mbar, tk_mbar,
     |  qo2p_aur, qop_aur, qn2p_aur ! auroral ionization for O2+, O+, N2+
!
! falfa1,...,falfa3 are 2-d only to exactly match tgcm24 results.
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  falfa1, falfa2, fcusp, fdrizl, falfa_sp, falfa3
!     real,dimension(lon0:lon1) ::
!    |  falfa1, falfa2, fcusp, fdrizl, falfa_sp, falfa3
      real :: qia(5) ! low energy proton source (not in use, 1/02)
!
      qia(:) = 0.
!
!     write(6,"('aurora_ions: lat=',i2,' alfa1=',/,(6e12.4))")
!    |  lat,alfa1
!     write(6,"('aurora_ions: lat=',i2,' alfa2=',/,(6e12.4))")
!    |  lat,alfa2
!     write(6,"('aurora_ions: lat=',i2,' cusp =',/,(6e12.4))")
!    |  lat,cusp 
!     write(6,"('aurora_ions: lat=',i2,' drizl=',/,(6e12.4))")
!    |  lat,drizl 
!
      do i=lon0,lon1
        do k=lev0,lev1
          p0ez(k,i) = (p0*expz(k)/(grav*4.e-6))**0.606  ! s6
          xalfa1(k,i) = p0ez(k,i)/alfa1(i)  ! s7 = s6 / t6
          xalfa2(k,i) = p0ez(k,i)/alfa2(i)  ! s8 = s6 / t7
          xalfa3(k,i) = p0ez(k,i)/alfa3(i)
          xcusp (k,i) = p0ez(k,i)/alfac
          xdrizl(k,i) = p0ez(k,i)/alfad
!
          p0ez_mbar(k,i) = p0*expz(k)/(boltz*tn(k,i)*
     |      (o2(k,i)*rmassinv_o2+
     |       o1(k,i)*rmassinv_o1+
     |       (1.-o2(k,i)-o1(k,i))*rmassinv_n2))/avo
          tk_mbar(k,i) = gask*tn(k,i)/(grav*barm(k,i))*p0ez_mbar(k,i)
          xalfa_sp(k,i) = ((tk_mbar(k,i)/0.00271)**0.58140)/alfa_sp ! s12
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1
!
!     call addfld('P0EZ'  ,' ',' ',p0ez  ,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('XALFA1',' ',' ',xalfa1,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('XALFA2',' ',' ',xalfa2,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('XALFA3',' ',' ',xalfa3,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('XCUSP' ,' ',' ',xcusp ,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('XDRIZL',' ',' ',xdrizl,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Init (whole array operations):
      flux1_ion = 1.e-20
      flux2_ion = 1.e-20
      alfa1_ion = 1.e-20
      alfa2_ion = 1.e-20
      cusp_ion  = 1.e-20
      drizl_ion = 1.e-20
      alfasp_bion = 1.e-20
      alfa3_ion   = 1.e-20
!
! Sub aion calculates auroral electrons: 
! (aion call with xalfa2,alfa2_ion not in tgcm24 (alfa2_ion (s4) is set 1.e-20)
!
      zlev = -7.5 ! time-gcm only
      klev = int((zlev-zpbot)/dzp)+1
      call aion(xalfa1,alfa1_ion,klev,lev0,lev1,lon0,lon1,lat) ! s7,s3
! 
      zlev = -5.0
      klev = int((zlev-zpbot)/dzp)+1
      call aion(xcusp ,cusp_ion ,klev,lev0,lev1,lon0,lon1,lat) ! s9,s5
!
      zlev = -7.5
      klev = int((zlev-zpbot)/dzp)+1
      call aion(xdrizl,drizl_ion,klev,lev0,lev1,lon0,lon1,lat) ! s10,s6
!
! Solar protons:
      klev = 1
      if (flx_sp > 1.e-19)
     |  call bion(xalfa_sp,alfasp_bion,klev,lev0,lev1,lon0,lon1,lat)
!
! High energy electrons:
      if (e30 > 1.e-19) 
     |  call aion(xalfa3, alfa3_ion,klev,lev0,lev1,lon0,lon1,lat)
!
! Fields needed (k,i) for qsum:
!     call addfld('ALFA1' ,' ',' ',alfa1_ion,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('ALFA2' ,' ',' ',alfa2_ion,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('ALFA3' ,' ',' ',alfa3_ion,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('XALFA2',' ',' ',xalfa2   ,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('FALFA2',' ',' ',falfa2   ,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('CUSP'  ,' ',' ',cusp_ion ,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('DRIZL' ,' ',' ',drizl_ion,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('ALFASPB',' ',' ',alfasp_bion,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          falfa1(k,i) = alfa1(i)*flux1(i)  ! s7=t6*t8
          falfa2(k,i) = alfa2(i)*flux2(i)  ! s8=t7*t9
          fcusp (k,i) = cusp(i)*alfac*fc   ! s9=t5*alfac*fc
          fdrizl(k,i) = drizl(i)*alfad*fd  ! s10
          falfa3(k,i) = alfa3(i)*flux3(i)  ! s13 (high energy electrons)
!
! Note if namelist parameter solar_protons was set, then flx_sp is
! time-dependent from sub get_solar_protons.
! 
          falfa_sp(k,i) = drizl(i)*p0ez_mbar(k,i)*flx_sp*1.e6/  ! s12
     |      (tk_mbar(k,i)*35.)
        enddo
      enddo
!
! Pressure column, bottom to top:
      do i=lon0,lon1
        do k=lev0,lev1-1
          barm_t(k,i) = grav*0.5*(barm(k,i)+barm(k+1,i))/(35.e-3*gask* ! s11
     |      tn(k,i))

!
! Small diffs ~ tgcm24 in qsum:
!         qsum(k,i) =                  ! s1
!    |      falfa1(i)*alfa1_ion(k,i)+  ! s7*s3
!    |      falfa2(i)*alfa2_ion(k,i)+  ! s8*s4 (s4 = alfa2_ion = 1.e-20)
!    |      fcusp (i)*cusp_ion (k,i)+  ! s9*s5
!    |      fdrizl(i)*drizl_ion(k,i)   ! s10*s6
!
! Include solar protons if add_sproton is set, and high energy electrons
! if add_helectron is set:
!         if (add_helectron)                                     ! true in tgcm24
!    |      qsum(k,i) = qsum(k,i) + falfa3(i)*alfa3_ion(k,i)     ! s13*s15
!         if (add_sproton)                                       ! true in tgcm24
!    |      qsum(k,i) = qsum(k,i) + falfa_sp(i)*alfasp_bion(k,i) ! s12*s14
!         qsum(k,i) = qsum(k,i)*barm_t(k,i)                      ! s1 = s1*s11

! Try making qsum exactly as in tgcm24 to avoid small diffs:
!       S1(I,1)=S7(I,1)*S3(I,1)+S8(I,1)*S4(I,1)
!       S1(I,1)=S1(I,1)+S13(I,1)*S15(I,1)
!       S1(I,1)=S1(I,1)+S9(I,1)*S5(I,1)
!       S1(I,1)=(S1(I,1)+S10(I,1)*S6(I,1))*S11(I,1)+S12(I,1)*S14(I,1)

!
! 8/29/03: qsum diffs with tgcm24: max 4.25e-20 from bottom boundary up to 
!          about zp -12:
! 5/24/05: Changed xalfa2 to falfa2 in first qsum statement:
!
          qsum(k,i)= falfa1(k,i)*alfa1_ion(k,i)+falfa2(k,i)*
     |      alfa2_ion(k,i)
          qsum(k,i) = qsum(k,i)+falfa3(k,i)*alfa3_ion(k,i)
          qsum(k,i) = qsum(k,i)+fcusp(k,i)*cusp_ion(k,i)
          qsum(k,i) = (qsum(k,i)+fdrizl(k,i)*drizl_ion(k,i))*barm_t(k,i)
     |      +falfa_sp(k,i)*alfasp_bion(k,i)
!
! Denominator of equations (13-16) in Roble,1987.
          denom(k,i) =                                                   ! s3
     |      0.92*(1.-o2(k,i)-o1(k,i))*rmassinv_n2 +
     |      1.5             *o2(k,i) *rmassinv_o2 + 
     |      0.56            *o1(k,i) *rmassinv_o1
!
! Production of O2+ (equation (15) in Roble,1987):
          qo2p_aur(k,i) = qsum(k,i)*o2(k,i)/(rmass_o2*denom(k,i))+qia(2) ! s13
!
! Production of O+ (equation (16) in Roble,1987):
          qop_aur(k,i) = qsum(k,i)*                                      ! s14
     |      (0.5 *o2(k,i)*rmassinv_o2+
     |       0.56*o1(k,i)*rmassinv_o1)/denom(k,i) + qia(3)
!
! Production of N2+ (equation (13) in Roble,1987)
          xn2 = (1.-o2(k,i)-o1(k,i))
!         if (xn2 <= 1.e-8) xn2 = 1.e-8 ! tiegcm
          if (xn2 <= 1.e-6) xn2 = 1.e-6 ! timegcm
          qn2p_aur(k,i) = qsum(k,i)*0.7*xn2/(rmass_n2*denom(k,i))+qia(1) ! s15
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1

!     real,intent(in),dimension(lon0:lon1) :: 
!    |  drizl,cusp,alfa1,alfa2,alfa3,flux1,flux2,flux3
!     real,dimension(lev0:lev1,lon0:lon1) ::
!    |  falfa1, falfa2, fcusp, fdrizl, falfa_sp, falfa3

!     write(6,"('aurora: lat=',i3,' flux1='  ,/,(3e28.20))") lat,flux1
!     write(6,"('aurora: lat=',i3,' alfa1='  ,/,(3e28.20))") lat,alfa1
!     write(6,"('aurora: lat=',i3,' alfa2='  ,/,(3e28.20))") lat,alfa2
!     write(6,"('aurora: lat=',i3,' falfa1=' ,/,(3e28.20))") 
!    |  lat,falfa1(1,:)
!     write(6,"('aurora: lat=',i3,' falfa2=' ,/,(3e28.20))") 
!    |  lat,falfa2(1,:)
!     write(6,"('aurora: lat=',i3,' falfa3=' ,/,(3e28.20))") 
!    |  lat,falfa3(1,:)
!     write(6,"('aurora: lat=',i3,' falfasp=',/,(3e28.20))")
!    |  lat,falfa_sp(1,:)
!     write(6,"('aurora: lat=',i3,' fcusp='  ,/,(3e28.20))") 
!    |  lat,fcusp(1,:)
!     write(6,"('aurora: lat=',i3,' fdrizl=' ,/,(3e28.20))") 
!    |  lat,fdrizl(1,:)

!     call addfld('BARMT'   ,' ',' ',barm_t  ,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('QSUM'    ,' ',' ',qsum    ,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('DENOM'   ,' ',' ',denom   ,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('QO2P_AUR',' ',' ',qo2p_aur,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('QOP_AUR' ,' ',' ',qop_aur ,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('QN2P_AUR',' ',' ',qn2p_aur,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('QO2P_PRE',' ',' ',qo2p(:,:,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('QOP_PRE' ,' ',' ',qop (:,:,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('QN2P_PRE',' ',' ',qn2p(:,:,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Add auroral to total ionization rates for o2+, o+, and n2+
! (see qrj_module)
!
      do i=lon0,lon1
        do k=lev0+1,lev1-1
          qo2p(k,i,lat) = qo2p(k,i,lat)+sqrt(qo2p_aur(k,i)*
     |      qo2p_aur(k-1,i))
          qop (k,i,lat) = qop (k,i,lat)+sqrt(qop_aur (k,i)*
     |      qop_aur(k-1,i))
          qn2p(k,i,lat) = qn2p(k,i,lat)+sqrt(qn2p_aur(k,i)*
     |      qn2p_aur(k-1,i))
          qnp (k,i,lat) = qnp (k,i,lat)+
     |      .22/.7 * sqrt(qn2p_aur(k,i)*qn2p_aur(k-1,i))
          disn2p(k,i,lat) = disn2p(k,i,lat)+sqrt(qn2p_aur(k,i)*
     |      qn2p_aur(k-1,i))
        enddo ! k=lev0+2,lev1-1
!
! Extrapolate lower boundary:
        qo2p(lev0,i,lat) = qo2p(lev0,i,lat)+1.5*qo2p_aur(lev0,i)-
     |    0.5*qo2p_aur(lev0+1,i)        
        qop (lev0,i,lat) = qop (lev0,i,lat)+1.5*qop_aur (lev0,i)-
     |    0.5*qop_aur (lev0+1,i)        
        qn2p(lev0,i,lat) = qn2p(lev0,i,lat)+1.5*qn2p_aur(lev0,i)-
     |    0.5*qn2p_aur(lev0+1,i)        
        qnp (lev0,i,lat) = qnp (lev0,i,lat)+ .22/.7 * 
     |    (1.5*qn2p_aur(lev0,i)-0.5*qn2p_aur(lev0+1,i))
        disn2p(1,i,lat) = disn2p(lev0,i,lat)+1.5*qn2p_aur(lev0,i)-
     |    0.5*qn2p_aur(lev0+1,i)        
!
! Extrapolate upper boundary:
        qo2p(lev1,i,lat)= qo2p(lev1,i,lat)+1.5*qo2p_aur(lev1-1,i)-
     |    0.5*qo2p_aur(lev1-2,i)
        qop (lev1,i,lat)= qop (lev1,i,lat)+1.5*qop_aur (lev1-1,i)-
     |    0.5*qop_aur (lev1-2,i)
        qn2p(lev1,i,lat)= qn2p(lev1,i,lat)+1.5*qn2p_aur(lev1-1,i)-
     |    0.5*qn2p_aur(lev1-2,i)
        qnp (lev1,i,lat)= qnp (lev1,i,lat)+ .22/.7 * 
     |    (1.5*qn2p_aur(lev1-1,i)-0.5*qn2p_aur(lev1-2,i))
        disn2p(lev1,i,lat) = disn2p(lev1,i,lat)+
     |    1.5*qn2p_aur(lev1-1,i)-0.5*qn2p_aur(lev1-2,i)
      enddo ! i=lon0,lon1

!     call addfld('QO2P',' ',' ',qo2p(lev0:lev1,lon0:lon1,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('QOP' ,' ',' ',qop (lev0:lev1,lon0:lon1,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('QN2P',' ',' ',qn2p(lev0:lev1,lon0:lon1,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('QNP' ,' ',' ',qnp (lev0:lev1,lon0:lon1,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('DISN2P',' ',' ',disn2p(lev0:lev1,lon0:lon1,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)

      end subroutine aurora_ions
!-----------------------------------------------------------------------
      subroutine aion(si,so,klev,lev0,lev1,lon0,lon1,lat)
      implicit none
!
! Calculates integrated f(x) needed for total auroral ionization.
! See equations (10-12) in Roble,1987.
!
! Args:
      integer,intent(in) :: klev,lev0,lev1,lon0,lon1,lat
      real,intent(in)  :: si(lev0:lev1,lon0:lon1)
      real,intent(out) :: so(lev0:lev1,lon0:lon1)
!
! Local:
      integer :: k,i
      real :: xlog
!
! Coefficients for equation (12) of Roble,1987, giving total 
! ionization rate (revised since 1987):
      real,parameter :: cc(8) =
     | (/3.2333134511131 ,  2.5658873458085 ,  2.2540957232641 ,
     |   0.72971983372673,  1.1069072431948 ,  1.7134937681128 ,
     |   1.8835442312993 ,  0.86472135072090/)
!
! Use the identity x**y = exp(y*ln(x)) for performance 
! (fewer (1/2) trancendental functions are required).
!
      do i=lon0,lon1
        do k=klev,lev1-1
          xlog = log(si(k,i))
          so(k,i) = cc(1)*exp(cc(2)*xlog-cc(3)*exp(cc(4)*xlog))+
     |              cc(5)*exp(cc(6)*xlog-cc(7)*exp(cc(8)*xlog))
        enddo ! k=klev,lev1
      enddo ! i=lon0,lon1
      end subroutine aion
!-----------------------------------------------------------------------
      subroutine bion(si,so,klev,lev0,lev1,lon0,lon1,lat)
      implicit none
!
! Calculates integrated f(x) needed for total auroral ionization.
! See equations (10-12) in Roble,1987.
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat,klev
      real,intent(in)  :: si(lev0:lev1,lon0:lon1)
      real,intent(out) :: so(lev0:lev1,lon0:lon1)
!
! Local:
      integer :: k,i
      real :: xlog
      real,parameter :: cc(8) =
     |  (/0.12718, 4.9119, 1.8429, 0.99336, 0.52472, 1.5565,
     |    0.85732, 1.4116/)
!
! Use the identity x**y = exp(y*ln(x)) for performance 
! (fewer (1/2) trancendental functions are required).
!
      do i=lon0,lon1
        do k=klev,lev1-1
          xlog = log(si(k,i))
          so(k,i) = cc(1)*exp(cc(2)*xlog-cc(3)*exp(cc(4)*xlog))+
     |              cc(5)*exp(cc(6)*xlog-cc(7)*exp(cc(8)*xlog))
          enddo ! i=lon0,lon1
      enddo ! k=lev0,lev1
      end subroutine bion