/[MITgcm]/MITgcm/model/src/port_rand.F
ViewVC logotype

Diff of /MITgcm/model/src/port_rand.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1.2.1 by adcroft, Thu Jan 25 18:44:07 2001 UTC revision 1.8 by jmc, Fri Aug 9 15:00:46 2013 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
 C Portable random number generator  
   
4  #undef _USE_INTEGERS  #undef _USE_INTEGERS
5    
6  ! ==============================================================================  C--  File port_rand.F: Portable random number generator functions
7        real*8 function port_rand()  C--   Contents
8        implicit none  C--   o PORT_RAND
9        integer nff  C--   o PORT_RANARR
10        parameter(nff=55)  C--   o PORT_RAND_NORM
11    
12    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
13    CBOP
14    C     !ROUTINE: PORT_RAND
15    C     !INTERFACE:
16          real*8 FUNCTION PORT_RAND(seed)
17    
18    C     !DESCRIPTION:
19    C     Portable random number generator
20    C      seed >=0 :: initialise using this seed ; and return 0
21    C      seed < 0 :: if first call then initialise using the default seed (=mseed)
22    C                  and always return a random number
23    
24    C     !USES:
25          IMPLICIT NONE
26    
27    C     !INPUT PARAMETERS:
28  #ifdef _USE_INTEGERS  #ifdef _USE_INTEGERS
29        integer mbig,mseed,mZ        INTEGER seed
30  #else  #else
31        real*8 mbig,mseed,mz        real*8  seed
32  #endif  #endif
33    CEOP
34    
35    C     !LOCAL VARIABLES:
36          INTEGER nff,idum
37          PARAMETER(nff=55)
38          PARAMETER(idum=-2)
39        real*8 fac        real*8 fac
40  #ifdef _USE_INTEGERS  #ifdef _USE_INTEGERS
41        parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1./mbig)        INTEGER mbig,mseed,mZ
42  #else        PARAMETER (mbig=1000000000,mz=0,fac=1.d0/mbig)
43        parameter (mbig=4000000.,mseed=1618033.,mz=0.,fac=1./mbig)        INTEGER mj,mk,ma(nff)
44  #endif        DATA mseed/161803398/
       integer i,ii,inext,inextp,k,idum  
       parameter(idum=-2)  
 #ifdef _USE_INTEGERS  
       integer mj,mk,ma(nff)  
