Spectral To Gauss¶
This block performs the following:
multiply
: the inverse of the spherical harmonic transform, from spectral space to a Gaussian mesh.multiplyAD
: the adjoint of the inverse of the spherical harmonic transform.leftInverseMultiply
: the direct spherical harmonic transform from a Gaussian mesh to spectral space.
(Additional notes for developers are included as footnotes.)
Requirements¶
This SABER block uses ECMWF libraries Atlas, ectrans and fiat and requires them to work.
Both the SABER block and the ECMWF libraries have been written to run in MPI and OpenMP contexts.
Example yaml¶
saber outer blocks:
 (...)
 saber block name: spectral to gauss
active variables:
 unbalanced_pressure
 moisture_control_variable
 streamfunction
 velocity_potential
 eastward_wind
outer inverse tolerance: 1e8
Overview of the multiply
method¶
This transforms spectral fields onto a Gaussian mesh. Note that the resolution of the Gaussian mesh is preset implicitly since it can be inferred from the preceding outer blocks or the model geometry (if there are no preceding outer blocks) 1.
The highest total wavenumber represented in the spectral space is currently fixed to \(2 N  1\) where \(N\) is the Gaussian resolution. Formally the spherical harmonic transform can support total wavenumbers up to \(2 N\). However, only \(2 N  1\) are exact. Also the transformation from spectral fields to horizontal winds require a maximum of \(2 N  1\), due to using a recursive formula for calculating the derivative of Legendre polynomials.
For active variables that are neither eastward_wind nor northward_wind, the standard scalar inverse spherical harmonic transform is used 2. Scientific details are explained in more detail in section Analytical representation 2.
For the horizontal wind components the spectral Helmholtz decomposition formula is used in conjunction with the standard scalar inverse spherical harmonic transform. In that case the active variables needs to include either streamfunction and velocity_potential or vorticity and divergence. The spectral fields of those variables need to be available from the the SABER block that is either the central block or one step closer to the central block 3.
Further scientific details are given in section Analytical representation 3.
Overview of the multiplyAD
method¶
As this is the adjoint of the multiply method, it is somewhat similar in nature, but the steps are reversed. It transforms the data from a gridpoint Gaussian mesh to spectral space and uses the adjoint of the scalar inverse spherical harmonic transform 4.
When the outer active variables include eastward_wind, northward_wind the adjoint of the inverse spherical harmonic transform is used with the adjoint of a spectral Helmholtz decomposition transformation 5.
Overview of the leftInverseMultiply
method¶
This uses the direct spherical harmonic transform that converts fields on a Gaussian mesh onto a spectral space 6.
When the outer active variables include eastward_wind, northward_wind the inverse of the spectral Helmholtz decomposition is used with the direct spherical harmonic transform 7.
Further details are given in section Analytical representation 4.
Analytical representation 0: Total, meridional wavenumber and normalised Legendre Polynomials¶
The spherical harmonic formulation considers a triangular truncation of zonal and meridional complex spectral coefficients.
Imagine a Gaussian mesh with \(Resolution=2\) e.g. F2. Then the spectral truncation is \(N = 2 * Resolution1\) which equal \(3\) in this case.
The triangular truncation can be visualised by the diagram below:
m=0 
m=1 
m=2 
m=3 


n=3 
(6)a 
(12)a+(13)bi 
(16)a+(17)bi 
(18)+(19)bi 
n=2 
(4)a 
(10)a+(11)bi 
(14)a+(15)bi 

n=1 
(2)a 
(8)a+(9)bi 

