diff --git a/docs/src/fieldexchanger.md b/docs/src/fieldexchanger.md
index f37d08a41a..663770cc25 100644
--- a/docs/src/fieldexchanger.md
+++ b/docs/src/fieldexchanger.md
@@ -9,7 +9,7 @@ The `FieldExchanger` needs to populate the coupler with
 The component models are updated by broadcasting the coupler fields, via the `update_model_sims!` function. For an update, this function requires that `update_field!` is defined for the particular variable and component model. Currently, we support the:
 - `AtmosModelSimulation`: `albedo`, `surface_temperature`
     - if calculating fluxes in the atmospheric model: `roughness_momentum`, `roughness_buoyancy`, `beta`
-- `SurfaceModelSimulation`: `air_density`, `turbulent_energy_flux`, `turbulent_moisture_flux`, `radiative_energy_flux`, `liquid_precipitation`, `snow_precipitation`
+- `SurfaceModelSimulation`: `air_density`, `turbulent_energy_flux`, `turbulent_moisture_flux`, `radiative_energy_flux_sfc`, `liquid_precipitation`, `snow_precipitation`
 
 If an `update_field!` function is not defined for a particular component model, it will be ignored.
 
diff --git a/docs/src/interfacer.md b/docs/src/interfacer.md
index 133d7e35aa..bdbad74258 100644
--- a/docs/src/interfacer.md
+++ b/docs/src/interfacer.md
@@ -1,39 +1,182 @@
 # Interfacer
+This module contains functions that define the interface for coupling
+component models, as well as stub objects that contain prescribed fields.
+Here we explain each type of component model, and the functions that must
+be implemented to use a component model with ClimaCoupler.jl
 
-This module contains functions for defining the interface for coupling component models, as well as stub objects that contain prescribed fields.
-
-## Coupled Simulation
-- `CoupledSimulation` (`cs`) stores info for ESM run. We require that each `cs` contains four (`atmos_sim`, `land_sim`, `ocean_sim` and `ice_sim`) components. While this requirement will not be eventually needed, for the time being, if a simulation surface type is not needed for a given run, it should be initialized with `SurfaceStub` with a zero `area_fracion`. The `atmos_sim` should always be specified.
-
-## Component model simulations
-- all Simulations that are not the `CoupledSimulation` fall under `ComponentModelSimulation`
-- the current version requires that there is:
-    - one `AtmosModelSimulation`
-    - one or more `SurfaceModelSimulation`s, which require the following adapter functions:
-        ```
-        get_field(sim::SurfaceModelSimulation, ::Val{:area_fraction}) = ...
-        get_field(sim::SurfaceModelSimulation, ::Val{:surface_temperature}) = ...
-        get_field(sim::SurfaceModelSimulation, ::Val{:albedo}) = ...
-        get_field(sim::SurfaceModelSimulation, ::Val{:roughness_momentum}) = ...
-        get_field(sim::SurfaceModelSimulation, ::Val{:roughness_buoyancy}) = ...
-        get_field(sim::SurfaceModelSimulation, ::Val{:beta}) = ...
-        update_field!(sim::SurfaceModelSimulation, ::Val{:area_fraction}, field::Fields.Field) = ...
-        update_field!(sim::SurfaceModelSimulation, ::Val{:surface_temperature}, field::Fields.Field) = ...
-        ```
-        - these adapter functions, to be defined in the component models' init files (preferably in their own repositories), allow the coupler to operate without having to assume particular data structures of the underlying component models. This allows easy swapping of model components, as well as a stable source code with coupler-specific unit tests.
-
-## Prescribed conditions
-- `SurfaceStub` is a `SurfaceModelSimulation`, but it only contains required data in `<surface_stub>.cache`, e.g., for the calculation of surface fluxes through a prescribed surface state.  The above adapter functions are already predefined for `SurfaceStub`, with the cache variables specified as:
+## Coupled simulation
+A `CoupledSimulation` stores info for ESM run and contains each of the
+component model simulations. We currently require that each `CoupledSimulation`
+contains four components: `atmos_sim`, `land_sim`, `ocean_sim` and `ice_sim`.
+If a simulation surface type is not needed for a given run, it should be
+initialized with `SurfaceStub` with a zero `area_fracion`.
+The `atmos_sim` should always be specified.
+
+
+## Component simulations
+Individual component model simulations fall under `ComponentModelSimulation`,
+which together combine to make the `CoupledSimulation`.
+We have two types of `ComponentModelSimulations`: `AtmosModelSimulation` and
+`SurfaceModelSimulation`. The two have different requirements,
+which are detailed below. `SurfaceModelSimulation` is further divided into
+`SeaIceModelSimulation`, `LandModelSimulation`, and `OceanModelSimulation`,
+representing the 3 currently-supported options for surface models.
+
+### ComponentModelSimulation - required functions
+A component model simulation should be implemented as a struct that is a concrete subtype
+of a `ComponentModelSimulation`. This struct should contain all of the
+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
+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.
+- `init`: construct and return an instance of the `ComponentModelSimulation`,
+and perform all initialization. This function should return a simulation that
+is ready to be stepped in the coupled simulation. The interface for this
+function varies across component models.
+- `name(::ComponentModelSimulation)`: return a string containing the name of
+this `ComponentModelSimulation`, which is used for printing information about
+component models and writing to checkpoint files.
+- `step!(::ComponentModelSimulation, t)`: A function to update the
+simulation in-place with values calculate for time `t`. For the
+models we currently have implemented, this is a simple wrapper around
+the `step!` function implemented in SciMLBase.jl.
+- `reinit!(::ComponentModelSimulation)`: A function to restart a simulation
+after solving of the simulation has been paused or interrupted. Like
+`step!`, this is currently a simple wrapper around the `reinit!` function
+of SciMLBase.jl.
+- `get_model_state_vector(::ComponentModelSimulation)`: A function that
+returns the state vector of the simulation at its current state. This
+is used for checkpointing the simulation.
+
+### ComponentModelSimulation - optional functions
+- `update_sim!(::ComponentModelSimulation, csf, turbulent_fluxes)`: A
+function to update each of the fields of the component model simulation
+that are updated by the coupler. ClimaCoupler.jl provides defaults of
+this function for both `AtmosModelSimulation` and
+`SurfaceModelSimulation` that update each of the fields expected by
+the coupler. This function can optionally be extended to include
+additional field updates as desired.
+
+### AtmosModelSimulation - required functions
+In addition to the functions required for a general
+`ComponentModelSimulation`, an `AtmosModelSimulation` requires the
+following functions to retrieve and update its fields.
+- `get_field(::AtmosModelSimulation. ::Val{property})`: This getter
+function returns the value of the field property for the simulation
+in its current state. For an `AtmosModelSimulation`, it must be extended
+for the following properties:
+
+| Coupler name      | Description | Units |
+|-------------------|-------------|-------|
+| `air_density`       | air density at the surface 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 |
+| `height_sfc`        | height at the surface (only required when using `PartitionedStateFluxes`) | m |
+| `liquid_precipitation` | liquid precipitation at the surface | kg m^-2 s^-1 |
+| `radiative_energy_flux_sfc` | net radiative flux at the surface | W m^-2 |
+| `radiative_energy_flux_toa` | net radiative flux at the top of the atmosphere | W m^-2 |
+| `snow_precipitation` | snow precipitation at the surface | kg m^-2 s^-1 |
+| `turbulent_energy_flux` | aerodynamic turbulent surface fluxes of energy (sensible and latent heat) | W m^-2 |
+| `turbulent_moisture_flux` | aerodynamic turbulent surface fluxes of energy (evaporation) | kg m^-2 s^-1 |
+| `thermo_state_int`  | thermodynamic state at the first internal model level | |
+| `uv_int`            | horizontal wind velocity vector at the first internal model level | m s^-1 |
+| `water`             | globally integrated water | kg |
+
+
+
+- `update_field!(::AtmosModelSimulation. ::Val{property}, field)`:
+A function to update the value of property in the component model
+simulation, using the values in `field`. This update should
+be done in place. If this function isn't extended for a property,
+that property will remain constant throughout the simulation
+and a warning will be raised.
+This function is expected to be extended for the
+following properties:
+
+| Coupler name      | Description | Units |
+|-------------------|-------------|-------|
+| `co2`              | global mean co2 | ppm |
+| `surface_albedo`   | bulk surface albedo over the whole surface space | |
+| `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 |
+
+
+### SurfaceModelSimulation - required functions
+Analogously to the `AtmosModelSimulation`, a `SurfaceModelSimulation`
+requires additional functions to those required for a general `ComponentModelSimulation`.
+- `get_field(::SurfaceModelSimulation. ::Val{property})`: This getter
+function returns the value of the field property for the simulation at
+the current time. For a `SurfaceModelSimulation`, it must be extended
+for the following properties:
+
+| Coupler name      | Description | Units |
+|-------------------|-------------|-------|
+| `air_density`       | surface air density | kg m^-3 |
+| `area_fraction`     | fraction of the simulation grid surface area this model covers | |
+| `beta`              | factor that scales evaporation based on its estimated level of saturation | |
+| `roughness_buoyancy` | aerodynamic roughness length for buoyancy | m |
+| `roughness_momentum` | aerodynamic roughness length for momentum | m |
+| `surface_albedo`    | bulk surface albedo | |
+| `surface_humidity`  | surface humidity | kg kg^-1 |
+| `surface_temperature` | surface temperature | K |
+
+- `update_field!(::SurfaceModelSimulation. ::Val{property}, field)`:
+A function to update the value of property in the component model
+simulation, using the values in `field` passed from the coupler
+This update should be done in place. If this function
+isn't extended for a property,
+that property will remain constant throughout the simulation
+and a warning will be raised.
+This function is expected to be extended for the
+following properties:
+
+| Coupler name      | Description | Units |
+|-------------------|-------------|-------|
+| `air_density`       | surface air density | kg m^-3 |
+| `area_fraction`     | fraction of the simulation grid surface area this model covers | |
+| `liquid_precipitation` | liquid precipitation at the surface | kg m^-2 s^-1 |
+| `radiative_energy_flux_sfc` | net radiative flux at the surface | W m^-2 |
+| `snow_precipitation` | snow precipitation at the surface | kg m^-2 s^-1 |
+| `turbulent_energy_flux` | aerodynamic turbulent surface fluxes of energy (sensible and latent heat) | W m^-2 |
+| `turbulent_moisture_flux` | aerodynamic turbulent surface fluxes of energy (evaporation) | kg m^-2 s^-1 |
+
+### SurfaceModelSimulation - optional functions
+- `update_turbulent_fluxes_point!(::ComponentModelSimulation, fields::NamedTuple, colidx)`:
+This function updates the turbulent fluxes of the component model simulation
+at this point in horizontal space. The values are updated using the energy
+and moisture turbulent fluxes stored in fields which are calculated by the
+coupler. Note that this function is only required when using the
+`PartitionedStateFluxes` option of ClimaCoupler.jl. See our `FluxCalculator`
+module docs for more information.
+
+### Prescribed surface conditions - SurfaceStub
+- `SurfaceStub` is a `SurfaceModelSimulation`, but it only contains
+required data in `<surface_stub>.cache`, e.g., for the calculation
+of surface fluxes through a prescribed surface state. The above
+adapter functions are already predefined for `SurfaceStub`, with
+the cache variables specified as:
 ```
