2.3 Stability & Convection

#!/usr/bin/env python3
"""
stability_convection.py - Atmospheric stability and CAPE calculation

Run: python3 stability_convection.py
Requires: pip install numpy matplotlib
"""
import numpy as np
import matplotlib.pyplot as plt

# Constants
g = 9.81        # m/s²
R_d = 287.0     # J/(kg·K)
c_p = 1005.0    # J/(kg·K)
L_v = 2.5e6     # J/kg

# Lapse rates
Gamma_d = g / c_p * 1000  # dry adiabatic lapse rate (K/km) ≈ 9.8
Gamma_m = 6.0             # typical moist adiabatic (K/km)

def saturation_vapor_pressure(T_celsius):
    """Tetens formula (hPa)"""
    return 6.11 * np.exp(17.27 * T_celsius / (T_celsius + 237.3))

def parcel_ascent(T_surface, p_surface, Td_surface, z_max=15000, dz=100):
    """
    Lift a parcel and calculate temperatures
    Returns: z, T_parcel, T_env, LCL, LFC, EL
    """
    z = np.arange(0, z_max, dz)
    T_parcel = np.zeros_like(z, dtype=float)
    T_env = np.zeros_like(z, dtype=float)

    # Initial conditions
    T_parcel[0] = T_surface
    T_env[0] = T_surface
    Td = Td_surface

    # Environmental lapse rate (conditionally unstable)
    Gamma_env = 7.0  # K/km

    lcl_z = None
    lfc_z = None
    el_z = None
    saturated = False

    for i in range(1, len(z)):
        # Environmental temperature
        T_env[i] = T_surface - Gamma_env * z[i] / 1000

        # Parcel temperature
        if not saturated:
            # Dry adiabatic ascent
            T_parcel[i] = T_parcel[i-1] - Gamma_d * dz / 1000
            # Check for saturation (rough approximation)
            # LCL: T_parcel ≈ Td
            Td_at_z = Td_surface - 1.8 * z[i] / 1000  # dew point lapse
            if T_parcel[i] <= Td_at_z + 273.15 and lcl_z is None:
                lcl_z = z[i]
                saturated = True
        else:
            # Moist adiabatic ascent
            T_parcel[i] = T_parcel[i-1] - Gamma_m * dz / 1000

        # Find LFC (Level of Free Convection)
        if saturated and T_parcel[i] > T_env[i] and lfc_z is None:
            lfc_z = z[i]

        # Find EL (Equilibrium Level)
        if lfc_z is not None and T_parcel[i] < T_env[i] and el_z is None:
            el_z = z[i]

    return z, T_parcel, T_env, lcl_z, lfc_z, el_z

# Surface conditions
T_sfc = 30 + 273.15  # K
p_sfc = 1013.25      # hPa
Td_sfc = 20          # °C

z, T_parcel, T_env, lcl, lfc, el = parcel_ascent(T_sfc, p_sfc, Td_sfc)

# Calculate CAPE and CIN
cape = 0
cin = 0
for i in range(1, len(z)):
    dz = z[i] - z[i-1]
    buoyancy = g * (T_parcel[i] - T_env[i]) / T_env[i]

    if lfc and z[i] > lfc and (el is None or z[i] < el):
        if buoyancy > 0:
            cape += buoyancy * dz
    elif z[i] < (lfc or 1e10):
        if buoyancy < 0:
            cin += abs(buoyancy) * dz

# Plot
fig, axes = plt.subplots(1, 2, figsize=(14, 8))

# Skew-T style plot
ax1 = axes[0]
ax1.plot(T_env - 273.15, z/1000, 'r-', linewidth=2, label='Environment')
ax1.plot(T_parcel - 273.15, z/1000, 'b-', linewidth=2, label='Parcel')

# Fill CAPE and CIN regions
for i in range(1, len(z)):
    if lfc and z[i] > lfc and (el is None or z[i] < el):
        if T_parcel[i] > T_env[i]:
            ax1.fill_betweenx([z[i-1]/1000, z[i]/1000],
                             [T_env[i-1]-273.15, T_env[i]-273.15],
                             [T_parcel[i-1]-273.15, T_parcel[i]-273.15],
                             color='red', alpha=0.3)

if lcl:
    ax1.axhline(y=lcl/1000, color='green', linestyle='--', label=f'LCL: {lcl/1000:.1f} km')
if lfc:
    ax1.axhline(y=lfc/1000, color='orange', linestyle='--', label=f'LFC: {lfc/1000:.1f} km')
if el:
    ax1.axhline(y=el/1000, color='purple', linestyle='--', label=f'EL: {el/1000:.1f} km')

ax1.set_xlabel('Temperature (°C)', fontsize=12)
ax1.set_ylabel('Height (km)', fontsize=12)
ax1.set_title('Parcel Ascent and Stability', fontsize=14)
ax1.legend(loc='upper right')
ax1.grid(True, alpha=0.3)
ax1.set_xlim(-60, 40)
ax1.set_ylim(0, 15)

# Stability regimes
ax2 = axes[1]
gamma_env_range = np.linspace(4, 12, 100)
colors = []
for g in gamma_env_range:
    if g < Gamma_m:
        colors.append('green')   # Absolutely stable
    elif g > Gamma_d:
        colors.append('red')     # Absolutely unstable
    else:
        colors.append('yellow')  # Conditionally unstable

ax2.scatter(gamma_env_range, np.ones_like(gamma_env_range), c=colors, s=100)
ax2.axvline(x=Gamma_m, color='green', linestyle='--', linewidth=2, label=f'Γ_m = {Gamma_m} K/km')
ax2.axvline(x=Gamma_d, color='red', linestyle='--', linewidth=2, label=f'Γ_d = {Gamma_d:.1f} K/km')
ax2.set_xlabel('Environmental Lapse Rate (K/km)', fontsize=12)
ax2.set_title('Stability Regimes', fontsize=14)
ax2.set_xlim(4, 12)
ax2.set_yticks([])
ax2.legend()

# Add text labels
ax2.text(5, 1.1, 'Absolutely\nStable', ha='center', fontsize=10, color='green')
ax2.text(7.9, 1.1, 'Conditionally\nUnstable', ha='center', fontsize=10, color='orange')
ax2.text(11, 1.1, 'Absolutely\nUnstable', ha='center', fontsize=10, color='red')

plt.tight_layout()
plt.savefig('stability_convection.png', dpi=150, bbox_inches='tight')
plt.show()

print(f"Results:")
print(f"  CAPE: {cape:.0f} J/kg")
print(f"  CIN: {cin:.0f} J/kg")
print(f"  LCL: {lcl/1000 if lcl else 'N/A':.1f} km")
print(f"  LFC: {lfc/1000 if lfc else 'N/A':.1f} km")
print(f"  EL: {el/1000 if el else 'N/A':.1f} km")