diff --git a/include/evalhyd/detail/probabilist/evaluator.hpp b/include/evalhyd/detail/probabilist/evaluator.hpp
index 3f94fd6d85dcdd94fa5da2224bad5dc914a8397f..619104e2a5b88dbd20f7cc6134044a5c5b92ffb5 100644
--- a/include/evalhyd/detail/probabilist/evaluator.hpp
+++ b/include/evalhyd/detail/probabilist/evaluator.hpp
@@ -36,6 +36,7 @@ namespace evalhyd
             // members for optional input data
             const XD2& _q_thr;
             const xt::xtensor<double, 1>& _c_lvl;
+            const xt::xtensor<double, 1>& _q_lvl;
             xtl::xoptional<const std::string, bool> _events;
             xt::xtensor<bool, 4> t_msk;
             const std::vector<xt::xkeep_slice<int>>& b_exp;
@@ -336,7 +337,7 @@ namespace evalhyd
                 if (!itv_bnds.has_value())
                 {
                     itv_bnds = elements::calc_itv_bnds(
-                            q_prd, get_c_lvl(),
+                            q_prd, get_c_lvl(), _q_lvl,
                             n_sit, n_ldt, n_itv, n_tim
                     );
                 }
@@ -393,7 +394,7 @@ namespace evalhyd
                 if (!qs.has_value())
                 {
                     qs = intermediate::calc_qs(
-                            q_obs, get_q_qnt(), n_mbr
+                            q_obs, get_q_qnt(), n_mbr, _q_lvl
                     );
                 }
                 return qs.value();
@@ -482,12 +483,13 @@ namespace evalhyd
                       const XD4& prd,
                       const XD2& thr,
                       const xt::xtensor<double, 1>& lvl,
+                      const xt::xtensor<double, 1>& qlvl,
                       xtl::xoptional<const std::string&, bool> events,
                       const XB4& msk,
                       const std::vector<xt::xkeep_slice<int>>& exp,
                       const long int seed) :
                     q_obs{obs}, q_prd{prd},
