diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl
index b6a662d4a0..e50f721836 100644
--- a/src/Oceananigans.jl
+++ b/src/Oceananigans.jl
@@ -114,8 +114,10 @@ export
     BinaryOutputWriter,
     NetCDFOutputWriter,
     JLD2OutputWriter,
+    Checkpointer,
     write_output,
     read_output,
+    restore_from_checkpoint,
 
     # Model diagnostics
     Diagnostic,
@@ -147,6 +149,7 @@ using
 using
     Adapt,
     FFTW,
+    Distributed,
     StaticArrays,
     OffsetArrays,
     JLD2,
diff --git a/src/output_writers.jl b/src/output_writers.jl
index 7a2fd969bb..0b50aad430 100644
--- a/src/output_writers.jl
+++ b/src/output_writers.jl
@@ -1,8 +1,218 @@
-using Distributed
-using NetCDF, JLD2
-
 ext(fw::OutputWriter) = throw("Extension for $(typeof(fw)) is not implemented.")
-time_to_write(clock, out::OutputWriter) = (clock.iteration % out.frequency) == 0
+
+function time_to_write(clock::Clock, diag::OutputWriter)
+    if :interval in propertynames(diag) && diag.interval != nothing
+        if clock.time >= diag.previous + diag.interval
+            diag.previous = clock.time - rem(clock.time, diag.interval)
+            return true
+        else
+            return false
+        end
+    elseif :frequency in propertynames(diag) && diag.frequency != nothing
+        return clock.iteration % diag.frequency == 0
+    else
+        error("$(typeof(diag)) must have a frequency or interval specified!")
+    end
+end
+
+function validate_interval(frequency, interval)
+    isnothing(frequency) && isnothing(interval) && @error "Must specify a frequency or interval!"
+    return
+end
+
+# When saving stuff to disk like a JLD2 file, `saveproperty!` is used, which
+# converts Julia objects to language-agnostic objects.
+saveproperty!(file, location, p::Number)        = file[location] = p
+saveproperty!(file, location, p::AbstractRange) = file[location] = collect(p)
+saveproperty!(file, location, p::AbstractArray) = file[location] = Array(p)
+saveproperty!(file, location, p::Field)         = file[location] = Array(p.data.parent)
+saveproperty!(file, location, p::Function) = @warn "Cannot save Function property into $location"
+saveproperty!(file, location, p) = [saveproperty!(file, location * "/$subp", getproperty(p, subp)) for subp in propertynames(p)]
+
+# Special saveproperty! so boundary conditions are easily readable outside julia.    
+function saveproperty!(file, location, cbcs::CoordinateBoundaryConditions)
+    for endpoint in propertynames(cbcs)
+        endpoint_bc = getproperty(cbcs, endpoint)
+        if isa(endpoint_bc.condition, Function)
+            @warn "$field.$coord.$endpoint boundary is of type Function and cannot be saved to disk!"
+            file["boundary_conditions/$field/$coord/$endpoint/type"] = string(bctype(endpoint_bc))
+            file["boundary_conditions/$field/$coord/$endpoint/condition"] = missing
+        else
+            file["boundary_conditions/$field/$coord/$endpoint/type"] = string(bctype(endpoint_bc))
+            file["boundary_conditions/$field/$coord/$endpoint/condition"] = endpoint_bc.condition
+        end
+    end
+end
+
+saveproperties!(file, structure, ps) = [saveproperty!(file, "$p", getproperty(structure, p)) for p in ps]
+
+# When checkpointing, `serializeproperty!` is used, which serializes objects
+# unless they need to be converted (basically CuArrays only).
+serializeproperty!(file, location, p) = file[location] = p
+serializeproperty!(file, location, p::Union{AbstractArray, Field}) = saveproperty!(file, location, p)
+serializeproperty!(file, location, p::Function) = @warn "Cannot serialize Function property into $location"
+
+serializeproperties!(file, structure, ps) = [serializeproperty!(file, "$p", getproperty(structure, p)) for p in ps]
+
+# Don't check arrays because we don't need that noise.
+has_reference(T, ::AbstractArray{<:Number}) = false
+
+# These two conditions are true, but should not necessary.
+has_reference(::Type{Function}, ::Field) = false
+has_reference(::Type{T}, ::NTuple{N, <:T}) where {N, T} = true
+
+function has_reference(has_type, obj)
+    if typeof(obj) <: has_type
+        return true
+    elseif applicable(iterate, obj) && length(obj) > 1
+        return any([has_reference(has_type, elem) for elem in obj])
+    elseif applicable(propertynames, obj) && length(propertynames(obj)) > 0
+        return any([has_reference(has_type, getproperty(obj, p)) for p in propertynames(obj)])
+    else
+        return typeof(obj) <: has_type
+    end
+end
+
+
+####
+####  JLD2 output writer
+####
+
+mutable struct JLD2OutputWriter{F, I, O, IF, IN} <: OutputWriter
+        filepath :: String
+         outputs :: O
+        interval :: I
+       frequency :: F
+            init :: IF
+       including :: IN
+        previous :: Float64
+            part :: Int
+    max_filesize :: Float64
+           async :: Bool
+           force :: Bool
+         verbose :: Bool
+end
+
+noinit(args...) = nothing
+
+"""
+    JLD2OutputWriter(model, outputs; interval=nothing, frequency=nothing, dir=".",
+                     prefix="", init=noinit, including=[:grid, :eos, :constants, :closure],
+                     part=1, max_filesize=Inf, force=false, async=false, verbose=false)
+
+Construct an `OutputWriter` that writes `label, func` pairs in the dictionary `outputs` to
+a single JLD2 file, where `label` is a symbol that labels the output and `func` is a
+function of the form `func(model)` that returns the data to be saved.
+
+# Keyword arguments
+- `frequency::Int`:    Save output every `n` model iterations.
+- `interval::Int`:     Save output every `t` units of model clock time.
+- `dir::String`:       Directory to save output to. Default: "." (current working
+                       directory).
+- `prefix::String`:    Descriptive filename prefixed to all output files. Default: "".
+- `init::Function`:    A function of the form `init(file, model)` that runs when a JLD2
+                       output file is initialized. Default: `noinit(args...) = nothing`.
+- `including::Array`:  List of model properties to save with every file. By default, the
+                       grid, equation of state, planetary constants, and the turbulence
+                       closure parameters are saved.
+- `part::Int`:         The starting part number used if `max_filesize` is finite.
+                       Default: 1.
+- `max_filesize::Int`: The writer will stop writing to the output file once the file size
+                       exceeds `max_filesize`, and write to a new one with a consistent
+                       naming scheme ending in `part1`, `part2`, etc. Defaults to `Inf`.
+- `force::Bool`:       Remove existing files if their filenames conflict. Default: `false`.
+- `async::Bool`:       Write output asynchronously. Default: `false`.
+- `verbose::Bool`:     Log what the output writer is doing with statistics on compute/write
+                       times and file sizes. Default: `false`.
+"""
+function JLD2OutputWriter(model, outputs; interval=nothing, frequency=nothing, dir=".", prefix="",
+                          init=noinit, including=[:grid, :eos, :constants, :closure],
+                          part=1, max_filesize=Inf, force=false, async=false, verbose=false)
+
+    validate_interval(interval, frequency)
+
+    mkpath(dir)
+    filepath = joinpath(dir, prefix * ".jld2")
+    force && isfile(filepath) && rm(filepath, force=true)
+
+    jldopen(filepath, "a+") do file
+        init(file, model)
+        saveproperties!(file, model, including)
+    end
+
+    return JLD2OutputWriter(filepath, outputs, interval, frequency, init, including,
+                            0.0, part, max_filesize, async, force, verbose)
+end
+
+function write_output(model, fw::JLD2OutputWriter)
+    verbose = fw.verbose
+    verbose && @info @sprintf("Calculating JLD2 output %s...", keys(fw.outputs))
+    tc = Base.@elapsed data = Dict((name, f(model)) for (name, f) in fw.outputs)
+    verbose && @info "Calculation time: $(prettytime(tc))"
+
+    iter = model.clock.iteration
+    time = model.clock.time
+
+    filesize(fw.filepath) >= fw.max_filesize && start_next_file(model, fw)
+
+    path = fw.filepath
+    verbose && @info "Writing JLD2 output $(keys(fw.outputs))..."
+    t0, sz = time_ns(), filesize(path)
+
+    if fw.async
+        @async remotecall(jld2output!, 2, path, iter, time, data)
+    else
+        jld2output!(path, iter, time, data)
+    end
+
+    t1, newsz = time_ns(), filesize(path)
+    verbose && @info "Writing done: time=$(prettytime((t1-t0)/1e9)), size=$(pretty_filesize(newsz)), Δsize=$(pretty_filesize(newsz-sz))"
+
+    return
+end
+
+function start_next_file(model::Model, fw::JLD2OutputWriter)
+    verbose = fw.verbose
+    sz = filesize(fw.filepath)
+    verbose && @info "Filesize $(pretty_filesize(sz)) has exceeded maximum file size $(pretty_filesize(fw.max_filesize))."
+
+    if fw.part == 1
+        part1_path = replace(fw.filepath, r".jld2$" => "_part1.jld2")
+        verbose && @info "Renaming first part: $(fw.filepath) -> $part1_path"
+        mv(fw.filepath, part1_path, force=fw.force)
+        fw.filepath = part1_path
+    end
+
+    fw.part += 1
+    fw.filepath = replace(fw.filepath, r"part\d+.jld2$" => "part" * string(fw.part) * ".jld2")
+    fw.force && isfile(fw.filepath) && rm(fw.filepath, force=true)
+    verbose && @info "Now writing to: $(fw.filepath)"
+
+    jldopen(fw.filepath, "a+") do file
+        fw.init(file, model)
+        saveproperties!(file, model, fw.including)
+    end
+end
+
+function jld2output!(path, iter, time, data)
+    jldopen(path, "r+") do file
+        file["timeseries/t/$iter"] = time
+        for (name, datum) in data
+            file["timeseries/$name/$iter"] = datum
+        end
+    end
+    return
+end
+
+function time_to_write(clock, out::JLD2OutputWriter{<:Nothing})
+    if clock.time > out.previous + out.interval
+        out.previous = clock.time - rem(clock.time, out.interval)
+        return true
+    else
+        return false
+    end
+end
+
 
 ####
 #### Binary output writer
