c
c------------------------------------------------------------------
c
      subroutine mke6300(tn,xo2,xo,xn2,ht,te,xne,xo2p,id,fout,ie6300,
     +  iyd,ut,glat,glon,f107d,f107a)
c
      real tn(id),xo2(id),xo(id),xn2(id),ht(id),te(id),xne(id),
     +  xo2p(id),fout(id),sr63(id),qn2p(id),qop(id),glyb(id)
      parameter (b1d=1.2, rk3=8.e-12, a1d=0.0093, a6300=0.0071)
      data ncalls/0/
c
c 3/10/94:
c Given tn,te,o2,o,n2,o2+ (all dimensioned id), return volume
c   emission rate at 6300A in fout. (from doptuv.f in this lib)
c
      ncalls = ncalls+1
c
c 6/95: call solred and add SR63 if ie6300 > 0:
c     subroutine solred(ifrst,stl,zpht,xno,xno2,xnn2,tn,f107,f107a,
c    +  day,glat,glong,kmx,sr63,qn2p,qop,glyb)
c
      if (ie6300.gt.0) then
        slt = fslt(dum,ut,glon,1)
        day = float(iyd - iyd/1000*1000)
        call solred(ncalls,slt,ht,xo,xo2,xn2,tn,f107d,f107a,
     +    day,glat,glon,id,sr63,qn2p,qop,glyb)
      endif
      do k=1,id                                                   
        rk1 = 2.e-11 * exp(107.8/tn(k))                                 
        rk2 = 2.9e-11 * exp(67.5/tn(k))                                 
        if (te(k).ge.1200.) rkdr = 1.6e-7 * (te(k)/300.)**(-0.55)
        if (te(k).lt.1200.) rkdr = 1.95e-7 * (te(k)/300.)**(-0.7)
        fout(k) = (b1d*rkdr*xo2p(k)*xne(k)) /               
     +            (a1d+rk1*xn2(k)+rk2*xo2(k)+rk3*xo(k))*a6300  
        if (ie6300.gt.0) fout(k) = fout(k) + sr63(k)
      enddo
      return
      end
c
c------------------------------------------------------------------
c
      subroutine mke5577(tn,xo2,xo,xn2,ht,te,xne,xo2p,xo21d,
     +  id,fout,ie5577,iyd,ut,glat,glon,f107d,f107a)
c
c Calculate greenline emission:
c
c 6/95: Up to five components may be requested:
c ie5577(1) > 0 -> o1 recombination (original)
c   (need tn,xo2,xo,xn2)
c ie5577(2) > 0 -> dissociative recombination of o2+
c   (need te,xo2p,xne,xo21d)
c
! Comments for ie5577(3) and ie5577(4) changed 4/12/12 as per
! Stan Solomon. (They were formerly "photoelectron impact" and
! "airglow" respectively). See also comments below.
!
c ie5577(3) > 0 -> N2(A)+O energy transfer
c   (need tn,ht,xo2,xo,xn2,xo21d and solred)
c ie5577(4) > 0 -> photoelectron impact
c   (need tn,ht,xo2,xo,xn2,xo21d and solred)
c
c ie5577(5) > 0 -> photo dissoc of o2 by solar lyman-beta
c   (glyb returned by solred)
c If solred is called, also need iyd, ut, glat, glon, f107d, f107a
c Solred returns sr63 (see e6300), qn2p, and qop
c
      real tn(id),xo2(id),xo(id),xn2(id),fout(id),
     +     ht(id),te(id),xne(id),xo2p(id),xo21d(id),
     +     sr63(id),qn2p(id),qop(id),glyb(id)
      integer ie5577(5)
      parameter (a5=1.18, a6=1.35, cp=15., cpp=211.)
