source: flexpart.git/src/random_mod.f90 @ 7999df47

devrelease-10univie
Last change on this file since 7999df47 was 5f9d14a, checked in by Espen Sollum ATMOS <eso@…>, 5 years ago

Updated wet depo scheme

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