dgam.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.
!######################################################################
!######################################################################
!
! THIS ROUTINE COMPUTE THE INCOMPLETE AND THE COMPLEMENTARY
! INCOMPLETE GAMMA FUNCTION REAL*8 IN BRIEF, THIS
! FUNCTION HELP TO COMPUTE THE DIFFERENT FRACTION PRESENT IN THE
! TOPMODEL FRAMEWORK.
! MORE EXPLANATION ON THE SUBROUTINE USE ARE GIVEN IN THE NEXT
! COMMENTARY
!
! THIS ROUTINE WAS FOUND ON http://www.netlib.org WEB SITE.
!
! REFERENCE - W. GAUTSCHI, :: A COMPUTATIONAL PROCEDURE FOR INCOMPLETE
! GAMMA FUNCTIONS, ACM TRANS. MATH. SOFTWARE., (1979) 482-489
!
!######################################################################
!######################################################################
!
!     LET  GAMMA(A)  DENOTE THE GAMMA FUNCTION AND  GAM(A,X)  THE
! (COMPLEMENTARY) INCOMPLETE GAMMA FUNCTION,
!
!    GAM(A,X)=INTEGRAL FROM T=X TO T=INFINITY OF EXP(-T)*T**(A-1).
!
! LET  GAMSTAR(A,X)  DENOTE TRICOMI:S FORM OF THE INCOMPLETE GAMMA
! FUNCTION, WHICH FOR A.GT.0. IS DEFINED BY
!
!  GAMSTAR(A,X)=(X**(-A)/GAMMA(A))*INTEGRAL FROM T=0 TO T=X OF
!    EXP(-T)*T**(A-1),
!
! AND FOR A.LE.0. BY ANALYTIC CONTINUATION. FOR THE PURPOSE OF
! THIS SUBROUTINE, THESE FUNCTIONS ARE NORMALIZED AS FOLLOWS&
!
! GAM(A,X)/GAMMA(A),  IF A.GT.0.,
!     G(A,X)=
! EXP(X)*X**(-A)*GAM(A,X),  IF A.LE.0.,
!
!     GSTAR(A,X)=(X**A)*GAMSTAR(A,X)
!   =(1/GAMMA(A))*INTEGRAL FROM T=0 TO T=X OF EXP(-T)*T**(A-1)
!
! THE PROGRAM BELOW ATTEMPTS TO EVALUATE  G(A,X)  AND  GSTAR(A,X),
! BOTH TO AN ACCURACY OF ACC SIGNIFICANT DECIMAL DIGITS, FOR ARBI-
! TRARY REAL VALUES OF  A  AND NONNEGATIVE VALUES OF  X. THE SUB-
! ROUTINE AUTOMATICALLY CHECKS FOR UNDERFLOW AND OVERFLOW CONDI-
! TIONS AND RETURNS APPROPRIATE WARNINGS THROUGH THE OUTPUT PARA-
! METERS  IFLG, IFLGST. A RESULT THAT UNDERFLOWS IS RETURNED WITH
! THE VALUE 0., ONE THAT OVERFLOWS WITH THE VALUE OF THE LARGEST
! MACHINE-REPRESENTABLE NUMBER.
!
!     NEAR LINES IN THE (A,X)-PLANE, A.LT.0., ALONG WHICH  GSTAR
! VANISHES, THE ACCURACY SPECIFIED WILL BE ATTAINED ONLY IN TERMS
! OF ABSOLUTE ERROR, NOT RELATIVE ERROR. THERE ARE OTHER (RARE)
! INSTANCES IN WHICH THE ACCURACY ATTAINED IS SOMEWHAT LESS THAN
! THE ACCURACY SPECIFIED. THE DISCREPANCY, HOWEVER, SHOULD NEVER
! EXCEED ONE OR TWO (DECIMAL) ORDERS OF ACCURACY# NO INDICATION
! OF THIS IS GIVEN THROUGH ERROR FLAGS.
!
!     PARAMETER LIST&
!
!  A - - THE FIRST ARGUMENT OF G AND GSTAR. TYPE REAL.
!  X - - THE SECOND ARGUMENT OF G AND GSTAR. TYPE REAL.
!ACC - - THE NUMBER OF CORRECT SIGNIFICANT DECIMAL DIGITS
!  DESIRED IN THE RESULTS. TYPE REAL.
!  G - - AN OUTPUT VARIABLE RETURNING THE VALUE OF G(A,X).
!  TYPE REAL.
!    GSTAR - - AN OUTPUT VARIABLE RETURNING THE VALUE OF
!  GSTAR(A,X). TYPE REAL.
!     IFLG - - AN ERROR FLAG INDICATING A NUMBER OF ERROR CONDI-
!  TIONS IN G UPON EXECUTION. TYPE INTEGER.
!   IFLGST - - AN ERROR FLAG INDICATING A NUMBER OF ERROR CONDI-
!  TIONS IN GSTAR UPON EXECUTION. TYPE INTEGER.
!  THE VALUES OF IFLG AND IFLGST HAVE THE FOLLOWING
!  MEANINGS&
!  0 - NO ERROR CONDITION.
!  1 - ILLEGAL NEGATIVE ARGUMENT X. THE ROUTINE EXITS
!WITH THE VALUES ZERO FOR G AND GSTAR.
!  2 - INFINITELY LARGE RESULT AT X=0. THE ROUTINE
!RETURNS THE LARGEST MACHINE-REPRESENTABLE NUMBER.
!  3 - (ONLY FOR IFLGST) GSTAR IS INDETERMINATE AT
!A=0. AND X=0. THE ROUTINE RETURNS THE VALUE 1.,
!WHICH IS THE LIMIT VALUE AS X TENDS TO ZERO FOR
!FIXED A=0.
!  4 - THE RESULT UNDERFLOWS. IT IS SET EQUAL TO 0.
!  5 - THE RESULT OVERFLOWS. IT IS SET EQUAL TO THE
!LARGEST MACHINE-REPRESENTABLE NUMBER, WITH
!PROPER SIGN.
!  6 - CONVERGENCE FAILS WITHIN 600 ITERATIONS, EITHER
!IN TAYLOR:S SERIES OR IN LEGENDRE:S CONTINUED
!FRACTION. REASON UNKNOWN. THE COMPUTATION IS
!ABORTED AND THE ROUTINE RETURNS WITH ZERO
!VALUES FOR G AND GSTAR.
!
!     ALL MACHINE-DEPENDENT PARAMETERS ARE COLLECTED IN THE FIRST
! DATA DECLARATION. THEY ARE AS FOLLOWS&
!
! IN THE PROGRAM BELOW THESE PARAMETERS ARE SET TO CORRESPOND TO
! THE MACHINE CHARACTERISTICS OF THE CDC 6500 COMPUTER.
!
!     THE SECOND DATA DECLARATION CONTAINS THE SINGLE PRECISION
! VALUE OF ALOG(10.). THE NEXT DATA DECLARATION CONTAINS THE SUCCES-
! SIVE COEFFICIENTS IN THE MACLAURIN EXPANSION OF (1/GAMMA(A+1))-1.
! THEY ARE GIVEN TO AS MANY DECIMAL PLACES AS IS NECESSARY TO ACHIEVE
! MACHINE PRECISION (ON THE CDC 6500 COMPUTER) IN THE RANGE
! ABS(A).LE..5. MORE ACCURATE VALUES OF THESE COEFFICIENTS (TO
! 31 DECIMAL PLACES) CAN BE FOUND IN TABLE 5 OF J.W.WRENCH,JR.,
! CONCERNING TWO SERIES FOR THE GAMMA FUNCTION, MATH. COMPUT.
! 22, 1968, 617-626.
!
!     THE SUBROUTINE CALLS ON A FUNCTION SUBROUTINE, NAMED  ALGA,
! WHICH IS TO PROVIDE SINGLE PRECISION VALUES OF THE LOGARITHM OF
! THE GAMMA FUNCTION FOR ARGUMENTS LARGER THAN OR EQUAL TO .5.
! A POSSIBLE VERSION OF SUCH A FUNCTION SUBROUTINE IS APPENDED
! TO THE PRESENT SUBROUTINE. IT IS TAYLORED TO THE ACCURACY RE-
! QUIREMENTS OF THE CDC 6500 COMPUTER, AND USES RATIONAL APPROXI-
! MATIONS DUE TO CODY AND HILLSTROM (MATH. COMPUT. 21, 1967, 198-
! 203).
!
!     REFERENCE - W. GAUTSCHI, ::A COMPUTATIONAL PROCEDURE FOR
! INCOMPLETE GAMMA FUNCTIONS, ACM TRANS. MATH. SOFTWARE.
!
!----------------------------------------------------
!!    MODIFICATIONS
!!    -------------
!!J.Escobar10/06/2013: replace DOUBLE PRECISION by REAL to handle problem for promotion of real on IBM SP
!----------------------------------------------------
!################################################################
!
SUBROUTINE dgam(PA, PX, PACC, PG, PGSTAR, KFLG, KFLGST)
!
!################################################################

USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb

IMPLICIT NONE

REAL, INTENT(IN) :: PA, PX, PACC
REAL, INTENT(OUT) :: PG, PGSTAR
INTEGER, INTENT(OUT) :: KFLG, KFLGST
!
REAL, PARAMETER :: XALPREC=digits(1.)*log(2.), &
        xtop=log10(huge(1.)), xbot=log10(tiny(1.))
!
REAL :: ZALX, ZALPHA
REAL :: ZAINF, ZEPS, ZEPS1, ZES, ZSGA
REAL :: ZAE, ZAA, ZAP1, ZAM1, ZAEP1
REAL :: ZAEM1, ZFMA, ZAEPS, ZSGAE, ZAAEPS, ZALGP1
REAL :: ZSGGA, ZALGEP1, ZSGGS, ZALGS
REAL :: ZALG, ZSUMM, ZGA, ZY, ZTERM, ZU, ZP, ZQ, ZR, ZV, ZT, ZH
REAL :: ZSGT, ZA1X, ZRHO, ZXPA
REAL :: ZXMA, ZS, ZTEMP
!
INTEGER :: JI, JK, JK1, JMA
LOGICAL :: GRET
!
!RJ REAL,SAVE::AL10=2.3025850929940456840179914547
REAL,DIMENSION(29),PARAMETER :: XC=(/ &
        0.57721566490153286060651209008, &
       -.65587807152025388107701951515, &
      -4.200263503409523552900393488e-2,  &
        0.1665386113822914895017007951,  &
      -4.21977345555443367482083013e-2,   &
      -9.6219715278769735621149217e-3,    &
       7.2189432466630995423950103e-3,    &
      -1.165167591859065112113971e-3,     &
      -2.15241674114950972815730e-4,&
       1.2805028238811618615320e-4, &
      -2.013485478078823865569e-5,  &
      -1.2504934821426706573e-6,    &
       1.1330272319816958824e-6,    &
      -2.056338416977607103e-7,     &
       6.1160951044814158e-9, &
       5.0020076444692229e-9, &
      -1.181274570487020e-9,  &
       1.04342671169110e-10,  &
       7.782263439905e-12,    &
      -3.696805618642e-12,    &
       5.1003702875e-13,&
      -2.058326054e-14, &
      -5.34812254e-15,  &
       1.2267786e-15,   &
      -1.181259e-16,    &
       1.187e-18, &
       1.412e-18, &
      -2.30e-19,  &
       1.7e-20/)

REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!
IF (lhook) CALL dr_hook('DGAM',0,zhook_handle)
!
pg     = 0.0
pgstar = 0.0
!
kflg   = 0
kflgst = 0
!
IF ( px<0. ) THEN
  ! 290 deb
  kflg = 1
  kflgst = 1
  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
  RETURN
  ! 290 end
ENDIF
!
! INITIALIZATION
!
zainf   = huge(1.)
zaa = abs(pa)
!
IF ( px==0. ) THEN
  !
  ! EVALUATION OF GSTAR(A,0.) AND G(A,0.)
  !
  IF (pa<0.) THEN
    kflgst = 2
    pgstar = zainf
    pg     = 1./zaa
  ELSEIF (pa==0.) THEN
    kflgst = 3
    kflg   = 2    
    pgstar = 1.
    pg     = zainf
  ELSEIF (pa>0.) THEN
    pg = 1.
  ENDIF
  !
  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
  RETURN
  !
ENDIF
!
! PX > 0.
!
! INITIALIZATION
!
ji = 0
!
zeps = 0.5*10.**(-pacc)
zeps1 = zeps/100.
!
zalx = log(px)
!
zalpha = px + 0.25
IF ( px<0.25 ) zalpha = log(0.5)/zalx
!
zsga = 1.
IF ( pa<0. ) zsga = -zsga
!
zae   = pa
zap1  = pa + 1.
zaep1 = zap1
!
jma    = int(0.5-pa)
zfma   = float(jma)
zaeps  = pa + zfma
zaaeps = abs(zaeps)
!
zsgae = 1.
IF (zaeps<0.) zsgae = -zsgae
!
! EVALUATION OF THE LOGARITHM OF THE ABSOLUTE VALUE OF
! GAMMA(A+1.) AND DETERMINATION OF THE SIGN OF GAMMA(A+1.)
!
zalgp1 = 0.
zsgga  = 1.
!
IF (jma>0 .AND. zaeps/=0.) THEN

  zsgga = zsgae
  IF ( jma==2*(jma/2) ) zsgga = -zsgga
  zalgp1 = dlga(zaeps+1.) - log(zaaeps)

  IF (jma/=1) THEN
    zalgp1 = zalgp1 + dlga(1.-zaeps) - dlga(zfma-zaeps)
  ENDIF
 
