10.4 Ocean Data Analysis

Modern oceanography generates massive datasets in formats like netCDF and HDF5, served via OPeNDAP. Analysis techniques span classical statistics (EOF/PCA, spectral analysis, optimal interpolation), time series methods, and increasingly machine learning—all essential for extracting physical insight from observations, reanalyses, and model output.

EOF / PCA Analysis

Empirical Orthogonal Functions (EOFs) decompose a space-time dataset into orthogonal spatial patterns and their temporal coefficients (Principal Components). Given a data matrix$\mathbf{X}$ (time × space), the anomaly covariance matrix is:

$$\mathbf{C} = \frac{1}{n-1} \mathbf{X}'^T \mathbf{X}', \quad \mathbf{C} \mathbf{e}_k = \lambda_k \mathbf{e}_k$$

$\mathbf{X}' = \mathbf{X} - \bar{\mathbf{X}}$ (anomalies),$\mathbf{e}_k$ = EOF pattern, $\lambda_k$ = eigenvalue (variance explained)

The principal component time series is: $\text{PC}_k(t) = \mathbf{X}'(t) \cdot \mathbf{e}_k$. The fraction of total variance explained by mode $k$ is:

$$f_k = \frac{\lambda_k}{\sum_j \lambda_j} \times 100\%$$

EOF1 of tropical Pacific SST captures ENSO (~50% of variance)

Spectral & Wavelet Analysis

Fourier spectral analysis identifies periodicities in ocean time series (tidal harmonics, seasonal cycles, ENSO). The discrete Fourier transform of a time series $x(t_n)$ is:

$$X(f_k) = \sum_{n=0}^{N-1} x(t_n) \, e^{-2\pi i f_k t_n}, \quad S(f_k) = \frac{2|\hat{X}(f_k)|^2}{N \Delta t}$$

$S(f_k)$ = power spectral density; peaks indicate dominant periodicities

For non-stationary signals (e.g., ENSO whose period and amplitude vary over time), wavelet analysis provides time-frequency localization. The continuous wavelet transform using a Morlet mother wavelet $\psi_0(\eta) = \pi^{-1/4} e^{i\omega_0 \eta} e^{-\eta^2/2}$ is:

$$W(s, \tau) = \frac{1}{\sqrt{s}} \sum_{n=0}^{N-1} x(t_n) \, \psi^*\!\left(\frac{t_n - \tau}{s}\right)$$

$s$ = scale (related to period), $\tau$ = time shift

Optimal Interpolation (Objective Analysis)

Optimal interpolation (OI) maps irregularly-spaced ocean observations onto a regular grid by minimizing the expected mean-squared error. The analysis at grid point $\mathbf{x}_0$ is:

$$\hat{\phi}(\mathbf{x}_0) = \bar{\phi} + \mathbf{c}^T (\mathbf{C} + \epsilon^2 \mathbf{I})^{-1} (\mathbf{d} - \bar{\phi})$$

$\mathbf{c}$ = covariance between grid point and observations,$\mathbf{C}$ = observation-observation covariance matrix,$\epsilon^2$ = noise variance

The background covariance is typically modeled as a Gaussian function of distance:$C(r) = \sigma^2 \exp(-r^2 / (2L^2))$, where $L$ is the decorrelation length scale (~100–300 km for mesoscale SST). The analysis error variance is:

$$\sigma_a^2 = \sigma^2 - \mathbf{c}^T (\mathbf{C} + \epsilon^2 \mathbf{I})^{-1} \mathbf{c}$$

Error is reduced where observations are dense and nearby

Derivation: EOF/PCA from the Covariance Matrix

Step 1: Form the Anomaly Matrix

Given a data matrix $\mathbf{X}$ of dimensions $n \times p$ ($n$ time steps, $p$ spatial locations), remove the temporal mean at each location to form the anomaly matrix:

$$\mathbf{X}' = \mathbf{X} - \mathbf{1}\bar{\mathbf{x}}^T, \quad \bar{x}_j = \frac{1}{n}\sum_{i=1}^n X_{ij}$$

Step 2: Compute the Covariance Matrix

The $p \times p$ spatial covariance matrix is:

$$\mathbf{C} = \frac{1}{n-1}\mathbf{X}'^T\mathbf{X}'$$

$C_{jk}$ gives the covariance between locations $j$ and $k$. The matrix is symmetric and positive semi-definite.

Step 3: Eigenvalue Problem

Find the eigenvectors $\mathbf{e}_k$ (EOF patterns) and eigenvalues $\lambda_k$ of $\mathbf{C}$:

$$\mathbf{C}\mathbf{e}_k = \lambda_k \mathbf{e}_k, \quad \mathbf{e}_k^T \mathbf{e}_l = \delta_{kl}$$

The EOFs are orthonormal: each represents an independent spatial pattern of variability.

Step 4: Optimality Property

EOF1 is the unit vector that maximizes the projected variance. Using the method of Lagrange multipliers with constraint $\|\mathbf{e}\| = 1$:

$$\max_{\mathbf{e}} \; \mathbf{e}^T \mathbf{C} \mathbf{e} \quad \text{s.t.} \quad \mathbf{e}^T\mathbf{e} = 1$$

$$\nabla_{\mathbf{e}}\left[\mathbf{e}^T\mathbf{C}\mathbf{e} - \lambda(\mathbf{e}^T\mathbf{e} - 1)\right] = 2\mathbf{C}\mathbf{e} - 2\lambda\mathbf{e} = 0 \quad \Rightarrow \quad \mathbf{C}\mathbf{e} = \lambda\mathbf{e}$$

Step 5: Principal Component Time Series

Project the anomaly data onto EOF pattern $k$ to obtain the $k$-th principal component (PC) time series:

$$\text{PC}_k(t) = \mathbf{X}'(t,:) \cdot \mathbf{e}_k = \sum_{j=1}^{p} X'_{tj} \, e_{k,j}$$

The PCs are uncorrelated: $\text{Cov}(\text{PC}_k, \text{PC}_l) = \lambda_k \delta_{kl}$. The variance of PC$_k$ equals $\lambda_k$.

Step 6: Variance Explained and Truncation

The fraction of total variance captured by mode $k$ is:

$$f_k = \frac{\lambda_k}{\sum_{j=1}^{p}\lambda_j} = \frac{\lambda_k}{\text{tr}(\mathbf{C})}$$

