    module random_normal
    use precision
    contains
!-------------------------------
    subroutine rnor(z,n,var)
!  RNOR - generate N normally distributed random numbers of mean 0 and
!  variance VAR.  Place them into Z on return.
!  Adapted from code by
!     Alan Miller
!     CSIRO Division of Mathematical & Information Sciences
!     Private Bag 10, Clayton South MDC
!     Clayton 3169, Victoria, Australia
!     Phone: (+61) 3 9545-8016      Fax: (+61) 3 9545-8080
!     e-mail: amiller @ bigpond.net.au
!  which in turn came from Fortran 77 code:
!      ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
!      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
!      VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.
!
!  The function random_normal() returns a normally distributed pseudo-random
!  number with zero mean and unit variance.
!
!  The algorithm uses the ratio of uniforms method of A.J. Kinderman
!  and J.F. Monahan augmented with quadratic bounding curves.
!
    integer,intent(in)::n
    real(R),intent(in)::var
    real(R),intent(out)::z(n)
!
!  Local variables
!
    real(R),parameter::S=0.449871_R,T=-0.386595_R,A=0.19600_R,B=0.25472_R, &
       R1=0.27597_R,R2=0.27846_R,R3=1.7156_R
    real(R)::u(n),v(n),x,y,q
    integer::j,left,have
!
!  First generate N(0,1) deviates.
!  For each deviate, generate P = (u,v) uniform in rectangle 
!  enclosing acceptance region.  Accept P if it is inside the ellipse
!  determined by the quadratic form Q.
!
    have=0
    do while(have.lt.n)
       left=n-have
       call random_number(u(1:left))
       call random_number(v(1:left))
       do j=1,left
          x=u(j)-S
          v(j)=R3*(v(j)-0.5)
          y=abs(v(j))-T
          q=x**2+y*(A*y-B*x)
          if(q.lt.R1.or.q.le.R2.and.v(j)**2.lt.-4.0*log(u(j))*u(j)**2) then
             have=have+1
             z(have)=v(j)/u(j)
          endif
       enddo
    enddo
    z=z*sqrt(var)  !  Convert to N(0,var)
    return
    end subroutine rnor
    end module random_normal
