gammas.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     #########
      FUNCTION gammas(PX)  RESULT(PGAMMA)
!     ##################################
!
!
!!****  *GAMMAS * -  Gamma  function
!!
!!
!!    PURPOSE
!!    -------
!       The purpose of this function is to compute the Generalized gamma
!    function of its argument.
!
!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!      NONE
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!      None
!!
!!    REFERENCE
!!    ---------
!!      Press, Teukolsky, Vetterling and Flannery: Numerical Recipes, 206-207
!!
!!    AUTHOR
!!    ------
!!      Jean-Pierre Pinty *LA/OMP*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original     7/11/95
!
!*       0. DECLARATIONS
!           ------------
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*       0.1 declarations of arguments and result
!
REAL, INTENT(IN)                     :: px
REAL                                 :: pgamma
!
!*       0.2 declarations of local variables
!
INTEGER                              :: jj ! Loop index
REAL                                 :: zser,zstp,ztmp,zx,zy,zcoef(6)
REAL                                 :: zpi
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('GAMMAS',0,zhook_handle)
zcoef(1) = 76.18009172947146
zcoef(2) =-86.50532032941677
zcoef(3) = 24.01409824083091
zcoef(4) = -1.231739572450155
zcoef(5) =  0.1208650973866179e-2
zcoef(6) = -0.5395239384953e-5
zstp     =  2.5066282746310005
!
zpi = 3.141592654
IF (px.LT.0.) THEN
  zx = 1.- px
ELSE 
  zx = px
END IF
zy = zx
ztmp =  zx + 5.5
ztmp = (zx + 0.5)*alog(ztmp) - ztmp
zser = 1.000000000190015
!
DO jj = 1 , 6
  zy = zy + 1.0
  zser = zser + zcoef(jj)/zy
END DO
!
IF (px.LT.0.) THEN
  pgamma = zpi/sin(zpi*px)/exp( ztmp + alog( zstp*zser/zx ) )
ELSE
  pgamma = exp( ztmp + alog( zstp*zser/zx ) )
END IF
IF (lhook) CALL dr_hook('GAMMAS',1,zhook_handle)
!
END FUNCTION gammas