source: branches/flexpart91_hasod/src_parallel/random_mod.f90 @ 8

Last change on this file since 8 was 8, checked in by hasod, 11 years ago

Added parallel version of Flexpart91

File size: 5.2 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
22module random_mod
23
24  use com_mod, only: maxpart
25  ! read: maxpart
26  ! write:
27
28  implicit none
29
30  private
31
32  integer, parameter :: ran1_ntab=32
33  integer :: ran1_iv(ran1_ntab) = 0, ran1_iy = 0
34
35  integer :: gasdev_iset = 0
36  real    :: gasdev_gset
37
38  integer :: ran3_iff = 0
39  integer :: ran3_inext,ran3_inextp
40  integer :: ran3_ma(55)
41  real    :: ran3_initialize(maxpart), ran3_advance(maxpart), ran3_conv(maxpart)
42
43  integer :: idummy_ran3_initialize = -7, idummy_ran3_advance = -7, idummy_ran3_conv = -7
44  ! idummy_ran3_initialize, idummy_ran3_advance, idummy_ran3_conv  Initial value for calculation
45  !                        of random numbers (see timemanager_computations, convmix_tiedtke,
46  !                        and ran3)
47
48  public :: ran1, gasdev, gasdev1, ran3
49  public :: ran3_initialize, ran3_advance, ran3_conv
50  public :: idummy_ran3_initialize, idummy_ran3_advance, idummy_ran3_conv
51
52contains
53  !  Taken from Press et al., Numerical Recipes
54
55  function ran1(idum)
56    implicit none
57    integer :: idum
58    real :: ran1
59    integer,parameter :: ia=16807, im=2147483647, iq=127773, ir=2836
60    integer,parameter :: ndiv=1+(im-1)/ran1_ntab
61    real,parameter :: am=1./real(im), eps=1.2e-7,rnmx=1.-eps
62    integer :: j,k
63    if (idum.le.0.or.ran1_iy.eq.0) then
64      idum=max(-idum,1)
65      do j=ran1_ntab+8,1,-1
66        k=idum/iq
67        idum=ia*(idum-k*iq)-ir*k
68        if (idum.lt.0) idum=idum+im
69        if (j.le.ran1_ntab) ran1_iv(j)=idum
70      enddo
71      ran1_iy=ran1_iv(1)
72    endif
73    k=idum/iq
74    idum=ia*(idum-k*iq)-ir*k
75    if (idum.lt.0) idum=idum+im
76    j=1+ran1_iy/ndiv
77    ran1_iy=ran1_iv(j)
78    ran1_iv(j)=idum
79    ran1=min(am*real(ran1_iy),rnmx)
80    return
81  end function ran1
82
83  function gasdev(idum)
84    implicit none
85    integer :: idum
86    real :: gasdev, fac, v1, v2, r
87    if (gasdev_iset.eq.0) then
881     v1=2.*ran3(idum)-1.
89      v2=2.*ran3(idum)-1.
90      r=v1**2+v2**2
91      if(r.ge.1.0 .or. r.eq.0.0) go to 1
92      fac=sqrt(-2.*log(r)/r)
93      gasdev_gset=v1*fac
94      gasdev=v2*fac
95      gasdev_iset=1
96    else
97      gasdev=gasdev_gset
98      gasdev_iset=0
99    endif
100    return
101  end function gasdev
102
103  subroutine gasdev1(idum,random1,random2)
104    implicit none
105    integer :: idum
106    real :: random1, random2, v1, v2, r, fac
1071     v1=2.*ran3(idum)-1.
108    v2=2.*ran3(idum)-1.
109    r=v1**2+v2**2
110    if(r.ge.1.0 .or. r.eq.0.0) go to 1
111    fac=sqrt(-2.*log(r)/r)
112    random1=v1*fac
113    random2=v2*fac
114    ! Limit the random numbers to lie within the interval -3 and +3
115    !**************************************************************
116    if (random1.lt.-3.) random1=-3.
117    if (random2.lt.-3.) random2=-3.
118    if (random1.gt.3.) random1=3.
119    if (random2.gt.3.) random2=3.
120    return
121  end subroutine gasdev1
122
123  function ran3(idum)
124    implicit none
125    integer :: idum
126    real :: ran3
127    integer, parameter :: mbig=1000000000, mseed=161803398, mz=0
128    real, parameter :: fac=1./real(mbig)
129    integer :: i,ii,k
130    integer :: mj,mk
131
132    if(idum.lt.0.or.ran3_iff.eq.0)then
133      ran3_iff=1
134      mj=mseed-iabs(idum)
135      mj=mod(mj,mbig)
136      ran3_ma(55)=mj
137      mk=1
138      do i=1,54
139        ii=mod(21*i,55)
140        ran3_ma(ii)=mk
141        mk=mj-mk
142        if(mk.lt.mz)mk=mk+mbig
143        mj=ran3_ma(ii)
144      end do
145      do k=1,4
146        do i=1,55
147          ran3_ma(i)=ran3_ma(i)-ran3_ma(1+mod(i+30,55))
148          if(ran3_ma(i).lt.mz)ran3_ma(i)=ran3_ma(i)+mbig
149        end do
150      end do
151      ran3_inext=0
152      ran3_inextp=31
153      idum=1
154    endif
155    ran3_inext=ran3_inext+1
156    if(ran3_inext.eq.56)ran3_inext=1
157    ran3_inextp=ran3_inextp+1
158    if(ran3_inextp.eq.56)ran3_inextp=1
159    mj=ran3_ma(ran3_inext)-ran3_ma(ran3_inextp)
160    if(mj.lt.mz)mj=mj+mbig
161    ran3_ma(ran3_inext)=mj
162    ran3=real(mj)*fac
163    return
164  end function ran3
165  !  (C) Copr. 1986-92 Numerical Recipes Software us.
166
167end module random_mod
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG