Diffusive fluxes and advective fluxes have different units

[simone-silvestri]

At the moment the divergence of the diffusive fluxes is calculated as

@inline function ∇_dot_qᶜ(i, j, k, grid, closure::AbstractTurbulenceClosure, diffusivities, tracer_index, args...)
    disc = time_discretization(closure)
    return 1/Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_qᶠᶜᶜ, _diffusive_flux_x, disc, closure, diffusivities, tracer_index, args...) +
                                    δyᵃᶜᵃ(i, j, k, grid, Ay_qᶜᶠᶜ, _diffusive_flux_y, disc, closure, diffusivities, tracer_index, args...) +
                                    δzᵃᵃᶜ(i, j, k, grid, Az_qᶜᶜᶠ, _diffusive_flux_z, disc, closure, diffusivities, tracer_index, args...))
end

with _diffusive_flux equal to $\kappa \partial c$ so it has the units of $[c] m / s$
but the divergence of the advective fluxes is

@inline function div_Uc(i, j, k, grid, advection, U, c)
    return 1/Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, _advective_tracer_flux_x, advection, U.u, c) +
                                    δyᵃᶜᵃ(i, j, k, grid, _advective_tracer_flux_y, advection, U.v, c) +
                                    δzᵃᵃᶜ(i, j, k, grid, _advective_tracer_flux_z, advection, U.w, c))

with _advective_tracer_flux equal to $A\cdot u c$ with units of $[c] m^3 / s$

I propose we use the same convention for diffusive and advective fluxes.

[glwagner]

True, perhaps we should modify advective_*_flux so that it doesn't include the areas, is that your suggestion?