+get_field(sim::SurfaceStub, ::Val{:air_density}) = sim.cache.ρ_sfc
 get_field(sim::SurfaceStub, ::Val{:area_fraction}) = sim.cache.area_fraction
-get_field(sim::SurfaceStub, ::Val{:surface_temperature}) = sim.cache.T_sfc
-get_field(sim::SurfaceStub, ::Val{:albedo}) = sim.cache.α
-get_field(sim::SurfaceStub, ::Val{:roughness_momentum}) = sim.cache.z0m
-get_field(sim::SurfaceStub, ::Val{:roughness_buoyancy}) = sim.cache.z0b
 get_field(sim::SurfaceStub, ::Val{:beta}) = sim.cache.beta
+get_field(sim::SurfaceStub, ::Val{:energy}) = nothing
+get_field(sim::SurfaceStub, ::Val{:roughness_buoyancy}) = sim.cache.z0b
+get_field(sim::SurfaceStub, ::Val{:roughness_momentum}) = sim.cache.z0m
+get_field(sim::SurfaceStub, ::Val{:surface_albedo}) = sim.cache.α
+get_field(sim::SurfaceStub, ::Val{:surface_humidity}) = TD.q_vap_saturation_generic.(sim.cache.thermo_params, sim.cache.T_sfc, sim.cache.ρ_sfc, sim.cache.phase)
+get_field(sim::SurfaceStub, ::Val{:surface_temperature}) = sim.cache.T_sfc
+get_field(sim::SurfaceStub, ::Val{:water}) = nothing
 ```
-with the corresponding `update_field!` functions
+and with the corresponding `update_field!` functions
 ```
+function update_field!(sim::SurfaceStub, ::Val{:air_density}, field)
+    sim.cache.ρ_sfc .= field
+end
 function update_field!(sim::SurfaceStub, ::Val{:area_fraction}, field::Fields.Field)
     sim.cache.area_fraction .= field
 end
