phys/module_cam_error_function.F

   1 #define WRF_PORT
   2 #define MODAL_AERO
   3 ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
   4
   5 module error_function
   6
   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.
   9
  10 implicit none
  11 private
  12 save
  13
  14 ! Public functions
  15 public :: erf, erfc, erfcx
  16
  17 interface erf
  18    module procedure erf_r4
  19    module procedure derf
  20 end interface
  21
  22 interface erfc
  23    module procedure erfc_r4
  24    module procedure derfc
  25 end interface
  26
  27 interface erfcx
  28    module procedure erfcx_r4
  29    module procedure derfcx
  30 end interface
  31
  32 ! Private variables
  33 integer, parameter :: r4 = selected_real_kind(6)  ! 4 byte real
  34 integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real
  35
  36 contains
  37
  38 !------------------------------------------------------------------
  39 !
  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.
  48 !
  49 !------------------------------------------------------------------
  50 !
  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:
  56 !
  57 !                   Y=ERF(X)     (or   Y=DERF(X)),
  58 !
  59 !                   Y=ERFC(X)    (or   Y=DERFC(X)),
  60 !   and
  61 !                   Y=ERFCX(X)   (or   Y=DERFCX(X)).
  62 !
  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
  66 !   statement
  67 !
  68 !          CALL CALERF(ARG,RESULT,JINT)
  69 !
  70 !   where the parameter usage is as follows
  71 !
  72 !      Function                     Parameters for CALERF
  73 !       call              ARG                  Result          JINT
  74 !
  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
  78 !
  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.
  87 !
  88 !*******************************************************************
  89 !*******************************************************************
  90 !
  91 ! Explanation of machine-dependent constants
  92 !
  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
  97 !            2*exp(x*x) = XINF.
  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
 107 !            1/[2*sqrt(XSMALL)]
 108 !   XMAX   = largest acceptable argument to ERFCX; the minimum
 109 !            of XINF and 1/[sqrt(pi)*XMIN].
 110 !
 111 !   Approximate values for some important machines are:
 112 !
 113 !                          XMIN       XINF        XNEG     XSMALL
 114 !
 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
 117 !  IEEE (IBM/XT,
 118 !    SUN, etc.)  (S.P.)  1.18E-38    3.40E+38     -9.382  5.96E-8
 119 !  IEEE (IBM/XT,
 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
 125 !
 126 !
 127 !                          XBIG       XHUGE       XMAX
 128 !
 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
 131 !  IEEE (IBM/XT,
 132 !    SUN, etc.)  (S.P.)   9.194      2.90E+3     4.79E+37
 133 !  IEEE (IBM/XT,
 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
 139 !
 140 !*******************************************************************
 141 !*******************************************************************
 142 !
 143 ! Error returns
 144 !
 145 !  The program returns  ERFC = 0      for  ARG .GE. XBIG;
 146 !
 147 !                       ERFCX = XINF  for  ARG .LT. XNEG;
 148 !      and
 149 !                       ERFCX = 0     for  ARG .GE. XMAX.
 150 !
 151 !
 152 ! Intrinsic functions required are:
 153 !
 154 !     ABS, AINT, EXP
 155 !
 156 !
 157 !  Author: W. J. Cody
 158 !          Mathematics and Computer Science Division
 159 !          Argonne National Laboratory
 160 !          Argonne, IL 60439
 161 !
 162 !  Latest modification: March 19, 1990
 163 !
 164 !------------------------------------------------------------------
 165
 166 SUBROUTINE CALERF_r8(ARG, RESULT, JINT)
 167
 168    !------------------------------------------------------------------
 169    !  This version uses 8-byte reals
 170    !------------------------------------------------------------------
 171    integer, parameter :: rk = r8
 172
 173    ! arguments
 174    real(rk), intent(in)  :: arg
 175    integer,  intent(in)  :: jint
 176    real(rk), intent(out) :: result
 177
 178    ! local variables
 179    INTEGER :: I
 180
 181    real(rk) :: X, Y, YSQ, XNUM, XDEN, DEL
 182
 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
 194
 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
 204
 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
 214
 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 /)
 223
 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 /)
 236
 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 /)
 246
 247    !------------------------------------------------------------------
 248    X = ARG
 249    Y = ABS(X)
 250    IF (Y .LE. THRESH) THEN
 251       !------------------------------------------------------------------
 252       !  Evaluate  erf  for  |X| <= 0.46875
 253       !------------------------------------------------------------------
 254       YSQ = ZERO
 255       IF (Y .GT. XSMALL) YSQ = Y * Y
 256       XNUM = A(5)*YSQ
 257       XDEN = YSQ
 258       DO I = 1, 3
 259          XNUM = (XNUM + A(I)) * YSQ
 260          XDEN = (XDEN + B(I)) * YSQ
 261       end do
 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
 265       GO TO 80
 266    ELSE IF (Y .LE. FOUR) THEN
 267       !------------------------------------------------------------------
 268       !  Evaluate  erfc  for 0.46875 <= |X| <= 4.0
 269       !------------------------------------------------------------------
 270       XNUM = C(9)*Y
 271       XDEN = Y
 272       DO I = 1, 7
 273          XNUM = (XNUM + C(I)) * Y
 274          XDEN = (XDEN + D(I)) * Y
 275       end do
 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
 281       END IF
 282    ELSE
 283       !------------------------------------------------------------------
 284       !  Evaluate  erfc  for |X| > 4.0
 285       !------------------------------------------------------------------
 286       RESULT = ZERO
 287       IF (Y .GE. XBIG) THEN
 288          IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 30
 289          IF (Y .GE. XHUGE) THEN
 290             RESULT = SQRPI / Y
 291             GO TO 30
 292          END IF
 293       END IF
 294       YSQ = ONE / (Y * Y)
 295       XNUM = P(6)*YSQ
 296       XDEN = YSQ
 297       DO I = 1, 4
 298          XNUM = (XNUM + P(I)) * YSQ
 299          XDEN = (XDEN + Q(I)) * YSQ
 300       end do
 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
 307       END IF
 308    END IF
 309 30 continue
 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
 318    ELSE
 319       IF (X .LT. ZERO) THEN
 320          IF (X .LT. XNEG) THEN
 321             RESULT = XINF
 322          ELSE
 323             YSQ = AINT(X*SIXTEN)/SIXTEN
 324             DEL = (X-YSQ)*(X+YSQ)
 325             Y = EXP(YSQ*YSQ) * EXP(DEL)
 326             RESULT = (Y+Y) - RESULT
 327          END IF
 328       END IF
 329    END IF
 330 80 continue
 331 end SUBROUTINE CALERF_r8
 332
 333 !------------------------------------------------------------------------------------------
 334
 335 SUBROUTINE CALERF_r4(ARG, RESULT, JINT)
 336
 337    !------------------------------------------------------------------
 338    !  This version uses 4-byte reals
 339    !------------------------------------------------------------------
 340    integer, parameter :: rk = r4
 341
 342    ! arguments
 343    real(rk), intent(in)  :: arg
 344    integer,  intent(in)  :: jint
 345    real(rk), intent(out) :: result
 346
 347    ! local variables
 348    INTEGER :: I
 349
 350    real(rk) :: X, Y, YSQ, XNUM, XDEN, DEL
 351
 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
 363
 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
 373
 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 /)
 382
 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 /)
 395
 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 /)
 405
 406    !------------------------------------------------------------------
 407    X = ARG
 408    Y = ABS(X)
 409    IF (Y .LE. THRESH) THEN
 410       !------------------------------------------------------------------
 411       !  Evaluate  erf  for  |X| <= 0.46875
 412       !------------------------------------------------------------------
 413       YSQ = ZERO
 414       IF (Y .GT. XSMALL) YSQ = Y * Y
 415       XNUM = A(5)*YSQ
 416       XDEN = YSQ
 417       DO I = 1, 3
 418          XNUM = (XNUM + A(I)) * YSQ
 419          XDEN = (XDEN + B(I)) * YSQ
 420       end do
 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
 424       GO TO 80
 425    ELSE IF (Y .LE. FOUR) THEN
 426       !------------------------------------------------------------------
 427       !  Evaluate  erfc  for 0.46875 <= |X| <= 4.0
 428       !------------------------------------------------------------------
 429       XNUM = C(9)*Y
 430       XDEN = Y
 431       DO I = 1, 7
 432          XNUM = (XNUM + C(I)) * Y
 433          XDEN = (XDEN + D(I)) * Y
 434       end do
 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
 440       END IF
 441    ELSE
 442       !------------------------------------------------------------------
 443       !  Evaluate  erfc  for |X| > 4.0
 444       !------------------------------------------------------------------
 445       RESULT = ZERO
 446       IF (Y .GE. XBIG) THEN
 447          IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 30
 448          IF (Y .GE. XHUGE) THEN
 449             RESULT = SQRPI / Y
 450             GO TO 30
 451          END IF
 452       END IF
 453       YSQ = ONE / (Y * Y)
 454       XNUM = P(6)*YSQ
 455       XDEN = YSQ
 456       DO I = 1, 4
 457          XNUM = (XNUM + P(I)) * YSQ
 458          XDEN = (XDEN + Q(I)) * YSQ
 459       end do
 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
 466       END IF
 467    END IF
 468 30 continue
 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
 477    ELSE
 478       IF (X .LT. ZERO) THEN
 479          IF (X .LT. XNEG) THEN
 480             RESULT = XINF
 481          ELSE
 482             YSQ = AINT(X*SIXTEN)/SIXTEN
 483             DEL = (X-YSQ)*(X+YSQ)
 484             Y = EXP(YSQ*YSQ) * EXP(DEL)
 485             RESULT = (Y+Y) - RESULT
 486          END IF
 487       END IF
 488    END IF
 489 80 continue
 490 end SUBROUTINE CALERF_r4
 491
 492 !------------------------------------------------------------------------------------------
 493
 494 FUNCTION DERF(X)
 495 !--------------------------------------------------------------------
 496 !
 497 ! This subprogram computes approximate values for erf(x).
 498 !   (see comments heading CALERF).
 499 !
 500 !   Author/date: W. J. Cody, January 8, 1985
 501 !
 502 !--------------------------------------------------------------------
 503    integer, parameter :: rk = r8 ! 8 byte real
 504
 505    ! argument
 506    real(rk), intent(in) :: X
 507
 508    ! return value
 509    real(rk) :: DERF
 510
 511    ! local variables
 512    INTEGER :: JINT = 0
 513    !------------------------------------------------------------------
 514
 515    CALL CALERF_r8(X, DERF, JINT)
 516 END FUNCTION DERF
 517
 518 !------------------------------------------------------------------------------------------
 519
 520 FUNCTION ERF_r4(X)
 521 !--------------------------------------------------------------------
 522 !
 523 ! This subprogram computes approximate values for erf(x).
 524 !   (see comments heading CALERF).
 525 !
 526 !   Author/date: W. J. Cody, January 8, 1985
 527 !
 528 !--------------------------------------------------------------------
 529    integer, parameter :: rk = r4 ! 4 byte real
 530
 531    ! argument
 532    real(rk), intent(in) :: X
 533
 534    ! return value
 535    real(rk) :: ERF_r4
 536
 537    ! local variables
 538    INTEGER :: JINT = 0
 539    !------------------------------------------------------------------
 540
 541    CALL CALERF_r4(X, ERF_r4, JINT)
 542 END FUNCTION ERF_r4
 543
 544 !------------------------------------------------------------------------------------------
 545
 546 FUNCTION DERFC(X)
 547 !--------------------------------------------------------------------
 548 !
 549 ! This subprogram computes approximate values for erfc(x).
 550 !   (see comments heading CALERF).
 551 !
 552 !   Author/date: W. J. Cody, January 8, 1985
 553 !
 554 !--------------------------------------------------------------------
 555    integer, parameter :: rk = r8 ! 8 byte real
 556
 557    ! argument
 558    real(rk), intent(in) :: X
 559
 560    ! return value
 561    real(rk) :: DERFC
 562
 563    ! local variables
 564    INTEGER :: JINT = 1
 565    !------------------------------------------------------------------
 566
 567    CALL CALERF_r8(X, DERFC, JINT)
 568 END FUNCTION DERFC
 569
 570 !------------------------------------------------------------------------------------------
 571
 572 FUNCTION ERFC_r4(X)
 573 !--------------------------------------------------------------------
 574 !
 575 ! This subprogram computes approximate values for erfc(x).
 576 !   (see comments heading CALERF).
 577 !
 578 !   Author/date: W. J. Cody, January 8, 1985
 579 !
 580 !--------------------------------------------------------------------
 581    integer, parameter :: rk = r4 ! 4 byte real
 582
 583    ! argument
 584    real(rk), intent(in) :: X
 585
 586    ! return value
 587    real(rk) :: ERFC_r4
 588
 589    ! local variables
 590    INTEGER :: JINT = 1
 591    !------------------------------------------------------------------
 592
 593    CALL CALERF_r4(X, ERFC_r4, JINT)
 594 END FUNCTION ERFC_r4
 595
 596 !------------------------------------------------------------------------------------------
 597
 598 FUNCTION DERFCX(X)
 599 !--------------------------------------------------------------------
 600 !
 601 ! This subprogram computes approximate values for exp(x*x) * erfc(x).
 602 !   (see comments heading CALERF).
 603 !
 604 !   Author/date: W. J. Cody, March 30, 1987
 605 !
 606 !--------------------------------------------------------------------
 607    integer, parameter :: rk = r8 ! 8 byte real
 608
 609    ! argument
 610    real(rk), intent(in) :: X
 611
 612    ! return value
 613    real(rk) :: DERFCX
 614
 615    ! local variables
 616    INTEGER :: JINT = 2
 617    !------------------------------------------------------------------
 618
 619    CALL CALERF_r8(X, DERFCX, JINT)
 620 END FUNCTION DERFCX
 621
 622 !------------------------------------------------------------------------------------------
 623
 624 FUNCTION ERFCX_R4(X)
 625 !--------------------------------------------------------------------
 626 !
 627 ! This subprogram computes approximate values for exp(x*x) * erfc(x).
 628 !   (see comments heading CALERF).
 629 !
 630 !   Author/date: W. J. Cody, March 30, 1987
 631 !
 632 !--------------------------------------------------------------------
 633    integer, parameter :: rk = r4 ! 8 byte real
 634
 635    ! argument
 636    real(rk), intent(in) :: X
 637
 638    ! return value
 639    real(rk) :: ERFCX_R4
 640
 641    ! local variables
 642    INTEGER :: JINT = 2
 643    !------------------------------------------------------------------
 644
 645    CALL CALERF_r4(X, ERFCX_R4, JINT)
 646 END FUNCTION ERFCX_R4
 647
 648 !------------------------------------------------------------------------------------------
 649
 650 end module error_function