
      module mo_aurora
!-----------------------------------------------------------------------
!
! 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 mo_aurora, mo_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 shr_kind_mod,  only: r8 => shr_kind_r8
      use mo_constants,  only: pi, &
                               avo => avogadro, &
                               boltz_cgs, &
                               gask => rgas_cgs
      use cam_logfile,   only: iulog
      use spmd_utils,    only: masterproc

      implicit none

      interface aurora
         module procedure aurora_prod
         module procedure aurora_hrate
      end interface

      save

      integer, parameter  :: isouth = 1
      integer, parameter  :: inorth = 2

      ! g = 8.7 m/s^2? Because this is 400 km up?
      real(r8), parameter :: grav   = 870._r8          ! (cm/s^2)

      integer  :: lev1 = 1
      real(r8) :: twopi
      real(r8) :: rmass_o1
      real(r8) :: rmass_o2
      real(r8) :: rmass_n2
      real(r8) :: rmassinv_o1
      real(r8) :: rmassinv_o2
      real(r8) :: rmassinv_n2
      real(r8) :: dtr

!-----------------------------------------------------------------------
! 	... 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)
!-----------------------------------------------------------------------
      real(r8), parameter :: alfad = 2.0_r8, &
                             ed    = 0.5_r8     
      real(r8) :: 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)
!-----------------------------------------------------------------------
      real(r8), parameter :: alfac = 1.0_r8, &
                             ec    = 0.5_r8
      real(r8) :: 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(r8) :: &
        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)
! Add solar protons to ionization if add_sproton is true (time-gcm only)
!-----------------------------------------------------------------------
      logical :: add_sproton = .false.
      real(r8), parameter :: &
        alfa_sp = 10._r8, &
        e_sp    = 1.e-20_r8
      real(r8) :: flx_sp

!-----------------------------------------------------------------------
! 	... high energy electron parameters in the auroral oval (time-gcm only):
!-----------------------------------------------------------------------
      logical :: add_helectron = .false.
      real(r8), parameter :: alfa30 = 40._r8, &                  ! Characteristic energy of auroral electrons (Mev)
                             e30    = .05_r8                     ! Column energy of auroral electrons (ergs/cm**2/s)

!-----------------------------------------------------------------------
! 	... additional auroral parameters
!-----------------------------------------------------------------------
      real(r8) :: &
        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
      real(r8) :: &
        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))
        rrad(2)         ! radius of auroral circle in radians
      real(r8) :: ctpoten    ! cross-cap potential (kV)
      real(r8) :: byimf      ! BY component of IMF (nT)


      private
      public :: aurora_inti, aurora_timestep_init, aurora
      public :: aurora_register

      logical :: has_ions = .false.
      integer :: indxAIPRS = -1

      contains

        
      !----------------------------------------------------------------------
      !----------------------------------------------------------------------
      subroutine aurora_register
        use ppgrid,       only : pver,pcols
        use physics_buffer, only : pbuf_add_field, dtype_r8

        ! add ionization rates to phys buffer for waccmx ionosphere module

        call pbuf_add_field('AurIPRateSum' , 'physpkg', dtype_r8, (/pcols,pver/), indxAIPRS)     ! Sum of ion auroral production rates for O2

      endsubroutine aurora_register

      subroutine aurora_inti
!-----------------------------------------------------------------------
! 	... initialize aurora module
!-----------------------------------------------------------------------

      use ppgrid,       only : pver
      use pmgrid,       only : plev, plevp
      use constituents, only : cnst_get_ind, cnst_mw
      use chem_mods,    only : adv_mass
      use ref_pres,     only : pref_mid
      use mo_chem_utls, only : get_spc_ndx
      use cam_history,  only : addfld

      implicit none

!-----------------------------------------------------------------------
! 	... local variables
!-----------------------------------------------------------------------
      integer             :: k, m
      real(r8), parameter :: e       = 1.e-10_r8
      real(r8), parameter :: convert = 3.1211e8_r8

      real(r8) :: plb
      real(r8) :: alfa_1, alfa_2, alfa21, alfa22
      real(r8) :: e21, e22

      integer :: op_ndx,o2p_ndx,np_ndx,n2p_ndx,e_ndx

      op_ndx   = get_spc_ndx( 'Op' )
      o2p_ndx  = get_spc_ndx( 'O2p' )
      np_ndx   = get_spc_ndx( 'Np' )
      n2p_ndx  = get_spc_ndx( 'N2p' )
      e_ndx    = get_spc_ndx( 'e' )

      has_ions = op_ndx > 0 .and. o2p_ndx > 0 .and. np_ndx > 0 .and. n2p_ndx > 0 .and. e_ndx > 0

      if (.not. has_ions) return

!-----------------------------------------------------------------------
!	... initialize module variables
!-----------------------------------------------------------------------
      twopi  = 2._r8*pi
      dtr    = pi/180._r8

!-----------------------------------------------------------------------
!	... set molecular weights
!-----------------------------------------------------------------------
      call cnst_get_ind( 'O2', m )
      rmass_o2    = cnst_mw(m)
      rmassinv_o2 = 1._r8/rmass_o2
      call cnst_get_ind( 'O', m )
      rmass_o1    = cnst_mw(m)
      rmassinv_o1 = 1._r8/rmass_o1
      call cnst_get_ind( 'N', m )
      rmass_n2    = 2._r8*cnst_mw(m)
      rmassinv_n2 = 1._r8/rmass_n2

      offa(isouth)   = 4.3_r8*dtr
      offa(inorth)   = 3.7_r8*dtr
      phid(isouth)   = 0._r8
      phid(inorth)   = 0._r8
      alfa_1         = 2._r8
      alfa_2         = 2._r8

