diff --git a/docs/src/fluxcalculator.md b/docs/src/fluxcalculator.md
index f77af1d570..7c1f1b0eb5 100644
--- a/docs/src/fluxcalculator.md
+++ b/docs/src/fluxcalculator.md
@@ -25,6 +25,7 @@ Fluxes over a heterogeneous surface (e.g., from a gridpoint where atmospheric ce
     ClimaCoupler.FluxCalculator.update_turbulent_fluxes_point!
     ClimaCoupler.FluxCalculator.extrapolate_ρ_to_sfc
     ClimaCoupler.FluxCalculator.surface_thermo_state
+    ClimaCoupler.FluxCalculator.water_albedo_from_atmosphere!
 ```
 
 ## FieldExchanger Internal Functions
diff --git a/docs/src/interfacer.md b/docs/src/interfacer.md
index bdf4b44798..1e237d5cdd 100644
--- a/docs/src/interfacer.md
+++ b/docs/src/interfacer.md
@@ -101,11 +101,12 @@ following properties:
 | Coupler name      | Description | Units |
 |-------------------|-------------|-------|
 | `co2`              | global mean co2 | ppm |
-| `surface_albedo`   | bulk surface albedo over the whole surface space | |
+| `surface_direct_albedo`   | bulk direct surface albedo over the whole surface space | |
+| `surface_diffuse_albedo`   | bulk diffuse 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 |
 
-- `calculate_surface_air_density(atmos_sim::Interfacer.AtmosModelSimulation, T_S::Fields.Field)`:
+- `calculate_surface_air_density(atmos_sim::Interfacer.AtmosModelSimulation, T_S::ClimaCore.Fields.Field)`:
 A function to return the air density of the atmosphere simulation
 extrapolated to the surface, with units of [kg m^-3].
 
@@ -124,7 +125,8 @@ for the following properties:
 | `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_direct albedo`    | bulk direct surface albedo | |
+| `surface_diffuse albedo`    | bulk diffuse surface albedo | |
 | `surface_humidity`  | surface humidity | kg kg^-1 |
 | `surface_temperature` | surface temperature | K |
 
@@ -170,7 +172,7 @@ 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, ::Union{Val{:surface_direct_albedo}, Val{:surface_diffuse_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
@@ -180,10 +182,10 @@ 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)
+function update_field!(sim::SurfaceStub, ::Val{:area_fraction}, field::ClimaCore.Fields.Field)
     sim.cache.area_fraction .= field
 end
-function update_field!(sim::SurfaceStub, ::Val{:surface_temperature}, field::Fields.Field)
+function update_field!(sim::SurfaceStub, ::Val{:surface_temperature}, field::ClimaCore.Fields.Field)
     sim.cache.T_sfc .= field
 end
 ```
diff --git a/experiments/AMIP/Manifest.toml b/experiments/AMIP/Manifest.toml
index e1cbe1ea60..eb90b68d1c 100644
--- a/experiments/AMIP/Manifest.toml
+++ b/experiments/AMIP/Manifest.toml
@@ -150,9 +150,9 @@ version = "0.4.2"
 
 [[deps.BandedMatrices]]
 deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "PrecompileTools"]
-git-tree-sha1 = "fe7ed7b33416db994f706f66e31bbd1b378f71b7"
+git-tree-sha1 = "c946c5014cf4cdbfacacb363b110e7bffba3e742"
 uuid = "aae01518-5342-5314-be14-df237901396f"
-version = "1.6.0"
+version = "1.6.1"
 weakdeps = ["SparseArrays"]
 
     [deps.BandedMatrices.extensions]
@@ -397,9 +397,9 @@ weakdeps = ["ClimaParams"]
 
 [[deps.ClimaParams]]
 deps = ["DocStringExtensions", "TOML", "Test"]
-git-tree-sha1 = "ec67949db856e01df4cbf7d6ddafefeda02f93ee"
+git-tree-sha1 = "1a3d2455fff201bcf130bbd5a71ac16fc3c21fd1"
 uuid = "5c42b081-d73a-476f-9059-fd94b934656c"
-version = "0.10.3"
+version = "0.10.4"
 
 [[deps.ClimaTimeSteppers]]
 deps = ["ClimaComms", "Colors", "DataStructures", "DiffEqBase", "DiffEqCallbacks", "KernelAbstractions", "Krylov", "LinearAlgebra", "LinearOperators", "NVTX", "SciMLBase", "StaticArrays"]
@@ -500,7 +500,7 @@ weakdeps = ["Dates", "LinearAlgebra"]
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.1.0+0"
+version = "1.0.5+1"
 
 [[deps.CompositionsBase]]
 git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad"
@@ -850,9 +850,9 @@ version = "0.3.1"
 
 [[deps.FileIO]]
 deps = ["Pkg", "Requires", "UUIDs"]
-git-tree-sha1 = "c5c28c245101bd59154f649e19b038d15901b5dc"
+git-tree-sha1 = "82d8afa92ecf4b52d78d869f038ebfb881267322"
 uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-version = "1.16.2"
+version = "1.16.3"
 
 [[deps.FilePaths]]
 deps = ["FilePathsBase", "MacroTools", "Reexport", "Requires"]
@@ -905,9 +905,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 = "5a626d6ef24ae0a8590c22dc12096fb65eb66325"
+git-tree-sha1 = "502d0232ec6430d40b6e5b57637333f32192592e"
 uuid = "587475ba-b771-5e3f-ad9e-33799f191a9c"
-version = "0.14.13"
+version = "0.14.14"
 
     [deps.Flux.extensions]
     FluxAMDGPUExt = "AMDGPU"
