model/src/port_rand.F

C $Header: /u/gcmpack/MITgcm/model/src/port_rand.F,v 1.4 2004/07/19 15:13:07 molod Exp $
C $Name:  $

C Portable random number generator

#undef _USE_INTEGERS

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C     !ROUTINE: port_rand
C     !INTERFACE:
      real*8 function port_rand(seed)

C     !DESCRIPTION:
C     Portable random number generator
      
C     !USES:
      implicit none
CEOP

C     !LOCAL VARIABLES:
      integer nff
      parameter(nff=55)
#ifdef _USE_INTEGERS
      integer mbig,mseed,mZ,seed
#else
      real*8 mbig,mseed,mz,seed
#endif
      real*8 fac
#ifdef _USE_INTEGERS
      parameter (mbig=1000000000,mz=0,fac=1./mbig)
      data mseed/161803398/
#else
      parameter (mbig=4000000.,mz=0.,fac=1./mbig)
      data mseed/1618033./
#endif
      integer i,ii,inext,inextp,k,idum
      parameter(idum=-2)
#ifdef _USE_INTEGERS
      integer mj,mk,ma(nff)
#else
      real*8 mj,mk,ma(nff)
#endif
      logical firstCall
      save firstCall,inext,inextp,ma
      data firstCall /.true./

      mseed = seed
      if(firstCall)then
        firstCall=.false.
        mj=mseed-iabs(idum)
        mj=mod(mj,mbig)
        ma(nff)=mj
        mk=1
        do i=1,nff-1
          ii=mod(21*i,nff)
          ma(ii)=mk
          mk=mj-mk
          if(mk.lt.mz)mk=mk+mbig
          mj=ma(ii)
        enddo
        do k=1,4
          do i=1,nff
            ma(i)=ma(i)-ma(1+mod(i+30,nff))
            if(ma(i).lt.MZ)ma(i)=ma(i)+mbig
          enddo
        enddo
        inext=0
        inextp=31
      endif
      inext=mod(inext,nff)+1
      inextp=mod(inextp,nff)+1
      mj=ma(inext)-ma(inextp)
      if(mj.lt.MZ)mj=mj+mbig
      ma(inext)=mj
      port_rand=mj*fac
      return

      end

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      subroutine port_ranarr(n,arr)
      implicit none
      integer n,i
      real arr(n)
      real*8 port_rand,seed
      seed=1618033.0d0
      do i=1,n
       arr(i)=port_rand(seed)
      enddo

      return
      end

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C     !ROUTINE: port_rand_norm
C     !INTERFACE:
      real*8 function port_rand_norm()

C     !DESCRIPTION:
C     This function generates a normally distributed random number with
C     the so called polar algorithm. The algorithm actually generates 2
C     numbers, but only 1 is returned for maximum compatibility with old
C     code.  The most obvious way to improve this function would be to
C     make sure that the second number is not wasted.

C     Changed: 2004.09.06 antti.westerlund@fimr.fi

C     !USES:
      implicit none
CEOP

C     !LOCAL VARIABLES:
      real*8 port_rand,seed
      real*8 x1, x2, xs, t           
      integer i

      seed=1618033.0d0
C     first generate 2 equally distributed random numbers (-1,1) 
      DO WHILE (1 .eq. 1)
         x1=2.0*port_rand(seed)-1.0
         x2=2.0*port_rand(seed)-1.0
         xs=x1**2+x2**2
         if(xs .lt. 1.0 .and. xs .ne. 0.0) then
            goto 100      
         end if
      END DO
 100  continue

      t = sqrt(-2.0*log(xs)/xs)
      port_rand_norm = t*x1
C
C     also t*x2 would be a gaussian random number and could be returned

      return
      end