!-----------------------------------------------------------------------
! Values from 10/05/94 HPI estimates (50% or more higher than old estimates):
!     alfa_1 = amin1(1.5,1.25+0.05*plevel)
!     alfa_2 = 1.2 + 0.095*plevel
!-----------------------------------------------------------------------
      alfa0  = 0.5_r8*(alfa_1 + alfa_2)
      ralfa  = (alfa_2 - alfa_1) / (alfa_1 + alfa_2 + e)
      alfa21 = 0.075_r8
      alfa22 = 0.075_r8
      alfa20 = 0.5_r8 * (alfa21 + alfa22)
      ralfa2 = (alfa22 - alfa21) / (alfa21 + alfa22 + e)
      e21    = 1.e-80_r8
      e22    = 1.e-80_r8
      e20    = 0.5_r8 * (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:
!-----------------------------------------------------------------------
      flx_sp = e_sp/1.602e-6_r8

!-----------------------------------------------------------------------
! 	... set auroral lower bndy index
!-----------------------------------------------------------------------
      plb = 5.e-4_r8*exp( 7._r8 ) * .1_r8             ! Pa
      do k = 1,pver
	 if( pref_mid(k) >= plb ) then
	    lev1 = k-1
	    exit
	 end if
      end do

      if (masterproc) write(iulog,*) ' '
      if (masterproc) write(iulog,*) 'aurora_inti: aurora will go down to lev,p = ',lev1,pref_mid(lev1)
      if (masterproc) write(iulog,*) ' '

!-----------------------------------------------------------------------
! Report to stdout:
!-----------------------------------------------------------------------
#ifdef AURORA_DIAGS
        write(iulog,"(/,'aurora_cons:')")
        write(iulog,"('  cusp:    alfac=',f8.3,' ec=',f8.3,' fc=',e10.4)") &
          alfac,ec,fc
        write(iulog,"('  drizzle: alfad=',f8.3,' ed=',f8.3,' fd=',e10.4)") &
          alfad,ed,fd
        write(iulog,"('  half-widths = h1,h2=',2f10.3)") h1,h2
        write(iulog,"('  energy flux = e1,e2=',2f10.3)") e1,e2
        write(iulog,"('  add_sproton = ',l1)") add_sproton
        if( add_sproton ) then
           write(iulog,"('  solar protons: alfa_sp=',f8.3, &
          ' e_sp = ',e10.4,' flx_sp = ',e10.4)") alfa_sp,e_sp,flx_sp
        end if
        if( add_helectron ) then
          write(iulog,"('  high-energy electrons: alfa30 = ',f8.3, &
          ' e30 = ',f8.3)") alfa30,e30
        end if
        write(iulog,"(' ')")
#endif

        call addfld( 'QSUM', (/ 'lev' /), 'I', '/s', 'total ion production' )

      end subroutine aurora_inti

      subroutine aurora_timestep_init
!-----------------------------------------------------------------------
! 	... per timestep initialization
!-----------------------------------------------------------------------

      use mag_parms,  only : get_mag_parms
      use spmd_utils, only : masterproc

!-----------------------------------------------------------------------
!	... local variables
!-----------------------------------------------------------------------
      real(r8) :: power, plevel
      real(r8) :: roth, rote, rcp, rhp
      real(r8) :: arad

      if (.not. has_ions) return

!-----------------------------------------------------------------------
!	... get hemispheric power
!-----------------------------------------------------------------------
      call get_mag_parms( by = byimf, hpower = power, ctpoten = ctpoten )
#ifdef AURORA_DIAGS
      if( masterproc ) then
         write(iulog,*) '----------------------------------------'
         write(iulog,*) 'aurora_timestep_init: by,power,ctpoten = ',byimf,power,ctpoten
         write(iulog,*) '----------------------------------------'
      end if
#endif

      if( power >= .01_r8 ) then
         plevel = 2.09_r8*log( power )
      else
         plevel = 0._r8
      end if

! fvitt -- moved the calc of h1, h2, rh, and h0 from aurora_inti
!          This was done to for bit-for-bit restarts.
!          These aurora oval dimension quantities power dependent and should be updated.
!-----------------------------------------------------------------------
! h1 = Gaussian half-width of the noon auroral oval in degrees
! h2 = Gaussian half-width of the midnight auroral oval in degrees
!-----------------------------------------------------------------------
!      h1 = 3._r8
!      h2 = 10._r8
! 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_r8, 0.83_r8 + 0.33_r8*plevel)
      h2 = 2.5_r8+0.025_r8*max(power,55._r8)+0.01_r8*min(0._r8,power-55._r8)
!-----------------------------------------------------------------------
! Values from corrections to Emery et al Parameterization report:
!     h1 = amin1(2.35, 0.83 + 0.33*plevel)
!     h2 = 2.87 + 0.15*plevel
!-----------------------------------------------------------------------

      rh = (h2 - h1) / (h1 + h2)
      h0     = 0.5_r8 * (h1 + h2) * dtr

      theta0(isouth) = (-3.80_r8 + 8.48_r8*(ctpoten**.1875_r8))*dtr
      theta0(inorth) = theta0(isouth)
      dskofa(isouth) = (-1.26_r8 + 0.15_r8 * byimf)*dtr
      dskofa(inorth) = dskofa(isouth)
      roth           = (12.18_r8 - 0.89_r8 * plevel)
      rote           = ( 2.62_r8 - 0.55_r8 * plevel)
      rroth          = roth * dtr
      rrote          = rote * dtr

!-----------------------------------------------------------------------
! 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 = (0.5_r8 + 0.15_r8 * power)
!      e2 = (1.5_r8 + 0.25_r8 * power)
! modified by LQIAN, 2008
! 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_r8, -2.15_r8 + 0.62_r8*plevel)
      e2=1._r8+0.11_r8*power
!-----------------------------------------------------------------------
! Values from corrections to Emery et al Parameterization report:
!-----------------------------------------------------------------------
      e0  = 0.5_r8 * (e1 + e2)
      ree = (e2 - e1) / (e1 + e2)

      rhp          = 14.20_r8 + 0.96_r8*plevel
      rcp          = -0.43_r8 + 9.69_r8 * (ctpoten**.1875_r8)
      arad         = max( rcp,rhp )
      rrad(isouth) = arad*dtr 
      rrad(inorth) = arad*dtr

      end subroutine aurora_timestep_init

      subroutine aurora_prod( tn, o2, o1, mbar, rlats, &
                              qo2p, qop, qn2p, qnp, pmid, &
                              lchnk, calday,  ncol, rlons, pbuf )
!-----------------------------------------------------------------------
! 	... auroral parameterization driver
!-----------------------------------------------------------------------

      use mo_apex,     only : alatm, alonm                      ! magnetic latitude,longitude grid (radians)
      use ppgrid,      only : pcols, pver
      use cam_history, only : outfld
      use physics_buffer,only: physics_buffer_desc

      implicit none

!-----------------------------------------------------------------------
! 	... dummy arguments
!-----------------------------------------------------------------------
      integer, intent(in) ::  &
       ncol, &                           ! column count
       lchnk                             ! chunk index
      real(r8), intent(in) :: &
        calday                           ! calendar day of year
      real(r8), intent(in) :: &
        tn(pcols,pver), &                ! neutral gas temperature (K)
        o2(ncol,pver), &                 ! O2 concentration (kg/kg)
        o1(ncol,pver), &                 ! O concentration (kg/kg)
        mbar(ncol,pver)                  ! mean molecular weight (g/mole)
      real(r8), intent(in) :: &
        pmid(pcols,pver)                 ! midpoint pressure (Pa)
      real(r8), intent(in) :: &
        rlats(ncol), &                   ! column latitudes (radians)
        rlons(ncol)
      real(r8), intent(out) :: &
        qo2p(ncol,pver), &               ! o2+ production
        qop(ncol,pver), &                ! o+ production
        qn2p(ncol,pver), &               ! n2+ production
        qnp(ncol,pver)                   ! n+ production

      type(physics_buffer_desc),pointer :: pbuf(:)