ELSEIF ( jma<=0 ) THEN
  !10 deb
  zalgp1 = dlga(zap1)
  ! 10 end
ENDIF
!
! 20 deb
zalgep1 = zalgp1
!
!Map: 
! PA<=ZALPHA
!   !
!   PX>1.5 
!     PA <0., ==0., >0. 
!     RETURN
!     !
!   PA<-0.5
!     PX<0.25 .AND. ZAE>ZALPHA
!   ELSE
!   !
!   PA <0.5
!   ELSE
!   !
!   CODE120
!   PA <0., ==0., >0.
!   RETURN
!   !
! PA>ZALPHA
!
! 60 deb
IF ( pa<=zalpha ) THEN
  !
  IF ( px>1.5 ) THEN 
    !
    ! EVALUATION OF G(A,X) FOR X.GT.1.5 AND A.LE.ALPHA(X) BY
    ! MEANS OF THE LEGENDRE CONTINUED FRACTION      
    !
    !240 deb
    pgstar = 1.
    zxpa = px + 1. - pa
    zxma = px - 1. - pa
    zp = 0.
    zq = zxpa*zxma
    zr = 4.*zxpa
    zs = -pa + 1.
    zterm = 1.
    zsumm = 1.
    zrho = 0.
    !
    DO jk = 1,600
      !
      zp = zp + zs
      zq = zq + zr
      zr = zr + 8.
      zs = zs + 2.
      zt = zp*(1.+zrho)
      zrho = zt/(zq-zt)
      zterm = zrho*zterm
      zsumm = zsumm + zterm
      !
      IF ( abs(zterm)<=zeps*zsumm ) EXIT
      !
      IF ( jk==600 ) THEN
        ! 330 deb
        kflg = 6
        IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
        RETURN
        ! 330 end
      ENDIF
      !
    ENDDO
    !
    IF ( pa<0. ) THEN
      !
      ! 260 deb
      pg = zsumm/zxpa
      zalg = log(pg)
      ! 260 end
      !
      ! EVALUATION OF GSTAR(A,X) IN TERMS OF G(A,X)
      !  200 deb (200 est la fin de 180)
      CALL code200(jma, px, pa, zaeps, zsga, zsgga, zaa, zalx, &
                   zalg, zalgp1, zainf, pgstar, kflgst, gret)
      !  
      IF ( gret ) THEN
        IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
        RETURN
      ENDIF
      !      
    ELSEIF ( pa==0. ) THEN
      !
      ! 270 deb
      pg = zsumm/zxpa
      ! 270 end
      !
    ELSEIF ( pa>0. ) THEN
      !
      ! 280 deb
      zalg = pa*zalx - px + log(pa*zsumm/zxpa) - zalgp1
      !
      IF ( zalg<=xbot ) THEN
        ! 300 deb
        kflg = 4
        IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
        RETURN
        ! 300 end
      ENDIF
      !
      pg = exp(zalg)
      pgstar = 1. - pg
      ! 280 end
      ! 
    ENDIF 
    !
    IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
    RETURN    
    !
  ENDIF ! end  PX>1.5 
  !
  !
  IF ( pa<-0.5 ) THEN
    !
    ! RECURSIVE EVALUATION OF G(A,X) FOR X.LE.1.5 AND A.LT.-.5
    !
    ! 170 deb
    ji = 1
    zae = zaeps
    zaep1 = zaeps + 1.
    ! 170 end
    !
    IF ( px<0.25 .AND. zae>zalpha ) THEN
      !
      ! 210 deb
      ji = 2
      zalgep1 = dlga(zaep1)
      !
      ! EVALUATION OF GSTAR(A,X) FOR A.GT.ALPHA(X) BY TAYLOR
      ! EXPANSION
      !
      CALL code220(px,zae,zalx,zalgep1,zeps,pg,pgstar,kflgst,gret)
      !
      IF ( gret ) THEN
        IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
        RETURN
      ENDIF
      !
      pg = pg*exp(zalgep1)/zae
      ! 210 end (220 230)
      !
      ! 180 deb
      CALL code180(jma, px, pa, zae, zaeps, zsga, zsgga, zaa, zalx, &
                   zalgp1, zainf, pg, pgstar, kflgst, gret)
      !
      IF ( gret ) THEN
        IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
        RETURN
      ENDIF
      !
    ENDIF ! PX<0.25 .AND. ZAE>ZALPHA
    !
  ELSE ! end PA < -0.5
    !
    ! DIRECT EVALUATION OF G(A,X) AND GSTAR(A,X) FOR X.LE.1.5
    ! AND -.5.LE.A.LE.ALPHA(X)
    !
    pgstar = 1.
    !
  ENDIF
  !
  IF ( pa<0.5 ) THEN 
    !
    ! 70 deb
    zsumm = xc(29)
    !
    DO jk = 1,28
      jk1 = 29 - jk
      zsumm = zae*zsumm + xc(jk1)
    ENDDO
    !
    zga = -zsumm / (1.+zae*zsumm)
    zy  = zae*zalx
    !
    IF ( abs(zy)<1. ) THEN
      !
      zsumm = 1.
      zterm = 1.
      !
      DO jk = 1,600
        !
        zterm = zy * zterm / jk
        zsumm = zsumm + zterm
        !
        IF ( abs(zterm)<=zeps1*zsumm ) EXIT
        !
        IF ( jk==600) THEN
          kflg = 6
          IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
          RETURN
        ENDIF
        !
      ENDDO
      !
      zu = zga - zsumm*zalx
      ! 70 end
      !
    ELSE
      !
      ! 100 deb
      zu = zga - (exp(zy)-1.)/zae
      ! 100 end
      !
    ENDIF
    !
  ELSE
    !
    ! 110 deb
    ztemp = dlga(pa*1.)
    zu = dexp(ztemp) - (px**pa)/pa
    ! 110 end
    !
  ENDIF
  !
  ! 120 deb
  zp = zae*px
  zq = zaep1
  zr = zae + 3.
  zterm = 1.
  zsumm = 1.
  !
  DO jk = 1,600
    !
    zp = zp + px
    zq = zq + zr
    zr = zr + 2.
    zterm = -zp * zterm / zq
    zsumm = zsumm + zterm
    !
    IF ( abs(zterm)<=zeps1*zsumm ) EXIT
    !
    IF ( jk==600) THEN
      kflg = 6
      IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
      RETURN
      !
    ENDIF
    !
  ENDDO
  !
  zv = (px**zaep1) * zsumm / zaep1
  pg = zu + zv
  !
  IF (ji==1) THEN
    !
    ! 180 deb
    CALL code180(jma, px, pa, zae, zaeps, zsga, zsgga, zaa, zalx, &
                 zalgp1, zainf, pg, pgstar, kflgst, gret)
    !
    IF ( gret ) THEN
      IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
      RETURN
    ENDIF
    !
  ENDIF
  !
  IF (pa<0.) THEN
    !
    zt = exp(px) * px**(-pa)
    pg = zt * pg
    pgstar = 1. - pa * pg * dexp(-zalgp1) / zt
    !
  ELSEIF (pa==0.) THEN
    !
    pg = exp(px)*pg
    !
  ELSEIF (pa>0.) THEN
    !
    pg = pa*pg*exp(-zalgp1)
    pgstar = 1. - pg
    !
  ENDIF
  !
  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
  RETURN  
  ! 120 end
  !
