diff --git a/tests/test_probabilist.cpp b/tests/test_probabilist.cpp
index e25233e2af6069df1d6cb61ce672b1c2b4130c83..1e15bf9c372c1320004af1bf60e3aee1086dde04 100644
--- a/tests/test_probabilist.cpp
+++ b/tests/test_probabilist.cpp
@@ -29,7 +29,8 @@ std::vector<std::string> all_metrics_p = {
         "BS", "BSS", "BS_CRD", "BS_LBD",
         "QS", "CRPS",
         "POD", "POFD", "FAR", "CSI", "ROCSS",
-        "RANK_DIAG", "DS", "AS"
+        "RANK_DIAG", "DS", "AS",
+        "CR", "AW", "AWN", "AWI", "WS", "WSS"
 };
 
 std::tuple<xt::xtensor<double, 1>, xt::xtensor<double, 2>> load_data_p()
@@ -406,7 +407,6 @@ TEST(ProbabilistTests, TestContingency)
     );
 }
 
-
 TEST(ProbabilistTests, TestRanks)
 {
     // read in data
@@ -424,6 +424,7 @@ TEST(ProbabilistTests, TestRanks)
                     {"RANK_DIAG", "DS", "AS"},
                     xt::xtensor<double, 2>({}),
                     "high",  // events
+                    {},  // c_lvl
                     xt::xtensor<bool, 4>({}),  // t_msk
                     {},  // m_cdt
                     xtl::missing<const std::unordered_map<std::string, int>>(),  // bootstrap
@@ -497,6 +498,52 @@ TEST(ProbabilistTests, TestRanks)
     );
 }
 
+TEST(ProbabilistTests, TestIntervals)
+{
+    // read in data
+    xt::xtensor<double, 1> observed;
+    xt::xtensor<double, 2> predicted;
+    std::tie(observed, predicted) = load_data_p();
+
+    // compute scores
+    std::vector<xt::xarray<double>> metrics =
+            evalhyd::evalp(
+                    // shape: (sites [1], time [t])
+                    xt::eval(xt::view(observed, xt::newaxis(), xt::all())),
+                    // shape: (sites [1], lead times [1], members [m], time [t])
+                    xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
+                    {"CR", "AW", "AWN", "AWI", "WS", "WSS"},
+                    xt::xtensor<double, 2>({}),
+                    "",  // events
+                    {30., 80.}  // c_lvl
+            );
+
+    // check results
+    // coverage ratios
+    xt::xtensor<double, 5> cr = {{{{{ 0.006431, 0.03537 }}}}};
+    EXPECT_TRUE(xt::all(xt::isclose(metrics[0], cr, 1e-05, 1e-06, true)));
+
+    // average widths
+    xt::xtensor<double, 5> aw = {{{{{ 9.27492, 31.321543}}}}};
+    EXPECT_TRUE(xt::all(xt::isclose(metrics[1], aw, 1e-05, 1e-06, true)));
+
+    // average widths normalised
+    xt::xtensor<double, 5> awn = {{{{{ 0.007383, 0.024931}}}}};
+    EXPECT_TRUE(xt::all(xt::isclose(metrics[2], awn, 1e-05, 1e-06, true)));
+
+    // average widths indices
+    xt::xtensor<double, 5> awi = {{{{{ 0.982112, 0.988095}}}}};
+    EXPECT_TRUE(xt::all(xt::isclose(metrics[3], awi, 1e-05, 1e-06, true)));
+
+    // Winkler scores
+    xt::xtensor<double, 5> ws = {{{{{ 764.447175, 2578.138264}}}}};
+    EXPECT_TRUE(xt::all(xt::isclose(metrics[4], ws, 1e-05, 1e-06, true)));
+
+    // Winkler skill scores
+    xt::xtensor<double, 5> wss = {{{{{ 0.662189, 0.436039}}}}};
+    EXPECT_TRUE(xt::all(xt::isclose(metrics[5], wss, 1e-05, 1e-06, true)));
+}
+
 TEST(ProbabilistTests, TestMasks)
 {
     // read in data
@@ -512,6 +559,7 @@ TEST(ProbabilistTests, TestMasks)
 
     // compute scores using masks to subset whole record
     xt::xtensor<double, 2> thresholds = {{690, 534, 445}};
+    std::vector<double> confidence_levels = {30., 80.};
 
     std::vector<xt::xarray<double>> metrics_masked =
             evalhyd::evalp(
@@ -522,6 +570,7 @@ TEST(ProbabilistTests, TestMasks)
                     all_metrics_p,
                     thresholds,
                     "high",
+                    confidence_levels,
                     // shape: (sites [1], lead times [1], subsets [1], time [t])
                     masks
             );
@@ -535,7 +584,8 @@ TEST(ProbabilistTests, TestMasks)
                     xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::range(20, _))),
                     all_metrics_p,
                     thresholds,
