Solver API#

pyCRO.solver.CRO_simulate(par, IC=[0, 0], N=2400, NE=100, NM='EH', dt=0.1, saveat=1.0, savemethod='sampling', freq='MS', use_cftime=False, start=None, EF=None, noise_custom=None, verbose=False)#

Run CRO ensemble simulation (same as RO_solver) and return results as an xarray Dataset.

This function wraps the low-level RO solver (RO_solver) and converts the output into a structured xarray.Dataset, including time and ensemble dimensions.

The CRO model simulates coupled SST (T) and thermocline (h) dynamics under stochastic forcing and optional external forcing.

Parameters:

  • par (dict) –

    CRO parameter dictionary (output of par_load), containing:

    • Linear and nonlinear deterministic parameters (R, F1, F2, epsilon, etc.)

    • Noise parameters (sigma_T, sigma_h, B)

    • Noise configuration (n_T, n_h, n_g, m_T, m_h)

  • IC (array-like of shape (2,), optional) –

    Initial conditions:

    • IC[0] : initial SST anomaly (T0)

    • IC[1] : initial thermocline anomaly (h0)

    Default is [0, 0].

  • N (float, optional) –

    Total simulation length in months.

    Default = 12 * 200 (200 years).

  • NE (int, optional) –

    Number of ensemble members.

    Default = 100.

  • NM ({"EM", "EH"}, optional) –

    Numerical integration scheme:

    • ”EM” : Euler–Maruyama method

    • ”EH” : Euler–Heun method (default)

  • dt (float, optional) – Time step (months). Default = 0.1.

  • saveat (float, optional) –

    Output saving interval (months). Must be a multiple of dt.

    Default = 1.0.

  • savemethod ({"sampling", "mean"}, optional) –

    Method for temporal aggregation:

    • ”sampling” : subsample at saveat interval (default)

    • ”mean” : block-average over each interval

  • freq (str, optional) – Time step frequency. Default is “MS” (month start).

  • use_cftime (bool, optional) – If True, return time axis as CFTime objects; otherwise use standard datetime objects. Default is False.

  • EF (dict, optional) –

    External forcing dictionary with keys:

    • ”E_T” : SST forcing time series

    • ”E_h” : thermocline forcing time series

    If None, zero forcing is used.

  • noise_custom (None, int, or ndarray, optional) –

    Custom stochastic forcing:

    • None : internally generated Gaussian noise

    • int : random seed for reproducibility

    • ndarray : pre-generated noise with shape (NT-1, 4, NE)

  • verbose (bool, optional) – If True, prints simulation configuration and progress.

Returns:

CRO simulation output with:

  • Nino34 : SST anomaly (T) [time × member]

  • WWV : thermocline depth anomaly (h) [time × member]

Coordinates: - time : monthly time index starting from year 0001 - member : ensemble member index

Return type:

xarray.Dataset

Notes

  • Time axis is generated using xarray.date_range(…, freq=”MS”) and uses CF-compatible time format.

  • Ensemble dimension corresponds to independent stochastic realizations.

  • Output is already aligned with climate-model-style diagnostics.

Examples

>>> par = pyCRO.par_load("ORAS5", "Linear-White-Additive")
>>> ds = pyCRO.CRO_simulate(par, NE=10, N=120)
>>> ds

Access variables:

>>> ds["T"]
>>> ds["h"]
pyCRO.solver.RO_analytic_solver(par, IC, N, NE, dt=0.1, saveat=1.0, savemethod='sampling', noise_custom=None)#

Analytical solution of the Recharge Oscillator (RO) model with deterministic and stochastic forcing. ONLY for linear RO model with white noise (annual mean parameters only)

Parameters:

  • par (dict) –

    Dictionary of model parameters, each as at least a one-element array:

    • 'R'float

      Damping parameter.

    • 'F1'float

      Feedback parameter relating thermocline depth to SST.

    • 'epsilon'float

      Thermocline damping parameter.

    • 'F2'float

      Feedback parameter relating SST to thermocline.

    • 'sigma_T'float

      Noise amplitude for SST.

    • 'sigma_h'float

      Noise amplitude for thermocline depth.

  • IC (tuple of float) –

    Initial condition (To, ho):

    • To : initial SST anomaly.

    • ho : initial thermocline depth anomaly.

  • N (float) – Total simulation length (time units).

  • NE (int) – Number of ensemble members.

  • dt (float, optional) – Numerical integration time step. Default is 0.1.

  • saveat (float, optional) – Output saving interval. Must be divisible by dt. Default is 1.0.

  • savemethod ({'sampling', 'mean'}, optional) –

    Method for saving results:

    • 'sampling' : take samples every saveat.

    • 'mean' : average values over each saveat interval.

  • noise_custom ({None, int, ndarray}, optional) –

    Specification of stochastic noise:

    • None (default): new Gaussian noise is generated for each ensemble.

    • int : seed for reproducible noise (same across ensembles).

    • ndarray of shape (NT-1, 4) : user-provided noise, repeated across ensembles.

    • ndarray of shape (NT-1, 4, NE) : user-provided noise for each ensemble.

