module rayscat

  use vtype_module

  implicit none

  private

  public :: CalcRayScatCoef


  !----------------------------------------------------------------------------

contains

  !----------------------------------------------------------------------------

  subroutine CalcRayScatCoef( &
    & MolNum, Wavelength, &
    & RayScatCoef &
    & )

    use constants0

    integer     , intent(in ) :: MolNum
    real(dp)    , intent(in ) :: Wavelength
    real(dp)    , intent(out) :: RayScatCoef


    ! Local variables
    !

    ! Pierrehumbert (2010), p.332, Table 5.2
    real(DP), parameter :: RayScatCoefRatio00Air =  4.4459d0
    real(DP), parameter :: RayScatCoefRatio01H2O =  3.3690d0
    real(DP), parameter :: RayScatCoefRatio02CO2 = 10.5611d0
    real(DP), parameter :: RayScatCoefRatio06CH4 = 10.1509d0
    real(DP), parameter :: RayScatCoefRatio07O2  =  3.8634d0
    real(DP), parameter :: RayScatCoefRatio22N2  =  4.6035d0
    real(DP), parameter :: RayScatCoefRatio11NH3 =  7.3427d0
    real(DP), parameter :: RayScatCoefRatio45H2  =  1.0d0
    real(DP), parameter :: RayScatCoefRatioXXHe  =  0.0641d0

    real(DP) :: ConvFactor