ENDIF
!
! EVALUATION OF GSTAR(A,X) FOR A.GT.ALPHA(X) BY TAYLOR
! EXPANSION
!
! PA > ZALPHA : 220 deb
  CALL code220(px,zae,zalx,zalgep1,zeps,pg,pgstar,kflgst,gret)
!
IF ( gret .OR. ji/=2) THEN
  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
  RETURN
ENDIF
!
pg = pg*exp(zalgep1)/zae
! 220 230 end
!
! 180 deb
! 180 deb
 CALL code180(jma, px, pa, zae, zaeps, zsga, zsgga, zaa, zalx, &
              zalgp1, zainf, pg, pgstar, kflgst, gret)
!
CONTAINS
!
!
SUBROUTINE code220(PX,PAE,PALX,PALGEP1,PEPS,PG,PGSTAR,KFLGST,ORET)
!
IMPLICIT NONE
!
REAL, INTENT(IN) :: PX
REAL, INTENT(IN) :: PAE
REAL, INTENT(IN) :: PALX
REAL, INTENT(IN) :: PALGEP1
REAL, INTENT(IN) :: PEPS
REAL, INTENT(OUT) :: PG
REAL, INTENT(INOUT) :: PGSTAR
INTEGER, INTENT(INOUT) :: KFLGST
LOGICAL, INTENT(OUT) :: ORET
!
REAL :: ZTERM, ZSUMM, ZALGS
INTEGER :: JK
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (lhook) CALL dr_hook('DGAM:CODE220',0,zhook_handle)
!
! 220 deb
pg = 1.
zterm = 1.
zsumm = 1.
!
DO jk = 1,600
  !
  zterm = px*zterm/(pae+jk)
  zsumm = zsumm + zterm
  IF (abs(zterm)<=peps*zsumm) EXIT
  !
  IF ( jk==600 ) THEN
    kflgst = 6
    IF (lhook) CALL dr_hook('DGAM:CODE220',1,zhook_handle)
    RETURN
  ENDIF
  !
