diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl
index 6538434c9..449affd31 100644
--- a/experiments/arctic_simulation.jl
+++ b/experiments/arctic_simulation.jl
@@ -81,7 +81,7 @@ dynamics = SeaIceMomentumEquation(grid;
                                   rheology = ElastoViscoPlasticRheology(),
                                   solver = SplitExplicitSolver(120))
 
-sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition, dynamics, advection=WENO(order=7)) 
+sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition) #, dynamics, advection=WENO(order=7)) 
 
 set!(sea_ice.model, h=Metadata(:sea_ice_thickness;     dataset), 
                     ℵ=Metadata(:sea_ice_concentration; dataset))
diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl
index 94a1a951e..f0f766663 100644
--- a/experiments/flux_climatology/flux_climatology.jl
+++ b/experiments/flux_climatology/flux_climatology.jl
@@ -115,32 +115,6 @@ end
 ##### A prescribed ocean...
 #####
 
-struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing}
-    architecture :: A       
-    grid :: G        
-    clock :: Clock{C}
-    velocities :: U
-    tracers :: T
-    timeseries :: F
-end
-
-function PrescribedOcean(timeseries; 
-                         grid, 
-                         clock=Clock{Float64}(time = 0)) 
-
-    τx = Field{Face, Center, Nothing}(grid)
-    τy = Field{Center, Face, Nothing}(grid)
-    Jᵀ = Field{Center, Center, Nothing}(grid)
-    Jˢ = Field{Center, Center, Nothing}(grid)
-
-    u = XFaceField(grid,  boundary_conditions=FieldBoundaryConditions(grid, (Face,   Center, Center), top = FluxBoundaryCondition(τx)))
-    v = YFaceField(grid,  boundary_conditions=FieldBoundaryConditions(grid, (Center, Face,   Center), top = FluxBoundaryCondition(τy)))
-    T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ)))
-    S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ)))
-
-    PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries)
-end
-
 # ...with prescribed velocity and tracer fields
 version = ECCO4Monthly()
 dates   = all_ECCO_dates(version)[1:24]
@@ -157,35 +131,6 @@ grid = ECCO.ECCO_immersed_grid(arch)
 ocean_model = PrescribedOcean((; u, v, T, S); grid)
 ocean = Simulation(ocean_model, Δt=3hour, stop_time=365days)
 
-#####
-##### Need to extend a couple of methods
-#####
-
-function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true)
-    tick!(model.clock, Δt)
-    time = Time(model.clock.time)
-
-    possible_fts = merge(model.velocities, model.tracers)
-    time_series_tuple = extract_field_time_series(possible_fts)
-
-    for fts in time_series_tuple
-        update_field_time_series!(fts, time)
-    end
-
-    parent(model.velocities.u) .= parent(model.timeseries.u[time])
-    parent(model.velocities.v) .= parent(model.timeseries.v[time])
-    parent(model.tracers.T)    .= parent(model.timeseries.T[time])
-    parent(model.tracers.S)    .= parent(model.timeseries.S[time])
-
-    return nothing
-end
-
-update_state!(::PrescribedOcean) = nothing
-timestepper(::PrescribedOcean) = nothing
-
-reference_density(ocean::Simulation{<:PrescribedOcean}) = 1025.6
-heat_capacity(ocean::Simulation{<:PrescribedOcean}) = 3995.6
-
 #####
 ##### A prescribed atmosphere...
 #####
