      MODULE numbers
      !storage of current generator values
      INTEGER, dimension(:), allocatable:: randx, idum2
      REAL, dimension(:), allocatable:: storage
      INTEGER :: pnodes
      END MODULE numbers
      

      MODULE global_constants
      
      !Floating Point Number initialization

      SAVE
      INTEGER, PARAMETER :: IM1 = 2147483563
      INTEGER, PARAMETER :: IM2 = 2147483399
      INTEGER, PARAMETER :: IMM1 = 2147483562
      INTEGER, PARAMETER :: IA1 = 40014
      INTEGER, PARAMETER :: IA2 = 40692
      INTEGER, PARAMETER :: IQ1 = 53668
      INTEGER, PARAMETER :: IQ2 = 52774
      INTEGER, PARAMETER :: IR1 = 12211
      INTEGER, PARAMETER :: IR2 = 3791
      INTEGER, PARAMETER :: MAXSEED = 900000000
!Integer Constants Set

      REAL, PARAMETER :: EPS = 1.2E-7
      REAL, PARAMETER :: RNMX = (1.0 - EPS)
      REAL, PARAMETER :: AM = (1.0 / IM1)
!Real Constants Set

      END MODULE global_constants


  

!Subroutine for lecuyer's method.
      SUBROUTINE lecuyer(array)
      USE numbers
      !set initial seeds.
      REAL array(:)
      INTEGER :: g,j
!HPF$ INDEPENDENT
      DO g = 1, pnodes
         array(g) = prand(g)
       END DO
      end subroutine lecuyer
      
           
      INTEGER FUNCTION array_size()
      
      USE numbers
      array_size = pnodes
      
      END FUNCTION array_size

      SUBROUTINE dump_seeds(array)
      
      USE numbers
      INTEGER :: array(pnodes,2)
      INTEGER :: i

      DO i = 1, pnodes
         array(i,1) = randx(i) 
         array(i,2) = idum2(i)
        
      END DO
      print *,array
      END SUBROUTINE dump_seeds

      SUBROUTINE put_seeds(array)
      
      USE numbers
      INTEGER :: array(pnodes,2)
      INTEGER :: i
      
      DO i = 1, pnodes
         randx(i) = array(i,1)
         idum2(i) = array(i,2)
      END DO
      END SUBROUTINE put_seeds

      REAL FUNCTION prand(g)
      USE global_constants
      USE numbers
      IMPLICIT NONE 
      INTEGER :: g, iy
      REAL :: temp, k

!set initial variables.
      k = (randx(g) / IQ1)
      randx(g) = IA1*(randx(g) - k * IQ1) -k * IR1
      IF (randx(g) < 0) THEN 
         randx(g) = randx(g) + IM1
      END IF
 
      k = (idum2(g) / IQ2)
      idum2(g) = IA2*(idum2(g) -k * IQ2) -k * IR2
      IF (idum2(g) < 0) THEN 
         idum2(g) = idum2(g) + IM2
      END IF

      iy = randx(g) - idum2(g)
      IF (iy < 1) THEN
         iy = iy + IMM1
      END IF
      temp = AM * iy
      prand = temp
      
      END FUNCTION prand

      SUBROUTINE set_init_seed(seed, nodes)
      USE numbers
      INTEGER :: seed,nodes,i,x
!HPF$ PROCESSORS p(NUMBER_OF_PROCESSORS())
!HPF$ ALIGN randx(:) with idum2(:)
!HPF$ DISTRIBUTE  (BLOCK) ONTO p:: idum2
      pnodes = nodes
      allocate(randx(pnodes))
      allocate(idum2(pnodes))

      randx(1) = seed
      idum2(1) = seed
      DO i = 2, pnodes
         !Each Partition has 8.58 x 10^9 numbers in it.
         call advance_state(i, 1)
      END DO

      END SUBROUTINE set_init_seed
      INTEGER FUNCTION seeder(seed)
      USE global_constants
      INTEGER :: ij, kl, i, j, k, l, ii, jj, m, s, t, seed

      
      if (seed .lt. 0) then
         seed = -seed
      end if
      seed = MOD(seed, MAXSEED)
      t = seed
      ij = seed / 30082
      kl = seed - (30082 * ij)
      i = ((ij / 177)/177) + 2
      j = (ij / 177) + 2
      k = ((kl/169) / 178) + 1
      l = kl / 169

      do ii = 1, 2
         s = 0
         do jj = 1, 32
            m = (((i * j) / 179) * k) / 179
            i = j
            j = k
            k = m
            l = (53*l+1) / 169
            if (((l*m)/64) .le. 32) then
               s = s + t
               t = 2 * t
            end if
         end do
      end do
      seeder = t
      end function seeder


      SUBROUTINE set_seed(seed1,seed2,g)

      USE numbers
      INTEGER :: seed1, seed2, g
      randx(g) = seed1
      idum2(g) = seed2
      END SUBROUTINE set_seed


      INTEGER FUNCTION multimod(a,s,m)
      !from L'Ecuyer and Cote'
      INTEGER, PARAMETER :: h = 32768
      !partition size (2^d)

      INTEGER :: a,s,m,q, r, k, p, sign
      
      q = (m / a)
      r = m - (a * q)
      k = s / q
      p = a*(s - (k * q))
     
      sign = -1
      IF (a .lt. h) THEN
         IF (r .ne. 0) THEN
            IF (k .ne. 0) THEN
              DO
                 a = r
                 s = k
                 q = m / a
                 r = m - (a * q)
                 k = s / q
                 p = p + sign*a*(s - (k * q))
                 sign = -sign
                 IF ((sign * p) .lt. 0) THEN
                    p = p - (sign*m)
                 END IF
                 IF (a .lt. h) THEN
                    IF (r .ne. 0) THEN
                       IF (k .ne. 0) THEN
                          EXIT
                       END IF
                    END IF
                 END IF
              END DO
           END IF
        END IF
      END IF
      IF (k .ne. 0) THEN
         p = p + (sign*k*r)
      END IF
      IF (p .lt. 0) THEN
         DO
            p = p + m
            IF (p .gt. 0) EXIT
         END DO 
      END IF
      IF (p .eq. 0) THEN
         p = 1
      END IF
      multimod = p

      END FUNCTION multimod


      subroutine advance_state(g,k)
      
      USE numbers
      USE global_constants

      integer g,k,b1,b2,i,j
      b1 = randx(g - 1)
      b2 = idum2(g - 1)
      b1 = multimod(IA1, b1, IM1)
      b2 = multimod(IA2, b2, IM2)
      do i  = 1, 64      
         b1 = multimod(b1, b1, IM1)
         b2 = multimod(b1, b2, IM2)
      end do
      
      call  set_seed(b1,b2,g)
      end subroutine advance_state


      PROGRAM rand_gen
!      subroutine random
      USE numbers
      INTERFACE
         SUBROUTINE lecuyer(array)
           real, dimension(:) :: array
         END SUBROUTINE lecuyer
      END INTERFACE
      REAL :: temp(0:7),t1,t2
      INTEGER :: i,j,dummy,dummy2,test(8)


!      print *, ' Test random'
      t1 = dtime(dummy)
      call set_init_seed(1234,8)

      DO i = 1, 4
         call lecuyer(temp)
         do j = 0, 7
!            if (temp(j) .lt. 0.001) then
               print *,temp(j)
!            end if
         end do
         print *, ' '
         print *, ' '

      END DO
      t2 = dtime(dummy)
      END PROGRAM rand_gen
!      end subroutine










