gaussrnd Function

function gaussrnd()
  real(8)       :: gaussrnd
  real(8)       :: fac,v1,v2,r
  real(8),save :: gset
  integer,save :: iset=0
  if (iset==0) then ! Create a new RN
     r=100.0
     do while (r>1.0)
        v1 = 2.0*grnd()-1.0
        v2 = 2.0*grnd()-1.0
        r = v1*v1+v2*v2
     end do
     fac = sqrt(-2.0*log(r)/r)
     gset = v1*fac
     gaussrnd = v2*fac
     iset = 1
  else ! Use the 2nd NR from the previous call
     gaussrnd = gset
     iset = 0
  endif
  return
end function gaussrnd