Interface for closures (and other model terms) that introduce auxiliary prognostic variables

[glwagner]

@simone-silvestri were discussing implementing some new closures such as K-epsilon / "Generic Length Scale", and the OSMOSIS scheme used by the NEMO ocean model. Both these closures, and CATKEVerticalDiffusivity, introduce new prognostic variables (such as turbulent kinetic energy, dissipation, or mixed layer depth) that must be evolved alongside the model state.

For CATKE we have hacked in a turbulent kinetic energy equation that gets used when i) the tracer is named :e and ii) one of the closures is CATKEVerticalDiffusivity:

https://github.com/CliMA/Oceananigans.jl/blob/9e620637670cdd4defa32c07c51fd5574520de73/src/Models/HydrostaticFreeSurfaceModels/calculate_hydrostatic_free_surface_tendencies.jl#L99-L105

this works but the design is not sustainable. There are a few issues: first, we can't continue to abuse model.tracers to add new prognostic variables. Some prognostic variables shouldn't be called "tracers" --- like mixed layer depth. Also, we probably want a design in which model.tracers is unaffected by the closure being used. This will make it easier to switch between different closures in, for example, a large-scale ocean simulation. The second issue is that we need to design an interface so that the "prognostic dynamics" of a closure's auxiliary prognostic variable can be implemented and confined to the file where the closure is defined (ie, HydrostaticFreeSurfaceModel shouldn't know anything about CATKE in particular).

One solution is to add a new model property called auxiliary_prognostic_fields. This might behave similarly to model.diffusivity_fields; eg every closure can populate this object in the model constructor. Then we'll add tendency fields for those terms, and time-step them forward. The "tendency kernel function" that get's called will also have be provided by the closure.

Just to make things more logical, we could also rename diffusivity_fields to auxiliary_diagnostic_fields when we do this.

This would all have to come with an API for accessing the data. Right now we have fields(model). We could divide this into two functions: prognostic_fields(model) and diagnostic_fields(model), and then also fields(model) which combines the two. One decision there is whether to "flatten" the lists, or to maintain the hierarchical structure that we currently have in diffusivity_fields.

I also think we need an interface for extracting mixing coefficients / diffusivities (rather than requiring that users spelunk into the model properties to find the information they want). We have a prototype with viscosity and diffusivity functions but it hasn't been used much and isn't really tested. So that would have to be expanded upon as well.

[glwagner]

I'm closing this issue because I'm judging that it's not of current, timely relevance to Oceananigans development. If you would like to make it a higher priority or if you think the issue was closed in error please feel free to re-open.