Clean up some diagnostics and ensure that diagnostics return AbstractOperations

src/Diagnostics/Diagnostics.jl

@@ -1,11 +1,10 @@
 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
+export meridional_heat_transport
+
 
 using KernelAbstractions: @index, @kernel
 

src/Diagnostics/interface_fluxes.jl

@@ -70,6 +70,8 @@ end
 ###########################
 
 """
+<<<<<<< HEAD
+=======
     frazil_heat_flux(esm::EarthSystemModel)
 
 Return the two-dimensional frazil heat flux (W m⁻²) in a coupled `esm`.
@@ -82,15 +84,17 @@ end
 frazil_heat_flux(esm::NoSeaIceOceanInterfaceModel) = ZeroField()
 
 """
+>>>>>>> main
     net_ocean_heat_flux(esm::EarthSystemModel)
 
 Return the net heat flux (W m⁻²) at the ocean's surface in a coupled `esm`.
 """
 function net_ocean_heat_flux(esm::EarthSystemModel)
+    Jᵀ = esm.ocean.model.tracers.T.boundary_conditions.top.condition # temperature flux
     ρᵒᶜ = 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)
+    frazil_heat_flux = esm.interfaces.sea_ice_ocean_interface.fluxes.frazil_heat
+    return ρᵒᶜ * cᵒᶜ * Jᵀ + frazil_heat_flux
 end
 
 """
@@ -100,9 +104,14 @@ Return the sea ice-ocean heat flux (W m⁻²) at the sea ice-ocean interface
 in a coupled `esm`.
 """
 function sea_ice_ocean_heat_flux(esm::EarthSystemModel)
+<<<<<<< HEAD
+    frazil_heat_flux = esm.interfaces.sea_ice_ocean_interface.fluxes.frazil_heat
+    return esm.interfaces.sea_ice_ocean_interface.fluxes.interface_heat + frazil_heat_flux
+=======
     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)
+>>>>>>> main
 end
 
 sea_ice_ocean_heat_flux(esm::NoSeaIceOceanInterfaceModel) = ZeroField()
@@ -113,6 +122,10 @@ sea_ice_ocean_heat_flux(esm::NoSeaIceOceanInterfaceModel) = ZeroField()
 Return the atmosphere-ocean heat flux (W m⁻²) at the atmosphere-ocean
 interface in a coupled `esm`.
 """
+<<<<<<< HEAD
+atmosphere_ocean_heat_flux(esm::EarthSystemModel) =
+    net_ocean_heat_flux(esm) - sea_ice_ocean_heat_flux(esm)
+=======
 function atmosphere_ocean_heat_flux(esm::EarthSystemModel)
     ρᵒᶜ = esm.interfaces.ocean_properties.reference_density
     cᵒᶜ = esm.interfaces.ocean_properties.heat_capacity
@@ -160,6 +173,7 @@ function atmosphere_ocean_salinity_flux(esm::EarthSystemModel)
         net_ocean_salinity_flux(esm) - sea_ice_ocean_salinity_flux(esm)
     return Field(atmosphere_ocean_salinity_flux)
 end
+>>>>>>> main
 
 
 ###########################
@@ -171,25 +185,17 @@ end
 
 Return the net freshwater mass flux (kg m⁻² s⁻¹) at the ocean's surface in a coupled `esm`.
 """
-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)
-end
+net_ocean_freshwater_flux(esm::EarthSystemModel) =
+    sea_ice_ocean_freshwater_flux(esm) + atmosphere_ocean_freshwater_flux(esm)
 
 """
     sea_ice_ocean_freshwater_flux(esm::EarthSystemModel)
 
 Return the sea ice-ocean freshwater mass flux (kg m⁻² s⁻¹) at the sea ice-ocean interface
 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)
-end
+sea_ice_ocean_freshwater_flux(esm::EarthSystemModel) =
+    esm.interfaces.sea_ice_ocean_interface.fluxes.freshwater_flux
 
 sea_ice_ocean_freshwater_flux(esm::NoSeaIceOceanInterfaceModel; kwargs...) = ZeroField()
 
