source: trunk/src/erf.f90 @ 28

Last change on this file since 28 was 4, checked in by mlanger, 11 years ago
  • Property svn:executable set to *
File size: 4.9 KB
Line 
1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22! To be used, if the non-standard Fortran function erf does not exist on
23! your machine
24!
25!aus:  Numerical Recipes (FORTRAN) / Chapter 6.
26!
27!6.1  FUNCTION GAMMLN
28!6.2  FUNCTION GAMMP   <6.2:GSER/6.2:GCF/6.1:GAMMLN>
29!6.2  FUNCTION GAMMQ   <6.2:GSER/6.2:GCF/6.1:GAMMLN>
30!6.2  SUBROUTINE GSER    <6.1:GAMMLN>
31!6.2  SUBROUTINE GCF     <6.1:GAMMLN>
32!6.2  FUNCTION ERF     <6.2:GAMMP/6.2:GSER/6.2:GCF/6.1:GAMMLN>
33!6.2  FUNCTION ERFC    <6.2.:GAMMP/6.2:GAMMQ/6.2:GSER/
34!                       6.2:GCF/6.1:GAMMLN>
35!6.2  FUNCTION ERFCC
36
37function gammln(xx)
38
39  use par_mod, only: dp
40
41  implicit none
42
43  integer :: j
44  real :: x,tmp,ser,xx,gammln
45  real :: cof(6) = (/ &
46       76.18009173_dp, -86.50532033_dp, 24.01409822_dp, &
47       -1.231739516_dp, .120858003e-2_dp, -.536382e-5_dp    /)
48  real :: stp = 2.50662827465_dp
49  real :: half = 0.5_dp, one = 1.0_dp, fpf = 5.5_dp
50
51  x=xx-one
52  tmp=x+fpf
53  tmp=(x+half)*log(tmp)-tmp
54  ser=one
55  do j=1,6
56    x=x+one
57    ser=ser+cof(j)/x
58  end do
59  gammln=tmp+log(stp*ser)
60end function gammln
61
62function gammp(a,x)
63
64  implicit none
65
66  real :: a, x, gln, gamser, gammp, gammcf
67
68  if(x .lt. 0. .or. a .le. 0.) then
69     print*, 'gammp'
70     stop
71  end if
72  if(x.lt.a+1.)then
73    call gser(gamser,a,x,gln)
74    gammp=gamser
75  else
76    call gcf(gammcf,a,x,gln)
77    gammp=1.-gammcf
78  endif
79end function gammp
80
81function gammq(a,x)
82
83  implicit none
84
85  real :: a, x, gln, gamser, gammq, gammcf
86
87  if(x.lt.0..or.a.le.0.) then
88     print*, 'gammq'
89     stop
90  end if
91  if(x.lt.a+1.)then
92    call gser(gamser,a,x,gln)
93    gammq=1.-gamser
94  else
95    call gcf(gammcf,a,x,gln)
96    gammq=gammcf
97  endif
98end function gammq
99
100subroutine gser(gamser,a,x,gln)
101
102  implicit none
103
104  integer :: n
105  real :: gamser, a, x, gln, ap, summ, del
106  real, external :: gammln
107
108  integer,parameter :: itmax=100
109  real,parameter    :: eps=3.e-7
110
111  gln=gammln(a)
112  if(x.le.0.)then
113    if(x.lt.0.) then
114       print*, 'gser'
115       stop
116    end if
117    gamser=0.
118    return
119  endif
120  ap=a
121  summ=1./a
122  del=summ
123  do n=1,itmax
124    ap=ap+1.
125    del=del*x/ap
126    summ=summ+del
127    if(abs(del).lt.abs(summ)*eps)go to 1
128  end do
129  print*, 'gser: a too large, itmax too small'
130  stop
1311   gamser=summ*exp(-x+a*log(x)-gln)
132end subroutine gser
133
134subroutine gcf(gammcf,a,x,gln)
135
136  implicit none
137
138  integer :: n
139  real :: gammcf, a, x, gln, gold, a0, a1, b0, b1, fac, an, anf, ana, g
140  real, external :: gammln
141
142  integer,parameter :: itmax=100
143  real,parameter    :: eps=3.e-7
144
145  gln=gammln(a)
146  gold=0.
147  a0=1.
148  a1=x
149  b0=0.
150  b1=1.
151  fac=1.
152  do n=1,itmax
153    an=real(n)
154    ana=an-a
155    a0=(a1+a0*ana)*fac
156    b0=(b1+b0*ana)*fac
157    anf=an*fac
158    a1=x*a0+anf*a1
159    b1=x*b0+anf*b1
160    if(a1.ne.0.)then
161      fac=1./a1
162      g=b1*fac
163      if(abs((g-gold)/g).lt.eps)go to 1
164      gold=g
165    endif
166  end do
167  print*, 'gcf: a too large, itmax too small'
168  stop
1691   gammcf=exp(-x+a*alog(x)-gln)*g
170end subroutine gcf
171
172function erf(x)
173
174  implicit none
175
176  real :: x, erf
177  real, external :: gammp
178
179  if(x.lt.0.)then
180    erf=-gammp(.5,x**2)
181  else
182    erf=gammp(.5,x**2)
183  endif
184end function erf
185
186function erfc(x)
187
188  implicit none
189
190  real :: x, erfc
191  real, external :: gammp, gammq
192
193  if(x.lt.0.)then
194    erfc=1.+gammp(.5,x**2)
195  else
196    erfc=gammq(.5,x**2)
197  endif
198end function erfc
199
200function erfcc(x)
201
202  implicit none
203
204  real :: x, z, t, erfcc
205
206  z=abs(x)
207  t=1./(1.+0.5*z)
208  erfcc=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+ &
209       t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+ &
210       t*(1.48851587+t*(-.82215223+t*.17087277)))))))))
211  if (x.lt.0.) erfcc=2.-erfcc
212end function erfcc
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG