source: lmdz_wrf/WRFV3/phys/module_cam_error_function.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 23.9 KB
Line 
1!------------------------------------------------------------------------
2! Based on error_function.F90 from CAM
3! Ported to WRF by William.Gustafson@pnl.gov, Dec. 2009
4!------------------------------------------------------------------------
5
6module error_function
7
8! This module provides generic interfaces for functions that evaluate
9! erf(x), erfc(x), and exp(x*x)*erfc(x) in either single or double precision.
10
11implicit none
12private
13save
14
15! Public functions
16public :: erf, erfc, erfcx
17
18interface erf
19   module procedure erf_r4
20   module procedure derf
21end interface
22
23interface erfc
24   module procedure erfc_r4
25   module procedure derfc
26end interface
27
28interface erfcx
29   module procedure erfcx_r4
30   module procedure derfcx
31end interface
32
33! Private variables
34integer, parameter :: r4 = selected_real_kind(6)  ! 4 byte real
35integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real
36
37contains
38
39!------------------------------------------------------------------
40!
41! 6 December 2006 -- B. Eaton
42! The following comments are from the original version of CALERF.
43! The only changes in implementing this module are that the function
44! names previously used for the single precision versions have been
45! adopted for the new generic interfaces.  To support these interfaces
46! there is now both a single precision version (calerf_r4) and a
47! double precision version (calerf_r8) of CALERF below.  These versions
48! are hardcoded to use IEEE arithmetic.
49!
50!------------------------------------------------------------------
51!
52! This packet evaluates  erf(x),  erfc(x),  and  exp(x*x)*erfc(x)
53!   for a real argument  x.  It contains three FUNCTION type
54!   subprograms: ERF, ERFC, and ERFCX (or DERF, DERFC, and DERFCX),
55!   and one SUBROUTINE type subprogram, CALERF.  The calling
56!   statements for the primary entries are:
57!
58!                   Y=ERF(X)     (or   Y=DERF(X)),
59!
60!                   Y=ERFC(X)    (or   Y=DERFC(X)),
61!   and
62!                   Y=ERFCX(X)   (or   Y=DERFCX(X)).
63!
64!   The routine  CALERF  is intended for internal packet use only,
65!   all computations within the packet being concentrated in this
66!   routine.  The function subprograms invoke  CALERF  with the
67!   statement
68!
69!          CALL CALERF(ARG,RESULT,JINT)
70!
71!   where the parameter usage is as follows
72!
73!      Function                     Parameters for CALERF
74!       call              ARG                  Result          JINT
75!
76!     ERF(ARG)      ANY REAL ARGUMENT         ERF(ARG)          0
77!     ERFC(ARG)     ABS(ARG) .LT. XBIG        ERFC(ARG)         1
78!     ERFCX(ARG)    XNEG .LT. ARG .LT. XMAX   ERFCX(ARG)        2
79!
80!   The main computation evaluates near-minimax approximations
81!   from "Rational Chebyshev approximations for the error function"
82!   by W. J. Cody, Math. Comp., 1969, PP. 631-638.  This
83!   transportable program uses rational functions that theoretically
84!   approximate  erf(x)  and  erfc(x)  to at least 18 significant
85!   decimal digits.  The accuracy achieved depends on the arithmetic
86!   system, the compiler, the intrinsic functions, and proper
87!   selection of the machine-dependent constants.
88!
89!*******************************************************************
90!*******************************************************************
91!
92! Explanation of machine-dependent constants
93!
94!   XMIN   = the smallest positive floating-point number.
95!   XINF   = the largest positive finite floating-point number.
96!   XNEG   = the largest negative argument acceptable to ERFCX;
97!            the negative of the solution to the equation
98!            2*exp(x*x) = XINF.
99!   XSMALL = argument below which erf(x) may be represented by
100!            2*x/sqrt(pi)  and above which  x*x  will not underflow.
101!            A conservative value is the largest machine number X
102!            such that   1.0 + X = 1.0   to machine precision.
103!   XBIG   = largest argument acceptable to ERFC;  solution to
104!            the equation:  W(x) * (1-0.5/x**2) = XMIN,  where
105!            W(x) = exp(-x*x)/[x*sqrt(pi)].
106!   XHUGE  = argument above which  1.0 - 1/(2*x*x) = 1.0  to
107!            machine precision.  A conservative value is
108!            1/[2*sqrt(XSMALL)]
109!   XMAX   = largest acceptable argument to ERFCX; the minimum
110!            of XINF and 1/[sqrt(pi)*XMIN].
111!
112!   Approximate values for some important machines are:
113!
114!                          XMIN       XINF        XNEG     XSMALL
115!
116!  CDC 7600      (S.P.)  3.13E-294   1.26E+322   -27.220  7.11E-15
117!  CRAY-1        (S.P.)  4.58E-2467  5.45E+2465  -75.345  7.11E-15
118!  IEEE (IBM/XT,
119!    SUN, etc.)  (S.P.)  1.18E-38    3.40E+38     -9.382  5.96E-8
120!  IEEE (IBM/XT,
121!    SUN, etc.)  (D.P.)  2.23D-308   1.79D+308   -26.628  1.11D-16
122!  IBM 195       (D.P.)  5.40D-79    7.23E+75    -13.190  1.39D-17
123!  UNIVAC 1108   (D.P.)  2.78D-309   8.98D+307   -26.615  1.73D-18
124!  VAX D-Format  (D.P.)  2.94D-39    1.70D+38     -9.345  1.39D-17
125!  VAX G-Format  (D.P.)  5.56D-309   8.98D+307   -26.615  1.11D-16
126!
127!
128!                          XBIG       XHUGE       XMAX
129!
130!  CDC 7600      (S.P.)  25.922      8.39E+6     1.80X+293
131!  CRAY-1        (S.P.)  75.326      8.39E+6     5.45E+2465
132!  IEEE (IBM/XT,
133!    SUN, etc.)  (S.P.)   9.194      2.90E+3     4.79E+37
134!  IEEE (IBM/XT,
135!    SUN, etc.)  (D.P.)  26.543      6.71D+7     2.53D+307
136!  IBM 195       (D.P.)  13.306      1.90D+8     7.23E+75
137!  UNIVAC 1108   (D.P.)  26.582      5.37D+8     8.98D+307
138!  VAX D-Format  (D.P.)   9.269      1.90D+8     1.70D+38
139!  VAX G-Format  (D.P.)  26.569      6.71D+7     8.98D+307
140!
141!*******************************************************************
142!*******************************************************************
143!
144! Error returns
145!
146!  The program returns  ERFC = 0      for  ARG .GE. XBIG;
147!
148!                       ERFCX = XINF  for  ARG .LT. XNEG;
149!      and
150!                       ERFCX = 0     for  ARG .GE. XMAX.
151!
152!
153! Intrinsic functions required are:
154!
155!     ABS, AINT, EXP
156!
157!
158!  Author: W. J. Cody
159!          Mathematics and Computer Science Division
160!          Argonne National Laboratory
161!          Argonne, IL 60439
162!
163!  Latest modification: March 19, 1990
164!
165!------------------------------------------------------------------
166
167SUBROUTINE CALERF_r8(ARG, RESULT, JINT)
168
169   !------------------------------------------------------------------
170   !  This version uses 8-byte reals
171   !------------------------------------------------------------------
172   integer, parameter :: rk = r8
173
174   ! arguments
175   real(rk), intent(in)  :: arg
176   integer,  intent(in)  :: jint
177   real(rk), intent(out) :: result
178
179   ! local variables
180   INTEGER :: I
181
182   real(rk) :: X, Y, YSQ, XNUM, XDEN, DEL
183
184   !------------------------------------------------------------------
185   !  Mathematical constants
186   !------------------------------------------------------------------
187   real(rk), parameter :: ZERO   = 0.0E0_rk
188   real(rk), parameter :: FOUR   = 4.0E0_rk
189   real(rk), parameter :: ONE    = 1.0E0_rk
190   real(rk), parameter :: HALF   = 0.5E0_rk
191   real(rk), parameter :: TWO    = 2.0E0_rk
192   real(rk), parameter :: SQRPI  = 5.6418958354775628695E-1_rk
193   real(rk), parameter :: THRESH = 0.46875E0_rk
194   real(rk), parameter :: SIXTEN = 16.0E0_rk
195
196!------------------------------------------------------------------
197!  Machine-dependent constants: IEEE single precision values
198!------------------------------------------------------------------
199!S      real, parameter :: XINF   =  3.40E+38
200!S      real, parameter :: XNEG   = -9.382E0
201!S      real, parameter :: XSMALL =  5.96E-8
202!S      real, parameter :: XBIG   =  9.194E0
203!S      real, parameter :: XHUGE  =  2.90E3
204!S      real, parameter :: XMAX   =  4.79E37
205
206   !------------------------------------------------------------------
207   !  Machine-dependent constants: IEEE double precision values
208   !------------------------------------------------------------------
209   real(rk), parameter :: XINF   =   1.79E308_r8
210   real(rk), parameter :: XNEG   = -26.628E0_r8
211   real(rk), parameter :: XSMALL =   1.11E-16_r8
212   real(rk), parameter :: XBIG   =  26.543E0_r8
213   real(rk), parameter :: XHUGE  =   6.71E7_r8
214   real(rk), parameter :: XMAX   =   2.53E307_r8
215
216   !------------------------------------------------------------------
217   !  Coefficients for approximation to  erf  in first interval
218   !------------------------------------------------------------------
219   real(rk), parameter :: A(5) = (/ 3.16112374387056560E00_rk, 1.13864154151050156E02_rk, &
220                                    3.77485237685302021E02_rk, 3.20937758913846947E03_rk, &
221                                    1.85777706184603153E-1_rk /)
222   real(rk), parameter :: B(4) = (/ 2.36012909523441209E01_rk, 2.44024637934444173E02_rk, &
223                                    1.28261652607737228E03_rk, 2.84423683343917062E03_rk /)
224
225   !------------------------------------------------------------------
226   !  Coefficients for approximation to  erfc  in second interval
227   !------------------------------------------------------------------
228   real(rk), parameter :: C(9) = (/ 5.64188496988670089E-1_rk, 8.88314979438837594E00_rk, &
229                                    6.61191906371416295E01_rk, 2.98635138197400131E02_rk, &
230                                    8.81952221241769090E02_rk, 1.71204761263407058E03_rk, &
231                                    2.05107837782607147E03_rk, 1.23033935479799725E03_rk, &
232                                    2.15311535474403846E-8_rk /)
233   real(rk), parameter :: D(8) = (/ 1.57449261107098347E01_rk, 1.17693950891312499E02_rk, &
234                                    5.37181101862009858E02_rk, 1.62138957456669019E03_rk, &
235                                    3.29079923573345963E03_rk, 4.36261909014324716E03_rk, &
236                                    3.43936767414372164E03_rk, 1.23033935480374942E03_rk /)
237
238   !------------------------------------------------------------------
239   !  Coefficients for approximation to  erfc  in third interval
240   !------------------------------------------------------------------
241   real(rk), parameter :: P(6) = (/ 3.05326634961232344E-1_rk, 3.60344899949804439E-1_rk, &
242                                    1.25781726111229246E-1_rk, 1.60837851487422766E-2_rk, &
243                                    6.58749161529837803E-4_rk, 1.63153871373020978E-2_rk /)
244   real(rk), parameter :: Q(5) = (/ 2.56852019228982242E00_rk, 1.87295284992346047E00_rk, &
245                                    5.27905102951428412E-1_rk, 6.05183413124413191E-2_rk, &
246                                    2.33520497626869185E-3_rk /)
247
248   !------------------------------------------------------------------
249   X = ARG
250   Y = ABS(X)
251   IF (Y .LE. THRESH) THEN
252      !------------------------------------------------------------------
253      !  Evaluate  erf  for  |X| <= 0.46875
254      !------------------------------------------------------------------
255      YSQ = ZERO
256      IF (Y .GT. XSMALL) YSQ = Y * Y
257      XNUM = A(5)*YSQ
258      XDEN = YSQ
259      DO I = 1, 3
260         XNUM = (XNUM + A(I)) * YSQ
261         XDEN = (XDEN + B(I)) * YSQ
262      end do
263      RESULT = X * (XNUM + A(4)) / (XDEN + B(4))
264      IF (JINT .NE. 0) RESULT = ONE - RESULT
265      IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT
266      GO TO 80
267   ELSE IF (Y .LE. FOUR) THEN
268      !------------------------------------------------------------------
269      !  Evaluate  erfc  for 0.46875 <= |X| <= 4.0
270      !------------------------------------------------------------------
271      XNUM = C(9)*Y
272      XDEN = Y
273      DO I = 1, 7
274         XNUM = (XNUM + C(I)) * Y
275         XDEN = (XDEN + D(I)) * Y
276      end do
277      RESULT = (XNUM + C(8)) / (XDEN + D(8))
278      IF (JINT .NE. 2) THEN
279         YSQ = AINT(Y*SIXTEN)/SIXTEN
280         DEL = (Y-YSQ)*(Y+YSQ)
281         RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
282      END IF
283   ELSE
284      !------------------------------------------------------------------
285      !  Evaluate  erfc  for |X| > 4.0
286      !------------------------------------------------------------------
287      RESULT = ZERO
288      IF (Y .GE. XBIG) THEN
289         IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 30
290         IF (Y .GE. XHUGE) THEN
291            RESULT = SQRPI / Y
292            GO TO 30
293         END IF
294      END IF
295      YSQ = ONE / (Y * Y)
296      XNUM = P(6)*YSQ
297      XDEN = YSQ
298      DO I = 1, 4
299         XNUM = (XNUM + P(I)) * YSQ
300         XDEN = (XDEN + Q(I)) * YSQ
301      end do
302      RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5))
303      RESULT = (SQRPI -  RESULT) / Y
304      IF (JINT .NE. 2) THEN
305         YSQ = AINT(Y*SIXTEN)/SIXTEN
306         DEL = (Y-YSQ)*(Y+YSQ)
307         RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
308      END IF
309   END IF
31030 continue
311   !------------------------------------------------------------------
312   !  Fix up for negative argument, erf, etc.
313   !------------------------------------------------------------------
314   IF (JINT .EQ. 0) THEN
315      RESULT = (HALF - RESULT) + HALF
316      IF (X .LT. ZERO) RESULT = -RESULT
317   ELSE IF (JINT .EQ. 1) THEN
318      IF (X .LT. ZERO) RESULT = TWO - RESULT
319   ELSE
320      IF (X .LT. ZERO) THEN
321         IF (X .LT. XNEG) THEN
322            RESULT = XINF
323         ELSE
324            YSQ = AINT(X*SIXTEN)/SIXTEN
325            DEL = (X-YSQ)*(X+YSQ)
326            Y = EXP(YSQ*YSQ) * EXP(DEL)
327            RESULT = (Y+Y) - RESULT
328         END IF
329      END IF
330   END IF
33180 continue
332end SUBROUTINE CALERF_r8
333
334!------------------------------------------------------------------------------------------
335
336SUBROUTINE CALERF_r4(ARG, RESULT, JINT)
337
338   !------------------------------------------------------------------
339   !  This version uses 4-byte reals
340   !------------------------------------------------------------------
341   integer, parameter :: rk = r4
342
343   ! arguments
344   real(rk), intent(in)  :: arg
345   integer,  intent(in)  :: jint
346   real(rk), intent(out) :: result
347
348   ! local variables
349   INTEGER :: I
350
351   real(rk) :: X, Y, YSQ, XNUM, XDEN, DEL
352
353   !------------------------------------------------------------------
354   !  Mathematical constants
355   !------------------------------------------------------------------
356   real(rk), parameter :: ZERO   = 0.0E0_rk
357   real(rk), parameter :: FOUR   = 4.0E0_rk
358   real(rk), parameter :: ONE    = 1.0E0_rk
359   real(rk), parameter :: HALF   = 0.5E0_rk
360   real(rk), parameter :: TWO    = 2.0E0_rk
361   real(rk), parameter :: SQRPI  = 5.6418958354775628695E-1_rk
362   real(rk), parameter :: THRESH = 0.46875E0_rk
363   real(rk), parameter :: SIXTEN = 16.0E0_rk
364
365   !------------------------------------------------------------------
366   !  Machine-dependent constants: IEEE single precision values
367   !------------------------------------------------------------------
368   real(rk), parameter :: XINF   =  3.40E+38_r4
369   real(rk), parameter :: XNEG   = -9.382E0_r4
370   real(rk), parameter :: XSMALL =  5.96E-8_r4
371   real(rk), parameter :: XBIG   =  9.194E0_r4
372   real(rk), parameter :: XHUGE  =  2.90E3_r4
373   real(rk), parameter :: XMAX   =  4.79E37_r4
374
375   !------------------------------------------------------------------
376   !  Coefficients for approximation to  erf  in first interval
377   !------------------------------------------------------------------
378   real(rk), parameter :: A(5) = (/ 3.16112374387056560E00_rk, 1.13864154151050156E02_rk, &
379                                    3.77485237685302021E02_rk, 3.20937758913846947E03_rk, &
380                                    1.85777706184603153E-1_rk /)
381   real(rk), parameter :: B(4) = (/ 2.36012909523441209E01_rk, 2.44024637934444173E02_rk, &
382                                    1.28261652607737228E03_rk, 2.84423683343917062E03_rk /)
383
384   !------------------------------------------------------------------
385   !  Coefficients for approximation to  erfc  in second interval
386   !------------------------------------------------------------------
387   real(rk), parameter :: C(9) = (/ 5.64188496988670089E-1_rk, 8.88314979438837594E00_rk, &
388                                    6.61191906371416295E01_rk, 2.98635138197400131E02_rk, &
389                                    8.81952221241769090E02_rk, 1.71204761263407058E03_rk, &
390                                    2.05107837782607147E03_rk, 1.23033935479799725E03_rk, &
391                                    2.15311535474403846E-8_rk /)
392   real(rk), parameter :: D(8) = (/ 1.57449261107098347E01_rk, 1.17693950891312499E02_rk, &
393                                    5.37181101862009858E02_rk, 1.62138957456669019E03_rk, &
394                                    3.29079923573345963E03_rk, 4.36261909014324716E03_rk, &
395                                    3.43936767414372164E03_rk, 1.23033935480374942E03_rk /)
396
397   !------------------------------------------------------------------
398   !  Coefficients for approximation to  erfc  in third interval
399   !------------------------------------------------------------------
400   real(rk), parameter :: P(6) = (/ 3.05326634961232344E-1_rk, 3.60344899949804439E-1_rk, &
401                                    1.25781726111229246E-1_rk, 1.60837851487422766E-2_rk, &
402                                    6.58749161529837803E-4_rk, 1.63153871373020978E-2_rk /)
403   real(rk), parameter :: Q(5) = (/ 2.56852019228982242E00_rk, 1.87295284992346047E00_rk, &
404                                    5.27905102951428412E-1_rk, 6.05183413124413191E-2_rk, &
405                                    2.33520497626869185E-3_rk /)
406
407   !------------------------------------------------------------------
408   X = ARG
409   Y = ABS(X)
410   IF (Y .LE. THRESH) THEN
411      !------------------------------------------------------------------
412      !  Evaluate  erf  for  |X| <= 0.46875
413      !------------------------------------------------------------------
414      YSQ = ZERO
415      IF (Y .GT. XSMALL) YSQ = Y * Y
416      XNUM = A(5)*YSQ
417      XDEN = YSQ
418      DO I = 1, 3
419         XNUM = (XNUM + A(I)) * YSQ
420         XDEN = (XDEN + B(I)) * YSQ
421      end do
422      RESULT = X * (XNUM + A(4)) / (XDEN + B(4))
423      IF (JINT .NE. 0) RESULT = ONE - RESULT
424      IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT
425      GO TO 80
426   ELSE IF (Y .LE. FOUR) THEN
427      !------------------------------------------------------------------
428      !  Evaluate  erfc  for 0.46875 <= |X| <= 4.0
429      !------------------------------------------------------------------
430      XNUM = C(9)*Y
431      XDEN = Y
432      DO I = 1, 7
433         XNUM = (XNUM + C(I)) * Y
434         XDEN = (XDEN + D(I)) * Y
435      end do
436      RESULT = (XNUM + C(8)) / (XDEN + D(8))
437      IF (JINT .NE. 2) THEN
438         YSQ = AINT(Y*SIXTEN)/SIXTEN
439         DEL = (Y-YSQ)*(Y+YSQ)
440         RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
441      END IF
442   ELSE
443      !------------------------------------------------------------------
444      !  Evaluate  erfc  for |X| > 4.0
445      !------------------------------------------------------------------
446      RESULT = ZERO
447      IF (Y .GE. XBIG) THEN
448         IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 30
449         IF (Y .GE. XHUGE) THEN
450            RESULT = SQRPI / Y
451            GO TO 30
452         END IF
453      END IF
454      YSQ = ONE / (Y * Y)
455      XNUM = P(6)*YSQ
456      XDEN = YSQ
457      DO I = 1, 4
458         XNUM = (XNUM + P(I)) * YSQ
459         XDEN = (XDEN + Q(I)) * YSQ
460      end do
461      RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5))
462      RESULT = (SQRPI -  RESULT) / Y
463      IF (JINT .NE. 2) THEN
464         YSQ = AINT(Y*SIXTEN)/SIXTEN
465         DEL = (Y-YSQ)*(Y+YSQ)
466         RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
467      END IF
468   END IF
46930 continue
470   !------------------------------------------------------------------
471   !  Fix up for negative argument, erf, etc.
472   !------------------------------------------------------------------
473   IF (JINT .EQ. 0) THEN
474      RESULT = (HALF - RESULT) + HALF
475      IF (X .LT. ZERO) RESULT = -RESULT
476   ELSE IF (JINT .EQ. 1) THEN
477      IF (X .LT. ZERO) RESULT = TWO - RESULT
478   ELSE
479      IF (X .LT. ZERO) THEN
480         IF (X .LT. XNEG) THEN
481            RESULT = XINF
482         ELSE
483            YSQ = AINT(X*SIXTEN)/SIXTEN
484            DEL = (X-YSQ)*(X+YSQ)
485            Y = EXP(YSQ*YSQ) * EXP(DEL)
486            RESULT = (Y+Y) - RESULT
487         END IF
488      END IF
489   END IF
49080 continue
491end SUBROUTINE CALERF_r4
492
493!------------------------------------------------------------------------------------------
494
495FUNCTION DERF(X)
496!--------------------------------------------------------------------
497!
498! This subprogram computes approximate values for erf(x).
499!   (see comments heading CALERF).
500!
501!   Author/date: W. J. Cody, January 8, 1985
502!
503!--------------------------------------------------------------------
504   integer, parameter :: rk = r8 ! 8 byte real
505
506   ! argument
507   real(rk), intent(in) :: X
508
509   ! return value
510   real(rk) :: DERF
511
512   ! local variables
513   INTEGER :: JINT = 0
514   !------------------------------------------------------------------
515
516   CALL CALERF_r8(X, DERF, JINT)
517END FUNCTION DERF
518
519!------------------------------------------------------------------------------------------
520
521FUNCTION ERF_r4(X)
522!--------------------------------------------------------------------
523!
524! This subprogram computes approximate values for erf(x).
525!   (see comments heading CALERF).
526!
527!   Author/date: W. J. Cody, January 8, 1985
528!
529!--------------------------------------------------------------------
530   integer, parameter :: rk = r4 ! 4 byte real
531
532   ! argument
533   real(rk), intent(in) :: X
534
535   ! return value
536   real(rk) :: ERF_r4
537
538   ! local variables
539   INTEGER :: JINT = 0
540   !------------------------------------------------------------------
541
542   CALL CALERF_r4(X, ERF_r4, JINT)
543END FUNCTION ERF_r4
544
545!------------------------------------------------------------------------------------------
546
547FUNCTION DERFC(X)
548!--------------------------------------------------------------------
549!
550! This subprogram computes approximate values for erfc(x).
551!   (see comments heading CALERF).
552!
553!   Author/date: W. J. Cody, January 8, 1985
554!
555!--------------------------------------------------------------------
556   integer, parameter :: rk = r8 ! 8 byte real
557
558   ! argument
559   real(rk), intent(in) :: X
560
561   ! return value
562   real(rk) :: DERFC
563
564   ! local variables
565   INTEGER :: JINT = 1
566   !------------------------------------------------------------------
567
568   CALL CALERF_r8(X, DERFC, JINT)
569END FUNCTION DERFC
570
571!------------------------------------------------------------------------------------------
572
573FUNCTION ERFC_r4(X)
574!--------------------------------------------------------------------
575!
576! This subprogram computes approximate values for erfc(x).
577!   (see comments heading CALERF).
578!
579!   Author/date: W. J. Cody, January 8, 1985
580!
581!--------------------------------------------------------------------
582   integer, parameter :: rk = r4 ! 4 byte real
583
584   ! argument
585   real(rk), intent(in) :: X
586
587   ! return value
588   real(rk) :: ERFC_r4
589
590   ! local variables
591   INTEGER :: JINT = 1
592   !------------------------------------------------------------------
593
594   CALL CALERF_r4(X, ERFC_r4, JINT)
595END FUNCTION ERFC_r4
596
597!------------------------------------------------------------------------------------------
598
599FUNCTION DERFCX(X)
600!--------------------------------------------------------------------
601!
602! This subprogram computes approximate values for exp(x*x) * erfc(x).
603!   (see comments heading CALERF).
604!
605!   Author/date: W. J. Cody, March 30, 1987
606!
607!--------------------------------------------------------------------
608   integer, parameter :: rk = r8 ! 8 byte real
609
610   ! argument
611   real(rk), intent(in) :: X
612
613   ! return value
614   real(rk) :: DERFCX
615
616   ! local variables
617   INTEGER :: JINT = 2
618   !------------------------------------------------------------------
619
620   CALL CALERF_r8(X, DERFCX, JINT)
621END FUNCTION DERFCX
622
623!------------------------------------------------------------------------------------------
624
625FUNCTION ERFCX_R4(X)
626!--------------------------------------------------------------------
627!
628! This subprogram computes approximate values for exp(x*x) * erfc(x).
629!   (see comments heading CALERF).
630!
631!   Author/date: W. J. Cody, March 30, 1987
632!
633!--------------------------------------------------------------------
634   integer, parameter :: rk = r4 ! 8 byte real
635
636   ! argument
637   real(rk), intent(in) :: X
638
639   ! return value
640   real(rk) :: ERFCX_R4
641
642   ! local variables
643   INTEGER :: JINT = 2
644   !------------------------------------------------------------------
645
646   CALL CALERF_r4(X, ERFCX_R4, JINT)
647END FUNCTION ERFCX_R4
648
649!------------------------------------------------------------------------------------------
650
651end module error_function
Note: See TracBrowser for help on using the repository browser.