diff --git a/unit_testing/test_tidegauge_methods.py b/unit_testing/test_tidegauge_methods.py
index 0a5466a6..acf0ff4a 100644
--- a/unit_testing/test_tidegauge_methods.py
+++ b/unit_testing/test_tidegauge_methods.py
@@ -132,6 +132,48 @@ def test_read_gesla_formats(self):
         self.assertTrue(check1, "check1")
         self.assertTrue(check2, "check2")
 
+    def test_read_bodc_and_compare_to_model(self):
+        sci = coast.Gridded(files.fn_nemo_dat, files.fn_nemo_dom, config=files.fn_config_t_grid)
+        sci.dataset["landmask"] = sci.dataset.bottom_level == 0
+        tganalysis = coast.TidegaugeAnalysis()
+
+        with self.subTest("Read BODC file"):
+            date0 = datetime.datetime(2020, 10, 10)
+            date1 = datetime.datetime(2020, 10, 12)
+            liv = coast.Tidegauge()
+            liv.read_bodc(files.fn_tidegauge_bodc, date_start=date0, date_end=date1)
+            # TEST: Check attribute dictionary and length of sea_level.
+            check1 = len(liv.dataset.ssh.isel(id_dim=0)) == 193
+            check2 = liv.dataset.id_name.values[0] == "liverpool,_gladstone_dock"
+            self.assertTrue(check1, "check1")
+            self.assertTrue(check2, "check2")
+
+        with self.subTest("Plot Single BODC"):
+            f, a = liv.plot_on_map()
+            f.savefig(files.dn_fig + "tidegauge_bodc_map.png")
+            plt.close("all")
+
+        with self.subTest("Obs operator"):
+            liv_interp = liv.obs_operator(sci, time_interp="linear")
+
+            # TEST: Check that the resulting interp_sossheig variable is of the same
+            # length as sea_level and that it is populated.
+            interp_len = liv_interp.dataset.ssh.shape[0]
+            orig_len = liv.dataset.ssh.shape[0]
+            check1 = interp_len == orig_len
+            check2 = False in np.isnan(liv_interp.dataset.ssh)
+            self.assertTrue(check1, "check1")
+            self.assertTrue(check2, "check2")
+
+        with self.subTest("Tidegauge BODC CRPS"):
+            crps = tganalysis.crps(liv.dataset.ssh, sci.dataset.ssh)
+
+            # TEST: Check length of crps and that it contains values
+            check1 = crps.dataset.crps.shape[0] == liv.dataset.ssh.shape[0]
+            check2 = False in np.isnan(crps.dataset.crps)
+            self.assertTrue(check1, "check1")
+            self.assertTrue(check2, "check2")
+
     def test_read_gesla_and_compare_to_model(self):
         sci = coast.Gridded(files.fn_nemo_dat, files.fn_nemo_dom, config=files.fn_config_t_grid)
         sci.dataset["landmask"] = sci.dataset.bottom_level == 0
@@ -151,7 +193,7 @@ def test_read_gesla_and_compare_to_model(self):
 
         with self.subTest("Plot Single GESLA"):
             f, a = lowestoft.plot_on_map()
-            f.savefig(files.dn_fig + "tidegauge_map.png")
+            f.savefig(files.dn_fig + "tidegauge_gesla_map.png")
             plt.close("all")
 
         with self.subTest("Obs operator"):
@@ -166,7 +208,7 @@ def test_read_gesla_and_compare_to_model(self):
             self.assertTrue(check1, "check1")
             self.assertTrue(check2, "check2")
 
-        with self.subTest("Tidegauge CRPS"):
+        with self.subTest("Tidegauge GESLA CRPS"):
             crps = tganalysis.crps(lowestoft.dataset.ssh, sci.dataset.ssh)
 
             # TEST: Check length of crps and that it contains values
@@ -325,7 +367,7 @@ def test_tidegauge_cubic_spline_extrema(self):
             # it assumes it is a GESLA  type object
             tg = coast.Tidegauge()
             # specify the data read as a High Low Water dataset
-            tg.read_bodc(files.fn_tidegauge2, date_start, date_end)
+            tg.read_bodc(files.fn_tidegauge_bodc, date_start, date_end)
 
             # Use cubic spline fitting method
             extrema_cubc = tganalysis.find_high_and_low_water(tg.dataset.ssh, method="cubic")
diff --git a/unit_testing/unit_test_contents.txt b/unit_testing/unit_test_contents.txt
index 4038fa3b..a133a8dd 100755
--- a/unit_testing/unit_test_contents.txt
+++ b/unit_testing/unit_test_contents.txt
@@ -49,15 +49,16 @@
 
 9. test_tidegauge_methods
       a. load_multiple_tidegauge
-      b. read_gesla_and_compare_to_model
-      c. read_gesla_formats
-      d. tidegauge_cubic_spline_extrema
-      e. tidegauge_finding_extrema
-      f. tidegauge_for_tabulated_data
-      g. tidegauge_resample_and_apply_doodsonx0
-      h. time_mean
-      i. time_slice
-      j. time_std
+      b. read_bodc_and_compare_to_model
+      c. read_gesla_and_compare_to_model
+      d. read_gesla_formats
+      e. tidegauge_cubic_spline_extrema
+      f. tidegauge_finding_extrema
+      g. tidegauge_for_tabulated_data
+      h. tidegauge_resample_and_apply_doodsonx0
+      i. time_mean
+      j. time_slice
+      k. time_std
 
 10. test_eof_methods
       a. compute_heofs
diff --git a/unit_testing/unit_test_files.py b/unit_testing/unit_test_files.py
index afb3f9fc..c76df89f 100644
--- a/unit_testing/unit_test_files.py
+++ b/unit_testing/unit_test_files.py
@@ -35,6 +35,7 @@
 file_names_amm7 = path.join(dn_files, "nemo_data_T_grid*.nc")
 fn_altimetry = path.join(dn_files, "coast_example_altimetry_data.nc")
 fn_tidegauge = path.join(dn_files, "tide_gauges/lowestoft-p024-uk-bodc")
+fn_tidegauge_bodc = path.join(dn_files, "tide_gauges/LIV200710.txt")
 fn_tidegauge2 = path.join(dn_files, "tide_gauges/LIV2010.txt")
 fn_tidegauge_v3 = path.join(dn_files, "tide_gauges/lowestoft-p024-uk-bodc")
 fn_tidegauge_v5 = path.join(dn_files, "tide_gauges/liverpool_gladstone_dock-liv-gbr-bodc")