diff --git a/Project.toml b/Project.toml
index 295aaf46a..998a6f71a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -49,7 +49,7 @@ JLD2 = "0.4, 0.5"
 KernelAbstractions = "0.9"
 MPI = "0.20"
 NCDatasets = "0.12, 0.13, 0.14"
-Oceananigans = "0.96 - 0.99"
+Oceananigans = "0.96.1 - 0.99"
 OffsetArrays = "1.14"
 PrecompileTools = "1"
 Reactant = "0.2.45"
diff --git a/examples/checkpointer_mwe.jl b/examples/checkpointer_mwe.jl
new file mode 100755
index 000000000..a0007fa0b
--- /dev/null
+++ b/examples/checkpointer_mwe.jl
@@ -0,0 +1,72 @@
+using ClimaOcean
+using Oceananigans
+using CFTime
+using Dates
+using Printf
+
+arch = CPU()
+
+Nx = 120
+Ny = 50
+Nz = 20
+
+grid = LatitudeLongitudeGrid(arch;
+                             size = (Nx, Ny, Nz),
+                             halo = (7, 7, 7),
+                             z = (-6000, 0),
+                             latitude  = (-75, 75),
+                             longitude = (0, 360))
+
+
+function testbed_coupled_simulation(grid; stop_iteration=8)
+    ocean = ocean_simulation(grid)
+        
+    radiation = Radiation(arch)
+    
+    atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(4))
+    
+    coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation)
+    
+    return Simulation(coupled_model; Δt=10, stop_iteration)
+end
+
+simulation = testbed_coupled_simulation(grid; stop_iteration=8)
+
+wall_time = Ref(time_ns())
+
+function progress(sim)
+    ocean = sim.model.ocean
+    atmosphere = sim.model.atmosphere
+
+    step_time = 1e-9 * (time_ns() - wall_time[])
+
+    msg = @sprintf("iteration: %d, sim time: %s, atmos time: %s, ocean time: %s, wall time: %s",
+                   iteration(sim), sim.model.clock.time, atmosphere.clock.time, ocean.model.clock.time, prettytime(step_time))
+
+    @info msg
+
+    wall_time[] = time_ns()
+end
+
+simulation.callbacks[:progress] = Callback(progress, IterationInterval(1))
+
+simulation.output_writers[:checkpoint] = Checkpointer(ocean.model;
+                                                      schedule = IterationInterval(3),
+                                                      prefix = "checkpointer",
+                                                      dir = ".",
+                                                      verbose = true,
+                                                      overwrite_existing = true)
+
+run!(simulation)
+
+new_simulation = testbed_coupled_simulation(grid; stop_iteration=14)
+
+new_simulation.output_writers[:checkpoint] = Checkpointer(ocean.model;
+                                                      schedule = IterationInterval(3),
+                                                      prefix = "checkpointer",
+                                                      dir = ".",
+                                                      verbose = true,
+                                                      overwrite_existing = true)
+
+
+run!(new_simulation, pickup=true)
diff --git a/examples/one_degree_simulation.jl b/examples/one_degree_simulation.jl
index b0483bb5c..fa787d87a 100644
--- a/examples/one_degree_simulation.jl
+++ b/examples/one_degree_simulation.jl
@@ -157,6 +157,14 @@ ocean.output_writers[:surface] = JLD2Writer(ocean.model, outputs;
                                             indices = (:, :, grid.Nz),
                                             overwrite_existing = true)
 
+# We also save a series of Checkpointers that will enable us to restart the simulation
+# from a later point; ideal for longer runs.
+
+simulation.output_writers[:checkpoint] = Checkpointer(ocean.model;
+                                                      schedule = TimeInterval(5days),
+                                                      prefix = "one_degree_checkpointer",
+                                                      overwrite_existing = true)
+
 # ### Ready to run
 
 # We are ready to press the big red button and run the simulation.
@@ -170,6 +178,16 @@ simulation.Δt = 20minutes
 simulation.stop_time = 360days
 run!(simulation)
 
