initial ClimaSeaIce changes

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

@@ -225,11 +225,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,129 @@
+"""
+    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.
+
+Note: `map_interpolate` does not support interpolation from `Field`s defined on
+`OrthogononalSphericalShellGrids` such as the `TripolarGrid`.
+
+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_vector!(vector, 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_vector!(vector, grid, u_cc, v_cc)
+    arch = OC.Architectures.architecture(grid)
+
+    # Rotate vector components onto the grid
+    OC.Utils.launch!(arch, grid, :xy, _rotate_vector!, u_cc, v_cc, grid)
+
+    # Fill halo regions with the rotated vector components so we can use them to interpolate
+    OC.fill_halo_regions!(u_cc)
+    OC.fill_halo_regions!(v_cc)
+
+    # Interpolate the vector components to face/center and center/face respectively
+    OC.Utils.launch!(
+        arch,
+        grid,
+        :xy,
+        _interpolate_vector!,
+        vector.u,
+        vector.v,
+        grid,
+        u_cc,
+        v_cc,
+    )
+    return nothing
+end
+
+"""
+    _rotate_vector!(τx, τy, grid)
+
+Rotate the velocities from the extrinsic coordinate system to the intrinsic
+coordinate system.
+"""
+@kernel function _rotate_vector!(τ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_vector!(τ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_vector!(τ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