From 446f3cb1b80a77689d2e010297286f8c5b8553df Mon Sep 17 00:00:00 2001
From: Julia Sloan <jsloan@caltech.edu>
Date: Thu, 7 Mar 2024 15:51:02 -0800
Subject: [PATCH 1/2] fix Manifest

---
 experiments/AMIP/Manifest.toml | 52 ++++++++++++++--------------------
 1 file changed, 22 insertions(+), 30 deletions(-)

diff --git a/experiments/AMIP/Manifest.toml b/experiments/AMIP/Manifest.toml
index e32fc08dc1..337548e393 100644
--- a/experiments/AMIP/Manifest.toml
+++ b/experiments/AMIP/Manifest.toml
@@ -1,6 +1,6 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.0"
+julia_version = "1.10.1"
 manifest_format = "2.0"
 project_hash = "8d2ff3613412ce745a56f29971ef6211131f932f"
 
@@ -91,9 +91,9 @@ version = "7.7.1"
 
 [[deps.ArrayLayouts]]
 deps = ["FillArrays", "LinearAlgebra"]
-git-tree-sha1 = "64d582bcb9c93ac741234789eeb4f16812413efb"
+git-tree-sha1 = "e46675dbc095ddfdf2b5fba247d5a25f34e1f8a2"
 uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
-version = "1.6.0"
+version = "1.6.1"
 weakdeps = ["SparseArrays"]
 
     [deps.ArrayLayouts.extensions]
@@ -188,12 +188,6 @@ git-tree-sha1 = "f1dff6729bc61f4d49e140da1af55dcd1ac97b2f"
 uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 version = "1.5.0"
 
-[[deps.BinaryProvider]]
-deps = ["Libdl", "Logging", "SHA"]
-git-tree-sha1 = "ecdec412a9abc8db54c0efc5548c64dfce072058"
-uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
-version = "0.5.10"
-
 [[deps.BitFlags]]
 git-tree-sha1 = "2dc09997850d68179b69dafb58ae806167a32b1b"
 uuid = "d1d4a3ce-64b1-5f1a-9ba4-7e7e69966f35"
@@ -502,7 +496,7 @@ weakdeps = ["Dates", "LinearAlgebra"]
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.0.5+1"
+version = "1.1.0+0"
 
 [[deps.CompositionsBase]]
 git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad"
@@ -579,9 +573,9 @@ version = "1.6.1"
 
 [[deps.DataStructures]]
 deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
-git-tree-sha1 = "1fb174f0d48fe7d142e1109a10636bc1d14f5ac2"
+git-tree-sha1 = "0f4b5d62a88d8f59003e43c25a8a90de9eb76317"
 uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-version = "0.18.17"
+version = "0.18.18"
 
 [[deps.DataValueInterfaces]]
 git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
@@ -744,9 +738,9 @@ uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 version = "1.0.4"
 
 [[deps.EnzymeCore]]
-git-tree-sha1 = "59c44d8fbc651c0395d8a6eda64b05ce316f58b4"
+git-tree-sha1 = "496c5455d6a61c2a6f2233ce07c1fcdbe4995ab6"
 uuid = "f151be2c-9106-41f4-ab19-57ee4f262869"
-version = "0.6.5"
+version = "0.7.0"
 weakdeps = ["Adapt"]
 
     [deps.EnzymeCore.extensions]
@@ -792,16 +786,16 @@ uuid = "411431e0-e8b7-467b-b5e0-f676ba4f2910"
 version = "0.1.2"
 
 [[deps.FFMPEG]]
-deps = ["BinaryProvider", "Libdl"]
-git-tree-sha1 = "9143266ba77d3313a4cf61d8333a1970e8c5d8b6"
+deps = ["FFMPEG_jll"]
+git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8"
 uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
-version = "0.2.4"
+version = "0.4.1"
 
 [[deps.FFMPEG_jll]]
 deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
-git-tree-sha1 = "ab3f7e1819dba9434a3a5126510c8fda3a4e7000"
+git-tree-sha1 = "466d45dc38e15794ec7d5d63ec03d776a9aff36e"
 uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
-version = "6.1.1+0"
+version = "4.4.4+1"
 
 [[deps.FFTW]]
 deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
@@ -907,9 +901,9 @@ version = "0.8.4"
 
 [[deps.Flux]]
 deps = ["Adapt", "ChainRulesCore", "Compat", "Functors", "LinearAlgebra", "MLUtils", "MacroTools", "NNlib", "OneHotArrays", "Optimisers", "Preferences", "ProgressLogging", "Random", "Reexport", "SparseArrays", "SpecialFunctions", "Statistics", "Zygote"]
-git-tree-sha1 = "fd7b23aa8013a7528563d429f6eaf406f60364ed"
+git-tree-sha1 = "5a626d6ef24ae0a8590c22dc12096fb65eb66325"
 uuid = "587475ba-b771-5e3f-ad9e-33799f191a9c"
-version = "0.14.12"
+version = "0.14.13"
 
     [deps.Flux.extensions]
     FluxAMDGPUExt = "AMDGPU"
@@ -1386,9 +1380,9 @@ version = "0.2.4"
 
 [[deps.KernelAbstractions]]
 deps = ["Adapt", "Atomix", "InteractiveUtils", "LinearAlgebra", "MacroTools", "PrecompileTools", "Requires", "SparseArrays", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"]
-git-tree-sha1 = "c7753cc3febe006708ce6798482004241f7d890b"
+git-tree-sha1 = "ed7167240f40e62d97c1f5f7735dea6de3cc5c49"
 uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
-version = "0.9.17"
+version = "0.9.18"
 weakdeps = ["EnzymeCore"]
 
     [deps.KernelAbstractions.extensions]
@@ -1831,20 +1825,18 @@ version = "4.5.1"
 
 [[deps.NNlib]]
 deps = ["Adapt", "Atomix", "ChainRulesCore", "GPUArraysCore", "KernelAbstractions", "LinearAlgebra", "Pkg", "Random", "Requires", "Statistics"]
-git-tree-sha1 = "877f15c331337d54cf24c797d5bcb2e48ce21221"
+git-tree-sha1 = "6e4e90c2e2ef091ef50b91af65fa4bb09c3d0728"
 uuid = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
-version = "0.9.12"
+version = "0.9.6"
 
     [deps.NNlib.extensions]
     NNlibAMDGPUExt = "AMDGPU"
     NNlibCUDACUDNNExt = ["CUDA", "cuDNN"]
     NNlibCUDAExt = "CUDA"
-    NNlibEnzymeCoreExt = "EnzymeCore"
 
     [deps.NNlib.weakdeps]
     AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e"
     CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
-    EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869"
     cuDNN = "02a925ec-e4fe-4b08-9a7e-0d78e3d38ccd"
 
 [[deps.NVTX]]
@@ -1928,7 +1920,7 @@ version = "0.3.24+0"
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.23+2"
+version = "0.3.23+4"
 
 [[deps.OpenEXR]]
 deps = ["Colors", "FileIO", "OpenEXR_jll"]
@@ -2095,9 +2087,9 @@ version = "3.1.0"
 
 [[deps.PlotUtils]]
 deps = ["ColorSchemes", "Colors", "Dates", "PrecompileTools", "Printf", "Random", "Reexport", "Statistics"]
-git-tree-sha1 = "862942baf5663da528f66d24996eb6da85218e76"
+git-tree-sha1 = "7b1a9df27f072ac4c9c7cbe5efb198489258d1f5"
 uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
-version = "1.4.0"
+version = "1.4.1"
 
 [[deps.Plots]]
 deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "JLFzf", "JSON", "LaTeXStrings", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "Pkg", "PlotThemes", "PlotUtils", "PrecompileTools", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "RelocatableFolders", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs", "UnicodeFun", "UnitfulLatexify", "Unzip"]

From 1e66ec00e51e4ac181da545041e7631ef856e658 Mon Sep 17 00:00:00 2001
From: Julia Sloan <jsloan@caltech.edu>
Date: Tue, 27 Feb 2024 16:55:19 -0800
Subject: [PATCH 2/2] restructure component models

---
 docs/src/interfacer.md                        |   7 +-
 .../{climaatmos_init.jl => climaatmos.jl}     | 255 ++++++-------
 .../AMIP/components/land/bucket_init.jl       | 156 --------
 .../{bucket_utils.jl => climaland_bucket.jl}  | 323 +++++++++++++----
 ...nman_seaice_init.jl => eisenman_seaice.jl} | 341 +++++++++---------
 ...prescr_seaice_init.jl => prescr_seaice.jl} | 125 +++----
 .../{slab_ocean_init.jl => slab_ocean.jl}     |  95 ++---
 experiments/AMIP/coupler_driver.jl            |  11 +-
 test/component_model_tests/bucket_tests.jl    |   2 +-
 .../component_model_tests/climaatmos_tests.jl |   2 +-
 .../eisenman_seaice_tests.jl                  |   2 +-
 .../prescr_seaice_tests.jl                    |   2 +-
 .../component_model_tests/slab_ocean_tests.jl |   2 +-
 13 files changed, 681 insertions(+), 642 deletions(-)
 rename experiments/AMIP/components/atmosphere/{climaatmos_init.jl => climaatmos.jl} (97%)
 delete mode 100644 experiments/AMIP/components/land/bucket_init.jl
 rename experiments/AMIP/components/land/{bucket_utils.jl => climaland_bucket.jl} (52%)
 rename experiments/AMIP/components/ocean/{eisenman_seaice_init.jl => eisenman_seaice.jl} (98%)
 rename experiments/AMIP/components/ocean/{prescr_seaice_init.jl => prescr_seaice.jl} (98%)
 rename experiments/AMIP/components/ocean/{slab_ocean_init.jl => slab_ocean.jl} (97%)

diff --git a/docs/src/interfacer.md b/docs/src/interfacer.md
index f63b519db7..b9e76f11ee 100644
--- a/docs/src/interfacer.md
+++ b/docs/src/interfacer.md
@@ -29,7 +29,7 @@ information needed to run that simulation.
 
 Each `ComponentModelSimulation` must extend the following functions to be able
 to use our coupler. For some existing models, these are defined within
-ClimaCoupler.jl in that model’s `init.jl` file, but it is preferable
+ClimaCoupler.jl in that model’s file in `experiments/AMIP/components/`, but it is preferable
 for these to be defined in a model’s own repository. Note that the dispatch
 `::ComponentModelSimulation` in the function definitions given below should
 be replaced with the particular component model extending these functions.
@@ -72,7 +72,7 @@ for the following properties:
 
 | Coupler name      | Description | Units |
 |-------------------|-------------|-------|
-| `air_density`       | air density at the surface of the atmosphere | kg m^-3 |
+| `air_density`       | air density of the atmosphere | kg m^-3 |
 | `air_temperature`   | air temperature at the surface of the atmosphere | K |
 | `energy`            | globally integrated energy | J |
 | `height_int`        | height at the first internal model level | m |
@@ -105,6 +105,9 @@ following properties:
 | `surface_temperature` | temperature over the combined surface space | K |
 | `turbulent_fluxes` | turbulent fluxes (note: only required when using `PartitionedStateFluxes` option - see our `FluxCalculator` module docs for more information) | W m^-2 |
 
+- `calculate_surface_air_density(atmos_sim::Interfacer.AtmosModelSimulation, T_S::Fields.Field)`:
+A function to return the air density of the atmosphere simulation
+extrapolated to the surface, with units of [kg m^-3].
 
 ### SurfaceModelSimulation - required functions
 Analogously to the `AtmosModelSimulation`, a `SurfaceModelSimulation`
diff --git a/experiments/AMIP/components/atmosphere/climaatmos_init.jl b/experiments/AMIP/components/atmosphere/climaatmos.jl
similarity index 97%
rename from experiments/AMIP/components/atmosphere/climaatmos_init.jl
rename to experiments/AMIP/components/atmosphere/climaatmos.jl
index e07532852b..08ba2e717a 100644
--- a/experiments/AMIP/components/atmosphere/climaatmos_init.jl
+++ b/experiments/AMIP/components/atmosphere/climaatmos.jl
@@ -22,7 +22,9 @@ import ClimaCoupler.Utilities: swap_space!
 
 include("climaatmos_extra_diags.jl")
 
-# the clima atmos `integrator` is now defined
+###
+### Functions required by ClimaCoupler.jl for an AtmosModelSimulation
+###
 struct ClimaAtmosSimulation{P, Y, D, I} <: AtmosModelSimulation
     params::P
     Y_init::Y
@@ -31,52 +33,6 @@ struct ClimaAtmosSimulation{P, Y, D, I} <: AtmosModelSimulation
 end
 name(::ClimaAtmosSimulation) = "ClimaAtmosSimulation"
 
-"""
-    get_atmos_config(coupler_dict::Dict)
-
-Returns the specified atmospheric configuration (`atmos_config_dict`) overwitten by arguments
-in the coupler dictionary (`config_dict`).
-"""
-function get_atmos_config(coupler_dict)
-    atmos_config_file = coupler_dict["atmos_config_file"]
-    # override default or specified configs with coupler arguments, and set the correct atmos config_file
-    if isnothing(atmos_config_file)
-        @info "Using Atmos default configuration"
-        atmos_config = merge(CA.default_config_dict(), coupler_dict, Dict("config_file" => atmos_config_file))
-    else
-        @info "Using Atmos configuration from $atmos_config_file"
-        atmos_config = merge(
-            CA.override_default_config(joinpath(pkgdir(CA), atmos_config_file)),
-            coupler_dict,
-            Dict("config_file" => atmos_config_file),
-        )
-    end
-
-    # use coupler toml if atmos is not defined
-    atmos_toml_file = atmos_config["toml"]
-    coupler_toml_file = coupler_dict["coupler_toml_file"]
-    default_toml_file = "toml/default_coarse.toml"
-
-    toml_file = !isempty(atmos_toml_file) ? joinpath(pkgdir(CA), atmos_toml_file[1]) : nothing
-    toml_file = !isnothing(coupler_toml_file) ? joinpath(pkgdir(ClimaCoupler), coupler_toml_file) : toml_file
-    toml_file = isnothing(toml_file) ? joinpath(pkgdir(ClimaCoupler), default_toml_file) : toml_file
-
-    if !isnothing(toml_file)
-        @info "Overwriting Atmos parameters from $toml_file"
-        atmos_config = merge(atmos_config, Dict("toml" => [toml_file]))
-    end
-
-    # specify atmos output directory to be inside the coupler output directory
-    atmos_output_dir = joinpath(
-        coupler_dict["coupler_output_dir"],
-        joinpath(coupler_dict["mode_name"], coupler_dict["run_name"]),
-        "clima_atmos",
-    )
-    atmos_config = merge(atmos_config, Dict("output_dir" => atmos_output_dir))
-
-    return atmos_config
-end
-
 function atmos_init(::Type{FT}, atmos_config_dict::Dict) where {FT}
     # By passing `parsed_args` to `AtmosConfig`, `parsed_args` overwrites the default atmos config
     atmos_config = CA.AtmosConfig(atmos_config_dict)
@@ -119,17 +75,78 @@ function atmos_init(::Type{FT}, atmos_config_dict::Dict) where {FT}
     return sim
 end
 
+"""
+    get_model_prog_state(sim::ClimaAtmosSimulation)
+
+Extension of Checkpointer.get_model_prog_state to get the model state.
+"""
+function get_model_prog_state(sim::ClimaAtmosSimulation)
+    return sim.integrator.u
+end
+
+"""
+    get_field(atmos_sim::ClimaAtmosSimulation, ::Val{:radiative_energy_flux_toa})
+
+Extension of Interfacer.get_field to get the net TOA radiation, which is a sum of the
+upward and downward longwave and shortwave radiation.
+"""
+function get_field(atmos_sim::ClimaAtmosSimulation, ::Val{:radiative_energy_flux_toa})
+    FT = eltype(atmos_sim.integrator.u)
+
+    if atmos_sim.integrator.p.radiation.radiation_model != nothing
+        face_space = axes(atmos_sim.integrator.u.f)
+        nz_faces = length(ClimaCore.Spaces.vertical_topology(face_space).mesh.faces)
+
+        (; face_lw_flux_dn, face_lw_flux_up, face_sw_flux_dn, face_sw_flux_up) =
+            atmos_sim.integrator.p.radiation.radiation_model
+
+        LWd_TOA = ClimaCore.Fields.level(CA.RRTMGPI.array2field(FT.(face_lw_flux_dn), face_space), nz_faces - half)
+        LWu_TOA = ClimaCore.Fields.level(CA.RRTMGPI.array2field(FT.(face_lw_flux_up), face_space), nz_faces - half)
+        SWd_TOA = ClimaCore.Fields.level(CA.RRTMGPI.array2field(FT.(face_sw_flux_dn), face_space), nz_faces - half)
+        SWu_TOA = ClimaCore.Fields.level(CA.RRTMGPI.array2field(FT.(face_sw_flux_up), face_space), nz_faces - half)
+
+        return @. -(LWd_TOA + SWd_TOA - LWu_TOA - SWu_TOA)
+    else
+        return FT(0)
+    end
+end
+
+function get_field(atmos_sim::ClimaAtmosSimulation, ::Val{:energy})
+    thermo_params = get_thermo_params(atmos_sim)
+
+    ᶜS_ρq_tot = atmos_sim.integrator.p.precipitation.ᶜS_ρq_tot
+    ᶜts = atmos_sim.integrator.p.precomputed.ᶜts
+    ᶜΦ = atmos_sim.integrator.p.core.ᶜΦ
+
+    # return total energy and (if Microphysics0Moment) the energy lost due to precipitation removal
+    if atmos_sim.integrator.p.atmos.precip_model isa CA.Microphysics0Moment
+        ᶜS_ρq_tot = atmos_sim.integrator.p.precipitation.ᶜS_ρq_tot
+        ᶜts = atmos_sim.integrator.p.precomputed.ᶜts
+        ᶜΦ = atmos_sim.integrator.p.core.ᶜΦ
+        return atmos_sim.integrator.u.c.ρe_tot .-
+               ᶜS_ρq_tot .* CA.e_tot_0M_precipitation_sources_helper.(Ref(thermo_params), ᶜts, ᶜΦ) .*
+               atmos_sim.integrator.dt
+    else
+        return atmos_sim.integrator.u.c.ρe_tot
+    end
+end
+
 # extensions required by the Interfacer
+get_field(sim::ClimaAtmosSimulation, ::Val{:air_density}) =
+    TD.air_density.(thermo_params, sim.integrator.p.precomputed.ᶜts)
+get_field(sim::ClimaAtmosSimulation, ::Val{:air_temperature}) =
+    TD.air_temperature.(thermo_params, sim.integrator.p.precomputed.ᶜts)
+get_field(sim::ClimaAtmosSimulation, ::Val{:liquid_precipitation}) = sim.integrator.p.precipitation.col_integrated_rain
 get_field(sim::ClimaAtmosSimulation, ::Val{:radiative_energy_flux_sfc}) =
     ClimaCore.Fields.level(sim.integrator.p.radiation.ᶠradiation_flux, half)
-get_field(sim::ClimaAtmosSimulation, ::Val{:liquid_precipitation}) = sim.integrator.p.precipitation.col_integrated_rain # kg/m^2/s
-get_field(sim::ClimaAtmosSimulation, ::Val{:snow_precipitation}) = sim.integrator.p.precipitation.col_integrated_snow  # kg/m^2/s
+get_field(sim::ClimaAtmosSimulation, ::Val{:snow_precipitation}) = sim.integrator.p.precipitation.col_integrated_snow
 get_field(sim::ClimaAtmosSimulation, ::Val{:turbulent_energy_flux}) =
     ClimaCore.Geometry.WVector.(sim.integrator.p.precomputed.sfc_conditions.ρ_flux_h_tot)
 get_field(sim::ClimaAtmosSimulation, ::Val{:turbulent_moisture_flux}) =
     ClimaCore.Geometry.WVector.(sim.integrator.p.precomputed.sfc_conditions.ρ_flux_q_tot)
 get_field(sim::ClimaAtmosSimulation, ::Val{:thermo_state_int}) =
     ClimaCore.Spaces.level(sim.integrator.p.precomputed.ᶜts, 1)
+get_field(atmos_sim::ClimaAtmosSimulation, ::Val{:water}) = atmos_sim.integrator.u.c.ρq_tot
 
 # extensions required by FluxCalculator (partitioned fluxes)
 get_field(sim::ClimaAtmosSimulation, ::Val{:height_int}) =
@@ -140,12 +157,6 @@ function get_field(sim::ClimaAtmosSimulation, ::Val{:uv_int})
     uₕ_int = ClimaCore.Geometry.UVVector.(ClimaCore.Spaces.level(sim.integrator.u.c.uₕ, 1))
     return @. StaticArrays.SVector(uₕ_int.components.data.:1, uₕ_int.components.data.:2)
 end
-get_field(sim::ClimaAtmosSimulation, ::Val{:air_density}) =
-    TD.air_density.(thermo_params, sim.integrator.p.precomputed.ᶜts)
-get_field(sim::ClimaAtmosSimulation, ::Val{:air_temperature}) =
-    TD.air_temperature.(thermo_params, sim.integrator.p.precomputed.ᶜts)
-
-get_surface_params(sim::ClimaAtmosSimulation) = CAP.surface_fluxes_params(sim.integrator.p.params)
 
 function update_field!(atmos_sim::ClimaAtmosSimulation, ::Val{:co2}, field)
     if atmos_sim.integrator.p.atmos.radiation_mode isa CA.RRTMGPI.GrayRadiation
@@ -167,7 +178,6 @@ function update_field!(sim::ClimaAtmosSimulation, ::Val{:surface_albedo}, field)
         reshape(CA.RRTMGPI.field2array(field), 1, length(parent(field)))
 end
 
-# get_surface_params required by FluxCalculator (partitioned fluxes)
 function update_field!(sim::ClimaAtmosSimulation, ::Val{:turbulent_fluxes}, fields)
     (; F_turb_energy, F_turb_moisture, F_turb_ρτxz, F_turb_ρτyz) = fields
 
@@ -197,7 +207,6 @@ end
 
 # extensions required by FieldExchanger
 step!(sim::ClimaAtmosSimulation, t) = step!(sim.integrator, t - sim.integrator.t, true)
-
 reinit!(sim::ClimaAtmosSimulation) = reinit!(sim.integrator)
 
 function update_sim!(atmos_sim::ClimaAtmosSimulation, csf, turbulent_fluxes)
@@ -209,6 +218,69 @@ function update_sim!(atmos_sim::ClimaAtmosSimulation, csf, turbulent_fluxes)
     end
 end
 
+"""
+    calculate_surface_air_density(atmos_sim::ClimaAtmosSimulation, T_S::ClimaCore.Fields.Field)
+
+Extension for this function to calculate surface density.
+"""
+function calculate_surface_air_density(atmos_sim::ClimaAtmosSimulation, T_S::ClimaCore.Fields.Field)
+    thermo_params = get_thermo_params(atmos_sim)
+    ts_int = get_field(atmos_sim, Val(:thermo_state_int))
+    extrapolate_ρ_to_sfc.(Ref(thermo_params), ts_int, swap_space!(ones(axes(ts_int.ρ)), T_S))
+end
+
+# get_surface_params required by FluxCalculator (partitioned fluxes)
+get_surface_params(sim::ClimaAtmosSimulation) = CAP.surface_fluxes_params(sim.integrator.p.params)
+
+###
+### ClimaAtmos.jl model-specific functions (not explicitly required by ClimaCoupler.jl)
+###
+"""
+    get_atmos_config(coupler_dict::Dict)
+
+Returns the specified atmospheric configuration (`atmos_config_dict`) overwitten by arguments
+in the coupler dictionary (`config_dict`).
+"""
+function get_atmos_config(coupler_dict)
+    atmos_config_file = coupler_dict["atmos_config_file"]
+    # override default or specified configs with coupler arguments, and set the correct atmos config_file
+    if isnothing(atmos_config_file)
+        @info "Using Atmos default configuration"
+        atmos_config = merge(CA.default_config_dict(), coupler_dict, Dict("config_file" => atmos_config_file))
+    else
+        @info "Using Atmos configuration from $atmos_config_file"
+        atmos_config = merge(
+            CA.override_default_config(joinpath(pkgdir(CA), atmos_config_file)),
+            coupler_dict,
+            Dict("config_file" => atmos_config_file),
+        )
+    end
+
+    # use coupler toml if atmos is not defined
+    atmos_toml_file = atmos_config["toml"]
+    coupler_toml_file = coupler_dict["coupler_toml_file"]
+    default_toml_file = "toml/default_coarse.toml"
+
+    toml_file = !isempty(atmos_toml_file) ? joinpath(pkgdir(CA), atmos_toml_file[1]) : nothing
+    toml_file = !isnothing(coupler_toml_file) ? joinpath(pkgdir(ClimaCoupler), coupler_toml_file) : toml_file
+    toml_file = isnothing(toml_file) ? joinpath(pkgdir(ClimaCoupler), default_toml_file) : toml_file
+
+    if !isnothing(toml_file)
+        @info "Overwriting Atmos parameters from $toml_file"
+        atmos_config = merge(atmos_config, Dict("toml" => [toml_file]))
+    end
+
+    # specify atmos output directory to be inside the coupler output directory
+    atmos_output_dir = joinpath(
+        coupler_dict["coupler_output_dir"],
+        joinpath(coupler_dict["mode_name"], coupler_dict["run_name"]),
+        "clima_atmos",
+    )
+    atmos_config = merge(atmos_config, Dict("output_dir" => atmos_output_dir))
+
+    return atmos_config
+end
+
 
 # flux calculation borrowed from atmos
 """
@@ -296,75 +368,6 @@ Returns the thermodynamic parameters from the atmospheric model simulation objec
 """
 get_thermo_params(sim::ClimaAtmosSimulation) = CAP.thermodynamics_params(sim.integrator.p.params)
 