c
c O(1s)/O+ ratios and corresponding hts, obtained from Marianna Shepard
c 6/95. This is used in po1s4 (component 4) calculation below.
c
      parameter(no1op=48)
      real o1op(no1op),o1opht(no1op)
      data o1op
     + /0.1377984,   0.1344115,   0.1306053,   0.1256085,   0.1185497,   
     +  0.1175393,   0.1165371,   0.1155426,   0.1145558,   0.1133328,   
     +  0.1121214,   0.1109214,   0.1092618,   0.1076284,   0.1057982,   
     +  0.1037934,   0.1016451,   0.09900434,  0.09598991,  0.09234733,  
     +  0.08779908,  0.08215878,  0.07523289,  0.07433964,  0.07339706,  
     +  0.07248987,  0.07148170,  0.07047267,  0.06947681,  0.0683660,
     +  0.06724124,  0.06598753,  0.06464776,  0.06320941,  0.0616136,
     +  0.05986323,  0.05803178,  0.05598638,  0.05314200,  0.05069219,  
     +  0.04738312,  0.04701184,  0.04661626,  0.04618508,  0.04574585,  
     +  0.04538802,  0.04496423,  0.04443988/
      data o1opht
     + /100.2, 101.3, 102.6, 104.4, 107.1, 107.5, 107.9, 108.3, 
     +  108.7, 109.2, 109.7, 110.2, 110.9, 111.6, 112.4, 113.3,
     +  114.3, 115.6, 117.2, 119.3, 122.2, 126.6, 134.1, 135.2, 
     +  136.4, 137.6, 139.0, 140.5, 142.1, 144.0, 146.0, 148.3,
     +  150.9, 153.9, 157.4, 161.6, 166.6, 173.0, 181.2, 192.6,
     +  210.1, 212.4, 214.9, 217.5, 220.4, 223.4, 226.8, 230.4/
      data ncalls/0/
c 
      ncalls = ncalls+1
c
c Call solred if necessary:
c
      if (ie5577(3).gt.0.or.ie5577(4).gt.0.or.ie5577(5).gt.0) then
        slt = fslt(dum,ut,glon,1)
        day = float(iyd - iyd/1000*1000)

!       write(6,"('mke5577 call solred: qn2p=',/,(6e12.4))") qn2p
!       write(6,"('mke5577 call solred: qop=',/,(6e12.4))") qop
!       write(6,"('mke5577 call solred: f107d,a=',2e12.4)")
!    |    f107d,f107a

        call solred(ncalls,slt,ht,xo,xo2,xn2,tn,f107d,f107a,
     +    day,glat,glon,id,sr63,qn2p,qop,glyb)

!       write(6,"('mke5577 after solred: glyb=',/,(6e12.4))") glyb

      endif
c
c Column loop:
      do k=1,id
