CubedSphere grid construction + tests

[navidcy]

The cubed sphere grid from the MIT as read from DataDeps has the faces ordered:

1. + x direction
2. + y direction
3. + z direction (Arctic)
4. - x direction
5. - y direction
6. - z direction (Antarctica)

So to avoid confusion for (users and developers!) I thought I make sure that the ConformalCubedSphere follows the same convention.

This PR does that + it adds tests to check whether the ConformalCubedSphereGrid metrics from the grid read from file and the grid constructed in Oceananigans.jl are the same.

After this is merged we'll move on to #2867 to constructing ConformalCubedSphereGrid using MultiRegion (an possibly nuke the CubedSpheres module completely). We will have a ConformalCubedSphereGrid (which is the same as the cs32 grid read from file) to bench our MultiRegion implementation with!

(This PR also moves some of the spherical trigonometry utility function to grid_utils.jl.)

Closes #1586

[navidcy]

this is ready for review :)

[glwagner]

We need a sketch with lat-lon lines overlaying the cube so we understand the "+z", etc terminology!

[navidcy]

Yeap, I plan to add some figures in the docs. For now, is the figures at https://clima.github.io/CubedSphere.jl/stable/conformal_cubed_sphere/ helpful? +z is the face that looks in the positive z direction.

[navidcy]

@glwagner asked why there was a warning saying that the metrics near the edge of the OrthogonalSphericalShellGrid might be wrong. I should have clarified there that I meant to write that "The area metrics might be wrong near the boundaries".

I explain here:

To compute the areas near the edge we actually need to know the node locations of the halos, and to do that we need to implement proper halo filling. Why is that? Take for For example, consider a face near the edge, then:

Δxᶠᶠᵃ[1, j, k] = xᶜᶠᵃ[1, j, k] - xᶜᶠᵃ[0, j, k]

but we can't compute the lat-lon position of xᶜᶠᵃ[0, j, k] from the conformal map because it belongs to a different panel. However, from symmetry argument, one can say that this would be just twice the distance from the first cell center to the edge of the panel:

Δxᶠᶠᵃ[1, j, k] = 2 * (xᶜᶠᵃ[1, j, k] - xᶠᶠᵃ[1, j, k])

So that's that. The distance metrics are actually correct, even near the edges.

How do we compute now areas? An approximation is

    @. Azᶜᶜᵃ = Δxᶜᶜᵃ * Δyᶜᶜᵃ
    @. Azᶠᶜᵃ = Δxᶠᶜᵃ * Δyᶠᶜᵃ
    @. Azᶜᶠᵃ = Δxᶜᶠᵃ * Δyᶜᶠᵃ
    @. Azᶠᶠᵃ = Δxᶠᶠᵃ * Δyᶠᶠᵃ

But that's only an approximation because these quadrilaterals are not symmetrical. Anywhere in the interior of the grid we can improve on this by computing the exact area of the spherical quadrilateral given that we know the location of each of its vertices. But for a quadrilateral that touches the edge some of its vertices would belong to the next panels, i.e., to get their location we'd need proper halo filling. E.g., Azᶠᶜᵃ[i, j] has vertices

   (φᶜᶠᵃ[i-1,  j ], λᶜᶠᵃ[i-1,  j ])
   (φᶜᶠᵃ[ i ,  j ], λᶜᶠᵃ[ i ,  j ])
   (φᶜᶠᵃ[ i , j+1], λᶜᶠᵃ[ i , j+1])
   (φᶜᶠᵃ[i-1, j+1], λᶜᶠᵃ[i-1, j+1])

thus Azᶠᶜᵃ[1, j] would require the location of (φᶜᶠᵃ[0, j+1], λᶜᶠᵃ[0, j+1]) and (φᶜᶠᵃ[0, j ], λᶜᶠᵃ[0, j ]).

But now that I'm thinking about it, again from symmetry and also due to the fact that these are right next to the edge, we would expect that the area of quadrilateral with vertices at

   (φᶜᶠᵃ[ 0 ,  j ], λᶜᶠᵃ[ 0 ,  j ])
   (φᶜᶠᵃ[ 1 ,  j ], λᶜᶠᵃ[ 1 ,  j ])
   (φᶜᶠᵃ[ 1 , j+1], λᶜᶠᵃ[ 1 , j+1])
   (φᶜᶠᵃ[ 0 , j+1], λᶜᶠᵃ[ 0 , j+1])

should be double the area of the quadrilateral with vertices

   (φᶠᶠᵃ[ 1 ,  j ], λᶠᶠᵃ[ 1 ,  j ])
   (φᶜᶠᵃ[ 1 ,  j ], λᶜᶠᵃ[ 1 ,  j ])
   (φᶜᶠᵃ[ 1 , j+1], λᶜᶠᵃ[ 1 , j+1])
   (φᶠᶠᵃ[ 1 , j+1], λᶠᶠᵃ[ 1 , j+1])

I'll fix this and remove the warning!