diff --git a/src/Diagnostics/Diagnostics.jl b/src/Diagnostics/Diagnostics.jl
index f1d65ae0..179aded1 100644
--- a/src/Diagnostics/Diagnostics.jl
+++ b/src/Diagnostics/Diagnostics.jl
@@ -2,9 +2,7 @@ module Diagnostics
 
 export MixedLayerDepthField, MixedLayerDepthOperand
 export meridional_heat_transport
-export frazil_temperature_flux, net_ocean_temperature_flux, sea_ice_ocean_temperature_flux, atmosphere_ocean_temperature_flux,
-       frazil_heat_flux, net_ocean_heat_flux, sea_ice_ocean_heat_flux, atmosphere_ocean_heat_flux,
-       net_ocean_salinity_flux, sea_ice_ocean_salinity_flux, atmosphere_ocean_salinity_flux,
+export net_ocean_heat_flux, sea_ice_ocean_heat_flux, atmosphere_ocean_heat_flux,
        net_ocean_freshwater_flux, sea_ice_ocean_freshwater_flux, atmosphere_ocean_freshwater_flux
 
 using KernelAbstractions: @index, @kernel
diff --git a/src/Diagnostics/interface_fluxes.jl b/src/Diagnostics/interface_fluxes.jl
index ef2ba76c..327af9c7 100644
--- a/src/Diagnostics/interface_fluxes.jl
+++ b/src/Diagnostics/interface_fluxes.jl
@@ -19,11 +19,10 @@ Return the two-dimensional frazil temperature flux (K m s⁻¹) in a coupled `es
 function frazil_temperature_flux(esm::EarthSystemModel)
     ρᵒᶜ = esm.interfaces.ocean_properties.reference_density
     cᵒᶜ = esm.interfaces.ocean_properties.heat_capacity
-    frazil_temperature_flux = 1 / (ρᵒᶜ * cᵒᶜ) * frazil_heat_flux(esm)
-    return Field(frazil_temperature_flux)
+    return 1 / (ρᵒᶜ * cᵒᶜ) * frazil_heat_flux(esm)
 end
 
-frazil_temperature_flux(esm::NoSeaIceOceanInterfaceModel) = ZeroField()
+frazil_temperature_flux(::NoSeaIceOceanInterfaceModel) = ZeroField()
 
 """
     net_ocean_temperature_flux(esm::EarthSystemModel)
@@ -32,8 +31,7 @@ Return the net temperature flux (K m s⁻¹) at the ocean's surface in a coupled
 """
 function net_ocean_temperature_flux(esm::EarthSystemModel)
     Jᵀ = flux_field(esm.ocean.model.tracers.T.boundary_conditions.top.condition)
-    net_ocean_temperature_flux = Jᵀ + frazil_temperature_flux(esm)
-    return Field(net_ocean_temperature_flux)
+    return Jᵀ + frazil_temperature_flux(esm)
 end
 
 
@@ -46,11 +44,10 @@ in a coupled `esm`.
 function sea_ice_ocean_temperature_flux(esm::EarthSystemModel)
     ρᵒᶜ = esm.interfaces.ocean_properties.reference_density
     cᵒᶜ = esm.interfaces.ocean_properties.heat_capacity
-    sea_ice_ocean_temperature_flux = 1 / (ρᵒᶜ * cᵒᶜ) * sea_ice_ocean_heat_flux(esm)
-    return Field(sea_ice_ocean_temperature_flux)
+    return 1 / (ρᵒᶜ * cᵒᶜ) * sea_ice_ocean_heat_flux(esm)
 end
 
-sea_ice_ocean_temperature_flux(esm::NoSeaIceOceanInterfaceModel) = ZeroField()
+sea_ice_ocean_temperature_flux(::NoSeaIceOceanInterfaceModel) = ZeroField()
 
 """
     atmosphere_ocean_temperature_flux(esm::EarthSystemModel)
@@ -58,11 +55,8 @@ sea_ice_ocean_temperature_flux(esm::NoSeaIceOceanInterfaceModel) = ZeroField()
 Return the atmosphere-ocean temperature flux (K m s⁻¹) at the atmosphere-ocean
 interface in a coupled `esm`.
 """
-function atmosphere_ocean_temperature_flux(esm::EarthSystemModel)
-    atmosphere_ocean_temperature_flux =
+atmosphere_ocean_temperature_flux(esm::EarthSystemModel) =
         net_ocean_temperature_flux(esm) - sea_ice_ocean_temperature_flux(esm)
-    return Field(atmosphere_ocean_temperature_flux)
-end
 
 
 ###########################
@@ -74,12 +68,10 @@ end
 
 Return the two-dimensional frazil heat flux (W m⁻²) in a coupled `esm`.
 """
-function frazil_heat_flux(esm::EarthSystemModel)
-    frazil_heat_flux = esm.interfaces.sea_ice_ocean_interface.fluxes.frazil_heat
-    return frazil_heat_flux
-end
+frazil_heat_flux(esm::EarthSystemModel) =
+    esm.interfaces.sea_ice_ocean_interface.fluxes.frazil_heat
 
-frazil_heat_flux(esm::NoSeaIceOceanInterfaceModel) = ZeroField()
+frazil_heat_flux(::NoSeaIceOceanInterfaceModel) = ZeroField()
 
 """
     net_ocean_heat_flux(esm::EarthSystemModel)
@@ -89,8 +81,7 @@ Return the net heat flux (W m⁻²) at the ocean's surface in a coupled `esm`.
 function net_ocean_heat_flux(esm::EarthSystemModel)
     ρᵒᶜ = esm.interfaces.ocean_properties.reference_density
     cᵒᶜ = esm.interfaces.ocean_properties.heat_capacity
-    net_ocean_heat_flux = ρᵒᶜ * cᵒᶜ * net_ocean_temperature_flux(esm)
-    return Field(net_ocean_heat_flux)
+    return ρᵒᶜ * cᵒᶜ * net_ocean_temperature_flux(esm)
 end
 
 """
@@ -101,11 +92,10 @@ in a coupled `esm`.
 """
 function sea_ice_ocean_heat_flux(esm::EarthSystemModel)
     sea_ice_ocean_fluxes = esm.interfaces.sea_ice_ocean_interface.fluxes
-    sea_ice_ocean_heat_flux = sea_ice_ocean_fluxes.interface_heat + frazil_heat_flux(esm)
-    return Field(sea_ice_ocean_heat_flux)
+    return sea_ice_ocean_fluxes.interface_heat + frazil_heat_flux(esm)
 end
 
-sea_ice_ocean_heat_flux(esm::NoSeaIceOceanInterfaceModel) = ZeroField()
+sea_ice_ocean_heat_flux(::NoSeaIceOceanInterfaceModel) = ZeroField()
 
 """
     atmosphere_ocean_heat_flux(esm::EarthSystemModel)
@@ -116,8 +106,7 @@ interface in a coupled `esm`.
 function atmosphere_ocean_heat_flux(esm::EarthSystemModel)
     ρᵒᶜ = esm.interfaces.ocean_properties.reference_density
     cᵒᶜ = esm.interfaces.ocean_properties.heat_capacity
-    atmosphere_ocean_heat_flux = ρᵒᶜ * cᵒᶜ * atmosphere_ocean_temperature_flux(esm)
-    return Field(atmosphere_ocean_heat_flux)
+    return ρᵒᶜ * cᵒᶜ * atmosphere_ocean_temperature_flux(esm)
 end
 
 
@@ -130,10 +119,8 @@ end
 
 Return the net salinity flux (g/kg m s⁻¹) at the ocean's surface in a coupled `esm`.
 """
-function net_ocean_salinity_flux(esm::EarthSystemModel)
-    Jˢ = flux_field(esm.ocean.model.tracers.S.boundary_conditions.top.condition)
-    return Jˢ
-end
+net_ocean_salinity_flux(esm::EarthSystemModel) =
+    flux_field(esm.ocean.model.tracers.S.boundary_conditions.top.condition)
 
 
 """
@@ -142,12 +129,10 @@ end
 Return the sea ice-ocean salinity flux (g/kg m s⁻¹) at the sea ice-ocean interface
 in a coupled `esm`.
 """
-function sea_ice_ocean_salinity_flux(esm::EarthSystemModel)
-    sea_ice_ocean_salinity_flux = esm.interfaces.sea_ice_ocean_interface.fluxes.salt
-    return sea_ice_ocean_salinity_flux
-end
+sea_ice_ocean_salinity_flux(esm::EarthSystemModel) =
+    esm.interfaces.sea_ice_ocean_interface.fluxes.salt
 
-sea_ice_ocean_salinity_flux(esm::NoSeaIceOceanInterfaceModel) = ZeroField()
+sea_ice_ocean_salinity_flux(::NoSeaIceOceanInterfaceModel) = ZeroField()
 
 """
     atmosphere_ocean_salinity_flux(esm::EarthSystemModel)
@@ -155,11 +140,8 @@ sea_ice_ocean_salinity_flux(esm::NoSeaIceOceanInterfaceModel) = ZeroField()
 Return the atmosphere-ocean salinity flux (g/kg m s⁻¹) at the atmosphere-ocean
 interface in a coupled `esm`.
 """
-function atmosphere_ocean_salinity_flux(esm::EarthSystemModel)
-    atmosphere_ocean_salinity_flux =
-        net_ocean_salinity_flux(esm) - sea_ice_ocean_salinity_flux(esm)
-    return Field(atmosphere_ocean_salinity_flux)
-end
+atmosphere_ocean_salinity_flux(esm::EarthSystemModel) =
+    net_ocean_salinity_flux(esm) - sea_ice_ocean_salinity_flux(esm)
 
 
 ###########################
@@ -173,9 +155,8 @@ Return the net freshwater mass flux (kg m⁻² s⁻¹) at the ocean's surface in
 """
 function net_ocean_freshwater_flux(esm::EarthSystemModel; reference_salinity = 35)
     ρᵒᶜ = esm.interfaces.ocean_properties.reference_density
-    S₀ = convert(typeof(ρᵒᶜ), reference_salinity)
-    net_ocean_frashwater_flux = - ρᵒᶜ / S₀ * net_ocean_salinity_flux(esm)
-    return Field(net_ocean_frashwater_flux)
+    Sᵒᶜ = convert(typeof(ρᵒᶜ), reference_salinity)
+    return - ρᵒᶜ / Sᵒᶜ * net_ocean_salinity_flux(esm)
 end
 
 """
@@ -186,12 +167,11 @@ in a coupled `esm`.
 """
 function sea_ice_ocean_freshwater_flux(esm::EarthSystemModel; reference_salinity = 35)
     ρᵒᶜ = esm.interfaces.ocean_properties.reference_density
-    S₀ = convert(typeof(ρᵒᶜ), reference_salinity)
-    sea_ice_ocean_freshwater_flux = - ρᵒᶜ / S₀ * sea_ice_ocean_salinity_flux(esm)
-    return Field(sea_ice_ocean_freshwater_flux)
+    Sᵒᶜ = convert(typeof(ρᵒᶜ), reference_salinity)
+    return - ρᵒᶜ / Sᵒᶜ * sea_ice_ocean_salinity_flux(esm)
 end
 
-sea_ice_ocean_freshwater_flux(esm::NoSeaIceOceanInterfaceModel; kwargs...) = ZeroField()
+sea_ice_ocean_freshwater_flux(::NoSeaIceOceanInterfaceModel; kwargs...) = ZeroField()
 
 """
     atmosphere_ocean_freshwater_flux(esm::EarthSystemModel)
@@ -201,7 +181,6 @@ interface in a coupled `esm`.
 """
 function atmosphere_ocean_freshwater_flux(esm::EarthSystemModel; reference_salinity = 35)
     ρᵒᶜ = esm.interfaces.ocean_properties.reference_density
-    S₀ = convert(typeof(ρᵒᶜ), reference_salinity)
-    atmosphere_ocean_freshwater_flux = - ρᵒᶜ / S₀ * atmosphere_ocean_salinity_flux(esm)
-    return Field(atmosphere_ocean_freshwater_flux)
+    Sᵒᶜ = convert(typeof(ρᵒᶜ), reference_salinity)
+    return - ρᵒᶜ / Sᵒᶜ * atmosphere_ocean_salinity_flux(esm)
 end
diff --git a/src/NumericalEarth.jl b/src/NumericalEarth.jl
index a44dc928..cd915b56 100644
--- a/src/NumericalEarth.jl
+++ b/src/NumericalEarth.jl
@@ -94,9 +94,7 @@ export
     default_sea_ice,
     sea_ice_dynamics,
     initialize!,
-    frazil_temperature_flux, net_ocean_temperature_flux, sea_ice_ocean_temperature_flux, atmosphere_ocean_temperature_flux,
-    frazil_heat_flux, net_ocean_heat_flux, sea_ice_ocean_heat_flux, atmosphere_ocean_heat_flux,
-    net_ocean_salinity_flux, sea_ice_ocean_salinity_flux, atmosphere_ocean_salinity_flux,
+    net_ocean_heat_flux, sea_ice_ocean_heat_flux, atmosphere_ocean_heat_flux,
     net_ocean_freshwater_flux, sea_ice_ocean_freshwater_flux, atmosphere_ocean_freshwater_flux,
     meridional_heat_transport,
     location,
diff --git a/test/test_diagnostics_2.jl b/test/test_diagnostics_2.jl
index f37a7150..48bd53e6 100644
--- a/test/test_diagnostics_2.jl
+++ b/test/test_diagnostics_2.jl
@@ -1,5 +1,7 @@
 include("runtests_setup.jl")
 
+using NumericalEarth.Diagnostics: Diagnostics
+
 for arch in test_architectures
     A = typeof(arch)
     @info "Testing InterfaceFluxOutputs [$A]..."
@@ -15,8 +17,8 @@ for arch in test_architectures
     interface_heat_flux_value = 0.3
     sea_ice_ocean_salt_flux_value = 0.9
 
-    model_configurations = ((name = "OceanSeaIceModel", with_sea_ice = true),
-                            (name = "OceanOnlyModel", with_sea_ice = false))
+    model_configurations = ((name = "OceanOnlyModel", with_sea_ice = false),
+                            (name = "OceanSeaIceModel", with_sea_ice = true))
 
     for config in model_configurations
         @testset "InterfaceFluxOutputs on $(config.name) [$A]" begin
@@ -54,21 +56,24 @@ for arch in test_architectures
             cᵒᶜ = esm.interfaces.ocean_properties.heat_capacity
             Sᵒᶜ = 35.0
 
-            frazil_temperature = frazil_temperature_flux(esm)
-            net_ocean_temperature = net_ocean_temperature_flux(esm)
-            sea_ice_ocean_temperature = sea_ice_ocean_temperature_flux(esm)
-            atmosphere_ocean_temperature = atmosphere_ocean_temperature_flux(esm)
-            frazil_heat = frazil_heat_flux(esm)
+            # exported diagnostics
             net_ocean_heat = net_ocean_heat_flux(esm)
             sea_ice_ocean_heat = sea_ice_ocean_heat_flux(esm)
             atmosphere_ocean_heat = atmosphere_ocean_heat_flux(esm)
-            net_ocean_salinity = net_ocean_salinity_flux(esm)
-            sea_ice_ocean_salinity = sea_ice_ocean_salinity_flux(esm)
-            atmosphere_ocean_salinity = atmosphere_ocean_salinity_flux(esm)
             net_ocean_freshwater = net_ocean_freshwater_flux(esm; reference_salinity = Sᵒᶜ)
             sea_ice_ocean_freshwater = sea_ice_ocean_freshwater_flux(esm; reference_salinity = Sᵒᶜ)
             atmosphere_ocean_freshwater = atmosphere_ocean_freshwater_flux(esm; reference_salinity = Sᵒᶜ)
 
+            # internal diagnostics
+            frazil_temperature = Diagnostics.frazil_temperature_flux(esm)
+            net_ocean_temperature = Diagnostics.net_ocean_temperature_flux(esm)
+            sea_ice_ocean_temperature = Diagnostics.sea_ice_ocean_temperature_flux(esm)
+            atmosphere_ocean_temperature = Diagnostics.atmosphere_ocean_temperature_flux(esm)
+            frazil_heat = Diagnostics.frazil_heat_flux(esm)
+            net_ocean_salinity = Diagnostics.net_ocean_salinity_flux(esm)
+            sea_ice_ocean_salinity = Diagnostics.sea_ice_ocean_salinity_flux(esm)
+            atmosphere_ocean_salinity = Diagnostics.atmosphere_ocean_salinity_flux(esm)
+
             for f in (frazil_temperature, net_ocean_temperature, sea_ice_ocean_temperature,
                       atmosphere_ocean_temperature, frazil_heat, net_ocean_heat, sea_ice_ocean_heat,
                       atmosphere_ocean_heat, net_ocean_salinity, sea_ice_ocean_salinity,