!  *********************************************************************
!  *                                                                   *
!  *                          function randu                           *
!  *                                                                   *
!  *********************************************************************
!  Single Precision Version 2.0
!  Adapted from Press etal. "Numerical Recipes in Fortran", 2nd Ed.,
!  RAN1 program, pg 271
!  by Gordon A. Fenton, TUNS, Oct 14, 1996
!
!  PURPOSE  returns a pseudo-random number uniformly distributed on the
!           interval (0,1)
!
!  This routine uses a pair of multiplicative congruential generators to
!  produce a sequence of random numbers uniformly distributed on the interval
!  (0,1) which exceeds Park and Miller's ``Minimal Standard.'' The generator
!  has a period of about 2.3E+18 and uses a shuffling algorithm due to Bays
!  and Durham to remove low-order serial correlations. This routine is
!  about 2.0 times slower than the minimal routine `RAN0' presented in
!  Numerical Recipes. The endpoints of the interval (0,1) are excluded.
!
!  Arguments to the function are as follows;
!
!      jseed    integer seed. If jseed is zero, the next pseudo random number
!               in the sequence is returned. Otherwise if jseed is positive,
!               it is used to initialize the generator and the first value
!               in the sequence is returned. (input)
!
!  Notes:
!    1) if you always call this routine without initializing it (using
!       a positive seed) then you'll always get the same sequence of
!       pseudo-random numbers.
!
!  REVISION HISTORY:
!  2.0	implemented RAN2 from 2nd Edition of Numerical Recipes
!-----------------------------------------------------------------------------
      real function randu2(jseed)
      integer jseed, idumm,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
      real AM,EPS,RNMX
      PARAMETER (IM1  = 2147483563, &
                IM2  = 2147483399,  &
                AM   = 1./IM1,      &
                IMM1 = IM1-1,       &
                IA1  = 40014,       &
                IA2  = 40692,       &
                IQ1  = 53668,       &
                IQ2  = 52774,       &
                IR1  = 12211,       &
                IR2  = 3791,        & 
                NTAB = 32,          &
                NDIV = 1+IMM1/NTAB, &
                EPS  = 1.2e-7,      &
                RNMX = 1.-EPS)      
      INTEGER idum22,j,k,ivv(NTAB),iyy
      SAVE ivv,iyy,idumm,idum22,ifirst2
      DATA idum22/123456789/, ivv/NTAB*0/, iyy/0/, ifirst2/0/

      if( (jseed .gt. 0) .or. (ifirst2 .eq. 0) ) then
         ifirst2 = 1
         idumm   = max0(jseed,1)
         idum22  = idumm
         do 10 j = NTAB+8, 1, -1
            k    = idumm/IQ1
            idumm = IA1*(idumm-k*IQ1) - k*IR1
            if( idumm .lt. 0 ) idumm  = idumm + IM1
            if( j .le. NTAB ) ivv(j) = idumm
  10     continue
         iyy = ivv(1)
      endif

      k    = idumm/IQ1
      idumm = IA1*(idumm-k*IQ1) - k*IR1
      if( idumm .lt. 0 ) idumm = idumm + IM1
      k     = idum22/IQ2
      idum22 = IA2*(idum22-k*IQ2) - k*IR2
      if( idum22 .lt. 0 ) idum22 = idum22 + IM2
      j     = 1 + iyy/NDIV
      iyy    = ivv(j) - idum22
      ivv(j) = idumm
      if( iyy .lt. 1 ) iyy = iyy + IMM1
      randu2 = amin1( AM*iyy, RNMX )
      return
      end