-"""
-    calculate_surface_air_density(atmos_sim::ClimaAtmosSimulation, T_S::ClimaCore.Fields.Field)
-
-Extension for this  to to calculate surface density.
-"""
-function calculate_surface_air_density(atmos_sim::ClimaAtmosSimulation, T_S::ClimaCore.Fields.Field)
-    thermo_params = get_thermo_params(atmos_sim)
-    ts_int = get_field(atmos_sim, Val(:thermo_state_int))
-    extrapolate_ρ_to_sfc.(Ref(thermo_params), ts_int, swap_space!(ones(axes(ts_int.ρ)), T_S))
-end
-
-"""
-    get_model_prog_state(sim::ClimaAtmosSimulation)
-
-Extension of Checkpointer.get_model_prog_state to get the model state.
-"""
-function get_model_prog_state(sim::ClimaAtmosSimulation)
-    return sim.integrator.u
-end
-
-"""
-    get_field(atmos_sim::ClimaAtmosSimulation, ::Val{:radiative_energy_flux_toa})
-
-Extension of Interfacer.get_field to get the net TOA radiation, which is a sum of the
-upward and downward longwave and shortwave radiation.
-"""
-function get_field(atmos_sim::ClimaAtmosSimulation, ::Val{:radiative_energy_flux_toa})
-    FT = eltype(atmos_sim.integrator.u)
-    # save radiation source
-    if atmos_sim.integrator.p.radiation.radiation_model != nothing
-        face_space = axes(atmos_sim.integrator.u.f)
-        nz_faces = length(ClimaCore.Spaces.vertical_topology(face_space).mesh.faces)
-
-        (; face_lw_flux_dn, face_lw_flux_up, face_sw_flux_dn, face_sw_flux_up) =
-            atmos_sim.integrator.p.radiation.radiation_model
-
-        LWd_TOA = ClimaCore.Fields.level(CA.RRTMGPI.array2field(FT.(face_lw_flux_dn), face_space), nz_faces - half)
-        LWu_TOA = ClimaCore.Fields.level(CA.RRTMGPI.array2field(FT.(face_lw_flux_up), face_space), nz_faces - half)
-        SWd_TOA = ClimaCore.Fields.level(CA.RRTMGPI.array2field(FT.(face_sw_flux_dn), face_space), nz_faces - half)
-        SWu_TOA = ClimaCore.Fields.level(CA.RRTMGPI.array2field(FT.(face_sw_flux_up), face_space), nz_faces - half)
-
-        return @. -(LWd_TOA + SWd_TOA - LWu_TOA - SWu_TOA) # [W/m^2]
-    else
-        return FT(0)
-    end
-end
-
-function get_field(atmos_sim::ClimaAtmosSimulation, ::Val{:energy})
-    thermo_params = get_thermo_params(atmos_sim)
-
-    ᶜS_ρq_tot = atmos_sim.integrator.p.precipitation.ᶜS_ρq_tot
-    ᶜts = atmos_sim.integrator.p.precomputed.ᶜts
-    ᶜΦ = atmos_sim.integrator.p.core.ᶜΦ
-
-    # return total energy and (if Microphysics0Moment) the energy lost due to precipitation removal
-    if atmos_sim.integrator.p.atmos.precip_model isa CA.Microphysics0Moment
-        ᶜS_ρq_tot = atmos_sim.integrator.p.precipitation.ᶜS_ρq_tot
-        ᶜts = atmos_sim.integrator.p.precomputed.ᶜts
-        ᶜΦ = atmos_sim.integrator.p.core.ᶜΦ
-        return atmos_sim.integrator.u.c.ρe_tot .-
-               ᶜS_ρq_tot .* CA.e_tot_0M_precipitation_sources_helper.(Ref(thermo_params), ᶜts, ᶜΦ) .*
-               atmos_sim.integrator.dt
-    else
-        return atmos_sim.integrator.u.c.ρe_tot
-    end
-end
-
-get_field(atmos_sim::ClimaAtmosSimulation, ::Val{:water}) = atmos_sim.integrator.u.c.ρq_tot
-
 """
     dss_state!(sim::ClimaAtmosSimulation)
 
diff --git a/experiments/AMIP/components/land/bucket_init.jl b/experiments/AMIP/components/land/bucket_init.jl
deleted file mode 100644
index a118d5cdf0..0000000000
--- a/experiments/AMIP/components/land/bucket_init.jl
+++ /dev/null
@@ -1,156 +0,0 @@
-# slab_rhs!
-using Dates: DateTime
-using Statistics: mean
-import SciMLBase: ODEProblem, init
-
-using ClimaCore
-import ClimaTimeSteppers as CTS
-import Thermodynamics as TD
-import CLIMAParameters
-import ClimaLand
-using ClimaLand.Bucket: BucketModel, BucketModelParameters, AbstractAtmosphericDrivers, AbstractRadiativeDrivers
-import ClimaLand.Bucket: BulkAlbedoTemporal, BulkAlbedoStatic, BulkAlbedoFunction
-using ClimaLand:
-    make_exp_tendency,
-    initialize,
-    make_set_initial_cache,
-    surface_evaporative_scaling,
-    CoupledRadiativeFluxes,
-    CoupledAtmosphere
-import ClimaLand.Parameters as LP
-
-import ClimaCoupler.Interfacer: LandModelSimulation, get_field, update_field!, name
-import ClimaCoupler.FieldExchanger: step!, reinit!
-import ClimaCoupler.FluxCalculator: update_turbulent_fluxes_point!, surface_thermo_state
-
-"""
-    BucketSimulation{M, Y, D, I}
-
-The bucket model simulation object.
-"""
-struct BucketSimulation{M, Y, D, I, A} <: LandModelSimulation
-    model::M
-    Y_init::Y
-    domain::D
-    integrator::I
-    area_fraction::A
-end
-name(::BucketSimulation) = "BucketSimulation"
-
-include("./bucket_utils.jl")
-
-"""
-    temp_anomaly_aquaplanet(coord)
-
-Introduce a temperature IC anomaly for the aquaplanet case.
-The values for this case follow the moist Held-Suarez setup of Thatcher &
-Jablonowski (2016, eq. 6), consistent with ClimaAtmos aquaplanet.
-"""
-temp_anomaly_aquaplanet(coord) = 29 * exp(-coord.lat^2 / (2 * 26^2))
-
-"""
-    temp_anomaly_amip(coord)
-
-Introduce a temperature IC anomaly for the AMIP case.
-The values used in this case have been tuned to align with observed temperature
-and result in stable simulations.
-"""
-temp_anomaly_amip(coord) = 40 * cosd(coord.lat)^4
-
-"""
-    bucket_init
-
-Initializes the bucket model variables.
-"""
-function bucket_init(
-    ::Type{FT},
-    tspan::Tuple{Float64, Float64},
-    config::String,
-    albedo_type::String,
-    land_temperature_anomaly::String,
-    regrid_dirpath::String;
-    space,
-    dt::Float64,
-    saveat::Float64,
-    area_fraction,
-    stepper = CTS.RK4(),
-    date_ref::DateTime,
-    t_start::Float64,
-) where {FT}
-    if config != "sphere"
-        println(
-            "Currently only spherical shell domains are supported; single column set-up will be addressed in future PR.",
-        )
-        @assert config == "sphere"
-    end
-
-    earth_param_set = LP.LandParameters(FT)
-
-    α_snow = FT(0.8) # snow albedo
-    if albedo_type == "map_static" # Read in albedo from static data file (default type)
-        # By default, this uses a file containing bareground albedo without a time component. Snow albedo is specified separately.
-        albedo = BulkAlbedoStatic{FT}(regrid_dirpath, space, α_snow = α_snow)
-    elseif albedo_type == "map_temporal" # Read in albedo from data file containing data over time
-        # By default, this uses a file containing linearly-interpolated monthly data of total albedo, generated by CESM2's land model (CLM).
-        albedo = BulkAlbedoTemporal{FT}(regrid_dirpath, date_ref, t_start, space)
-    elseif albedo_type == "function" # Use prescribed function of lat/lon for surface albedo
-        function α_bareground(coordinate_point)
-            (; lat, long) = coordinate_point
-            return typeof(lat)(0.38)
-        end
-        albedo = BulkAlbedoFunction{FT}(α_snow, α_bareground, space)
-    else
-        error("invalid albedo type $albedo_type")
-    end
-
-    σS_c = FT(0.2)
-    W_f = FT(10)
-    d_soil = FT(3.5) # soil depth
-    z_0m = FT(1e-3) # roughness length for momentum over smooth bare soil
-    z_0b = FT(1e-3) # roughness length for tracers over smooth bare soil
-    κ_soil = FT(0.7)
-    ρc_soil = FT(2e8)
-    t_crit = FT(dt) # This is the timescale on which snow exponentially damps to zero, in the case where all
-    # the snow would melt in time t_crit. It prevents us from having to specially time step in cases where
-    # all the snow melts in a single timestep.
-    params = BucketModelParameters(κ_soil, ρc_soil, albedo, σS_c, W_f, z_0m, z_0b, t_crit, earth_param_set)
-    n_vertical_elements = 7
-    # Note that this does not take into account topography of the surface, which is OK for this land model.
-    # But it must be taken into account when computing surface fluxes, for Δz.
-    domain = make_land_domain(space, (-d_soil, FT(0.0)), n_vertical_elements)
-    args = (params, CoupledAtmosphere{FT}(), CoupledRadiativeFluxes{FT}(), domain)
-    model = BucketModel{FT, typeof.(args)...}(args...)
-
-    # Initial conditions with no moisture
-    Y, p, coords = initialize(model)
-
-    # Get temperature anomaly function
-    T_functions = Dict("aquaplanet" => temp_anomaly_aquaplanet, "amip" => temp_anomaly_amip)
-    haskey(T_functions, land_temperature_anomaly) ||
-        error("land temp anomaly function $land_temperature_anomaly not supported")
-    temp_anomaly = T_functions[land_temperature_anomaly]
-
-    # Set temperature IC including anomaly, based on atmospheric setup
-    T_sfc_0 = FT(271.0)
-    @. Y.bucket.T = T_sfc_0 + temp_anomaly(coords.subsurface)
-
-    Y.bucket.W .= 6.5
-    Y.bucket.Ws .= 0.0
-    Y.bucket.σS .= 0.0
-
-    # Set initial aux variable values
-    set_initial_cache! = make_set_initial_cache(model)
-    set_initial_cache!(p, Y, tspan[1])
-
-    exp_tendency! = make_exp_tendency(model)
-    ode_algo = CTS.ExplicitAlgorithm(stepper)
-    bucket_ode_function = CTS.ClimaODEFunction(T_exp! = exp_tendency!, dss! = ClimaLand.dss!)
-    prob = ODEProblem(bucket_ode_function, Y, tspan, p)
-    integrator = init(prob, ode_algo; dt = dt, saveat = saveat, adaptive = false)
-
-    sim = BucketSimulation(model, Y, (; domain = domain, soil_depth = d_soil), integrator, area_fraction)
-
-    # DSS state to ensure we have continuous fields
-    dss_state!(sim)
-    return sim
-end
diff --git a/experiments/AMIP/components/land/bucket_utils.jl b/experiments/AMIP/components/land/climaland_bucket.jl
similarity index 52%
rename from experiments/AMIP/components/land/bucket_utils.jl
rename to experiments/AMIP/components/land/climaland_bucket.jl
index 77f9bed646..d9cc5ba672 100644
--- a/experiments/AMIP/components/land/bucket_utils.jl
+++ b/experiments/AMIP/components/land/climaland_bucket.jl
@@ -1,81 +1,213 @@
+# slab_rhs!
+using ClimaCore
+import ClimaTimeSteppers as CTS
+import Thermodynamics as TD
+using Dates: DateTime
+using ClimaComms: AbstractCommsContext
+import CLIMAParameters
+
+import ClimaLand
+using ClimaLand.Bucket: BucketModel, BucketModelParameters, AbstractAtmosphericDrivers, AbstractRadiativeDrivers
+import ClimaLand.Bucket: BulkAlbedoTemporal, BulkAlbedoStatic, BulkAlbedoFunction
+using ClimaLand:
+    make_exp_tendency,
+    initialize,
+    make_set_initial_cache,
+    surface_evaporative_scaling,
+    CoupledRadiativeFluxes,
+    CoupledAtmosphere
+import ClimaLand.Parameters as LP
+
+import ClimaCoupler.Interfacer: LandModelSimulation, get_field, update_field!, name
+import ClimaCoupler.FieldExchanger: step!, reinit!
+import ClimaCoupler.FluxCalculator: update_turbulent_fluxes_point!, surface_thermo_state
+
+###
+### Functions required by ClimaCoupler.jl for a SurfaceModelSimulation
+###
 """
