diff --git a/.github/workflows/unit_tests.yml b/.github/workflows/unit_tests.yml
index 8a83305a..2a41ae92 100644
--- a/.github/workflows/unit_tests.yml
+++ b/.github/workflows/unit_tests.yml
@@ -64,6 +64,7 @@ jobs:
         pytest test_utils.py
         pytest test_validate_mv_python.py
         pytest test_future_warnings.py
+        pytest test_sl1l2.py
         coverage run -m pytest test_agg_eclv.py test_agg_stats_and_boot.py test_agg_stats_with_groups.py test_calc_difficulty_index.py test_convert_lon_indices.py test_event_equalize.py test_event_equalize_against_values.py test_lon_360_to_180.py test_statistics.py test_tost_paired.py test_utils.py test_future_warnings.py
         coverage html
         
diff --git a/metcalcpy/util/sl1l2_statistics.py b/metcalcpy/util/sl1l2_statistics.py
index 9324a0dd..12508b45 100644
--- a/metcalcpy/util/sl1l2_statistics.py
+++ b/metcalcpy/util/sl1l2_statistics.py
@@ -524,10 +524,13 @@ def calculate_bcmse(input_data, columns_names, aggregation=False):
     """
     warnings.filterwarnings('error')
     try:
+
         mse = calculate_mse(input_data, columns_names, aggregation)
         me = calculate_me(input_data, columns_names, aggregation)
         result = mse - me ** 2
         result = round_half_up(result, PRECISION)
+        if result < 0:
+            return 0.
     except (TypeError, Warning):
         result = None
     warnings.filterwarnings('ignore')
diff --git a/test/test_sl1l2.py b/test/test_sl1l2.py
new file mode 100644
index 00000000..2cbce584
--- /dev/null
+++ b/test/test_sl1l2.py
@@ -0,0 +1,22 @@
+import numpy as np
+import pytest
+
+from metcalcpy.util import sl1l2_statistics as sl1l2
+
+def test_calculate_bcmse():
+    # Test that negative BCMSE values are no longer returned.
+    input_data_list = []
+
+    # These data produce negative BCMSE values. Test that the modified code no longer returns negative values
+    # for the BCMSE.
+    input_data_list.append(np.array([[4.37978400e+01, 4.70115800e+01, 1.91825108e+03, 2.21008843e+03, 2.05900571e+03, 1.00000000e+00]]))
+    input_data_list.append(np.array([[8.66233900e+01, 4.83037900e+01, 7.50361146e+03, 2.33325660e+03, 4.18423840e+03, 1.00000000e+00]]))
+    input_data_list.append(np.array([[3.68089000e+01, 1.64253370e+02, 1.35489535e+03, 2.69791703e+04, 6.04598647e+03, 1.00000000e+00]]))
+    columns_names = np.array(['fbar', 'obar', 'ffbar', 'oobar', 'fobar', 'total'], dtype='<U5')
+
+    for input in input_data_list:
+      result = sl1l2.calculate_bcmse(input, columns_names, aggregation=False)
+      assert result >= 0.
+
+
+