diff --git a/experiments/single_column_sea_ice.jl b/experiments/single_column_sea_ice.jl
new file mode 100644
index 000000000..63cc91d1d
--- /dev/null
+++ b/experiments/single_column_sea_ice.jl
@@ -0,0 +1,114 @@
+using ClimaOcean
+using ClimaSeaIce
+using Oceananigans
+using Oceananigans.Grids
+using Oceananigans.Units
+using ClimaOcean.OceanSimulations
+using ClimaOcean.ECCO
+using ClimaOcean.DataWrangling
+using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium
+using Printf
+
+# A single column grid with one point in the Arctic ocean
+grid = LatitudeLongitudeGrid(size = 1, 
+                             latitude = 87.0,
+                             longitude = 33.0,
+                             z = (-10, 0),
+                             topology = (Flat, Flat, Bounded))
+
+#####
+##### A Prescribed Ocean model
+#####
+
+dataset = ECCO4Monthly()
+dates   = all_dates(dataset)[1:24]
+
+T = ECCOFieldTimeSeries(:temperature, grid; dates, inpainting=nothing, time_indices_in_memory=24)
+S = ECCOFieldTimeSeries(:salinity, grid;    dates, inpainting=nothing, time_indices_in_memory=24)
+
+ocean_model = PrescribedOceanModel((; T, S); grid)
+ocean = Simulation(ocean_model, Δt=30minutes, stop_time=730days)
+
+##### 
+##### A Prescribed Atmosphere model
+#####
+
+atmosphere = JRA55PrescribedAtmosphere(latitude = 87.0,
+                                       longitude = 33.0,
+                                       backend = InMemory())
+
+#####
+##### A Prognostic Sea-ice model with no dynamics
+#####
+
+# Remember to pass the SSS as a bottom bc to the sea ice!
+SSS = view(ocean.model.tracers.S, :, :, 1)
+bottom_heat_boundary_condition = IceWaterThermalEquilibrium(SSS)
+
+sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition) 
+
+set!(sea_ice.model, h=Metadata(:sea_ice_thickness;     dataset, dates=first(dates)), 
+                    ℵ=Metadata(:sea_ice_concentration; dataset, dates=first(dates)))
+
+#####
+##### Radiation of the Ocean and the Sea ice
+#####
+
+# to mimick snow, as in Semtner (1975) we use a sea ice albedo that is 0.64 for top temperatures 
+# above -0.1 ᵒC and 0.75 for top temperatures below -0.1 ᵒC. We also disable radiation for the ocean.
+radiation = Radiation(sea_ice_albedo=0.7, ocean_albedo=1, ocean_emissivity=0)
+
+#####
+##### Arctic coupled model
+#####
+
+arctic = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation)
+arctic = Simulation(arctic, Δt=30minutes, stop_time=730days)
+
+# Sea-ice variables
+h = sea_ice.model.ice_thickness
+ℵ = sea_ice.model.ice_concentration
+
+# Fluxes
+Tu = arctic.model.interfaces.atmosphere_sea_ice_interface.temperature
+Qˡ = arctic.model.interfaces.atmosphere_sea_ice_interface.fluxes.latent_heat
+Qˢ = arctic.model.interfaces.atmosphere_sea_ice_interface.fluxes.sensible_heat
+Qⁱ = arctic.model.interfaces.sea_ice_ocean_interface.fluxes.interface_heat
+Qᶠ = arctic.model.interfaces.sea_ice_ocean_interface.fluxes.frazil_heat
+Qᵗ = arctic.model.interfaces.net_fluxes.sea_ice_top.heat
+Qᴮ = arctic.model.interfaces.net_fluxes.sea_ice_bottom.heat
+
+# Output writers
+arctic.output_writers[:vars] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ),
+                                                 filename = "sea_ice_quantities.jld2",
+                                                 schedule = IterationInterval(12),
+                                                 overwrite_existing=true)
+
+wall_time = Ref(time_ns())
+
+using Statistics
+
+function progress(sim)
+    sea_ice = sim.model.sea_ice
+    h = first(sea_ice.model.ice_thickness)
+    ℵ = first(sea_ice.model.ice_concentration)
+    T = first(sim.model.interfaces.atmosphere_sea_ice_interface.temperature)
+
+    step_time = 1e-9 * (time_ns() - wall_time[])
+
+    msg1 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt))
+    msg2 = @sprintf("h: %.2e m, ℵ: %.2e , T: %.2e ᴼC", h, ℵ, T)
+    msg3 = @sprintf("wall time: %s \n", prettytime(step_time))
+
+    @info msg1 * msg2 * msg3
+
+    wall_time[] = time_ns()
+    return nothing
+end
+
+# And add it as a callback to the simulation.
+add_callback!(arctic, progress, IterationInterval(10))
+
+run!(arctic)
+
+h = FieldTimeSeries("sea_ice_quantities.jld2", "h")
\ No newline at end of file
diff --git a/src/ClimaOcean.jl b/src/ClimaOcean.jl
index 2a90d73e8..96a543673 100644
--- a/src/ClimaOcean.jl
+++ b/src/ClimaOcean.jl
@@ -41,6 +41,16 @@ using DataDeps
 using Oceananigans.OutputReaders: GPUAdaptedFieldTimeSeries, FieldTimeSeries
 using Oceananigans.Grids: node
 
