diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index 2978c58ed..a47d51a1d 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -1,25 +1,19 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.2"
+julia_version = "1.10.3"
 manifest_format = "2.0"
-project_hash = "6c87fe520a2f2573c1d64ad2168076eae19c6786"
+project_hash = "a503b4f951cf10fc76e1520a620cb8fdc7c63464"
 
 [[deps.ADTypes]]
-git-tree-sha1 = "fcdb00b4d412b80ab08e39978e3bdef579e5e224"
+git-tree-sha1 = "daf26bbdec60d9ca1c0003b70f389d821ddb4224"
 uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
-version = "1.0.0"
+version = "1.2.1"
 weakdeps = ["ChainRulesCore", "EnzymeCore"]
 
     [deps.ADTypes.extensions]
     ADTypesChainRulesCoreExt = "ChainRulesCore"
     ADTypesEnzymeCoreExt = "EnzymeCore"
 
-[[deps.AMD]]
-deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse_jll"]
-git-tree-sha1 = "45a1272e3f809d36431e57ab22703c6896b8908f"
-uuid = "14f7f29c-3bd6-536c-9a0b-7339e30b5a3e"
-version = "0.5.3"
-
 [[deps.ANSIColoredPrinters]]
 git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
 uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9"
@@ -73,10 +67,10 @@ weakdeps = ["StaticArrays"]
     AdaptStaticArraysExt = "StaticArrays"
 
 [[deps.AliasTables]]
-deps = ["Random"]
-git-tree-sha1 = "07591db28451b3e45f4c0088a2d5e986ae5aa92d"
+deps = ["PtrArrays", "Random"]
+git-tree-sha1 = "9876e1e164b144ca45e9e3198d0b689cadfed9ff"
 uuid = "66dad0bd-aa9a-41b7-9441-69ab47430ed8"
-version = "1.1.1"
+version = "1.1.3"
 
 [[deps.ArgTools]]
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
@@ -112,9 +106,9 @@ version = "7.10.0"
 
 [[deps.ArrayLayouts]]
 deps = ["FillArrays", "LinearAlgebra"]
-git-tree-sha1 = "33207a8be6267bc389d0701e97a9bce6a4de68eb"
+git-tree-sha1 = "29649b61e0313db0a7ad5ecf41210e4e85aea234"
 uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
-version = "1.9.2"
+version = "1.9.3"
 weakdeps = ["SparseArrays"]
 
     [deps.ArrayLayouts.extensions]
@@ -129,12 +123,6 @@ git-tree-sha1 = "c06a868224ecba914baa6942988e2f2aade419be"
 uuid = "a9b6321e-bd34-4604-b9c9-b65b8de01458"
 version = "0.1.0"
 
-[[deps.BFloat16s]]
-deps = ["LinearAlgebra", "Printf", "Random", "Test"]
-git-tree-sha1 = "2c7cc21e8678eff479978a0a2ef5ce2f51b63dff"
-uuid = "ab4f0b2a-ad5b-11e8-123f-65d77653426b"
-version = "0.5.0"
-
 [[deps.BandedMatrices]]
 deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "PrecompileTools"]
 git-tree-sha1 = "30b7ea34abc4fe816eb1a5f434a43da804836163"
@@ -195,44 +183,15 @@ version = "0.5.0"
 
 [[deps.CPUSummary]]
 deps = ["CpuId", "IfElse", "PrecompileTools", "Static"]
-git-tree-sha1 = "601f7e7b3d36f18790e2caf83a882d88e9b71ff1"
+git-tree-sha1 = "585a387a490f1c4bd88be67eea15b93da5e85db7"
 uuid = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9"
-version = "0.2.4"
-
-[[deps.CUDA]]
-deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CUDA_Driver_jll", "CUDA_Runtime_Discovery", "CUDA_Runtime_jll", "Crayons", "DataFrames", "ExprTools", "GPUArrays", "GPUCompiler", "KernelAbstractions", "LLVM", "LLVMLoopInfo", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "NVTX", "Preferences", "PrettyTables", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "StaticArrays", "Statistics"]
-git-tree-sha1 = "4e33522a036b39fc6f5cb7447ae3b28eb8fbe99b"
-uuid = "052768ef-5323-5732-b1bb-66c8b64840ba"
-version = "5.3.3"
-weakdeps = ["ChainRulesCore", "SpecialFunctions"]
-
-    [deps.CUDA.extensions]
-    ChainRulesCoreExt = "ChainRulesCore"
-    SpecialFunctionsExt = "SpecialFunctions"
-
-[[deps.CUDA_Driver_jll]]
-deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
-git-tree-sha1 = "dc172b558adbf17952001e15cf0d6364e6d78c2f"
-uuid = "4ee394cb-3365-5eb0-8335-949819d2adfc"
-version = "0.8.1+0"
-
-[[deps.CUDA_Runtime_Discovery]]
-deps = ["Libdl"]
-git-tree-sha1 = "38f830504358e9972d2a0c3e5d51cb865e0733df"
-uuid = "1af6417a-86b4-443c-805f-a4643ffb695f"
-version = "0.2.4"
-
-[[deps.CUDA_Runtime_jll]]
-deps = ["Artifacts", "CUDA_Driver_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"]
-git-tree-sha1 = "4ca7d6d92075906c2ce871ea8bba971fff20d00c"
-uuid = "76a88914-d11a-5bdc-97e0-2f5a05c973a2"
-version = "0.12.1+0"
+version = "0.2.5"
 
 [[deps.Cairo_jll]]
 deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
-git-tree-sha1 = "a4c43f59baa34011e303e76f5c8c91bf58415aaf"
+git-tree-sha1 = "a2f1c8c668c8e3cb4cca4e57a8efdb09067bb3fd"
 uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
-version = "1.18.0+1"
+version = "1.18.0+2"
 
 [[deps.Calculus]]
 deps = ["LinearAlgebra"]
@@ -251,32 +210,53 @@ weakdeps = ["SparseArrays"]
     ChainRulesCoreSparseArraysExt = "SparseArrays"
 
 [[deps.ClimaComms]]
-deps = ["CUDA", "MPI"]
-git-tree-sha1 = "ef5d206be51fdf62cd0cbd63058e237128652cf7"
+git-tree-sha1 = "f19a9d42ef27affa8d57a5702dcb8a318077aa86"
 uuid = "3a4d1b5c-c61d-41fd-a00a-5873ba7a1b0d"
-version = "0.5.8"
+version = "0.6.0"
+
+    [deps.ClimaComms.extensions]
+    ClimaCommsCUDAExt = "CUDA"
+    ClimaCommsMPIExt = "MPI"
+
+    [deps.ClimaComms.weakdeps]
+    CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+    MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 
 [[deps.ClimaCore]]
-deps = ["Adapt", "BandedMatrices", "BlockArrays", "CUDA", "ClimaComms", "CubedSphere", "DataStructures", "DocStringExtensions", "ForwardDiff", "GaussQuadrature", "GilbertCurves", "HDF5", "InteractiveUtils", "IntervalSets", "KrylovKit", "LinearAlgebra", "MultiBroadcastFusion", "PkgVersion", "RecursiveArrayTools", "RootSolvers", "SparseArrays", "Static", "StaticArrays", "Statistics", "Unrolled"]
-git-tree-sha1 = "ce60e793254676e4ff0121997bae6fc4b78ebd09"
+deps = ["Adapt", "BandedMatrices", "BlockArrays", "ClimaComms", "CubedSphere", "DataStructures", "DocStringExtensions", "ForwardDiff", "GaussQuadrature", "GilbertCurves", "HDF5", "InteractiveUtils", "IntervalSets", "KrylovKit", "LinearAlgebra", "MultiBroadcastFusion", "NVTX", "PkgVersion", "RecursiveArrayTools", "RootSolvers", "SparseArrays", "Static", "StaticArrays", "Statistics", "Unrolled"]
+git-tree-sha1 = "684a90d619335c8160cb0cfae7c602701649e4e9"
 uuid = "d414da3d-4745-48bb-8d80-42e94e092884"
-version = "0.13.4"
-weakdeps = ["Krylov"]
+version = "0.14.5"
 
     [deps.ClimaCore.extensions]
+    ClimaCoreCUDAExt = "CUDA"
     KrylovExt = "Krylov"
 
+    [deps.ClimaCore.weakdeps]
+    CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+    Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
+
 [[deps.ClimaCorePlots]]
 deps = ["ClimaCore", "RecipesBase", "StaticArrays", "TriplotBase"]
-git-tree-sha1 = "ded3e0f3e7069f7c807f7b56caff232921bc2f5f"
+git-tree-sha1 = "4fc0c2a0845fc34a8d6c0e0a276d62a0d4e50286"
 uuid = "cf7c7e5a-b407-4c48-9047-11a94a308626"
-version = "0.2.8"
+version = "0.2.10"
 
 [[deps.ClimaTimeSteppers]]
 deps = ["ClimaComms", "Colors", "DataStructures", "DiffEqBase", "DiffEqCallbacks", "KernelAbstractions", "Krylov", "LinearAlgebra", "LinearOperators", "NVTX", "SciMLBase", "StaticArrays"]
 path = ".."
 uuid = "595c0a79-7f3d-439a-bc5a-b232dc3bde79"
-version = "0.7.20"
+version = "0.7.30"
+
+    [deps.ClimaTimeSteppers.extensions]
+    ClimaTimeSteppersBenchmarkToolsExt = ["CUDA", "BenchmarkTools", "OrderedCollections", "StatsBase", "PrettyTables"]
+
+    [deps.ClimaTimeSteppers.weakdeps]
+    BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+    CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+    OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+    PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+    StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [[deps.CloseOpenIntervals]]
 deps = ["Static", "StaticArrayInterface"]
@@ -292,9 +272,9 @@ version = "0.7.4"
 
 [[deps.ColorSchemes]]
 deps = ["ColorTypes", "ColorVectorSpace", "Colors", "FixedPointNumbers", "PrecompileTools", "Random"]
-git-tree-sha1 = "67c1f244b991cad9b0aa4b7540fb758c2488b129"
+git-tree-sha1 = "4b270d6465eb21ae89b732182c20dc165f8bf9f2"
 uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
-version = "3.24.0"
+version = "3.25.0"
 
 [[deps.ColorTypes]]
 deps = ["FixedPointNumbers", "Random"]
@@ -314,9 +294,9 @@ weakdeps = ["SpecialFunctions"]
 
 [[deps.Colors]]
 deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
-git-tree-sha1 = "fc08e5930ee9a4e03f84bfb5211cb54e7769758a"
+git-tree-sha1 = "362a287c3aa50601b0bc359053d5c2468f0e7ce0"
 uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
-version = "0.12.10"
+version = "0.12.11"
 
 [[deps.CommonSolve]]
 git-tree-sha1 = "0eee5eb66b1cf62cd6ad1b460238e60e4b09400c"
@@ -331,9 +311,9 @@ version = "0.3.0"
 
 [[deps.Compat]]
 deps = ["TOML", "UUIDs"]
-git-tree-sha1 = "c955881e3c981181362ae4088b35995446298b80"
+git-tree-sha1 = "b1c55339b7c6c350ee89f2c1604299660525b248"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.14.0"
+version = "4.15.0"
 weakdeps = ["Dates", "LinearAlgebra"]
 
     [deps.Compat.extensions]
@@ -342,7 +322,7 @@ weakdeps = ["Dates", "LinearAlgebra"]
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.1.0+0"
+version = "1.1.1+0"
 
 [[deps.CompositionsBase]]
 git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad"
@@ -402,12 +382,6 @@ git-tree-sha1 = "abe83f3a2f1b857aac70ef8b269080af17764bbe"
 uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
 version = "1.16.0"
 
