%function plot_OHC(num)
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;

    filename = 'ETAT_NA.nc';
    filename = 'ETAT_FAFHEAT_NA.nc';
%    filename = 'ETAT_FAFSTRESS_NA.nc';
%    filename = 'ETAT_FAFWATER_NA.nc';
    title1 = 'EOF1 (SLA)'

latitude  = double(ncread(filename,'YT_OCEAN138_267'));
longitude = double(ncread(filename,'XT_OCEAN181_310'));
time      = double(ncread(filename,'TIME'));
yfrac = 1:length(time);

longitude(longitude>180) = longitude(longitude>180)-360;

disp('Load Data')
A_apress = ncread(filename,'ETA_T');

disp('NORTH ATLANTIC')
[nx,ny,nt] = size(A_apress);
%apress = permute(A_apress,[3 2 1]);      
%clear A_apress;


apresm = A_apress*nan;
disp('INTERP')
for jj=1:12
    ii = jj:12:nt;
    aux = mean(A_apress(:,:,ii),3);
    apresm(:,:,ii) = aux(:,:,ones(1,length(ii)));
   
end 
apressa = A_apress-apresm;
apressa = permute(apressa,[3 2 1]);

%PACIFIC AREA
apressa(:,1:30,1:30)=nan;
oo = std(apressa);
apressa(:,oo==0)=nan;

 
clear apresm A_apress

disp('Anomalies')
mask = ~isnan(squeeze(apressa(1,:,:)));
W = sw_area(longitude,latitude);
W=W'./max(W(:));
apressa=apressa.*W(ones(1,nt),:,:);
apresd = map2mat(mask,apressa)*100;

%PREWHITEN
if 1
for ii = 1:size(apresd,2)
      aux = apresd(:,ii);
      if ~isnan(aux)
      apresd(:,ii) = prewhiten(aux);
      end
end
end
%ssta  = map2mat(mask,ssta);
%[X,Y] = meshgrid(longitude,latitude);
%apresmi = apresmi.*sqrt(cos(shiftdim(Y,-1)*pi/180));
%stop

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

apresd(isnan(apresd))=0;

