diff --git a/experiments/AMIP/modular/components/ocean/eisenman_seaice.jl b/experiments/AMIP/modular/components/ocean/eisenman_seaice.jl
index fe16cbe88a..d1d7e54994 100644
--- a/experiments/AMIP/modular/components/ocean/eisenman_seaice.jl
+++ b/experiments/AMIP/modular/components/ocean/eisenman_seaice.jl
@@ -140,6 +140,11 @@ function solve_eisenman_model!(Y, Ya, p, thermo_params, Δt)
     return Y, Ya
 end
 
+"""
+    ∑tendencies(dY, Y, cache, _)
+
+Calculate the tendencies for the Eisenman-Zhang sea ice model.
+"""
 function ∑tendencies(dY, Y, cache, _)
     FT = eltype(dY)
     Δt = cache.Δt
@@ -189,10 +194,6 @@ function eisenman_seaice_init(
     # initiate prognostic variables
     Y, Ya = state_init(params_ice, space)
 
-
-    # problem = OrdinaryDiffEq.ODEProblem(∑tendencies, Y, tspan, (; Ya = Ya, Δt = dt, params_ice = params_ice, area_fraction = area_fraction))
-    # integrator = OrdinaryDiffEq.init(problem, stepper, dt = dt, saveat = saveat)
-
     ode_algo = CTS.ExplicitAlgorithm(stepper)
     ode_function = CTS.ClimaODEFunction(T_exp! = ∑tendencies, dss! = weighted_dss_slab!)
 
diff --git a/experiments/AMIP/modular/user_io/viz_explorer.jl b/experiments/AMIP/modular/user_io/viz_explorer.jl
index cfe77dcce2..765d56e052 100644
--- a/experiments/AMIP/modular/user_io/viz_explorer.jl
+++ b/experiments/AMIP/modular/user_io/viz_explorer.jl
@@ -33,7 +33,7 @@ function plot_anim(cs, out_dir = ".")
     combined_field = zeros(boundary_space)
     sol_slab = slab_land_sim.integrator.sol
 
-    if mode_name[1:10] == "slabplanet"
+    if mode_name == "slabplanet"
         sol_slab_ocean = slab_ocean_sim.integrator.sol
         anim = Plots.@animate for (bucketu, oceanu) in zip(sol_slab.u, sol_slab_ocean.u)
             land_T_sfc = get_land_temp_from_state(cs.model_sims.land_sim, bucketu)
@@ -44,6 +44,17 @@ function plot_anim(cs, out_dir = ".")
             )
             Plots.plot(combined_field)
         end
+    elseif mode_name == "slabplanet_eisenman"
+        slab_ice_sim = slab_ice_sim.integrator.sol
+        anim = Plots.@animate for (bucketu, iceu) in zip(sol_slab.u, slab_ice_sim.u)
+            land_T_sfc = get_land_temp_from_state(cs.model_sims.land_sim, bucketu)
+            combine_surfaces_from_sol!(
+                combined_field,
+                cs.surface_fractions,
+                (; land = land_T_sfc, ocean = FT(0), ice = iceu.T_sfc),
+            )
+            Plots.plot(combined_field)
+        end
 
     elseif mode_name == "amip"
         sol_slab_ice = slab_ice_sim.integrator.sol
diff --git a/src/FluxCalculator.jl b/src/FluxCalculator.jl
index 59ce31eeef..bde9ce62ad 100644
--- a/src/FluxCalculator.jl
+++ b/src/FluxCalculator.jl
@@ -309,9 +309,9 @@ function surface_thermo_state(
     colidx::Fields.ColumnIndex,
     ΔT_sfc = 0,
 )
-    FT = eltype(thermo_state_int)
+    FT = eltype(thermo_state_int[1])
     @warn("Simulation " * Interfacer.name(sim) * " uses the default thermo (saturated) surface state", maxlog = 10)
-    T_sfc = Interfacer.get_field(sim, Val(:surface_temperature), colidx) .+ FT(ΔT_sfc) #
+    T_sfc = Interfacer.get_field(sim, Val(:surface_temperature), colidx) .+ FT(ΔT_sfc)
     ρ_sfc = extrapolate_ρ_to_sfc.(thermo_params, thermo_state_int, T_sfc) # ideally the # calculate elsewhere, here just getter...
     q_sfc = TD.q_vap_saturation_generic.(thermo_params, T_sfc, ρ_sfc, TD.Liquid()) # default: saturated liquid surface
     @. TD.PhaseEquil_ρTq.(thermo_params, ρ_sfc, T_sfc, q_sfc)
@@ -352,7 +352,13 @@ function get_surface_fluxes_point!(inputs, surface_params::SF.Parameters.Surface
     # moisture
     F_turb_moisture = @. SF.evaporation(surface_params, inputs, outputs.Ch)
 
-    return F_turb_ρτxz, F_turb_ρτyz, F_shf, F_lhf, F_turb_moisture
+    return (;
+        F_turb_ρτxz = F_turb_ρτxz,
+        F_turb_ρτyz = F_turb_ρτyz,
+        F_shf = F_shf,
+        F_lhf = F_lhf,
+        F_turb_moisture = F_turb_moisture,
+    )
 end
 
 """
diff --git a/test/flux_calculator_tests.jl b/test/flux_calculator_tests.jl
index 2f93f9f9b7..460f4c0534 100644
--- a/test/flux_calculator_tests.jl
+++ b/test/flux_calculator_tests.jl
@@ -307,7 +307,7 @@ end
 @testset "surface_thermo_state" begin
     boundary_space = TestHelper.create_space(FT)
     _ones = Fields.ones(boundary_space)
-    sim = DummySurfaceSimulation3(
+    surface_sim = DummySurfaceSimulation3(
         [],
         [],
         [],
@@ -317,5 +317,6 @@ end
     thermo_params = get_thermo_params(atmos_sim)
     thermo_state_int = Interfacer.get_field(atmos_sim, Val(:thermo_state_int))
     colidx = Fields.ColumnIndex{2}((1, 1), 73) # arbitrary index
-    @test surface_thermo_state(sim, thermo_params, thermo_state_int[colidx], colidx).ρ == thermo_state_int[colidx].ρ
+    @test surface_thermo_state(surface_sim, thermo_params, thermo_state_int[colidx], colidx).ρ ==
+          thermo_state_int[colidx].ρ
 end