!!$    call CalcRayScatCoefL2002Air( &
!!$      & Wavelength, &
!!$      & RayScatCoef &
!!$      & )
    call CalcRayScatCoefHT1974Air( &
      & Wavelength, &
      & RayScatCoef &
      & )

    select case ( MolNum )
    case (  0 ) ! Air
      ConvFactor = 1.0d0
    case (  1 ) ! H2O
      ConvFactor = RayScatCoefRatio01H2O / RayScatCoefRatio00Air
    case (  2 ) ! CO2
      ConvFactor = RayScatCoefRatio02CO2 / RayScatCoefRatio00Air
      !
      call CalcRayScatCoefHT1974CO2( &
        & Wavelength, &
        & RayScatCoef &
        & )
      ConvFactor = 1.0d0
    case (  6 ) ! CH4
      ConvFactor = RayScatCoefRatio06CH4 / RayScatCoefRatio00Air
    case (  7 ) ! O2
      ConvFactor = RayScatCoefRatio07O2  / RayScatCoefRatio00Air
    case ( 11 ) ! NH3
      ConvFactor = RayScatCoefRatio11NH3 / RayScatCoefRatio00Air
    case ( 22 ) ! N2
      ConvFactor = RayScatCoefRatio22N2  / RayScatCoefRatio00Air
    case ( 45 ) ! H2
      ConvFactor = RayScatCoefRatio45H2  / RayScatCoefRatio00Air
    case ( 99 ) ! He
      ConvFactor = RayScatCoefRatioXXHe  / RayScatCoefRatio00Air
    case default
      ConvFactor = 1.0d0
    end select

    RayScatCoef = ConvFactor * RayScatCoef


  end subroutine CalcRayScatCoef

  !----------------------------------------------------------------------------

  subroutine CalcRayScatCoefL2002Air( &
    & Wavelength, &
    & RayScatCoef &
    & )

    use constants0

    real(dp)    , intent(in ) :: Wavelength
    real(dp)    , intent(out) :: RayScatCoef


    ! Local variables
    !
    real(dp)            :: RefIndexRealMinusOne
    real(dp), parameter :: Delta = 0.035d0
    real(dp)            :: AnIsoProp
    real(dp), parameter :: NumDenSTP = 101325.0d0 / ( Boltz * 273.15d0 )


    ! Liou (2002), p.93
    RefIndexRealMinusOne = &
      &   6432.8d0 &
      & + 2949810.0d0 / ( 146.0d0 - ( WaveLength*1.0d6 )**(-2) ) &
      & +   25540.0d0 / (  41.0d0 - ( WaveLength*1.0d6 )**(-2) )
    RefIndexRealMinusOne = RefIndexRealMinusOne * 1.0d-8
    AnIsoProp = ( 6.0d0 + 3.0d0 * Delta ) / ( 6.0d0 - 7.0d0 * Delta )
    RayScatCoef = &
      & 8.0d0 * PI**3 &
      &   * ( RefIndexRealMinusOne * ( RefIndexRealMinusOne + 2.0d0 ) )**2 &
      &   / ( 3.0d0 * WaveLength**4 * NumDenSTP**2 ) &
      &   * AnIsoProp


  end subroutine CalcRayScatCoefL2002Air

  !----------------------------------------------------------------------------

  subroutine CalcRayScatCoefI2008CO2( &
    & Wavelength, &
    & RayScatCoef &
    & )

    use constants0

    real(dp)    , intent(in ) :: Wavelength
    real(dp)    , intent(out) :: RayScatCoef


    ! Local variables
    !
    real(dp)            :: WavelengthCM
    real(dp)            :: RefIndexRealMinusOne
    real(dp)            :: AnIsoProp
    real(dp), parameter :: NumDenSTP = 101325.0d0 / ( Boltz * 273.15d0 )


    WaveLengthCM = WaveLength * 1.0d2

    ! Ityaksov et al. (2008) Chem Phys. Lett., 462, 31-34.
    RefIndexRealMinusOne = &
      &   1.1427d6 * WaveLengthCM**2 &
      & * ( &
      &       5799.3d0   / ( 16.6d9  * WaveLengthCM**2 - 1.0d0 ) &
      &     +  120.0d0   / (  7.96d9 * WaveLengthCM**2 - 1.0d0 ) &
      &     +    5.33d0  / (  5.63d9 * WaveLengthCM**2 - 1.0d0 ) &
      &     +    4.32d0  / (  4.6d9  * WaveLengthCM**2 - 1.0d0 ) &
      &     +    1.22d-5 / (  5.85d6 * WaveLengthCM**2 - 1.0d0 ) &
      &   )
    AnIsoProp = 1.14d0 + 25.3d-12 / WaveLengthCM**2
    RayScatCoef = &
      & 24.0d0 * PI**3 &
      &   / ( WaveLengthCM**4 * (NumDenSTP*1d-6)**2 ) &
      &   * (   ( RefIndexRealMinusOne*(RefIndexRealMinusOne+2.0d0)         )      &
      &       / ( RefIndexRealMinusOne*(RefIndexRealMinusOne+2.0d0) + 3.0d0 ) )**2 &
      &   * AnIsoProp &
      &   * 1.0d-4


  end subroutine CalcRayScatCoefI2008CO2

  !----------------------------------------------------------------------------
  ! Hansen and Travis (1974)
  !----------------------------------------------------------------------------

  subroutine CalcRayScatCoefHT1974Air( &
    & Wavelength, &
    & RayScatCoef &
    & )

    use constants0

    real(dp)    , intent(in ) :: Wavelength
    real(dp)    , intent(out) :: RayScatCoef


    ! Local variables
    !
    real(dp), parameter :: Press0    = 1013.25d2
    real(dp), parameter :: Grav      = 9.8d0
    real(dp), parameter :: MeanMolWt = 28.971d0
    real(dp)            :: WavelengthMicron
    real(dp)            :: OptDep


    WavelengthMicron = Wavelength * 1.0d6

    ! Hansen and Travis (1974) (2.30)
    OptDep = &
      &     0.008569d0 * WavelengthMicron**(-4) &
      & * ( 1.0d0 + 0.0113d0  * WavelengthMicron**(-2) &
      &           + 0.00013d0 * WavelengthMicron**(-4) )
    RayScatCoef = OptDep / ( Press0 / Grav / ( MeanMolWt * AtomMassConst ) )


  end subroutine CalcRayScatCoefHT1974Air

  !----------------------------------------------------------------------------
  ! Hansen and Travis (1974)
  !----------------------------------------------------------------------------

  subroutine CalcRayScatCoefHT1974CO2( &
    & Wavelength, &
    & RayScatCoef &
    & )

    use constants0

    real(dp)    , intent(in ) :: Wavelength
    real(dp)    , intent(out) :: RayScatCoef


    ! Local variables
    !
    real(dp), parameter :: Press0    = 93.0d0 * 1013.25d2
    real(dp), parameter :: Grav      = 8.7d0
    real(dp), parameter :: MeanMolWt = 44.0d0
    real(dp)            :: WavelengthMicron
    real(dp)            :: OptDep


    WavelengthMicron = Wavelength * 1.0d6

    ! Hansen and Travis (1974) (2.32)
    OptDep = &
      &     1.527d0 * WavelengthMicron**(-4) &
      & * ( 1.0d0 + 0.013d0 * WavelengthMicron**(-2) )
    RayScatCoef = OptDep / ( Press0 / Grav / ( MeanMolWt * AtomMassConst ) )


  end subroutine CalcRayScatCoefHT1974CO2

  !----------------------------------------------------------------------------

end module rayscat
