Update QBO amplitude calculation

[chengzhuzhang]

Change the NFFT0, frequency, and interpolation of period calculations

Description

@jjbenedict and @whannah1 identified an odd behavior that the QBO amplitude peak (bottom figure) gets weaker with longer years of simulation being analysed.

@chengzhuzhang tested with multiple year lengths (1-10, 1-20, …1-100) and can reproduce this behavior.
QBO_diags_tests.pptx

This PR is contributed by @justin-richling and Jack Chen ([email protected]) for coming up with a fix.

• Closes #<ISSUE_NUMBER_HERE>

Checklist

• My code follows the style guidelines of this project
• I have performed a self-review of my own code
• My changes generate no new warnings
• Any dependent changes have been merged and published in downstream modules

If applicable:

• New and existing unit tests pass with my changes (locally and CI/CD build)
• I have added tests that prove my fix is effective or that my feature works
• I have commented my code, particularly in hard-to-understand areas
• I have made corresponding changes to the documentation
• I have noted that this is a breaking change for a major release (fix or feature that would cause existing functionality to not work as expected)

[chengzhuzhang]

@justin-richling Thanks for the PR! I still won't be able to figure out why it is needed to add the default=str argument when json.dump (json.dump(json_dict, outfile, default=str)). But it seems to be a fine workaround for saving the files in json format.

I ran a test with the new algorithm and compared to the original one, the power spectrum density changed a lot, and I'm not sure if it is the intended behavior. Could you provide an explanation of the code change here? Thank you!

Before:

After:

Command line to reproduce on Perlmutter:e3sm_diags qbo --no_viewer --reference_data_path '/global/cfs/cdirs/e3sm/diagnostics/observations/Atm/time-series' --test_data_path '/global/cfs/cdirs/e3sm/chengzhu/tutorial2024/v3.LR.historical_0101/post/atm/180x360_aave/ts/monthly/15yr' --results_dir '/global/cfs/cdirs/e3sm/www/chengzhu/tests/qbofix' --case_id 'QBO-ERA-Interim' --ref_timeseries_input --test_timeseries_input --run_type 'model_vs_obs' --sets 'qbo' --variables 'U' --seasons 'ANN' 'DJF' 'MAM' 'JJA' 'SON' --regions 'global' --regrid_tool 'esmf' --regrid_method 'conservative' --multiprocessing --num_workers '24' --main_title 'QBO index, amplitude, and power spectral density for U' --backend 'mpl' --output_format 'png' --canvas_size_w '1212' --canvas_size_h '1628' --figsize '8.5' '11.0' --dpi '150' --arrows --test_name 'extendedOutput.v3.LR.historical_0101' --short_test_name 'v3.LR.historical_0101' --test_colormap 'cet_rainbow.rgb' --ref_name 'ERA-Interim' --reference_name 'ERA-Interim' --reference_colormap 'cet_rainbow.rgb' --diff_title 'Model - Observations' --diff_colormap 'diverging_bwr.rgb' --granulate 'variables' 'plevs' 'regions' --selectors 'sets' 'seasons' --test_start_yr '2000' --test_end_yr '2014' --ref_start_yr '2000' --ref_end_yr '2014' --start_yr '2000' --end_yr '2014'

[justin-richling]

Sure thing, @chengzhuzhang!

From Jack:

Fast Fourier transform (FFT) is used for the QBO spectra calculator in which power spectra in frequency space is carried out in a discrete form.

The original code used an interpolation approach to obtain power spectra in the period space (inverse of frequency).

However, this approach could not produce consistent results under simulations of different durations.
For a longer simulation, more frequencies will be resolved through FFT which could lower the peak amplitude in the frequency space.

The revised code uses a bin approach which sums all the power within each period bin (2 months).
Consequently, this approach produces robust power spectra.

[jjbenedict]

@justin-richling -- so it looks like in this 15-year demonstration case the higher frequency activity is better resolved and the lower frequency activity is poorly resolved. I assume amplitude line in the bottom panel would become smoother if a longer time window were used, correct?

[justin-richling]

@jjbenedict Correct, for something like a typical 30-year run (ERA Interim), the result will be smoother.

[chengzhuzhang]

Thank you for the explanation @justin-richling!
My understanding is that the original one does provide smoother power spectrum, but the amplitude of the peak can get weaker when more frequency is resolved given longer time series. The new one could be difficult to use for simulations shorter than i.e. 30 years. @jjbenedict Could you provide some recommendation how we should move forward?

[jjbenedict]

@chengzhuzhang I talked with @justin-richling and Jack earlier this week. Justin is actively working with Walter Hannah to implement the wavelet analysis (Walter had earlier written a python script to compute the wavelets). We're thinking it might be good to include in E3SM Diags both the traditional Fourier method (as implemented now, and appropriate for time windows of 30 yr or longer) -and- the wavelet method (which may be better for shorter time windows). We'll have to do some tests with the wavelet method to see how effective it is at shorter time windows, but the expectation is that is will be smoother than Jack's existing Fourier method. Maintaining the Fourier method would also allow greater consistency with QBO plots from previous E3SM versions. We realize the urgency of implementing the wavelet analysis, given that v3 simulations are now being run... we'll do our best to get this set up and tested as quickly as possible.

FYI -- I will be on vacation July 6-22 (though interrupted by a meeting July 16-17, when I will be checking email).