C----------------------------------------------------------------------C
C                                                                      C
C  Lagged Fibonnacci random number generator 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  MUST be initialized with rmarin() before any calls to ranmar().     C
C                                                                      C
C  For the parallel version, initialize with randinit().               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  Note that the number of processors (nprocs) and the processor       C
C  number (procnum) are needed for the subroutine randinit.            C
C  These must be passed as parameters.                                 C
C                                                                      C
C----------------------------------------------------------------------C

	SUBROUTINE randinit(seed,nprocs,procnum)
        IMPLICIT NONE
	INTEGER seed
	INTEGER procnum, nprocs
	INTEGER MULT,ADD,M 
	PARAMETER ( MULT = 23159 )
	PARAMETER ( ADD	 = 6779 )
	PARAMETER ( M = 2147483647 )    ! 2**31 - 1
	INTEGER AAA,BBB,randx
	INTEGER proc
C
C ----  Initialize seeds using a linear congruential random number 
C ----  generator, so that every processor gets a different seed.
C
	IF (seed .LT. 0) seed = - seed
	seed = MOD(seed, M)
 
	AAA = 1
	BBB = 0
	DO proc=0,nprocs-1
	  AAA = MOD(MULT * AAA, M)
	  BBB = MOD(MULT * BBB + ADD, M)
	  IF (proc .EQ. procnum) randx = MOD(AAA*seed + BBB, M)
	ENDDO

C
C ----  Initialize ranmar.
C
	CALL rmarin(randx)

	RETURN
	END



C----------------------------------------------------------------------C
C                                                                      C
C  This is the initialization routine RMARIN for the random number     C
C  generator RANMAR.                                                   C
C  It uses a small lag Fibonacci generator to set up the lag table.    C
C                                                                      C
C  NOTE: The seed variables can have values between: 0 <= ij <= 31328  C
C                                                    0 <= kl <= 30081  C
C                                                                      C
C----------------------------------------------------------------------C

      SUBROUTINE rmarin(seed)
      IMPLICIT NONE
      INTEGER seed
      INTEGER MaxSeed
      PARAMETER (MaxSeed = 900000000)
      INTEGER i,j,k,l, ij,kl, ii,jj, m
      REAL s,t
      REAL u(97), c, cd, cm
      INTEGER i97, j97, ivec
      COMMON /raset1/ u, c, cd, cm, i97, j97, ivec

      seed  = MOD(seed, MaxSeed) 
      IF (seed.EQ.0) seed = 1
      ij = seed / 30082
      kl = seed - 30082 * ij

      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 ranmar()
      IMPLICIT NONE
      REAL u(97), c, cd, cm, uni, ranmar
      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
      ranmar = uni
      RETURN
      END


C----------------------------------------------------------------------C
C                                                                      C
C   Returns a random integer between 1 and max inclusive               C
C                                                                      C
C----------------------------------------------------------------------C
 
        integer function ran(max)
        integer max
        real ranmar
        
        ran=(ranmar()*max)+1
        return
        end

