% KOE Analysis

% Compute circulation statstics and response of KOE to SSHa NPGO.
Figure_KOE_circulation



zoo=koe_zooplantlon;
inp=find(zoo.sbww >0.5);
in=find(ofes.month== 5);
v1=mean(U(:,:,in(inp)),3);
v2=mean(V(:,:,in(inp)),3);
amp=sqrt(v1.^2 + v2.^2);

inm=find(zoo.sbww < -0.5);
in=find(ofes.month== 5);
v1=mean(U(:,:,in(inm)),3);
v2=mean(V(:,:,in(inm)),3);
amp=sqrt(v1.^2 + v2.^2);


%% Long forward integrations
files=rnt_getfilenames('.','LONG_');
ctl=rnt_timectl(files,'time')

for ntime=1:length(ctl.time)
    PSTORE=rnt_loadvar(ctl, ntime, 'P');
    clf; rnc_map(PSTORE, grd);
    caxis([0 110]); colorbar 'h'
    title (datestr(ctl.time(ntime)));
    pause(0.1)
end


%% Load the forward passive tracer data
filesTR=rnt_getfilenames('/nas/edl/OFES_UV_KOE/matlib/pstore_nc','SB');
ctl=rnt_timectl(filesTR, 'time', 'd');
[year,month,day]=dates_datenum(ctl.time);
ctl.year=year;
ctl.month=month;
ctl.day=day;
ctl_sb=ctl;

clear P
ctl=ctl_sb;
k=0; clear year
for yr=1966:1999
    in=find(ctl.month >5  & ctl.year ==yr);
    k=k+1;
    P(:,:,k)=mean(rnt_loadvar(ctl, in, 'P'),3);
    year(k)=yr;
end

%
%yr=1980:1999;
Pa=P*0;
yr=1966:1999;
in=find(yr <1980);
Pa(:,:,in)=P(:,:,in)-repmat(mean(P(:,:,in),3),[1 1 length(in)]);
Pa(:,:,in)=Pa(:,:,in)./(repmat(stdNaN(Pa(:,:,in),3),[1 1 length(in)])+0.1);


in=find(yr >1979);
Pa(:,:,in)=P(:,:,in)-repmat(mean(P(:,:,in),3),[1 1 length(in)]);
Pa(:,:,in)=Pa(:,:,in)./(repmat(stdNaN(Pa(:,:,in),3),[1 1 length(in)])+0.1);

% Transition Region 36-40N 143-150W
xlim1=[143 150]-360; ylim1=[36 40];
% Oyashio Region 40-43N 143-150W
xlim2=[143 150]-360; ylim2=[40 44];

[Z,Zyear]=koe_zooplantlonM;
[eof,pc,vexp]=rnt_doEof(Pa, grd.mask);

C=Z*pc(:,1:10);

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

C=Z*Z';
[U,L]=eig(C);
Zpc=U'*Z;






% mfig(4); clf;
% rnc_map(mean(P,3), grd);
% caxis([-50 50]); gradsmap4; colorbar 'h'

[eof,pc,vexp]=rnt_doEof(Pa, grd.mask);
K=10;
ii=find(zoo.year>1579);
y=zoo.sbww(ii);
E=pc(:,1:K);
x=inv(E'*E)*(E'*y);
yhat=E*x;
figure(4)
red_signi(zoo.year, zoo.sbww, zoo.year(ii), yhat, 1, 1500, 1, 'r','b',0);
pmap=eof(:,:,1)*0;
for k=1:K
    pmap=pmap+eof(:,:,k)*x(k);
end
mfig(1); clf;
rnc_map(pmap, grd);

% TR box
xlim=xlim1; ylim=ylim1;
x=[xlim(1) xlim(1) xlim(2) xlim(2) xlim(1)];
y=[ylim(1) ylim(2) ylim(2) ylim(1) ylim(1)];
plot(x,y,'color','k', 'linewidth', 2);
% OY box
xlim=xlim2; ylim=ylim2;
x=[xlim(1) xlim(1) xlim(2) xlim(2) xlim(1)];
y=[ylim(1) ylim(2) ylim(2) ylim(1) ylim(1)];
plot(x,y,'color','k', 'linewidth', 2);

year=1966:1999;
o=rnc_corr3(year, Pa, grd.mask, zoo.year, ptr);
clf;rnc_map(o.corr, grd);
caxis([-50 50]); gradsmap4; colorbar 'h'

in=find(zoo.sbww >0.5);
mfig(4);clf;rnc_map(mean(Pa(:,:,in),3), grd);
caxis([-5 5]); gradsmap4; colorbar 'h'

in=find(zoo.sbww <-0.5);
mfig(5);clf;rnc_map(mean(Pa(:,:,in),3), grd);
caxis([-5 5]); gradsmap4; colorbar 'h'

