KE Spectra - best practices
This page documents how to compute kinetic energy spectra from EAM atmosphere simulations.
Summary of best practices:
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 vorticity, divergence and PS on the GLL grid (not the PG2) grid. PS is needed if interpolating to pressure levels
in EAM, use: ‘VOR:I','DIV:I','DYN_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 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 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 seen by Nastrom and Gage 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 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 spectra for higher resolutions now that we are able to run at 6km and 3km resolutions.
Computing KE Spectra
Interpolate to the correct lat/lon grid.
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 (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
To compute spectral from native grid EAM output, we first interpolate to a (N+1)x2N lat/lon cap grid. The KE spectra turns out to be sensitive to both the resolution of the lat/lon grid, and the algorithm used for interpolation. To examine this sensitivity, we look at spectra below from an NE=256 (13km) Aqua planet simulation. This grid has 3072 degrees of freedom along the Equator, meaning it can potentially resolve spherical harmonics up to degree (or wavenumber) 1536 (corresponding to 2dx).
To look at the impacts of the interpolation algorithm, we compute the spectra from a single snapshot using bilinear, integrated bilinear and TR’s highorder (integrated native SE shape functions), and two different resolutions. We plot E(k)/E_0(k), where we normalize all spectra by E_0, the spectra compute using the high order algorithm on the highest resolution 3073x6144 lat/lon cap grid.
Three key points from the data below:
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 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 (GLL 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 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.
Note that in the actual spectra plot, the differences are quite small, with the lower resolution spectra matching very well at low wave numbers, and only slightly below at high wave numbers (zoomed in version on right). Only the yellow curve (1025x2048 high order) has visible roll off different that the highest resolution result.
Time Averaging
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 were used when averaging E(k). Left plot is with daily snapshots, right plot is with 3h snapshots. With 1 month of averaging, the data is still a touch noisy, 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
To generate a TR “highorder” mapping file for GLL data using NCO and TempestRemap:
ncremap -G ttl="lat-lon-cap-grid"#latlon=1537,3072#lat_typ=fv#lon_typ=grn_ctr -g latlon_scrip.nc
GenerateCSMesh --alt --res 256 --out_format Netcdf4 --file atmgrid.g
GenerateOverlapMesh --b atmgrid.g --a latlon_scrip.nc --out overlap_mesh.g
GenerateOfflineMap --in_mesh atmgrid.g --out_mesh latlon_scrip.nc --ov_mesh overlap_mesh.g --in_type cgll --in_np 4 --out_type fv --out_double --out_format Netcdf4 --correct_areas --out_map map_ne256np4_to_1537x3072_highorder.nc
Remap a native grid EAM output file:
ncremap -m map_ne256np4_to_1537x3072_highorder.nc -i casename.cam.h2.0001-02-02-00000.nc -o casename.cam.h2.0001-02-02-00000.latlon.nc
Compute spectra.