n=0 
(0)a 
The bracketed number refers to the index of the spectral coefficient. The value is represented in the form of a complex number of the form \(a + bi\). Note that the imaginary values for m=0 are zero and as such are not included. Note also that the values for negative m are held implicitly are they are the complex conjugate of the values for positive m.
Spherical harmonics \(Y^m_n (\theta, \lambda)\), where \(\theta\) is the latitude and \(\lambda\) is the longitude, are calculated from the product of
complex exponent \(e^{i m \lambda}\)
an associated Legendre polynomial of degree \(n\) and order \(m\) i.e. \(P^m_n (\sin \theta)\),
and a normalisation coefficient \(s^m_n\)
i.e.
Henceforth in this note we will consider the normalised associated Legendre polynomial, denoted by \(\tilde{P^m_n} = s^m_n P^m_n\).
A few notes:
The normalisation is defined as in Belousov 1962, p.5 eq.8 such that
and the normalisation coefficient is effectively
The CondonShortley phase term (Condon, 1970) is not included in the normalised associated Legendre Polynomial \(\tilde{P^m_n}\).
The normalised associated Legendre Polynomials for \(m=0\) and \(m=1\) are calculated using the NewtonRaphson method. The other normalised Legendre Polynomials are calculated recursively (Wedi, 2013).
Analytical representation 1: The direct spherical harmonic transform¶
The direct complex spherical harmonic transform is given by (Wedi et al, 2013)
where the scalar field on the Gaussian mesh is represented by \(f(\lambda, \theta)\).
The discrete analogue of the above equation involves an FFT for the inner square bracket as documented in (Temperton, 1983). Let us assume that the output of this FFT is \(\xi^m (\theta)\).
For each \(m\) a Gaussian quadrature is employed to get the spectral coefficients
where \(w_k\) is a Gaussian weight at the \(k^{\text{th}}\) Gaussian latitudes. Gaussian latitudes are where the roots of the normalised associated Legendre polyomials are 0 for \(P^0_N\).
The Gaussian weights are calculated as in Swatztrauber (2002), namely from
The Gaussian weights are equal to the ratio of surface area segments that are associated with each Gaussian latitude to the surface area of the sphere. So the sum of the Gaussian weights across all Gaussian latitudes will sum to 1. This is different to the standard textbook set of weights which are typically twice as large (but then use a different normalisation for the normalised associated Legendre polynomials).
Analytical representation 2: The inverse spherical harmonic transform¶
The inverse transform is defined as (Wedi, 2013)
Note that the actual code does not use this formulation, but takes numerous shortcuts to the same result.
Analytical representation 3: The inverse spherical harmonic transform for horizontal winds¶
This section explains the transformations that convert the horizontal spectral vorticity and divergence to horizontal wind components on the Gaussian mesh 8.
A spectral equivalent version of the 2D Helmholtz equation is used for this.
First the spectral vorticity and divergence coefficients, \(\zeta^m_n\) and \(D^m_n\), are converted into spectral streamfunction and velocity potential coefficients, \(\psi^m_n\) and \(\chi^m_n\) using the inverse Laplacian scaling.
where \(\psi^0_0 = \chi^0_0 = 0\).
Then a spectral version of the 2D Helmholtz equation is used to get spectral coefficients U and V. The equations use a Legendre difference relationship to represent the derivative of the Legendre polynomial.
A grid point version of \(U, V\) is given by spherical harmonic synthesis.
The actual wind components are calculated by dividing by \(\cos \theta\).
Analytical representation 4: The direct spherical harmonic transform for horizontal winds¶
This section explains the transformation from horizontal wind components on the Gaussian mesh to the spectral vorticity and divergence fields 9.
The input to this transformation comprises the zonal and meridional wind components which are stored in a single Atlas vector field on a Gaussian mesh. The output is the twodimensional vorticity and divergence stored as spectral coefficients.
In general the steps are the inverse of Analytical representation 3.
The trick that is used to solve the former is to scale the grid point wind components (u,v) into two scalar fields (U,V)
Each \(U\) and \(V\) is decomposed into spectral coefficients and treated as if they were scalars.
Then the inverse of the spectral Helmholtz decomposition (described in Section 3) is performed. This involves solving a tridiagonal banded system of equations for each meridional wavenumber \(m\). The spectral horizontal vorticity and divergence are calculated by using Laplacian scaling
where \(D^0_0 = \zeta^0_0 = 0\).
References¶
Belousov, S. L. (1962), Tables of normalized associated Legendre polynomials, Mathematical tables, vol. 18, Pergamon Press. https://www.google.de/books/edition/Tables_of_Normalized_Associated_Legendre/u_viBQAAQBAJ?hl=en&gbpv=1&dq=belousov+legendre+polynomials&printsec=frontcover
Condon, E. U.; Shortley, G. H. (1970), The Theory of Atomic Spectra, Cambridge, England: Cambridge University Press, OCLC 5388084; Chapter 3.
Swarztrauber, P. N. (2002) On computing the points and weights for GaussLegendre quadrature, SIAM J. Sci. Comput. Vol. 24 (3) pp. 945954, https://epubs.siam.org/doi/abs/10.1137/S1064827500379690
Temperton, C., 1991: On Scalar and Vector Transform Methods for Global Spectral Models. Mon. Wea. Rev., 119, 1303–1307, https://doi.org/10.1175/1520049311951303.1.
Wedi, N. P., M. Hamrud, and G. Mozdzynski, 2013: A Fast Spherical Harmonics Transform for Global NWP and Climate Models. Mon. Wea. Rev., 141, 3450–3461, https://doi.org/10.1175/MWRD1300016.1.
Footnotes (for developers)
 1
