1 | !*********************************************************************** |
---|
2 | !* Copyright 2012,2013 * |
---|
3 | !* Jerome Brioude, Delia Arnold, Jerome Fast, |
---|
4 | !* Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa * |
---|
5 | !* Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * |
---|
6 | !* * |
---|
7 | !* This file is part of FLEXPART WRF * |
---|
8 | !* * |
---|
9 | !* FLEXPART is free software: you can redistribute it and/or modify * |
---|
10 | !* it under the terms of the GNU General Public License as published by* |
---|
11 | !* the Free Software Foundation, either version 3 of the License, or * |
---|
12 | !* (at your option) any later version. * |
---|
13 | !* * |
---|
14 | !* FLEXPART is distributed in the hope that it will be useful, * |
---|
15 | !* but WITHOUT ANY WARRANTY; without even the implied warranty of * |
---|
16 | !* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * |
---|
17 | !* GNU General Public License for more details. * |
---|
18 | !* * |
---|
19 | !* You should have received a copy of the GNU General Public License * |
---|
20 | !* along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * |
---|
21 | !*********************************************************************** |
---|
22 | ! |
---|
23 | ! J. Brioude: the version in fortran 9 gives some error. better use f77 |
---|
24 | ! |
---|
25 | c To be used, if the non-standard Fortran function erf does not exist on |
---|
26 | c your machine |
---|
27 | C |
---|
28 | C aus: Numerical Recipes (FORTRAN) / Chapter 6. |
---|
29 | C |
---|
30 | C 6.1 FUNCTION GAMMLN |
---|
31 | C 6.2 FUNCTION GAMMP <6.2:GSER/6.2:GCF/6.1:GAMMLN> |
---|
32 | C 6.2 FUNCTION GAMMQ <6.2:GSER/6.2:GCF/6.1:GAMMLN> |
---|
33 | C 6.2 SUBROUTINE GSER <6.1:GAMMLN> |
---|
34 | C 6.2 SUBROUTINE GCF <6.1:GAMMLN> |
---|
35 | C 6.2 FUNCTION ERF <6.2:GAMMP/6.2:GSER/6.2:GCF/6.1:GAMMLN> |
---|
36 | C 6.2 FUNCTION ERFC <6.2.:GAMMP/6.2:GAMMQ/6.2:GSER/ |
---|
37 | C 6.2:GCF/6.1:GAMMLN> |
---|
38 | C 6.2 FUNCTION ERFCC |
---|
39 | C |
---|
40 | FUNCTION GAMMLN(XX) |
---|
41 | REAL*8 COF(6),STP,HALF,ONE,FPF,X,TMP,SER |
---|
42 | DATA COF,STP/76.18009173D0,-86.50532033D0,24.01409822D0, |
---|
43 | * -1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/ |
---|
44 | DATA HALF,ONE,FPF/0.5D0,1.0D0,5.5D0/ |
---|
45 | X=XX-ONE |
---|
46 | TMP=X+FPF |
---|
47 | TMP=(X+HALF)*LOG(TMP)-TMP |
---|
48 | SER=ONE |
---|
49 | DO 11 J=1,6 |
---|
50 | X=X+ONE |
---|
51 | SER=SER+COF(J)/X |
---|
52 | 11 CONTINUE |
---|
53 | GAMMLN=TMP+LOG(STP*SER) |
---|
54 | RETURN |
---|
55 | END |
---|
56 | C |
---|
57 | FUNCTION GAMMP(A,X) |
---|
58 | IF(X.LT.0..OR.A.LE.0.) STOP 'GAMMP' |
---|
59 | IF(X.LT.A+1.)THEN |
---|
60 | CALL GSER(GAMSER,A,X,GLN) |
---|
61 | GAMMP=GAMSER |
---|
62 | ELSE |
---|
63 | CALL GCF(GAMMCF,A,X,GLN) |
---|
64 | GAMMP=1.-GAMMCF |
---|
65 | ENDIF |
---|
66 | RETURN |
---|
67 | END |
---|
68 | C |
---|
69 | FUNCTION GAMMQ(A,X) |
---|
70 | IF(X.LT.0..OR.A.LE.0.) STOP 'GAMMQ' |
---|
71 | IF(X.LT.A+1.)THEN |
---|
72 | CALL GSER(GAMSER,A,X,GLN) |
---|
73 | GAMMQ=1.-GAMSER |
---|
74 | ELSE |
---|
75 | CALL GCF(GAMMCF,A,X,GLN) |
---|
76 | GAMMQ=GAMMCF |
---|
77 | ENDIF |
---|
78 | RETURN |
---|
79 | END |
---|
80 | C |
---|
81 | SUBROUTINE GSER(GAMSER,A,X,GLN) |
---|
82 | PARAMETER (ITMAX=100,EPS=3.E-7) |
---|
83 | GLN=GAMMLN(A) |
---|
84 | IF(X.LE.0.)THEN |
---|
85 | ! IF(X.LT.0.) PAUSE 'GSER' |
---|
86 | GAMSER=0. |
---|
87 | RETURN |
---|
88 | ENDIF |
---|
89 | AP=A |
---|
90 | SUM=1./A |
---|
91 | DEL=SUM |
---|
92 | DO 11 N=1,ITMAX |
---|
93 | AP=AP+1. |
---|
94 | DEL=DEL*X/AP |
---|
95 | SUM=SUM+DEL |
---|
96 | IF(ABS(DEL).LT.ABS(SUM)*EPS)GO TO 1 |
---|
97 | 11 CONTINUE |
---|
98 | STOP 'GSER: A too large, ITMAX too small' |
---|
99 | 1 GAMSER=SUM*EXP(-X+A*LOG(X)-GLN) |
---|
100 | RETURN |
---|
101 | END |
---|
102 | C |
---|
103 | SUBROUTINE GCF(GAMMCF,A,X,GLN) |
---|
104 | PARAMETER (ITMAX=100,EPS=3.E-7) |
---|
105 | GLN=GAMMLN(A) |
---|
106 | GOLD=0. |
---|
107 | A0=1. |
---|
108 | A1=X |
---|
109 | B0=0. |
---|
110 | B1=1. |
---|
111 | FAC=1. |
---|
112 | DO 11 N=1,ITMAX |
---|
113 | AN=FLOAT(N) |
---|
114 | ANA=AN-A |
---|
115 | A0=(A1+A0*ANA)*FAC |
---|
116 | B0=(B1+B0*ANA)*FAC |
---|
117 | ANF=AN*FAC |
---|
118 | A1=X*A0+ANF*A1 |
---|
119 | B1=X*B0+ANF*B1 |
---|
120 | IF(A1.NE.0.)THEN |
---|
121 | FAC=1./A1 |
---|
122 | G=B1*FAC |
---|
123 | IF(ABS((G-GOLD)/G).LT.EPS)GO TO 1 |
---|
124 | GOLD=G |
---|
125 | ENDIF |
---|
126 | 11 CONTINUE |
---|
127 | STOP 'GCF: A too large, ITMAX too small' |
---|
128 | 1 GAMMCF=EXP(-X+A*ALOG(X)-GLN)*G |
---|
129 | RETURN |
---|
130 | END |
---|
131 | C |
---|
132 | FUNCTION erf(X) |
---|
133 | IF(X.LT.0.)THEN |
---|
134 | erf=-GAMMP(.5,X**2) |
---|
135 | ELSE |
---|
136 | erf=GAMMP(.5,X**2) |
---|
137 | ENDIF |
---|
138 | RETURN |
---|
139 | END |
---|
140 | C |
---|
141 | FUNCTION ERFC(X) |
---|
142 | IF(X.LT.0.)THEN |
---|
143 | ERFC=1.+GAMMP(.5,X**2) |
---|
144 | ELSE |
---|
145 | ERFC=GAMMQ(.5,X**2) |
---|
146 | ENDIF |
---|
147 | RETURN |
---|
148 | END |
---|
149 | C |
---|
150 | FUNCTION ERFCC(X) |
---|
151 | Z=ABS(X) |
---|
152 | T=1./(1.+0.5*Z) |
---|
153 | ERFCC=T*EXP(-Z*Z-1.26551223+T*(1.00002368+T*(.37409196+ |
---|
154 | * T*(.09678418+T*(-.18628806+T*(.27886807+T*(-1.13520398+ |
---|
155 | * T*(1.48851587+T*(-.82215223+T*.17087277))))))))) |
---|
156 | IF (X.LT.0.) ERFCC=2.-ERFCC |
---|
157 | RETURN |
---|
158 | END |
---|
159 | C |
---|