include file: ERF_LSP 4
C*LL IN-LINE FUNCTION ERF_LSP------------------------------------------ ERF_LSP.2
!LL ERF_LSP.3
!LL Purpose: ERF_LSP.4
!LL Vectorizable calculation of (x**3 / 3) - x + (PI**0.5 * ERF(x)/2 ERF_LSP.5
!LL needed for LSPFOR routine evaluation of autoconversion rate. ERF_LSP.6
!LL Need to declare REAL ERF_LSP, XERF, TERF and to calculate ERF_LSP.7
!LL SQRT_PI = SQRT(PI) in calling routine. ERF_LSP.8
!LL ERF_LSP.9
!LL Author: A Bushell ERF_LSP.10
!LL ERF_LSP.11
!LL Model Modification history from model version 4.0: ERF_LSP.12
!LL version date ERF_LSP.13
!LL ERF_LSP.14
!LL Documentation: ERF_LSP.15
!LL For abs(XERF) <= 1.161, use the lowest order terms of the integrated ERF_LSP.16
!LL exp(-XERF**2) Maclaurin expansion. For a ERF_LSP.17
!LL mantissa precision error 1E-13, this will be ERF_LSP.18
!LL EXACT for abs(XERF) < 0.03 with max error at ERF_LSP.19
!LL 1.161 of about 0.01%. ERF_LSP.20
!LL For abs(XERF) > 1.161, Press, Teukolsky, Vettering & Flannery ERF_LSP.21
!LL NUMERICAL RECIPES IN FORTRAN eq 6.2.6 gives ERF_LSP.22
!LL a continued fraction representation. The 1st ERF_LSP.23
!LL 4 terms give ERF to high accuracy for ERF_LSP.24
!LL abs(XERF) > 2 & max error < 0.02% at 1.161. ERF_LSP.25
!LL ERF_LSP.26
!LL NOTE: ERF = (2/SQRT(Pi)) * (XERF - ((XERF**3)/3) + ERF_LSP) ERF_LSP.27
!LL BUT I have a provisional ERF COMDECK if you are interested. ACB. ERF_LSP.28
C*---------------------------------------------------------------------- ERF_LSP.29
! ERF_LSP.30
TERF = XERF * XERF ERF_LSP.31
! ERF_LSP.32
IF (ABS(XERF) .LE. 1.161) THEN ERF_LSP.33
!------Hardwired version of small x approximation----------------------- ERF_LSP.34
ERF_LSP = XERF * TERF * TERF * 0.1 ERF_LSP.35
& * (1. - ((5. / 21.) * TERF * (1. - ((7. / 36.) * TERF ERF_LSP.36
& * (1. - ((9. / 55.) * TERF * (1. - ((11. / 78.) * TERF )))) ERF_LSP.37
& )))) ERF_LSP.38
! ERF_LSP.39
ELSEIF (ABS(XERF) .GT. 1.161 .AND. ABS(XERF) .LT. 6.) THEN ERF_LSP.40
!------Continued fraction approximation at large xerf------------------- ERF_LSP.41
ERF_LSP = SIGN((SQRT_PI * 0.5),XERF) - (XERF * (1. - (TERF / 3.) ERF_LSP.42
& + (EXP(- TERF) / (2. * (TERF + 0.5 - (0.5 / (TERF + 2.5 ERF_LSP.43
& - (3. / (TERF + 4.5 - (7.5 / (TERF + 6.5)))))))) ))) ERF_LSP.44
! ERF_LSP.45
ELSE ERF_LSP.46
ERF_LSP = SIGN((SQRT_PI * 0.5),XERF) + (XERF * ((TERF / 3.) - 1.)) ERF_LSP.47
! ERF_LSP.48
ENDIF ERF_LSP.49
C*---------------------------------------------------------------------- ERF_LSP.50