SURFEX v7.3: /home/dasprezs/EXPORT_v7_3/src/SURFEX/dgam.F Source File

Go to the documentation of this file.
00001 C######################################################################
00002 C######################################################################
00003 C
00004 C       THIS ROUTINE COMPUTE THE INCOMPLETE AND THE COMPLEMENTARY 
00005 C INCOMPLETE GAMMA FUNCTION DOUBLE PRECISION. IN BRIEF, THIS 
00006 C FUNCTION HELP TO COMPUTE THE DIFFERENT FRACTION PRESENT IN THE 
00007 C TOPMODEL FRAMEWORK. 
00008 C       MORE EXPLANATION ON THE SUBROUTINE USE ARE GIVEN IN THE NEXT 
00009 C COMMENTARY 
00010 C
00011 C THIS ROUTINE WAS FOUND ON http://www.netlib.org WEB SITE.
00012 C      
00013 C REFERENCE - W. GAUTSCHI, :: A COMPUTATIONAL PROCEDURE FOR INCOMPLETE 
00014 C GAMMA FUNCTIONS, ACM TRANS. MATH. SOFTWARE., (1979) 482-489
00015 C
00016 C######################################################################
00017 C######################################################################
00018 C
00019 C     LET  GAMMA(A)  DENOTE THE GAMMA FUNCTION AND  GAM(A,X)  THE
00020 C (COMPLEMENTARY) INCOMPLETE GAMMA FUNCTION,
00021 C
00022 C    GAM(A,X)=INTEGRAL FROM T=X TO T=INFINITY OF EXP(-T)*T**(A-1).
00023 C
00024 C LET  GAMSTAR(A,X)  DENOTE TRICOMI:S FORM OF THE INCOMPLETE GAMMA
00025 C FUNCTION, WHICH FOR A.GT.0. IS DEFINED BY
00026 C
00027 C  GAMSTAR(A,X)=(X**(-A)/GAMMA(A))*INTEGRAL FROM T=0 TO T=X OF
00028 C                EXP(-T)*T**(A-1),
00029 C
00030 C AND FOR A.LE.0. BY ANALYTIC CONTINUATION. FOR THE PURPOSE OF
00031 C THIS SUBROUTINE, THESE FUNCTIONS ARE NORMALIZED AS FOLLOWS&
00032 C
00033 C             GAM(A,X)/GAMMA(A),  IF A.GT.0.,
00034 C     G(A,X)=
00035 C             EXP(X)*X**(-A)*GAM(A,X),  IF A.LE.0.,
00036 C
00037 C     GSTAR(A,X)=(X**A)*GAMSTAR(A,X)
00038 C               =(1/GAMMA(A))*INTEGRAL FROM T=0 TO T=X OF EXP(-T)*T**(A-1)
00039 C
00040 C THE PROGRAM BELOW ATTEMPTS TO EVALUATE  G(A,X)  AND  GSTAR(A,X),
00041 C BOTH TO AN ACCURACY OF ACC SIGNIFICANT DECIMAL DIGITS, FOR ARBI-
00042 C TRARY REAL VALUES OF  A  AND NONNEGATIVE VALUES OF  X. THE SUB-
00043 C ROUTINE AUTOMATICALLY CHECKS FOR UNDERFLOW AND OVERFLOW CONDI-
00044 C TIONS AND RETURNS APPROPRIATE WARNINGS THROUGH THE OUTPUT PARA-
00045 C METERS  IFLG, IFLGST. A RESULT THAT UNDERFLOWS IS RETURNED WITH
00046 C THE VALUE 0., ONE THAT OVERFLOWS WITH THE VALUE OF THE LARGEST
00047 C MACHINE-REPRESENTABLE NUMBER.
00048 C
00049 C     NEAR LINES IN THE (A,X)-PLANE, A.LT.0., ALONG WHICH  GSTAR
00050 C VANISHES, THE ACCURACY SPECIFIED WILL BE ATTAINED ONLY IN TERMS
00051 C OF ABSOLUTE ERROR, NOT RELATIVE ERROR. THERE ARE OTHER (RARE)
00052 C INSTANCES IN WHICH THE ACCURACY ATTAINED IS SOMEWHAT LESS THAN
00053 C THE ACCURACY SPECIFIED. THE DISCREPANCY, HOWEVER, SHOULD NEVER
00054 C EXCEED ONE OR TWO (DECIMAL) ORDERS OF ACCURACY# NO INDICATION
00055 C OF THIS IS GIVEN THROUGH ERROR FLAGS.
00056 C
00057 C     PARAMETER LIST&
00058 C
00059 C        A - - THE FIRST ARGUMENT OF G AND GSTAR. TYPE REAL.
00060 C        X - - THE SECOND ARGUMENT OF G AND GSTAR. TYPE REAL.
00061 C      ACC - - THE NUMBER OF CORRECT SIGNIFICANT DECIMAL DIGITS
00062 C              DESIRED IN THE RESULTS. TYPE REAL.
00063 C        G - - AN OUTPUT VARIABLE RETURNING THE VALUE OF G(A,X).
00064 C              TYPE REAL.
00065 C    GSTAR - - AN OUTPUT VARIABLE RETURNING THE VALUE OF
00066 C              GSTAR(A,X). TYPE REAL.
00067 C     IFLG - - AN ERROR FLAG INDICATING A NUMBER OF ERROR CONDI-
00068 C              TIONS IN G UPON EXECUTION. TYPE INTEGER.
00069 C   IFLGST - - AN ERROR FLAG INDICATING A NUMBER OF ERROR CONDI-
00070 C              TIONS IN GSTAR UPON EXECUTION. TYPE INTEGER.
00071 C              THE VALUES OF IFLG AND IFLGST HAVE THE FOLLOWING
00072 C              MEANINGS&
00073 C              0 - NO ERROR CONDITION.
00074 C              1 - ILLEGAL NEGATIVE ARGUMENT X. THE ROUTINE EXITS
00075 C                  WITH THE VALUES ZERO FOR G AND GSTAR.
00076 C              2 - INFINITELY LARGE RESULT AT X=0. THE ROUTINE
00077 C                  RETURNS THE LARGEST MACHINE-REPRESENTABLE NUMBER.
00078 C              3 - (ONLY FOR IFLGST) GSTAR IS INDETERMINATE AT
00079 C                  A=0. AND X=0. THE ROUTINE RETURNS THE VALUE 1.,
00080 C                  WHICH IS THE LIMIT VALUE AS X TENDS TO ZERO FOR
00081 C                  FIXED A=0.
00082 C              4 - THE RESULT UNDERFLOWS. IT IS SET EQUAL TO 0.
00083 C              5 - THE RESULT OVERFLOWS. IT IS SET EQUAL TO THE
00084 C                  LARGEST MACHINE-REPRESENTABLE NUMBER, WITH
00085 C                  PROPER SIGN.
00086 C              6 - CONVERGENCE FAILS WITHIN 600 ITERATIONS, EITHER
00087 C                  IN TAYLOR:S SERIES OR IN LEGENDRE:S CONTINUED
00088 C                  FRACTION. REASON UNKNOWN. THE COMPUTATION IS
00089 C                  ABORTED AND THE ROUTINE RETURNS WITH ZERO
00090 C                  VALUES FOR G AND GSTAR.
00091 C
00092 C     ALL MACHINE-DEPENDENT PARAMETERS ARE COLLECTED IN THE FIRST
00093 C DATA DECLARATION. THEY ARE AS FOLLOWS&
00094 C
00095 C IN THE PROGRAM BELOW THESE PARAMETERS ARE SET TO CORRESPOND TO
00096 C THE MACHINE CHARACTERISTICS OF THE CDC 6500 COMPUTER.
00097 C
00098 C     THE SECOND DATA DECLARATION CONTAINS THE SINGLE PRECISION
00099 C VALUE OF ALOG(10.). THE NEXT DATA DECLARATION CONTAINS THE SUCCES-
00100 C SIVE COEFFICIENTS IN THE MACLAURIN EXPANSION OF (1/GAMMA(A+1))-1.
00101 C THEY ARE GIVEN TO AS MANY DECIMAL PLACES AS IS NECESSARY TO ACHIEVE
00102 C MACHINE PRECISION (ON THE CDC 6500 COMPUTER) IN THE RANGE
00103 C ABS(A).LE..5. MORE ACCURATE VALUES OF THESE COEFFICIENTS (TO
00104 C 31 DECIMAL PLACES) CAN BE FOUND IN TABLE 5 OF J.W.WRENCH,JR.,
00105 C CONCERNING TWO SERIES FOR THE GAMMA FUNCTION, MATH. COMPUT.
00106 C 22, 1968, 617-626.
00107 C
00108 C     THE SUBROUTINE CALLS ON A FUNCTION SUBROUTINE, NAMED  ALGA,
00109 C WHICH IS TO PROVIDE SINGLE PRECISION VALUES OF THE LOGARITHM OF
00110 C THE GAMMA FUNCTION FOR ARGUMENTS LARGER THAN OR EQUAL TO .5.
00111 C A POSSIBLE VERSION OF SUCH A FUNCTION SUBROUTINE IS APPENDED
00112 C TO THE PRESENT SUBROUTINE. IT IS TAYLORED TO THE ACCURACY RE-
00113 C QUIREMENTS OF THE CDC 6500 COMPUTER, AND USES RATIONAL APPROXI-
00114 C MATIONS DUE TO CODY AND HILLSTROM (MATH. COMPUT. 21, 1967, 198-
00115 C 203).
00116 C
00117 C     REFERENCE - W. GAUTSCHI, ::A COMPUTATIONAL PROCEDURE FOR
00118 C INCOMPLETE GAMMA FUNCTIONS, ACM TRANS. MATH. SOFTWARE.
00119 C
00120 !################################################################
00121 !
00122       SUBROUTINE DGAM(A, X, ACC, G, GSTAR, IFLG, IFLGST)
00123 !
00124 !################################################################
00125 
00126       USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
00127       USE PARKIND1  ,ONLY : JPRB
00128 
00129       IMPLICIT NONE
00130 
00131       REAL(KIND=JPRB) :: ZHOOK_HANDLE
00132 
00133       REAL A, X, ACC, G, GSTAR
00134       INTEGER IFLG, IFLGST, I, K, K1, MA
00135 
00136       DOUBLE PRECISION C, AL10, ALX, ALPHA, ALPREC,
00137      * TOP, BOT, AINF, EPS, EPS1, ES, SGA, AE, AA, AP1, AM1, AEP1,
00138      * AEM1, FMA, AEPS, SGAE, AAEPS, ALGP1, SGGA, ALGEP1, SGGS, ALGS,
00139      * ALG, SUMM, GA, Y, TERM, U, P, Q, R, V, T, H, SGT, A1X, RHO, XPA,
00140      * XMA, S, TEMP
00141       DIMENSION C(29)
00142       DATA AL10 /2.3025850929940456840179914547D0/
00143       DATA C /.57721566490153286060651209008D0,-.65587807152025388107701
00144      * 951515D0,-4.200263503409523552900393488D-2,.166538611382291489501
00145      * 7007951D0,-4.21977345555443367482083013D-2,-9.6219715278769735621
00146      * 149217D-3,7.2189432466630995423950103D-3,-1.165167591859065112113
00147      * 971D-3,-2.15241674114950972815730D-4,1.2805028238811618615320D-4,
00148      * -2.013485478078823865569D-5,-1.2504934821426706573D-6,
00149      * 1.1330272319816958824D-6,-2.056338416977607103D-7,
00150      * 6.1160951044814158D-9,5.0020076444692229D-9,-1.181274570487020D-9
00151      * ,1.04342671169110D-10,7.782263439905D-12,-3.696805618642D-12,
00152      * 5.1003702875D-13,-2.058326054D-14,-5.34812254D-15,1.2267786D-15,
00153      * -1.181259D-16,1.187D-18,1.412D-18,-2.30D-19,1.7D-20/
00154 
00155       IF (LHOOK) CALL DR_HOOK('DGAM',0,ZHOOK_HANDLE)
00156 
00157       G = 0.D0
00158       GSTAR = 0.D0
00159       IF (X.LT.0.D0) GO TO 290
00160 C
00161 C INITIALIZATION
00162 C
00163       IFLG = 0
00164       IFLGST = 0
00165       I = 0
00166       IF (X.GT.0.D0) ALX = DLOG(X*1.D0)
00167       ALPHA = X + .25D0
00168       IF (X.LT..25D0 .AND. X.GT.0.D0) ALPHA = DLOG(.5D0)/ALX
00169       ALPREC = DBLE(DIGITS(1.D0))*DLOG(2.D0)
00170       TOP    = DLOG10(HUGE(1.D0))
00171       BOT    = DLOG10(TINY(1.D0))
00172       AINF   = HUGE(1.D0)
00173       EPS = .5D0*10.D0**(-ACC)
00174       EPS1 = EPS/1.D2
00175       SGA = 1.D0
00176       IF (A.LT.0.D0) SGA = -SGA
00177       AE = A
00178       AA = DABS(A*1.D0)
00179       AP1 = A + 1.D0
00180       AEP1 = AP1
00181       MA = SNGL(.5D0-A)
00182       FMA = DBLE(FLOAT(MA))
00183       AEPS = A + FMA
00184       SGAE = 1.D0
00185       IF (AEPS.LT.0.D0) SGAE = -SGAE
00186       AAEPS = DABS(AEPS)
00187       ALGP1 = 0.D0
00188 C
00189 C EVALUATION OF THE LOGARITHM OF THE ABSOLUTE VALUE OF
00190 C GAMMA(A+1.) AND DETERMINATION OF THE SIGN OF GAMMA(A+1.)
00191 C
00192       SGGA = 1.D0
00193       IF (MA.LE.0) GO TO 10
00194       IF (AEPS.EQ.0.D0) GO TO 20
00195       SGGA = SGAE
00196       IF (MA.EQ.2*(MA/2)) SGGA = -SGGA
00197       ALGP1 = DLGA(AEPS+1.D0) - DLOG(AAEPS)
00198       IF (MA.EQ.1) GO TO 20
00199       ALGP1 = ALGP1 + DLGA(1.D0-AEPS) - DLGA(FMA-AEPS)
00200       GO TO 20
00201    10 ALGP1 = DLGA(AP1)
00202    20 ALGEP1 = ALGP1
00203       IF (X.GT.0.D0) GO TO 60
00204 C
00205 C EVALUATION OF GSTAR(A,0.) AND G(A,0.)
00206 C
00207       IF (A) 30, 40, 50
00208    30 IFLGST = 2
00209       GSTAR = AINF
00210       G = 1.D0/AA
00211       IF (LHOOK) CALL DR_HOOK('DGAM',1,ZHOOK_HANDLE)
00212       RETURN
00213    40 IFLGST = 3
00214       GSTAR = 1.D0
00215       IFLG = 2
00216       G = AINF
00217       IF (LHOOK) CALL DR_HOOK('DGAM',1,ZHOOK_HANDLE)
00218       RETURN
00219    50 G = 1.D0
00220       IF (LHOOK) CALL DR_HOOK('DGAM',1,ZHOOK_HANDLE)
00221       RETURN
00222    60 IF (A.GT.ALPHA) GO TO 220
00223       IF (X.GT.1.5D0) GO TO 240
00224       IF (A.LT.-.5D0) GO TO 170
00225 C
00226 C DIRECT EVALUATION OF G(A,X) AND GSTAR(A,X) FOR X.LE.1.5
00227 C AND -.5.LE.A.LE.ALPHA(X)
00228 C
00229       GSTAR = 1.D0
00230       IF (A.GE..5D0) GO TO 110
00231    70 SUMM = C(29)
00232       DO 80 K=1,28
00233         K1 = 29 - K
00234         SUMM = AE*SUMM + C(K1)
00235    80 CONTINUE
00236       GA = -SUMM/(1.D0+AE*SUMM)
00237       Y = AE*ALX
00238       IF (DABS(Y).GE.1.D0) GO TO 100
00239       SUMM = 1.D0
00240       TERM = 1.D0
00241       K = 1
00242    90 K = K + 1
00243       IF (K.GT.600) GO TO 330
00244       TERM = Y*TERM/DBLE(FLOAT(K))
00245       SUMM = SUMM + TERM
00246       IF (DABS(TERM).GT.EPS1*SUMM) GO TO 90
00247       U = GA - SUMM*ALX
00248       GO TO 120
00249   100 U = GA - (DEXP(Y)-1.D0)/AE
00250       GO TO 120
00251   110 TEMP=DLGA(A*1.D0)
00252       U = DEXP(TEMP) - (X**A)/A
00253   120 P = AE*X
00254       Q = AEP1
00255       R = AE + 3.D0
00256       TERM = 1.D0
00257       SUMM = 1.D0
00258       K = 1
00259   130 K = K + 1
00260       IF (K.GT.600) GO TO 330
00261       P = P + X
00262       Q = Q + R
00263       R = R + 2.D0
00264       TERM = -P*TERM/Q
00265       SUMM = SUMM + TERM
00266       IF (DABS(TERM).GT.EPS1*SUMM) GO TO 130
00267       V = (X**AEP1)*SUMM/AEP1
00268       G = U + V
00269       IF (I.EQ.1) GO TO 180
00270       IF (A) 140, 150, 160
00271   140 T = DEXP(X*1.D0)*X**(-A)
00272       G = T*G
00273       GSTAR = 1.D0 - A*G*DEXP(-ALGP1)/T
00274       IF (LHOOK) CALL DR_HOOK('DGAM',1,ZHOOK_HANDLE)
00275       RETURN
00276   150 G = DEXP(X*1.D0)*G
00277       IF (LHOOK) CALL DR_HOOK('DGAM',1,ZHOOK_HANDLE)
00278       RETURN
00279   160 G = A*G*DEXP(-ALGP1)
00280       GSTAR = 1.D0 - G
00281       IF (LHOOK) CALL DR_HOOK('DGAM',1,ZHOOK_HANDLE)
00282       RETURN
00283 C
00284 C RECURSIVE EVALUATION OF G(A,X) FOR X.LE.1.5 AND A.LT.-.5
00285 C
00286   170 I = 1
00287       AE = AEPS
00288       AEP1 = AEPS + 1.D0
00289       IF (X.LT..25D0 .AND. AE.GT.ALPHA) GO TO 210
00290       GO TO 70
00291   180 G = G*DEXP(X*1.D0)*X**(-AE)
00292       DO 190 K=1,MA
00293         G = (1.D0-X*G)/(DBLE(FLOAT(K))-AE)
00294   190 CONTINUE
00295       ALG = DLOG(G*1.D0)
00296 C
00297 C EVALUATION OF GSTAR(A,X) IN TERMS OF G(A,X)
00298 C
00299   200 GSTAR = 1.D0
00300       IF (MA.GE.0 .AND. AEPS.EQ.0.D0 .AND. LHOOK) 
00301      *      CALL DR_HOOK('DGAM',1,ZHOOK_HANDLE)
00302       IF (MA.GE.0 .AND. AEPS.EQ.0.D0) RETURN
00303       SGT = SGA*SGGA
00304       T = DLOG(AA) - X + A*ALX + ALG - ALGP1
00305       IF (T.LT.-ALPREC .AND. LHOOK) CALL DR_HOOK('DGAM',1,ZHOOK_HANDLE)
00306       IF (T.LT.-ALPREC) RETURN
00307       IF (T.GE.TOP) GO TO 320
00308       GSTAR = 1.D0 - SGT*DEXP(T)
00309       IF (LHOOK) CALL DR_HOOK('DGAM',1,ZHOOK_HANDLE)
00310       RETURN
00311   210 I = 2
00312       ALGEP1 = DLGA(AEP1)
00313 C
00314 C EVALUATION OF GSTAR(A,X) FOR A.GT.ALPHA(X) BY TAYLOR
00315 C EXPANSION
00316 C
00317   220 G = 1.D0
00318       TERM = 1.D0
00319       SUMM = 1.D0
00320       K = 0
00321   230 K = K + 1
00322       IF (K.GT.600) GO TO 340
00323       TERM = X*TERM/(AE+DBLE(FLOAT(K)))
00324       SUMM = SUMM + TERM
00325       IF (DABS(TERM).GT.EPS*SUMM) GO TO 230
00326       ALGS = AE*ALX - X + DLOG(SUMM) - ALGEP1
00327       IF (ALGS.LE.BOT) GO TO 310
00328       GSTAR = DEXP(ALGS)
00329       G = 1.D0 - GSTAR
00330       IF (I.NE.2 .AND. LHOOK) CALL DR_HOOK('DGAM',1,ZHOOK_HANDLE)
00331       IF (I.NE.2) RETURN
00332       G = G*DEXP(ALGEP1)/AE
00333       GO TO 180
00334 C
00335 C EVALUATION OF G(A,X) FOR X.GT.1.5 AND A.LE.ALPHA(X) BY
00336 C MEANS OF THE LEGENDRE CONTINUED FRACTION
00337 C
00338   240 GSTAR = 1.D0
00339       XPA = X + 1.D0 - A
00340       XMA = X - 1.D0 - A
00341       P = 0.D0
00342       Q = XPA*XMA
00343       R = 4.D0*XPA
00344       S = -A + 1.D0
00345       TERM = 1.D0
00346       SUMM = 1.D0
00347       RHO = 0.D0
00348       K = 1
00349   250 K = K + 1
00350       IF (K.GT.600) GO TO 330
00351       P = P + S
00352       Q = Q + R
00353       R = R + 8.D0
00354       S = S + 2.D0
00355       T = P*(1.D0+RHO)
00356       RHO = T/(Q-T)
00357       TERM = RHO*TERM
00358       SUMM = SUMM + TERM
00359       IF (DABS(TERM).GT.EPS*SUMM) GO TO 250
00360       IF (A) 260, 270, 280
00361   260 G = SUMM/XPA
00362       ALG = DLOG(G*1.D0)
00363       GO TO 200
00364   270 G = SUMM/XPA
00365       IF (LHOOK) CALL DR_HOOK('DGAM',1,ZHOOK_HANDLE)
00366       RETURN
00367   280 ALG = A*ALX - X + DLOG(A*SUMM/XPA) - ALGP1
00368       IF (ALG.LE.BOT) GO TO 300
00369       G = DEXP(ALG)
00370       GSTAR = 1.D0 - G
00371       IF (LHOOK) CALL DR_HOOK('DGAM',1,ZHOOK_HANDLE)
00372       RETURN
00373   290 IFLG = 1
00374       IFLGST = 1
00375       IF (LHOOK) CALL DR_HOOK('DGAM',1,ZHOOK_HANDLE)
00376       RETURN
00377   300 IFLG = 4
00378       IF (LHOOK) CALL DR_HOOK('DGAM',1,ZHOOK_HANDLE)
00379       RETURN
00380   310 IFLGST = 4
00381       IF (LHOOK) CALL DR_HOOK('DGAM',1,ZHOOK_HANDLE)
00382       RETURN
00383   320 IFLGST = 5
00384       GSTAR = -SGT*AINF
00385       IF (LHOOK) CALL DR_HOOK('DGAM',1,ZHOOK_HANDLE)      
00386       RETURN
00387   330 IFLG = 6
00388       IF (LHOOK) CALL DR_HOOK('DGAM',1,ZHOOK_HANDLE)
00389       RETURN
00390   340 IFLGST = 6
00391       IF (LHOOK) CALL DR_HOOK('DGAM',1,ZHOOK_HANDLE)
00392       RETURN
00393      
00394       CONTAINS 
00395 
00396       FUNCTION DLGA(DX) RESULT(DLGAR)
00397 
00398       USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
00399       USE PARKIND1  ,ONLY : JPRB
00400 
00401       REAL(KIND=JPRB) :: ZHOOK_HANDLE
00402 
00403       DOUBLE PRECISION DLGAR
00404       DOUBLE PRECISION DBNUM, DBDEN, DX, DC, DP, DY, DT, DS
00405       DIMENSION DBNUM(8), DBDEN(8)
00406       DATA DBNUM /-3.617D3,1.D0,-6.91D2,1.D0,-1.D0,1.D0,-1.D0,1.D0/,
00407      * DBDEN /1.224D5,1.56D2,3.6036D5,1.188D3,1.68D3,1.26D3,3.6D2,1.2D1/
00408 
00409       IF (LHOOK) CALL DR_HOOK('DLGA',0,ZHOOK_HANDLE)
00410 
00411       DC = .5D0*DLOG(8.D0*DATAN(1.D0))
00412       DP = 1.D0
00413       DY = DX
00414       Y = SNGL(DY)
00415 C
00416 C THE CONDITIONAL CLAUSE IN THE NEXT STATEMENT EXPRESSES THE
00417 C INEQUALITY Y.GT.EXP(.121189*DPREC+.053905), WHERE DPREC IS THE
00418 C NUMBER OF DECIMAL DIGITS CARRIED IN DOUBLE PRECISION FLOATING-POINT
00419 C ARITHMETIC.
00420 C
00421    10 IF (Y.GT.35.027) GO TO 20
00422       DP = DY*DP
00423       DY = DY + 1.D0
00424       Y = SNGL(DY)
00425       GO TO 10
00426    20 DT = 1.D0/(DY*DY)
00427       DS = 4.3867D4/2.44188D5
00428       DO 30 I=1,8
00429         DS = DT*DS + DBNUM(I)/DBDEN(I)
00430    30 CONTINUE
00431       DLGAR = (DY-.5D0)*DLOG(DY) - DY + DC + DS/DY - DLOG(DP)
00432       IF (LHOOK) CALL DR_HOOK('DLGA',1,ZHOOK_HANDLE)
00433       RETURN
00434       END FUNCTION DLGA
00435 
00436       END SUBROUTINE DGAM