!-----------------------------------------------------------------------
! 	... local variables
!-----------------------------------------------------------------------
      integer  :: i, k
      integer  :: hemis(ncol)
      real(r8) :: r2d
      real(r8) :: ofda, cosofa, sinofa, aslona
      real(r8) :: sunlons(ncol)                     ! sun's mag longitudes
      real(r8) :: dlat_aur(ncol)
      real(r8) :: dlon_aur(ncol)
      real(r8) :: colat(ncol)
      real(r8) :: sinlat(ncol)
      real(r8) :: coslat(ncol)
      real(r8) :: coslon(ncol)
      real(r8) :: sinlon(ncol)
      real(r8) :: alon(ncol)
      real(r8) :: cusp(ncol)
      real(r8) :: alfa(ncol)
      real(r8) :: alfa2(ncol)
      real(r8) :: alfa3(ncol)
      real(r8) :: flux(ncol)
      real(r8) :: flux2(ncol)
      real(r8) :: flux3(ncol)
      real(r8) :: drizl(ncol)
      real(r8) :: qteaur(ncol)                         ! for electron temperature
      logical  :: do_aurora(ncol)

      if (.not. has_ions) return

!-----------------------------------------------------------------------
! 	... initialize ion production
!-----------------------------------------------------------------------
      do k = 1,pver
        qo2p(:,k) = 0._r8
        qop(:,k)  = 0._r8
        qn2p(:,k) = 0._r8
        qnp(:,k)  = 0._r8
      end do

      r2d = 180._r8/pi

!-----------------------------------------------------------------------
! 	... output mag lons, lats
!-----------------------------------------------------------------------
      call outfld( 'ALONM', r2d*alonm(:ncol,lchnk), pcols, lchnk )
      call outfld( 'ALATM', r2d*alatm(:ncol,lchnk), pcols, lchnk )

!-----------------------------------------------------------------------
! 	... check latitudes, and return if all below 32.5 deg
!-----------------------------------------------------------------------
      do_aurora(:) = abs( rlats(:) ) > pi/6._r8
      if( all( .not. do_aurora(:) ) ) then
         return
      end if

!-----------------------------------------------------------------------
! 	... set sun location
!-----------------------------------------------------------------------
      call sunloc( calday, sunlons, lchnk, ncol )

      do i = 1,ncol
        if( do_aurora(i) ) then
          dlat_aur(i) = alatm(i,lchnk)
          dlon_aur(i) = alonm(i,lchnk) - sunlons(i)
          if( dlon_aur(i) > pi ) then
             dlon_aur(i) = dlon_aur(i) - twopi
          else if( dlon_aur(i) < -pi ) then
             dlon_aur(i) = dlon_aur(i) + twopi
          end if
          if( dlat_aur(i) > 0._r8 ) then
            hemis(i) = 2
          else
            hemis(i) = 1
          end if
!-----------------------------------------------------------------------
! 	... find auroral circle coordinates
!-----------------------------------------------------------------------
            ofda      = sqrt( offa(hemis(i))**2 + dskofa(hemis(i))**2)
            cosofa    = cos( ofda )
            sinofa    = sin( ofda )
            aslona    = asin( dskofa(hemis(i))/ofda )
            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)   = mod( atan2( sinlon(i)*coslat(i),sinlat(i)*sinofa &
                                    + cosofa*coslat(i)*coslon(i) ) - aslona + 3._r8*pi,twopi) - pi
        end if
      end do
#ifdef AURORA_DIAGS
      write(iulog,*) '-----------------------------------------------------'
      write(iulog,*) 'aurora: diagnostics for lchnk = ',lchnk
      write(iulog,*) '        geo lats'
      write(iulog,'(1p,5g15.7)') r2d*rlats(:ncol)
      write(iulog,*) '        geo lons'
      write(iulog,'(1p,5g15.7)') r2d*rlons(:ncol)
      write(iulog,*) '        mag lats'
      write(iulog,'(1p,5g15.7)') r2d*dlat_aur(:ncol)
      write(iulog,*) '        mag lons'
      write(iulog,'(1p,5g15.7)') r2d*alonm(:ncol,lchnk)
      write(iulog,*) '        mag table lons'
      write(iulog,'(1p,5g15.7)') r2d*dlon_aur(:ncol)
      write(iulog,*) '        sunlons'
      write(iulog,'(1p,5g15.7)') r2d*sunlons(:ncol)
      write(iulog,*) '     min,max mag lons = ',r2d*minval(alonm(:ncol,lchnk)),r2d*maxval(alonm(:ncol,lchnk))
      write(iulog,*) '-----------------------------------------------------'
#endif

!-----------------------------------------------------------------------
! 	... make cusp
!-----------------------------------------------------------------------
      call aurora_cusp( cusp, do_aurora, hemis, colat, alon, ncol )

!-----------------------------------------------------------------------
! 	... make alfa, flux, and drizzle
!-----------------------------------------------------------------------
      call aurora_heat( flux, flux2, flux3, alfa, alfa2, &
                        alfa3, qteaur, drizl, do_aurora, hemis, &
                        alon, colat, ncol )

!-----------------------------------------------------------------------
! 	... auroral additions to ionization rates
!-----------------------------------------------------------------------
      call aurora_ions( drizl, cusp, alfa, alfa2, alfa3, &
                        flux, flux2, flux3, tn, o2, &
                        o1, mbar, qo2p, qop, qn2p, &
                        qnp, pmid, do_aurora, ncol, lchnk, pbuf )

      end subroutine aurora_prod

      subroutine aurora_hrate( tn, o2, o1, mbar, rlats, &
                               aur_hrate, cpair, pmid, lchnk, calday, &
                               ncol, rlons )
!-----------------------------------------------------------------------
! 	... auroral parameterization driver
!-----------------------------------------------------------------------

      use mo_apex, only : alatm, alonm                      ! magnetic latitude,longitude grid (radians)
      use ppgrid,  only : pcols, pver

      implicit none

!-----------------------------------------------------------------------
! 	... dummy arguments
!-----------------------------------------------------------------------
      integer, intent(in) ::  &
       ncol, &                           ! column count
       lchnk                             ! chunk index
      real(r8), intent(in) :: &
        calday                           ! calendar day of year
      real(r8), intent(in) :: &
        tn(pcols,pver), &                ! neutral gas temperature (K)
        o2(ncol,pver), &                 ! O2 concentration (kg/kg)
        o1(ncol,pver), &                 ! O concentration (kg/kg)
        mbar(ncol,pver)                  ! mean molecular weight (g/mole)
      real(r8), intent(in) :: &
        cpair(ncol,pver)                 ! specific heat capacity (J/K/kg)
      real(r8), intent(in) :: &
        pmid(pcols,pver)                 ! midpoint pressure (Pa)
      real(r8), intent(in) :: &
        rlats(ncol), &                   ! column latitudes (radians)
        rlons(ncol)
      real(r8), intent(out) :: &
        aur_hrate(ncol,pver)             ! auroral heating rate

