Support setting sea ice thermodynamics to nothing

src/EarthSystemModels/InterfaceComputations/component_interfaces.jl

@@ -254,6 +254,28 @@ end
 ##### Sea Ice-Ocean Interface
 #####
 
+struct DefaultSeaIceOceanSalinityFlux end
+const default_sea_ice_ocean_salinity_flux = DefaultSeaIceOceanSalinityFlux()
+
+with_sea_ice_ocean_salinity_flux(::Nothing, salinity_flux) = nothing
+
+function with_sea_ice_ocean_salinity_flux(flux::IceBathHeatFlux, salinity_flux)
+    return IceBathHeatFlux(flux.heat_transfer_coefficient,
+                           flux.friction_velocity;
+                           sea_ice_heat_flux = flux.sea_ice_heat_flux,
+                           salinity_flux)
+end
+
+function with_sea_ice_ocean_salinity_flux(flux::ThreeEquationHeatFlux, salinity_flux)
+    return ThreeEquationHeatFlux(flux.conductive_flux,
+                                 flux.internal_temperature,
+                                 flux.heat_transfer_coefficient,
+                                 flux.salt_transfer_coefficient,
+                                 flux.friction_velocity;
+                                 sea_ice_heat_flux = flux.sea_ice_heat_flux,
+                                 salinity_flux)
+end
+
 sea_ice_ocean_interface(grid, ::Nothing, ocean,     flux_formulation; kwargs...) = nothing
 sea_ice_ocean_interface(grid, ::Nothing, ::Nothing, flux_formulation; kwargs...) = nothing
 sea_ice_ocean_interface(grid, sea_ice,   ::Nothing, flux_formulation; kwargs...) = nothing
@@ -280,7 +302,9 @@ Arguments
 - `ocean`: ocean simulation
 - `flux_formulation`: heat flux formulation (`IceBathHeatFlux` or `ThreeEquationHeatFlux`)
 """
-function sea_ice_ocean_interface(grid, sea_ice, ocean, flux_formulation)
+function sea_ice_ocean_interface(grid, sea_ice, ocean, flux_formulation;
+                                 sea_ice_ocean_salinity_flux = default_sea_ice_ocean_salinity_flux)
+    flux_formulation = with_sea_ice_ocean_salinity_flux(flux_formulation, sea_ice_ocean_salinity_flux)
     io_fluxes = SeaIceOceanFluxes(grid)
 
     # For default flux formulations, interface temperature and salinity point to ocean surface
@@ -290,7 +314,9 @@ function sea_ice_ocean_interface(grid, sea_ice, ocean, flux_formulation)
     return SeaIceOceanInterface(io_fluxes, flux_formulation, Tⁱⁿᵗ, Sⁱⁿᵗ)
 end
 
-function sea_ice_ocean_interface(grid, sea_ice, ocean, flux_formulation::ThreeEquationHeatFlux)
+function sea_ice_ocean_interface(grid, sea_ice, ocean, flux_formulation::ThreeEquationHeatFlux;
+                                 sea_ice_ocean_salinity_flux = default_sea_ice_ocean_salinity_flux)
+    flux_formulation = with_sea_ice_ocean_salinity_flux(flux_formulation, sea_ice_ocean_salinity_flux)
     io_fluxes = SeaIceOceanFluxes(grid)
 
     # Interface temperature and salinity are computed fields
@@ -321,12 +347,14 @@ Construct component interfaces for atmosphere-ocean-sea ice coupling.
 Keyword Arguments
 =================
 
-- `sea_ice_ocean_heat_flux`: formulation for sea ice-ocean heat flux. Options are:
+- `sea_ice_ocean_thermodynamics`: formulation for sea ice-ocean thermodynamics. Options are:
   - `IceBathHeatFlux()`: bulk heat flux with interface at freezing point
   - `ThreeEquationHeatFlux()`: coupled heat/salt/freezing point system (default)
+  - `nothing`: omit sea ice-ocean thermodynamic exchange
+- `sea_ice_ocean_heat_flux`: alias for `sea_ice_ocean_thermodynamics`
+- `sea_ice_ocean_salinity_flux`: formulation for sea ice-ocean salinity flux. Set to `nothing` to omit it.
 
-- `radiation`: radiation component. Default: `Radiation()`.
-+ `radiation`: radiation component. Default: `nothing`.
+- `radiation`: radiation component. Default: `nothing`.
 - `freshwater_density`: reference density of freshwater. Default: `default_freshwater_density`.
 - `atmosphere_ocean_fluxes`: flux formulation for atmosphere-ocean interface. Default: `SimilarityTheoryFluxes()`.
 - `atmosphere_sea_ice_fluxes`: flux formulation for atmosphere-sea ice interface. Default: `SimilarityTheoryFluxes()`.
@@ -348,7 +376,9 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing;
                              freshwater_density = default_freshwater_density,
                              atmosphere_ocean_fluxes = SimilarityTheoryFluxes(eltype(exchange_grid)),
                              atmosphere_sea_ice_fluxes = atmosphere_sea_ice_similarity_theory(eltype(exchange_grid)),
-                             sea_ice_ocean_heat_flux = ThreeEquationHeatFlux(sea_ice),
+                             sea_ice_ocean_thermodynamics = ThreeEquationHeatFlux(sea_ice),
+                             sea_ice_ocean_heat_flux = sea_ice_ocean_thermodynamics,
+                             sea_ice_ocean_salinity_flux = default_sea_ice_ocean_salinity_flux,
                              atmosphere_ocean_interface_temperature = BulkTemperature(),
                              atmosphere_ocean_velocity_difference = RelativeVelocity(),
                              atmosphere_ocean_interface_specific_humidity = default_ao_specific_humidity(ocean),
@@ -399,7 +429,8 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing;
                                               atmosphere_ocean_velocity_difference,
                                               atmosphere_ocean_interface_specific_humidity)
 
-    io_interface = sea_ice_ocean_interface(exchange_grid, sea_ice, ocean, sea_ice_ocean_heat_flux)
+    io_interface = sea_ice_ocean_interface(exchange_grid, sea_ice, ocean, sea_ice_ocean_heat_flux;
+                                           sea_ice_ocean_salinity_flux)
 
     ai_interface = atmosphere_sea_ice_interface(exchange_grid,
                                                 atmosphere,

src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_fluxes.jl

@@ -57,6 +57,7 @@ function compute_sea_ice_ocean_fluxes!(interface, ocean, sea_ice, ocean_properti
     if !isnothing(dynamics)
         kernel_parameters = interface_kernel_parameters(grid)
         τₛ = dynamics.external_momentum_stresses.bottom
+
         launch!(arch, grid, kernel_parameters, _compute_sea_ice_ocean_stress!,
                 fluxes, grid, clock, hˢⁱ, ℵ, uˢⁱ, vˢⁱ, τₛ)
     else
@@ -98,43 +99,13 @@ end
     end
 end
 
-@kernel function _compute_sea_ice_ocean_fluxes!(flux_formulation,
-                                                fluxes,
-                                                interface_temperature,
-                                                interface_salinity,
-                                                grid,
-                                                clock,
-                                                ice_thickness,
-                                                ice_consolidation_thickness,
-                                                ice_concentration,
-                                                ice_salinity,
-                                                ocean_temperature,
-                                                ocean_salinity,
-                                                sea_ice_u_velocity,
-                                                sea_ice_v_velocity,
-                                                sea_ice_ocean_stresses,
-                                                liquidus,
-                                                ocean_properties,
-                                                latent_heat,
-                                                Δt)
+@inline compute_frazil_heat_flux!(::Nothing, i, j, grid, Tᵒᶜ, Sᵒᶜ,
+                                  liquidus, ρᵒᶜ, cᵒᶜ, Δt) = zero(grid)
 
-    i, j = @index(Global, NTuple)
+@inline function compute_frazil_heat_flux!(sea_ice_heat_flux, i, j, grid, Tᵒᶜ, Sᵒᶜ,
+                                           liquidus, ρᵒᶜ, cᵒᶜ, Δt)
 
     Nz = size(grid, 3)
-    𝒬ᶠʳᶻ = fluxes.frazil_heat
-    𝒬ⁱⁿᵗ = fluxes.interface_heat
-    Jˢ = fluxes.salt
-    τˣ = fluxes.x_momentum
-    τʸ = fluxes.y_momentum
-    T★ = interface_temperature
-    S★ = interface_salinity
-    Tᵒᶜ = ocean_temperature
-    Sᵒᶜ = ocean_salinity
-    hc = ice_consolidation_thickness
-    ℰ  = latent_heat
-
-    ρᵒᶜ = ocean_properties.reference_density
-    cᵒᶜ = ocean_properties.heat_capacity
 
     # =============================================
     # Part 1: Frazil ice formation (all formulations)
@@ -169,48 +140,166 @@ end
         δ𝒬ᶠʳᶻ -= δE * Δz / Δt
     end
 
-    # Store frazil heat flux
-    @inbounds 𝒬ᶠʳᶻ[i, j, 1] = δ𝒬ᶠʳᶻ
+    return δ𝒬ᶠʳᶻ
+end
 
-    # Freezing rate
-    qᶠ = δ𝒬ᶠʳᶻ / ℰ
+@inline function compute_sea_ice_heat_flux(::Nothing,
+                                           flux_formulation,
+                                           ocean_surface_state,
+                                           ice_state,
+                                           liquidus,
+                                           ocean_properties,
+                                           latent_heat,
+                                           u★,
+                                           grid)
+
+    return zero(grid), zero(grid), zero(grid), zero(grid)
+end
 
-    @inbounds begin
-        Tᴺ  = Tᵒᶜ[i, j, Nz]
-        Sᴺ  = Sᵒᶜ[i, j, Nz]
-        Sˢⁱ = ice_salinity[i, j, 1]
-        hˢⁱ = ice_thickness[i, j, 1]
-        ℵᵢ  = ice_concentration[i, j, 1]
-        hc  = ice_consolidation_thickness[i, j, 1]
-    end
+@inline function compute_sea_ice_heat_flux(sea_ice_heat_flux,
+                                           flux_formulation,
+                                           ocean_surface_state,
+                                           ice_state,
+                                           liquidus,
+                                           ocean_properties,
+                                           latent_heat,
+                                           u★,
+                                           grid)
+
+    return compute_interface_heat_flux(flux_formulation,
+                                       ocean_surface_state,
+                                       ice_state,
+                                       liquidus,
+                                       ocean_properties,
+                                       latent_heat,
+                                       u★)
+end
+
+@inline compute_salinity_flux(::Nothing, qᵐ, qᶠ, ρᵒᶜ, Sᴺ, Sˢⁱ, grid) = zero(grid)
+
+@inline compute_salinity_flux(salinity_flux, qᵐ, qᶠ, ρᵒᶜ, Sᴺ, Sˢⁱ, grid) =
+    (qᵐ + qᶠ) / ρᵒᶜ * (Sᴺ - Sˢⁱ)
+
+@kernel function _compute_sea_ice_ocean_fluxes!(flux_formulation,
+                                                fluxes,
+                                                interface_temperature,
+                                                interface_salinity,
+                                                grid,
+                                                clock,
+                                                ice_thickness,
+                                                ice_consolidation_thickness,
+                                                ice_concentration,
+                                                ice_salinity,
+                                                ocean_temperature,
+                                                ocean_salinity,
+                                                sea_ice_u_velocity,
+                                                sea_ice_v_velocity,
+                                                sea_ice_ocean_stresses,
+                                                liquidus,
+                                                ocean_properties,
+                                                latent_heat,
+                                                Δt)
 
-    # Extract internal temperature (for ConductiveFluxTEF, zero otherwise)
-    Tˢⁱ = extract_internal_temperature(flux_formulation, i, j)
+    i, j = @index(Global, NTuple)
 
-    # Package states
-    ocean_surface_state = (; T = Tᴺ, S = Sᴺ)
-    ice_state = (; S = Sˢⁱ, h = hˢⁱ, hc = hc, ℵ = ℵᵢ, T = Tˢⁱ)
+    𝒬ᶠʳᶻ = fluxes.frazil_heat
+    𝒬ⁱⁿᵗ = fluxes.interface_heat
+    Jˢ = fluxes.salt
+    T★ = interface_temperature
+    S★ = interface_salinity
 
-    # Compute friction velocity
-    u★ = get_friction_velocity(flux_formulation.friction_velocity, i, j, grid, τˣ, τʸ, ρᵒᶜ)
+    if isnothing(flux_formulation)
 
-    # =============================================
-    # Part 3: Interface heat flux (formulation-specific)
-    # =============================================
-    # 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★)
+        @inbounds begin
+            𝒬ᶠʳᶻ[i, j, 1] = zero(grid)
+            𝒬ⁱⁿᵗ[i, j, 1] = zero(grid)
+            Jˢ[i, j, 1] = zero(grid)
+            T★[i, j, 1] = zero(grid)
+            S★[i, j, 1] = zero(grid)
+        end
 
-    # Store interface values and heat flux
-    @inbounds 𝒬ⁱⁿᵗ[i, j, 1] = 𝒬ⁱᵒ
-    store_interface_state!(flux_formulation, T★, S★, i, j, Tᵦ, Sᵦ)
+    else
 
-    # =============================================
-    # Part 4: Salt flux
-    # =============================================
-    # 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ˢⁱ)
-end
+        Nz = size(grid, 3)
+
+        τˣ = fluxes.x_momentum
+        τʸ = fluxes.y_momentum
+        Tᵒᶜ = ocean_temperature
+        Sᵒᶜ = ocean_salinity
+        ℰ = latent_heat
+
+        ρᵒᶜ = ocean_properties.reference_density
+        cᵒᶜ = ocean_properties.heat_capacity
+
+        sea_ice_heat_flux = flux_formulation.sea_ice_heat_flux
+        salinity_flux = flux_formulation.salinity_flux
+
+        # Store frazil heat flux
+        δ𝒬ᶠʳᶻ = compute_frazil_heat_flux!(sea_ice_heat_flux,
+                                           i, j, grid, Tᵒᶜ, Sᵒᶜ,
+                                           liquidus, ρᵒᶜ, cᵒᶜ, Δt)
+
+        @inbounds 𝒬ᶠʳᶻ[i, j, 1] = δ𝒬ᶠʳᶻ
+
+        # Freezing rate
+        qᶠ = δ𝒬ᶠʳᶻ / ℰ
+
+        @inbounds begin
+            Tᴺ  = Tᵒᶜ[i, j, Nz]
+            Sᴺ  = Sᵒᶜ[i, j, Nz]
+            Sˢⁱ = ice_salinity[i, j, 1]
+            hˢⁱ = ice_thickness[i, j, 1]
+            ℵᵢ  = ice_concentration[i, j, 1]
+            hc  = ice_consolidation_thickness[i, j, 1]
+        end
+
+        # Extract internal temperature (for ConductiveFluxTEF, zero otherwise)
+        Tˢⁱ = extract_internal_temperature(flux_formulation, i, j)
+
+        # Package states
+        ocean_surface_state = (; T = Tᴺ, S = Sᴺ)
+        ice_state = (; S = Sˢⁱ, h = hˢⁱ, hc = hc, ℵ = ℵᵢ, T = Tˢⁱ)
+
+        # Compute friction velocity
+        u★ = get_friction_velocity(flux_formulation.friction_velocity,
+                                    i, j, grid, τˣ, τʸ, ρᵒᶜ)
+
+        # =============================================
+        # Part 3: Interface heat flux (formulation-specific)
+        # =============================================
+        # Returns interfacial heat flux, melt rate qᵐ, and interface T, S
+
+        𝒬ⁱᵒ, qᵐ, Tᵦ, Sᵦ = compute_sea_ice_heat_flux(sea_ice_heat_flux,
+                                                       flux_formulation,
+                                                       ocean_surface_state,
+                                                       ice_state,
+                                                       liquidus,
+                                                       ocean_properties,
+                                                       ℰ,
+                                                       u★,
+                                                       grid)
+
+        # Store interface values and heat flux
+        @inbounds 𝒬ⁱⁿᵗ[i, j, 1] = 𝒬ⁱᵒ
+
+        if isnothing(sea_ice_heat_flux)
+            @inbounds begin
+                T★[i, j, 1] = zero(grid)
+                S★[i, j, 1] = zero(grid)
+            end
+        else
+            store_interface_state!(flux_formulation, T★, S★, i, j, Tᵦ, Sᵦ)
+        end
+
+        # =============================================
+        # Part 4: Salt flux
+        # =============================================
+        # 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] = compute_salinity_flux(salinity_flux,
+                                                       qᵐ, qᶠ, ρᵒᶜ,
+                                                       Sᴺ, Sˢⁱ, grid)
+    end
+end