FUNCTION rnorm() RESULT( fn_val )

!   This subroutine was taken from: http://jblevins.org/mirror/amiller/rnorm.f90

!   Generate a random normal deviate using the polar method.
!   Reference: Marsaglia,G. & Bray,T.A. 'A convenient method for generating
!              normal variables', Siam Rev., vol.6, 260-264, 1964.

REAL(dp)  :: fn_val

! Local variables

REAL(dp)            :: u, sum
REAL(dp), SAVE      :: v, sln
LOGICAL, SAVE       :: second = .FALSE.
REAL(dp), PARAMETER :: one = 1, vsmall = TINY( one )

IF (second) THEN
! If second, use the second random number generated on last call

  second = .false.
  fn_val = v*sln

ELSE
! First call; generate a pair of random normals

  second = .true.
  DO
    CALL RANDOM_NUMBER( u )
    CALL RANDOM_NUMBER( v )
    u = 2*u - one
    v = 2*v - one
    sum = u*u + v*v + vsmall         ! vsmall added to prevent LOG(zero) / zero
    IF(sum < one) EXIT
  END DO
  sln = SQRT(- 2*LOG(sum) / sum)
  fn_val = u*sln
END IF

END FUNCTION rnorm