Add SeawaterDensity to Oceananigans.Models

[jbisits]

See #3316.

This PR adds SeawaterDensity which returns a KernelFunctionOperation to compute either the in-situ seawater density or a potential density a some reference geopotential height. To calculate the seawater density SeawaterPolynomials.ρ is used. This means that at this stage you can only compute the seawater density if the equation of state is a BoussinesqEquationOfState because a reference density is required. I was not sure how/if to do this for other SeawaterBuoyancy models so have not tried to implement yet.

I made density_model.jl a new script so that all the new code could be easily viewed and moved to another location if necessary.

[jbisits]

I will fix these tests later today (Sydney time..)!

[glwagner]

I do love tests but I think it might make sense to simplify / reduce the lines added for testing compared to what's in the PR now. Typically we find that maintaining tests is a significant fraction of developer time so it's best to make them as minimal and simple as possible. Awesome PR, thank you!

[jbisits]

@glwagner I have simplified the test script (good idea as it was clunky for me to sieve through and I wrote it..) and the new tests pass on both CPU() and GPU(). The error from the docs build does not look related to this PR as far as I can tell! The docs built on a previous version and I do not think I changed any code related to the docstring.

This PR is already likely enough but one thing that could be done is to move src/Buoyancy/buoyancy_field.jl

# TODO: move to Models
buoyancy(model) = buoyancy(model.buoyancy, model.grid, model.tracers)
buoyancy(b, grid, tracers) = KernelFunctionOperation{Center, Center, Center}(buoyancy_perturbationᶜᶜᶜ, grid, b.model, tracers)
BuoyancyField(model) = Field(buoyancy(model))

to Oceananigans.Models and alter this so instead of returning BuoyancyField it returns the KernelFunctionOperation

sewater_buoyancy_perturbation(model) = sewater_buoyancy_perturbation(model.buoyancy, model.grid, model.tracers)
sewater_buoyancy_perturbation(b, grid, tracers) = KernelFunctionOperation{Center, Center, Center}(buoyancy_perturbationᶜᶜᶜ, grid, b.model, tracers)
SewaterBuoyancyPerturbation(model) = sewater_buoyancy_perturbation(model)

to match the behaviour of SeawaterDensity.
Likely better to do another day but if you think it worth doing in this PR let me know.