Implement Wheeler-Kiladis diagnostic

Author: ApexCapybara

.zenodo.json

@@ -260,6 +260,11 @@
       "affiliation": "Met Office, UK",
       "name": "Little, Bill"
     },
+    {
+      "affiliation": "DLR, Germany",
+      "name": "Londoño Castillo, Santiago",
+      "orcid": "0000-0002-8437-5031"
+    },
     {
       "affiliation": "BSC, Spain",
       "name": "Loosveldt-Tomas, Saskia"

CITATION.cff

@@ -452,6 +452,11 @@ authors:
     affiliation: "Met Office, UK"
     family-names: Ellis
     given-names: Hannah
+  -
+    affiliation: "DLR, Germany"
+    family-names: Londoño Castillo
+    given-names: Santiago
+    orcid: "https://orcid.org/0000-0002-8437-5031"
 
 
 cff-version: 1.2.0

doc/sphinx/source/recipes/figures/wheeler_kiladis/ERA5_pr_wk_normalized_example.png

@@ -66,6 +66,7 @@ Atmosphere
    recipe_aod_aeronet_assess
    recipe_surface_trace_gas
    recipe_weathertyping
+   recipe_wheeler_kiladis_diagnostic
 
 Climate metrics
 ^^^^^^^^^^^^^^^

doc/sphinx/source/recipes/figures/wheeler_kiladis/MPI-ESM1-2-LR_rlut_wk_normalized_example.png

@@ -0,0 +1,131 @@
+.. _recipe_wheeler_kiladis_diagnostic:
+
+Wheeler-Kiladis tropical wave spectra
+=====================================
+
+Overview
+--------
+
+This recipe computes normalized Wheeler-Kiladis wavenumber-frequency
+spectra for tropical precipitation and outgoing longwave radiation.
+
+The diagnostic is intended to evaluate tropical wave variability,
+including Kelvin waves, equatorial Rossby waves, and broad intraseasonal
+variability such as the Madden-Julian Oscillation.
+
+The recipe uses standard ESMValTool preprocessors to convert units,
+compute daily means, extract the tropical belt between 15S and 15N, and
+regrid the data to a regular latitude-longitude grid. The diagnostic then
+computes a cosine-latitude-weighted equatorial mean, removes leap days,
+removes the time mean, linear trend, and annual harmonics, optionally
+removes the zonal mean, computes two-dimensional Fourier spectra in time
+and longitude, estimates a smoothed background spectrum, and plots
+normalized spectra with optional shallow-water dispersion curves.
+
+The diagnostic saves raw, background, and normalized Wheeler-Kiladis
+spectra as NetCDF files and records provenance information for both data
+outputs and figures.
+
+Available recipes and diagnostics
+---------------------------------
+
+Recipes are stored in ``esmvaltool/recipes/``:
+
+    * ``recipe_wheeler_kiladis_diagnostic.yml``
+
+Diagnostics are stored in ``esmvaltool/diag_scripts/wheeler_kiladis/``:
+
+    * ``wheeler_kiladis.py``: ESMValTool wrapper for running the
+      Wheeler-Kiladis diagnostic.
+    * ``spectra.py``: utilities for preparing equatorial fields and
+      computing wavenumber-frequency spectra.
+    * ``plot.py``: plotting utilities for normalized Wheeler-Kiladis
+      spectra and theoretical dispersion curves.
+
+User settings in recipe
+-----------------------
+
+#. Script ``wheeler_kiladis/wheeler_kiladis.py``
+
+   *Optional settings for script*
+
+   * ``annual_harmonics``: Number of annual harmonics removed during
+     seasonal-cycle correction. Default: 3.
+
+   * ``remove_zonal_mean``: If ``True``, remove the zonal mean before
+     the wavenumber-frequency transform. This emphasizes propagating
+     equatorial disturbances. Default: ``True``.
+
+   * ``segment_length``: Length of each spectral segment in days.
+     Default: 180.
+
+   * ``segment_overlap``: Overlap between consecutive spectral segments
+     in days. Default: 90.
+
+   * ``sampling_frequency_per_day``: Sampling frequency of the input
+     data in samples per day. For daily data, this should be 1.0.
+     Default: 1.0.
+
+   * ``sigma_freq``: Gaussian smoothing width in the frequency direction
+     used to estimate the background spectrum. Default: 4.0.
+
+   * ``sigma_wn``: Gaussian smoothing width in the zonal-wavenumber
+     direction used to estimate the background spectrum. Default: 4.0.
+
+   * ``max_wavenumber``: Maximum zonal wavenumber shown in the plot.
+     Default: 15.
+
+   * ``max_frequency``: Maximum frequency shown in the plot in cycles
+     per day. Default: 0.5.
+
+   * ``period_ticks``: Periods in days shown as labels on the frequency
+     axis. Default: ``[2, 3, 5, 10, 20, 30, 60, 100]``.
+
+   * ``equivalent_depths``: Equivalent depths in metres used to draw
+     theoretical shallow-water dispersion curves. Default:
+     ``[8, 12, 25, 50]``.
+
+   * ``show_dispersion``: If ``True``, overlay theoretical Kelvin and
+     equatorial Rossby wave dispersion curves. Default: ``True``.
+
+   * ``mask_zero_wavenumber``: If ``True``, mask the
+     zonal-wavenumber-zero column in the normalized spectrum plot. This
+     is useful when the zonal mean has been removed. Default: ``True``.
+
+Variables
+---------
+
+* pr (atmos, daily, longitude, latitude, time)
+* rlut (atmos, daily, longitude, latitude, time)
+
+Example plots
+-------------
+
+.. _wk_era5_pr_example:
+.. figure:: /recipes/figures/wheeler_kiladis/ERA5_pr_wk_normalized_example.png
+   :align: center
+
+   Normalized Wheeler-Kiladis spectrum for ERA5 precipitation.
+
+.. _wk_persiann_cdr_pr_example:
+.. figure:: /recipes/figures/wheeler_kiladis/PERSIANN-CDR_pr_wk_normalized_example.png
+   :align: center
+
+   Normalized Wheeler-Kiladis spectrum for PERSIANN-CDR precipitation.
+
+.. _wk_mpi_esm1_2_lr_rlut_example:
+.. figure:: /recipes/figures/wheeler_kiladis/MPI-ESM1-2-LR_rlut_wk_normalized_example.png
+   :align: center
+
+   Normalized Wheeler-Kiladis spectrum for MPI-ESM1-2-LR outgoing
+   longwave radiation.
+
+References
+----------
+
+* Wheeler, M. C. and Kiladis, G. N. (1999): Convectively Coupled
+  Equatorial Waves: Analysis of Clouds and Temperature in the
+  Wavenumber-Frequency Domain.
+* Kiladis, G. N. et al. (2009): Convectively Coupled Equatorial Waves.
+* Hannah et al. (2020): Evaluation of tropical variability using
+  Wheeler-Kiladis-type spectra.

doc/sphinx/source/recipes/figures/wheeler_kiladis/PERSIANN-CDR_pr_wk_normalized_example.png

@@ -402,6 +402,11 @@ authors:
     name: Lledo, Llorenc
     institute: BSC, Spain
     orcid:
+  londono-castillo_santiago:
+    name: Londono Castillo, Santiago
+    institute: DLR, Germany
+    orcid: https://orcid.org/0000-0002-8437-5031
+    github: ApexCapybara
   loosveldt-tomas_saskia:
     name: Loosveldt-Tomas, Saskia
     institute: BSC, Spain

doc/sphinx/source/recipes/index.rst

@@ -0,0 +1 @@
+"""Wheeler-Kiladis diagnostic."""

doc/sphinx/source/recipes/recipe_wheeler_kiladis_diagnostic.rst

@@ -0,0 +1,186 @@
+"""Plot utilities for Wheeler-Kiladis diagnostic."""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+def kelvin_curve(wavenumber, depth_m):
+    """
+    Kelvin wave dispersion curve in cycles per day.
+
+    Parameters
+    ----------
+    wavenumber : array-like
+        Dimensionless zonal wavenumber.
+    depth_m : float
+        Equivalent depth [m].
+    """
+    gravity = 9.80665  # m s^-2
+    earth_radius = 6.371e6  # m
+
+    phase_speed = np.sqrt(gravity * depth_m)  # m s^-1
+    zonal_wavenumber = np.asarray(wavenumber) / earth_radius  # m^-1
+
+    angular_frequency = phase_speed * zonal_wavenumber  # rad s^-1
+    frequency = angular_frequency / (2.0 * np.pi) * 86400.0  # cycles day^-1
+
+    return frequency
+
+
+def equatorial_rossby_curve(wavenumber, depth_m, meridional_mode=1):
+    """
+    Compute the dispersion relation for equatorial Rossby waves.
+
+    Parameters
+    ----------
+    wavenumber : array-like
+        Dimensionless zonal wavenumber.
+    depth_m : float
+        Equivalent depth [m].
+    meridional_mode : int, optional
+        Meridional mode number of the equatorial Rossby wave.
+        Default is 1.
+
+    Returns
+    -------
+    numpy.ndarray
+        Wave frequency [cycles day^-1] corresponding to the input
+        zonal wavenumbers.
+    """
+    earth_rotation = 7.2921159e-5  # rad s^-1
+    earth_radius = 6.371e6  # m
+    gravity = 9.80665  # m s^-2
+
+    beta = 2.0 * earth_rotation / earth_radius  # m^-1 s^-1
+    phase_speed = np.sqrt(gravity * depth_m)  # m s^-1
+
+    wavenumber = np.asarray(wavenumber)  # dimensionless zonal wavenumber
+    zonal_wavenumber = wavenumber / earth_radius  # m^-1
+
+    denominator = (
+        zonal_wavenumber**2
+        + (2 * meridional_mode + 1) * beta / phase_speed  # m^-2
+    )
+
+    angular_frequency = -beta * zonal_wavenumber / denominator  # rad s^-1
+    frequency = angular_frequency / (2.0 * np.pi) * 86400.0  # cycles day^-1
+
+    return frequency
+
+
+def plot_wk_spectrum(
+    spectrum,
+    title,
+    output_file,
+    max_wn=15,
+    max_freq=0.5,
+    levels=None,
+    period_ticks=(2, 3, 5, 10, 20, 30, 60, 100),
+    depths=(8, 12, 25, 50),
+    show_dispersion=True,
+    caption=None,
+):
+    """Plot normalized Wheeler-Kiladis spectrum."""
+    if levels is None:
+        levels = np.linspace(0.8, 4.0, 12)
+
+    plot_data = spectrum.sel(
+        wavenumber=slice(-max_wn, max_wn),
+        frequency=slice(0, max_freq),
+    )
+
+    fig, ax = plt.subplots(figsize=(12, 11.0))
+
+    plot_data.plot.contourf(
+        ax=ax,
+        x="wavenumber",
+        y="frequency",
+        levels=levels,
+        cmap="Spectral_r",
+        extend="both",
+        cbar_kwargs={"label": "normalized power"},
+    )
+
+    ax.axvline(0, color="k", linewidth=0.8)
+
+    period_ticks = np.asarray(period_ticks)
+    freq_ticks = 1.0 / period_ticks
+    valid_ticks = (freq_ticks > 0) & (freq_ticks <= max_freq)
+
+    ax.set_yticks(freq_ticks[valid_ticks])
+    ax.set_yticklabels([f"{period:g}" for period in period_ticks[valid_ticks]])
+
+    ax.set_ylabel("Frequency [cycles day$^{-1}$]")
+    ax.set_xlabel("Zonal wavenumber")
+
+    ax.set_xlim(-max_wn, max_wn)
+    ax.set_ylim(0, max_freq)
+
+    ax.set_title(title)
+
+    for period in [30, 100]:
+        frequency = 1.0 / period
+
+        if frequency <= max_freq:
+            ax.axhline(
+                frequency,
+                color="k",
+                linestyle="--",
+                linewidth=0.8,
+                alpha=0.7,
+            )
+
+    if show_dispersion:
+        kelvin_wn = np.linspace(0.05, max_wn, 800)
+        rossby_wn = -np.linspace(0.05, max_wn, 800)
+
+        for depth in depths:
+            kelvin_freq = kelvin_curve(kelvin_wn, depth)
+            rossby_freq = np.abs(
+                equatorial_rossby_curve(
+                    rossby_wn,
+                    depth,
+                    meridional_mode=1,
+                )
+            )
+
+            kelvin_valid = (
+                np.isfinite(kelvin_freq)
+                & (kelvin_freq >= 0)
+                & (kelvin_freq <= max_freq)
+            )
+
+            rossby_valid = (
+                np.isfinite(rossby_freq)
+                & (rossby_freq >= 0)
+                & (rossby_freq <= max_freq)
+            )
+
+            ax.plot(
+                kelvin_wn[kelvin_valid],
+                kelvin_freq[kelvin_valid],
+                color="k",
+                linewidth=1.0,
+            )
+
+            ax.plot(
+                rossby_wn[rossby_valid],
+                rossby_freq[rossby_valid],
+                color="k",
+                linewidth=1.0,
+            )
+    if caption is not None:
+        fig.text(
+            0.5,
+            0.01,
+            caption,
+            ha="center",
+            va="bottom",
+            fontsize=8,
+            wrap=True,
+        )
+    output_file = str(output_file)
+    fig.savefig(output_file, bbox_inches="tight", dpi=150)
+    plt.close(fig)
+
+    return output_file