Integrate sea ice in the flux computation

[simone-silvestri]

This PR is a tentative implementation for atmosphere-sea-ice computations to be able to integrate sea-ice in the OceanSeaIceModel.

The objective is trying to reuse as much as possible of the infrastructure in place for atmosphere - ocean fluxes

I think everything is roughly in place, we just need to change the surface properties and implement a way to limit the surface temperature to the melting temperature in the surface temperature computation

[glwagner]

idealized_atmosphere.mp4

[glwagner]

@simone-silvestri can you hold commits on this branch for a bit? I'm working on it. Just got a bunch of conflicts.

[glwagner]

As part of this PR we have some important design decisions to make regarding how fluxes are computed and stored. These decisions are not numerical in nature (we know what to do for the physics) but more questions of the user API and internal code design.

I'll start by noting that the way fluxes are currently computed is not correct. This PR will require 3 types of fluxes: atmosphere-ocean, atmosphere-sea-ice, and sea-ice-ocean. (A fourth flux, land-ocean, is implicitly included in atmosphere-ocean for now but we may need to break this out eventually too.)

Right now we do not explicitly store these three types of fluxes. Instead we have an object that represents the total flux into the ocean, which combines the atmosphere and sea-ice fluxes. For example, when we compute atmospheric fluxes, we store the result in the total ocean surface flux fields, then multiply those by $(1 - \aleph)$ where $\aleph$ is the sea ice concentration. (The sea ice concentration $\aleph$ will play a prominent role in this discussion.) See:

ClimaOcean.jl/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_sea_ice_fluxes.jl

Lines 23 to 27 in e1d0066

	ℵ = ice_concentration[i, j, kᴺ]
	τx[i, j, kᴺ] = τx[i, j, kᴺ] * (1 - ℵ)
	τy[i, j, kᴺ] = τy[i, j, kᴺ] * (1 - ℵ)
	Jᵀ[i, j, kᴺ] = Jᵀ[i, j, kᴺ] * (1 - ℵ)
	Jˢ[i, j, kᴺ] = Jˢ[i, j, kᴺ] * (1 - ℵ)

The next step is to add the sea-ice contribution to the total ocean flux, multiplied by $\aleph$. This actually isn't implemented yet (though we do compute sea-ice-ocean fluxes, we end up overwriting the salinity flux which is wrong of course):

ClimaOcean.jl/src/OceanSeaIceModels/CrossRealmFluxes/sea_ice_ocean_fluxes.jl

Line 54 in e1d0066

Qˢ[i, j, 1] = Δh / Δt * (Sᵢ[i, j, 1] - Sₒ[i, j, Nz])

In proceeding we have a choice to make. We can proceed along the current path, and fix the sea-ice-ocean fluxes to incorporate concentration correctly. This will necessarily assume that the atmosphere-ocean fluxes have been previously computed, and therefore add the sea-ice contribution to them (rather than overwriting the values as we do now, which assumes that the previous values are from a prior time step for example). Note that currently, the sea ice does not require explicit fluxes at the ice-ocean interface, because the ice-ocean interface is assumed to be in thermal equilibrium at the melting temperature. Therefore, we do not functionally require storing the ice-ocean fluxes; merely adding them to the total ocean flux is sufficient. (This could change in the future though.) One advantage of this design is that we do not store unnecessary data in memory.

Another possibility is to explicitly store each interface flux separately. This requires more memory allocation (12 extra 2D fields, I think, in the worst case scenario). For a model that's high resolution in the vertical this is fairly trivial extra storage. But it might represent a lot of memory for a low-vertical resolution in the vertical, nearly tantamount to an extra 3D field. But perhaps even that is not a huge price to pay...

To spell this out, consider that the computational model for the coupled system is effectively like this:

    atmosphere
------------------
     sea ice     |
------------------
      ocean      |
------------------

Consider the flux of some quantity denoted by $J$. There are three fluxes: $J_{ao}$ for atmosphere-ocean, $J_{ai}$ for atmosphere-sea-ice, and $J_{io}$ for sea-ice-ocean. The total flux into the top of the ocean $J_o$ is computed using the sea ice concentration $\aleph$:

$$ J_o = \aleph J_{io} + (1 - \aleph) J_{ao} $$

The total flux into the top of the sea ice is just $J_i^{top} = J_{ai}$ and at the bottom of the sea ice is just $J_i^{bot} = J_{io}$.

Our current strategy is to avoid allocating memory for $J_{ao}$ explicitly; instead we only have $J_o$. So our algorithm is

• Compute $J_{ao}$ and store in $J_o$.
• Multiply $J_o$ by $1 - \aleph$ (see above).
• Add $\aleph J_{io}$ to $J_o$?

Another option though is to explicitly compute $J_{ao}$ and $J_{io}$ and then construct $J_o$ from the two. Which I think makes more sense than what we are doing now. @simone-silvestri what do you think?

Finally note that there are at least 4 fluxes of the form $J$: $\tau_x$, $\tau_y$, $Q$, $J^s$. There can also be additional fluxes for tracers.

[simone-silvestri]

I also think explicitly computing $J_{ao}$ and $J_{io}$ and later constructing $J_o$ might be better. This is because being able to explore the two fluxes individually will make it easier to debug/analyze.
Also, we need a $J_i$ for which $J_{ao}$ does not participate, so it is nice to have the individual cross-component fluxes to be able to later discriminate which contributes to which flux. Finally, for coupling purposes, only $J_{ao}$ is required (actually only the turbulent component of $J_{ao}$, the radiative component might need to be handled differently), so it is good to hang on to each individual turbulent component.