@@ -91,7 +301,7 @@ function write_output(model::Model, fw::NetCDFOutputWriter)
         write_output_netcdf(fw, fields, model.clock.iteration)
     end
 
-    return nothing
+    return
 end
 
 function write_output_netcdf(fw::NetCDFOutputWriter, fields, iteration)
@@ -156,7 +366,7 @@ function write_output_netcdf(fw::NetCDFOutputWriter, fields, iteration)
 
     ncclose(filepath)
 
-    return nothing
+    return
 end
 
 function read_output(fw::NetCDFOutputWriter, field_name, iter)
@@ -167,220 +377,101 @@ function read_output(fw::NetCDFOutputWriter, field_name, iter)
 end
 
 ####
-####  JLD2 output writer
+#### Checkpointer
 ####
 
-mutable struct JLD2OutputWriter{F, I, O, IF, IN} <: OutputWriter
-        filepath :: String
-         outputs :: O
-        interval :: I
-       frequency :: F
-            init :: IF
-       including :: IN
-        previous :: Float64
-            part :: Int
-    max_filesize :: Float64
-           async :: Bool
-           force :: Bool
-         verbose :: Bool
+mutable struct Checkpointer{I, T, P, A} <: OutputWriter
+         frequency :: I
+          interval :: T
+               dir :: String
+            prefix :: String
+        properties :: P
+             force :: Bool
+    has_array_refs :: A
 end
 
