C
C  Lagged Fibonacci random number generator RANMAR().
C  Must be initialized with RMARIN() before use.
C
C  See F. James, Comp. Phys. Comm. 60, 329 (1990), or
C  G. Marsaglia, A. Zaman and W.-W. Tsang, Stat. Prob. Lett. 9, 35 (1990).
C

C
C  Need to initialize the seeds for every processor differently, so take 
C  a single seed and generate independent numbers on each processor using
C  a DIFFERENT random number generator, in this case, prand.
C

	SUBROUTINE rand_init(seed)
	INCLUDE "exp_wolff.h"
	INTEGER seed
	INTEGER MULT,ADD 
	PARAMETER ( MULT = 23159 )
	PARAMETER ( ADD	 = 6779 )
	INTEGER AAA,BBB
	INTEGER irand,seed1,seed2,proc,randx
	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,nproc-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 Return a random int in [0, 2^15).
C
	FUNCTION irand()
	INCLUDE "exp_wolff.h"
	INTEGER irand
	INTEGER AAA,BBB,randx
	COMMON /ranset/ AAA,BBB,randx
	randx = MOD(AAA*randx+BBB, MASK)
	irand = randx
	RETURN
	END


C **************************************************************************


      SUBROUTINE rmarin(ij,kl)
C This is the initialization routine for the random number generator
C     RANMAR()
C NOTE: The seed variables can have values between:    0 <= IJ <= 31328
C                                                      0 <= KL <= 30081
      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  Lagged Fibonacci random number generator RANMAR().

      FUNCTION ranmar()
      IMPLICIT NONE
      REAL rand_u(97), rand_c, rand_cd, rand_cm, rand_uni, ranmar
      INTEGER rand_i97, rand_j97, rand_ivec
      COMMON /raset1/ rand_u, rand_c, rand_cd, rand_cm, rand_i97,
     $                rand_j97, rand_ivec 

      rand_uni = rand_u(rand_i97) - rand_u(rand_j97)
      IF( rand_uni .LT. 0.0 ) rand_uni = rand_uni + 1.0
      rand_u(rand_i97) = rand_uni
      rand_i97 = rand_i97 - 1
      IF(rand_i97 .EQ. 0) rand_i97 = 97
      rand_j97 = rand_j97 - 1
      IF(rand_j97 .EQ. 0) rand_j97 = 97
      rand_c = rand_c - rand_cd
      IF( rand_c .LT. 0.0 ) rand_c = rand_c + rand_cm
      rand_uni = rand_uni - rand_c
      IF( rand_uni .LT. 0.0 ) rand_uni = rand_uni + 1.0
      ranmar = rand_uni
      RETURN
      END




