module spehox

!----------------------------------------------------------------------- 
!
! BOP
!
! !MODULE: spedata
!
! !DESCRIPTION
! Determines the HOx production factor assoctioned with
! solar proton ionization.
!
! !USES
  use shr_kind_mod, only: r8 => shr_kind_r8,r4 => shr_kind_r4
  use cam_abortutils,   only: endrun

  implicit none

  private  ! all unless made public
  save

  public hox_prod_factor

! !REVISION HISTORY:
!   17 Nov 2005  Francis Vitt     Creation
!
! EOP
!----------------------------------------------------------------------- 
! $Id: spehox.F90,v 1.1.2.1 2006/05/03 20:53:09 stacy Exp $
! $Author: stacy $
!----------------------------------------------------------------------- 
! HOx production per ion pair (cm^-3 s-1) from Figure 2 
! of Solomon et al. (1981)
! 
! Data source: 
!           Dr. Charles H. Jackman, Code 916                       
!           NASA/Goddard Space Flight Center                       
!                Greenbelt, MD  20771                              
!  ph:301-614-6053  fx:301-614-5903  Charles.H.Jackman@nasa.gov    
! 
!   Alt(km)     10        100      1000      10000    100000
!  
!    40.0      2.00      2.00      2.00      1.99      1.99
!    42.5      2.00      2.00      1.99      1.99      1.99
!    45.0      2.00      2.00      1.99      1.99      1.99
!    47.5      2.00      2.00      1.99      1.99      1.98
!    50.0      2.00      1.99      1.99      1.98      1.98
!    52.5      2.00      1.99      1.99      1.98      1.95
!    55.0      2.00      1.99      1.98      1.97      1.93
!    57.5      2.00      1.99      1.98      1.95      1.89
!    60.0      1.99      1.98      1.97      1.94      1.85
!    62.5      1.99      1.98      1.96      1.90      1.81
!    65.0      1.99      1.98      1.94      1.87      1.77
!    67.5      1.98      1.96      1.91      1.82      1.72
!    70.0      1.98      1.94      1.87      1.77      1.64
!    72.5      1.96      1.90      1.80      1.70      1.50
!    75.0      1.93      1.84      1.73      1.60      1.30
!    77.5      1.84      1.72      1.60      1.40      0.93
!    80.0      1.60      1.40      1.20      0.95      0.40
!    82.5      0.80      0.60      0.40      0.15      0.00
!    85.0      0.30      0.15      0.10      0.00      0.00
!    87.5      0.00      0.00      0.00      0.00      0.00
!    90.0      0.00      0.00      0.00      0.00      0.00
!----------------------------------------------------------------------- 

  integer, parameter :: nalts = 21
  integer, parameter :: nprods = 5
  real(r8) :: alts(nalts)
  real(r8) :: log_ion_prod(nprods)
  real(r8) :: factor_tbl(nalts,nprods)

  data alts(1:21) / 40.0_r8, 42.5_r8, 45.0_r8, 47.5_r8, 50.0_r8, 52.5_r8, 55.0_r8, 57.5_r8, &
                    60.0_r8, 62.5_r8, 65.0_r8, 67.5_r8, 70.0_r8, 72.5_r8, 75.0_r8, 77.5_r8, &
                    80.0_r8, 82.5_r8, 85.0_r8, 87.5_r8, 90.0_r8 / 

  data log_ion_prod(1:5)  /  1._r8,    2._r8,    3._r8,    4._r8,    5._r8 /

  data factor_tbl( 1,1:5) /  2.00_r8,  2.00_r8,  2.00_r8,  1.99_r8,  1.99_r8 /
  data factor_tbl( 2,1:5) /  2.00_r8,  2.00_r8,  2.00_r8,  1.99_r8,  1.99_r8 /
  data factor_tbl( 3,1:5) /  2.00_r8,  2.00_r8,  2.00_r8,  1.99_r8,  1.99_r8 /
  data factor_tbl( 4,1:5) /  2.00_r8,  2.00_r8,  2.00_r8,  1.99_r8,  1.99_r8 /
  data factor_tbl( 5,1:5) /  2.00_r8,  2.00_r8,  2.00_r8,  1.99_r8,  1.99_r8 /
  data factor_tbl( 6,1:5) /  2.00_r8,  1.99_r8,  1.99_r8,  1.98_r8,  1.95_r8 /
  data factor_tbl( 7,1:5) /  2.00_r8,  1.99_r8,  1.98_r8,  1.97_r8,  1.93_r8 /
  data factor_tbl( 8,1:5) /  2.00_r8,  1.99_r8,  1.98_r8,  1.95_r8,  1.89_r8 /
  data factor_tbl( 9,1:5) /  1.99_r8,  1.98_r8,  1.97_r8,  1.94_r8,  1.85_r8 /
  data factor_tbl(10,1:5) /  1.99_r8,  1.98_r8,  1.96_r8,  1.90_r8,  1.81_r8 /
  data factor_tbl(11,1:5) /  1.99_r8,  1.98_r8,  1.94_r8,  1.87_r8,  1.77_r8 /
  data factor_tbl(12,1:5) /  1.98_r8,  1.96_r8,  1.91_r8,  1.82_r8,  1.72_r8 /
  data factor_tbl(13,1:5) /  1.98_r8,  1.94_r8,  1.87_r8,  1.77_r8,  1.64_r8 /
  data factor_tbl(14,1:5) /  1.96_r8,  1.90_r8,  1.80_r8,  1.70_r8,  1.50_r8 /
  data factor_tbl(15,1:5) /  1.93_r8,  1.84_r8,  1.73_r8,  1.60_r8,  1.30_r8 /
  data factor_tbl(16,1:5) /  1.84_r8,  1.72_r8,  1.60_r8,  1.40_r8,  0.93_r8 /
  data factor_tbl(17,1:5) /  1.60_r8,  1.40_r8,  1.20_r8,  0.95_r8,  0.40_r8 /
  data factor_tbl(18,1:5) /  0.80_r8,  0.60_r8,  0.40_r8,  0.15_r8,  0.00_r8 /
  data factor_tbl(19,1:5) /  0.30_r8,  0.15_r8,  0.10_r8,  0.00_r8,  0.00_r8 /
  data factor_tbl(20,1:5) /  0.00_r8,  0.00_r8,  0.00_r8,  0.00_r8,  0.00_r8 /
  data factor_tbl(21,1:5) /  0.00_r8,  0.00_r8,  0.00_r8,  0.00_r8,  0.00_r8 /

