Spherical Harmonic Transform

The following sections outline the implementation of the spherical harmonic transform (in short spectral transform) between the coefficients of the spherical harmonics (the spectral space) and the grid space which can be any of the Implemented grids as defined by RingGrids. This includes the classical full Gaussian grid, a regular longitude-latitude grid called the full Clenshaw grid (FullClenshawGrid), ECMWF's octahedral Gaussian grid[^Malardel2016], and HEALPix grids[^Gorski2004]. SpeedyWeather.jl's spectral transform module SpeedyTransforms is grid-flexible and can be used with any of these, see Grids.

SpeedyTransforms is a module too!

SpeedyTransform is the underlying module that SpeedyWeather imports to transform between spectral and grid-point space, which also implements Derivatives in spherical coordinates. You can use this module independently of SpeedyWeather for spectral transforms, see SpeedyTransforms.

Inspiration

The spectral transform implemented by SpeedyWeather.jl follows largely Justin Willmert's CMB.jl and SphericalHarmonicTransforms.jl package and makes use of AssociatedLegendrePolynomials.jl and FFTW.jl for the Fourier transform. Justin described his work in a Blog series [^Willmert2020].

Spherical harmonics

The spherical harmonics of degree and order over the longitude   and latitudes   , are

with being the pre-normalized associated Legendre polynomials, and are the complex exponentials (the Fourier modes). Together they form a set of orthogonal basis functions on the sphere. For an interactive visualisation of the spherical harmonics, see here.

Latitudes versus colatitudes

The implementation of the spectral transforms in SpeedyWeather.jl uses colatitudes   (0 at the north pole) but the dynamical core uses latitudes    ( at the north pole). Note: We may also use latitudes in the spherical harmonic transform in the future for consistency.

Synthesis (spectral to grid)

The synthesis (or inverse transform) takes the spectral coefficients and transforms them to grid-point values (for the sake of simplicity first regarded as continuous). The synthesis is a linear combination of the spherical harmonics with non-zero coefficients.

We obtain an approximation with a finite set of by truncating the series in both degree and order somehow. Most commonly, a triangular truncation is applied, such that all degrees after   are discarded. Triangular because the retained array of the coefficients looks like a triangle. Other truncations like rhomboidal have been studied[^Daley78] but are rarely used since. Choosing also constrains and determines the (horizontal) spectral resolution. In SpeedyWeather.jl this resolution as chosen as trunc when creating the SpectralGrid.

For being a real-valued there is a symmetry