@@ -41,15 +184,3 @@ function update_field!(sim::SurfaceStub, ::Val{:surface_temperature}, field::Fie
     sim.cache.T_sfc .= field
 end
 ```
-
-## Interfacer API
-
-```@docs
-    ClimaCoupler.Interfacer.ComponentModelSimulation
-    ClimaCoupler.Interfacer.AtmosModelSimulation
-    ClimaCoupler.Interfacer.SurfaceModelSimulation
-    ClimaCoupler.Interfacer.SurfaceStub
-    ClimaCoupler.Interfacer.name
-    ClimaCoupler.Interfacer.get_field
-    ClimaCoupler.Interfacer.update_field!
-```
diff --git a/experiments/AMIP/components/atmosphere/climaatmos_init.jl b/experiments/AMIP/components/atmosphere/climaatmos_init.jl
index 336d4a8fee..a3b4e217e8 100644
--- a/experiments/AMIP/components/atmosphere/climaatmos_init.jl
+++ b/experiments/AMIP/components/atmosphere/climaatmos_init.jl
@@ -95,7 +95,7 @@ function atmos_init(::Type{FT}, atmos_config_dict::Dict) where {FT}
 end
 
 # extensions required by the Interfacer
-get_field(sim::ClimaAtmosSimulation, ::Val{:radiative_energy_flux}) =
+get_field(sim::ClimaAtmosSimulation, ::Val{:radiative_energy_flux_sfc}) =
     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
@@ -118,13 +118,10 @@ 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{:cv_m}) = TD.cv_m.(thermo_params, sim.integrator.p.precomputed.ᶜts)
-get_field(sim::ClimaAtmosSimulation, ::Val{:gas_constant_air}) =
-    TD.gas_constant_air.(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_gm}, field)
+function update_field!(atmos_sim::ClimaAtmosSimulation, ::Val{:co2}, field)
     if atmos_sim.integrator.p.atmos.radiation_mode isa CA.RRTMGPI.GrayRadiation
         @warn("Gray radiation model initialized, skipping CO2 update", maxlog = 1)
         return
@@ -137,7 +134,7 @@ function update_field!(sim::ClimaAtmosSimulation, ::Val{:surface_temperature}, c
     sim.integrator.p.radiation.radiation_model.surface_temperature .= CA.RRTMGPI.field2array(csf.T_S)
 end
 
-function update_field!(sim::ClimaAtmosSimulation, ::Val{:albedo}, field)
+function update_field!(sim::ClimaAtmosSimulation, ::Val{:surface_albedo}, field)
     sim.integrator.p.radiation.radiation_model.diffuse_sw_surface_albedo .=
         reshape(CA.RRTMGPI.field2array(field), 1, length(parent(field)))
     sim.integrator.p.radiation.radiation_model.direct_sw_surface_albedo .=
@@ -178,7 +175,7 @@ step!(sim::ClimaAtmosSimulation, t) = step!(sim.integrator, t - sim.integrator.t
 reinit!(sim::ClimaAtmosSimulation) = reinit!(sim.integrator)
 
 function update_sim!(atmos_sim::ClimaAtmosSimulation, csf, turbulent_fluxes)
-    update_field!(atmos_sim, Val(:albedo), csf.albedo)
+    update_field!(atmos_sim, Val(:surface_albedo), csf.surface_albedo)
     update_field!(atmos_sim, Val(:surface_temperature), csf)
 
     if turbulent_fluxes isa PartitionedStateFluxes
@@ -294,12 +291,12 @@ function get_model_state_vector(sim::ClimaAtmosSimulation)
 end
 
 """
-    get_field(atmos_sim::ClimaAtmosSimulation, ::Val{:F_radiative_TOA})
+    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{:F_radiative_TOA})
+function get_field(atmos_sim::ClimaAtmosSimulation, ::Val{:radiative_energy_flux_toa})
     radiation = atmos_sim.integrator.p.radiation.radiation_model
     FT = eltype(atmos_sim.integrator.u)
     # save radiation source
diff --git a/experiments/AMIP/components/land/bucket_utils.jl b/experiments/AMIP/components/land/bucket_utils.jl
index 58b94a5613..46be1dd87d 100644
--- a/experiments/AMIP/components/land/bucket_utils.jl
+++ b/experiments/AMIP/components/land/bucket_utils.jl
@@ -44,7 +44,7 @@ get_field(sim::BucketSimulation, ::Val{:roughness_momentum}) = sim.model.paramet
 get_field(sim::BucketSimulation, ::Val{:roughness_buoyancy}) = sim.model.parameters.z_0b
 get_field(sim::BucketSimulation, ::Val{:beta}) =
     ClimaLand.surface_evaporative_scaling(sim.model, sim.integrator.u, sim.integrator.p)
-get_field(sim::BucketSimulation, ::Val{:albedo}) =
+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
@@ -56,7 +56,7 @@ function update_field!(sim::BucketSimulation, ::Val{:turbulent_moisture_flux}, f
     ρ_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
-function update_field!(sim::BucketSimulation, ::Val{:radiative_energy_flux}, field)
+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{:liquid_precipitation}, field)
diff --git a/experiments/AMIP/components/ocean/eisenman_seaice_init.jl b/experiments/AMIP/components/ocean/eisenman_seaice_init.jl
index 8a1652e875..e9e46ecaad 100644
--- a/experiments/AMIP/components/ocean/eisenman_seaice_init.jl
+++ b/experiments/AMIP/components/ocean/eisenman_seaice_init.jl
@@ -1,6 +1,7 @@
 import ClimaCoupler: FluxCalculator
 import ClimaTimeSteppers as CTS
-import ClimaCoupler.FluxCalculator: update_turbulent_fluxes_point!, differentiate_turbulent_fluxes!
+import ClimaCoupler.FluxCalculator:
+    update_turbulent_fluxes_point!, differentiate_turbulent_fluxes!, surface_thermo_state
 import ClimaCoupler.Interfacer: get_field, update_field!
 using ClimaCore.Fields: getindex
 import SciMLBase: step!, reinit!, init, ODEProblem
@@ -247,7 +248,7 @@ 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{:albedo}) =
+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
@@ -259,7 +260,7 @@ 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}, field)
+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)
@@ -293,6 +294,31 @@ Extension of differentiate_turbulent_fluxes! from FluxCalculator to get the turb
 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
diff --git a/experiments/AMIP/components/ocean/prescr_seaice_init.jl b/experiments/AMIP/components/ocean/prescr_seaice_init.jl
index 926c92c1ea..805d733931 100644
--- a/experiments/AMIP/components/ocean/prescr_seaice_init.jl
+++ b/experiments/AMIP/components/ocean/prescr_seaice_init.jl
@@ -124,10 +124,10 @@ end
 
 # file-specific
 """
-    clean_sic(SIC, _info)
+    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.
 """
-clean_sic(SIC, _info) = swap_space!(zeros(axes(_info.land_fraction)), SIC) ./ float_type_bcf(_info)(100.0)
+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} =
@@ -139,7 +139,7 @@ get_field(sim::PrescribedIceSimulation, ::Val{:surface_humidity}) = sim.integrat
 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{:beta}) = convert(eltype(sim.integrator.u), 1.0)
-get_field(sim::PrescribedIceSimulation, ::Val{:albedo}) = sim.integrator.p.params.α
+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
 
@@ -150,7 +150,7 @@ 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}, field)
+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)
diff --git a/experiments/AMIP/components/ocean/slab_ocean_init.jl b/experiments/AMIP/components/ocean/slab_ocean_init.jl
index 1ab7dc66a7..c18e180f9c 100644
--- a/experiments/AMIP/components/ocean/slab_ocean_init.jl
+++ b/experiments/AMIP/components/ocean/slab_ocean_init.jl
@@ -126,10 +126,10 @@ end
 
 # file specific
 """
-    clean_sst(SST::FT, _info)
+    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)
 """
-clean_sst(SST, _info) = (swap_space!(zeros(axes(_info.land_fraction)), SST) .+ float_type_bcf(_info)(273.15))
+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
@@ -137,7 +137,7 @@ get_field(sim::SlabOceanSimulation, ::Val{:surface_humidity}) = sim.integrator.p
 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{:beta}) = convert(eltype(sim.integrator.u), 1.0)
-get_field(sim::SlabOceanSimulation, ::Val{:albedo}) = sim.integrator.p.params.α
+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
 
@@ -148,7 +148,7 @@ end
 function update_field!(sim::SlabOceanSimulation, ::Val{:turbulent_energy_flux}, field)
     parent(sim.integrator.p.F_turb_energy) .= parent(field)
 end
-function update_field!(sim::SlabOceanSimulation, ::Val{:radiative_energy_flux}, field)
+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)
diff --git a/experiments/AMIP/coupler_driver.jl b/experiments/AMIP/coupler_driver.jl
index 5b38393939..f7047a36d0 100644
--- a/experiments/AMIP/coupler_driver.jl
+++ b/experiments/AMIP/coupler_driver.jl
@@ -276,7 +276,7 @@ if mode_name == "amip"
         boundary_space,
         comms_ctx,
         interpolate_daily = true,
-        scaling_function = clean_sst, ## convert to Kelvin
+        scaling_function = scale_sst, ## convert to Kelvin
         land_fraction = land_fraction,
         date0 = date0,
         mono = mono_surface,
@@ -305,7 +305,7 @@ if mode_name == "amip"
         boundary_space,
         comms_ctx,
         interpolate_daily = true,
-        scaling_function = clean_sic, ## convert to fraction
+        scaling_function = scale_sic, ## convert to fraction
         land_fraction = land_fraction,
         date0 = date0,
         mono = mono_surface,
@@ -339,7 +339,7 @@ if mode_name == "amip"
 
     update_midmonth_data!(date0, CO2_info)
     CO2_init = interpolate_midmonth_to_daily(date0, CO2_info)
-    update_field!(atmos_sim, Val(:co2_gm), CO2_init)
+    update_field!(atmos_sim, Val(:co2), CO2_init)
 
     mode_specifics = (; name = mode_name, SST_info = SST_info, SIC_info = SIC_info, CO2_info = CO2_info)
 
@@ -447,7 +447,7 @@ coupler_field_names = (
     :z0b_S,
     :ρ_sfc,
     :q_sfc,
-    :albedo,
+    :surface_albedo,
     :beta,
     :F_turb_energy,
     :F_turb_moisture,
@@ -456,7 +456,7 @@ coupler_field_names = (
     :F_radiative,
     :P_liq,
     :P_snow,
-    :F_radiative_TOA,
+    :radiative_energy_flux_toa,
     :P_net,
 )
 coupler_fields =
@@ -569,7 +569,7 @@ step!(ice_sim, Δt_cpl)
 # 3) coupler re-imports updated surface fields and calculates turbulent fluxes, while updating atmos sfc_conditions
 if turbulent_fluxes isa CombinedStateFluxes
     # calculate fluxes using combined surface states on the atmos grid
-    import_combined_surface_fields!(cs.fields, cs.model_sims, cs.boundary_space, turbulent_fluxes) # i.e. T_sfc, albedo, z0, beta, q_sfc
+    import_combined_surface_fields!(cs.fields, cs.model_sims, cs.boundary_space, turbulent_fluxes) # i.e. T_sfc, surface_albedo, z0, beta, q_sfc
     combined_turbulent_fluxes!(cs.model_sims, cs.fields, turbulent_fluxes) # this updates the atmos thermo state, sfc_ts
 elseif turbulent_fluxes isa PartitionedStateFluxes
     # calculate turbulent fluxes in surface models and save the weighted average in coupler fields
@@ -631,7 +631,7 @@ function solve_coupler!(cs)
                 update_midmonth_data!(cs.dates.date[1], cs.mode.CO2_info)
             end
             CO2_current = interpolate_midmonth_to_daily(cs.dates.date[1], cs.mode.CO2_info)
-            update_field!(atmos_sim, Val(:co2_gm), CO2_current)
+            update_field!(atmos_sim, Val(:co2), CO2_current)
 
             ## calculate and accumulate diagnostics at each timestep
             ClimaComms.barrier(comms_ctx)
@@ -654,7 +654,7 @@ function solve_coupler!(cs)
         step_model_sims!(cs.model_sims, t)
 
         ## exchange combined fields and (if specified) calculate fluxes using combined states
-        import_combined_surface_fields!(cs.fields, cs.model_sims, cs.boundary_space, turbulent_fluxes) # i.e. T_sfc, albedo, z0, beta
+        import_combined_surface_fields!(cs.fields, cs.model_sims, cs.boundary_space, turbulent_fluxes) # i.e. T_sfc, surface_albedo, z0, beta
         if turbulent_fluxes isa CombinedStateFluxes
             combined_turbulent_fluxes!(cs.model_sims, cs.fields, turbulent_fluxes) # this updates the surface thermo state, sfc_ts, in ClimaAtmos (but also unnecessarily calculates fluxes)
         elseif turbulent_fluxes isa PartitionedStateFluxes
diff --git a/experiments/AMIP/user_io/debug_plots.jl b/experiments/AMIP/user_io/debug_plots.jl
index 243c35c78a..85f49b9b29 100644
--- a/experiments/AMIP/user_io/debug_plots.jl
+++ b/experiments/AMIP/user_io/debug_plots.jl
@@ -27,8 +27,19 @@ If `cs_fields_ref` is provided (e.g., using a copy of cs.fields from the initial
 plot the anomalies of the fields with respect to `cs_fields_ref`.
 """
 function debug(cs_fields::NamedTuple, dir, cs_fields_ref = nothing)
-    field_names =
-        (:albedo, :F_radiative, :F_turb_energy, :F_turb_moisture, :P_liq, :T_S, :ρ_sfc, :q_sfc, :beta, :z0b_S, :z0m_S)
+    field_names = (
+        :surface_albedo,
+        :F_radiative,
+        :F_turb_energy,
+        :F_turb_moisture,
+        :P_liq,
+        :T_S,
+        :ρ_sfc,
+        :q_sfc,
+        :beta,
+        :z0b_S,
+        :z0m_S,
+    )
     all_plots = []
     for field_name in field_names
         field = getproperty(cs_fields, field_name)
diff --git a/src/ConservationChecker.jl b/src/ConservationChecker.jl
index b6a2845775..5c2fdaab39 100644
--- a/src/ConservationChecker.jl
+++ b/src/ConservationChecker.jl
@@ -97,13 +97,15 @@ function check_conservation!(
         sim_name = Symbol(Interfacer.name(sim))
         if sim isa Interfacer.AtmosModelSimulation
             # save radiation source
-            parent(coupler_sim.fields.F_radiative_TOA) .= parent(Interfacer.get_field(sim, Val(:F_radiative_TOA)))
+            parent(coupler_sim.fields.radiative_energy_flux_toa) .=
+                parent(Interfacer.get_field(sim, Val(:radiative_energy_flux_toa)))
 
             if isempty(ccs.toa_net_source)
-                radiation_sources_accum = sum(coupler_sim.fields.F_radiative_TOA .* FT(coupler_sim.Δt_cpl)) # ∫ J / m^2 dA
+                radiation_sources_accum = sum(coupler_sim.fields.radiative_energy_flux_toa .* FT(coupler_sim.Δt_cpl)) # ∫ J / m^2 dA
             else
                 radiation_sources_accum =
-                    sum(coupler_sim.fields.F_radiative_TOA .* FT(coupler_sim.Δt_cpl)) .+ ccs.toa_net_source[end] # ∫ J / m^2 dA
+                    sum(coupler_sim.fields.radiative_energy_flux_toa .* FT(coupler_sim.Δt_cpl)) .+
+                    ccs.toa_net_source[end] # ∫ J / m^2 dA
             end
             push!(ccs.toa_net_source, radiation_sources_accum)
 
diff --git a/src/FieldExchanger.jl b/src/FieldExchanger.jl
index a40b39ae92..8101d6e80c 100644
--- a/src/FieldExchanger.jl
+++ b/src/FieldExchanger.jl
@@ -16,7 +16,7 @@ using ClimaCoupler: Interfacer, FluxCalculator, Regridder, Utilities
     import_atmos_fields!(csf, model_sims, boundary_space, turbulent_fluxes)
 
 Updates the coupler with the atmospheric fluxes. The `Interfacer.get_field` functions
-(`:turbulent_energy_flux`, `:turbulent_moisture_flux`, `:radiative_energy_flux`, `:liquid_precipitation`, `:snow_precipitation`)
+(`:turbulent_energy_flux`, `:turbulent_moisture_flux`, `:radiative_energy_flux_sfc`, `:liquid_precipitation`, `:snow_precipitation`)
 have to be defined for the amtospheric component model type.
 
 # Arguments
@@ -38,7 +38,7 @@ function import_atmos_fields!(csf, model_sims, boundary_space, turbulent_fluxes)
     Regridder.dummmy_remap!(csf.ρ_sfc, FluxCalculator.calculate_surface_air_density(atmos_sim, csf.T_S))
 
     # radiative fluxes
-    Regridder.dummmy_remap!(csf.F_radiative, Interfacer.get_field(atmos_sim, Val(:radiative_energy_flux)))
+    Regridder.dummmy_remap!(csf.F_radiative, Interfacer.get_field(atmos_sim, Val(:radiative_energy_flux_sfc)))
 
     # precipitation
     Regridder.dummmy_remap!(csf.P_liq, Interfacer.get_field(atmos_sim, Val(:liquid_precipitation)))
@@ -49,7 +49,7 @@ end
     import_combined_surface_fields!(csf, model_sims, boundary_space, turbulent_fluxes)
 
 Updates the coupler with the surface properties. The `Interfacer.get_field` functions for
-(`:surface_temperature`, `:albedo`, `:roughness_momentum`, `:roughness_buoyancy`, `:beta`)
+(`:surface_temperature`, `:surface_albedo`, `:roughness_momentum`, `:roughness_buoyancy`, `:beta`)
 need to be specified for each surface model.
 
 # Arguments
@@ -67,8 +67,8 @@ function import_combined_surface_fields!(csf, model_sims, boundary_space, turbul
     Regridder.combine_surfaces!(combined_field, model_sims, Val(:surface_temperature))
     Regridder.dummmy_remap!(csf.T_S, combined_field)
 
-    Regridder.combine_surfaces!(combined_field, model_sims, Val(:albedo))
-    Regridder.dummmy_remap!(csf.albedo, combined_field)
+    Regridder.combine_surfaces!(combined_field, model_sims, Val(:surface_albedo))
+    Regridder.dummmy_remap!(csf.surface_albedo, combined_field)
 
     if turbulent_fluxes isa FluxCalculator.CombinedStateFluxes
         Regridder.combine_surfaces!(combined_field, model_sims, Val(:roughness_momentum))
@@ -97,7 +97,7 @@ Updates the surface fields for temperature, roughness length, albedo, and specif
 """
 function update_sim!(atmos_sim::Interfacer.AtmosModelSimulation, csf, turbulent_fluxes)
 
-    Interfacer.update_field!(atmos_sim, Val(:albedo), csf.albedo)
+    Interfacer.update_field!(atmos_sim, Val(:surface_albedo), csf.surface_albedo)
     Interfacer.update_field!(atmos_sim, Val(:surface_temperature), csf.T_S)
 
     if turbulent_fluxes isa FluxCalculator.CombinedStateFluxes
@@ -137,7 +137,7 @@ function update_sim!(sim::Interfacer.SurfaceModelSimulation, csf, turbulent_flux
     end
 
     # radiative fluxes
-    Interfacer.update_field!(sim, Val(:radiative_energy_flux), FT.(mask .* csf.F_radiative))
+    Interfacer.update_field!(sim, Val(:radiative_energy_flux_sfc), FT.(mask .* csf.F_radiative))
 
     # precipitation
     Interfacer.update_field!(sim, Val(:liquid_precipitation), csf.P_liq)
diff --git a/src/FluxCalculator.jl b/src/FluxCalculator.jl
index 461c7ff2d1..c45e6b9217 100644
--- a/src/FluxCalculator.jl
+++ b/src/FluxCalculator.jl
@@ -379,38 +379,12 @@ end
 
 update_turbulent_fluxes_point!(sim::Interfacer.SurfaceStub, fields::NamedTuple, colidx::Fields.ColumnIndex) = nothing
 
-differentiate_turbulent_fluxes!(::Interfacer.SurfaceModelSimulation, args) = nothing
-
 """
-    differentiate_turbulent_fluxes(sim::Interfacer.SurfaceModelSimulation, thermo_params, input_args, fluxes, δT_sfc = 0.1)
+    differentiate_turbulent_fluxes!(sim::Interfacer.SurfaceModelSimulation, args)
 
-Differentiates the turbulent fluxes in the surface model simulation `sim` with respect to the surface temperature,
-using δT_sfc as the perturbation.
+This function provides a placeholder for differentiating fluxes with respect to
+surface temperature in surface energy balance calculations.
 """
-function differentiate_turbulent_fluxes!(
-    sim::Interfacer.SurfaceModelSimulation,
-    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
+differentiate_turbulent_fluxes!(::Interfacer.SurfaceModelSimulation, args) = nothing
 
 end # module
diff --git a/src/Interfacer.jl b/src/Interfacer.jl
index 343113dbd1..cb854ff06a 100644
--- a/src/Interfacer.jl
+++ b/src/Interfacer.jl
@@ -101,6 +101,56 @@ abstract type SeaIceModelSimulation <: SurfaceModelSimulation end
 abstract type LandModelSimulation <: SurfaceModelSimulation end
 abstract type OceanModelSimulation <: SurfaceModelSimulation end
 
+"""
+    get_field(sim::AtmosModelSimulation, val::Val)
+
+A getter function that should not allocate. Here we implement a default that
+will raise an error if `get_field` isn't defined for all required fields of
+an atmosphere component model.
+"""
+get_field(
+    sim::AtmosModelSimulation,
+    val::Union{
+        Val{:air_density},
+        Val{:air_temperature},
+        Val{:energy},
+        Val{:radiative_energy_flux_toa},
+        Val{:height_int},
+        Val{:height_sfc},
+        Val{:liquid_precipitation},
+        Val{:radiative_energy_flux_sfc},
+        Val{:snow_precipitation},
+        Val{:turblent_energy_flux},
+        Val{:turbulent_moisture_flux},
+        Val{:thermo_state_int},
+        Val{:uv_int},
+        Val{:water},
+    },
+) = get_field_error(sim, val)
+
+"""
+    get_field(sim::SurfaceModelSimulation, val::Val)
+
+A getter function that should not allocate. Here we implement a default that
+will raise an error if `get_field` isn't defined for all required fields of
+a surface component model.
+"""
+get_field(
+    sim::SurfaceModelSimulation,
+    val::Union{
+        Val{:air_density},
+        Val{:area_fraction},
+        Val{:beta},
+        Val{:roughness_buoyancy},
+        Val{:roughness_momentum},
+        Val{:surface_albedo},
+        Val{:surface_humidity},
+        Val{:surface_temperature},
+    },
+) = get_field_error(sim, val)
+
+get_field_error(sim, val::Val{X}) where {X} = error("undefined field `$X` for " * name(sim))
+
 
 """
     SurfaceStub
@@ -111,25 +161,28 @@ struct SurfaceStub{I} <: SurfaceModelSimulation
     cache::I
 end
 
+"""
+    stub_init(cache)
+
+Initialization function for SurfaceStub simulation type.
+"""
+stub_init(cache) = SurfaceStub(cache)
+
 """
     get_field(::SurfaceStub, ::Val)
 
 A getter function, that should not allocate. If undefined, it returns a descriptive error.
 """
+get_field(sim::SurfaceStub, ::Val{:air_density}) = sim.cache.ρ_sfc
 get_field(sim::SurfaceStub, ::Val{:area_fraction}) = sim.cache.area_fraction
-get_field(sim::SurfaceStub, ::Val{:surface_temperature}) = sim.cache.T_sfc
-get_field(sim::SurfaceStub, ::Val{:albedo}) = sim.cache.α
-get_field(sim::SurfaceStub, ::Val{:roughness_momentum}) = sim.cache.z0m
-get_field(sim::SurfaceStub, ::Val{:roughness_buoyancy}) = sim.cache.z0b
 get_field(sim::SurfaceStub, ::Val{:beta}) = sim.cache.beta
+get_field(sim::SurfaceStub, ::Val{:energy}) = nothing
+get_field(sim::SurfaceStub, ::Val{:roughness_buoyancy}) = sim.cache.z0b
+get_field(sim::SurfaceStub, ::Val{:roughness_momentum}) = sim.cache.z0m
+get_field(sim::SurfaceStub, ::Val{:surface_albedo}) = sim.cache.α
 get_field(sim::SurfaceStub, ::Val{:surface_humidity}) =
     TD.q_vap_saturation_generic.(sim.cache.thermo_params, sim.cache.T_sfc, sim.cache.ρ_sfc, sim.cache.phase)
-
-function get_field(sim::ComponentModelSimulation, val::Val)
-    error("undefined field $val for " * name(sim))
-end
-
-get_field(sim::SurfaceStub, ::Val{:energy}) = nothing
+get_field(sim::SurfaceStub, ::Val{:surface_temperature}) = sim.cache.T_sfc
 get_field(sim::SurfaceStub, ::Val{:water}) = nothing
 
 """
@@ -146,21 +199,41 @@ function get_field(sim::ComponentModelSimulation, val::Val, colidx::Fields.Colum
 end
 
 """
-    update_field!(::ComponentModelSimulation, ::Val, _...)
+    update_field!(::AtmosModelSimulation, ::Val, _...)
 
-No update in unspecified in the particular component model simulation.
+Default functions for updating fields at each timestep in an atmosphere
+component model simulation. This should be extended by component models.
+If it isn't extended, the field won't be updated and a warning will be raised.
 """
-update_field!(sim::ComponentModelSimulation, val::Val, _...) = nothing
+update_field!(
+    sim::AtmosModelSimulation,
+    val::Union{Val{:co2}, Val{:surface_albedo}, Val{:surface_temperature}, Val{:turbulent_fluxes}},
+    _,
+) = update_field_warning(sim, val)
 
-# TODO:
-# function update_field!(sim::ComponentModelSimulation, val::Val, _...)
-#     warning = Warning("undefined `update!` for $val in " * name(sim) * ": skipping")
-#     @warn(warning.message, maxlog=10)
-#     return warning
-# end
-# struct Warning
-#     message::String
-# end
+"""
+    update_field!(::SurfaceModelSimulation, ::Val, _...)
+
+Default functions for updating fields at each timestep in an atmosphere
+component model simulation. This should be extended by component models.
+If it isn't extended, the field won't be updated and a warning will be raised.
+"""
+update_field!(
+    sim::SurfaceModelSimulation,
+    val::Union{
+        Val{:air_density},
+        Val{:area_fraction},
+        Val{:liquid_precipitation},
+        Val{:radiative_energy_flux_sfc},
+        Val{:snow_precipitation},
+        Val{:turbulent_energy_flux},
+        Val{:turbulent_moisture_flux},
+    },
+    _,
+) = update_field_warning(sim, val)
+
+update_field_warning(sim, val::Val{X}) where {X} =
+    @warn("`update_field!` is not extended for the `$X` field of " * name(sim) * ": skipping update.", maxlog = 1)
 
 """
     update_field!(sim::SurfaceStub, ::Val{:area_fraction}, field::Fields.Field)
diff --git a/test/conservation_checker_tests.jl b/test/conservation_checker_tests.jl
index 77e7285fd6..327392cf58 100644
--- a/test/conservation_checker_tests.jl
+++ b/test/conservation_checker_tests.jl
@@ -25,7 +25,7 @@ struct TestAtmos{I} <: Interfacer.AtmosModelSimulation
     i::I
 end
 name(s::TestAtmos) = "TestAtmos"
-get_field(s::TestAtmos, ::Val{:F_radiative_TOA}) = ones(s.i.space) .* 200
+get_field(s::TestAtmos, ::Val{:radiative_energy_flux_toa}) = ones(s.i.space) .* 200
 get_field(s::TestAtmos, ::Val{:water}) = ones(s.i.space) .* 1
 function get_field(s::TestAtmos, ::Val{:energy})
     FT = Domains.float_type(Meshes.domain(s.i.space.grid.topology.mesh))
@@ -74,13 +74,13 @@ for FT in (Float32, Float64)
 
         # coupler fields
         cf = (;
-            F_radiative_TOA = Fields.ones(space),
+            radiative_energy_flux_toa = Fields.ones(space),
             P_net = Fields.zeros(space),
             P_liq = Fields.zeros(space),
             P_snow = Fields.zeros(space),
             F_turb_moisture = Fields.zeros(space),
         )
-        @. cf.F_radiative_TOA = 200
+        @. cf.radiative_energy_flux_toa = 200
         @. cf.P_liq = -100
 
         # init
@@ -103,7 +103,7 @@ for FT in (Float32, Float64)
         )
 
         # set non-zero radiation and precipitation
-        F_r = cf.F_radiative_TOA
+        F_r = cf.radiative_energy_flux_toa
         P = cf.P_liq
         Δt = cs.Δt_cpl
 
@@ -152,13 +152,13 @@ for FT in (Float32, Float64)
 
         # coupler fields
         cf = (;
-            F_radiative_TOA = Fields.ones(space),
+            radiative_energy_flux_toa = Fields.ones(space),
             P_net = Fields.zeros(space),
             P_liq = Fields.zeros(space),
             P_snow = Fields.zeros(space),
             F_turb_moisture = Fields.zeros(space),
         )
-        @. cf.F_radiative_TOA = 200
+        @. cf.radiative_energy_flux_toa = 200
         @. cf.P_liq = -100
 
         # init
diff --git a/test/debug/debug_amip_plots.jl b/test/debug/debug_amip_plots.jl
index b15c6c0d17..fdda9cb14a 100644
--- a/test/debug/debug_amip_plots.jl
+++ b/test/debug/debug_amip_plots.jl
@@ -38,8 +38,19 @@ plot_field_names(sim::SurfaceStub) = (:stub_field,)
 @testset "import_atmos_fields!" begin
 
     boundary_space = TestHelper.create_space(FT)
-    coupler_names =
-        (:albedo, :F_radiative, :F_turb_energy, :F_turb_moisture, :P_liq, :T_S, :ρ_sfc, :q_sfc, :beta, :z0b_S, :z0m_S)
+    coupler_names = (
+        :surface_albedo,
+        :F_radiative,
+        :F_turb_energy,
+        :F_turb_moisture,
+        :P_liq,
+        :T_S,
+        :ρ_sfc,
+        :q_sfc,
+        :beta,
+        :z0b_S,
+        :z0m_S,
+    )
     atmos_names = (:atmos_field,)
     surface_names = (:surface_field,)
     stub_names = (:stub_field,)
diff --git a/test/field_exchanger_tests.jl b/test/field_exchanger_tests.jl
index 5eb57521c6..5ebba43db5 100644
--- a/test/field_exchanger_tests.jl
+++ b/test/field_exchanger_tests.jl
@@ -21,7 +21,7 @@ end
 
 get_field(sim::DummySimulation, ::Val{:turbulent_energy_flux}) = sim.cache.turbulent_energy_flux
 get_field(sim::DummySimulation, ::Val{:turbulent_moisture_flux}) = sim.cache.turbulent_moisture_flux
-get_field(sim::DummySimulation, ::Val{:radiative_energy_flux}) = sim.cache.radiative_energy_flux
+get_field(sim::DummySimulation, ::Val{:radiative_energy_flux_sfc}) = sim.cache.radiative_energy_flux_sfc
 get_field(sim::DummySimulation, ::Val{:liquid_precipitation}) = sim.cache.liquid_precipitation
 get_field(sim::DummySimulation, ::Val{:snow_precipitation}) = sim.cache.snow_precipitation
 
@@ -41,7 +41,7 @@ end
 
 get_field(sim::Union{TestSurfaceSimulation1, TestSurfaceSimulation2}, ::Val{:surface_temperature}) =
     sim.cache_field .* eltype(sim.cache_field)(1.0)
-get_field(sim::Union{TestSurfaceSimulation1, TestSurfaceSimulation2}, ::Val{:albedo}) =
+get_field(sim::Union{TestSurfaceSimulation1, TestSurfaceSimulation2}, ::Val{:surface_albedo}) =
     sim.cache_field .* eltype(sim.cache_field)(1.0)
 get_field(sim::Union{TestSurfaceSimulation1, TestSurfaceSimulation2}, ::Val{:roughness_momentum}) =
     sim.cache_field .* eltype(sim.cache_field)(1.0)
@@ -63,13 +63,16 @@ get_field(sim::Union{TestSurfaceSimulation2, TestSurfaceSimulation2}, ::Val{:sur
 struct TestAtmosSimulation{C} <: AtmosModelSimulation
     cache::C
 end
-function update_field!(sim::TestAtmosSimulation, ::Val{:albedo}, field)
-    parent(sim.cache.albedo) .= parent(field)
+function update_field!(sim::TestAtmosSimulation, ::Val{:surface_albedo}, field)
+    parent(sim.cache.surface_albedo) .= parent(field)
 end
 function update_field!(sim::TestAtmosSimulation, ::Val{:roughness_momentum}, field)
     parent(sim.cache.roughness_momentum) .= parent(field)
 end
 
+update_field!(sim::TestAtmosSimulation, ::Val{:roughness_buoyancy}, field) = nothing
+update_field!(sim::TestAtmosSimulation, ::Val{:beta}, field) = nothing
+
 #surface sim
 struct TestSurfaceSimulationLand{C} <: SurfaceModelSimulation
     cache::C
@@ -92,7 +95,7 @@ for FT in (Float32, Float64)
         atmos_names = (
             :turbulent_energy_flux,
             :turbulent_moisture_flux,
-            :radiative_energy_flux,
+            :radiative_energy_flux_sfc,
             :liquid_precipitation,
             :snow_precipitation,
         )
@@ -116,10 +119,10 @@ for FT in (Float32, Float64)
     @testset "import_combined_surface_fields! for FT=$FT" begin
         # coupler cache setup
         boundary_space = TestHelper.create_space(FT)
-        coupler_names = (:T_S, :z0m_S, :z0b_S, :albedo, :beta, :q_sfc)
+        coupler_names = (:T_S, :z0m_S, :z0b_S, :surface_albedo, :beta, :q_sfc)
 
         # coupler cache setup
-        exchanged_fields = (:surface_temperature, :albedo, :roughness_momentum, :roughness_buoyancy, :beta)
+        exchanged_fields = (:surface_temperature, :surface_albedo, :roughness_momentum, :roughness_buoyancy, :beta)
 
         sims = (; a = TestSurfaceSimulation1(ones(boundary_space)), b = TestSurfaceSimulation2(ones(boundary_space)))
 
@@ -131,7 +134,7 @@ for FT in (Float32, Float64)
             coupler_fields = NamedTuple{coupler_names}(ntuple(i -> Fields.zeros(boundary_space), length(coupler_names)))
             import_combined_surface_fields!(coupler_fields, sims, boundary_space, t)
             @test parent(coupler_fields.T_S)[1] == results[1]
-            @test parent(coupler_fields.albedo)[1] == results[1]
+            @test parent(coupler_fields.surface_albedo)[1] == results[1]
             @test parent(coupler_fields.z0m_S)[1] == results[i]
             @test parent(coupler_fields.z0b_S)[1] == results[i]
             @test parent(coupler_fields.beta)[1] == results[i]
@@ -146,7 +149,7 @@ for FT in (Float32, Float64)
             :T_S,
             :z0m_S,
             :z0b_S,
-            :albedo,
+            :surface_albedo,
             :beta,
             :F_turb_energy,
             :F_turb_moisture,
@@ -159,13 +162,13 @@ for FT in (Float32, Float64)
 
         # model cache setup
 
-        atmos_names = (:surface_temperature, :albedo, :roughness_momentum, :roughness_buoyancy, :beta)
+        atmos_names = (:surface_temperature, :surface_albedo, :roughness_momentum, :roughness_buoyancy, :beta)
         atmos_fields = NamedTuple{atmos_names}(ntuple(i -> Fields.zeros(boundary_space), length(atmos_names)))
 
         land_names = (
             :turbulent_energy_flux,
             :turbulent_moisture_flux,
-            :radiative_energy_flux,
+            :radiative_energy_flux_sfc,
             :liquid_precipitation,
             :snow_precipitation,
             :ρ_sfc,
@@ -189,7 +192,7 @@ for FT in (Float32, Float64)
             update_model_sims!(model_sims, coupler_fields, t)
 
             # test atmos
-            @test parent(model_sims.atmos_sim.cache.albedo)[1] == results[2]
+            @test parent(model_sims.atmos_sim.cache.surface_albedo)[1] == results[2]
             if t isa CombinedStateFluxes
                 @test parent(model_sims.atmos_sim.cache.roughness_momentum)[1] == results[2]
             else
@@ -206,7 +209,7 @@ for FT in (Float32, Float64)
             @test parent(model_sims.land_sim.cache.turbulent_moisture_flux)[1] == results[2]
 
             # unspecified variables
-            @test parent(model_sims.land_sim.cache.radiative_energy_flux)[1] == results[1]
+            @test parent(model_sims.land_sim.cache.radiative_energy_flux_sfc)[1] == results[1]
             @test parent(model_sims.land_sim.cache.liquid_precipitation)[1] == results[1]
             @test parent(model_sims.land_sim.cache.snow_precipitation)[1] == results[1]
 
diff --git a/test/flux_calculator_tests.jl b/test/flux_calculator_tests.jl
index f95a7ae1a1..417ee6c745 100644
--- a/test/flux_calculator_tests.jl
+++ b/test/flux_calculator_tests.jl
@@ -17,8 +17,7 @@ import ClimaCoupler.FluxCalculator:
     surface_inputs,
     get_surface_fluxes_point!,
     get_scheme_properties,
-    surface_thermo_state,
-    differentiate_turbulent_fluxes!
+    surface_thermo_state
 import ClimaCoupler: Interfacer
 
 import CLIMAParameters as CP
@@ -97,7 +96,7 @@ Interfacer.get_field(sim::TestOcean, ::Val{:beta}) = sim.integrator.p.beta
 Interfacer.get_field(sim::TestOcean, ::Val{:area_fraction}) = sim.integrator.p.area_fraction
 Interfacer.get_field(sim::TestOcean, ::Val{:heat_transfer_coefficient}) = sim.integrator.p.Ch
 Interfacer.get_field(sim::TestOcean, ::Val{:drag_coefficient}) = sim.integrator.p.Cd
-Interfacer.get_field(sim::TestOcean, ::Val{:albedo}) = sim.integrator.p.α
+Interfacer.get_field(sim::TestOcean, ::Val{:surface_albedo}) = sim.integrator.p.α
 
 function surface_thermo_state(
     sim::TestOcean,
@@ -130,9 +129,6 @@ Interfacer.get_field(sim::DummySurfaceSimulation3, ::Val{:heat_transfer_coeffici
 Interfacer.get_field(sim::DummySurfaceSimulation3, ::Val{:drag_coefficient}) = sim.integrator.p.Cd
 Interfacer.get_field(sim::DummySurfaceSimulation3, ::Val{:beta}) = sim.integrator.p.beta
 
-function Interfacer.update_field!(sim::DummySurfaceSimulation3, ::Val{:∂F_turb_energy∂T_sfc}, field, colidx)
-    sim.integrator.p.∂F_turb_energy∂T_sfc[colidx] .= field
-end
 function surface_thermo_state(
     sim::DummySurfaceSimulation3,
     thermo_params::ThermodynamicsParameters,
@@ -224,7 +220,7 @@ for FT in (Float32, Float64)
 
             coupler_cache_names = (
                 :T_S,
-                :albedo,
+                :surface_albedo,
                 :F_R_sfc,
                 :F_R_toa,
                 :P_liq,
@@ -326,63 +322,4 @@ for FT in (Float32, Float64)
         @test surface_thermo_state(surface_sim, thermo_params, thermo_state_int[colidx], colidx).ρ ==
               thermo_state_int[colidx].ρ
     end
-
-    @testset "differentiate_turbulent_fluxes! for FT=$FT" begin
-        boundary_space = TestHelper.create_space(FT)
-        _ones = Fields.ones(boundary_space)
-        surface_sim = DummySurfaceSimulation3(
-            [],
-            [],
-            [],
-            (;
-                T = _ones .* FT(300),
-                ρ = _ones .* FT(1.2),
-                p = (;
-                    q = _ones .* FT(0.00),
-                    area_fraction = _ones,
-                    Ch = FT(0.001),
-                    Cd = FT(0.001),
-                    beta = _ones,
-                    ∂F_turb_energy∂T_sfc = _ones .* 0,
-                    q_sfc = _ones .* 0,
-                ),
-            ),
-        )
-        atmos_sim =
-            TestAtmos((; FT = FT), [], [], (; T = _ones .* FT(300), ρ = _ones .* FT(1.2), q = _ones .* FT(0.00)))
-        thermo_params = get_thermo_params(atmos_sim)
-        colidx = Fields.ColumnIndex{2}((1, 1), 73) # arbitrary index
-
-        thermo_state_int = Interfacer.get_field(atmos_sim, Val(:thermo_state_int))[colidx]
-        surface_scheme = BulkScheme()
-        surface_params = get_surface_params(atmos_sim)
-        uₕ_int = Geometry.UVVector.(Geometry.Covariant12Vector.(_ones .* FT(1), _ones .* FT(1)))[colidx]
-        z_int = _ones[colidx]
-        z_sfc = (_ones .* FT(0))[colidx]
-        thermo_state_sfc = surface_thermo_state(surface_sim, thermo_params, thermo_state_int[colidx], colidx)
-        scheme_properties = get_scheme_properties(surface_scheme, surface_sim, colidx)
-        input_args = (;
-            thermo_state_sfc,
-            thermo_state_int,
-            uₕ_int,
-            z_int,
-            z_sfc,
-            surface_params,
-            surface_scheme,
-            scheme_properties,
-            colidx,
-        )
-
-        inputs = surface_inputs(surface_scheme, input_args)
-        fluxes = get_surface_fluxes_point!(inputs, surface_params)
-
-        dFdTs = differentiate_turbulent_fluxes!(surface_sim, thermo_params, input_args, fluxes, δT_sfc = 1)
-
-        sf_out = SF.surface_conditions.(surface_params, inputs)
-
-        cp_m = TD.cp_m.(thermo_params, thermo_state_int)
-        dFdTs_analytical = @. thermo_state_sfc.ρ * sf_out.Ch * SF.windspeed.(inputs) * cp_m
-
-        @test all(isapprox(parent(dFdTs), parent(dFdTs_analytical), atol = 0.1))
-    end
 end
diff --git a/test/interfacer_tests.jl b/test/interfacer_tests.jl
index 9ca7705443..0be4ec2f0c 100644
--- a/test/interfacer_tests.jl
+++ b/test/interfacer_tests.jl
@@ -26,6 +26,7 @@ end
 struct DummySimulation3{S} <: LandModelSimulation
     space::S
 end
+name(::DummySimulation3) = "DummySimulation3"
 
 get_field(sim::SurfaceModelSimulation, ::Val{:var}) = ones(sim.space)
 get_field(sim::SurfaceModelSimulation, ::Val{:surface_temperature}) = ones(sim.space) .* 300
@@ -92,7 +93,7 @@ for FT in (Float32, Float64)
         ))
         @test get_field(stub, Val(:area_fraction)) == FT(1)
         @test get_field(stub, Val(:surface_temperature)) == FT(280)
-        @test get_field(stub, Val(:albedo)) == 3
+        @test get_field(stub, Val(:surface_albedo)) == 3
         @test get_field(stub, Val(:roughness_momentum)) == 4
         @test get_field(stub, Val(:roughness_buoyancy)) == 5
         @test get_field(stub, Val(:beta)) == 6
@@ -116,3 +117,61 @@ for FT in (Float32, Float64)
         @test parent(get_field(stub, Val(:surface_temperature)))[1] == FT(2)
     end
 end
+
+@testset "update_field! warnings" begin
+    FT = Float64
+    space = TestHelper.create_space(FT)
+    dummy_field = Fields.ones(space)
+    sim = DummySimulation3(space)
+
+    # Test that update_field! gives correct warnings for unextended fields
+    @test_logs (
+        :warn,
+        "`update_field!` is not extended for the `air_density` field of " * name(sim) * ": skipping update.",
+    ) update_field!(sim, Val(:air_density), dummy_field)
+    @test_logs (
+        :warn,
+        "`update_field!` is not extended for the `area_fraction` field of " * name(sim) * ": skipping update.",
+    ) update_field!(sim, Val(:area_fraction), dummy_field)
+    @test_logs (:warn, "`update_field!` is not extended for the `co2` field of " * name(sim) * ": skipping update.") update_field!(
+        sim,
+        Val(:co2),
+        dummy_field,
+    )
+    @test_logs (
+        :warn,
+        "`update_field!` is not extended for the `liquid_precipitation` field of " * name(sim) * ": skipping update.",
+    ) update_field!(sim, Val(:liquid_precipitation), dummy_field)
+    @test_logs (
+        :warn,
+        "`update_field!` is not extended for the `radiative_energy_flux_sfc` field of " *
+        name(sim) *
+        ": skipping update.",
+    ) update_field!(sim, Val(:radiative_energy_flux_sfc), dummy_field)
+    @test_logs (
+        :warn,
+        "`update_field!` is not extended for the `snow_precipitation` field of " * name(sim) * ": skipping update.",
+    ) update_field!(sim, Val(:snow_precipitation), dummy_field)
+    @test_logs (
+        :warn,
+        "`update_field!` is not extended for the `surface_albedo` field of " * name(sim) * ": skipping update.",
+    ) update_field!(sim, Val(:surface_albedo), dummy_field)
+    @test_logs (
+        :warn,
+        "`update_field!` is not extended for the `surface_temperature` field of " * name(sim) * ": skipping update.",
+    ) update_field!(sim, Val(:surface_temperature), dummy_field)
+    @test_logs (
+        :warn,
+        "`update_field!` is not extended for the `turbulent_energy_flux` field of " * name(sim) * ": skipping update.",
+    ) update_field!(sim, Val(:turbulent_energy_flux), dummy_field)
+    @test_logs (
+        :warn,
+        "`update_field!` is not extended for the `turbulent_fluxes` field of " * name(sim) * ": skipping update.",
+    ) update_field!(sim, Val(:turbulent_fluxes), dummy_field)
+    @test_logs (
+        :warn,
+        "`update_field!` is not extended for the `turbulent_moisture_flux` field of " *
+        name(sim) *
+        ": skipping update.",
+    ) update_field!(sim, Val(:turbulent_moisture_flux), dummy_field)
+end