...
Spherical harmonic transform available in NCL and pyNGL for lat/lon data.
Interpolate EAM native grid output to a (N+1)x2N lat/lon cap grid with TempestRemap’s “highorder” (native SE shape function) algorithm.
For cubed sphere grid with resolution NE, take N>=NE*6 (degrees of freedom, pole to pole)
Need instantaneous output of (U,V) or (vor,div) vorticity, divergence and PS on the GLL grid (not the PG2) grid. PS is needed if interpolating to pressure levels
in EAM, U and V may be on the PG2 grid in EAMv2. To get GLL output, use U_DYN and V_DYN. VOR and DIV output is always on the GLL grid. fincl3 = , use: ‘VOR:I','DIV:I','DYN_PS:I' or: ‘U_DYN:I’,'V_DYN:I','PS:I'
Need ‘DYN_PS’ instead of ‘PS’ to get PS on the GLL grid instead of the PG2 grid.
For smooth results, average over at least 1 month (post spinup) with 12h snapshots. Running longer rather than with higher frequency I/O results in smoother spectra in Aqua planet. Are there seasonal effects in Real planet?
At NE256, interpolating and computing the spectra for each snapshot takes about 5min and 30GB of memory. At NE1024 this will probably require ~500GB of memory and could take 320min per snapshot
At NE1024 needed to use homme_tool for interpolation (running on 1024 KNL nodes: 20min init time + 5min per 3D snapshot). NCL spherical harmonic transform code required 100GB of memory, taking about 3h per transform.
KE Spectra
KE spectra is useful for understanding and tuning the dissipation mechanisms in the model. For a global model, it can be computed via a vector spherical harmonic transform of the instantaneous velocity field (or two spherical harmonic transforms of the instantaneous vorticity and divergence scalars) on a spherical surface. For all the data below, we interpolate to 250mb (I think because this is the most energetic part of the atmosphere - check this?). Spherical harmonics can be thought of as polynomials in Cartesian coordinates (x,y,z) restricted to the sphere. For each degree k there will be 2k+1 polynomials (spherical harmonics) of total degree k. To compute the KE spectra, E(k), we sum the power (coefficient squared) over all degree k spherical harmonics.
The atmosphere is close to geostrophic balance, resulting in a kinetic energy spectrum that scales like E(k) ~ k^-3 as k ranges from 3000km to 800km. But at 600km (wave number 134), the atmosphere transitions into a k^−5/3 scaling that holds down to a few km. The transition, or breakdown of the geostrophic balance, is still not fully explained. The −5/3 scaling suggests a transition to three-dimensional isotropic turbulence, but this explanation can be quickly ruled out. This transition was first observed seen by Nastrom and Gage , who from aircraft observations. They suggested a combination of the known enstrophy generation from baroclinic instability with an unknown small scale energy source. Later investigators have speculated on where this source could come from: there are pure dynamical explanations and explanations which require convection and moist physics (such as thunderstorms).
Example KE spectra. Plotted is the “compensated” spectra, E(k)*k^5/3, to better illustrate the possible transition to ^-5/3. E(k) was computed from many flow snapshots and then E(k) is averaged over all the results to get the nice curves shown below. Unfortunately, one cannot compute E(k) from time averaged flow fields.
...
In CAM4 at 1/8 degree (14km) resolution the model had a nice Nastrom-Gage transition. This has slowly deteriorated, perhaps due to increased dissipation being added to the model as it the physics has evolved. It may also have been that CAM4 had this transition for the wrong reasons. The breakdown of the spectra into divergence and vorticity components (right) shows that most of the differences are due to small scale divergence - associated with wave activity. CAM5.3 is EAM V0. It will be interesting to see the specta spectra for higher resolutions now that we are able to run at 6km and 3km resolutions.
...
Computing the KE spectra requires performing spherical harmonic transforms. This can be done via the long standing SPHEREPACK Fortran package ( https://www2.cisl.ucar.edu/resources/legacy/spherepack ), which has an easy to use interface in NCL (and pyNGL?as of 2020/7, not available in pyNGL).
SPHEREPACK requires gridded data on either a Gaussian grid or an equal angle lat/lon “cap” grid which includes points at the poles. It relies on a combination of associated Legendre transforms and FFTs. For efficiency the grid should have dimensions N x 2N (Gaussian), or (N+1) x 2N (equal angle cap) with N the product of powers of small primes. The most common lat/lon data is an equal angle offset grid, which avoids pole points (see. CAM-FV Grid Overview ). Data on this offset grid would have to first be interpolated to a cap grid before computing spherical harmonic transforms.
Choice of variables.
SPHEREPACK contains both vector and scalar transforms. KE spectra can be computed from either (U,V) using the vector transform, or (VOR,DIV) using the scalar transforms. We recommend using VOR & DIV above for two reasons:
Interpolating (U,V) from a cubed-sphere grid to a lat/lon grid needs to be done carefully at the poles. The common approach of treating U and V components as scalars and interpolating them individually will be very inaccurate at the poles. This is due to the fact that even for a smooth velocity vector, the individual components U and V will be discontinuous at the poles. Although it’s not clear how much this error would impact the KE spectra. For correct interpolation of (U,V) near the poles, one needs to represent the velocity vector in a more well behaved coordinate system, interpolating these coefficients, and then transform back to a spherical coordinate representation. (The E3SM coupler uses Cartesian coordinates when mapping vectors between components).
Using (VOR,DIV) also allows one to separate the spectra into the rotational and compressible parts. This can be done using (U,V) but requires twice as many spherical harmonic transforms.
Choice of interpolation algorithm and resolution
...
Computing E(k) at lower resolution or with certain algorithms introduces a roll off (reduced energy) and high frequencies. Since the behavior of the spectra, especially the rolloff starting around wave number 200 is of most interest, this should be avoided.
Using “PG2” output (cyan curve below) is not recommended - this is most likely due to the downscaling downsampling of the velocity when transforming from the dynamics to the PG2 physics grid, removing energy from small scales.
“intbilin” algorithm is not recommended (brow and green curves below). This algorithm is very nice for analysis, but does remove energy from small scales. “bilin” is a little better (orange, blue)
Best results are obtained by the highorder algorithm (red and straight line referenceGLL 2048HO and GLL 3072HO).
...
We thus recommend spectra be computed by first transforming to a lat/lon grid using TR’s high-order algorithm.
We next look at the effects of interpolation grid resolution: Again plotting E(k)/E_0(k). Observations:
The HO algorithm is quite good enough (no significant rolloff at wave number 200) at a resolution of 1537x3072 with slight rolloff visible at 1025x2048.
The bilin algorithm requires more than twice as much resolution to match the results of highorder.
...
The spectra above from single snaphots have a lot of noise which can be removed by time averaging. We next analyze how much averaging is necessary from an aqua planet NE256 simulation with 1 month of daily snapshots. The legend indicates how many snapshots (spaced 1 day apart) were used when averaging E(k). With 28 snapshots . Left plot is with daily snapshots, right plot is with 3h snapshots. With 1 month of averaging, the data is still a touch noisy at the lower frequencies, so we recommend at least 60 snapshots. Need to test the minimum spacing: 60 snapshots from hourly data is probably not as good as 60 snapshots from daily data. In the plot , with either daily or 3h output. The 3h data is only slightly less noisy than daily (comparing 233 snapshots of 3h data to 28 snapshots of daily data). 2 months of 12h data would probably be sufficient. In the plots below, each curves are all shifted downward so they would not be plotted on top of each other.
...
Software needed
...