-    make_land_domain(
-        atmos_boundary_space::ClimaCore.Spaces.SpectralElementSpace2D,
-        zlim::Tuple{FT, FT},
-        nelements_vert::Int,) where {FT}
+    BucketSimulation{M, Y, D, I}
 
-Creates the land model domain from the horizontal space of the atmosphere, and information
-about the number of elements and extent of the vertical domain.
+The bucket model simulation object.
 """
-function make_land_domain(
-    atmos_boundary_space::ClimaCore.Spaces.SpectralElementSpace2D,
-    zlim::Tuple{FT, FT},
-    nelements_vert::Int,
+struct BucketSimulation{M, Y, D, I, A} <: LandModelSimulation
+    model::M
+    Y_init::Y
+    domain::D
+    integrator::I
+    area_fraction::A
+end
+name(::BucketSimulation) = "BucketSimulation"
+
+
+"""
+    bucket_init
+
+Initializes the bucket model variables.
+"""
+function bucket_init(
+    ::Type{FT},
+    tspan::Tuple{Float64, Float64},
+    config::String,
+    albedo_type::String,
+    land_temperature_anomaly::String,
+    regrid_dirpath::String;
+    space,
+    dt::Float64,
+    saveat::Float64,
+    area_fraction,
+    stepper = CTS.RK4(),
+    date_ref::DateTime,
+    t_start::Float64,
 ) where {FT}
-    @assert zlim[1] < zlim[2]
-    depth = zlim[2] - zlim[1]
+    if config != "sphere"
+        println(
+            "Currently only spherical shell domains are supported; single column set-up will be addressed in future PR.",
+        )
+        @assert config == "sphere"
+    end
 
-    mesh = ClimaCore.Spaces.topology(atmos_boundary_space).mesh
+    earth_param_set = LP.LandParameters(FT)
 
-    radius = mesh.domain.radius
-    nelements_horz = mesh.ne
-    npolynomial =
-        ClimaCore.Spaces.Quadratures.polynomial_degree(ClimaCore.Spaces.quadrature_style(atmos_boundary_space))
-    nelements = (nelements_horz, nelements_vert)
-    vertdomain = ClimaCore.Domains.IntervalDomain(
-        ClimaCore.Geometry.ZPoint(FT(zlim[1])),
-        ClimaCore.Geometry.ZPoint(FT(zlim[2]));
-        boundary_tags = (:bottom, :top),
-    )
+    α_snow = FT(0.8) # snow albedo
+    if albedo_type == "map_static" # Read in albedo from static data file (default type)
+        # By default, this uses a file containing bareground albedo without a time component. Snow albedo is specified separately.
+        albedo = BulkAlbedoStatic{FT}(regrid_dirpath, space, α_snow = α_snow)
+    elseif albedo_type == "map_temporal" # Read in albedo from data file containing data over time
+        # By default, this uses a file containing linearly-interpolated monthly data of total albedo, generated by CESM2's land model (CLM).
+        albedo = BulkAlbedoTemporal{FT}(regrid_dirpath, date_ref, t_start, space)
+    elseif albedo_type == "function" # Use prescribed function of lat/lon for surface albedo
+        function α_bareground(coordinate_point)
+            (; lat, long) = coordinate_point
+            return typeof(lat)(0.38)
+        end
+        albedo = BulkAlbedoFunction{FT}(α_snow, α_bareground, space)
+    else
+        error("invalid albedo type $albedo_type")
+    end
 
-    vertmesh = ClimaCore.Meshes.IntervalMesh(vertdomain, ClimaCore.Meshes.Uniform(), nelems = nelements[2])
-    verttopology = ClimaCore.Topologies.IntervalTopology(vertmesh)
-    vert_center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(verttopology)
-    subsurface_space = ClimaCore.Spaces.ExtrudedFiniteDifferenceSpace(atmos_boundary_space, vert_center_space)
-    space = (; surface = atmos_boundary_space, subsurface = subsurface_space)
+    σS_c = FT(0.2)
+    W_f = FT(10)
+    d_soil = FT(3.5) # soil depth
+    z_0m = FT(1e-3) # roughness length for momentum over smooth bare soil
+    z_0b = FT(1e-3) # roughness length for tracers over smooth bare soil
+    κ_soil = FT(0.7)
+    ρc_soil = FT(2e8)
+    t_crit = FT(dt) # This is the timescale on which snow exponentially damps to zero, in the case where all
+    # the snow would melt in time t_crit. It prevents us from having to specially time step in cases where
+    # all the snow melts in a single timestep.
+    params = BucketModelParameters(κ_soil, ρc_soil, albedo, σS_c, W_f, z_0m, z_0b, t_crit, earth_param_set)
+    n_vertical_elements = 7
+    # Note that this does not take into account topography of the surface, which is OK for this land model.
+    # But it must be taken into account when computing surface fluxes, for Δz.
+    domain = make_land_domain(space, (-d_soil, FT(0.0)), n_vertical_elements)
+    args = (params, CoupledAtmosphere{FT}(), CoupledRadiativeFluxes{FT}(), domain)
+    model = BucketModel{FT, typeof.(args)...}(args...)
 
-    return ClimaLand.Domains.SphericalShell{FT}(radius, depth, nothing, nelements, npolynomial, space)
+    # Initial conditions with no moisture
+    Y, p, coords = initialize(model)
+
+    # Get temperature anomaly function
+    T_functions = Dict("aquaplanet" => temp_anomaly_aquaplanet, "amip" => temp_anomaly_amip)
+    haskey(T_functions, land_temperature_anomaly) ||
+        error("land temp anomaly function $land_temperature_anomaly not supported")
+    temp_anomaly = T_functions[land_temperature_anomaly]
+
+    # Set temperature IC including anomaly, based on atmospheric setup
+    T_sfc_0 = FT(271.0)
+    @. Y.bucket.T = T_sfc_0 + temp_anomaly(coords.subsurface)
+
+    Y.bucket.W .= 6.5
+    Y.bucket.Ws .= 0.0
+    Y.bucket.σS .= 0.0
+
+    # Set initial aux variable values
+    set_initial_cache! = make_set_initial_cache(model)
+    set_initial_cache!(p, Y, tspan[1])
+
+    exp_tendency! = make_exp_tendency(model)
+    ode_algo = CTS.ExplicitAlgorithm(stepper)
+    bucket_ode_function = CTS.ClimaODEFunction(T_exp! = exp_tendency!, dss! = ClimaLand.dss!)
+    prob = ODEProblem(bucket_ode_function, Y, tspan, p)
+    integrator = init(prob, ode_algo; dt = dt, saveat = saveat, adaptive = false)
+
+    sim = BucketSimulation(model, Y, (; domain = domain, soil_depth = d_soil), integrator, area_fraction)
+
+    # DSS state to ensure we have continuous fields
+    dss_state!(sim)
+    return sim
 end
 
 # extensions required by Interfacer
-get_field(sim::BucketSimulation, ::Val{:surface_temperature}) =
-    ClimaLand.surface_temperature(sim.model, sim.integrator.u, sim.integrator.p, sim.integrator.t)
-get_field(sim::BucketSimulation, ::Val{:surface_humidity}) =
-    ClimaLand.surface_specific_humidity(sim.model, sim.integrator.u, sim.integrator.p, sim.integrator.t)
-get_field(sim::BucketSimulation, ::Val{:roughness_momentum}) = sim.model.parameters.z_0m
-get_field(sim::BucketSimulation, ::Val{:roughness_buoyancy}) = sim.model.parameters.z_0b
+get_field(sim::BucketSimulation, ::Val{:air_density}) = sim.integrator.p.bucket.ρ_sfc
+get_field(sim::BucketSimulation, ::Val{:area_fraction}) = sim.area_fraction
 get_field(sim::BucketSimulation, ::Val{:beta}) =
     ClimaLand.surface_evaporative_scaling(sim.model, sim.integrator.u, sim.integrator.p)
+get_field(sim::BucketSimulation, ::Val{:roughness_buoyancy}) = sim.model.parameters.z_0b
+get_field(sim::BucketSimulation, ::Val{:roughness_momentum}) = sim.model.parameters.z_0m
 get_field(sim::BucketSimulation, ::Val{:surface_albedo}) =
     ClimaLand.surface_albedo(sim.model, sim.integrator.u, sim.integrator.p)
-get_field(sim::BucketSimulation, ::Val{:area_fraction}) = sim.area_fraction
-get_field(sim::BucketSimulation, ::Val{:air_density}) = sim.integrator.p.bucket.ρ_sfc
+get_field(sim::BucketSimulation, ::Val{:surface_humidity}) =
+    ClimaLand.surface_specific_humidity(sim.model, sim.integrator.u, sim.integrator.p, sim.integrator.t)
+get_field(sim::BucketSimulation, ::Val{:surface_temperature}) =
+    ClimaLand.surface_temperature(sim.model, sim.integrator.u, sim.integrator.p, sim.integrator.t)
 
-function update_field!(sim::BucketSimulation, ::Val{:turbulent_energy_flux}, field)
-    parent(sim.integrator.p.bucket.turbulent_fluxes.shf) .= parent(field)
+"""
+    get_field(bucket_sim::BucketSimulation, ::Val{:energy})
+
+Extension of Interfacer.get_field that provides the total energy contained in the bucket, including the latent heat due to snow melt.
+"""
+function get_field(bucket_sim::BucketSimulation, ::Val{:energy})
+    # required by ConservationChecker
+    e_per_area = zeros(axes(bucket_sim.integrator.u.bucket.W))
+    ClimaCore.Fields.bycolumn(axes(bucket_sim.integrator.u.bucket.T)) do colidx
+        e_per_area[colidx] .=
+            bucket_sim.model.parameters.ρc_soil .* mean(bucket_sim.integrator.u.bucket.T[colidx]) .*
+            bucket_sim.domain.soil_depth
+    end
+
+    e_per_area .+=
+        -LP.LH_f0(bucket_sim.model.parameters.earth_param_set) .*
+        LP.ρ_cloud_liq(bucket_sim.model.parameters.earth_param_set) .* bucket_sim.integrator.u.bucket.σS
+    return e_per_area
 end
-function update_field!(sim::BucketSimulation, ::Val{:turbulent_moisture_flux}, field)
-    ρ_liq = (LP.ρ_cloud_liq(sim.model.parameters.earth_param_set))
-    parent(sim.integrator.p.bucket.turbulent_fluxes.vapor_flux) .= parent(field ./ ρ_liq) # TODO: account for sublimation
+
+"""
+    get_field(bucket_sim::BucketSimulation, ::Val{:water})
+
+Extension of Interfacer.get_field that provides the total water contained in the bucket, including the liquid water in snow.
+"""
+function get_field(bucket_sim::BucketSimulation, ::Val{:water})
+    ρ_cloud_liq = ClimaLand.LP.ρ_cloud_liq(bucket_sim.model.parameters.earth_param_set)
+    return
+    @. (bucket_sim.integrator.u.bucket.σS + bucket_sim.integrator.u.bucket.W + bucket_sim.integrator.u.bucket.Ws) *
+       ρ_cloud_liq  # kg water / m2
 end
-function update_field!(sim::BucketSimulation, ::Val{:radiative_energy_flux_sfc}, field)
-    parent(sim.integrator.p.bucket.R_n) .= parent(field)
+
+function update_field!(sim::BucketSimulation, ::Val{:air_density}, field)
+    parent(sim.integrator.p.bucket.ρ_sfc) .= parent(field)
 end
 function update_field!(sim::BucketSimulation, ::Val{:liquid_precipitation}, field)
     ρ_liq = (LP.ρ_cloud_liq(sim.model.parameters.earth_param_set))
     parent(sim.integrator.p.drivers.P_liq) .= parent(field ./ ρ_liq)
 end
+function update_field!(sim::BucketSimulation, ::Val{:radiative_energy_flux_sfc}, field)
+    parent(sim.integrator.p.bucket.R_n) .= parent(field)
+end
+function update_field!(sim::BucketSimulation, ::Val{:turbulent_energy_flux}, field)
+    parent(sim.integrator.p.bucket.turbulent_fluxes.shf) .= parent(field)
+end
 function update_field!(sim::BucketSimulation, ::Val{:snow_precipitation}, field)
     ρ_ice = (LP.ρ_cloud_ice(sim.model.parameters.earth_param_set))
     parent(sim.integrator.p.drivers.P_snow) .= parent(field ./ ρ_ice)
 end
-
-function update_field!(sim::BucketSimulation, ::Val{:air_density}, field)
-    parent(sim.integrator.p.bucket.ρ_sfc) .= parent(field)
+function update_field!(sim::BucketSimulation, ::Val{:turbulent_moisture_flux}, field)
+    ρ_liq = (LP.ρ_cloud_liq(sim.model.parameters.earth_param_set))
+    parent(sim.integrator.p.bucket.turbulent_fluxes.vapor_flux) .= parent(field ./ ρ_liq) # TODO: account for sublimation
 end
 
+# extensions required by FieldExchanger
 step!(sim::BucketSimulation, t) = step!(sim.integrator, t - sim.integrator.t, true)
-
 reinit!(sim::BucketSimulation) = reinit!(sim.integrator)
 
 # extensions required by FluxCalculator (partitioned fluxes)
@@ -114,36 +246,87 @@ function get_model_prog_state(sim::BucketSimulation)
     return sim.integrator.u.bucket
 end
 
+###
+### ClimaLand.jl bucket model-specific functions (not explicitly required by ClimaCoupler.jl)
+###
+
+# TODO remove this function after ClimaLand v0.8.1 update
+function ClimaLand.turbulent_fluxes(atmos::CoupledAtmosphere, model::BucketModel, Y, p, t)
+    # coupler has done its thing behind the scenes already
+    model_name = ClimaLand.name(model)
+    model_cache = getproperty(p, model_name)
+    return model_cache.turbulent_fluxes
+end
+
+
+function ClimaLand.initialize_drivers(a::CoupledAtmosphere{FT}, coords) where {FT}
+    keys = (:P_liq, :P_snow)
+    types = ([FT for k in keys]...,)
+    domain_names = ([:surface for k in keys]...,)
+    model_name = :drivers
+    # intialize_vars packages the variables as a named tuple,
+    # as part of a named tuple with `model_name` as the key.
+    # Here we just want the variable named tuple itself
+    vars = ClimaLand.initialize_vars(keys, types, domain_names, coords, model_name)
+    return vars.drivers
+end
+
+
 """
