use SIC instead of area fraction [skip ci]

Author: juliasloan25

experiments/ClimaEarth/components/ocean/clima_seaice.jl

@@ -262,6 +262,8 @@ end
 Compute the fluxes between the ocean and sea ice, storing them in the `ocean_ice_fluxes`
 fields of the ocean and sea ice simulations.
 
+This function assumes both simulations share the same grid, so no remapping is done.
+
 !!! note
     This function must be called after the turbulent fluxes have been updated in both
     simulations. Here only the contributions from the sea ice/ocean interactions
@@ -275,11 +277,14 @@ function FluxCalculator.ocean_seaice_fluxes!(
     ocean_properties = ocean_sim.ocean_properties
     ice_concentration = Interfacer.get_field(ice_sim, Val(:ice_concentration))
 
+    # Update the sea ice concentration in the ocean simulation
+    ocean_sim.ice_concentration .= ice_concentration
+
     # Compute the fluxes and store them in the both simulations
     ocean_properties = (;
         reference_density = ocean_properties.ocean_reference_density,
         heat_capacity = ocean_properties.ocean_heat_capacity,
-    ) # TODO rename in constructor
+    ) # TODO rename in ocean constructor
     CO.OceanSeaIceModels.InterfaceComputations.compute_sea_ice_ocean_fluxes!(
         ice_sim.ocean_ice_fluxes,
         ocean_sim.ocean,

experiments/ClimaEarth/components/ocean/oceananigans.jl

@@ -10,7 +10,7 @@ using KernelAbstractions: @kernel, @index, @inbounds
 include("climaocean_helpers.jl")
 
 """
-    OceananigansSimulation{SIM, A, OPROP, REMAP}
+    OceananigansSimulation{SIM, A, OPROP, REMAP, SIC}
 
 The ClimaCoupler simulation object used to run with Oceananigans.
 This type is used by the coupler to indicate that this simulation
@@ -21,12 +21,14 @@ It contains the following objects:
 - `area_fraction::A`: A ClimaCore Field representing the surface area fraction of this component model on the exchange grid.
 - `ocean_properties::OPROP`: A NamedTuple of ocean properties and parameters
 - `remapping::REMAP`: Objects needed to remap from the exchange (spectral) grid to Oceananigans spaces.
+- `ice_concentration::SIC`: An Oceananigans Field representing the sea ice concentration on the ocean/sea ice grid.
 """
-struct OceananigansSimulation{SIM, A, OPROP, REMAP} <: Interfacer.OceanModelSimulation
+struct OceananigansSimulation{SIM, A, OPROP, REMAP, SIC} <: Interfacer.OceanModelSimulation
     ocean::SIM
     area_fraction::A
     ocean_properties::OPROP
     remapping::REMAP
+    ice_concentration::SIC
 end
 
 """
@@ -178,8 +180,11 @@ function OceananigansSimulation(
         ocean.output_writers[:diagnostics] = netcdf_writer
     end
 
-    sim = OceananigansSimulation(ocean, area_fraction, ocean_properties, remapping)
-    return sim
+    # Initialize with 0 ice concentration; this will be updated in `resolve_ocean_ice_fractions!`
+    # if the ocean is coupled to a non-prescribed sea ice model.
+    ice_concentration = OC.Field{OC.Center, OC.Center, Nothing}(grid)
+
+    return OceananigansSimulation(ocean, area_fraction, ocean_properties, remapping, ice_concentration)
 end
 
 """
@@ -193,6 +198,9 @@ degrees latitude, and make sure the ocean fraction is 0 there.
 
 The land fraction is expected to be set to 1 at the poles before calling this function,
 and doesn't need to be set again since its fraction is static.
+
+This function also updates the ice concentration field in the ocean simulation
+so that it can be used for weighting flux updates.
 """
 function FieldExchanger.resolve_ocean_ice_fractions!(
     ocean_sim::OceananigansSimulation,
@@ -215,6 +223,9 @@ function FieldExchanger.resolve_ocean_ice_fractions!(
         @. ice_fraction = ifelse.(polar_mask == FT(1), FT(1) - land_fraction, ice_fraction)
         @. ocean_fraction = ifelse.(polar_mask == FT(1), FT(0), ocean_fraction)
     end
+
+    # Update the ice concentration field in the ocean simulation
+    ocean_sim.ice_concentration .= Interfacer.get_field(ice_sim, Val(:ice_concentration))
     return nothing
 end
 
@@ -268,14 +279,7 @@ function FluxCalculator.update_turbulent_fluxes!(sim::OceananigansSimulation, fi
     # TODO clarify where we need to add and where we set fluxes directly
     (; F_lh, F_sh, F_turb_ρτxz, F_turb_ρτyz, F_turb_moisture) = fields
     grid = sim.ocean.model.grid
-
-    # Remap the area fraction from the boundary space to the Oceananigans grid
-    CC.Remapping.interpolate!(
-        sim.remapping.scratch_arr3,
-        sim.remapping.remapper_cc,
-        sim.area_fraction,
-    )
-    area_fraction = sim.remapping.scratch_arr3
+    ice_concentration = sim.ice_concentration
 
     # Remap momentum fluxes onto reduced 2D Center, Center fields using scratch arrays and fields
     CC.Remapping.interpolate!(
@@ -295,6 +299,11 @@ function FluxCalculator.update_turbulent_fluxes!(sim::OceananigansSimulation, fi
     F_turb_ρτxz_cc = sim.remapping.scratch_cc1
     F_turb_ρτyz_cc = sim.remapping.scratch_cc2
 
+    # Weight by (1 - sea ice concentration)
+    # TODO does this work with OC fields?
+    OC.interior(F_turb_ρτxz_cc, :, :, 1) .= OC.interior(F_turb_ρτxz_cc, :, :, 1) .* (1.0 .- ice_concentration)
+    OC.interior(F_turb_ρτyz_cc, :, :, 1) .= OC.interior(F_turb_ρτyz_cc, :, :, 1) .* (1.0 .- ice_concentration)
+
     # Set the momentum flux BCs at the correct locations using the remapped scratch fields
     oc_flux_u = surface_flux(sim.ocean.model.velocities.u)
     oc_flux_v = surface_flux(sim.ocean.model.velocities.v)
@@ -303,7 +312,7 @@ function FluxCalculator.update_turbulent_fluxes!(sim::OceananigansSimulation, fi
         grid,
         F_turb_ρτxz_cc,
         F_turb_ρτyz_cc,
-    ) # TODO multiply by area_fraction?
+    )
 
     (; ocean_reference_density, ocean_heat_capacity, ocean_fresh_water_density) =
         sim.ocean_properties
@@ -322,7 +331,7 @@ function FluxCalculator.update_turbulent_fluxes!(sim::OceananigansSimulation, fi
     oc_flux_T = surface_flux(sim.ocean.model.tracers.T)
     OC.interior(oc_flux_T, :, :, 1) .=
         OC.interior(oc_flux_T, :, :, 1) .+
-        area_fraction .* (remapped_F_lh .+ remapped_F_sh) ./
+        (1.0 .- ice_concetration) .* (remapped_F_lh .+ remapped_F_sh) ./
         (ocean_reference_density * ocean_heat_capacity)
 
     # Add the part of the salinity flux that comes from the moisture flux, we also need to
@@ -337,7 +346,7 @@ function FluxCalculator.update_turbulent_fluxes!(sim::OceananigansSimulation, fi
     surface_salinity = OC.interior(sim.ocean.model.tracers.S, :, :, 1)
     OC.interior(oc_flux_S, :, :, 1) .=
         OC.interior(oc_flux_S, :, :, 1) .-
-        area_fraction .* surface_salinity .* moisture_fresh_water_flux
+        (1.0 .- ice_concetration) .* surface_salinity .* moisture_fresh_water_flux
     return nothing
 end
 
@@ -369,14 +378,7 @@ so a sign change is needed when we convert from precipitation to salinity flux.
 function FieldExchanger.update_sim!(sim::OceananigansSimulation, csf)
     (; ocean_reference_density, ocean_heat_capacity, ocean_fresh_water_density) =
         sim.ocean_properties
-    # TODO use SIC instead?
-    # Remap the area fraction from the boundary space to the Oceananigans grid
-    CC.Remapping.interpolate!(
-        sim.remapping.scratch_arr3,
-        sim.remapping.remapper_cc,
-        sim.area_fraction,
-    )
-    area_fraction = sim.remapping.scratch_arr3
+    ice_concentration = sim.ice_concentration
 
     # Remap radiative flux onto scratch array; rename for clarity
     CC.Remapping.interpolate!(
@@ -401,7 +403,7 @@ function FieldExchanger.update_sim!(sim::OceananigansSimulation, csf)
     α = Interfacer.get_field(sim, Val(:surface_direct_albedo)) # scalar
     ϵ = Interfacer.get_field(sim, Val(:emissivity)) # scalar
     OC.interior(oc_flux_T, :, :, 1) .=
-        area_fraction .* (
+        (1.0 .- ice_concetration) .* (
             -(1 - α) .* remapped_SW_d .-
             ϵ * (
                 remapped_LW_d .-
@@ -413,20 +415,21 @@ function FieldExchanger.update_sim!(sim::OceananigansSimulation, csf)
     CC.Remapping.interpolate!(
         sim.remapping.scratch_arr1,
         sim.remapping.remapper_cc,
-        sim.area_fraction .* csf.P_liq,
+        csf.P_liq,
     )
     CC.Remapping.interpolate!(
         sim.remapping.scratch_arr2,
         sim.remapping.remapper_cc,
-        sim.area_fraction .* csf.P_snow,
+        csf.P_snow,
     )
     remapped_P_liq = sim.remapping.scratch_arr1
     remapped_P_snow = sim.remapping.scratch_arr2
 
     # Virtual salt flux
+    # TODO why is this split on multiple lines?
     oc_flux_S = surface_flux(sim.ocean.model.tracers.S)
     precipitating_fresh_water_flux =
-        area_fraction .* (remapped_P_liq .+ remapped_P_snow) ./ ocean_fresh_water_density
+        (1.0 .- ice_concetration) .* (remapped_P_liq .+ remapped_P_snow) ./ ocean_fresh_water_density
     surface_salinity_flux =
         OC.interior(sim.ocean.model.tracers.S, :, :, 1) .* precipitating_fresh_water_flux
     OC.interior(oc_flux_S, :, :, 1) .=