!-----------------------------------------------------------------------
! 	... local variables
!-----------------------------------------------------------------------
      real(r8), parameter :: aur_therm     = 807._r8
      real(r8), parameter :: jkcal         = 4184._r8
      real(r8), parameter :: aur_heat_eff  = .05_r8
      real(r8), parameter :: aur_hconst    = 1.e3_r8*jkcal*aur_therm*aur_heat_eff

      integer  :: i, k
      integer  :: hemis(ncol)
      real(r8) :: r2d
      real(r8) :: ofda, cosofa, sinofa, aslona
      real(r8) :: sunlons(ncol)                     ! sun's mag longitudes
      real(r8) :: dlat_aur(ncol)
      real(r8) :: dlon_aur(ncol)
      real(r8) :: colat(ncol)
      real(r8) :: sinlat(ncol)
      real(r8) :: coslat(ncol)
      real(r8) :: coslon(ncol)
      real(r8) :: sinlon(ncol)
      real(r8) :: alon(ncol)
      real(r8) :: cusp(ncol)
      real(r8) :: alfa(ncol)
      real(r8) :: alfa2(ncol)
      real(r8) :: alfa3(ncol)
      real(r8) :: flux(ncol)
      real(r8) :: flux2(ncol)
      real(r8) :: flux3(ncol)
      real(r8) :: drizl(ncol)
      real(r8) :: qteaur(ncol)                         ! for electron temperature
      real(r8) :: qsum(ncol,pver)                      ! total ion production (1/s)
      logical  :: do_aurora(ncol)

!-----------------------------------------------------------------------
! 	... initialize ion production
!-----------------------------------------------------------------------
      do k = 1,pver
        aur_hrate(:,k) = 0._r8
      end do

      if (.not. has_ions) return

      r2d = 180._r8/pi

!-----------------------------------------------------------------------
! 	... check latitudes, and return if all below 32.5 deg
!-----------------------------------------------------------------------
      do_aurora(:) = abs( rlats(:) ) > pi/6._r8
      if( all( .not. do_aurora(:) ) ) then
         return
      end if

!-----------------------------------------------------------------------
! 	... set sun location
!-----------------------------------------------------------------------
      call sunloc( calday, sunlons, lchnk, ncol )

      do i = 1,ncol
        if( do_aurora(i) ) then
          dlat_aur(i) = alatm(i,lchnk)
          dlon_aur(i) = alonm(i,lchnk) - sunlons(i)
          if( dlon_aur(i) > pi ) then
             dlon_aur(i) = dlon_aur(i) - twopi
          else if( dlon_aur(i) < -pi ) then
             dlon_aur(i) = dlon_aur(i) + twopi
          end if
          if( dlat_aur(i) > 0._r8 ) then
            hemis(i) = 2
          else
            hemis(i) = 1
          end if
!-----------------------------------------------------------------------
! 	... find auroral circle coordinates
!-----------------------------------------------------------------------
            ofda      = sqrt( offa(hemis(i))**2 + dskofa(hemis(i))**2)
            cosofa    = cos( ofda )
            sinofa    = sin( ofda )
            aslona    = asin( dskofa(hemis(i))/ofda )
            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)   = mod( atan2( sinlon(i)*coslat(i),sinlat(i)*sinofa &
                                    + cosofa*coslat(i)*coslon(i) ) - aslona + 3._r8*pi,twopi) - pi
        end if
      end do
#ifdef AURORA_DIAGS
      write(iulog,*) '-----------------------------------------------------'
      write(iulog,*) 'aurora: diagnostics for lchnk = ',lchnk
      write(iulog,*) '        geo lats'
      write(iulog,'(1p,5g15.7)') r2d*rlats(:ncol)
      write(iulog,*) '        geo lons'
      write(iulog,'(1p,5g15.7)') r2d*rlons(:ncol)
      write(iulog,*) '        mag lats'
      write(iulog,'(1p,5g15.7)') r2d*dlat_aur(:ncol)
      write(iulog,*) '        mag lons'
      write(iulog,'(1p,5g15.7)') r2d*alonm(:ncol,lchnk)
      write(iulog,*) '        mag table lons'
      write(iulog,'(1p,5g15.7)') r2d*dlon_aur(:ncol)
      write(iulog,*) '        sunlons'
      write(iulog,'(1p,5g15.7)') r2d*sunlons(:ncol)
      write(iulog,*) '     min,max mag lons = ',r2d*minval(alonm(:ncol,lchnk)),r2d*maxval(alonm(:ncol,lchnk))
      write(iulog,*) '-----------------------------------------------------'
#endif

!-----------------------------------------------------------------------
! 	... make cusp
!-----------------------------------------------------------------------
      call aurora_cusp( cusp, do_aurora, hemis, colat, alon, ncol )

!-----------------------------------------------------------------------
! 	... make alfa, flux, and drizzle
!-----------------------------------------------------------------------
      call aurora_heat( flux, flux2, flux3, alfa, alfa2, &
                        alfa3, qteaur, drizl, do_aurora, hemis, &
                        alon, colat, ncol )

!-----------------------------------------------------------------------
! 	... auroral additions to ionization rates
!-----------------------------------------------------------------------
      call total_ion_prod( drizl, cusp, alfa, alfa2, alfa3, &
                           flux, flux2, flux3, tn, o2, &
                           o1, mbar, qsum, pmid, do_aurora, &
                           ncol, lchnk )
!-----------------------------------------------------------------------
! 	... form auroral heating rate
!-----------------------------------------------------------------------
      do k = 1,pver
         aur_hrate(:,k) = aur_hconst * qsum(:,k) / (cpair(:,k) * mbar(:,k))
      end do

      end subroutine aurora_hrate

      subroutine aurora_cusp( cusp, do_aurora, hemis, colat, alon, ncol )
!-----------------------------------------------------------------------
! 	... calculate horizontal variation of polar cusp heating
!-----------------------------------------------------------------------

      implicit none

!-----------------------------------------------------------------------
! 	... dummy arguments
!-----------------------------------------------------------------------
      integer, intent(in)     :: ncol
      integer, intent(in)     :: hemis(ncol)
      real(r8), intent(in)    :: colat(ncol)
      real(r8), intent(in)    :: alon(ncol)
      real(r8), intent(out)   :: cusp(ncol)
      logical, intent(in)     :: do_aurora(ncol)

