SURFEX  V8_0
Surfex V8_0 release
 All Classes Files Functions Variables
dgam.F90
Go to the documentation of this file.
1 !SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2 !SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
3 !SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
4 !SFX_LIC for details. version 1.
5 !######################################################################
6 !######################################################################
7 !
8 ! THIS ROUTINE COMPUTE THE INCOMPLETE AND THE COMPLEMENTARY
9 ! INCOMPLETE GAMMA FUNCTION REAL*8 IN BRIEF, THIS
10 ! FUNCTION HELP TO COMPUTE THE DIFFERENT FRACTION PRESENT IN THE
11 ! TOPMODEL FRAMEWORK.
12 ! MORE EXPLANATION ON THE SUBROUTINE USE ARE GIVEN IN THE NEXT
13 ! COMMENTARY
14 !
15 ! THIS ROUTINE WAS FOUND ON http://www.netlib.org WEB SITE.
16 !
17 ! REFERENCE - W. GAUTSCHI, :: A COMPUTATIONAL PROCEDURE FOR INCOMPLETE
18 ! GAMMA FUNCTIONS, ACM TRANS. MATH. SOFTWARE., (1979) 482-489
19 !
20 !######################################################################
21 !######################################################################
22 !
23 ! LET GAMMA(A) DENOTE THE GAMMA FUNCTION AND GAM(A,X) THE
24 ! (COMPLEMENTARY) INCOMPLETE GAMMA FUNCTION,
25 !
26 ! GAM(A,X)=INTEGRAL FROM T=X TO T=INFINITY OF EXP(-T)*T**(A-1).
27 !
28 ! LET GAMSTAR(A,X) DENOTE TRICOMI:S FORM OF THE INCOMPLETE GAMMA
29 ! FUNCTION, WHICH FOR A.GT.0. IS DEFINED BY
30 !
31 ! GAMSTAR(A,X)=(X**(-A)/GAMMA(A))*INTEGRAL FROM T=0 TO T=X OF
32 ! EXP(-T)*T**(A-1),
33 !
34 ! AND FOR A.LE.0. BY ANALYTIC CONTINUATION. FOR THE PURPOSE OF
35 ! THIS SUBROUTINE, THESE FUNCTIONS ARE NORMALIZED AS FOLLOWS&
36 !
37 ! GAM(A,X)/GAMMA(A), IF A.GT.0.,
38 ! G(A,X)=
39 ! EXP(X)*X**(-A)*GAM(A,X), IF A.LE.0.,
40 !
41 ! GSTAR(A,X)=(X**A)*GAMSTAR(A,X)
42 ! =(1/GAMMA(A))*INTEGRAL FROM T=0 TO T=X OF EXP(-T)*T**(A-1)
43 !
44 ! THE PROGRAM BELOW ATTEMPTS TO EVALUATE G(A,X) AND GSTAR(A,X),
45 ! BOTH TO AN ACCURACY OF ACC SIGNIFICANT DECIMAL DIGITS, FOR ARBI-
46 ! TRARY REAL VALUES OF A AND NONNEGATIVE VALUES OF X. THE SUB-
47 ! ROUTINE AUTOMATICALLY CHECKS FOR UNDERFLOW AND OVERFLOW CONDI-
48 ! TIONS AND RETURNS APPROPRIATE WARNINGS THROUGH THE OUTPUT PARA-
49 ! METERS IFLG, IFLGST. A RESULT THAT UNDERFLOWS IS RETURNED WITH
50 ! THE VALUE 0., ONE THAT OVERFLOWS WITH THE VALUE OF THE LARGEST
51 ! MACHINE-REPRESENTABLE NUMBER.
52 !
53 ! NEAR LINES IN THE (A,X)-PLANE, A.LT.0., ALONG WHICH GSTAR
54 ! VANISHES, THE ACCURACY SPECIFIED WILL BE ATTAINED ONLY IN TERMS
55 ! OF ABSOLUTE ERROR, NOT RELATIVE ERROR. THERE ARE OTHER (RARE)
56 ! INSTANCES IN WHICH THE ACCURACY ATTAINED IS SOMEWHAT LESS THAN
57 ! THE ACCURACY SPECIFIED. THE DISCREPANCY, HOWEVER, SHOULD NEVER
58 ! EXCEED ONE OR TWO (DECIMAL) ORDERS OF ACCURACY# NO INDICATION
59 ! OF THIS IS GIVEN THROUGH ERROR FLAGS.
60 !
61 ! PARAMETER LIST&
62 !
63 ! A - - THE FIRST ARGUMENT OF G AND GSTAR. TYPE REAL.
64 ! X - - THE SECOND ARGUMENT OF G AND GSTAR. TYPE REAL.
65 ! ACC - - THE NUMBER OF CORRECT SIGNIFICANT DECIMAL DIGITS
66 ! DESIRED IN THE RESULTS. TYPE REAL.
67 ! G - - AN OUTPUT VARIABLE RETURNING THE VALUE OF G(A,X).
68 ! TYPE REAL.
69 ! GSTAR - - AN OUTPUT VARIABLE RETURNING THE VALUE OF
70 ! GSTAR(A,X). TYPE REAL.
71 ! IFLG - - AN ERROR FLAG INDICATING A NUMBER OF ERROR CONDI-
72 ! TIONS IN G UPON EXECUTION. TYPE INTEGER.
73 ! IFLGST - - AN ERROR FLAG INDICATING A NUMBER OF ERROR CONDI-
74 ! TIONS IN GSTAR UPON EXECUTION. TYPE INTEGER.
75 ! THE VALUES OF IFLG AND IFLGST HAVE THE FOLLOWING
76 ! MEANINGS&
77 ! 0 - NO ERROR CONDITION.
78 ! 1 - ILLEGAL NEGATIVE ARGUMENT X. THE ROUTINE EXITS
79 ! WITH THE VALUES ZERO FOR G AND GSTAR.
80 ! 2 - INFINITELY LARGE RESULT AT X=0. THE ROUTINE
81 ! RETURNS THE LARGEST MACHINE-REPRESENTABLE NUMBER.
82 ! 3 - (ONLY FOR IFLGST) GSTAR IS INDETERMINATE AT
83 ! A=0. AND X=0. THE ROUTINE RETURNS THE VALUE 1.,
84 ! WHICH IS THE LIMIT VALUE AS X TENDS TO ZERO FOR
85 ! FIXED A=0.
86 ! 4 - THE RESULT UNDERFLOWS. IT IS SET EQUAL TO 0.
87 ! 5 - THE RESULT OVERFLOWS. IT IS SET EQUAL TO THE
88 ! LARGEST MACHINE-REPRESENTABLE NUMBER, WITH
89 ! PROPER SIGN.
90 ! 6 - CONVERGENCE FAILS WITHIN 600 ITERATIONS, EITHER
91 ! IN TAYLOR:S SERIES OR IN LEGENDRE:S CONTINUED
92 ! FRACTION. REASON UNKNOWN. THE COMPUTATION IS
93 ! ABORTED AND THE ROUTINE RETURNS WITH ZERO
94 ! VALUES FOR G AND GSTAR.
95 !
96 ! ALL MACHINE-DEPENDENT PARAMETERS ARE COLLECTED IN THE FIRST
97 ! DATA DECLARATION. THEY ARE AS FOLLOWS&
98 !
99 ! IN THE PROGRAM BELOW THESE PARAMETERS ARE SET TO CORRESPOND TO
100 ! THE MACHINE CHARACTERISTICS OF THE CDC 6500 COMPUTER.
101 !
102 ! THE SECOND DATA DECLARATION CONTAINS THE SINGLE PRECISION
103 ! VALUE OF ALOG(10.). THE NEXT DATA DECLARATION CONTAINS THE SUCCES-
104 ! SIVE COEFFICIENTS IN THE MACLAURIN EXPANSION OF (1/GAMMA(A+1))-1.
105 ! THEY ARE GIVEN TO AS MANY DECIMAL PLACES AS IS NECESSARY TO ACHIEVE
106 ! MACHINE PRECISION (ON THE CDC 6500 COMPUTER) IN THE RANGE
107 ! ABS(A).LE..5. MORE ACCURATE VALUES OF THESE COEFFICIENTS (TO
108 ! 31 DECIMAL PLACES) CAN BE FOUND IN TABLE 5 OF J.W.WRENCH,JR.,
109 ! CONCERNING TWO SERIES FOR THE GAMMA FUNCTION, MATH. COMPUT.
110 ! 22, 1968, 617-626.
111 !
112 ! THE SUBROUTINE CALLS ON A FUNCTION SUBROUTINE, NAMED ALGA,
113 ! WHICH IS TO PROVIDE SINGLE PRECISION VALUES OF THE LOGARITHM OF
114 ! THE GAMMA FUNCTION FOR ARGUMENTS LARGER THAN OR EQUAL TO .5.
115 ! A POSSIBLE VERSION OF SUCH A FUNCTION SUBROUTINE IS APPENDED
116 ! TO THE PRESENT SUBROUTINE. IT IS TAYLORED TO THE ACCURACY RE-
117 ! QUIREMENTS OF THE CDC 6500 COMPUTER, AND USES RATIONAL APPROXI-
118 ! MATIONS DUE TO CODY AND HILLSTROM (MATH. COMPUT. 21, 1967, 198-
119 ! 203).
120 !
121 ! REFERENCE - W. GAUTSCHI, ::A COMPUTATIONAL PROCEDURE FOR
122 ! INCOMPLETE GAMMA FUNCTIONS, ACM TRANS. MATH. SOFTWARE.
123 !
124 !----------------------------------------------------
125 !! MODIFICATIONS
126 !! -------------
127 !! J.Escobar 10/06/2013: replace DOUBLE PRECISION by REAL to handle problem for promotion of real on IBM SP
128 !----------------------------------------------------
129 !################################################################
130 !
131  SUBROUTINE dgam(A, X, ACC, G, GSTAR, IFLG, IFLGST)
132 !
133 !################################################################
134 
135  USE yomhook ,ONLY : lhook, dr_hook
136  USE parkind1 ,ONLY : jprb
137 
138  IMPLICIT NONE
139 
140  REAL(KIND=JPRB) :: zhook_handle
141 
142  INTEGER,PARAMETER :: jdbl=8
143  REAL :: a, x, acc, g, gstar
144  INTEGER :: iflg, iflgst, i, k, k1, ma
145 
146  REAL(KIND=JDBL) :: alx, alpha, alprec
147  REAL(KIND=JDBL) :: top, bot, ainf, eps, eps1, es, sga
148  REAL(KIND=JDBL) :: ae, aa, ap1, am1, aep1
149  REAL(KIND=JDBL) :: aem1, fma, aeps, sgae, aaeps, algp1
150  REAL(KIND=JDBL) :: sgga, algep1, sggs, algs
151  REAL(KIND=JDBL) :: alg, summ, ga, y, term, u, p, q, r, v, t, h
152  REAL(KIND=JDBL) :: sgt, a1x, rho, xpa
153  REAL(KIND=JDBL) :: xma, s, temp
154 
155 !RJ REAL(KIND=JDBL),SAVE::AL10=2.3025850929940456840179914547e0_JDBL
156  REAL(KIND=JDBL),DIMENSION(29),PARAMETER :: c=(/ &
157  & .57721566490153286060651209008e0_JDBL, &
158  & -.65587807152025388107701951515e0_JDBL, &
159  & -4.200263503409523552900393488e-2_JDBL, &
160  & .1665386113822914895017007951e0_JDBL, &
161  & -4.21977345555443367482083013e-2_JDBL, &
162  & -9.6219715278769735621149217e-3_JDBL, &
163  & 7.2189432466630995423950103e-3_JDBL, &
164  & -1.165167591859065112113971e-3_JDBL, &
165  & -2.15241674114950972815730e-4_JDBL, &
166  & 1.2805028238811618615320e-4_JDBL, &
167  & -2.013485478078823865569e-5_JDBL, &
168  & -1.2504934821426706573e-6_JDBL, &
169  & 1.1330272319816958824e-6_JDBL, &
170  & -2.056338416977607103e-7_JDBL, &
171  & 6.1160951044814158e-9_JDBL, &
172  & 5.0020076444692229e-9_JDBL, &
173  & -1.181274570487020e-9_JDBL, &
174  & 1.04342671169110e-10_JDBL, &
175  & 7.782263439905e-12_JDBL, &
176  & -3.696805618642e-12_JDBL, &
177  & 5.1003702875e-13_JDBL, &
178  & -2.058326054e-14_JDBL, &
179  & -5.34812254e-15_JDBL, &
180  & 1.2267786e-15_JDBL, &
181  & -1.181259e-16_JDBL, &
182  & 1.187e-18_JDBL, &
183  & 1.412e-18_JDBL, &
184  & -2.30e-19_JDBL, &
185  & 1.7e-20_JDBL/)
186 
187  IF (lhook) CALL dr_hook('DGAM',0,zhook_handle)
188 
189  g = 0.e0_jdbl
190  gstar = 0.e0_jdbl
191  IF (x<0.e0_jdbl) go to 290
192 !
193 ! INITIALIZATION
194 !
195  iflg = 0
196  iflgst = 0
197  i = 0
198  IF (x>0.e0_jdbl) alx = log(x*1.e0_jdbl)
199  alpha = x + .25e0_jdbl
200  IF (x<0.25e0_jdbl .AND. x>0.e0_jdbl) alpha = log(0.5e0_jdbl)/alx
201  alprec = digits(1.e0_jdbl)*log(2.e0_jdbl)
202  top = log10(huge(1.e0_jdbl))
203  bot = log10(tiny(1.e0_jdbl))
204  ainf = huge(1.e0_jdbl)
205  eps = 0.5e0_jdbl*10.e0_jdbl**(-acc)
206  eps1 = eps/1.e2_jdbl
207  sga = 1.e0_jdbl
208  IF (a<0.e0_jdbl) sga = -sga
209  ae = a
210  aa = abs(a*1.e0_jdbl)
211  ap1 = a + 1.e0_jdbl
212  aep1 = ap1
213  ma = int(0.5e0_jdbl-a)
214  fma = ma
215  aeps = a + fma
216  sgae = 1.e0_jdbl
217  IF (aeps<0.e0_jdbl) sgae = -sgae
218  aaeps = abs(aeps)
219  algp1 = 0.e0_jdbl
220 !
221 ! EVALUATION OF THE LOGARITHM OF THE ABSOLUTE VALUE OF
222 ! GAMMA(A+1.) AND DETERMINATION OF THE SIGN OF GAMMA(A+1.)
223 !
224  sgga = 1.e0_jdbl
225  IF (ma<=0) go to 10
226  IF (aeps==0.e0_jdbl) go to 20
227  sgga = sgae
228  IF (ma==2*(ma/2)) sgga = -sgga
229  algp1 = dlga(aeps+1.e0_jdbl) - log(aaeps)
230  IF (ma==1) go to 20
231  algp1 = algp1 + dlga(1.e0_jdbl-aeps) - dlga(fma-aeps)
232  go to 20
233  10 algp1 = dlga(ap1)
234  20 algep1 = algp1
235  IF (x>0.e0_jdbl) go to 60
236 !
237 ! EVALUATION OF GSTAR(A,0.) AND G(A,0.)
238 !
239  IF (a) 30, 40, 50
240  30 iflgst = 2
241  gstar = ainf
242  g = 1.e0_jdbl/aa
243  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
244  RETURN
245  40 iflgst = 3
246  gstar = 1.e0_jdbl
247  iflg = 2
248  g = ainf
249  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
250  RETURN
251  50 g = 1.e0_jdbl
252  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
253  RETURN
254  60 IF (a>alpha) go to 220
255  IF (x>1.5e0_jdbl) go to 240
256  IF (a<-.5e0_jdbl) go to 170
257 !
258 ! DIRECT EVALUATION OF G(A,X) AND GSTAR(A,X) FOR X.LE.1.5
259 ! AND -.5.LE.A.LE.ALPHA(X)
260 !
261  gstar = 1.e0_jdbl
262  IF (a>=0.5e0_jdbl) go to 110
263  70 summ = c(29)
264  DO 80 k=1,28
265  k1 = 29 - k
266  summ = ae*summ + c(k1)
267  80 CONTINUE
268  ga = -summ/(1.e0_jdbl+ae*summ)
269  y = ae*alx
270  IF (abs(y)>=1.e0_jdbl) go to 100
271  summ = 1.e0_jdbl
272  term = 1.e0_jdbl
273  k = 1
274  90 k = k + 1
275  IF (k>600) go to 330
276  term = y*term/k
277  summ = summ + term
278  IF (abs(term)>eps1*summ) go to 90
279  u = ga - summ*alx
280  go to 120
281  100 u = ga - (exp(y)-1.e0_jdbl)/ae
282  go to 120
283  110 temp=dlga(a*1.e0_jdbl)
284  u = dexp(temp) - (x**a)/a
285  120 p = ae*x
286  q = aep1
287  r = ae + 3.e0_jdbl
288  term = 1.e0_jdbl
289  summ = 1.e0_jdbl
290  k = 1
291  130 k = k + 1
292  IF (k>600) go to 330
293  p = p + x
294  q = q + r
295  r = r + 2.e0_jdbl
296  term = -p*term/q
297  summ = summ + term
298  IF (abs(term)>eps1*summ) go to 130
299  v = (x**aep1)*summ/aep1
300  g = u + v
301  IF (i==1) go to 180
302  IF (a) 140, 150, 160
303  140 t = exp(x*1.e0_jdbl)*x**(-a)
304  g = t*g
305  gstar = 1.e0_jdbl - a*g*dexp(-algp1)/t
306  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
307  RETURN
308  150 g = exp(x*1.e0_jdbl)*g
309  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
310  RETURN
311  160 g = a*g*exp(-algp1)
312  gstar = 1.e0_jdbl - g
313  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
314  RETURN
315 !
316 ! RECURSIVE EVALUATION OF G(A,X) FOR X.LE.1.5 AND A.LT.-.5
317 !
318  170 i = 1
319  ae = aeps
320  aep1 = aeps + 1.e0_jdbl
321  IF (x<0.25e0_jdbl .AND. ae>alpha) go to 210
322  go to 70
323  180 g = g*exp(x*1.e0_jdbl)*x**(-ae)
324  DO 190 k=1,ma
325  g = (1.e0_jdbl-x*g)/(k-ae)
326  190 CONTINUE
327  alg = log(g*1.e0_jdbl)
328 !
329 ! EVALUATION OF GSTAR(A,X) IN TERMS OF G(A,X)
330 !
331  200 gstar = 1.e0_jdbl
332  IF (ma>=0 .AND. aeps==0.e0_jdbl .AND. lhook) &
333  & CALL dr_hook('DGAM',1,zhook_handle)
334  IF (ma>=0 .AND. aeps==0.0_jdbl) RETURN
335  sgt = sga*sgga
336  t = log(aa) - x + a*alx + alg - algp1
337  IF (t<-alprec .AND. lhook) CALL dr_hook('DGAM',1,zhook_handle)
338  IF (t<-alprec) RETURN
339  IF (t>=top) go to 320
340  gstar = 1.e0_jdbl - sgt*exp(t)
341  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
342  RETURN
343  210 i = 2
344  algep1 = dlga(aep1)
345 !
346 ! EVALUATION OF GSTAR(A,X) FOR A.GT.ALPHA(X) BY TAYLOR
347 ! EXPANSION
348 !
349  220 g = 1.e0_jdbl
350  term = 1.e0_jdbl
351  summ = 1.e0_jdbl
352  k = 0
353  230 k = k + 1
354  IF (k>600) go to 340
355  term = x*term/(ae+k)
356  summ = summ + term
357  IF (abs(term)>eps*summ) go to 230
358  algs = ae*alx - x + log(summ) - algep1
359  IF (algs<=bot) go to 310
360  gstar = exp(algs)
361  g = 1.e0_jdbl - gstar
362  IF (i/=2 .AND. lhook) CALL dr_hook('DGAM',1,zhook_handle)
363  IF (i/=2) RETURN
364  g = g*exp(algep1)/ae
365  go to 180
366 !
367 ! EVALUATION OF G(A,X) FOR X.GT.1.5 AND A.LE.ALPHA(X) BY
368 ! MEANS OF THE LEGENDRE CONTINUED FRACTION
369 !
370  240 gstar = 1.e0_jdbl
371  xpa = x + 1.e0_jdbl - a
372  xma = x - 1.e0_jdbl - a
373  p = 0.e0_jdbl
374  q = xpa*xma
375  r = 4.e0_jdbl*xpa
376  s = -a + 1.e0_jdbl
377  term = 1.e0_jdbl
378  summ = 1.e0_jdbl
379  rho = 0.e0_jdbl
380  k = 1
381  250 k = k + 1
382  IF (k>600) go to 330
383  p = p + s
384  q = q + r
385  r = r + 8.e0_jdbl
386  s = s + 2.e0_jdbl
387  t = p*(1.e0_jdbl+rho)
388  rho = t/(q-t)
389  term = rho*term
390  summ = summ + term
391  IF (abs(term)>eps*summ) go to 250
392  IF (a) 260, 270, 280
393  260 g = summ/xpa
394  alg = log(g*1.e0_jdbl)
395  go to 200
396  270 g = summ/xpa
397  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
398  RETURN
399  280 alg = a*alx - x + log(a*summ/xpa) - algp1
400  IF (alg<=bot) go to 300
401  g = exp(alg)
402  gstar = 1.e0_jdbl - g
403  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
404  RETURN
405  290 iflg = 1
406  iflgst = 1
407  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
408  RETURN
409  300 iflg = 4
410  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
411  RETURN
412  310 iflgst = 4
413  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
414  RETURN
415  320 iflgst = 5
416  gstar = -sgt*ainf
417  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
418  RETURN
419  330 iflg = 6
420  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
421  RETURN
422  340 iflgst = 6
423  IF (lhook) CALL dr_hook('DGAM',1,zhook_handle)
424  RETURN
425 
426  CONTAINS
427 
428  FUNCTION dlga(DX) RESULT(DLGAR)
429 
430  USE yomhook ,ONLY : lhook, dr_hook
431  USE parkind1 ,ONLY : jprb
432 
433  IMPLICIT NONE
434 
435  REAL(KIND=JPRB) :: zhook_handle
436 
437  REAL(KIND=JDBL) :: dlgar
438  REAL(KIND=JDBL) :: dx, dc, dp, dy, dt, ds
439  REAL(KIND=JDBL),DIMENSION(8),PARAMETER :: dbnum=(/ &
440  & -3.617e3_JDBL,1.e0_JDBL,-6.91e2_JDBL,1.e0_JDBL, &
441  & -1.e0_JDBL, 1.e0_JDBL,-1.e0_JDBL, 1.e0_JDBL/)
442  REAL(KIND=JDBL),DIMENSION(8),PARAMETER :: dbden=(/ &
443  & 1.224e5_JDBL,1.56e2_JDBL,3.6036e5_JDBL,1.188e3_JDBL, &
444  & 1.68e3_JDBL, 1.26e3_JDBL,3.6e2_JDBL, 1.2e1_JDBL/)
445 
446  IF (lhook) CALL dr_hook('DLGA',0,zhook_handle)
447 
448  dc = 0.5e0_jdbl*log(8.e0_jdbl*atan(1.e0_jdbl))
449  dp = 1.e0_jdbl
450  dy = dx
451  y = dy
452 !
453 ! THE CONDITIONAL CLAUSE IN THE NEXT STATEMENT EXPRESSES THE
454 ! INEQUALITY Y.GT.EXP(.121189*DPREC+.053905), WHERE DPREC IS THE
455 ! NUMBER OF DECIMAL DIGITS CARRIED IN REAL*8 FLOATING-POINT
456 ! ARITHMETIC.
457 !
458 !RJ-remark this magic number is for 53 significant bit prec only +/-!!!!
459  10 IF (y>35.027_jdbl) go to 20
460  dp = dy*dp
461  dy = dy + 1.e0_jdbl
462  y = dy
463  go to 10
464  20 dt = 1.e0_jdbl/(dy*dy)
465  ds = 4.3867e4_jdbl/2.44188e5_jdbl
466  DO 30 i=1,8
467  ds = dt*ds + dbnum(i)/dbden(i)
468  30 CONTINUE
469  dlgar = (dy-0.5e0_jdbl)*log(dy) - dy + dc + ds/dy - log(dp)
470  IF (lhook) CALL dr_hook('DLGA',1,zhook_handle)
471  RETURN
472  END FUNCTION dlga
473 
474  END SUBROUTINE dgam
dlga
real(kind=jdbl) function dlga(DX)
Definition: dgam.F90:428
dgam
subroutine dgam(A, X, ACC, G, GSTAR, IFLG, IFLGST)
Definition: dgam.F90:131
Generated on Fri May 13 2016 09:31:03 for SURFEX by   doxygen 1.8.5