I’m trying to integrate following equation

The equation express the propagation of wave, where WE represents Ekman pumping and eta_e represents disturbance from eastern boundary. y is specified a certain latitude.

Though the algorithm is to integrate along time and longitude (x,t), the actual WE is integrated diagonally.
I’m wondering if there’s a better way to integrate to reduce computation time.

Following is part of my codes, I can provide more if need.

import numpy as np
import xarray as xr
from scipy.interpolate import interp1d
import pandas as pd
import matplotlib.pyplot as plt
import cftime
from datetime import timedelta
import datetime

def calc_ip_wsc(wsc, tshift, t_index):
    original_time = wsc.time[t_index].values
    shifted_time  = original_time - timedelta(days=tshift)
    if shifted_time < wsc.time[0].values or shifted_time > wsc.time[-1].values:
        return np.nan

    time_before   = wsc.time.sel(time=slice(None, shifted_time)).max().values
    time_after    = wsc.time.sel(time=slice(shifted_time, None)).min().values

    # Select data around the boundary times
    data_before = wsc.sel(time=time_before).values
    return data_before
    data_after = wsc.sel(time=time_after).values

    # Compute interpolation weights
    time_diff_before = (shifted_time - time_before).total_seconds()/86400
    time_diff_after  = (time_after - shifted_time).total_seconds()/86400
    total_diff       = (time_after - time_before).total_seconds()/86400

    if total_diff == 0:
        return data_before  # Edge case where the times are the same

    interpolated_value = (data_before * (time_diff_after / total_diff) +
                          data_after * (time_diff_before / total_diff))

    return interpolated_value

def calc_eta_e(eta_e, wsc):
    eta_e_x = np.full_like(wsc, np.nan)
    eta_e_x = eta_e[:,None] * np.exp(-(wsc.lon - wsc.lon[-1])[None,:] * 110e3 * epsilon/ CR)
    return eta_e_x

def calc_eta_m(wsc, eta_e_x):
    eta_m = np.full_like(wsc, np.nan)
    for i in range(len(time)):
        print(i,flush=True)
        for j in range(len(wsc.lon)):
            RHS1 = np.full(len(wsc.lon) - j, np.nan)
            for k in range(len(RHS1)):
                tshift =int((wsc.lon[j+k] - wsc.lon[j]).values * 110e3 / CR / 86400)
                WEp = calc_ip_wsc(wsc[:,j+k], tshift, i)
                if np.isnan(WEp):
                    RHS1[k] = 0
                else:
                    RHS1[k] = 1/CR * WEp * np.exp(-epsilon/CR * (wsc.lon[j+k] - wsc.lon[j]) * 110e3)

            RHS1_r = RHS1[::-1]
            lon_r  = wsc.lon[j:].values*110e3
            x_r    = lon_r[::-1]
            integral1 = np.trapz(RHS1, x=x_r)

            RHS2   = eta_e_x[i, j:]
            RHS2_r = RHS2[::-1]
            integral2 = np.trapz(RHS2, x=x_r)
            eta_m[i,j] = integral1 + integral2

    return eta_m

Thanks!

I think it would be better if I use cumtrapz for second integration.