@@ -199,9 +205,7 @@ sea_ice_ocean_freshwater_flux(esm::NoSeaIceOceanInterfaceModel; kwargs...) = Zer
 Return the atmosphere-ocean freshwater mass flux (kg m⁻² s⁻¹) at the atmosphere-ocean
 interface in a coupled `esm`.
 """
-function atmosphere_ocean_freshwater_flux(esm::EarthSystemModel; reference_salinity = 35)
+function atmosphere_ocean_freshwater_flux(esm::EarthSystemModel)
     ρᵒᶜ = 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)
+    return ρᵒᶜ * esm.interfaces.atmosphere_ocean_interface.fluxes.freshwater_flux
 end

src/EarthSystemModels/InterfaceComputations/component_interfaces.jl

@@ -32,7 +32,8 @@ Container for sea ice-ocean interface data including fluxes, formulation, and in
 Fields
 ======
 
-- `fluxes::J`: `SeaIceOceanFluxes` containing interface_heat, frazil_heat, salt, x_momentum, y_momentum
+- `fluxes::J`: `SeaIceOceanFluxes` containing `interface_heat`, `frazil_heat`, `salt`, `x_momentum`, `y_momentum`,
+                `and freshwater_flux`.
 - `flux_formulation::F`: heat flux formulation (`IceBathHeatFlux` or `ThreeEquationHeatFlux`)
 - `temperature::T`: interface temperature field (ocean surface view or computed field)
 - `salinity::S`: interface salinity field (ocean surface view or computed field)
@@ -205,9 +206,45 @@ function atmosphere_ocean_interface(grid,
                                     velocity_formulation,
                                     specific_humidity_formulation)
 
+<<<<<<< HEAD
+    water_vapor           = Field{Center, Center, Nothing}(grid)
+    latent_heat           = Field{Center, Center, Nothing}(grid)
+    sensible_heat         = Field{Center, Center, Nothing}(grid)
+    x_momentum            = Field{Center, Center, Nothing}(grid)
+    y_momentum            = Field{Center, Center, Nothing}(grid)
+    friction_velocity     = Field{Center, Center, Nothing}(grid)
+    temperature_scale     = Field{Center, Center, Nothing}(grid)
+    water_vapor_scale     = Field{Center, Center, Nothing}(grid)
+    freshwater_flux       = Field{Center, Center, Nothing}(grid)
+    upwelling_longwave    = Field{Center, Center, Nothing}(grid)
+    downwelling_longwave  = Field{Center, Center, Nothing}(grid)
+    downwelling_shortwave = Field{Center, Center, Nothing}(grid)
+
+    ao_fluxes = (; latent_heat,
+                   sensible_heat,
+                   water_vapor,
+                   x_momentum,
+                   y_momentum,
+                   friction_velocity,
+                   temperature_scale,
+                   water_vapor_scale,
+                   freshwater_flux,
+                   upwelling_longwave,
+                   downwelling_longwave,
+                   downwelling_shortwave)
+
+    σ = radiation.stefan_boltzmann_constant
+    αₐₒ = radiation.reflection.ocean
+    ϵₐₒ = radiation.emission.ocean
+    radiation = (σ=σ, α=αₐₒ, ϵ=ϵₐₒ)
+
+    ao_properties = InterfaceProperties(radiation,
+                                        specific_humidity_formulation,
+=======
     ao_fluxes = AtmosphereOceanFluxes(grid)
 
     ao_properties = InterfaceProperties(specific_humidity_formulation,
+>>>>>>> main
                                         temperature_formulation,
                                         velocity_formulation)
 
@@ -281,7 +318,24 @@ Arguments
 - `flux_formulation`: heat flux formulation (`IceBathHeatFlux` or `ThreeEquationHeatFlux`)
 """
 function sea_ice_ocean_interface(grid, sea_ice, ocean, flux_formulation)
+<<<<<<< HEAD
+
+    io_bottom_heat_flux = Field{Center, Center, Nothing}(grid)
+    io_frazil_heat_flux = Field{Center, Center, Nothing}(grid)
+    io_salt_flux = Field{Center, Center, Nothing}(grid)
+    io_freshwater_flux = Field{Center, Center, Nothing}(grid)
+    x_momentum = Field{Face, Center, Nothing}(grid)
+    y_momentum = Field{Center, Face, Nothing}(grid)
+
+    io_fluxes = (interface_heat = io_bottom_heat_flux,
+                 frazil_heat = io_frazil_heat_flux,
+                 salt = io_salt_flux,
+                 freshwater_flux = io_freshwater_flux,
+                 x_momentum = x_momentum,
+                 y_momentum = y_momentum)
+=======
     io_fluxes = SeaIceOceanFluxes(grid)
+>>>>>>> main
 
     # For default flux formulations, interface temperature and salinity point to ocean surface
     Tⁱⁿᵗ = ocean_surface_temperature(ocean)
@@ -291,7 +345,24 @@ function sea_ice_ocean_interface(grid, sea_ice, ocean, flux_formulation)
 end
 
 function sea_ice_ocean_interface(grid, sea_ice, ocean, flux_formulation::ThreeEquationHeatFlux)
+<<<<<<< HEAD
+
+    io_bottom_heat_flux = Field{Center, Center, Nothing}(grid)
+    io_frazil_heat_flux = Field{Center, Center, Nothing}(grid)
+    io_salt_flux = Field{Center, Center, Nothing}(grid)
+    io_freshwater_flux = Field{Center, Center, Nothing}(grid)
+    x_momentum = Field{Face, Center, Nothing}(grid)
+    y_momentum = Field{Center, Face, Nothing}(grid)
+
+    io_fluxes = (interface_heat = io_bottom_heat_flux,
+                 frazil_heat = io_frazil_heat_flux,
+                 salt = io_salt_flux,
+                 freshwater_flux = io_freshwater_flux,
+                 x_momentum = x_momentum,
+                 y_momentum = y_momentum)
+=======
     io_fluxes = SeaIceOceanFluxes(grid)
