1 2! KGEN-generated Fortran source file 3! 4! Filename : shr_spfn_mod.F90 5! Generated at: 2015-03-31 09:44:41 6! KGEN version: 0.4.5 7 8 9 10 MODULE shr_spfn_mod 11 ! Module for common mathematical functions 12 ! This #ifdef is to allow the module to be compiled with no dependencies, 13 ! even on shr_kind_mod. 14 USE shr_kind_mod, ONLY: r8 => shr_kind_r8 15 USE shr_const_mod, ONLY: pi => shr_const_pi 16 IMPLICIT NONE 17 PRIVATE 18 ! Error functions 19 20 21 22 ! Gamma functions 23 ! Note that we lack an implementation of log_gamma, but we do have an 24 ! implementation of the upper incomplete gamma function, which is not in 25 ! Fortran 2008. 26 ! Note also that this gamma function is only for double precision. We 27 ! haven't needed an r4 version yet. 28 PUBLIC shr_spfn_gamma 29 30 INTERFACE shr_spfn_gamma 31 MODULE PROCEDURE shr_spfn_gamma_r8 32 END INTERFACE shr_spfn_gamma 33 ! Mathematical constants 34 ! sqrt(pi) 35 ! Define machine-specific constants needed in this module. 36 ! These were used by the original gamma and calerf functions to guarantee 37 ! safety against overflow, and precision, on many different machines. 38 ! By defining the constants in this way, we assume that 1/xmin is 39 ! representable (i.e. does not overflow the real type). This assumption was 40 ! not in the original code, but is valid for IEEE single and double 41 ! precision. 42 ! Double precision 43 !--------------------------------------------------------------------- 44 ! Machine epsilon 45 REAL(KIND=r8), parameter :: epsr8 = epsilon(1._r8) 46 ! "Huge" value is returned when actual value would be infinite. 47 REAL(KIND=r8), parameter :: xinfr8 = huge(1._r8) 48 ! Smallest normal value. 49 REAL(KIND=r8), parameter :: xminr8 = tiny(1._r8) 50 ! Largest number that, when added to 1., yields 1. 51 ! Largest argument for which erfcx > 0. 52 ! Single precision 53 !--------------------------------------------------------------------- 54 ! Machine epsilon 55 ! "Huge" value is returned when actual value would be infinite. 56 ! Smallest normal value. 57 ! Largest number that, when added to 1., yields 1. 58 ! Largest argument for which erfcx > 0. 59 ! For gamma/igamma 60 ! Approximate value of largest acceptable argument to gamma, 61 ! for IEEE double-precision. 62 REAL(KIND=r8), parameter :: xbig_gamma = 171.624_r8 63 CONTAINS 64 65 ! write subroutines 66 ! No subroutines 67 ! No module extern variables 68 ! Wrapper functions for erf 69 70 71 72 73 74 75 ! Wrapper functions for erfc 76 77 78 79 80 81 82 ! Wrapper functions for erfc_scaled 83 84 85 86 elemental FUNCTION shr_spfn_gamma_r8(x) RESULT ( res ) 87 REAL(KIND=r8), intent(in) :: x 88 REAL(KIND=r8) :: res 89 ! No intrinsic 90 res = shr_spfn_gamma_nonintrinsic_r8(x) 91 END FUNCTION shr_spfn_gamma_r8 92 !------------------------------------------------------------------ 93 ! 94 ! 6 December 2006 -- B. Eaton 95 ! The following comments are from the original version of CALERF. 96 ! The only changes in implementing this module are that the function 97 ! names previously used for the single precision versions have been 98 ! adopted for the new generic interfaces. To support these interfaces 99 ! there is now both a single precision version (calerf_r4) and a 100 ! double precision version (calerf_r8) of CALERF below. These versions 101 ! are hardcoded to use IEEE arithmetic. 102 ! 103 !------------------------------------------------------------------ 104 ! 105 ! This packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x) 106 ! for a real argument x. It contains three FUNCTION type 107 ! subprograms: ERF, ERFC, and ERFCX (or ERF_R8, ERFC_R8, and ERFCX_R8), 108 ! and one SUBROUTINE type subprogram, CALERF. The calling 109 ! statements for the primary entries are: 110 ! 111 ! Y=ERF(X) (or Y=ERF_R8(X)), 112 ! 113 ! Y=ERFC(X) (or Y=ERFC_R8(X)), 114 ! and 115 ! Y=ERFCX(X) (or Y=ERFCX_R8(X)). 116 ! 117 ! The routine CALERF is intended for internal packet use only, 118 ! all computations within the packet being concentrated in this 119 ! routine. The function subprograms invoke CALERF with the 120 ! statement 121 ! 122 ! CALL CALERF(ARG,RESULT,JINT) 123 ! 124 ! where the parameter usage is as follows 125 ! 126 ! Function Parameters for CALERF 127 ! call ARG Result JINT 128 ! 129 ! ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0 130 ! ERFC(ARG) ABS(ARG) .LT. XBIG ERFC(ARG) 1 131 ! ERFCX(ARG) XNEG .LT. ARG .LT. XMAX ERFCX(ARG) 2 132 ! 133 ! The main computation evaluates near-minimax approximations 134 ! from "Rational Chebyshev approximations for the error function" 135 ! by W. J. Cody, Math. Comp., 1969, PP. 631-638. This 136 ! transportable program uses rational functions that theoretically 137 ! approximate erf(x) and erfc(x) to at least 18 significant 138 ! decimal digits. The accuracy achieved depends on the arithmetic 139 ! system, the compiler, the intrinsic functions, and proper 140 ! selection of the machine-dependent constants. 141 ! 142 !******************************************************************* 143 !******************************************************************* 144 ! 145 ! Explanation of machine-dependent constants 146 ! 147 ! XMIN = the smallest positive floating-point number. 148 ! XINF = the largest positive finite floating-point number. 149 ! XNEG = the largest negative argument acceptable to ERFCX; 150 ! the negative of the solution to the equation 151 ! 2*exp(x*x) = XINF. 152 ! XSMALL = argument below which erf(x) may be represented by 153 ! 2*x/sqrt(pi) and above which x*x will not underflow. 154 ! A conservative value is the largest machine number X 155 ! such that 1.0 + X = 1.0 to machine precision. 156 ! XBIG = largest argument acceptable to ERFC; solution to 157 ! the equation: W(x) * (1-0.5/x**2) = XMIN, where 158 ! W(x) = exp(-x*x)/[x*sqrt(pi)]. 159 ! XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to 160 ! machine precision. A conservative value is 161 ! 1/[2*sqrt(XSMALL)] 162 ! XMAX = largest acceptable argument to ERFCX; the minimum 163 ! of XINF and 1/[sqrt(pi)*XMIN]. 164 ! 165 ! Approximate values for some important machines are: 166 ! 167 ! XMIN XINF XNEG XSMALL 168 ! 169 ! CDC 7600 (S.P.) 3.13E-294 1.26E+322 -27.220 7.11E-15 170 ! CRAY-1 (S.P.) 4.58E-2467 5.45E+2465 -75.345 7.11E-15 171 ! IEEE (IBM/XT, 172 ! SUN, etc.) (S.P.) 1.18E-38 3.40E+38 -9.382 5.96E-8 173 ! IEEE (IBM/XT, 174 ! SUN, etc.) (D.P.) 2.23D-308 1.79D+308 -26.628 1.11D-16 175 ! IBM 195 (D.P.) 5.40D-79 7.23E+75 -13.190 1.39D-17 176 ! UNIVAC 1108 (D.P.) 2.78D-309 8.98D+307 -26.615 1.73D-18 177 ! VAX D-Format (D.P.) 2.94D-39 1.70D+38 -9.345 1.39D-17 178 ! VAX G-Format (D.P.) 5.56D-309 8.98D+307 -26.615 1.11D-16 179 ! 180 ! 181 ! XBIG XHUGE XMAX 182 ! 183 ! CDC 7600 (S.P.) 25.922 8.39E+6 1.80X+293 184 ! CRAY-1 (S.P.) 75.326 8.39E+6 5.45E+2465 185 ! IEEE (IBM/XT, 186 ! SUN, etc.) (S.P.) 9.194 2.90E+3 4.79E+37 187 ! IEEE (IBM/XT, 188 ! SUN, etc.) (D.P.) 26.543 6.71D+7 2.53D+307 189 ! IBM 195 (D.P.) 13.306 1.90D+8 7.23E+75 190 ! UNIVAC 1108 (D.P.) 26.582 5.37D+8 8.98D+307 191 ! VAX D-Format (D.P.) 9.269 1.90D+8 1.70D+38 192 ! VAX G-Format (D.P.) 26.569 6.71D+7 8.98D+307 193 ! 194 !******************************************************************* 195 !******************************************************************* 196 ! 197 ! Error returns 198 ! 199 ! The program returns ERFC = 0 for ARG .GE. XBIG; 200 ! 201 ! ERFCX = XINF for ARG .LT. XNEG; 202 ! and 203 ! ERFCX = 0 for ARG .GE. XMAX. 204 ! 205 ! 206 ! Intrinsic functions required are: 207 ! 208 ! ABS, AINT, EXP 209 ! 210 ! 211 ! Author: W. J. Cody 212 ! Mathematics and Computer Science Division 213 ! Argonne National Laboratory 214 ! Argonne, IL 60439 215 ! 216 ! Latest modification: March 19, 1990 217 ! 218 !------------------------------------------------------------------ 219 220 !------------------------------------------------------------------------------------------ 221 222 !------------------------------------------------------------------------------------------ 223 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 224 225 pure FUNCTION shr_spfn_gamma_nonintrinsic_r8(x) RESULT ( gamma ) 226 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 227 ! 228 ! 7 Feb 2013 -- S. Santos 229 ! The following comments are from the original version. Changes have 230 ! been made to update syntax and allow inclusion into this module. 231 ! 232 !---------------------------------------------------------------------- 233 ! 234 ! THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X. 235 ! COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1. 236 ! THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA 237 ! FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS 238 ! FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED. 239 ! THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2. 240 ! THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE 241 ! COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE 242 ! MACHINE-DEPENDENT CONSTANTS. 243 ! 244 ! 245 !******************************************************************* 246 !******************************************************************* 247 ! 248 ! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS 249 ! 250 ! BETA - RADIX FOR THE FLOATING-POINT REPRESENTATION 251 ! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS 252 ! XBIG - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE 253 ! IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION 254 ! GAMMA(XBIG) = BETA**MAXEXP 255 ! XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER; 256 ! APPROXIMATELY BETA**MAXEXP 257 ! EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT 258 ! 1.0+EPS .GT. 1.0 259 ! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT 260 ! 1/XMININ IS MACHINE REPRESENTABLE 261 ! 262 ! APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE: 263 ! 264 ! BETA MAXEXP XBIG 265 ! 266 ! CRAY-1 (S.P.) 2 8191 966.961 267 ! CYBER 180/855 268 ! UNDER NOS (S.P.) 2 1070 177.803 269 ! IEEE (IBM/XT, 270 ! SUN, ETC.) (S.P.) 2 128 35.040 271 ! IEEE (IBM/XT, 272 ! SUN, ETC.) (D.P.) 2 1024 171.624 273 ! IBM 3033 (D.P.) 16 63 57.574 274 ! VAX D-FORMAT (D.P.) 2 127 34.844 275 ! VAX G-FORMAT (D.P.) 2 1023 171.489 276 ! 277 ! XINF EPS XMININ 278 ! 279 ! CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466 280 ! CYBER 180/855 281 ! UNDER NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294 282 ! IEEE (IBM/XT, 283 ! SUN, ETC.) (S.P.) 3.40E+38 1.19E-7 1.18E-38 284 ! IEEE (IBM/XT, 285 ! SUN, ETC.) (D.P.) 1.79D+308 2.22D-16 2.23D-308 286 ! IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76 287 ! VAX D-FORMAT (D.P.) 1.70D+38 1.39D-17 5.88D-39 288 ! VAX G-FORMAT (D.P.) 8.98D+307 1.11D-16 1.12D-308 289 ! 290 !******************************************************************* 291 !******************************************************************* 292 ! 293 ! ERROR RETURNS 294 ! 295 ! THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR 296 ! WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED 297 ! TO BE FREE OF UNDERFLOW AND OVERFLOW. 298 ! 299 ! 300 ! INTRINSIC FUNCTIONS REQUIRED ARE: 301 ! 302 ! INT, DBLE, EXP, LOG, REAL, SIN 303 ! 304 ! 305 ! REFERENCES: AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL 306 ! FUNCTIONS W. J. CODY, LECTURE NOTES IN MATHEMATICS, 307 ! 506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON 308 ! (ED.), SPRINGER VERLAG, BERLIN, 1976. 309 ! 310 ! COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND 311 ! SONS, NEW YORK, 1968. 312 ! 313 ! LATEST MODIFICATION: OCTOBER 12, 1989 314 ! 315 ! AUTHORS: W. J. CODY AND L. STOLTZ 316 ! APPLIED MATHEMATICS DIVISION 317 ! ARGONNE NATIONAL LABORATORY 318 ! ARGONNE, IL 60439 319 ! 320 !---------------------------------------------------------------------- 321 REAL(KIND=r8), intent(in) :: x 322 REAL(KIND=r8) :: gamma 323 REAL(KIND=r8) :: fact 324 REAL(KIND=r8) :: sum 325 REAL(KIND=r8) :: y 326 REAL(KIND=r8) :: y1 327 REAL(KIND=r8) :: res 328 REAL(KIND=r8) :: z 329 REAL(KIND=r8) :: xnum 330 REAL(KIND=r8) :: xden 331 REAL(KIND=r8) :: ysq 332 INTEGER :: n 333 INTEGER :: i 334 LOGICAL :: negative_odd 335 ! log(2*pi)/2 336 REAL(KIND=r8), parameter :: logsqrt2pi = 0.9189385332046727417803297e0_r8 337 !---------------------------------------------------------------------- 338 ! NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX 339 ! APPROXIMATION OVER (1,2). 340 !---------------------------------------------------------------------- 341 REAL(KIND=r8), parameter :: p(8) = (/-1.71618513886549492533811e+0_r8, 2.47656508055759199108314e+1_r8, & 342 -3.79804256470945635097577e+2_r8, 6.29331155312818442661052e+2_r8, 8.66966202790413211295064e+2_r8,& 343 -3.14512729688483675254357e+4_r8, -3.61444134186911729807069e+4_r8, 6.64561438202405440627855e+4_r8 /) 344 REAL(KIND=r8), parameter :: q(8) = (/-3.08402300119738975254353e+1_r8, 3.15350626979604161529144e+2_r8, & 345 -1.01515636749021914166146e+3_r8,-3.10777167157231109440444e+3_r8, 2.25381184209801510330112e+4_r8, & 346 4.75584627752788110767815e+3_r8, -1.34659959864969306392456e+5_r8,-1.15132259675553483497211e+5_r8 /) 347 !---------------------------------------------------------------------- 348 ! COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF). 349 !---------------------------------------------------------------------- 350 REAL(KIND=r8), parameter :: c(7) = (/-1.910444077728e-03_r8, 8.4171387781295e-04_r8, & 351 -5.952379913043012e-04_r8, 7.93650793500350248e-04_r8, -2.777777777777681622553e-03_r8, & 352 8.333333333333333331554247e-02_r8, 5.7083835261e-03_r8 /) 353 negative_odd = .false. 354 fact = 1._r8 355 n = 0 356 y = x 357 IF (y <= 0._r8) THEN 358 !---------------------------------------------------------------------- 359 ! ARGUMENT IS NEGATIVE 360 !---------------------------------------------------------------------- 361 y = -x 362 y1 = aint(y) 363 res = y - y1 364 IF (res /= 0._r8) THEN 365 negative_odd = (y1 /= aint(y1*0.5_r8)*2._r8) 366 fact = -pi/sin(pi*res) 367 y = y + 1._r8 368 ELSE 369 gamma = xinfr8 370 RETURN 371 END IF 372 END IF 373 !---------------------------------------------------------------------- 374 ! ARGUMENT IS POSITIVE 375 !---------------------------------------------------------------------- 376 IF (y < epsr8) THEN 377 !---------------------------------------------------------------------- 378 ! ARGUMENT .LT. EPS 379 !---------------------------------------------------------------------- 380 IF (y >= xminr8) THEN 381 res = 1._r8/y 382 ELSE 383 gamma = xinfr8 384 RETURN 385 END IF 386 ELSE IF (y < 12._r8) THEN 387 y1 = y 388 IF (y < 1._r8) THEN 389 !---------------------------------------------------------------------- 390 ! 0.0 .LT. ARGUMENT .LT. 1.0 391 !---------------------------------------------------------------------- 392 z = y 393 y = y + 1._r8 394 ELSE 395 !---------------------------------------------------------------------- 396 ! 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY 397 !---------------------------------------------------------------------- 398 n = int(y) - 1 399 y = y - real(n, r8) 400 z = y - 1._r8 401 END IF 402 !---------------------------------------------------------------------- 403 ! EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0 404 !---------------------------------------------------------------------- 405 xnum = 0._r8 406 xden = 1._r8 407 DO i=1,8 408 xnum = (xnum+p(i))*z 409 xden = xden*z + q(i) 410 END DO 411 res = xnum/xden + 1._r8 412 IF (y1 < y) THEN 413 !---------------------------------------------------------------------- 414 ! ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0 415 !---------------------------------------------------------------------- 416 res = res/y1 417 ELSE IF (y1 > y) THEN 418 !---------------------------------------------------------------------- 419 ! ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0 420 !---------------------------------------------------------------------- 421 DO i = 1,n 422 res = res*y 423 y = y + 1._r8 424 END DO 425 END IF 426 ELSE 427 !---------------------------------------------------------------------- 428 ! EVALUATE FOR ARGUMENT .GE. 12.0, 429 !---------------------------------------------------------------------- 430 IF (y <= xbig_gamma) THEN 431 ysq = y*y 432 sum = c(7) 433 DO i=1,6 434 sum = sum/ysq + c(i) 435 END DO 436 sum = sum/y - y + logsqrt2pi 437 sum = sum + (y-0.5_r8)*log(y) 438 res = exp(sum) 439 ELSE 440 gamma = xinfr8 441 RETURN 442 END IF 443 END IF 444 !---------------------------------------------------------------------- 445 ! FINAL ADJUSTMENTS AND RETURN 446 !---------------------------------------------------------------------- 447 IF (negative_odd) res = -res 448 IF (fact /= 1._r8) res = fact/res 449 gamma = res 450 ! ---------- LAST LINE OF GAMMA ---------- 451 END FUNCTION shr_spfn_gamma_nonintrinsic_r8 452 !! Incomplete Gamma function 453 !! 454 !! @author Tianyi Fan 455 !! @version August-2010 456 457 END MODULE shr_spfn_mod 458