struct Buoyancy{M, G}
                   model :: M
    vertical_unit_vector :: G
end

struct ZDirection end

"""
    Buoyancy(; model, vertical_unit_vector=ZDirection())

Uses a given buoyancy `model` to create buoyancy in a model. The optional keyword argument 
`vertical_unit_vector` can be used to specify the direction opposite to the gravitational
acceleration (which we take here to mean the "vertical" direction).

Example
=======

```julia
using Oceananigans

grid = RegularRectilinearGrid(size=(1, 8, 8), extent=(1, 1000, 100))
θ = 45 # degrees
g̃ = (0, sind(θ), cosd(θ))

buoyancy = Buoyancy(model=BuoyancyTracer(), vertical_unit_vector=g̃)

model = NonhydrostaticModel(grid=grid, buoyancy=buoyancy, tracers=:b)
```
"""
function Buoyancy(; model, vertical_unit_vector=ZDirection())
    vertical_unit_vector = validate_vertical_unit_vector(vertical_unit_vector)
    return Buoyancy(model, vertical_unit_vector)
end

validate_vertical_unit_vector(ĝ::ZDirection) = ĝ

function validate_vertical_unit_vector(ĝ)
    length(ĝ) == 3 || throw(ArgumentError("`vertical_unit_vector` must have length 3"))

    gx, gy, gz = ĝ

    gx^2 + gy^2 + gz^2 ≈ 1 ||
        throw(ArgumentError("`vertical_unit_vector` must be a unit vector with g[1]² + g[2]² + g[3]² ≈ 1"))

    return tuple(ĝ...)
end

@inline ĝ_x(buoyancy) = @inbounds buoyancy.vertical_unit_vector[1]
@inline ĝ_y(buoyancy) = @inbounds buoyancy.vertical_unit_vector[2]
@inline ĝ_z(buoyancy) = @inbounds buoyancy.vertical_unit_vector[3]

@inline ĝ_x(::Buoyancy{M, ZDirection}) where M = 0
@inline ĝ_y(::Buoyancy{M, ZDirection}) where M = 0
@inline ĝ_z(::Buoyancy{M, ZDirection}) where M = 1

#####
##### For convenience
#####

@inline required_tracers(bm::Buoyancy) = required_tracers(bm.model)

@inline get_temperature_and_salinity(bm::Buoyancy, C) = get_temperature_and_salinity(bm.model, C)

@inline ∂x_b(i, j, k, grid, b::Buoyancy, C) = ∂x_b(i, j, k, grid, b.model, C)
@inline ∂y_b(i, j, k, grid, b::Buoyancy, C) = ∂y_b(i, j, k, grid, b.model, C)
@inline ∂z_b(i, j, k, grid, b::Buoyancy, C) = ∂z_b(i, j, k, grid, b.model, C)

@inline top_buoyancy_flux(i, j, grid, b::Buoyancy, args...) = top_buoyancy_flux(i, j, grid, b.model, args...)

regularize_buoyancy(b) = b
regularize_buoyancy(b::AbstractBuoyancyModel) = Buoyancy(model=b)