ENDDO
!
zalgs = pae*palx - px + log(zsumm) - palgep1
!
IF ( zalgs<=xbot ) THEN
  kflgst = 4
  IF (lhook) CALL dr_hook('DGAM:CODE220',1,zhook_handle)
  RETURN
ENDIF
!
pgstar = exp(zalgs)
pg = 1. - pgstar
!
IF (lhook) CALL dr_hook('DGAM:CODE220',1,zhook_handle)
!
END SUBROUTINE code220
!
!
SUBROUTINE code200(KMA, PX, PA, PAEPS, PSGA, PSGGA, PAA, PALX, &
                   PALG, PALGP1, PAINF, PGSTAR, KFLGST, ORET)
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: KMA
REAL, INTENT(IN) :: PX
REAL, INTENT(IN) :: PA
REAL, INTENT(IN) :: PAEPS
REAL, INTENT(IN) :: PSGA
REAL, INTENT(IN) :: PSGGA
REAL, INTENT(IN) :: PAA
REAL, INTENT(IN) :: PALX
REAL, INTENT(IN) :: PALG
REAL, INTENT(IN) :: PALGP1
REAL, INTENT(IN) :: PAINF
REAL, INTENT(OUT) :: PGSTAR
INTEGER, INTENT(INOUT) :: KFLGST
LOGICAL, INTENT(OUT) :: ORET
!
REAL :: ZSGT, ZT
INTEGER :: JK
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (lhook) CALL dr_hook('DGAM:CODE200',0,zhook_handle)
!
oret = .false.
!
! EVALUATION OF GSTAR(A,X) IN TERMS OF G(A,X)
!
pgstar = 1.
!
IF ( kma>=0 .AND. paeps==0. ) THEN
  oret = .true.        
  IF ( lhook ) CALL dr_hook('DGAM:CODE200',1,zhook_handle)
  RETURN
