Generalize ErtelPotentialVorticity and RichardsonNumber to work when model.buoyancy is not BuoyancyTracer

[tomchor]

Easy fix. I wasn't aware that the function Oceananigans.BuoyancyFormulations.buoyancy existed when I first coded EPV.

[tomchor]

@glwagner I'd like your opinion on something. In a few of the CI runs, I get errors like these with the LagrangianAveragedDynamicSmag:

Where the two numbers are small, but differently small. Generally one of them is order 1e-16 and other is 1e-30. It doesn't happen every time, but often enough that I've noticed. That's what made the last CI run fail btw. Specifically this is the test that failed this previous time:

Oceanostics.jl/test/runtests.jl

Lines 363 to 371 in 56f7f8f

	if model isa NonhydrostaticModel
	χ = TracerVarianceTendency(model, :b)
	χ_field = compute!(Field(χ))
	@test χ isa AbstractOperation
	@test χ_field isa Field
	@compute ∂ₜc² = Field(Average(TracerVarianceTendency(model, :b)))
	@test ≈(Array(interior(ε̄ₚ, 1, 1, 1)), -Array(interior(∂ₜc², 1, 1, 1)), rtol=1e-10, atol=2*eps())
	end

Which compares the KE dissipation from the closure model with the total tendency. Since all that's there is the closure model and no other tendencies, they should be the same. I know that eps(Float64) ~ 1e-16 so I'm wondering if this is just natural uncertainty due to limited machine precision, of if there's something else at play here (since the tolerance is already 2*eps() for that test).

Do you have thoughts?

EDIT: FYI I just restarted that run. So if things are passing now, it's because of the normal variability in these tests.