-    get_field(bucket_sim::BucketSimulation, ::Val{:energy})
+    temp_anomaly_aquaplanet(coord)
 
-Extension of Interfacer.get_field that provides the total energy contained in the bucket, including the latent heat due to snow melt.
+Introduce a temperature IC anomaly for the aquaplanet case.
+The values for this case follow the moist Held-Suarez setup of Thatcher &
+Jablonowski (2016, eq. 6), consistent with ClimaAtmos aquaplanet.
 """
-function get_field(bucket_sim::BucketSimulation, ::Val{:energy})
-    # required by ConservationChecker
-    e_per_area = zeros(axes(bucket_sim.integrator.u.bucket.W))
-    ClimaCore.Fields.bycolumn(axes(bucket_sim.integrator.u.bucket.T)) do colidx
-        e_per_area[colidx] .=
-            bucket_sim.model.parameters.ρc_soil .* mean(bucket_sim.integrator.u.bucket.T[colidx]) .*
-            bucket_sim.domain.soil_depth
-    end
+temp_anomaly_aquaplanet(coord) = 29 * exp(-coord.lat^2 / (2 * 26^2))
 
-    e_per_area .+=
-        -LP.LH_f0(bucket_sim.model.parameters.earth_param_set) .*
-        LP.ρ_cloud_liq(bucket_sim.model.parameters.earth_param_set) .* bucket_sim.integrator.u.bucket.σS
-    return e_per_area
-end
+"""
+    temp_anomaly_amip(coord)
 
+Introduce a temperature IC anomaly for the AMIP case.
+The values used in this case have been tuned to align with observed temperature
+and result in stable simulations.
 """
-    get_field(bucket_sim::BucketSimulation, ::Val{:water})
+temp_anomaly_amip(coord) = 40 * cosd(coord.lat)^4
 
-Extension of Interfacer.get_field that provides the total water contained in the bucket, including the liquid water in snow.
 """
-function get_field(bucket_sim::BucketSimulation, ::Val{:water})
-    ρ_cloud_liq = ClimaLand.LP.ρ_cloud_liq(bucket_sim.model.parameters.earth_param_set)
-    return
-    @. (bucket_sim.integrator.u.bucket.σS + bucket_sim.integrator.u.bucket.W + bucket_sim.integrator.u.bucket.Ws) *
-       ρ_cloud_liq  # kg water / m2
+    make_land_domain(
+        atmos_boundary_space::ClimaCore.Spaces.SpectralElementSpace2D,
+        zlim::Tuple{FT, FT},
+        nelements_vert::Int,) where {FT}
+
+Creates the land model domain from the horizontal space of the atmosphere, and information
+about the number of elements and extent of the vertical domain.
+"""
+function make_land_domain(
+    atmos_boundary_space::ClimaCore.Spaces.SpectralElementSpace2D,
+    zlim::Tuple{FT, FT},
+    nelements_vert::Int,
+) where {FT}
+    @assert zlim[1] < zlim[2]
+    depth = zlim[2] - zlim[1]
+
+    mesh = ClimaCore.Spaces.topology(atmos_boundary_space).mesh
+
+    radius = mesh.domain.radius
+    nelements_horz = mesh.ne
+    npolynomial =
+        ClimaCore.Spaces.Quadratures.polynomial_degree(ClimaCore.Spaces.quadrature_style(atmos_boundary_space))
+    nelements = (nelements_horz, nelements_vert)
+    vertdomain = ClimaCore.Domains.IntervalDomain(
+        ClimaCore.Geometry.ZPoint(FT(zlim[1])),
+        ClimaCore.Geometry.ZPoint(FT(zlim[2]));
+        boundary_tags = (:bottom, :top),
+    )
+
+    vertmesh = ClimaCore.Meshes.IntervalMesh(vertdomain, ClimaCore.Meshes.Uniform(), nelems = nelements[2])
+    verttopology = ClimaCore.Topologies.IntervalTopology(vertmesh)
+    vert_center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(verttopology)
+    subsurface_space = ClimaCore.Spaces.ExtrudedFiniteDifferenceSpace(atmos_boundary_space, vert_center_space)
+    space = (; surface = atmos_boundary_space, subsurface = subsurface_space)
+
+    return ClimaLand.Domains.SphericalShell{FT}(radius, depth, nothing, nelements, npolynomial, space)
 end
 
 """
diff --git a/experiments/AMIP/components/ocean/eisenman_seaice_init.jl b/experiments/AMIP/components/ocean/eisenman_seaice.jl
similarity index 98%
rename from experiments/AMIP/components/ocean/eisenman_seaice_init.jl
rename to experiments/AMIP/components/ocean/eisenman_seaice.jl
index 6957196311..0bfe3bd678 100644
--- a/experiments/AMIP/components/ocean/eisenman_seaice_init.jl
+++ b/experiments/AMIP/components/ocean/eisenman_seaice.jl
@@ -10,6 +10,9 @@ import ClimaCoupler.FluxCalculator:
 import ClimaCoupler.Interfacer: get_field, update_field!
 import ClimaCoupler.FieldExchanger: step!, reinit!
 
