diff --git a/Project.toml b/Project.toml
index bdb05b22e73..90b19a4131d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,7 @@ ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 ArtifactWrappers = "a14bc488-3040-4b00-9dc1-f6467924858a"
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 AtmosphericProfilesLibrary = "86bc3604-9858-485a-bdbe-831ec50de11d"
+BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 ClimaComms = "3a4d1b5c-c61d-41fd-a00a-5873ba7a1b0d"
 ClimaCore = "d414da3d-4745-48bb-8d80-42e94e092884"
 ClimaDiagnostics = "1ecacbb8-0713-4841-9a07-eb5aa8a2d53f"
@@ -19,6 +20,7 @@ CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Insolation = "e98cc03f-d57e-4e3c-b70c-8d51efe9e0d8"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
@@ -41,6 +43,7 @@ ArgParse = "1"
 ArtifactWrappers = "0.2"
 Artifacts = "1"
 AtmosphericProfilesLibrary = "0.1.7"
+BlockArrays = "1"
 ClimaComms = "0.6.4"
 ClimaCore = "0.14.12"
 ClimaDiagnostics = "0.2.4"
@@ -51,6 +54,7 @@ CloudMicrophysics = "0.22.3"
 Dates = "1"
 DiffEqBase = "6.145"
 FastGaussQuadrature = "0.5, 1"
+ForwardDiff = "0.10.36"
 Insolation = "0.9.2"
 Interpolations = "0.15.1"
 LazyArtifacts = "1"
diff --git a/config/default_configs/default_config.yml b/config/default_configs/default_config.yml
index 56abd972b97..a2531eb718a 100644
--- a/config/default_configs/default_config.yml
+++ b/config/default_configs/default_config.yml
@@ -86,6 +86,15 @@ eisenstat_walker_forcing_alpha:
 jvp_step_adjustment:
   help: "Amount by which the step size of the forward difference approximation of the Jacobian-vector product in the Krylov method should be scaled (only used if `use_krylov_method` is `true`)"
   value: 1.0
+use_exact_jacobian:
+  help: "Whether to solve the linear system using the exact Jacobian, which is computed with forward-mode automatic differentiation (default is `false`)"
+  value: true # false
+n_steps_update_exact_jacobian:
+  help: "Number of timesteps between updates to the exact Jacobian (default is 0, which corresponds to updating the exact Jacobian on each Newton iteration)"
+  value: 10 # 0
+debug_approximate_jacobian:
+  help: "Whether to check the accuracy of the approximate Jacobian against the exact Jacobian every `n_steps_update_exact_jacobian` timesteps (default is `false`)"
+  value: false
 # Radiation
 rad:
   help: "Radiation model [`nothing` (default), `gray`, `clearsky`, `allsky`, `allskywithclear`]"
diff --git a/config/model_configs/prognostic_edmfx_bomex_implicit_column.yml b/config/model_configs/prognostic_edmfx_bomex_implicit_column.yml
index a7e29317ff6..faea36866cf 100644
--- a/config/model_configs/prognostic_edmfx_bomex_implicit_column.yml
+++ b/config/model_configs/prognostic_edmfx_bomex_implicit_column.yml
@@ -7,7 +7,7 @@ turbconv: "prognostic_edmfx"
 implicit_diffusion: true
 implicit_sgs_advection: true
 approximate_linear_solve_iters: 2
-max_newton_iters_ode: 3
+max_newton_iters_ode: 10
 edmfx_upwinding: first_order
 edmfx_entr_model: "Generalized"
 edmfx_detr_model: "Generalized"
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index 498ac7b3244..35067d7da73 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -295,7 +295,7 @@ weakdeps = ["SparseArrays"]
     ChainRulesCoreSparseArraysExt = "SparseArrays"
 
 [[deps.ClimaAtmos]]