1	edhill	1.5	C $Header: /u/gcmpack/MITgcm/model/src/port_rand.F,v 1.4 2004/07/19 15:13:07 molod Exp $
2	cnh	1.3	C $Name: $
3	adcroft	1.2
4			C Portable random number generator
5
6			#undef _USE_INTEGERS
7
8	edhill	1.5	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
9	cnh	1.3	CBOP
10			C !ROUTINE: port_rand
11			C !INTERFACE:
12	molod	1.4	real*8 function port_rand(seed)
13	cnh	1.3
14	edhill	1.5	C !DESCRIPTION:
15	cnh	1.3	C Portable random number generator
16	edhill	1.5
17			C !USES:
18			implicit none
19			CEOP
20	cnh	1.3
21			C !LOCAL VARIABLES:
22	adcroft	1.2	integer nff
23			parameter(nff=55)
24			#ifdef _USE_INTEGERS
25	molod	1.4	integer mbig,mseed,mZ,seed
26	adcroft	1.2	#else
27	molod	1.4	real*8 mbig,mseed,mz,seed
28	adcroft	1.2	#endif
29			real*8 fac
30			#ifdef _USE_INTEGERS
31	molod	1.4	parameter (mbig=1000000000,mz=0,fac=1./mbig)
32			data mseed/161803398/
33	adcroft	1.2	#else
34	molod	1.4	parameter (mbig=4000000.,mz=0.,fac=1./mbig)
35			data mseed/1618033./
36	adcroft	1.2	#endif
37			integer i,ii,inext,inextp,k,idum
38			parameter(idum=-2)
39			#ifdef _USE_INTEGERS
40			integer mj,mk,ma(nff)
41			#else
42			real*8 mj,mk,ma(nff)
43			#endif
44			logical firstCall
45			save firstCall,inext,inextp,ma
46			data firstCall /.true./
47	edhill	1.5
48	molod	1.4	mseed = seed
49	adcroft	1.2	if(firstCall)then
50			firstCall=.false.
51			mj=mseed-iabs(idum)
52			mj=mod(mj,mbig)
53			ma(nff)=mj
54			mk=1
55			do i=1,nff-1
56			ii=mod(21*i,nff)
57			ma(ii)=mk
58			mk=mj-mk
59			if(mk.lt.mz)mk=mk+mbig
60			mj=ma(ii)
61			enddo
62			do k=1,4
63			do i=1,nff
64			ma(i)=ma(i)-ma(1+mod(i+30,nff))
65			if(ma(i).lt.MZ)ma(i)=ma(i)+mbig
66			enddo
67			enddo
68			inext=0
69			inextp=31
70			endif
71			inext=mod(inext,nff)+1
72			inextp=mod(inextp,nff)+1
73			mj=ma(inext)-ma(inextp)
74			if(mj.lt.MZ)mj=mj+mbig
75			ma(inext)=mj
76			port_rand=mj*fac
77			return
78	edhill	1.5
79	adcroft	1.2	end
80
81	edhill	1.5	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
82
83	adcroft	1.2	subroutine port_ranarr(n,arr)
84			implicit none
85			integer n,i
86			real arr(n)
87	molod	1.4	real*8 port_rand,seed
88			seed=1618033.0d0
89	adcroft	1.2	do i=1,n
90	molod	1.4	arr(i)=port_rand(seed)
91	adcroft	1.2	enddo
92
93			return
94			end
95	edhill	1.5
96			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
97			CBOP
98			C !ROUTINE: port_rand_norm
99			C !INTERFACE:
100			real*8 function port_rand_norm()
101
102			C !DESCRIPTION:
103			C This function generates a normally distributed random number with
104			C the so called polar algorithm. The algorithm actually generates 2
105			C numbers, but only 1 is returned for maximum compatibility with old
106			C code. The most obvious way to improve this function would be to
107			C make sure that the second number is not wasted.
108
109			C Changed: 2004.09.06 antti.westerlund@fimr.fi
110
111			C !USES:
112			implicit none
113			CEOP
114
115			C !LOCAL VARIABLES:
116			real*8 port_rand,seed
117			real*8 x1, x2, xs, t
118			integer i
119
120			seed=1618033.0d0
121			C first generate 2 equally distributed random numbers (-1,1)
122			DO WHILE (1 .eq. 1)
123			x1=2.0*port_rand(seed)-1.0
124			x2=2.0*port_rand(seed)-1.0
125			xs=x12+x22
126			if(xs .lt. 1.0 .and. xs .ne. 0.0) then
127			goto 100
128			end if
129			END DO
130			100 continue
131
132			t = sqrt(-2.0*log(xs)/xs)
133			port_rand_norm = t*x1
134			C
135			C also t*x2 would be a gaussian random number and could be returned
136
137			return
138			end