The Gaussian mesh geometry is passed to the saber block from the
outerGeometryData
object when the ‘gauss to spectral` SABER block is instantiated. 2
The SABER method
multiplyScalarField
is used for the inverse spherical harmonic transform. Within this method the Atlas methodinvtrans
is called which is an interface to the underlying ECMWF spherical harmonic transform code. 3
In the case where we need to create the eastward and northward components of the wind the following steps are taken in SABER method
multiplyVectorFields
:identify whether inner active variables contain the
streamfunction
andvelocity_potential
spectral fields or horizontalvorticity
anddivergence
fields.if the inner active variables contain
streamfunction
andvelocity_potential
spectral fields then convert them into horizontalvorticity
anddivergence
spectral fields.apply an altas method, called
invtrans_vordiv2wind
, that convertsvorticity
anddivergence
spectral fields into an Atlas vector field holding the horizontal wind components on the Gaussian grid.since Atlas and
ectrans
happens to store the two components as a single field, it is split into separate eastward_wind and northward_wind fields for each longitude, latitude pair.
 4
The adjoint of the scalar inverse spherical harmonic transform is done in the SABER method
multiplyScalarFieldAdj
and within this method the Atlas methodinvtrans_adj
is called. 5
The adjoint of the scalar inverse spherical harmonic transform for the horizontal wind components is done in the SABER method
multiplyVectorFieldsAdj
. The steps in this method are:copy the data from the two fields
eastward_wind
andnorthward_wind
into a single vector Atlas field.apply an Atlas method, called
invtrans_vortdiv2wind_adj
that converts the single horizontal vector field intovorticity
anddivergence
spectral fields.if the inner active variables contain
streamfunction
andvelocity_potential
spectral fields convert the horizontalvorticity
anddivergence
spectral fields into them. The scaling is the same that is used in the equivalentmultiply
method as it is selfadjoint.
 6
The scalar transform is done in SABER method
invertMultiplyScalarFields
using the Atlas methoddirtrans
. 7
When the outer active variables include
eastward_wind
,northward_wind
we do the steps ininvertMultiplyVectorFields
:copy the data from the two fields
eastward_wind
andnorthward_wind
into a single vector Atlas field.apply an altas method, called
dirtrans_wind2vortdiv
that converts the single horizontal vector field intovorticity
anddivergence
spectral fields.for inner active variables
streamfunction
andvelocity_potential
, the associated spectral fields are converted into horizontal vorticity and divergence spectral fields. The scaling is the reciprocal of what is used in the equivalentmultiply
method.
 8
The Atlas routine
invtrans_vortdiv2wind
is the interface to the ECWMF code which calculates the grid point horizontal wind components from the horizontal spectral vorticity and the divergence. 9
The Atlas routine
dirtrans_wind2vortdiv
is the interface to the ECMWF code which calculates the horizontal spectral vorticity and the divergence from the horizontal wind components on the Gaussian mesh.