!     generate gaussian random numbers - Abramovitz & Stegun
      SUBROUTINE GaussRand (nranst,rand1)
      IMPLICIT none
      REAL*8, PARAMETER:: A0 = 2.30753d+00, A1 = 0.27061d+00,  &
                          B1 = 0.99229d+00, B2 = 0.04481d+00
      INTEGER, INTENT(IN):: nranst
      REAL*8, INTENT(OUT):: rand1(nranst)
      REAL*8:: ri,t
      INTEGER:: i
      Call FvRand(rand1,1,nranst)
      DO i=1,nranst
        ri = rand1(i)
        t  = SQRT(-2*LOG(MIN(ri,1.0d0-ri)))
        rand1(i) = SIGN(1.0d0,ri-0.5d0)*(t-(A0+A1*t)/(1.0d0+(B1+B2*t)*t))
      END DO
      RETURN
      END SUBROUTINE GaussRand


!     generate an array of random numbers
      SUBROUTINE FvRand(rnd,incr,nel)
      IMPLICIT none
      INTEGER, INTENT(IN):: incr, nel
      INTEGER:: i
      REAL*8, INTENT(OUT):: rnd(nel)
      REAL*8:: rn
      DO i=1,nel,incr
         rnd(i)=rn()
      END DO
      RETURN
      END SUBROUTINE FvRand


      FUNCTION RN()
      IMPLICIT none
      REAL*8:: RN, ST, U1, U2
      INTEGER:: ISEED
      COMMON/RNDN/ ISEED(64), U1, U2
      ST=U1+U2
      U2=U1
 64   CONTINUE
      IF (ST-4) 8,9,9
 9    ST=ST-4
      GO TO 64
 8    U1=ST
      RN=.25*ST
      RETURN
      END