initial ClimaSeaIce infrastructure changes

Author: juliasloan25

docs/src/fieldexchanger.md

@@ -12,7 +12,7 @@ The specific fields that are exchanged depend on the requirements of the compone
 The fields imported from the atmosphere to the coupler are specified by extending `FieldExchanger.import_atmos_fields!`
 The default `import_atmos_fields!` imports radiative fluxes, liquid precipitation, and snow precipitation.
 
-The fields of a component model that get updated by the coupler are specified by extending `FieldExchanger.update_sim!`
+The fields of a component model that get updated by the coupler are specified by extending `FieldExchanger.update_sim!`.
 The default `update_sim!` for an atmosphere model updates the direct and diffuse surface albedos,
 the surface temperature, and the turbulent fluxes.
 The default `update_sim!` for a surface model updates the air density, radiative fluxes,
@@ -21,18 +21,24 @@ These updates are done via the `update_field!` function, which must be extended
 particular variable and component model.
 If an `update_field!` function is not defined for a particular component model, it will be ignored.
 
+Note that turbulent fluxes are not updated in `update_sim!`, but rather via
+`FluxCalculator.update_turbulent_fluxes!`, where fluxes are computed between
+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!
 ```

docs/src/fluxcalculator.md

@@ -11,24 +11,23 @@ turbulent fluxes and ancillary quantities, such as the Obukhov length, using
 function is called at the end of each coupling step.
 
 All the quantities computed in `turbulent_fluxes!` are calculated
-separately for each surface type using the
+separately for each surface model using the
 [`FluxCalculator.compute_surface_fluxes!`](@ref) function. This function can be
 extended by component models if they need specific type of flux calculation, and
 a default is provided for models that can use the standard flux calculation.
-The final result of `turbulent_fluxes!` is an area-weighted sum of
-all the contributions of the various surfaces.
 
 The default method of [`FluxCalculator.compute_surface_fluxes!`](@ref), in turn,
 calls [`FluxCalculator.get_surface_fluxes`](@ref). This function uses a thermal
 state obtained by using the model surface temperature, extrapolates atmospheric
 density adiabatically to the surface, and with the surface humidity (if
 available, if not, assuming a saturation specific humidity for liquid phase).
-`compute_surface_fluxes!` also updates the component internal fluxes fields with
+`compute_surface_fluxes!` also updates the component internal fluxes fields via
 [`FluxCalculator.update_turbulent_fluxes!`](@ref), and adds the area-weighted
 contribution from this component model to the `CoupledSimulation` fluxes fields.
-Any extension of this function for a particular surface model is also expected
-to update both the component models' internal fluxes and the CoupledSimulation
-object's fluxes fields.
+
+Any extension of `FluxCalculator.compute_surface_fluxes!` for a particular
+surface model is also expected to update both the component models' internal
+fluxes and the CoupledSimulation object's fluxes fields.
 
 [`FluxCalculator.compute_surface_fluxes!`](@ref) sets:
 - the flux of momenta, `F_turb_ρτxz`, `F_turb_ρτyz`;

docs/src/interfacer.md

@@ -223,11 +223,26 @@ for the following properties:
 | `surface_direct albedo`  | bulk direct surface albedo                                     |         |
 | `surface_diffuse albedo` | bulk diffuse surface albedo                                    |         |
 | `surface_temperature`    | surface temperature                                            | K       |
+| `ice_concentration`      | sea ice concentration (*sea ice models only*)                  |         |
 
 !!! note
     `area_fraction` is expected to be defined on the boundary space of the simulation,
     while all other fields will likely be on the simulation's own space.
 
+!!! note "Sea ice concentration vs. area fraction"
+Sea ice models are expected to provide both `area_fraction` and `ice_concentration`.
+This may seem redundant, but there are subtle differences between the two.
+`ice_concentration` is internal to the ice model and may be determined
+via a prognostic variable, prescribed data, etc. `area_fraction` is defined
+at the coupler level and may follow some constraints that don't apply to `ice_concentration`.
+For example, we require that surface model area fractions sum to 1 at all points;
+this constraint is enforced for `area_fraction`, but not for `ice_concentration`.
+Additionally, since `area_fraction` is a coupler-defined concept, it is defined on
+the coupler boundary space, whereas `ice_concentration` exists on the model's space.
+Generally, `ice_concentration` and `area_fraction` should largely agree,
+with differences only arising from `area_fraction` corrections and
+the fields existing on different spaces.
+
 - `update_field!(::SurfaceModelSimulation, ::Val{property}, field)`:
 A function to update the value of property in the component model
 simulation, using the values in `field` passed from the coupler

experiments/ClimaEarth/components/ocean/climaocean_helpers.jl

@@ -0,0 +1,126 @@
+"""
+    to_node(pt::CC.Geometry.LatLongPoint)
+
+Transform `LatLongPoint` into a tuple (long, lat, 0), where the 0 is needed because we only
+care about the surface.
+"""
+@inline to_node(pt::CC.Geometry.LatLongPoint) = pt.long, pt.lat, zero(pt.lat)
+# This next one is needed if we have "LevelGrid"
+@inline to_node(pt::CC.Geometry.LatLongZPoint) = pt.long, pt.lat, zero(pt.lat)
+
+"""
+    map_interpolate(points, oc_field::OC.Field)
+
+Interpolate the given 3D field onto the target points.
+
+If the underlying grid does not contain a given point, return 0 instead.
+
+TODO: Use a non-allocating version of this function (simply replace `map` with `map!`)
+"""
+function map_interpolate(points, oc_field::OC.Field)
+    loc = map(L -> L(), OC.Fields.location(oc_field))
+    grid = oc_field.grid
+    data = oc_field.data
+
+    # TODO: There has to be a better way
+    min_lat, max_lat = extrema(OC.φnodes(grid, OC.Center(), OC.Center(), OC.Center()))
+
+    map(points) do pt
+        FT = eltype(pt)
+
+        # The oceananigans grid does not cover the entire globe, so we should not
+        # interpolate outside of its latitude bounds. Instead we return 0
+        min_lat < pt.lat < max_lat || return FT(0)
+
+        fᵢ = OC.Fields.interpolate(to_node(pt), data, loc, grid)
+        convert(FT, fᵢ)::FT
+    end
+end
+
+"""
+    surface_flux(f::OC.AbstractField)
+
+Extract the top boundary conditions for the given field.
+"""
+function surface_flux(f::OC.AbstractField)
+    top_bc = f.boundary_conditions.top
+    if top_bc isa OC.BoundaryCondition{<:OC.BoundaryConditions.Flux}
+        return top_bc.condition
+    else
+        return nothing
+    end
+end
+
+function Interfacer.remap(field::OC.Field, target_space)
+    return map_interpolate(CC.Fields.coordinate_field(target_space), field)
+end
+
+function Interfacer.remap(operation::OC.AbstractOperations.AbstractOperation, target_space)
+    evaluated_field = OC.Field(operation)
+    OC.compute!(evaluated_field)
+    return Interfacer.remap(evaluated_field, target_space)
+end
+
+"""
+    set_from_extrinsic_vectors!(vectors, grid, u_cc, v_cc)
+
+Given the extrinsic vector components `u_cc` and `v_cc` as `Center, Center`
+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 = OC.Architectures.architecture(grid)
+
+    # Rotate vectors onto the 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)
+    OC.fill_halo_regions!(v_cc)
+
+    # Interpolate the vectors to face/center and center/face respectively
+    OC.Utils.launch!(
+        arch,
+        grid,
+        :xy,
+        _interpolate_vectors!,
+        vectors.u,
+        vectors.v,
+        grid,
+        u_cc,
+        v_cc,
+    )
+    return nothing
+end
+
+"""
+    _rotate_vectors!(τx, τy, grid)
+
+Rotate the velocities from the extrinsic coordinate system to the intrinsic
+coordinate system.
+"""
+@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
+    τxr, τyr = OC.Operators.intrinsic_vector(i, j, 1, grid, τx, τy)
+    @inbounds begin
+        τx[i, j, 1] = τxr
+        τy[i, j, 1] = τyr
+    end
+end
+
+"""
+    _interpolate_vectors!(τx, τy, grid, τx_cc, τy_cc)
+
+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_vectors!(τx, τy, grid, τx_cc, τy_cc)
+    # Use `k = 1` to index into the reduced Fields
+    i, j = @index(Global, NTuple)
+    @inbounds begin
+        τ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

@@ -7,6 +7,8 @@ import Thermodynamics as TD
 import ClimaOcean.EN4: download_dataset
 using KernelAbstractions: @kernel, @index, @inbounds
 
+include("climaocean_helpers.jl")
+
 """
     OceananigansSimulation{SIM, A, OPROP, REMAP}
 
@@ -207,70 +209,6 @@ end
 Interfacer.step!(sim::OceananigansSimulation, t) =
     OC.time_step!(sim.ocean, float(t) - sim.ocean.model.clock.time)
 
-# We always want the surface, so we always set zero(pt.lat) for z
-"""
-    to_node(pt::CA.ClimaCore.Geometry.LatLongPoint)
-
-Transform `LatLongPoint` into a tuple (long, lat, 0), where the 0 is needed because we only
-care about the surface.
-"""
-@inline to_node(pt::CA.ClimaCore.Geometry.LatLongPoint) = pt.long, pt.lat, zero(pt.lat)
-# This next one is needed if we have "LevelGrid"
-@inline to_node(pt::CA.ClimaCore.Geometry.LatLongZPoint) = pt.long, pt.lat, zero(pt.lat)
-
-"""
-    map_interpolate(points, oc_field::OC.Field)
-
-Interpolate the given 3D field onto the target points.
-
-If the underlying grid does not contain a given point, return 0 instead.
-
-TODO: Use a non-allocating version of this function (simply replace `map` with `map!`)
-"""
-function map_interpolate(points, oc_field::OC.Field)
-    loc = map(L -> L(), OC.Fields.location(oc_field))
-    grid = oc_field.grid
-    data = oc_field.data
-
-    # TODO: There has to be a better way
-    min_lat, max_lat = extrema(OC.φnodes(grid, OC.Center(), OC.Center(), OC.Center()))
-
-    map(points) do pt
-        FT = eltype(pt)
-
-        # The oceananigans grid does not cover the entire globe, so we should not
-        # interpolate outside of its latitude bounds. Instead we return 0
-        min_lat < pt.lat < max_lat || return FT(0)
-
-        fᵢ = OC.Fields.interpolate(to_node(pt), data, loc, grid)
-        convert(FT, fᵢ)::FT
-    end
-end
-
-"""
-    surface_flux(f::OC.AbstractField)
-
-Extract the top boundary conditions for the given field.
-"""
-function surface_flux(f::OC.AbstractField)
-    top_bc = f.boundary_conditions.top
-    if top_bc isa OC.BoundaryCondition{<:OC.BoundaryConditions.Flux}
-        return top_bc.condition
-    else
-        return nothing
-    end
-end
-
-function Interfacer.remap(field::OC.Field, target_space)
-    return map_interpolate(CC.Fields.coordinate_field(target_space), field)
-end
-
-function Interfacer.remap(operation::OC.AbstractOperations.AbstractOperation, target_space)
-    evaluated_field = OC.Field(operation)
-    OC.compute!(evaluated_field)
-    return Interfacer.remap(evaluated_field, target_space)
-end
-
 Interfacer.get_field(sim::OceananigansSimulation, ::Val{:area_fraction}) = sim.area_fraction
 
 # TODO: Better values for this
@@ -375,70 +313,6 @@ function FluxCalculator.update_turbulent_fluxes!(sim::OceananigansSimulation, fi
     return nothing
 end
 
-"""
-    set_from_extrinsic_vectors!(vectors, grid, u_cc, v_cc)
-
-Given the extrinsic vector components `u_cc` and `v_cc` as `Center, Center`
-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
-
-    # Rotate vectors onto the grid
-    OC.Utils.launch!(arch, grid, :xy, _rotate_velocities!, 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)
-    OC.fill_halo_regions!(v_cc)
-
-    # Interpolate the vectors to face/center and center/face respectively
-    OC.Utils.launch!(
-        arch,
-        grid,
-        :xy,
-        _interpolate_velocities!,
-        vectors.u,
-        vectors.v,
-        grid,
-        u_cc,
-        v_cc,
-    )
-    return nothing
-end
-
-"""
-    _rotate_velocities!(u, v, grid)
-
-Rotate the velocities from the extrinsic coordinate system to the intrinsic
-coordinate system.
-"""
-@kernel function _rotate_velocities!(u, v, 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)
-    @inbounds begin
-        u[i, j, 1] = ur
-        v[i, j, 1] = vr
-    end
-end
-
-"""
-    _interpolate_velocities!(u, v, grid, u_cc, v_cc)
-
-Interpolate the input velocities `u_cc` and `v_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)
-    # 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)
-    end
-end
-
 function Interfacer.update_field!(sim::OceananigansSimulation, ::Val{:area_fraction}, field)
     sim.area_fraction .= field
     return nothing

experiments/ClimaEarth/components/ocean/prescr_seaice.jl

@@ -165,6 +165,8 @@ end
 # extensions required by Interfacer
 Interfacer.get_field(sim::PrescribedIceSimulation, ::Val{:area_fraction}) =
     sim.integrator.p.area_fraction
+Interfacer.get_field(sim::PrescribedIceSimulation, ::Val{:ice_concentration}) =
+    sim.integrator.p.area_fraction
 Interfacer.get_field(sim::PrescribedIceSimulation, ::Val{:roughness_buoyancy}) =
     sim.integrator.p.params.z0b
 Interfacer.get_field(sim::PrescribedIceSimulation, ::Val{:roughness_momentum}) =