-                    _q_thr{thr}, _c_lvl{lvl}, _events{events},
+                    _q_thr{thr}, _c_lvl{lvl}, _q_lvl{qlvl}, _events{events},
                     t_msk(msk), b_exp(exp),
                     random_seed{seed}
             {
diff --git a/include/evalhyd/detail/probabilist/intervals.hpp b/include/evalhyd/detail/probabilist/intervals.hpp
index adbd487b592ab4903210072e9ae01836ba41e0ac..21316f2657fe5658865370bad66df97de1c00019 100644
--- a/include/evalhyd/detail/probabilist/intervals.hpp
+++ b/include/evalhyd/detail/probabilist/intervals.hpp
@@ -44,6 +44,7 @@ namespace evalhyd
             inline xt::xtensor<double, 5> calc_itv_bnds(
                     const XD4& q_prd,
                     const xt::xtensor<double, 1>& c_lvl,
+                    const xt::xtensor<double, 1>& q_lvl,
                     std::size_t n_sit,
                     std::size_t n_ldt,
                     std::size_t n_itv,
@@ -60,15 +61,38 @@ namespace evalhyd
                 xt::col(quantiles, 0) = 0.5 - c_lvl / 200.;
                 xt::col(quantiles, 1) = 0.5 + c_lvl / 200.;
 
-                // compute predictive interval bounds from quantiles
-                for (std::size_t i = 0; i < n_itv; i++)
+                // get or compute predictive interval bounds from quantiles
+                if (q_lvl.size() > 0)
+                {   
+                    for (std::size_t i = 0; i < n_itv; i++)
+                    {
+                    auto res = xt::where(xt::equal(xt::view(quantiles, i), q_lvl));
+                    if (res.size() != 2) 
+                    {
+                        throw std::runtime_error(
+                                "interval-based metric requested, "
+                                "but *c_lvl* not matching *q_lvl"
+                        );
+                    } else 
+                    {
+                        xt::view(itv_bnds, xt::all(), xt::all(), i, 0, xt::all()) =
+                                xt::view(q_prd, res[0]);
+                        xt::view(itv_bnds, xt::all(), xt::all(), i, 1, xt::all()) = 
+                                xt::view(q_prd, res[1]);
+                    }
+                    }
+                }
+                else
                 {
+                    for (std::size_t i = 0; i < n_itv; i++)
+                    {
                     auto q =  xt::quantile(q_prd, xt::view(quantiles, i), 2);
 
                     xt::view(itv_bnds, xt::all(), xt::all(), i, 0, xt::all()) =
                             xt::view(q, 0);
                     xt::view(itv_bnds, xt::all(), xt::all(), i, 1, xt::all()) =
                             xt::view(q, 1);
+                    }
                 }
 
                 return itv_bnds;
diff --git a/include/evalhyd/detail/probabilist/quantiles.hpp b/include/evalhyd/detail/probabilist/quantiles.hpp
index 88bd5284b5c3983a2b562ce8e88e37bd92d8b95a..bb13751ea119c9a4538269c29d9b9abfb1e7d62f 100644
--- a/include/evalhyd/detail/probabilist/quantiles.hpp
+++ b/include/evalhyd/detail/probabilist/quantiles.hpp
@@ -60,13 +60,22 @@ namespace evalhyd
             inline xt::xtensor<double, 4> calc_qs(
                     const XD2 &q_obs,
                     const xt::xtensor<double, 4>& q_qnt,
-                    std::size_t n_mbr
+                    std::size_t n_mbr,
+                    const xt::xtensor<double, 1>& q_lvl
             )
             {
-                // compute the quantile order $alpha$
-                xt::xtensor<double, 1> alpha =
+                xt::xtensor<double, 1> alpha;
+                // get or compute the quantile order $alpha$
+                if (q_lvl.size() > 0)
+                {
+                    alpha = xt::sort(q_lvl);
+                }
+                else
+                {
+                    alpha =
                         xt::arange<double>(1., double(n_mbr + 1))
                         / double(n_mbr + 1);
+                }
 
                 // calculate the difference
                 xt::xtensor<double, 4> diff =
diff --git a/include/evalhyd/evalp.hpp b/include/evalhyd/evalp.hpp
index f7f4340f210a2027d49d30d6fcf0a8505c0e4110..0039677616398742211b14cb5ce05e63feea3e08 100644
--- a/include/evalhyd/evalp.hpp
+++ b/include/evalhyd/evalp.hpp
@@ -181,6 +181,7 @@ namespace evalhyd
             xtl::xoptional<const std::string, bool> events =
                     xtl::missing<const std::string>(),
             const std::vector<double>& c_lvl = {},
+            const std::vector<double>& q_lvl = {},
             const xt::xexpression<XB4>& t_msk = XB4({}),
             const xt::xexpression<XS2>& m_cdt = XS2({}),
             xtl::xoptional<const std::unordered_map<std::string, int>, bool> bootstrap =
@@ -233,6 +234,7 @@ namespace evalhyd
 
         // adapt vector to tensor
         const xt::xtensor<double, 1> c_lvl_ = xt::adapt(c_lvl);
+        const xt::xtensor<double, 1> q_lvl_ = xt::adapt(q_lvl);
 
         // check that the metrics/diagnostics to be computed are valid
         utils::check_metrics(
@@ -418,7 +420,7 @@ namespace evalhyd
 
         // instantiate determinist evaluator
         probabilist::Evaluator<XD2, XD4, XB4> evaluator(
-                q_obs_, q_prd_, q_thr_, c_lvl_, events,
+                q_obs_, q_prd_, q_thr_, c_lvl_, q_lvl_, events,
                 t_msk_.size() > 0 ? t_msk_: (m_cdt_.size() > 0 ? c_msk : t_msk_),
                 b_exp,
                 random_seed