addpath /WSmounts/phodnet/share/mgoes/matlab2/Time_transform/       
addpath /data/mgoes/matlab/timeplt
addpath /home/mgoes/matlab/EOF_rotated
addpath /data/mgoes/matlab/m_map/
addpath /data/mgoes/matlab/mystuff
addpath /data/mgoes/matlab/pcatool
addpath /data/mgoes/matlab/Seawater
addpath /data/mgoes/matlab/drtoolbox/

num = 1;

    filename1 = 'PSL_coupled_passive.nc';
    title1 = 'EOF1 (SLP)'
    filename2 = 'PSL_coupled_passive.nc';
    title2 = 'EOF1 (SLA)'

[apresd,mask1,lon1,lat1,yfrac] = load_PSL;
apresd = smoothdata(apresd,'movmean',9);
%apresd = medfilt1(apresd,11,[],1);
[sstd,mask,lon2,lat2,yfrac]   = load_SLA;
sstd = smoothdata(sstd,'movmean',3);

%DELETE FIRST 48 months
apresd = apresd(49:end,:);
sstd = sstd(49:end,:);
yfrac= yfrac(49:end);

nx = length(lon1);
ny = length(lat1);
nt = length(yfrac);

%[nt,ny,nx] = size(ssta);

%Calculate the svds
disp('SVDs')
eplot = 1;
neofs = 3;   %number of eofs

[SVD3,SVD4,PC3,PC4,EXPVAR,Lambda] = calsvd(sstd,apresd,neofs);
SVD3 = mat2map(mask,SVD3);
%SVD4 = reshape(SVD4',neofs, ny,nx);
SVD4 = mat2map(mask1,SVD4);
%EXPVAR = L/sum(L);

load /data/mgoes/SACZ/mapa.mat
mapa2 = interp1(0:length(mapa)-1,mapa,linspace(0,length(mapa)-1,17));

disp('PLOTS')
for mm=1:neofs
figure(mm);clf
set(gcf,'PaperPosition',[1 1 5 7]); 
set(gcf,'DefaultLineLineWidth',1.5)
wysiwyg;

subplot('position',[.15 .4 .8 .5])
m_proj('miller','long',[-100 40],'lat',[0 75])
m_pcolor(lon2,lat2,(squeeze(SVD3(mm,:,:))));
shading flat
colormap(jet(17))
m_grid('box','fancy')
m_coast('patch',[.5 .5 .5]);
hold on
m_contour(lon1,lat1,squeeze(SVD4(mm,:,:)),[-.1:.01:-.01],'--k');
m_contour(lon1,lat1,squeeze(SVD4(mm,:,:)),[.01:.01:1],'k');
[C,LAGS] = xcorr(PC3(mm,:),PC4(mm,:),50,'coef');
[maxc,ic]=max(C);
title(sprintf('OHC SVD%i:  R= %3.2f,%s%2.0f%s',mm,maxc,'lag = ',LAGS(ic),' mo'))%num))
%title(sprintf('OHC SVD%i',num))

subplot('position',[.15 .2 .8 .2])
plot(yfrac,PC3(mm,:)/std(PC3(mm,:)))
hold on
plot(yfrac,PC4(mm,:)/std(PC4(mm,:)))
ylim([-5 5])
legend('OHC','SLP','orientation','horiz')
xlabel('year')
ylabel('PC')
axis tight
%xlim([1850 2020])
%print('-djpeg','-r300',sprintf('OHCUVIC_SVD%i%s',num,'.jpg'))

end

