import xarray as xr
import numpy as np
import matplotlib.pyplot as plt

#ds = xr.open_dataset('../00/04008.atmos_daily.nc', decode_times=False)
ds = xr.open_dataset('04008.atmos_daily.nc', decode_times=False)
p_meridiani = ds.sel(lat=0.0, lon=358.0, method='nearest')

# 1. Convertiamo i giorni cumulativi in Ls
# Nel MGCM Ames, 'areo' spesso contiene i gradi totali percorsi.
# Usiamo il modulo 360 per riportarlo nel range di un anno.
days_cum = p_meridiani.areo.values[:, 0]
ls_normalized = days_cum % 360 

# 2. Estrazione componenti per lo stress
u = p_meridiani.ukd.values
v = p_meridiani.vkd.values
tau = p_meridiani.stress.values
wind_speed = np.sqrt(u**2 + v**2)
tau_y = np.where(wind_speed > 0, tau * (v / wind_speed), 0)

# 3. PLOT
plt.figure(figsize=(10, 5))
plt.plot(ls_normalized, tau_y, label='Stress Meridionale (Tau_y)')
plt.axhline(0, color='black', lw=1)
plt.title(f"Stress su Meridiani - Finestra temporale: Ls {ls_normalized.min():.1f} a {ls_normalized.max():.1f}")
plt.xlabel("Ls (0-360)")
plt.ylabel("Tau_y [Pa]")
plt.grid(True)
plt.savefig('stresssumeridiani03.png')

print(f"Il tuo file copre solo da Ls {ls_normalized.min():.2f} a {ls_normalized.max():.2f}")