+###
+### Functions required by ClimaCoupler.jl for a SurfaceModelSimulation
+###
 """
     EisenmanIceSimulation{P, Y, D, I}
 
@@ -45,6 +48,174 @@ Base.@kwdef struct EisenmanOceanParameters{FT <: AbstractFloat}
     α::FT = 0.1                     # albedo
 end
 
+"""
+    eisenman_seaice_init(::Type{FT}, tspan; space = nothing, area_fraction = nothing, thermo_params = nothing, stepper = CTS.RK4(), dt = 0.02, saveat = 1.0e10)
+
+Initialize the Eisenman-Zhang sea ice model and simulation.
+"""
+function eisenman_seaice_init(
+    ::Type{FT},
+    tspan;
+    space = nothing,
+    area_fraction = nothing,
+    thermo_params = nothing,
+    stepper = CTS.RK4(),
+    dt = 0.02,
+    saveat = 1.0e10,
+) where {FT}
+
+    params_ice = EisenmanIceParameters{FT}()
+    params_ocean = EisenmanOceanParameters{FT}()
+    params = (; p_i = params_ice, p_o = params_ocean)
+
+    # initiate prognostic variables
+    Y, Ya = state_init(params_ice, space)
+
+    ode_algo = CTS.ExplicitAlgorithm(stepper)
+    ode_function = CTS.ClimaODEFunction(T_exp! = ∑tendencies, dss! = weighted_dss_slab!)
+
+    cache = (;
+        Ya = Ya,
+        Δt = dt,
+        params = params,
+        area_fraction = area_fraction,
+        ice_area_fraction = zeros(space),
+        thermo_params = thermo_params,
+        dss_buffer = ClimaCore.Spaces.create_dss_buffer(ClimaCore.Fields.zeros(space)),
+    )
+    problem = ODEProblem(ode_function, Y, Float64.(tspan), cache)
+    integrator = init(problem, ode_algo, dt = Float64(dt), saveat = Float64(saveat), adaptive = false)
+
+    sim = EisenmanIceSimulation(params, Y, space, integrator)
+    @warn name(sim) *
+          " assumes gray radiation, no snow coverage, and PartitionedStateFluxes for the surface flux calculation."
+    return sim
+end
+
+# extensions required by Interfacer
+get_field(sim::EisenmanIceSimulation, ::Val{:air_density}) = sim.integrator.p.Ya.ρ_sfc
+get_field(sim::EisenmanIceSimulation, ::Val{:area_fraction}) = sim.integrator.p.area_fraction
+get_field(sim::EisenmanIceSimulation, ::Val{:beta}) = convert(eltype(sim.integrator.u), 1.0)
+get_field(sim::EisenmanIceSimulation, ::Val{:roughness_buoyancy}) =
+    @. sim.integrator.p.params.p_i.z0b * (sim.integrator.p.ice_area_fraction) +
+       sim.integrator.p.params.p_o.z0b .* (1 - sim.integrator.p.ice_area_fraction)
+get_field(sim::EisenmanIceSimulation, ::Val{:roughness_momentum}) =
+    @. sim.integrator.p.params.p_i.z0m * (sim.integrator.p.ice_area_fraction) +
+       sim.integrator.p.params.p_o.z0m .* (1 - sim.integrator.p.ice_area_fraction)
+get_field(sim::EisenmanIceSimulation, ::Val{:surface_albedo}) =
+    @. sim.integrator.p.params.p_i.α * (sim.integrator.p.ice_area_fraction) +
+       sim.integrator.p.params.p_o.α .* (1 - sim.integrator.p.ice_area_fraction)
+get_field(sim::EisenmanIceSimulation, ::Val{:surface_humidity}) = sim.integrator.u.q_sfc
+get_field(sim::EisenmanIceSimulation, ::Val{:surface_temperature}) = sim.integrator.u.T_sfc
+get_field(sim::EisenmanIceSimulation, ::Val{:water}) = nothing
+
+"""
+    get_field(sim::EisenmanIceSimulation, ::Val{:energy})
+
+Extension of Interfacer.get_field to get the energy of the ocean.
+It is the sum of the heat content of the mixed layer, the heat content of the ice, the heat flux from the ocean below ice.
+"""
+function get_field(sim::EisenmanIceSimulation, ::Val{:energy})
+    p_i = sim.integrator.p.params.p_i
+    p_o = sim.integrator.p.params.p_o
+    C0_base = p_i.C0_base
+    T_base = p_i.T_base
+    L_ice = p_i.L_ice
+    T_freeze = p_i.T_freeze
+    k_ice = p_i.k_ice
+    ocean_qflux = sim.integrator.p.Ya.ocean_qflux
+
+    cache = sim.integrator.p
+    Δt = cache.Δt
+    e_base = cache.Ya.e_base
+    ocean_qflux = cache.Ya.ocean_qflux
+
+    FT = eltype(sim.integrator.u)
+
+    hρc_ml = p_o.h * p_o.ρ * p_o.c
+
+    e_ml = @. p_o.h * p_o.ρ * p_o.c * sim.integrator.u.T_ml # heat
+    e_ice = @. p_i.L_ice * sim.integrator.u.h_ice # phase
+    e_qflux = @. ocean_qflux * FT(sim.integrator.t)
+
+    return @. e_ml + e_ice + e_qflux + e_base
+end
+
+function update_field!(sim::EisenmanIceSimulation, ::Val{:air_density}, field)
+    parent(sim.integrator.p.Ya.ρ_sfc) .= parent(field)
+end
+function update_field!(sim::EisenmanIceSimulation, ::Val{:area_fraction}, field::ClimaCore.Fields.Field)
+    sim.integrator.p.area_fraction .= field
+end
+function update_field!(sim::EisenmanIceSimulation, ::Val{:∂F_turb_energy∂T_sfc}, field, colidx)
+    sim.integrator.p.Ya.∂F_turb_energy∂T_sfc[colidx] .= field
+end
+function update_field!(sim::EisenmanIceSimulation, ::Val{:radiative_energy_flux_sfc}, field)
+    parent(sim.integrator.p.Ya.F_rad) .= parent(field)
+end
+function update_field!(sim::EisenmanIceSimulation, ::Val{:turbulent_energy_flux}, field)
+    parent(sim.integrator.p.Ya.F_turb) .= parent(field)
+end
+
+# extensions required by FieldExchanger
+step!(sim::EisenmanIceSimulation, t) = step!(sim.integrator, t - sim.integrator.t, true)
+reinit!(sim::EisenmanIceSimulation) = reinit!(sim.integrator)
+
+# extensions required by FluxCalculator (partitioned fluxes)
+function update_turbulent_fluxes_point!(
+    sim::EisenmanIceSimulation,
+    fields::NamedTuple,
+    colidx::ClimaCore.Fields.ColumnIndex,
+)
+    (; F_turb_energy) = fields
+    @. sim.integrator.p.Ya.F_turb[colidx] = F_turb_energy
+end
+
+"""
+    get_model_prog_state(sim::EisenmanIceSimulation)
+
+Extension of Checkpointer.get_model_prog_state to get the model state.
+"""
+function get_model_prog_state(sim::EisenmanIceSimulation)
+    return sim.integrator.u
+end
+
+"""
+    differentiate_turbulent_fluxes!(sim::EisenmanIceSimulation, args)
+
+Extension of differentiate_turbulent_fluxes! from FluxCalculator to get the turbulent fluxes.
+"""
+differentiate_turbulent_fluxes!(sim::EisenmanIceSimulation, args) =
+    differentiate_turbulent_fluxes!(sim::EisenmanIceSimulation, args..., ΔT_sfc = 0.1)
+
+"""
+    differentiate_turbulent_fluxes(sim::Interfacer.SurfaceModelSimulation, thermo_params, input_args, fluxes, δT_sfc = 0.1)
+
+Differentiates the turbulent fluxes in the surface model simulation `sim` with respect to the surface temperature,
+using δT_sfc as the perturbation.
+"""
+function differentiate_turbulent_fluxes!(sim::EisenmanIceSimulation, thermo_params, input_args, fluxes; δT_sfc = 0.1)
+    (; thermo_state_int, surface_params, surface_scheme, colidx) = input_args
+    thermo_state_sfc_dT = surface_thermo_state(sim, thermo_params, thermo_state_int, colidx, δT_sfc = δT_sfc)
+    input_args = merge(input_args, (; thermo_state_sfc = thermo_state_sfc_dT))
+
+    # set inputs based on whether the surface_scheme is `MoninObukhovScheme` or `BulkScheme`
+    inputs = surface_inputs(surface_scheme, input_args)
+
+    # calculate the surface fluxes
+    _, _, F_shf_δT_sfc, F_lhf_δT_sfc, _ = get_surface_fluxes_point!(inputs, surface_params)
+
+    (; F_shf, F_lhf) = fluxes
+
+    # calculate the derivative
+    ∂F_turb_energy∂T_sfc = @. (F_shf_δT_sfc + F_lhf_δT_sfc - F_shf - F_lhf) / δT_sfc
+
+    Interfacer.update_field!(sim, Val(:∂F_turb_energy∂T_sfc), ∂F_turb_energy∂T_sfc, colidx)
+end
+
+###
+### Eisenman-Zhang sea ice model-specific functions (not explicitly required by ClimaCoupler.jl)
+###
 """
     state_init(p::EisenmanIceParameters, space)
 
@@ -195,174 +366,4 @@ function ∑tendencies(dY, Y, cache, _)
 
     # update ice area fraction (binary mask for now)
     cache.ice_area_fraction .= Regridder.binary_mask.(Y.h_ice)
-
 end
-
-"""
-    eisenman_seaice_init(::Type{FT}, tspan; space = nothing, area_fraction = nothing, thermo_params = nothing, stepper = CTS.RK4(), dt = 0.02, saveat = 1.0e10)
-
-Initialize the Eisenman-Zhang sea ice model and simulation.
-"""
-function eisenman_seaice_init(
-    ::Type{FT},
-    tspan;
-    space = nothing,
-    area_fraction = nothing,
-    thermo_params = nothing,
-    stepper = CTS.RK4(),
-    dt = 0.02,
-    saveat = 1.0e10,
-) where {FT}
-
-    params_ice = EisenmanIceParameters{FT}()
-    params_ocean = EisenmanOceanParameters{FT}()
-    params = (; p_i = params_ice, p_o = params_ocean)
-
-    # initiate prognostic variables
-    Y, Ya = state_init(params_ice, space)
-
-    ode_algo = CTS.ExplicitAlgorithm(stepper)
-    ode_function = CTS.ClimaODEFunction(T_exp! = ∑tendencies, dss! = weighted_dss_slab!)
-
-    cache = (;
-        Ya = Ya,
-        Δt = dt,
-        params = params,
-        area_fraction = area_fraction,
-        ice_area_fraction = zeros(space),
-        thermo_params = thermo_params,
-        dss_buffer = ClimaCore.Spaces.create_dss_buffer(ClimaCore.Fields.zeros(space)),
-    )
-    problem = ODEProblem(ode_function, Y, Float64.(tspan), cache)
-    integrator = init(problem, ode_algo, dt = Float64(dt), saveat = Float64(saveat), adaptive = false)
-
-    sim = EisenmanIceSimulation(params, Y, space, integrator)
-    @warn name(sim) *
-          " assumes gray radiation, no snow coverage, and PartitionedStateFluxes for the surface flux calculation."
-    return sim
-end
-
-# extensions required by Interfacer
-get_field(sim::EisenmanIceSimulation, ::Val{:surface_temperature}) = sim.integrator.u.T_sfc
-get_field(sim::EisenmanIceSimulation, ::Val{:surface_humidity}) = sim.integrator.u.q_sfc
-get_field(sim::EisenmanIceSimulation, ::Val{:roughness_momentum}) =
-    @. sim.integrator.p.params.p_i.z0m * (sim.integrator.p.ice_area_fraction) +
-       sim.integrator.p.params.p_o.z0m .* (1 - sim.integrator.p.ice_area_fraction)
-get_field(sim::EisenmanIceSimulation, ::Val{:roughness_buoyancy}) =
-    @. sim.integrator.p.params.p_i.z0b * (sim.integrator.p.ice_area_fraction) +
-       sim.integrator.p.params.p_o.z0b .* (1 - sim.integrator.p.ice_area_fraction)
-get_field(sim::EisenmanIceSimulation, ::Val{:beta}) = convert(eltype(sim.integrator.u), 1.0)
-get_field(sim::EisenmanIceSimulation, ::Val{:surface_albedo}) =
-    @. sim.integrator.p.params.p_i.α * (sim.integrator.p.ice_area_fraction) +
-       sim.integrator.p.params.p_o.α .* (1 - sim.integrator.p.ice_area_fraction)
-get_field(sim::EisenmanIceSimulation, ::Val{:area_fraction}) = sim.integrator.p.area_fraction
-get_field(sim::EisenmanIceSimulation, ::Val{:air_density}) = sim.integrator.p.Ya.ρ_sfc
-
-function update_field!(sim::EisenmanIceSimulation, ::Val{:area_fraction}, field::ClimaCore.Fields.Field)
-    sim.integrator.p.area_fraction .= field
-end
-function update_field!(sim::EisenmanIceSimulation, ::Val{:turbulent_energy_flux}, field)
-    parent(sim.integrator.p.Ya.F_turb) .= parent(field)
-end
-function update_field!(sim::EisenmanIceSimulation, ::Val{:radiative_energy_flux_sfc}, field)
-    parent(sim.integrator.p.Ya.F_rad) .= parent(field)
-end
-function update_field!(sim::EisenmanIceSimulation, ::Val{:air_density}, field)
-    parent(sim.integrator.p.Ya.ρ_sfc) .= parent(field)
-end
-
-# extensions required by FieldExchanger
-step!(sim::EisenmanIceSimulation, t) = step!(sim.integrator, t - sim.integrator.t, true)
-reinit!(sim::EisenmanIceSimulation) = reinit!(sim.integrator)
-
-# extensions required by FluxCalculator (partitioned fluxes)
-function update_turbulent_fluxes_point!(
-    sim::EisenmanIceSimulation,
-    fields::NamedTuple,
-    colidx::ClimaCore.Fields.ColumnIndex,
-)
-    (; F_turb_energy) = fields
-    @. sim.integrator.p.Ya.F_turb[colidx] = F_turb_energy
-end
-
-"""
-    get_model_prog_state(sim::EisenmanIceSimulation)
-
-Extension of Checkpointer.get_model_prog_state to get the model state.
-"""
-function get_model_prog_state(sim::EisenmanIceSimulation)
-    return sim.integrator.u
-end
-
-"""
-    differentiate_turbulent_fluxes!(sim::EisenmanIceSimulation, args)
-
-Extension of differentiate_turbulent_fluxes! from FluxCalculator to get the turbulent fluxes.
-"""
-differentiate_turbulent_fluxes!(sim::EisenmanIceSimulation, args) =
-    differentiate_turbulent_fluxes!(sim::EisenmanIceSimulation, args..., ΔT_sfc = 0.1)
-
-"""
-    differentiate_turbulent_fluxes(sim::Interfacer.SurfaceModelSimulation, thermo_params, input_args, fluxes, δT_sfc = 0.1)
-
-Differentiates the turbulent fluxes in the surface model simulation `sim` with respect to the surface temperature,
-using δT_sfc as the perturbation.
-"""
-function differentiate_turbulent_fluxes!(sim::EisenmanIceSimulation, thermo_params, input_args, fluxes; δT_sfc = 0.1)
-    (; thermo_state_int, surface_params, surface_scheme, colidx) = input_args
-    thermo_state_sfc_dT = surface_thermo_state(sim, thermo_params, thermo_state_int, colidx, δT_sfc = δT_sfc)
-    input_args = merge(input_args, (; thermo_state_sfc = thermo_state_sfc_dT))
-
-    # set inputs based on whether the surface_scheme is `MoninObukhovScheme` or `BulkScheme`
-    inputs = surface_inputs(surface_scheme, input_args)
-
-    # calculate the surface fluxes
-    _, _, F_shf_δT_sfc, F_lhf_δT_sfc, _ = get_surface_fluxes_point!(inputs, surface_params)
-
-    (; F_shf, F_lhf) = fluxes
-
-    # calculate the derivative
-    ∂F_turb_energy∂T_sfc = @. (F_shf_δT_sfc + F_lhf_δT_sfc - F_shf - F_lhf) / δT_sfc
-
-    Interfacer.update_field!(sim, Val(:∂F_turb_energy∂T_sfc), ∂F_turb_energy∂T_sfc, colidx)
-end
-
-function update_field!(sim::EisenmanIceSimulation, ::Val{:∂F_turb_energy∂T_sfc}, field, colidx)
-    sim.integrator.p.Ya.∂F_turb_energy∂T_sfc[colidx] .= field
-end
-
-
-"""
-    get_field(sim::EisenmanIceSimulation, ::Val{:energy})
-
-Extension of Interfacer.get_field to get the energy of the ocean.
-It is the sum of the heat content of the mixed layer, the heat content of the ice, the heat flux from the ocean below ice.
-"""
-function get_field(sim::EisenmanIceSimulation, ::Val{:energy})
-    p_i = sim.integrator.p.params.p_i
-    p_o = sim.integrator.p.params.p_o
-    C0_base = p_i.C0_base
-    T_base = p_i.T_base
-    L_ice = p_i.L_ice
-    T_freeze = p_i.T_freeze
-    k_ice = p_i.k_ice
-    ocean_qflux = sim.integrator.p.Ya.ocean_qflux
-
-    cache = sim.integrator.p
-    Δt = cache.Δt
-    e_base = cache.Ya.e_base
-    ocean_qflux = cache.Ya.ocean_qflux
-
-    FT = eltype(sim.integrator.u)
-
-    hρc_ml = p_o.h * p_o.ρ * p_o.c
-
-    e_ml = @. p_o.h * p_o.ρ * p_o.c * sim.integrator.u.T_ml # heat
-    e_ice = @. p_i.L_ice * sim.integrator.u.h_ice # phase
-    e_qflux = @. ocean_qflux * FT(sim.integrator.t)
-
-    return @. e_ml + e_ice + e_qflux + e_base
-
-end
-
-get_field(sim::EisenmanIceSimulation, ::Val{:water}) = nothing
diff --git a/experiments/AMIP/components/ocean/prescr_seaice_init.jl b/experiments/AMIP/components/ocean/prescr_seaice.jl
similarity index 98%
rename from experiments/AMIP/components/ocean/prescr_seaice_init.jl
rename to experiments/AMIP/components/ocean/prescr_seaice.jl
index f7c51cc141..88d499b248 100644
--- a/experiments/AMIP/components/ocean/prescr_seaice_init.jl
+++ b/experiments/AMIP/components/ocean/prescr_seaice.jl
@@ -13,6 +13,9 @@ import ClimaCoupler.BCReader: float_type_bcf
 
 include("../slab_utils.jl")
 
+###
+### Functions required by ClimaCoupler.jl for a SurfaceModelSimulation
+###
 """
     PrescribedIceSimulation{P, Y, D, I}
 
@@ -61,38 +64,6 @@ function slab_ice_space_init(::Type{FT}, space, p) where {FT}
     return Y
 end
 
-"""
-    ice_rhs!(du, u, p, _)
-
-Rhs method in the form as required by `ClimeTimeSteppers`, with the tendency vector `du`,
-the state vector `u` and the parameter vector, `p`, as input arguments.
-
-This sea-ice energy formulation follows [Holloway and Manabe 1971](https://journals.ametsoc.org/view/journals/mwre/99/5/1520-0493_1971_099_0335_socbag_2_3_co_2.xml?tab_body=pdf),
-where sea-ice concentrations and thicknes are prescribed, and the model solves for temperature (curbed at the freezing point).
-"""
-function ice_rhs!(du, u, p, _)
-    dY = du
-    Y = u
-    FT = eltype(dY)
-
-    params = p.params
-    F_turb_energy = p.F_turb_energy
-    F_radiative = p.F_radiative
-    area_fraction = p.area_fraction
-    T_freeze = params.T_freeze
-
-    F_conductive = @. params.k_ice / (params.h) * (params.T_base - Y.T_sfc) # fluxes are defined to be positive when upward
-    rhs = @. (-F_turb_energy - F_radiative + F_conductive) / (params.h * params.ρ * params.c)
-
-    # do not count tendencies that lead to temperatures above freezing, and mask out no-ice areas
-    area_mask = Regridder.binary_mask.(area_fraction)
-    threshold = zero(FT)
-    unphysical = @. Regridder.binary_mask.(T_freeze - (Y.T_sfc + FT(rhs) * FT(p.dt)), threshold) .* area_mask
-    parent(dY.T_sfc) .= parent(rhs .* unphysical)
-
-    @. p.q_sfc = TD.q_vap_saturation_generic.(p.thermo_params, Y.T_sfc, p.ρ_sfc, TD.Ice())
-end
-
 """
     ice_init(::Type{FT}; tspan, dt, saveat, space, ice_fraction, stepper = CTS.RK4()) where {FT}
 
@@ -128,39 +99,37 @@ function ice_init(::Type{FT}; tspan, saveat, dt, space, area_fraction, thermo_pa
     return sim
 end
 
-# file-specific
-"""
-    scale_sic(SIC, _info)
-Ensures that the space of the SIC struct matches that of the mask, and converts the units from area % to area fraction.
-"""
-scale_sic(SIC, _info) = swap_space!(zeros(axes(_info.land_fraction)), SIC) ./ float_type_bcf(_info)(100.0)
-
-# setting that SIC < 0.5 is counted as ocean if binary remapping.
-get_ice_fraction(h_ice::FT, mono::Bool, threshold = 0.5) where {FT} =
-    mono ? h_ice : Regridder.binary_mask(h_ice, threshold)
-
 # extensions required by Interfacer
-get_field(sim::PrescribedIceSimulation, ::Val{:surface_temperature}) = sim.integrator.u.T_sfc
-get_field(sim::PrescribedIceSimulation, ::Val{:surface_humidity}) = sim.integrator.p.q_sfc
-get_field(sim::PrescribedIceSimulation, ::Val{:roughness_momentum}) = sim.integrator.p.params.z0m
-get_field(sim::PrescribedIceSimulation, ::Val{:roughness_buoyancy}) = sim.integrator.p.params.z0b
+get_field(sim::PrescribedIceSimulation, ::Val{:air_density}) = sim.integrator.p.ρ_sfc
+get_field(sim::PrescribedIceSimulation, ::Val{:area_fraction}) = sim.integrator.p.area_fraction
 get_field(sim::PrescribedIceSimulation, ::Val{:beta}) = convert(eltype(sim.integrator.u), 1.0)
+get_field(sim::PrescribedIceSimulation, ::Val{:roughness_buoyancy}) = sim.integrator.p.params.z0b
+get_field(sim::PrescribedIceSimulation, ::Val{:roughness_momentum}) = sim.integrator.p.params.z0m
 get_field(sim::PrescribedIceSimulation, ::Val{:surface_albedo}) = sim.integrator.p.params.α
-get_field(sim::PrescribedIceSimulation, ::Val{:area_fraction}) = sim.integrator.p.area_fraction
-get_field(sim::PrescribedIceSimulation, ::Val{:air_density}) = sim.integrator.p.ρ_sfc
+get_field(sim::PrescribedIceSimulation, ::Val{:surface_humidity}) = sim.integrator.p.q_sfc
+get_field(sim::PrescribedIceSimulation, ::Val{:surface_temperature}) = sim.integrator.u.T_sfc
+get_field(sim::PrescribedIceSimulation, ::Val{:water}) = nothing
+
+"""
+    get_field(sim::PrescribedIceSimulation, ::Val{:energy})
+
+Extension of Interfacer.get_field to get the energy of the ocean.
+It multiplies the the slab temperature by the heat capacity, density, and depth.
+"""
+get_field(sim::PrescribedIceSimulation, ::Val{:energy}) =
+    sim.integrator.p.params.ρ .* sim.integrator.p.params.c .* sim.integrator.u.T_sfc .* sim.integrator.p.params.h
 
+function update_field!(sim::PrescribedIceSimulation, ::Val{:air_density}, field)
+    parent(sim.integrator.p.ρ_sfc) .= parent(field)
+end
 function update_field!(sim::PrescribedIceSimulation, ::Val{:area_fraction}, field::ClimaCore.Fields.Field)
     sim.integrator.p.area_fraction .= field
 end
-
-function update_field!(sim::PrescribedIceSimulation, ::Val{:turbulent_energy_flux}, field)
-    parent(sim.integrator.p.F_turb_energy) .= parent(field)
-end
 function update_field!(sim::PrescribedIceSimulation, ::Val{:radiative_energy_flux_sfc}, field)
     parent(sim.integrator.p.F_radiative) .= parent(field)
 end
-function update_field!(sim::PrescribedIceSimulation, ::Val{:air_density}, field)
-    parent(sim.integrator.p.ρ_sfc) .= parent(field)
+function update_field!(sim::PrescribedIceSimulation, ::Val{:turbulent_energy_flux}, field)
+    parent(sim.integrator.p.F_turb_energy) .= parent(field)
 end
 
 # extensions required by FieldExchanger
@@ -186,16 +155,50 @@ function get_model_prog_state(sim::PrescribedIceSimulation)
     return sim.integrator.u
 end
 
+###
+### Sea ice model-specific functions (not explicitly required by ClimaCoupler.jl)
+###
 """
-    get_field(sim::PrescribedIceSimulation, ::Val{:energy})
+    scale_sic(SIC, _info)
+Ensures that the space of the SIC struct matches that of the mask, and converts the units from area % to area fraction.
+"""
+scale_sic(SIC, _info) = swap_space!(zeros(axes(_info.land_fraction)), SIC) ./ float_type_bcf(_info)(100.0)
+
+# setting that SIC < 0.5 is counted as ocean if binary remapping.
+get_ice_fraction(h_ice::FT, mono::Bool, threshold = 0.5) where {FT} =
+    mono ? h_ice : Regridder.binary_mask(h_ice, threshold)
 
-Extension of Interfacer.get_field to get the energy of the ocean.
-It multiplies the the slab temperature by the heat capacity, density, and depth.
 """
-get_field(sim::PrescribedIceSimulation, ::Val{:energy}) =
-    sim.integrator.p.params.ρ .* sim.integrator.p.params.c .* sim.integrator.u.T_sfc .* sim.integrator.p.params.h
+    ice_rhs!(du, u, p, _)
 
-get_field(sim::PrescribedIceSimulation, ::Val{:water}) = nothing
+Rhs method in the form as required by `ClimeTimeSteppers`, with the tendency vector `du`,
+the state vector `u` and the parameter vector, `p`, as input arguments.
+
+This sea-ice energy formulation follows [Holloway and Manabe 1971](https://journals.ametsoc.org/view/journals/mwre/99/5/1520-0493_1971_099_0335_socbag_2_3_co_2.xml?tab_body=pdf),
+where sea-ice concentrations and thicknes are prescribed, and the model solves for temperature (curbed at the freezing point).
+"""
+function ice_rhs!(du, u, p, _)
+    dY = du
+    Y = u
+    FT = eltype(dY)
+
+    params = p.params
+    F_turb_energy = p.F_turb_energy
+    F_radiative = p.F_radiative
+    area_fraction = p.area_fraction
+    T_freeze = params.T_freeze
+
+    F_conductive = @. params.k_ice / (params.h) * (params.T_base - Y.T_sfc) # fluxes are defined to be positive when upward
+    rhs = @. (-F_turb_energy - F_radiative + F_conductive) / (params.h * params.ρ * params.c)
+
+    # do not count tendencies that lead to temperatures above freezing, and mask out no-ice areas
+    area_mask = Regridder.binary_mask.(area_fraction)
+    threshold = zero(FT)
+    unphysical = @. Regridder.binary_mask.(T_freeze - (Y.T_sfc + FT(rhs) * FT(p.dt)), threshold) .* area_mask
+    parent(dY.T_sfc) .= parent(rhs .* unphysical)
+
+    @. p.q_sfc = TD.q_vap_saturation_generic.(p.thermo_params, Y.T_sfc, p.ρ_sfc, TD.Ice())
+end
 
 """
     dss_state!(sim::PrescribedIceSimulation)
diff --git a/experiments/AMIP/components/ocean/slab_ocean_init.jl b/experiments/AMIP/components/ocean/slab_ocean.jl
similarity index 97%
rename from experiments/AMIP/components/ocean/slab_ocean_init.jl
rename to experiments/AMIP/components/ocean/slab_ocean.jl
index 601ce1918c..2a708bb52b 100644
--- a/experiments/AMIP/components/ocean/slab_ocean_init.jl
+++ b/experiments/AMIP/components/ocean/slab_ocean.jl
@@ -10,6 +10,9 @@ import ClimaCoupler.BCReader: float_type_bcf
 
 include("../slab_utils.jl")
 
+###
+### Functions required by ClimaCoupler.jl for a SurfaceModelSimulation
+###
 """
     SlabOceanSimulation{P, Y, D, I}
 
@@ -40,18 +43,6 @@ end
 
 name(::SlabOceanSimulation) = "SlabOceanSimulation"
 
-"""
-    temp_anomaly(coord)
-
-Calculates a an anomaly to be added to the initial condition for temperature.
-This default case includes only an anomaly at the tropics.
-"""
-function temp_anomaly(coord)
-    # include tropics anomaly
-    anom = FT(29 * exp(-coord.lat^2 / (2 * 26^2)))
-    return anom
-end
-
 """
     slab_ocean_space_init(space, params)
 
@@ -72,15 +63,6 @@ function slab_ocean_space_init(space, params)
     return Y, space
 end
 
-# ode
-function slab_ocean_rhs!(dY, Y, cache, t)
-    p, F_turb_energy, F_radiative, area_fraction = cache
-    FT = eltype(Y.T_sfc)
-    rhs = @. -(F_turb_energy + F_radiative) / (p.h * p.ρ * p.c)
-    @. dY.T_sfc = rhs * Regridder.binary_mask(area_fraction) * p.evolving_switch
-    @. cache.q_sfc = TD.q_vap_saturation_generic.(cache.thermo_params, Y.T_sfc, cache.ρ_sfc, TD.Liquid())
-end
-
 """
     ocean_init(::Type{FT}; tspan, dt, saveat, space, area_fraction, stepper = CTS.RK4()) where {FT}
 
@@ -127,40 +109,41 @@ function ocean_init(
     return sim
 end
 
-# file specific
-"""
-    scale_sst(SST::FT, _info)
-Ensures that the space of the SST struct matches that of the land_fraction, and converts the units to Kelvin (N.B.: this is dataset specific)
-"""
-scale_sst(SST, _info) = (swap_space!(zeros(axes(_info.land_fraction)), SST) .+ float_type_bcf(_info)(273.15))
-
 # extensions required by Interfacer
-get_field(sim::SlabOceanSimulation, ::Val{:surface_temperature}) = sim.integrator.u.T_sfc
-get_field(sim::SlabOceanSimulation, ::Val{:surface_humidity}) = sim.integrator.p.q_sfc
-get_field(sim::SlabOceanSimulation, ::Val{:roughness_momentum}) = sim.integrator.p.params.z0m
-get_field(sim::SlabOceanSimulation, ::Val{:roughness_buoyancy}) = sim.integrator.p.params.z0b
+get_field(sim::SlabOceanSimulation, ::Val{:air_density}) = sim.integrator.p.ρ_sfc
+get_field(sim::SlabOceanSimulation, ::Val{:area_fraction}) = sim.integrator.p.area_fraction
 get_field(sim::SlabOceanSimulation, ::Val{:beta}) = convert(eltype(sim.integrator.u), 1.0)
+get_field(sim::SlabOceanSimulation, ::Val{:roughness_buoyancy}) = sim.integrator.p.params.z0b
+get_field(sim::SlabOceanSimulation, ::Val{:roughness_momentum}) = sim.integrator.p.params.z0m
 get_field(sim::SlabOceanSimulation, ::Val{:surface_albedo}) = sim.integrator.p.params.α
-get_field(sim::SlabOceanSimulation, ::Val{:area_fraction}) = sim.integrator.p.area_fraction
-get_field(sim::SlabOceanSimulation, ::Val{:air_density}) = sim.integrator.p.ρ_sfc
+get_field(sim::SlabOceanSimulation, ::Val{:surface_humidity}) = sim.integrator.p.q_sfc
+get_field(sim::SlabOceanSimulation, ::Val{:surface_temperature}) = sim.integrator.u.T_sfc
+get_field(sim::SlabOceanSimulation, ::Val{:water}) = nothing
+
+"""
+    get_field(sim::SlabOceanSimulation, ::Val{:energy})
+
+Extension of Interfacer.get_field to get the energy of the ocean.
+It multiplies the the slab temperature by the heat capacity, density, and depth.
+"""
+get_field(sim::SlabOceanSimulation, ::Val{:energy}) =
+    sim.integrator.p.params.ρ .* sim.integrator.p.params.c .* sim.integrator.u.T_sfc .* sim.integrator.p.params.h
 
 function update_field!(sim::SlabOceanSimulation, ::Val{:area_fraction}, field::ClimaCore.Fields.Field)
     sim.integrator.p.area_fraction .= field
 end
-
-function update_field!(sim::SlabOceanSimulation, ::Val{:turbulent_energy_flux}, field)
-    parent(sim.integrator.p.F_turb_energy) .= parent(field)
+function update_field!(sim::SlabOceanSimulation, ::Val{:air_density}, field)
+    parent(sim.integrator.p.ρ_sfc) .= parent(field)
 end
 function update_field!(sim::SlabOceanSimulation, ::Val{:radiative_energy_flux_sfc}, field)
     parent(sim.integrator.p.F_radiative) .= parent(field)
 end
-function update_field!(sim::SlabOceanSimulation, ::Val{:air_density}, field)
-    parent(sim.integrator.p.ρ_sfc) .= parent(field)
+function update_field!(sim::SlabOceanSimulation, ::Val{:turbulent_energy_flux}, field)
+    parent(sim.integrator.p.F_turb_energy) .= parent(field)
 end
 
 # extensions required by FieldExchanger
 step!(sim::SlabOceanSimulation, t) = step!(sim.integrator, t - sim.integrator.t, true)
-
 reinit!(sim::SlabOceanSimulation) = reinit!(sim.integrator)
 
 # extensions required by FluxCalculator (partitioned fluxes)
@@ -182,16 +165,35 @@ function get_model_prog_state(sim::SlabOceanSimulation)
     return sim.integrator.u
 end
 
+###
+### Slab ocean model-specific functions (not explicitly required by ClimaCoupler.jl)
+###
 """
-    get_field(sim::SlabOceanSimulation, ::Val{:energy})
+    scale_sst(SST::FT, _info)
+Ensures that the space of the SST struct matches that of the land_fraction, and converts the units to Kelvin (N.B.: this is dataset specific)
+"""
+scale_sst(SST, _info) = (swap_space!(zeros(axes(_info.land_fraction)), SST) .+ float_type_bcf(_info)(273.15))
+
+# ode
+function slab_ocean_rhs!(dY, Y, cache, t)
+    p, F_turb_energy, F_radiative, area_fraction = cache
+    FT = eltype(Y.T_sfc)
+    rhs = @. -(F_turb_energy + F_radiative) / (p.h * p.ρ * p.c)
+    @. dY.T_sfc = rhs * Regridder.binary_mask(area_fraction) * p.evolving_switch
+    @. cache.q_sfc = TD.q_vap_saturation_generic.(cache.thermo_params, Y.T_sfc, cache.ρ_sfc, TD.Liquid())
+end
 
-Extension of Interfacer.get_field to get the energy of the ocean.
-It multiplies the the slab temperature by the heat capacity, density, and depth.
 """
-get_field(sim::SlabOceanSimulation, ::Val{:energy}) =
-    sim.integrator.p.params.ρ .* sim.integrator.p.params.c .* sim.integrator.u.T_sfc .* sim.integrator.p.params.h
+    temp_anomaly(coord)
 
-get_field(sim::SlabOceanSimulation, ::Val{:water}) = nothing
+Calculates a an anomaly to be added to the initial condition for temperature.
+This default case includes only an anomaly at the tropics.
+"""
+function temp_anomaly(coord)
+    # include tropics anomaly
+    anom = FT(29 * exp(-coord.lat^2 / (2 * 26^2)))
+    return anom
+end
 
 """
     dss_state!(sim::SlabOceanSimulation)
@@ -207,4 +209,5 @@ function dss_state!(sim::SlabOceanSimulation)
         buffer = get_dss_buffer(axes(field), p)
         ClimaCore.Spaces.weighted_dss!(field, buffer)
     end
+
 end
diff --git a/experiments/AMIP/coupler_driver.jl b/experiments/AMIP/coupler_driver.jl
index b0a833d96f..5a89f5401e 100644
--- a/experiments/AMIP/coupler_driver.jl
+++ b/experiments/AMIP/coupler_driver.jl
@@ -75,12 +75,11 @@ contain any internals of the ClimaCoupler source code, except extensions to the
 =#
 
 ## helpers for component models
-include("components/atmosphere/climaatmos_init.jl")
-include("components/land/bucket_init.jl")
-include("components/land/bucket_utils.jl")
-include("components/ocean/slab_ocean_init.jl")
-include("components/ocean/prescr_seaice_init.jl")
-include("components/ocean/eisenman_seaice_init.jl")
+include("components/atmosphere/climaatmos.jl")
+include("components/land/climaland_bucket.jl")
+include("components/ocean/slab_ocean.jl")
+include("components/ocean/prescr_seaice.jl")
+include("components/ocean/eisenman_seaice.jl")
 
 ## helpers for user-specified IO
 include("user_io/user_diagnostics.jl")
diff --git a/test/component_model_tests/bucket_tests.jl b/test/component_model_tests/bucket_tests.jl
index 0791a8565a..e4f5537553 100644
--- a/test/component_model_tests/bucket_tests.jl
+++ b/test/component_model_tests/bucket_tests.jl
@@ -3,7 +3,7 @@ import ClimaCoupler
 using ClimaCoupler.TestHelper: create_space
 using ClimaCore: Fields, Spaces
 
-include(pkgdir(ClimaCoupler, "experiments/AMIP/components/land/bucket_init.jl"))
+include(pkgdir(ClimaCoupler, "experiments/AMIP/components/land/climaland_bucket.jl"))
 
 for FT in (Float32, Float64)
     @testset "dss_state! BucketSimulation for FT=$FT" begin
diff --git a/test/component_model_tests/climaatmos_tests.jl b/test/component_model_tests/climaatmos_tests.jl
index 2a3cce6272..ca1ad2a4d3 100644
--- a/test/component_model_tests/climaatmos_tests.jl
+++ b/test/component_model_tests/climaatmos_tests.jl
@@ -4,7 +4,7 @@ using ClimaCoupler.Interfacer: AtmosModelSimulation
 using ClimaCoupler.TestHelper: create_space
 using ClimaCore: Fields, Spaces
 
-include(pkgdir(ClimaCoupler, "experiments/AMIP/components/atmosphere/climaatmos_init.jl"))
+include(pkgdir(ClimaCoupler, "experiments/AMIP/components/atmosphere/climaatmos.jl"))
 
 for FT in (Float32, Float64)
     @testset "dss_state! ClimaAtmosSimulation for FT=$FT" begin
diff --git a/test/component_model_tests/eisenman_seaice_tests.jl b/test/component_model_tests/eisenman_seaice_tests.jl
index ca7805023c..5141982626 100644
--- a/test/component_model_tests/eisenman_seaice_tests.jl
+++ b/test/component_model_tests/eisenman_seaice_tests.jl
@@ -12,7 +12,7 @@ import CLIMAParameters as CP
 import Thermodynamics as TD
 import Thermodynamics.Parameters as TP
 
-include("../../experiments/AMIP/components/ocean/eisenman_seaice_init.jl")
+include("../../experiments/AMIP/components/ocean/eisenman_seaice.jl")
 
 
 for FT in (Float32, Float64)
diff --git a/test/component_model_tests/prescr_seaice_tests.jl b/test/component_model_tests/prescr_seaice_tests.jl
index 3c7ab7bf33..5a4c8abcfd 100644
--- a/test/component_model_tests/prescr_seaice_tests.jl
+++ b/test/component_model_tests/prescr_seaice_tests.jl
@@ -7,7 +7,7 @@ using ClimaCore: Fields, Spaces
 import CLIMAParameters as CP
 import Thermodynamics.Parameters as TDP
 
-include(pkgdir(ClimaCoupler, "experiments/AMIP/components/ocean/prescr_seaice_init.jl"))
+include(pkgdir(ClimaCoupler, "experiments/AMIP/components/ocean/prescr_seaice.jl"))
 
 for FT in (Float32, Float64)
     @testset "test sea-ice energy slab for FT=$FT" begin
diff --git a/test/component_model_tests/slab_ocean_tests.jl b/test/component_model_tests/slab_ocean_tests.jl
index 87c48bc62b..7a02f54584 100644
--- a/test/component_model_tests/slab_ocean_tests.jl
+++ b/test/component_model_tests/slab_ocean_tests.jl
@@ -7,7 +7,7 @@ using ClimaCore: Fields, Spaces
 import CLIMAParameters as CP
 import Thermodynamics.Parameters as TDP
 
-include(pkgdir(ClimaCoupler, "experiments/AMIP/components/ocean/slab_ocean_init.jl"))
+include(pkgdir(ClimaCoupler, "experiments/AMIP/components/ocean/slab_ocean.jl"))
 
 for FT in (Float32, Float64)
     @testset "dss_state! SlabOceanSimulation for FT=$FT" begin