diff --git a/docs/make.jl b/docs/make.jl
index 732dd020..585b2c02 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,14 +1,18 @@
 using Documenter
 using Documenter: doctest
+using CalibrateAtmos
 using Base.CoreLogging
 using DocumenterCitations
 
 disable_logging(Base.CoreLogging.Info) # Hide doctest's `@info` printing
+bib = CitationBibliography(joinpath(@__DIR__, "bibliography.bib"))
+
+doctest(CalibrateAtmos; plugins = [bib])
 disable_logging(Base.CoreLogging.BelowMinLevel) # Re-enable all logging
 
 makedocs(
-    CitationBibliography(joinpath(@__DIR__, "bibliography.bib")),
-    modules = Vector{Module}(),
+    plugins = [bib],
+    modules = [CalibrateAtmos],
     sitename = "CalibrateAtmos.jl",
     authors = "Clima",
     strict = true,
diff --git a/docs/src/quickstart.md b/docs/src/quickstart.md
new file mode 100644
index 00000000..e69de29b
diff --git a/plot/latitude_contour_plots.jl b/plot/latitude_contour_plots.jl
index 4ccd99a4..99117287 100644
--- a/plot/latitude_contour_plots.jl
+++ b/plot/latitude_contour_plots.jl
@@ -5,14 +5,7 @@ import Statistics: mean
 
 # Adapted from contours_and_plots.jl in ClimaAtmos
 
-function time_from_filename(file)
-    arr = split(basename(file), ".")
-    day = parse(Float64, replace(arr[1], "day" => ""))
-    sec = parse(Float64, arr[2])
-    return day * (60 * 60 * 24) + sec
-end
-
-function read_hdf5_file(file_path)
+function read_diagnostics(file_path)
     reader = InputOutput.HDF5Reader(
         file_path,
         ClimaComms.SingletonCommsContext(ClimaComms.CPUSingleThreaded()),
@@ -22,6 +15,17 @@ function read_hdf5_file(file_path)
     return diagnostics
 end
 
+"""
+    time_from_filename(file_path)
+Returns the time in seconds from the name of the given file.
+"""
+function time_from_filename(file_path)
+    arr = split(basename(file_path), ".")
+    day = parse(Float64, replace(arr[1], "day" => ""))
+    sec = parse(Float64, arr[2])
+    return day * (60 * 60 * 24) + sec
+end
+
 horizontal_space(diagnostics) =
     Spaces.horizontal_space(axes(diagnostics.temperature))
 
@@ -33,6 +37,11 @@ get_column_1(diagnostics) =
     isnothing(diagnostics) ? (nothing, nothing) :
     column_view_diagnostics(diagnostics, ((1, 1), 1))
 
+
+"""
+    column_view_diagnostics(diagnostics, column)
+Retrieves the diagnostics from the given column.
+"""
 function column_view_diagnostics(diagnostics, column)
     ((i, j), h) = column
     is_extruded_field(object) =
@@ -62,6 +71,10 @@ function column_view_diagnostics(diagnostics, column)
     return column_zs_and_values, col_string
 end
 
+"""
+    column_at_coord_getter(latitude, longitude)
+Retrieves the diagnostics from the column closest to the given latitude and longitude.
+"""
 column_at_coord_getter(latitude, longitude) =
     diagnostics -> begin
         isnothing(diagnostics) && return (nothing, nothing)
@@ -96,7 +109,7 @@ function latitude_contour_plot(
     longitude = 0,
     out_path = nothing,
 )
-    diagnostics = read_hdf5_file(file_path)
+    diagnostics = read_diagnostics(file_path)
     latitudes = -90:5:90
     cols = map(latitudes) do lat
         column_at_coord_getter(lat, longitude)(diagnostics)[1]
@@ -123,3 +136,15 @@ function latitude_contour_plot(
     Makie.save(out_path, figure)
     return nothing
 end
+
+"""
+    all_variables(file_path)
+Returns the propertynames of the diagnostics of the given file path. 
+Assumes all columns have the same propertynames.
+"""
+function all_variables(file_path)
+    diagnostics = read_diagnostics(file_path)
+    latitude = -90
+    col = column_at_coord_getter(latitude, 0)(diagnostics)[1]
+    return propertynames(col[2])
+end