load KOE_GRD.mat
% 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];



%% Compute UV used for tracer experiments OFES 3day
yrstr=1966;
yrend=1999;

k=0;
for yr=yrstr:yrend
    for imon=1:12
        k=k+1;
        ofes.time(k)=datenum(yr,imon, 15);
        ofes.month(k)=imon;
        ofes.year(k)=yr;
    end
end

[I,J]=size(grd.lon);
U=zeros(I,J,length(ofes.time));
V=U;
k=0;
for yr=yrstr:yrend
    for imon=1:12
        in=find(ctl.month==imon & ctl.year==yr);
        u=rnt_loadvar(ctl, in, 'u');
        v=rnt_loadvar(ctl, in, 'v');
        k=k+1;
        U(:,:,k)=mean(u,3);
        V(:,:,k)=mean(v,3);
        disp(datestr(ctl.time(in(1))));
    end
end

[Ua,Um]=RemoveSeas(U,ofes.month, 3);
[Va,Vm]=RemoveSeas(V,ofes.month, 3);


%% OFES SSHa analysis 1 month
%/ocean/drive/edl/NEPD/process_data/SalinityMODE/KOE
addpath /drive/edl/NEPD/OFES
xlim=[143 180]; ylim=[32 34];
% Load SSHa index
[ktime, kssha32]=OFES_KOE_ssha(xlim, ylim);

% Load KOE data in box
xlim(1)=min(grd.lon(:));  xlim(2)=max(grd.lon(:));
ylim(1)=min(grd.lat(:));  ylim(2)=max(grd.lat(:));
[ktime, kssha, ktrans, ktrans_lon, ksshap, koe]=OFES_KOE_ssha(xlim+360, ylim);


% KOE modes Taguchi
xlim=[143 180]; ylim=[30 45]
[ktime, kssha, ktrans, ktrans_lon, ksshap, koe2]=OFES_KOE_ssha(xlim, ylim);
ssha2=mean(koe2.ssha,1);

%ssha=dtrend2d(ssha);
[eof,pc,vexp]=rnt_doEof(ssha2, ssha2(:,:,1)*0+1);


%% Make figures
mfig;
p=[1.1603      3.2141      6.1749      4.5745];
set(gcf, 'paperposition',p);

% Mean circulation of KOE
amp=sqrt(mean(Um,3).^2 + mean(Vm,3).^2).*grd.mask;
v1=mean(Um,3).*grd.mask; v2=mean(Vm,3).*grd.mask;

clf
rnc_map(amp, grd); 
colormap (getpmap(7)); caxis([0 50]); colorbar 'h'
d=6; I=1:d:399; J=1:d:249;
lon=grd.lon(I,J).*grd.mask(I,J); lat=grd.lat(I,J).*grd.mask(I,J);
quiver(lon, lat, v1(I,J), v2(I,J), 4,'k')
hold on
contour(koe.lon, koe.lat, mean(koe.ssh,3), 'color', 'w', 'linewidth', 1.1);

clf
rnc_map(mean(koe.ssh,3), koe); 
hold on
%contour(koe.lon, koe.lat, mean(koe.ssh,3), 'color', [.8 .8 .8], 'linewidth', 1.1);
contour(koe.lon, koe.lat, mean(koe.ssh,3), 'color', 'w', 'linewidth', 1.5);

d=6; I=1:d:399; J=1:d:249;
lon=grd.lon(I,J).*grd.mask(I,J); lat=grd.lat(I,J).*grd.mask(I,J);
quiver(lon, lat, v1(I,J), v2(I,J), 4,'k')


% 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);

%SBOX
xlim=[-228.5265 -221.0984]; ylim=[28.7471   34.1533];
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','r', 'linewidth', 2);


ou=rnc_corr3(ofes.time, Ua, grd.mask, ktime, -kssha32);
ov=rnc_corr3(ofes.time, Va, grd.mask, ktime, -kssha32);
ossh=rnc_corr3(ktime, koe.ssha, koe.mask, ktime, -kssha32);

clf
rnc_map(ossh.r1, koe); 

caxis([-10 10]);  gradsmap4; colorbar 'h'
d=6; I=1:d:399; J=1:d:249;
lon=grd.lon(I,J).*grd.mask(I,J); lat=grd.lat(I,J).*grd.mask(I,J);
quiver(lon, lat, ou.r1(I,J), ov.r1(I,J), 2,'k')


zoo=koe_zooplantlon
p=interp1(datenum(zoo.year, 7,15), zoo.all(:,6), ktime);
p(isnan(p))=0;
ozoo=rnc_corr3(ktime, koe.ssha, koe.mask, ktime, p);

clf; rnc_map(ozoo.corr, koe);
caxis([-.4 .4]); gradsmap4; colorbar 'h'

%% Long forward integrations
files=rnt_getfilenames('.','LONG_KOE.nc');
ctl=rnt_timectl(files,'time')
PKOE=grd.lon*0;

