Multi-layer Quasi-Geostrophic (Baroclinic) Model

### Formulation
2-layer (or multi-layer) QG equations to capture baroclinic instability.

#### Governing Equations (classic 2-layer):
$$
\begin{align*}
\frac{\partial q_1}{\partial t} + \mathbf{u}_1\cdot\nabla q_1 &= D_1[q_1] + F_1 \\
\frac{\partial q_2}{\partial t} + \mathbf{u}_2\cdot\nabla q_2 &= D_2[q_2] + F_2
\end{align*}
$$
where $q_1$, $q_2$ are PV in each layer, $\mathbf{u}_1$, $\mathbf{u}_2$ from their respective streamfunctions.

#### Forcing
- Wind forcing, bottom drag, thermal damping (as desired).

#### Starter code snippet
```python
# Layer setup (assuming layer variables)
q1 = ...; q2 = ...
psi1 = ...; psi2 = ...

for t in timeline:
dq1dt = -advect(q1, u1, v1) + D1[q1] + F1
dq2dt = -advect(q2, u2, v2) + D2[q2] + F2
# Coupling: update psi from q (inversion)
```

#### Dependencies
- [nlocked by] PR: Multi-layer solver, vertical stacking operators

#### Reference: Vallis, Atmospheric & Oceanic Fluid Dynamics.