@@ -1127,9 +1127,9 @@ version = "1.3.14+0"
 
 [[deps.GridLayoutBase]]
 deps = ["GeometryBasics", "InteractiveUtils", "Observables"]
-git-tree-sha1 = "af13a277efd8a6e716d79ef635d5342ccb75be61"
+git-tree-sha1 = "6f93a83ca11346771a93bbde2bdad2f65b61498f"
 uuid = "3955a311-db13-416c-9275-1d80ed98e5e9"
-version = "0.10.0"
+version = "0.10.2"
 
 [[deps.Grisu]]
 git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
@@ -1154,9 +1154,9 @@ version = "1.14.3+1"
 
 [[deps.HTTP]]
 deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"]
-git-tree-sha1 = "db864f2d91f68a5912937af80327d288ea1f3aee"
+git-tree-sha1 = "995f762e0182ebc50548c434c171a5bb6635f8e4"
 uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
-version = "1.10.3"
+version = "1.10.4"
 
 [[deps.HarfBuzz_jll]]
 deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"]
@@ -1430,9 +1430,9 @@ version = "3.0.0+1"
 
 [[deps.LLVM]]
 deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Preferences", "Printf", "Requires", "Unicode"]
-git-tree-sha1 = "ddab4d40513bce53c8e3157825e245224f74fae7"
+git-tree-sha1 = "7c6650580b4c3169d9905858160db895bff6d2e2"
 uuid = "929cbde3-209d-540e-8aea-75f648917ca0"
-version = "6.6.0"
+version = "6.6.1"
 weakdeps = ["BFloat16s"]
 
     [deps.LLVM.extensions]
@@ -1926,7 +1926,7 @@ version = "0.3.24+0"
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.23+4"
+version = "0.3.23+2"
 
 [[deps.OpenEXR]]
 deps = ["Colors", "FileIO", "OpenEXR_jll"]
@@ -2405,10 +2405,10 @@ uuid = "76ed43ae-9a5d-5a62-8c75-30186b810ce8"
 version = "3.45.0+0"
 
 [[deps.SciMLBase]]
-deps = ["ADTypes", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables"]
-git-tree-sha1 = "3a281a9fce9cd62b849d7f16e412933a5fe755cb"
+deps = ["ADTypes", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "SciMLStructures", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables"]
+git-tree-sha1 = "dce2d07e3519cb7817f2d89a7978c13fef30be87"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-version = "2.29.0"
+version = "2.30.0"
 
     [deps.SciMLBase.extensions]
     SciMLBaseChainRulesCoreExt = "ChainRulesCore"
@@ -2435,6 +2435,11 @@ git-tree-sha1 = "10499f619ef6e890f3f4a38914481cc868689cd5"
 uuid = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 version = "0.3.8"
 
+[[deps.SciMLStructures]]
+git-tree-sha1 = "5833c10ce83d690c124beedfe5f621b50b02ba4d"
+uuid = "53ae85a6-f571-4167-b2af-e1d143709226"
+version = "1.1.0"
+
 [[deps.Scratch]]
 deps = ["Dates"]
 git-tree-sha1 = "3bac05bc7e74a75fd9cba4295cde4045d9fe2386"
@@ -2791,9 +2796,9 @@ uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 version = "0.5.23"
 
 [[deps.TranscodingStreams]]
-git-tree-sha1 = "3caa21522e7efac1ba21834a03734c57b4611c7e"
+git-tree-sha1 = "a09c933bebed12501890d8e92946bbab6a1690f1"
 uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
-version = "0.10.4"
+version = "0.10.5"
 weakdeps = ["Random", "Test"]
 
     [deps.TranscodingStreams.extensions]
diff --git a/experiments/AMIP/components/atmosphere/climaatmos.jl b/experiments/AMIP/components/atmosphere/climaatmos.jl
index 08ba2e717a..27f64a5f40 100644
--- a/experiments/AMIP/components/atmosphere/climaatmos.jl
+++ b/experiments/AMIP/components/atmosphere/climaatmos.jl
@@ -1,6 +1,7 @@
 # atmos_init: for ClimaAtmos pre-AMIP interface
 using StaticArrays
 using Statistics: mean
+using LinearAlgebra: norm
 
 import ClimaAtmos as CA
 import ClimaAtmos: CT1, CT2, CT12, CT3, C3, C12, unit_basis_vector_data, ⊗
@@ -14,7 +15,8 @@ import ClimaCoupler.FluxCalculator:
     calculate_surface_air_density,
     PartitionedStateFluxes,
     extrapolate_ρ_to_sfc,
-    get_surface_params
+    get_surface_params,
+    water_albedo_from_atmosphere!
 import ClimaCoupler.Interfacer: get_field, update_field!, name
 import ClimaCoupler.Checkpointer: get_model_prog_state
 import ClimaCoupler.FieldExchanger: update_sim!, step!, reinit!
@@ -35,6 +37,7 @@ name(::ClimaAtmosSimulation) = "ClimaAtmosSimulation"
 
 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_dict["surface_albedo"] = "CouplerAlbedo"
     atmos_config = CA.AtmosConfig(atmos_config_dict)
     simulation = CA.get_simulation(atmos_config)
     (; integrator) = simulation
@@ -171,13 +174,16 @@ 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{:surface_albedo}, field)
-    sim.integrator.p.radiation.radiation_model.diffuse_sw_surface_albedo .=
-        reshape(CA.RRTMGPI.field2array(field), 1, length(parent(field)))
+function update_field!(sim::ClimaAtmosSimulation, ::Val{:surface_direct_albedo}, field)
     sim.integrator.p.radiation.radiation_model.direct_sw_surface_albedo .=
         reshape(CA.RRTMGPI.field2array(field), 1, length(parent(field)))
 end
 
+function update_field!(sim::ClimaAtmosSimulation, ::Val{:surface_diffuse_albedo}, field)
+    sim.integrator.p.radiation.radiation_model.diffuse_sw_surface_albedo .=
+        reshape(CA.RRTMGPI.field2array(field), 1, length(parent(field)))
+end
+
 function update_field!(sim::ClimaAtmosSimulation, ::Val{:turbulent_fluxes}, fields)
     (; F_turb_energy, F_turb_moisture, F_turb_ρτxz, F_turb_ρτyz) = fields
 
@@ -210,7 +216,8 @@ 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(:surface_albedo), csf.surface_albedo)
+    update_field!(atmos_sim, Val(:surface_direct_albedo), csf.surface_direct_albedo)
+    update_field!(atmos_sim, Val(:surface_diffuse_albedo), csf.surface_diffuse_albedo)
     update_field!(atmos_sim, Val(:surface_temperature), csf)
 
     if turbulent_fluxes isa PartitionedStateFluxes