45  #else  #else
46          real*8 mbig,mseed,mz
47          PARAMETER (mbig=4000000.,mz=0.,fac=1.d0/mbig)
48        real*8 mj,mk,ma(nff)        real*8 mj,mk,ma(nff)
49          DATA mseed/1618033./
50  #endif  #endif
51        logical firstCall        LOGICAL firstCall
52        save firstCall,inext,inextp,ma        INTEGER i,ii,inext,inextp,k
53        data firstCall /.true./        DATA firstCall /.true./
54  ! ------------------------------------------------------------------------------        SAVE firstCall,inext,inextp,ma
55        if(firstCall)then  
56    C-    Initialise the random number generator
57          IF (firstCall .OR. seed.GE.mz) THEN
58            IF (seed.GE.mz) mseed = seed
59          firstCall=.false.          firstCall=.false.
60          mj=mseed-iabs(idum)          mj=mseed-iabs(idum)
61          mj=mod(mj,mbig)          mj=mod(mj,mbig)
62          ma(nff)=mj          ma(nff)=mj
63          mk=1          mk=1
64          do i=1,nff-1          DO i=1,nff-1
65            ii=mod(21*i,nff)            ii=mod(21*i,nff)
66            ma(ii)=mk            ma(ii)=mk
67            mk=mj-mk            mk=mj-mk
68            if(mk.lt.mz)mk=mk+mbig            IF (mk.LT.mz) mk=mk+mbig
69            mj=ma(ii)            mj=ma(ii)
70          enddo          ENDDO
71          do k=1,4          DO k=1,4
72            do i=1,nff            DO i=1,nff
73              ma(i)=ma(i)-ma(1+mod(i+30,nff))              ma(i)=ma(i)-ma(1+mod(i+30,nff))
74              if(ma(i).lt.MZ)ma(i)=ma(i)+mbig              IF (ma(i).LT.mz) ma(i)=ma(i)+mbig
75            enddo            ENDDO
76          enddo          ENDDO
77          inext=0          inext=0
78          inextp=31          inextp=31
79        endif        ENDIF
80        inext=mod(inext,nff)+1  
81        inextp=mod(inextp,nff)+1  C-    Compute a random number (only if seed < 0)
82        mj=ma(inext)-ma(inextp)        IF (seed.GE.mz) THEN
83        if(mj.lt.MZ)mj=mj+mbig          port_rand=0.d0
84        ma(inext)=mj        ELSE
85        port_rand=mj*fac          inext=mod(inext,nff)+1
86        return          inextp=mod(inextp,nff)+1
87  ! ------------------------------------------------------------------------------          mj=ma(inext)-ma(inextp)
88        end          IF (mj.LT.mz) mj=mj+mbig
89  ! ==============================================================================          ma(inext)=mj
90            port_rand=mj*fac
91  ! ==============================================================================        ENDIF
92        subroutine port_ranarr(n,arr)  
93        implicit none        RETURN
94        integer n,i        END
95        real arr(n)  
96    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
97    
98    CBOP
99    C     !ROUTINE: PORT_RANARR
100    C     !INTERFACE:
101          SUBROUTINE PORT_RANARR(n,arr)
102    
103    C     !DESCRIPTION:
104    C     Portable random number array generator
105    
106    C     !USES:
107          IMPLICIT NONE
108    
109    C     !INPUT PARAMETERS:
110          INTEGER n
111          real*8 arr(n)
112    CEOP
113    
114    C     !LOCAL VARIABLES:
115          INTEGER i
116        real*8 port_rand        real*8 port_rand
117  ! ------------------------------------------------------------------------------  #ifdef _USE_INTEGERS
118        do i=1,n        INTEGER seed
119         arr(i)=port_rand()        seed=-1
120        enddo  #else
121          real*8  seed
122        return        seed=-1.d0
123  ! ------------------------------------------------------------------------------  #endif
124        end  c     seed=1618033.0d0
125  ! ==============================================================================        DO i=1,n
126           arr(i)=port_rand(seed)
127          ENDDO
128    
129          RETURN
130          END
131    
132    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
133    CBOP
134    C     !ROUTINE: PORT_RAND_NORM
135    C     !INTERFACE:
136          real*8 FUNCTION PORT_RAND_NORM()
137    
138    C     !DESCRIPTION:
139    C     This function generates a normally distributed random number with
140    C     the so called polar algorithm. The algorithm actually generates 2
141    C     numbers, but only 1 is returned for maximum compatibility with old
142    C     code.  The most obvious way to improve this function would be to
143    C     make sure that the second number is not wasted.
144    
145    C     Changed: 2004.09.06 antti.westerlund@fimr.fi
146    
147    C     !USES:
148          IMPLICIT NONE
149    CEOP
150    
151    C     !LOCAL VARIABLES:
152          real*8 port_rand
153          real*8 x1, x2, xs, t
154    
155    #ifdef _USE_INTEGERS
156          INTEGER seed
157          seed=-1
158    #else
159          real*8  seed
160          seed=-1.d0
161    #endif
162    c     seed=1618033.0d0
163    
164    C     first generate 2 equally distributed random numbers (-1,1)
165          xs = 0.0
166          DO WHILE ( xs.GE.1.0 .OR. xs.EQ.0.0 )
167             x1=2.0*port_rand(seed)-1.0
168             x2=2.0*port_rand(seed)-1.0
169             xs=x1**2+x2**2
170          ENDDO
171    
172          t = SQRT(-2.0*LOG(xs)/xs)
173          port_rand_norm = t*x1
174    
175    C     also t*x2 would be a gaussian random number and could be returned
176    
177          RETURN
178          END
Legend:
Removed from v.1.1.2.1  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22