source: trunk/src/random.f90 @ 28

Last change on this file since 28 was 4, checked in by mlanger, 11 years ago
File size: 4.0 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!  Taken from Press et al., Numerical Recipes
23
24function ran1(idum)
25
26  implicit none
27
28  integer :: idum
29  real    :: ran1
30  integer,parameter :: ia=16807, im=2147483647, iq=127773, ir=2836
31  integer,parameter :: ntab=32, ndiv=1+(im-1)/ntab
32  real,parameter    :: am=1./im, eps=1.2e-7, rnmx=1.-eps
33  integer :: j, k
34  integer :: iv(ntab) = (/ (0,j=1,ntab) /)
35  integer :: iy=0
36
37  if (idum.le.0.or.iy.eq.0) then
38    idum=max(-idum,1)
39    do j=ntab+8,1,-1
40      k=idum/iq
41      idum=ia*(idum-k*iq)-ir*k
42      if (idum.lt.0) idum=idum+im
43      if (j.le.ntab) iv(j)=idum
44    enddo
45    iy=iv(1)
46  endif
47  k=idum/iq
48  idum=ia*(idum-k*iq)-ir*k
49  if (idum.lt.0) idum=idum+im
50  j=1+iy/ndiv
51  iy=iv(j)
52  iv(j)=idum
53  ran1=min(am*iy,rnmx)
54end function ran1
55
56
57function gasdev(idum)
58
59  implicit none
60
61  integer :: idum
62  real    :: gasdev, fac, r, v1, v2
63  integer :: iset = 0
64  real    :: gset = 0.
65  real, external :: ran3
66
67  if (iset.eq.0) then
681   v1=2.*ran3(idum)-1.
69    v2=2.*ran3(idum)-1.
70    r=v1**2+v2**2
71    if(r.ge.1.0 .or. r.eq.0.0) go to 1
72    fac=sqrt(-2.*log(r)/r)
73    gset=v1*fac
74    gasdev=v2*fac
75    iset=1
76  else
77    gasdev=gset
78    iset=0
79  endif
80end function gasdev
81
82
83subroutine gasdev1(idum,random1,random2)
84
85  implicit none
86
87  integer :: idum
88  real :: random1, random2, fac, v1, v2, r
89  real, external :: ran3
90
911   v1=2.*ran3(idum)-1.
92  v2=2.*ran3(idum)-1.
93  r=v1**2+v2**2
94  if(r.ge.1.0 .or. r.eq.0.0) go to 1
95  fac=sqrt(-2.*log(r)/r)
96  random1=v1*fac
97  random2=v2*fac
98  ! Limit the random numbers to lie within the interval -3 and +3
99  !**************************************************************
100  if (random1.lt.-3.) random1=-3.
101  if (random2.lt.-3.) random2=-3.
102  if (random1.gt.3.) random1=3.
103  if (random2.gt.3.) random2=3.
104end subroutine gasdev1
105
106
107function ran3(idum)
108
109  implicit none
110
111  integer :: idum
112  real :: ran3
113
114  integer,parameter :: mbig=1000000000, mseed=161803398, mz=0
115  real,parameter    :: fac=1./mbig
116  integer :: i,ii,inext,inextp,k
117  integer :: mj,mk,ma(55)
118
119  save inext,inextp,ma
120  integer :: iff = 0
121
122  if(idum.lt.0.or.iff.eq.0)then
123    iff=1
124    mj=mseed-iabs(idum)
125    mj=mod(mj,mbig)
126    ma(55)=mj
127    mk=1
128    do i=1,54
129      ii=mod(21*i,55)
130      ma(ii)=mk
131      mk=mj-mk
132      if(mk.lt.mz)mk=mk+mbig
133      mj=ma(ii)
134    end do
135    do k=1,4
136      do i=1,55
137        ma(i)=ma(i)-ma(1+mod(i+30,55))
138        if(ma(i).lt.mz)ma(i)=ma(i)+mbig
139      end do
140    end do
141    inext=0
142    inextp=31
143    idum=1
144  endif
145  inext=inext+1
146  if(inext.eq.56)inext=1
147  inextp=inextp+1
148  if(inextp.eq.56)inextp=1
149  mj=ma(inext)-ma(inextp)
150  if(mj.lt.mz)mj=mj+mbig
151  ma(inext)=mj
152  ran3=mj*fac
153end function ran3
154!  (C) Copr. 1986-92 Numerical Recipes Software US.
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG