clear
% /dods/matlib/rnt/util/lprf.m
mfig;
p=[1.0267      3.6669      6.4404      3.6536];

[tmp1,tmp2]=rout_npgo;
npgo.time=tmp1;
npgo.index=tmp2;
%================================================================
%% Load Background data for SPEEDY
%================================================================
addpath /sdd/PALEO/speedy2/matlib

lat=[ -87.16  -83.47  -79.78  -76.07  -72.36  -68.65  -64.94 ...
-61.23  -57.52  -53.81  -50.10  -46.39  -42.68  -38.97  -35.26  -31.54  ...
-27.83  -24.12  -20.41  -16.70  -12.99   -9.28   -5.57   -1.86    1.86  ...
  5.57    9.28   12.99   16.70   20.41   24.12   27.83   31.54   35.26  ...
 38.97   42.68   46.39   50.10   53.81   57.52   61.23   64.94   68.65  ...
 72.36   76.07   79.78   83.47   87.16]';
lon=0:3.75:3.75*95;
[LON,LAT]=meshgrid(lon,lat); LON=LON'; LAT=LAT'; lon=LON-360; lat=LAT;


%================================================================
%% NPO classical computation
%================================================================
%lat_range=[20 70]; w_type=1; runmean=3;
%atm=Compute_NPO_AL(lat_range, w_type, runmean,'ncep');
%atms=Compute_NPO_AL(lat_range, w_type, runmean,'speedy');

%================================================================
% Load ENSEMBLE SPEEDY AGCM
%================================================================ 
% Select the region of interest
pac.xlim=[-326  -64.0553];
pac.ylim=[-50 65];
i=find(lon(:,1) >= pac.xlim(1) & lon(:,1) <= pac.xlim(2));
j=find(lat(1,:) >= pac.ylim(1) & lat(1,:) <= pac.ylim(2));

slpc.lon=lon(i,j);
slpc.lat=lat(i,j);
slpc.mask=lon(i,j)*0+1;

load /nas/edl/SPEEDY/NPO_ensemble/SLPA_ENSEMBLE_2008.mat
in=[1:14 16:17 19:26 28:30]; in=1:40;
slpc.ano=SLPA(i,j,:,in);

% change date for this ensemble
clear k; k=0;
for yr=1950:2008
for imon=1:12
  k=k+1;
  slpc.month(k)=imon;
  slpc.year(k)=yr;
  slpc.datenum(k)=datenum(yr,imon, 15);
end
end
slpc.time=slpc.datenum;


%================================================================
%% Compute NPO Southern Pole from SPEDDY ensemble
%================================================================

agcm.time=slpc.time;
agcm.npo=-sq(mean(slpc.ano(43:45,19,25:end,:),1));
agcm.al =-sq(mean(slpc.ano(43:45,23,25:end,:),1));

agcm.times=slpc.time(25:end);
agcm.npos=mean(agcm.npo,2);
agcm.als=mean(agcm.al,2);
[year,month]=dates_datenum(agcm.times);
agcm.years=year; agcm.months=month;

sstc.ano=SSTA(:,:,:,25:end);

for i=1:684
 agcm.std(i)=std(agcm.npo(i,:));
end

% Compute average SST pattern associated with NPOhi
in=find(sstn.month==2);
in2=find(agcm.months==2);
CORR=zeros(134,61,40); REGR=CORR;
for i=1:40
i
tic
 slpano=ndnanfilter(agcm.npo(:,i), 'rectwin', 3);
 o=rnc_corr2(sstn.year(in), sstn.ano3(:,:,in), agcm.years(in2), slpano(in2), 1);
 CORR(:,:,i)=o.corr;
 REGR(:,:,i)=o.r1;
toc
end

in=find(slpc.month==2);
in1=find(slpn.month==2);
in2=find(ba.month==2);
CORRcpi=zeros(69,31,40); REGRcpi=CORRcpi;
for i=1:40
i
tic
 tmp=rnt_filter3D(slpc.ano(:,:,:,i), slpc.mask,3,'rectwin');
 tmp(isnan(tmp))=0;
 o=rnc_corr2(slpc.year(in), tmp(:,:,in), ba.year(in2), ba.cpi(in2), 1);
 CORRcpi(:,:,i)=o.corr;
 REGRcpi(:,:,i)=o.r1;
toc
end
% obs
o=rnc_corr2(slpn.year(in), slpn.ano3(:,:,in), ba.year(in2), ba.cpi(in2), 1);

for i=1:40
  plot(agcm.times, detrend(agcm.npo(:,i)),'color',[0.5 0.5 0.5],'linewidth',0.5)
  hold on
end
datetick
 plot(agcm.times, detrend(agcm.npos),'color','k','linewidth',1.3)
set(gca,'ylim',[-4.5 4.5]); set(gca,'xlim', [npgo_hs.time(1) npgo_hs.time(end)]); 

%================================================================
%   Hindcast the NPGO
%================================================================   
mfig;
p=[2.0135      1.7468       4.467      7.4938];
set(gcf, 'paperposition',p);
npgo_hs=AR1_model(agcm.times, agcm.npos, 1/10); 
red_signi(npgo_hs.time, lowpass(npgo_hs.sig,1/(12*8)), npgo.time, lowpass(npgo.index,1/(12*8)), 30, 1500);
tmp2=npgo_hs.sig;

npgo_hs=AR1_model(agcm.times, agcm.npos, 1/5); 
red_signi(npgo_hs.time, lowpass(npgo_hs.sig,1/(12*8)), pdo.datenum(2:end-3), lowpass(pdo.index(2:end-3),1/(12*8)), 30, 1500);
nnpgo=interp1(npgo.time, npgo.index, agcm.times);

% compute the ensemble reconstruction of NPGO
for i=1:40
   npgo_hs=AR1_model(agcm.times, agcm.npo(:,i)-agcm.npos, 1/10);
   sig(:,i)=npgo_hs.sig;
   tmp1=corrcoef(agcm.npos, agcm.npo(:,i));
   corrnp(i)=tmp1(1,2);
   n=1/(12*8);
   tmp1=corrcoef(sig(:,i), tmp2);
   corrnpl(i)=tmp1(1,2);
%   tmp1=corrcoef(lowpass(sig(:,i),n), lowpass(tmp2,n));
%   corrnpgo(i)=tmp1(1,2);
   
end 
[N,bin]=histc(corrnpl,[0.1:0.1:1]);
bar([0.1:0.1:1], N/sum(N)*100, 'histc');
%/dods/matlib/rnt/util/lprf.m

red_signi(npgo_hs.time, mean(sig,2), npgo.time, npgo.index, 30, 1500);  
red_signi(npgo_hs.time, lowpass(mean(sig,2),1/(12*8)), npgo.time, lowpass(npgo.index,1/(12*8)), 30, 1500);

% Compute the avergae spectra
%/dods/matlib/obsolete/spectra/myspec.m

pp=reshape(agcm.npo, [684*40 1]);
pp=pp-repmat(agcm.npos, [40 1]);
[sigh2, freq , period, sigh2_AR1, sigh2_AR1_90, phase ] = myspec(pp,1/12,40);
%/dods/matlib/obsolete/spectra/plot_myspec.m   
plot_myspec (sigh2, sigh2_AR1, sigh2_AR1_90, freq)

clf
plot(npgo.time, -detrend(lowpass(npgo.index,1/(12*8))),'color','b','linewidth',1.3); hold on
plot(npgo_hs.time, detrend(lowpass(npgo_hs.sig,1/(12*8))),'color','k','linewidth',1.3)
datetick;
set(gca,'ylim',[-3 3]); set(gca,'xlim', [npgo_hs.time(1) npgo_hs.time(end)]); 

clf
plot(npgo.time, -detrend((npgo.index)),'color','b','linewidth',1.3); hold on
plot(npgo_hs.time, detrend((npgo_hs.sig)),'color','k','linewidth',1.3)
datetick;
set(gca,'ylim',[-3.5 3.5]); set(gca,'xlim', [npgo_hs.time(1) npgo_hs.time(end)]); 

clf
plot(agcm.times, agcm.npos,'color','k','linewidth',1.3); hold on
set(gca,'ylim',[-3.5 3.5]); set(gca,'xlim', [npgo_hs.time(1) npgo_hs.time(end)]); 

npgo_hn=AR1_model(n.npo.time, n.npo.slpi, 1/10);
red_signi(npgo_hn.time, npgo_hn.sig, npgo.time, npgo.index, 30, 1500);

