scale_specfunc.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_specfunc
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: sf_gamma
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  ! Lanczos approximation
  ! g = 7.0, N = 8
  function sf_gamma(x)
    implicit none
 
    real(rp), intent(in) :: x
    real(rp)             :: sf_gamma
 
    real(dp), parameter :: lrtp = 0.9189385332046727_dp ! log(sqrt(2*PI))
    real(dp)            :: ag, dx, work
    !---------------------------------------------------------------------------
 
    dx = real(x-1.0_rp, kind=dp)
 
    ag = 0.999999999999809932276847004735_dp                     &
       + 676.520368121885098567009190444_dp    / ( dx + 1.0_dp ) &
       - 1259.13921672240287047156078755_dp    / ( dx + 2.0_dp ) &
       + 771.323428777653078848652825889_dp    / ( dx + 3.0_dp ) &
       - 176.615029162140599065845513540_dp    / ( dx + 4.0_dp ) &
       + 12.5073432786869048144589368533_dp    / ( dx + 5.0_dp ) &
       - 0.138571095265720116895547069851_dp   / ( dx + 6.0_dp ) &
       + 9.98436957801957085956266899569e-6_dp / ( dx + 7.0_dp ) &
       + 1.50563273514931155833835577539e-7_dp / ( dx + 8.0_dp )
 
    work = exp( lrtp + (dx+0.5_dp)*log(dx+7.5_dp) - (dx+7.5_dp) + log(ag) )
    sf_gamma = real(work, kind=rp)
 
    return
  end function sf_gamma
 
end module scale_specfunc