!
! See D.J.Strickland, et.al., 
!  "Analytical representation of g-factors for rapid, accurate 
!   calculation of excitation rates in the dayside thermosphere", 
!   JGR Vol.102 No.A7, pp.14485-14498
! See http://www.cpi.com (products link) for access to coefficients
!  and idl routines for reading coeffs and calculating g-factors.
!
      module gfactors
      implicit none
      integer,parameter :: 
     +  nstate=21,	! # excited states
     +  ncoeff=10,	! # coefficients
     +  nsza=8,		! # solar zen angles
     +  nf107=3		! # f10.7 flux
!
! Names of excited states:
! (note 1085SOL was changed to 1085S)
      character(len=8) :: enames(nstate)=
     + (/'1085E   ','1085S   ','1134N   ','1134N2  ','1200N   ',
     +   '1200N2  ','1304    ','1356    ','1493N   ','1493N2  ',
     +   '5577    ','6300    ','834E    ','834SOL  ','989     ',
     +   'N22PG   ','N2C4    ','N2LBH   ','N2VK    ','OPLS2D  ',
     +   'OPLS2P  '/)
      character(len=8) :: exsp(nstate) = 
     + (/'N2      ','N2      ','N2      ','N2      ','N2      ',
     +   'N2      ','O       ','O       ','N2      ','N2      ',
     +   'O       ','O       ','O       ','O       ','O       ',
     +   'N2      ','N2      ','N2      ','N2      ','O       ',
     +   'O       '/)
!
! Solar zenith angles and f10.7 fluxes for which coefficients
! are provided:
      real ::
     +  sza(nsza)  = (/0.,45.,60.,70.,80.,85.,87.,90./),	
     +  f107(nf107)= (/75.,150.,250./)
!
! Strickland provides ncoeff (10) coefficients for each excited
!   state (enames(nstate)) at each  zenith angle (sza(nsza)), and
!   each f10.7 (f107(nf107)).
! These coeffs were taken from data files obtained from the 
!   above URL for Computational Physics, Inc. 
! (see ~foster/gfac for reformatting procedure to write the file
!  gfac_coeffs.dat):
!
      real :: coeffs(ncoeff,nsza,nf107,nstate)
      include "gfac_coeffs.dat"
      contains
!-------------------------------------------------------------------
      subroutine exrate(fname,fout,rho,xn2,xo,z,imx,kmx,glat,h,hdr)
      use hist,only: history,tgcmheader
      use proc,only: gcmlon,spval
!
! Calculate excitation rates for field fname at current latitude.
! (called by mkderived)
!
! Args:
      integer,intent(in) :: imx,kmx
      character(len=*),intent(in) :: fname
      real,intent(out) :: fout(imx,kmx)	! vol excitation rate output
      real,intent(in) :: 
     +  rho(imx,kmx),	! total density (o2+o+n2 (cm3))
     +  xn2(imx,kmx),	! n2 (cm3)
     +  xo(imx,kmx),	! o  (cm3)
     +  z(imx,kmx),	! z  (km)
     +  glat
      type(history),intent(in) :: h
      type(tgcmheader),intent(in) :: hdr
!
! Locals:
      character(len=8) :: ename		! state name
      integer :: i,ii,k,iday,ixstate
      real :: fmin,fmax,za,slt,gcoeffs(ncoeff),banks,gaus1,gaus2,
     +  tden,tdenlog,gfac
      real :: a,r,c,d,e,f,g,t,u,v	! coefficient names
!
! Externals:
      real,external :: getsza,fslt,bilin,quadrat
      integer,external :: ixfindc
!
! Processor field names for excitation states are: "E-" followed by 
! the Strickland state name:
!
      if (fname(1:2)/='E-') then
        write(6,"('>>> exrate: field ',a,' not an excitation state',
     +    ' (should begin with ''E-'')')") fname
        return
      endif
      write(ename,"(a)") fname(3:len_trim(fname))
      ixstate = ixfindc(enames,nstate,ename)
      if (ixstate==0) then
        write(6,"('>>> exrate: cannot find field ',a,' in enames=',
     +    /(5(a,' ')))") ename,enames
        return
      endif
      iday = h%iyd-h%iyd/1000*1000
!
! Longitude loop (sza will vary):
      lonloop: do i=1,imx
        slt = fslt(slt,h%ut,gcmlon(i),1)	! current local time
        za = getsza(iday,slt,glat,gcmlon(i))	! current solar zenith angle
        if (za >= 90.) then	! night
          fout(i,:) = 1.e-20
!         fout(i,:) = spval
          cycle lonloop
        endif
!
! Get coefficients for excited state at f107d, solar zenith angle za,
!   and current local time (bilinearly interpolate coeffs).
! real function bilin(f,xa,ya,nx,ny,x,y,spv,iprnt)
!
        do ii=1,ncoeff
          gcoeffs(ii) = bilin(coeffs(ii,:,:,ixstate),sza,f107,nsza,
     +      nf107,za,hdr%f107d,spval,1)
        enddo 
        a = gcoeffs(1) 
        r = gcoeffs(2) 
        c = gcoeffs(3) 
        d = gcoeffs(4) 
        e = gcoeffs(5) 
        f = gcoeffs(6) 
        g = gcoeffs(7) 
        t = gcoeffs(8) 
        u = gcoeffs(9) 
        v = gcoeffs(10) 
!
! Calculate g-factors (see Strickland, et.al., equation (9), p.14493)
!
        do k=1,kmx
          if (k < kmx) then
            tden = quadrat(rho(i,k:kmx),z(i,k:kmx),kmx-k+1)
          else
            tden = rho(i,k)
          endif
          tdenlog = log(tden)
          banks = a*(2.-(1.+r)*exp(-(c*tden)))*exp(-d*tden)
          gaus1 = e*exp(-((tdenlog-log(f))/log(g))**2)
          gaus2 = t*exp(-((tdenlog-log(u))/log(v))**2)
          gfac  = banks+gaus1+gaus2
!
! Volume excitation rate is g-factor times species density:
! (see Strickland, et.al., equation (5), p.14489):
!
          if (trim(exsp(ixstate))=='N2') then
            fout(i,k) = gfac*xn2(i,k)
          else
            fout(i,k) = gfac*xo(i,k)
          endif
        enddo
      enddo lonloop
      end subroutine exrate
      end module gfactors