diff --git a/doc/release/upcoming_changes/21623.new_feature.rst b/doc/release/upcoming_changes/21623.new_feature.rst
new file mode 100644
index 000000000000..fafb2e8d6907
--- /dev/null
+++ b/doc/release/upcoming_changes/21623.new_feature.rst
@@ -0,0 +1,6 @@
+New parameter ``equal_nans`` added to `np.unique`
+-----------------------------------------------------------------------------------
+
+`np.unique` was changed in 1.21 to treat all ``NaN`` values as equal and return
+a single ``NaN``. Setting ``equal_nans=False`` will restore pre-1.21 behavior
+to treat ``NaNs`` as unique. Defaults to ``True``.
diff --git a/numpy/lib/arraysetops.py b/numpy/lib/arraysetops.py
index d44e1a983ebf..6d36fdcbddc9 100644
--- a/numpy/lib/arraysetops.py
+++ b/numpy/lib/arraysetops.py
@@ -131,13 +131,13 @@ def _unpack_tuple(x):
 
 
 def _unique_dispatcher(ar, return_index=None, return_inverse=None,
-                       return_counts=None, axis=None):
+                       return_counts=None, axis=None, *, equal_nans=None):
     return (ar,)
 
 
 @array_function_dispatch(_unique_dispatcher)
 def unique(ar, return_index=False, return_inverse=False,
-           return_counts=False, axis=None):
+           return_counts=False, axis=None, *, equal_nans=True):
     """
     Find the unique elements of an array.
 
@@ -162,8 +162,10 @@ def unique(ar, return_index=False, return_inverse=False,
     return_counts : bool, optional
         If True, also return the number of times each unique item appears
         in `ar`.
+    equal_nans : bool, optional
+        If True, collapses multiple NaN values in return array into 1
 
-        .. versionadded:: 1.9.0
+        .. versionchanged: 1.24
 
     axis : int or None, optional
         The axis to operate on. If None, `ar` will be flattened. If an integer,
@@ -269,7 +271,8 @@ def unique(ar, return_index=False, return_inverse=False,
     """
     ar = np.asanyarray(ar)
     if axis is None:
-        ret = _unique1d(ar, return_index, return_inverse, return_counts)
+        ret = _unique1d(ar, return_index, return_inverse, return_counts, 
+        equal_nans = equal_nans)
         return _unpack_tuple(ret)
 
     # axis was specified and not None
@@ -312,13 +315,13 @@ def reshape_uniq(uniq):
         return uniq
 
     output = _unique1d(consolidated, return_index,
-                       return_inverse, return_counts)
+                       return_inverse, return_counts, equal_nans = equal_nans)
     output = (reshape_uniq(output[0]),) + output[1:]
     return _unpack_tuple(output)
 
 
 def _unique1d(ar, return_index=False, return_inverse=False,
-              return_counts=False):
+              return_counts=False, *, equal_nans=True):
     """
     Find the unique elements of an array, ignoring shape.
     """
@@ -334,7 +337,8 @@ def _unique1d(ar, return_index=False, return_inverse=False,
         aux = ar
     mask = np.empty(aux.shape, dtype=np.bool_)
     mask[:1] = True
-    if aux.shape[0] > 0 and aux.dtype.kind in "cfmM" and np.isnan(aux[-1]):
+    if (equal_nans and aux.shape[0] > 0 and aux.dtype.kind in "cfmM" and 
+            np.isnan(aux[-1])):
         if aux.dtype.kind == "c":  # for complex all NaNs are considered equivalent
             aux_firstnan = np.searchsorted(np.isnan(aux), True, side='left')
         else:
diff --git a/numpy/lib/tests/test_arraysetops.py b/numpy/lib/tests/test_arraysetops.py
index 13385cd2409d..f97fea310dbb 100644
--- a/numpy/lib/tests/test_arraysetops.py
+++ b/numpy/lib/tests/test_arraysetops.py
@@ -765,3 +765,11 @@ def _run_axis_tests(self, dtype):
         assert_array_equal(uniq[:, inv], data)
         msg = "Unique's return_counts=True failed with axis=1"
         assert_array_equal(cnt, np.array([2, 1, 1]), msg)
+
+    def test_unique_nanequals(self):
+        # issue 20326
+        a = np.array([1, 1, np.nan, np.nan, np.nan])
+        unq = np.unique(a)
+        not_unq = np.unique(a, equal_nans = False)
+        assert_array_equal(unq, np.array([1, np.nan]))
+        assert_array_equal(not_unq, np.array([1, np.nan, np.nan, np.nan]))