+# Does not really matter if there is no model
+reference_density(::Nothing) = 0
+heat_capacity(::Nothing) = 0
+
+reference_density(unsupported) =
+    throw(ArgumentError("Cannot extract reference density from $(typeof(unsupported))"))
+
+heat_capacity(unsupported) =
+    throw(ArgumentError("Cannot deduce the heat capacity from $(typeof(unsupported))"))
+
 const SomeKindOfFieldTimeSeries = Union{FieldTimeSeries,
                                         GPUAdaptedFieldTimeSeries}
 
diff --git a/src/DataWrangling/ECCO/ECCO_restoring.jl b/src/DataWrangling/ECCO/ECCO_restoring.jl
index 23483d032..f76d195d1 100644
--- a/src/DataWrangling/ECCO/ECCO_restoring.jl
+++ b/src/DataWrangling/ECCO/ECCO_restoring.jl
@@ -138,15 +138,14 @@ function ECCOFieldTimeSeries(metadata::ECCOMetadata, grid::AbstractGrid;
     return fts	
 end
 
-function ECCOFieldTimeSeries(variable_name::Symbol; 
+function ECCOFieldTimeSeries(variable_name::Symbol, arch_or_grid=CPU();
                              dataset = ECCO4Monthly(),
-                             architecture = CPU(),
                              dates = all_dates(dataset, variable_name),
                              dir = download_ECCO_cache,
                              kw...)
 
     metadata = Metadata(variable_name, dates, dataset, dir)
-    return ECCOFieldTimeSeries(metadata, architecture; kw...)
+    return ECCOFieldTimeSeries(metadata, arch_or_grid; kw...)
 end
 
 # Variable names for restorable data
diff --git a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl
index 20d3c5d32..7391dfd70 100644
--- a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl
+++ b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl
@@ -174,6 +174,7 @@ function compute_net_sea_ice_fluxes!(coupled_model)
     kernel_parameters = interface_kernel_parameters(grid)
 
     sea_ice_surface_temperature = coupled_model.interfaces.atmosphere_sea_ice_interface.temperature
+    ice_concentration = sea_ice_concentration(sea_ice)
 
     launch!(arch, grid, kernel_parameters, 
             _assemble_net_sea_ice_fluxes!,
@@ -186,6 +187,7 @@ function compute_net_sea_ice_fluxes!(coupled_model)
             freshwater_flux,
             sea_ice_surface_temperature,
             downwelling_radiation,
+            ice_concentration,
             sea_ice_properties,
             atmos_sea_ice_properties)
 
@@ -201,6 +203,7 @@ end
                                                freshwater_flux, # Where do we add this one?
                                                surface_temperature,
                                                downwelling_radiation,
+                                               ice_concentration,
                                                sea_ice_properties,
                                                atmos_sea_ice_properties)
 
@@ -211,13 +214,13 @@ end
     @inbounds begin
         Ts = surface_temperature[i, j, kᴺ]
         Ts = convert_to_kelvin(sea_ice_properties.temperature_units, Ts)
-
+        ℵi = ice_concentration[i, j, 1]
         Qs = downwelling_radiation.Qs[i, j, 1]
         Qℓ = downwelling_radiation.Qℓ[i, j, 1]
         Qc = atmosphere_sea_ice_fluxes.sensible_heat[i, j, 1] # sensible or "conductive" heat flux
         Qv = atmosphere_sea_ice_fluxes.latent_heat[i, j, 1]   # latent heat flux
         Qf = sea_ice_ocean_fluxes.frazil_heat[i, j, 1]        # frazil heat flux
-        Qi = sea_ice_ocean_fluxes.interface_heat[i, j, 1]   # interfacial heat flux
+        Qi = sea_ice_ocean_fluxes.interface_heat[i, j, 1]     # interfacial heat flux
     end
 
     ρτx = atmosphere_sea_ice_fluxes.x_momentum  # zonal momentum flux                      
@@ -230,8 +233,8 @@ end
     Qu = upwelling_radiation(Ts, σ, ϵ)
     Qd = net_downwelling_radiation(i, j, grid, time, α, ϵ, Qs, Qℓ)
 
-    ΣQt = Qd + Qu + Qc + Qv
-    ΣQb = Qf + Qi
+    ΣQt = (Qd + Qu + Qc + Qv) * ℵi # We need to multiply these times the concentration?
+    ΣQb = Qf + Qi 
 
     # Mask fluxes over land for convenience
     inactive = inactive_node(i, j, kᴺ, grid, c, c, c)
diff --git a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl
index 97627ceab..aff140335 100644
--- a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl
+++ b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl
@@ -1,5 +1,6 @@
 using Oceananigans.Operators: intrinsic_vector
 using Oceananigans.Grids: inactive_node
+using ClimaOcean.OceanSimulations
 using ClimaOcean.OceanSeaIceModels.PrescribedAtmospheres: thermodynamics_parameters,
                                                           surface_layer_height,
                                                           boundary_layer_height
@@ -70,7 +71,6 @@ end
 
     i, j = @index(Global, NTuple)
     kᴺ   = size(grid, 3) # index of the top ocean cell
-    time = Time(clock.time)
 
     @inbounds begin
         uₐ = atmosphere_state.u[i, j, 1]
diff --git a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl
index c09e60f32..0bc5f1c72 100644
--- a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl
+++ b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl
@@ -116,7 +116,7 @@ end
                               h_bℓ = atmosphere_state.h_bℓ)
 
     downwelling_radiation = (; Qs, Qℓ)
