diff --git a/.dev/clima_formatter_options.jl b/.dev/clima_formatter_options.jl
index 1db65e8c..68b624c9 100644
--- a/.dev/clima_formatter_options.jl
+++ b/.dev/clima_formatter_options.jl
@@ -5,4 +5,5 @@ clima_formatter_options = (
     whitespace_typedefs = true,
     whitespace_ops_in_indices = true,
     remove_extra_newlines = false,
+    ignore = ["examples"],
 )
diff --git a/.gitignore b/.gitignore
index 245dc283..05df12a5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -20,6 +20,9 @@
 !docs/src/assets/*.svg
 docs/build/
 docs/site/
+docs/src/literated/
 
 # Deps
 Manifest.toml
+
+*.jld2
diff --git a/docs/Project.toml b/docs/Project.toml
index 38eeb2a3..94e0973f 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,13 +1,17 @@
 [deps]
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 ExprTools = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
+Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
-Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 
 [compat]
+CairoMakie = "0.11"
+Documenter = "1"
 julia = "1.7"
diff --git a/docs/make.jl b/docs/make.jl
index 2e044ca2..be57774e 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,16 +1,53 @@
-using Thermodynamics, Documenter
-using DocumenterCitations
+using Thermodynamics
+using Documenter, DocumenterCitations, Literate, Printf
 
 # https://github.com/jheinen/GR.jl/issues/278#issuecomment-587090846
 ENV["GKSwstype"] = "nul"
 
 bib = CitationBibliography(joinpath(@__DIR__, "bibliography.bib"))
 
+const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
+const OUTPUT_DIR = joinpath(@__DIR__, "src/literated")
+
+example_scripts = ["density_from_temperature_pressure_humidity.jl"]
+
+if isdir(OUTPUT_DIR)
+    rm(OUTPUT_DIR)
+end
+
+mkdir(OUTPUT_DIR)
+
+cp(
+    joinpath(EXAMPLES_DIR, "JRA55_atmospheric_state_Jan_1_1991.jld2"),
+    joinpath(OUTPUT_DIR, "JRA55_atmospheric_state_Jan_1_1991.jld2"),
+)
+
+for example in example_scripts
+    example_filepath = joinpath(EXAMPLES_DIR, example)
+
+    withenv("JULIA_DEBUG" => "Literate") do
+        start_time = time_ns()
+        Literate.markdown(
+            example_filepath,
+            OUTPUT_DIR;
+            flavor = Literate.DocumenterFlavor(),
+            execute = true,
+        )
+        elapsed = 1e-9 * (time_ns() - start_time)
+        @info @sprintf("%s example took %s seconds to build.", example, elapsed)
+    end
+end
+
+example_pages = [
+    "Density from temperature, pressure, and humidity" => "literated/density_from_temperature_pressure_humidity.md",
+]
+
 pages = Any[
     "Home" => "index.md",
     "Installation" => "Installation.md",
     "API" => "API.md",
     "How-to-guide" => "HowToGuide.md",
+    "Examples" => example_pages,
     "Tested profiles" => "TestedProfiles.md",
     "Temperature profiles" => "TemperatureProfiles.md",
     "Developer docs" => "DevDocs.md",
@@ -28,10 +65,13 @@ mathengine = MathJax(
     ),
 )
 
+const MiB = 2^20
+
 format = Documenter.HTML(
     prettyurls = get(ENV, "CI", nothing) == "true",
     mathengine = mathengine,
     collapselevel = 1,
+    size_threshold = 2MiB,
 )
 
 makedocs(;
@@ -45,6 +85,27 @@ makedocs(;
     pages = pages,
 )
 
+@info "Clean up temporary .jld2 and .nc output created by doctests or literated examples..."
+
+"""
+    recursive_find(directory, pattern)
+
+Return list of filepaths within `directory` that contains the `pattern::Regex`.
+"""
+recursive_find(directory, pattern) =
+    mapreduce(vcat, walkdir(directory)) do (root, dirs, files)
+        joinpath.(root, filter(contains(pattern), files))
+    end
+
+files = []
+for pattern in [r"\.jld2", r"\.nc"]
+    global files = vcat(files, recursive_find(@__DIR__, pattern))
+end
+
+for file in files
+    rm(file)
+end
+
 deploydocs(
     repo = "github.com/CliMA/Thermodynamics.jl.git",
     target = "build",
diff --git a/examples/JRA55_atmospheric_state_Jan_1_1991.jld2 b/examples/JRA55_atmospheric_state_Jan_1_1991.jld2
new file mode 100644
index 00000000..48c22f78
Binary files /dev/null and b/examples/JRA55_atmospheric_state_Jan_1_1991.jld2 differ
diff --git a/examples/density_from_temperature_pressure_humidity.jl b/examples/density_from_temperature_pressure_humidity.jl
new file mode 100644
index 00000000..30ae4fde
--- /dev/null
+++ b/examples/density_from_temperature_pressure_humidity.jl
@@ -0,0 +1,150 @@
+# # Defining a simple parameter set and using it to compute density
+#
+# This script shows how to define a simple parameter set, and then using it to
+# compute density as a function of pressure, temperature, and humidity.
+
+# # Define parameters for computing the density of air
+#
+# First, we build a parameter set suitable for computing the density of
+# _moist_ air --- that is, a mixture of dry air, water vapor, liquid droplets,
+# and ice crystals (the latter two are called "condensates").
+
+using Thermodynamics
+using Thermodynamics.Parameters: AbstractThermodynamicsParameters
+
+struct ConstitutiveParameters{FT} <: AbstractThermodynamicsParameters{FT}
+    gas_constant::FT
+    dry_air_molar_mass::FT
+    water_molar_mass::FT
+end
+
+"""
+    ConstitutiveParameters(FT; gas_constant       = 8.3144598,
+                               dry_air_molar_mass = 0.02897,
+                               water_molar_mass   = 0.018015)
+
+Construct a set of parameters that define the density of moist air,
+```math
+ρ = p / Rᵐ(q) T,
+```
+where ``p`` is pressure, ``T`` is temperature, ``q`` defines the partition
+of total mass into vapor, liqiud, and ice mass fractions, and
+``Rᵐ`` is the effective specific gas constant for the mixture,
+```math
+Rᵐ(q) = 
+```
+where 
+For more information see [reference docs].
+"""
+function ConstitutiveParameters(
+    FT = Float64;
+    gas_constant = 8.3144598,
+    dry_air_molar_mass = 0.02897,
+    water_molar_mass = 0.018015,
+)
+
+    return ConstitutiveParameters(
+        convert(FT, gas_constant),
+        convert(FT, dry_air_molar_mass),
+        convert(FT, water_molar_mass),
+    )
+end
+
+# Next, we define functions that return:
+#   1. The specific gas constant for dry air
+#   2. The specific gas constant for water vapor
+#   3. The ratio between the dry air molar mass and water molar mass
+
+const VTP = ConstitutiveParameters
+using Thermodynamics: Parameters
+
+Parameters.R_d(p::VTP) = p.gas_constant / p.dry_air_molar_mass
+Parameters.R_v(p::VTP) = p.gas_constant / p.water_molar_mass
+Parameters.molmass_ratio(p::VTP) = p.dry_air_molar_mass / p.water_molar_mass
+
+# # The density of dry air
+
+import Thermodynamics as AtmosphericThermodynamics
+
+# To compute the density of dry air, we first build a default
+# parameter set,
+
+parameters = ConstitutiveParameters()
+
+# Next, we construct a phase partitioning representing dry air --- the trivial
+# case where the all mass ratios are 0.
+#
+# Note that the syntax for `PhasePartition` is
+#
+# ```
+# PhasePartition(q_total, q_liquid, q_ice)
+# ```
+# where the q's are mass fractions of total water components,
+# liquid droplet condensate, and ice crystal condensate.
+
+q_dry = AtmosphericThermodynamics.PhasePartition(0.0)
+
+# Finally, we define the pressure and temperature, which constitute
+# the "state" of our atmosphere,
+
+p₀ = 101325.0 # sea level pressure in Pascals (here taken to be mean sea level pressure)
+T₀ = 273.15   # temperature in Kelvin
+
+# Note that the above must be defined with the same float point precision as
+# used for the parameters and PhasePartition.
+# We're now ready to compute the density of dry air,
+
+ρ = air_density(parameters, T₀, p₀, q_dry)
+
+@show ρ
+
+# Note that the above must be defined with the same float point precision as
+# used for the parameters and PhasePartition.
+# We're now ready to compute the density of dry air,
+
+using JLD2
+
+# Next, we load an atmospheric state correpsonding to atmospheric surface
+# variables substampled from the JRA55 dataset, from the date Jan 1, 1991:
+
+@load "JRA55_atmospheric_state_Jan_1_1991.jld2" q T p
+
+# The variables `q`, `T`, and `p` correspond to the total specific humidity
+# (a mass fraction), temperature (Kelvin), and sea level pressure (Pa)
+#
+# We use `q` to build a vector of PhasePartition,
+
+qp = PhasePartition.(q)
+
+# And then compute the density using the same parameters as before:
+
+ρ = air_density.(parameters, T, p, qp)
+
+# Finally, we plot the density as a function of temperature and specific humidity,
+
+using CairoMakie
+
+## Pressure range, centered around the mean sea level pressure defined above
+pmax = maximum(abs, p)
+dp = 3 / 4 * (pmax - p₀)
+prange = (p₀ - dp, p₀ + dp)
+pmap = :balance
+
+## Compute temperature range
+Tmin = minimum(T)
+Tmax = maximum(T)
+Trange = (Tmin, Tmax)
+Tmap = :viridis
+
+fig = Figure(size = (1000, 450))
+
+axρ = Axis(fig[2, 1], xlabel = "Temperature (K) ", ylabel = "Density (kg m⁻³)")
+axq = Axis(fig[2, 2], xlabel = "Specific humidity", ylabel = "Density (kg m⁻³)")
+
+scatter!(axρ, T[:], ρ[:], color = p[:], colorrange = prange, colormap = pmap, alpha = 0.1)
+scatter!(axq, q[:], ρ[:], color = T[:], colorrange = Trange, colormap = Tmap, alpha = 0.1)
+
+Colorbar(fig[1, 1], label = "Pressure (Pa)", vertical = false, colorrange = prange, colormap = pmap)
+Colorbar(fig[1, 2], label = "Temperature (K)", vertical = false, colorrange = Trange, colormap = Tmap)
+
+current_figure()
diff --git a/src/Parameters.jl b/src/Parameters.jl
index fd51b602..f6921fbd 100644
--- a/src/Parameters.jl
+++ b/src/Parameters.jl
@@ -2,6 +2,9 @@ module Parameters
 
 export ThermodynamicsParameters
 
+abstract type AbstractThermodynamicsParameters{FT} end
+const ATP = AbstractThermodynamicsParameters
+
 """
     ThermodynamicsParameters
 
@@ -21,7 +24,8 @@ param_set = TP.ThermodynamicsParameters(toml_dict)
 
 ```
 """
-Base.@kwdef struct ThermodynamicsParameters{FT}
+Base.@kwdef struct ThermodynamicsParameters{FT} <:
+                   AbstractThermodynamicsParameters{FT}
     T_0::FT
     MSLP::FT
     p_ref_theta::FT
@@ -50,14 +54,12 @@ Base.@kwdef struct ThermodynamicsParameters{FT}
     pow_icenuc::FT
 end
 
-const ATP = ThermodynamicsParameters
-
 Base.broadcastable(ps::ATP) = tuple(ps)
 Base.eltype(::ThermodynamicsParameters{FT}) where {FT} = FT
 
 # wrappers
-for fn in fieldnames(ATP)
-    @eval @inline $(fn)(ps::ATP) = ps.$(fn)
+for fn in fieldnames(ThermodynamicsParameters)
+    @eval $(fn)(ps::ThermodynamicsParameters) = ps.$(fn)
 end
 
 # Derived parameters
diff --git a/src/TemperatureProfiles.jl b/src/TemperatureProfiles.jl
index dbe4ba8b..7dd1a889 100644
--- a/src/TemperatureProfiles.jl
+++ b/src/TemperatureProfiles.jl
@@ -8,7 +8,7 @@ export TemperatureProfile,
 
 import ..Parameters
 const TP = Parameters
-const APS = TP.ThermodynamicsParameters
+const APS = TP.AbstractThermodynamicsParameters
 
 """
     TemperatureProfile
diff --git a/src/Thermodynamics.jl b/src/Thermodynamics.jl
index 4cef0f87..fd2c6b07 100644
--- a/src/Thermodynamics.jl
+++ b/src/Thermodynamics.jl
@@ -54,7 +54,7 @@ const KA = KernelAbstractions
 include("Parameters.jl")
 import .Parameters
 const TP = Parameters
-const APS = TP.ThermodynamicsParameters
+const APS = TP.AbstractThermodynamicsParameters
 
 # Allow users to skip error on non-convergence
 # by importing:
diff --git a/src/relations.jl b/src/relations.jl
index 8d68762d..8a180895 100644
--- a/src/relations.jl
+++ b/src/relations.jl
@@ -1022,12 +1022,11 @@ end
 
 Compute the saturation specific humidity over a plane surface of condensate, given
 
- - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
- - `T` temperature
- - `ρ` (moist-)air density
-and, optionally,
- - `Liquid()` indicating condensate is liquid
- - `Ice()` indicating condensate is ice
+    - `param_set`: an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
+    - `T`: temperature
+    - `ρ`: air density
+    - (optional) `Liquid()`: indicating condensate is liquid
+    - (optional) `Ice()`: indicating condensate is ice
 """
 @inline function q_vap_saturation_generic(
     param_set::APS,
@@ -1046,12 +1045,11 @@ q_vap_saturation_generic(param_set::APS, T, ρ, phase::Phase) =
 
 Compute the saturation specific humidity, given
 
- - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
- - `T` temperature
- - `ρ` (moist-)air density
- - `phase_type` a thermodynamic state type
-and, optionally,
- - `q` [`PhasePartition`](@ref)
+- `param_set`: an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
+- `T`: temperature
+- `ρ`: air density
+- `phase_type`: a thermodynamic state type
+- (optional) `q`: [`PhasePartition`](@ref)
 
 If the `PhasePartition` `q` is given, the saturation specific humidity is that of a
 mixture of liquid and ice, computed in a thermodynamically consistent way from the