k=0;
clear month year time
for yr=ctl.year(1):ctl.year(end)
    yr
    for imon=1:12
        in=find(ctl.year == yr & ctl.month ==imon);
        p=rnt_loadvar(ctl,in, 'P');
        k=k+1; PKOE(:,:,k)=mean(p,3);
        month(k)=imon;
        year(k)=yr;
        time(k)=datenum(yr,imon, 15);
    end
end
ir=find( ~isnan(sq(PKOE(1,1,:))) );
pkoe.tracer=PKOE(:,:,ir);
pkoe.year=year(ir);
pkoe.month=month(ir);
pkoe.time=time(ir);


files=rnt_getfilenames('.','LONG_SB.nc');
ctl=rnt_timectl(files,'time')
PKOE=grd.lon*0;

k=0;
clear month year time
for yr=ctl.year(1):ctl.year(end)
    yr
    for imon=1:12
        in=find(ctl.year == yr & ctl.month ==imon);
        p=rnt_loadvar(ctl,in, 'P');
        k=k+1; PKOE(:,:,k)=mean(p,3);
        month(k)=imon;
        year(k)=yr;
        time(k)=datenum(yr,imon, 15);
    end
end
ir=2:408-7;
psb.tracer=PKOE(:,:,ir);
psb.year=year(ir);
psb.month=month(ir);
psb.time=time(ir);

files=rnt_getfilenames('.','LONG_NB.nc');
ctl=rnt_timectl(files,'time')
PKOE=grd.lon*0;

k=0;
clear month year time
for yr=ctl.year(1):ctl.year(end)
    yr
    for imon=1:12
        in=find(ctl.year == yr & ctl.month ==imon);
        p=rnt_loadvar(ctl,in, 'P');
        k=k+1; PKOE(:,:,k)=mean(p,3);
        month(k)=imon;
        year(k)=yr;
        time(k)=datenum(yr,imon, 15);
    end
end
ir=find( ~isnan(sq(PKOE(1,1,:))) );
pnb.tracer=PKOE(:,:,ir);
pnb.year=year(ir);
pnb.month=month(ir);
pnb.time=time(ir);

%% Long Adjoint integrations
files=rnt_getfilenames('.','LONG_ADJTR.nc');
ctl=rnt_timectl(files,'time')
PKOE=grd.lon*0;

k=0;
clear month year time
for yr=ctl.year(end):ctl.year(1)
    yr
    for imon=1:12
        in=find(ctl.year == yr & ctl.month ==imon);
        p=rnt_loadvar(ctl,in, 'P');
        k=k+1; PKOE(:,:,k)=mean(p,3);
        month(k)=imon;
        year(k)=yr;
        time(k)=datenum(yr,imon, 15);
    end
end
ir=find( ~isnan(sq(PKOE(1,1,:))) );
ptr.tracer=PKOE(:,:,ir);
ptr.year=year(ir);
ptr.month=month(ir);
ptr.time=time(ir);

% Mean tracer
pp=ptr;
exp_name='TR';
clf;rnc_map(mean(pp.tracer,3), grd);
caxis([-100 100]); gradsmap4; colorbar 'h'
mask1=koe.mask;
i=find(koe.lon(:,1) < -220);
j=find(koe.lat(1,:) > 35);
mask1(i,j)=nan;
hold on
contour(koe.lon, koe.lat, mean(koe.ssh,3).*koe.mask.*mask1, 'color', 'k', 'linewidth', 1.5);
sbox
lpr ppmean.eps

pa=RemoveSeas(pp.tracer, pp.month, 3);
opa=rnc_corr3(pp.time, pa, grd.mask, ktime, -kssha32);
time=pp.time;


opas=rnc_corr3(time, pa, grd.mask, ktime, -nn(pc(:,1)));
mask1=abs(opas.r1);
mask1(mask1 < 0.1)=nan;
mask1(mask1>0)=1;
clf; rnc_map(opas.corr.*mask1, grd);
caxis([-.5 .5]); gradsmap4; colorbar 'h'
sbox
lpr ps1.eps

clf; rnc_map(opas.r1, grd);
caxis([-7 7]); gradsmap4; colorbar 'h'
sbox
lpr ps2.eps

opas=rnc_corr3(psb.time, pa, grd.mask, ktime, nn(pc(:,1)));
mask1=abs(opas.r1);
mask1(mask1 < 1)=nan;
mask1(mask1>0)=1;
clf; rnc_map(opas.corr.*mask1, grd);
caxis([-.5 .5]); gradsmap4; colorbar 'h'
sbox
lpr ps3.eps

clf; rnc_map(opas.r1, grd);
caxis([-7 7]); gradsmap4; colorbar 'h'
sbox
lpr ps4.eps


% hold on
% contour(koe.lon, koe.lat, mean(koe.ssh,3), 'color', 'k', 'linewidth', 1.1);
d=6; I=1:d:399; J=1:d:249;
lon=grd.lon(I,J).*grd.mask(I,J); lat=grd.lat(I,J).*grd.mask(I,J);
quiver(lon, lat, ou.r1(I,J), ov.r1(I,J), 2,'k')


% 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);

%SBOX
xlim=[-228.5265 -221.0984]; ylim=[28.7471   34.1533];
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','r', 'linewidth', 2);





for ntime=700: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