subroutine gauss2d(x, y)
      implicit none
     
      real x, y, rnd(2), rndm, dum, r1, r2, w
      real twopi/6.2831853/

         call rnorml(rnd, 2)

 1       r1 = 2. * rndm(dum) - 1.0
         r2 = 2. * rndm(dum) - 1.0
         w = r1 * r1 + r2 * r2
         if( w.ge.1.0 ) go to 1
         w = sqrt(( -2.0 * alog(w)) / w )
         x = r1 * w
         y = r2 * w

         return
      end