diff --git a/MANIFEST.in b/MANIFEST.in
index b4dcb420..77b292de 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -7,6 +7,7 @@ exclude camb/camblib.so
include camb/PArthENoPE_880.2_marcucci.dat
include camb/PArthENoPE_880.2_standard.dat
include camb/PRIMAT_Yp_DH_Error.dat
+include camb/PRIMAT_Yp_DH_ErrorMC_2021.dat
exclude camb/HighLExtrapTemplate_lenspotentialCls.dat
exclude fortran/sigma8.f90
exclude fortran/writefits.f90
diff --git a/camb/bbn.py b/camb/bbn.py
index 84d03803..4dfac0b8 100644
--- a/camb/bbn.py
+++ b/camb/bbn.py
@@ -74,8 +74,8 @@ class BBN_table_interpolator(BBNPredictor):
BBN predictor based on interpolation from a numerical table calculated by a BBN code.
Tables are supplied for `Parthenope `_ 2017 (PArthENoPE_880.2_standard.dat, default),
- similar but with Marucci rates (PArthENoPE_880.2_marcucci.dat), and
- `PRIMAT `_ (PRIMAT_Yp_DH_Error.dat).
+ similar but with Marucci rates (PArthENoPE_880.2_marcucci.dat),
+ `PRIMAT `_ (PRIMAT_Yp_DH_Error.dat, PRIMAT_Yp_DH_ErrorMC_2021.dat).
:param interpolation_table: filename of interpolation table to use.
:param function_of: two variables that determine the interpolation grid (x,y) in the table,
@@ -230,3 +230,5 @@ def get_predictor(predictor_name=None):
print(BBN_fitting_parthenope().DH(0.02463, -0.6))
print(BBN_table_interpolator('PArthENoPE_880.2_marcucci.dat').DH(0.02463, -0.6))
print(BBN_table_interpolator('PRIMAT_Yp_DH_Error.dat').DH(0.02463, -0.6))
+ print(BBN_table_interpolator('PRIMAT_Yp_DH_ErrorMC_2021.dat').DH(0.02463, -0.6))
+
diff --git a/setup.py b/setup.py
index 5d4b8895..0bef88e3 100644
--- a/setup.py
+++ b/setup.py
@@ -287,7 +287,7 @@ def run(self):
platforms="any",
package_data={'camb': [DLLNAME, 'HighLExtrapTemplate_lenspotentialCls.dat',
'PArthENoPE_880.2_marcucci.dat', 'PArthENoPE_880.2_standard.dat',
- 'PRIMAT_Yp_DH_Error.dat']},
+ 'PRIMAT_Yp_DH_Error.dat', 'PRIMAT_Yp_DH_ErrorMC_2021.dat']},
test_suite='camb.tests',
entry_points={
'console_scripts': [