meaning that the coefficients at and are the same, but the sign of the real and imaginary component can be flipped, as denoted with the   and the complex conjugate . As we are only dealing with real-valued fields anyway, we therefore never have to store the negative orders and end up with a lower triangular matrix of size     or technically   where is the truncation trunc. One is added here because the degree and order use 0-based indexing but sizes (and so is Julia's indexing) are 1-based.

For correctness we want to mention here that vector quantities require one more degree due to the recurrence relation in the Meridional derivative. Hence for practical reasons all spectral fields are represented as a lower triangular matrix of size    . And the scalar quantities would just not make use of that last degree, and its entries would be simply zero. We will, however, for the following sections ignore this and only discuss it again in Meridional derivative.

Another consequence of the symmetry mentioned above is that the zonal harmonics, meaning have no imaginary component. Because these harmonics are zonally constant, a non-zero imaginary component would rotate them around the Earth's axis, which, well, doesn't actually change a real-valued field.

Following the notation of [^Willmert2020] we can therefore write the truncated synthesis as

The   factor using the Kronecker is used here because of the symmetry we have to count both the   order pairs (hence the ) except for the zonal harmonics which do not have a pair.

Another symmetry arises from the fact that the spherical harmonics are either symmetric or anti-symmetric around the Equator. There is an even/odd combination of degrees and orders so that the sign flips like a checkerboard

This means that one only has to compute the Legendre polynomials for one hemisphere and the other one follows with this equality.

Analysis (grid to spectral)

Starting in grid-point space we can transform a field into the spectral space of the spherical harmonics by

Note that this notation again uses colatitudes , for latitudes the becomes a and the bounds have to be changed accordingly to  . A discretization with grid points at location , indexed by can be written as [^Willmert2020]

The hat on just means that it is an approximation, or an estimate of the true  . We can essentially make use of the same symmetries as already discussed in Synthesis. Splitting into the Fourier modes and the Legendre polynomials (which are defined over   so the argument maps them to colatitudes) we have

So the term in brackets can be separated out as long as the latitude is constant, which motivates us to restrict the spectral transform to grids with iso-latitude rings, see Grids. Furthermore, this term can be written as a fast Fourier transform, if the are equally spaced on the latitude ring . Note that the in-ring index can depend on the ring index , so that one can have reduced grids, which have fewer grid points towards the poles, for example. Also the Legendre polynomials only have to be computed for the colatitudes (and in fact only one hemisphere, due to the north-south symmetry discussed in the Synthesis). It is therefore practical and efficient to design a spectral transform implementation for ring grids, but there is no need to hardcode a specific grid.

Spectral packing

Spectral packing is the way how the coefficients of the spherical harmonics of a given spectral field are stored in an array. SpeedyWeather.jl uses the conventional spectral packing of degree and order as illustrated in the following image (Cyp, CC BY-SA 3.0, via Wikimedia Commons)

Every row represents an order  , starting from   at the top. Every column represents an order  , starting from   on the left. The coefficients of these spherical harmonics are directly mapped into a matrix as

which is consistently extended for higher degrees and orders. Consequently, all spectral fields are lower-triangular matrices with complex entries. The upper triangle excluding the diagonal are zero. Note that internally vector fields include an additional degree, such that    (see Derivatives in spherical coordinates for more information). The harmonics with (the first column) are also called zonal harmonics as they are constant with longitude . The harmonics with (the main diagonal) are also called sectoral harmonics as they essentially split the sphere into sectors in longitude without a zero-crossing in latitude.

For correctness it is mentioned here that SpeedyWeather.jl uses a LowerTriangularMatrix type to store the spherical harmonic coefficients. By doing so, the upper triangle is actually not explicitly stored and the data technically unravelled into a vector, but this is hidden as much as possible from the user. For more details see LowerTriangularArrays.

Array indices

For a spectral field a note that due to Julia's 1-based indexing the coefficient is obtained via a[l+1, m+1]. Alternatively, we may index over 1-based l, m but a comment is usually added for clarification.

Fortran SPEEDY does not use the same spectral packing as SpeedyWeather.jl. The alternative packing therein uses   and    as summarized in the following table.

degree	order
0	0	0	0
1	0	0	1
1	1	1	0
2	0	0	2
2	1	1	1
2	2	2	0
3	0	0	3
...	...	...	...

This alternative packing uses the top-left triangle of a coefficient matrix, and the degrees and orders from above are stored at the following indices

This spectral packing is not used in SpeedyWeather.jl but illustrated here for completeness and comparison with Fortran SPEEDY.

SpeedyWeather.jl uses triangular truncation such that only spherical harmonics with   and   are explicitly represented. This is usually described as , with   (although in vector quantities require one more degree in the recursion relation of meridional gradients). For example, T31 is the spectral resolution with   . Note that the degree and order are mathematically 0-based, such that the corresponding coefficient matrix is of size 32x32.

Available horizontal resolutions

Technically, SpeedyWeather.jl supports arbitrarily chosen resolution parameter trunc when creating the SpectralGrid that refers to the highest non-zero degree that is resolved in spectral space. SpeedyWeather.jl will always try to choose an easily-Fourier transformable[^FFT] size of the grid, but as we use FFTW.jl there is quite some flexibility without performance sacrifice. However, this has traditionally lead to typical resolutions that we also use for testing we therefore recommend to use. They are as follows with more details below

trunc	nlon	nlat
31 (default)	96	48	400 km
42	128	64	312 km
63	192	96	216 km
85	256	128	165 km
127	384	192	112 km
170	512	256	85 km
255	768	384	58 km
341	1024	512	43 km
511	1536	768	29 km
682	2048	1024	22 km
1024	3072	1536	14 km
1365	4092	2048	11 km

Some remarks on this table

• This assumes the default quadratic truncation, you can always adapt the grid resolution via the dealiasing option, see Matching spectral and grid resolution

• nlat refers to the total number of latitude rings, see Grids. With non-Gaussian grids, nlat will be one one less, e.g. 47 instead of 48 rings.

• nlon is the number of longitude points on the Full Gaussian Grid, for other grids there will be at most these number of points around the Equator.

• is the horizontal resolution. For a spectral model there are many ways of estimating this[^Randall2021]. We use here the square root of the average area a grid cell covers, see Effective grid resolution

Effective grid resolution

There are many ways to estimate the effective grid resolution of spectral models[^Randall2021]. Some of them are based on the wavelength a given spectral resolution allows to represent, others on the total number of real variables per area. However, as many atmospheric models do represent a considerable amount of physics on the grid (see Parameterizations) there is also a good argument to include the actual grid resolution into this estimate and not just the spectral resolution. We therefore use the average grid cell area to estimate the resolution

with number of grid points over a sphere with radius . However, we have to acknowledge that this usually gives higher resolution compared to other methods of estimating the effective resolution, see [^Randall2021] for a discussion. You may therefore need to be careful to make claims that, e.g. trunc=85 can resolve the atmospheric dynamics at a scale of 165km.

Derivatives in spherical coordinates

Horizontal gradients in spherical coordinates are defined for a scalar field and the latitudes and longitudes as

However, the divergence of a vector field   includes additional scalings

and similar for the curl

The radius of the sphere (i.e. Earth) is . The zonal gradient scales with as the longitudes converge towards the poles (note that describes latitudes here, definitions using colatitudes replace the with a .)

Starting with a spectral field of vorticity and divergence one can obtain stream function and velocity potential by inverting the Laplace operator :

The velocities are then obtained from    following the definition from above and   

How the operators    can be implemented with spherical harmonics is presented in the following sections. However, note that the actually implemented operators differ slightly in their scaling with respect to the radius and the cosine of latitude . For further details see Gradient operators which describes those as implemented in the SpeedyTransforms module. Also note that the equations in SpeedyWeather.jl are scaled with the radius (see Radius scaling) which turns most operators into non-dimensional operators on the unit sphere anyway.

Zonal derivative

The zonal derivative of a scalar field in spectral space is the zonal derivative of all its respective spherical harmonics (now we use for longitudes to avoid confusion with the Legendre polynomials )

So for every spectral harmonic, is obtained from via a multiplication with . Unscaling the -factor is done after transforming the spectral coefficients into grid-point space. As discussed in Radius scaling, SpeedyWeather.jl scales the stream function as   such that the division by radius in the gradients can be omitted. The zonal derivative becomes therefore effectively for each spherical harmonic a scaling with its (imaginary) order . The spherical harmonics are essentially just a Fourier transform in zonal direction and the derivative a multiplication with the respective wave number times imaginary .

Meridional derivative

The meridional derivative of the spherical harmonics is a derivative of the Legendre polynomials for which the following recursion relation applies[^Randall2021], [^Durran2010], [^GFDL], [^Orszag70]

with recursion factors

In the following we use the example of obtaining the zonal velocity from the stream function , which is through the negative meridional gradient. For the meridional derivative itself the leading minus sign has to be omitted. Starting with the spectral expansion

we multiply with to obtain

at which point the recursion from above can be applied. Collecting terms proportional to then yields

To obtain the coefficient of each spherical harmonic of the meridional gradient of a spectral field, two coefficients at   and   have to be combined. This means that the coefficient of a gradient is a linear combination of the coefficients of one higher and one lower degree . As the coefficient with   are zero, the sectoral harmonics ( ) of the gradients are obtained from the first off-diagonal only. However, the   harmonics of the gradients require the   as well as the   harmonics. As a consequence vector quantities like velocity components require one more degree than scalar quantities like vorticity[^Bourke72]. However, for easier compatibility all spectral fields in SpeedyWeather.jl use one more degree , but scalar quantities should not make use of it. Equivalently, the last degree is set to zero before the time integration, which only advances scalar quantities.

In SpeedyWeather.jl, vector quantities like use therefore one more meridional mode than scalar quantities such as vorticity or stream function . The meridional derivative in SpeedyWeather.jl also omits the -scaling as explained for the Zonal derivative and in Radius scaling.

Divergence and curl in spherical harmonics

The meridional gradient as described above can be applied to scalars, such as and in the conversion to velocities   , however, the operators curl   and divergence   in spherical coordinates involve a scaling before the meridional gradient is applied. How to translate this to spectral coefficients has to be derived separately[^Randall2021], [^Durran2010].

The spectral transform of vorticity is

Given that   , we therefore have to evaluate a meridional integral of the form

which can be solved through integration by parts. As   at    only the integral

remains. Inserting the recurrence relation from the Meridional derivative turns this into

Now we expand but only the harmonic will project onto. Let    we then have in total

And the divergence is similar, but   . We have moved the scaling with the radius directly into as further described in Radius scaling.

Laplacian

The spectral Laplacian is easily applied to the coefficients of a spectral field as the spherical harmonics are eigenfunctions of the Laplace operator in spherical coordinates with eigenvalues   divided by the radius squared , i.e. becomes in spectral space. For example, vorticity and stream function are related by   in the barotropic vorticity model. Hence, in spectral space this is equivalent for every spectral mode of degree and order to

This can be easily inverted to obtain the stream function from vorticity instead. In order to avoid division by zero, we set here, given that the stream function is only defined up to a constant anyway.

See also Horizontal diffusion and Normalization of diffusion.

U, V from vorticity and divergence

After having discussed the zonal and meridional derivatives with spherical harmonics as well as the Laplace operator, we can derive the conversion from vorticity and divergence (which are prognostic variables) to   . Both are linear operations that act either solely on a given harmonic (the zonal gradient and the Laplace operator) or are linear combinations between one lower and one higher degree (the meridional gradient). It is therefore computationally more efficient to compute directly from instead of calculating stream function and velocity potential first. In total we have

We have moved the scaling with the radius directly into as further described in Radius scaling.

References

[^Malardel2016]: Malardel S, Wedi N, Deconinck W, Diamantakis M, Kühnlein C, Mozdzynski G, Hamrud M, Smolarkiewicz P. A new grid for the IFS. ECMWF newsletter. 2016; 146(23-28):321. doi: 10.21957/zwdu9u5i

[^Gorski2004]: Górski, Hivon, Banday, Wandelt, Hansen, Reinecke, Bartelmann, 2004. HEALPix: A FRAMEWORK FOR HIGH-RESOLUTION DISCRETIZATION AND FAST ANALYSIS OF DATA DISTRIBUTED ON THE SPHERE, The Astrophysical Journal. doi:10.1086/427976

[^Willmert2020]: Justin Willmert, 2020. justinwillmert.com

• Introduction to Associated Legendre Polynomials (Legendre.jl Series, Part I)

• Calculating Legendre Polynomials (Legendre.jl Series, Part II)

• Pre-normalizing Legendre Polynomials (Legendre.jl Series, Part III)

• Maintaining numerical accuracy in the Legendre recurrences (Legendre.jl Series, Part IV)

• Introducing Legendre.jl (Legendre.jl Series, Part V)

• Numerical Accuracy of the Spherical Harmonic Recurrence Coefficient (Legendre.jl Series Addendum)

• Notes on Calculating the Spherical Harmonics

• More Notes on Calculating the Spherical Harmonics: Analysis of maps to harmonic coefficients

[^Daley78]: Roger Daley & Yvon Bourassa (1978) Rhomboidal versus triangular spherical harmonic truncation: Some verification statistics, Atmosphere-Ocean, 16:2, 187-196, DOI: 10.1080/07055900.1978.9649026

[^Randall2021]: David Randall, 2021. An Introduction to Numerical Modeling of the Atmosphere, Chapter 22.

[^Durran2010]: Dale Durran, 2010. Numerical Methods for Fluid Dynamics, Springer. In particular section 6.2, 6.4.

[^GFDL]: Geophysical Fluid Dynamics Laboratory, The barotropic vorticity equation.

[^FFT]: Depending on the implementation of the Fast Fourier Transform (Cooley-Tukey algorithm, or or the Bluestein algorithm) easily Fourier-transformable can mean different things: Vectors of the length that is a power of two, i.e.   is certainly easily Fourier-transformable, but for most FFT implementations so are   with some positive integers. In fact, FFTW uses algorithms even for prime sizes.

[^Bourke72]: Bourke, W. An Efficient, One-Level, Primitive-Equation Spectral Model. Mon. Wea. Rev. 100, 683–689 (1972). doi:10.1175/1520-0493(1972)100<0683:AEOPSM>2.3.CO;2

[^Orszag70]: Orszag, S. A., 1970: Transform Method for the Calculation of Vector-Coupled Sums: Application to the Spectral Form of the Vorticity Equation. J. Atmos. Sci., 27, 890–895, 10.1175/1520-0469(1970)027<0890:TMFTCO>2.0.CO;2.