[WIP] Implement "MEKE" parameterization with prognostic mesoscale eddy kinetic energy

[glwagner]

This PR implements the "MEKE" mesoscale parameterization proposed by Jansen et al (2015) (see also Kong and Jansen (2020)).

To implement this PR we add a property to HydrostaticFreeSurfaceModel called auxiliary_prognostic_fields which provides a container for closure-specific prognostic fields that need to be stepped forward alongside the rest of the model's prognostic state.

Another change is that we introduce AbstractSkewSymmetricDiffusivity that provides an interface for implementing mesoscale closures with skew and symmetric diffusivities.

The "MEKE" parameterization has a two-dimensional prognostic eddy kinetic energy variable that this feature supports.

It's WIP now, but a few notes are:

• We need a slightly more descriptive name than MEKE or MesoscaleEddyKineticEnergy I think (though it does pronounce well "mee-key")... maybe PrognosticMEKEDiffusivity or something? We may want to distinguish from (or better, combine with?) a similar parameterization with 2D prognostic MEKE called "GEOMETRIC"

• There's probably a way to improve the auxiliary_prognostic_fields design... in particular, I'm wondering if we should use a more hierarchical structure for auxiliary fields that also encompasses diffusivity_fields, something like

struct AuxiliaryFields
    user_defined
    diagnostic_closure_fields
    prognostic_closure_fields
end

The downside is that eddy diffusivities for LES are then buried in model.auxiliary_fields.diagnostic_closure_fields. BUT we can also design an interface for extracting these like eddy_diffusivity(model). Curious what people think about that (@tomchor you've had opinions). Basically it's easier to separate the user API (here, functions that extract properties) from the model struct design (which is motivated more by internal considerations). But it's tradeoffs all the way down so good to discuss.

I also did a little clean up and changed the internal function TurbulenceClosures.DiffusivityFields to TurbulenceClosures.diffusivity_fields. The code is inconsistent about the use of TitleCase (is it a constructor / struct? is it a function?) and there's a bit of clean up to do...

Resolves #2422

[glwagner]

@NoraLoose

[NoraLoose]

Hi @glwagner, thanks for starting this PR! Sorry for chiming in so late - I'm finally getting started with Oceananigans.jl (but am already loving it 🎉).

I don't know the code well enough (yet!) to comment on your code-related questions above. But I can give some general input concerning the MEKE scheme.

1. References: The most up-to-date references for MEKE are Jansen et al. (2019) and Kong and Jansen (2020); in the latter mostly Appendices A b,c.
2. Naming convention: The prognostic MEKE equation is used to inform the GM coefficient via kappa_GM = c * sqrt(2 MEKE) * L_mix with a mixing length L_mix, but also to energetically constrain the backscatter in the momentum equation. See this schematic, modified from Figure 1, Jansen et al. (2019):
   
   So how about PrognosticMEKEDiffusivity and PrognosticMEKEBackscatter? Probably good to keep these 2 separated to be able to switch on/off different terms in the MEKE budget. For instance, as of now, MOM6 OM4 only uses the PrognosticMEKEDiffusivity (and no backscatter), but the CPT and others are looking into the energetically-constrained backscatter term intensively these days!
3. MEKE vs. GEOMETRIC: As you mention above, the two schemes are very similar. There are maybe three notable differences:
   a) energy type: MEKE solves a prognostic equation for EKE, whereas GEOMETRIC solves a prognostic equation for EKE + EPE. To inform GM, the two schemes then use the variable 2 * EKE and EKE + EPE, respectively. Bachman (2017) found that in practice, kappa_GM is not sensitive to whether one solves for EKE (as MEKE does) or EKE + EPE (as GEOMETRIC does). So we could probably start by following the MEKE approach, and just solve for EKE.
   b) GM/energy source term: In MEKE, it is left to the user how to specify the mixing length scale L_mix in kappa_GM = c * sqrt(2 MEKE) * L_mix and in the inferred MEKE source term. (Options are: the deformation scale, Rhines scale, grid scale, frictional halting scale, Eady length scale, or any combination of those, see here.) In GEOMETRIC, the mixing length scale is essentially chosen as the Eady length scale. So if you use MEKE with the Eady length scale L_e (see Appendix Ab in Kong and Jansen, 2021), you recover the GEOMETRIC source term.
   c) sub-grid dissipation: MEKE assumes only bottom drag, GEOMETRIC uses one "catch-all" linear damping term as the energy sink term in the sub-grid energy budget. As highlighted in Mak et al. (2022), this sub-grid dissipation term is probably the big elephant in the room, which the solution is very sensitive to but noone knows what the term should look like.

[glwagner]

@NoraLoose my turn to apologize for a very slow response --- thank you for those explanations, they are very helpful! I think this PR is going a bit stale so we may have to restart it, but when that happens I'll port your comments over to the new PR. Probably we need to wait for #2477, which among other things implements a more stable discretization for GM+Redi a la Griffies et al 1998.

[glwagner]

It's way stale!