wrfv2_fire/phys/module_cam_error_function.F

   1 !------------------------------------------------------------------------
   2 ! Based on error_function.F90 from CAM
   3 ! Ported to WRF by William.Gustafson@pnl.gov, Dec. 2009
   4 !------------------------------------------------------------------------
   5
   6 module 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
  11 implicit none
  12 private
  13 save
  14
  15 ! Public functions
  16 public :: erf, erfc, erfcx
  17
  18 interface erf
  19    module procedure erf_r4
  20    module procedure derf
  21 end interface
  22
  23 interface erfc
  24    module procedure erfc_r4
  25    module procedure derfc
  26 end interface
  27
  28 interface erfcx
  29    module procedure erfcx_r4
  30    module procedure derfcx
  31 end interface
  32
  33 ! Private variables
  34 integer, parameter :: r4 = selected_real_kind(6)  ! 4 byte real
  35 integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real
  36
  37 contains
  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
 167 SUBROUTINE 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
 310 30 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
 331 80 continue
 332 end SUBROUTINE CALERF_r8
 333
 334 !------------------------------------------------------------------------------------------
 335
 336 SUBROUTINE 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
 469 30 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
 490 80 continue
 491 end SUBROUTINE CALERF_r4
 492
 493 !------------------------------------------------------------------------------------------
 494
 495 FUNCTION 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)
 517 END FUNCTION DERF
 518
 519 !------------------------------------------------------------------------------------------
 520
 521 FUNCTION 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)
 543 END FUNCTION ERF_r4
 544
 545 !------------------------------------------------------------------------------------------
 546
 547 FUNCTION 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)
 569 END FUNCTION DERFC
 570
 571 !------------------------------------------------------------------------------------------
 572
 573 FUNCTION 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)
 595 END FUNCTION ERFC_r4
 596
 597 !------------------------------------------------------------------------------------------
 598
 599 FUNCTION 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)
 621 END FUNCTION DERFCX
 622
 623 !------------------------------------------------------------------------------------------
 624
 625 FUNCTION 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)
 647 END FUNCTION ERFCX_R4
 648
 649 !------------------------------------------------------------------------------------------
 650
 651 end module error_function