c
c po1s1: o1 recombination (need tn,xo2,xo,xn2)
c The greenline volume emission rates due to the recombination source, which 
c dominates in the 85 to 110 km region at night.
c   (from Ian McDade, (mcdade@windic.yorku.ca) Mar 7, 1994 mail 
c    communication to Roble)
c
c                    A5 k1 [O]^3 {[N2]+[O2]}
c           V1S = -----------------------------      (photons cm-3 sec-1)
c                 {A6 + k5[O2]} {C'[O2] + C"[O]} 
c
c where [O], [O2] and [N2] are number densities cm-3
c A5 = 1.18 sec-1
c A6 = 1.35 sec-1
c k1 = 4.7E-33 (300/T)^2  cm6 sec-1
c k5 = 4.0E-12 exp(-865/T) cm3 sec-1
c C' = 15 dimensionless
c C" = 211 dimensionless
c
c po1s1: o1 recombination
c   (need tn,o,o2,n2)
c
        po1s1 = 0.
        if (ie5577(1).gt.0) then
          rk1 = 4.7e-33*(300./tn(k))**2
          rk5 = 4.0e-12*exp(-865./tn(k))
          po1s1 = (a5 * rk1 * xo(k)**3 * (xn2(k)+xo2(k))) /
     +      ((a6 + rk5*xo2(k)) * (cp*xo2(k) + cpp*xo(k)))
        endif
c
c po1s2: dissociative recombination of o2+
c   (need te,xo2p,xne,xo21d)
c
        po1s2 = 0.
        if (ie5577(2).gt.0) then
          if (te(k).ge.1200.) rkdr = 1.6e-7 * (te(k)/300.)**(-0.55)
          if (te(k).lt.1200.) rkdr = 1.95e-7 * (te(k)/300.)**(-0.7)
          qyld = 0.12	! average
!         qyld = 0.04	! 6/28: experimental lowest value
          po1s2 = qyld*rkdr*xo2p(k)*xne(k)
        endif
c
c po1s3: N2(A)+O energy transfer (this comment changed 4/12/12)
c   (need tn,ht,xo2,xo,xn2,xo21d and qn2p from solred)
c
        po1s3 = 0.
        if (ie5577(3).gt.0) then
          xn2p = 0.3896*qn2p(k) / 
     +      (2.8e-11*xo(k)+2.3e-11*xo2(k)+0.77)
          po1s3 = 0.75*2.8e-11*xo(k)*xn2p
        endif
c
c po1s4: photoelectron impact (from Marianna Shepard)
c   (need tn,ht,xo2,xo,xn2,xo21d and qop from solred)
c o1oprat is O(1s)/O+ ratio corresponding to height, as given by
c   Marianna Shepard (she ran msis90 to get these numbers):
c
        po1s4 = 0.
        if (ie5577(4).gt.0) then
          if (ht(k).lt.o1opht(1)) then
            o1oprat = 0.
          elseif (ht(k).gt.o1opht(no1op)) then
            o1oprat = o1op(no1op)
          else
            do kk=1,no1op-1
              if (ht(k).ge.o1opht(kk).and.ht(k).le.o1opht(kk+1)) then
                o1oprat = o1op(kk)
                if (o1opht(kk+1)-ht(k).lt.ht(k)-o1opht(kk)) 
     +            o1oprat = o1op(kk+1)
              endif
            enddo 
          endif
          po1s4 = o1oprat*qop(k)
        endif
c
c po1s5: photo dissociation of o2 by solar lyman-beta:
c   (glby returned by solred)
c
        po1s5 = 0.
        if (ie5577(5).gt.0) po1s5 = glyb(k)
c
c Sum the components:
c
        if (ie5577(2).gt.0.or.ie5577(3).gt.0.or.ie5577(4).gt.0.or.
     +      ie5577(5).gt.0) then
          rk12 = 4.e-12*exp(-865.*tn(k))
          fout(k) = (po1s2 + po1s3 + po1s4 + po1s5)*1.18 /
     +      (1.35+rk12*xo2(k)+2.e-14*xo(k)+2.6e-10*xo21d(k)) + po1s1

!         write(6,"('mke5577: k=',i3,' po1s2,s3,s4,s5=',4e12.4)")
!    |      k,po1s2,po1s3,po1s4,po1s5

        elseif (ie5577(1).gt.0) then
          fout(k) = po1s1
        endif
      enddo
      return
      end
c
c------------------------------------------------------------------
c
      subroutine mkeoh83(tn,xo2,xo,xn2,id,fout)	! without ohrad
c
c 6/95: Return Oh emission v(8,3) band
c
      real tn(id),xo2(id),xo(id),xn2(id),fout(id)
c     real xh(id),xo3(id),xho2(id),xoh(id),rho(id)
      parameter (f8=0.29)
c
      do k=1,id
        pk6n2 = 5.70e-34*(300./tn(k))**2.62
        pk6o2 = 5.96e-34*(300./tn(k))**2.37
        fout(k) = f8*xo(k)*xo2(k)*(pk6n2*xn2(k)+pk6o2*xo2(k)) /
     +            (260.+2.e-11*xo2(k))
      enddo
      return
      end
c
c------------------------------------------------------------------
c
      subroutine mkeco215(tn,xo,xco2,id,fout)
c
c Return 15 micron co2 emission (from John Wise):
c
c N (15 mic) =  (5.94E-26) Sqrt[Tk] exp[-960/Tk] [O] [CO2]
c               -----------------------------------------
c                    (4 Pi) (1.28 + 3.5E-13 Sqrt[Tk] [O])
c
c The 15 micron term is only the O-CO2 collisional component
c   but it accounts for at least 80% of the radiance above 110 km
c
      real tn(id),xo(id),xco2(id),fout(id)
      parameter (pi=3.14156)
c
      do k=1,id
        fout(k) = (5.94e-26*sqrt(tn(k))*exp(-960./tn(k))*xo(k)*xco2(k))/ 
     +            (4.*pi*(1.28+3.5e-13*sqrt(tn(k))*xo(k)))
      enddo
      return
      end
c
c------------------------------------------------------------------
c
      subroutine mkeno53(tn,xo,xno,id,fout)
c
c Return 5.3 micron NO emission (from John Wise):
c
c N (5.3 mic) = (2.63E-22) exp[-2715/Tk] [O] [NO]
c               ---------------------------------
c                 (4 Pi) (10.78 + 6.5E-11 [O]) 
c
      real tn(id),xo(id),xno(id),fout(id)
      parameter (pi=3.14156)
c
      do k=1,id
        fout(k) = (2.63e-22*exp(-2715./tn(k))*xo(k)*xno(k)) /
     +            (4.*pi*(10.78+6.5e-11*xo(k))) 
      enddo
      return
      end
c-----------------------------------------------------------------
      subroutine mkeo200(tn, h, xo2, xo, xn2, xo1d, id, sza,
     +  ieo200, fout)
c------------------------------------------------------------------
c
c   Calculate Atmospheric Band O2(b1-Sigma -> X3-Sigma)(0-0) emission:
c
c   Up to three components can be requested:
c
c   ieo200(1) > 0 -> Production via Barth-type 2-step process:
c                    O + O + M -> O2* + M
c                    O2* + O2 -> O2(b1-Sigma) + O2
c
c   ieo200(2) > 0 -> O2(b1-Sigma) production via resonance scattering
c                    O2 + hv -> O2(b1-Sigma)
c
c   ieo200(3) > 0 -> O2(b1-Sigma) production via collisional excitation
c                    with O(1D) 
c
c------------------------------------------------------------------
c   Change history:
c
c 8/29/00:  Added dayglow mechanisms - D.Marsh (marsh@ucar.edu)
c
c------------------------------------------------------------------

      use proc, only : p0
      use hist
      implicit none

      integer, intent(in) :: id 
      type(history), intent(in) :: h
      real, dimension(id), intent(in) :: tn
      real, dimension(id), intent(in) :: xo2, xo, xn2, xo1d
      real, intent(in) :: sza
      integer, dimension(3), intent(in) :: ieo200
      real, dimension(id), intent(out) :: fout

c... parameters 
      real, parameter :: a1s =.085
      real, parameter :: fc =.93
      real, parameter :: phi =.95                 ! Green, 2000
      real, parameter :: k2a=3.9e-17, k2b=2.1e-15 ! JPL-97

      real, parameter :: a1=.079, a2=.083
      real, parameter :: rk2a=4.e-17, rk2b=2.2e-15, rk2c=8.e-14
      real, parameter :: dp=6.6, dpp=19.

c... locals
      integer :: k
      real :: rk1, dlev, zp
      real, dimension(id) :: k0, q, pro1d, press, p
      real :: tbar(id-1), o2b(id-1)
      real :: gf(id-1)

      fout = 0.0
 
      dlev = (h%zpt-h%zpb)/float(h%nzp-1)
      do k=1,h%nzp
        zp = h%zpb+(k-1)*dlev
	press(k) = p0*exp(-zp)*1.e-3 ! pressure in mbars
      enddo ! k

c... ieo200(1) > 0 -> Production via Barth-type 2-step process
c   
c   3/9/94: from Ian McDade, (mcdade@windic.yorku.ca) Mar 7, 1994 mail 
c   communication to Roble
c
c   The nighttime O2 Atmospheric (0-0) band emission rates
c   (photons cm-3 sec-1) in the region may be modelled using:
c
c                     A1 k1 [O]^2 {[N2]+[O2]} [O2]
c   Vbarth= --------------------------------------------------
c           {A2 + k2a[O2] + k2b[N2] + k2c[O]} {D'[O2] + D"[O]} 
c
c   where [O], [O2] and [N2] are number densities cm-3
c   A1 = 0.079 sec-1
c   A2 = 0.083 sec-1
c   k1 = 4.7E-33 (300/T)^2  cm6 sec-1
c   k2a = 4.0E-17 cm3 sec-1
c   k2b = 2.2E-15 cm3 sec-1
c   k2c = 8.0E-14 cm3 sec-1
c   D' = 6.6 dimensionless
c   D" = 19  dimensionless
c
c   ref: Planet. Space Sci., 34, 789-800, 1986.
c
      if (ieo200(1).gt.0) then
!        write(6,*) 'ieo200(1)'
        do k=1,id
          rk1 = 4.7e-33*(300./tn(k))**2

          fout(k) = a1*rk1*xo(k)**2 * 
     +      ((xn2(k)+xo2(k)) / (dp*xo2(k)+dpp*xo(k))) *
     +      ( xo2(k) / (a2 + rk2a*xo2(k) + rk2b*xn2(k) + rk2c*xo(k)) )

c          fout(k) = (a1*rk1*xo(k)**2 * (xn2(k)+xo2(k)) * xo2(k)) /
c     +      ((a2 + rk2a*xo2(k) + rk2b*xn2(k) + rk2c*xo(k)) *
c     +      (dp*xo2(k)+dpp*xo(k)))

        enddo
      endif

c... ieo200(2) > 0 -> O2(b1-Sigma) production via resonance scattering
c
c    code supplied by M.Mlynczak
c
    
      gf = 0.0
      if (ieo200(2).gt.0) then
!        write(6,*) 'ieo200(2)'

        ! calculate average temperature and optical mass
        call l13( press, tn, id, sza, tbar, o2b )

        ! calculate g-factor
        call o2egg(tbar, o2b, id, gf)
        gf = gf * 6.192e-9

      endif
 
c... ieo200(3) > 0 -> O2(b1-Sigma) production via collisional excitation
c                    with O(1D) 
c         phi = 0.95 (efficiency for production of b1-sigma)
c         (ref: Green, Shi and Barker, J Phys Chem 104, 6218, 2000)
c         k0 =  3.2(-11)exp(70/T) cm3 s-1 (rate for O(1D) + O2 from JPL97)

      pro1d = 0.0
      if (ieo200(3).gt.0) then

!        write(6,*) 'ieo200(3)'
        k0 =  3.2e-11 * exp(70./tn)
        pro1d = phi * k0 * xo1d 

      endif
 
c... calculate total band VER (photons cm-3 sec-1):
c
c    Vat = Fc Q { g + pro1d } * xo2 + Vbarth
c   
c    where, Fc =  0.93  ( Franck-Condon factor for 0-0 )
c    and Q is the ratio of radiative losses to all losses
c
c          Q =        A1s
c              -----------------------
c              A1s + k2a[O2] + k2b[N2] 

      if (ieo200(2).gt.0 .or. ieo200(3).gt.0) then 
        q = a1s / ( a1s + k2a*xo2 + k2b*xn2 )
        fout = fc * q * ( pro1d + gf ) * xo2 + fout
      endif

      return
      end 


c----------------------------------------------------------------------
        subroutine l13( pmb, t, nz, sza, tbar, o2b )
c----------------------------------------------------------------------
c
c  program to calculate average temperature and optical mass amounts for
c  input into the line-by-line program.  
c 
c----------------------------------------------------------------------

	implicit none

	integer, intent(in) :: nz
	real, dimension(nz), intent(in) :: pmb, t
	real, intent(in) :: sza

	real, dimension(nz), intent(out) :: tbar, o2b

c... locals
	real, parameter :: kb=1.381e-23, pi=3.1415926
	real, dimension(nz) :: p
	real :: cos_sza
	integer :: i
        
!	write(6,*) 'l13'
c... convert pressure to pascals 
	p = pmb*100.0

c... calculate optical mass exactly (see notes 2-18-93)
        if (sza .lt. 89.9) then
	  cos_sza = cos(pi*sza/180.)
        else
	  cos_sza = 0.00001
        endif
	do i=1,nz-1
	  o2b(i)=0.21*2.1188e+22*(p(i)-p(i+1))/100.
	  o2b(i)=o2b(i)/cos_sza
	enddo
	o2b(nz) = 0.21*2.1188e+22*p(nz)/(100*cos_sza)

c... calculate mean temperatures
	do i=1,nz-1
	  tbar(i)=(t(i)+t(i+1))/2.
	enddo
	tbar(nz) = t(nz)

	return

	end subroutine l13
	
c----------------------------------------------------------------------
	subroutine o2egg(tbar, o2b, nz, gf)
c----------------------------------------------------------------------
	use egadb2, only : n, utab, etab

	implicit none

	integer, intent(in) :: nz
	real, dimension(nz), intent(in) :: tbar, o2b
	real, intent(out) :: gf(nz)

c... locals
c... parameters 
	real :: tk(7) = (/150.,175.,200.,225.,250.,275.,300./)
c       real, parameter :: sb=1.946e-22      ! HITRAN 86 bandstrength
	real, parameter :: sb=2.244e-22      ! HITRAN 96 bandstrength

	real :: tb, u, us, utot 
	real :: e1, e2, a1, w1, w2, u1, u2 
	real :: eint(n)
	integer :: i, ij, jlo, jt

c... do the emission calculations
	jlo=1
	utot=0.
	e1=0.
	e2=0.

c... loop over altitude
	do i=nz,1,-1
	  if (tbar(i) > 300) then
	    tb = 300. 
	  else
	    tb = tbar(i)
	  endif

	  u = o2b(i)

c... use liner interpolation to generate emissivity table at tb
	  jt=1
	  call locate(tk,7,tb,jt)
	  a1=(tb-tk(jt))/25.
	  do ij=1,n
	    eint(ij)=etab(ij,jt)+a1*(etab(ij,jt+1)-etab(ij,jt))
	  enddo

c... calculate psuedo mass
	  if(e1.eq.0.) then
	    us=0.
	  else
c... locate u from e1 and interpolate u
  	    jlo=1
	    call locate(eint,n,e1,jlo)
	    u1=utab(jlo)
	    u2=utab(jlo+1)
	    w1=eint(jlo)
	    w2=eint(jlo+1)
	    a1=(e1-w1)/(w2-w1)
	    us=u1+a1*(u2-u1)
!	    write(6,*) 'us', us 
	  endif

	  utot=u+us

c... now get e2 from utot
	  jlo=1
	  call locate(utab,n,utot,jlo)
	  a1=(utot-utab(jlo))/(utab(jlo+1)-utab(jlo))
	  e2=a1*(eint(jlo+1)-eint(jlo))+eint(jlo)
	  gf(i)=(e2-e1)/u/sb
!	  write(6,*) 'gf', gf(i)
	  e1=e2
	enddo

	return

	end subroutine o2egg

c----------------------------------------------------------------------
      subroutine locate(xx,n,x,j)
c----------------------------------------------------------------------
      dimension xx(n)
      jl=0
      ju=n+1
10    if(ju-jl.gt.1)then
        jm=(ju+jl)/2
        if((xx(n).gt.xx(1)).eqv.(x.gt.xx(jm)))then
          jl=jm
        else
          ju=jm
        endif
      go to 10
      endif
      j=jl
      return

      end