The data can be reconstructed from $m \ll p$ modes: $\mathbf{X}' \approx \sum_{k=1}^{m} \text{PC}_k \, \mathbf{e}_k^T$. The North et al. (1982) rule of thumb tests whether eigenvalues are well-separated: $\delta\lambda_k \approx \lambda_k \sqrt{2/n}$.

Step 7: SVD Equivalence (Practical Computation)

For large datasets, the SVD of the anomaly matrix is more efficient than eigendecomposition of $\mathbf{C}$:

$$\mathbf{X}' = \mathbf{U}\mathbf{\Sigma}\mathbf{V}^T, \quad \mathbf{e}_k = \mathbf{V}_{:,k}, \quad \text{PC}_k = \mathbf{U}_{:,k}\sigma_k, \quad \lambda_k = \frac{\sigma_k^2}{n-1}$$

Derivation: Spectral Analysis & the Periodogram

Step 1: Discrete Fourier Transform

Given a uniformly sampled time series $x_n = x(n\Delta t)$ for $n = 0, 1, \ldots, N-1$, the discrete Fourier transform at frequency $f_k = k/(N\Delta t)$ is:

$$\hat{X}(f_k) = \sum_{n=0}^{N-1} x_n \, e^{-2\pi i k n / N}$$

Step 2: Parseval's Theorem (Energy Conservation)

The total energy is preserved between time and frequency domains:

$$\sum_{n=0}^{N-1} |x_n|^2 = \frac{1}{N}\sum_{k=0}^{N-1}|\hat{X}(f_k)|^2$$

Step 3: Define the Periodogram

The periodogram is the sample estimate of the power spectral density (PSD). It distributes the total variance across frequency bins:

$$S(f_k) = \frac{2\Delta t}{N}|\hat{X}(f_k)|^2 = \frac{2}{N\Delta t}\left|\sum_{n=0}^{N-1}x_n \, e^{-2\pi i f_k n\Delta t}\right|^2$$

The factor of 2 accounts for folding negative frequencies onto positive ones (one-sided spectrum). Units are $[\text{unit}]^2/\text{Hz}$.

Step 4: Frequency Resolution and Nyquist Limit

The fundamental frequency resolution is $\Delta f = 1/(N\Delta t)$ (the Rayleigh resolution). The maximum resolvable frequency (Nyquist) is:

$$f_{\text{Nyquist}} = \frac{1}{2\Delta t}$$

For monthly ocean data with $\Delta t = 1$ month: $f_{\text{Ny}} = 0.5$ cycle/month, and $N = 300$ months gives $\Delta f = 1/300$ cycle/month.

Step 5: Statistical Properties and Smoothing