+# After we had run long-enough to produce at least one checkpointer file, if we want to restart the
+# simulation we call
+#
+# ```julia
+# run!(simulation, pickup=true)
+# ```
+#
+# Doing so, simulation will be picked up from the latest checkpointer.
+
+
 # ### A pretty movie
 #
 # We load the saved output and make a pretty movie of the simulation. First we plot a snapshot:
diff --git a/src/ClimaOcean.jl b/src/ClimaOcean.jl
index 586fb678f..d7c77ff4b 100644
--- a/src/ClimaOcean.jl
+++ b/src/ClimaOcean.jl
@@ -108,4 +108,3 @@ using PrecompileTools: @setup_workload, @compile_workload
 end
 
 end # module
-
diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl
index 73e1203ba..746636c72 100644
--- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl
+++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl
@@ -418,4 +418,3 @@ function sea_ice_similarity_theory(sea_ice::SeaIceSimulation)
     interface_temperature_type = SkinTemperature(internal_flux)
     return SimilarityTheoryFluxes(; interface_temperature_type)
 end
-
diff --git a/src/OceanSeaIceModels/OceanSeaIceModels.jl b/src/OceanSeaIceModels/OceanSeaIceModels.jl
index 54d43ead0..170b805d7 100644
--- a/src/OceanSeaIceModels/OceanSeaIceModels.jl
+++ b/src/OceanSeaIceModels/OceanSeaIceModels.jl
@@ -26,6 +26,11 @@ using ClimaOcean: stateindex
 using KernelAbstractions: @kernel, @index
 using KernelAbstractions.Extras.LoopInfo: @unroll
 
+import Oceananigans.OutputWriters: checkpointer_address,
+                                   required_checkpointed_properties,
+                                   default_checkpointed_properties,
+                                   prognostic_fields
+
 function downwelling_radiation end
 function freshwater_flux end
 function reference_density end
@@ -69,6 +74,7 @@ import .InterfaceComputations:
     compute_sea_ice_ocean_fluxes!
 
 include("ocean_sea_ice_model.jl")
+include("output_writers.jl")
 include("freezing_limited_ocean_temperature.jl")
 include("time_step_ocean_sea_ice_model.jl")
 
diff --git a/src/OceanSeaIceModels/PrescribedAtmospheres.jl b/src/OceanSeaIceModels/PrescribedAtmospheres.jl
index 967317222..9a85fdaac 100644
--- a/src/OceanSeaIceModels/PrescribedAtmospheres.jl
+++ b/src/OceanSeaIceModels/PrescribedAtmospheres.jl
@@ -7,9 +7,15 @@ using Oceananigans.OutputReaders: FieldTimeSeries, update_field_time_series!, ex
 using Oceananigans.TimeSteppers: Clock, tick!
 
 using Adapt
+using JLD2
 using Thermodynamics.Parameters: AbstractThermodynamicsParameters
 
-import Oceananigans.TimeSteppers: time_step!
+import Oceananigans: initialize!, prognostic_fields
+import Oceananigans.Fields: set!
+import Oceananigans.TimeSteppers: time_step!, set_clock!
+import Oceananigans.OutputWriters: checkpointer_address,
+                                   required_checkpointed_properties,
+                                   default_checkpointed_properties
 
 import Thermodynamics.Parameters:
     gas_constant,   #
@@ -22,7 +28,7 @@ import Thermodynamics.Parameters:
     T_0,            # Enthalpy reference temperature
     LH_v0,          # Vaporization enthalpy at the reference temperature
     LH_s0,          # Sublimation enthalpy at the reference temperature
-    LH_f0,          # Fusionn enthalpy at the reference temperature
+    LH_f0,          # Fusion enthalpy at the reference temperature
     cp_d,           # Heat capacity of dry air at constant pressure
     cp_v,           # Isobaric specific heat capacity of gaseous water vapor
     cp_l,           # Isobaric specific heat capacity of liquid water