-deps = ["Adapt", "ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaParams", "ClimaTimeSteppers", "ClimaUtilities", "CloudMicrophysics", "Dates", "DiffEqBase", "FastGaussQuadrature", "Insolation", "Interpolations", "LazyArtifacts", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "RRTMGP", "SciMLBase", "StaticArrays", "Statistics", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities", "YAML"]
+deps = ["Adapt", "ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "BlockArrays", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaParams", "ClimaTimeSteppers", "ClimaUtilities", "CloudMicrophysics", "Dates", "DiffEqBase", "FastGaussQuadrature", "ForwardDiff", "Insolation", "Interpolations", "LazyArtifacts", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "RRTMGP", "SciMLBase", "StaticArrays", "Statistics", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities", "YAML"]
 path = ".."
 uuid = "b2c96348-7fb7-4fe0-8da9-78d88439e717"
 version = "0.27.5"
diff --git a/examples/Manifest.toml b/examples/Manifest.toml
index d247f216bf8..710f48519aa 100644
--- a/examples/Manifest.toml
+++ b/examples/Manifest.toml
@@ -319,7 +319,7 @@ version = "0.5.7"
     Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 
 [[deps.ClimaAtmos]]
-deps = ["Adapt", "ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaParams", "ClimaTimeSteppers", "ClimaUtilities", "CloudMicrophysics", "Dates", "DiffEqBase", "FastGaussQuadrature", "Insolation", "Interpolations", "LazyArtifacts", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "RRTMGP", "SciMLBase", "StaticArrays", "Statistics", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities", "YAML"]
+deps = ["Adapt", "ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "BlockArrays", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaParams", "ClimaTimeSteppers", "ClimaUtilities", "CloudMicrophysics", "Dates", "DiffEqBase", "FastGaussQuadrature", "ForwardDiff", "Insolation", "Interpolations", "LazyArtifacts", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "RRTMGP", "SciMLBase", "StaticArrays", "Statistics", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities", "YAML"]
 path = ".."
 uuid = "b2c96348-7fb7-4fe0-8da9-78d88439e717"
 version = "0.27.5"
diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl
index 8777ee59eaf..bb2971a2c56 100644
--- a/examples/hybrid/driver.jl
+++ b/examples/hybrid/driver.jl
@@ -145,10 +145,17 @@ if ClimaComms.iamroot(config.comms_ctx)
     include(
         joinpath(pkgdir(CA), "regression_tests", "self_reference_or_path.jl"),
     )
+
+    # TODO: Improve design of Diagnostics to make Jacobian plotting less clunky.
+    (; dict_writer) = simulation.output_writers
+    Yₜ_end = similar(integrator.u)
+    CA.implicit_tendency!(Yₜ_end, integrator.u, integrator.p, integrator.t)
+
     @info "Plotting"
     path = self_reference_or_path() # __build__ path (not job path)
     if path == :self_reference
-        make_plots(Val(Symbol(reference_job_id)), simulation.output_dir)
+        paths = [simulation.output_dir]
+        make_plots(Val(Symbol(reference_job_id)), paths, dict_writer, Yₜ_end)
     else
         main_job_path = joinpath(path, reference_job_id)
         nc_dir = joinpath(main_job_path, "nc_files")
@@ -170,11 +177,11 @@ if ClimaComms.iamroot(config.comms_ctx)
         end
 
         paths = if isempty(readdir(nc_dir))
-            simulation.output_dir
+            [simulation.output_dir]
         else
             [simulation.output_dir, nc_dir]
         end
-        make_plots(Val(Symbol(reference_job_id)), paths)
+        make_plots(Val(Symbol(reference_job_id)), paths, dict_writer, Yₜ_end)
     end
     @info "Plotting done"
 
diff --git a/perf/Manifest.toml b/perf/Manifest.toml
index 052f4d32a10..f72b3260929 100644
--- a/perf/Manifest.toml
+++ b/perf/Manifest.toml
@@ -330,7 +330,7 @@ version = "0.5.6"
     GeoMakie = "db073c08-6b98-4ee5-b6a4-5efafb3259c6"
 
 [[deps.ClimaAtmos]]
-deps = ["Adapt", "ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaParams", "ClimaTimeSteppers", "ClimaUtilities", "CloudMicrophysics", "Dates", "DiffEqBase", "FastGaussQuadrature", "Insolation", "Interpolations", "LazyArtifacts", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "RRTMGP", "SciMLBase", "StaticArrays", "Statistics", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities", "YAML"]
+deps = ["Adapt", "ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "BlockArrays", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaParams", "ClimaTimeSteppers", "ClimaUtilities", "CloudMicrophysics", "Dates", "DiffEqBase", "FastGaussQuadrature", "ForwardDiff", "Insolation", "Interpolations", "LazyArtifacts", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "RRTMGP", "SciMLBase", "StaticArrays", "Statistics", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities", "YAML"]
 path = ".."
 uuid = "b2c96348-7fb7-4fe0-8da9-78d88439e717"
 version = "0.27.5"
diff --git a/post_processing/ci_plots.jl b/post_processing/ci_plots.jl
index b058c2d71a6..4256d32e8af 100644
--- a/post_processing/ci_plots.jl
+++ b/post_processing/ci_plots.jl
@@ -49,6 +49,7 @@ import ClimaCoreSpectra: power_spectrum_2d
 
 using Poppler_jll: pdfunite
 import Base.Filesystem
+import LinearAlgebra: diag
 
 const days = 86400
 
@@ -71,19 +72,27 @@ function make_plots(sim, _)
 end
 
 """
-    make_plots(sim, simulation_path::AbstractString)
-    make_plots(sim, simulation_paths::Iterable{AbstractString})
+    make_plots(sim, simulation_paths, dict_writer, Yₜ_end)
 
-Plot the corresponding sets for the given `sim`.
+Plot the available data for the given `sim`.
+    
+When `simulation_paths` contains one string, the data at that path is plotted.
 
-When `simulation_path` is a string, use the data in the path.
+When `simulations_paths` contains multiple strings, the data at those paths is
+compared side-by-side.
 
-When `simulation_paths` is a collection of strings, use all those paths and making
-side-by-side comparisons.
+All of the data stored in the `dict_writer` is also plotted by default.
 """
-function make_plots(sim, simulation_path::AbstractString)
-    paths = [simulation_path]
-    make_plots(sim, paths)
+function make_plots(sim, simulation_paths, dict_writer, Yₜ_end)
+    for (time_series_name, time_series) in dict_writer.dict
+        filename = if startswith(time_series_name, "ejac")
+            "exact_jacobian"
+        else
+            error("Unrecognized time series name $time_series_name")
+        end
+        plot_jacobian(simulation_paths[1], time_series, Yₜ_end, filename)
+    end
+    make_plots(sim, simulation_paths)
 end
 
 # The contour plot functions in ClimaAnalysis work by finding the nearest slice available.
@@ -469,6 +478,372 @@ function plot_spectrum_with_line!(grid_loc, spectrum; exponent = -3.0)
     return nothing
 end
 
+function plot_jacobian(output_path, jacobian_time_series, Yₜ_end, filename)
+    # Take the transpose of each Jacobian matrix so that its rows are plotted
+    # along the y-axis and its columns are plotted along the x-axis.
+    time_series =
+        Dict(map(((t, matrix),) -> t => matrix', [jacobian_time_series...]))
+
+    n_max_plot_times = 6
+    n_available_times = length(time_series)
+    available_times = sort(collect(keys(time_series)))
+    plot_times = if n_available_times < n_max_plot_times
+        available_times
+    else
+        n_times_to_skip = fld(n_available_times, n_max_plot_times)
+        available_times[1:n_times_to_skip:(n_max_plot_times * n_times_to_skip)]
+    end
+
+    field_names = collect(CA.scalar_field_names(Yₜ_end))
+    first_tick_positions = map(field_names) do field_name
+        tick_position = 1
+        for prev_field_name in field_names
+            field_name == prev_field_name && return tick_position
+            prev_field = CA.MatrixFields.get_field(Yₜ_end, prev_field_name)
+            if prev_field isa CA.Fields.SpectralElementField
+                tick_position += 1
+            else
+                tick_position += CA.Spaces.nlevels(axes(prev_field))
+            end
+        end
+    end
+    last_field = CA.MatrixFields.get_field(Yₜ_end, field_names[end])
+    last_index = if last_field isa CA.Fields.SpectralElementField
+        first_tick_positions[end]
+    else
+        first_tick_positions[end] + CA.Spaces.nlevels(axes(last_field)) - 1
+    end
+    last_tick_positions = [(first_tick_positions[2:end] .- 1)..., last_index]
+
+    center_tick_positions = (first_tick_positions .+ last_tick_positions) ./ 2
+    boundary_tick_positions = [0.5, (last_tick_positions .+ 0.5)...]
+    index_ranges = map(UnitRange, first_tick_positions, last_tick_positions)
+
+    limit_padding = length(index_ranges[1]) / 20
+    limit_positions =
+        extrema(boundary_tick_positions) .+ (-limit_padding, limit_padding)
+
+    tick_labels = map(field_names) do field_name
+        replace(
+            string(field_name),
+            "@name(" => "",
+            ".components.data" => "",
+            ":(" => "",
+            ')' => "",
+        )
+    end
+    axis_kwargs = (;
+        limits = (limit_positions, limit_positions),
+        titlegap = 100,
+        yreversed = true,
+        xlabel = "∂Y index",
+        ylabel = "∂Yₜ index",
+        xlabelpadding = 60,
+        ylabelpadding = 60,
+        xticks = (center_tick_positions, tick_labels),
+        yticks = (center_tick_positions, tick_labels),
+        xticksvisible = false,
+        yticksvisible = false,
+        xticklabelrotation = pi / 4,
+        xticklabelpad = 50,
+        yticklabelpad = 80,
+        xminorticks = boundary_tick_positions,
+        yminorticks = boundary_tick_positions,
+        xminorticksvisible = true,
+        yminorticksvisible = true,
+        xminorticksize = 50,
+        yminorticksize = 50,
+        xminortickwidth = 10,
+        yminortickwidth = 10,
+        xgridvisible = false,
+        ygridvisible = false,
+        xminorgridvisible = true,
+        yminorgridvisible = true,
+        xminorgridwidth = 10,
+        yminorgridwidth = 10,
+        spinewidth = 10,
+    ) # Flip the y-axis so that the diagonal runs from top-left to bottom-right.
+    colorbar_kwargs = (;
+        size = 150,
+        labelpadding = 60,
+        ticklabelpad = 50,
+        ticksize = 50,
+        tickwidth = 10,
+        spinewidth = 10,
+    )
+    figure_kwargs = (; figure_padding = 300, fontsize = 150)
+    plot_size = 5000
+    figure_size = (2.2, cld(length(plot_times), 2)) .* plot_size
+
+    typical_tendency_values = Dict(
+        CA.MatrixFields.@name(c.ρ) => 1e-8,
+        CA.MatrixFields.@name(c.uₕ.components.data.:1) => 1e-6,
+        CA.MatrixFields.@name(c.uₕ.components.data.:2) => 1e-6,
+        CA.MatrixFields.@name(c.ρe_tot) => 1e-1,
+        CA.MatrixFields.@name(c.ρq_tot) => 1e-7,
+        CA.MatrixFields.@name(c.ρq_liq) => 1e-8,
+        CA.MatrixFields.@name(c.ρq_ice) => 1e-8,
+        CA.MatrixFields.@name(c.ρq_rai) => 1e-8,
+        CA.MatrixFields.@name(c.ρq_sno) => 1e-8,
+        CA.MatrixFields.@name(c.sgs⁰.ρatke) => 1e-3,
+        CA.MatrixFields.@name(c.sgsʲs.:(1).ρa) => 1e-5,
+        CA.MatrixFields.@name(c.sgsʲs.:(1).mse) => 1e-1,
+        CA.MatrixFields.@name(c.sgsʲs.:(1).q_tot) => 1e-7,
+        CA.MatrixFields.@name(f.u₃.components.data.:1) => 1e-2,
+        CA.MatrixFields.@name(f.sgsʲs.:(1).u₃.components.data.:1) => 1e-2,
+        CA.MatrixFields.@name(sfc.T) => 1e-3,
+        CA.MatrixFields.@name(sfc.water) => 1e-7,
+    )
+    rescaling_vector = ClimaComms.allowscalar(
+        Vector,
+        ClimaComms.device(Yₜ_end),
+        CA.Fields.column(Yₜ_end, 1, 1, 1),
+    )
+    for (field_name, index_range) in zip(field_names, index_ranges)
+        block_max = CA.smooth_maximum(rescaling_vector[index_range])
+        rescaling_vector[index_range] .=
+            block_max < 1e-10 ? typical_tendency_values[field_name] : block_max
+    end
+    rescaling_matrix = (inv.(rescaling_vector) * rescaling_vector')'
+    # Take the transpose of the rescaling matrix so that its rows and columns
+    # are consistent with the Jacobian matrices.
+
+    sign_or_nan(x) = iszero(x) ? typeof(x)(NaN) : sign(x)
+    logabs_or_nan(x) = iszero(x) ? typeof(x)(NaN) : log10(abs(x))
+    entry_sign(matrix) = sign_or_nan.(matrix)
+    entry_logabs(matrix) = logabs_or_nan.(matrix)
+    entry_logabs_rescaled(matrix) = entry_logabs(rescaling_matrix .* matrix)
+    function block_logabs_rescaled(matrix)
+        block_max_matrix = similar(matrix)
+        for row_index_range in index_ranges, col_index_range in index_ranges
+            block_max_matrix[row_index_range, col_index_range] .=
+                CA.smooth_maximum(matrix[row_index_range, col_index_range])
+        end
+        return entry_logabs_rescaled(block_max_matrix)
+    end
+    function row_logabs_rescaled(matrix)
+        row_max_matrix = similar(matrix)
+        for row_index_range in index_ranges, col in axes(matrix, 1)
+            row_max_matrix[row_index_range, col] .=
+                CA.smooth_maximum(matrix[row_index_range, col])
+        end
+        return entry_logabs_rescaled(row_max_matrix)
+    end
+    function block_bandwidth(matrix)
+        block_bandwidth_matrix = similar(matrix, Int)
+        for row_index_range in index_ranges, col_index_range in index_ranges
+            block = matrix[row_index_range, col_index_range]
+            bands = (1 - length(row_index_range)):(length(col_index_range) - 1)
+            block_bandwidth_matrix[row_index_range, col_index_range] .=
+                all(iszero, block) ? 0 :
+                findlast(band -> any(!iszero, diag(block, band)), bands) -
+                findfirst(band -> any(!iszero, diag(block, band)), bands) + 1
+        end
+        return block_bandwidth_matrix
+    end
+
+    bandwidth_matrices = map(t -> block_bandwidth(time_series[t]), plot_times)
+    sign_matrices = map(t -> entry_sign(time_series[t]), plot_times)
+    max_bandwidth = maximum(maximum, bandwidth_matrices)
+
+    categorical(colormap, n) = CairoMakie.cgrad(colormap, n; categorical = true)
+    main_colormap = categorical(:tol_iridescent, 21)
+    bandwidth_colors = [CairoMakie.RGB(1, 1, 1), main_colormap[1:(end - 1)]...]
+    bandwidth_colormap = categorical(bandwidth_colors, max_bandwidth + 1)
+    sign_colormap = categorical(:RdBu_5, 2)
+    rescaling_colors =
+        setindex!(CairoMakie.to_colormap(:RdBu_11), CairoMakie.RGB(1, 1, 1), 6)
+    rescaling_colormap = categorical(rescaling_colors, 21)
+
+    first_column_field = CA.Fields.column(Yₜ_end.c, 1, 1, 1)
+    coord_field =
+        CA.Fields.coordinate_field(CA.Fields.level(first_column_field, 1))
+    coord =
+        ClimaComms.allowscalar(getindex, ClimaComms.device(Yₜ_end), coord_field)
+    column_loc_str = if coord isa CA.Geometry.XZPoint
+        "x = $(coord.x) m"
+    elseif coord isa CA.Geometry.XYZPoint
+        "x = $(coord.x) m, y = $(coord.y) m"
+    elseif coord isa CA.Geometry.LatLongZPoint
+        "lat = $(coord.lat)°, long = $(coord.long)°"
+    else
+        error("Unrecognized coordinate type $(typeof(coord))")
+    end
+
+    page_file_paths = map(n -> joinpath(output_path, "$filename$n.pdf"), 1:6)
+    full_file_path = joinpath(output_path, "$filename.pdf")
+
+    figure = CairoMakie.Figure(; size = figure_size, figure_kwargs...)
+    for (index, (t, bandwidth_matrix, sign_matrix)) in
+        enumerate(zip(plot_times, bandwidth_matrices, sign_matrices))
+        grid_position = figure[cld(index, 2), (index - 1) % 2 + 1]
+        grid_layout = CairoMakie.GridLayout(grid_position)
+        title = "∂Yₜ/∂Y at $column_loc_str, t = $(CA.time_and_units_str(t))"
+        axis = CairoMakie.Axis(grid_layout[1, 1]; title, axis_kwargs...)
+        bandwidth_plot = CairoMakie.heatmap!(
+            axis,
+            boundary_tick_positions,
+            boundary_tick_positions,
+            bandwidth_matrix[first_tick_positions, first_tick_positions];
+            colormap = bandwidth_colormap,
+            colorrange = (-0.5, max_bandwidth + 0.5),
+        )
+        CairoMakie.translate!(bandwidth_plot, 0, 0, -100) # Put plot under grid.
+        sign_plot = CairoMakie.heatmap!(
+            axis,
+            sign_matrix;
+            colormap = sign_colormap,
+            colorrange = (-1, 1),
+        )
+        grid_sublayout = CairoMakie.GridLayout(grid_layout[1, 2])
+        CairoMakie.Colorbar(
+            grid_sublayout[1, 1],
+            bandwidth_plot;
+            label = "Bandwidth of matrix block",
+            colorbar_kwargs...,
+        )
+        CairoMakie.Colorbar(
+            grid_sublayout[2, 1],
+            sign_plot;
+            label = "Sign of matrix entry",
+            ticks = ([-0.5, 0.5], Makie.latexstring.(["-", "+"])),
+            colorbar_kwargs...,
+            labelpadding = 30,
+        )
+        CairoMakie.colgap!(grid_layout, CairoMakie.Relative(0.05))
+        CairoMakie.rowgap!(grid_sublayout, CairoMakie.Relative(0.05))
+    end
+    CairoMakie.colsize!(figure.layout, 1, CairoMakie.Aspect(1, 1))
+    CairoMakie.colsize!(figure.layout, 2, CairoMakie.Aspect(1, 1))
+    CairoMakie.colgap!(figure.layout, CairoMakie.Relative(0.08))
+    CairoMakie.rowgap!(figure.layout, CairoMakie.Relative(0.05))
+    CairoMakie.save(page_file_paths[1], figure)
+
+    tickformat = Base.Fix1(map, x -> Makie.latexstring("10^{$(round(Int, x))}"))
+    for (index, (logabs_transform, value_type_str, is_rescaled)) in enumerate((
+        (block_logabs_rescaled, "block", true),
+        (row_logabs_rescaled, "block row", true),
+        (entry_logabs_rescaled, "entry", true),
+        (entry_logabs, "entry", false),
+    ))
+        matrices = map(t -> logabs_transform(time_series[t]), plot_times)
+        max_logabs =
+            maximum(matrix -> maximum(filter(!isnan, matrix)), matrices)
+        figure = CairoMakie.Figure(; size = figure_size, figure_kwargs...)
+        for (index, (t, matrix)) in enumerate(zip(plot_times, matrices))
+            grid_position = figure[cld(index, 2), (index - 1) % 2 + 1]
+            grid_layout = CairoMakie.GridLayout(grid_position)
+            title = "∂Yₜ/∂Y at $column_loc_str, t = $(CA.time_and_units_str(t))"
+            axis = CairoMakie.Axis(grid_layout[1, 1]; title, axis_kwargs...)
+            plot_args = if value_type_str == "block"
+                submatrix = matrix[first_tick_positions, first_tick_positions]
+                (boundary_tick_positions, boundary_tick_positions, submatrix)
+            elseif value_type_str == "block row"
+                submatrix = matrix[first_tick_positions, :]
+                (boundary_tick_positions, 1:last_index, submatrix)
+            else
+                (matrix,)
+            end
+            plot = CairoMakie.heatmap!(
+                axis,
+                plot_args...;
+                colormap = main_colormap,
+                colorrange = (max_logabs - 6, max_logabs),
+                lowclip = main_colormap[1],
+            ) # Only color in the top 6 orders of magnitude.
+            CairoMakie.translate!(plot, 0, 0, -100) # Put plot under grid.
+            is_max = value_type_str == "entry"
+            label = "Absolute value of $(is_max ? " max of" : "") matrix \
+                     $value_type_str$(is_rescaled ? "  [s⁻¹]" : "")"
+            CairoMakie.Colorbar(
+                grid_layout[1, 2],
+                plot;
+                label,
+                tickformat,
+                colorbar_kwargs...,
+            )
+            CairoMakie.colgap!(grid_layout, CairoMakie.Relative(0.05))
+        end
+        CairoMakie.colsize!(figure.layout, 1, CairoMakie.Aspect(1, 1))
+        CairoMakie.colsize!(figure.layout, 2, CairoMakie.Aspect(1, 1))
+        CairoMakie.colgap!(figure.layout, CairoMakie.Relative(0.08))
+        CairoMakie.rowgap!(figure.layout, CairoMakie.Relative(0.05))
+        CairoMakie.save(page_file_paths[index + 1], figure)
+    end
+
+    figure = CairoMakie.Figure(; size = (2.2, 1) .* plot_size, figure_kwargs...)
+
+    grid_position = figure[1, 1]
+    title = "Rescaling vector for Yₜ"
+    ylimits = (minimum(rescaling_vector) / 2, maximum(rescaling_vector) * 2)
+    vector_axis_kwargs = (;
+        limits = (limit_positions, ylimits),
+        titlegap = 100,
+        yscale = log10,
+        xlabel = "Yₜ index",
+        ylabel = "Vector entry",
+        xlabelpadding = 60,
+        ylabelpadding = 60,
+        xticks = (center_tick_positions, tick_labels),
+        xticksvisible = false,
+        xticklabelrotation = pi / 4,
+        xticklabelpad = 50,
+        yticklabelpad = 50,
+        xminorticks = boundary_tick_positions,
+        xminorticksvisible = true,
+        xminorticksize = 50,
+        yticksize = 50,
+        xminortickwidth = 10,
+        ytickwidth = 10,
+        xgridvisible = false,
+        xminorgridvisible = true,
+        xminorgridwidth = 10,
+        ygridwidth = 10,
+        spinewidth = 10,
+    )
+    axis = CairoMakie.Axis(grid_position; title, vector_axis_kwargs...)
+    CairoMakie.lines!(
+        axis,
+        1:last_index,
+        rescaling_vector;
+        color = rescaling_colormap[end],
+        linewidth = 40,
+    )
+
+    matrix = logabs_or_nan.(rescaling_matrix)
+    grid_layout = CairoMakie.GridLayout(figure[1, 2])
+    title = "Rescaling matrix for ∂Yₜ/∂Y"
+    axis = CairoMakie.Axis(grid_layout[1, 1]; title, axis_kwargs...)
+    plot = CairoMakie.heatmap!(
+        axis,
+        boundary_tick_positions,
+        boundary_tick_positions,
+        matrix[first_tick_positions, first_tick_positions];
+        colormap = rescaling_colormap,
+    )
+    CairoMakie.translate!(plot, 0, 0, -100) # Put plot under grid.
+    CairoMakie.Colorbar(
+        grid_layout[1, 2],
+        plot;
+        label = "Matrix entry",
+        tickformat,
+        colorbar_kwargs...,
+    )
+    CairoMakie.colgap!(grid_layout, CairoMakie.Relative(0.05))
+
+    CairoMakie.colsize!(figure.layout, 1, CairoMakie.Aspect(1, 1))
+    CairoMakie.colsize!(figure.layout, 2, CairoMakie.Aspect(1, 1))
+    CairoMakie.colgap!(figure.layout, CairoMakie.Relative(0.15))
+    CairoMakie.save(page_file_paths[end], figure)
+
+    pdfunite() do unite
+        run(Cmd([unite, page_file_paths..., full_file_path]))
+    end
+    for page_file_path in page_file_paths
+        Filesystem.rm(page_file_path; force = true)
+    end
+end
+
 ColumnPlots = Union{
     Val{:single_column_hydrostatic_balance_ft64},
     Val{:single_column_radiative_equilibrium_gray},
diff --git a/regression_tests/ref_counter.jl b/regression_tests/ref_counter.jl
index b763073bad7..c01ebd1d981 100644
--- a/regression_tests/ref_counter.jl
+++ b/regression_tests/ref_counter.jl
@@ -1,4 +1,4 @@
-179
+180
 
 #=
 179:
diff --git a/src/ClimaAtmos.jl b/src/ClimaAtmos.jl
index d0c6cbda5ba..1523c2aa9dd 100644
--- a/src/ClimaAtmos.jl
+++ b/src/ClimaAtmos.jl
@@ -48,6 +48,10 @@ include(joinpath("prognostic_equations", "zero_tendency.jl"))
 
 include(joinpath("prognostic_equations", "implicit", "implicit_tendency.jl"))
 include(joinpath("prognostic_equations", "implicit", "implicit_solver.jl"))
+include(joinpath("prognostic_equations", "implicit", "approx_jacobian.jl"))
+include(joinpath("prognostic_equations", "implicit", "exact_jacobian.jl"))
+include(joinpath("prognostic_equations", "implicit", "debug_jacobian.jl"))
+include(joinpath("prognostic_equations", "implicit", "dual_fixes.jl"))
 
 include(joinpath("prognostic_equations", "remaining_tendency.jl"))
 include(joinpath("prognostic_equations", "forcing", "large_scale_advection.jl")) # TODO: should this be in tendencies/?
diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl
index 77ee03de1be..999dc2eea03 100644
--- a/src/cache/diagnostic_edmf_precomputed_quantities.jl
+++ b/src/cache/diagnostic_edmf_precomputed_quantities.jl
@@ -564,7 +564,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
                     ρ_prev_level
                 ) / local_geometry_level.J / ρ_level
 
-            @. detrʲ_prev_level = detrainment(
+            @. detrʲ_prev_level = detrainment_from_thermo_state(
                 params,
                 z_prev_level,
                 z_sfc_halflevel,
diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl
index c109e80df46..1903c5c5aa6 100644
--- a/src/cache/precomputed_quantities.jl
+++ b/src/cache/precomputed_quantities.jl
@@ -4,6 +4,43 @@
 import Thermodynamics as TD
 import ClimaCore: Spaces, Fields
 
+function implicit_precomputed_quantities(Y, atmos)
+    FT = eltype(Y)
+    TST = thermo_state_type(atmos.moisture_model, FT)
+    n = n_mass_flux_subdomains(atmos.turbconv_model)
+    gs_quantities = (;
+        ᶜspecific = specific_gs.(Y.c),
+        ᶠuₕ³ = similar(Y.f, CT3{FT}),
+        ᶜu = similar(Y.c, C123{FT}),
+        ᶠu³ = similar(Y.f, CT3{FT}),
+        ᶜK = similar(Y.c, FT),
+        ᶜts = similar(Y.c, TST),
+        ᶜp = similar(Y.c, FT),
+        ᶜh_tot = similar(Y.c, FT),
+    )
+    advective_sgs_quantities =
+        atmos.turbconv_model isa PrognosticEDMFX ?
+        (;
+            ᶜtke⁰ = similar(Y.c, FT),
+            ᶜρa⁰ = similar(Y.c, FT),
+            ᶜmse⁰ = similar(Y.c, FT),
+            ᶜq_tot⁰ = similar(Y.c, FT),
+            ᶠu₃⁰ = similar(Y.f, C3{FT}),
+            ᶜu⁰ = similar(Y.c, C123{FT}),
+            ᶠu³⁰ = similar(Y.f, CT3{FT}),
+            ᶜK⁰ = similar(Y.c, FT),
+            ᶜts⁰ = similar(Y.c, TST),
+            ᶜρ⁰ = similar(Y.c, FT),
+            ᶜuʲs = similar(Y.c, NTuple{n, C123{FT}}),
+            ᶠu³ʲs = similar(Y.f, NTuple{n, CT3{FT}}),
+            ᶜKʲs = similar(Y.c, NTuple{n, FT}),
+            ᶠKᵥʲs = similar(Y.f, NTuple{n, FT}),
+            ᶜtsʲs = similar(Y.c, NTuple{n, TST}),
+            ᶜρʲs = similar(Y.c, NTuple{n, FT}),
+        ) : (;)
+    return (; gs_quantities..., advective_sgs_quantities...)
+end
+
 """
     precomputed_quantities(Y, atmos)
 
@@ -42,16 +79,8 @@ function precomputed_quantities(Y, atmos)
     TST = thermo_state_type(atmos.moisture_model, FT)
     SCT = SurfaceConditions.surface_conditions_type(atmos, FT)
     cspace = axes(Y.c)
-    fspace = axes(Y.f)
     n = n_mass_flux_subdomains(atmos.turbconv_model)
     gs_quantities = (;
-        ᶜspecific = specific_gs.(Y.c),
-        ᶜu = similar(Y.c, C123{FT}),
-        ᶠu³ = similar(Y.f, CT3{FT}),
-        ᶜK = similar(Y.c, FT),
-        ᶜts = similar(Y.c, TST),
-        ᶜp = similar(Y.c, FT),
-        ᶜh_tot = similar(Y.c, FT),
         ᶜmixing_length = similar(Y.c, FT),
         ᶜlinear_buoygrad = similar(Y.c, FT),
         ᶜstrain_rate_norm = similar(Y.c, FT),
@@ -76,42 +105,26 @@ function precomputed_quantities(Y, atmos)
     advective_sgs_quantities =
         atmos.turbconv_model isa PrognosticEDMFX ?
         (;
-            ᶜtke⁰ = similar(Y.c, FT),
-            ᶜρa⁰ = similar(Y.c, FT),
-            ᶠu₃⁰ = similar(Y.f, C3{FT}),
-            ᶜu⁰ = similar(Y.c, C123{FT}),
-            ᶠu³⁰ = similar(Y.f, CT3{FT}),
-            ᶜK⁰ = similar(Y.c, FT),
-            ᶜmse⁰ = similar(Y.c, FT),
-            ᶜq_tot⁰ = similar(Y.c, FT),
-            ᶜts⁰ = similar(Y.c, TST),
-            ᶜρ⁰ = similar(Y.c, FT),
             ᶜmixing_length_tuple = similar(Y.c, MixingLength{FT}),
             ᶜK_u = similar(Y.c, FT),
             ᶜK_h = similar(Y.c, FT),
             ρatke_flux = similar(Fields.level(Y.f, half), C3{FT}),
-            ᶜuʲs = similar(Y.c, NTuple{n, C123{FT}}),
-            ᶠu³ʲs = similar(Y.f, NTuple{n, CT3{FT}}),
-            ᶜKʲs = similar(Y.c, NTuple{n, FT}),
-            ᶠKᵥʲs = similar(Y.f, NTuple{n, FT}),
-            ᶜtsʲs = similar(Y.c, NTuple{n, TST}),
             bdmr_l = similar(Y.c, BidiagonalMatrixRow{FT}),
             bdmr_r = similar(Y.c, BidiagonalMatrixRow{FT}),
             bdmr = similar(Y.c, BidiagonalMatrixRow{FT}),
-            ᶜρʲs = similar(Y.c, NTuple{n, FT}),
             ᶜentrʲs = similar(Y.c, NTuple{n, FT}),
             ᶜdetrʲs = similar(Y.c, NTuple{n, FT}),
             ᶜturb_entrʲs = similar(Y.c, NTuple{n, FT}),
             ᶠnh_pressure₃ʲs = similar(Y.f, NTuple{n, C3{FT}}),
-            ᶜgradᵥ_θ_virt⁰ = Fields.Field(C3{FT}, cspace),
-            ᶜgradᵥ_q_tot⁰ = Fields.Field(C3{FT}, cspace),
-            ᶜgradᵥ_θ_liq_ice⁰ = Fields.Field(C3{FT}, cspace),
+            ᶜgradᵥ_θ_virt⁰ = similar(Y.c, C3{FT}),
+            ᶜgradᵥ_q_tot⁰ = similar(Y.c, C3{FT}),
+            ᶜgradᵥ_θ_liq_ice⁰ = similar(Y.c, C3{FT}),
             precipitation_sgs_quantities...,
         ) : (;)
     sgs_quantities = (;
-        ᶜgradᵥ_θ_virt = Fields.Field(C3{FT}, cspace),
-        ᶜgradᵥ_q_tot = Fields.Field(C3{FT}, cspace),
-        ᶜgradᵥ_θ_liq_ice = Fields.Field(C3{FT}, cspace),
+        ᶜgradᵥ_θ_virt = similar(Y.c, C3{FT}),
+        ᶜgradᵥ_q_tot = similar(Y.c, C3{FT}),
+        ᶜgradᵥ_θ_liq_ice = similar(Y.c, C3{FT}),
     )
 
     diagnostic_sgs_quantities =
@@ -157,6 +170,7 @@ function precomputed_quantities(Y, atmos)
             ᶜqₛ = similar(Y.c, FT),
         ) : (;)
     return (;
+        implicit_precomputed_quantities(Y, atmos)...,
         gs_quantities...,
         sgs_quantities...,
         advective_sgs_quantities...,
@@ -328,14 +342,14 @@ end
 function eddy_diffusivity_coefficient(
     z::FT,
     z₀,
-    f_b::FT,
-    h::FT,
+    f_b,
+    h,
     uₐ,
-    C_E::FT,
-    Ri::FT,
-    Ri_a::FT,
-    Ri_c::FT,
-    κ::FT,
+    C_E,
+    Ri,
+    Ri_a,
+    Ri_c,
+    κ,
 ) where {FT}
     # Equations (17), (18)
     if z <= f_b * h
@@ -414,12 +428,12 @@ end
 
 function compute_surface_layer_diffusivity(
     z::FT,
-    z₀::FT,
-    κ::FT,
-    C_E::FT,
-    Ri::FT,
-    Ri_a::FT,
-    Ri_c::FT,
+    z₀,
+    κ,
+    C_E,
+    Ri,
+    Ri_a,
+    Ri_c,
     norm_uₐ,
 ) where {FT}
     # Equations (19), (20)
@@ -454,15 +468,19 @@ function instead of recomputing the value yourself. Otherwise, it will be
 difficult to ensure that the duplicated computations are consistent.
 """
 NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
+    set_implicit_precomputed_quantities!(Y, p, t)
+    set_explicit_precomputed_quantities!(Y, p, t)
+end
+function set_implicit_precomputed_quantities!(Y, p, t)
     (; moisture_model, turbconv_model, vert_diff, precip_model, cloud_model) =
         p.atmos
     (; call_cloud_diagnostics_per_stage) = p.atmos
     thermo_params = CAP.thermodynamics_params(p.params)
+    turbconv_params = CAP.turbconv_params(p.params)
     n = n_mass_flux_subdomains(turbconv_model)
     thermo_args = (thermo_params, moisture_model)
     (; ᶜΦ) = p.core
-    (; ᶜspecific, ᶜu, ᶠu³, ᶜK, ᶜts, ᶜp) = p.precomputed
-    ᶠuₕ³ = p.scratch.ᶠtemp_CT3
+    (; ᶜspecific, ᶠuₕ³, ᶜu, ᶠu³, ᶜK, ᶜts, ᶜp, ᶜh_tot) = p.precomputed
 
     @. ᶜspecific = specific_gs(Y.c)
     set_ᶠuₕ³!(ᶠuₕ³, Y)
@@ -490,6 +508,63 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
     end
     @. ᶜts = ts_gs(thermo_args..., ᶜspecific, ᶜK, ᶜΦ, Y.c.ρ)
     @. ᶜp = TD.air_pressure(thermo_params, ᶜts)
+    @. ᶜh_tot = TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜspecific.e_tot)
+
+    if turbconv_model isa PrognosticEDMFX
+        @assert !(moisture_model isa DryModel)
+        (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶠKᵥʲs, ᶜtsʲs, ᶜρʲs) = p.precomputed
+        (; ᶜtke⁰, ᶜρa⁰, ᶜmse⁰, ᶜq_tot⁰, ᶠu₃⁰, ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶜts⁰, ᶜρ⁰) =
+            p.precomputed
+
+        for j in 1:n
+            ᶠu₃ʲ = Y.f.sgsʲs.:($j).u₃
+            ᶜmseʲ = Y.c.sgsʲs.:($j).mse
+            ᶜq_totʲ = Y.c.sgsʲs.:($j).q_tot
+
+            ᶜuʲ = ᶜuʲs.:($j)
+            ᶠu³ʲ = ᶠu³ʲs.:($j)
+            ᶜKʲ = ᶜKʲs.:($j)
+            ᶠKᵥʲ = ᶠKᵥʲs.:($j)
+            ᶜtsʲ = ᶜtsʲs.:($j)
+            ᶜρʲ = ᶜρʲs.:($j)
+
+            set_velocity_quantities!(ᶜuʲ, ᶠu³ʲ, ᶜKʲ, ᶠu₃ʲ, Y.c.uₕ, ᶠuₕ³)
+            @. ᶠKᵥʲ = (adjoint(CT3(ᶠu₃ʲ)) * ᶠu₃ʲ) / 2
+            @. ᶜtsʲ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜmseʲ - ᶜΦ, ᶜq_totʲ)
+            @. ᶜρʲ = TD.air_density(thermo_params, ᶜtsʲ)
+        end
+
+        @. ᶜρa⁰ = ρa⁰(Y.c)
+        @. ᶜtke⁰ = divide_by_ρa(Y.c.sgs⁰.ρatke, ᶜρa⁰, 0, Y.c.ρ, turbconv_model)
+        @. ᶜmse⁰ = divide_by_ρa(
+            Y.c.ρ * (ᶜh_tot - ᶜK) - ρamse⁺(Y.c.sgsʲs),
+            ᶜρa⁰,
+            Y.c.ρ * (ᶜh_tot - ᶜK),
+            Y.c.ρ,
+            turbconv_model,
+        )
+        @. ᶜq_tot⁰ = divide_by_ρa(
+            Y.c.ρq_tot - ρaq_tot⁺(Y.c.sgsʲs),
+            ᶜρa⁰,
+            Y.c.ρq_tot,
+            Y.c.ρ,
+            turbconv_model,
+        )
+        set_sgs_ᶠu₃!(u₃⁰, ᶠu₃⁰, Y, turbconv_model)
+        set_velocity_quantities!(ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶠu₃⁰, Y.c.uₕ, ᶠuₕ³)
+        # @. ᶜK⁰ += ᶜtke⁰
+        @. ᶜts⁰ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜmse⁰ - ᶜΦ, ᶜq_tot⁰)
+        @. ᶜρ⁰ = TD.air_density(thermo_params, ᶜts⁰)
+    end
+
+    return nothing
+end
+function set_explicit_precomputed_quantities!(Y, p, t)
+    (; moisture_model, turbconv_model, vert_diff, precip_model, cloud_model) =
+        p.atmos
+    (; call_cloud_diagnostics_per_stage) = p.atmos
+    thermo_params = CAP.thermodynamics_params(p.params)
+    (; ᶠuₕ³, ᶜts, ᶜp) = p.precomputed
 
     if turbconv_model isa AbstractEDMF
         @. p.precomputed.ᶜgradᵥ_θ_virt =
@@ -500,9 +575,6 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
             ᶜgradᵥ(ᶠinterp(TD.liquid_ice_pottemp(thermo_params, ᶜts)))
     end
 
-    (; ᶜh_tot) = p.precomputed
-    @. ᶜh_tot = TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜspecific.e_tot)
-
     if !isnothing(p.sfc_setup)
         SurfaceConditions.update_surface_conditions!(Y, p, t)
     end
@@ -519,12 +591,11 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
 
     if turbconv_model isa PrognosticEDMFX
         set_prognostic_edmf_precomputed_quantities_draft_and_bc!(Y, p, ᶠuₕ³, t)
-        set_prognostic_edmf_precomputed_quantities_environment!(Y, p, ᶠuₕ³, t)
         set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
         set_prognostic_edmf_precomputed_quantities_precipitation!(
             Y,
             p,
-            p.atmos.precip_model,
+            precip_model,
         )
     end
 
@@ -537,7 +608,7 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
             Y,
             p,
             t,
-            p.atmos.precip_model,
+            precip_model,
         )
     end
 
@@ -546,7 +617,7 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
         interior_uₕ = Fields.level(Y.c.uₕ, 1)
         ᶜΔz_surface = Fields.Δz_field(interior_uₕ)
         @. ᶜK_h = eddy_diffusivity_coefficient(
-            p.atmos.vert_diff.C_E,
+            vert_diff.C_E,
             norm(interior_uₕ),
             ᶜΔz_surface / 2,
             ᶜp,
diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl
index eb620888bb9..4733d532644 100644
--- a/src/cache/prognostic_edmf_precomputed_quantities.jl
+++ b/src/cache/prognostic_edmf_precomputed_quantities.jl
@@ -5,50 +5,6 @@ import NVTX
 import Thermodynamics as TD
 import ClimaCore: Spaces, Fields
 
-"""
-    set_prognostic_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
-
-Updates the edmf environment precomputed quantities stored in `p` for edmfx.
-"""
-NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_environment!(
-    Y,
-    p,
-    ᶠuₕ³,
-    t,
-)
-    @assert !(p.atmos.moisture_model isa DryModel)
-
-    thermo_params = CAP.thermodynamics_params(p.params)
-    (; turbconv_model) = p.atmos
-    (; ᶜΦ,) = p.core
-    (; ᶜp, ᶜh_tot, ᶜK) = p.precomputed
-    (; ᶜtke⁰, ᶜρa⁰, ᶠu₃⁰, ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶜts⁰, ᶜρ⁰, ᶜmse⁰, ᶜq_tot⁰) =
-        p.precomputed
-
-    @. ᶜρa⁰ = ρa⁰(Y.c)
-    @. ᶜtke⁰ = divide_by_ρa(Y.c.sgs⁰.ρatke, ᶜρa⁰, 0, Y.c.ρ, turbconv_model)
-    @. ᶜmse⁰ = divide_by_ρa(
-        Y.c.ρ * (ᶜh_tot - ᶜK) - ρamse⁺(Y.c.sgsʲs),
-        ᶜρa⁰,
-        Y.c.ρ * (ᶜh_tot - ᶜK),
-        Y.c.ρ,
-        turbconv_model,
-    )
-    @. ᶜq_tot⁰ = divide_by_ρa(
-        Y.c.ρq_tot - ρaq_tot⁺(Y.c.sgsʲs),
-        ᶜρa⁰,
-        Y.c.ρq_tot,
-        Y.c.ρ,
-        turbconv_model,
-    )
-    set_sgs_ᶠu₃!(u₃⁰, ᶠu₃⁰, Y, turbconv_model)
-    set_velocity_quantities!(ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶠu₃⁰, Y.c.uₕ, ᶠuₕ³)
-    # @. ᶜK⁰ += ᶜtke⁰
-    @. ᶜts⁰ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜmse⁰ - ᶜΦ, ᶜq_tot⁰)
-    @. ᶜρ⁰ = TD.air_density(thermo_params, ᶜts⁰)
-    return nothing
-end
-
 """
     set_prognostic_edmf_precomputed_quantities_draft_and_bc!(Y, p, ᶠuₕ³, t)
 
@@ -62,8 +18,6 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_draft_and_bc!
     t,
 )
     (; moisture_model, turbconv_model) = p.atmos
-    #EDMFX BCs only support total energy as state variable
-    @assert !(moisture_model isa DryModel)
 
     FT = Spaces.undertype(axes(Y.c))
     n = n_mass_flux_subdomains(turbconv_model)
@@ -78,21 +32,10 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_draft_and_bc!
     (; ustar, obukhov_length, buoyancy_flux) = p.precomputed.sfc_conditions
 
     for j in 1:n
-        ᶜuʲ = ᶜuʲs.:($j)
-        ᶠu³ʲ = ᶠu³ʲs.:($j)
-        ᶜKʲ = ᶜKʲs.:($j)
-        ᶠKᵥʲ = ᶠKᵥʲs.:($j)
-        ᶠu₃ʲ = Y.f.sgsʲs.:($j).u₃
         ᶜtsʲ = ᶜtsʲs.:($j)
-        ᶜρʲ = ᶜρʲs.:($j)
         ᶜmseʲ = Y.c.sgsʲs.:($j).mse
         ᶜq_totʲ = Y.c.sgsʲs.:($j).q_tot
 
-        set_velocity_quantities!(ᶜuʲ, ᶠu³ʲ, ᶜKʲ, ᶠu₃ʲ, Y.c.uₕ, ᶠuₕ³)
-        @. ᶠKᵥʲ = (adjoint(CT3(ᶠu₃ʲ)) * ᶠu₃ʲ) / 2
-        @. ᶜtsʲ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜmseʲ - ᶜΦ, ᶜq_totʲ)
-        @. ᶜρʲ = TD.air_density(thermo_params, ᶜtsʲ)
-
         # EDMFX boundary condition:
 
         # We need field_values everywhere because we are mixing
diff --git a/src/cache/temporary_quantities.jl b/src/cache/temporary_quantities.jl
index 4482e861475..2943f524051 100644
--- a/src/cache/temporary_quantities.jl
+++ b/src/cache/temporary_quantities.jl
@@ -10,6 +10,7 @@ function temporary_quantities(Y, atmos)
 
     FT = Spaces.undertype(center_space)
     CTh = CTh_vector_type(Y.c)
+    num_maxtrix_rows = length(first(column_iterator(Y)))
     return (;
         ᶠtemp_scalar = Fields.Field(FT, face_space), # ᶠp, ᶠρK_E
         ᶜtemp_scalar = Fields.Field(FT, center_space), # ᶜ1
@@ -68,5 +69,6 @@ function temporary_quantities(Y, atmos)
         ᶜdiffusion_h_matrix = similar(Y.c, TridiagonalMatrixRow{FT}),
         ᶜdiffusion_u_matrix = similar(Y.c, TridiagonalMatrixRow{FT}),
         ᶠtridiagonal_matrix_c3 = similar(Y.f, TridiagonalMatrixRow{C3{FT}}),
+        temp_matrix = Array{FT}(undef, num_maxtrix_rows, num_maxtrix_rows),
     )
 end
diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl
index 56bbb7076a8..e9cf9b021fb 100644
--- a/src/callbacks/callbacks.jl
+++ b/src/callbacks/callbacks.jl
@@ -404,7 +404,7 @@ function print_walltime_estimate(integrator)
             ""
         end
         estimated_finish_date =
-            Dates.now() + compound_period(wall_time_remaining, Dates.Second)
+            Dates.now() + compound_period(wall_time_remaining)
         @info "Progress" simulation_time = simulation_time n_steps_completed =
             n_steps wall_time_per_step = wall_time_ave_per_step_str wall_time_total =
             wall_time_total wall_time_remaining = wall_time_remaining_str wall_time_spent =
@@ -452,6 +452,32 @@ function gc_func(integrator)
     return nothing
 end
 
+function update_exact_jacobian!(integrator)
+    integrator.alg isa CTS.RosenbrockAlgorithm ||
+        integrator.alg.name isa CTS.IMEXARKAlgorithmName ||
+        return
+    (; tableau) = integrator.alg
+    tableau_coefficients = LinearAlgebra.diag(
+        integrator.alg isa CTS.RosenbrockAlgorithm ? tableau.Γ : tableau.a_imp,
+    )
+    γs = unique(filter(!iszero, tableau_coefficients))
+    length(γs) == 1 ||
+        error("The exact Jacobian must be updated on every Newton iteration, \
+               rather than on every timestep (or every N steps), because the \
+               specified IMEX algorithm has implicit stages with distinct \
+               tableau coefficients (i.e., it is not an SDIRK algorithm).")
+    j =
+        integrator.alg isa CTS.RosenbrockAlgorithm ? integrator.cache.W :
+        integrator.cache.newtons_method_cache.j
+    update_exact_jacobian!(
+        j,
+        integrator.u,
+        integrator.p,
+        integrator.dt * γs[1],
+        integrator.t,
+    )
+end
+
 """
     maybe_graceful_exit(integrator)
 
diff --git a/src/callbacks/get_callbacks.jl b/src/callbacks/get_callbacks.jl
index 45ac6115301..b700d224aa0 100644
--- a/src/callbacks/get_callbacks.jl
+++ b/src/callbacks/get_callbacks.jl
@@ -23,6 +23,7 @@ function get_diagnostics(parsed_args, atmos_model, Y, p, dt)
         "average" => ((+), CAD.average_pre_output_hook!),
     )
 
+    dict_writer = CAD.DictWriter()
     hdf5_writer = CAD.HDF5Writer(p.output_dir)
 
     if !isnothing(parsed_args["netcdf_interpolation_num_points"])
@@ -45,11 +46,13 @@ function get_diagnostics(parsed_args, atmos_model, Y, p, dt)
         z_sampling_method,
         sync_schedule = CAD.EveryStepSchedule(),
     )
-    writers = (hdf5_writer, netcdf_writer)
 
-    # The default writer is HDF5
+    writers = (; dict_writer, hdf5_writer, netcdf_writer)
+
+    # The default writer is NetCDF
     ALLOWED_WRITERS = Dict(
         "nothing" => netcdf_writer,
+        "dict" => dict_writer,
         "h5" => hdf5_writer,
         "hdf5" => hdf5_writer,
         "nc" => netcdf_writer,
@@ -157,7 +160,30 @@ function get_diagnostics(parsed_args, atmos_model, Y, p, dt)
             diagnostics...,
         ]
     end
-    diagnostics = collect(diagnostics)
+    if parsed_args["output_default_diagnostics"] && (
+        parsed_args["use_exact_jacobian"] ||
+        parsed_args["debug_approximate_jacobian"]
+    )
+        period_seconds =
+            min(
+                max(
+                    1,
+                    parsed_args["n_steps_update_exact_jacobian"],
+                    fld(p.t_end / 99, p.dt), # Save no more than 100 matrices.
+                ),
+                fld(p.t_end, p.dt), # Save at least 2 matrices.
+            ) * p.dt
+        schedule = CAD.EveryDtSchedule(period_seconds)
+        exact_jacobian_diagnostic = CAD.ScheduledDiagnostic(;
+            variable = CAD.get_diagnostic_variable("ejac"),
+            output_schedule_func = schedule,
+            compute_schedule_func = deepcopy(schedule),
+            reduction_time_func = nothing,
+            pre_output_hook! = nothing,
+            output_writer = dict_writer,
+        )
+        diagnostics = [exact_jacobian_diagnostic, diagnostics...]
+    end
 
     for writer in writers
         writer_str = nameof(typeof(writer))
@@ -271,6 +297,19 @@ function get_callbacks(config, sim_info, atmos, params, Y, p, t_start)
             (callbacks..., call_every_dt(cloud_fraction_model_callback!, dt_cf))
     end
 
+    if (
+        parsed_args["use_exact_jacobian"] ||
+        parsed_args["debug_approximate_jacobian"]
+    ) && parsed_args["n_steps_update_exact_jacobian"] != 0
+        callbacks = (
+            callbacks...,
+            call_every_n_steps(
+                update_exact_jacobian!,
+                parsed_args["n_steps_update_exact_jacobian"],
+            ),
+        )
+    end
+
     if atmos.radiation_mode isa RRTMGPI.AbstractRRTMGPMode
         dt_rad = FT(time_to_seconds(parsed_args["dt_rad"]))
         callbacks =
diff --git a/src/diagnostics/Diagnostics.jl b/src/diagnostics/Diagnostics.jl
index 2c8a5c19281..4c8e8f008e6 100644
--- a/src/diagnostics/Diagnostics.jl
+++ b/src/diagnostics/Diagnostics.jl
@@ -56,6 +56,7 @@ import ClimaDiagnostics.Schedules:
     EveryStepSchedule, EveryDtSchedule, EveryCalendarDtSchedule
 
 import ClimaDiagnostics.Writers:
+    DictWriter,
     HDF5Writer,
     NetCDFWriter,
     write_field!,
diff --git a/src/diagnostics/core_diagnostics.jl b/src/diagnostics/core_diagnostics.jl
index ab76d8a4a80..c00fba3d45f 100644
--- a/src/diagnostics/core_diagnostics.jl
+++ b/src/diagnostics/core_diagnostics.jl
@@ -1163,3 +1163,19 @@ add_diagnostic_variable!(
     comments = "Mass of water vapor per mass of air",
     compute! = compute_husv!,
 )
+
+###
+# Implicit solver
+###
+add_diagnostic_variable!(
+    short_name = "ejac",
+    long_name = "Exact Jacobian matrix of first column",
+    standard_name = "exact_jacobian",
+    units = "s^-1",
+    comments = "Exact Jacobian matrix for column (1, 1, 1) of tendency",
+    compute! = (out, state, cache, time) -> if isnothing(out)
+        copy(cache.scratch.temp_matrix)
+    else
+        out .= cache.scratch.temp_matrix
+    end,
+)
diff --git a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl
index ac118ef08df..6880ffd090d 100644
--- a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl
+++ b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl
@@ -108,14 +108,8 @@ end
 Returns the total energy source term multiplier from precipitation formation
 for the 0-moment scheme
 """
-function e_tot_0M_precipitation_sources_helper(
-    thp,
-    ts,
-    Φ::FT,
-) where {FT <: Real}
-
-    λ::FT = TD.liquid_fraction(thp, ts)
-
+function e_tot_0M_precipitation_sources_helper(thp, ts, Φ)
+    λ = TD.liquid_fraction(thp, ts)
     return λ * Iₗ(thp, ts) + (1 - λ) * Iᵢ(thp, ts) + Φ
 end
 
diff --git a/src/prognostic_equations/edmfx_boundary_condition.jl b/src/prognostic_equations/edmfx_boundary_condition.jl
index 20810c9df9a..5d02d96799e 100644
--- a/src/prognostic_equations/edmfx_boundary_condition.jl
+++ b/src/prognostic_equations/edmfx_boundary_condition.jl
@@ -4,9 +4,9 @@
 
 function sgs_scalar_first_interior_bc(
     ᶜz_int::FT,
-    ᶜρ_int::FT,
-    ᶜaʲ_int::FT,
-    ᶜscalar_int::FT,
+    ᶜρ_int,
+    ᶜaʲ_int,
+    ᶜscalar_int,
     sfc_buoyancy_flux,
     sfc_ρ_flux_scalar,
     ustar,
@@ -32,9 +32,9 @@ end
 
 function get_first_interior_variance(
     flux,
-    ustar::FT,
+    ustar,
     z::FT,
-    oblength::FT,
+    oblength,
     local_geometry,
 ) where {FT}
     c_star = -flux / ustar
@@ -61,13 +61,10 @@ function guass_quantile(p::FT) where {FT <: Real}
     return sqrt(2) * approximate_inverf(2p - 1)
 end
 
-function percentile_bounds_mean_norm(
-    low_percentile::FT,
-    high_percentile::FT,
-) where {FT <: Real}
+function percentile_bounds_mean_norm(low_percentile, high_percentile)
     gauss_int(x) = -exp(-x * x / 2) / sqrt(2 * pi)
     xp_high = guass_quantile(high_percentile)
     xp_low = guass_quantile(low_percentile)
     return (gauss_int(xp_high) - gauss_int(xp_low)) /
-           max(high_percentile - low_percentile, eps(FT))
+           max(high_percentile - low_percentile, eps(typeof(high_percentile)))
 end
diff --git a/src/prognostic_equations/edmfx_closures.jl b/src/prognostic_equations/edmfx_closures.jl
index 018385d1894..2d15310b02b 100644
--- a/src/prognostic_equations/edmfx_closures.jl
+++ b/src/prognostic_equations/edmfx_closures.jl
@@ -153,17 +153,17 @@ Smagorinsky length scale.
 """
 function mixing_length(
     params,
-    ustar::FT,
+    ustar,
     ᶜz::FT,
-    z_sfc::FT,
-    ᶜdz::FT,
-    sfc_tke::FT,
-    ᶜlinear_buoygrad::FT,
-    ᶜtke::FT,
-    obukhov_length::FT,
-    ᶜstrain_rate_norm::FT,
-    ᶜPr::FT,
-    ᶜtke_exch::FT,
+    z_sfc,
+    ᶜdz,
+    sfc_tke,
+    ᶜlinear_buoygrad,
+    ᶜtke,
+    obukhov_length,
+    ᶜstrain_rate_norm,
+    ᶜPr,
+    ᶜtke_exch,
 ) where {FT}
 
     turbconv_params = CAP.turbconv_params(params)
@@ -234,7 +234,7 @@ function mixing_length(
     l_smin = lamb_smooth_minimum(l, smin_ub, smin_rm)
     l_limited = max(l_smag, min(l_smin, l_z))
 
-    return MixingLength{FT}(l_limited, l_W, l_TKE, l_N)
+    return MixingLength(l_limited, l_W, l_TKE, l_N)
 end
 
 """
diff --git a/src/prognostic_equations/edmfx_entr_detr.jl b/src/prognostic_equations/edmfx_entr_detr.jl
index 0d37e891e0f..2b3148e12e4 100644
--- a/src/prognostic_equations/edmfx_entr_detr.jl
+++ b/src/prognostic_equations/edmfx_entr_detr.jl
@@ -19,17 +19,17 @@
 function entrainment(
     params,
     ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    ᶜtke⁰::FT,
+    z_sfc,
+    ᶜp,
+    ᶜρ,
+    ᶜaʲ,
+    ᶜwʲ,
+    ᶜRHʲ,
+    ᶜbuoyʲ,
+    ᶜw⁰,
+    ᶜRH⁰,
+    ᶜbuoy⁰,
+    ᶜtke⁰,
     ::NoEntrainment,
 ) where {FT}
     return FT(0)
@@ -38,17 +38,17 @@ end
 function entrainment(
     params,
     ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    ᶜtke⁰::FT,
+    z_sfc,
+    ᶜp,
+    ᶜρ,
+    ᶜaʲ,
+    ᶜwʲ,
+    ᶜRHʲ,
+    ᶜbuoyʲ,
+    ᶜw⁰,
+    ᶜRH⁰,
+    ᶜbuoy⁰,
+    ᶜtke⁰,
     ::PiGroupsEntrainment,
 ) where {FT}
     if ᶜaʲ <= FT(0)
@@ -86,17 +86,17 @@ end
 function entrainment(
     params,
     ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    ᶜtke⁰::FT,
+    z_sfc,
+    ᶜp,
+    ᶜρ,
+    ᶜaʲ,
+    ᶜwʲ,
+    ᶜRHʲ,
+    ᶜbuoyʲ,
+    ᶜw⁰,
+    ᶜRH⁰,
+    ᶜbuoy⁰,
+    ᶜtke⁰,
     ::GeneralizedEntrainment,
 ) where {FT}
     turbconv_params = CAP.turbconv_params(params)
@@ -117,7 +117,7 @@ function entrainment(
     return max(entr, 0)
 end
 
-function detrainment(
+function detrainment_from_thermo_state(
     params,
     z_prev_level::FT,
     z_sfc_halflevel,
@@ -177,21 +177,21 @@ end
 function detrainment(
     params,
     ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    ᶜρaʲ::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    ᶜentr::FT,
-    ᶜvert_div::FT,
-    ᶜmassflux_vert_div::FT,
-    ᶜtke⁰::FT,
+    z_sfc,
+    ᶜp,
+    ᶜρ,
+    ᶜρaʲ,
+    ᶜaʲ,
+    ᶜwʲ,
+    ᶜRHʲ,
+    ᶜbuoyʲ,
+    ᶜw⁰,
+    ᶜRH⁰,
+    ᶜbuoy⁰,
+    ᶜentr,
+    ᶜvert_div,
+    ᶜmassflux_vert_div,
+    ᶜtke⁰,
     ::NoDetrainment,
 ) where {FT}
     return FT(0)
@@ -200,21 +200,21 @@ end
 function detrainment(
     params,
     ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    ᶜρaʲ::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    ᶜentr::FT,
-    ᶜvert_div::FT,
-    ᶜmassflux_vert_div::FT,
-    ᶜtke⁰::FT,
+    z_sfc,
+    ᶜp,
+    ᶜρ,
+    ᶜρaʲ,
+    ᶜaʲ,
+    ᶜwʲ,
+    ᶜRHʲ,
+    ᶜbuoyʲ,
+    ᶜw⁰,
+    ᶜRH⁰,
+    ᶜbuoy⁰,
+    ᶜentr,
+    ᶜvert_div,
+    ᶜmassflux_vert_div,
+    ᶜtke⁰,
     ::PiGroupsDetrainment,
 ) where {FT}
 
@@ -251,21 +251,21 @@ end
 function detrainment(
     params,
     ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    ᶜρaʲ::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    ᶜentr::FT,
-    ᶜvert_div::FT,
-    ᶜmassflux_vert_div::FT,
-    ᶜtke⁰::FT,
+    z_sfc,
+    ᶜp,
+    ᶜρ,
+    ᶜρaʲ,
+    ᶜaʲ,
+    ᶜwʲ,
+    ᶜRHʲ,
+    ᶜbuoyʲ,
+    ᶜw⁰,
+    ᶜRH⁰,
+    ᶜbuoy⁰,
+    ᶜentr,
+    ᶜvert_div,
+    ᶜmassflux_vert_div,
+    ᶜtke⁰,
     ::GeneralizedDetrainment,
 ) where {FT}
     turbconv_params = CAP.turbconv_params(params)
@@ -301,20 +301,20 @@ end
 function detrainment(
     params,
     ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    ᶜρaʲ::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    ᶜentr::FT,
-    ᶜvert_div::FT,
-    ᶜmassflux_vert_div::FT,
+    z_sfc,
+    ᶜp,
+    ᶜρ,
+    ᶜρaʲ,
+    ᶜaʲ,
+    ᶜwʲ,
+    ᶜRHʲ,
+    ᶜbuoyʲ,
+    ᶜw⁰,
+    ᶜRH⁰,
+    ᶜbuoy⁰,
+    ᶜentr,
+    ᶜvert_div,
+    ᶜmassflux_vert_div,
     ::ConstantAreaDetrainment,
 ) where {FT}
     detr = ᶜentr - ᶜvert_div
diff --git a/src/prognostic_equations/gm_sgs_closures.jl b/src/prognostic_equations/gm_sgs_closures.jl
index 13730270ed5..c6af7fea766 100644
--- a/src/prognostic_equations/gm_sgs_closures.jl
+++ b/src/prognostic_equations/gm_sgs_closures.jl
@@ -19,7 +19,7 @@ function smagorinsky_lilly_length(c_smag, N_eff, dz, Pr, ϵ_st)
     return N_eff > FT(0) ?
            c_smag *
            dz *
-           max(0, 1 - N_eff^2 / Pr / 2 / max(ϵ_st, eps(FT)))^(1 / 4) :
+           max(FT(0), 1 - N_eff^2 / Pr / 2 / max(ϵ_st, eps(FT)))^(1 / FT(4)) :
            c_smag * dz
 end
 
diff --git a/src/prognostic_equations/implicit/approx_jacobian.jl b/src/prognostic_equations/implicit/approx_jacobian.jl
new file mode 100644
index 00000000000..079f13cd25a
--- /dev/null
+++ b/src/prognostic_equations/implicit/approx_jacobian.jl
@@ -0,0 +1,683 @@
+abstract type DerivativeFlag end
+struct UseDerivative <: DerivativeFlag end
+struct IgnoreDerivative <: DerivativeFlag end
+use_derivative(::UseDerivative) = true
+use_derivative(::IgnoreDerivative) = false
+
+"""
+    ApproxJacobian(Y, p; approximate_solve_iters, diffusion_flag, topography_flag)
+
+An approximation of the exact `ImplicitEquationJacobian`, which is updated using
+analytically derived tendency derivatives and inverted using a specialized
+nested linear solver.
+
+# Keyword Arguments
+
+- `diffusion_flag::DerivativeFlag`: whether the derivative of the
+  diffusion tendency with respect to the quantities being diffused should be
+  computed or approximated as 0; must be either `UseDerivative()` or
+  `Ignoreerivative()` instead of a `Bool` to ensure type-stability
+- `topography_flag::DerivativeFlag`: whether the derivative of vertical
+  contravariant velocity with respect to horizontal covariant velocity should be
+  computed or approximated as 0; must be either `UseDerivative()` or
+  `IgnoreDerivative()` instead of a `Bool` to ensure type-stability
+- `sgs_advection_flag::DerivativeFlag`: whether the derivative of the
+  subgrid-scale advection tendency with respect to the updraft quantities should
+  be computed or approximated as 0; must be either `UseDerivative()` or
+  `IgnoreDerivative()` instead of a `Bool` to ensure type-stability
+- `approximate_solve_iters::Int`: number of iterations to take for the
+  approximate linear solve required when `diffusion_flag = UseDerivative()`
+"""
+@kwdef struct ApproxJacobian{F1, F2, F3} <: JacobianAlgorithm
+    diffusion_flag::F1 = nothing
+    topography_flag::F2 = nothing
+    sgs_advection_flag::F3 = nothing
+    approximate_solve_iters::Int = 1
+end
+
+function jacobian_cache(alg::ApproxJacobian, Y, p)
+    diffusion_flag = if isnothing(alg.diffusion_flag)
+        p.atmos.diff_mode == Implicit() ? UseDerivative() : IgnoreDerivative()
+    else
+        alg.diffusion_flag
+    end
+    topography_flag = if isnothing(alg.topography_flag)
+        has_topography(axes(Y.c)) ? UseDerivative() : IgnoreDerivative()
+    else
+        alg.topography_flag
+    end
+    sgs_advection_flag = if isnothing(alg.sgs_advection_flag)
+        p.atmos.sgs_adv_mode == Implicit() ? UseDerivative() : IgnoreDerivative()
+    else
+        alg.sgs_advection_flag
+    end
+
+    FT = Spaces.undertype(axes(Y.c))
+    CTh = CTh_vector_type(axes(Y.c))
+
+    DiagonalRow = DiagonalMatrixRow{FT}
+    TridiagonalRow = TridiagonalMatrixRow{FT}
+    BidiagonalRow_C3 = BidiagonalMatrixRow{C3{FT}}
+    TridiagonalRow_ACTh = TridiagonalMatrixRow{Adjoint{FT, CTh{FT}}}
+    BidiagonalRow_ACT3 = BidiagonalMatrixRow{Adjoint{FT, CT3{FT}}}
+    BidiagonalRow_C3xACTh =
+        BidiagonalMatrixRow{typeof(zero(C3{FT}) * zero(CTh{FT})')}
+    DiagonalRow_C3xACT3 =
+        DiagonalMatrixRow{typeof(zero(C3{FT}) * zero(CT3{FT})')}
+    TridiagonalRow_C3xACT3 =
+        TridiagonalMatrixRow{typeof(zero(C3{FT}) * zero(CT3{FT})')}
+
+    is_in_Y(name) = MatrixFields.has_field(Y, name)
+
+    ρq_tot_if_available = is_in_Y(@name(c.ρq_tot)) ? (@name(c.ρq_tot),) : ()
+    ρatke_if_available =
+        is_in_Y(@name(c.sgs⁰.ρatke)) ? (@name(c.sgs⁰.ρatke),) : ()
+    sfc_if_available = is_in_Y(@name(sfc)) ? (@name(sfc),) : ()
+
+    tracer_names = (
+        @name(c.ρq_tot),
+        @name(c.ρq_liq),
+        @name(c.ρq_ice),
+        @name(c.ρq_rai),
+        @name(c.ρq_sno),
+    )
+    available_tracer_names = MatrixFields.unrolled_filter(is_in_Y, tracer_names)
+
+    # Note: We have to use FT(-1) * I instead of -I because inv(-1) == -1.0,
+    # which means that multiplying inv(-1) by a Float32 will yield a Float64.
+    identity_blocks = MatrixFields.unrolled_map(
+        name -> (name, name) => FT(-1) * I,
+        (@name(c.ρ), sfc_if_available...),
+    )
+
+    active_scalar_names = (@name(c.ρ), @name(c.ρe_tot), ρq_tot_if_available...)
+    advection_blocks = (
+        (
+            use_derivative(topography_flag) ?
+            MatrixFields.unrolled_map(
+                name ->
+                    (name, @name(c.uₕ)) =>
+                        similar(Y.c, TridiagonalRow_ACTh),
+                active_scalar_names,
+            ) : ()
+        )...,
+        MatrixFields.unrolled_map(
+            name -> (name, @name(f.u₃)) => similar(Y.c, BidiagonalRow_ACT3),
+            active_scalar_names,
+        )...,
+        MatrixFields.unrolled_map(
+            name -> (@name(f.u₃), name) => similar(Y.f, BidiagonalRow_C3),
+            active_scalar_names,
+        )...,
+        (@name(f.u₃), @name(c.uₕ)) => similar(Y.f, BidiagonalRow_C3xACTh),
+        (@name(f.u₃), @name(f.u₃)) => similar(Y.f, TridiagonalRow_C3xACT3),
+    )
+
+    diffused_scalar_names =
+        (@name(c.ρe_tot), available_tracer_names..., ρatke_if_available...)
+    diffusion_blocks = if use_derivative(diffusion_flag)
+        (
+            MatrixFields.unrolled_map(
+                name -> (name, @name(c.ρ)) => similar(Y.c, TridiagonalRow),
+                diffused_scalar_names,
+            )...,
+            MatrixFields.unrolled_map(
+                name -> (name, name) => similar(Y.c, TridiagonalRow),
+                diffused_scalar_names,
+            )...,
+            (
+                is_in_Y(@name(c.ρq_tot)) ?
+                (
+                    (@name(c.ρe_tot), @name(c.ρq_tot)) =>
+                        similar(Y.c, TridiagonalRow),
+                ) : ()
+            )...,
+            (@name(c.uₕ), @name(c.uₕ)) =>
+                !isnothing(p.atmos.turbconv_model) ||
+                    diffuse_momentum(p.atmos.vert_diff) ?
+                similar(Y.c, TridiagonalRow) : FT(-1) * I,
+        )
+    else
+        MatrixFields.unrolled_map(
+            name -> (name, name) => FT(-1) * I,
+            (diffused_scalar_names..., @name(c.uₕ)),
+        )
+    end
+
+    sgs_advection_blocks = if p.atmos.turbconv_model isa PrognosticEDMFX
+        @assert n_prognostic_mass_flux_subdomains(p.atmos.turbconv_model) == 1
+        sgs_scalar_names = (
+            @name(c.sgsʲs.:(1).q_tot),
+            @name(c.sgsʲs.:(1).mse),
+            @name(c.sgsʲs.:(1).ρa)
+        )
+        if use_derivative(sgs_advection_flag)
+            (
+                MatrixFields.unrolled_map(
+                    name -> (name, name) => similar(Y.c, TridiagonalRow),
+                    sgs_scalar_names,
+                )...,
+                (@name(c.sgsʲs.:(1).mse), @name(c.ρ)) =>
+                    similar(Y.c, DiagonalRow),
+                (@name(c.sgsʲs.:(1).mse), @name(c.sgsʲs.:(1).q_tot)) =>
+                    similar(Y.c, DiagonalRow),
+                (@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).q_tot)) =>
+                    similar(Y.c, TridiagonalRow),
+                (@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).mse)) =>
+                    similar(Y.c, TridiagonalRow),
+                (@name(f.sgsʲs.:(1).u₃), @name(c.ρ)) =>
+                    similar(Y.f, BidiagonalRow_C3),
+                (@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).q_tot)) =>
+                    similar(Y.f, BidiagonalRow_C3),
+                (@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).mse)) =>
+                    similar(Y.f, BidiagonalRow_C3),
+                (@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)) =>
+                    similar(Y.f, TridiagonalRow_C3xACT3),
+            )
+        else
+            (
+                MatrixFields.unrolled_map(
+                    name -> (name, name) => FT(-1) * I,
+                    sgs_scalar_names,
+                )...,
+                (@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)) =>
+                    !isnothing(p.atmos.rayleigh_sponge) ?
+                    similar(Y.f, DiagonalRow_C3xACT3) : FT(-1) * I,
+            )
+        end
+    else
+        ()
+    end
+
+    matrix = MatrixFields.FieldMatrix(
+        identity_blocks...,
+        sgs_advection_blocks...,
+        advection_blocks...,
+        diffusion_blocks...,
+    )
+
+    sgs_names_if_available = if p.atmos.turbconv_model isa PrognosticEDMFX
+        (
+            @name(c.sgsʲs.:(1).q_tot),
+            @name(c.sgsʲs.:(1).mse),
+            @name(c.sgsʲs.:(1).ρa),
+            @name(f.sgsʲs.:(1).u₃),
+        )
+    else
+        ()
+    end
+    names₁_group₁ = (@name(c.ρ), sfc_if_available...)
+    names₁_group₂ = (available_tracer_names..., ρatke_if_available...)
+    names₁_group₃ = (@name(c.ρe_tot),)
+    names₁ = (
+        names₁_group₁...,
+        sgs_names_if_available...,
+        names₁_group₂...,
+        names₁_group₃...,
+    )
+
+    alg₂ = MatrixFields.BlockLowerTriangularSolve(@name(c.uₕ))
+    alg =
+        if use_derivative(diffusion_flag) || use_derivative(sgs_advection_flag)
+            alg₁_subalg₂ =
+                if p.atmos.turbconv_model isa PrognosticEDMFX &&
+                   use_derivative(sgs_advection_flag)
+                    diff_subalg =
+                        use_derivative(diffusion_flag) ?
+                        (;
+                            alg₂ = MatrixFields.BlockLowerTriangularSolve(
+                                names₁_group₂...,
+                            )
+                        ) : (;)
+                    (;
+                        alg₂ = MatrixFields.BlockLowerTriangularSolve(
+                            @name(c.sgsʲs.:(1).q_tot);
+                            alg₂ = MatrixFields.BlockLowerTriangularSolve(
+                                @name(c.sgsʲs.:(1).mse);
+                                alg₂ = MatrixFields.BlockLowerTriangularSolve(
+                                    @name(c.sgsʲs.:(1).ρa),
+                                    @name(f.sgsʲs.:(1).u₃);
+                                    diff_subalg...,
+                                ),
+                            ),
+                        )
+                    )
+                else
+                    is_in_Y(@name(c.ρq_tot)) ?
+                    (;
+                        alg₂ = MatrixFields.BlockLowerTriangularSolve(
+                            names₁_group₂...,
+                        )
+                    ) : (;)
+                end
+            alg₁ = MatrixFields.BlockLowerTriangularSolve(
+                names₁_group₁...;
+                alg₁_subalg₂...,
+            )
+            MatrixFields.ApproximateBlockArrowheadIterativeSolve(
+                names₁...;
+                alg₁,
+                alg₂,
+                P_alg₁ = MatrixFields.MainDiagonalPreconditioner(),
+                n_iters = approximate_solve_iters,
+            )
+        else
+            MatrixFields.BlockArrowheadSolve(names₁...; alg₂)
+        end
+
+    solver = MatrixFields.FieldMatrixSolver(alg, matrix, Y)
+
+    return (;
+        diffusion_flag,
+        topography_flag,
+        sgs_advection_flag,
+        matrix,
+        solver,
+    )
+end
+
+always_update_exact_jacobian(::ApproxJacobian) = false
+
+factorize_exact_jacobian!(::ApproxJacobian, _, _, _, _, _) = nothing
+
+function approximate_jacobian!(::ApproxJacobian, cache, Y, p, dtγ, t)
+    (; diffusion_flag, topography_flag, sgs_advection_flag, matrix) = cache
+    (; params) = p
+    (; edmfx_upwinding) = p.atmos.numerics
+    (; ᶜΦ, ᶠgradᵥ_ᶜΦ) = p.core
+    (; ᶜspecific, ᶜK, ᶜts, ᶜp) = p.precomputed
+    (;
+        ∂ᶜK_∂ᶜuₕ,
+        ∂ᶜK_∂ᶠu₃,
+        ᶠp_grad_matrix,
+        ᶜadvection_matrix,
+        ᶜdiffusion_h_matrix,
+        ᶜdiffusion_u_matrix,
+        ᶠbidiagonal_matrix_ct3,
+        ᶠbidiagonal_matrix_ct3_2,
+        ᶠtridiagonal_matrix_c3,
+    ) = p.scratch
+
+    FT = Spaces.undertype(axes(Y.c))
+    CTh = CTh_vector_type(axes(Y.c))
+    one_C3xACT3 = C3(FT(1)) * CT3(FT(1))'
+
+    cv_d = FT(CAP.cv_d(params))
+    Δcv_v = FT(CAP.cv_v(params)) - cv_d
+    T_tri = FT(CAP.T_triple(params))
+    e_int_v0 = T_tri * Δcv_v - FT(CAP.e_int_v0(params))
+    thermo_params = CAP.thermodynamics_params(params)
+
+    ᶜρ = Y.c.ρ
+    ᶜuₕ = Y.c.uₕ
+    ᶠu₃ = Y.f.u₃
+    ᶜJ = Fields.local_geometry_field(Y.c).J
+    ᶜgⁱʲ = Fields.local_geometry_field(Y.c).gⁱʲ
+    ᶠgⁱʲ = Fields.local_geometry_field(Y.f).gⁱʲ
+
+    ᶜkappa_m = p.scratch.ᶜtemp_scalar
+    @. ᶜkappa_m =
+        TD.gas_constant_air(thermo_params, ᶜts) / TD.cv_m(thermo_params, ᶜts)
+
+    if use_derivative(topography_flag)
+        @. ∂ᶜK_∂ᶜuₕ = DiagonalMatrixRow(
+            adjoint(CTh(ᶜuₕ)) + adjoint(ᶜinterp(ᶠu₃)) * g³ʰ(ᶜgⁱʲ),
+        )
+    else
+        @. ∂ᶜK_∂ᶜuₕ = DiagonalMatrixRow(adjoint(CTh(ᶜuₕ)))
+    end
+    @. ∂ᶜK_∂ᶠu₃ =
+        ᶜinterp_matrix() ⋅ DiagonalMatrixRow(adjoint(CT3(ᶠu₃))) +
+        DiagonalMatrixRow(adjoint(CT3(ᶜuₕ))) ⋅ ᶜinterp_matrix()
+
+    @. ᶠp_grad_matrix = DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ)) ⋅ ᶠgradᵥ_matrix()
+
+    @. ᶜadvection_matrix =
+        -(ᶜadvdivᵥ_matrix()) ⋅ DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρ))
+
+    if use_derivative(topography_flag)
+        ∂ᶜρ_err_∂ᶜuₕ = matrix[@name(c.ρ), @name(c.uₕ)]
+        @. ∂ᶜρ_err_∂ᶜuₕ =
+            dtγ * ᶜadvection_matrix ⋅ ᶠwinterp_matrix(ᶜJ * ᶜρ) ⋅
+            DiagonalMatrixRow(g³ʰ(ᶜgⁱʲ))
+    end
+    ∂ᶜρ_err_∂ᶠu₃ = matrix[@name(c.ρ), @name(f.u₃)]
+    @. ∂ᶜρ_err_∂ᶠu₃ = dtγ * ᶜadvection_matrix ⋅ DiagonalMatrixRow(g³³(ᶠgⁱʲ))
+
+    tracer_info = (
+        (@name(c.ρe_tot), @name(ᶜh_tot)),
+        (@name(c.ρq_tot), @name(ᶜspecific.q_tot)),
+    )
+    MatrixFields.unrolled_foreach(tracer_info) do (ρχ_name, χ_name)
+        MatrixFields.has_field(Y, ρχ_name) || return
+        ᶜχ = MatrixFields.get_field(p.precomputed, χ_name)
+        if use_derivative(topography_flag)
+            ∂ᶜρχ_err_∂ᶜuₕ = matrix[ρχ_name, @name(c.uₕ)]
+        end
+        ∂ᶜρχ_err_∂ᶠu₃ = matrix[ρχ_name, @name(f.u₃)]
+        use_derivative(topography_flag) && @. ∂ᶜρχ_err_∂ᶜuₕ =
+            dtγ * ᶜadvection_matrix ⋅ DiagonalMatrixRow(ᶠinterp(ᶜχ)) ⋅
+            ᶠwinterp_matrix(ᶜJ * ᶜρ) ⋅ DiagonalMatrixRow(g³ʰ(ᶜgⁱʲ))
+        @. ∂ᶜρχ_err_∂ᶠu₃ =
+            dtγ * ᶜadvection_matrix ⋅ DiagonalMatrixRow(ᶠinterp(ᶜχ) * g³³(ᶠgⁱʲ))
+    end
+
+    ∂ᶠu₃_err_∂ᶜρ = matrix[@name(f.u₃), @name(c.ρ)]
+    ∂ᶠu₃_err_∂ᶜρe_tot = matrix[@name(f.u₃), @name(c.ρe_tot)]
+    @. ∂ᶠu₃_err_∂ᶜρ =
+        dtγ * (
+            ᶠp_grad_matrix ⋅
+            DiagonalMatrixRow(ᶜkappa_m * (T_tri * cv_d - ᶜK - ᶜΦ)) +
+            DiagonalMatrixRow(ᶠgradᵥ(ᶜp) / abs2(ᶠinterp(ᶜρ))) ⋅
+            ᶠinterp_matrix()
+        )
+    @. ∂ᶠu₃_err_∂ᶜρe_tot = dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow(ᶜkappa_m)
+    if MatrixFields.has_field(Y, @name(c.ρq_tot))
+        ∂ᶠu₃_err_∂ᶜρq_tot = matrix[@name(f.u₃), @name(c.ρq_tot)]
+        @. ∂ᶠu₃_err_∂ᶜρq_tot =
+            dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow(ᶜkappa_m * e_int_v0)
+    end
+
+    ∂ᶠu₃_err_∂ᶜuₕ = matrix[@name(f.u₃), @name(c.uₕ)]
+    ∂ᶠu₃_err_∂ᶠu₃ = matrix[@name(f.u₃), @name(f.u₃)]
+    I_u₃ = DiagonalMatrixRow(one_C3xACT3)
+    @. ∂ᶠu₃_err_∂ᶜuₕ =
+        dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow(-(ᶜkappa_m) * ᶜρ) ⋅ ∂ᶜK_∂ᶜuₕ
+    if p.atmos.rayleigh_sponge isa RayleighSponge
+        @. ∂ᶠu₃_err_∂ᶠu₃ =
+            dtγ * (
+                ᶠp_grad_matrix ⋅ DiagonalMatrixRow(-(ᶜkappa_m) * ᶜρ) ⋅
+                ∂ᶜK_∂ᶠu₃ + DiagonalMatrixRow(-p.ᶠβ_rayleigh_w * (one_C3xACT3,))
+            ) - (I_u₃,)
+    else
+        @. ∂ᶠu₃_err_∂ᶠu₃ =
+            dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow(-(ᶜkappa_m) * ᶜρ) ⋅
+            ∂ᶜK_∂ᶠu₃ - (I_u₃,)
+    end
+
+    if use_derivative(diffusion_flag)
+        (; ᶜK_h, ᶜK_u) = p
+        @. ᶜdiffusion_h_matrix =
+            ᶜadvdivᵥ_matrix() ⋅ DiagonalMatrixRow(ᶠinterp(ᶜρ) * ᶠinterp(ᶜK_h)) ⋅
+            ᶠgradᵥ_matrix()
+        if (
+            MatrixFields.has_field(Y, @name(c.sgs⁰.ρatke)) ||
+            !isnothing(p.atmos.turbconv_model) ||
+            diffuse_momentum(p.atmos.vert_diff)
+        )
+            @. ᶜdiffusion_u_matrix =
+                ᶜadvdivᵥ_matrix() ⋅
+                DiagonalMatrixRow(ᶠinterp(ᶜρ) * ᶠinterp(ᶜK_u)) ⋅ ᶠgradᵥ_matrix()
+        end
+
+        ∂ᶜρe_tot_err_∂ᶜρ = matrix[@name(c.ρe_tot), @name(c.ρ)]
+        ∂ᶜρe_tot_err_∂ᶜρe_tot = matrix[@name(c.ρe_tot), @name(c.ρe_tot)]
+        @. ∂ᶜρe_tot_err_∂ᶜρ =
+            dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow(
+                (
+                    -(1 + ᶜkappa_m) * ᶜspecific.e_tot -
+                    ᶜkappa_m * e_int_v0 * ᶜspecific.q_tot
+                ) / ᶜρ,
+            )
+        @. ∂ᶜρe_tot_err_∂ᶜρe_tot =
+            dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow((1 + ᶜkappa_m) / ᶜρ) -
+            (I,)
+        if MatrixFields.has_field(Y, @name(c.ρq_tot))
+            ∂ᶜρe_tot_err_∂ᶜρq_tot = matrix[@name(c.ρe_tot), @name(c.ρq_tot)]
+            @. ∂ᶜρe_tot_err_∂ᶜρq_tot =
+                dtγ * ᶜdiffusion_h_matrix ⋅
+                DiagonalMatrixRow(ᶜkappa_m * e_int_v0 / ᶜρ)
+        end
+
+        tracer_info = (
+            (@name(c.ρq_tot), @name(q_tot)),
+            (@name(c.ρq_liq), @name(q_liq)),
+            (@name(c.ρq_ice), @name(q_ice)),
+            (@name(c.ρq_rai), @name(q_rai)),
+            (@name(c.ρq_sno), @name(q_sno)),
+        )
+        MatrixFields.unrolled_foreach(tracer_info) do (ρq_name, q_name)
+            MatrixFields.has_field(Y, ρq_name) || return
+            ᶜq = MatrixFields.get_field(ᶜspecific, q_name)
+            ∂ᶜρq_err_∂ᶜρ = matrix[ρq_name, @name(c.ρ)]
+            ∂ᶜρq_err_∂ᶜρq = matrix[ρq_name, ρq_name]
+            @. ∂ᶜρq_err_∂ᶜρ =
+                dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow(-(ᶜq) / ᶜρ)
+            @. ∂ᶜρq_err_∂ᶜρq =
+                dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow(1 / ᶜρ) - (I,)
+        end
+
+        if MatrixFields.has_field(Y, @name(c.sgs⁰.ρatke))
+            turbconv_params = CAP.turbconv_params(params)
+            c_d = CAP.tke_diss_coeff(turbconv_params)
+            (; ᶜtke⁰, ᶜmixing_length, dt) = p
+            ᶜρa⁰ = p.atmos.turbconv_model isa PrognosticEDMFX ? p.ᶜρa⁰ : ᶜρ
+            ᶜρatke⁰ = Y.c.sgs⁰.ρatke
+
+            @inline dissipation_rate(tke⁰, mixing_length) =
+                tke⁰ >= 0 ? c_d * sqrt(tke⁰) / max(mixing_length, 1) : 1 / dt
+            @inline ∂dissipation_rate_∂tke⁰(tke⁰, mixing_length) =
+                tke⁰ > 0 ? c_d / (2 * max(mixing_length, 1) * sqrt(tke⁰)) :
+                typeof(tke⁰)(0)
+
+            ᶜdissipation_matrix_diagonal = p.scratch.ᶜtemp_scalar
+            @. ᶜdissipation_matrix_diagonal =
+                ᶜρatke⁰ * ∂dissipation_rate_∂tke⁰(ᶜtke⁰, ᶜmixing_length)
+
+            ∂ᶜρatke⁰_err_∂ᶜρ = matrix[@name(c.sgs⁰.ρatke), @name(c.ρ)]
+            ∂ᶜρatke⁰_err_∂ᶜρatke⁰ =
+                matrix[@name(c.sgs⁰.ρatke), @name(c.sgs⁰.ρatke)]
+            @. ∂ᶜρatke⁰_err_∂ᶜρ =
+                dtγ * (
+                    ᶜdiffusion_u_matrix -
+                    DiagonalMatrixRow(ᶜdissipation_matrix_diagonal)
+                ) ⋅ DiagonalMatrixRow(-(ᶜtke⁰) / ᶜρa⁰)
+            @. ∂ᶜρatke⁰_err_∂ᶜρatke⁰ =
+                dtγ * (
+                    (
+                        ᶜdiffusion_u_matrix -
+                        DiagonalMatrixRow(ᶜdissipation_matrix_diagonal)
+                    ) ⋅ DiagonalMatrixRow(1 / ᶜρa⁰) -
+                    DiagonalMatrixRow(dissipation_rate(ᶜtke⁰, ᶜmixing_length))
+                ) - (I,)
+        end
+
+        if (
+            !isnothing(p.atmos.turbconv_model) ||
+            diffuse_momentum(p.atmos.vert_diff)
+        )
+            ∂ᶜuₕ_err_∂ᶜuₕ = matrix[@name(c.uₕ), @name(c.uₕ)]
+            @. ∂ᶜuₕ_err_∂ᶜuₕ =
+                dtγ * DiagonalMatrixRow(1 / ᶜρ) ⋅ ᶜdiffusion_u_matrix - (I,)
+        end
+
+        ᶠlg = Fields.local_geometry_field(Y.f)
+        precip_info =
+            ((@name(c.ρq_rai), @name(ᶜwᵣ)), (@name(c.ρq_sno), @name(ᶜwₛ)))
+        MatrixFields.unrolled_foreach(precip_info) do (ρqₚ_name, wₚ_name)
+            MatrixFields.has_field(Y, ρqₚ_name) || return
+            ∂ᶜρqₚ_err_∂ᶜρqₚ = matrix[ρqₚ_name, ρqₚ_name]
+            ᶜwₚ = MatrixFields.get_field(p.precomputed, wₚ_name)
+            ᶠtmp = p.scratch.ᶠtemp_CT3
+            @. ᶠtmp = CT3(unit_basis_vector_data(CT3, ᶠlg)) * ᶠwinterp(ᶜJ, ᶜρ)
+            @. ∂ᶜρqₚ_err_∂ᶜρqₚ +=
+                dtγ * -(ᶜprecipdivᵥ_matrix()) ⋅ DiagonalMatrixRow(ᶠtmp) ⋅
+                ᶠright_bias_matrix() ⋅ DiagonalMatrixRow(-(ᶜwₚ) / ᶜρ)
+        end
+    end
+
+    if p.atmos.turbconv_model isa PrognosticEDMFX
+        if use_derivative(sgs_advection_flag)
+            (; ᶜgradᵥ_ᶠΦ, ᶜρʲs, ᶠu³ʲs, ᶜtsʲs) = p
+            (; bdmr_l, bdmr_r, bdmr) = p
+            is_third_order = edmfx_upwinding == Val(:third_order)
+            ᶠupwind = is_third_order ? ᶠupwind3 : ᶠupwind1
+            ᶠset_upwind_bcs = Operators.SetBoundaryOperator(;
+                top = Operators.SetValue(zero(CT3{FT})),
+                bottom = Operators.SetValue(zero(CT3{FT})),
+            ) # Need to wrap ᶠupwind in this for well-defined boundaries.
+            UpwindMatrixRowType =
+                is_third_order ? QuaddiagonalMatrixRow : BidiagonalMatrixRow
+            ᶠupwind_matrix = is_third_order ? ᶠupwind3_matrix : ᶠupwind1_matrix
+            ᶠset_upwind_matrix_bcs = Operators.SetBoundaryOperator(;
+                top = Operators.SetValue(zero(UpwindMatrixRowType{CT3{FT}})),
+                bottom = Operators.SetValue(zero(UpwindMatrixRowType{CT3{FT}})),
+            ) # Need to wrap ᶠupwind_matrix in this for well-defined boundaries.
+            ᶜkappa_mʲ = p.scratch.ᶜtemp_scalar
+            @. ᶜkappa_mʲ =
+                TD.gas_constant_air(thermo_params, ᶜtsʲs.:(1)) /
+                TD.cv_m(thermo_params, ᶜtsʲs.:(1))
+            ∂ᶜq_totʲ_err_∂ᶜq_totʲ =
+                matrix[@name(c.sgsʲs.:(1).q_tot), @name(c.sgsʲs.:(1).q_tot)]
+            @. ∂ᶜq_totʲ_err_∂ᶜq_totʲ =
+                dtγ * (
+                    DiagonalMatrixRow(ᶜadvdivᵥ(ᶠu³ʲs.:(1))) -
+                    ᶜadvdivᵥ_matrix() ⋅
+                    ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1)))
+                ) - (I,)
+
+            ∂ᶜmseʲ_err_∂ᶜq_totʲ =
+                matrix[@name(c.sgsʲs.:(1).mse), @name(c.sgsʲs.:(1).q_tot)]
+            @. ∂ᶜmseʲ_err_∂ᶜq_totʲ =
+                dtγ * (
+                    -DiagonalMatrixRow(
+                        adjoint(ᶜinterp(ᶠu³ʲs.:(1))) *
+                        ᶜgradᵥ_ᶠΦ *
+                        Y.c.ρ *
+                        ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp) * e_int_v0,
+                    )
+                )
+            ∂ᶜmseʲ_err_∂ᶜρ = matrix[@name(c.sgsʲs.:(1).mse), @name(c.ρ)]
+            @. ∂ᶜmseʲ_err_∂ᶜρ =
+                dtγ * (
+                    -DiagonalMatrixRow(
+                        adjoint(ᶜinterp(ᶠu³ʲs.:(1))) * ᶜgradᵥ_ᶠΦ / ᶜρʲs.:(1),
+                    )
+                )
+            ∂ᶜmseʲ_err_∂ᶜmseʲ =
+                matrix[@name(c.sgsʲs.:(1).mse), @name(c.sgsʲs.:(1).mse)]
+            @. ∂ᶜmseʲ_err_∂ᶜmseʲ =
+                dtγ * (
+                    DiagonalMatrixRow(ᶜadvdivᵥ(ᶠu³ʲs.:(1))) -
+                    ᶜadvdivᵥ_matrix() ⋅
+                    ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) -
+                    DiagonalMatrixRow(
+                        adjoint(ᶜinterp(ᶠu³ʲs.:(1))) *
+                        ᶜgradᵥ_ᶠΦ *
+                        Y.c.ρ *
+                        ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp),
+                    )
+                ) - (I,)
+
+            ∂ᶜρaʲ_err_∂ᶜq_totʲ =
+                matrix[@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).q_tot)]
+            @. ᶠbidiagonal_matrix_ct3 =
+                DiagonalMatrixRow(
+                    ᶠset_upwind_bcs(
+                        ᶠupwind(
+                            ᶠu³ʲs.:(1),
+                            draft_area(Y.c.sgsʲs.:(1).ρa, ᶜρʲs.:(1)),
+                        ),
+                    ),
+                ) ⋅ ᶠwinterp_matrix(ᶜJ) ⋅ DiagonalMatrixRow(
+                    ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp) *
+                    e_int_v0,
+                )
+            @. ᶠbidiagonal_matrix_ct3_2 =
+                DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρʲs.:(1))) ⋅
+                ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) ⋅
+                DiagonalMatrixRow(
+                    Y.c.sgsʲs.:(1).ρa * ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp) *
+                    e_int_v0,
+                )
+
+            @. ∂ᶜρaʲ_err_∂ᶜq_totʲ =
+                dtγ * ᶜadvdivᵥ_matrix() ⋅
+                (ᶠbidiagonal_matrix_ct3 - ᶠbidiagonal_matrix_ct3_2)
+            ∂ᶜρaʲ_err_∂ᶜmseʲ =
+                matrix[@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).mse)]
+            @. ᶠbidiagonal_matrix_ct3 =
+                DiagonalMatrixRow(
+                    ᶠset_upwind_bcs(
+                        ᶠupwind(
+                            ᶠu³ʲs.:(1),
+                            draft_area(Y.c.sgsʲs.:(1).ρa, ᶜρʲs.:(1)),
+                        ),
+                    ),
+                ) ⋅ ᶠwinterp_matrix(ᶜJ) ⋅ DiagonalMatrixRow(
+                    ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp),
+                )
+            @. ᶠbidiagonal_matrix_ct3_2 =
+                DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρʲs.:(1))) ⋅
+                ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) ⋅
+                DiagonalMatrixRow(
+                    Y.c.sgsʲs.:(1).ρa * ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp),
+                )
+            @. ∂ᶜρaʲ_err_∂ᶜmseʲ =
+                dtγ * ᶜadvdivᵥ_matrix() ⋅
+                (ᶠbidiagonal_matrix_ct3 - ᶠbidiagonal_matrix_ct3_2)
+            ∂ᶜρaʲ_err_∂ᶜρaʲ =
+                matrix[@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).ρa)]
+            @. ᶜadvection_matrix =
+                -(ᶜadvdivᵥ_matrix()) ⋅
+                DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρʲs.:(1)))
+            @. ∂ᶜρaʲ_err_∂ᶜρaʲ =
+                dtγ * ᶜadvection_matrix ⋅
+                ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) ⋅
+                DiagonalMatrixRow(1 / ᶜρʲs.:(1)) - (I,)
+
+            ∂ᶠu₃ʲ_err_∂ᶜρ = matrix[@name(f.sgsʲs.:(1).u₃), @name(c.ρ)]
+            @. ∂ᶠu₃ʲ_err_∂ᶜρ =
+                dtγ * DiagonalMatrixRow(ᶠgradᵥ_ᶜΦ / ᶠinterp(ᶜρʲs.:(1))) ⋅
+                ᶠinterp_matrix()
+            ∂ᶠu₃ʲ_err_∂ᶜq_totʲ =
+                matrix[@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).q_tot)]
+            @. ∂ᶠu₃ʲ_err_∂ᶜq_totʲ =
+                dtγ * DiagonalMatrixRow(
+                    ᶠgradᵥ_ᶜΦ * ᶠinterp(Y.c.ρ) / (ᶠinterp(ᶜρʲs.:(1)))^2,
+                ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow(
+                    ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp) *
+                    e_int_v0,
+                )
+            ∂ᶠu₃ʲ_err_∂ᶜmseʲ =
+                matrix[@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).mse)]
+            @. ∂ᶠu₃ʲ_err_∂ᶜmseʲ =
+                dtγ * DiagonalMatrixRow(
+                    ᶠgradᵥ_ᶜΦ * ᶠinterp(Y.c.ρ) / (ᶠinterp(ᶜρʲs.:(1)))^2,
+                ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow(
+                    ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp),
+                )
+            ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ =
+                matrix[@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)]
+            ᶜu₃ʲ = p.scratch.ᶜtemp_C3
+            @. ᶜu₃ʲ = ᶜinterp(Y.f.sgsʲs.:(1).u₃)
+
+            @. bdmr_l = convert(BidiagonalMatrixRow{FT}, ᶜleft_bias_matrix())
+            @. bdmr_r = convert(BidiagonalMatrixRow{FT}, ᶜright_bias_matrix())
+            @. bdmr = ifelse(ᶜu₃ʲ.components.data.:1 > 0, bdmr_l, bdmr_r)
+
+            @. ᶠtridiagonal_matrix_c3 = -(ᶠgradᵥ_matrix()) ⋅ bdmr
+            if p.atmos.rayleigh_sponge isa RayleighSponge
+                @. ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ =
+                    dtγ * (
+                        ᶠtridiagonal_matrix_c3 ⋅
+                        DiagonalMatrixRow(adjoint(CT3(Y.f.sgsʲs.:(1).u₃))) -
+                        DiagonalMatrixRow(p.ᶠβ_rayleigh_w * (one_C3xACT3,))
+                    ) - (I_u₃,)
+            else
+                @. ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ =
+                    dtγ * ᶠtridiagonal_matrix_c3 ⋅
+                    DiagonalMatrixRow(adjoint(CT3(Y.f.sgsʲs.:(1).u₃))) - (I_u₃,)
+            end
+        elseif p.atmos.rayleigh_sponge isa RayleighSponge
+            ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ =
+                matrix[@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)]
+            @. ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ =
+                dtγ * -DiagonalMatrixRow(p.ᶠβ_rayleigh_w * (one_C3xACT3,)) -
+                (I_u₃,)
+        end
+    end
+end
+
+invert_jacobian!(::ApproxJacobian, cache, x, b) =
+    MatrixFields.field_matrix_solve!(cache.solver, x, cache.matrix, b)
diff --git a/src/prognostic_equations/implicit/debug_jacobian.jl b/src/prognostic_equations/implicit/debug_jacobian.jl
new file mode 100644
index 00000000000..4be4a130189
--- /dev/null
+++ b/src/prognostic_equations/implicit/debug_jacobian.jl
@@ -0,0 +1,108 @@
+struct DebugJacobian{E <: ExactJacobian, A <: ApproxJacobian} <:
+       JacobianAlgorithm
+    exact_jacobian_mode::E
+    approx_jacobian_mode::A
+    use_exact_jacobian::Bool
+end
+function DebugJacobian(
+    approx_jacobian_mode;
+    use_exact_jacobian = true,
+    exact_jacobian_mode_kwargs...,
+)
+    exact_jacobian_mode = ExactJacobian(;
+        preserve_unfactorized_jacobian = true,
+        exact_jacobian_mode_kwargs...,
+    )
+    return DebugJacobian(
+        exact_jacobian_mode,
+        approx_jacobian_mode,
+        use_exact_jacobian,
+    )
+end
+
+jacobian_cache(alg::DebugJacobian, Y, p) = (;
+    exact_jacobian_cache = jacobian_cache(alg.exact_jacobian_mode, Y, p),
+    approx_jacobian_cache = jacobian_cache(alg.approx_jacobian_mode, Y, p),
+    similar_to_x = similar(Y),
+    x_error = similar(Y),
+    x_rel_error = similar(Y),
+    compute_error_ref = Ref{Bool}(),
+    t_ref = Ref{typeof(p.t_end)}(),
+)
+
+always_update_exact_jacobian(alg::DebugJacobian) =
+    always_update_exact_jacobian(alg.exact_jacobian_mode)
+
+function factorize_exact_jacobian!(alg::DebugJacobian, cache, Y, p, dtγ, t)
+    factorize_exact_jacobian!(
+        alg.exact_jacobian_mode,
+        cache.exact_jacobian_cache,
+        Y,
+        p,
+        dtγ,
+        t,
+    )
+    cache.compute_error_ref[] = true
+    cache.t_ref[] = t
+end
+
+approximate_jacobian!(alg::DebugJacobian, cache, Y, p, dtγ, t) =
+    (!alg.use_exact_jacobian || cache.compute_error_ref[]) &&
+    approximate_jacobian!(
+        alg.approx_jacobian_mode,
+        cache.approx_jacobian_cache,
+        Y,
+        p,
+        dtγ,
+        t,
+    )
+
+import Dates: @sprintf # TODO: Add Printf as a dependency of ClimaAtmos.
+
+function invert_jacobian!(alg::DebugJacobian, cache, x, b)
+    (; exact_jacobian_mode, approx_jacobian_mode, use_exact_jacobian) = alg
+    (; exact_jacobian_cache, approx_jacobian_cache) = cache
+    (; similar_to_x, x_error, x_rel_error, compute_error_ref, t_ref) = cache
+    if compute_error_ref[]
+        x_exact = use_exact_jacobian ? x : similar_to_x
+        x_approx = use_exact_jacobian ? similar_to_x : x
+        x_exact .= x_approx .= NaN # TODO: Remove this sanity check
+        invert_jacobian!(exact_jacobian_mode, exact_jacobian_cache, x_exact, b)
+        invert_jacobian!(
+            approx_jacobian_mode,
+            approx_jacobian_cache,
+            x_approx,
+            b,
+        )
+        compute_error_ref[] = false
+
+        @. x_error = abs(x_approx - x_exact)
+        @. x_rel_error = ifelse(
+            x_approx == x_exact == 0,
+            0,
+            abs((x_approx - x_exact) / x_exact),
+        )
+
+        rms(f, x) = sqrt(mean(abs2 ∘ f, x))
+        rms(x) = rms(identity, x)
+
+        field_vector_rms(x) = rms(rms ∘ Fields.backing_array, Fields._values(x))
+        field_vector_max(x) =
+            maximum(maximum ∘ Fields.backing_array, Fields._values(x))
+
+        rel_rms_error = field_vector_rms(x_error) / field_vector_rms(x_exact)
+        rel_max_error = field_vector_max(x_error) / field_vector_max(x_exact)
+        rms_rel_error = field_vector_rms(x_rel_error)
+        max_rel_error = field_vector_max(x_rel_error)
+        @info "Error of approximate implicit solve at t = $(t_ref[]):\n  \
+               relative RMS error = $(@sprintf("%.3e", rel_rms_error))\n  \
+               relative max error = $(@sprintf("%.3e", rel_max_error))\n  \
+               RMS relative error = $(@sprintf("%.3e", rms_rel_error))\n  \
+               max relative error = $(@sprintf("%.3e", max_rel_error))"
+        @assert rel_rms_error < 1 # TODO: Make this test more rigorous
+    elseif use_exact_jacobian
+        invert_jacobian!(exact_jacobian_mode, exact_jacobian_cache, x, b)
+    else
+        invert_jacobian!(approx_jacobian_mode, approx_jacobian_cache, x, b)
+    end
+end
diff --git a/src/prognostic_equations/implicit/dual_fixes.jl b/src/prognostic_equations/implicit/dual_fixes.jl
new file mode 100644
index 00000000000..2dffd032892
--- /dev/null
+++ b/src/prognostic_equations/implicit/dual_fixes.jl
@@ -0,0 +1,182 @@
+# TODO: Move these to Thermodynamics.jl
+TD.air_density(param_set, T, p, q) =
+    TD.air_density(param_set, TD.promote_phase_partition(T, p, q)...)
+TD.air_pressure(param_set, T, ρ, q) =
+    TD.air_pressure(param_set, TD.promote_phase_partition(T, ρ, q)...)
+TD.internal_energy(param_set, T, q) =
+    TD.internal_energy(param_set, TD.promote_phase_partition(T, T, q)[2:end]...)
+TD.internal_energy_sat(param_set, T, ρ, q_tot, type) =
+    TD.internal_energy_sat(param_set, promote(T, ρ, q_tot)..., type)
+TD.PhasePartition_equil(param_set, T, ρ, q_tot, type) =
+    TD.PhasePartition_equil(param_set, promote(T, ρ, q_tot)..., type)
+TD.PhasePartition_equil_given_p(param_set, T, p, q_tot, type) =
+    TD.PhasePartition_equil_given_p(param_set, promote(T, p, q_tot)..., type)
+TD.specific_enthalpy_sat(param_set, T, ρ, q_tot, type) =
+    TD.specific_enthalpy_sat(param_set, promote(T, ρ, q_tot)..., type)
+TD.saturation_vapor_pressure(param_set, T, LH_0, Δcp) =
+    TD.saturation_vapor_pressure(param_set, promote(T, LH_0, Δcp)...)
+
+# TODO: Move these to SurfaceFluxes.jl
+import SurfaceFluxes as SF
+SF.FluxesAndFrictionVelocity(int_values, surf_values, args...) =
+    SF.FluxesAndFrictionVelocity(int_values, surf_values, promote(args...)...)
+
+# TODO: Move this to ClimaCore.jl
+ClimaComms.device(fv::Fields.FieldVector) =
+    ClimaComms.device(first(Fields._values(fv)))
+
+# Ensure that we can broadcast between FieldVectors and other BlockVectors, but
+# without going through the type-unstable BlockArrays broadcasting code.
+# TODO: Replace the FieldVector broadcasting code in ClimaCore.jl with this
+struct NewFieldVectorStyle <: Base.Broadcast.AbstractArrayStyle{1} end
+Base.Broadcast.BroadcastStyle(::Type{<:Fields.FieldVector}) =
+    NewFieldVectorStyle()
+Base.Broadcast.BroadcastStyle(
+    fs::NewFieldVectorStyle,
+    as::Base.Broadcast.DefaultArrayStyle{N},
+) where {N} = N in (0, 1) ? fs : as # prevents an ambiguity with Base.Broadcast
+Base.Broadcast.BroadcastStyle(
+    fs::NewFieldVectorStyle,
+    as::Base.Broadcast.AbstractArrayStyle{N},
+) where {N} = N in (0, 1) ? fs : as
+Base.@propagate_inbounds function BlockArrays.viewblock(
+    fv::Fields.FieldVector,
+    block::BlockArrays.Block{1},
+)
+    array = Fields.backing_array(Fields._values(fv)[block.n...])
+    # vec(array) allocates; see https://github.com/JuliaLang/julia/issues/36313
+    return Base.ReshapedArray(array, (length(array),), ())
+end
+
+first_field_vector(x::Fields.FieldVector, args...) = x
+@inline first_field_vector(x::Base.Broadcast.Broadcasted, args...) =
+    first_field_vector(x.args..., args...)
+@inline first_field_vector(_, args...) = first_field_vector(args...)
+check_broadcast_axes(fv, x::Base.Broadcast.Broadcasted) =
+    check_broadcast_args_axes(fv, x.args...)
+@inline check_broadcast_axes(fv, x) =
+    isempty(axes(x)) || # scalar 
+    axes(x) == axes(fv) || # AbstractBlockVector
+    x isa AbstractVector && length(x) == length(fv) || # fallback
+    error("Mismatched broadcast axes: $(axes(x)) != $(axes(fv))")
+check_broadcast_args_axes(_) = true
+@inline check_broadcast_args_axes(fv, arg, args...) =
+    check_broadcast_axes(fv, arg) && check_broadcast_args_axes(fv, args...)
+
+block_view(x, _, _) = x
+Base.@propagate_inbounds block_view(
+    x::BlockArrays.AbstractBlockVector,
+    ::Val{index},
+    _,
+) where {index} = BlockArrays.viewblock(x, BlockArrays.Block(index))
+Base.@propagate_inbounds block_view(x::AbstractVector, val, fv) =
+    block_view(BlockArrays.PseudoBlockVector(x, axes(fv)), val, fv)
+Base.@propagate_inbounds block_view(x::Base.Broadcast.Broadcasted, val, fv) =
+    Base.Broadcast.broadcasted(x.f, block_view_args(val, fv, x.args...)...)
+block_view_args(_, _) = ()
+Base.@propagate_inbounds block_view_args(val, fv, arg, args...) =
+    (block_view(arg, val, fv), block_view_args(val, fv, args...)...)
+
+# Bypass the method for materialize! defined in BlockArrays, which results in
+# type instabilities and allocations.
+Base.materialize!(
+    dest::Fields.FieldVector,
+    bc::Base.Broadcast.Broadcasted{NewFieldVectorStyle},
+) = materialize_field_vector!(dest, bc)
+Base.materialize!(
+    dest::AbstractVector,
+    bc::Base.Broadcast.Broadcasted{NewFieldVectorStyle},
+) = materialize_field_vector!(dest, bc)
+Base.materialize!(
+    dest::Fields.FieldVector,
+    bc::Base.Broadcast.Broadcasted{<:Base.Broadcast.AbstractArrayStyle{0}},
+) = materialize_field_vector!(dest, bc)
+Base.materialize!(
+    dest::Fields.FieldVector,
+    bc::Base.Broadcast.Broadcasted{<:Base.Broadcast.AbstractArrayStyle{1}},
+) = materialize_field_vector!(dest, bc)
+Base.materialize!(
+    dest::Fields.FieldVector,
+    bc::Base.Broadcast.Broadcasted{<:BlockArrays.AbstractBlockStyle{1}},
+) = materialize_field_vector!(dest, bc) # prevents an ambiguity with BlockArrays
+function materialize_field_vector!(dest, bc)
+    fv = first_field_vector(dest, bc)
+    check_broadcast_args_axes(fv, dest, bc)
+    block_index_vals = ntuple(Val, BlockArrays.blocksize(fv, 1))
+    MatrixFields.unrolled_foreach(block_index_vals) do val
+        dest_view = block_view(dest, val, fv)
+        bc_view = block_view(bc, val, fv)
+        @inbounds Base.materialize!(dest_view, bc_view)
+    end
+    return dest
+end
+Base.@propagate_inbounds function Base.getindex(
+    fv::Fields.FieldVector,
+    i::Integer,
+)
+    @warn "scalar indexing into FieldVector (this can be very slow)" maxlog = 1
+    getindex(fv, BlockArrays.findblockindex(axes(fv, 1), i))
+end
+Base.@propagate_inbounds function Base.setindex!(
+    fv::Fields.FieldVector,
+    val,
+    i::Integer,
+)
+    @warn "scalar indexing into FieldVector (this can be very slow)" maxlog = 1
+    setindex!(fv, val, BlockArrays.findblockindex(axes(fv, 1), i))
+end
+
+# TODO: Move all of the following to ClimaCore.jl
+
+filtered_names(f::F, x) where {F} = filtered_names_at_name(f, x, @name())
+function filtered_names_at_name(f::F, x, name) where {F}
+    field = MatrixFields.get_field(x, name)
+    f(field) && return (name,)
+    internal_names = MatrixFields.top_level_names(field)
+    isempty(internal_names) && return ()
+    tuples_of_names = MatrixFields.unrolled_map(internal_names) do internal_name
+        Base.@_inline_meta
+        child_name = MatrixFields.append_internal_name(name, internal_name)
+        filtered_names_at_name(f, x, child_name)
+    end
+    return MatrixFields.unrolled_flatten(tuples_of_names)
+end
+if hasfield(Method, :recursion_relation)
+    dont_limit = (args...) -> true
+    for m in methods(filtered_names_at_name)
+        m.recursion_relation = dont_limit
+    end
+end
+
+scalar_field_names(fv) =
+    filtered_names(x -> x isa Fields.Field && eltype(x) == eltype(fv), fv)
+
+function scalar_level_iterator(field_vector)
+    Iterators.flatmap(scalar_field_names(field_vector)) do name
+        field = MatrixFields.get_field(field_vector, name)
+        if field isa Fields.SpectralElementField
+            (field,)
+        else
+            Iterators.map(1:Spaces.nlevels(axes(field))) do v
+                Fields.level(field, v - 1 + Operators.left_idx(axes(field)))
+            end
+        end
+    end
+end
+
+function column_iterator(field_vector)
+    example_field = field_vector.c # TODO: Generalize this.
+    horz_space = Spaces.horizontal_space(axes(example_field))
+    qs = 1:Quadratures.degrees_of_freedom(Spaces.quadrature_style(horz_space))
+    hs = Spaces.eachslabindex(horz_space)
+    return if Fields.field_values(example_field) isa DataLayouts.VIFH
+        Iterators.map(Iterators.product(qs, hs)) do (i, h)
+            field_vector[Fields.ColumnIndex((i,), h)]
+        end
+    else
+        @assert Fields.field_values(example_field) isa DataLayouts.VIJFH
+        Iterators.map(Iterators.product(qs, qs, hs)) do (i, j, h)
+            field_vector[Fields.ColumnIndex((i, j), h)]
+        end
+    end
+end
diff --git a/src/prognostic_equations/implicit/exact_jacobian.jl b/src/prognostic_equations/implicit/exact_jacobian.jl
new file mode 100644
index 00000000000..2e565313ad8
--- /dev/null
+++ b/src/prognostic_equations/implicit/exact_jacobian.jl
@@ -0,0 +1,196 @@
+@kwdef struct ExactJacobian <: JacobianAlgorithm
+    always_update_exact_jacobian::Bool = true
+    preserve_unfactorized_jacobian::Bool = false # only true for DebugJacobian
+end
+
+function jacobian_cache(alg::ExactJacobian, Y, p)
+    FT = eltype(Y)
+    DA = ClimaComms.array_type(Y)
+    Y_columns = column_iterator(Y)
+    n_columns = length(Y_columns)
+    n_εs = length(first(Y_columns))
+    dual_type = ForwardDiff.Dual{ExactJacobian, FT, 1}
+
+    similar_with_dual_fields(named_tuple) =
+        Fields._values(similar(Fields.FieldVector(; named_tuple...), dual_type))
+
+    Y_dual = similar(Y, dual_type)
+    Y_dual_copy = similar(Y_dual)
+    Y_err_dual = similar(Y_dual)
+    # precomputed_dual =
+    #     similar_with_dual_fields(precomputed_quantities(Y, p.atmos))
+    implicit_precomputed_dual =
+        similar_with_dual_fields(implicit_precomputed_quantities(Y, p.atmos))
+    scratch_dual = similar_with_dual_fields(p.scratch)
+
+    column_partials_vector = DA{ForwardDiff.Partials{1, FT}}(undef, n_εs)
+    column_vectors = DA{FT}(undef, n_columns, n_εs)
+    column_matrices = DA{FT}(undef, n_columns, n_εs, n_εs)
+    column_factorized_matrices =
+        alg.preserve_unfactorized_jacobian ? copy(column_matrices) :
+        column_matrices
+    lu_cache = DA{FT}(undef, n_columns)
+
+    return (;
+        Y_dual,
+        Y_dual_copy,
+        Y_err_dual,
+        implicit_precomputed_dual, # precomputed_dual,
+        scratch_dual,
+        column_partials_vector,
+        column_vectors,
+        column_matrices,
+        column_factorized_matrices,
+        lu_cache,
+    )
+end
+
+always_update_exact_jacobian(alg::ExactJacobian) =
+    alg.always_update_exact_jacobian
+
+# function factorize_exact_jacobian!(alg::ExactJacobian, cache, Y, p, dtγ, t)
+#     ...
+#     Y_dual .= Y
+#     zeros_tuple = ntuple(_ -> 0, n_εs)
+#     for (ε_index, Y_dual_scalar_level) in enumerate(Y_dual_scalar_levels)
+#         ε_partials = Base.setindex(zeros_tuple, 1, ε_index)
+#         ε = ForwardDiff.Dual{ExactJacobian}(0, ε_partials...)
+#         parent(Y_dual_scalar_level) .+= ε
+#     end # add a unique infinitesimal ε to the scalar values on each level of Y
+#     Y_dual_copy .= Y_dual # need a copy because Y_dual gets changed at boundary
+#     set_implicit_precomputed_quantities!(Y_dual, p_dual, t) # changes Y_dual
+#     implicit_tendency!(Y_err_dual, Y_dual, p_dual, t) # the tendency Jacobian
+#     @. Y_err_dual = dtγ * Y_err_dual - Y_dual_copy # the residual Jacobian
+#     for (col_index, Y_err_dual_column) in enumerate(Y_err_dual_columns)
+#         column_partials_vector = cache.column_partials_vectors[col_index]
+#         column_matrix = cache.column_matrices[col_index]
+#         column_factorized_matrix = cache.column_factorized_matrices[col_index]
+#         column_partials_vector .= ForwardDiff.partials.(Y_err_dual_column)
+#         column_matrix .= reinterpret(reshape, FT, column_partials_vector)'
+#     end # TODO: Parallelize this loop
+#     ...
+# end
+function factorize_exact_jacobian!(alg::ExactJacobian, cache, Y, p, dtγ, t)
+    (; preserve_unfactorized_jacobian) = alg
+    (; Y_dual, Y_dual_copy, Y_err_dual) = cache
+    # (; precomputed_dual, scratch_dual, column_partials_vector) = cache
+    (; implicit_precomputed_dual, scratch_dual, column_partials_vector) = cache
+    (; column_matrices, column_factorized_matrices, lu_cache) = cache
+
+    set_precomputed_quantities!(Y, p, t)
+    p_dual_args = ntuple(Val(fieldcount(typeof(p)))) do cache_field_index
+        cache_field_name = fieldname(typeof(p), cache_field_index)
+        if cache_field_name == :precomputed
+            (; p.precomputed..., implicit_precomputed_dual...) # precomputed_dual
+        elseif cache_field_name == :scratch
+            scratch_dual
+        else
+            getfield(p, cache_field_index)
+        end
+    end
+    p_dual = AtmosCache(p_dual_args...)
+
+    FT = eltype(Y)
+    Y_dual_scalar_levels = scalar_level_iterator(Y_dual)
+    Y_err_dual_columns = column_iterator(Y_err_dual)
+    ε = ForwardDiff.Dual{ExactJacobian}(0, 1)
+
+    Y_dual .= Y
+    for (ε_index, Y_dual_scalar_level) in enumerate(Y_dual_scalar_levels)
+        parent(Y_dual_scalar_level) .+= ε # add ε to this level of Y_dual
+        Y_dual_copy .= Y_dual # make a copy because Y_dual gets changed at boundary
+        # set_precomputed_quantities!(Y_dual, p_dual, t)
+        set_implicit_precomputed_quantities!(Y_dual, p_dual, t) # changes Y_dual
+        implicit_tendency!(Y_err_dual, Y_dual, p_dual, t) # get tendency derivative
+        @. Y_err_dual = dtγ * Y_err_dual - Y_dual_copy # get residual derivative
+        parent(Y_dual_scalar_level) .-= ε # remove ε from this level of Y_dual
+
+        for (col_index, Y_err_dual_column) in enumerate(Y_err_dual_columns)
+            column_partials_vector .= ForwardDiff.partials.(Y_err_dual_column)
+            view(column_matrices, col_index, :, ε_index) .=
+                reinterpret(reshape, FT, column_partials_vector)
+        end # TODO: Parallelize this loop.
+    end # TODO: Make this loop use batches of up to 32 values of ε at a time.
+
+    first_column_matrix = view(cache.column_matrices, 1, :, :)
+    first_column_matrix_cpu =
+        ClimaComms.allowscalar(Array, ClimaComms.device(Y), first_column_matrix)
+    p.scratch.temp_matrix .= (first_column_matrix_cpu + I) / dtγ
+
+    if preserve_unfactorized_jacobian
+        column_factorized_matrices .= column_matrices
+    end
+    parallel_lu_factorize!(column_factorized_matrices, lu_cache)
+end
+
+approximate_jacobian!(::ExactJacobian, _, _, _, _, _) = nothing
+
+function invert_jacobian!(::ExactJacobian, cache, x, b)
+    (; column_vectors, column_factorized_matrices, lu_cache) = cache
+    for (col_index, b_column) in enumerate(column_iterator(b))
+        view(column_vectors, col_index, :) .= b_column
+    end # TODO: Parallelize this loop.
+    parallel_lu_solve!(column_vectors, column_factorized_matrices, lu_cache)
+    for (col_index, x_column) in enumerate(column_iterator(x))
+        x_column .= view(column_vectors, col_index, :)
+    end # TODO: Parallelize this loop.
+end
+
+# Set the derivative of `sqrt(x)` to `iszero(x) ? zero(x) : inv(2 * sqrt(x))` in
+# order to properly handle derivatives of `x * sqrt(x)`. Without this change,
+# the derivative of `x * sqrt(x)` is `NaN` when `x` is zero. This method
+# specializes on the tag `ExactJacobian` because not specializing on any tag
+# overwrites the generic method for `Dual` in `ForwardDiff` and breaks
+# precompilation, while specializing on the default tag `Nothing` causes the
+# type piracy Aqua test to fail.
+@inline function Base.sqrt(d::ForwardDiff.Dual{ExactJacobian})
+    tag = Val{ExactJacobian}()
+    x = ForwardDiff.value(d)
+    partials = ForwardDiff.partials(d)
+    val = sqrt(x)
+    deriv = iszero(x) ? zero(x) : inv(2 * val)
+    return ForwardDiff.dual_definition_retval(tag, val, deriv, partials)
+end
+
+function parallel_lu_factorize!(As, temporary_vector)
+    @assert ndims(As) == 3
+    n = size(As, 2)
+    @assert size(As, 3) == n
+    @inbounds for k in 1:n
+        all(!isnan, view(As, :, k, k)) || error("LU error: NaN on diagonal")
+        all(!iszero, view(As, :, k, k)) || error("LU error: 0 on diagonal")
+        temporary_vector .= inv.(view(As, :, k, k))
+        for i in (k + 1):n
+            view(As, :, i, k) .*= temporary_vector
+        end
+        for j in (k + 1):n
+            for i in (k + 1):n
+                view(As, :, i, j) .-= view(As, :, i, k) .* view(As, :, k, j)
+            end
+        end
+    end
+end
+
+function parallel_lu_solve!(xs, As, temporary_vector)
+    @assert ndims(xs) == 2 && ndims(As) == 3
+    n = size(As, 2)
+    @assert size(xs, 2) == n && size(As, 3) == n
+    @inbounds begin
+        for i in 2:n
+            temporary_vector .= zero(eltype(xs))
+            for j in 1:(i - 1)
+                temporary_vector .+= view(As, :, i, j) .* view(xs, :, j)
+            end
+            view(xs, :, i) .-= temporary_vector
+        end
+        view(xs, :, n) ./= view(As, :, n, n)
+        for i in (n - 1):-1:1
+            temporary_vector .= zero(eltype(xs))
+            for j in (i + 1):n
+                temporary_vector .+= view(As, :, i, j) .* view(xs, :, j)
+            end
+            view(xs, :, i) .-= temporary_vector
+            view(xs, :, i) ./= view(As, :, i, i)
+        end
+    end
+end
diff --git a/src/prognostic_equations/implicit/implicit_solver.jl b/src/prognostic_equations/implicit/implicit_solver.jl
index d8f721e9679..c818aeb0fd2 100644
--- a/src/prognostic_equations/implicit/implicit_solver.jl
+++ b/src/prognostic_equations/implicit/implicit_solver.jl
@@ -1,28 +1,25 @@
-import LinearAlgebra: I, Adjoint, ldiv!
+import BlockArrays, ForwardDiff
+import LinearAlgebra: I, Adjoint, LU, lu, lu!, ldiv!
 import ClimaCore.MatrixFields: @name
-using ClimaCore.MatrixFields
-
-abstract type DerivativeFlag end
-struct UseDerivative <: DerivativeFlag end
-struct IgnoreDerivative <: DerivativeFlag end
-use_derivative(::UseDerivative) = true
-use_derivative(::IgnoreDerivative) = false
 
 """
-    ImplicitEquationJacobian(
-        Y, atmos; 
-        approximate_solve_iters, diffusion_flag, topography_flag, sgs_advection_flag, transform_flag
-    )
+    JacobianAlgorithm
 
-A wrapper for the matrix ``∂E/∂Y``, where ``E(Y)`` is the "error" of the
-implicit step with the state ``Y``.
+A description of how to compute the matrix ``∂R/∂Y``, where ``R(Y)`` denotes the
+residual of an implicit step with the state ``Y``. Concrete implementations of
+this abstract type should define 4 methods:
+ - `jacobian_cache(alg::JacobianAlgorithm, Y, p)`
+ - `always_update_exact_jacobian(alg::JacobianAlgorithm)`
+ - `factorize_exact_jacobian!(alg::JacobianAlgorithm, cache, Y, p, dtγ, t)`
+ - `approximate_jacobian!(alg::JacobianAlgorithm, cache, Y, p, dtγ, t)`
+ - `invert_jacobian!(alg::JacobianAlgorithm, cache, x, b)`
 
 # Background
 
 When we use an implicit or split implicit-explicit (IMEX) timestepping scheme,
-we end up with a nonlinear equation of the form ``E(Y) = 0``, where
+we end up with a nonlinear equation of the form ``R(Y) = 0``, where
 ```math
-    E(Y) = Y_{imp}(Y) - Y = \\hat{Y} + Δt * T_{imp}(Y) - Y.
+    R(Y) = Y_{imp}(Y) - Y = \\hat{Y} + Δt * T_{imp}(Y) - Y.
 ```
 In this expression, ``Y_{imp}(Y)`` denotes the state at some time ``t + Δt``.
 This can be expressed as the sum of ``\\hat{Y}``, the contribution from the
@@ -38,827 +35,90 @@ divided into several sub-steps or "stages", where the duration of stage ``i`` is
 ``Δt * γ_i`` for some constant ``γ_i`` between 0 and 1.
 
 In order to solve this equation using Newton's method, we must specify the
-derivative ``∂E/∂Y``. Since ``\\hat{Y}`` does not depend on ``Y`` (it is only a
+derivative ``∂R/∂Y``. Since ``\\hat{Y}`` does not depend on ``Y`` (it is only a
 function of the state at or before time ``t``), this derivative is
 ```math
-    E'(Y) = Δt * T_{imp}'(Y) - I.
+    R'(Y) = Δt * T_{imp}'(Y) - I.
 ```
-In addition, we must specify how to divide ``E(Y)`` by this derivative, i.e.,
+In addition, we must specify how to divide ``R(Y)`` by this derivative, i.e.,
 how to solve the linear equation
 ```math
-    E'(Y) * ΔY = E(Y).
+    R'(Y) * ΔY = R(Y).
 ```
 
 Note: This equation comes from assuming that there is some ``ΔY`` such that
-``E(Y - ΔY) = 0`` and making the first-order approximation
+``R(Y - ΔY) = 0`` and making the first-order approximation
 ```math
-    E(Y - ΔY) \\approx E(Y) - E'(Y) * ΔY.
+    R(Y - ΔY) \\approx R(Y) - R'(Y) * ΔY.
 ```
 
 After initializing ``Y`` to ``Y[0] = \\hat{Y}``, Newton's method executes the
 following steps:
-- Compute the derivative ``E'(Y[0])``.
-- Compute the implicit tendency ``T_{imp}(Y[0])`` and use it to get ``E(Y[0])``.
-- Solve the linear equation ``E'(Y[0]) * ΔY[0] = E(Y[0])`` for ``ΔY[0]``.
+- Compute the derivative ``R'(Y[0])``.
+- Compute the implicit tendency ``T_{imp}(Y[0])`` and use it to get ``R(Y[0])``.
+- Solve the linear equation ``R'(Y[0]) * ΔY[0] = R(Y[0])`` for ``ΔY[0]``.
 - Update ``Y`` to ``Y[1] = Y[0] - ΔY[0]``.
 If the number of Newton iterations is limited to 1, this new value of ``Y`` is
 taken to be the solution of the implicit equation. Otherwise, this sequence of
 steps is repeated, i.e., ``ΔY[1]`` is computed and used to update ``Y`` to
 ``Y[2] = Y[1] - ΔY[1]``, then ``ΔY[2]`` is computed and used to update ``Y`` to
 ``Y[3] = Y[2] - ΔY[2]``, and so on. The iterative process is terminated either
-when the error ``E(Y)`` is sufficiently close to 0 (according to the convergence
-condition passed to Newton's method), or when the maximum number of iterations
-is reached.
-
-# Arguments
-
-- `Y::FieldVector`: the state of the simulation
-- `atmos::AtmosModel`: the model configuration
-- `approximate_solve_iters::Int`: number of iterations to take for the
-  approximate linear solve required when `diffusion_flag = UseDerivative()`
-- `diffusion_flag::DerivativeFlag`: whether the derivative of the
-  diffusion tendency with respect to the quantities being diffused should be
-  computed or approximated as 0; must be either `UseDerivative()` or
-  `IgnoreDerivative()` instead of a `Bool` to ensure type-stability
-- `topography_flag::DerivativeFlag`: whether the derivative of vertical
-  contravariant velocity with respect to horizontal covariant velocity should be
-  computed or approximated as 0; must be either `UseDerivative()` or
-  `IgnoreDerivative()` instead of a `Bool` to ensure type-stability
-- `sgs_advection_flag::DerivativeFlag`: whether the derivative of the
-  subgrid-scale advection tendency with respect to the updraft quantities should be
-  computed or approximated as 0; must be either `UseDerivative()` or
-  `IgnoreDerivative()` instead of a `Bool` to ensure type-stability
-- `transform_flag::Bool`: whether the error should be transformed from ``E(Y)``
-  to ``E(Y)/Δt``, which is required for non-Rosenbrock timestepping schemes from
-  OrdinaryDiffEq.jl
+when the residual ``R(Y)`` is sufficiently close to 0 (according to the
+convergence condition passed to Newton's method), or when the maximum number of
+iterations is reached.
 """
-struct ImplicitEquationJacobian{
-    M <: MatrixFields.FieldMatrix,
-    S <: MatrixFields.FieldMatrixSolver,
-    F1 <: DerivativeFlag,
-    F2 <: DerivativeFlag,
-    F3 <: DerivativeFlag,
-    T <: Fields.FieldVector,
-    R <: Base.RefValue,
-}
-    # stores the matrix E'(Y) = Δt * T_imp'(Y) - I
-    matrix::M
-
-    # solves the linear equation E'(Y) * ΔY = E(Y) for ΔY
-    solver::S
+abstract type JacobianAlgorithm end
 
-    # flags that determine how E'(Y) is approximated
-    diffusion_flag::F1
-    topography_flag::F2
-    sgs_advection_flag::F3
-
-    # required by Krylov.jl to evaluate ldiv! with AbstractVector inputs
-    temp_b::T
-    temp_x::T
-
-    # required by OrdinaryDiffEq.jl to run non-Rosenbrock timestepping schemes
-    transform_flag::Bool
-    dtγ_ref::R
+struct ImplicitEquationJacobian{A <: JacobianAlgorithm, C}
+    alg::A
+    cache::C
 end
-
-function ImplicitEquationJacobian(
-    Y,
-    atmos;
-    approximate_solve_iters = 1,
-    diffusion_flag = atmos.diff_mode == Implicit() ? UseDerivative() :
-                     IgnoreDerivative(),
-    topography_flag = has_topography(axes(Y.c)) ? UseDerivative() :
-                      IgnoreDerivative(),
-    sgs_advection_flag = atmos.sgs_adv_mode == Implicit() ? UseDerivative() :
-                         IgnoreDerivative(),
-    transform_flag = false,
-)
-    FT = Spaces.undertype(axes(Y.c))
-    CTh = CTh_vector_type(axes(Y.c))
-
-    DiagonalRow = DiagonalMatrixRow{FT}
-    TridiagonalRow = TridiagonalMatrixRow{FT}
-    BidiagonalRow_C3 = BidiagonalMatrixRow{C3{FT}}
-    TridiagonalRow_ACTh = TridiagonalMatrixRow{Adjoint{FT, CTh{FT}}}
-    BidiagonalRow_ACT3 = BidiagonalMatrixRow{Adjoint{FT, CT3{FT}}}
-    BidiagonalRow_C3xACTh =
-        BidiagonalMatrixRow{typeof(zero(C3{FT}) * zero(CTh{FT})')}
-    DiagonalRow_C3xACT3 =
-        DiagonalMatrixRow{typeof(zero(C3{FT}) * zero(CT3{FT})')}
-    TridiagonalRow_C3xACT3 =
-        TridiagonalMatrixRow{typeof(zero(C3{FT}) * zero(CT3{FT})')}
-
-    is_in_Y(name) = MatrixFields.has_field(Y, name)
-
-    ρq_tot_if_available = is_in_Y(@name(c.ρq_tot)) ? (@name(c.ρq_tot),) : ()
-    ρatke_if_available =
-        is_in_Y(@name(c.sgs⁰.ρatke)) ? (@name(c.sgs⁰.ρatke),) : ()
-    sfc_if_available = is_in_Y(@name(sfc)) ? (@name(sfc),) : ()
-
-    tracer_names = (
-        @name(c.ρq_tot),
-        @name(c.ρq_liq),
-        @name(c.ρq_ice),
-        @name(c.ρq_rai),
-        @name(c.ρq_sno),
-    )
-    available_tracer_names = MatrixFields.unrolled_filter(is_in_Y, tracer_names)
-
-    # Note: We have to use FT(-1) * I instead of -I because inv(-1) == -1.0,
-    # which means that multiplying inv(-1) by a Float32 will yield a Float64.
-    identity_blocks = MatrixFields.unrolled_map(
-        name -> (name, name) => FT(-1) * I,
-        (@name(c.ρ), sfc_if_available...),
+function ImplicitEquationJacobian(alg, Y, p)
+    cache = (;
+        jacobian_cache(alg, Y, p)...,
+        x_krylov = similar(Y),
+        b_krylov = similar(Y),
     )
+    return ImplicitEquationJacobian(alg, cache)
+end
 
-    active_scalar_names = (@name(c.ρ), @name(c.ρe_tot), ρq_tot_if_available...)
-    advection_blocks = (
-        (
-            use_derivative(topography_flag) ?
-            MatrixFields.unrolled_map(
-                name ->
-                    (name, @name(c.uₕ)) =>
-                        similar(Y.c, TridiagonalRow_ACTh),
-                active_scalar_names,
-            ) : ()
-        )...,
-        MatrixFields.unrolled_map(
-            name -> (name, @name(f.u₃)) => similar(Y.c, BidiagonalRow_ACT3),
-            active_scalar_names,
-        )...,
-        MatrixFields.unrolled_map(
-            name -> (@name(f.u₃), name) => similar(Y.f, BidiagonalRow_C3),
-            active_scalar_names,
-        )...,
-        (@name(f.u₃), @name(c.uₕ)) => similar(Y.f, BidiagonalRow_C3xACTh),
-        (@name(f.u₃), @name(f.u₃)) => similar(Y.f, TridiagonalRow_C3xACT3),
-    )
-
-    diffused_scalar_names =
-        (@name(c.ρe_tot), available_tracer_names..., ρatke_if_available...)
-    diffusion_blocks = if use_derivative(diffusion_flag)
-        (
-            MatrixFields.unrolled_map(
-                name -> (name, @name(c.ρ)) => similar(Y.c, TridiagonalRow),
-                diffused_scalar_names,
-            )...,
-            MatrixFields.unrolled_map(
-                name -> (name, name) => similar(Y.c, TridiagonalRow),
-                diffused_scalar_names,
-            )...,
-            (
-                is_in_Y(@name(c.ρq_tot)) ?
-                (
-                    (@name(c.ρe_tot), @name(c.ρq_tot)) =>
-                        similar(Y.c, TridiagonalRow),
-                ) : ()
-            )...,
-            (@name(c.uₕ), @name(c.uₕ)) =>
-                !isnothing(atmos.turbconv_model) ||
-                    diffuse_momentum(atmos.vert_diff) ?
-                similar(Y.c, TridiagonalRow) : FT(-1) * I,
-        )
-    else
-        MatrixFields.unrolled_map(
-            name -> (name, name) => FT(-1) * I,
-            (diffused_scalar_names..., @name(c.uₕ)),
-        )
-    end
-
-    sgs_advection_blocks = if atmos.turbconv_model isa PrognosticEDMFX
-        @assert n_prognostic_mass_flux_subdomains(atmos.turbconv_model) == 1
-        sgs_scalar_names = (
-            @name(c.sgsʲs.:(1).q_tot),
-            @name(c.sgsʲs.:(1).mse),
-            @name(c.sgsʲs.:(1).ρa)
-        )
-        if use_derivative(sgs_advection_flag)
-            (
-                MatrixFields.unrolled_map(
-                    name -> (name, name) => similar(Y.c, TridiagonalRow),
-                    sgs_scalar_names,
-                )...,
-                (@name(c.sgsʲs.:(1).mse), @name(c.ρ)) =>
-                    similar(Y.c, DiagonalRow),
-                (@name(c.sgsʲs.:(1).mse), @name(c.sgsʲs.:(1).q_tot)) =>
-                    similar(Y.c, DiagonalRow),
-                (@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).q_tot)) =>
-                    similar(Y.c, TridiagonalRow),
-                (@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).mse)) =>
-                    similar(Y.c, TridiagonalRow),
-                (@name(f.sgsʲs.:(1).u₃), @name(c.ρ)) =>
-                    similar(Y.f, BidiagonalRow_C3),
-                (@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).q_tot)) =>
-                    similar(Y.f, BidiagonalRow_C3),
-                (@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).mse)) =>
-                    similar(Y.f, BidiagonalRow_C3),
-                (@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)) =>
-                    similar(Y.f, TridiagonalRow_C3xACT3),
-            )
-        else
-            (
-                MatrixFields.unrolled_map(
-                    name -> (name, name) => FT(-1) * I,
-                    sgs_scalar_names,
-                )...,
-                (@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)) =>
-                    !isnothing(atmos.rayleigh_sponge) ?
-                    similar(Y.f, DiagonalRow_C3xACT3) : FT(-1) * I,
-            )
-        end
-    else
-        ()
-    end
+# We only use A, but ClimaTimeSteppers.jl currently requires us to pass
+# jac_prototype and then calls similar(jac_prototype) to obtain A. Since
+# this was only done for compatibility with OrdinaryDiffEq.jl, which we are no
+# longer using, we should fix this in a new version of ClimaTimeSteppers.jl.
+Base.similar(A::ImplicitEquationJacobian) = A
 
-    matrix = MatrixFields.FieldMatrix(
-        identity_blocks...,
-        sgs_advection_blocks...,
-        advection_blocks...,
-        diffusion_blocks...,
-    )
+# This is called from a callback when always_update_exact_jacobian is false.
+NVTX.@annotate function update_exact_jacobian!(A, Y, p, dtγ, t)
+    FT = eltype(Y)
+    factorize_exact_jacobian!(A.alg, A.cache, Y, p, FT(dtγ), t)
+end
 
-    sgs_names_if_available = if atmos.turbconv_model isa PrognosticEDMFX
-        (
-            @name(c.sgsʲs.:(1).q_tot),
-            @name(c.sgsʲs.:(1).mse),
-            @name(c.sgsʲs.:(1).ρa),
-            @name(f.sgsʲs.:(1).u₃),
-        )
-    else
-        ()
+# This is passed to ClimaTimeSteppers.jl and called on each Newton iteration.
+NVTX.@annotate function update_jacobian!(A, Y, p, dtγ, t)
+    FT = eltype(Y)
+    if always_update_exact_jacobian(A.alg)
+        factorize_exact_jacobian!(A.alg, A.cache, Y, p, FT(dtγ), t)
     end
-    names₁_group₁ = (@name(c.ρ), sfc_if_available...)
-    names₁_group₂ = (available_tracer_names..., ρatke_if_available...)
-    names₁_group₃ = (@name(c.ρe_tot),)
-    names₁ = (
-        names₁_group₁...,
-        sgs_names_if_available...,
-        names₁_group₂...,
-        names₁_group₃...,
-    )
-
-    alg₂ = MatrixFields.BlockLowerTriangularSolve(@name(c.uₕ))
-    alg =
-        if use_derivative(diffusion_flag) || use_derivative(sgs_advection_flag)
-            alg₁_subalg₂ =
-                if atmos.turbconv_model isa PrognosticEDMFX &&
-                   use_derivative(sgs_advection_flag)
-                    diff_subalg =
-                        use_derivative(diffusion_flag) ?
-                        (;
-                            alg₂ = MatrixFields.BlockLowerTriangularSolve(
-                                names₁_group₂...,
-                            )
-                        ) : (;)
-                    (;
-                        alg₂ = MatrixFields.BlockLowerTriangularSolve(
-                            @name(c.sgsʲs.:(1).q_tot);
-                            alg₂ = MatrixFields.BlockLowerTriangularSolve(
-                                @name(c.sgsʲs.:(1).mse);
-                                alg₂ = MatrixFields.BlockLowerTriangularSolve(
-                                    @name(c.sgsʲs.:(1).ρa),
-                                    @name(f.sgsʲs.:(1).u₃);
-                                    diff_subalg...,
-                                ),
-                            ),
-                        )
-                    )
-                else
-                    is_in_Y(@name(c.ρq_tot)) ?
-                    (;
-                        alg₂ = MatrixFields.BlockLowerTriangularSolve(
-                            names₁_group₂...,
-                        )
-                    ) : (;)
-                end
-            alg₁ = MatrixFields.BlockLowerTriangularSolve(
-                names₁_group₁...;
-                alg₁_subalg₂...,
-            )
-            MatrixFields.ApproximateBlockArrowheadIterativeSolve(
-                names₁...;
-                alg₁,
-                alg₂,
-                P_alg₁ = MatrixFields.MainDiagonalPreconditioner(),
-                n_iters = approximate_solve_iters,
-            )
-        else
-            MatrixFields.BlockArrowheadSolve(names₁...; alg₂)
-        end
-
-    return ImplicitEquationJacobian(
-        matrix,
-        MatrixFields.FieldMatrixSolver(alg, matrix, Y),
-        diffusion_flag,
-        topography_flag,
-        sgs_advection_flag,
-        similar(Y),
-        similar(Y),
-        transform_flag,
-        Ref{FT}(),
-    )
+    approximate_jacobian!(A.alg, A.cache, Y, p, FT(dtγ), t)
 end
 
-# We only use A, but ClimaTimeSteppers.jl require us to
-# pass jac_prototype and then call similar(jac_prototype) to
-# obtain A. This is a workaround to avoid unnecessary allocations.
-Base.similar(A::ImplicitEquationJacobian) = A
-
-# This method specifies how to solve the equation E'(Y) * ΔY = E(Y) for ΔY.
-NVTX.@annotate function ldiv!(
+# This is called by ClimaTimeSteppers.jl on each Newton iteration.
+NVTX.@annotate ldiv!(
     x::Fields.FieldVector,
     A::ImplicitEquationJacobian,
     b::Fields.FieldVector,
-)
-    MatrixFields.field_matrix_solve!(A.solver, x, A.matrix, b)
-    if A.transform_flag
-        @. x *= -A.dtγ_ref[]
-    end
-end
+) = invert_jacobian!(A.alg, A.cache, x, b)
 
-# This method for ldiv! is called by Krylov.jl from inside ClimaTimeSteppers.jl.
-# See https://github.com/JuliaSmoothOptimizers/Krylov.jl/issues/605 for a
-# related issue that requires the same workaround.
-NVTX.@annotate function ldiv!(
+# This is called by Krylov.jl from inside ClimaTimeSteppers.jl. See
+# https://github.com/JuliaSmoothOptimizers/Krylov.jl/issues/605 for a related
+# issue that requires the same workaround.
+function ldiv!(
     x::AbstractVector,
     A::ImplicitEquationJacobian,
     b::AbstractVector,
 )
-    A.temp_b .= b
-    ldiv!(A.temp_x, A, A.temp_b)
-    x .= A.temp_x
-end
-
-# This function is used by DiffEqBase.jl instead of ldiv!.
-linsolve!(::Type{Val{:init}}, f, u0; kwargs...) = _linsolve!
-_linsolve!(x, A, b, update_matrix = false; kwargs...) = ldiv!(x, A, b)
-
-# This method specifies how to compute E'(Y), which is referred to as "Wfact" in
-# DiffEqBase.jl.
-NVTX.@annotate function Wfact!(A, Y, p, dtγ, t)
-    # Remove unnecessary values from p to avoid allocations in bycolumn.
-    p′ = (;
-        p.precomputed.ᶜspecific,
-        p.precomputed.ᶜK,
-        p.precomputed.ᶜts,
-        p.precomputed.ᶜp,
-        (
-            p.atmos.precip_model isa Microphysics1Moment ?
-            (; p.precomputed.ᶜwᵣ, p.precomputed.ᶜwₛ) : (;)
-        )...,
-        p.precomputed.ᶜh_tot,
-        (
-            use_derivative(A.diffusion_flag) ?
-            (; p.precomputed.ᶜK_u, p.precomputed.ᶜK_h) : (;)
-        )...,
-        (
-            use_derivative(A.diffusion_flag) &&
-            p.atmos.turbconv_model isa AbstractEDMF ?
-            (; p.precomputed.ᶜtke⁰, p.precomputed.ᶜmixing_length) : (;)
-        )...,
-        (
-            use_derivative(A.diffusion_flag) &&
-            p.atmos.turbconv_model isa PrognosticEDMFX ?
-            (; p.precomputed.ᶜρa⁰) : (;)
-        )...,
-        (
-            use_derivative(A.sgs_advection_flag) &&
-            p.atmos.turbconv_model isa PrognosticEDMFX ?
-            (;
-                p.core.ᶜgradᵥ_ᶠΦ,
-                p.precomputed.ᶜρʲs,
-                p.precomputed.ᶠu³ʲs,
-                p.precomputed.ᶜtsʲs,
-                p.precomputed.bdmr_l,
-                p.precomputed.bdmr_r,
-                p.precomputed.bdmr,
-            ) : (;)
-        )...,
-        p.core.ᶜΦ,
-        p.core.ᶠgradᵥ_ᶜΦ,
-        p.scratch.ᶜtemp_scalar,
-        p.scratch.ᶜtemp_C3,
-        p.scratch.ᶠtemp_CT3,
-        p.scratch.∂ᶜK_∂ᶜuₕ,
-        p.scratch.∂ᶜK_∂ᶠu₃,
-        p.scratch.ᶠp_grad_matrix,
-        p.scratch.ᶜadvection_matrix,
-        p.scratch.ᶜdiffusion_h_matrix,
-        p.scratch.ᶜdiffusion_u_matrix,
-        p.scratch.ᶠbidiagonal_matrix_ct3,
-        p.scratch.ᶠbidiagonal_matrix_ct3_2,
-        p.scratch.ᶠtridiagonal_matrix_c3,
-        p.dt,
-        p.params,
-        p.atmos,
-        (
-            p.atmos.rayleigh_sponge isa RayleighSponge ?
-            (; p.rayleigh_sponge.ᶠβ_rayleigh_w) : (;)
-        )...,
-    )
-
-    # Convert dtγ from a Float64 to an FT.
-    FT = Spaces.undertype(axes(Y.c))
-    dtγ′ = FT(dtγ)
-
-    A.dtγ_ref[] = dtγ′
-    update_implicit_equation_jacobian!(A, Y, p′, dtγ′)
-end
-
-function update_implicit_equation_jacobian!(A, Y, p, dtγ)
-    (; matrix, diffusion_flag, sgs_advection_flag, topography_flag) = A
-    (; ᶜspecific, ᶜK, ᶜts, ᶜp, ᶜΦ, ᶠgradᵥ_ᶜΦ) = p
-    (;
-        ᶜtemp_C3,
-        ∂ᶜK_∂ᶜuₕ,
-        ∂ᶜK_∂ᶠu₃,
-        ᶠp_grad_matrix,
-        ᶜadvection_matrix,
-        ᶠbidiagonal_matrix_ct3,
-        ᶠbidiagonal_matrix_ct3_2,
-        ᶠtridiagonal_matrix_c3,
-    ) = p
-    (; ᶜdiffusion_h_matrix, ᶜdiffusion_u_matrix, params) = p
-    (; edmfx_upwinding) = p.atmos.numerics
-
-    FT = Spaces.undertype(axes(Y.c))
-    CTh = CTh_vector_type(axes(Y.c))
-    one_C3xACT3 = C3(FT(1)) * CT3(FT(1))'
-
-    cv_d = FT(CAP.cv_d(params))
-    Δcv_v = FT(CAP.cv_v(params)) - cv_d
-    T_tri = FT(CAP.T_triple(params))
-    e_int_v0 = T_tri * Δcv_v - FT(CAP.e_int_v0(params))
-    thermo_params = CAP.thermodynamics_params(params)
-
-    ᶜρ = Y.c.ρ
-    ᶜuₕ = Y.c.uₕ
-    ᶠu₃ = Y.f.u₃
-    ᶜJ = Fields.local_geometry_field(Y.c).J
-    ᶜgⁱʲ = Fields.local_geometry_field(Y.c).gⁱʲ
-    ᶠgⁱʲ = Fields.local_geometry_field(Y.f).gⁱʲ
-
-    ᶜkappa_m = p.ᶜtemp_scalar
-    @. ᶜkappa_m =
-        TD.gas_constant_air(thermo_params, ᶜts) / TD.cv_m(thermo_params, ᶜts)
-
-    if use_derivative(topography_flag)
-        @. ∂ᶜK_∂ᶜuₕ = DiagonalMatrixRow(
-            adjoint(CTh(ᶜuₕ)) + adjoint(ᶜinterp(ᶠu₃)) * g³ʰ(ᶜgⁱʲ),
-        )
-    else
-        @. ∂ᶜK_∂ᶜuₕ = DiagonalMatrixRow(adjoint(CTh(ᶜuₕ)))
-    end
-    @. ∂ᶜK_∂ᶠu₃ =
-        ᶜinterp_matrix() ⋅ DiagonalMatrixRow(adjoint(CT3(ᶠu₃))) +
-        DiagonalMatrixRow(adjoint(CT3(ᶜuₕ))) ⋅ ᶜinterp_matrix()
-
-    @. ᶠp_grad_matrix = DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ)) ⋅ ᶠgradᵥ_matrix()
-
-    @. ᶜadvection_matrix =
-        -(ᶜadvdivᵥ_matrix()) ⋅ DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρ))
-
-    if use_derivative(topography_flag)
-        ∂ᶜρ_err_∂ᶜuₕ = matrix[@name(c.ρ), @name(c.uₕ)]
-        @. ∂ᶜρ_err_∂ᶜuₕ =
-            dtγ * ᶜadvection_matrix ⋅ ᶠwinterp_matrix(ᶜJ * ᶜρ) ⋅
-            DiagonalMatrixRow(g³ʰ(ᶜgⁱʲ))
-    end
-    ∂ᶜρ_err_∂ᶠu₃ = matrix[@name(c.ρ), @name(f.u₃)]
-    @. ∂ᶜρ_err_∂ᶠu₃ = dtγ * ᶜadvection_matrix ⋅ DiagonalMatrixRow(g³³(ᶠgⁱʲ))
-
-    tracer_info = (
-        (@name(c.ρe_tot), @name(ᶜh_tot)),
-        (@name(c.ρq_tot), @name(ᶜspecific.q_tot)),
-    )
-    MatrixFields.unrolled_foreach(tracer_info) do (ρχ_name, χ_name)
-        MatrixFields.has_field(Y, ρχ_name) || return
-        ᶜχ = MatrixFields.get_field(p, χ_name)
-        if use_derivative(topography_flag)
-            ∂ᶜρχ_err_∂ᶜuₕ = matrix[ρχ_name, @name(c.uₕ)]
-        end
-        ∂ᶜρχ_err_∂ᶠu₃ = matrix[ρχ_name, @name(f.u₃)]
-        use_derivative(topography_flag) && @. ∂ᶜρχ_err_∂ᶜuₕ =
-            dtγ * ᶜadvection_matrix ⋅ DiagonalMatrixRow(ᶠinterp(ᶜχ)) ⋅
-            ᶠwinterp_matrix(ᶜJ * ᶜρ) ⋅ DiagonalMatrixRow(g³ʰ(ᶜgⁱʲ))
-        @. ∂ᶜρχ_err_∂ᶠu₃ =
-            dtγ * ᶜadvection_matrix ⋅ DiagonalMatrixRow(ᶠinterp(ᶜχ) * g³³(ᶠgⁱʲ))
-    end
-
-    ∂ᶠu₃_err_∂ᶜρ = matrix[@name(f.u₃), @name(c.ρ)]
-    ∂ᶠu₃_err_∂ᶜρe_tot = matrix[@name(f.u₃), @name(c.ρe_tot)]
-    @. ∂ᶠu₃_err_∂ᶜρ =
-        dtγ * (
-            ᶠp_grad_matrix ⋅
-            DiagonalMatrixRow(ᶜkappa_m * (T_tri * cv_d - ᶜK - ᶜΦ)) +
-            DiagonalMatrixRow(ᶠgradᵥ(ᶜp) / abs2(ᶠinterp(ᶜρ))) ⋅
-            ᶠinterp_matrix()
-        )
-    @. ∂ᶠu₃_err_∂ᶜρe_tot = dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow(ᶜkappa_m)
-    if MatrixFields.has_field(Y, @name(c.ρq_tot))
-        ∂ᶠu₃_err_∂ᶜρq_tot = matrix[@name(f.u₃), @name(c.ρq_tot)]
-        @. ∂ᶠu₃_err_∂ᶜρq_tot =
-            dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow(ᶜkappa_m * e_int_v0)
-    end
-
-    ∂ᶠu₃_err_∂ᶜuₕ = matrix[@name(f.u₃), @name(c.uₕ)]
-    ∂ᶠu₃_err_∂ᶠu₃ = matrix[@name(f.u₃), @name(f.u₃)]
-    I_u₃ = DiagonalMatrixRow(one_C3xACT3)
-    @. ∂ᶠu₃_err_∂ᶜuₕ =
-        dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow(-(ᶜkappa_m) * ᶜρ) ⋅ ∂ᶜK_∂ᶜuₕ
-    if p.atmos.rayleigh_sponge isa RayleighSponge
-        @. ∂ᶠu₃_err_∂ᶠu₃ =
-            dtγ * (
-                ᶠp_grad_matrix ⋅ DiagonalMatrixRow(-(ᶜkappa_m) * ᶜρ) ⋅
-                ∂ᶜK_∂ᶠu₃ + DiagonalMatrixRow(-p.ᶠβ_rayleigh_w * (one_C3xACT3,))
-            ) - (I_u₃,)
-    else
-        @. ∂ᶠu₃_err_∂ᶠu₃ =
-            dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow(-(ᶜkappa_m) * ᶜρ) ⋅
-            ∂ᶜK_∂ᶠu₃ - (I_u₃,)
-    end
-
-    if use_derivative(diffusion_flag)
-        (; ᶜK_h, ᶜK_u) = p
-        @. ᶜdiffusion_h_matrix =
-            ᶜadvdivᵥ_matrix() ⋅ DiagonalMatrixRow(ᶠinterp(ᶜρ) * ᶠinterp(ᶜK_h)) ⋅
-            ᶠgradᵥ_matrix()
-        if (
-            MatrixFields.has_field(Y, @name(c.sgs⁰.ρatke)) ||
-            !isnothing(p.atmos.turbconv_model) ||
-            diffuse_momentum(p.atmos.vert_diff)
-        )
-            @. ᶜdiffusion_u_matrix =
-                ᶜadvdivᵥ_matrix() ⋅
-                DiagonalMatrixRow(ᶠinterp(ᶜρ) * ᶠinterp(ᶜK_u)) ⋅ ᶠgradᵥ_matrix()
-        end
-
-        ∂ᶜρe_tot_err_∂ᶜρ = matrix[@name(c.ρe_tot), @name(c.ρ)]
-        ∂ᶜρe_tot_err_∂ᶜρe_tot = matrix[@name(c.ρe_tot), @name(c.ρe_tot)]
-        @. ∂ᶜρe_tot_err_∂ᶜρ =
-            dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow(
-                (
-                    -(1 + ᶜkappa_m) * ᶜspecific.e_tot -
-                    ᶜkappa_m * e_int_v0 * ᶜspecific.q_tot
-                ) / ᶜρ,
-            )
-        @. ∂ᶜρe_tot_err_∂ᶜρe_tot =
-            dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow((1 + ᶜkappa_m) / ᶜρ) -
-            (I,)
-        if MatrixFields.has_field(Y, @name(c.ρq_tot))
-            ∂ᶜρe_tot_err_∂ᶜρq_tot = matrix[@name(c.ρe_tot), @name(c.ρq_tot)]
-            @. ∂ᶜρe_tot_err_∂ᶜρq_tot =
-                dtγ * ᶜdiffusion_h_matrix ⋅
-                DiagonalMatrixRow(ᶜkappa_m * e_int_v0 / ᶜρ)
-        end
-
-        tracer_info = (
-            (@name(c.ρq_tot), @name(q_tot)),
-            (@name(c.ρq_liq), @name(q_liq)),
-            (@name(c.ρq_ice), @name(q_ice)),
-            (@name(c.ρq_rai), @name(q_rai)),
-            (@name(c.ρq_sno), @name(q_sno)),
-        )
-        MatrixFields.unrolled_foreach(tracer_info) do (ρq_name, q_name)
-            MatrixFields.has_field(Y, ρq_name) || return
-            ᶜq = MatrixFields.get_field(ᶜspecific, q_name)
-            ∂ᶜρq_err_∂ᶜρ = matrix[ρq_name, @name(c.ρ)]
-            ∂ᶜρq_err_∂ᶜρq = matrix[ρq_name, ρq_name]
-            @. ∂ᶜρq_err_∂ᶜρ =
-                dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow(-(ᶜq) / ᶜρ)
-            @. ∂ᶜρq_err_∂ᶜρq =
-                dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow(1 / ᶜρ) - (I,)
-        end
-
-        if MatrixFields.has_field(Y, @name(c.sgs⁰.ρatke))
-            turbconv_params = CAP.turbconv_params(params)
-            c_d = CAP.tke_diss_coeff(turbconv_params)
-            (; ᶜtke⁰, ᶜmixing_length, dt) = p
-            ᶜρa⁰ = p.atmos.turbconv_model isa PrognosticEDMFX ? p.ᶜρa⁰ : ᶜρ
-            ᶜρatke⁰ = Y.c.sgs⁰.ρatke
-
-            @inline dissipation_rate(tke⁰, mixing_length) =
-                tke⁰ >= 0 ? c_d * sqrt(tke⁰) / max(mixing_length, 1) : 1 / dt
-            @inline ∂dissipation_rate_∂tke⁰(tke⁰, mixing_length) =
-                tke⁰ > 0 ? c_d / (2 * max(mixing_length, 1) * sqrt(tke⁰)) :
-                typeof(tke⁰)(0)
-
-            ᶜdissipation_matrix_diagonal = p.ᶜtemp_scalar
-            @. ᶜdissipation_matrix_diagonal =
-                ᶜρatke⁰ * ∂dissipation_rate_∂tke⁰(ᶜtke⁰, ᶜmixing_length)
-
-            ∂ᶜρatke⁰_err_∂ᶜρ = matrix[@name(c.sgs⁰.ρatke), @name(c.ρ)]
-            ∂ᶜρatke⁰_err_∂ᶜρatke⁰ =
-                matrix[@name(c.sgs⁰.ρatke), @name(c.sgs⁰.ρatke)]
-            @. ∂ᶜρatke⁰_err_∂ᶜρ =
-                dtγ * (
-                    ᶜdiffusion_u_matrix -
-                    DiagonalMatrixRow(ᶜdissipation_matrix_diagonal)
-                ) ⋅ DiagonalMatrixRow(-(ᶜtke⁰) / ᶜρa⁰)
-            @. ∂ᶜρatke⁰_err_∂ᶜρatke⁰ =
-                dtγ * (
-                    (
-                        ᶜdiffusion_u_matrix -
-                        DiagonalMatrixRow(ᶜdissipation_matrix_diagonal)
-                    ) ⋅ DiagonalMatrixRow(1 / ᶜρa⁰) -
-                    DiagonalMatrixRow(dissipation_rate(ᶜtke⁰, ᶜmixing_length))
-                ) - (I,)
-        end
-
-        if (
-            !isnothing(p.atmos.turbconv_model) ||
-            diffuse_momentum(p.atmos.vert_diff)
-        )
-            ∂ᶜuₕ_err_∂ᶜuₕ = matrix[@name(c.uₕ), @name(c.uₕ)]
-            @. ∂ᶜuₕ_err_∂ᶜuₕ =
-                dtγ * DiagonalMatrixRow(1 / ᶜρ) ⋅ ᶜdiffusion_u_matrix - (I,)
-        end
-
-        ᶠlg = Fields.local_geometry_field(Y.f)
-        precip_info =
-            ((@name(c.ρq_rai), @name(ᶜwᵣ)), (@name(c.ρq_sno), @name(ᶜwₛ)))
-        MatrixFields.unrolled_foreach(precip_info) do (ρqₚ_name, wₚ_name)
-            MatrixFields.has_field(Y, ρqₚ_name) || return
-            ∂ᶜρqₚ_err_∂ᶜρqₚ = matrix[ρqₚ_name, ρqₚ_name]
-            ᶜwₚ = MatrixFields.get_field(p, wₚ_name)
-            ᶠtmp = p.ᶠtemp_CT3
-            @. ᶠtmp = CT3(unit_basis_vector_data(CT3, ᶠlg)) * ᶠwinterp(ᶜJ, ᶜρ)
-            @. ∂ᶜρqₚ_err_∂ᶜρqₚ +=
-                dtγ * -(ᶜprecipdivᵥ_matrix()) ⋅ DiagonalMatrixRow(ᶠtmp) ⋅
-                ᶠright_bias_matrix() ⋅ DiagonalMatrixRow(-(ᶜwₚ) / ᶜρ)
-        end
-    end
-
-    if p.atmos.turbconv_model isa PrognosticEDMFX
-        if use_derivative(sgs_advection_flag)
-            (; ᶜgradᵥ_ᶠΦ, ᶜρʲs, ᶠu³ʲs, ᶜtsʲs) = p
-            (; bdmr_l, bdmr_r, bdmr) = p
-            is_third_order = edmfx_upwinding == Val(:third_order)
-            ᶠupwind = is_third_order ? ᶠupwind3 : ᶠupwind1
-            ᶠset_upwind_bcs = Operators.SetBoundaryOperator(;
-                top = Operators.SetValue(zero(CT3{FT})),
-                bottom = Operators.SetValue(zero(CT3{FT})),
-            ) # Need to wrap ᶠupwind in this for well-defined boundaries.
-            UpwindMatrixRowType =
-                is_third_order ? QuaddiagonalMatrixRow : BidiagonalMatrixRow
-            ᶠupwind_matrix = is_third_order ? ᶠupwind3_matrix : ᶠupwind1_matrix
-            ᶠset_upwind_matrix_bcs = Operators.SetBoundaryOperator(;
-                top = Operators.SetValue(zero(UpwindMatrixRowType{CT3{FT}})),
-                bottom = Operators.SetValue(zero(UpwindMatrixRowType{CT3{FT}})),
-            ) # Need to wrap ᶠupwind_matrix in this for well-defined boundaries.
-            ᶜkappa_mʲ = p.ᶜtemp_scalar
-            @. ᶜkappa_mʲ =
-                TD.gas_constant_air(thermo_params, ᶜtsʲs.:(1)) /
-                TD.cv_m(thermo_params, ᶜtsʲs.:(1))
-            ∂ᶜq_totʲ_err_∂ᶜq_totʲ =
-                matrix[@name(c.sgsʲs.:(1).q_tot), @name(c.sgsʲs.:(1).q_tot)]
-            @. ∂ᶜq_totʲ_err_∂ᶜq_totʲ =
-                dtγ * (
-                    DiagonalMatrixRow(ᶜadvdivᵥ(ᶠu³ʲs.:(1))) -
-                    ᶜadvdivᵥ_matrix() ⋅
-                    ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1)))
-                ) - (I,)
-
-            ∂ᶜmseʲ_err_∂ᶜq_totʲ =
-                matrix[@name(c.sgsʲs.:(1).mse), @name(c.sgsʲs.:(1).q_tot)]
-            @. ∂ᶜmseʲ_err_∂ᶜq_totʲ =
-                dtγ * (
-                    -DiagonalMatrixRow(
-                        adjoint(ᶜinterp(ᶠu³ʲs.:(1))) *
-                        ᶜgradᵥ_ᶠΦ *
-                        Y.c.ρ *
-                        ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp) * e_int_v0,
-                    )
-                )
-            ∂ᶜmseʲ_err_∂ᶜρ = matrix[@name(c.sgsʲs.:(1).mse), @name(c.ρ)]
-            @. ∂ᶜmseʲ_err_∂ᶜρ =
-                dtγ * (
-                    -DiagonalMatrixRow(
-                        adjoint(ᶜinterp(ᶠu³ʲs.:(1))) * ᶜgradᵥ_ᶠΦ / ᶜρʲs.:(1),
-                    )
-                )
-            ∂ᶜmseʲ_err_∂ᶜmseʲ =
-                matrix[@name(c.sgsʲs.:(1).mse), @name(c.sgsʲs.:(1).mse)]
-            @. ∂ᶜmseʲ_err_∂ᶜmseʲ =
-                dtγ * (
-                    DiagonalMatrixRow(ᶜadvdivᵥ(ᶠu³ʲs.:(1))) -
-                    ᶜadvdivᵥ_matrix() ⋅
-                    ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) -
-                    DiagonalMatrixRow(
-                        adjoint(ᶜinterp(ᶠu³ʲs.:(1))) *
-                        ᶜgradᵥ_ᶠΦ *
-                        Y.c.ρ *
-                        ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp),
-                    )
-                ) - (I,)
-
-            ∂ᶜρaʲ_err_∂ᶜq_totʲ =
-                matrix[@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).q_tot)]
-            @. ᶠbidiagonal_matrix_ct3 =
-                DiagonalMatrixRow(
-                    ᶠset_upwind_bcs(
-                        ᶠupwind(
-                            ᶠu³ʲs.:(1),
-                            draft_area(Y.c.sgsʲs.:(1).ρa, ᶜρʲs.:(1)),
-                        ),
-                    ),
-                ) ⋅ ᶠwinterp_matrix(ᶜJ) ⋅ DiagonalMatrixRow(
-                    ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp) *
-                    e_int_v0,
-                )
-            @. ᶠbidiagonal_matrix_ct3_2 =
-                DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρʲs.:(1))) ⋅
-                ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) ⋅
-                DiagonalMatrixRow(
-                    Y.c.sgsʲs.:(1).ρa * ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp) *
-                    e_int_v0,
-                )
-
-            @. ∂ᶜρaʲ_err_∂ᶜq_totʲ =
-                dtγ * ᶜadvdivᵥ_matrix() ⋅
-                (ᶠbidiagonal_matrix_ct3 - ᶠbidiagonal_matrix_ct3_2)
-            ∂ᶜρaʲ_err_∂ᶜmseʲ =
-                matrix[@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).mse)]
-            @. ᶠbidiagonal_matrix_ct3 =
-                DiagonalMatrixRow(
-                    ᶠset_upwind_bcs(
-                        ᶠupwind(
-                            ᶠu³ʲs.:(1),
-                            draft_area(Y.c.sgsʲs.:(1).ρa, ᶜρʲs.:(1)),
-                        ),
-                    ),
-                ) ⋅ ᶠwinterp_matrix(ᶜJ) ⋅ DiagonalMatrixRow(
-                    ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp),
-                )
-            @. ᶠbidiagonal_matrix_ct3_2 =
-                DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρʲs.:(1))) ⋅
-                ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) ⋅
-                DiagonalMatrixRow(
-                    Y.c.sgsʲs.:(1).ρa * ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp),
-                )
-            @. ∂ᶜρaʲ_err_∂ᶜmseʲ =
-                dtγ * ᶜadvdivᵥ_matrix() ⋅
-                (ᶠbidiagonal_matrix_ct3 - ᶠbidiagonal_matrix_ct3_2)
-            ∂ᶜρaʲ_err_∂ᶜρaʲ =
-                matrix[@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).ρa)]
-            @. ᶜadvection_matrix =
-                -(ᶜadvdivᵥ_matrix()) ⋅
-                DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρʲs.:(1)))
-            @. ∂ᶜρaʲ_err_∂ᶜρaʲ =
-                dtγ * ᶜadvection_matrix ⋅
-                ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) ⋅
-                DiagonalMatrixRow(1 / ᶜρʲs.:(1)) - (I,)
-
-            ∂ᶠu₃ʲ_err_∂ᶜρ = matrix[@name(f.sgsʲs.:(1).u₃), @name(c.ρ)]
-            @. ∂ᶠu₃ʲ_err_∂ᶜρ =
-                dtγ * DiagonalMatrixRow(ᶠgradᵥ_ᶜΦ / ᶠinterp(ᶜρʲs.:(1))) ⋅
-                ᶠinterp_matrix()
-            ∂ᶠu₃ʲ_err_∂ᶜq_totʲ =
-                matrix[@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).q_tot)]
-            @. ∂ᶠu₃ʲ_err_∂ᶜq_totʲ =
-                dtγ * DiagonalMatrixRow(
-                    ᶠgradᵥ_ᶜΦ * ᶠinterp(Y.c.ρ) / (ᶠinterp(ᶜρʲs.:(1)))^2,
-                ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow(
-                    ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp) *
-                    e_int_v0,
-                )
-            ∂ᶠu₃ʲ_err_∂ᶜmseʲ =
-                matrix[@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).mse)]
-            @. ∂ᶠu₃ʲ_err_∂ᶜmseʲ =
-                dtγ * DiagonalMatrixRow(
-                    ᶠgradᵥ_ᶜΦ * ᶠinterp(Y.c.ρ) / (ᶠinterp(ᶜρʲs.:(1)))^2,
-                ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow(
-                    ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp),
-                )
-            ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ =
-                matrix[@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)]
-            ᶜu₃ʲ = ᶜtemp_C3
-            @. ᶜu₃ʲ = ᶜinterp(Y.f.sgsʲs.:(1).u₃)
-
-            @. bdmr_l = convert(BidiagonalMatrixRow{FT}, ᶜleft_bias_matrix())
-            @. bdmr_r = convert(BidiagonalMatrixRow{FT}, ᶜright_bias_matrix())
-            @. bdmr = ifelse(ᶜu₃ʲ.components.data.:1 > 0, bdmr_l, bdmr_r)
-
-            @. ᶠtridiagonal_matrix_c3 = -(ᶠgradᵥ_matrix()) ⋅ bdmr
-            if p.atmos.rayleigh_sponge isa RayleighSponge
-                @. ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ =
-                    dtγ * (
-                        ᶠtridiagonal_matrix_c3 ⋅
-                        DiagonalMatrixRow(adjoint(CT3(Y.f.sgsʲs.:(1).u₃))) -
-                        DiagonalMatrixRow(p.ᶠβ_rayleigh_w * (one_C3xACT3,))
-                    ) - (I_u₃,)
-            else
-                @. ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ =
-                    dtγ * ᶠtridiagonal_matrix_c3 ⋅
-                    DiagonalMatrixRow(adjoint(CT3(Y.f.sgsʲs.:(1).u₃))) - (I_u₃,)
-            end
-        elseif p.atmos.rayleigh_sponge isa RayleighSponge
-            ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ =
-                matrix[@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)]
-            @. ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ =
-                dtγ * -DiagonalMatrixRow(p.ᶠβ_rayleigh_w * (one_C3xACT3,)) -
-                (I_u₃,)
-        end
-    end
-
+    A.cache.b_krylov .= b
+    ldiv!(A.cache.x_krylov, A, A.cache.b_krylov)
+    x .= A.cache.x_krylov
 end
diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl
index 20f437a9119..5562303ee21 100644
--- a/src/solver/type_getters.jl
+++ b/src/solver/type_getters.jl
@@ -9,8 +9,8 @@ import ClimaAtmos as CA
 import LinearAlgebra
 import ClimaCore.Fields
 import ClimaTimeSteppers as CTS
-
 import ClimaDiagnostics
+import DiffEqBase: KeywordArgSilent
 
 function get_atmos(config::AtmosConfig, params)
     (; turbconv_params) = params
@@ -369,56 +369,29 @@ function get_surface_setup(parsed_args)
     return getproperty(SurfaceConditions, Symbol(parsed_args["surface_setup"]))()
 end
 
-is_explicit_CTS_algo_type(alg_or_tableau) =
-    alg_or_tableau <: CTS.ERKAlgorithmName
-
-is_imex_CTS_algo_type(alg_or_tableau) =
-    alg_or_tableau <: CTS.IMEXARKAlgorithmName
-
-is_implicit_type(alg_or_tableau) = is_imex_CTS_algo_type(alg_or_tableau)
-
-is_imex_CTS_algo(::CTS.IMEXAlgorithm) = true
-is_imex_CTS_algo(::CTS.RosenbrockAlgorithm) = true
-is_imex_CTS_algo(::SciMLBase.AbstractODEAlgorithm) = false
-
-is_implicit(ode_algo) = is_imex_CTS_algo(ode_algo)
-
-is_rosenbrock(::SciMLBase.AbstractODEAlgorithm) = false
-use_transform(ode_algo) =
-    !(is_imex_CTS_algo(ode_algo) || is_rosenbrock(ode_algo))
-
-additional_integrator_kwargs(::SciMLBase.AbstractODEAlgorithm) = (;
-    adaptive = false,
-    progress = isinteractive(),
-    progress_steps = isinteractive() ? 1 : 1000,
-)
-import DiffEqBase
-additional_integrator_kwargs(::CTS.DistributedODEAlgorithm) = (;
-    kwargshandle = DiffEqBase.KeywordArgSilent, # allow custom kwargs
-    adjustfinal = true,
-    # TODO: enable progress bars in ClimaTimeSteppers
-)
-
-is_cts_algo(::SciMLBase.AbstractODEAlgorithm) = false
-is_cts_algo(::CTS.DistributedODEAlgorithm) = true
-
-function jac_kwargs(ode_algo, Y, atmos, parsed_args)
-    if is_implicit(ode_algo)
-        A = ImplicitEquationJacobian(
-            Y,
-            atmos;
-            approximate_solve_iters = parsed_args["approximate_linear_solve_iters"],
-            transform_flag = use_transform(ode_algo),
-        )
-        if use_transform(ode_algo)
-            return (; jac_prototype = A, Wfact_t = Wfact!)
+jac_kwargs(ode_algo, Y, p, parsed_args) =
+    if ode_algo isa Union{CTS.IMEXAlgorithm, CTS.RosenbrockAlgorithm}
+        use_exact_jacobian = parsed_args["use_exact_jacobian"]
+        always_update_exact_jacobian =
+            parsed_args["n_steps_update_exact_jacobian"] == 0
+        approximate_solve_iters = parsed_args["approximate_linear_solve_iters"]
+        jacobian_algorithm = if parsed_args["debug_approximate_jacobian"]
+            DebugJacobian(
+                ApproxJacobian(; approximate_solve_iters);
+                use_exact_jacobian,
+                always_update_exact_jacobian,
+            )
         else
-            return (; jac_prototype = A, Wfact = Wfact!)
+            use_exact_jacobian ?
+            ExactJacobian(; always_update_exact_jacobian) :
+            ApproxJacobian(; approximate_solve_iters)
         end
+        @info "Jacobian algorithm: $(dump_string(jacobian_algorithm))"
+        A = ImplicitEquationJacobian(jacobian_algorithm, Y, p)
+        (; jac_prototype = A, Wfact = update_jacobian!)
     else
-        return NamedTuple()
+        (;)
     end
-end
 
 #=
     ode_configuration(Y, parsed_args)
@@ -427,17 +400,14 @@ Returns the ode algorithm
 =#
 function ode_configuration(::Type{FT}, parsed_args) where {FT}
     ode_name = parsed_args["ode_algo"]
-    alg_or_tableau = getproperty(CTS, Symbol(ode_name))
-    @info "Using ODE config: `$alg_or_tableau`"
-    if ode_name == "SSPKnoth"
-        return CTS.RosenbrockAlgorithm(CTS.tableau(CTS.SSPKnoth()))
-    end
-
-    if is_explicit_CTS_algo_type(alg_or_tableau)
-        return CTS.ExplicitAlgorithm(alg_or_tableau())
-    elseif !is_implicit_type(alg_or_tableau)
-        return alg_or_tableau()
-    elseif is_imex_CTS_algo_type(alg_or_tableau)
+    ode_algo_name = getproperty(CTS, Symbol(ode_name))
+    @info "Using ODE config: `$ode_algo_name`"
+    return if ode_algo_name <: CTS.RosenbrockAlgorithmName
+        CTS.RosenbrockAlgorithm(CTS.tableau(ode_algo_name()))
+    elseif ode_algo_name <: CTS.ERKAlgorithmName
+        CTS.ExplicitAlgorithm(ode_algo_name())
+    else
+        @assert ode_algo_name <: CTS.IMEXARKAlgorithmName
         newtons_method = CTS.NewtonsMethod(;
             max_iters = parsed_args["max_newton_iters_ode"],
             krylov_method = if parsed_args["use_krylov_method"]
@@ -466,9 +436,7 @@ function ode_configuration(::Type{FT}, parsed_args) where {FT}
                 nothing
             end,
         )
-        return CTS.IMEXAlgorithm(alg_or_tableau(), newtons_method)
-    else
-        return alg_or_tableau(; linsolve = linsolve!)
+        CTS.IMEXAlgorithm(ode_algo_name(), newtons_method)
     end
 end
 
@@ -522,29 +490,20 @@ function get_sim_info(config::AtmosConfig)
 end
 
 function args_integrator(parsed_args, Y, p, tspan, ode_algo, callback)
-    (; atmos, dt) = p
+    (; dt) = p
     dt_save_to_sol = time_to_seconds(parsed_args["dt_save_to_sol"])
 
     s = @timed_str begin
         func = if parsed_args["split_ode"]
-            implicit_func = SciMLBase.ODEFunction(
-                implicit_tendency!;
-                jac_kwargs(ode_algo, Y, atmos, parsed_args)...,
-                tgrad = (∂Y∂t, Y, p, t) -> (∂Y∂t .= 0),
-            )
-            if is_cts_algo(ode_algo)
-                CTS.ClimaODEFunction(;
-                    T_exp_T_lim! = remaining_tendency!,
-                    T_imp! = implicit_func,
-                    # Can we just pass implicit_tendency! and jac_prototype etc.?
-                    lim! = limiters_func!,
-                    dss!,
-                    post_explicit! = set_precomputed_quantities!,
-                    post_implicit! = set_precomputed_quantities!,
-                )
-            else
-                SciMLBase.SplitFunction(implicit_func, remaining_tendency!)
-            end
+            jacobian = jac_kwargs(ode_algo, Y, p, parsed_args)
+            CTS.ClimaODEFunction(;
+                T_exp_T_lim! = remaining_tendency!,
+                T_imp! = SciMLBase.ODEFunction(implicit_tendency!; jacobian...),
+                lim! = limiters_func!,
+                dss!,
+                post_explicit! = set_precomputed_quantities!,
+                post_implicit! = set_precomputed_quantities!,
+            ) # TODO: Separate the Wfact! and jac_prototype kwargs from T_imp!.
         else
             remaining_tendency! # should be total_tendency!
         end
@@ -560,7 +519,13 @@ function args_integrator(parsed_args, Y, p, tspan, ode_algo, callback)
     end # ensure that tspan[2] is always saved
     @info "dt_save_to_sol: $dt_save_to_sol, length(saveat): $(length(saveat))"
     args = (problem, ode_algo)
-    kwargs = (; saveat, callback, dt, additional_integrator_kwargs(ode_algo)...)
+    kwargs = (;
+        saveat,
+        callback,
+        dt,
+        kwargshandle = KeywordArgSilent, # allow custom kwargs
+        adjustfinal = true,
+    )
     return (args, kwargs)
 end
 
@@ -719,6 +684,13 @@ function get_simulation(config::AtmosConfig)
     end
     @info "init integrator: $s"
 
+    if (
+        config.parsed_args["use_exact_jacobian"] ||
+        config.parsed_args["debug_approximate_jacobian"]
+    ) && config.parsed_args["n_steps_update_exact_jacobian"] == 0
+        update_exact_jacobian!(integrator)
+    end # TODO: Is this necessary for initialization?
+
     if config.parsed_args["enable_diagnostics"]
         s = @timed_str begin
             integrator = ClimaDiagnostics.IntegratorWithDiagnostics(
diff --git a/src/solver/types.jl b/src/solver/types.jl
index 81a6deacdf3..1e54fa9062a 100644
--- a/src/solver/types.jl
+++ b/src/solver/types.jl
@@ -175,6 +175,7 @@ struct MixingLength{FT}
     tke::FT
     buoy::FT
 end
+MixingLength(args...) = MixingLength{promote_type(typeof.(args)...)}(args...)
 
 abstract type AbstractEDMF end
 
diff --git a/src/utils/abbreviations.jl b/src/utils/abbreviations.jl
index 8b686bc2e9d..6b5dea700df 100644
--- a/src/utils/abbreviations.jl
+++ b/src/utils/abbreviations.jl
@@ -1,4 +1,5 @@
-using ClimaCore: Geometry, Operators, MatrixFields
+using ClimaCore: Geometry, Operators
+using ClimaCore.MatrixFields
 
 # Alternatively, we could use Vec₁₂₃, Vec³, etc., if that is more readable.
 const C1 = Geometry.Covariant1Vector
diff --git a/src/utils/utilities.jl b/src/utils/utilities.jl
index c3ac7a4d4a2..79f82fcb2b1 100644
--- a/src/utils/utilities.jl
+++ b/src/utils/utilities.jl
@@ -4,6 +4,7 @@
 import ClimaComms
 import ClimaCore: Spaces, Topologies, Fields, Geometry
 import LinearAlgebra: norm_sqr
+import Statistics: quantile
 
 is_energy_var(symbol) = symbol in (:ρe_tot, :ρae_tot)
 is_momentum_var(symbol) = symbol in (:uₕ, :ρuₕ, :u₃, :ρw)
@@ -283,28 +284,40 @@ end
 import Dates
 
 """
-    time_and_units_str(x::Real)
-
-Returns a truncated string of time and units,
-given a time `x` in Seconds.
-"""
-time_and_units_str(x::Real) =
-    trunc_time(string(compound_period(x, Dates.Second)))
+    time_and_units_str(seconds)
+
+Converts the given number of seconds into a string that contains a numerical
+value (rounded to the nearest 2 decimal digits) and its corresponding units.
+"""
+function time_and_units_str(seconds)
+    # Convert seconds from Float32 to Float64 to avoid round-off errors.
+    precise_seconds = seconds isa Float32 ? Float64(seconds) : seconds
+    full_period = compound_period(precise_seconds)
+    isempty(Dates.periods(full_period)) && return "0 Seconds"
+    whole_period = Dates.periods(full_period)[1]
+    whole_period_value = Dates.value(whole_period)
+    remaining_period = full_period - whole_period
+    no_remaining_digits = isempty(Dates.periods(remaining_period))
+    period_units =
+        string(typeof(whole_period).name.name) *
+        (whole_period_value == 1 && no_remaining_digits ? "" : "s")
+    no_remaining_digits && return "$whole_period_value $period_units"
+    remaining_period_value =
+        Dates.value(convert(Dates.Nanosecond, remaining_period)) /
+        Dates.value(convert(Dates.Nanosecond, typeof(whole_period)(1)))
+    remaining_period_value < 0.01 && return "$whole_period_value $period_units"
+    remaining_digits = lpad(round(Int, remaining_period_value * 100), 2, '0')
+    return "$whole_period_value.$remaining_digits $period_units"
+end
 
 """
-    compound_period(x::Real, ::Type{T}) where {T <: Dates.Period}
+    compound_period(seconds)
 
-A canonicalized `Dates.CompoundPeriod` given a real value
-`x`, and its units via the period type `T`.
+A `Dates.CompoundPeriod` that represents the given number of seconds (rounded to
+the nearest nanosecond).
 """
-function compound_period(x::Real, ::Type{T}) where {T <: Dates.Period}
-    nf = Dates.value(convert(Dates.Nanosecond, T(1)))
-    ns = Dates.Nanosecond(ceil(x * nf))
-    return Dates.canonicalize(Dates.CompoundPeriod(ns))
-end
-
-trunc_time(s::String) = count(',', s) > 1 ? join(split(s, ",")[1:2], ",") : s
-
+compound_period(seconds) =
+    Dates.canonicalize(Dates.Nanosecond(round(Int, seconds * 10^9)))
 
 function prettymemory(b)
     if b < 1024
@@ -332,6 +345,19 @@ macro timed_str(ex)
     end
 end
 
+"""
+    dump_string(x)
+
+Returns a string that contains the output of `dump(x)`.
+"""
+function dump_string(x)
+    buffer = IOBuffer()
+    dump(buffer, x)
+    result = String(take!(buffer))
+    close(buffer)
+    return result
+end
+
 struct AllNothing end
 const all_nothing = AllNothing()
 Base.getproperty(::AllNothing, ::Symbol) = nothing
@@ -352,6 +378,19 @@ function horizontal_integral_at_boundary(f::Fields.Field)
     sum(f ./ Fields.Δz_field(axes(f)) .* 2) # TODO: is there a way to ensure this is derived from face z? The 2d topology doesn't contain this info
 end
 
+"""
+    smooth_maximum(array, [percentile_rank])
+
+Filters out all non-negligible values in the given `array` and computes the
+specified percentile of the result.
+"""
+function smooth_maximum(array, percentile_rank = 1 - 0.1^ndims(array))
+    minimum_filtered_value = maximum(abs, array) / 1e6
+    filtered_array_values = filter(>(minimum_filtered_value), map(abs, array))
+    isempty(filtered_array_values) && return zero(eltype(array))
+    return quantile(filtered_array_values, percentile_rank)
+end
+
 """
     gaussian_smooth(arr, sigma = 1)
 
diff --git a/src/utils/variable_manipulations.jl b/src/utils/variable_manipulations.jl
index 721bc78b0e4..9a2dcd80360 100644
--- a/src/utils/variable_manipulations.jl
+++ b/src/utils/variable_manipulations.jl
@@ -62,9 +62,9 @@ function as
        (1 + tanh(2 * atanh(1 - 2 * (1 - a)^(-1 / log2(1 - a_half))))) / 2`.
 """
 sgs_weight_function(a, a_half) =
-    if a < 0
+    if a <= 0 # autodiff generates NaNs when a is 0
         zero(a)
-    elseif a > 1
+    elseif a > min(1, 42 * a_half) # autodiff generates NaNs when a is large
         one(a)
     else
         (1 + tanh(2 * atanh(1 - 2 * (1 - a)^(-1 / log2(1 - a_half))))) / 2
diff --git a/test/dependencies.jl b/test/dependencies.jl
index 92e507ee529..267fac7c684 100644
--- a/test/dependencies.jl
+++ b/test/dependencies.jl
@@ -38,6 +38,7 @@ known_dependencies = Set([
     "ArtifactWrappers",
     "Artifacts",
     "AtmosphericProfilesLibrary",
+    "BlockArrays",
     "ClimaComms",
     "ClimaCore",
     "ClimaDiagnostics",
@@ -48,6 +49,7 @@ known_dependencies = Set([
     "Dates",
     "DiffEqBase",
     "FastGaussQuadrature",
+    "ForwardDiff",
     "Insolation",
     "Interpolations",
     "LazyArtifacts",