/[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.4 by molod, Mon Jul 19 15:13:07 2004 UTC revision 1.5 by edhill, Tue Sep 7 15:16:58 2004 UTC
# Line 5  C Portable random number generator Line 5  C Portable random number generator
5    
6  #undef _USE_INTEGERS  #undef _USE_INTEGERS
7    
8  ! ==============================================================================  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
9  CBOP  CBOP
10  C     !ROUTINE: port_rand  C     !ROUTINE: port_rand
11  C     !INTERFACE:  C     !INTERFACE:
12        real*8 function port_rand(seed)        real*8 function port_rand(seed)
13    
14  C     !DESCRIPTION: \bv  C     !DESCRIPTION:
15  C     Portable random number generator  C     Portable random number generator
16  C     \ev        
17    C     !USES:
18          implicit none
19    CEOP
20    
21  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
       implicit none  
22        integer nff        integer nff
23        parameter(nff=55)        parameter(nff=55)
24  #ifdef _USE_INTEGERS  #ifdef _USE_INTEGERS
# Line 42  C     !LOCAL VARIABLES: Line 44  C     !LOCAL VARIABLES:
44        logical firstCall        logical firstCall
45        save firstCall,inext,inextp,ma        save firstCall,inext,inextp,ma
46        data firstCall /.true./        data firstCall /.true./
47  CEOP  
 ! ------------------------------------------------------------------------------  
48        mseed = seed        mseed = seed
49        if(firstCall)then        if(firstCall)then
50          firstCall=.false.          firstCall=.false.
# Line 74  CEOP Line 75  CEOP
75        ma(inext)=mj        ma(inext)=mj
76        port_rand=mj*fac        port_rand=mj*fac
77        return        return
78  ! ------------------------------------------------------------------------------  
79        end        end
 ! ==============================================================================  
80    
81  ! ==============================================================================  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
82    
83        subroutine port_ranarr(n,arr)        subroutine port_ranarr(n,arr)
84        implicit none        implicit none
85        integer n,i        integer n,i
86        real arr(n)        real arr(n)
87        real*8 port_rand,seed        real*8 port_rand,seed
88        seed=1618033.0d0        seed=1618033.0d0
 ! ------------------------------------------------------------------------------  
89        do i=1,n        do i=1,n
90         arr(i)=port_rand(seed)         arr(i)=port_rand(seed)
91        enddo        enddo
92    
93        return        return
 ! ------------------------------------------------------------------------------  
94        end        end
95  ! ==============================================================================  
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=x1**2+x2**2
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
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22