@@ -30,7 +36,7 @@ import Thermodynamics.Parameters:
     cv_v,           # Heat capacity of dry air at constant volume
     cv_l,           # Isobaric specific heat capacity of liquid water
     cv_i,           # Isobaric specific heat capacity of liquid water
-    e_int_v0,       # what? someting about reference internal energy of water vapor
+    e_int_v0,       # what? something about reference internal energy of water vapor
     T_freeze,       # Freezing temperature of _pure_ water
     T_triple,       # Triple point temperature of _pure_ water
     press_triple,   # Triple point pressure of pure water
@@ -75,7 +81,7 @@ Construct a set of parameters that define the density of moist air,
 ```
 
 where ``p`` is pressure, ``T`` is temperature, ``q`` defines the partition
-of total mass into vapor, liqiud, and ice mass fractions, and
+of total mass into vapor, liquid, and ice mass fractions, and
 ``Rᵐ`` is the effective specific gas constant for the mixture,
 
 ```math
@@ -318,6 +324,29 @@ function Base.show(io::IO, pa::PrescribedAtmosphere)
     print(io, "└── boundary_layer_height: ", prettysummary(pa.boundary_layer_height))
 end
 
+checkpointer_address(::PrescribedAtmosphere)             = "PrescribedAtmosphere"
+default_checkpointed_properties(::PrescribedAtmosphere)  = [:clock]
+required_checkpointed_properties(::PrescribedAtmosphere) = [:clock]
+prognostic_fields(::PrescribedAtmosphere)                = ()
+
+function set!(model::PrescribedAtmosphere, checkpoint_file_path::AbstractString)
+    addr = checkpointer_address(model)
+
+    checkpointed_clock = nothing
+    jldopen(checkpoint_file_path, "r") do file
+        checkpointed_clock = file["$addr/clock"]
+    end
+
+    # Update model clock
+    set_clock!(model, checkpointed_clock)
+
+    @info "hi, I set the clock of the atmosphere"
+    return nothing
+end
+
+set_clock!(atmos::PrescribedAtmosphere, new_clock) =
+    set_clock!(atmos.clock, new_clock)
+
 function default_atmosphere_velocities(grid, times)
     ua = FieldTimeSeries{Center, Center, Nothing}(grid, times)
     va = FieldTimeSeries{Center, Center, Nothing}(grid, times)
@@ -327,7 +356,7 @@ end
 function default_atmosphere_tracers(grid, times)
     Ta = FieldTimeSeries{Center, Center, Nothing}(grid, times)
     qa = FieldTimeSeries{Center, Center, Nothing}(grid, times)
-    parent(Ta) .= 273.15 + 20
+    parent(Ta) .= 273.15 + 20 # ᵒK
     return (T=Ta, q=qa)
 end
 
@@ -347,13 +376,20 @@ end
 This is approximately equal to the mean sea-level atmospheric pressure on Earth. """
 function default_atmosphere_pressure(grid, times)
     pa = FieldTimeSeries{Center, Center, Nothing}(grid, times)
-    parent(pa) .= 101325
+    parent(pa) .= 101325 # Pa
     return pa
 end
 
 @inline function time_step!(atmos::PrescribedAtmosphere, Δt)
     tick!(atmos.clock, Δt)
+    update_state!(atmos)
+
+    return nothing
+end
+
+initialize!(atmos::PrescribedAtmosphere) = update_state!(atmos)
 
+function update_state!(atmos::PrescribedAtmosphere)
     time = Time(atmos.clock.time)
     ftses = extract_field_time_series(atmos)
 
@@ -422,7 +458,7 @@ struct TwoBandDownwellingRadiation{SW, LW}
 end
 
 """
-    TwoBandDownwellingRadiation(shortwave=nothing, longwave=nothing)
+    TwoBandDownwellingRadiation(; shortwave=nothing, longwave=nothing)
 
 Return a two-band model for downwelling radiation (split into a shortwave band
 and a longwave band) that passes through the atmosphere and arrives at the surface of ocean
