Commit 6f3a20a6 authored by Espen Sollum's avatar Espen Sollum
Browse files

Deleted random.f90 (now using random_mod.f90 instead)

parent 2f8635b3
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART is free software: you can redistribute it and/or modify *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or *
! (at your option) any later version. *
! *
! FLEXPART is distributed in the hope that it will be useful, *
! but WITHOUT ANY WARRANTY; without even the implied warranty of *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
! GNU General Public License for more details. *
! *
! You should have received a copy of the GNU General Public License *
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
! Taken from Press et al., Numerical Recipes
function ran1(idum)
implicit none
integer :: idum
real :: ran1
integer,parameter :: ia=16807, im=2147483647, iq=127773, ir=2836
integer,parameter :: ntab=32, ndiv=1+(im-1)/ntab
real,parameter :: am=1./im, eps=1.2e-7, rnmx=1.-eps
integer :: j, k
integer :: iv(ntab) = (/ (0,j=1,ntab) /)
integer :: iy=0
if (idum.le.0.or.iy.eq.0) then
idum=max(-idum,1)
do j=ntab+8,1,-1
k=idum/iq
idum=ia*(idum-k*iq)-ir*k
if (idum.lt.0) idum=idum+im
if (j.le.ntab) iv(j)=idum
enddo
iy=iv(1)
endif
k=idum/iq
idum=ia*(idum-k*iq)-ir*k
if (idum.lt.0) idum=idum+im
j=1+iy/ndiv
iy=iv(j)
iv(j)=idum
ran1=min(am*iy,rnmx)
end function ran1
function gasdev(idum)
implicit none
integer :: idum
real :: gasdev, fac, r, v1, v2
integer :: iset = 0
real :: gset = 0.
real, external :: ran3
if (iset.eq.0) then
1 v1=2.*ran3(idum)-1.
v2=2.*ran3(idum)-1.
r=v1**2+v2**2
if(r.ge.1.0 .or. r.eq.0.0) go to 1
fac=sqrt(-2.*log(r)/r)
gset=v1*fac
gasdev=v2*fac
iset=1
else
gasdev=gset
iset=0
endif
end function gasdev
subroutine gasdev1(idum,random1,random2)
implicit none
integer :: idum
real :: random1, random2, fac, v1, v2, r
real, external :: ran3
1 v1=2.*ran3(idum)-1.
v2=2.*ran3(idum)-1.
r=v1**2+v2**2
if(r.ge.1.0 .or. r.eq.0.0) go to 1
fac=sqrt(-2.*log(r)/r)
random1=v1*fac
random2=v2*fac
! Limit the random numbers to lie within the interval -3 and +3
!**************************************************************
if (random1.lt.-3.) random1=-3.
if (random2.lt.-3.) random2=-3.
if (random1.gt.3.) random1=3.
if (random2.gt.3.) random2=3.
end subroutine gasdev1
function ran3(idum)
implicit none
integer :: idum
real :: ran3
integer,parameter :: mbig=1000000000, mseed=161803398, mz=0
real,parameter :: fac=1./mbig
integer :: i,ii,inext,inextp,k
integer :: mj,mk,ma(55)
save inext,inextp,ma
integer :: iff = 0
if(idum.lt.0.or.iff.eq.0)then
iff=1
mj=mseed-iabs(idum)
mj=mod(mj,mbig)
ma(55)=mj
mk=1
do i=1,54
ii=mod(21*i,55)
ma(ii)=mk
mk=mj-mk
if(mk.lt.mz)mk=mk+mbig
mj=ma(ii)
end do
do k=1,4
do i=1,55
ma(i)=ma(i)-ma(1+mod(i+30,55))
if(ma(i).lt.mz)ma(i)=ma(i)+mbig
end do
end do
inext=0
inextp=31
idum=1
endif
inext=inext+1
if(inext.eq.56)inext=1
inextp=inextp+1
if(inextp.eq.56)inextp=1
mj=ma(inext)-ma(inextp)
if(mj.lt.mz)mj=mj+mbig
ma(inext)=mj
ran3=mj*fac
end function ran3
! (C) Copr. 1986-92 Numerical Recipes Software US.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment