model/src/port_rand.F

fb481a83c2 Alis*0001 #undef _USE_INTEGERS
                0002
f04f2001af Jean*0003 C--  File port_rand.F: Portable random number generator functions
                0004 C--   Contents
                0005 C--   o PORT_RAND
                0006 C--   o PORT_RANARR
                0007 C--   o PORT_RAND_NORM
                0008
019d02e9b2 Ed H*0009 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
9366854e02 Chri*0010 CBOP
f04f2001af Jean*0011 C     !ROUTINE: PORT_RAND
9366854e02 Chri*0012 C     !INTERFACE:
f04f2001af Jean*0013       real*8 FUNCTION PORT_RAND(seed)
9366854e02 Chri*0014
019d02e9b2 Ed H*0015 C     !DESCRIPTION:
9366854e02 Chri*0016 C     Portable random number generator
29a56a7b59 Jean*0017 C      seed >=0 :: initialise using this seed ; and return 0
                0018 C      seed < 0 :: if first call then initialise using the default seed (=mseed)
                0019 C                  and always return a random number
f04f2001af Jean*0020
019d02e9b2 Ed H*0021 C     !USES:
f04f2001af Jean*0022       IMPLICIT NONE
                0023
29a56a7b59 Jean*0024 C     !INPUT PARAMETERS:
                0025 #ifdef _USE_INTEGERS
f04f2001af Jean*0026       INTEGER seed
29a56a7b59 Jean*0027 #else
                0028       real*8  seed
                0029 #endif
019d02e9b2 Ed H*0030 CEOP
9366854e02 Chri*0031
                0032 C     !LOCAL VARIABLES:
f04f2001af Jean*0033       INTEGER nff,idum
                0034       PARAMETER(nff=55)
                0035       PARAMETER(idum=-2)
29a56a7b59 Jean*0036       real*8 fac
fb481a83c2 Alis*0037 #ifdef _USE_INTEGERS
f04f2001af Jean*0038       INTEGER mbig,mseed,mZ
                0039       PARAMETER (mbig=1000000000,mz=0,fac=1.d0/mbig)
                0040       INTEGER mj,mk,ma(nff)
                0041       DATA mseed/161803398/
fb481a83c2 Alis*0042 #else
29a56a7b59 Jean*0043       real*8 mbig,mseed,mz
f04f2001af Jean*0044       PARAMETER (mbig=4000000.,mz=0.,fac=1.d0/mbig)
fb481a83c2 Alis*0045       real*8 mj,mk,ma(nff)
f04f2001af Jean*0046       DATA mseed/1618033./
fb481a83c2 Alis*0047 #endif
f04f2001af Jean*0048       LOGICAL firstCall
                0049       INTEGER i,ii,inext,inextp,k
                0050       DATA firstCall /.true./
                0051       SAVE firstCall,inext,inextp,ma
019d02e9b2 Ed H*0052
29a56a7b59 Jean*0053 C-    Initialise the random number generator
f04f2001af Jean*0054       IF (firstCall .OR. seed.GE.mz) THEN
                0055         IF (seed.GE.mz) mseed = seed
fb481a83c2 Alis*0056         firstCall=.false.
                0057         mj=mseed-iabs(idum)
                0058         mj=mod(mj,mbig)
                0059         ma(nff)=mj
                0060         mk=1
f04f2001af Jean*0061         DO i=1,nff-1
fb481a83c2 Alis*0062           ii=mod(21*i,nff)
                0063           ma(ii)=mk
                0064           mk=mj-mk
f04f2001af Jean*0065           IF (mk.LT.mz) mk=mk+mbig
fb481a83c2 Alis*0066           mj=ma(ii)
f04f2001af Jean*0067         ENDDO
                0068         DO k=1,4
                0069           DO i=1,nff
fb481a83c2 Alis*0070             ma(i)=ma(i)-ma(1+mod(i+30,nff))
f04f2001af Jean*0071             IF (ma(i).LT.mz) ma(i)=ma(i)+mbig
                0072           ENDDO
                0073         ENDDO
fb481a83c2 Alis*0074         inext=0
                0075         inextp=31
f04f2001af Jean*0076       ENDIF
019d02e9b2 Ed H*0077
29a56a7b59 Jean*0078 C-    Compute a random number (only if seed < 0)
f04f2001af Jean*0079       IF (seed.GE.mz) THEN
29a56a7b59 Jean*0080         port_rand=0.d0
f04f2001af Jean*0081       ELSE
29a56a7b59 Jean*0082         inext=mod(inext,nff)+1
                0083         inextp=mod(inextp,nff)+1
                0084         mj=ma(inext)-ma(inextp)
f04f2001af Jean*0085         IF (mj.LT.mz) mj=mj+mbig
29a56a7b59 Jean*0086         ma(inext)=mj
                0087         port_rand=mj*fac
f04f2001af Jean*0088       ENDIF
29a56a7b59 Jean*0089
f04f2001af Jean*0090       RETURN
                0091       END
fb481a83c2 Alis*0092
019d02e9b2 Ed H*0093 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0094
f04f2001af Jean*0095 CBOP
                0096 C     !ROUTINE: PORT_RANARR
                0097 C     !INTERFACE:
                0098       SUBROUTINE PORT_RANARR(n,arr)
                0099
                0100 C     !DESCRIPTION:
                0101 C     Portable random number array generator
                0102
                0103 C     !USES:
                0104       IMPLICIT NONE
                0105
                0106 C     !INPUT PARAMETERS:
                0107       INTEGER n
                0108       real*8 arr(n)
                0109 CEOP
                0110
                0111 C     !LOCAL VARIABLES:
                0112       INTEGER i
29a56a7b59 Jean*0113       real*8 port_rand
                0114 #ifdef _USE_INTEGERS
f04f2001af Jean*0115       INTEGER seed
29a56a7b59 Jean*0116       seed=-1
                0117 #else
                0118       real*8  seed
                0119       seed=-1.d0
                0120 #endif
                0121 c     seed=1618033.0d0
f04f2001af Jean*0122       DO i=1,n
3cd17a95e0 Andr*0123        arr(i)=port_rand(seed)
f04f2001af Jean*0124       ENDDO
fb481a83c2 Alis*0125
f04f2001af Jean*0126       RETURN
                0127       END
019d02e9b2 Ed H*0128
                0129 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0130 CBOP
f04f2001af Jean*0131 C     !ROUTINE: PORT_RAND_NORM
019d02e9b2 Ed H*0132 C     !INTERFACE:
f04f2001af Jean*0133       real*8 FUNCTION PORT_RAND_NORM()
019d02e9b2 Ed H*0134
                0135 C     !DESCRIPTION:
                0136 C     This function generates a normally distributed random number with
                0137 C     the so called polar algorithm. The algorithm actually generates 2
                0138 C     numbers, but only 1 is returned for maximum compatibility with old
                0139 C     code.  The most obvious way to improve this function would be to
                0140 C     make sure that the second number is not wasted.
                0141
                0142 C     Changed: 2004.09.06 antti.westerlund@fimr.fi
                0143
                0144 C     !USES:
f04f2001af Jean*0145       IMPLICIT NONE
019d02e9b2 Ed H*0146 CEOP
                0147
                0148 C     !LOCAL VARIABLES:
29a56a7b59 Jean*0149       real*8 port_rand
f04f2001af Jean*0150       real*8 x1, x2, xs, t
019d02e9b2 Ed H*0151
29a56a7b59 Jean*0152 #ifdef _USE_INTEGERS
f04f2001af Jean*0153       INTEGER seed
29a56a7b59 Jean*0154       seed=-1
                0155 #else
                0156       real*8  seed
                0157       seed=-1.d0
                0158 #endif
                0159 c     seed=1618033.0d0
                0160
f04f2001af Jean*0161 C     first generate 2 equally distributed random numbers (-1,1)
231007c5f0 Jean*0162       xs = 0.0
                0163       DO WHILE ( xs.GE.1.0 .OR. xs.EQ.0.0 )
019d02e9b2 Ed H*0164          x1=2.0*port_rand(seed)-1.0
                0165          x2=2.0*port_rand(seed)-1.0
                0166          xs=x1**2+x2**2
f04f2001af Jean*0167       ENDDO
019d02e9b2 Ed H*0168
f04f2001af Jean*0169       t = SQRT(-2.0*LOG(xs)/xs)
019d02e9b2 Ed H*0170       port_rand_norm = t*x1
231007c5f0 Jean*0171
019d02e9b2 Ed H*0172 C     also t*x2 would be a gaussian random number and could be returned
                0173
f04f2001af Jean*0174       RETURN
                0175       END