!-----------------------------------------------------------------------
! 	... local variables
!-----------------------------------------------------------------------
      real(r8), parameter :: s5  =.08726646_r8, &
                             s20 =.34906585_r8

      where( do_aurora(:) )
         cusp(:) = (exp( -((theta0(hemis(:)) - colat(:))/s5)**2 ) &
                      + exp( -((pi - theta0(hemis(:)) - colat(:))/s5)**2) ) &
                        *exp( -(atan2( sin(alon(:) - phid(hemis(:))), cos(alon(:) - phid(hemis(:))) )/s20)**2 )
      elsewhere
         cusp(:) = 0._r8
      endwhere

      end subroutine aurora_cusp 

      subroutine aurora_heat( flux, flux2, flux3, alfa, alfa2, &
                              alfa3, qteaur, drizl, do_aurora, hemis, &
                              alon, colat, ncol )
!-----------------------------------------------------------------------
! 	... calculate alfa, flux, and drizzle
!-----------------------------------------------------------------------

      implicit none

!-----------------------------------------------------------------------
! 	... dummy arguments
!-----------------------------------------------------------------------
      integer, intent(in)     :: ncol
      integer, intent(in)     :: hemis(ncol)
      real(r8), intent(in)    :: colat(ncol)
      real(r8), intent(in)    :: alon(ncol)
      real(r8), intent(inout) :: flux(ncol)
      real(r8), intent(inout) :: flux2(ncol)
      real(r8), intent(inout) :: flux3(ncol)
      real(r8), intent(inout) :: drizl(ncol)
      real(r8), intent(inout) :: qteaur(ncol)
      real(r8), intent(inout) :: alfa(ncol)
      real(r8), intent(inout) :: alfa2(ncol)
      real(r8), intent(inout) :: alfa3(ncol)
      logical, intent(in)     :: do_aurora(ncol)

!-----------------------------------------------------------------------
! 	... local variables
!-----------------------------------------------------------------------
      real(r8), dimension(ncol) :: &
        coslamda, &                             ! cos(angle from throat)
        halfwidth, &                            ! oval half-width
        wrk, &                                  ! temp wrk array
        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
!-----------------------------------------------------------------------
      where( do_aurora(:) )
         coslamda(:) = cos( atan2( sin( alon(:) - rrote ),cos( alon(:) - rrote ) ) )
!-----------------------------------------------------------------------
! 	... auroral oval half-width (equation (1) in Roble,1987):
!-----------------------------------------------------------------------
         halfwidth(:) = h0*(1._r8 - rh*cos( atan2( sin(alon(:) - rroth),cos( alon(:) - rroth ) ) ) )
         dtheta(:)    = colat(:) - rrad(hemis(:))
      endwhere
!-----------------------------------------------------------------------
! 	... characteristic energy (equation (2) in Roble,1987):
!-----------------------------------------------------------------------
      if( alfa0 > .01_r8 ) then
         where( do_aurora(:) )
            alfa(:) =  alfa0*(1._r8 - ralfa*coslamda(:))
         endwhere
      else
         alfa(:) =  0._r8
      end if

      where( do_aurora(:) )
!-----------------------------------------------------------------------
! 	... flux, drizzle, alfa2, flux2
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
! 	... energy flux (equation (3) in Roble,1987):
!-----------------------------------------------------------------------
         wrk(:)   = exp( -(dtheta(:)/halfwidth(:))**2 )
         flux(:)  = e0*(1._r8 - ree*coslamda(:))*wrk(:) / (2._r8*alfa(:)*1.602e-9_r8)
         drizl(:) = exp( -((dtheta(:) + abs(dtheta(:)))/(2._r8*h0))**2 )
         alfa2(:) = alfa20*(1._r8 - ralfa2*coslamda(:))
         flux2(:) = e20*(1._r8 - re2*coslamda(:))*wrk(:) / (2._r8*alfa2(:)*1.602e-9_r8)
!-----------------------------------------------------------------------
! 	... alfa3, flux3 for high energy electrons:
!-----------------------------------------------------------------------
         alfa3(:) = alfa30
         flux3(:) = e30*wrk(:) / 1.602e-6_r8
!-----------------------------------------------------------------------
! 	... for electron temperature (used in settei):  
!-----------------------------------------------------------------------
         qteaur(:) = -7.e8_r8*wrk(:)
      endwhere

      end subroutine aurora_heat

      subroutine aurora_ions( drizl, cusp, alfa1, alfa2, alfa3, &
                              flux1, flux2, flux3, tn, o2, &
                              o1, mbar, qo2p, qop, qn2p, &
                              qnp, pmid, do_aurora, ncol, lchnk, pbuf )
!-----------------------------------------------------------------------
! 	... calculate auroral additions to ionization rates
!-----------------------------------------------------------------------

      use ppgrid,      only : pcols, pver
      use cam_history, only : outfld

      use physics_buffer,only: physics_buffer_desc, pbuf_set_field

      implicit none

!-----------------------------------------------------------------------
! 	... dummy arguments
!-----------------------------------------------------------------------
      integer, intent(in) :: ncol
      integer, intent(in) :: lchnk
      real(r8), intent(in), dimension(ncol) :: &
                             drizl, &
                             cusp, &
                             alfa1, &
                             alfa2, &
                             alfa3, &
                             flux1, &
                             flux2, &
                             flux3
      real(r8), dimension(pcols,pver), intent(in) :: &
                             tn, &                     ! midpoint neutral temperature (K)
                             pmid                      ! midpoint pressure (Pa)
      real(r8), dimension(ncol,pver), intent(in) :: &
                             o2, &                     ! midpoint o2 concentration (kg/kg)
                             o1, &                     ! midpoint o  concentration (kg/kg)
                             mbar                      ! mean molecular mass (g/mole)
      real(r8), dimension(ncol,pver), intent(inout) :: &
                             qo2p, &                   ! o2p prod from aurora (molecules/cm^3/s)
                             qop, &                    ! op prod from aurora (molecules/cm^3/s)
                             qn2p, &                   ! n2p prod from aurora (molecules/cm^3/s)
                             qnp                       ! np prod from aurora (molecules/cm^3/s)
      logical, intent(in) :: do_aurora(ncol)

      type(physics_buffer_desc),pointer :: pbuf(:)

!-----------------------------------------------------------------------
! 	... local variables
!-----------------------------------------------------------------------
      real(r8), parameter :: const0        = 1.e-20_r8

      integer  :: i, k
      real(r8), dimension(ncol) :: &
        p0ez, &
        press, &                                   ! pressure at interface levels (dyne/cm^2)
        tempi, &                                   ! temperature at interface levels (K)
        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, &
        barm, &
        falfa1, &
        falfa2, &
        fcusp, &
        fdrizl, &
        falfa_sp, &
        xn2, &
        falfa3
      real(r8), dimension(ncol) :: &
        qo2p_aur, &
        qop_aur, &
        qn2p_aur                                   ! auroral ionization for O2+, O+, N2+
      real(r8) :: qia(5)                           ! low energy proton source (not in use, 1/02)
      real(r8) :: wrk(ncol,pver)

      qia(:) = 0._r8
      wrk(:,:) = 0._r8

