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

      double precision function drnormal(ix)

C     A Box-Muller transform (by George Edward Pelham Box and Mervin 
C     Edgar Muller 1958) is a method of generating pairs of independent 
C     standard normally distributed (zero expectation, unit variance) 
C     random numbers, given a source of uniformly distributed random 
C     numbers.

C     SCALAR ARGUMENTS
      double precision ix

C     LOCAL SCALARS
      double precision t1,t2

C     PARAMETERS
      double precision two,twopi
      parameter ( two   = 2.0d+0                           )
      parameter ( twopi = 6.2831853071795864769252867663d0 )

C     EXTERNAL FUNCTIONS
      double precision drand
      external drand

      t1 = drand(ix)

      t2 = drand(ix)

      drnormal = sqrt( - two * log( t1 ) ) * cos( twopi * t2 )

      return

      end

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

      double precision function drand(ix)

C     This is the random number generator of Schrage:
C
C     L. Schrage, A more portable Fortran random number generator, ACM
C     Transactions on Mathematical Software 5 (1979), 132-138.

      double precision ix

      double precision a,p,b15,b16,xhi,xalo,leftlo,fhi,k
      data a/16807.d0/,b15/32768.d0/,b16/65536.d0/,p/2147483647.d0/

      xhi= ix/b16
      xhi= xhi - dmod(xhi,1.d0)
      xalo= (ix-xhi*b16)*a
      leftlo= xalo/b16
      leftlo= leftlo - dmod(leftlo,1.d0)
      fhi= xhi*a + leftlo
      k= fhi/b15
      k= k - dmod(k,1.d0)
      ix= (((xalo-leftlo*b16)-p)+(fhi-k*b15)*b16)+k
      if (ix.lt.0) ix= ix + p
      drand= ix*4.656612875d-10

      return

      end