-    local_interior_state = (u=uᵢ, v=vᵢ, T=Tᵢ, S=Sᵢ, h=hᵢ)
+    local_interior_state = (u=uᵢ, v=vᵢ, T=Tᵢ, S=Sᵢ, h=hᵢ, ℵ=ℵᵢ)
 
     # Estimate initial interface state
     u★ = convert(FT, 1e-4)
diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl
index 284a541ef..0336dd3d0 100644
--- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl
+++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl
@@ -245,7 +245,7 @@ function default_ao_specific_humidity(ocean)
 end
 
 """
-    ComponentInterfaces(atmosphere, ocean, sea_ice=nothing;
+    ComponentInterfaces(atmosphere, ocean, sea_ice;
                         radiation = Radiation(),
                         freshwater_density = 1000,
                         atmosphere_ocean_flux_formulation = SimilarityTheoryFluxes(),
@@ -262,11 +262,11 @@ end
                         sea_ice_reference_density = reference_density(sea_ice),
                         sea_ice_heat_capacity = heat_capacity(sea_ice))
 """
-function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing;
+function ComponentInterfaces(atmosphere, ocean, sea_ice;
                              radiation = Radiation(),
                              freshwater_density = 1000,
                              atmosphere_ocean_flux_formulation = SimilarityTheoryFluxes(eltype(ocean.model.grid)),
-                             atmosphere_sea_ice_flux_formulation = SimilarityTheoryFluxes(eltype(ocean.model.grid)),
+                             atmosphere_sea_ice_flux_formulation = SimilarityTheoryFluxes(eltype(ocean.model.grid), solver_maxiter=10000, stability_functions=atmosphere_sea_ice_stability_functions()),
                              atmosphere_ocean_interface_temperature = BulkTemperature(),
                              atmosphere_ocean_velocity_difference = RelativeVelocity(),
                              atmosphere_ocean_interface_specific_humidity = default_ao_specific_humidity(ocean),
diff --git a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl
index 108c9413e..8a693b54f 100644
--- a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl
+++ b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl
@@ -212,12 +212,13 @@ end
     F = st.internal_flux
     k = F.conductivity
     h = Ψᵢ.h
+    ℵ = Ψᵢ.ℵ
 
     # Bottom temperature at the melting temperature
     Tᵢ  = ClimaSeaIce.SeaIceThermodynamics.melting_temperature(ℙᵢ.liquidus, Ψᵢ.S)
     Tᵢ  = convert_to_kelvin(ℙᵢ.temperature_units, Tᵢ)
     Tₛ⁻ = Ψₛ.T
-
+    
     #=
     σ = ℙₛ.radiation.σ
     ϵ = ℙₛ.radiation.ϵ
@@ -225,7 +226,7 @@ end
     Tₛ = (Tᵢ - h / k * (Qₐ + 4α * Tₛ⁻^4)) / (1 + 4α * h * Tₛ⁻^3 / k)
     =#
 
-    T★ = Tᵢ - Qₐ * h / k
+    T★ = Tᵢ - Qₐ * h / k * ℵ # We need to multiply the flux by the concentration?
 
     # Fix a NaN
     T★ = ifelse(isnan(T★), Tₛ⁻, T★)
diff --git a/src/OceanSeaIceModels/InterfaceComputations/radiation.jl b/src/OceanSeaIceModels/InterfaceComputations/radiation.jl
index faf17481c..551b8c4ed 100644
--- a/src/OceanSeaIceModels/InterfaceComputations/radiation.jl
+++ b/src/OceanSeaIceModels/InterfaceComputations/radiation.jl
@@ -67,7 +67,7 @@ function Base.show(io::IO, r::Radiation)
     print(io, "├── stefan_boltzmann_constant: ", prettysummary(σ), '\n')
     print(io, "├── emission: ", summary(r.emission), '\n')
     print(io, "│   ├── ocean: ", prettysummary(r.emission.ocean), '\n')
-    print(io, "│   └── sea_ice: ", prettysummary(r.emission.ocean), '\n')
+    print(io, "│   └── sea_ice: ", prettysummary(r.emission.sea_ice), '\n')
     print(io, "└── reflection: ", summary(r.reflection), '\n')
     print(io, "    ├── ocean: ", prettysummary(r.reflection.ocean), '\n')
     print(io, "    └── sea_ice: ", prettysummary(r.reflection.sea_ice))
diff --git a/src/OceanSeaIceModels/OceanSeaIceModels.jl b/src/OceanSeaIceModels/OceanSeaIceModels.jl
index 8041f497d..8316a157f 100644
--- a/src/OceanSeaIceModels/OceanSeaIceModels.jl
+++ b/src/OceanSeaIceModels/OceanSeaIceModels.jl
@@ -22,14 +22,13 @@ using ClimaSeaIce: SeaIceModel
 using ClimaSeaIce.SeaIceThermodynamics: melting_temperature
 
 using ClimaOcean: stateindex
-
 using KernelAbstractions: @kernel, @index
 using KernelAbstractions.Extras.LoopInfo: @unroll
 
+import ClimaOcean: reference_density, heat_capacity
+
 function downwelling_radiation end
 function freshwater_flux end
-function reference_density end
-function heat_capacity end
 
 const default_gravitational_acceleration = 9.80665
 const default_freshwater_density = 1000
diff --git a/src/OceanSeaIceModels/ocean_sea_ice_model.jl b/src/OceanSeaIceModels/ocean_sea_ice_model.jl
index d2344dabb..81386b13f 100644
--- a/src/OceanSeaIceModels/ocean_sea_ice_model.jl
+++ b/src/OceanSeaIceModels/ocean_sea_ice_model.jl
@@ -73,25 +73,6 @@ function initialize!(model::OSIM)
     return nothing
 end
 
-reference_density(unsupported) =
-    throw(ArgumentError("Cannot extract reference density from $(typeof(unsupported))"))
-
-heat_capacity(unsupported) =
-    throw(ArgumentError("Cannot deduce the heat capacity from $(typeof(unsupported))"))
-
-reference_density(ocean::Simulation) = reference_density(ocean.model.buoyancy.formulation)
-reference_density(buoyancy_formulation::SeawaterBuoyancy) = reference_density(buoyancy_formulation.equation_of_state)
-reference_density(eos::TEOS10EquationOfState) = eos.reference_density
-reference_density(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_density
-
-heat_capacity(ocean::Simulation) = heat_capacity(ocean.model.buoyancy.formulation)
-heat_capacity(buoyancy_formulation::SeawaterBuoyancy) = heat_capacity(buoyancy_formulation.equation_of_state)
-heat_capacity(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_heat_capacity
-
-# Does not really matter if there is no model
-reference_density(::Nothing) = 0
-heat_capacity(::Nothing) = 0
-
 function heat_capacity(::TEOS10EquationOfState{FT}) where FT
     cₚ⁰ = SeawaterPolynomials.TEOS10.teos10_reference_heat_capacity
     return convert(FT, cₚ⁰)
@@ -160,19 +141,15 @@ function default_nan_checker(model::OceanSeaIceModel)
     return nan_checker
 end
 
-@kernel function _above_freezing_ocean_temperature!(T, grid, S, ℵ, liquidus)
+@kernel function _above_freezing_ocean_temperature!(T, grid, S, liquidus)
     i, j = @index(Global, NTuple)
     Nz = size(grid, 3)
 
     @inbounds begin
-        for k in 1:Nz-1
+        for k in 1:Nz
             Tm = melting_temperature(liquidus, S[i, j, k])
             T[i, j, k] = max(T[i, j, k], Tm)
         end
-
-        ℵi = ℵ[i, j, 1]
-        Tm = melting_temperature(liquidus, S[i, j, Nz])
-        T[i, j, Nz] = ifelse(ℵi > 0, Tm, T[i, j, Nz])
     end
 end
 
@@ -187,7 +164,7 @@ function above_freezing_ocean_temperature!(ocean, sea_ice::SeaIceSimulation)
 
     grid = ocean.model.grid
     arch = architecture(grid)
-    launch!(arch, grid, :xy, _above_freezing_ocean_temperature!, T, grid, S, ℵ, liquidus)
+    launch!(arch, grid, :xy, _above_freezing_ocean_temperature!, T, grid, S, liquidus)
 
     return nothing
 end
diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl
index 32423c690..49d48aa6e 100644
--- a/src/OceanSimulations/OceanSimulations.jl
+++ b/src/OceanSimulations/OceanSimulations.jl
@@ -1,6 +1,6 @@
 module OceanSimulations
 
-export ocean_simulation
+export ocean_simulation, PrescribedOceanModel
 
 using Oceananigans
 using Oceananigans.Units
@@ -21,6 +21,15 @@ using Oceananigans.BuoyancyFormulations: g_Earth
 using Oceananigans.Coriolis: Ω_Earth
 using Oceananigans.Operators
 
+import ClimaOcean: reference_density, heat_capacity
+
+reference_density(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) = reference_density(ocean.model.buoyancy.formulation)
+reference_density(buoyancy_formulation::SeawaterBuoyancy) = reference_density(buoyancy_formulation.equation_of_state)
+reference_density(eos::TEOS10EquationOfState) = eos.reference_density
+
+heat_capacity(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) = heat_capacity(ocean.model.buoyancy.formulation)
+heat_capacity(buoyancy_formulation::SeawaterBuoyancy) = heat_capacity(buoyancy_formulation.equation_of_state)
+
 struct Default{V}
     value :: V
 end
@@ -41,5 +50,6 @@ default_or_override(override, alternative_default=nothing) = override
 
 include("barotropic_potential_forcing.jl")
 include("ocean_simulation.jl")
+include("prescribed_ocean_model.jl")
 
 end # module
diff --git a/src/OceanSimulations/prescribed_ocean_model.jl b/src/OceanSimulations/prescribed_ocean_model.jl
new file mode 100644
index 000000000..e1ff82296
--- /dev/null
+++ b/src/OceanSimulations/prescribed_ocean_model.jl
@@ -0,0 +1,106 @@
+using Oceananigans.Models: AbstractModel
+using Oceananigans.Fields: ZeroField
+using Oceananigans.OutputReaders: extract_field_time_series, update_field_time_series!
+
+import Oceananigans.TimeSteppers: time_step!, update_state!, reset!, tick!
+import Oceananigans.Models: timestepper, update_model_field_time_series!
+
+import ClimaOcean: reference_density, heat_capacity
+import Oceananigans.Architectures: on_architecture
+
+#####
+##### A prescribed ocean...
+#####
+
+struct PrescribedOceanModel{G, C, U, T, F, Arch} <: AbstractModel{Nothing, Arch}
+    architecture :: Arch      
+    grid :: G        
+    clock :: Clock{C}
+    velocities :: U
+    tracers :: T
+    timeseries :: F
+end
+
+"""
+    PrescribedOceanModel(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0))
+
+Create a prescribed ocean model to be used in combination with ClimaOcean's `OceanSeaIceModel` 
+on a `grid` with a `clock`.
+
+Arguments
+=========
+- `timeseries`: A `NamedTuple` containing time series data for various fields. The named tuple can be empty 
+                (meaning that the ocean does not evolve in time) or include any combination of the 
+                following fields: `u`, `v`, `T`, `S`. All elements provided must be of type `FieldTimeSeries` 
+                and reside on the provided `grid`.
+"""
+function PrescribedOceanModel(timeseries=NamedTuple(); 
+                         grid, 
+                         clock=Clock{Float64}(time = 0)) 
+
+    # Make sure all elements of the timeseries are on the same grid
+    # If we decide to pass a timeseries
+    if !isempty(timeseries)
+        for k in keys(timeseries)
+            f = timeseries[k]
+            isa(f, FieldTimeSeries) ||
+                throw(ArgumentError("All variables in the `timeseries` argument must be `FieldTimeSeries`"))
+            f.grid == grid ||
+                throw(ArgumentError("All variables in the timeseries reside on the provided grid"))
+        end
+    end
+
+    τx = Field{Face, Center, Nothing}(grid)
+    τy = Field{Center, Face, Nothing}(grid)
+    Jᵀ = Field{Center, Center, Nothing}(grid)
+    Jˢ = Field{Center, Center, Nothing}(grid)
+
+    u = XFaceField(grid,  boundary_conditions=FieldBoundaryConditions(grid, (Face,   Center, Center), top = FluxBoundaryCondition(τx)))
+    v = YFaceField(grid,  boundary_conditions=FieldBoundaryConditions(grid, (Center, Face,   Center), top = FluxBoundaryCondition(τy)))
+    T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ)))
+    S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ)))
+
+    return PrescribedOceanModel(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries)
+end
+
+#####
+##### Need to extend a couple of methods
+#####
+
+function time_step!(model::PrescribedOceanModel, Δt; callbacks=[], euler=true)
+    tick!(model.clock, Δt)
+    time = Time(model.clock.time)
+
+    possible_fts = merge(model.velocities, model.tracers)
+    time_series_tuple = extract_field_time_series(possible_fts)
+
+    for fts in time_series_tuple
+        update_field_time_series!(fts, time)
+    end
+
+    # Time stepping the model!
+
+    if haskey(model.timeseries, :u)  
+        arent(model.velocities.u) .= parent(model.timeseries.u[time])
+    end
+
+    if haskey(model.timeseries, :v)  
+        parent(model.velocities.v) .= parent(model.timeseries.v[time])
+    end
+     
+    if haskey(model.timeseries, :T)  
+        parent(model.tracers.T) .= parent(model.timeseries.T[time])
+    end
+
+    if haskey(model.timeseries, :S)  
+        parent(model.tracers.S) .= parent(model.timeseries.S[time])
+    end
+
+    return nothing
+end
+
+update_state!(::PrescribedOceanModel) = nothing
+timestepper(::PrescribedOceanModel) = nothing
+
+reference_density(ocean::Simulation{<:PrescribedOceanModel}) = 1025.6
+heat_capacity(ocean::Simulation{<:PrescribedOceanModel}) = 3995.6
diff --git a/src/SeaIceSimulations.jl b/src/SeaIceSimulations.jl
index 2b7dc17b8..59576aa0c 100644
--- a/src/SeaIceSimulations.jl
+++ b/src/SeaIceSimulations.jl
@@ -20,6 +20,11 @@ using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium
 
 using ClimaOcean.OceanSimulations: Default
 