disp('EOFs')
eplot = 1;
%Calculate the eofs
neofs = 3;   %number of eofs
%EOFs = nan*ones(length(inan),neofs);
%SVD1 = EOFs;SVD2=EOFs;
%norma = 9.7235e+08; % (This is the sum of all eigenvalues)
[L, EOFs1, EC1 ] = EOF( apresd, 100, true);
EOFs1 = EOFs1(:,1:neofs);EC1=EC1(:,1:neofs);
EOFs1 = mat2map(mask,EOFs1');EOFs1 = permute(EOFs1,[2 3 1]);
EXPVAR = L/sum(L);

%[SVD1,SVD2,PC1,PC2,EXPVAR,Lambda] = calsvd(hilbert(prmsl),hilbert(apresd),neofs);
%[CEOFs,PC,EXPVAR] = calCeof(prmsl,neofs,2);
%rec1 = transpose(PC(1,:))*conj(CEOFs(1,:));
%rec1 = mat2map(mask,rec1);
%SVD1 = mat2map(mask,CEOFs);
%PC1 = real(PC(1,:));
%[SVD1,SVD2,PC1,PC2,EXPVAR,Lambda] = calsvd(hilbert(prmsl),(apresd),neofs);
%[SVD3,SVD4,PC3,PC4,EXPVAR,Lambda] = calsvd(prmsl,apresd,neofs);
%SVD1 = mat2map(mask,SVD1);
%SVD2 = reshape(SVD2',ny,nx,neofs);
%SVD3 = mat2map(mask,SVD3);
%SVD4 = reshape(SVD4',ny,nx,neofs);
%EOFs = reshape(EOFs,ny,nx,neofs);
%SVD1 = reshape(SVD1',ny,nx,neofs);

% oo1 = transpose(PC1([1],:))*conj(SVD1([1],:));
% oo2 = transpose(PC2([1],:))*conj(SVD2([1],:));
% rec1 = mat2map(mask,oo1);
% rec2 = mat2map(mask,oo2);
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',[.2 .45 .7 .5])
m_proj('miller','long',[-100 40],'lat',[0 75])
m_pcolor(longitude,latitude,EOFs1(:,:,mm));
shading flat
%colormap(jet(17))
colormap(mapa2)
m_grid('box','fancy')
m_coast('patch',[.5 .5 .5]);
title(sprintf('%s%i = %3.2f%s','EOF',mm,EXPVAR(mm)*100))
hold on
subplot('position',[.2 .33 .7 .15])
plot(yfrac,EC1(:,mm))
axis tight
ylim([-120 120])
set(gca,'ytick',-100:50:100)
xlabel('year')
ylabel('PC')
%print('-djpeg','-r300',sprintf('OHC_UVIC%i_EOF%i%s',num,mm,'.jpg'))

end

return
% mm=1
% figure;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(longitude,latitude,EOFs2(:,:,mm));
% shading flat
% colormap(jet(17))
% m_grid('box','fancy')
% m_coast('patch',[.5 .5 .5]);
% title(sprintf('SLP EOF%i',mm))%num))
% subplot('position',[.15 .2 .8 .2])
% plot(yfrac,EC2(:,mm))
% xlabel('year')
% ylabel('PC')
% xlim([1850 2020])
% %print('-djpeg','-r300',sprintf('OHCUVIC_EOF%i%s',num,'.jpg'))

figure;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(longitude,latitude,imag(squeeze(SVD1(1,:,:,1))));
shading flat
colormap(jet(17))
m_grid('box','fancy')
m_coast('patch',[.5 .5 .5]);
hold on
m_contour(longitude,latitude,imag(SVD2(:,:,1)),[-.1:.01:-.01],'--k');
m_contour(longitude,latitude,imag(SVD2(:,:,1)),[.01:.01:1],'k');
title(sprintf('OHC CSVDI%i',num))
subplot('position',[.15 .2 .8 .2])
plot(yfrac,imag(PC1(1,:)))
hold on
plot(yfrac,imag(PC2(1,:)))
legend('OHC','SLP')
xlabel('year')
ylabel('PC')
%xlim([1850 2020])
axis tight
%print('-djpeg','-r300',sprintf('OHCUVIC_SVD%i%s',num,'.jpg'))

figure;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(longitude,latitude,real(squeeze(SVD1(1,:,:,1))));
shading flat
colormap(jet(17))
m_grid('box','fancy')
m_coast('patch',[.5 .5 .5]);
hold on
m_contour(longitude,latitude,real(SVD2(:,:,1)),[-.1:.01:-.01],'--k');
m_contour(longitude,latitude,real(SVD2(:,:,1)),[.01:.01:1],'k');
title(sprintf('OHC CSVDR%i',num))
subplot('position',[.15 .2 .8 .2])
plot(yfrac,real(PC1(1,:)))
hold on
plot(yfrac,real(PC2(1,:)))
legend('OHC','SLP')
xlabel('year')
ylabel('PC')
%xlim([1850 2020])
axis tight
%print('-djpeg','-r300',sprintf('OHCUVIC_SVD%i%s',num,'.jpg'))

figure;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(longitude,latitude,(squeeze(SVD3(1,:,:,1))));
shading flat
colormap(jet(17))
m_grid('box','fancy')
m_coast('patch',[.5 .5 .5]);
hold on
m_contour(longitude,latitude,(SVD4(:,:,1)),[-.1:.01:-.01],'--k');
m_contour(longitude,latitude,(SVD4(:,:,1)),[.01:.01:1],'k');
[C,LAGS] = xcorr(PC3(1,:),PC4(1,:),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(1,:)))
hold on
plot(yfrac,(PC4(1,:)))
ylim([-10 10])
legend('OHC','SLP','orientation','horiz')
xlabel('year')
ylabel('PC')
axis tight
%xlim([1850 2020])
%print('-djpeg','-r300',sprintf('OHCUVIC_SVD%i%s',num,'.jpg'))
