Biharmonic viscosity on Cubed-Sphere grid

[jm-c]

It seems that there are 2 potential issues with implementation of biharmonic viscosity when used on cubed-sphere grid:

1. in abstract_scalar_biharmonic_diffusivity_closure.jl, we compute "ζ★ᶠᶠᶜ" (starting at line 86) directly without the special cases for the cubed-sphere corners (as we do for the vorticity), so I think this will not be right at the cube corners. An easy fix could be to use the vorticity function instead.
2. the second issue is related to how we evaluate the horizontal laplacian of u and v, i.e., "∇²hfcᶜ" and "∇²hᶜᶠᶜ", that are used for biharmonic viscosity. Currently, the way these operators are defined in Operators/laplacian_operators.jl is more like a flux-form expression that does not work at the cubed-sphere corners. A way to fix this would be to use more like a vector-invariant analog where the laplacian is evaluated from the vorticity and divergence of the horizontal velocity component.

[glwagner]

Oceananigans.jl/src/TurbulenceClosures/abstract_scalar_biharmonic_diffusivity_closure.jl

Lines 86 to 95 in 3f62872

	@inline function ζ★ᶠᶠᶜ(i, j, k, grid, u, v)
	# These closures seem to be needed to help the compiler infer types
	# (either of u and v or of the function arguments)
	@inline Δy_∇²v(i, j, k, grid, v) = Δy_qᶜᶠᶜ(i, j, k, grid, biharmonic_mask_y, ∇²hᶜᶠᶜ, v)
	@inline Δx_∇²u(i, j, k, grid, u) = Δx_qᶠᶜᶜ(i, j, k, grid, biharmonic_mask_x, ∇²hᶠᶜᶜ, u)
	return 1 / Azᶠᶠᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, Δy_∇²v, v) -
	δyᵃᶠᵃ(i, j, k, grid, Δx_∇²u, u))
	end

[glwagner]

Oceananigans.jl/src/Operators/laplacian_operators.jl

Lines 10 to 18 in 3f62872

	@inline function ∇²hᶠᶜᶜ(i, j, k, grid, u)
	return 1/Vᶠᶜᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, Ax_∂xᶜᶜᶜ, u) +
	δyᵃᶜᵃ(i, j, k, grid, Ay_∂yᶠᶠᶜ, u))
	end
	@inline function ∇²hᶜᶠᶜ(i, j, k, grid, v)
	return 1/Vᶜᶠᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_∂xᶠᶠᶜ, v) +
	δyᵃᶠᵃ(i, j, k, grid, Ay_∂yᶜᶜᶜ, v))
	end