npo.mean=rnc_corr2(time_ens, mean(slpc.ano(:,:,25:end,:),4) , time_ens, l,30);
npo_sst.mean=rnc_corr2(sstn.time, sstn.ano , time_ens, l,30);
npo.mean12=rnc_corr2(time_ens, mean(SLPA(:,:,25:end,:),4) , time_ens, lowpassa(l,12),30);
npo_sst.mean12=rnc_corr2(sstn.time, sstn.ano , time_ens, lowpassa(l,12),30);

clf;rnc_map(npo_sst.mean12.corr, sstc);
set(gca, 'xlim', [-270 -60], 'ylim', [-20 20]); 
caxis([-.5 .5]); colorbar 'h'

% compute correlation for each month
[yr1,m1]=dates_datenum(sstn.time);
[yr2,m2]=dates_datenum(time_ens);

for imon=1:12
imon
in1=find(m1 == imon); in2=find(m2 == imon);
seas(imon)=rnc_corr2(yr1(in1), sstc.ssta(:,:,in1) , yr2(in2), l(in2),1);
end


% cpw region 
i=38:48; j=12:15;
for imon=1:12
clf;rnc_map(seas(imon).corr, sstc);
set(gca, 'xlim', [-270 -60], 'ylim', [-13 13]); 
caxis([-.7 .7]); colorbar 'h'
title (imon);
cind(imon)=sq(mean(mean( seas(imon).corr(i,j),1),2));
pause
end

for imon=1:12
imon1=imon-2; 
lag=0
if imon1 <= 0
   imon1=12+imon1; lag=+1
end
in1=find(m1 == imon1); in2=find(m2 == imon);
seas_lag2(imon)=rnc_corr2(yr1(in1)+lag, sstc.ssta(:,:,in1) , yr2(in2), l(in2),1);
cind_lag2(imon)=sq(mean(mean( seas_lag2(imon).corr(i,j),1),2));
end


in1=find(m1 == 12); in2=find(m2 == 2);
tmp=rnc_corr2(yr1(in1)+1, sstc.ssta(:,:,in1) , yr2(in2), l(in2),1);


slpi.ind=[];
slpi.ssta=zeros(69,31,15876);
slpi.ind_time=[];
slpi.sst_time=[];
for i=1:27
i
itrange=[1:588]+588*(i-1);
 slpi.ind = [slpi.ind ; slpi_ens(:,i)];
 slpi.ind_time=[slpi.ind_time ; time_ens];
 slpi.sst_time=[slpi.sst_time ; sstc.datenum(25:end)];
 slpi.ssta(:,:,itrange)=sstc.ssta(:,:,25:end);
end
[yr1,m1]=dates_datenum(slpi.sst_time);
[yr2,m2]=dates_datenum(slpi.ind_time);

in1=find(m1 == 12); in2=find(m2 == 2);
yr=1:49*27;
tmp=rnc_corr2(yr+1, slpi.ssta(:,:,in1) , yr, slpi.ind(in2),1);

% cpw region 
i=38:48; j=12:15;
for imon=1:12
imon
imon1=imon-2; 
lag=0;
if imon1 <= 0
   imon1=12+imon1; lag=+1;
end
in1=find(m1 == imon1); in2=find(m2 == imon);
seas_lag2(imon)=rnc_corr2(yr+lag, slpi.ssta(:,:,in1) , yr, slpi.ind(in2),1);
cind_lag2(imon)=sq(mean(mean( seas_lag2(imon).corr(i,j),1),2));
end



for ie=1:27
i
tic
 npo.ens(ie)=rnc_corr2(time_ens, SLPA(:,:,25:end,ie), time_ens, slpi_ens(:,ie) ,30);
 npo_sst.ens(ie)=rnc_corr2(sstc.time, sstc.ssta, time_ens, slpi_ens(:,ie) ,30);

 ll=lowpassa(slpi_ens(:,ie),12);
 npo.ens12(ie)=rnc_corr2(time_ens, SLPA(:,:,25:end,ie), time_ens, ll ,30);
 npo_sst.ens12(ie)=rnc_corr2(sstc.time, sstc.ssta, time_ens, ll ,30);
toc
end

SLP_CORR=zeros(69,31);
SLP_R1=SLP_CORR; SST_CORR=SLP_CORR; SST_R1=SLP_CORR;
num=27;
for ie=1:num
  SLP_CORR=SLP_CORR + npo.ens(ie).corr/num;
  SLP_R1=SLP_R1 + npo.ens(ie).r1/num;  
  SST_CORR=SST_CORR + npo_sst.ens(ie).corr/num;
  SST_R1=SST_R1 + npo_sst.ens(ie).r1/num;
end



for i=1:27
 [c99, c95, c12, signi]=red_signi(time_ens, lowpassa(l,2), n.time, lowpassa(time_ens(:,i),2), 30, 1500);


% SLP/SST maps form regression with PC2 (NPGO/NPO) modesd
k=0;
clear CORR_SLP CORR_SST REGR_SST REGR_SLP
for imon=+12:-2:-12
imon
tic
k=k+1;
o=rnc_corr(sstc,'ano', slpc.time-imon/12*365, slpi,30);
CORR_SST(:,:,k)=o.corr;
REGR_SST(:,:,k)=o.r1;
o=rnc_corr(slpc,'ano', slpc.time-imon/12*365, slpi,30);
CORR_SLP(:,:,k)=o.corr;
REGR_SLP(:,:,k)=o.r1;
MON(k)=-imon;
YEAR(k)=-imon/12;
toc
end


MON=-12:2:12
for i=1:2:11
clf;rnc_map(CORR_SST(:,:,i), sstc);
title(MON(i)); caxis([-.3 .3]); colorbar off
str=['psst_',num2str(i),'.eps'];
pause(0.2);%lpr(str);
pause
end

MON=-12:2:12
for i=1:13
clf;rnc_map(CORR_SLP(:,:,i), slp);
title(MON(i)); caxis([-.5 .5]); colorbar off
str=['slp_',num2str(i),'.eps'];
pause;lpr(str);
end

% SLP/SST maps form regression with PC2 (NPGO/NPO) modesd
k=0;
clear  CORR_SST REGR_SST CORR_SSTn REGR_SSTn 
for imon=+12:-1:-12
imon
tic
k=k+1;
o=rnc_corr(sstp,'ano', atms.time-imon/12*365, -p,30);
CORR_SST(:,:,k)=o.corr;
REGR_SST(:,:,k)=o.r1;

o=rnc_corr(sstn,'ano', atm.time-imon/12*365, atm.npo,30);
CORR_SSTn(:,:,k)=o.corr;
REGR_SSTn(:,:,k)=o.r1;
MON(k)=-imon;
YEAR(k)=-imon/12;
toc
end

pac.xlim=[-326  -64.0553];
pac.ylim=[-40 65];
trange=[datenum(1950, 1,1) datenum(2000, 12, 31)];
% Load SLP
p=pac;
tmp=rnc_Extract_SPEEDY_SLP(p.xlim,p.ylim,trange,0); 
slp = rnc_forcdStats(tmp,'SLP');
p=sq(slp.ano(44,17,:));


% MOVIE
for k=1:25
figure(1); clf
rnc_map(CORR_SLP(:,:,k),slp);
caxis([-0.4 0.4]); gradsmap4; colorbar 'h';
title ([num2str(YEAR(k)), ' YEAR LAG']);
figure(2); clf
rnc_map(CORR_SST(:,:,k),sstp);
caxis([-0.3 0.3]); colorbar 'h';
title ([num2str(YEAR(k)), ' YEAR LAG']);

figure(3); clf
rnc_map(REGR_SLP(:,:,k),slp);
caxis([-1 1]); gradsmap4; colorbar 'h';
title ([num2str(YEAR(k)), ' YEAR LAG']);
figure(4); clf
rnc_map(REGR_SST(:,:,k),sstp);
caxis([-0.5 0.5]); colorbar 'h';
title ([num2str(YEAR(k)), ' YEAR LAG']);

k
MON(k)
pause
end


% MOVIE
for k=1:25
clf;
rnc_map(CORR_SST(:,:,k),sstp);
caxis([-0.5 0.5]); colorbar 'h';
title ([num2str(YEAR(k)), ' YEAR LAG']);


k
MON(k)
pause
end




xlim=[-250 -66];
ylim=[-20 20];
i=find (lon(:,1) > xlim(1) & lon(:,1) < xlim(2));
j=find (lat(1,:) > ylim(1) & lat(1,:) < ylim(2));