-noinit(args...) = nothing
-
-"""
-    JLD2OutputWriter(model, outputs; interval=nothing, frequency=nothing, dir=".",
-                     prefix="", init=noinit, including=[:grid, :eos, :constants, :closure],
-                     part=1, max_filesize=Inf, force=false, async=false, verbose=false)
+function Checkpointer(model; frequency=nothing, interval=nothing, dir=".", prefix="checkpoint", force=false,
+                      properties = [:arch, :boundary_conditions, :grid, :clock, :eos, :constants, :closure,
+                                    :velocities, :tracers, :timestepper])
 
-Construct an `OutputWriter` that writes `label, func` pairs in the dictionary `outputs` to
-a single JLD2 file, where `label` is a symbol that labels the output and `func` is a
-function of the form `func(model)` that returns the data to be saved.
+    validate_interval(frequency, interval)
 
-# Keyword arguments
-- `frequency::Int`:    Save output every `n` model iterations.
-- `interval::Int`:     Save output every `t` units of model clock time.
-- `dir::String`:       Directory to save output to. Default: "." (current working
-                       directory).
-- `prefix::String`:    Descriptive filename prefixed to all output files. Default: "".
-- `init::Function`:    A function of the form `init(file, model)` that runs when a JLD2
-                       output file is initialized. Default: `noinit(args...) = nothing`.
-- `including::Array`:  List of model properties to save with every file. By default, the
-                       grid, equation of state, planetary constants, and the turbulence
-                       closure parameters are saved.
-- `part::Int`:         The starting part number used if `max_filesize` is finite.
-                       Default: 1.
-- `max_filesize::Int`: The writer will stop writing to the output file once the file size
-                       exceeds `max_filesize`, and write to a new one with a consistent
-                       naming scheme ending in `part1`, `part2`, etc. Defaults to `Inf`.
-- `force::Bool`:       Remove existing files if their filenames conflict. Default: `false`.
-- `async::Bool`:       Write output asynchronously. Default: `false`.
-- `verbose::Bool`:     Log what the output writer is doing with statistics on compute/write
-                       times and file sizes. Default: `false`.
-"""
-function JLD2OutputWriter(model, outputs; interval=nothing, frequency=nothing, dir=".", prefix="",
-                          init=noinit, including=[:grid, :eos, :constants, :closure],
-                          part=1, max_filesize=Inf, force=false, async=false, verbose=false)
+    # Grid needs to be checkpointed for restoring to work.
+    :grid ∉ properties && push!(properties, :grid)
 