@@ -384,3 +391,37 @@ function dss_state!(sim::ClimaAtmosSimulation)
         ClimaCore.Spaces.weighted_dss!(field, buffer)
     end
 end
+
+"""
+    water_albedo_from_atmosphere!(atmos_sim::ClimaAtmosSimulation, direct_albedo::ClimaCore.Fields.Field, diffuse_albedo::ClimaCore.Fields.Field)
+
+Extension to calculate the water surface albedo from wind speed and insolation. It can be used for prescribed ocean and lakes.
+"""
+function water_albedo_from_atmosphere!(
+    atmos_sim::ClimaAtmosSimulation,
+    direct_albedo::ClimaCore.Fields.Field,
+    diffuse_albedo::ClimaCore.Fields.Field,
+)
+
+    Y = atmos_sim.integrator.u
+    p = atmos_sim.integrator.p
+    t = atmos_sim.integrator.t
+
+    radiation_model = atmos_sim.integrator.p.radiation.radiation_model
+    FT = eltype(Y)
+    λ = FT(0) # spectral wavelength (not used for now)
+
+    # update for current zenith angle
+    CA.set_insolation_variables!(Y, p, t)
+
+    bottom_coords = ClimaCore.Fields.coordinate_field(ClimaCore.Spaces.level(Y.c, 1))
+    μ = CA.RRTMGPI.array2field(radiation_model.cos_zenith, axes(bottom_coords))
+    FT = eltype(atmos_sim.integrator.u)
+    α_model = CA.RegressionFunctionAlbedo{FT}()
+
+
+    # set the direct and diffuse surface albedos
+    direct_albedo .= CA.surface_albedo_direct(α_model).(λ, μ, norm.(ClimaCore.Fields.level(Y.c.uₕ, 1)))
+    diffuse_albedo .= CA.surface_albedo_diffuse(α_model).(λ, μ, norm.(ClimaCore.Fields.level(Y.c.uₕ, 1)))
+
+end
diff --git a/experiments/AMIP/components/land/climaland_bucket.jl b/experiments/AMIP/components/land/climaland_bucket.jl
index 4b4d9cd64f..916e20b57a 100644
--- a/experiments/AMIP/components/land/climaland_bucket.jl
+++ b/experiments/AMIP/components/land/climaland_bucket.jl
@@ -145,7 +145,9 @@ 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}) =
+get_field(sim::BucketSimulation, ::Val{:surface_direct_albedo}) =
+    ClimaLand.surface_albedo(sim.model, sim.integrator.u, sim.integrator.p)
+get_field(sim::BucketSimulation, ::Val{:surface_diffuse_albedo}) =
     ClimaLand.surface_albedo(sim.model, sim.integrator.u, sim.integrator.p)
 get_field(sim::BucketSimulation, ::Val{:surface_humidity}) =
     ClimaLand.surface_specific_humidity(sim.model, sim.integrator.u, sim.integrator.p, sim.integrator.t)
diff --git a/experiments/AMIP/components/ocean/eisenman_seaice.jl b/experiments/AMIP/components/ocean/eisenman_seaice.jl
index 0bfe3bd678..75ac89f91e 100644
--- a/experiments/AMIP/components/ocean/eisenman_seaice.jl
+++ b/experiments/AMIP/components/ocean/eisenman_seaice.jl
@@ -102,7 +102,7 @@ get_field(sim::EisenmanIceSimulation, ::Val{:roughness_buoyancy}) =
 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}) =
