# Zonal-symmetric model

This section documents the zonal-symmetric model: its assumptions, required inputs, typical outputs, and an example usage. For theoretical background information, see Paper ... (see References)

---

## Overview and outputs

- **Assumptions:** Zonal symmetry, meridional variations only. The grid is equidistant from 90° N to −90° N. For diagnostic purposes, a stationary zonal periodic pseudo-orographic forcing in the jet region is assumed. The geographical latitude of the forcing is specified as jetlat. In addition, a sponge is used to dampen wave activity that leaves the jet region.
- **Core outputs:**
  - **Basic state plots:** Zonal wind \(u_0(\phi)\), basic-state PV \(q_0(\phi)\), meridional PV gradient \(q_{0y}(\phi)\), and stationary wave number \(\hat{K}_s\).
  - **Resonance plots:** Maximum stationary perturbation streamfunction amplitude as a function of zonal wavenumber.
  - **Solution for one zonal wavenumber (e.g. resonant wavenumber):** spatial distribution of perturbation streamfunction, meridional enstrophy profile, enstrophy tendencies, and meridional Eliassen–Palm (EP) flux.
- **Further reading:** For theory, variable definitions, and parameterization details, please refer to the paper PAPER2.

---

## File structure

- **model_methods.py:** Numerical routines for the model and routines for plotting..
  *Note:* Users do not need to modify this file but can inspect it to understand the implementation.
- **zonal_symmetric_model.py:** User-facing script to define the grid, the basic state wind profile, jet latitude and jet range, and to run diagnostics.

---

## Grid and wind profile requirements

- **Equidistant meridional grid:** Only the meridional resolution is specified; the grid spans from 90°N to −90°N.
- **Wind profile definition:** Provide the basic-state zonal wind \(u_0(\phi)\) as a NumPy array with one value per grid point from 90°N to −90°N. An idealized wind profile is supplied as standard. A Gaussian jet stream is superimposed on a solid-body rotation.
- **Important:**
  ⚠️ **The wind profile must be strictly zonal-symmetric (only depends on latitude) and defined at every grid point on the equidistant meridional grid.**
- **Jet range and sponge:** Within the jet range, sponge damping is off; outside this range, damping increases toward the poles.

---

## Example

```python
import numpy as np
from model_methods import resonanceplot, solve_for_one_wavenumber

"""
Define grid:
The model requires an equidistant grid between the North Pole (90°N)
and the South Pole (−90°N). Only the meridional resolution (Δlat) is specified.
"""
dellat = 0.5  # meridional grid spacing in degrees

"""
Define basic state wind profile u0:
Provide a NumPy array with one value per grid point between 90°N and −90°N.
Example: idealized solid-body rotation with a superimposed Gaussian jet (code in github repository).
"""
# Example placeholder (replace with actual array of length N = 180/dellat + 1):
u0 = np.array([...])

"""
Specify jet characteristics:
- jetlat: latitude of the jet maximum (e.g., 45 for 45°N) where the pseudo-orographic forcing is placed.
- jetrange: latitudinal interval where sponge damping is off
  Outside this range, damping increases toward the poles.
"""
jetlat = 45
jetrange = [0, 60]

# 1) Create basic-state plots for a single zonal wavenumber (e.g., s = 4)
solve_for_one_wavenumber(
    wavenumber_s=4,
    dellat=dellat,
    u0=u0,
    jetlat=jetlat,
    jetrange=jetrange,
    plot_basic_state_plots=True,
    plot_solution_plots=False
)

# 2) Create resonance plot (streamfunction amplitude for several wavenumbers)
s_values, ampl, phase, sfampl, sfphase, Q_res, S_res = resonanceplot(
    dellat=dellat,
    u0=u0,
    jetlat=jetlat,
    jetrange=jetrange,
    s_delta=0.25
)

# 3) Plot spatial distribution, EP flux, and enstrophy diagnostics at resonance
solve_for_one_wavenumber(
    S_res,
    dellat,
    u0,
    jetlat,
    jetrange,
    plot_basic_state_plots=False,
    plot_solution_plots=True
)