the code runs

Author: juliasloan25

docs/src/fieldexchanger.md

@@ -28,15 +28,17 @@ the atmosphere and each surface model.
 ## FieldExchanger API
 
 ```@docs
+    ClimaCoupler.FieldExchanger.exchange!
     ClimaCoupler.FieldExchanger.update_sim!
     ClimaCoupler.FieldExchanger.step_model_sims!
     ClimaCoupler.FieldExchanger.update_surface_fractions!
-    ClimaCoupler.FieldExchanger.exchange!
     ClimaCoupler.FieldExchanger.set_caches!
 ```
 
 ## FieldExchanger Internal Functions
 
 ```@docs
     ClimaCoupler.FieldExchanger.combine_surfaces!
+    ClimaCoupler.FieldExchanger.resolve_ocean_ice_fractions!
+    ClimaCoupler.FieldExchanger.import_atmos_fields!
 ```

experiments/ClimaEarth/Manifest-v1.11.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.11.6"
 manifest_format = "2.0"
-project_hash = "bf902fbe4cdfe69d6bdc0c99b12a9c00c5aa28ee"
+project_hash = "205aaadcd07952f66e0fa3e53c74e23c6fca069d"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "27cecae79e5cc9935255f90c53bb831cc3c870d7"

experiments/ClimaEarth/components/ocean/clima_seaice.jl

@@ -1,4 +1,5 @@
 import Oceananigans as OC
+import ClimaSeaIce as CSI
 import ClimaOcean as CO
 import ClimaCoupler: Checkpointer, FieldExchanger, FluxCalculator, Interfacer, Utilities
 import ClimaComms
@@ -51,12 +52,12 @@ can be found in the documentation for `ClimaOcean.ocean_simulation`.
 """
 function ClimaSeaIceSimulation(land_fraction, ocean; output_dir)
     # Initialize the sea ice with the same grid as the ocean
-    grid = ocean.ocean.model.grid # TODO can't use lat/lon grid at poles for ice
-    arch = grid.architecture
+    grid = ocean.ocean.model.grid # TODO can't use lat/lon grid at poles for ice, need to fill with bucket
+    arch = OC.Architectures.architecture(grid)
     advection = ocean.ocean.model.advection.T
-    top_heat_boundary_condition = CO.MeltingConstrainedFluxBalance()
+    top_heat_boundary_condition = CSI.MeltingConstrainedFluxBalance()
 
-    # TODO use branch ss-js/top-heat-bc for this constructor
+    # TODO use ClimaOcean 0.8.6
     ice = CO.sea_ice_simulation(grid, ocean.ocean; advection, top_heat_boundary_condition)
 
     melting_speed = 1e-4
@@ -68,7 +69,7 @@ function ClimaSeaIceSimulation(land_fraction, ocean; output_dir)
     if arch isa OC.CPU || pkgversion(OC) >= v"0.96.22"
         # Save all tracers and velocities to a NetCDF file at daily frequency
         outputs = OC.prognostic_fields(ice.model)
-        jld_writer = OC.JLDWriter(
+        jld_writer = OC.JLD2Writer(
             ice.model,
             outputs;
             schedule = OC.TimeInterval(86400), # Daily output
@@ -80,7 +81,19 @@ function ClimaSeaIceSimulation(land_fraction, ocean; output_dir)
     end
 
     # Allocate space for the sea ice-ocean (io) fluxes
-    ocean_ice_fluxes = ocean.ocean_ice_fluxes
+    io_bottom_heat_flux = OC.Field{OC.Center, OC.Center, Nothing}(grid)
+    io_frazil_heat_flux = OC.Field{OC.Center, OC.Center, Nothing}(grid)
+    io_salt_flux = OC.Field{OC.Center, OC.Center, Nothing}(grid)
+    x_momentum = OC.Field{OC.Face, OC.Center, Nothing}(grid)
+    y_momentum = OC.Field{OC.Center, OC.Face, Nothing}(grid)
+
+    ocean_ice_fluxes = (
+        interface_heat = io_bottom_heat_flux,
+        frazil_heat = io_frazil_heat_flux,
+        salt = io_salt_flux,
+        x_momentum = x_momentum,
+        y_momentum = y_momentum,
+    )
 
     # Get the initial area fraction from the fractional ice concentration
     boundary_space = axes(ocean.area_fraction)
@@ -110,7 +123,7 @@ Interfacer.step!(sim::ClimaSeaIceSimulation, t) =
     OC.time_step!(sim.ice, float(t) - sim.ice.model.clock.time)
 
 Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:area_fraction}) = sim.area_fraction
-Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val(:ice_concentration)) =
+Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:ice_concentration}) =
     sim.ice.model.ice_concentration
 
 # At the moment, we return always Float32. This is because we always want to run
@@ -122,6 +135,7 @@ Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:roughness_buoyancy}) =
 Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:roughness_momentum}) =
     Float32(5.8e-5)
 Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:beta}) = Float32(1)
+Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:emissivity}) = Float32(1)
 Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:surface_direct_albedo}) =
     Float32(0.7)
 Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:surface_diffuse_albedo}) =
@@ -130,7 +144,7 @@ Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:surface_diffuse_albedo})
 # Approximate the sea ice surface temperature as the temperature computed from the
 #  fluxes at the previous timestep.
 Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:surface_temperature}) =
-    273.15 .+ OC.interior(sim.ice.model.ice_thermodynamics.top_surface_temperature, :, :, 1)
+    273.15 + sim.ice.model.ice_thermodynamics.top_surface_temperature
 
 """
     FluxCalculator.update_turbulent_fluxes!(sim::ClimaSeaIceSimulation, fields)
