!
      module o2src
      use cons_module,only: check_exp
      implicit none
!-----------------------------------------------------------------------
!
! Module to implement parameterization of photolysis and energy 
!   deposition/heating of O2 by solar radiation in the Schumann-Runge 
!   Continuum (do2src)
!
! Code history:
!   Source for SRC received by Ray Roble from Sam Yee, June 2000.
!   Standardized for waccm and tgcm by Ben Foster, Oct 2000.
!
! Reference:
!   DeMajistre, R., Jeng-Hwa Yee, and Xun Zhu, Parameterizations of 
!     Oxygen Photoysis and Heating Rates due to Solar Energy Absorption
!     in the Schumann-Runge Continuum, GRL(?), in press.
!
!-----------------------------------------------------------------------
      real,external :: expo ! util.F
      contains
      subroutine mkdo2src(sco2,f107d,do2src,nlev)
!
! Calculate O2 photolysis (dissociation) rate in Shumann-Runge Continuum:
!
! Input arguments:
      integer,intent(in) :: nlev        ! number of vertical levels 
      real,intent(in)    :: sco2(nlev)  ! O2 slant column density (cm-2)
      real,intent(in)    :: f107d       ! 10.7 cm daily solar flux
!
! Output arguments:
      real,intent(out) :: do2src(nlev)  ! O2 photolysis rate (s-1) in SRC
!
! Local:
      integer :: ic,jc,k
      real :: pjinf,w
      real,parameter :: a = 1.031e-17
      real,parameter :: pji(3) = (/2.118e-6, 5.658e-9, -9.179e-12/)
      real,parameter :: alpha(4,3) = reshape(source =
     |  (/ 6.255e-1,  2.577e-1,  1.082e-1, 4.206e-3,
     |     2.488e-4, -1.659e-4, -7.647e-5, 7.783e-6,
     |    -4.871e-7,  3.577e-7,  1.281e-7, 1.752e-8/),
     |     shape = (/4,3/))
!
      pjinf=0.
      do ic=1,3
        pjinf=pjinf+pji(ic)*(f107d**(ic-1))
      enddo
!
      do k=1,nlev
        do2src(k)=0.
        do ic=1,4
          w=0.
          do jc=1,3
            w = w+alpha(ic,jc)*(f107d**(jc-1))
          enddo
          if (.not.check_exp) then
            do2src(k) = do2src(k)+w*exp(-sco2(k)*a*(4.**(1-ic)))
          else
            do2src(k) = do2src(k)+w*expo(-sco2(k)*a*(4.**(1-ic)),0)
          endif
        enddo
        do2src(k) = do2src(k)*pjinf
      enddo ! k=1,nlev
      return
      end subroutine mkdo2src
!-----------------------------------------------------------------------
      subroutine mkho2src(sco2,xno2,rho,cp,f107d,ho2src,nlev,mks)
!
! Calculate energy deposition rate (heating) from O2 photo dissociation
! in the Schumann-Runge Continuum.
!
! Input arguments:
      integer,intent(in) :: nlev         ! number of vertical levels
      real,intent(in)    :: sco2(nlev)   ! O2 slant column density (cm-2)
      real,intent(in)    :: xno2(nlev)   ! O2 number density (cm-3)
      real,intent(in)    :: rho(nlev)    ! total mass density (gm/cm-3)
      real,intent(in)    :: cp(nlev)     ! specific heat (mks > 0)
      real,intent(in)    :: f107d        ! 10.7 cm daily solar flux
      integer,intent(in) :: mks          ! units flag for ho2src
!
! Output arguments:
!   If mks >  0, ho2src is returned in deg K/sec (mks)
!   If mks <= 0, ho2src is returned in ergs/gm-1/s-1 (cgs)
!
      real,intent(out) :: ho2src(nlev)   ! energy deposition (heating)
!
! Local:
      integer :: ic,jc,k
      real :: einf,u
      real,parameter :: b = 1.094e-17
      real,parameter :: ei(3) = (/3.179e-18, 9.979e-21, -1.609e-23/)
      real,parameter :: beta(4,3) = reshape(source =
     |  (/7.932e-1,  1.702e-1,  2.467e-1,  1.440e-3,
     |    1.291e-4, -1.540e-4,  3.988e-5, -8.463e-6,
     |   -3.110e-7,  3.815e-7, -1.092e-7,  2.474e-8/),
     |    shape = (/4,3/))
!
      einf=0.
      do ic=1,3
        einf = einf+ei(ic)*(f107d**(ic-1))
      enddo
!
      do k=1,nlev
        ho2src(k) = 0.
        do ic=1,4
          u=0.
          do jc=1,3
            u = u+beta(ic,jc)*(f107d**(jc-1))
          enddo
          if (.not.check_exp) then
            ho2src(k) = ho2src(k)+u*exp(-sco2(k)*B*(4.**(1-ic)))
          else
            ho2src(k) = ho2src(k)+u*expo(-sco2(k)*B*(4.**(1-ic)))
          endif
        enddo
        ho2src(k) = ho2src(k)*einf
      enddo ! k=1,nlev
!
! cgs units: ergs/gm-1/s-1
      do k=1,nlev
        ho2src(k) = xno2(k) * ho2src(k) / rho(k)
!
! mks units: deg K/s
        if (mks > 0) ho2src(k) = ho2src(k) / cp(k)
      enddo
      return
      end subroutine mkho2src
!-----------------------------------------------------------------------
      end module o2src