
C----------------------------------------------------------------------C
C                                                                      C
C  Parallel version of the lagged Fibonnacci random number generator   C
C  RANMAR.                                                             C
C                                                                      C
C  See F. James, Comp. Phys. Comm. 60, 329 (1990), or                  C
C  G. Marsaglia et al., Stat. Prob. Lett. 9, 35 (1990).                C
C                                                                      C
C  Need to initialize the seeds for every processor differently, so    C
C  take a single seed and generate independent numbers on each         C
C  processor using a DIFFERENT random number generator, in this case,  C
C  the parallel linear congruential generator prand (see Fox et al.,   C
C  Solving Problems on Concurrent Processors).                         C
C                                                                      C
C  MUST be initialized with randinit() before any calls to rand().     C
C                                                                      C
C----------------------------------------------------------------------C

	SUBROUTINE randinit(seed)
	INCLUDE "common.h"
	INTEGER seed
	INTEGER MULT,ADD 
	INTEGER AAA,BBB,randx
	PARAMETER ( MULT = 23159 )
	PARAMETER ( ADD	 = 6779 )
	INTEGER irand,seed1,seed2,proc
	COMMON /ranset/ AAA,BBB,randx

C ----  Initialize seeds for the random number generator irand, 
C ----  so that every processor gets a different seed.
	AAA = 1
	BBB = 0

	seed = MOD(seed, MASK)

	DO proc=0,nprocs-1
	  AAA = MOD(MULT * AAA, MASK)
	  BBB = MOD(MULT * BBB + ADD, MASK)
	  IF (proc .EQ. procnum) randx = MOD(AAA*seed + BBB, MASK)
	ENDDO

C ----  Initialize the random number generator ranmar using irand to
C ----  generate the seeds, so that every processor gets different seeds.
	seed1 = irand()
	seed2 = irand()
	CALL rmarin(seed1,seed2)

	RETURN
	END




C----------------------------------------------------------------------C
C                                                                      C
C  Linear congruential generator to return a random int in [0, 2^15).  C
C                                                                      C
C----------------------------------------------------------------------C

	FUNCTION irand()
	INCLUDE "common.h"
	INTEGER irand
	INTEGER AAA,BBB,randx
	COMMON /ranset/ AAA,BBB,randx
	randx = MOD(AAA*randx+BBB, MASK)
	irand = randx
	RETURN
	END



C----------------------------------------------------------------------C
C                                                                      C
C This is the initialization routine RMARIN for the random number      C
C     generator RANMAR                                                 C
C                                                                      C
C NOTE: The seed variables can have values between:  0 <= IJ <= 31328  C
C                                                    0 <= KL <= 30081  C
C----------------------------------------------------------------------C

      SUBROUTINE rmarin(ij,kl)
      REAL u(97), c, cd, cm
      INTEGER i97, j97, ivec
      COMMON /raset1/ u, c, cd, cm, i97, j97, ivec
      IF( ij .LT. 0  .OR.  ij .GT. 31328 ) THEN
        ij = MOD(ij,31328+1)
      ENDIF
      IF( kl .LT. 0  .OR.  kl .GT. 30081 ) THEN
        kl = MOD(kl,30081+1)
      ENDIF
      i = mod(ij/177, 177) + 2
      j = mod(ij    , 177) + 2
      k = mod(kl/169, 178) + 1
      l = mod(kl,     169)
      DO 2 ii = 1, 97
        s = 0.0
        t = 0.5
        DO 3 jj = 1, 24
          m = mod(mod(i*j, 179)*k, 179)
          i = j
          j = k
          k = m
          l = mod(53*l+1, 169)
          IF (mod(l*m, 64) .ge. 32) then
            s = s + t
          ENDIF
          t = 0.5 * t
    3   CONTINUE
        u(ii) = s
    2 CONTINUE
      c = 362436.0 / 16777216.0
      cd = 7654321.0 / 16777216.0
      cm = 16777213.0 /16777216.0
      i97 = 97
      j97 = 33
      RETURN
      END



C----------------------------------------------------------------------C
C                                                                      C
C  Lagged Fibonacci random number generator RANMAR().                  C
C                                                                      C
C----------------------------------------------------------------------C

      FUNCTION rand()
      IMPLICIT NONE
      REAL u(97), c, cd, cm, uni, rand
      INTEGER i97, j97, ivec
      COMMON /raset1/ u, c, cd, cm, i97, j97, ivec

      uni = u(i97) - u(j97)
      IF( uni .LT. 0.0 ) uni = uni + 1.0
      u(i97) = uni
      i97 = i97 - 1
      IF(i97 .EQ. 0) i97 = 97
      j97 = j97 - 1
      IF(j97 .EQ. 0) j97 = 97
      c = c - cd
      IF( c .LT. 0.0 ) c = c + cm
      uni = uni - c
      IF( uni .LT. 0.0 ) uni = uni + 1.0
      rand = uni
      RETURN
      END