ENDIF
!
zsgt = psga*psgga
zt = log(paa) - px + pa*palx + palg - palgp1
!
IF ( zt<-xalprec .OR. zt>=xtop ) THEN
  !
  IF ( zt>=xtop ) THEN
    kflgst = 5
    pgstar = -zsgt*painf
  ENDIF
  !
  oret = .true.        
  IF ( lhook ) CALL dr_hook('DGAM:CODE200',1,zhook_handle)
  RETURN
  !
ENDIF
!
pgstar = 1. - zsgt*exp(zt)
!
IF (lhook) CALL dr_hook('DGAM:CODE200',1,zhook_handle)
!
END SUBROUTINE code200
!
!
SUBROUTINE code180(KMA, PX, PA, PAE, PAEPS, PSGA, PSGGA, PAA, PALX, &
                   PALGP1, PAINF, PG, PGSTAR, KFLGST, ORET)
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: KMA
REAL, INTENT(IN) :: PX
REAL, INTENT(IN) :: PA
REAL, INTENT(IN) :: PAE
REAL, INTENT(IN) :: PAEPS
REAL, INTENT(IN) :: PSGA
REAL, INTENT(IN) :: PSGGA
REAL, INTENT(IN) :: PAA
REAL, INTENT(IN) :: PALX
REAL, INTENT(IN) :: PALGP1
REAL, INTENT(IN) :: PAINF
REAL, INTENT(INOUT) :: PG
REAL, INTENT(OUT) :: PGSTAR
INTEGER, INTENT(INOUT) :: KFLGST
LOGICAL, INTENT(OUT) :: ORET
!
REAL :: ZALG
INTEGER :: JK
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (lhook) CALL dr_hook('DGAM:CODE180',0,zhook_handle)
!
pg = pg*exp(px)*px**(-pae)
!
DO jk = 1,kma
  pg = (1.-px*pg)/(jk-pae)
ENDDO
!
zalg = log(pg)
!
! EVALUATION OF GSTAR(A,X) IN TERMS OF G(A,X)
!
 CALL code200(kma, px, pa, paeps, psga, psgga, paa, palx, &
              zalg, palgp1, painf, &
              pgstar, kflgst, oret)
!
IF (lhook) CALL dr_hook('DGAM:CODE180',1,zhook_handle)
!
END SUBROUTINE code180
!
!
FUNCTION dlga(PDX) RESULT(PDLGAR)

USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb

IMPLICIT NONE

REAL :: PDLGAR
!
REAL, INTENT(IN) :: PDX
!
REAL :: ZDC, ZDP, ZDY, ZDT, ZDS
!
REAL,DIMENSION(8),PARAMETER :: &
        XDBNUM=(/-3617., 1., -691., 1., -1., 1., -1., 1./)
REAL,DIMENSION(8),PARAMETER :: &
        XDBDEN=(/122400., 156., 360360., 1188., 1680., 1260., 360. , 12./)
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('DGAMM:DLGA',0,zhook_handle)

zdc = 0.5*log(8.*atan(1.))
zdp = 1.
zdy = pdx
zy  = zdy
!
! THE CONDITIONAL CLAUSE IN THE NEXT STATEMENT EXPRESSES THE
! INEQUALITY Y.GT.EXP(.121189*DPREC+.053905), WHERE DPREC IS THE
! NUMBER OF DECIMAL DIGITS CARRIED IN REAL*8 FLOATING-POINT
! ARITHMETIC.
!
!RJ-remark this magic number is for 53 significant bit prec only +/-!!!!
DO WHILE ( zy<=35.027 )
  zdp = zdy*zdp
  zdy = zdy + 1.
  zy  = zdy
ENDDO
!
zdt = 1./(zdy*zdy)
zds = 43867./244188.
!
DO ji=1,8
  zds = zdt*zds + xdbnum(ji)/xdbden(ji)
ENDDO
!
pdlgar = (zdy-0.5)*log(zdy) - zdy + zdc + zds/zdy - log(zdp)
!
IF (lhook) CALL dr_hook('DGAMM:DLGA',1,zhook_handle)
!
END FUNCTION dlga

END SUBROUTINE dgam