diff --git a/src/OceanSeaIceModels/ocean_sea_ice_model.jl b/src/OceanSeaIceModels/ocean_sea_ice_model.jl
index f7fc690c4..c8ea29d06 100644
--- a/src/OceanSeaIceModels/ocean_sea_ice_model.jl
+++ b/src/OceanSeaIceModels/ocean_sea_ice_model.jl
@@ -3,21 +3,25 @@ using Oceananigans.TimeSteppers: Clock
 using Oceananigans: SeawaterBuoyancy
 using ClimaSeaIce.SeaIceThermodynamics: melting_temperature
 using KernelAbstractions: @kernel, @index
-
 using SeawaterPolynomials: TEOS10EquationOfState
 
 import Thermodynamics as AtmosphericThermodynamics
 
 # Simulations interface
-import Oceananigans: fields, prognostic_fields
+import Oceananigans: fields,
+                     prognostic_fields
 import Oceananigans.Architectures: architecture
 import Oceananigans.Fields: set!
-import Oceananigans.Models: timestepper, NaNChecker, default_nan_checker
-import Oceananigans.OutputWriters: default_included_properties
-import Oceananigans.Simulations: reset!, initialize!, iteration
-import Oceananigans.TimeSteppers: time_step!, update_state!, time
-import Oceananigans.Utils: prettytime
 import Oceananigans.Models: timestepper, NaNChecker, default_nan_checker, initialization_update_state!
+import Oceananigans.OutputWriters: checkpointer_address,
+                                   required_checkpointed_properties,
+                                   default_checkpointed_properties
+
+import Oceananigans.Simulations: reset!, initialize!, iteration, run!
+import Oceananigans.TimeSteppers: time_step!, update_state!, time, set_clock!
+import Oceananigans.Utils: prettytime
+
+import .PrescribedAtmospheres: set_clock!
 
 mutable struct OceanSeaIceModel{I, A, O, F, C, Arch} <: AbstractModel{Nothing, Arch}
     architecture :: Arch
