O2.3.11 Add land/sea mask to ClimaCore domain

[kmdeck]

Purpose

For global runs, we currently set up a spherical shell domain (extruded FD) and create fields that are defined everywhere on the domain. The issues arise in three areas:

• Interpolation: land parameters are defined over land only. When we interpolate these to the model grid, we need to be aware of the mask so as not to get unphysical land parameters over the ocean. For example, the van genuchten n parameter must be > 1. If over the ocean, this is zero, a cell with a large fraction of ocean will have an unphysical value. This same type of issue arises when we interpolate/remap the state to a lat/lon grid for plotting.
• Solving the equations over the ocean: we need to set up initial conditions in regions where there is no land that wont break anything, and we evaluate tendencies over fake land in the ocean, which is inefficient, and these cant break anything
• Other: there are sometimes multiple masks (coming from different native resolutions of parameter files, e.g.).

A related issue is that ultimately we would like to support lateral flow in the soil. This requires boundary conditions in the horizontal.

The benefits to what we are doing now is that (I think) regridding in the coupler is easy because land and atmos have the exact same space at the surface

Priorities

• In the short term (i.e. done by this fall, for our paper), treating the land as a collection of independent columns is fine. We would not evaluate the model over ocean and not have to worry about setting values of parameters/IC over the ocean. That is, it would be nice to have a custom ClimaCore domain which is a set of independent columns which does not span the entire globe. Then we alleviate basically all of the issues above.

• In the long term (end of next year?) supporting boundary conditions and lateral flow is a goal, but this is not essential and requires a lot more work. Note that there are two levels to this:

• • regional simulations (hybrid box - no mask required, but boundary conditions in the horizontal are required). Here we can learn from Parflow, e.g.,
• • global simulations (mask + horizontal BC) - can we learn from the ocean model/do what they do?

Proposed Work Summer 2024

ClimaLand:

• short term: (what is in place now as of Aug 2024): set ocean values to mean over land, Ignore NaNs over the ocean that result as a function of the simulation. (on CPU, we might get breaking errors, while on GPU, running with NaNs does not break the sim)

ClimaCore:

• Define a new type of domain which is a set of independent columns not defined over the entire globe
• no DSS should be required and this should make the land model more efficient
• Interpolation to the model grid needs to be mask aware
• How would the flux exchange need to adapt if the land model had a very different grid?

Proposed Work 2025

• Enable boundary conditions for the horizontal using a mask. These BC would be applied where the mask transitions from 1->0 (land -> ocean). The prognostic variable fields would still be defined globally.

Notes from a while ago - need to revisit When we have lateral flow we want to make use of the current topology/connectivity of points and how this is set up in ClimaCore. This is because it is easier to use a mask in boundary condition setting in the horizontal than to create a totally unstructured mesh.

I checked with Andre and he confirmed that this is what the ocean model does. there are prognostic variables over land but the tendencies are never evaluated. Only indices which are in the mask = 1 region are loaded onto the GPU.

on the ClimaLand side: This solution means that we need to use the mask on the land side because we will still have fields defined over the ocean.

We will need to:

1. ingest the land/sea mask and store somewhere where it is accessible to the tendency
2. initial state vector = 0 everywhere or we can set to NaN everywhere (?).
3. Set initial conditions using the mask (NaN over ocean)
4. Parameter data should be masked already (NaN over the ocean). When regridding, we'll get back Nans in the ocean for the field. Interpolating to model grid needs to be mask aware
5. Write a macro which can be used in any broadcasted function when mask = 1, like:

#=
This assumes that `mask` is in local scope
=#
macro masked(expr)
    :(@. ifelse(mask == 0, 0, expr))
end

function tendency_func!(Y, p, t)
    (; mask) = p
    (; x, z) = p # fields
    (; ϕ, ψ) = Y # fields
   # @. ϕ = some_tendency(ψ, x, z) # current
  #  @. ϕ = ifelse(mask == 0, 0, some_tendency(ψ, x, z)) # what we want
    @masked ϕ = some_tendency(ψ, x, z) # cleaner way to implement
end

6. Repeat for all broadcasted expressions in the tendency, cache, etc.
7. use the mask in plotting (maybe not if NaN in ocean)

• is memory an issue since we are stored fields everywhere, when land is only 30% of the Earth's surface? We should be computed limited on GPU and not memory limited (this is true for the ocean model).

@sriharshakandala @juliasloan25 @Sbozzolo

### Tasks
- [ ] Use ETOPO 2022 global relief model data to generate element-level `land-sea` mask
- [ ] Add `element mask` to ClimaCore spaces
- [ ] Block FD/SE/tendency operations for masked elements in ClimaCore