-    interval === nothing && frequency === nothing &&
-        error("Either interval or frequency must be passed to the JLD2OutputWriter!")
+    has_array_refs = Symbol[]
 
-    mkpath(dir)
-    filepath = joinpath(dir, prefix * ".jld2")
-    force && isfile(filepath) && rm(filepath, force=true)
+    for p in properties
+        isa(p, Symbol) || @error "Property $p to be checkpointed must be a Symbol."
+        p ∉ propertynames(model) && @error "Cannot checkpoint $p, it is not a model property!"
+        
+        if has_reference(Function, getproperty(model, p))
+            @warn "model.$p contains a function somewhere in its hierarchy and will not be checkpointed."
+            filter!(e -> e != p, properties)
+        end
 
-    jldopen(filepath, "a+") do file
-        init(file, model)
-        savesubstructs!(file, model, including)
+        has_reference(Field, getproperty(model, p)) && push!(has_array_refs, p)
     end
 
-    return JLD2OutputWriter(filepath, outputs, interval, frequency, init, including,
-                            0.0, part, max_filesize, async, force, verbose)
+    mkpath(dir)
+    return Checkpointer(frequency, interval, dir, prefix, properties, force, has_array_refs)
 end
 
-function savesubstruct!(file, model, name, flds=propertynames(getproperty(model, name)))
-    for fld in flds
-        file["$name/$fld"] = getproperty(getproperty(model, name), fld)
+function write_output(model, c::Checkpointer)
+    filepath = joinpath(c.dir, c.prefix * string(model.clock.iteration) * ".jld2")
+    jldopen(filepath, "w") do file
+        file["checkpointed_properties"] = c.properties
+        file["has_array_refs"] = c.has_array_refs
+
+        serializeproperties!(file, model, filter(e -> !(e ∈ c.has_array_refs), c.properties))
+        saveproperties!(file, model, filter(e -> e ∈ c.has_array_refs, c.properties))
     end
