@testset "surface radiative flux consistency (AMIP + integrated land)" begin
    # Build AMIP configuration with integrated land
    config_file = joinpath(pkgdir(ClimaCoupler), "config/ci_configs/amip_default.yml")
    config_dict = get_coupler_config_dict(config_file)
    config_dict["land_model"] = "integrated"

    # Make sure radiation is computed during the first step
    config_dict["dt_rad"] = config_dict["dt"]

    # Construct coupled simulation and run one coupling step
    cs = CoupledSimulation(config_dict)
    step!(cs)
    boundary_space = Interfacer.boundary_space(cs)

    # Unpack component models
    (; atmos_sim, land_sim, ocean_sim, ice_sim) = cs.model_sims

    # Atmosphere: radiative flux on the surface interface
    # Convention: positive downward to the surface
    atmos_flux = CC.Spaces.level(
        atmos_sim.integrator.p.radiation.ᶠradiation_flux.components.data.:1,
        CC.Utilities.half,
    )

    # Integrated land: compute net radiation similar to bucket
    # Convention note: integrated land doesn't store R_n directly, so we compute it
    # using the same formula as bucket (positive upward).
    p = land_sim.integrator.p
    Y = land_sim.integrator.u
    FT = eltype(Y)
    land_fraction = Interfacer.get_field(land_sim, Val(:area_fraction))

    # Compute net radiative flux (R_n) similar to bucket
    # TODO: get sigma from parameters
    σ = FT(5.67e-8)
    α = Interfacer.get_field(land_sim, Val(:surface_direct_albedo))
    ϵ = Interfacer.get_field(land_sim, Val(:emissivity))
    T_sfc = Interfacer.get_field(land_sim, Val(:surface_temperature))

    # Compute net radiation on land surface space, then remap to boundary space
    land_flux = @. -(1 - α) * p.drivers.SW_d -
        ϵ * (p.drivers.LW_d - σ * T_sfc ^ 4)
    land_flux = Interfacer.remap(land_flux, boundary_space)
    @. land_flux = ifelse(land_fraction ≈ 0, zero(land_flux), land_flux)

    err_land = @. atmos_flux - land_flux
    @. err_land = ifelse(land_fraction ≈ 0, zero(err_land), err_land)
    @show "Integrated land flux error: $(maximum(abs.(err_land)))"
    @test maximum(abs.(err_land)) < 5

    # Prescribed ice: radiative fluxes aren't stored; compute from cache and compare
    p = ice_sim.integrator.p
    Y = ice_sim.integrator.u
    FT = eltype(Y)

    # Radiative flux toward surface (positive downward)
    # TODO: get sigma from parameters
    σ = FT(5.67e-8)
    (; k_ice, h, T_base, ρ, c, α, ϵ) = p.params
    ice_rad_flux =
        (1 .- α) .* p.SW_d .+
        ϵ .* (p.LW_d .- σ .* Interfacer.get_field(ice_sim, Val(:surface_temperature)) .^ 4)
    @. ice_rad_flux = ifelse(p.area_fraction ≈ 0, zero(ice_rad_flux), ice_rad_flux)
    ice_fraction = Interfacer.get_field(ice_sim, Val(:area_fraction))

    # Prescribed ocean: SST is prescribed, but for this test we can still compute
    # the radiative flux seen by the ocean surface using the same formula.
    α = Interfacer.get_field(ocean_sim, Val(:surface_direct_albedo))
    ϵ = Interfacer.get_field(ocean_sim, Val(:emissivity))
    ocean_rad_flux =
        (1 .- α) .* cs.fields.SW_d .+
        ϵ .* (
            cs.fields.LW_d .-
            σ .* Interfacer.get_field(ocean_sim, Val(:surface_temperature)) .^ 4
        )
    ocean_fraction = Interfacer.get_field(ocean_sim, Val(:area_fraction))
    @. ocean_rad_flux = ifelse(ocean_fraction ≈ 0, zero(ocean_rad_flux), ocean_rad_flux)

    # Combine component fluxes by area-weighted sum (incl. land sign convention):
    combined_fluxes =
        .-land_fraction .* land_flux .+ ice_fraction .* ice_rad_flux .+
        ocean_fraction .* ocean_rad_flux
    err_fluxes = atmos_flux .+ combined_fluxes
    @show "Combined fluxes error: $(maximum(abs.(err_fluxes)))"
    @test maximum(abs.(err_fluxes)) < 8
end