The raw periodogram is an inconsistent estimator: its variance does not decrease with more data. Each periodogram ordinate follows a $\chi^2_2$ distribution. Smoothing (Welch's method, Bartlett averaging, or multitaper) reduces variance:

$$\hat{S}_{\text{smooth}}(f) = \frac{1}{M}\sum_{m=1}^{M} S_m(f), \quad \text{d.o.f.} = 2M$$

Confidence intervals for the smoothed PSD: $\frac{\nu \hat{S}(f)}{\chi^2_{\nu,\alpha/2}} \leq S(f) \leq \frac{\nu \hat{S}(f)}{\chi^2_{\nu,1-\alpha/2}}$ where $\nu = 2M$ degrees of freedom.

Step 6: Red Noise Background and Significance Testing

Ocean time series typically have red spectra (more variance at low frequencies). The AR(1) red noise null spectrum for a process with lag-1 autocorrelation $r_1$ is:

$$S_{\text{red}}(f) = \frac{S_0 (1 - r_1^2)}{1 - 2r_1\cos(2\pi f\Delta t) + r_1^2}$$

Spectral peaks exceeding the 95% confidence level above this background are considered statistically significant periodicities.

Derivation: Objective Analysis / Optimal Interpolation

Step 1: Statement of the Problem

We have $m$ observations $d_i = \phi(\mathbf{x}_i) + \epsilon_i$ of a field $\phi$ at irregular locations, where $\epsilon_i$ is measurement noise with variance $\epsilon^2$. We seek the best linear unbiased estimate at an arbitrary point $\mathbf{x}_0$:

$$\hat{\phi}(\mathbf{x}_0) = \bar{\phi} + \sum_{i=1}^{m} w_i (d_i - \bar{\phi})$$

Step 2: Minimize the Expected Mean-Squared Error

Define the analysis error as $e = \hat{\phi}(\mathbf{x}_0) - \phi(\mathbf{x}_0)$. Minimize the expected squared error $\langle e^2 \rangle$ with respect to the weights $w_i$:

$$\frac{\partial}{\partial w_j}\left\langle\left[\sum_i w_i(d_i - \bar{\phi}) - (\phi_0 - \bar{\phi})\right]^2\right\rangle = 0$$

Step 3: Expand Using Covariances

Define the signal covariance between two points: $C(\mathbf{x}_i, \mathbf{x}_j) = \langle\phi'(\mathbf{x}_i)\phi'(\mathbf{x}_j)\rangle$. The observation covariance matrix includes noise: $[\mathbf{C}_{obs}]_{ij} = C(\mathbf{x}_i, \mathbf{x}_j) + \epsilon^2\delta_{ij}$. Taking the derivative yields:

$$\sum_{i=1}^{m}\left[C(\mathbf{x}_j, \mathbf{x}_i) + \epsilon^2\delta_{ji}\right]w_i = C(\mathbf{x}_j, \mathbf{x}_0) \quad \forall j$$

Step 4: Matrix Form of the Solution

In matrix notation, $(\mathbf{C} + \epsilon^2\mathbf{I})\mathbf{w} = \mathbf{c}$, giving:

$$\mathbf{w} = (\mathbf{C} + \epsilon^2\mathbf{I})^{-1}\mathbf{c}$$

$$\hat{\phi}(\mathbf{x}_0) = \bar{\phi} + \mathbf{c}^T(\mathbf{C} + \epsilon^2\mathbf{I})^{-1}(\mathbf{d} - \bar{\phi})$$

Step 5: Analysis Error Variance

Substituting the optimal weights back into the error expression gives the minimum achievable error variance:

$$\sigma_a^2(\mathbf{x}_0) = \sigma^2 - \mathbf{c}^T(\mathbf{C} + \epsilon^2\mathbf{I})^{-1}\mathbf{c}$$

When $\mathbf{x}_0$ coincides with an observation, the error reduces to $\sigma_a^2 \to \epsilon^2 \sigma^2/(\sigma^2 + \epsilon^2)$. Far from any observation, $\mathbf{c} \to 0$ and $\sigma_a^2 \to \sigma^2$ (background variance).

Step 6: Choice of Covariance Model

The background covariance function must be specified a priori. Common choices for isotropic, homogeneous fields:

$$C(r) = \sigma^2 \exp\!\left(-\frac{r^2}{2L^2}\right) \quad \text{(Gaussian)}, \qquad C(r) = \sigma^2\left(1 + \frac{r}{L}\right)e^{-r/L} \quad \text{(SOAR)}$$

The decorrelation length $L$ controls the smoothness of the analysis. For SST: $L \sim 100\text{--}300$ km; for SSH: $L \sim 100\text{--}200$ km. Anisotropic models use different length scales along and across mean flow.

Step 7: Connection to Kriging and Kalman Filtering

OI is mathematically identical to simple Kriging in geostatistics and to the analysis step of the Kalman filter. The Kalman gain matrix is:

$$\mathbf{K} = \mathbf{P}^f\mathbf{H}^T(\mathbf{H}\mathbf{P}^f\mathbf{H}^T + \mathbf{R})^{-1}$$

where $\mathbf{P}^f$ is the forecast error covariance, $\mathbf{H}$ is the observation operator, and $\mathbf{R}$ is the observation error covariance. In OI, $\mathbf{P}^f$ is replaced by the static background covariance.

Machine Learning in Oceanography

Machine learning has transformed oceanographic data analysis, enabling pattern recognition in massive satellite datasets, gap-filling of sparse in-situ observations, and emulation of computationally expensive numerical models. The key techniques span supervised learning (regression, classification), unsupervised learning (clustering, dimensionality reduction), and physics-informed approaches that combine data-driven methods with domain knowledge.

Neural Networks for SST Prediction

Sea surface temperature forecasting is one of the most successful applications of deep learning in oceanography. Traditional approaches include persistence (tomorrow's SST = today's SST), linear inverse models (LIM), and analog forecasting. Neural networks outperform these methods, particularly for ENSO prediction at lead times of 6–18 months.

LSTM (Long Short-Term Memory): Recurrent neural network designed for sequential data. The cell state carries information over long time lags, making LSTMs ideal for capturing ENSO's multi-year memory. The gating mechanism (forget gate, input gate, output gate) selectively retains or discards information:

$$f_t = \sigma(W_f \cdot [h_{t-1}, x_t] + b_f), \quad i_t = \sigma(W_i \cdot [h_{t-1}, x_t] + b_i)$$

$$C_t = f_t \odot C_{t-1} + i_t \odot \tanh(W_C \cdot [h_{t-1}, x_t] + b_C)$$

ConvLSTM: Extends LSTM by replacing matrix multiplications with convolutions, preserving spatial structure. Input is a spatiotemporal field (e.g., monthly SST maps). The convolution kernels capture spatial patterns (warm pool, cold tongue) while the LSTM architecture captures temporal evolution. Ham et al. (2019) demonstrated ConvLSTM ENSO forecasts with correlation skill >0.5 at 17-month lead time, surpassing all dynamical models.

Transfer Learning: Pre-train on CMIP6 climate model output (abundant data), then fine-tune on short observational record (1982–present). This overcomes the limited training data problem that plagues purely data-driven ocean forecasting. The model learns physics-like relationships from simulations and adapts to real-world biases.

Water Mass Classification

Traditional water mass identification relies on expert analysis of T-S diagrams. Machine learning automates this process and reveals structure invisible to manual inspection:

K-means Clustering: Partitions T-S-depth-nutrient profiles into K clusters by minimizing within-cluster variance. The objective function is:

$$J = \sum_{k=1}^{K} \sum_{i \in C_k} \|\mathbf{x}_i - \boldsymbol{\mu}_k\|^2$$

Applied to Argo float profiles, K-means identifies water masses (e.g., AAIW, NADW, AABW) and their boundaries more objectively than manual methods. The silhouette score or BIC helps determine optimal K.

Gaussian Mixture Models (GMM): Soft clustering — each profile has a probability of belonging to each water mass. The EM algorithm fits a mixture of K Gaussians. Captures overlapping water masses and mixing zones that hard clustering misses. Probability maps show transition zones between water masses.

Self-Organizing Maps (SOM): Neural network that projects high-dimensional T/S/nutrient profiles onto a 2D map while preserving topological relationships. Adjacent neurons represent similar water masses. SOMs have been used to classify biogeochemical provinces in the Southern Ocean and to track water mass transformation along Argo float trajectories.

Physics-Informed Neural Networks (PINNs)

PINNs embed physical laws directly into the neural network loss function, combining the flexibility of neural networks with the constraints of known physics:

$$\mathcal{L} = \underbrace{\frac{1}{N_d}\sum_{i=1}^{N_d}|u_\theta(\mathbf{x}_i) - u_i^{\text{obs}}|^2}_{\text{data loss}} + \underbrace{\lambda \frac{1}{N_r}\sum_{j=1}^{N_r}|\mathcal{N}[u_\theta](\mathbf{x}_j)|^2}_{\text{physics loss}}$$

where $\mathcal{N}$ is a differential operator encoding the PDE (e.g., advection-diffusion, shallow water equations) and $\lambda$ balances data fidelity against physical consistency.

Ocean applications: Inferring subsurface temperature from surface observations and sparse Argo profiles; reconstructing velocity fields from altimetry while satisfying geostrophic balance; estimating turbulent diffusivity from tracer observations while obeying the advection-diffusion equation.

Neural Operator approaches: Fourier Neural Operators (FNO) and DeepONet learn the solution operator of PDEs, enabling 1000× faster ocean model emulation. Trained on NEMO or MOM6 output, they can predict ocean state evolution in milliseconds instead of hours.

Ocean Reanalyses and Data Assimilation

Ocean reanalyses combine numerical ocean models with observations through data assimilation algorithms to produce the best estimate of the 3D ocean state over time:

4D-Var (Variational): Minimizes the cost function $J(\mathbf{x}_0) = \frac{1}{2}(\mathbf{x}_0 - \mathbf{x}_b)^T \mathbf{B}^{-1}(\mathbf{x}_0 - \mathbf{x}_b) + \frac{1}{2}\sum_k (\mathbf{y}_k - H_k[\mathbf{x}_k])^T \mathbf{R}_k^{-1}(\mathbf{y}_k - H_k[\mathbf{x}_k])$ over a time window, requiring the adjoint model for gradient computation. Used in ECMWF's ORAS5.

Ensemble Kalman Filter (EnKF): Propagates an ensemble of model states forward, then updates each member using observations. The ensemble spread estimates the forecast error covariance. Used in GLORYS12 (Mercator Océan). Advantage: no adjoint model needed; handles nonlinearity better.

Key products: ORAS5 (ECMWF, 1958–present, 0.25°), GLORYS12 (Mercator, 1993–present, 1/12°), SODA (UMD, 1980–present, 0.25°). These provide global 3D temperature, salinity, currents, and sea level fields used for climate monitoring, initialization of seasonal forecasts, and validation of ocean models.

Python: Water Mass Classification & SST Prediction

Python

K-means clustering, GMM, and LSTM-like SST forecasting

script.py241 lines

#!/usr/bin/env python3
"""ml_oceanography.py - Water mass classification & SST prediction"""
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# ============================================================
# 1. WATER MASS CLASSIFICATION (K-means on T-S profiles)
# ============================================================
# Generate synthetic T-S profiles for 5 water masses
water_masses = {
    'AAIW': {'T': (3.5, 0.8), 'S': (34.3, 0.15)},   # Antarctic Intermediate
    'NADW': {'T': (2.5, 0.6), 'S': (34.9, 0.1)},     # North Atlantic Deep
    'AABW': {'T': (0.5, 0.3), 'S': (34.65, 0.08)},   # Antarctic Bottom
    'MW':   {'T': (11.0, 1.0), 'S': (35.7, 0.2)},    # Mediterranean
    'SACW': {'T': (16.0, 2.0), 'S': (35.5, 0.3)},    # South Atlantic Central
}

n_per_mass = 80
T_all, S_all, labels_true = [], [], []
for idx, (name, params) in enumerate(water_masses.items()):
    T = np.random.normal(params['T'][0], params['T'][1], n_per_mass)
    S = np.random.normal(params['S'][0], params['S'][1], n_per_mass)
    T_all.extend(T); S_all.extend(S)
    labels_true.extend([idx] * n_per_mass)

T_all = np.array(T_all)
S_all = np.array(S_all)
labels_true = np.array(labels_true)

# --- K-means from scratch ---
def kmeans(X, K, max_iter=100):
    """K-means clustering algorithm."""
    n = len(X)
    # Random initialization
    idx = np.random.choice(n, K, replace=False)
    centroids = X[idx].copy()
    labels = np.zeros(n, dtype=int)

for iteration in range(max_iter):
        # Assign to nearest centroid
        for i in range(n):
            dists = np.sum((centroids - X[i])**2, axis=1)
            labels[i] = np.argmin(dists)
        # Update centroids
        new_centroids = np.zeros_like(centroids)
        for k in range(K):
            mask = labels == k
            if np.any(mask):
                new_centroids[k] = X[mask].mean(axis=0)
            else:
                new_centroids[k] = centroids[k]
        if np.allclose(centroids, new_centroids):
            break
        centroids = new_centroids

# Compute inertia
    inertia = sum(np.sum((X[labels == k] - centroids[k])**2) for k in range(K))
    return labels, centroids, inertia

# Normalize features
T_norm = (T_all - T_all.mean()) / T_all.std()
S_norm = (S_all - S_all.mean()) / S_all.std()
X = np.column_stack([T_norm, S_norm])

# Run K-means with K=5
labels_km, centroids, inertia = kmeans(X, K=5)

# Elbow method: try K = 2..8
inertias = []
for K in range(2, 9):
    _, _, iner = kmeans(X, K)
    inertias.append(iner)

fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# T-S diagram with true labels
colors = ['#e74c3c', '#3498db', '#2ecc71', '#f39c12', '#9b59b6']
wm_names = list(water_masses.keys())
for idx in range(5):
    mask = labels_true == idx
    axes[0].scatter(S_all[mask], T_all[mask], c=colors[idx],
                    s=15, alpha=0.6, label=wm_names[idx])
axes[0].set_xlabel('Salinity (PSU)'); axes[0].set_ylabel('Temperature (C)')
axes[0].set_title('T-S Diagram: True Water Masses')
axes[0].legend(fontsize=8); axes[0].grid(True, alpha=0.3)

# K-means result
for k in range(5):
    mask = labels_km == k
    axes[1].scatter(S_all[mask], T_all[mask], s=15, alpha=0.6,
                    label=f'Cluster {k+1}')
# Denormalize centroids for plotting
c_S = centroids[:, 1] * S_all.std() + S_all.mean()
c_T = centroids[:, 0] * T_all.std() + T_all.mean()
axes[1].scatter(c_S, c_T, c='red', marker='X', s=200, edgecolors='black',
                linewidths=2, zorder=5, label='Centroids')
axes[1].set_xlabel('Salinity (PSU)'); axes[1].set_ylabel('Temperature (C)')
axes[1].set_title('K-means Clustering (K=5)')
axes[1].legend(fontsize=8); axes[1].grid(True, alpha=0.3)

# Elbow plot
axes[2].plot(range(2, 9), inertias, 'bo-', lw=2)
axes[2].axvline(5, color='red', ls='--', label='Optimal K=5')
axes[2].set_xlabel('Number of clusters K')
axes[2].set_ylabel('Inertia (within-cluster SS)')
axes[2].set_title('Elbow Method')
axes[2].legend(); axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('water_mass_classification.png', dpi=150)
plt.show()

# ============================================================
# 2. SST PREDICTION (Simple neural network from scratch)
# ============================================================
# Generate synthetic Nino 3.4 monthly SST anomaly
t = np.arange(600)  # 50 years monthly
sst = (1.5 * np.sin(2*np.pi*t/42) *    # quasi-periodic ENSO (~3.5 yr)
       (1 + 0.4*np.sin(2*np.pi*t/144)) + # amplitude modulation
       0.3 * np.sin(2*np.pi*t/12) +      # annual cycle
       0.4 * np.random.randn(len(t)))     # noise

# Create input/output pairs: predict SST(t+lead) from SST(t-11:t)
window = 12  # 12 months lookback
lead = 6     # 6-month forecast
X_seq, y_seq = [], []
for i in range(window, len(sst) - lead):
    X_seq.append(sst[i-window:i])
    y_seq.append(sst[i + lead])
X_seq = np.array(X_seq)
y_seq = np.array(y_seq)

# Normalize
X_mean, X_std = X_seq.mean(), X_seq.std()
X_seq = (X_seq - X_mean) / X_std
y_seq_norm = (y_seq - X_mean) / X_std

# Train/test split (80/20)
split = int(0.8 * len(X_seq))
X_train, X_test = X_seq[:split], X_seq[split:]
y_train, y_test = y_seq_norm[:split], y_seq_norm[split:]

# --- Simple 2-layer neural network ---
class SimpleNN:
    def __init__(self, n_in, n_hidden, n_out, lr=0.001):
        self.W1 = np.random.randn(n_in, n_hidden) * 0.1
        self.b1 = np.zeros(n_hidden)
        self.W2 = np.random.randn(n_hidden, n_out) * 0.1
        self.b2 = np.zeros(n_out)
        self.lr = lr

def forward(self, X):
        self.z1 = X @ self.W1 + self.b1
        self.a1 = np.tanh(self.z1)
        self.z2 = self.a1 @ self.W2 + self.b2
        return self.z2.flatten()

def backward(self, X, y, y_pred):
        n = len(y)
        dz2 = (y_pred - y).reshape(-1, 1) / n
        dW2 = self.a1.T @ dz2
        db2 = dz2.sum(axis=0)
        da1 = dz2 @ self.W2.T
        dz1 = da1 * (1 - self.a1**2)
        dW1 = X.T @ dz1
        db1 = dz1.sum(axis=0)
        # Update
        self.W1 -= self.lr * dW1
        self.b1 -= self.lr * db1
        self.W2 -= self.lr * dW2
        self.b2 -= self.lr * db2

def train(self, X, y, epochs=200, batch_size=32):
        losses = []
        for ep in range(epochs):
            # Mini-batch SGD
            idx = np.random.permutation(len(X))
            for start in range(0, len(X), batch_size):
                batch = idx[start:start+batch_size]
                y_pred = self.forward(X[batch])
                self.backward(X[batch], y[batch], y_pred)
            # Epoch loss
            y_full = self.forward(X)
            loss = np.mean((y_full - y)**2)
            losses.append(loss)
        return losses

nn = SimpleNN(window, 32, 1, lr=0.005)
losses = nn.train(X_train, y_train, epochs=300)

# Predictions
y_pred_train = nn.forward(X_train) * X_std + X_mean
y_pred_test = nn.forward(X_test) * X_std + X_mean
y_true_test = y_test * X_std + X_mean

# Persistence forecast (baseline)
y_persist = X_test[:, -1] * X_std + X_mean  # last known value

# Correlation skill
corr_nn = np.corrcoef(y_pred_test, y_true_test)[0, 1]
corr_persist = np.corrcoef(y_persist, y_true_test)[0, 1]
rmse_nn = np.sqrt(np.mean((y_pred_test - y_true_test)**2))
rmse_persist = np.sqrt(np.mean((y_persist - y_true_test)**2))

fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# Training loss
axes[0].semilogy(losses, 'b-', lw=1)
axes[0].set_xlabel('Epoch'); axes[0].set_ylabel('MSE Loss')
axes[0].set_title('Training Loss'); axes[0].grid(True, alpha=0.3)

# Test predictions
t_test = np.arange(len(y_true_test))
axes[1].plot(t_test, y_true_test, 'k-', lw=1.5, label='Observed')
axes[1].plot(t_test, y_pred_test, 'r-', lw=1.5, alpha=0.8,
             label=f'NN (r={corr_nn:.3f})')
axes[1].plot(t_test, y_persist, 'g--', lw=1, alpha=0.6,
             label=f'Persist (r={corr_persist:.3f})')
axes[1].set_xlabel('Month'); axes[1].set_ylabel('SST anomaly (K)')
axes[1].set_title(f'6-Month SST Forecast (RMSE: NN={rmse_nn:.2f}, Persist={rmse_persist:.2f})')
axes[1].legend(); axes[1].grid(True, alpha=0.3)

# Scatter plot
axes[2].scatter(y_true_test, y_pred_test, s=10, alpha=0.5, c='blue')
lims = [min(y_true_test.min(), y_pred_test.min()),
        max(y_true_test.max(), y_pred_test.max())]
axes[2].plot(lims, lims, 'r--', lw=2)
axes[2].set_xlabel('Observed SST anomaly')
axes[2].set_ylabel('Predicted SST anomaly')
axes[2].set_title(f'Scatter (r = {corr_nn:.3f})')
axes[2].set_aspect('equal'); axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('sst_prediction.png', dpi=150)
plt.show()

print(f"Water mass classification: K=5, inertia={inertia:.1f}")
print(f"SST NN forecast: r={corr_nn:.3f}, RMSE={rmse_nn:.3f}")
print(f"Persistence:     r={corr_persist:.3f}, RMSE={rmse_persist:.3f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran: Physics-Informed Temperature Interpolation

This Fortran program demonstrates a physics-informed approach to ocean temperature interpolation. Instead of purely statistical interpolation, it solves the steady-state advection-diffusion equation constrained by sparse observations — the core idea behind PINNs, implemented as a relaxation solver.

Fortran: Physics-Informed Temperature Interpolation

Fortran

Advection-diffusion constrained interpolation

program.f90135 lines

program physics_informed_interp
  ! Physics-informed ocean temperature interpolation
  ! Solves steady advection-diffusion nudged toward sparse observations
  ! T_xx + T_yy - Pe*(u*T_x + v*T_y) = -lambda*(T - T_obs) at obs points
  implicit none
  integer, parameter :: nx = 50, ny = 50, n_obs = 25
  real(8) :: T(nx,ny), T_new(nx,ny), T_obs(n_obs)
  real(8) :: x_obs(n_obs), y_obs(n_obs)
  integer :: ix_obs(n_obs), iy_obs(n_obs)
  real(8) :: dx, dy, Lx, Ly, kappa, Pe
  real(8) :: u_vel, v_vel, lam, residual
  real(8) :: T_xx, T_yy, T_x, T_y, nudge
  integer :: i, j, k, iter, max_iter
  logical :: is_obs

! Domain: 1000 km x 1000 km ocean region
  Lx = 1000.0d0; Ly = 1000.0d0
  dx = Lx / dble(nx-1)
  dy = Ly / dble(ny-1)
  kappa = 500.0d0   ! diffusivity (m^2/s equivalent in grid units)
  Pe = 0.5d0        ! Peclet number (advection/diffusion ratio)
  u_vel = 0.3d0     ! eastward current
  v_vel = -0.1d0    ! southward current
  lam = 10.0d0      ! observation nudging strength
  max_iter = 5000

! --- Generate synthetic observations ---
  ! True field: warm in west, cold in east with a front
  do k = 1, n_obs
    x_obs(k) = Lx * mod(dble(k*37), 47.0d0) / 47.0d0
    y_obs(k) = Ly * mod(dble(k*53), 43.0d0) / 43.0d0
    ix_obs(k) = max(2, min(nx-1, nint(x_obs(k)/dx) + 1))
    iy_obs(k) = max(2, min(ny-1, nint(y_obs(k)/dy) + 1))
    ! True temperature: warm pool + thermal front
    T_obs(k) = 28.0d0 - 8.0d0 * x_obs(k)/Lx &
               + 3.0d0 * sin(3.14159d0 * y_obs(k)/Ly) &
               + 0.5d0 * sin(2.0d0*3.14159d0*x_obs(k)/Lx)
    ! Add measurement noise
    T_obs(k) = T_obs(k) + 0.3d0 * sin(dble(k*17))
  end do

! Initial guess: mean of observations
  T = sum(T_obs) / dble(n_obs)

! Boundary conditions
  do j = 1, ny
    T(1,j)  = 28.0d0  ! warm western boundary
    T(nx,j) = 20.0d0  ! cold eastern boundary
  end do
  do i = 1, nx
    T(i,1)  = T(i,1); T(i,ny) = T(i,ny)  ! Neumann (zero gradient)
  end do

write(*,'(A)') '=== Physics-Informed Temperature Interpolation ==='
  write(*,'(A,I4,A,I4)') 'Grid: ', nx, ' x ', ny
  write(*,'(A,I4)') 'Observations: ', n_obs
  write(*,'(A,F4.1)') 'Peclet number: ', Pe
  write(*,*)
  write(*,'(A)') 'Iter    Max_Residual    Mean_T     Min_T      Max_T'

! --- Iterative relaxation (SOR with physics) ---
  do iter = 1, max_iter
    residual = 0.0d0
    T_new = T

do j = 2, ny-1
      do i = 2, nx-1
        ! Diffusion (Laplacian)
        T_xx = (T(i+1,j) - 2.0d0*T(i,j) + T(i-1,j)) / (dx*dx)
        T_yy = (T(i,j+1) - 2.0d0*T(i,j) + T(i,j-1)) / (dy*dy)

! Advection (upwind)
        if (u_vel > 0.0d0) then
          T_x = (T(i,j) - T(i-1,j)) / dx
        else
          T_x = (T(i+1,j) - T(i,j)) / dx
        end if
        if (v_vel > 0.0d0) then
          T_y = (T(i,j) - T(i,j-1)) / dy
        else
          T_y = (T(i,j+1) - T(i,j)) / dy
        end if

! Observation nudging
        nudge = 0.0d0
        do k = 1, n_obs
          if (ix_obs(k) == i .and. iy_obs(k) == j) then
            nudge = lam * (T_obs(k) - T(i,j))
            exit
          end if
        end do

! Update: kappa*Laplacian - Pe*(u*Tx + v*Ty) + nudge = 0
        T_new(i,j) = T(i,j) + 0.2d0 * ( &
            kappa * (T_xx + T_yy) - Pe*(u_vel*T_x + v_vel*T_y) + nudge )

residual = max(residual, abs(T_new(i,j) - T(i,j)))
      end do
    end do

! Neumann BC top/bottom
    T_new(:,1) = T_new(:,2)
    T_new(:,ny) = T_new(:,ny-1)

T = T_new

if (mod(iter, 500) == 0 .or. iter == 1) then
      write(*,'(I5, E16.6, 3F11.3)') iter, residual, &
            sum(T)/(nx*ny), minval(T), maxval(T)
    end if

if (residual < 1.0d-8) exit
  end do

! --- Output final field cross-sections ---
  write(*,*)
  write(*,'(A)') 'Zonal cross-section at mid-latitude (y=500km):'
  write(*,'(A)') '  x(km)    T_interp(C)    T_true(C)'
  do i = 1, nx, 5
    write(*,'(F8.1, 2F14.3)') dble(i-1)*dx, T(i, ny/2), &
        28.0d0 - 8.0d0*dble(i-1)*dx/Lx + &
        3.0d0*sin(3.14159d0*0.5d0) + &
        0.5d0*sin(2.0d0*3.14159d0*dble(i-1)*dx/Lx)

end do

! Compute RMSE at observation points
  residual = 0.0d0
  do k = 1, n_obs
    residual = residual + (T(ix_obs(k), iy_obs(k)) - T_obs(k))**2
  end do
  write(*,*)
  write(*,'(A,F8.4)') 'RMSE at obs points: ', sqrt(residual/dble(n_obs))
  write(*,'(A)') 'Physics constraint ensures smooth, dynamically consistent field'
end program physics_informed_interp

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Python: EOF Analysis, Wavelet Transform & Optimal Interpolation

Python

!/usr/bin/env python3

script.py150 lines

#!/usr/bin/env python3
"""data_analysis.py - EOF, wavelet, optimal interpolation, ML classification"""
import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg

# === 1. EOF analysis of synthetic SST data ===
np.random.seed(42)
nt, nx = 300, 60
t = np.arange(nt)
x = np.arange(nx)

# Mode 1: ENSO-like (quasi-periodic ~36 months)
e1_spatial = np.sin(np.pi * x / nx)
e1_temporal = np.sin(2*np.pi*t/36) + 0.2*np.random.randn(nt)

# Mode 2: Annual cycle
e2_spatial = np.cos(2*np.pi*x/nx)
e2_temporal = np.cos(2*np.pi*t/12)

# Mode 3: Trend
e3_spatial = np.ones(nx) * 0.5
e3_temporal = t / nt

data = (2*np.outer(e1_temporal, e1_spatial) +
        1.2*np.outer(e2_temporal, e2_spatial) +
        0.5*np.outer(e3_temporal, e3_spatial) +
        0.4*np.random.randn(nt, nx))

# Remove mean, compute EOFs
data_anom = data - data.mean(axis=0)
cov = np.cov(data_anom.T)
eigenvalues, eigenvectors = linalg.eigh(cov)
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
pcs = data_anom @ eigenvectors
var_exp = eigenvalues / eigenvalues.sum() * 100

fig, axes = plt.subplots(2, 3, figsize=(16, 10))

# EOF patterns
for k in range(3):
    axes[0, k].plot(x, eigenvectors[:, k], 'b-', lw=2)
    axes[0, k].set_title(f'EOF{k+1} ({var_exp[k]:.1f}%)')
    axes[0, k].set_xlabel('Space'); axes[0, k].grid(True, alpha=0.3)

# PC time series
for k in range(3):
    axes[1, k].plot(t, pcs[:, k], lw=0.8)
    axes[1, k].set_title(f'PC{k+1}')
    axes[1, k].set_xlabel('Time (months)'); axes[1, k].grid(True, alpha=0.3)

plt.suptitle('EOF Analysis of SST-like Data', fontsize=14)
plt.tight_layout()
plt.savefig('eof_analysis.png', dpi=150)
plt.show()

# === 2. Wavelet analysis of ENSO-like signal ===
fig, axes = plt.subplots(2, 1, figsize=(14, 8))

# Time series
enso = np.sin(2*np.pi*t/36) * (1 + 0.5*np.sin(2*np.pi*t/120)) + \
       0.3*np.random.randn(nt)
axes[0].plot(t, enso, 'b-', lw=1)
axes[0].set_xlabel('Time (months)'); axes[0].set_ylabel('ENSO Index')
axes[0].set_title('Synthetic ENSO Signal'); axes[0].grid(True, alpha=0.3)

# Wavelet power spectrum using manual Morlet CWT
widths = np.arange(2, 80)
def morlet_cwt(signal_data, scales, w=6):
    """Manual CWT with Morlet wavelet."""
    N = len(signal_data)
    result = np.zeros((len(scales), N), dtype=complex)
    for i, s in enumerate(scales):
        # Morlet wavelet in frequency domain
        freqs = np.fft.fftfreq(N)
        wavelet_ft = np.exp(-0.5 * ((2*np.pi*freqs*s - w)**2)) * (np.pi**(-0.25))
        signal_ft = np.fft.fft(signal_data)
        result[i] = np.fft.ifft(signal_ft * wavelet_ft)
    return result
cwt_matrix = morlet_cwt(enso, widths, w=6)
power = np.abs(cwt_matrix)**2
periods = widths * 2  # approximate period in months

axes[1].contourf(t, periods, power, levels=20, cmap='jet')
axes[1].set_xlabel('Time (months)'); axes[1].set_ylabel('Period (months)')
axes[1].set_title('Wavelet Power Spectrum')
axes[1].set_ylim(0, 100)
axes[1].axhline(36, color='white', ls='--', alpha=0.7, label='ENSO period')
axes[1].legend()

plt.tight_layout()
plt.savefig('wavelet_enso.png', dpi=150)
plt.show()

# === 3. Optimal interpolation of scattered data ===
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Generate scattered observations
n_obs = 40
x_obs = np.random.uniform(0, 100, n_obs)
y_obs = np.random.uniform(0, 100, n_obs)
truth = lambda x, y: 2*np.sin(0.05*x)*np.cos(0.04*y) + 0.5
z_obs = truth(x_obs, y_obs) + np.random.normal(0, 0.2, n_obs)

# OI parameters
L = 20.0          # decorrelation length (km)
eps2 = 0.1        # noise-to-signal ratio

# Grid
xg = np.linspace(0, 100, 50)
yg = np.linspace(0, 100, 50)
Xg, Yg = np.meshgrid(xg, yg)
Z_oi = np.zeros_like(Xg)

# Build observation covariance matrix
C_obs = np.zeros((n_obs, n_obs))
for ii in range(n_obs):
    for jj in range(n_obs):
        d = np.sqrt((x_obs[ii]-x_obs[jj])**2 + (y_obs[ii]-y_obs[jj])**2)
        C_obs[ii, jj] = np.exp(-d**2 / (2*L**2))
C_obs += eps2 * np.eye(n_obs)
C_inv = np.linalg.inv(C_obs)

# Interpolate to each grid point
for ix in range(len(xg)):
    for iy in range(len(yg)):
        c_vec = np.array([
            np.exp(-((xg[ix]-x_obs[k])**2 + (yg[iy]-y_obs[k])**2) / (2*L**2))
            for k in range(n_obs)])
        Z_oi[iy, ix] = np.mean(z_obs) + c_vec @ C_inv @ (z_obs - np.mean(z_obs))

cf = axes[0].contourf(Xg, Yg, Z_oi, levels=15, cmap='RdBu_r')
axes[0].scatter(x_obs, y_obs, c=z_obs, edgecolors='k', s=30, cmap='RdBu_r')
plt.colorbar(cf, ax=axes[0])
axes[0].set_title('Optimal Interpolation Result')

# Truth
Z_true = truth(Xg, Yg)
cf2 = axes[1].contourf(Xg, Yg, Z_true, levels=15, cmap='RdBu_r')
plt.colorbar(cf2, ax=axes[1])
axes[1].set_title('True Field')

plt.tight_layout()
plt.savefig('optimal_interp.png', dpi=150)
plt.show()

print(f"EOF1: {var_exp[0]:.1f}%  EOF2: {var_exp[1]:.1f}%  EOF3: {var_exp[2]:.1f}%")
print(f"OI RMSE: {np.sqrt(np.mean((Z_oi - Z_true)**2)):.3f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran: Optimal Interpolation for Irregularly-Spaced Ocean Data

This program implements the Bretherton, Davis, and Fandry (1976) optimal interpolation scheme for mapping scattered ocean observations (e.g., ship tracks, Argo floats) onto a regular grid with uncertainty estimates.

Fortran: Optimal Interpolation for Irregularly-Spaced Ocean Data

Fortran

Optimal interpolation (objective analysis) for ocean data

program.f90148 lines

program optimal_interpolation
  ! Optimal interpolation (objective analysis) for ocean data
  ! Maps irregularly-spaced observations to regular grid
  implicit none

integer, parameter :: n_obs = 50       ! number of observations
  integer, parameter :: nx = 30, ny = 30 ! analysis grid
  real(8) :: x_obs(n_obs), y_obs(n_obs), z_obs(n_obs)
  real(8) :: x_grid(nx), y_grid(ny)
  real(8) :: z_analysis(nx, ny), z_error(nx, ny)
  real(8) :: C_obs(n_obs, n_obs), C_inv(n_obs, n_obs)
  real(8) :: c_vec(n_obs), z_anom(n_obs)
  real(8) :: L_decorr, eps2, sigma2
  real(8) :: d, z_mean, err_var
  integer :: i, j, k, m
  integer :: info

! Parameters
  L_decorr = 15.0d0     ! decorrelation length (grid units)
  eps2 = 0.05d0         ! noise-to-signal variance ratio
  sigma2 = 1.0d0        ! signal variance

! --- Generate synthetic scattered observations ---
  ! (In practice, read from file)
  z_mean = 0.0d0
  do i = 1, n_obs
    x_obs(i) = 100.0d0 * mod(dble(i * 37), 97.0d0) / 97.0d0
    y_obs(i) = 100.0d0 * mod(dble(i * 53), 89.0d0) / 89.0d0
    z_obs(i) = 2.0d0 * sin(0.05d0 * x_obs(i)) * &
               cos(0.04d0 * y_obs(i)) + 0.3d0 * sin(0.1d0 * x_obs(i))
    z_mean = z_mean + z_obs(i)
  end do
  z_mean = z_mean / dble(n_obs)

! Anomalies
  do i = 1, n_obs
    z_anom(i) = z_obs(i) - z_mean
  end do

! --- Build observation covariance matrix ---
  do i = 1, n_obs
    do j = 1, n_obs
      d = sqrt((x_obs(i) - x_obs(j))**2 + (y_obs(i) - y_obs(j))**2)
      C_obs(i, j) = sigma2 * exp(-d**2 / (2.0d0 * L_decorr**2))
    end do
    C_obs(i, i) = C_obs(i, i) + eps2 * sigma2   ! add noise
  end do

! --- Invert covariance matrix (Gauss-Jordan) ---
  C_inv = C_obs
  call gauss_jordan_inv(C_inv, n_obs, info)
  if (info /= 0) then
    write(*,'(A)') 'ERROR: Matrix inversion failed'
    stop
  end if

! --- Set up analysis grid ---
  do i = 1, nx
    x_grid(i) = 100.0d0 * dble(i - 1) / dble(nx - 1)
  end do
  do j = 1, ny
    y_grid(j) = 100.0d0 * dble(j - 1) / dble(ny - 1)
  end do

! --- Perform optimal interpolation ---
  write(*,'(A)') '=== Optimal Interpolation ==='
  write(*,'(A,I4,A,I4)') 'Grid: ', nx, ' x ', ny
  write(*,'(A,I4)') 'Observations: ', n_obs
  write(*,'(A,F6.1)') 'Decorrelation length: ', L_decorr

do i = 1, nx
    do j = 1, ny
      ! Covariance between grid point and observations
      do k = 1, n_obs
        d = sqrt((x_grid(i) - x_obs(k))**2 + (y_grid(j) - y_obs(k))**2)
        c_vec(k) = sigma2 * exp(-d**2 / (2.0d0 * L_decorr**2))
      end do

! Analysis value: phi_a = mean + c^T * C_inv * (obs - mean)
      z_analysis(i, j) = z_mean
      do k = 1, n_obs
        do m = 1, n_obs
          z_analysis(i, j) = z_analysis(i, j) + &
            c_vec(k) * C_inv(k, m) * z_anom(m)
        end do
      end do

! Analysis error variance
      err_var = sigma2
      do k = 1, n_obs
        do m = 1, n_obs
          err_var = err_var - c_vec(k) * C_inv(k, m) * c_vec(m)
        end do
      end do
      z_error(i, j) = sqrt(max(err_var, 0.0d0))
    end do
  end do

! --- Output sample results ---
  write(*,'(/,A)') 'Sample analysis values (every 5th point):'
  write(*,'(A)') '  x      y      Analysis   Error'
  do i = 1, nx, 5
    do j = 1, ny, 10
      write(*,'(2F7.1, F10.4, F9.4)') &
        x_grid(i), y_grid(j), z_analysis(i,j), z_error(i,j)
    end do
  end do

contains
  subroutine gauss_jordan_inv(A, n, info)
    integer, intent(in) :: n
    real(8), intent(inout) :: A(n,n)
    integer, intent(out) :: info
    real(8) :: Aug(n, 2*n), pivot, factor
    integer :: i, j, k, prow
    info = 0
    ! Build augmented matrix [A | I]
    Aug(:, 1:n) = A
    Aug(:, n+1:2*n) = 0.0d0
    do i = 1, n
      Aug(i, n+i) = 1.0d0
    end do
    ! Forward elimination with partial pivoting
    do k = 1, n
      prow = k
      do i = k+1, n
        if (abs(Aug(i,k)) > abs(Aug(prow,k))) prow = i
      end do
      if (prow /= k) then
        do j = 1, 2*n
          factor = Aug(k,j); Aug(k,j) = Aug(prow,j); Aug(prow,j) = factor
        end do
      end if
      pivot = Aug(k,k)
      if (abs(pivot) < 1.0d-15) then
        info = k; return
      end if
      Aug(k,:) = Aug(k,:) / pivot
      do i = 1, n
        if (i /= k) then
          factor = Aug(i,k)
          Aug(i,:) = Aug(i,:) - factor * Aug(k,:)
        end if
      end do
    end do
    A = Aug(:, n+1:2*n)
  end subroutine
end program optimal_interpolation

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← AUVs Future →