% ----------------------------------------------------------------------- %
%
% Extended Data Figure 2 of Archer et al. (2025)
%
% Matlab code shared on Zenodo.
%
% Paper:
% "Wide-swath satellite altimetry unveils global submesoscale ocean
% dynamics" by M. Archer, J. Wang, P. Klein, G. Dibarboure, and L-L Fu.
% Nature 2025.
%
% Requirements:
% --> Matlab: https://www.mathworks.com/products/matlab.html
% --> M_Map: www.eoas.ubc.ca/~rich/map.html
%
% Data links:
% --> SWOT L3: https://doi.org/10.24400/527896/a01-2023.018
%
% ----------------------------------------------------------------------- %
clear;clc

%% Load SWOT

% Path to Files
files = dir('/path/to/SWOT/L3/cycle13/all_passes/');
lenfiles = length(files)

ln = 50; ls = 20; le = 175; lw = 120; % Kuroshio lon/lat limits (north,south,east,west)

% >>>>>>>>>>>>>>>>> Load in data from folder

count=0;
for i = 1:lenfiles
    
    disp([num2str(i) ' out of ' num2str(lenfiles)])
    
    fname = [files(i).folder '/' files(i).name]; 
    
    % Load in Coordinates
    lon1 = ncread(fname,'longitude');lon1(lon1>180) = lon1(lon1>180)-360;
    lat1 = ncread(fname,'latitude');
    
    % Conditional
    cswath = 35; % get center of SWOT
    lonc = lon1(cswath,:); latc = lat1(cswath,:);
    indc = find(lonc >= lw & lonc <= le & latc >= ls & latc <= ln);    
    
    if length(indc) > 1
        
        count = count + 1;
                
        % Subset
        name{count} = files(i).name;
        time{count} = ncread(fname,'time',indc(1),length(indc))/60/60/24 + datenum(2000,1,1);
        timem(count) = nanmean(time{count});
        lon{count} = lon1(:,indc);
        lat{count} = lat1(:,indc);
        %        
        karin = ncread(fname,'ssha_unedited',[1 indc(1)],[69 length(indc)]);
        flag = ncread(fname,'quality_flag',[1 indc(1)],[69 length(indc)]);
        karin(flag~=0)=NaN; clear flag % applying strictest QC from L3 – but can choose lesser QC
        %
        ssha{count} = karin; clear karin

    end    
    
end
datestr(timem)

%% 5 eddies

% Index for the 5 eddies:
i = 28; alonglength = 945:970; alk = 35:69; % Pass 155
i = 37; alonglength = 910:940; alk = 1:35;  % Pass 196
i = 79; alonglength = 850:900; alk = 1:69;  % Pass 446
i = 82; alonglength = 790:820; alk = 1:35;  % Pass 461
i = 85; alonglength = 880:910; alk = 35:69; % Pass 474

% Load in swath
sshf = ssha{i};
lona = lon{i};
lata = lat{i};

% Cut swaths to eddy location
latc = lata(alk,alonglength);
lonc = lona(alk,alonglength);
eta = sshf(alk,alonglength); 

% Smooth field
sp=1;
etas = smooth2a(eta,sp,sp);

% Set f and g
f = gsw_f(latc);
g = 9.81*ones(size(f));

% Distances
[E,~] = lonlat2km(lonc(1:end-1,:),latc(1:end-1,:),lonc(2:end,:),latc(2:end,:));
[~,N] = lonlat2km(lonc(:,1:end-1),latc(:,1:end-1),lonc(:,2:end),latc(:,2:end));
DZ = cat(1,zeros(1,size(lonc,2)),cumsum(E,1)); % distance-zonal
DM = cat(2,zeros(1,size(lonc,1))',cumsum(N,2)); % distance-meridional
%
dx = vdiff(DZ,1).*1000; % km to m
dy = vdiff(DM,2).*1000; % km to m

% >>>>>>>> Get d(eta) in x and y
detax = vdiff(etas,1);
detay = vdiff(etas,2);

% >>>>>>>> Geostrophic velocity
vga = (g./f).*(detax./dx);
uga = -(g./f).*(detay./dy);
% Cut eta to same size
eta(isnan(uga))=NaN;eta(isnan(vga))=NaN;

% Plot figure
figure
pcolor(lonc,latc,eta),shading flat,hold on
quiver(lonc,latc,uga,vga,1,'k'),hold on
%
set(gcf,'color','w'),set(gca,'Layer','Top'), box on,
Y=get(gca,'ylim');set(gca,'dataaspectratio',[1 cos((Y(2)+Y(1))/2/180*pi) 1])