Returns:

  • T_anal_out (ndarray of shape (N_out, NE)) – SST anomalies including deterministic and stochastic forcing, after applying the saving scheme. The final row is NaN-padded for dimensional consistency.

  • h_anal_out (ndarray of shape (N_out, NE)) – Thermocline anomalies including deterministic and stochastic forcing, after applying the saving scheme. The final row is NaN-padded.

  • noise_array (ndarray of shape (NT-1, 4, NE)) – Realizations of Gaussian noise used in the simulation.

Notes

The analytical solution consists of:

  • Deterministic part: exponential-sinusoidal functions of time.

  • Stochastic part: weighted sums (discrete convolutions) of noise with exponential and sinusoidal kernels.

Ensemble members are evaluated simultaneously using vectorized NumPy broadcasting, avoiding explicit Python loops.

Examples

>>> par = {'R':[0.5], 'F1':[1.2], 'epsilon':[0.3], 'F2':[0.8],
...        'sigma_T':[0.2], 'sigma_h':[0.1]}
>>> IC = (0.1, -0.2)
>>> T, h, noise = RO_solver_analytic(par, IC, N=50, NE=10, dt=0.1, saveat=1.0)
>>> T.shape, h.shape
((51, 10), (51, 10))
pyCRO.solver.RO_solver(par, IC, N, NE, NM='EH', dt=0.1, saveat=1.0, savemethod='sampling', EF=None, noise_custom=None, verbose=True)#

Numerical solution of the Recharge Oscillator (RO) model with stochastic forcing.

This function integrates the Recharge Oscillator (RO) system numerically using either the Euler–Maruyama (EM) or Euler–Heun (EH) scheme. It includes deterministic dynamics, stochastic noise (white or red), and optional external forcing. Ensemble simulations are supported, with results returned at user-specified output intervals.

Parameters:

  • par (dict) –

    Dictionary of model parameters with the following keys (as 1-element arrays):

    • 'R'float

      Damping parameter.

    • 'F1'float

      Feedback parameter relating thermocline depth to SST.

    • 'epsilon'float

      Thermocline damping parameter.

    • 'F2'float

      Feedback parameter relating SST to thermocline.

    • 'sigma_T'float

      Noise amplitude for SST.

    • 'sigma_h'float

      Noise amplitude for thermocline depth.

    • 'B'float

      Multiplicative noise coefficient (used when n_g = 0 or 1).

    • 'n_T'int

      Noise type in SST equation (1 = white noise, 0 = red noise).

    • 'm_T'ndarray

      Memory kernel for red noise in SST equation (ignored if n_T=1).

    • 'n_h'int

      Noise type in thermocline equation (1 = white noise, 0 = red noise).

    • 'm_h'ndarray

      Memory kernel for red noise in thermocline equation (ignored if n_h=1).

    • 'n_g'int

      Noise structure in SST equation: 0 = multiplicative, 1 = multiplicative + Heaviside, 2 = additive.

  • IC (tuple of float) –

    Initial condition (T0, h0):

    • T0 : initial SST anomaly.

    • h0 : initial thermocline depth anomaly.

  • N (float) – Total simulation length (time units).

  • NE (int) – Number of ensemble members.

  • NM ({'EM', 'EH'}, optional) – Numerical integration method. Default is 'EH' (Euler–Heun).

  • dt (float, optional) – Numerical integration time step. Default is 0.1.

  • saveat (float, optional) – Output saving interval. Must be divisible by dt. Default is 1.0.

  • savemethod ({'sampling', 'mean'}, optional) –

    Method for saving results:

    • 'sampling' : store values every saveat steps (default).

    • 'mean' : average values within each saveat interval.

  • EF (dict, optional) –

    External forcing with the following keys (as 1D arrays of length 5):

    • 'E_T' : SST forcing coefficients.

    • 'E_h' : Thermocline forcing coefficients.

    If None (default), no external forcing is applied.

  • noise_custom ({None, int, ndarray}, optional) –

    Specification of stochastic noise:

    • None (default): generate new Gaussian noise for each ensemble.

    • int : random seed for reproducible noise (same across ensembles).

    • ndarray of shape (NT-1, 4, NE) : user-provided noise realizations.

  • verbose (bool, optional) – If True (default), print detailed setup and progress messages.

Returns:

  • T_out (ndarray of shape (N_out, NE)) – SST anomalies from the numerical integration, after applying the saving scheme.

  • h_out (ndarray of shape (N_out, NE)) – Thermocline anomalies from the numerical integration, after applying the saving scheme.

  • noise_out (ndarray of shape (NT-1, 4, NE)) – Realizations of Gaussian noise used in the simulation.

Notes

  • White noise forcing corresponds to n_T = 1 or n_h = 1.

  • Red noise forcing (n_T = 0 or n_h = 0) requires a nonzero memory kernel m_T or m_h.

  • For the Euler–Maruyama method (NM = 'EM') with multiplicative noise, an Ito-to-Stratonovich correction is applied automatically.

Examples

>>> par = {
...     'R':[0.5], 'F1':[1.2], 'epsilon':[0.3], 'F2':[0.8],
...     'sigma_T':[0.2], 'sigma_h':[0.1], 'B':[0.05],
...     'n_T':[1], 'm_T':[0], 'n_h':[1], 'm_h':[0], 'n_g':[0]
... }
>>> IC = (0.1, -0.2)
>>> T, h, noise = RO_solver(par, IC, N=50, NE=5, dt=0.1, saveat=1.0)
>>> T.shape, h.shape
((51, 5), (51, 5))