[i,j]=rgrd_FindIJ(grd.lon, grd.lat, xlim1, ylim1);
ptr=sq(meanNaN(meanNaN(Pa(i,j,:),1),2));
zoo=koe_zooplantlon;
red_signi(zoo.year, zoo.sbww, zoo.year, ptr, 1, 1500, 1, 'r','b',0);


%% Load the inverse passive tracer data
filesTR=rnt_getfilenames('/nas/edl/OFES_UV_KOE/matlib/pstore_nc','TR');
ctl=rnt_timectl(filesTR, 'time', 'd');
[year,month,day]=dates_datenum(ctl.time);
ctl.year=year;
ctl.month=month;
ctl.day=day;
ctl_tr=ctl;

filesOY=rnt_getfilenames('/nas/edl/OFES_UV_KOE/matlib/pstore_nc','OY');
ctl=rnt_timectl(filesOY, 'time', 'd');
[year,month,day]=dates_datenum(ctl.time);
ctl.year=year;
ctl.month=month;
ctl.day=day;
ctl_oy=ctl;

% Transition Region 36-40N 143-150W
xlim1=[143 150]-360; ylim1=[36 40];

% Oyashio Region 40-43N 143-150W
xlim2=[143 150]-360; ylim2=[40 44];

zoo=koe_zooplantlon;
red_signi(zoo.year, zoo.sbww, zoo.year, zoo.sbTmodS, 1, 1500, 1, 'r','b',0);

clear P
ctl=ctl_tr;
k=0; clear year
for yr=1966:1999
    in=find(ctl.month == 4 & ctl.year ==yr);
    k=k+1;
    P(:,:,k)=mean(rnt_loadvar(ctl, in, 'P'),3);
    year(k)=yr;
end

Pa=P*0;
yr=1966:1999;
in=find(yr <1980);
Pa(:,:,in)=P(:,:,in)-repmat(mean(P(:,:,in),3),[1 1 length(in)]);
%Pa(:,:,in)=Pa(:,:,in)./repmat(stdNaN(Pa(:,:,in),3),[1 1 length(in)]);

in=find(yr >1979);
Pa(:,:,in)=P(:,:,in)-repmat(mean(P(:,:,in),3),[1 1 length(in)]);
%Pa(:,:,in)=Pa(:,:,in)./repmat(stdNaN(Pa(:,:,in),3),[1 1 length(in)]);


mfig; clf;
rnc_map(mean(P,3), grd);

Pa=P-repmat(mean(P,3),[1 1 length(year)]);
[eof,pc,vexp]=rnt_doEof(Pa, grd.mask);

K=10;
y=zoo.sbww;
E=pc(:,1:K);
x=inv(E'*E)*(E'*y);
yhat=E*x;
figure(4)
red_signi(zoo.year, zoo.sbww, zoo.year, yhat, 1, 1500, 1, 'r','b',0);
red_signi(zoo.year, zoo.sbww, zoo.year, p, 1, 1500, 1, 'r','b',0);

pmap=eof(:,:,1)*0;
for k=1:K
    pmap=pmap+eof(:,:,k)*x(k);
end
mfig(1); clf;
rnc_map(pmap, grd);
mfig; clf;
rnc_map(eof(:,:,2), grd);


for k=1:10
    
    pause
end





% TR box
xlim=xlim1; ylim=ylim1;
x=[xlim(1) xlim(1) xlim(2) xlim(2) xlim(1)];
y=[ylim(1) ylim(2) ylim(2) ylim(1) ylim(1)];
plot(x,y,'color','k', 'linewidth', 2);
% OY box
xlim=xlim2; ylim=ylim2;
x=[xlim(1) xlim(1) xlim(2) xlim(2) xlim(1)];
y=[ylim(1) ylim(2) ylim(2) ylim(1) ylim(1)];
plot(x,y,'color','k', 'linewidth', 2);


in=find(zoo.sbsc > 0.8);
mfig; clf;
rnc_map(mean(P(:,:,in),3)-mean(P,3), grd);
caxis([-10 10]); colorbar 'h' ; gradsmap4

in=find(zoo.sbsc < -0.8);
mfig; clf;
rnc_map(mean(P(:,:,in),3)-mean(P,3), grd);
caxis([-10 10]); colorbar 'h' ; gradsmap4



Pa=P-repmat(mean(P,3),[1 1 length(year)]);
o=rnc_corr3(year, Pa, grd.mask, zoo.year, zoo.sbww);

mfig(1); clf;
rnc_map(o.corr, grd);




%% KOE veloctiy

k=0; clear year
for yr=1966:2000
    in=find(ctl.month == 4 & ctl.year ==yr);
    k=k+1;
    P(:,:,k)=mean(rnt_loadvar(ctl, in, 'v'),3);
    year(k)=yr;
end