-[[deps.DataFrames]]
-deps = ["Compat", "DataAPI", "DataStructures", "Future", "InlineStrings", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrecompileTools", "PrettyTables", "Printf", "REPL", "Random", "Reexport", "SentinelArrays", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
-git-tree-sha1 = "04c738083f29f86e62c8afc341f0967d8717bdb8"
-uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-version = "1.6.1"
-
 [[deps.DataStructures]]
 deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
 git-tree-sha1 = "1d0a14036acb104d9e89698bd408f63ab58cdc82"
@@ -431,11 +405,12 @@ version = "1.9.1"
 
 [[deps.DiffEqBase]]
 deps = ["ArrayInterface", "ConcreteStructs", "DataStructures", "DocStringExtensions", "EnumX", "EnzymeCore", "FastBroadcast", "FastClosures", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "Parameters", "PreallocationTools", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "Setfield", "SparseArrays", "Static", "StaticArraysCore", "Statistics", "Tricks", "TruncatedStacktraces"]
-git-tree-sha1 = "531c53fd0405716712a8b4960216c3b7b5ec89b9"
+git-tree-sha1 = "37d49a1f8eedfe68b7622075ff3ebe3dd0e8f327"
 uuid = "2b5f629d-d688-5b77-993f-72d75c75574e"
-version = "6.149.1"
+version = "6.151.2"
 
     [deps.DiffEqBase.extensions]
+    DiffEqBaseCUDAExt = "CUDA"
     DiffEqBaseChainRulesCoreExt = "ChainRulesCore"
     DiffEqBaseDistributionsExt = "Distributions"
     DiffEqBaseEnzymeExt = ["ChainRulesCore", "Enzyme"]
@@ -448,6 +423,7 @@ version = "6.149.1"
     DiffEqBaseUnitfulExt = "Unitful"
 
     [deps.DiffEqBase.weakdeps]
+    CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
     ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
     Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
     Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
@@ -498,9 +474,9 @@ uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 
 [[deps.Distributions]]
 deps = ["AliasTables", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns"]
-git-tree-sha1 = "22c595ca4146c07b16bcf9c8bea86f731f7109d2"
+git-tree-sha1 = "9c405847cc7ecda2dc921ccf18b47ca150d7317e"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
-version = "0.25.108"
+version = "0.25.109"
 
     [deps.Distributions.extensions]
     DistributionsChainRulesCoreExt = "ChainRulesCore"
@@ -520,9 +496,9 @@ version = "0.9.3"
 
 [[deps.Documenter]]
 deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "CodecZlib", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "TOML", "Test", "Unicode"]
-git-tree-sha1 = "f15a91e6e3919055efa4f206f942a73fedf5dfe6"
+git-tree-sha1 = "5461b2a67beb9089980e2f8f25145186b6d34f91"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "1.4.0"
+version = "1.4.1"
 
 [[deps.DocumenterCitations]]
 deps = ["AbstractTrees", "Bibliography", "Dates", "Documenter", "Logging", "Markdown", "MarkdownAST", "OrderedCollections", "Unicode"]
@@ -552,9 +528,9 @@ uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 version = "1.0.4"
 
 [[deps.EnzymeCore]]
-git-tree-sha1 = "18394bc78ac2814ff38fe5e0c9dc2cd171e2810c"
+git-tree-sha1 = "0910982db2490a20f81dc7db5d4bbea236c027b3"
 uuid = "f151be2c-9106-41f4-ab19-57ee4f262869"
-version = "0.7.2"
+version = "0.7.3"
 weakdeps = ["Adapt"]
 
     [deps.EnzymeCore.extensions]
@@ -574,9 +550,9 @@ version = "0.1.10"
 
 [[deps.Expat_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "4558ab818dcceaab612d1bb8c19cee87eda2b83c"
+git-tree-sha1 = "1c6317308b9dc757616f0b5cb379db10494443a7"
 uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
-version = "2.5.0+0"
+version = "2.6.2+0"
 
 [[deps.ExprTools]]
 git-tree-sha1 = "27415f162e6028e81c72b82ef756bf321213b6ec"
@@ -609,9 +585,9 @@ version = "3.3.10+0"
 
 [[deps.FastBroadcast]]
 deps = ["ArrayInterface", "LinearAlgebra", "Polyester", "Static", "StaticArrayInterface", "StrideArraysCore"]
-git-tree-sha1 = "a6e756a880fc419c8b41592010aebe6a5ce09136"
+git-tree-sha1 = "e17367f052035620d832499496080f792fa7ea47"
 uuid = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
-version = "0.2.8"
+version = "0.3.2"
 
 [[deps.FastClosures]]
 git-tree-sha1 = "acebe244d53ee1b461970f8910c235b259e772ef"
@@ -629,9 +605,9 @@ uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
 
 [[deps.FillArrays]]
 deps = ["LinearAlgebra"]
-git-tree-sha1 = "57f08d5665e76397e96b168f9acc12ab17c84a68"
+git-tree-sha1 = "0653c0a2396a6da5bc4766c43041ef5fd3efbe57"
 uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "1.10.2"
+version = "1.11.0"
 weakdeps = ["PDMats", "SparseArrays", "Statistics"]
 
     [deps.FillArrays.extensions]
@@ -657,15 +633,15 @@ version = "2.23.1"
 
 [[deps.FixedPointNumbers]]
 deps = ["Statistics"]
-git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
+git-tree-sha1 = "05882d6995ae5c12bb5f36dd2ed3f61c98cbb172"
 uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
-version = "0.8.4"
+version = "0.8.5"
 
 [[deps.Fontconfig_jll]]
-deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
-git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03"
+deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Zlib_jll"]
+git-tree-sha1 = "db16beca600632c95fc8aca29890d83788dd8b23"
 uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
-version = "2.13.93+0"
+version = "2.13.96+0"
 
 [[deps.Format]]
 git-tree-sha1 = "9c68794ef81b08086aeb32eeaf33531668d5f5fc"
@@ -684,15 +660,15 @@ weakdeps = ["StaticArrays"]
 
 [[deps.FreeType2_jll]]
 deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Zlib_jll"]
-git-tree-sha1 = "d8db6a5a2fe1381c1ea4ef2cab7c69c2de7f9ea0"
+git-tree-sha1 = "5c1d8ae0efc6c2e7b1fc502cbe25def8f661b7bc"
 uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
-version = "2.13.1+0"
+version = "2.13.2+0"
 
 [[deps.FriBidi_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "1ed150b39aebcc805c26b93a8d0122c940f64ce2"
 uuid = "559328eb-81f9-559d-9380-de523a88c83c"
-version = "1.0.10+0"
+version = "1.0.14+0"
 
 [[deps.FunctionWrappers]]
 git-tree-sha1 = "d62485945ce5ae9c0c48f124a84998d755bae00e"
@@ -721,35 +697,23 @@ git-tree-sha1 = "ff38ba61beff76b8f4acad8ab0c97ef73bb670cb"
 uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89"
 version = "3.3.9+0"
 
-[[deps.GPUArrays]]
-deps = ["Adapt", "GPUArraysCore", "LLVM", "LinearAlgebra", "Printf", "Random", "Reexport", "Serialization", "Statistics"]
-git-tree-sha1 = "68e8ff56a4a355a85d2784b94614491f8c900cde"
-uuid = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
-version = "10.1.0"
-
 [[deps.GPUArraysCore]]
 deps = ["Adapt"]
 git-tree-sha1 = "ec632f177c0d990e64d955ccc1b8c04c485a0950"
 uuid = "46192b85-c4d5-4398-a991-12ede77f4527"
 version = "0.1.6"
 
-[[deps.GPUCompiler]]
-deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "Scratch", "TimerOutputs", "UUIDs"]
-git-tree-sha1 = "1600477fba37c9fc067b9be21f5e8101f24a8865"
-uuid = "61eb1bfa-7361-4325-ad38-22787b887f55"
-version = "0.26.4"
-
 [[deps.GR]]
-deps = ["Artifacts", "Base64", "DelimitedFiles", "Downloads", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Preferences", "Printf", "Random", "Serialization", "Sockets", "TOML", "Tar", "Test", "UUIDs", "p7zip_jll"]
-git-tree-sha1 = "3437ade7073682993e092ca570ad68a2aba26983"
+deps = ["Artifacts", "Base64", "DelimitedFiles", "Downloads", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Preferences", "Printf", "Random", "Serialization", "Sockets", "TOML", "Tar", "Test", "p7zip_jll"]
+git-tree-sha1 = "ddda044ca260ee324c5fc07edb6d7cf3f0b9c350"
 uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
-version = "0.73.3"
+version = "0.73.5"
 
 [[deps.GR_jll]]
 deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "FreeType2_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Qt6Base_jll", "Zlib_jll", "libpng_jll"]
-git-tree-sha1 = "a96d5c713e6aa28c242b0d25c1347e258d6541ab"
+git-tree-sha1 = "278e5e0f820178e8a26df3184fcb2280717c79b1"
 uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
-version = "0.73.3+0"
+version = "0.73.5+0"
 
 [[deps.GaussQuadrature]]
 deps = ["SpecialFunctions"]
@@ -782,9 +746,9 @@ version = "2.44.0+2"
 
 [[deps.Glib_jll]]
 deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Zlib_jll"]
-git-tree-sha1 = "359a1ba2e320790ddbe4ee8b4d54a305c0ea2aff"
+git-tree-sha1 = "7c82e6a6cd34e9d935e9aa4051b66c6ff3af59ba"
 uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
-version = "2.80.0+0"
+version = "2.80.2+0"
 
 [[deps.Graphite2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -802,11 +766,13 @@ deps = ["Compat", "HDF5_jll", "Libdl", "MPIPreferences", "Mmap", "Preferences",
 git-tree-sha1 = "e856eef26cf5bf2b0f95f8f4fc37553c72c8641c"
 uuid = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 version = "0.17.2"
-weakdeps = ["MPI"]
 
     [deps.HDF5.extensions]
     MPIExt = "MPI"
 
+    [deps.HDF5.weakdeps]
+    MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
+
 [[deps.HDF5_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "LibCURL_jll", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "OpenSSL_jll", "TOML", "Zlib_jll", "libaec_jll"]
 git-tree-sha1 = "82a471768b513dc39e471540fdadc84ff80ff997"
@@ -815,9 +781,9 @@ version = "1.14.3+3"
 
 [[deps.HTTP]]
 deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"]
-git-tree-sha1 = "2c3ec1f90bb4a8f7beafb0cffea8a4c3f4e636ab"
+git-tree-sha1 = "d1d712be3164d61d1fb98e7ce9bcbc6cc06b45ed"
 uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
-version = "1.10.6"
+version = "1.10.8"
 
 [[deps.HarfBuzz_jll]]
 deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"]
@@ -848,12 +814,6 @@ git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
 uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 version = "0.1.1"
 
-[[deps.InlineStrings]]
-deps = ["Parsers"]
-git-tree-sha1 = "9cc2baf75c6d09f9da536ddf58eb2f29dedaf461"
-uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
-version = "1.4.0"
-
 [[deps.IntelOpenMP_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
 git-tree-sha1 = "be50fe8df3acbffa0274a744f1a99d29c45a57f4"
@@ -885,11 +845,6 @@ weakdeps = ["Dates"]
     [deps.InverseFunctions.extensions]
     DatesExt = "Dates"
 
-[[deps.InvertedIndices]]
-git-tree-sha1 = "0dc7b50b8d436461be01300fd8cd45aa0274b038"
-uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
-version = "1.3.0"
-
 [[deps.IrrationalConstants]]
 git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2"
 uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
@@ -938,9 +893,9 @@ version = "1.3.0"
 
 [[deps.JpegTurbo_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "3336abae9a713d2210bb57ab484b1e065edd7d23"
+git-tree-sha1 = "c84a835e1a09b289ffcd2271bf2a337bbdda6637"
 uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
-version = "3.0.2+0"
+version = "3.0.3+0"
 
 [[deps.JuliaNVTXCallbacks_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -950,9 +905,9 @@ version = "0.2.1+0"
 
 [[deps.KernelAbstractions]]
 deps = ["Adapt", "Atomix", "InteractiveUtils", "LinearAlgebra", "MacroTools", "PrecompileTools", "Requires", "SparseArrays", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"]
-git-tree-sha1 = "ed7167240f40e62d97c1f5f7735dea6de3cc5c49"
+git-tree-sha1 = "db02395e4c374030c53dc28f3c1d33dec35f7272"
 uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
-version = "0.9.18"
+version = "0.9.19"
 weakdeps = ["EnzymeCore"]
 
     [deps.KernelAbstractions.extensions]
@@ -960,9 +915,9 @@ weakdeps = ["EnzymeCore"]
 
 [[deps.Krylov]]
 deps = ["LinearAlgebra", "Printf", "SparseArrays"]
-git-tree-sha1 = "8a6837ec02fe5fb3def1abc907bb802ef11a0729"
+git-tree-sha1 = "267dad6b4b7b5d529c76d40ff48d33f7e94cb834"
 uuid = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
-version = "0.9.5"
+version = "0.9.6"
 
 [[deps.KrylovKit]]
 deps = ["ChainRulesCore", "GPUArraysCore", "LinearAlgebra", "Printf"]
@@ -971,16 +926,10 @@ uuid = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 version = "0.6.1"
 
 [[deps.LAME_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "170b660facf5df5de098d866564877e119141cbd"
 uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d"
-version = "3.100.1+0"
-
-[[deps.LDLFactorizations]]
-deps = ["AMD", "LinearAlgebra", "SparseArrays", "Test"]
-git-tree-sha1 = "70f582b446a1c3ad82cf87e62b878668beef9d13"
-uuid = "40e66cde-538c-5869-a4ad-c39174c6795b"
-version = "0.10.1"
+version = "3.100.2+0"
 
 [[deps.LERC_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -990,25 +939,22 @@ version = "3.0.0+1"
 
 [[deps.LLVM]]
 deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Preferences", "Printf", "Requires", "Unicode"]
-git-tree-sha1 = "839c82932db86740ae729779e610f07a1640be9a"
+git-tree-sha1 = "065c36f95709dd4a676dc6839a35d6fa6f192f24"
 uuid = "929cbde3-209d-540e-8aea-75f648917ca0"
-version = "6.6.3"
-weakdeps = ["BFloat16s"]
+version = "7.1.0"
 
     [deps.LLVM.extensions]
     BFloat16sExt = "BFloat16s"
 
+    [deps.LLVM.weakdeps]
+    BFloat16s = "ab4f0b2a-ad5b-11e8-123f-65d77653426b"
+
 [[deps.LLVMExtra_jll]]
 deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"]
 git-tree-sha1 = "88b916503aac4fb7f701bb625cd84ca5dd1677bc"
 uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab"
 version = "0.0.29+0"
 
-[[deps.LLVMLoopInfo]]
-git-tree-sha1 = "2e5c102cfc41f48ae4740c7eca7743cc7e7b75ea"
-uuid = "8b046642-f1f6-4319-8d3c-209ddc03c586"
-version = "1.0.0"
-
 [[deps.LLVMOpenMP_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
 git-tree-sha1 = "d986ce2d884d49126836ea94ed5bfb0f12679713"
@@ -1016,10 +962,10 @@ uuid = "1d63c593-3942-5779-bab2-d838dc0a180e"
 version = "15.0.7+0"
 
 [[deps.LZO_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "70c5da094887fd2cae843b8db33920bac4b6f07d"
 uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
-version = "2.10.1+0"
+version = "2.10.2+0"
 
 [[deps.LaTeXStrings]]
 git-tree-sha1 = "50901ebc375ed41dbf8058da26f9de442febbbec"
@@ -1101,10 +1047,10 @@ uuid = "7e76a0d4-f3c7-5321-8279-8d96eeed0f29"
 version = "1.6.0+0"
 
 [[deps.Libgpg_error_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "fbb1f2bef882392312feb1ede3615ddc1e9b99ed"
 uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
-version = "1.42.0+0"
+version = "1.49.0+0"
 
 [[deps.Libiconv_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -1114,9 +1060,9 @@ version = "1.17.0+0"
 
 [[deps.Libmount_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "4b683b19157282f50bfd5dcaa2efe5295814ea22"
+git-tree-sha1 = "0c4f9c4f1a50d8f35048fa0532dabbadf702f81e"
 uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
-version = "2.40.0+0"
+version = "2.40.1+0"
 
 [[deps.Libtiff_jll]]
 deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "LERC_jll", "Libdl", "XZ_jll", "Zlib_jll", "Zstd_jll"]
@@ -1126,9 +1072,9 @@ version = "4.5.1+1"
 
 [[deps.Libuuid_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "27fd5cc10be85658cacfe11bb81bee216af13eda"
+git-tree-sha1 = "5ee6203157c120d79034c748a2acba45b82b8807"
 uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
-version = "2.40.0+0"
+version = "2.40.1+0"
 
 [[deps.LineSearches]]
 deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"]
@@ -1141,14 +1087,26 @@ deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [[deps.LinearOperators]]
-deps = ["FastClosures", "LDLFactorizations", "LinearAlgebra", "Printf", "Requires", "SparseArrays", "TimerOutputs"]
-git-tree-sha1 = "f06df3a46255879cbccae1b5b6dcb16994c31be7"
+deps = ["FastClosures", "LinearAlgebra", "Printf", "Requires", "SparseArrays", "TimerOutputs"]
+git-tree-sha1 = "ae5d90280094348c32fda8bc8b5a88bb16514d43"
 uuid = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125"
-version = "2.7.0"
-weakdeps = ["ChainRulesCore"]
+version = "2.8.0"
 
     [deps.LinearOperators.extensions]
+    LinearOperatorsCUDAExt = "CUDA"
     LinearOperatorsChainRulesCoreExt = "ChainRulesCore"
+    LinearOperatorsLDLFactorizationsExt = "LDLFactorizations"
+
+    [deps.LinearOperators.weakdeps]
+    CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+    LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b"
+
+[[deps.Literate]]
+deps = ["Base64", "IOCapture", "JSON", "REPL"]
+git-tree-sha1 = "596df2daea9c27da81eee63ef2cf101baf10c24c"
+uuid = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+version = "2.18.0"
 
 [[deps.LogExpFunctions]]
 deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
@@ -1181,37 +1139,23 @@ git-tree-sha1 = "80b2833b56d466b3858d565adcd16a4a05f2089b"
 uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
 version = "2024.1.0+0"
 
-[[deps.MPI]]
-deps = ["Distributed", "DocStringExtensions", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "PkgVersion", "PrecompileTools", "Requires", "Serialization", "Sockets"]
-git-tree-sha1 = "4e3136db3735924f96632a5b40a5979f1f53fa07"
-uuid = "da04e1cc-30fd-572f-bb4f-1f8673147195"
-version = "0.20.19"
-
-    [deps.MPI.extensions]
-    AMDGPUExt = "AMDGPU"
-    CUDAExt = "CUDA"
-
-    [deps.MPI.weakdeps]
-    AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e"
-    CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
-
 [[deps.MPICH_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"]
-git-tree-sha1 = "d8a7bf80c88326ebc98b7d38437208c3a0f20725"
+git-tree-sha1 = "4099bb6809ac109bfc17d521dad33763bcf026b7"
 uuid = "7cb0a576-ebde-5e09-9194-50597f1243b4"
-version = "4.2.1+0"
+version = "4.2.1+1"
 
 [[deps.MPIPreferences]]
 deps = ["Libdl", "Preferences"]
-git-tree-sha1 = "8f6af051b9e8ec597fa09d8885ed79fd582f33c9"
+git-tree-sha1 = "c105fe467859e7f6e9a852cb15cb4301126fac07"
 uuid = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267"
-version = "0.1.10"
+version = "0.1.11"
 
 [[deps.MPItrampoline_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"]
-git-tree-sha1 = "3f884417b47a96d87e7c6219f8f7b30ce67f4f2c"
+git-tree-sha1 = "8c35d5420193841b2f367e658540e8d9e0601ed0"
 uuid = "f1f71cc9-e9ae-5b93-9b94-4fe0e1ad3748"
-version = "5.3.3+0"
+version = "5.4.0+0"
 
 [[deps.MacroTools]]
 deps = ["Markdown", "Random"]
@@ -1275,9 +1219,16 @@ uuid = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 version = "0.2.4"
 
 [[deps.MultiBroadcastFusion]]
-git-tree-sha1 = "cb2fb12bdd092eea393b4bd3252815cc8a395f7f"
+git-tree-sha1 = "19b2f184b5882538bee0d5355b152212110e0fd2"
 uuid = "c3c07f87-98de-43f2-a76f-835b330b2cbb"
-version = "0.2.0"
+version = "0.3.1"
+
+    [deps.MultiBroadcastFusion.extensions]
+    MultiBroadcastFusionCUDAExt = ["CUDA", "Adapt"]
+
+    [deps.MultiBroadcastFusion.weakdeps]
+    Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+    CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 
 [[deps.NLSolversBase]]
 deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"]
@@ -1394,9 +1345,9 @@ version = "1.3.0"
 
 [[deps.Pixman_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl"]
-git-tree-sha1 = "64779bc4c9784fee475689a1752ef4d5747c5e87"
+git-tree-sha1 = "35621f10a7531bc8fa58f74610b1bfb70a3cfc6b"
 uuid = "30392449-352a-5448-841d-b1acce4e97dc"
-version = "0.42.2+0"
+version = "0.43.4+0"
 
 [[deps.Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
@@ -1411,9 +1362,9 @@ version = "0.3.3"
 
 [[deps.PlotThemes]]
 deps = ["PlotUtils", "Statistics"]
-git-tree-sha1 = "1f03a2d339f42dca4a4da149c7e15e9b896ad899"
+git-tree-sha1 = "6e55c6841ce3411ccb3457ee52fc48cb698d6fb0"
 uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
-version = "3.1.0"
+version = "3.2.0"
 
 [[deps.PlotUtils]]
 deps = ["ColorSchemes", "Colors", "Dates", "PrecompileTools", "Printf", "Random", "Reexport", "Statistics"]
@@ -1443,9 +1394,9 @@ version = "1.40.4"
 
 [[deps.Polyester]]
 deps = ["ArrayInterface", "BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "ManualMemory", "PolyesterWeave", "Requires", "Static", "StaticArrayInterface", "StrideArraysCore", "ThreadingUtilities"]
-git-tree-sha1 = "2ba5f33cbb51a85ef58a850749492b08f9bf2193"
+git-tree-sha1 = "b3e2bae88cf07baf0a051fe09666b8ef97aefe93"
 uuid = "f517fe37-dbe3-4b94-8317-1923a5111588"
-version = "0.7.13"
+version = "0.7.14"
 
 [[deps.PolyesterWeave]]
 deps = ["BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "Static", "ThreadingUtilities"]
@@ -1453,12 +1404,6 @@ git-tree-sha1 = "240d7170f5ffdb285f9427b92333c3463bf65bf6"
 uuid = "1d0040c9-8b98-4ee7-8388-3f51789ca0ad"
 version = "0.2.1"
 
-[[deps.PooledArrays]]
-deps = ["DataAPI", "Future"]
-git-tree-sha1 = "36d8b4b899628fb92c2749eb488d884a926614d3"
-uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
-version = "1.4.3"
-
 [[deps.PreallocationTools]]
 deps = ["Adapt", "ArrayInterface", "ForwardDiff"]
 git-tree-sha1 = "a660e9daab5db07adf3dedfe09b435cc530d855e"
@@ -1485,9 +1430,9 @@ version = "1.4.3"
 
 [[deps.PrettyTables]]
 deps = ["Crayons", "LaTeXStrings", "Markdown", "PrecompileTools", "Printf", "Reexport", "StringManipulation", "Tables"]
-git-tree-sha1 = "88b895d13d53b5577fd53379d913b9ab9ac82660"
+git-tree-sha1 = "66b20dd35966a748321d3b2537c4584cf40387c7"
 uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
-version = "2.3.1"
+version = "2.3.2"
 
 [[deps.Printf]]
 deps = ["Unicode"]
@@ -1499,6 +1444,11 @@ git-tree-sha1 = "b437cdb0385ed38312d91d9c00c20f3798b30256"
 uuid = "49802e3a-d2f1-5c88-81d8-b72133a6f568"
 version = "1.5.1"
 
+[[deps.PtrArrays]]
+git-tree-sha1 = "f011fbb92c4d401059b2212c05c0601b70f8b759"
+uuid = "43287f4e-b6f4-7ad1-bb20-aadabca52c3d"
+version = "1.2.0"
+
 [[deps.Qt6Base_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Vulkan_Loader_jll", "Xorg_libSM_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_cursor_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "libinput_jll", "xkbcommon_jll"]
 git-tree-sha1 = "37b7bb7aabf9a085e0044307e1717436117f2b3b"
@@ -1519,18 +1469,6 @@ uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 deps = ["SHA"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
-[[deps.Random123]]
-deps = ["Random", "RandomNumbers"]
-git-tree-sha1 = "4743b43e5a9c4a2ede372de7061eed81795b12e7"
-uuid = "74087812-796a-5b5d-8853-05524746bad3"
-version = "1.7.0"
-
-[[deps.RandomNumbers]]
-deps = ["Random", "Requires"]
-git-tree-sha1 = "043da614cc7e95c703498a491e2c21f58a2b8111"
-uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143"
-version = "1.5.3"
-
 [[deps.RecipesBase]]
 deps = ["PrecompileTools"]
 git-tree-sha1 = "5c3d09cc4f31f5fc6af001c250bf1278733100ff"
@@ -1545,9 +1483,9 @@ version = "0.6.12"
 
 [[deps.RecursiveArrayTools]]
 deps = ["Adapt", "ArrayInterface", "DocStringExtensions", "GPUArraysCore", "IteratorInterfaceExtensions", "LinearAlgebra", "RecipesBase", "SparseArrays", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables"]
-git-tree-sha1 = "d8f131090f2e44b145084928856a561c83f43b27"
+git-tree-sha1 = "23a2f876de8c1f1fda6345b289cd7a3641fd0340"
 uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
-version = "3.13.0"
+version = "3.21.0"
 
     [deps.RecursiveArrayTools.extensions]
     RecursiveArrayToolsFastBroadcastExt = "FastBroadcast"
@@ -1597,10 +1535,10 @@ uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
 version = "0.7.1"
 
 [[deps.Rmath_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "6ed52fdd3382cf21947b15e8870ac0ddbff736da"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "d483cd324ce5cf5d61b77930f0bbd6cb61927d21"
 uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
-version = "0.4.0+0"
+version = "0.4.2+0"
 
 [[deps.RootSolvers]]
 deps = ["ForwardDiff"]
@@ -1625,9 +1563,9 @@ version = "0.1.0"
 
 [[deps.SciMLBase]]
 deps = ["ADTypes", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "SciMLStructures", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables"]
-git-tree-sha1 = "beb1f94b08c4976ed1db0ca01b9e6bac89706faf"
+git-tree-sha1 = "9f59654e2a85017ee27b0f59c7fac5a57aa10ced"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-version = "2.35.0"
+version = "2.39.0"
 
     [deps.SciMLBase.extensions]
     SciMLBaseChainRulesCoreExt = "ChainRulesCore"
@@ -1655,9 +1593,9 @@ uuid = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 version = "0.3.8"
 
 [[deps.SciMLStructures]]
-git-tree-sha1 = "5833c10ce83d690c124beedfe5f621b50b02ba4d"
+git-tree-sha1 = "d778a74df2f64059c38453b34abad1953b2b8722"
 uuid = "53ae85a6-f571-4167-b2af-e1d143709226"
-version = "1.1.0"
+version = "1.2.0"
 
 [[deps.Scratch]]
 deps = ["Dates"]
@@ -1665,12 +1603,6 @@ git-tree-sha1 = "3bac05bc7e74a75fd9cba4295cde4045d9fe2386"
 uuid = "6c6a2e73-6563-6170-7368-637461726353"
 version = "1.2.1"
 
-[[deps.SentinelArrays]]
-deps = ["Dates", "Random"]
-git-tree-sha1 = "0e7508ff27ba32f26cd459474ca2ede1bc10991f"
-uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
-version = "1.4.1"
-
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
@@ -1707,9 +1639,9 @@ version = "1.10.0"
 
 [[deps.SpecialFunctions]]
 deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
-git-tree-sha1 = "e2cfc4012a19088254b3950b85c3c1d8882d864d"
+git-tree-sha1 = "2f5d4697f21388cbe1ff299430dd169ef97d7e14"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "2.3.1"
+version = "2.4.0"
 weakdeps = ["ChainRulesCore"]
 
     [deps.SpecialFunctions.extensions]
@@ -1737,9 +1669,9 @@ version = "1.5.0"
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "bf074c045d3d5ffd956fa0a461da38a44685d6b2"
+git-tree-sha1 = "9ae599cd7529cfce7fea36cf00a62cfc56f0f37c"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.3"
+version = "1.9.4"
 weakdeps = ["ChainRulesCore", "Statistics"]
 
     [deps.StaticArrays.extensions]
@@ -1813,10 +1745,10 @@ uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c"
 version = "7.2.1+1"
 
 [[deps.SymbolicIndexingInterface]]
-deps = ["Accessors", "ArrayInterface", "MacroTools", "RuntimeGeneratedFunctions", "StaticArraysCore"]
-git-tree-sha1 = "7a7be02e16d11c17e2407bab80c2dd1410f774cb"
+deps = ["Accessors", "ArrayInterface", "RuntimeGeneratedFunctions", "StaticArraysCore"]
+git-tree-sha1 = "a5f6f138b740c9d93d76f0feddd3092e6ef002b7"
 uuid = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
-version = "0.3.17"
+version = "0.3.22"
 
 [[deps.TOML]]
 deps = ["Dates"]
@@ -1870,14 +1802,14 @@ version = "0.5.2"
 
 [[deps.TimerOutputs]]
 deps = ["ExprTools", "Printf"]
-git-tree-sha1 = "f548a9e9c490030e545f72074a41edfd0e5bcdd7"
+git-tree-sha1 = "5a13ae8a41237cff5ecf34f73eb1b8f42fff6531"
 uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
-version = "0.5.23"
+version = "0.5.24"
 
 [[deps.TranscodingStreams]]
-git-tree-sha1 = "71509f04d045ec714c4748c785a59045c3736349"
+git-tree-sha1 = "a947ea21087caba0a798c5e494d0bb78e3a1a3a0"
 uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
-version = "0.10.7"
+version = "0.10.9"
 weakdeps = ["Random", "Test"]
 
     [deps.TranscodingStreams.extensions]
@@ -1924,9 +1856,9 @@ version = "0.4.1"
 
 [[deps.Unitful]]
 deps = ["Dates", "LinearAlgebra", "Random"]
-git-tree-sha1 = "3c793be6df9dd77a0cf49d80984ef9ff996948fa"
+git-tree-sha1 = "dd260903fdabea27d9b6021689b3cd5401a57748"
 uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
-version = "1.19.0"
+version = "1.20.0"
 weakdeps = ["ConstructionBase", "InverseFunctions"]
 
     [deps.Unitful.extensions]
@@ -1952,9 +1884,9 @@ version = "0.2.1"
 
 [[deps.UnsafeAtomicsLLVM]]
 deps = ["LLVM", "UnsafeAtomics"]
-git-tree-sha1 = "323e3d0acf5e78a56dfae7bd8928c989b4f3083e"
+git-tree-sha1 = "d9f5962fecd5ccece07db1ff006fb0b5271bdfdd"
 uuid = "d80eeb9a-aca5-4d75-85e5-170c8b632249"
-version = "0.1.3"
+version = "0.1.4"
 
 [[deps.Unzip]]
 git-tree-sha1 = "ca0969166a028236229f63514992fc073799bb78"
@@ -1981,9 +1913,9 @@ version = "1.31.0+0"
 
 [[deps.XML2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"]
-git-tree-sha1 = "532e22cf7be8462035d092ff21fada7527e2c488"
+git-tree-sha1 = "52ff2af32e591541550bd753c0da8b9bc92bb9d9"
 uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
-version = "2.12.6+0"
+version = "2.12.7+0"
 
 [[deps.XSLT_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
@@ -1998,16 +1930,16 @@ uuid = "ffd25f8a-64ca-5728-b0f7-c24cf3aae800"
 version = "5.4.6+0"
 
 [[deps.Xorg_libICE_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "e5becd4411063bdcac16be8b66fc2f9f6f1e8fe5"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "326b4fea307b0b39892b3e85fa451692eda8d46c"
 uuid = "f67eecfb-183a-506d-b269-f58e52b52d7c"
-version = "1.0.10+1"
+version = "1.1.1+0"
 
 [[deps.Xorg_libSM_jll]]
-deps = ["Libdl", "Pkg", "Xorg_libICE_jll"]
-git-tree-sha1 = "4a9d9e4c180e1e8119b5ffc224a7b59d3a7f7e18"
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libICE_jll"]
+git-tree-sha1 = "3796722887072218eabafb494a13c963209754ce"
 uuid = "c834827a-8449-5923-a945-d239c165b7dd"
-version = "1.2.3+0"
+version = "1.2.4+0"
 
 [[deps.Xorg_libX11_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
@@ -2034,10 +1966,10 @@ uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
 version = "1.1.4+0"
 
 [[deps.Xorg_libXext_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
-git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3"
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll"]
+git-tree-sha1 = "d2d1a5c49fae4ba39983f63de6afcbea47194e85"
 uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
-version = "1.3.4+4"
+version = "1.3.6+0"
 
 [[deps.Xorg_libXfixes_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
@@ -2064,10 +1996,10 @@ uuid = "ec84b674-ba8e-5d96-8ba1-2a689ba10484"
 version = "1.5.2+4"
 
 [[deps.Xorg_libXrender_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
-git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96"
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll"]
+git-tree-sha1 = "47e45cd78224c53109495b3e324df0c37bb61fbe"
 uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
-version = "0.9.10+4"
+version = "0.9.11+0"
 
 [[deps.Xorg_libpthread_stubs_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -2143,9 +2075,9 @@ version = "1.5.0+0"
 
 [[deps.YAML]]
 deps = ["Base64", "Dates", "Printf", "StringEncodings"]
-git-tree-sha1 = "e6330e4b731a6af7959673621e91645eb1356884"
+git-tree-sha1 = "669a78c59a307fa3ebc0a0bffd7ae83bd2184361"
 uuid = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
-version = "0.4.9"
+version = "0.4.10"
 
 [[deps.Zlib_jll]]
 deps = ["Libdl"]
@@ -2183,10 +2115,10 @@ uuid = "477f73a3-ac25-53e9-8cc3-50b2fa2566f0"
 version = "1.1.2+0"
 
 [[deps.libaom_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "3a2ea60308f0996d26f1e5354e10c24e9ef905d4"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "1827acba325fdcdf1d2647fc8d5301dd9ba43a9d"
 uuid = "a4ae2306-e953-59d6-aa16-d00cac43593b"
-version = "3.4.0+0"
+version = "3.9.0+0"
 
 [[deps.libass_jll]]
 deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
diff --git a/docs/Project.toml b/docs/Project.toml
index c03d5b3e2..ead45860d 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -11,8 +11,10 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 NVTX = "5da4648a-3479-48b8-97b9-01cb529c0a1f"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
diff --git a/docs/make.jl b/docs/make.jl
index 259c8bbaf..1e16fdf72 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,5 +1,23 @@
 using Documenter, DocumenterCitations
 using ClimaTimeSteppers
+using Literate
+
+tutorial_basedir = "tutorials"
+tutorial_basedir_from_here = joinpath(@__DIR__, "src", tutorial_basedir)
+
+jl_files_in_basedir = filter(endswith(".jl"), readdir(tutorial_basedir_from_here))
+
+println("Building literate tutorials...")
+generated_tutorials = String[]
+for filename in jl_files_in_basedir
+    Literate.markdown(
+        joinpath(tutorial_basedir_from_here, filename),
+        tutorial_basedir_from_here;
+        execute = true,
+        flavor = Literate.CommonMarkFlavor(),
+    )
+    push!(generated_tutorials, joinpath(tutorial_basedir, replace(filename, ".jl" => ".md")))
+end
 
 # https://github.com/jheinen/GR.jl/issues/278#issuecomment-587090846
 ENV["GKSwstype"] = "nul"
@@ -12,12 +30,14 @@ pages = [
     "Algorithm Formulations" => [
         "ODE Solvers" => "algorithm_formulations/ode_solvers.md",
         "Newtons Method" => "algorithm_formulations/newtons_method.md",
+        "Rosenbrock Method" => "algorithm_formulations/rosenbrock.md",
         "Old LSRK Formulations" => "algorithm_formulations/lsrk.md",
         "Old MRRK Formulations" => "algorithm_formulations/mrrk.md",
     ],
     "Test problems" => [
         "test_problems/index.md",
     ],
+    "Tutorials" => generated_tutorials,
     "API docs" => [
         "ODE Solvers" => "api/ode_solvers.md",
         "Newtons Method" => "api/newtons_method.md",
diff --git a/docs/src/algorithm_formulations/rosenbrock.md b/docs/src/algorithm_formulations/rosenbrock.md
new file mode 100644
index 000000000..5b1b730dc
--- /dev/null
+++ b/docs/src/algorithm_formulations/rosenbrock.md
@@ -0,0 +1,122 @@
+# Rosenbrock methods
+
+In this page, we introduce Rosenbrock-type methods to solve ordinary
+differential equations. In doing do, we roughly follow Chapter IV.7 of "Solving
+Ordinary Differential Equations II" by Hairer and Wanner. (Beware, notation is
+not the same.). In this spirit, let us introduce the need for Rosenbrock-type
+methods in the same way as Hairer and Wanner, by quoting Rosenbrock himself:
+
+> When the functions are non-linear, implicit equations can in general be solved
+> only by iteration. This is a severe drawback, as it adds to the problem of
+> stability, that of convergence of the iterative process. An alternative, which
+> avoids this difficulty, is ...
+
+Rosenbrock method!
+
+Before reading this page, we recommend reading the page on ODE solvers first
+[TODO: Add link]
+
+## Introduction to the formalism
+
+Let us consider an ordinary differential equation of the form
+
+$$ \frac{d}{dt}u(t) = T(u(t), t)\,,$$
+
+where $u$ is the state, $T$ is the tendency,
+and $t$ the time. For the sake of simplicity, let us ignore the explicit time
+dependency in the tendency (we will get back to that at the end).
+
+The simplest way to introduce the Rosenbrock method is to start from a
+diagonally implicit Runge-Kutta scheme (see page on DIRK). In an implicit
+Runge-Kutta method with $s$ stages and tableau $a, b, c$, the updated value
+$u_1$ for a step of size $\Delta t$ is obtained starting from the known value $u_0$
+with
+
+$$ u_1 = u_0 + \sum_{i=1}^s b_i k_i\,, $$
+with
+$$ k_i = \Delta t  T ( u_0 + \sum_{j=1}^{i-1}\alpha_{ij}k_j + \alpha_{ii} k_i)\,. $$
+$\alpha_{ij}$, $b_i$ are carefully chosen coefficients.
+
+Rosenbrock's idea consists in linearizing $T$. This operation can be interpreted
+as taking one iteration for Netwon solver and yields
+
+$$ k_i = \Delta t T(g_i) + \Delta t T'(g_i) \alpha_{ii} k_i $$
+
+with $J(g_i)$ Jacobian of the tendency and with
+
+$$ g_i = u_0 + \sum_{j=1}^{i-1}\alpha_{ij} k_j  $$
+
+In Rosenbrock-type methods, the Jacobian $T'$ is replaced with linear combinations
+of the Jacobian evaluated at $u_0$, $J$:
+$$ T'(g_i) \alpha_{ii} k_i \approx J \sum_{j=1}^{i-1}\gamma_{ij} k_j + J \gamma_{ii} k_i\,,$$
+with $\gamma_{ij}$ other coefficients that can be chosen.
+
+Now, each stage consists of solving a system of linear equations in $k_i$ with
+matrix $I - \Delta t \gamma_{ii}$: $$ (I - J \Delta t \gamma_{ii}) k_i = \Delta
+t T(g_i) + J \sum_{j=1}^{i-1}\gamma_{ij} k_j$$ for each $i$ in $1, \dots, s$.
+Once $k_i$ are known, $u_1$ can is easily computed and the process repeated for
+the next step.
+
+In practice, there are computational advantages at implementing a slight
+variation of what presented here.
+
+Let us define $\tilde{k}_{i} = \sum_{j=1}^{i} \gamma_{ij} k_j$. If the matrix
+$\gamma_{ij}$ is invertible, we can move freely from $k_i$ to $\tilde{k}_i$. A
+convenient step is to also define $C$, as
+
+$$ C = diag(\gamma_{11}^{-1}, \gamma_{22}^{-1}, \dots, \gamma_{ss}^{-1}) - \Gamma^{-1} $$
+
+Which establishes that
+
+$$ k_i = \frac{\tilde{k}_i}{\gamma_{ii}} - \sum_{j=1}^{i -1} c_{ij} \tilde{k}_j$$
+Substituting this, we obtain the following equations
+
+$$  (J \Delta t \gamma_{ii} - 1) \tilde{k}_i = - \Delta
+t \gamma_{ii} T(g_i) - \gamma_{ii} J \sum_{j=1}^{i-1}c_{ij} \tilde{k}_j \,, $$
+
+$$ g_i =  u_0 + \sum_{j=1}^{i-1}a_{ij} \tilde{k}_j \,,$$
+
+$$ u_1 = u_0 + \sum_{j=1}^{s} m_j \tilde{k}_j\,,$$
+with
+$$ (a_{ij}) = (\alpha_{ij}) \Gamma^{-1}\,, $$ and $$ (m_i) = (b_i) \Gamma^{-1} $$
+
+Finally, small changes are required to add support for explicit time derivative:
+
+$$ (J \Delta t \gamma_{ii} - 1) \tilde{k}_i = - \Delta
+t \gamma_{ii} T( t_0 + \alpha_i \Delta t,  g_i) - \gamma_{ii} J \sum_{j=1}^{i-1}c_{ij} \tilde{k}_j - \Delta
+t \gamma_{ii} \gamma_i \Delta
+t \frac{\partial T}{\partial t}(t_0, u_0) $$
+
+where we defined
+
+$\alpha_{i} = \sum_{j=1}^{i-1}\alpha_{ij} $
+
+$ \gamma _{i} = \sum_{j=1}^{i}\gamma_{ij} $
+
+> :note: You may wonder, what is up with all these ugly minus signs? Why don't
+> we cancel them out? The reason for this is because we want to have the
+> prefactor $(J \Delta t \gamma_{ii} - 1)$. This is called `Wfact` in the
+> language of `DifferentialEquations.jl`. Most other algorithms specify this
+> quantity, so it is convenient to be consistent.
+
+## Implementation
+
+In `ClimaTimeSteppers.jl`, we implement the last equation presented in the
+previous section. Currently, the only Rosenbrock-type algorithm implemented is
+`SSPKnoth`. `SSPKnoth` is 3-stage, second-order accurate Rosenbrock with
+
+$$ \alpha = \begin{bmatrix}
+    0 & 0 & 0 \\
+    1 & 0 & 0 \\
+    \frac{1}{4} & \frac{1}{4} & 0 \\
+    \end{bmatrix} $$
+
+$$ \Gamma = \begin{bmatrix}
+    1 & 0 & 0 \\
+    1 & 1 & 0 \\
+    -\frac{3}{4} & \frac{3}{4} & 1 \\
+    \end{bmatrix} $$
+
+and $$ b = (\frac{1}{6}, \frac{1}{6}, \frac{2}{3}) $$.
+
+At each stage, the state is updated to $g_i$ and precomputed quantities are updated. Tendencies are computed and, unless the stage is the last, DSS is applied to the sum of all the tendencies. After the last stage, DSS is applied to the incremented state.
diff --git a/docs/src/api/ode_solvers.md b/docs/src/api/ode_solvers.md
index 44b72a2a2..c19659ff2 100644
--- a/docs/src/api/ode_solvers.md
+++ b/docs/src/api/ode_solvers.md
@@ -51,6 +51,7 @@ HOMMEM1
 ARK2GKC
 ARK437L2SA1
 ARK548L2SA2
+SSPKnoth
 ```
 
 ## Explicit Algorithm Names
diff --git a/docs/src/dev/report_gen.jl b/docs/src/dev/report_gen.jl
index 738f83fca..bfcef9be5 100644
--- a/docs/src/dev/report_gen.jl
+++ b/docs/src/dev/report_gen.jl
@@ -17,6 +17,10 @@ let # Convergence
     title = "All Algorithms"
     algorithm_names = map(T -> T(), all_subtypes(ClimaTimeSteppers.AbstractAlgorithmName))
     algorithm_names = filter(name -> !(name isa ARK437L2SA1 || name isa ARK548L2SA2), algorithm_names) # too high order
+
+    # NOTE: Some imperfections in the convergence order for SSPKnoth are to be
+    # expected because we are not using the exact Jacobian
+
     verify_convergence(title, algorithm_names, ark_analytic_nonlin_test_cts(Float64), 200)
     verify_convergence(title, algorithm_names, ark_analytic_sys_test_cts(Float64), 400)
     verify_convergence(title, algorithm_names, ark_analytic_test_cts(Float64), 40000; super_convergence = (ARS121(),))
@@ -29,6 +33,8 @@ let # Convergence
         num_steps_scaling_factor = 4,
         numerical_reference_algorithm_name = ARS343(),
     )
+    rosenbrock_schems = filter(name -> name isa ClimaTimeSteppers.RosenbrockAlgorithmName, algorithm_names)
+    verify_convergence(title, rosenbrock_schems, climacore_1Dheat_test_implicit_cts(Float64), 400)
     verify_convergence(
         title,
         algorithm_names,
diff --git a/docs/src/plotting_utils.jl b/docs/src/plotting_utils.jl
index 8b715d500..2119f72de 100644
--- a/docs/src/plotting_utils.jl
+++ b/docs/src/plotting_utils.jl
@@ -59,6 +59,8 @@ imex_convergence_orders(::ARK548L2SA2) = (5, 5, 5)
 imex_convergence_orders(::SSP22Heuns) = (2, 2, 2)
 imex_convergence_orders(::SSP33ShuOsher) = (3, 3, 3)
 imex_convergence_orders(::RK4) = (4, 4, 4)
+# SSPKnoth is not really an IMEX method
+imex_convergence_orders(::SSPKnoth) = (2, 2, 2)
 
 # Compute a confidence interval for the convergence order, returning the
 # estimated convergence order and its uncertainty.
@@ -113,6 +115,8 @@ function verify_convergence(
     default_dt = t_end / num_steps
 
     algorithm(algorithm_name::ClimaTimeSteppers.ERKAlgorithmName) = ExplicitAlgorithm(algorithm_name)
+    algorithm(algorithm_name::ClimaTimeSteppers.SSPKnoth) =
+        ClimaTimeSteppers.RosenbrockAlgorithm(ClimaTimeSteppers.tableau(ClimaTimeSteppers.SSPKnoth()))
     algorithm(algorithm_name::ClimaTimeSteppers.IMEXARKAlgorithmName) =
         IMEXAlgorithm(algorithm_name, NewtonsMethod(; max_iters = linear_implicit ? 1 : 2))
 
diff --git a/docs/src/tutorials/diffusion.jl b/docs/src/tutorials/diffusion.jl
new file mode 100644
index 000000000..2b202bc51
--- /dev/null
+++ b/docs/src/tutorials/diffusion.jl
@@ -0,0 +1,273 @@
+# # A full example on how to use ClimaTimesteppers with ClimaCore
+#
+# In this tutorial, we will solve a diffusion equation on a sphere using
+# ClimaCore, and a mixed implicit-explicit solver in ClimaTimesteppers.
+#
+# First, we will set up the ClimaCore Spaces. We will define the implicit and
+# explicit tendencies and Jacobian. Next, we will set up the ODE problem and
+# solve it.
+#
+# ## The setup
+#
+# In this example, we consider a 3D spherical shell. We set initial data up on
+# the lower face of this shell, and we will see it diffuse horizontally and
+# vertically. We will treat the horizontal diffusion explicitly and the
+# vertical one implicitly.
+#
+# Let us start by importing the required pieces
+import SciMLBase
+import LinearAlgebra
+import ClimaTimeSteppers
+import ClimaCore
+import Plots
+import ClimaCore.MatrixFields: @name, ⋅
+
+# We also define some units, mostly to make the code more explicit
+const meters = meter = 1.0
+const kilometers = kilometer = 1.0
+const seconds = second = 1.0
+
+# Now, we will set up the ClimaCore Spaces. There is some boilerplate involved.
+
+# Radius and height of the spherical shell:
+radius = 6000kilometers
+height = 1kilometers
+
+# Details of the computational grid:
+number_horizontal_elements = 10
+horizontal_polynomial_order = 3
+number_vertical_elements = 10
+
+# We prepare a face-centered vertical grid by first creating the Domain, then
+# the Mesh, and finally the Space
+vertdomain = ClimaCore.Domains.IntervalDomain(
+    ClimaCore.Geometry.ZPoint(0kilometers),
+    ClimaCore.Geometry.ZPoint(height);
+    boundary_names = (:bottom, :top),
+);
+vertmesh = ClimaCore.Meshes.IntervalMesh(vertdomain; nelems = number_vertical_elements);
+vertspace = ClimaCore.Spaces.FaceFiniteDifferenceSpace(vertmesh);
+
+# For the horizontal grid, we create a 2D spectral element grid on a
+# EquiangularCubedSphere mesh with Gauss-Lagrange-Lobatto quadrature points
+horzdomain = ClimaCore.Domains.SphereDomain(radius);
+horzmesh = ClimaCore.Meshes.EquiangularCubedSphere(horzdomain, number_horizontal_elements);
+horztopology = ClimaCore.Topologies.Topology2D(ClimaCore.ClimaComms.context(), horzmesh);
+horzquad = ClimaCore.Spaces.Quadratures.GLL{horizontal_polynomial_order + 1}();
+horzspace = ClimaCore.Spaces.SpectralElementSpace2D(horztopology, horzquad);
+# ClimaComms.context() specify where the simulation should be run (ie, on CPU/GPU).
+
+# Finally, we can combine `vertspace` and `horzspace` by extruding them (ie,
+# taking their "Cartesian product")
+space = ClimaCore.Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vertspace);
+
+# Now, we define the initial data. As initial data, we prepare Gaussian
+# perturbation defined only on the lower vertical face. This Gaussian
+# perturbation will diffuse horizontally and vertically
+σ = 15.0
+
+(; lat, long, z) = ClimaCore.Fields.coordinate_field(space);
+φ_gauss = @. exp(-(lat^2 + long^2) / σ^2) * (z < 0.005);
+
+# Note how we multiplied by `(z < 0.005)` to ensure that the Gaussian
+# perturbation is only on the lowest face.
+#
+# When working with `ClimaCore`, `ClimaTimesteppers` requires all the evolved
+# variables to be packed in `ClimaCore.FieldVector`s. In this case, we only
+# have one variable, which we will call `my_var`
+#
+# Let us pack this into a a `FieldVector`
+Y₀ = ClimaCore.Fields.FieldVector(; my_var = copy(φ_gauss));
+
+# We copy `φ_gauss` because it it will be modified by the integrator.
+
+# ## The equations and the integrator
+
+# The diffusion equation for $u$ is $\partial_t u = K \nabla^2 u$, with $K$
+# diffusivity. To solve this equation with `ClimaTimesteppers`, we have to
+# provide the tendencies (ie, the right-hand-sides). We want to solve the
+# diffusion equation horizontally and vertically. Horizontally, we will treat
+# the evolution explicitly (currently, there is no other option for our spectral
+# elements). On the other hand, we will solve the vertical diffusion implicitly.
+# So, this problem has two tendencies, the explicit and the implicit ones.
+#
+# So, let us start by defining the explicit tendency.
+
+# A tendency for `ClimaTimesteppers.jl` is a function that generally takes four
+# arguments `(∂ₜY, Y, p, t)`. `∂ₜY` is the right-hand side, which has to be
+# modified by the function, `Y` is the current state, `p` is the cache (an
+# optional collection of auxiliary variables) and `t` the time. Tendencies in
+# `ClimaTimesteppers.jl` are in-place. In our case, we want the divergence of
+# the gradient (the Laplacian) of `my_var`.
+
+diverg = ClimaCore.Operators.WeakDivergence();
+grad = ClimaCore.Operators.Gradient();
+
+K = 3.0
+
+# Notice how we picked `WeakDivergence` instead of normal divergence.
+# TODO: Explain why
+
+function T_exp!(∂ₜY, Y, _, _)
+    ∂ₜY.my_var .= K .* diverg.(grad.(Y.my_var))
+    return nothing
+end
+
+# We do not specify boundary conditions here because we are on a sphere, where
+# the only meaningful boundary conditions are periodic.
+
+# Next, we move to the implicit tendency. This is much more involved to set up.
+# First, we need again divergence and gradient. However, this time we want them
+# for the vertical direction
+diverg_vert = ClimaCore.Operators.DivergenceC2F(;
+    bottom = ClimaCore.Operators.SetDivergence(0.0),
+    top = ClimaCore.Operators.SetDivergence(0.0),
+);
+grad_vert = ClimaCore.Operators.GradientF2C();
+
+# We choose `DivergenceC2F` and `GradientF2C` to implement the Laplacian. `F2C`
+# and `C2F` stand for `face-to-center` and `center-to-face`. The choice is
+# motivated so that we have a chain of operations that brings us back to faces.
+# We are working with faces so that we can cleanly impose boundary conditions.
+# The boundary condition we are setting here is null divergence at the top and
+# bottom of the domain.
+
+function actual_T_imp!(∂ₜY, Y, _, _)
+    ∂ₜY.my_var .= K .* diverg_vert.(grad_vert.(Y.my_var))
+    return nothing
+end
+
+# For an implicit solver, we also need to provide the Jacobian of the implicit
+# tendency, or precisely, the `Wfact`. `Wfact` is `dt*γ*J - 1`, where `J` is the
+# Jacobian, `dt` the timestep, and γ a factor that enter the specific solver. In
+# some places, you might find `Wfact_t`, meaning, Wfact` transformed. The
+# relationship between `Wfact` and `Wfact_t` is `Wfact = - dt*γ Wfact_t`.
+
+# TODO: This next section needs to be written by someone who knows how things
+# work
+
+# Defining `Wfact` for ClimaCore simulations is a little involved. First, we
+# need to define a new type, which we call `_ImplicitJacobian`. This type
+# contains a representation of the Jacobian and a solver from ClimaCore to solve
+# the resulting system of linear equations.
+
+struct _ImplicitJacobian{M, S}
+    matrix::M
+    solver::S
+end
+
+# Alongside the struct, we need to provide a new method to `LinearAlgebra.ldiv!`
+# that lowers to the ClimaCore solver. This will almost always be like this:
+
+# Using `a` to solve `x = a \ b`
+function LinearAlgebra.ldiv!(x, a::_ImplicitJacobian, b)
+    ClimaCore.MatrixFields.field_matrix_solve!(a.solver, x, a.matrix, b)
+    return nothing
+end
+
+# It is also convenient to specify a new method for similar, which some solvers
+# use internally.
+Base.similar(implicitjac::_ImplicitJacobian) = implicitjac
+
+# Finally, we provide a constructor that takes the ClimaCore objects and returns
+# the _ImplicitJacobian
+function ImplicitJacobian(matrix, Y; alg = ClimaCore.MatrixFields.BlockDiagonalSolve())
+    return _ImplicitJacobian(matrix, ClimaCore.MatrixFields.FieldMatrixSolver(alg, matrix, Y))
+end
+
+# The Jacobian matrix has to be specified as a
+# `ClimaCore.MatrixFields.FieldMatrix`. This matrix takes pairs of variable
+# names `name1` and `name2` and returns the type for `∂name1\∂name2`. Given that
+# we are working with gradient/divergence operations, our operations are
+# tridiagonal, so have that
+jacobian_matrix = ClimaCore.MatrixFields.FieldMatrix(
+    (@name(my_var), @name(my_var)) => similar(φ_gauss, ClimaCore.MatrixFields.TridiagonalMatrixRow{Float64}),
+);
+
+# Similarly, we define `Wfact`, as `dtγ J - I`.
+# TODO: Why are things the way the are?
+div_matrix = ClimaCore.MatrixFields.operator_matrix(diverg_vert)
+grad_matrix = ClimaCore.MatrixFields.operator_matrix(grad_vert)
+function Wfact(W, Y, p, dtγ, t)
+    @. W.matrix[@name(my_var), @name(my_var)] = dtγ * div_matrix() ⋅ grad_matrix() - (LinearAlgebra.I,)
+    return nothing
+end
+
+# With all of this, we are ready to define the implicit tendency. Implicit
+# tendencies are `SciMLBase.ODEFunction`s and take in the actual tendency
+# (similar to `T_exp!`), the Jacobian and `Wfact`:
+T_imp! = SciMLBase.ODEFunction(actual_T_imp!; jac_prototype = ImplicitJacobian(jacobian_matrix, Y₀), Wfact = Wfact);
+
+# On this type of spaces, we need to apply DSS to ensure continuity
+function dss!(state, p, t)
+    ClimaCore.Spaces.weighted_dss!(state.my_var)
+end
+
+# Now, we have all the pieces to set up the integrator
+t0 = 0seconds
+t_end = 500seconds
+dt = 5seconds
+
+prob = SciMLBase.ODEProblem(ClimaTimeSteppers.ClimaODEFunction(; T_imp!, T_exp!, dss!), Y₀, (t0, t_end), nothing);
+# prob = SciMLBase.ODEProblem(ClimaTimeSteppers.ClimaODEFunction(; T_exp!, dss!), Y₀, (t0, t_end), nothing)
+
+# We use SSPKnoth for this example
+algo = ClimaTimeSteppers.RosenbrockAlgorithm(ClimaTimeSteppers.tableau(ClimaTimeSteppers.SSPKnoth()));
+
+# And here is the integrator, where we set `saveat = dt` to save a snapshot of
+# the solution at every timestep.
+integrator = SciMLBase.init(prob, algo; dt, saveat = dt);
+
+## Solution and visualization
+#
+# To visualize the solution, we use `ClimaCore.Remapper` and Plots
+
+# Let us a prepare a convenience function that remaps a ClimaCore.Fields onto a
+# 2D Cartesian grid at a target z.
+function remap(; target_z = 0.0, integrator = integrator)
+    longpts = range(-180.0, 180.0, 180)
+    latpts = range(-90.0, 90.0, 90)
+    hcoords = [ClimaCore.Geometry.LatLongPoint(lat, long) for long in longpts, lat in latpts]
+    zcoords = [ClimaCore.Geometry.ZPoint(target_z)]
+    field = integrator.u.my_var
+    space = axes(field)
+    remapper = ClimaCore.Remapping.Remapper(space, hcoords, zcoords)
+    return ClimaCore.Remapping.interpolate(remapper, field)[:, :, begin]
+end
+
+# First, let us have a look at the surface, where we have the initial Gaussian
+# perturbation
+Plots.heatmap(remap());
+Plots.savefig("diff-hm1.png")
+# ![](diff-hm1.png)
+
+# Now, let us double check that it is empty at higher elevation
+Plots.heatmap(remap(; target_z = 0.1kilometers));
+Plots.savefig("diff-hm2.png")
+# ![](diff-hm2.png)
+
+# The extrema for `my_var` at the beginning of the simulation are
+extrema(integrator.u)
+
+# Let us focus on the surface level
+extrema(ClimaCore.Fields.level(integrator.u.my_var, ClimaCore.Utilities.PlusHalf(0)))
+
+# And the first level
+extrema(ClimaCore.Fields.level(integrator.u.my_var, ClimaCore.Utilities.PlusHalf(1)))
+
+# Let us solve the equation
+SciMLBase.solve!(integrator);
+
+# Now, the extreme for `my_var` will have decreased, due to diffusion
+extrema(integrator.u)
+
+# Let us focus on the surface level
+extrema(ClimaCore.Fields.level(integrator.u.my_var, ClimaCore.Utilities.PlusHalf(0)))
+
+# And the first level
+extrema(ClimaCore.Fields.level(integrator.u.my_var, ClimaCore.Utilities.PlusHalf(1)))
+
+# And we will see some development on the layers that did not have data before
+Plots.heatmap(remap(; target_z = 0.1kilometers))
+Plots.savefig("diff-hm3.png")
+# ![](diff-hm3.png)
diff --git a/ext/ClimaTimeSteppersBenchmarkToolsExt.jl b/ext/ClimaTimeSteppersBenchmarkToolsExt.jl
index 62ed39202..d7dc95ba2 100644
--- a/ext/ClimaTimeSteppersBenchmarkToolsExt.jl
+++ b/ext/ClimaTimeSteppersBenchmarkToolsExt.jl
@@ -36,6 +36,18 @@ n_calls_per_step(::CTS.ARS343, max_newton_iters) = Dict(
     "post_implicit!" => 4,
     "step!" => 1,
 )
+function get_n_calls_per_step(i::CTS.RosenbrockAlgorithm)
+    return Dict(
+        "Wfact" => 3,
+        "ldiv!" => 3,
+        "T_imp!" => 3,
+        "T_exp_T_lim!" => 4,
+        "lim!" => 4,
+        "dss!" => 4,
+        "post_explicit!" => 3,
+        "step!" => 1,
+    )
+end
 
 function maybe_push!(trials₀, name, f!, args, kwargs, only)
     if isnothing(only) || name in only
diff --git a/ext/benchmark_utils.jl b/ext/benchmark_utils.jl
index 95c7df58f..227563672 100644
--- a/ext/benchmark_utils.jl
+++ b/ext/benchmark_utils.jl
@@ -83,7 +83,8 @@ function get_trial(f, args, name; device, with_cu_prof = :bprofile, trace = fals
     return trial
 end
 
-get_W(i::CTS.DistributedODEIntegrator) = i.cache.newtons_method_cache.j
+get_W(i::CTS.DistributedODEIntegrator) = hasproperty(i.cache, :W) ? i.cache.W : i.cache.newtons_method_cache.j
+get_W(i::CTS.RosenbrockAlgorithm) = i.cache.W
 get_W(i) = i.cache.W
 f_args(i, f::CTS.ForwardEulerODEFunction) = (copy(i.u), i.u, i.p, i.t, i.dt)
 f_args(i, f) = (similar(i.u), i.u, i.p, i.t)
diff --git a/src/solvers/rosenbrock.jl b/src/solvers/rosenbrock.jl
index ca871dbba..84ed542bd 100644
--- a/src/solvers/rosenbrock.jl
+++ b/src/solvers/rosenbrock.jl
@@ -1,93 +1,239 @@
 export SSPKnoth
+using StaticArrays
+import DiffEqBase
+import LinearAlgebra: ldiv!, diagm
+import LinearAlgebra
 
-abstract type RosenbrockAlgorithm <: DistributedODEAlgorithm end
+abstract type RosenbrockAlgorithmName <: AbstractAlgorithmName end
 
-struct RosenbrockTableau{N, RT, N²}
-    A::SMatrix{N, N, RT, N²}
-    C::SMatrix{N, N, RT, N²}
-    Γ::SMatrix{N, N, RT, N²}
-    m::SMatrix{N, 1, RT, N}
+"""
+    RosenbrockTableau{N, RT, N²}
+
+Contains everything that defines a Rosenbrock-type method.
+
+- N: number of stages,
+
+Refer to the documentation for the precise meaning of the symbols below.
+"""
+struct RosenbrockTableau{N}
+    """A = α Γ⁻¹"""
+    A::SMatrix{N, N}
+    """Tableau used for the time-dependent part"""
+    α::SMatrix{N, N}
+    """Stepping matrix. C = 1/diag(Γ) - Γ⁻¹"""
+    C::SMatrix{N, N}
+    """Substage contribution matrix"""
+    Γ::SMatrix{N, N}
+    """m = b Γ⁻¹, used to compute the increments k"""
+    m::SMatrix{N, 1}
 end
 
-struct RosenbrockCache{Nstages, RT, N², A}
-    tableau::RosenbrockTableau{Nstages, RT, N²}
+"""
+    RosenbrockAlgorithm(tableau)
+
+Constructs a Rosenbrock algorithm for solving ODEs.
+"""
+struct RosenbrockAlgorithm{T <: RosenbrockTableau} <: ClimaTimeSteppers.DistributedODEAlgorithm
+    tableau::T
+end
+
+"""
+    RosenbrockCache{N, A, WT}
+
+Contains everything that is needed to run a Rosenbrock-type method.
+
+- Nstages: number of stages,
+- A: type of the evolved state (e.g., a ClimaCore.FieldVector),
+- WT: type of the Jacobian (Wfact)
+"""
+struct RosenbrockCache{Nstages, A, WT}
+    """Preallocated space for the state"""
     U::A
+
+    """Preallocated space for the tendency"""
     fU::A
+
+    """Preallocated space for the implicit contribution to the tendency"""
+    fU_imp::A
+
+    """Preallocated space for the explicit contribution to the tendency"""
+    fU_exp::A
+
+    """Preallocated space for the limited contribution to the tendency"""
+    fU_lim::A
+
+    """Contributions to the state for each stage"""
     k::NTuple{Nstages, A}
-    W::Any
-    linsolve!::Any
+
+    """Preallocated space for the Wfact, dtγJ - 𝕀, or Wfact_t, 𝕀/dtγ - J, with J the Jacobian of the implicit tendency"""
+    W::WT
+
+    """Preallocated space for the explicit time derivative of the tendency"""
+    ∂Y∂t::A
 end
 
 function init_cache(prob::DiffEqBase.AbstractODEProblem, alg::RosenbrockAlgorithm; kwargs...)
-
-    tab = tableau(alg, eltype(prob.u0))
-    Nstages = length(tab.m)
+    Nstages = length(alg.tableau.m)
     U = zero(prob.u0)
     fU = zero(prob.u0)
+    fU_imp = zero(prob.u0)
+    fU_exp = zero(prob.u0)
+    fU_lim = zero(prob.u0)
+    ∂Y∂t = zero(prob.u0)
     k = ntuple(n -> similar(prob.u0), Nstages)
-    W = prob.f.jac_prototype
-    linsolve! = alg.linsolve(Val{:init}, W, prob.u0; kwargs...)
-
-    return RosenbrockCache(tab, U, fU, k, W, linsolve!)
+    if !isnothing(prob.f.T_imp!)
+        W = prob.f.T_imp!.jac_prototype
+    else
+        W = nothing
+    end
+    return RosenbrockCache{Nstages, typeof(U), typeof(W)}(U, fU, fU_imp, fU_exp, fU_lim, k, W, ∂Y∂t)
 end
 
+"""
+    step_u!(int, cache::RosenbrockCache{Nstages})
+
+Take one step with the Rosenbrock-method with the given `cache`.
 
-function step_u!(int, cache::RosenbrockCache{Nstages, RT}) where {Nstages, RT}
-    (; m, Γ, A, C) = cache.tableau
+Some choices are being made here. Most of these are empirically motivated and should be
+revisited on different problems.
+- We do not update dtγ across stages
+- We apply DSS to the sum of the explicit and implicit tendency at all the stages but the last
+- We apply DSS to incremented state (ie, after the final stage is applied)
+"""
+function step_u!(int, cache::RosenbrockCache{Nstages}) where {Nstages}
+    (; m, Γ, A, α, C) = int.alg.tableau
     (; u, p, t, dt) = int
-    (; W, U, fU, k, linsolve!) = cache
-    f! = int.sol.prob.f
-    Wfact_t! = int.sol.prob.f.Wfact_t
+    (; W, U, fU, fU_imp, fU_exp, fU_lim, k, ∂Y∂t) = cache
+    T_imp! = int.sol.prob.f.T_imp!
+    T_exp! = int.sol.prob.f.T_exp!
+    T_exp_lim! = int.sol.prob.f.T_exp_T_lim!
+    tgrad! = isnothing(T_imp!) ? nothing : T_imp!.tgrad
+
+    (; post_explicit!, post_implicit!, dss!) = int.sol.prob.f
+
+    has_Wfact_t = false
+
+    # TODO: This is only valid when Γ[i, i] is constant, otherwise we have to
+    # move this in the for loop
+    dtγ = dt * Γ[1, 1]
+
+    if !isnothing(T_imp!)
+        Wfact! = int.sol.prob.f.T_imp!.Wfact
+        Wfact_t! = int.sol.prob.f.T_imp!.Wfact_t
+        if !isnothing(Wfact!) && !isnothing(Wfact_t!)
+            error("Only one between Wfact and Wfact_t can be non-nothing")
+        end
+        if !isnothing(Wfact!)
+            Wfact!(W, u, p, dtγ, t)
+        elseif !isnothing(Wfact_t!)
+            has_Wfact_t = true
+            Wfact_t!(W, u, p, dtγ, t)
+        end
+    end
+
+    if !isnothing(tgrad!)
+        tgrad!(∂Y∂t, u, p, t)
+    end
 
-    # 1) compute jacobian factorization
-    γ = dt * Γ[1, 1]
-    Wfact_t!(W, u, p, γ, t)
     for i in 1:Nstages
+        # Reset tendency
+        fill!(fU, 0)
+
+        αi = sum(α[i, 1:(i - 1)])
+        γi = sum(Γ[i, 1:i])
+
         U .= u
         for j in 1:(i - 1)
             U .+= A[i, j] .* k[j]
         end
-        # TODO: there should be a time modification here (t + c * dt)
-        # if f does depend on time, would need to add tgrad term as well
-        f!(fU, U, p, t)
+
+        if !isnothing(post_implicit!)
+            # NOTE: post_implicit! is a misnomer
+            post_implicit!(U, p, t + αi * dt)
+        end
+
+        if !isnothing(T_imp!)
+            T_imp!(fU_imp, U, p, t + αi * dt)
+            fU .+= fU_imp
+        end
+
+        if !isnothing(T_exp!)
+            T_exp!(fU_exp, U, p, t + αi * dt)
+            fU .+= fU_exp
+        end
+
+        if !isnothing(T_exp_lim!)
+            T_exp_lim!(fU_exp, fU_lim, U, p, t + αi * dt)
+            fU .+= fU_exp
+            fU .+= fU_lim
+        end
+
+        if !isnothing(tgrad!)
+            fU .+= γi .* dt .* ∂Y∂t
+        end
+
+        # We dss the tendency at every stage but the last. At the last stage, we
+        # dss the incremented state
+        (i != Nstages) && dss!(fU, p, t + αi * dt)
+
         for j in 1:(i - 1)
             fU .+= (C[i, j] / dt) .* k[j]
         end
-        linsolve!(k[i], W, fU)
+
+        if !has_Wfact_t
+            fU .*= -dtγ
+        end
+
+        if !isnothing(T_imp!)
+            if W isa Matrix
+                ldiv!(k[i], lu(W), fU)
+            else
+                ldiv!(k[i], W, fU)
+            end
+        else
+            k[i] .= .-fU
+        end
     end
+
     for i in 1:Nstages
         u .+= m[i] .* k[i]
     end
-end
 
-struct SSPKnoth{L} <: RosenbrockAlgorithm
-    linsolve::L
+    dss!(u, p, t + dt)
+    return nothing
 end
-SSPKnoth(; linsolve) = SSPKnoth(linsolve)
 
+"""
+    SSPKnoth
 
-function tableau(::SSPKnoth, RT)
-    # ROS.transformed=true;
+`SSPKnoth` is a second-order Rosenbrock method developed by Osward Knoth.
+
+The coefficients are the same as in `CGDycore.jl`, except that for C we add the
+diagonal elements too. Note, however, that the elements on the diagonal of C do
+not really matter because C is only used in its lower triangular part. We add them
+mostly to match literature on the subject
+"""
+struct SSPKnoth <: RosenbrockAlgorithmName end
+
+function tableau(::SSPKnoth)
     N = 3
     N² = N * N
-    α = @SMatrix RT[
+    α = @SMatrix [
         0 0 0
         1 0 0
         1/4 1/4 0
     ]
-    # ROS.d=ROS.alpha*ones(ROS.nStage,1);
-    b = @SMatrix RT[1 / 6 1 / 6 2 / 3]
-    Γ = @SMatrix RT[
+    b = @SMatrix [1 / 6 1 / 6 2 / 3]
+    Γ = @SMatrix [
         1 0 0
         0 1 0
         -3/4 -3/4 1
     ]
     A = α / Γ
-    C = -inv(Γ)
+    invΓ = inv(Γ)
+    diag_invΓ = SMatrix{N, N}(diagm([invΓ[i, i] for i in 1:N]))
+    # C is diag(γ₁₁⁻¹, γ₂₂⁻¹, ...) - Γ⁻¹
+    C = diag_invΓ .- inv(Γ)
     m = b / Γ
-    return RosenbrockTableau{N, RT, N²}(A, C, Γ, m)
-    #   ROS.SSP.alpha=[1 0 0
-    #                  3/4 1/4 0
-    #                  1/3 0 2/3];
-
+    return RosenbrockTableau{N}(A, α, C, Γ, m)
 end
diff --git a/test/Project.toml b/test/Project.toml
index 4d132ceb0..c3911b92c 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -27,7 +27,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Aqua = "0.8"
 CUDA = "3, 4, 5"
 ClimaComms = "0.4, 0.5, 0.6"
-ClimaCore = "0.10, 0.11, 0.12"
+ClimaCore = "0.14"
 DataStructures = "0.18"
 DiffEqBase = "6"
 DiffEqCallbacks = "2"
diff --git a/test/callbacks.jl b/test/callbacks.jl
index 914c5a82f..6fb2b39f4 100644
--- a/test/callbacks.jl
+++ b/test/callbacks.jl
@@ -2,7 +2,7 @@ using Test
 using ClimaComms
 using ClimaTimeSteppers, DiffEqBase
 using ClimaTimeSteppers.Callbacks
-
+import MPI
 @static isdefined(ClimaComms, :device_type) && ClimaComms.@import_required_backends
 device = ClimaComms.device()
 comm_ctx = ClimaComms.MPICommsContext(device)
diff --git a/test/convergence.jl b/test/convergence.jl
index 217e10396..72d208154 100644
--- a/test/convergence.jl
+++ b/test/convergence.jl
@@ -115,3 +115,19 @@ function tabulate_convergence_orders_multirate()
 end
 
 tabulate_convergence_orders_multirate()
+
+function tabulate_convergence_orders_rosenbrock()
+    tabs = [SSPKnoth]
+    tabs = map(t -> t(), tabs)
+    test_cases = all_test_cases(Float64)
+    results = convergence_order_results(tabs, test_cases)
+
+    prob_names = map(t -> t.test_name, test_cases)
+    expected_orders = SciMLBase.alg_order.(tabs)
+    algs = algorithm.(tabs)
+
+    tabulate_convergence_orders(prob_names, algs, results, expected_orders; tabs)
+    return results
+end
+
+tabulate_convergence_orders_rosenbrock()
diff --git a/test/convergence_orders.jl b/test/convergence_orders.jl
index c5e8c3b3f..8248068dc 100644
--- a/test/convergence_orders.jl
+++ b/test/convergence_orders.jl
@@ -31,6 +31,7 @@ second_order_tableau() = [
     SSP222,
     SSP322,
     SSP332,
+    SSPKnoth,
 ]
 
 #####
diff --git a/test/problems.jl b/test/problems.jl
index ce15d5fe4..033435702 100644
--- a/test/problems.jl
+++ b/test/problems.jl
@@ -1,9 +1,12 @@
-using DiffEqBase, ClimaTimeSteppers, LinearAlgebra, StaticArrays
+using DiffEqBase, ClimaTimeSteppers, StaticArrays
+import LinearAlgebra
+import LinearAlgebra: norm, Diagonal, mul!
 using ClimaCore
 using ClimaComms
+import ClimaCore.MatrixFields: @name, ⋅
 import ClimaCore: Domains, Geometry, Meshes, Topologies, Spaces, Fields, Operators, Limiters
 
-@static isdefined(ClimaComms, :device_type) && ClimaComms.@import_required_backends
+@static pkgversion(ClimaComms) >= v"0.6" && ClimaComms.@import_required_backends
 
 import Krylov # Trigger ClimaCore/ext/KrylovExt
 
@@ -366,6 +369,20 @@ function ark_analytic_test_cts(::Type{FT}) where {FT}
     )
 end
 
+function analytic_linear_no_source_test_cts(::Type{FT}) where {FT}
+    λ = FT(-1)
+    ClimaIntegratorTestCase(;
+        test_name = "analytic_linear_no_source",
+        linear_implicit = true,
+        t_end = FT(10),
+        Y₀ = FT[1],
+        analytic_sol = (t) -> [exp(λ * t)],
+        tendency! = (Yₜ, Y, _, t) -> Yₜ .= λ .* Y,
+        implicit_tendency! = (Yₜ, Y, _, t) -> Yₜ .= λ .* Y,
+        Wfact! = (W, Y, _, dtγ, t) -> W .= dtγ * λ - 1,
+    )
+end
+
 # From Section 1.2 of "Example Programs for ARKode v4.4.0" by D. R. Reynolds
 function ark_analytic_nonlin_test_cts(::Type{FT}) where {FT}
     ClimaIntegratorTestCase(;
@@ -422,19 +439,20 @@ function onewaycouple_mri_test_cts(::Type{FT}) where {FT}
     )
 end
 
-# Used only for inference
-struct DummySchurComplementW{T}
-    temp1::T
-    temp2::T
+
+struct ImplicitJacobianSolver{M, S}
+    matrix::M
+    solver::S
+end
+Base.similar(implicitjac::ImplicitJacobianSolver) = implicitjac
+function ImplicitJacobian(matrix, Y::ClimaCore.Fields.FieldVector; alg = ClimaCore.MatrixFields.BlockDiagonalSolve())
+    solver = ClimaCore.MatrixFields.FieldMatrixSolver(alg, matrix, Y)
+    return ImplicitJacobianSolver(matrix, solver)
 end
-DummySchurComplementW(Y::T) where {T} = DummySchurComplementW{T}(similar(Y), similar(Y))
-Base.similar(w::DummySchurComplementW) = w
-function LinearAlgebra.ldiv!(x, A::DummySchurComplementW, b)
-    A.temp1 .= b
-    LinearAlgebra.ldiv!(A.temp2, A, A.temp1)
-    x .= A.temp2
+function LinearAlgebra.ldiv!(x, a::ImplicitJacobianSolver, b)
+    ClimaCore.MatrixFields.field_matrix_solve!(a.solver, x, a.matrix, b)
+    return nothing
 end
-LinearAlgebra.ldiv!(x::Fields.FieldVector, A::DummySchurComplementW, b::Fields.FieldVector) = nothing
 
 Wfact!(W, Y, p, dtγ, t) = nothing
 
@@ -493,9 +511,11 @@ function climacore_2Dheat_test_cts(::Type{FT}) where {FT}
     post_implicit! = (u, _, t) -> nothing
     post_explicit! = (u, _, t) -> nothing
 
+    jacobian = ClimaCore.MatrixFields.FieldMatrix((@name(u), @name(u)) => FT(-1) * LinearAlgebra.I)
+
     T_imp! = SciMLBase.ODEFunction(
         (Yₜ, u, _, t) -> nothing;
-        jac_prototype = DummySchurComplementW(init_state),
+        jac_prototype = ImplicitJacobian(jacobian, init_state),
         Wfact = Wfact!,
         tgrad = (∂Y∂t, Y, p, t) -> (∂Y∂t .= 0),
     )
@@ -512,6 +532,7 @@ function climacore_2Dheat_test_cts(::Type{FT}) where {FT}
         make_prob(split_tendency_func),
     )
 end
+
 function climacore_1Dheat_test_cts(::Type{FT}) where {FT}
     n_elem_z = 10
     n_z = 1
@@ -554,7 +575,72 @@ function climacore_1Dheat_test_cts(::Type{FT}) where {FT}
     )
 end
 
-# "Dynamical Core Model Intercomparison Project (DCMIP) Test Case Document" by 
+function climacore_1Dheat_test_implicit_cts(::Type{FT}) where {FT}
+    n_elem_z = 10000
+    n_z = 1
+    f_0 = FT(0.0) # denoted by f̂₀ above
+    Δλ = FT(1) # denoted by Δλ̂ above
+    t_end = FT(0.1) # denoted by t̂ above
+
+    domain = Domains.IntervalDomain(Geometry.ZPoint(FT(0)), Geometry.ZPoint(FT(1)), boundary_names = (:bottom, :top))
+    mesh = Meshes.IntervalMesh(domain, nelems = n_elem_z)
+    space = Spaces.FaceFiniteDifferenceSpace(mesh)
+    (; z) = Fields.coordinate_field(space)
+
+    λ = (2 * FT(π) * n_z)^2
+    φ_sin = @. sin(2 * FT(π) * n_z * z)
+
+    init_state = Fields.FieldVector(; u = φ_sin)
+
+    diverg = Operators.DivergenceC2F(; bottom = Operators.SetDivergence(FT(0)), top = Operators.SetDivergence(FT(0)))
+    grad = Operators.GradientF2C()
+
+    diverg_matrix = ClimaCore.MatrixFields.operator_matrix(diverg)
+    grad_matrix = ClimaCore.MatrixFields.operator_matrix(grad)
+
+    function Wfact(W, Y, p, dtγ, t)
+        name = @name(u)
+        # NOTE: We need MatrixFields.⋅, not LinearAlgebra.⋅
+        @. W.matrix[name, name] = diverg_matrix() ⋅ grad_matrix() - (LinearAlgebra.I,)
+        return nothing
+    end
+
+    function T_imp_func!(tendency, state, _, t)
+        @. tendency.u = diverg.(grad.(state.u)) + f_0 * exp(-(λ + Δλ) * t) * φ_sin
+    end
+
+    function tgrad(∂Y∂t, state, _, t)
+        @. ∂Y∂t.u = -f_0 * (λ + Δλ) * exp(-(λ + Δλ) * t) * φ_sin
+    end
+
+    jacobian = ClimaCore.MatrixFields.FieldMatrix(
+        (@name(u), @name(u)) => similar(φ_sin, ClimaCore.MatrixFields.TridiagonalMatrixRow{FT}),
+    )
+
+    jac_prototype = ImplicitJacobian(jacobian, init_state)
+
+    T_imp! = SciMLBase.ODEFunction(T_imp_func!; jac_prototype = jac_prototype, Wfact = Wfact, tgrad = tgrad)
+
+    function analytic_sol(t)
+        state = similar(init_state)
+        @. state.u = (1 + f_0 / Δλ * (1 - exp(-Δλ * t))) * exp(-λ * t) * φ_sin
+        return state
+    end
+
+    tendency_func = ClimaODEFunction(; T_imp!)
+    split_tendency_func = tendency_func
+    make_prob(func) = ODEProblem(func, init_state, (FT(0), t_end), nothing)
+    IntegratorTestCase(
+        "1D Heat Equation Implicit",
+        false,
+        t_end,
+        analytic_sol,
+        make_prob(tendency_func),
+        make_prob(split_tendency_func),
+    )
+end
+
+# "Dynamical Core Model Intercomparison Project (DCMIP) Test Case Document" by
 # Ullrich et al., Section 1.1
 # (http://www-personal.umich.edu/~cjablono/DCMIP-2012_TestCaseDocument_v1.7.pdf)
 # Implemented in flux form, with an optional limiter and hyperdiffusion.
diff --git a/test/utils.jl b/test/utils.jl
index 600670bef..7f50be0e1 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -5,4 +5,5 @@ problem(test_case, tab::CTS.IMEXARKAlgorithmName) = test_case.split_prob
 problem(test_case, tab) = test_case.prob
 
 algorithm(tab::CTS.IMEXARKAlgorithmName) = CTS.IMEXAlgorithm(tab, NewtonsMethod(; max_iters = 2))
+algorithm(tab::CTS.RosenbrockAlgorithmName) = CTS.RosenbrockAlgorithm(ClimaTimeSteppers.tableau(tab))
 algorithm(tab) = tab