-                    "high"
+                    "high",
+                    confidence_levels
             );
 
     // check results are identical
@@ -563,6 +613,7 @@ TEST(ProbabilistTests, TestMasks)
 TEST(ProbabilistTests, TestMaskingConditions)
 {
     xt::xtensor<double, 2> thresholds = {{690, 534, 445}};
+    std::vector<double> confidence_levels = {30., 80.};
 
     // read in data
     xt::xtensor<double, 1> observed_;
@@ -589,6 +640,7 @@ TEST(ProbabilistTests, TestMaskingConditions)
                     all_metrics_p,
                     thresholds,
                     "high",
+                    confidence_levels,
                     masks,
                     q_conditions
             );
@@ -600,7 +652,8 @@ TEST(ProbabilistTests, TestMaskingConditions)
                     xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
                     all_metrics_p,
                     thresholds,
-                    "high"
+                    "high",
+                    confidence_levels
             );
 
     // check results are identical
@@ -644,6 +697,7 @@ TEST(ProbabilistTests, TestMaskingConditions)
                     all_metrics_p,
                     thresholds,
                     "high",
+                    confidence_levels,
                     masks,
                     q_conditions_
             );
@@ -655,7 +709,8 @@ TEST(ProbabilistTests, TestMaskingConditions)
                     xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
                     all_metrics_p,
                     thresholds,
-                    "high"
+                    "high",
+                    confidence_levels
             );
 
     // check results are identical
@@ -696,6 +751,7 @@ TEST(ProbabilistTests, TestMaskingConditions)
                     all_metrics_p,
                     thresholds,
                     "high",
+                    confidence_levels,
                     masks,
                     t_conditions
             );
@@ -707,7 +763,8 @@ TEST(ProbabilistTests, TestMaskingConditions)
                     xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::range(0, 100))),
                     all_metrics_p,
                     thresholds,
-                    "high"
+                    "high",
+                    confidence_levels
             );
 
     // check results are identical
@@ -737,8 +794,8 @@ TEST(ProbabilistTests, TestMaskingConditions)
 
 TEST(ProbabilistTests, TestMissingData)
 {
-    xt::xtensor<double, 2> thresholds
-        {{ 4., 5. }};
+    xt::xtensor<double, 2> thresholds = {{ 4., 5. }};
+    std::vector<double> confidence_levels = {30., 80.};
 
     // compute metrics on series with NaN
     xt::xtensor<double, 4> forecast_nan {{
@@ -761,7 +818,8 @@ TEST(ProbabilistTests, TestMissingData)
                 forecast_nan,
                 all_metrics_p,
                 thresholds,
-                "high"
+                "high",
+                confidence_levels
         );
 
     // compute metrics on manually subset series (one leadtime at a time)
@@ -781,7 +839,8 @@ TEST(ProbabilistTests, TestMissingData)
                 forecast_pp1,
                 all_metrics_p,
                 thresholds,
-                "high"
+                "high",
+                confidence_levels
         );
 
     xt::xtensor<double, 4> forecast_pp2 {{
@@ -800,7 +859,8 @@ TEST(ProbabilistTests, TestMissingData)
                 forecast_pp2,
                 all_metrics_p,
                 thresholds,
-                "high"
+                "high",
+                confidence_levels
         );
 
     // check that numerical results are identical
@@ -824,8 +884,8 @@ TEST(ProbabilistTests, TestMissingData)
 
 TEST(ProbabilistTests, TestBootstrap)
 {
-    xt::xtensor<double, 2> thresholds
-            {{ 33.87, 55.67 }};
+    xt::xtensor<double, 2> thresholds = {{ 33.87, 55.67 }};
+    std::vector<double> confidence_levels = {30., 80.};
 
     // read in data
     std::ifstream ifs;
@@ -854,6 +914,7 @@ TEST(ProbabilistTests, TestBootstrap)
                     all_metrics_p,
                     thresholds,
                     "high",  // events
+                    confidence_levels,
                     xt::xtensor<bool, 4>({}),  // t_msk
                     {},  // m_cdt
                     bootstrap,
@@ -877,7 +938,8 @@ TEST(ProbabilistTests, TestBootstrap)
                     xt::eval(xt::view(predicted_x3, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
                     all_metrics_p,
                     thresholds,
-                    "high"
+                    "high",
+                    confidence_levels
             );
 
     // check results are identical