level_loop : &
      do k = 1,lev1
          where( do_aurora(:) )
             press(:ncol) = 10._r8*pmid(:ncol,k)              ! from Pa to dyne/cm^2
             tempi(:ncol) = tn(:ncol,k)
             barm(:)      = mbar(:,k)
             p0ez(:)      = (press(:)/(grav*4.e-6_r8))**.606_r8
             xalfa1(:)    = p0ez(:)/alfa1(:)
             xalfa2(:)    = p0ez(:)/alfa2(:)
             xcusp (:)    = p0ez(:)/alfac
             xdrizl(:)    = p0ez(:)/alfad

!-----------------------------------------------------------------------
! 	... initialize (whole array operations):
!-----------------------------------------------------------------------
             flux1_ion(:) = const0
             flux2_ion(:) = const0
             alfa1_ion(:) = const0
             alfa2_ion(:) = const0
             cusp_ion(:)  = const0
             drizl_ion(:) = const0
          endwhere
!-----------------------------------------------------------------------
! 	... auroral electrons
!-----------------------------------------------------------------------
          call aion( xalfa1, alfa1_ion, do_aurora, ncol )
          call aion( xalfa2, alfa2_ion, do_aurora, ncol )
          call aion( xcusp , cusp_ion, do_aurora, ncol  )
          call aion( xdrizl, drizl_ion, do_aurora, ncol )
          where( do_aurora(:) )
             falfa1(:) = alfa1(:)*flux1(:)  ! s7
             falfa2(:) = alfa2(:)*flux2(:)  ! s8
             fcusp (:) = cusp(:)*alfac*fc   ! s9
             fdrizl(:) = drizl(:)*alfad*fd  ! s10
             qsum(:)   = falfa1(:)*alfa1_ion(:) &    ! s7*s3
                       + falfa2(:)*alfa2_ion(:) &    ! s8*s4
                       + fcusp(:)*cusp_ion (:) &     ! s9*s5
                       + fdrizl(:)*drizl_ion(:)       ! s10*s6
          endwhere
!-----------------------------------------------------------------------
! 	... include solar protons if add_sproton is set, 
!           and high energy electrons if add_helectron is set
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
! 	... solar protons
!-----------------------------------------------------------------------
          if( add_sproton ) then
             if( flx_sp > 1.e-19_r8 ) then
                where( do_aurora(:) )
                   p0ez_mbar(:) = press(:)/(boltz_cgs*tempi(:)*barm(:))/avo
                   tk_mbar(:)   = gask*tempi(:)/(grav*barm(:))*p0ez_mbar(:)
                   xalfa_sp(:)  = ((tk_mbar(:)/.00271_r8)**.58140_r8)/alfa_sp
                endwhere
                call bion( xalfa_sp, alfasp_bion, do_aurora, ncol )
                where( do_aurora(:) )
                   falfa_sp(:) = drizl(:)*p0ez_mbar(:)*flx_sp*1.e6_r8/(tk_mbar(:)*35._r8)
                   qsum(:)     = qsum(:) + falfa_sp(:)*alfasp_bion(:)
                endwhere
             end if
          end if
!-----------------------------------------------------------------------
! 	... high energy electrons
!-----------------------------------------------------------------------
          if( add_helectron ) then
             if( e30 > 1.e-19_r8 ) then
                where( do_aurora(:) )
                   xalfa3(:) = p0ez(:)/alfa3(:)        ! alfa3(:)==alfa30
                endwhere
                call aion( xalfa3, alfa3_ion, do_aurora, ncol )
                where( do_aurora(:) )
                   falfa3(:) = alfa3(:)*flux3(:)  ! s13 (high energy electrons)
                   qsum(:)   = qsum(:) + falfa3(:)*alfa3_ion(:)
                endwhere
             end if
          end if

!-----------------------------------------------------------------------
! 	... form production
!-----------------------------------------------------------------------
          where( do_aurora(:) )
             barm_t(:) = grav*barm(:)/(35.e-3_r8*gask*tempi(:))
             qsum(:)   = qsum(:)*barm_t(:)               ! s1 = s1*s11
             wrk(:,k)  = qsum(:)
!-----------------------------------------------------------------------
! 	... denominator of equations (13-16) in Roble,1987.
!-----------------------------------------------------------------------
             xn2(:)   = max( (1._r8 - o2(:,k) - o1(:,k)),1.e-8_r8 )
             denom(:) = 0.92_r8*xn2(:)*rmassinv_n2 &
                      + 1.5_r8*o2(:,k) *rmassinv_o2 + 0.56_r8*o1(:,k) *rmassinv_o1
!-----------------------------------------------------------------------
! 	... production of O2+ (equation (15) in Roble,1987):
!-----------------------------------------------------------------------
             qo2p_aur(:) = qsum(:)*o2(:,k)/(rmass_o2*denom(:)) + qia(2)
!-----------------------------------------------------------------------
! 	... production of O+ (equation (16) in Roble,1987):
!-----------------------------------------------------------------------
             qop_aur(:) = qsum(:)*(.5_r8 *o2(:,k)*rmassinv_o2 &
                                   + .56_r8*o1(:,k)*rmassinv_o1)/denom(:) + qia(3)
!-----------------------------------------------------------------------
! 	... production of N2+ (equation (13) in Roble,1987)
!-----------------------------------------------------------------------
             qn2p_aur(:) = qsum(:)*.7_r8*xn2(:)/(rmass_n2*denom(:)) + qia(1)
             qo2p(:,k)   = qo2p(:,k) + qo2p_aur(:)
             qop(:,k)    = qop(:,k) + qop_aur(:)
             qn2p(:,k)   = qn2p(:,k) + qn2p_aur(:)
             qnp(:,k)    = qnp (:,k) + .22_r8/.7_r8 * qn2p_aur(:)
          endwhere
      end do level_loop

!--------------------------------------------------------------------------------------------
!  Save sum of auroral ion production rates to be accessed in ionosphere module in physics 
!  buffer using pointer for WACCM-X
!--------------------------------------------------------------------------------------------
      if (indxAIPRS>0) then
         call pbuf_set_field(pbuf, indxAIPRS,  wrk(1:ncol,1:pver))
      endif

      call outfld( 'QSUM', wrk, ncol, lchnk )

      end subroutine aurora_ions

      subroutine total_ion_prod( drizl, cusp, alfa1, alfa2, alfa3, &
                                 flux1, flux2, flux3, tn, o2, &
                                 o1, mbar, tpions, pmid, do_aurora, &
                                 ncol, lchnk )
