Geostrophic Functions

def bfrq(s, t, p, lat=None):
    """
    Brünt-Väisälä frequency squared and potential vorticity.
    
    Calculates the buoyancy frequency (also known as the Brünt-Väisälä frequency)
    which indicates ocean stability and stratification strength.
    
    Parameters:
    - s: array_like, salinity profile [psu (PSS-78)]
    - t: array_like, temperature profile [°C (ITS-90)]
    - p: array_like, pressure profile [db]
    - lat: array_like, latitude [decimal degrees], optional
    
    Returns:
    - tuple: (n2, q, p_ave)
        - n2: array_like, Brünt-Väisälä frequency squared [rad²/s²]
        - q: array_like, potential vorticity [(m·s)⁻¹]
        - p_ave: array_like, mid-point pressures [db]
    """

def svan(s, t, p=0):
    """
    Specific volume anomaly of seawater.
    
    Calculates the specific volume anomaly which is used in geostrophic
    velocity calculations and dynamic height computations.
    
    Parameters:
    - s: array_like, salinity [psu (PSS-78)]
    - t: array_like, temperature [°C (ITS-90)]
    - p: array_like, pressure [db], default 0
    
    Returns:
    - array_like: specific volume anomaly [m³/kg]
    """

def gpan(s, t, p):
    """
    Geopotential anomaly of seawater.
    
    Calculates geopotential anomaly which represents the work done against
    gravity in lifting a unit mass from a reference level.
    
    Parameters:
    - s: array_like, salinity [psu (PSS-78)]
    - t: array_like, temperature [°C (ITS-90)]
    - p: array_like, pressure [db]
    
    Returns:
    - array_like: geopotential anomaly [m²/s²]
    """

def gvel(ga, lat, lon):
    """
    Geostrophic velocity from geopotential anomaly gradients.
    
    Calculates geostrophic velocities from horizontal gradients of
    geopotential anomaly, accounting for Earth's rotation (Coriolis effect).
    
    Parameters:
    - ga: array_like, geopotential anomaly [m²/s²]
    - lat: array_like, latitude [decimal degrees]
    - lon: array_like, longitude [decimal degrees]
    
    Returns:
    - array_like: geostrophic velocity [m/s]
    """

import seawater as sw
import numpy as np
import matplotlib.pyplot as plt

# Create a stratified ocean profile
pressure = np.arange(0, 1000, 10)  # 0 to 1000 db
# Typical thermocline structure
temperature = 20 * np.exp(-pressure/200) + 2
salinity = 34.5 + 0.5 * (pressure/1000)
latitude = 30.0  # degrees N

# Calculate Brünt-Väisälä frequency
n2, q, p_mid = sw.bfrq(salinity, temperature, pressure, latitude)

# Convert to buoyancy frequency in cycles per hour
n_cph = np.sqrt(np.maximum(n2, 0)) * 3600 / (2 * np.pi)

# Plot stratification
plt.figure(figsize=(10, 6))
plt.subplot(1, 2, 1)
plt.plot(temperature, pressure, 'b-', label='Temperature')
plt.plot(salinity, pressure, 'r-', label='Salinity')
plt.gca().invert_yaxis()
plt.xlabel('Temperature (°C) / Salinity (psu)')
plt.ylabel('Pressure (db)')
plt.legend()
plt.title('Ocean Profile')

plt.subplot(1, 2, 2)
plt.plot(n_cph, p_mid, 'g-', linewidth=2)
plt.gca().invert_yaxis()
plt.xlabel('Buoyancy Frequency (cph)')
plt.ylabel('Pressure (db)')
plt.title('Ocean Stability')
plt.grid(True)

plt.tight_layout()
plt.show()

import seawater as sw
import numpy as np

# Create two ocean stations for geostrophic calculation
pressure = np.arange(0, 500, 25)  # pressure levels

# Station 1 (warmer, less saline)
s1 = np.full_like(pressure, 34.8) - 0.001 * pressure
t1 = 15.0 * np.exp(-pressure/150) + 4

# Station 2 (cooler, more saline)  
s2 = np.full_like(pressure, 35.0) - 0.0005 * pressure
t2 = 12.0 * np.exp(-pressure/180) + 3

# Calculate geopotential anomaly at each station
gpa1 = sw.gpan(s1, t1, pressure)
gpa2 = sw.gpan(s2, t2, pressure)

# Calculate geopotential difference (proxy for geostrophic velocity)
dgpa = gpa2 - gpa1
print("Geopotential anomaly difference:")
for i in range(0, len(pressure), 4):
    print(f"P: {pressure[i]:3.0f} db, Δφ: {dgpa[i]:8.3f} m²/s²")

# Calculate specific volume anomaly
sva1 = sw.svan(s1, t1, pressure)
sva2 = sw.svan(s2, t2, pressure)

print(f"\\nSpecific volume anomaly range:")
print(f"Station 1: {sva1.min():.2e} to {sva1.max():.2e} m³/kg")
print(f"Station 2: {sva2.min():.2e} to {sva2.max():.2e} m³/kg")

import seawater as sw
import numpy as np

# Typical subtropical ocean profile
depth = np.arange(0, 2000, 50)
lat = 25.0  # degrees N

# Convert depth to pressure
pressure = sw.pres(depth, lat)

# Create realistic T-S profile
temperature = 22 * np.exp(-depth/300) + 2.5
salinity = 36.0 - 1.0 * np.exp(-depth/500)

# Calculate dynamic properties
n2, q, p_mid = sw.bfrq(salinity, temperature, pressure, lat)
sva = sw.svan(salinity, temperature, pressure)
gpa = sw.gpan(salinity, temperature, pressure)

# Calculate dynamic height (integrated specific volume anomaly)
# This is a simplified calculation - proper dynamic height requires
# integration between pressure levels
dyn_height = np.cumsum(sva * np.gradient(pressure)) / 10  # convert to dynamic meters

print("Dynamic oceanography summary:")
print(f"Maximum N²: {np.max(n2):.2e} rad²/s²")
print(f"Surface specific volume anomaly: {sva[0]:.2e} m³/kg")
print(f"Deep specific volume anomaly: {sva[-1]:.2e} m³/kg")
print(f"Total dynamic height: {dyn_height[-1]:.3f} dynamic meters")

# Find thermocline depth (maximum stratification)
max_strat_idx = np.argmax(n2)
thermocline_depth = depth[max_strat_idx//2]  # approximate due to mid-point pressures
print(f"Thermocline depth: {thermocline_depth:.0f} m")