@@ -176,8 +190,8 @@ function FluxCalculator.update_turbulent_fluxes!(sim::ClimaSeaIceSimulation, fie
 
     # Set the momentum flux BCs at the correct locations using the remapped scratch fields
     # Note that this requires the sea ice model to always be run with dynamics turned on
-    si_flux_u = sim.ice.model.dynamics.external_stresses.top.u
-    si_flux_v = sim.ice.model.dynamics.external_stresses.top.v
+    si_flux_u = sim.ice.model.dynamics.external_momentum_stresses.top.u
+    si_flux_v = sim.ice.model.dynamics.external_momentum_stresses.top.v
     set_from_extrinsic_vectors!(
         (; u = si_flux_u, v = si_flux_v),
         grid,
@@ -247,23 +261,32 @@ 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.
+
+!!! 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
+    are added to the fluxes.
 """
 function FluxCalculator.ocean_seaice_fluxes!(
     ocean_sim::OceananigansSimulation,
     ice_sim::ClimaSeaIceSimulation,
 )
     melting_speed = ice_sim.melting_speed
     ocean_properties = ocean_sim.ocean_properties
+    ice_concentration = Interfacer.get_field(ice_sim, Val(:ice_concentration))
 
     # Compute the fluxes and store them in the both simulations
-    CO.compute_sea_ice_ocean_fluxes!(
+    ocean_properties = (;
+        reference_density = ocean_properties.ocean_reference_density,
+        heat_capacity = ocean_properties.ocean_heat_capacity,
+    ) # TODO rename in constructor
+    CO.OceanSeaIceModels.InterfaceComputations.compute_sea_ice_ocean_fluxes!(
         ice_sim.ocean_ice_fluxes,
         ocean_sim.ocean,
         ice_sim.ice,
         melting_speed,
         ocean_properties,
     )
-    ocean_sim.ocean_ice_fluxes = ice_sim.ocean_ice_fluxes
 
     ## Update the internals of the sea ice model
     # Set the bottom heat flux to the sum of the frazil and interface heat fluxes
@@ -274,21 +297,62 @@ function FluxCalculator.ocean_seaice_fluxes!(
     bottom_heat_flux .= Qf .+ Qi
 
     ## Update the internals of the ocean model
-    ρₒ⁻¹ = 1 / ocean_sim.ocean_properties.reference_density
-    cₒ = ocean_sim.ocean_properties.heat_capacity
+    ρₒ⁻¹ = 1 / ocean_sim.ocean_properties.ocean_reference_density
+    cₒ = ocean_sim.ocean_properties.ocean_heat_capacity
 
-    Jᵀio = Qio * ρₒ⁻¹ / cₒ
-    Jˢio = ice_sim.ocean_ice_fluxes.salt[i, j, 1] * ℵᵢ
+    # Compute fluxes for u, v, T, and S from momentum, heat, and freshwater fluxes
+    oc_flux_u = surface_flux(ocean_sim.ocean.model.velocities.u)
+    oc_flux_v = surface_flux(ocean_sim.ocean.model.velocities.v)
 
-    # ℑxᶠᵃᵃ: interpolate faces to centers
-    τxio = ρτxio[i, j, 1] * ρₒ⁻¹ * ℑxᶠᵃᵃ(i, j, 1, grid, ℵ)
-    τyio = ρτyio[i, j, 1] * ρₒ⁻¹ * ℑyᵃᶠᵃ(i, j, 1, grid, ℵ)
+    ρτxio = ice_sim.ocean_ice_fluxes.x_momentum # sea_ice - ocean zonal momentum flux
+    ρτyio = ice_sim.ocean_ice_fluxes.y_momentum # sea_ice - ocean meridional momentum flux
+
+    # Update the momentum flux contributions from ocean/sea ice fluxes
+    grid = ocean_sim.ocean.model.grid
+    arch = OC.Architectures.architecture(grid)
+    OC.Utils.launch!(
+        arch,
+        grid,
+        :xy,
+        _add_ocean_ice_stress!,
+        oc_flux_u,
+        oc_flux_v,
+        grid,
+        ρτxio,
+        ρτyio,
+        ρₒ⁻¹,
+        ice_concentration,
+    )
+
+    oc_flux_T = surface_flux(ocean_sim.ocean.model.tracers.T)
+    OC.interior(oc_flux_T, :, :, 1) .+=
+        OC.interior(ice_concentration, :, :, 1) .* OC.interior(Qi, :, :, 1) .* ρₒ⁻¹ ./ cₒ
 
-    # TODO finish this
+    oc_flux_S = surface_flux(ocean_sim.ocean.model.tracers.S)
+    OC.interior(oc_flux_S, :, :, 1) .+=
+        OC.interior(ice_concentration, :, :, 1) .*
+        OC.interior(ice_sim.ocean_ice_fluxes.salt, :, :, 1)
 
     return nothing
 end
 
+@kernel function _add_ocean_ice_stress!(
+    oc_flux_u,
+    oc_flux_v,
+    grid,
+    ρτxio,
+    ρτyio,
+    ρₒ⁻¹,
+    ice_concentration,
+)
+    i, j = @index(Global, NTuple)
+
+    # ℑxᶠᵃᵃ: interpolate faces to centers
+    oc_flux_u +=
+        ρτxio[i, j, 1] * ρₒ⁻¹ * OC.Operators.ℑxᶠᵃᵃ(i, j, 1, grid, ice_concentration)
+    oc_flux_v +=
+        ρτyio[i, j, 1] * ρₒ⁻¹ * OC.Operators.ℑyᵃᶠᵃ(i, j, 1, grid, ice_concentration)
+end
 
 """
     get_model_prog_state(sim::ClimaSeaIceSimulation)

experiments/ClimaEarth/components/ocean/climaocean_helpers.jl

@@ -71,10 +71,10 @@ fields, rotate them onto the target grid and remap to `Face, Center` and
 `Center, Face` fields, respectively.
 """
 function set_from_extrinsic_vectors!(vectors, grid, u_cc, v_cc)
-    arch = grid.architecture
+    arch = OC.Architectures.architecture(grid)
 
     # Rotate vectors onto the grid
-    OC.Utils.launch!(arch, grid, :xy, _rotate_velocities!, u_cc, v_cc, grid)
+    OC.Utils.launch!(arch, grid, :xy, _rotate_vectors!, u_cc, v_cc, grid)
 
     # Fill halo regions with the rotated vectors so we can use them to interpolate
     OC.fill_halo_regions!(u_cc)
@@ -85,7 +85,7 @@ function set_from_extrinsic_vectors!(vectors, grid, u_cc, v_cc)
         arch,
         grid,
         :xy,
-        _interpolate_velocities!,
+        _interpolate_vectors!,
         vectors.u,
         vectors.v,
         grid,
@@ -96,33 +96,33 @@ function set_from_extrinsic_vectors!(vectors, grid, u_cc, v_cc)
 end
 
 """
-    _rotate_velocities!(u, v, grid)
+    _rotate_vectors!(τx, τy, grid)
 
 Rotate the velocities from the extrinsic coordinate system to the intrinsic
 coordinate system.
 """
-@kernel function _rotate_velocities!(u, v, grid)
+@kernel function _rotate_vectors!(τx, τy, grid)
     # Use `k = 1` to index into the reduced Fields
     i, j = @index(Global, NTuple)
     # Rotate u, v from extrinsic to intrinsic coordinate system
-    ur, vr = OC.Operators.intrinsic_vector(i, j, 1, grid, u, v)
+    τxr, τyr = OC.Operators.intrinsic_vector(i, j, 1, grid, τx, τy)
     @inbounds begin
-        u[i, j, 1] = ur
-        v[i, j, 1] = vr
+        τx[i, j, 1] = τxr
+        τy[i, j, 1] = τyr
     end
 end
 
 """
-    _interpolate_velocities!(u, v, grid, u_cc, v_cc)
+    _interpolate_vectors!(τx, τy, grid, τx_cc, τy_cc)
 
-Interpolate the input velocities `u_cc` and `v_cc`, which are Center/Center
+Interpolate the input fluxes `τx_cc` and `τy_cc`, which are Center/Center
 Fields to Face/Center and Center/Face coordinates, respectively.
 """
-@kernel function _interpolate_velocities!(u, v, grid, u_cc, v_cc)
+@kernel function _interpolate_vectors!(τx, τy, grid, τx_cc, τy_cc)
     # Use `k = 1` to index into the reduced Fields
     i, j = @index(Global, NTuple)
     @inbounds begin
-        u[i, j, 1] = OC.Operators.ℑxyᶠᶜᵃ(i, j, 1, grid, u_cc)
-        v[i, j, 1] = OC.Operators.ℑxyᶜᶠᵃ(i, j, 1, grid, v_cc)
+        τx[i, j, 1] = OC.Operators.ℑxᶠᵃᵃ(i, j, 1, grid, τx_cc)
+        τy[i, j, 1] = OC.Operators.ℑyᵃᶠᵃ(i, j, 1, grid, τy_cc)
     end
 end

experiments/ClimaEarth/components/ocean/oceananigans.jl

@@ -21,15 +21,12 @@ 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.
-- `ocean_ice_fluxes::NT`: A NamedTuple of fluxes between the ocean and sea ice, computed at each coupling step.
 """
-struct OceananigansSimulation{SIM, A, OPROP, REMAP, NT, SIC} <:
-       Interfacer.OceanModelSimulation
+struct OceananigansSimulation{SIM, A, OPROP, REMAP} <: Interfacer.OceanModelSimulation
     ocean::SIM
     area_fraction::A
     ocean_properties::OPROP
     remapping::REMAP
-    ocean_ice_fluxes::NT
 end
 
 """
@@ -44,19 +41,14 @@ Specific details about the default model configuration
 can be found in the documentation for `ClimaOcean.ocean_simulation`.
 """
 function OceananigansSimulation(
-    ice_fraction,
+    area_fraction,
     start_date,
     stop_date;
     output_dir,
     comms_ctx = ClimaComms.context(),
 )
     arch = comms_ctx.device isa ClimaComms.CUDADevice ? OC.GPU() : OC.CPU()
 
-    # Initialize the area fraction, making sure no overlap between ice and ocean
-    boundary_space = axes(ice_fraction)
-    area_fraction = CC.Fields.zeros(boundary_space)
-    @. area_fraction = max(min(area_fraction, FT(1) - ice_fraction), FT(0))
-
     # Retrieve EN4 data (monthly)
     # (It requires username and password)
     dates = range(start_date, step = Dates.Month(1), stop = stop_date)
@@ -72,10 +64,6 @@ function OceananigansSimulation(
     z = OC.ExponentialDiscretization(Nz, -depth, 0; scale = 0.85 * depth)
 
     # Regular LatLong because we know how to do interpolation there
-
-    # TODO: When moving to TripolarGrid, note that we need to be careful about
-    # ensuring the coordinate systems align (ie, rotate vectors on the OC grid)
-
     underlying_grid = OC.LatitudeLongitudeGrid(
         arch;
         size = resolution_points,
@@ -88,7 +76,7 @@ function OceananigansSimulation(
     bottom_height = CO.regrid_bathymetry(
         underlying_grid;
         minimum_depth = 30,
-        interpolation_passes = 20,
+        interpolation_passes = 1, # TODO revert
         major_basins = 1,
     )
 
@@ -189,29 +177,7 @@ function OceananigansSimulation(
         ocean.output_writers[:diagnostics] = netcdf_writer
     end
 
-    # TODO make this depend on sea ice model type
-    # Allocate space for the sea ice-ocean (io) fluxes
-    io_bottom_heat_flux = OC.Field{OC.Center, OC.Center, Nothing}(grid)
-    io_frazil_heat_flux = OC.Field{OC.Center, OC.Center, Nothing}(grid)
-    io_salt_flux = OC.Field{OC.Center, OC.Center, Nothing}(grid)
-    x_momentum = OC.Field{OC.Face, OC.Center, Nothing}(grid)
-    y_momentum = OC.Field{OC.Center, OC.Face, Nothing}(grid)
-
-    ocean_ice_fluxes = (
-        interface_heat = io_bottom_heat_flux,
-        frazil_heat = io_frazil_heat_flux,
-        salt = io_salt_flux,
-        x_momentum = x_momentum,
-        y_momentum = y_momentum,
-    )
-
-    sim = OceananigansSimulation(
-        ocean,
-        area_fraction,
-        ocean_properties,
-        remapping,
-        ocean_ice_fluxes,
-    )
+    sim = OceananigansSimulation(ocean, area_fraction, ocean_properties, remapping)
     return sim
 end