!-----------------------------------------------------------------------
! 	... calculate auroral additions to ionization rates
!-----------------------------------------------------------------------

      use ppgrid,      only : pcols, pver
      use cam_history, only : outfld

      implicit none

!-----------------------------------------------------------------------
! 	... dummy arguments
!-----------------------------------------------------------------------
      integer, intent(in) :: ncol
      integer, intent(in) :: lchnk
      real(r8), intent(in), dimension(ncol) :: &
                             drizl, &
                             cusp, &
                             alfa1, &
                             alfa2, &
                             alfa3, &
                             flux1, &
                             flux2, &
                             flux3
      real(r8), dimension(pcols,pver), intent(in) :: &
                             tn, &                     ! midpoint neutral temperature (K)
                             pmid                      ! midpoint pressure (Pa)
      real(r8), dimension(ncol,pver), intent(in) :: &
                             o2, &                     ! midpoint o2 concentration (kg/kg)
                             o1, &                     ! midpoint o  concentration (kg/kg)
                             mbar                      ! mean molecular mass (g/mole)
      real(r8), dimension(ncol,pver), intent(inout) :: &
                             tpions                    ! total ion production (1/s)
      logical, intent(in) :: do_aurora(ncol)

!-----------------------------------------------------------------------
! 	... local variables
!-----------------------------------------------------------------------
      real(r8), parameter :: const0        = 1.e-20_r8

      integer  :: i, k
      real(r8), dimension(ncol) :: &
        p0ez, &
        press, &                                   ! pressure at interface levels (dyne/cm^2)
        tempi, &                                   ! temperature at interface levels (K)
        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, &
        barm, &
        falfa1, &
        falfa2, &
        fcusp, &
        fdrizl, &
        falfa_sp, &
        xn2, &
        falfa3
      real(r8), dimension(ncol) :: &
        qo2p_aur, &
        qop_aur, &
        qn2p_aur                                   ! auroral ionization for O2+, O+, N2+
      real(r8) :: qia(5)                           ! low energy proton source (not in use, 1/02)

      qia(:)      = 0._r8
      tpions(:,:) = 0._r8

level_loop : &
      do k = 1,lev1
          where( do_aurora(:) )
             press(:ncol) = 10._r8*pmid(:ncol,k)              ! from Pa to dyne/cm^2
             tempi(:ncol) = tn(:ncol,k)
             barm(:)      = mbar(:,k)
             p0ez(:)      = (press(:)/(grav*4.e-6_r8))**.606_r8
             xalfa1(:)    = p0ez(:)/alfa1(:)
             xalfa2(:)    = p0ez(:)/alfa2(:)
             xcusp (:)    = p0ez(:)/alfac
             xdrizl(:)    = p0ez(:)/alfad

!-----------------------------------------------------------------------
! 	... initiliaze (whole array operations):
!-----------------------------------------------------------------------
             flux1_ion(:) = const0
             flux2_ion(:) = const0
             alfa1_ion(:) = const0
             alfa2_ion(:) = const0
             cusp_ion(:)  = const0
             drizl_ion(:) = const0
          endwhere
!-----------------------------------------------------------------------
! 	... auroral electrons
!-----------------------------------------------------------------------
          call aion( xalfa1, alfa1_ion, do_aurora, ncol )
          call aion( xalfa2, alfa2_ion, do_aurora, ncol )
          call aion( xcusp , cusp_ion, do_aurora, ncol  )
          call aion( xdrizl, drizl_ion, do_aurora, ncol )
          where( do_aurora(:) )
             falfa1(:) = alfa1(:)*flux1(:)  ! s7
             falfa2(:) = alfa2(:)*flux2(:)  ! s8
             fcusp (:) = cusp(:)*alfac*fc   ! s9
             fdrizl(:) = drizl(:)*alfad*fd  ! s10
             qsum(:)   = falfa1(:)*alfa1_ion(:) &    ! s7*s3
                       + falfa2(:)*alfa2_ion(:) &    ! s8*s4
                       + fcusp(:)*cusp_ion (:) &     ! s9*s5
                       + drizl(:)*drizl_ion(:)       ! s10*s6
          endwhere
!-----------------------------------------------------------------------
! 	... include solar protons if add_sproton is set, 
!           and high energy electrons if add_helectron is set
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
! 	... solar protons
!-----------------------------------------------------------------------
          if( add_sproton ) then
             if( flx_sp > 1.e-19_r8 ) then
                where( do_aurora(:) )
                   p0ez_mbar(:) = press(:)/(boltz_cgs*tempi(:)*barm(:))/avo
                   tk_mbar(:)   = gask*tempi(:)/(grav*barm(:))*p0ez_mbar(:)
                   xalfa_sp(:)  = ((tk_mbar(:)/.00271_r8)**.58140_r8)/alfa_sp
                endwhere
                call bion( xalfa_sp, alfasp_bion, do_aurora, ncol )
                where( do_aurora(:) )
                   falfa_sp(:) = drizl(:)*p0ez_mbar(:)*flx_sp*1.e6_r8/(tk_mbar(:)*35._r8)
                   qsum(:)     = qsum(:) + falfa_sp(:)*alfasp_bion(:)
                endwhere
             end if
          end if
!-----------------------------------------------------------------------
! 	... high energy electrons
!-----------------------------------------------------------------------
          if( add_helectron ) then
             if( e30 > 1.e-19_r8 ) then
                where( do_aurora(:) )
                   xalfa3(:) = p0ez(:)/alfa3(:)        ! alfa3(:)==alfa30
                endwhere
                call aion( xalfa3, alfa3_ion, do_aurora, ncol )
                where( do_aurora(:) )
                   falfa3(:) = alfa3(:)*flux3(:)  ! s13 (high energy electrons)
                   qsum(:)   = qsum(:) + falfa3(:)*alfa3_ion(:)
                endwhere
             end if
          end if

!-----------------------------------------------------------------------
! 	... form production
!-----------------------------------------------------------------------
          where( do_aurora(:) )
             barm_t(:)   = grav*barm(:)/(35.e-3_r8*gask*tempi(:))
             tpions(:,k) = qsum(:)*barm_t(:)               ! s1 = s1*s11
          endwhere
      end do level_loop

      end subroutine total_ion_prod

      subroutine aion( si, so, do_aurora, ncol )
!-----------------------------------------------------------------------
! Calculates integrated f(x) needed for total auroral ionization.
! See equations (10-12) in Roble,1987.
! Coefficients for equation (12) of Roble,1987 are in variable cc 
! (revised since 1987):
! Uses the identity x**y = exp(y*ln(x)) for performance 
! (fewer (1/2) trancendental functions are required).
!------------------------------------------------------------------------

      implicit none