contains

!----------------------------------------------------------------------- 
! Returns the HOx production factor for the ionization profile
!----------------------------------------------------------------------- 
  function hox_prod_factor( ion_pairs, zmid )

    use ppgrid,  only : pcols, pver
    ! for each level interpolate the factor table

    implicit none

    real(r8),intent(in) :: ion_pairs(pver)
    real(r8),intent(in) :: zmid(pver)
    real(r8) :: hox_prod_factor(pver)

    integer :: k
    integer :: lastk


    lastk = 1

    ! start at the bottom since the table goes from bottom to top
    do k = pver,1,-1
       hox_prod_factor(k) = interp_factor_tbl( ion_pairs(k), zmid(k), lastk )
    enddo

  end function hox_prod_factor

!----------------------------------------------------------------------- 
! bilinear interpolates the above table of factors
!----------------------------------------------------------------------- 
  function interp_factor_tbl( ionp, z, lastk )

    implicit none

    real(r8),intent(in) :: ionp
    real(r8),intent(in) :: z
    integer, intent(inout) :: lastk
    real(r8) :: interp_factor_tbl

    integer :: i,j
    real(r8) :: logp
    real(r8) :: atlwgt1,atlwgt2
    real(r8) :: prodwgt1,prodwgt2
    real(r8) :: fact(nprods)

    if (ionp <= 0._r8 ) then 
       interp_factor_tbl = 0._r8
       return
    endif

    ! interpolate log10 of the ionization rate since the table is 
    ! on a log scale.
    logp = log10(ionp)

    if ( z <= alts(1) ) then 
       fact(:) = factor_tbl(1,:)
    else if ( z >= alts(nalts) ) then
       fact(:) = factor_tbl(nalts,:)
    else
       do i = lastk,nalts
          if ( z > alts(i) .and. z <= alts(i+1) ) then 
             atlwgt1 = (alts(i+1) - z)/(alts(i+1) - alts(i))
             atlwgt2 =   (z - alts(i))/(alts(i+1) - alts(i))
             fact(:) = atlwgt1*factor_tbl(i,:) + atlwgt2*factor_tbl(i+1,:)
             lastk = i
             exit
          endif
       enddo
    endif

    if ( logp <= log_ion_prod(1) ) then 
       interp_factor_tbl = fact(1)
    else if ( logp >= log_ion_prod(nprods) ) then
       interp_factor_tbl = fact(nprods)
    else 
       do i = 1,nprods
          if ( logp > log_ion_prod(i) .and. logp <= log_ion_prod(i+1) ) then 
             prodwgt1 = (log_ion_prod(i+1) - logp)/(log_ion_prod(i+1) - log_ion_prod(i))
             prodwgt2 =   (logp - log_ion_prod(i))/(log_ion_prod(i+1) - log_ion_prod(i))
             interp_factor_tbl = prodwgt1*fact(i) + prodwgt2*fact(i+1)
             exit
          endif
       enddo
    endif

  end function interp_factor_tbl

end module spehox