+>>>>>>> main
 
     # Interface temperature and salinity are computed fields
     Tⁱⁿᵗ = Field{Center, Center, Nothing}(grid)

src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_fluxes.jl

@@ -124,6 +124,7 @@ end
     𝒬ᶠʳᶻ = fluxes.frazil_heat
     𝒬ⁱⁿᵗ = fluxes.interface_heat
     Jˢ = fluxes.salt
+    ΣF = fluxes.freshwater_flux
     τˣ = fluxes.x_momentum
     τʸ = fluxes.y_momentum
     T★ = interface_temperature
@@ -199,8 +200,8 @@ end
     # =============================================
     # Returns interfacial heat flux, melt rate qᵐ, and interface T, S
     𝒬ⁱᵒ, qᵐ, Tᵦ, Sᵦ = compute_interface_heat_flux(flux_formulation,
-                                                   ocean_surface_state, ice_state,
-                                                   liquidus, ocean_properties, ℰ, u★)
+                                                  ocean_surface_state, ice_state,
+                                                  liquidus, ocean_properties, ℰ, u★)
 
     # Store interface values and heat flux
     @inbounds 𝒬ⁱⁿᵗ[i, j, 1] = 𝒬ⁱᵒ
@@ -212,5 +213,7 @@ end
     # Salt flux from melting/freezing:
     # - during ice melt   (qᵐ > 0), fresh meltwater dilutes the ocean
     # - during ice growth (qᶠ < 0), brine rejection adds salt to ocean
-    @inbounds Jˢ[i, j, 1] = (qᵐ + qᶠ) / ρᵒᶜ * (Sᴺ - Sˢⁱ)
+    ΣFio = qᵐ + qᶠ
+    @inbounds Jˢ[i, j, 1] = ΣFio / ρᵒᶜ * (Sᴺ - Sˢⁱ)
+    @inbounds ΣF[i, j, 1] = ΣFio
 end

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,

src/Oceans/assemble_net_ocean_fluxes.jl

@@ -88,14 +88,47 @@ Base.@propagate_inbounds get_land_freshwater_flux(i, j, flux) = flux[i, j, 1]
     # Turbulent contributions to surface heat flux (radiation added later)
     ΣQao = (𝒬ᵀ + 𝒬ᵛ) * (1 - ℵᵢ)
 
+<<<<<<< HEAD
+    # Compute the interior + surface absorbed shortwave radiation
+    ℐₜˢʷ = transmitted_shortwave_radiation(i, j, kᴺ, grid, time, α, ℐꜜˢʷ)
+
+    ℐₐˡʷ *= (1 - ℵᵢ)
+    ℐₜˢʷ *= (1 - ℵᵢ)
+
+    Qss = shortwave_radiative_forcing(i, j, grid, penetrating_radiation, ℐₜˢʷ, ocean_properties)
+
+    # Compute the total heat flux
+    ΣQao = (ℐꜛˡʷ + 𝒬ᵀ + 𝒬ᵛ) * (1 - ℵᵢ) + ℐₐˡʷ + Qss
+
+    # Convert from a mass flux to a volume flux (aka velocity)
+    # by dividing with the ocean reference density.
+    # Also switch the sign, for some reason we are given freshwater flux as positive down.
+=======
     # Freshwater flux to the ocean per unit cell area (volume flux, positive up = leaving ocean):
     # - rain and land runoff reach the ocean everywhere (rain runs through cracks in ice)
     # - snow only reaches the ocean through the open-water fraction (1 - ℵᵢ);
     # - evaporation acts only over the open-water fraction (1 - ℵᵢ)
     # The atmospheric mass-flux convention is positive down; Jᵛ is positive up.
+>>>>>>> main
     ρᵒᶜ⁻¹ = 1 / ocean_properties.reference_density
     ΣFao = - (Jʳⁿ + Jˡⁿ + (1 - ℵᵢ) * Jˢⁿ) * ρᵒᶜ⁻¹ + (1 - ℵᵢ) * Jᵛ * ρᵒᶜ⁻¹
 
+<<<<<<< HEAD
+    # Add the contribution from the turbulent water vapor flux, which has
+    # a different sign convention as the prescribed water mass fluxes (positive upwards)
+    Jᵛᵒᶜ = Jᵛ * ρᵒᶜ⁻¹
+    ΣFao += Jᵛᵒᶜ
+    @inbounds begin
+        # Write radiative components of the heat flux for diagnostic purposes
+        atmos_ocean_fluxes.upwelling_longwave[i, j, 1] = ℐꜛˡʷ
+        atmos_ocean_fluxes.downwelling_longwave[i, j, 1] = - ℐₐˡʷ
+        atmos_ocean_fluxes.downwelling_shortwave[i, j, 1] = - ℐₜˢʷ
+        atmos_ocean_fluxes.freshwater_flux[i, j, 1] = ΣFao
+    end
+
+    # Compute fluxes for u, v, T, and S from momentum, heat, and freshwater fluxes
+=======
+>>>>>>> main
     τˣ = net_ocean_fluxes.u
     τʸ = net_ocean_fluxes.v
     Jᵀ = net_ocean_fluxes.T