! Define flags for compilers supporting Fortran 2008 intrinsics ! HAVE_GAMMA_INTRINSICS: gamma and log_gamma ! HAVE_ERF_INTRINSICS: erf, erfc, and erfc_scaled ! erfc_scaled(x) = (exp(x**2)*erfc(x)) ! Use this flag for compilers that don't have real intrinsics, but link in ! a library for you. ! HAVE_ERF_EXTERNALS: erf and erfc ! These compilers have the intrinsics. ! Intel also has them (and Cray), but as of mid-2015, our implementation is ! actually faster, in part because they do not properly vectorize, so we ! pretend that the compiler version doesn't exist. #if defined CPRIBM || defined __GFORTRAN__ #define HAVE_GAMMA_INTRINSICS #define HAVE_ERF_INTRINSICS #endif ! PGI has external erf/derf and erfc/derfc, and will link them for you, but ! it does not consider them "intrinsics" right now. #if defined CPRPGI #define HAVE_ERF_EXTERNALS #endif ! As of 5.3.1, NAG does not have any of these. module shr_spfn_mod ! Module for common mathematical functions ! This #ifdef is to allow the module to be compiled with no dependencies, ! even on shr_kind_mod. #ifndef NO_CSM_SHARE use shr_kind_mod, only: & r4 => shr_kind_r4, & r8 => shr_kind_r8 use shr_const_mod, only: & pi => shr_const_pi #endif implicit none private save #ifdef NO_CSM_SHARE integer, parameter :: r4 = selected_real_kind(6) ! 4 byte real integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real real(r8), parameter :: pi = 3.1415926535897932384626434E0_r8 #endif ! Error functions public :: shr_spfn_erf public :: shr_spfn_erfc public :: shr_spfn_erfc_scaled interface shr_spfn_erf module procedure shr_spfn_erf_r4 module procedure shr_spfn_erf_r8 end interface shr_spfn_erf interface shr_spfn_erfc module procedure shr_spfn_erfc_r4 module procedure shr_spfn_erfc_r8 end interface shr_spfn_erfc interface shr_spfn_erfc_scaled module procedure shr_spfn_erfc_scaled_r4 module procedure shr_spfn_erfc_scaled_r8 end interface shr_spfn_erfc_scaled ! Gamma functions ! Note that we lack an implementation of log_gamma, but we do have an ! implementation of the upper incomplete gamma function, which is not in ! Fortran 2008. ! Note also that this gamma function is only for double precision. We ! haven't needed an r4 version yet. public :: shr_spfn_gamma public :: shr_spfn_igamma interface shr_spfn_gamma module procedure shr_spfn_gamma_r8 end interface shr_spfn_gamma ! Mathematical constants ! sqrt(pi) real(r8), parameter :: sqrtpi = 1.77245385090551602729_r8 ! Define machine-specific constants needed in this module. ! These were used by the original gamma and calerf functions to guarantee ! safety against overflow, and precision, on many different machines. ! By defining the constants in this way, we assume that 1/xmin is ! representable (i.e. does not overflow the real type). This assumption was ! not in the original code, but is valid for IEEE single and double ! precision. ! Double precision !--------------------------------------------------------------------- ! Machine epsilon real(r8), parameter :: epsr8 = epsilon(1._r8) ! "Huge" value is returned when actual value would be infinite. real(r8), parameter :: xinfr8 = huge(1._r8) ! Smallest normal value. real(r8), parameter :: xminr8 = tiny(1._r8) ! Largest number that, when added to 1., yields 1. real(r8), parameter :: xsmallr8 = epsr8/2._r8 ! Largest argument for which erfcx > 0. real(r8), parameter :: xmaxr8 = 1._r8/(sqrtpi*xminr8) ! Single precision !--------------------------------------------------------------------- ! Machine epsilon real(r4), parameter :: epsr4 = epsilon(1._r4) ! "Huge" value is returned when actual value would be infinite. real(r4), parameter :: xinfr4 = huge(1._r4) ! Smallest normal value. real(r4), parameter :: xminr4 = tiny(1._r4) ! Largest number that, when added to 1., yields 1. real(r4), parameter :: xsmallr4 = epsr4/2._r4 ! Largest argument for which erfcx > 0. real(r4), parameter :: xmaxr4 = 1._r4/(real(sqrtpi,r4)*xminr4) ! For gamma/igamma ! Approximate value of largest acceptable argument to gamma, ! for IEEE double-precision. real(r8), parameter :: xbig_gamma = 171.624_r8 contains ! Wrapper functions for erf function shr_spfn_erf_r4(x) result(res) real(r4), intent(in) :: x real(r4) :: res #ifdef HAVE_ERF_EXTERNALS ! If erf is provided as an external, provide ! explicit interface here. interface function erf(x) import :: r4 real(r4) :: x, erf end function erf end interface #endif #ifdef HAVE_ERF_INTRINSICS ! Call intrinsic erf. intrinsic erf res = erf(x) #else #ifdef HAVE_ERF_EXTERNALS ! Call compiler-provided external erf. res = erf(x) #else ! No compiler-provided erf, so call local version. call calerf_r4(x, res, 0) #endif #endif end function shr_spfn_erf_r4 function shr_spfn_erf_r8(x) result(res) real(r8), intent(in) :: x real(r8) :: res #ifdef HAVE_ERF_EXTERNALS ! If erf is provided as an external, provide ! explicit interface here. interface function derf(x) import :: r8 real(r8) :: x, derf end function derf end interface #endif #ifdef HAVE_ERF_INTRINSICS ! Call intrinsic erf. intrinsic erf res = erf(x) #else #ifdef HAVE_ERF_EXTERNALS ! Call compiler-provided external erf. res = derf(x) #else ! No compiler-provided erf, so call local version. call calerf_r8(x, res, 0) #endif #endif end function shr_spfn_erf_r8 ! Wrapper functions for erfc function shr_spfn_erfc_r4(x) result(res) real(r4), intent(in) :: x real(r4) :: res #ifdef HAVE_ERF_EXTERNALS ! If erfc is provided as an external, provide ! explicit interface here. interface function erfc(x) import :: r4 real(r4) :: x, erfc end function erfc end interface #endif #ifdef HAVE_ERF_INTRINSICS ! Call intrinsic erfc. intrinsic erfc res = erfc(x) #else #ifdef HAVE_ERF_EXTERNALS ! Call compiler-provided external erfc. res = erfc(x) #else ! No compiler-provided erfc, so call local version. call calerf_r4(x, res, 1) #endif #endif end function shr_spfn_erfc_r4 function shr_spfn_erfc_r8(x) result(res) real(r8), intent(in) :: x real(r8) :: res #ifdef HAVE_ERF_EXTERNALS ! If erfc is provided as an external, provide ! explicit interface here. interface function derfc(x) import :: r8 real(r8) :: x, derfc end function derfc end interface #endif #ifdef HAVE_ERF_INTRINSICS ! Call intrinsic erfc. intrinsic erfc res = erfc(x) #else #ifdef HAVE_ERF_EXTERNALS ! Call compiler-provided external erfc. res = derfc(x) #else ! No compiler-provided erfc, so call local version. call calerf_r8(x, res, 1) #endif #endif end function shr_spfn_erfc_r8 ! Wrapper functions for erfc_scaled function shr_spfn_erfc_scaled_r4(x) result(res) real(r4), intent(in) :: x real(r4) :: res #if defined HAVE_ERF_INTRINSICS ! Call intrinsic erfc_scaled. intrinsic erfc_scaled res = erfc_scaled(x) #else ! No intrinsic. call calerf_r4(x, res, 2) #endif end function shr_spfn_erfc_scaled_r4 function shr_spfn_erfc_scaled_r8(x) result(res) real(r8), intent(in) :: x real(r8) :: res #if defined HAVE_ERF_INTRINSICS ! Call intrinsic erfc_scaled. intrinsic erfc_scaled res = erfc_scaled(x) #else ! No intrinsic. call calerf_r8(x, res, 2) #endif end function shr_spfn_erfc_scaled_r8 elemental function shr_spfn_gamma_r8(x) result(res) real(r8), intent(in) :: x real(r8) :: res #if defined HAVE_GAMMA_INTRINSICS ! Call intrinsic gamma. intrinsic gamma res = gamma(x) #else ! No intrinsic res = shr_spfn_gamma_nonintrinsic_r8(x) #endif end function shr_spfn_gamma_r8 !------------------------------------------------------------------ ! ! 6 December 2006 -- B. Eaton ! The following comments are from the original version of CALERF. ! The only changes in implementing this module are that the function ! names previously used for the single precision versions have been ! adopted for the new generic interfaces. To support these interfaces ! there is now both a single precision version (calerf_r4) and a ! double precision version (calerf_r8) of CALERF below. These versions ! are hardcoded to use IEEE arithmetic. ! !------------------------------------------------------------------ ! ! This packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x) ! for a real argument x. It contains three FUNCTION type ! subprograms: ERF, ERFC, and ERFCX (or ERF_R8, ERFC_R8, and ERFCX_R8), ! and one SUBROUTINE type subprogram, CALERF. The calling ! statements for the primary entries are: ! ! Y=ERF(X) (or Y=ERF_R8(X)), ! ! Y=ERFC(X) (or Y=ERFC_R8(X)), ! and ! Y=ERFCX(X) (or Y=ERFCX_R8(X)). ! ! The routine CALERF is intended for internal packet use only, ! all computations within the packet being concentrated in this ! routine. The function subprograms invoke CALERF with the ! statement ! ! CALL CALERF(ARG,RESULT,JINT) ! ! where the parameter usage is as follows ! ! Function Parameters for CALERF ! call ARG Result JINT ! ! ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0 ! ERFC(ARG) ABS(ARG) .LT. XBIG ERFC(ARG) 1 ! ERFCX(ARG) XNEG .LT. ARG .LT. XMAX ERFCX(ARG) 2 ! ! The main computation evaluates near-minimax approximations ! from "Rational Chebyshev approximations for the error function" ! by W. J. Cody, Math. Comp., 1969, PP. 631-638. This ! transportable program uses rational functions that theoretically ! approximate erf(x) and erfc(x) to at least 18 significant ! decimal digits. The accuracy achieved depends on the arithmetic ! system, the compiler, the intrinsic functions, and proper ! selection of the machine-dependent constants. ! !******************************************************************* !******************************************************************* ! ! Explanation of machine-dependent constants ! ! XMIN = the smallest positive floating-point number. ! XINF = the largest positive finite floating-point number. ! XNEG = the largest negative argument acceptable to ERFCX; ! the negative of the solution to the equation ! 2*exp(x*x) = XINF. ! XSMALL = argument below which erf(x) may be represented by ! 2*x/sqrt(pi) and above which x*x will not underflow. ! A conservative value is the largest machine number X ! such that 1.0 + X = 1.0 to machine precision. ! XBIG = largest argument acceptable to ERFC; solution to ! the equation: W(x) * (1-0.5/x**2) = XMIN, where ! W(x) = exp(-x*x)/[x*sqrt(pi)]. ! XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to ! machine precision. A conservative value is ! 1/[2*sqrt(XSMALL)] ! XMAX = largest acceptable argument to ERFCX; the minimum ! of XINF and 1/[sqrt(pi)*XMIN]. ! ! Approximate values for some important machines are: ! ! XMIN XINF XNEG XSMALL ! ! CDC 7600 (S.P.) 3.13E-294 1.26E+322 -27.220 7.11E-15 ! CRAY-1 (S.P.) 4.58E-2467 5.45E+2465 -75.345 7.11E-15 ! IEEE (IBM/XT, ! SUN, etc.) (S.P.) 1.18E-38 3.40E+38 -9.382 5.96E-8 ! IEEE (IBM/XT, ! SUN, etc.) (D.P.) 2.23D-308 1.79D+308 -26.628 1.11D-16 ! IBM 195 (D.P.) 5.40D-79 7.23E+75 -13.190 1.39D-17 ! UNIVAC 1108 (D.P.) 2.78D-309 8.98D+307 -26.615 1.73D-18 ! VAX D-Format (D.P.) 2.94D-39 1.70D+38 -9.345 1.39D-17 ! VAX G-Format (D.P.) 5.56D-309 8.98D+307 -26.615 1.11D-16 ! ! ! XBIG XHUGE XMAX ! ! CDC 7600 (S.P.) 25.922 8.39E+6 1.80X+293 ! CRAY-1 (S.P.) 75.326 8.39E+6 5.45E+2465 ! IEEE (IBM/XT, ! SUN, etc.) (S.P.) 9.194 2.90E+3 4.79E+37 ! IEEE (IBM/XT, ! SUN, etc.) (D.P.) 26.543 6.71D+7 2.53D+307 ! IBM 195 (D.P.) 13.306 1.90D+8 7.23E+75 ! UNIVAC 1108 (D.P.) 26.582 5.37D+8 8.98D+307 ! VAX D-Format (D.P.) 9.269 1.90D+8 1.70D+38 ! VAX G-Format (D.P.) 26.569 6.71D+7 8.98D+307 ! !******************************************************************* !******************************************************************* ! ! Error returns ! ! The program returns ERFC = 0 for ARG .GE. XBIG; ! ! ERFCX = XINF for ARG .LT. XNEG; ! and ! ERFCX = 0 for ARG .GE. XMAX. ! ! ! Intrinsic functions required are: ! ! ABS, AINT, EXP ! ! ! Author: W. J. Cody ! Mathematics and Computer Science Division ! Argonne National Laboratory ! Argonne, IL 60439 ! ! Latest modification: March 19, 1990 ! !------------------------------------------------------------------ SUBROUTINE CALERF_r8(ARG, RESULT, JINT) !------------------------------------------------------------------ ! This version uses 8-byte reals !------------------------------------------------------------------ integer, parameter :: rk = r8 ! arguments real(rk), intent(in) :: arg integer, intent(in) :: jint real(rk), intent(out) :: result ! local variables INTEGER :: I real(rk) :: X, Y, YSQ, XNUM, XDEN, DEL !------------------------------------------------------------------ ! Mathematical constants !------------------------------------------------------------------ real(rk), parameter :: ZERO = 0.0E0_rk real(rk), parameter :: FOUR = 4.0E0_rk real(rk), parameter :: ONE = 1.0E0_rk real(rk), parameter :: HALF = 0.5E0_rk real(rk), parameter :: TWO = 2.0E0_rk ! 1/sqrt(pi) real(rk), parameter :: SQRPI = 5.6418958354775628695E-1_rk real(rk), parameter :: THRESH = 0.46875E0_rk real(rk), parameter :: SIXTEN = 16.0E0_rk !------------------------------------------------------------------ ! Machine-dependent constants: IEEE double precision values !------------------------------------------------------------------ real(rk), parameter :: XNEG = -26.628E0_r8 real(rk), parameter :: XBIG = 26.543E0_r8 real(rk), parameter :: XHUGE = 6.71E7_r8 !------------------------------------------------------------------ ! Coefficients for approximation to erf in first interval !------------------------------------------------------------------ real(rk), parameter :: A(5) = (/ 3.16112374387056560E00_rk, 1.13864154151050156E02_rk, & 3.77485237685302021E02_rk, 3.20937758913846947E03_rk, & 1.85777706184603153E-1_rk /) real(rk), parameter :: B(4) = (/ 2.36012909523441209E01_rk, 2.44024637934444173E02_rk, & 1.28261652607737228E03_rk, 2.84423683343917062E03_rk /) !------------------------------------------------------------------ ! Coefficients for approximation to erfc in second interval !------------------------------------------------------------------ real(rk), parameter :: C(9) = (/ 5.64188496988670089E-1_rk, 8.88314979438837594E00_rk, & 6.61191906371416295E01_rk, 2.98635138197400131E02_rk, & 8.81952221241769090E02_rk, 1.71204761263407058E03_rk, & 2.05107837782607147E03_rk, 1.23033935479799725E03_rk, & 2.15311535474403846E-8_rk /) real(rk), parameter :: D(8) = (/ 1.57449261107098347E01_rk, 1.17693950891312499E02_rk, & 5.37181101862009858E02_rk, 1.62138957456669019E03_rk, & 3.29079923573345963E03_rk, 4.36261909014324716E03_rk, & 3.43936767414372164E03_rk, 1.23033935480374942E03_rk /) !------------------------------------------------------------------ ! Coefficients for approximation to erfc in third interval !------------------------------------------------------------------ real(rk), parameter :: P(6) = (/ 3.05326634961232344E-1_rk, 3.60344899949804439E-1_rk, & 1.25781726111229246E-1_rk, 1.60837851487422766E-2_rk, & 6.58749161529837803E-4_rk, 1.63153871373020978E-2_rk /) real(rk), parameter :: Q(5) = (/ 2.56852019228982242E00_rk, 1.87295284992346047E00_rk, & 5.27905102951428412E-1_rk, 6.05183413124413191E-2_rk, & 2.33520497626869185E-3_rk /) !------------------------------------------------------------------ X = ARG Y = ABS(X) IF (Y .LE. THRESH) THEN !------------------------------------------------------------------ ! Evaluate erf for |X| <= 0.46875 !------------------------------------------------------------------ YSQ = ZERO IF (Y .GT. XSMALLR8) YSQ = Y * Y XNUM = A(5)*YSQ XDEN = YSQ DO I = 1, 3 XNUM = (XNUM + A(I)) * YSQ XDEN = (XDEN + B(I)) * YSQ end do RESULT = X * (XNUM + A(4)) / (XDEN + B(4)) IF (JINT .NE. 0) RESULT = ONE - RESULT IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT GO TO 80 ELSE IF (Y .LE. FOUR) THEN !------------------------------------------------------------------ ! Evaluate erfc for 0.46875 <= |X| <= 4.0 !------------------------------------------------------------------ XNUM = C(9)*Y XDEN = Y DO I = 1, 7 XNUM = (XNUM + C(I)) * Y XDEN = (XDEN + D(I)) * Y end do RESULT = (XNUM + C(8)) / (XDEN + D(8)) IF (JINT .NE. 2) THEN YSQ = AINT(Y*SIXTEN)/SIXTEN DEL = (Y-YSQ)*(Y+YSQ) RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT END IF ELSE !------------------------------------------------------------------ ! Evaluate erfc for |X| > 4.0 !------------------------------------------------------------------ RESULT = ZERO IF (Y .GE. XBIG) THEN IF ((JINT .NE. 2) .OR. (Y .GE. XMAXR8)) GO TO 30 IF (Y .GE. XHUGE) THEN RESULT = SQRPI / Y GO TO 30 END IF END IF YSQ = ONE / (Y * Y) XNUM = P(6)*YSQ XDEN = YSQ DO I = 1, 4 XNUM = (XNUM + P(I)) * YSQ XDEN = (XDEN + Q(I)) * YSQ end do RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5)) RESULT = (SQRPI - RESULT) / Y IF (JINT .NE. 2) THEN YSQ = AINT(Y*SIXTEN)/SIXTEN DEL = (Y-YSQ)*(Y+YSQ) RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT END IF END IF 30 continue !------------------------------------------------------------------ ! Fix up for negative argument, erf, etc. !------------------------------------------------------------------ IF (JINT .EQ. 0) THEN RESULT = (HALF - RESULT) + HALF IF (X .LT. ZERO) RESULT = -RESULT ELSE IF (JINT .EQ. 1) THEN IF (X .LT. ZERO) RESULT = TWO - RESULT ELSE IF (X .LT. ZERO) THEN IF (X .LT. XNEG) THEN RESULT = XINFR8 ELSE YSQ = AINT(X*SIXTEN)/SIXTEN DEL = (X-YSQ)*(X+YSQ) Y = EXP(YSQ*YSQ) * EXP(DEL) RESULT = (Y+Y) - RESULT END IF END IF END IF 80 continue end SUBROUTINE CALERF_r8 !------------------------------------------------------------------------------------------ SUBROUTINE CALERF_r4(ARG, RESULT, JINT) !------------------------------------------------------------------ ! This version uses 4-byte reals !------------------------------------------------------------------ integer, parameter :: rk = r4 ! arguments real(rk), intent(in) :: arg integer, intent(in) :: jint real(rk), intent(out) :: result ! local variables INTEGER :: I real(rk) :: X, Y, YSQ, XNUM, XDEN, DEL !------------------------------------------------------------------ ! Mathematical constants !------------------------------------------------------------------ real(rk), parameter :: ZERO = 0.0E0_rk real(rk), parameter :: FOUR = 4.0E0_rk real(rk), parameter :: ONE = 1.0E0_rk real(rk), parameter :: HALF = 0.5E0_rk real(rk), parameter :: TWO = 2.0E0_rk ! 1/sqrt(pi) real(rk), parameter :: SQRPI = 5.6418958354775628695E-1_rk real(rk), parameter :: THRESH = 0.46875E0_rk real(rk), parameter :: SIXTEN = 16.0E0_rk !------------------------------------------------------------------ ! Machine-dependent constants: IEEE single precision values !------------------------------------------------------------------ real(rk), parameter :: XNEG = -9.382E0_r4 real(rk), parameter :: XBIG = 9.194E0_r4 real(rk), parameter :: XHUGE = 2.90E3_r4 !------------------------------------------------------------------ ! Coefficients for approximation to erf in first interval !------------------------------------------------------------------ real(rk), parameter :: A(5) = (/ 3.16112374387056560E00_rk, 1.13864154151050156E02_rk, & 3.77485237685302021E02_rk, 3.20937758913846947E03_rk, & 1.85777706184603153E-1_rk /) real(rk), parameter :: B(4) = (/ 2.36012909523441209E01_rk, 2.44024637934444173E02_rk, & 1.28261652607737228E03_rk, 2.84423683343917062E03_rk /) !------------------------------------------------------------------ ! Coefficients for approximation to erfc in second interval !------------------------------------------------------------------ real(rk), parameter :: C(9) = (/ 5.64188496988670089E-1_rk, 8.88314979438837594E00_rk, & 6.61191906371416295E01_rk, 2.98635138197400131E02_rk, & 8.81952221241769090E02_rk, 1.71204761263407058E03_rk, & 2.05107837782607147E03_rk, 1.23033935479799725E03_rk, & 2.15311535474403846E-8_rk /) real(rk), parameter :: D(8) = (/ 1.57449261107098347E01_rk, 1.17693950891312499E02_rk, & 5.37181101862009858E02_rk, 1.62138957456669019E03_rk, & 3.29079923573345963E03_rk, 4.36261909014324716E03_rk, & 3.43936767414372164E03_rk, 1.23033935480374942E03_rk /) !------------------------------------------------------------------ ! Coefficients for approximation to erfc in third interval !------------------------------------------------------------------ real(rk), parameter :: P(6) = (/ 3.05326634961232344E-1_rk, 3.60344899949804439E-1_rk, & 1.25781726111229246E-1_rk, 1.60837851487422766E-2_rk, & 6.58749161529837803E-4_rk, 1.63153871373020978E-2_rk /) real(rk), parameter :: Q(5) = (/ 2.56852019228982242E00_rk, 1.87295284992346047E00_rk, & 5.27905102951428412E-1_rk, 6.05183413124413191E-2_rk, & 2.33520497626869185E-3_rk /) !------------------------------------------------------------------ X = ARG Y = ABS(X) IF (Y .LE. THRESH) THEN !------------------------------------------------------------------ ! Evaluate erf for |X| <= 0.46875 !------------------------------------------------------------------ YSQ = ZERO IF (Y .GT. XSMALLR4) YSQ = Y * Y XNUM = A(5)*YSQ XDEN = YSQ DO I = 1, 3 XNUM = (XNUM + A(I)) * YSQ XDEN = (XDEN + B(I)) * YSQ end do RESULT = X * (XNUM + A(4)) / (XDEN + B(4)) IF (JINT .NE. 0) RESULT = ONE - RESULT IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT GO TO 80 ELSE IF (Y .LE. FOUR) THEN !------------------------------------------------------------------ ! Evaluate erfc for 0.46875 <= |X| <= 4.0 !------------------------------------------------------------------ XNUM = C(9)*Y XDEN = Y DO I = 1, 7 XNUM = (XNUM + C(I)) * Y XDEN = (XDEN + D(I)) * Y end do RESULT = (XNUM + C(8)) / (XDEN + D(8)) IF (JINT .NE. 2) THEN YSQ = AINT(Y*SIXTEN)/SIXTEN DEL = (Y-YSQ)*(Y+YSQ) RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT END IF ELSE !------------------------------------------------------------------ ! Evaluate erfc for |X| > 4.0 !------------------------------------------------------------------ RESULT = ZERO IF (Y .GE. XBIG) THEN IF ((JINT .NE. 2) .OR. (Y .GE. XMAXR4)) GO TO 30 IF (Y .GE. XHUGE) THEN RESULT = SQRPI / Y GO TO 30 END IF END IF YSQ = ONE / (Y * Y) XNUM = P(6)*YSQ XDEN = YSQ DO I = 1, 4 XNUM = (XNUM + P(I)) * YSQ XDEN = (XDEN + Q(I)) * YSQ end do RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5)) RESULT = (SQRPI - RESULT) / Y IF (JINT .NE. 2) THEN YSQ = AINT(Y*SIXTEN)/SIXTEN DEL = (Y-YSQ)*(Y+YSQ) RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT END IF END IF 30 continue !------------------------------------------------------------------ ! Fix up for negative argument, erf, etc. !------------------------------------------------------------------ IF (JINT .EQ. 0) THEN RESULT = (HALF - RESULT) + HALF IF (X .LT. ZERO) RESULT = -RESULT ELSE IF (JINT .EQ. 1) THEN IF (X .LT. ZERO) RESULT = TWO - RESULT ELSE IF (X .LT. ZERO) THEN IF (X .LT. XNEG) THEN RESULT = XINFR4 ELSE YSQ = AINT(X*SIXTEN)/SIXTEN DEL = (X-YSQ)*(X+YSQ) Y = EXP(YSQ*YSQ) * EXP(DEL) RESULT = (Y+Y) - RESULT END IF END IF END IF 80 continue end SUBROUTINE CALERF_r4 !------------------------------------------------------------------------------------------ !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc pure function shr_spfn_gamma_nonintrinsic_r8(X) result(gamma) !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! ! 7 Feb 2013 -- S. Santos ! The following comments are from the original version. Changes have ! been made to update syntax and allow inclusion into this module. ! !---------------------------------------------------------------------- ! ! THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X. ! COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1. ! THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA ! FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS ! FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED. ! THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2. ! THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE ! COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE ! MACHINE-DEPENDENT CONSTANTS. ! ! !******************************************************************* !******************************************************************* ! ! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS ! ! BETA - RADIX FOR THE FLOATING-POINT REPRESENTATION ! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS ! XBIG - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE ! IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION ! GAMMA(XBIG) = BETA**MAXEXP ! XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER; ! APPROXIMATELY BETA**MAXEXP ! EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT ! 1.0+EPS .GT. 1.0 ! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT ! 1/XMININ IS MACHINE REPRESENTABLE ! ! APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE: ! ! BETA MAXEXP XBIG ! ! CRAY-1 (S.P.) 2 8191 966.961 ! CYBER 180/855 ! UNDER NOS (S.P.) 2 1070 177.803 ! IEEE (IBM/XT, ! SUN, ETC.) (S.P.) 2 128 35.040 ! IEEE (IBM/XT, ! SUN, ETC.) (D.P.) 2 1024 171.624 ! IBM 3033 (D.P.) 16 63 57.574 ! VAX D-FORMAT (D.P.) 2 127 34.844 ! VAX G-FORMAT (D.P.) 2 1023 171.489 ! ! XINF EPS XMININ ! ! CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466 ! CYBER 180/855 ! UNDER NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294 ! IEEE (IBM/XT, ! SUN, ETC.) (S.P.) 3.40E+38 1.19E-7 1.18E-38 ! IEEE (IBM/XT, ! SUN, ETC.) (D.P.) 1.79D+308 2.22D-16 2.23D-308 ! IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76 ! VAX D-FORMAT (D.P.) 1.70D+38 1.39D-17 5.88D-39 ! VAX G-FORMAT (D.P.) 8.98D+307 1.11D-16 1.12D-308 ! !******************************************************************* !******************************************************************* ! ! ERROR RETURNS ! ! THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR ! WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED ! TO BE FREE OF UNDERFLOW AND OVERFLOW. ! ! ! INTRINSIC FUNCTIONS REQUIRED ARE: ! ! INT, DBLE, EXP, LOG, REAL, SIN ! ! ! REFERENCES: AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL ! FUNCTIONS W. J. CODY, LECTURE NOTES IN MATHEMATICS, ! 506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON ! (ED.), SPRINGER VERLAG, BERLIN, 1976. ! ! COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND ! SONS, NEW YORK, 1968. ! ! LATEST MODIFICATION: OCTOBER 12, 1989 ! ! AUTHORS: W. J. CODY AND L. STOLTZ ! APPLIED MATHEMATICS DIVISION ! ARGONNE NATIONAL LABORATORY ! ARGONNE, IL 60439 ! !---------------------------------------------------------------------- real(r8), intent(in) :: x real(r8) :: gamma real(r8) :: fact, res, sum, xden, xnum, y, y1, ysq, z integer :: i, n logical :: negative_odd ! log(2*pi)/2 real(r8), parameter :: logsqrt2pi = 0.9189385332046727417803297E0_r8 !---------------------------------------------------------------------- ! NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX ! APPROXIMATION OVER (1,2). !---------------------------------------------------------------------- real(r8), parameter :: P(8) = & (/-1.71618513886549492533811E+0_r8, 2.47656508055759199108314E+1_r8, & -3.79804256470945635097577E+2_r8, 6.29331155312818442661052E+2_r8, & 8.66966202790413211295064E+2_r8,-3.14512729688483675254357E+4_r8, & -3.61444134186911729807069E+4_r8, 6.64561438202405440627855E+4_r8 /) real(r8), parameter :: Q(8) = & (/-3.08402300119738975254353E+1_r8, 3.15350626979604161529144E+2_r8, & -1.01515636749021914166146E+3_r8,-3.10777167157231109440444E+3_r8, & 2.25381184209801510330112E+4_r8, 4.75584627752788110767815E+3_r8, & -1.34659959864969306392456E+5_r8,-1.15132259675553483497211E+5_r8 /) !---------------------------------------------------------------------- ! COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF). !---------------------------------------------------------------------- real(r8), parameter :: C(7) = & (/-1.910444077728E-03_r8, 8.4171387781295E-04_r8, & -5.952379913043012E-04_r8, 7.93650793500350248E-04_r8, & -2.777777777777681622553E-03_r8, 8.333333333333333331554247E-02_r8, & 5.7083835261E-03_r8 /) negative_odd = .false. fact = 1._r8 n = 0 y = x if (y <= 0._r8) then !---------------------------------------------------------------------- ! ARGUMENT IS NEGATIVE !---------------------------------------------------------------------- y = -x y1 = aint(y) res = y - y1 if (res /= 0._r8) then negative_odd = (y1 /= aint(y1*0.5_r8)*2._r8) fact = -pi/sin(pi*res) y = y + 1._r8 else gamma = xinfr8 return end if end if !---------------------------------------------------------------------- ! ARGUMENT IS POSITIVE !---------------------------------------------------------------------- if (y < epsr8) then !---------------------------------------------------------------------- ! ARGUMENT .LT. EPS !---------------------------------------------------------------------- if (y >= xminr8) then res = 1._r8/y else gamma = xinfr8 return end if elseif (y < 12._r8) then y1 = y if (y < 1._r8) then !---------------------------------------------------------------------- ! 0.0 .LT. ARGUMENT .LT. 1.0 !---------------------------------------------------------------------- z = y y = y + 1._r8 else !---------------------------------------------------------------------- ! 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY !---------------------------------------------------------------------- n = int(y) - 1 y = y - real(n, r8) z = y - 1._r8 end if !---------------------------------------------------------------------- ! EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0 !---------------------------------------------------------------------- xnum = 0._r8 xden = 1._r8 do i=1,8 xnum = (xnum+P(i))*z xden = xden*z + Q(i) end do res = xnum/xden + 1._r8 if (y1 < y) then !---------------------------------------------------------------------- ! ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0 !---------------------------------------------------------------------- res = res/y1 elseif (y1 > y) then !---------------------------------------------------------------------- ! ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0 !---------------------------------------------------------------------- do i = 1,n res = res*y y = y + 1._r8 end do end if else !---------------------------------------------------------------------- ! EVALUATE FOR ARGUMENT .GE. 12.0, !---------------------------------------------------------------------- if (y <= xbig_gamma) then ysq = y*y sum = C(7) do i=1,6 sum = sum/ysq + C(i) end do sum = sum/y - y + logsqrt2pi sum = sum + (y-0.5_r8)*log(y) res = exp(sum) else gamma = xinfr8 return end if end if !---------------------------------------------------------------------- ! FINAL ADJUSTMENTS AND RETURN !---------------------------------------------------------------------- if (negative_odd) res = -res if (fact /= 1._r8) res = fact/res gamma = res ! ---------- LAST LINE OF GAMMA ---------- end function shr_spfn_gamma_nonintrinsic_r8 !! Incomplete Gamma function !! !! @author Tianyi Fan !! @version August-2010 real(r8) elemental function shr_spfn_igamma(a, x) ! Upper incomplete gamma function. ! Modified for inclusion in this module and made ! pure elemental, September 2012 real(r8), intent(in) :: a real(r8), intent(in) :: x ! local variable real(r8) :: xam, gin, s, r, t0 integer :: k if (x == 0.0_r8) then shr_spfn_igamma = shr_spfn_gamma(a) return end if xam = -x + a * log(x) if ((xam > 700.0_r8) .or. (a > xbig_gamma)) then ! Out of bounds ! Return "huge" value. shr_spfn_igamma = xinfr8 return else if (x <= (1.0_r8 + a)) then s = 1.0_r8 / a r = s do k = 1,60 r = r * x / (a+k) s = s + r if (abs(r/s) < 1.0e-15_r8) exit end do gin = exp(xam) * s shr_spfn_igamma = shr_spfn_gamma(a) - gin else t0 = 0.0_r8 do k = 60,1,-1 t0 = (k - a) / (1.0_r8 + k / (x + t0)) end do shr_spfn_igamma = exp(xam) / (x + t0) endif end function shr_spfn_igamma end module shr_spfn_mod