-    return nothing
 end
 
-savesubstructs!(file, model, names) = [savesubstruct!(file, model, name) for name in names]
+_arr(::CPU, a) = a
+_arr(::GPU, a) = CuArray(a)
 
-function write_output(model, fw::JLD2OutputWriter)
-    verbose = fw.verbose
-    verbose && @info @sprintf("Calculating JLD2 output %s...", keys(fw.outputs))
-    tc = Base.@elapsed data = Dict((name, f(model)) for (name, f) in fw.outputs)
-    verbose && @info "Calculation time: $(prettytime(tc))"
-
-    iter = model.clock.iteration
-    time = model.clock.time
-
-    filesize(fw.filepath) >= fw.max_filesize && start_next_file(model, fw)
-
-    path = fw.filepath
-    verbose && @info "Writing JLD2 output $(keys(fw.outputs))..."
-    t0, sz = time_ns(), filesize(path)
-
-    if fw.async
-        @async remotecall(jld2output!, 2, path, iter, time, data)
+function restore_fields!(model, file, arch, fieldset; location="$fieldset")
+    if fieldset == :timestepper
+        restore_fields!(model.timestepper, file, arch, :Gⁿ; location="timestepper/Gⁿ")
+        restore_fields!(model.timestepper, file, arch, :G⁻; location="timestepper/G⁻")
     else
-        jld2output!(path, iter, time, data)
+        for p in propertynames(getproperty(model, fieldset))
+            getproperty(getproperty(model, fieldset), p).data.parent .= _arr(arch, file[location * "/$p"])
+        end
     end
-
-    t1, newsz = time_ns(), filesize(path)
-    verbose && @info "Writing done: time=$(prettytime((t1-t0)/1e9)), size=$(pretty_filesize(newsz)), Δsize=$(pretty_filesize(newsz-sz))"
-
-    return nothing
 end
 
-function start_next_file(model::Model, fw::JLD2OutputWriter)
-    verbose = fw.verbose
-    sz = filesize(fw.filepath)
-    verbose && @info "Filesize $(pretty_filesize(sz)) has exceeded maximum file size $(pretty_filesize(fw.max_filesize))."
+#const array_containing_structs = (:velocities, :tracers, :timestepper)
 
-    if fw.part == 1
-        part1_path = replace(fw.filepath, r".jld2$" => "_part1.jld2")
-        verbose && @info "Renaming first part: $(fw.filepath) -> $part1_path"
-        mv(fw.filepath, part1_path, force=fw.force)
-        fw.filepath = part1_path
-    end
-
-    fw.part += 1
-    fw.filepath = replace(fw.filepath, r"part\d+.jld2$" => "part" * string(fw.part) * ".jld2")
-    fw.force && isfile(fw.filepath) && rm(fw.filepath, force=true)
-    verbose && @info "Now writing to: $(fw.filepath)"
+function restore_from_checkpoint(filepath; kwargs=Dict())
+    file = jldopen(filepath, "r")
 
-    jldopen(fw.filepath, "a+") do file
-        fw.init(file, model)
-        savesubstructs!(file, model, fw.including)
-    end
-end
+    cps = file["checkpointed_properties"]
+    has_array_refs = file["has_array_refs"]
 
-function jld2output!(path, iter, time, data)
-    jldopen(path, "r+") do file
-        file["timeseries/t/$iter"] = time
-        for (name, datum) in data
-            file["timeseries/$name/$iter"] = datum
+    # Restore model properties that were just serialized.
+    # We skip fields that contain structs and restore them after model creation.
+    for p in cps
+        if p ∉ has_array_refs
+            kwargs[p] = file["$p"]
         end
     end
