[16] | 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 |
---|