+import ClimaOcean: reference_density, heat_capacity
+
+reference_density(sea_ice::Simulation{<:SeaIceModel}) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_density
+heat_capacity(sea_ice::Simulation{<:SeaIceModel}) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_heat_capacity
+
 function sea_ice_simulation(grid;
                             Δt = 5minutes,
                             ice_salinity = 0, # psu
diff --git a/test/test_ocean_sea_ice_model.jl b/test/test_ocean_sea_ice_model.jl
index 55b6ac3f5..16d55fc9d 100644
--- a/test/test_ocean_sea_ice_model.jl
+++ b/test/test_ocean_sea_ice_model.jl
@@ -100,5 +100,24 @@ using ClimaSeaIce.Rheologies
             true
         end
     end
+
+    @testset "Test a prescribed ocean model" begin 
+        grid = LatitudeLongitudeGrid(size = (10, 10, 10), latitude = (20, 30), longitude = (40, 50), z = (-100, 0))
+        T = ECCOFieldTimeSeries(:temperature, grid; time_indices_in_memory = 2)
+        S = ECCOFieldTimeSeries(:salinity, grid; time_indices_in_memory = 2)
+        ocean_model = PrescribedOcean((; T, S); grid)
+
+        time_step!(ocean_model, T.times[2])
+
+        @test ocean_model.clock.time == T.times[2]
+        @test ocean_model.tracers.T.data == T[2].data
+        @test ocean_model.tracers.S.data == S[2].data
+
+        time_step!(ocean_model, 10days)
+
+        @test ocean_model.clock.time == 10days
+        @test T.backend.start == 2
+        @test ocean_model.tracers.T.data != T[3].data 
+    end
 end