-    return nothing
-end
-
-function time_to_write(clock, out::JLD2OutputWriter{<:Nothing})
-    if clock.time > out.previous + out.interval
-        out.previous = clock.time - rem(clock.time, out.interval)
-        return true
-    else
-        return false
-    end
-end
-
-####
-#### Output utils
-####
-
-"""
-    HorizontalAverages(arch, grid)
 
-Instantiate a struct of 1D arrays on `arch` and `grid`
-for storing the results of horizontal averages of 3D fields.
-"""
-struct HorizontalAverages{A}
-    U :: A
-    V :: A
-    T :: A
-    S :: A
-end
-
-function HorizontalAverages(arch::CPU, grid::Grid{FT}) where FT
-    U = zeros(FT, 1, 1, grid.Tz)
-    V = zeros(FT, 1, 1, grid.Tz)
-    T = zeros(FT, 1, 1, grid.Tz)
-    S = zeros(FT, 1, 1, grid.Tz)
-
-    HorizontalAverages(U, V, T, S)
-end
-
-function HorizontalAverages(arch::GPU, grid::Grid{FT}) where FT
-    U = CuArray{FT}(undef, 1, 1, grid.Tz)
-    V = CuArray{FT}(undef, 1, 1, grid.Tz)
-    T = CuArray{FT}(undef, 1, 1, grid.Tz)
-    S = CuArray{FT}(undef, 1, 1, grid.Tz)
-
-    HorizontalAverages(U, V, T, S)
-end
-
-HorizontalAverages(model) = HorizontalAverages(model.arch, model.grid)
+    # The Model constructor needs N and L.
+    kwargs[:N] = (kwargs[:grid].Nx, kwargs[:grid].Ny, kwargs[:grid].Nz)
+    kwargs[:L] = (kwargs[:grid].Lx, kwargs[:grid].Ly, kwargs[:grid].Lz)
 
-"""
-    VerticalPlanes(arch, grid)
-
-Instantiate a struct of 2D arrays on `arch` and `grid`
-for storing the results of y-averages of 3D fields.
-"""
-struct VerticalPlanes{A}
-    U :: A
-    V :: A
-    T :: A
-    S :: A
-end
-
-function VerticalPlanes(arch::CPU, grid::Grid{FT}) where FT
-    U = zeros(FT, grid.Tx, 1, grid.Tz)
-    V = zeros(FT, grid.Tx, 1, grid.Tz)
-    T = zeros(FT, grid.Tx, 1, grid.Tz)
-    S = zeros(FT, grid.Tx, 1, grid.Tz)
+    model = Model(; kwargs...)
 
-    VerticalPlanes(U, V, T, S)
-end
+    # Now restore fields.
+    for p in cps
+        if p in has_array_refs
+            restore_fields!(model, file, model.arch, p)
+        end
+    end
 
-function VerticalPlanes(arch::GPU, grid::Grid{FT}) where FT
-    U = CuArray{FT}(undef, grid.Tx, 1, grid.Tz)
-    V = CuArray{FT}(undef, grid.Tx, 1, grid.Tz)
-    T = CuArray{FT}(undef, grid.Tx, 1, grid.Tz)
-    S = CuArray{FT}(undef, grid.Tx, 1, grid.Tz)
+    close(file)
 
-    VerticalPlanes(U, V, T, S)
+    return model
 end
-
-VerticalPlanes(model) = VerticalPlanes(model.arch, model.grid)
diff --git a/test/test_output_writers.jl b/test/test_output_writers.jl
index 90b2e08f33..07699c6f7d 100644
--- a/test/test_output_writers.jl
+++ b/test/test_output_writers.jl
@@ -39,10 +39,10 @@ end
 
 function run_jld2_file_splitting_tests(arch)
     model = Model(N=(16, 16, 16), L=(1, 1, 1))
-    
+
     u(model) = Array(model.velocities.u.data.parent)
     fields = Dict(:u => u)
-    
+
     function fake_bc_init(file, model)
         file["boundary_conditions/fake"] = π
     end
@@ -51,7 +51,7 @@ function run_jld2_file_splitting_tests(arch)
                           init=fake_bc_init, including=[:grid],
                           max_filesize=200KiB, force=true)
     push!(model.output_writers, ow)
-    
+
     # 531 KiB of output will be written which should get split into 3 files.
     time_step!(model, 10, 1)
 
@@ -62,10 +62,10 @@ function run_jld2_file_splitting_tests(arch)
 
     for n in string.(1:3)
         filename = "test_part" * n * ".jld2"
