%% Load in the SST data
%Get info out of netcdf
file = 'thetao_Omon_GFDL-ESM2M_piControl_r1i1p1_049601-050012.nc';
info = ncinfo(file);

%Extract SST, lat, lon data from netcdf
SST = ncread(file,'thetao',[1 1 1 1], [inf inf 1 inf]); %only want top level of data
LAT = ncread(file,'lat');
LON = ncread(file,'lon');
LON(LON<0) = LON(LON<0)+360;

SST = squeeze(SST);

%Convert to columns
SST = reshape(SST, [size(SST,1)*size(SST,2) size(SST,3)]);

%Remove Seasonal Cycle
for i=1:12
    SSTclim(:,i) = nanmean(SST(:,i:12:end),2);
end
SST = SST-repmat(SSTclim, [1 size(SST,2)/12]);

%Keep only SST in Equatorial Pacific: 130E-280E, 15S-15N
SST_eqP = SST;
SST_eqP(LAT(:)>15,:) = NaN; SST_eqP(LAT(:)<-15,:) = NaN;
SST_eqP(LON(:)<130,:) = NaN; SST_eqP(LON(:)>280,:) = NaN;

%remove rows with NaN
SST_nonNaN = SST_eqP(~any(isnan(SST_eqP),2),:);

%% Load in the Surface Stress Data
%Get info out of netcdf
file = 'tauuo_Omon_GFDL-ESM2M_piControl_r1i1p1_049601-050012.nc';
info = ncinfo(file);

%Extract SST, lat, lon data from netcdf
TAU = ncread(file,'tauuo',[1 1 1], [inf inf inf]); 

TAU = squeeze(TAU);

%Convert to columns
TAU = reshape(TAU, [size(TAU,1)*size(TAU,2) size(TAU,3)]);

%Remove Seasonal Cycle
for i=1:12
    TAUclim(:,i) = nanmean(TAU(:,i:12:end),2);
end
TAU = TAU-repmat(TAUclim, [1 size(TAU,2)/12]);

%Keep only TAU in Equatorial Pacific: 130E-280E, 15S-15N
TAU_eqP = TAU;
TAU_eqP(LAT(:)>15,:) = NaN; TAU_eqP(LAT(:)<-15,:) = NaN;
TAU_eqP(LON(:)<130,:) = NaN; TAU_eqP(LON(:)>280,:) = NaN;
TAU_nonNaN = TAU_eqP(~any(isnan(TAU_eqP),2),:);

%% Calculate SVD of covariance Matrix
N = size(SST_nonNaN,2);
C = SST_nonNaN*TAU_nonNaN'/N;

[U,S,V] = svd(C);

%% Plot the results
%reshape the leading pattern of SST and tau to fit the lat/lon matrix
U1 = nan(size(SST,1),1);
U1(~any(isnan(SST_eqP),2))=U(:,1);
U1 = reshape(U1,size(LAT));

V1 = nan(size(TAU,1),1);
V1(~any(isnan(TAU_eqP),2))=V(:,1);
V1 = reshape(V1,size(LAT));

%Plot the Temperature and Tau patterns
load coast;
figure;
subplot(2,1,1);
axesm('miller','MapLonLimit',[110 300], 'MapLatLimit',[-20 20]);
pcolorm(double(LAT),double(LON),U1);
plotm(lat,long,'k');
tightmap;
colormap(b2r_noWhite(-.03,.03));
title('Temp Pattern');

subplot(2,1,2);
axesm('miller','MapLonLimit',[110 300], 'MapLatLimit',[-20 20]);
pcolorm(double(LAT),double(LON),V1);
plotm(lat,long,'k');
tightmap;
colormap(b2r_noWhite(-.03,.03));
title('Tau Pattern');

f = S(1,1)^2/(sum(sum(S).^2))