! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt ! SPDX-License-Identifier: GPL-3.0-or-later ! Taken from Press et al., Numerical Recipes module random_mod implicit none integer, parameter :: ran1_ntab=32 integer, allocatable :: ran1_iv(:,:), ran1_iy(:) integer, allocatable :: gasdev_iset(:) real, allocatable :: gasdev_gset(:) integer, allocatable :: ran3_iff(:) integer, allocatable :: ran3_inext(:),ran3_inextp(:) integer, allocatable :: ma(:,:) integer, allocatable :: iseed1(:), iseed2(:) contains subroutine alloc_random(num_threads) implicit none integer :: num_threads, i,stat allocate(ran1_iv(ran1_ntab,0:num_threads-1),ran1_iy(0:num_threads-1), stat=stat) if (stat.ne.0) error stop "Could not allocate ran1_iv" allocate(gasdev_iset(0:num_threads-1),gasdev_gset(0:num_threads-1), stat=stat) if (stat.ne.0) error stop "Could not allocate gasdev_iset" allocate(ran3_iff(0:num_threads-1),ran3_inext(0:num_threads-1),ran3_inextp(0:num_threads-1), stat=stat) if (stat.ne.0) error stop "Could not allocate ran3_iff" allocate(ma(55,0:num_threads-1), stat=stat) if (stat.ne.0) error stop "Could not allocate ma" allocate(iseed1(0:num_threads-1),iseed2(0:num_threads-1), stat=stat) if (stat.ne.0) error stop "Could not allocate iseed1" do i=0,num_threads-1 iseed1(i) = -7-i iseed2(i) = -88-i end do ran3_iff(0:num_threads-1)=0 ran1_iv(:,0:num_threads-1)=0 ran1_iy(0:num_threads-1)=0 gasdev_iset(0:num_threads-1)=0 gasdev_gset(0:num_threads-1)=0 end subroutine alloc_random subroutine dealloc_random() deallocate(ran1_iv,ran1_iy) deallocate(gasdev_iset,gasdev_gset) deallocate(ran3_iff,ran3_inext,ran3_inextp) deallocate(ma) deallocate(iseed1,iseed2) end subroutine dealloc_random function ran1(idum,ithread) implicit none integer :: idum,ithread real :: ran1 integer,parameter :: ia=16807, im=2147483647, iq=127773, ir=2836 integer,parameter :: ndiv=1+(im-1)/ran1_ntab real,parameter :: am=1./im, eps=1.2e-7, rnmx=1.-eps integer :: j, k if (idum.le.0.or.ran1_iy(ithread).eq.0) then idum=max(-idum,1) do j=ran1_ntab+8,1,-1 k=idum/iq idum=ia*(idum-k*iq)-ir*k if (idum.lt.0) idum=idum+im if (j.le.ran1_ntab) ran1_iv(j,ithread)=idum enddo ran1_iy(ithread)=ran1_iv(1,ithread) endif k=idum/iq idum=ia*(idum-k*iq)-ir*k if (idum.lt.0) idum=idum+im j=1+ran1_iy(ithread)/ndiv ran1_iy(ithread)=ran1_iv(j,ithread) ran1_iv(j,ithread)=idum ran1=min(am*ran1_iy(ithread),rnmx) end function ran1 function gasdev(idum,ithread) implicit none integer :: idum,ithread real :: gasdev, fac, r, v1, v2 if (gasdev_iset(ithread).eq.0) then 1 v1=2.*ran3(idum,ithread)-1. v2=2.*ran3(idum,ithread)-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) gasdev_gset(ithread)=v1*fac gasdev=v2*fac gasdev_iset(ithread)=1 else gasdev=gasdev_gset(ithread) gasdev_iset(ithread)=0 endif end function gasdev subroutine gasdev1(idum,random1,random2) implicit none integer :: idum real :: random1, random2, fac, v1, v2, r 1 v1=2.*ran3(idum,0)-1. v2=2.*ran3(idum,0)-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,ithread) implicit none integer :: idum,ithread real :: ran3 integer,parameter :: mbig=1000000000, mseed=161803398, mz=0 real,parameter :: fac=1./mbig integer :: i,ii,k integer :: mj,mk if(idum.lt.0 .or. ran3_iff(ithread).eq.0)then ran3_iff(ithread)=1 mj=mseed-iabs(idum) mj=mod(mj,mbig) ma(55,ithread)=mj mk=1 do i=1,54 ii=mod(21*i,55) ma(ii,ithread)=mk mk=mj-mk if(mk.lt.mz)mk=mk+mbig mj=ma(ii,ithread) end do do k=1,4 do i=1,55 ma(i,ithread)=ma(i,ithread)-ma(1+mod(i+30,55),ithread) if(ma(i,ithread).lt.mz) ma(i,ithread)=ma(i,ithread)+mbig end do end do ran3_inext(ithread)=0 ran3_inextp(ithread)=31 idum=1 endif ran3_inext(ithread)=ran3_inext(ithread)+1 if(ran3_inext(ithread).eq.56) ran3_inext(ithread)=1 ran3_inextp(ithread)=ran3_inextp(ithread)+1 if(ran3_inextp(ithread).eq.56) ran3_inextp(ithread)=1 mj=ma(ran3_inext(ithread),ithread)-ma(ran3_inextp(ithread),ithread) if(mj.lt.mz)mj=mj+mbig ma(ran3_inext(ithread),ithread)=mj ran3=mj*fac end function ran3 ! (C) Copr. 1986-92 Numerical Recipes Software US. end module random_mod