+get_field(sim::EisenmanIceSimulation, ::Union{Val{:surface_direct_albedo}, Val{:surface_diffuse_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
diff --git a/experiments/AMIP/components/ocean/prescr_seaice.jl b/experiments/AMIP/components/ocean/prescr_seaice.jl
index 88d499b248..47807fa5c6 100644
--- a/experiments/AMIP/components/ocean/prescr_seaice.jl
+++ b/experiments/AMIP/components/ocean/prescr_seaice.jl
@@ -105,7 +105,8 @@ get_field(sim::PrescribedIceSimulation, ::Val{:area_fraction}) = sim.integrator.
 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, ::Union{Val{:surface_direct_albedo}, Val{:surface_diffuse_albedo}}) =
+    sim.integrator.p.params.α
 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
diff --git a/experiments/AMIP/components/ocean/slab_ocean.jl b/experiments/AMIP/components/ocean/slab_ocean.jl
index 2a708bb52b..fa03a7471c 100644
--- a/experiments/AMIP/components/ocean/slab_ocean.jl
+++ b/experiments/AMIP/components/ocean/slab_ocean.jl
@@ -115,7 +115,8 @@ get_field(sim::SlabOceanSimulation, ::Val{:area_fraction}) = sim.integrator.p.ar
 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, ::Union{Val{:surface_direct_albedo}, Val{:surface_diffuse_albedo}}) =
+    sim.integrator.p.params.α
 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
diff --git a/experiments/AMIP/coupler_driver.jl b/experiments/AMIP/coupler_driver.jl
index 8fff9b38c8..fa47c2636c 100644
--- a/experiments/AMIP/coupler_driver.jl
+++ b/experiments/AMIP/coupler_driver.jl
@@ -58,7 +58,8 @@ using ClimaCoupler.FluxCalculator:
     CombinedStateFluxes,
     combined_turbulent_fluxes!,
     MoninObukhovScheme,
-    partitioned_turbulent_fluxes!
+    partitioned_turbulent_fluxes!,
+    water_albedo_from_atmosphere!
 using ClimaCoupler.Interfacer: CoupledSimulation, SurfaceStub, get_field, update_field!
 using ClimaCoupler.Regridder
 using ClimaCoupler.Regridder: update_surface_fractions!, combine_surfaces!, binary_mask
@@ -133,6 +134,7 @@ hourly_checkpoint = config_dict["hourly_checkpoint"]
 restart_dir = config_dict["restart_dir"]
 restart_t = Int(config_dict["restart_t"])
 evolving_ocean = config_dict["evolving_ocean"]
+dt_rad = config_dict["dt_rad"]
 
 #=
 ## Setup Communication Context
@@ -280,7 +282,8 @@ if mode_name == "amip"
         z0m = FT(1e-3),
         z0b = FT(1e-3),
         beta = FT(1),
-        α = FT(0.06),
+        α_direct = ClimaCore.Fields.ones(boundary_space) .* FT(0.06),
+        α_diffuse = ClimaCore.Fields.ones(boundary_space) .* FT(0.06),
         area_fraction = (FT(1) .- land_fraction),
         phase = TD.Liquid(),
         thermo_params = thermo_params,
@@ -374,7 +377,8 @@ elseif mode_name in ("slabplanet", "slabplanet_aqua", "slabplanet_terra")
         z0m = FT(0),
         z0b = FT(0),
         beta = FT(1),
-        α = FT(1),
+        α_direct = ClimaCore.Fields.ones(boundary_space) .* FT(1),
+        α_diffuse = ClimaCore.Fields.ones(boundary_space) .* FT(1),
         area_fraction = ClimaCore.Fields.zeros(boundary_space),
         phase = TD.Ice(),
         thermo_params = thermo_params,
@@ -438,7 +442,8 @@ coupler_field_names = (
     :z0b_S,
     :ρ_sfc,
     :q_sfc,
-    :surface_albedo,
+    :surface_direct_albedo,
+    :surface_diffuse_albedo,
     :beta,
     :F_turb_energy,
     :F_turb_moisture,
@@ -449,6 +454,8 @@ coupler_field_names = (
     :P_snow,
     :radiative_energy_flux_toa,
     :P_net,
+    :temp1,
+    :temp2,
 )
 coupler_fields =
     NamedTuple{coupler_field_names}(ntuple(i -> ClimaCore.Fields.zeros(boundary_space), length(coupler_field_names)))
@@ -515,14 +522,23 @@ the `ref_date` field specifies the reference (first) date for the callback, and
 The currently implemented callbacks are:
 - `checkpoint_cb`: generates a checkpoint of all model states at a specified interval. This is mainly used for restarting simulations.
 - `update_firstdayofmonth!_cb`: generates a callback to update the first day of the month for monthly message print (and other monthly operations).
+- `albedo_cb`: for the amip mode, the water albedo is time varying (since the reflectivity of water depends on insolation and wave characteristics, with the latter
+  being approximated from wind speed). It is updated at the same frequency as the atmospheric radiation.
+  NB: Eventually, we will call all of radiation from the coupler, in addition to the albedo calculation.
 =#
-
-## checkpoint_cb generates a checkpoint of all model states at a specified interval. This mainly used for restarting simulations.
 checkpoint_cb =
     HourlyCallback(dt = FT(480), func = checkpoint_sims, ref_date = [dates.date[1]], active = hourly_checkpoint) # 20 days
 update_firstdayofmonth!_cb =
     MonthlyCallback(dt = FT(1), func = update_firstdayofmonth!, ref_date = [dates.date1[1]], active = true)
-callbacks = (; checkpoint = checkpoint_cb, update_firstdayofmonth! = update_firstdayofmonth!_cb)
+dt_water_albedo = parse(FT, filter(x -> !occursin(x, "hours"), dt_rad))
+albedo_cb = HourlyCallback(
+    dt = dt_water_albedo,
+    func = water_albedo_from_atmosphere!,
+    ref_date = [dates.date[1]],
+    active = mode_name == "amip",
+)
+callbacks =
+    (; checkpoint = checkpoint_cb, update_firstdayofmonth! = update_firstdayofmonth!_cb, water_albedo = albedo_cb)
 
 #=
 ## Initialize turbulent fluxes
@@ -690,9 +706,12 @@ function solve_coupler!(cs)
 
         ## compute global energy
         !isnothing(cs.conservation_checks) ? check_conservation!(cs) : nothing
+        ClimaComms.barrier(comms_ctx)
+
+        ## update water albedo from wind at dt_water_albedo (this will be extended to a radiation callback from the coupler)
+        trigger_callback!(cs, cs.callbacks.water_albedo)
 
         ## run component models sequentially for one coupling timestep (Δt_cpl)
-        ClimaComms.barrier(comms_ctx)
         update_surface_fractions!(cs)
         update_model_sims!(cs.model_sims, cs.fields, cs.turbulent_fluxes)
 
diff --git a/experiments/AMIP/user_io/debug_plots.jl b/experiments/AMIP/user_io/debug_plots.jl
index c4f5f6b8c6..4cc8fc9d2b 100644
--- a/experiments/AMIP/user_io/debug_plots.jl
+++ b/experiments/AMIP/user_io/debug_plots.jl
@@ -29,7 +29,8 @@ plot the anomalies of the fields with respect to `cs_fields_ref`.
 """
 function debug(cs_fields::NamedTuple, dir, cs_fields_ref = nothing)
     field_names = (
-        :surface_albedo,
+        :surface_direct_albedo,
+        :surface_diffuse_albedo,
         :F_radiative,
         :F_turb_energy,
         :F_turb_moisture,
diff --git a/src/FieldExchanger.jl b/src/FieldExchanger.jl
index ec5ca5d922..c38558415f 100644
--- a/src/FieldExchanger.jl
+++ b/src/FieldExchanger.jl
@@ -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`, `:surface_albedo`, `:roughness_momentum`, `:roughness_buoyancy`, `:beta`)
+(`:surface_temperature`, `:surface_direct_albedo`, `:surface_diffuse_albedo`, `:roughness_momentum`, `:roughness_buoyancy`, `:beta`)
 need to be specified for each surface model.
 
 # Arguments
@@ -67,8 +67,11 @@ 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(:surface_albedo))
-    Regridder.dummmy_remap!(csf.surface_albedo, combined_field)
+    Regridder.combine_surfaces!(combined_field, model_sims, Val(:surface_direct_albedo))
+    Regridder.dummmy_remap!(csf.surface_direct_albedo, combined_field)
+
+    Regridder.combine_surfaces!(combined_field, model_sims, Val(:surface_diffuse_albedo))
+    Regridder.dummmy_remap!(csf.surface_diffuse_albedo, combined_field)
 
     if turbulent_fluxes isa FluxCalculator.CombinedStateFluxes
         Regridder.combine_surfaces!(combined_field, model_sims, Val(:roughness_momentum))
@@ -97,7 +100,8 @@ 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(:surface_albedo), csf.surface_albedo)
+    Interfacer.update_field!(atmos_sim, Val(:surface_direct_albedo), csf.surface_direct_albedo)
+    Interfacer.update_field!(atmos_sim, Val(:surface_diffuse_albedo), csf.surface_diffuse_albedo)
     Interfacer.update_field!(atmos_sim, Val(:surface_temperature), csf.T_S)
 
     if turbulent_fluxes isa FluxCalculator.CombinedStateFluxes
diff --git a/src/FluxCalculator.jl b/src/FluxCalculator.jl
index c45e6b9217..b8a924cc00 100644
--- a/src/FluxCalculator.jl
+++ b/src/FluxCalculator.jl
@@ -9,7 +9,6 @@ module FluxCalculator
 import SurfaceFluxes as SF
 import Thermodynamics as TD
 using StaticArrays
-
 using ClimaCoupler: Interfacer, Regridder
 using ClimaCore: Fields, Spaces
 export PartitionedStateFluxes,
@@ -23,7 +22,8 @@ export PartitionedStateFluxes,
     BulkScheme,
     partitioned_turbulent_fluxes!,
     get_surface_params,
-    update_turbulent_fluxes_point!
+    update_turbulent_fluxes_point!,
+    water_albedo_from_atmosphere!
 
 """
     TurbulentFluxPartition
@@ -387,4 +387,31 @@ surface temperature in surface energy balance calculations.
 """
 differentiate_turbulent_fluxes!(::Interfacer.SurfaceModelSimulation, args) = nothing
 
+"""
+    water_albedo_from_atmosphere!(cs::Interfacer.CoupledSimulation, _)
+
+Callback to calculate the water albedo from atmospheric state. This is a placeholder for the full radiation callback.
+"""
+function water_albedo_from_atmosphere!(cs::Interfacer.CoupledSimulation, _)
+    atmos_sim = cs.model_sims.atmos_sim
+    ocean_sim = cs.model_sims.ocean_sim
+    cf = cs.fields
+
+    # use temp fields
+    water_albedo_from_atmosphere!(atmos_sim, cf.temp1, cf.temp2)
+
+    Interfacer.update_field!(ocean_sim, Val(:surface_direct_albedo), cf.temp1)
+    Interfacer.update_field!(ocean_sim, Val(:surface_diffuse_albedo), cf.temp2)
+
+end
+
+"""
+    water_albedo_from_atmosphere!(atmos_sim::Interfacer.AtmosModelSimulation, ::Fields.Field, ::Fields.Field)
+
+Placeholder for the water albedo calculation from the atmosphere. It returns an error if not extended.
+"""
+function water_albedo_from_atmosphere!(atmos_sim::Interfacer.AtmosModelSimulation, ::Fields.Field, ::Fields.Field)
+    error("this function is required to be dispatched on" * Interfacer.name(atmos_sim) * ", but no method defined")
+end
+
 end # module
diff --git a/src/Interfacer.jl b/src/Interfacer.jl
index 18e59368f0..530b123d01 100644
--- a/src/Interfacer.jl
+++ b/src/Interfacer.jl
@@ -147,7 +147,8 @@ get_field(
         Val{:beta},
         Val{:roughness_buoyancy},
         Val{:roughness_momentum},
-        Val{:surface_albedo},
+        Val{:surface_direct_albedo},
+        Val{:surface_diffuse_albedo},
         Val{:surface_humidity},
         Val{:surface_temperature},
     },
@@ -190,7 +191,8 @@ 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_direct_albedo}) = sim.cache.α_direct
+get_field(sim::SurfaceStub, ::Val{:surface_diffuse_albedo}) = sim.cache.α_diffuse
 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
@@ -218,7 +220,13 @@ If it isn't extended, the field won't be updated and a warning will be raised.
 """
 update_field!(
     sim::AtmosModelSimulation,
-    val::Union{Val{:co2}, Val{:surface_albedo}, Val{:surface_temperature}, Val{:turbulent_fluxes}},
+    val::Union{
+        Val{:co2},
+        Val{:surface_direct_albedo},
+        Val{:surface_diffuse_albedo},
+        Val{:surface_temperature},
+        Val{:turbulent_fluxes},
+    },
     _,
 ) = update_field_warning(sim, val)
 
@@ -239,6 +247,8 @@ update_field!(
         Val{:snow_precipitation},
         Val{:turbulent_energy_flux},
         Val{:turbulent_moisture_flux},
+        Val{:surface_direct_albedo},
+        Val{:surface_diffuse_albedo},
     },
     _,
 ) = update_field_warning(sim, val)
@@ -260,6 +270,12 @@ end
 function update_field!(sim::SurfaceStub, ::Val{:air_density}, field)
     parent(sim.cache.ρ_sfc) .= parent(field)
 end
+function update_field!(sim::SurfaceStub, ::Val{:surface_direct_albedo}, field::Fields.Field)
+    sim.cache.α_direct .= field
+end
+function update_field!(sim::SurfaceStub, ::Val{:surface_diffuse_albedo}, field::Fields.Field)
+    sim.cache.α_diffuse .= field
+end
 
 """
     name(::ComponentModelSimulation)
diff --git a/test/Project.toml b/test/Project.toml
index cead746cae..1b19987fba 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -17,6 +17,7 @@ GeoMakie = "db073c08-6b98-4ee5-b6a4-5efafb3259c6"
 Insolation = "e98cc03f-d57e-4e3c-b70c-8d51efe9e0d8"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 MPIPreferences = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
diff --git a/test/debug/debug_amip_plots.jl b/test/debug/debug_amip_plots.jl
index f9ce30f5ad..8470c389ba 100644
--- a/test/debug/debug_amip_plots.jl
+++ b/test/debug/debug_amip_plots.jl
@@ -39,7 +39,8 @@ plot_field_names(sim::SurfaceStub) = (:stub_field,)
 
     boundary_space = TestHelper.create_space(FT)
     coupler_names = (
-        :surface_albedo,
+        :surface_direct_albedo,
+        :surface_diffuse_albedo,
         :F_radiative,
         :F_turb_energy,
         :F_turb_moisture,
diff --git a/test/field_exchanger_tests.jl b/test/field_exchanger_tests.jl
index 5ebba43db5..d3fca071a6 100644
--- a/test/field_exchanger_tests.jl
+++ b/test/field_exchanger_tests.jl
@@ -41,8 +41,10 @@ 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{:surface_albedo}) =
-    sim.cache_field .* eltype(sim.cache_field)(1.0)
+get_field(
+    sim::Union{TestSurfaceSimulation1, TestSurfaceSimulation2},
+    ::Union{Val{:surface_direct_albedo}, Val{:surface_diffuse_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)
 get_field(sim::Union{TestSurfaceSimulation1, TestSurfaceSimulation2}, ::Val{:roughness_buoyancy}) =
@@ -63,8 +65,11 @@ get_field(sim::Union{TestSurfaceSimulation2, TestSurfaceSimulation2}, ::Val{:sur
 struct TestAtmosSimulation{C} <: AtmosModelSimulation
     cache::C
 end
-function update_field!(sim::TestAtmosSimulation, ::Val{:surface_albedo}, field)
-    parent(sim.cache.surface_albedo) .= parent(field)
+function update_field!(sim::TestAtmosSimulation, ::Val{:surface_direct_albedo}, field)
+    parent(sim.cache.albedo_direct) .= parent(field)
+end
+function update_field!(sim::TestAtmosSimulation, ::Val{:surface_diffuse_albedo}, field)
+    parent(sim.cache.albedo_diffuse) .= parent(field)
 end
 function update_field!(sim::TestAtmosSimulation, ::Val{:roughness_momentum}, field)
     parent(sim.cache.roughness_momentum) .= parent(field)
@@ -119,10 +124,17 @@ 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, :surface_albedo, :beta, :q_sfc)
+        coupler_names = (:T_S, :z0m_S, :z0b_S, :surface_direct_albedo, :surface_diffuse_albedo, :beta, :q_sfc)
 
         # coupler cache setup
-        exchanged_fields = (:surface_temperature, :surface_albedo, :roughness_momentum, :roughness_buoyancy, :beta)
+        exchanged_fields = (
+            :surface_temperature,
+            :surface_direct_albedo,
+            :surface_diffuse_albedo,
+            :roughness_momentum,
+            :roughness_buoyancy,
+            :beta,
+        )
 
         sims = (; a = TestSurfaceSimulation1(ones(boundary_space)), b = TestSurfaceSimulation2(ones(boundary_space)))
 
@@ -134,7 +146,8 @@ 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.surface_albedo)[1] == results[1]
+            @test parent(coupler_fields.surface_direct_albedo)[1] == results[1]
+            @test parent(coupler_fields.surface_diffuse_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]
@@ -149,7 +162,8 @@ for FT in (Float32, Float64)
             :T_S,
             :z0m_S,
             :z0b_S,
-            :surface_albedo,
+            :surface_direct_albedo,
+            :surface_diffuse_albedo,
             :beta,
             :F_turb_energy,
             :F_turb_moisture,
@@ -162,7 +176,8 @@ for FT in (Float32, Float64)
 
         # model cache setup
 
-        atmos_names = (:surface_temperature, :surface_albedo, :roughness_momentum, :roughness_buoyancy, :beta)
+        atmos_names =
+            (:surface_temperature, :albedo_direct, :albedo_diffuse, :roughness_momentum, :roughness_buoyancy, :beta)
         atmos_fields = NamedTuple{atmos_names}(ntuple(i -> Fields.zeros(boundary_space), length(atmos_names)))
 
         land_names = (
@@ -181,18 +196,22 @@ for FT in (Float32, Float64)
             stub_sim = SurfaceStub((;
                 area_fraction = Fields.ones(boundary_space),
                 ρ_sfc = Fields.ones(boundary_space),
+                albedo_direct = Fields.ones(boundary_space),
+                albedo_diffuse = Fields.ones(boundary_space),
             )),
         )
+        coupler_fields.surface_diffuse_albedo .= FT(0.5)
 
         # test the sim update under CombinedStateFluxes (update all) and PartitionedStateFluxes (update all except turbulent fluxes)
         flux_types = (CombinedStateFluxes(), PartitionedStateFluxes())
-        results = [FT(0), FT(1)]
+        results = [FT(0), FT(1), FT(0.5)]
         for (i, t) in enumerate(flux_types)
             model_sims.atmos_sim.cache.roughness_momentum .= FT(0)
             update_model_sims!(model_sims, coupler_fields, t)
 
             # test atmos
-            @test parent(model_sims.atmos_sim.cache.surface_albedo)[1] == results[2]
+            @test parent(model_sims.atmos_sim.cache.albedo_direct)[1] == results[2]
+            @test parent(model_sims.atmos_sim.cache.albedo_diffuse)[1] == results[3]
             if t isa CombinedStateFluxes
                 @test parent(model_sims.atmos_sim.cache.roughness_momentum)[1] == results[2]
             else
@@ -213,6 +232,10 @@ for FT in (Float32, Float64)
             @test parent(model_sims.land_sim.cache.liquid_precipitation)[1] == results[1]
             @test parent(model_sims.land_sim.cache.snow_precipitation)[1] == results[1]
 
+            # test stub - albedo should be updated by update_sim!
+            @test parent(model_sims.stub_sim.cache.albedo_direct)[1] == results[2]
+            @test parent(model_sims.stub_sim.cache.albedo_diffuse)[1] == results[2]
+
         end
     end
     @testset "reinit_model_sims! for FT=$FT" begin
diff --git a/test/flux_calculator_tests.jl b/test/flux_calculator_tests.jl
index 7b00324250..078f9c0956 100644
--- a/test/flux_calculator_tests.jl
+++ b/test/flux_calculator_tests.jl
@@ -1,6 +1,6 @@
 
 using ClimaCore: Meshes, Domains, Topologies, Spaces, Fields, InputOutput, Geometry
-using ClimaCoupler: Utilities, Regridder, TestHelper
+using ClimaCoupler: Utilities, Regridder, TestHelper, Interfacer
 using Test
 import ClimaCoupler.FluxCalculator:
     atmos_turbulent_fluxes!,
@@ -17,16 +17,15 @@ import ClimaCoupler.FluxCalculator:
     surface_inputs,
     get_surface_fluxes_point!,
     get_scheme_properties,
-    surface_thermo_state
+    surface_thermo_state,
+    water_albedo_from_atmosphere!
 import ClimaCoupler: Interfacer
-
 import ClimaParams as CP
 import Thermodynamics as TD
 import Thermodynamics.Parameters.ThermodynamicsParameters
 import SurfaceFluxes as SF
 import SurfaceFluxes.Parameters.SurfaceFluxesParameters
 import SurfaceFluxes.UniversalFunctions as UF
-
 using StaticArrays
 
 # simple generic atmos model
@@ -50,6 +49,8 @@ struct TestAtmos{P, Y, D, I} <: Interfacer.AtmosModelSimulation
     domain::D
     integrator::I
 end
+struct TestAtmos2 <: Interfacer.AtmosModelSimulation end
+name(sim::TestAtmos2) = "TestAtmos2"
 
 Interfacer.get_field(sim::TestAtmos, ::Val{:height_int}) = sim.integrator.p.z
 Interfacer.get_field(sim::TestAtmos, ::Val{:height_sfc}) = sim.integrator.p.z_sfc
@@ -96,7 +97,8 @@ 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{:surface_albedo}) = sim.integrator.p.α
+Interfacer.get_field(sim::TestOcean, ::Union{Val{:surface_direct_albedo}, Val{:surface_diffuse_albedo}}) =
+    sim.integrator.p.α
 
 function surface_thermo_state(
     sim::TestOcean,
@@ -143,6 +145,10 @@ function surface_thermo_state(
     @. TD.PhaseEquil_ρTq.(thermo_params, ρ_sfc, T_sfc, q_sfc)
 end
 
+function water_albedo_from_atmosphere!(::TestAtmos, temp1::Fields.Field, temp2::Fields.Field)
+    temp1 .*= 2
+    temp2 .*= 3
+end
 
 for FT in (Float32, Float64)
     @testset "combined_turbulent_fluxes! for FT=$FT" begin
@@ -220,7 +226,8 @@ for FT in (Float32, Float64)
 
             coupler_cache_names = (
                 :T_S,
-                :surface_albedo,
+                :surface_direct_albedo,
+                :surface_diffuse_albedo,
                 :F_R_sfc,
                 :F_R_toa,
                 :P_liq,
@@ -322,4 +329,39 @@ for FT in (Float32, Float64)
         @test surface_thermo_state(surface_sim, thermo_params, thermo_state_int[colidx], colidx).ρ ==
               thermo_state_int[colidx].ρ
     end
+
+    @testset "water_albedo_from_atmosphere!" begin
+        boundary_space = TestHelper.create_space(FT)
+        ocean_sim = Interfacer.SurfaceStub((; α_direct = zeros(boundary_space), α_diffuse = zeros(boundary_space)))
+        atmos_sim = TestAtmos(1, 2, 3, 4)
+        coupler_fields = (; temp1 = ones(boundary_space), temp2 = ones(boundary_space))
+        model_sims = (; atmos_sim, ocean_sim)
+        cs = Interfacer.CoupledSimulation{FT}(
+            nothing, # comms_ctx
+            nothing, # dates
+            nothing, # boundary_space
+            coupler_fields, # fields
+            nothing, # parsed_args
+            nothing, # conservation_checks
+            (Int(0), Int(1)), # tspan
+            0, # t
+            0, # Δt_cpl
+            (;), # surface_masks
+            model_sims, # model_sims
+            (;), # mode
+            (), # diagnostics
+            (;), # callbacks
+            (;), # dirs
+            nothing, # turbulent_fluxes
+            nothing, # thermo_params
+        )
+        water_albedo_from_atmosphere!(cs, nothing)
+        @test sum(parent(cs.model_sims.ocean_sim.cache.α_direct) .- parent(ones(boundary_space)) .* 2) == 0
+        @test sum(parent(cs.model_sims.ocean_sim.cache.α_diffuse) .- parent(ones(boundary_space)) .* 3) == 0
+
+        atmos_sim2 = TestAtmos2()
+        @test_throws ErrorException(
+            "this function is required to be dispatched on" * Interfacer.name(atmos_sim2) * ", but no method defined",
+        ) water_albedo_from_atmosphere!(atmos_sim2, ones(boundary_space), ones(boundary_space))
+    end
 end
diff --git a/test/interfacer_tests.jl b/test/interfacer_tests.jl
index 635a173c7f..ef09cfdc4f 100644
--- a/test/interfacer_tests.jl
+++ b/test/interfacer_tests.jl
@@ -82,7 +82,8 @@ for FT in (Float32, Float64)
         stub = SurfaceStub((;
             area_fraction = FT(1),
             T_sfc = FT(280),
-            α = 3,
+            α_direct = 3,
+            α_diffuse = 3,
             z0m = 4,
             z0b = 5,
             beta = 6,
@@ -92,23 +93,37 @@ 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(:surface_albedo)) == 3
+        @test get_field(stub, Val(:surface_direct_albedo)) == 3
+        @test get_field(stub, Val(:surface_diffuse_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
+        @test get_field(stub, Val(:air_density)) == FT(1)
         @test ≈(get_field(stub, Val(:surface_humidity))[1], FT(0.0076), atol = FT(1e-4))
     end
 
     @testset "update_field! the SurfaceStub area_fraction for FT=$FT" begin
         boundary_space = TestHelper.create_space(FT)
 
-        stub = SurfaceStub((; area_fraction = zeros(boundary_space), T_sfc = zeros(boundary_space)))
+        stub = SurfaceStub((;
+            area_fraction = zeros(boundary_space),
+            T_sfc = zeros(boundary_space),
+            α_direct = zeros(boundary_space),
+            α_diffuse = zeros(boundary_space),
+            z0m = zeros(boundary_space),
+            z0b = zeros(boundary_space),
+            beta = zeros(boundary_space),
+        ))
 
         update_field!(stub, Val(:area_fraction), ones(boundary_space))
         update_field!(stub, Val(:surface_temperature), ones(boundary_space) .* 2)
+        update_field!(stub, Val(:surface_direct_albedo), ones(boundary_space) .* 3)
+        update_field!(stub, Val(:surface_diffuse_albedo), ones(boundary_space) .* 4)
 
         @test parent(get_field(stub, Val(:area_fraction)))[1] == FT(1)
         @test parent(get_field(stub, Val(:surface_temperature)))[1] == FT(2)
+        @test parent(get_field(stub, Val(:surface_direct_albedo)))[1] == FT(3)
+        @test parent(get_field(stub, Val(:surface_diffuse_albedo)))[1] == FT(4)
     end
 end
 
@@ -137,7 +152,8 @@ end
         :beta,
         :roughness_buoyancy,
         :roughness_momentum,
-        :surface_albedo,
+        :surface_direct_albedo,
+        :surface_diffuse_albedo,
         :surface_humidity,
         :surface_temperature,
     )
@@ -205,7 +221,7 @@ end
     sim = DummySimulation4(space)
 
     # Test that update_field! gives correct warnings for unextended fields
-    for value in (:co2, :surface_albedo, :surface_temperature, :turbulent_fluxes)
+    for value in (:co2, :surface_direct_albedo, :surface_diffuse_albedo, :surface_temperature, :turbulent_fluxes)
         val = Val(value)
         @test_logs (
             :warn,