-        jldopen(filename, "r") do file    
+        jldopen(filename, "r") do file
             # Test to make sure all files contain structs from `including`.
             @test file["grid/Nx"] == 16
-    
+
             # Test to make sure all files contain info from `init` function.
             @test file["boundary_conditions/fake"] == π
         end
@@ -75,18 +75,72 @@ function run_jld2_file_splitting_tests(arch)
     end
 end
 
+"""
+Run two coarse rising thermal bubble simulations and make sure that when
+restarting from a checkpoint, the restarted simulation matches the non-restarted
+simulation numerically.
+"""
+function run_thermal_bubble_checkpointer_tests(arch)
+    Nx, Ny, Nz = 16, 16, 16
+    Lx, Ly, Lz = 100, 100, 100
+    Δt = 6
+
+    true_model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), ν=4e-2, κ=4e-2, arch=arch)
+
+    # Add a cube-shaped warm temperature anomaly that takes up the middle 50%
+    # of the domain volume.
+    i1, i2 = round(Int, Nx/4), round(Int, 3Nx/4)
+    j1, j2 = round(Int, Ny/4), round(Int, 3Ny/4)
+    k1, k2 = round(Int, Nz/4), round(Int, 3Nz/4)
+    true_model.tracers.T.data[i1:i2, j1:j2, k1:k2] .+= 0.01
+
+    checkpointed_model = deepcopy(true_model)
+
+    time_step!(true_model, 9, Δt)
+
+    checkpointer = Checkpointer(checkpointed_model; frequency=5, force=true)
+    push!(checkpointed_model.output_writers, checkpointer)
+
+    # Checkpoint should be saved as "checkpoint5.jld" after the 5th iteration.
+    time_step!(checkpointed_model, 5, Δt)
+
+    # Remove all knowledge of the checkpointed model.
+    checkpointed_model = nothing
+
+    restored_model = restore_from_checkpoint("checkpoint5.jld2")
+
+    time_step!(restored_model, 4, Δt; init_with_euler=false)
+
+    # Now the true_model and restored_model should be identical.
+    @test all(restored_model.velocities.u.data      .≈ true_model.velocities.u.data)
+    @test all(restored_model.velocities.v.data      .≈ true_model.velocities.v.data)
+    @test all(restored_model.velocities.w.data      .≈ true_model.velocities.w.data)
+    @test all(restored_model.tracers.T.data         .≈ true_model.tracers.T.data)
+    @test all(restored_model.tracers.S.data         .≈ true_model.tracers.S.data)
+    @test all(restored_model.timestepper.Gⁿ.Gu.data .≈ true_model.timestepper.Gⁿ.Gu.data)
+    @test all(restored_model.timestepper.Gⁿ.Gv.data .≈ true_model.timestepper.Gⁿ.Gv.data)
+    @test all(restored_model.timestepper.Gⁿ.Gw.data .≈ true_model.timestepper.Gⁿ.Gw.data)
+    @test all(restored_model.timestepper.Gⁿ.GT.data .≈ true_model.timestepper.Gⁿ.GT.data)
+    @test all(restored_model.timestepper.Gⁿ.GS.data .≈ true_model.timestepper.Gⁿ.GS.data)
+end
+
 @testset "Output writers" begin
     println("Testing output writers...")
 
     for arch in archs
-        @testset "NetCDF [$(typeof(arch))]" begin
-            println("  Testing NetCDF output writer [$(typeof(arch))]...")
-            run_thermal_bubble_netcdf_tests(arch)
-        end
+         @testset "NetCDF [$(typeof(arch))]" begin
+             println("  Testing NetCDF output writer [$(typeof(arch))]...")
+             run_thermal_bubble_netcdf_tests(arch)
+         end
 
         @testset "JLD2 [$(typeof(arch))]" begin
             println("  Testing JLD2 output writer [$(typeof(arch))]...")
             run_jld2_file_splitting_tests(arch)
         end
+
+        @testset "Checkpointer [$(typeof(arch))]" begin
+            println("  Testing Checkpointer [$(typeof(arch))]...")
+            run_thermal_bubble_checkpointer_tests(arch)
+        end
     end
 end