diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md
index 860aec2dd..88e692ccc 100644
--- a/doc/specs/stdlib_stats_distribution_normal.md
+++ b/doc/specs/stdlib_stats_distribution_normal.md
@@ -30,16 +30,16 @@ Note: the algorithm used for generating normal random variates is fundamentally
 
 ### Class
 
-Function
+Elemental function (passing both `loc` and `scale`).
 
 ### Arguments
 
-`array_size`: optional argument has `intent(in)` and is a scalar of type `integer`.
-
 `loc`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`.
 
 `scale`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`.
 
+`array_size`: optional argument has `intent(in)` and is a scalar of type `integer`.
+
 `loc` and `scale` arguments must be of the same type.
 
 ### Return value
diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp
index e21b7af14..de4f4b275 100644
--- a/src/stdlib_stats_distribution_normal.fypp
+++ b/src/stdlib_stats_distribution_normal.fypp
@@ -160,6 +160,7 @@ contains
 
 
     #:for k1, t1 in REAL_KINDS_TYPES
+    impure elemental &
     function rvs_norm_${t1[0]}$${k1}$(loc, scale) result(res)
     !
     ! Normal random variate (loc, scale)
@@ -178,6 +179,7 @@ contains
 
 
     #:for k1, t1 in CMPLX_KINDS_TYPES
+    impure elemental &
     function rvs_norm_${t1[0]}$${k1}$(loc, scale) result(res)
     !
     ! Normally distributed complex. The real part and imaginary part are       &