3 ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
7 ! This module provides generic interfaces for functions that evaluate
8 ! erf(x), erfc(x), and exp(x*x)*erfc(x) in either single or double precision.
15 public :: erf, erfc, erfcx
18 module procedure erf_r4
23 module procedure erfc_r4
24 module procedure derfc
28 module procedure erfcx_r4
29 module procedure derfcx
33 integer, parameter :: r4 = selected_real_kind(6) ! 4 byte real
34 integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real
38 !------------------------------------------------------------------
40 ! 6 December 2006 -- B. Eaton
41 ! The following comments are from the original version of CALERF.
42 ! The only changes in implementing this module are that the function
43 ! names previously used for the single precision versions have been
44 ! adopted for the new generic interfaces. To support these interfaces
45 ! there is now both a single precision version (calerf_r4) and a
46 ! double precision version (calerf_r8) of CALERF below. These versions
47 ! are hardcoded to use IEEE arithmetic.
49 !------------------------------------------------------------------
51 ! This packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x)
52 ! for a real argument x. It contains three FUNCTION type
53 ! subprograms: ERF, ERFC, and ERFCX (or DERF, DERFC, and DERFCX),
54 ! and one SUBROUTINE type subprogram, CALERF. The calling
55 ! statements for the primary entries are:
57 ! Y=ERF(X) (or Y=DERF(X)),
59 ! Y=ERFC(X) (or Y=DERFC(X)),
61 ! Y=ERFCX(X) (or Y=DERFCX(X)).
63 ! The routine CALERF is intended for internal packet use only,
64 ! all computations within the packet being concentrated in this
65 ! routine. The function subprograms invoke CALERF with the
68 ! CALL CALERF(ARG,RESULT,JINT)
70 ! where the parameter usage is as follows
72 ! Function Parameters for CALERF
73 ! call ARG Result JINT
75 ! ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0
76 ! ERFC(ARG) ABS(ARG) .LT. XBIG ERFC(ARG) 1
77 ! ERFCX(ARG) XNEG .LT. ARG .LT. XMAX ERFCX(ARG) 2
79 ! The main computation evaluates near-minimax approximations
80 ! from "Rational Chebyshev approximations for the error function"
81 ! by W. J. Cody, Math. Comp., 1969, PP. 631-638. This
82 ! transportable program uses rational functions that theoretically
83 ! approximate erf(x) and erfc(x) to at least 18 significant
84 ! decimal digits. The accuracy achieved depends on the arithmetic
85 ! system, the compiler, the intrinsic functions, and proper
86 ! selection of the machine-dependent constants.
88 !*******************************************************************
89 !*******************************************************************
91 ! Explanation of machine-dependent constants
93 ! XMIN = the smallest positive floating-point number.
94 ! XINF = the largest positive finite floating-point number.
95 ! XNEG = the largest negative argument acceptable to ERFCX;
96 ! the negative of the solution to the equation
98 ! XSMALL = argument below which erf(x) may be represented by
99 ! 2*x/sqrt(pi) and above which x*x will not underflow.
100 ! A conservative value is the largest machine number X
101 ! such that 1.0 + X = 1.0 to machine precision.
102 ! XBIG = largest argument acceptable to ERFC; solution to
103 ! the equation: W(x) * (1-0.5/x**2) = XMIN, where
104 ! W(x) = exp(-x*x)/[x*sqrt(pi)].
105 ! XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to
106 ! machine precision. A conservative value is
108 ! XMAX = largest acceptable argument to ERFCX; the minimum
109 ! of XINF and 1/[sqrt(pi)*XMIN].
111 ! Approximate values for some important machines are:
113 ! XMIN XINF XNEG XSMALL
115 ! CDC 7600 (S.P.) 3.13E-294 1.26E+322 -27.220 7.11E-15
116 ! CRAY-1 (S.P.) 4.58E-2467 5.45E+2465 -75.345 7.11E-15
118 ! SUN, etc.) (S.P.) 1.18E-38 3.40E+38 -9.382 5.96E-8
120 ! SUN, etc.) (D.P.) 2.23D-308 1.79D+308 -26.628 1.11D-16
121 ! IBM 195 (D.P.) 5.40D-79 7.23E+75 -13.190 1.39D-17
122 ! UNIVAC 1108 (D.P.) 2.78D-309 8.98D+307 -26.615 1.73D-18
123 ! VAX D-Format (D.P.) 2.94D-39 1.70D+38 -9.345 1.39D-17
124 ! VAX G-Format (D.P.) 5.56D-309 8.98D+307 -26.615 1.11D-16
129 ! CDC 7600 (S.P.) 25.922 8.39E+6 1.80X+293
130 ! CRAY-1 (S.P.) 75.326 8.39E+6 5.45E+2465
132 ! SUN, etc.) (S.P.) 9.194 2.90E+3 4.79E+37
134 ! SUN, etc.) (D.P.) 26.543 6.71D+7 2.53D+307
135 ! IBM 195 (D.P.) 13.306 1.90D+8 7.23E+75
136 ! UNIVAC 1108 (D.P.) 26.582 5.37D+8 8.98D+307
137 ! VAX D-Format (D.P.) 9.269 1.90D+8 1.70D+38
138 ! VAX G-Format (D.P.) 26.569 6.71D+7 8.98D+307
140 !*******************************************************************
141 !*******************************************************************
145 ! The program returns ERFC = 0 for ARG .GE. XBIG;
147 ! ERFCX = XINF for ARG .LT. XNEG;
149 ! ERFCX = 0 for ARG .GE. XMAX.
152 ! Intrinsic functions required are:
158 ! Mathematics and Computer Science Division
159 ! Argonne National Laboratory
162 ! Latest modification: March 19, 1990
164 !------------------------------------------------------------------
166 SUBROUTINE CALERF_r8(ARG, RESULT, JINT)
168 !------------------------------------------------------------------
169 ! This version uses 8-byte reals
170 !------------------------------------------------------------------
171 integer, parameter :: rk = r8
174 real(rk), intent(in) :: arg
175 integer, intent(in) :: jint
176 real(rk), intent(out) :: result
181 real(rk) :: X, Y, YSQ, XNUM, XDEN, DEL
183 !------------------------------------------------------------------
184 ! Mathematical constants
185 !------------------------------------------------------------------
186 real(rk), parameter :: ZERO = 0.0E0_rk
187 real(rk), parameter :: FOUR = 4.0E0_rk
188 real(rk), parameter :: ONE = 1.0E0_rk
189 real(rk), parameter :: HALF = 0.5E0_rk
190 real(rk), parameter :: TWO = 2.0E0_rk
191 real(rk), parameter :: SQRPI = 5.6418958354775628695E-1_rk
192 real(rk), parameter :: THRESH = 0.46875E0_rk
193 real(rk), parameter :: SIXTEN = 16.0E0_rk
195 !------------------------------------------------------------------
196 ! Machine-dependent constants: IEEE single precision values
197 !------------------------------------------------------------------
198 !S real, parameter :: XINF = 3.40E+38
199 !S real, parameter :: XNEG = -9.382E0
200 !S real, parameter :: XSMALL = 5.96E-8
201 !S real, parameter :: XBIG = 9.194E0
202 !S real, parameter :: XHUGE = 2.90E3
203 !S real, parameter :: XMAX = 4.79E37
205 !------------------------------------------------------------------
206 ! Machine-dependent constants: IEEE double precision values
207 !------------------------------------------------------------------
208 real(rk), parameter :: XINF = 1.79E308_r8
209 real(rk), parameter :: XNEG = -26.628E0_r8
210 real(rk), parameter :: XSMALL = 1.11E-16_r8
211 real(rk), parameter :: XBIG = 26.543E0_r8
212 real(rk), parameter :: XHUGE = 6.71E7_r8
213 real(rk), parameter :: XMAX = 2.53E307_r8
215 !------------------------------------------------------------------
216 ! Coefficients for approximation to erf in first interval
217 !------------------------------------------------------------------
218 real(rk), parameter :: A(5) = (/ 3.16112374387056560E00_rk, 1.13864154151050156E02_rk, &
219 3.77485237685302021E02_rk, 3.20937758913846947E03_rk, &
220 1.85777706184603153E-1_rk /)
221 real(rk), parameter :: B(4) = (/ 2.36012909523441209E01_rk, 2.44024637934444173E02_rk, &
222 1.28261652607737228E03_rk, 2.84423683343917062E03_rk /)
224 !------------------------------------------------------------------
225 ! Coefficients for approximation to erfc in second interval
226 !------------------------------------------------------------------
227 real(rk), parameter :: C(9) = (/ 5.64188496988670089E-1_rk, 8.88314979438837594E00_rk, &
228 6.61191906371416295E01_rk, 2.98635138197400131E02_rk, &
229 8.81952221241769090E02_rk, 1.71204761263407058E03_rk, &
230 2.05107837782607147E03_rk, 1.23033935479799725E03_rk, &
231 2.15311535474403846E-8_rk /)
232 real(rk), parameter :: D(8) = (/ 1.57449261107098347E01_rk, 1.17693950891312499E02_rk, &
233 5.37181101862009858E02_rk, 1.62138957456669019E03_rk, &
234 3.29079923573345963E03_rk, 4.36261909014324716E03_rk, &
235 3.43936767414372164E03_rk, 1.23033935480374942E03_rk /)
237 !------------------------------------------------------------------
238 ! Coefficients for approximation to erfc in third interval
239 !------------------------------------------------------------------
240 real(rk), parameter :: P(6) = (/ 3.05326634961232344E-1_rk, 3.60344899949804439E-1_rk, &
241 1.25781726111229246E-1_rk, 1.60837851487422766E-2_rk, &
242 6.58749161529837803E-4_rk, 1.63153871373020978E-2_rk /)
243 real(rk), parameter :: Q(5) = (/ 2.56852019228982242E00_rk, 1.87295284992346047E00_rk, &
244 5.27905102951428412E-1_rk, 6.05183413124413191E-2_rk, &
245 2.33520497626869185E-3_rk /)
247 !------------------------------------------------------------------
250 IF (Y .LE. THRESH) THEN
251 !------------------------------------------------------------------
252 ! Evaluate erf for |X| <= 0.46875
253 !------------------------------------------------------------------
255 IF (Y .GT. XSMALL) YSQ = Y * Y
259 XNUM = (XNUM + A(I)) * YSQ
260 XDEN = (XDEN + B(I)) * YSQ
262 RESULT = X * (XNUM + A(4)) / (XDEN + B(4))
263 IF (JINT .NE. 0) RESULT = ONE - RESULT
264 IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT
266 ELSE IF (Y .LE. FOUR) THEN
267 !------------------------------------------------------------------
268 ! Evaluate erfc for 0.46875 <= |X| <= 4.0
269 !------------------------------------------------------------------
273 XNUM = (XNUM + C(I)) * Y
274 XDEN = (XDEN + D(I)) * Y
276 RESULT = (XNUM + C(8)) / (XDEN + D(8))
277 IF (JINT .NE. 2) THEN
278 YSQ = AINT(Y*SIXTEN)/SIXTEN
279 DEL = (Y-YSQ)*(Y+YSQ)
280 RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
283 !------------------------------------------------------------------
284 ! Evaluate erfc for |X| > 4.0
285 !------------------------------------------------------------------
287 IF (Y .GE. XBIG) THEN
288 IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 30
289 IF (Y .GE. XHUGE) THEN
298 XNUM = (XNUM + P(I)) * YSQ
299 XDEN = (XDEN + Q(I)) * YSQ
301 RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5))
302 RESULT = (SQRPI - RESULT) / Y
303 IF (JINT .NE. 2) THEN
304 YSQ = AINT(Y*SIXTEN)/SIXTEN
305 DEL = (Y-YSQ)*(Y+YSQ)
306 RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
310 !------------------------------------------------------------------
311 ! Fix up for negative argument, erf, etc.
312 !------------------------------------------------------------------
313 IF (JINT .EQ. 0) THEN
314 RESULT = (HALF - RESULT) + HALF
315 IF (X .LT. ZERO) RESULT = -RESULT
316 ELSE IF (JINT .EQ. 1) THEN
317 IF (X .LT. ZERO) RESULT = TWO - RESULT
319 IF (X .LT. ZERO) THEN
320 IF (X .LT. XNEG) THEN
323 YSQ = AINT(X*SIXTEN)/SIXTEN
324 DEL = (X-YSQ)*(X+YSQ)
325 Y = EXP(YSQ*YSQ) * EXP(DEL)
326 RESULT = (Y+Y) - RESULT
331 end SUBROUTINE CALERF_r8
333 !------------------------------------------------------------------------------------------
335 SUBROUTINE CALERF_r4(ARG, RESULT, JINT)
337 !------------------------------------------------------------------
338 ! This version uses 4-byte reals
339 !------------------------------------------------------------------
340 integer, parameter :: rk = r4
343 real(rk), intent(in) :: arg
344 integer, intent(in) :: jint
345 real(rk), intent(out) :: result
350 real(rk) :: X, Y, YSQ, XNUM, XDEN, DEL
352 !------------------------------------------------------------------
353 ! Mathematical constants
354 !------------------------------------------------------------------
355 real(rk), parameter :: ZERO = 0.0E0_rk
356 real(rk), parameter :: FOUR = 4.0E0_rk
357 real(rk), parameter :: ONE = 1.0E0_rk
358 real(rk), parameter :: HALF = 0.5E0_rk
359 real(rk), parameter :: TWO = 2.0E0_rk
360 real(rk), parameter :: SQRPI = 5.6418958354775628695E-1_rk
361 real(rk), parameter :: THRESH = 0.46875E0_rk
362 real(rk), parameter :: SIXTEN = 16.0E0_rk
364 !------------------------------------------------------------------
365 ! Machine-dependent constants: IEEE single precision values
366 !------------------------------------------------------------------
367 real(rk), parameter :: XINF = 3.40E+38_r4
368 real(rk), parameter :: XNEG = -9.382E0_r4
369 real(rk), parameter :: XSMALL = 5.96E-8_r4
370 real(rk), parameter :: XBIG = 9.194E0_r4
371 real(rk), parameter :: XHUGE = 2.90E3_r4
372 real(rk), parameter :: XMAX = 4.79E37_r4
374 !------------------------------------------------------------------
375 ! Coefficients for approximation to erf in first interval
376 !------------------------------------------------------------------
377 real(rk), parameter :: A(5) = (/ 3.16112374387056560E00_rk, 1.13864154151050156E02_rk, &
378 3.77485237685302021E02_rk, 3.20937758913846947E03_rk, &
379 1.85777706184603153E-1_rk /)
380 real(rk), parameter :: B(4) = (/ 2.36012909523441209E01_rk, 2.44024637934444173E02_rk, &
381 1.28261652607737228E03_rk, 2.84423683343917062E03_rk /)
383 !------------------------------------------------------------------
384 ! Coefficients for approximation to erfc in second interval
385 !------------------------------------------------------------------
386 real(rk), parameter :: C(9) = (/ 5.64188496988670089E-1_rk, 8.88314979438837594E00_rk, &
387 6.61191906371416295E01_rk, 2.98635138197400131E02_rk, &
388 8.81952221241769090E02_rk, 1.71204761263407058E03_rk, &
389 2.05107837782607147E03_rk, 1.23033935479799725E03_rk, &
390 2.15311535474403846E-8_rk /)
391 real(rk), parameter :: D(8) = (/ 1.57449261107098347E01_rk, 1.17693950891312499E02_rk, &
392 5.37181101862009858E02_rk, 1.62138957456669019E03_rk, &
393 3.29079923573345963E03_rk, 4.36261909014324716E03_rk, &
394 3.43936767414372164E03_rk, 1.23033935480374942E03_rk /)
396 !------------------------------------------------------------------
397 ! Coefficients for approximation to erfc in third interval
398 !------------------------------------------------------------------
399 real(rk), parameter :: P(6) = (/ 3.05326634961232344E-1_rk, 3.60344899949804439E-1_rk, &
400 1.25781726111229246E-1_rk, 1.60837851487422766E-2_rk, &
401 6.58749161529837803E-4_rk, 1.63153871373020978E-2_rk /)
402 real(rk), parameter :: Q(5) = (/ 2.56852019228982242E00_rk, 1.87295284992346047E00_rk, &
403 5.27905102951428412E-1_rk, 6.05183413124413191E-2_rk, &
404 2.33520497626869185E-3_rk /)
406 !------------------------------------------------------------------
409 IF (Y .LE. THRESH) THEN
410 !------------------------------------------------------------------
411 ! Evaluate erf for |X| <= 0.46875
412 !------------------------------------------------------------------
414 IF (Y .GT. XSMALL) YSQ = Y * Y
418 XNUM = (XNUM + A(I)) * YSQ
419 XDEN = (XDEN + B(I)) * YSQ
421 RESULT = X * (XNUM + A(4)) / (XDEN + B(4))
422 IF (JINT .NE. 0) RESULT = ONE - RESULT
423 IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT
425 ELSE IF (Y .LE. FOUR) THEN
426 !------------------------------------------------------------------
427 ! Evaluate erfc for 0.46875 <= |X| <= 4.0
428 !------------------------------------------------------------------
432 XNUM = (XNUM + C(I)) * Y
433 XDEN = (XDEN + D(I)) * Y
435 RESULT = (XNUM + C(8)) / (XDEN + D(8))
436 IF (JINT .NE. 2) THEN
437 YSQ = AINT(Y*SIXTEN)/SIXTEN
438 DEL = (Y-YSQ)*(Y+YSQ)
439 RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
442 !------------------------------------------------------------------
443 ! Evaluate erfc for |X| > 4.0
444 !------------------------------------------------------------------
446 IF (Y .GE. XBIG) THEN
447 IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 30
448 IF (Y .GE. XHUGE) THEN
457 XNUM = (XNUM + P(I)) * YSQ
458 XDEN = (XDEN + Q(I)) * YSQ
460 RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5))
461 RESULT = (SQRPI - RESULT) / Y
462 IF (JINT .NE. 2) THEN
463 YSQ = AINT(Y*SIXTEN)/SIXTEN
464 DEL = (Y-YSQ)*(Y+YSQ)
465 RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
469 !------------------------------------------------------------------
470 ! Fix up for negative argument, erf, etc.
471 !------------------------------------------------------------------
472 IF (JINT .EQ. 0) THEN
473 RESULT = (HALF - RESULT) + HALF
474 IF (X .LT. ZERO) RESULT = -RESULT
475 ELSE IF (JINT .EQ. 1) THEN
476 IF (X .LT. ZERO) RESULT = TWO - RESULT
478 IF (X .LT. ZERO) THEN
479 IF (X .LT. XNEG) THEN
482 YSQ = AINT(X*SIXTEN)/SIXTEN
483 DEL = (X-YSQ)*(X+YSQ)
484 Y = EXP(YSQ*YSQ) * EXP(DEL)
485 RESULT = (Y+Y) - RESULT
490 end SUBROUTINE CALERF_r4
492 !------------------------------------------------------------------------------------------
495 !--------------------------------------------------------------------
497 ! This subprogram computes approximate values for erf(x).
498 ! (see comments heading CALERF).
500 ! Author/date: W. J. Cody, January 8, 1985
502 !--------------------------------------------------------------------
503 integer, parameter :: rk = r8 ! 8 byte real
506 real(rk), intent(in) :: X
513 !------------------------------------------------------------------
515 CALL CALERF_r8(X, DERF, JINT)
518 !------------------------------------------------------------------------------------------
521 !--------------------------------------------------------------------
523 ! This subprogram computes approximate values for erf(x).
524 ! (see comments heading CALERF).
526 ! Author/date: W. J. Cody, January 8, 1985
528 !--------------------------------------------------------------------
529 integer, parameter :: rk = r4 ! 4 byte real
532 real(rk), intent(in) :: X
539 !------------------------------------------------------------------
541 CALL CALERF_r4(X, ERF_r4, JINT)
544 !------------------------------------------------------------------------------------------
547 !--------------------------------------------------------------------
549 ! This subprogram computes approximate values for erfc(x).
550 ! (see comments heading CALERF).
552 ! Author/date: W. J. Cody, January 8, 1985
554 !--------------------------------------------------------------------
555 integer, parameter :: rk = r8 ! 8 byte real
558 real(rk), intent(in) :: X
565 !------------------------------------------------------------------
567 CALL CALERF_r8(X, DERFC, JINT)
570 !------------------------------------------------------------------------------------------
573 !--------------------------------------------------------------------
575 ! This subprogram computes approximate values for erfc(x).
576 ! (see comments heading CALERF).
578 ! Author/date: W. J. Cody, January 8, 1985
580 !--------------------------------------------------------------------
581 integer, parameter :: rk = r4 ! 4 byte real
584 real(rk), intent(in) :: X
591 !------------------------------------------------------------------
593 CALL CALERF_r4(X, ERFC_r4, JINT)
596 !------------------------------------------------------------------------------------------
599 !--------------------------------------------------------------------
601 ! This subprogram computes approximate values for exp(x*x) * erfc(x).
602 ! (see comments heading CALERF).
604 ! Author/date: W. J. Cody, March 30, 1987
606 !--------------------------------------------------------------------
607 integer, parameter :: rk = r8 ! 8 byte real
610 real(rk), intent(in) :: X
617 !------------------------------------------------------------------
619 CALL CALERF_r8(X, DERFCX, JINT)
622 !------------------------------------------------------------------------------------------
625 !--------------------------------------------------------------------
627 ! This subprogram computes approximate values for exp(x*x) * erfc(x).
628 ! (see comments heading CALERF).
630 ! Author/date: W. J. Cody, March 30, 1987
632 !--------------------------------------------------------------------
633 integer, parameter :: rk = r4 ! 8 byte real
636 real(rk), intent(in) :: X
643 !------------------------------------------------------------------
645 CALL CALERF_r4(X, ERFCX_R4, JINT)
646 END FUNCTION ERFCX_R4
648 !------------------------------------------------------------------------------------------
650 end module error_function