!------------------------------------------------------------------------
! 	... dummy arguments
!------------------------------------------------------------------------
      integer,  intent(in)  :: ncol
      real(r8), intent(in)  :: si(ncol)
      real(r8), intent(out) :: so(ncol)
      logical,  intent(in)  :: do_aurora(ncol)

!------------------------------------------------------------------------
! 	... local variables
!------------------------------------------------------------------------
      real(r8), parameter :: cc(8) = &
       (/ 3.2333134511131_r8 ,  2.5658873458085_r8 ,  2.2540957232641_r8 , &
          0.72971983372673_r8,  1.1069072431948_r8 ,  1.7134937681128_r8 , &
          1.8835442312993_r8 ,  0.86472135072090_r8 /)

      real(r8) :: xlog(ncol)

      where( do_aurora(:) )
         xlog(:) = log( si(:) )
         so(:)   = cc(1)*exp( cc(2)*xlog(:) - cc(3)*exp( cc(4)*xlog(:) ) ) &
                   + cc(5)*exp( cc(6)*xlog(:) - cc(7)*exp( cc(8)*xlog(:) ) )
      elsewhere
         so(:) = 0._r8
      endwhere

      end subroutine aion

      subroutine bion( si, so, do_aurora, ncol )
!-----------------------------------------------------------------------
! Calculates integrated f(x) needed for total auroral ionization.
! See equations (10-12) in Roble,1987.
! Use the identity x**y = exp(y*ln(x)) for performance 
! (fewer (1/2) trancendental functions are required).
!-----------------------------------------------------------------------

      implicit none

!------------------------------------------------------------------------
! 	... dummy arguments
!------------------------------------------------------------------------
      integer, intent(in)   :: ncol
      real(r8), intent(in)  :: si(ncol)
      real(r8), intent(out) :: so(ncol)
      logical,  intent(in)  :: do_aurora(ncol)

!------------------------------------------------------------------------
! 	... local variables
!------------------------------------------------------------------------
      real(r8), parameter :: cc(8) = &
        (/ 0.12718_r8, 4.9119_r8, 1.8429_r8, 0.99336_r8, 0.52472_r8, &
           1.5565_r8,  .85732_r8, 1.4116_r8 /)

      real(r8) :: xlog(ncol)

      where( do_aurora(:) )
         xlog(:) = log( si(:) )
         so(:)   = cc(1)*exp( cc(2)*xlog(:) - cc(3)*exp( cc(4)*xlog(:) ) ) &
                   + cc(5)*exp( cc(6)*xlog(:) - cc(7)*exp( cc(8)*xlog(:) ) )
      elsewhere
         so(:) = 0._r8
      endwhere

      end subroutine bion

      subroutine sunloc( calday, sunlons, lchnk, ncol )
!-----------------------------------------------------------------------
! 	... calculate sun's longitude in dipole coordinates, defining sunlon
!-----------------------------------------------------------------------
      use dyn_grid,   only : get_dyn_grid_parm, get_horiz_grid_d
      use spmd_utils, only : masterproc
      use  phys_grid, only : get_lat_all_p, get_rlat_all_p, get_rlon_all_p
      use  mo_apex,   only : glonm                                              ! magnetic longitude grid (radians)

      implicit none

!-----------------------------------------------------------------------
! 	... dummy arguments
!-----------------------------------------------------------------------
      integer,  intent(in)   :: ncol
      integer,  intent(in)   :: lchnk
      real(r8), intent(in)   :: calday  ! calendar day of year
      real(r8), intent(out)  :: sunlons(ncol)

!-----------------------------------------------------------------------
! 	... local variables
!-----------------------------------------------------------------------
      integer  :: col
      integer  :: latndx(ncol)
      integer  :: sunlon_ndx
      integer  :: sunlon_ndxp1
      integer :: plon, plat
      real(r8), allocatable :: clon(:)
      real(r8) :: wght1, wght2, dellon, r2d
      real(r8) :: rlats(ncol), rlons(ncol)

      r2d = 180._r8/pi
!-----------------------------------------------------------------------
! 	... sun's geographic coordinates
!-----------------------------------------------------------------------
      plon = get_dyn_grid_parm('plon')
      plat = get_dyn_grid_parm('plat')
      allocate(clon(plon))
      call get_horiz_grid_d(plon,clon_d_out=clon)


      dellon       = .5_r8 - (calday - int(calday))
      sunlon_ndx   = mod( nint( dellon*plon ) - 1,plon ) + 1
      if( sunlon_ndx < 1 ) then
         sunlon_ndx = plon + sunlon_ndx
      end if
      sunlon_ndxp1 = mod( sunlon_ndx,plon ) + 1
      wght2        = min( 1._r8, max( (dellon*twopi - clon(sunlon_ndx))*plon/twopi,0._r8 ) )
      wght1        = 1._r8 - wght2
      deallocate(clon)

!-----------------------------------------------------------------------      
!        ... get chunck latitudes
!-----------------------------------------------------------------------      
      call get_lat_all_p( lchnk, ncol, latndx )
      call get_rlat_all_p( lchnk, ncol, rlats )
      call get_rlon_all_p( lchnk, ncol, rlons )

      do col = 1,ncol
!        sunlons(col) = wght1*glonm(sunlon_ndx,latndx(col)) + wght2*glonm(sunlon_ndxp1,latndx(col))
	 sunlons(col) = wght1*glonm(sunlon_ndx,plat/2) + wght2*glonm(sunlon_ndxp1,plat/2)
      end do

#ifdef AURORA_DIAGS
      write(iulog,*) '-----------------------------------------------------'
      write(iulog,*) 'sunloc: diagnostics for lchnk = ',lchnk
      write(iulog,*) '        dellon,sunlon_ndx,sunlon_ndxp1,wght1,wght2'
      write(iulog,'(1p,g15.7,2i6,2g15.7)') dellon,sunlon_ndx,sunlon_ndxp1,wght1,wght2
      write(iulog,*) '   phys lons'
      write(iulog,'(1p,5g15.7)') r2d*rlons(:ncol)
      write(iulog,*) '  phys lats'
!     write(iulog,'(10i5)') latndx(:ncol)
      write(iulog,'(1p,5g15.7)') r2d*rlats(:ncol)
      write(iulog,*) '    mag lons'
      write(iulog,'(1p,5g15.7)') r2d*glonm(sunlon_ndx,latndx(:ncol))
      write(iulog,*) '    mag lons'
      write(iulog,'(1p,5g15.7)') r2d*glonm(sunlon_ndxp1,latndx(:ncol))
      write(iulog,*) '-----------------------------------------------------'
      write(iulog,"('sunloc: sunlons=',/,(e12.4))") sunlons
#endif

      end subroutine sunloc

      end module mo_aurora