@@ -29,6 +33,7 @@ mutable struct OceanSeaIceModel{I, A, O, F, C, Arch} <: AbstractModel{Nothing, A
 end
 
 const OSIM = OceanSeaIceModel
+const OSIMPA = OceanSeaIceModel{<:Any, <:PrescribedAtmosphere}
 
 function Base.summary(model::OSIM)
     A = nameof(typeof(architecture(model)))
@@ -53,20 +58,23 @@ function Base.show(io::IO, cm::OSIM)
 end
 
 # Assumption: We have an ocean!
-architecture(model::OSIM)           = architecture(model.ocean.model)
-Base.eltype(model::OSIM)            = Base.eltype(model.ocean.model)
-prettytime(model::OSIM)             = prettytime(model.clock.time)
-iteration(model::OSIM)              = model.clock.iteration
-timestepper(::OSIM)                 = nothing
-default_included_properties(::OSIM) = tuple()
-prognostic_fields(cm::OSIM)         = nothing
-fields(::OSIM)                      = NamedTuple()
-default_clock(TT)                   = Oceananigans.TimeSteppers.Clock{TT}(0, 0, 1)
-
-function reset!(model::OSIM)
-    reset!(model.ocean)
-    return nothing
-end
+architecture(model::OSIM)                = architecture(model.ocean.model)
+Base.eltype(model::OSIM)                 = Base.eltype(model.ocean.model)
+prettytime(model::OSIM)                  = prettytime(model.clock.time)
+iteration(model::OSIM)                   = model.clock.iteration
+timestepper(::OSIM)                      = nothing
+default_included_properties(::OSIM)      = tuple()
+prognostic_fields(::OSIM)                = tuple()
+fields(::OSIM)                           = NamedTuple()
+default_clock(TT)                        = Oceananigans.TimeSteppers.Clock{TT}(0, 0, 1)
+time(model::OSIM)                        = model.clock.time
+required_checkpointed_properties(::OSIM) = [:clock]
+default_checkpointed_properties(::OSIM)  = [:clock]
+checkpointer_address(::OSIM)             = "OceanSeaIceModel"
+
+reset!(model::OSIM) = reset!(model.ocean)
+
+components(model::OSIM) = (model.atmosphere, model.ocean.model)
 
 # Make sure to initialize the exchanger here
 function initialization_update_state!(model::OSIM)
@@ -77,10 +85,25 @@ end
 
 function initialize!(model::OSIM)
     initialize!(model.ocean)
+    initialize!(model.atmosphere)
     initialize!(model.interfaces.exchanger, model.atmosphere)
     return nothing
 end
 
+function set_clock!(model::OSIM, clock)
+
+    set_clock!(model.clock, clock)
+
+    # set the component clocks
+    atmos, ocean = model.atmosphere, model.ocean
+    set_clock!(atmos, clock)
+    set_clock!(ocean, clock)
+
+    return nothing
+end
+
+set_clock!(sim::Simulation, new_clock) = set_clock!(sim.model, new_clock)
+
 reference_density(unsupported) =
     throw(ArgumentError("Cannot extract reference density from $(typeof(unsupported))"))
 
@@ -159,10 +182,8 @@ function OceanSeaIceModel(ocean, sea_ice=FreezingLimitedOceanTemperature(eltype(
     return ocean_sea_ice_model
 end
 
-time(coupled_model::OceanSeaIceModel) = coupled_model.clock.time
-
 # Check for NaNs in the first prognostic field (generalizes to prescribed velocities).
-function default_nan_checker(model::OceanSeaIceModel)
+function default_nan_checker(model::OSIM)
     u_ocean = model.ocean.model.velocities.u
     nan_checker = NaNChecker((; u_ocean))
     return nan_checker
diff --git a/src/OceanSeaIceModels/output_writers.jl b/src/OceanSeaIceModels/output_writers.jl
new file mode 100644
index 000000000..9144e2e39
--- /dev/null
+++ b/src/OceanSeaIceModels/output_writers.jl
@@ -0,0 +1,62 @@
+using JLD2
+
+using Oceananigans.Utils: pretty_filesize
+using Oceananigans.OutputWriters: checkpoint_path,
+                                  cleanup_checkpoints,
+                                  validate_checkpointed_properties
+
+import Oceananigans.Fields: set!
+import Oceananigans.OutputWriters: write_output!
+
+wall_time = Ref(time_ns())
+
+function set!(model::OSIM, checkpoint_file_path::AbstractString)
+    addr = checkpointer_address(model)
+
+    checkpointed_clock = nothing
+    jldopen(checkpoint_file_path, "r") do file
+        checkpointed_clock = file["$addr/clock"]
+    end
+
+    @show checkpointed_clock
+
+    # Update model clock
+    set_clock!(model, checkpointed_clock)
+
+    # deal with model components
+
+    for component in components(model)
+        set!(component, checkpoint_file_path)
+    end
+
+    return nothing
+end
+
+
+function write_output!(c::Checkpointer, model::OSIM)
+
+    filepath = checkpoint_path(model.clock.iteration, c)
+    c.verbose && @info "Checkpointing to file $filepath..."
+
+    t1 = time_ns()
+
+    # write OceanSeaIceModel
+    mode = "w"
+    write_output!(c, model, filepath, mode,
+                  properties = default_checkpointed_properties(model))
+
+    # write the OceanSeaIceModel components
+    mode = "a" # to append in the file already written above
+    for component in components(model)
+        write_output!(c, component, filepath, mode,
+                      properties = default_checkpointed_properties(component))
+    end
+
+    t2, sz = time_ns(), filesize(filepath)
+
+    c.verbose && @info "Checkpointing done: time=$(prettytime((t2 - t1) * 1e-9)), size=$(pretty_filesize(sz))"
+
+    c.cleanup && cleanup_checkpoints(c)
+
+    return nothing
+end
diff --git a/src/OceanSimulations/ocean_simulation.jl b/src/OceanSimulations/ocean_simulation.jl
index 4d0379884..ac6a2c109 100644
--- a/src/OceanSimulations/ocean_simulation.jl
+++ b/src/OceanSimulations/ocean_simulation.jl
@@ -229,4 +229,3 @@ end
 
 hasclosure(closure, ClosureType) = closure isa ClosureType
 hasclosure(closure_tuple::Tuple, ClosureType) = any(hasclosure(c, ClosureType) for c in closure_tuple)
-
diff --git a/test/runtests.jl b/test/runtests.jl
index 6376a22ac..ede1be784 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,6 +5,7 @@ using CUDA
 
 test_group = get(ENV, "TEST_GROUP", :all)
 test_group = Symbol(test_group)
+test_group = :ocean_sea_ice_model
 
 using ClimaOcean.DataWrangling: download_dataset
 
diff --git a/test/runtests_setup.jl b/test/runtests_setup.jl
index 448de8da5..2fde40b76 100644
--- a/test/runtests_setup.jl
+++ b/test/runtests_setup.jl
@@ -1,6 +1,9 @@
 using ClimaOcean
 using Oceananigans
+
 using CUDA
+using CFTime
+using Dates
 using Test
 
 using ClimaOcean.Bathymetry: download_bathymetry_cache
@@ -12,9 +15,6 @@ using ClimaOcean.JRA55
 using Oceananigans.Architectures: architecture, on_architecture
 using Oceananigans.OutputReaders: interpolate!
 
-using CFTime
-using Dates
-
 using CUDA: @allowscalar
 
 gpu_test = parse(Bool, get(ENV, "GPU_TEST", "false"))
diff --git a/test/test_ocean_sea_ice_model.jl b/test/test_ocean_sea_ice_model.jl
index 24b5fa75e..a53462487 100644
--- a/test/test_ocean_sea_ice_model.jl
+++ b/test/test_ocean_sea_ice_model.jl
@@ -1,8 +1,11 @@
 include("runtests_setup.jl")
 
 using CUDA
+using Oceananigans.Units
+using Oceananigans: prognostic_fields
+using Oceananigans.TimeSteppers: QuasiAdamsBashforth2TimeStepper
 using Oceananigans.OrthogonalSphericalShellGrids
-using ClimaOcean.OceanSeaIceModels: above_freezing_ocean_temperature!
+using ClimaOcean.OceanSeaIceModels: above_freezing_ocean_temperature!, OceanSeaIceModel
 using ClimaSeaIce.SeaIceThermodynamics: melting_temperature
 using ClimaSeaIce.SeaIceMomentumEquations
 using ClimaSeaIce.Rheologies
@@ -57,9 +60,9 @@ using ClimaSeaIce.Rheologies
                                 z = (-6000, 0))
 
             bottom_height = regrid_bathymetry(grid;
-                                              minimum_depth = 10,
-                                              interpolation_passes = 20,
-                                              major_basins = 1)
+                                            minimum_depth = 10,
+                                            interpolation_passes = 20,
+                                            major_basins = 1)
 
             grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true)
 
@@ -86,10 +89,10 @@ using ClimaSeaIce.Rheologies
             τva = Field{Center, Face, Nothing}(grid)
 
             dynamics = SeaIceMomentumEquation(grid;
-                                              coriolis = ocean.model.coriolis,
-                                              top_momentum_stress = (u=τua, v=τva),
-                                              rheology = ElastoViscoPlasticRheology(),
-                                              solver = SplitExplicitSolver(120))
+                                            coriolis = ocean.model.coriolis,
+                                            top_momentum_stress = (u=τua, v=τva),
+                                            rheology = ElastoViscoPlasticRheology(),
+                                            solver = SplitExplicitSolver(120))
 
             sea_ice  = sea_ice_simulation(grid; dynamics, advection=WENO(order=7))
             liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus
@@ -115,3 +118,210 @@ using ClimaSeaIce.Rheologies
         end
     end
 end
+
+"""
+    test_ocean_sea_ice_model_checkpointer(grid)
+
+Return two ocean-sea ice coupled models to be used for testing purposes.
+"""
+function test_ocean_sea_ice_model_checkpointer(grid)
+    ocean = ocean_simulation(grid, closure=nothing, free_surface = ExplicitFreeSurface())
+    atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(4))
+
+    true_coupled_model = OceanSeaIceModel(ocean; atmosphere)
+    # true_coupled_model = OceanSeaIceModel(ocean; atmosphere=nothing, radiation=nothing)
+
+    test_coupled_model = deepcopy(true_coupled_model)
+
+    return true_coupled_model, test_coupled_model
+end
+
+function test_model_equality(test_model::OceanSeaIceModel, true_model::OceanSeaIceModel)
+    test_clock_equality(test_model.clock, true_model.clock)
+
+    # test equality for components
+    test_model_equality(test_model.atmosphere, true_model.atmosphere)
+    test_model_equality(test_model.ocean.model, true_model.ocean.model)
+
+    return nothing
+end
+
+test_model_equality(test_model::PrescribedAtmosphere, true_model::PrescribedAtmosphere) =
+    test_clock_equality(test_model.clock, true_model.clock)
+
+function test_clock_equality(clock1::Clock, clock2::Clock)
+    @test clock1.iteration == clock2.iteration
+    @test clock1.time == clock2.time
+    @test clock1.last_Δt == clock2.last_Δt
+    return nothing
+end
+
+function test_model_equality(test_model::HydrostaticFreeSurfaceModel, true_model::HydrostaticFreeSurfaceModel)
+
+    test_clock_equality(test_model.clock, true_model.clock)
+
+    CUDA.@allowscalar begin
+        test_model_fields = prognostic_fields(test_model)
+        true_model_fields = prognostic_fields(true_model)
+        field_names = keys(test_model_fields)
+
+        for name in field_names
+            @show name
+            @show all(test_model_fields[name].data .≈ true_model_fields[name].data)
+
+            if test_model.timestepper isa QuasiAdamsBashforth2TimeStepper
+                if name ∈ keys(test_model.timestepper.Gⁿ)
+                    @show name
+                    # @show test_model.timestepper.Gⁿ[name]
+                    # @show true_model.timestepper.Gⁿ[name]
+                    @show all(test_model.timestepper.Gⁿ[name].data .≈ true_model.timestepper.Gⁿ[name].data)
+                    @show all(test_model.timestepper.G⁻[name].data .≈ true_model.timestepper.G⁻[name].data)
+                end
+            end
+        end
+    end
+
+    return nothing
+end
+
+function run_checkpointer_tests(true_model, test_model, Δt)
+
+    true_simulation = Simulation(true_model; Δt, stop_iteration = 4)
+    true_simulation.output_writers[:checkpoint] =
+        Checkpointer(true_model, schedule=IterationInterval(4), overwrite_existing=true)
+
+    run!(true_simulation)
+
+    checkpointed_model = deepcopy(true_simulation.model)
+
+    true_simulation.stop_iteration = 7
+    run!(true_simulation)
+
+    #####
+    ##### Test `set!(model, checkpoint_file)`
+    #####
+
+    @info "Testing `set!(model, checkpoint_file)`"
+
+    set!(test_model, "checkpoint_iteration4.jld2")
+
+    test_model_equality(test_model, checkpointed_model)
+
+    #=
+    #####
+    ##### Test pickup from explicit checkpoint path
+    #####
+
+    @info "Testing pickup from explicit checkpoint path"
+
+    test_simulation = Simulation(test_model; Δt, stop_iteration=7)
+    run!(test_simulation, pickup="checkpoint_iteration0.jld2")
+
+    test_model_equality(test_model, true_model)
+
+    run!(test_simulation, pickup="checkpoint_iteration4.jld2")
+
+    @info "Testing model equality when running with pickup=checkpoint_iteration4.jld2."
+
+    test_model_equality(test_model, true_model)
+
+    @test test_simulation.model.clock.iteration == true_simulation.model.clock.iteration
+    @test test_simulation.model.clock.time == true_simulation.model.clock.time
+    test_model_equality(test_simulation.model, true_simulation.model)
+
+    #####
+    ##### Test `run!(sim, pickup=true)
+    #####
+
+    @info "Testing run!(sim, pickup=true)"
+
+    # Pickup using existing checkpointer
+    test_simulation.output_writers[:checkpoint] =
+        Checkpointer(test_model, schedule=IterationInterval(4), overwrite_existing=true)
+
+    run!(test_simulation, pickup=true)
+    @info "    Testing model equality when running with pickup=true."
+
+    @test test_simulation.model.clock.iteration == true_simulation.model.clock.iteration
+    @test test_simulation.model.clock.time == true_simulation.model.clock.time
+    test_model_equality(test_simulation.model, true_simulation.model)
+
+    run!(test_simulation, pickup=0)
+    @info "    Testing model equality when running with pickup=0."
+
+    @test test_simulation.model.clock.iteration == true_simulation.model.clock.iteration
+    @test test_simulation.model.clock.time == true_simulation.model.clock.time
+    test_model_equality(test_simulation.model, true_simulation.model)
+
+    run!(test_simulation, pickup=4)
+    @info "    Testing model equality when running with pickup=4."
+
+    @test test_simulation.model.clock.iteration == true_simulation.model.clock.iteration
+    @test test_simulation.model.clock.time == true_simulation.model.clock.time
+    test_model_equality(test_simulation.model, true_simulation.model)
+
+    =#
+
+    for file in Tuple("checkpoint_iteration$i.jld2" for i in 0:4:10)
+        rm(file, force=true)
+    end
+end
+
+#=
+@testset "Checkpointer tests" begin
+    for arch in test_architectures
+
+        Δt = 1.2hours
+
+        Nx, Ny, Nz = 40, 25, 10
+
+        grid = LatitudeLongitudeGrid(arch;
+                                     size = (Nx, Ny, Nz),
+                                     halo = (7, 7, 7),
+                                     z = (-6000, 0),
+                                     latitude  = (-75, 75),
+                                     longitude = (0, 360))
+
+        true_coupled_model, test_coupled_model = test_ocean_sea_ice_model_checkpointer(grid)
+
+        run_checkpointer_tests(true_coupled_model, test_coupled_model, Δt)
+
+        #=
+        true_simulation = testbed_coupled_simulation(true_coupled_model; Δt, checkpoint_iteration_interval,
+                                                     stop_iteration=intermediate_iteration)
+
+        test_simulation = testbed_coupled_simulation(test_coupled_model; Δt, checkpoint_iteration_interval
+                                                     stop_iteration = final_iteration)
+
+        run!(true_simulation) # up to intermediate_iteration
+
+        checkpointed_model = deepcopy(true_simulation.model)
+
+        set!(test_model, "checkpoint_iteration5.jld2")
+
+        @test test_model.clock.iteration == checkpointed_model.clock.iteration
+        @test test_model.clock.time == checkpointed_model.clock.time
+
+        run!(new_simulation, pickup=true)
+
+        # ensure all clocks all at same time and iteration
+        for clock in (simulation.model.clock,
+                      simulation.model.atmosphere.clock,
+                      simulation.model.ocean.model.clock)
+
+            @test clock.iteration ≈ intermediate_iteration
+            @test clock.time ≈ intermediate_iteration * Δt
+        end
+
+        for clock in (new_simulation.model.clock,
+                      new_simulation.model.atmosphere.clock,
+                      new_simulation.model.ocean.model.clock)
+
+            @test clock.iteration ≈ final_iteration
+            @test clock.time ≈ final_iteration * Δt
+        end
+        =#
+    end
+
+end
+=#