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/


  name = 'ALL'; %'STRESS' %'HEAT'
[mocd,mask1,lon1,lat1,yfrac1] = load_MOC(name);
%[sstd,mask,lon2,lat2,yfrac]   = load_SLA;
%[sstd,mask,lon2,lat2,yfrac] = load_PSL;
%[sstd,mask,lon2,lat2,yfrac] = load_coupled_SLA;
[sstd,mask,lon2,lat2,yfrac] = load_SLA(name);
%[sstd,mask,lon2,lat2,yfrac] = load_coupled_PSL;

%SMOOTH WITH MOVING AVERAGE
mocd = smoothdata(mocd,'movmean',13);
sstd = smoothdata(sstd,'movmean',13);

%DELETE FIRST 5 years months
it = 2 * 12+1:length(yfrac);
mocd = mocd(it,:);
sstd = sstd(it,:);
yfrac= yfrac(it);

%nx1 = length(lon1);
%ny1 = length(lat1);
nx = length(lon2);
ny = length(lat2);
nt = length(yfrac);

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

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

[SVD3,SVD4,PC3,PC4,EXPVAR,Lambda] = calsvd(sstd,mocd,neofs);
%keyboard
SVD3 = mat2map(mask,SVD3);  %SLA
%SVD3 = reshape(SVD3,neofs,ny,nx);
SVD4 = mat2map(mask1,SVD4); %MOC

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 7 5]); 
set(gcf,'DefaultLineLineWidth',1.5)
wysiwyg;

subplot('position',[.12 .4 .5 .5])
m_proj('miller','long',[-100 40],'lat',[0 75])
m_pcolor(lon2,lat2,(squeeze(SVD3(mm,:,:))));
shading flat
caxis([-0.04 0.04])
%colormap(jet(17))
colormap(mapa2)
m_grid('box','fancy')
m_coast('patch',[.5 .5 .5]);
[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))

subplot('position',[.7 .4 .25 .5])
pcolor(lon1,lat1,squeeze(SVD4(mm,:,:))');%,[-.1:.01:-.01],'--k');
shading flat
ylim([0 75])

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

end

